Enhancement of the double Higgs production via leptoquarks at the LHC

Measurements of single Higgs production and its decays are in good agreement with the Standard Model. There is still room for large modifications in double Higgs production at LHC, though these effects may be correlated with large corrections to other observables, in particular single Higgs production. In this work we address the issue of enhancing double Higgs production in the presence of scalar leptoquarks while satisfying all experimental constraints. We show at leading order that including more than one species of leptoquarks, large cubic interactions with the Higgs can lead to sizable enhancement of di-Higgs production cross section at LHC, while at the same time keeping other Higgs observables and precision measurements under control. For masses above 800 GeV these corrections are in general below 30%, whereas in a viable scenario in which one of the leptoquarks can be light, specifically in the mass range $400-600$ GeV, we show that it is possible to roughly double the SM cross section for di-Higgs production, implying that possible first hints of it may be probed at the high luminosity LHC at $\mathcal{L}\sim 2$ ab$^{-1}$.


Introduction
A particle which highly resembles the Higgs boson of the Standard Model (SM) of particle physics has been discovered at the LHC in 2012 [1,2] and since then, many of its couplings to the rest of the SM spectrum have been constrained by studying its single production and dominant decays [3,4], pointing towards its SM-like nature. In fact no sign of new physics has been discovered at the LHC by direct or indirect probes so far, but nonetheless there still remains room for the possibility that some of the Higgs interactions may be affected by new physics.
We consider in this work an example of the last situation, how the double Higgs production process at the LHC can possibly be enhanced in a model of beyond the SM physics (BSM) that contains scalar leptoquarks (LQ) which interact with the Higgs and SM fermions. Leptoquarks have been also considered in connection with the phenomenology of the Higgs boson in Refs. [5][6][7][8]. It is well known that within the SM and in the present experiments, the multiple Higgs production process is the only way to access the interactions provided by the potential responsible for the spontaneous breaking of the Electroweak (EW) gauge symmetry. However, double Higgs production cross section in the SM case will only be probed at the high luminosity LHC (L ∼ 3000 fb −1 ) [9][10][11][12]. On the other hand, LQs can run in the loops of gluon fusion diagrams since they are scalars that carry color charge and may also posses large cubic and quartic interactions that involve the Higgs, allowing an enhancement in the double Higgs production [13].
LQs have also recently been advocated as some of the best candidates to explain leptonic anomalies, as the anomalous magnetic moment of the muon [14,15] 1 and decays of B-mesons [17][18][19][20][21][22][23][24]. However, due to the leptoquark nature, they are usually involved in many flavor changing transitions which constrain couplings, mostly to 1st and 2nd generations. There are also strong constraints for LQs coming from their direct QCD production and subsequent decay into SM fermions as well indirect constraints to some of the possible cubic and quartic interactions allowed in the scalar sector due to T -parameter contributions and potential contributions to the Zbb-coupling [25]. Since, as we will show in the following sections, one needs large values of the cubic interactions involving Higgs and LQs in order to enhance double Higgs production, there are perturbativity limits on how large these can be, as well as possible implications in a destabilization of the electroweak vacuum, both of which we address. Finally there can be a strong tension in enhancing double Higgs production and at the same time keeping under control the contributions to single Higgs production, a delicate problem on which a central part of this work is dedicated to.
We take into account the bounds described in the previous paragraph and, by performing a numerical study which in no way is intended to be thorough, we show that it is indeed possible in our BSM model to obtain cross sections for double Higgs production which can be slightly larger than 2 times the SM one, while at the same time satisfying all experimental constraints. This provides for an opportunity of observing double Higgs production at LHC at a lower luminosity than the ones required for the SM case, possibly reaching discovery levels at L < 3000 fb −1 .
The paper is organized as follows. In Sec. 2 we introduce the model of LQs that we will use in order to enhance the double Higgs production at the LHC, considering stability constraints, working out the physical states and their mixing, and writing out the LQs interactions with fermions. Sec. 3 deals with the proper calculation of the double Higgs production at the LHC in our model, taking into account as well the contributions to single Higgs production which constrain the enhancement. In this section we also introduce a very useful limit known as the low energy theorem (LET), in which the interactions between the Higgs and gluons take a simple form that allow us to understand better the physical behavior in certain regions of the parameter space. In Sec. 4 we focus on the constraints on the model from possible oblique corrections to electroweak precision tests (EWPT), LQs contributions to the well measured value of the Zbb coupling, flavor changing transitions and LQ direct searches. We show and analyze our numerical results in Sec. 5 in the case in which one of the LQ remains light, whereas Sec. 6 is devoted to the heavy LQs scenario. In Sec. 7 we briefly comment on how some of our LQs could address the recently B-anomalies measurements. Finally we give our conclusions in Sec. 8.

Model
The aim of our work consists in constructing and studying, from a phenomenological bottom-up approach, a model beyond the SM that includes only a small number of LQs which provide an enhancement in the double Higgs production at the LHC. In that spirit, we consider a theory in which we add to the SM content the following three scalar LQs following the notation used in [25] R 2 ∼ (3, 2) 1/6 , S 1 ∼ (3, 1) 1/3 ,S 1 ∼ (3, 1) −2/3 , where the numbers in parenthesis represent the SU (3) C × SU (2) L SM gauge group transformation properties, and the subindex the LQ hypercharge. We denote generically the scalar fields as: φ = H,R 2 , S 1 ,S 1 , and use greek: α, β, . . . , and latin: i, j, . . . indices for color and weak isospin in the fundamental representation, respectively. The reason behind the choice of these particular LQs becomes clearer once we write the most general renormalizable potential, and in particular the possible cubic and quartic interactions between the LQs and the SM Higgs [25,26] with γk 2 H * αβγ ij k + h.c. , (2.5) and where i 1 ...i d is the Levi-Civita tensor in d dimensions (we use the convention 12...d = +1), σ i are the Pauli matrices and t a the Gell-Mann color matrices. Note how in particular the terms in Eq. (2.4) and some of Eq. (2.5) (those involving two Higgs fields H) lead to cubic and quartic interactions, respectively, that provide contributions to the 1-loop gluon fusion diagrams for double Higgs production. In particular the cubic interactions, before EW symmetry breaking, contribute to the double cross section but not to the single one, showing a hint on the enhancement of double production in the presence of large cubic couplings. For the LHC phenomenology that we are interested in, we can neglect the quartic terms, except those with couplings: λ H , λ 1,2,3 and λ 1 , where λ H is the quartic Higgs self-coupling. These terms, when the Higgs vacuum expectation value (vev) is considered, contribute to the LQ masses and they also can be relevant to single and double Higgs production. It is clear from this potential that in principle all scalars involved could possibly obtain a vev in regions of the parameter space. In order to have proper Electroweak symmetry breaking (EWSB) taking place and no color/EM charge breaking minima stemming from the LQs acquiring a vev, we will in the following show the conditions that must be demanded to the couplings that enter in the potential [27][28][29].
Given the large amount of parameters that enter in the potential and the fact that we are interested in regions in which the LQs acquire vanishing vevs, we do an analysis of the extrema of the potential in which we neglect the quartic contributions, given that they will be only important in regions of large field values. We do require though that all of the quartic couplings are real, λ H and λ 1,2,3 positive, and the condition λ 1 > |λ 1 | in order to guarantee a potential bounded from below. Thus the chances of obtaining a non-vanishing value for R 2 , S 1 and S 1 depend mostly on the parameters µ 1 , µ 2 , m 2 R 2 , m 2 S 1 and m 2 S 1 . The extrema condition then simplifies to which in the unitary gauge for H and separating the up and down components ofR 2 leads to (2.10) The system described in Eqs. (2.9) and (2.10) have a vanishing solution for the LQs vevs whenever These mass matrices have to be simultaneously diagonalized. We assume that all parameters are real. We can find the eigenvalues and eigenvectors of both mass matrices which represent the physical fields analytically with the eigenvalues as A simple relationship can be found between the mixing angles and the mass eigenvalues Through these sets of equations and from a practical point of view, we are able to write the Lagrangian parameters in terms of the physical mass eigenvalues and mixing angles. This turns out to be very helpful when calculating numerically the Higgs production at the LHC. In Appendix A we derive the couplings between the Higgs and LQs in the physical basis, and in Appendix B we present a Naive Dimensional Analysis (NDA) of the size of the couplings to be under control in the perturbative expansions. Finally, we write the Lagrangian for the LQs and the SM fermions in the interaction basis. In the absence of ν R We have not included generation indices and the z-coupling is antisymmetric. The presence of ν R would introduce new interactions, opening more decay channels.

Higgs phenomenology
The discovery of a scalar particle compatible with the Standard Model Higgs boson at the Large Hadron Collider (LHC) has motivated searches for new physics using the Higgs boson as a probe. In particular, Higgs pair production is not only interesting as a probe of the trilinear Higgs self-coupling, but because its rate could be sensible to the presence of new physics effects. Several mechanisms give rise to the production of pairs of neutral Higgs bosons in hadrons collisions. Multiple Higgs states can be excited through Higgs-strahlung off W bosons or through W W fusion in proton-proton collisions but, the largest production rates are provided by gluon fusion. Since Higgs bosons are colorless particles, their pair production from the collision of two gluons is necessarily mediated through loop contributions which could be permeable to new physics effects. In this section, we explore the Higgs production in gluon fusion under the influence of scalar LQs. First, we focus on the dissection of the ultraviolet complete calculation of the double Higgs production cross section, at leading order on α S . Then, we review it but from the perspective of the Higgs LET framework, comparing both results to check the limit of this approximation. Finally, we take advantage of the simplicity of the results obtained through the LET and evaluate the correlation between the double Higgs production cross section and the single Higgs production, which is certainly more constrained by the experimental data. We also consider Higgs decays into a pair of photons and into a Z boson and a photon, where in both cases the LQs contribute at one loop level. However, we observe that their impact on the phenomenology is mild.
All the results shown throughout the paper are given at leading order. Notice, however, that when considering the phenomenological impact of the double Higgs production cross section in Secs. 5 and 6 we take ratios against the SM prediction at leading order. Therefore, we expect that most high order corrections (in particular QCD ones) drop in the ratio. 2

Double Higgs production
The Lorentz invariant differential cross section of gluon fusion into Higgs pairs could be split into two main pieces as follows [48][49][50][51] whereŝ represents the invariant mass of the Higgs pair finally produced andt is the transferred squared momentum from one of the gluons in the initial state to one of the Higgs bosons states. The form factors, F and G, hold the contributions linked to the same or the opposite polarization configurations of the incoming gluons, respectively. Each of these form factors could be generically written as the sum of several terms related with various loop configurations Higher order QCD corrections to the double Higgs production have been performed under different approximations, for example in [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].
The first two terms of each form factor are provided by the Standard Model, through the contributions mediated by triangle and box loops of a top quark. Taking into account that this form factors are proportional to the mass of the quark involved in the loop, the remaining flavors bring negligible contributions. Representative diagrams are shown in Fig. 1. The remaining terms arise from the loops induced by the presence of the mass eigenstate LQs, χ u 1 , χ u 2 , χ d 1 and χ d 2 , presented in Sec. 2, and their interaction with the Higgs boson. The generic diagrams that contribute to these terms are represented in Fig. 2. It is worth mentioning that the new trilinear effective couplings enable the possibility of blending different types of physical LQs in the loops, see Fig. 2c and 2d. This feature is not present in the Standard Model diagrams due to the diagonal nature of the Higgs couplings to quarks. The dependence of F , F , G , G and the new physics form factors on the Mandelstam variablesŝ,t andû can be found in the Appendix C.  In the upper frame of Fig. 3, we can see how the differential cross section of Eq. (3.1) looks as a function of the center-of-mass energy of the gluon system √ŝ , after being in-tegrated int. The red dashed line represents the purely Standard Model result. At the leading order, the shape of the cross section is basically dominated by the squared module of the box diagrams and the destructive interference between the former and the triangle diagram. The blue solid line is produced by the inclusion of a pair of scalar LQs in the loop, obtained from the diagonalization of the mass matrices of Eqs. (2.12) and (2.13) with arbitrary parameters that satisfy all the experimental constraints considered in this work. The Standard Model cross section is reshaped by the presence of two resonance peaks at twice the masses of the LQs, taken around 0.5 TeV and 2 TeV. The interference between the new physics form factors and the Standard Model contribution is responsible for the final shape of the resonance peaks and the enhancement of the cross-section at low values ofŝ, see the lower panel of Fig. 3. As it can also be appreciated in this figure, the gluon-gluon luminosity dL gg /dŝ increases extremely fast asŝ decreases. Consequently, the enhancement produced by the interference between the Standard Model and the new physics contributions results especially relevant to boost the double Higgs production in proton-proton collisions. The total cross section for Higgs pair production can be derived by integrating overt and performing the convolution of the parton cross section with the gluon-gluon luminosity as follows where s denotes the center-of-mass energy of the proton to proton collision that gives rise to the gluon fusion process and where we made use of the relationŝ = τ s. The limits To compute the gluon luminosity dL gg /dτ we made use of the set of parton distribution functions (PDF) PDF4LHC15 nnlo mc, provided by the PDF4LHC working group [52] at next-to-next-toleading order.

Higgs low energy theorem
The LET allows to obtain the interactions with photons and gluons generated at loop level by heavy particles, in the limit where their masses are large compared with the other scales [53,54]. We are interested in the couplings of one and two Higgs with two gluons, involved in single and double Higgs production. In the following we include the top quark within the heavy particles, as well as the LQs, although it is well known that the heavy top mass limit disagrees with the full calculation by 20% for di-Higgs production. Similar discussions of LETs with colored scalars in the context of supersymmetry have been considered, for example, in Ref. [55].
For the interactions of the Higgs with gluons we use the following effective Lagrangian .

(3.5)
Replacing v → v + h and expanding m i in powers of h where i runs over all the particles, β i = 2/3, 1/6 for Dirac fermions and complex scalars, respectively, are the coefficients of the β function of QCD, and s labels the species of particles: up-type LQs, down-type LQs and the top (we neglect the effect of the other SM fields).
For an eigenmass LQ χ, the coupling |χ| 2 h can be obtained from: where g t h arises from the top quark contribution: g t h = 2. 3 For our model with LQs, we split the contributions of the up-and down-type states .
These effective couplings can be computed straightforwardly from the mass matrix. Thus, we have a simple explicit expression for these couplings in terms of the parameters of the model Similar Eqs. can be obtained for the contributions to g hh : g LQu hh and g LQd hh , by using Eq. (3.8).
It results quite illuminating to rewrite the expressions for g h and g hh in terms of the physical masses. For that purpose it is useful to consider that Considering only the LQ contributions, concentrating initially in the calculation of g h and using the expressions for the LQ physical masses, it is easy to calculate how the cubic and quartic interactions contribute from the up and down sectors quartic . Notice how the second term inside the parenthesis in both g u h, quartic and g d h, quartic tend to be smaller than the first term. Comparing against the sign from the cubic contributions, this implies that the quartic tend to contribute with opposite sign. This opposing behavior is seen even in the numerical scans in the full UV theory. Furthermore, for µ 2 2 m 2 , the cubic contributions dominate over the quartic ones, g h, cubic g h, quartic . In order to obtain a simple expression and assuming that the cubic terms dominate, we can also calculate the contributions to g LQ (3.19)

Double Higgs production
The calculation of double Higgs production cross section using LETs is greatly simplified. The amplitude is simply given by the sum of two terms, one with the hhgg vertex, and another one with the hgg vertex, a Higgs propagator and the Higgs trilinear coupling. To each kind of term there are contributions from the top and LQs. The hadronic cross section is [56] where the center of mass energy of the collider s is related withŝ by:ŝ = τ s and (3.22) Q and µ are the factorization and renormalization scales that we take as the center of mass energy of the system of two Higgs states: Q = µ = √ŝ . We use the MSTW2008 parton distribution functions for the calculations with LET [57].
Since g h and g hh are independent of momentum, using the PDF of gluons we can integrate Eq. (3.20) and obtain an explicit expression for σ hh in terms of the parameters of the model as hh are integrals of the PDF, independent of the parameters of the NP sector, whereas g h and g hh contain all the information of the NP, Eqs. (3.7) and (3.8).
In Fig. 4 we show the results of the full calculation and the calculation with LET, there are several interesting features to comment. On the left we show the results for a scenario where the mass of one of the LQs is taken as m light LQ 400 GeV, whereas the rest of them are larger than 800 GeV. Although the LET correctly predicts the order of magnitude of the cross section, there is a dispersion for large values of σ hh , that corresponds to m light LQ being near the lowest value. This dispersion can be expected, since for light masses the assumptions for the LET are not valid anymore. There is another correction, that is responsible for the departures of the points from a straight line, that is: for low LQ masses the interference term is sizable, but the LET fails to reproduce it properly. Given that the top mass is not heavy enough compared with other scales of the process, the LET overestimates the interference by a factor 2. On the right we show the case where all the LQ masses are larger than 800 GeV, such that they are heavier than the Higgs mass and the energy flowing through the loop. Unlike the previous case, we have divided by two the interference term, such that the vertical axis corresponds to: σ LET|SM hh . The red line is the identity, with a shift to compensate the top contribution that, in the LET, underestimates the cross section by 20% approximately.
, that is: we correct the interference term by a factor 1/2, the red line is the identity, shifted by a constant value to compensate the underestimation of the SM cross section using the LET.
To show that the origin of the factor two in the interference is the heavy top mass approximation, we have checked that, by making the top mass larger, both calculations converge. In Fig. 5 we show the ratio of interference term using the LET and using the full calculation: , as function of the top mass, for a point of the parameter space with: m light LQ = 1.3 TeV. This ratio is of order 1/2 for the SM top mass, and converges to 1 for large m t and LQ masses. In the figure, the curve does not converge to 1 because the LQ masses are finite, but we have checked that as we increase their values it goes to 1.
In Fig. 6 we show only the LQ (SM-subtracted) contribution to the cross section normalized with respect to the SM one, , as function of the physical mass of the lightest LQ state. We have taken λ 1,2,3 = λ 1 = µ 1 = 0 and different values for µ 2 , as shown by the different curves. We have varied mS 1 , whereas the other quadratic coefficients are fixed to 2 TeV. For this choice of parameters the mass of the light mass eigenstate is driven by mS 1 . As expected, for lighter mass and large µ 2 , the cross  The Fig. 7 shows the ratio σ LET hh /σ LET|SM hh vs. µ 2 , using the full LET result for the double Higgs production in purple and using the cubic approximations to g h and g hh in blue, Eqs. (3.14) and (3.18), respectively, replaced in Eq. (3.23), considering in both cases only contributions from the up-type sector. Note that the approximation does a good job in mimicking the full trend of the LET result, however it tends to overestimate its value. This is expected since in the approximation we are neglecting the quartic contributions which, as mentioned previously, provide an opposite behavior with respect to the cubic ones, decreasing the overall value of the cross section.
to LO as a function of µ 2 . The purple dots correspond to the full LET computation and the blue ones to the inclusion of µ 2 terms only.

Higgs couplings
The LQs interactions with the Higgs can significantly affect the gluon fusion Higgs production and this impact can be quantified through the ratio between the total and the SM cross sections 5 which is expressed in terms of the top (A 1/2 ) and LQs (A 0 ) one-loop functions where τ i = m 2 h /(4m 2 i ) for i = t, χ, and the one-loop functions are normalized as A 1/2 (0) = 4/3 and A 0 (0) = 1/3. Defining κ g = 1 + δκ g we separate the LQs contributions The numerical factor in the right hand side results from evaluating all the LQ masses at 800 GeV 6 , meaning that the LQs contributions reproduce the ones obtained through the 5 At leading order the LQs effects on the gluon fusion Higgs cross section can be also measured as the ratio between the modified and the SM gluonic partial decay widths of the Higgs (|κg| 2 = Γ(h → gg)/Γ(h → gg) SM ) following the relation among the gluon fusion cross section and the partial decay width into gluons [58]. 6 This is a good approximation even if not all the LQs have a mass of 800 GeV. The extreme case appears in the scenario analyzed in Sec. 5 where the mass of the lightest LQ (∼ 400 GeV) produces a correction lower than 1% to Eq. (3.26).
LET discussed in the previous section, see Eq. (3.9). The sum is taken over the four LQs physical states. Accordingly, C (k) ii stand for the trilinear scalar couplings in the physical basis. Following Eqs. (2.4) and (2.5), we can see that C (k) ii depends after mixing on the cubic and quartic (once one of the Higgs is evaluated at the vev) scalar couplings in the potential. 7 It is simple to derive an expression that relates the LQ contributions to the single Higgsdigluon coupling, δκ g , and the double Higgs cross section, σ hh , as described in Eq. (3.23) using the LET and considering only a cubic contribution from the up-sector. Noticing that δκ g = (1/8)g LQ h , one can write g h and g hh (note that these include the top quark contribution, see Sec. 3.2) in terms of δκ g as Replacing these expressions in Eq. (3.23) for σ hh we get the dependence of the cross section on δκ g . Furthermore, the LQs couplings with the Higgs can produce important modifications to the partial decay width of the Higgs into a pair of photons. In this case, these variations are captured by the decay width normalized to its SM value 29) which is determined by the W one-loop function, A 1 , as well as by A 1/2 and A 0 for the top and LQs contributions, respectively, 30) where N c is the number of colors of the LQs and Q χ their electric charges, and A 1 is normalized as A 1 (0) = −7. Analogously to Eq. (3.26), we define in this case κ γ = 1 + δκ γ The numerical factor in the right hand side results again now from evaluating all the LQ masses at 800 GeV. From Eqs. (3.26) and (3.31) we see that a combination of large trilinear couplings and small LQ masses produces larger deviations and, consequently, it tends to put these observables under more tension. Both κ g and κ γ have been measured at the LHC [59]. We will analyze in the next sections the constraints imposed by the measurements on the parameter space of the model and show the expected correlation between κ g and κ γ [25].
Finally, another possible important LQ effect might appear in the partial decay width of the Higgs into a Z and a photon, Γ(h → Zγ). In this case, however, we expect a contribution of the same order as in Γ(h → γγ) [25] with smaller experimental restrictions [60]. As we show in Sec. 5.2, Γ(h → γγ) introduces no more than slight constraints on the parameter space in regions where the double Higgs production is enhanced, therefore, Γ(h → Zγ) ends up having a negligible impact on our analysis.

Constraints
In this section we analyze the restrictions imposed on the model by different observables. In particular, we concentrate on the constraints arising from the oblique corrections to electroweak precision tests, LQs contributions to the Zbb coupling, flavor changing transitions and LQ direct searches. We also study the requirements for baryon and lepton number conservation and define the scenarios of interest for the phenomenological analysis.

Oblique corrections
The LQs give contributions to the self-energies of the EW gauge bosons at one-loop level [25,61,62]. The leading constraints to these quantities are captured by the S and T parameters [63,64]: S = 0.05 ± 0.11 and T = 0.09 ± 0.13.
There are two sources for the T -parameter: the splitting of the doublet induced by λ 1 and the mixing induced by µ 1,2 . Defining and using Eqs. (2.14) and (2.15), we obtain Notice that for m χ ∼ O(0.5 − 2) TeV and splittings between the components of a doublet of O(1), T can take values 1, thus T is expected to give interesting constraints on the parameter space of the model.
In order to understand how λ 1 and µ 1,2 generate contributions to T , one can consider the limit µ 1,2 v m 2 where we show the first corrections generated by λ 1 and µ 1,2 in the first line, and in the second line a mixed correction. We have checked by comparison with the full result at a numerical level that, in the proper region of parameter space, this approximation gives a good estimate of T . Similar expressions can be derived for the S-parameter. Defining we obtain We do not expect large corrections to S. Considering an example, for a single doublet: 10) gives ∆S ∈ (−0.06, +0.06). For our model and region of parameter space S does not lead to important constraints.

Z-couplings
LQs also give corrections at one-loop level on vertex interactions of SM fermions and EW gauge bosons [65,66]. Precision measurements at LEP give very stringent constraints on the modifications of the Z-couplings [67]. Focussing on Zbb, Ref. [25] has shown that for m LQ 0.5 TeV, the couplings to fermions are required to be 10. Since in our approach those couplings will be taken small, the constraints from Z-couplings can be easily avoided.

Direct searches
Bounds from direct searches of LQs at colliders depend on the final state. For pair production, that is dominated by QCD, these bounds are usually presented in terms of the BRs.
Final states with top quarks and e or µ have stronger bounds, with m LQ 1.6 TeV for BRs to light leptons 1 [68]. LQs dominantly coupled to light quarks and leptons are also excluded up to m LQ 1.6 TeV [69,70]. For decays to quarks and leptons of the third generation: m LQ 800 − 900 GeV [71]. For the decay to a pair of down type quarks, [72] excludes masses below 410 GeV and 610 GeV for decays to qq and bq, respectively.
The above searches assume that the LQs decay promptly and no longer apply if they are long lived. In that case, bounds arising from searches of long lived squarks must be considered instead. Different final states involving decays into down-type quarks and a lepton [73,74] or two down-type quarks [74][75][76] have been explored. All the exclusion limits reported in these searches restrict the long-lived particle mass to be well above 1 TeV. For example, in [74] top squark masses up to 1.6 TeV are excluded for mean proper decay lengths between 3 and 300 mm.

Baryon and Lepton numbers
The interactions of LQs can violate baryon and lepton number, inducing processes forbidden in the SM, as proton decay and neutron-antineutron oscillations. In fact, if one includes all the renormalizable interactions conserving the gauge symmetries of the SM and the Poincaré symmetry, then B and L are violated. Below we discuss different cases where these symmetries can be conserved, as well as situations in which they are violated, but proton decay can be sufficiently suppressed. In the rest of the paper we will concentrate in the first and second cases of the B and L conserving scenarios.

B and L conservation
Taking different combinations of vanishing couplings in Eq. (2.18), the LQs can be assigned B and L charges such that these numbers are conserved.
1. Taking y 2 = z S 1 q L = z S 1 q R = µ 1 = 0, one can asign: B(R 2 ) = −B(S 1 ) = 2B(S 1 ) = −2/3, L(R 2 ) = L(S 1 ) = 0 and L(S 1 ) = −1. Notice that in this caseR 2 andS 1 do not interact with leptons at the level of dimension four operators. For phenomenology this is a very interesting situation, since it allows one of the cubic interactions, whereas R 2 andS 1 decay to dijets and S 1 decays to q . These assignments and a lighter LQ, with a mass ∼ 400 − 600 GeV, define our scenario 1, dubbed "Scenario with a light LQ", and it will be analyzed in Sec. 5. In this case S 1 and S 1 decay to q after mixing. This is also a very interesting situation, since it allows both cubic interactions, as well as decay of all the LQs. These assignments and LQs with masses 800 GeV, define our scenario 2, dubbed "Scenario with heavy LQs", and it will be analyzed in Sec. 6.

Taking vanishing interactions with fermions
3. Taking zS 1 d = z S 1 q L = z S 1 q R = µ 1 = 0, one can assign: B(R 2 ) = −B(S 1 ) = −B(S 1 ) = 1/3 and L(R 2 ) = −L(S 1 ) = L(S 1 ) = −1. Diquark interactions are not present andS 1 decays after mixing with the up-component ofR 2 to q . Regarding the LQ-Higgs interactions at leading order, this scenario is already included in the "Scenario with heavy LQs". However, it is disfavored from the phenomenological perspective since it is constrained by the same bounds on the LQ masses arising from the direct searches and the condition µ 1 = 0 implies a smaller contribution to the double Higgs production.

B and L violation
In the general case dimension six operators as dduν can be induced at 1-loop level, as well as dimension nine operators as ddudue at tree level [77]. They can produce n → π 0 ν, π + e, nn oscillations and most dangerous: p → π + π + e − . Proton decay is proportional to |y 2 zS 1 d | 2 /m 4 LQ , thus it can be suppressed by taking y 2 = 0 for d and s, leaving only coupling with b, or taking y 2 small enough to satisfy the constraints.

Flavor changing transitions
Since the LQs interact with the SM fermions, they can induce flavor transitions and CP violation. At tree level they can contribute to processes involving quarks and leptons, as semileptonic and leptonic decays of mesons, decay of heavy leptons as the τ and transitions in nuclei. At loop level they can also contribute to processes as meson mixing, dipole moments and rare decays. The absence of conclusive evidence of physics BSM in the flavor sector strongly constrains the flavor structure and the size of the couplings to fermions. Below we give estimates for the most constraining processes, that usually involve light fermions, and we refer to [25] for an extended list of processes as well as the details of the calculations.
Semileptonic interactions can be generated by tree level exchange of LQs, their coefficients being estimated to be of order C 4f ∼ y 2 /m 2 LQ , with y a generic coupling to fermions. Leptonic decays of pseudoscalar mesons are some of the most sensitive processes receiving contributions from these interactions. The ratio R π µ/e = Γ(π → eν)/Γ(π → µν) can be predicted with good precision in the SM due to cancellation of uncertainties, leading after comparison with the experimental result to the estimates [78]: |Re(y * S 1 e,11 y S 1 ,11 )| 10 −7 (m LQ /TeV) 2 and |Re(y * S 1 e,12 y S 1 ,12 )| 2 × 10 −3 (m LQ /TeV) 2 . Four-fermion operators involving only quarks are generated at one-loop level, with Wilson coeficients that roughly can be estimated as: C 4f ∼ y 4 /(4πm LQ ) 2 [79]. To get an idea on the size of the allowed couplings, in the absence of flavor symmetries, one can consider one of the most stringent constraints that arises from K in the Kaon system, where for complex couplings C 4f 7 × 10 −9 TeV −2 [80]. Assuming that d and s quarks have couplings of similar size, this bound leads to y 6 × 10 −3 , when m LQ 500 GeV. 8 These bounds are relaxed by one order of magnitude if the couplings are real, and at least by one order of magnitude for mesons involving heavy quarks.
The electromagnetic dipole moment (EDM) of the electron has a very tight constraint. LQs with couplings to Left and Right fermions can induce corrections to the EDM at oneloop level [81], leading to strong bounds on the imaginary component of the couplings to electrons, as can be the case of S 1 . Following [25] leads to: Im[y S 1 e,31 (V t y S 1 ) * 31 ] 10 −12 , where the contribution is dominated by the top quark in the loop due to the requirement of a chiral flip. Neglecting Cabibbo suppressed terms and assuming anarchic complex couplings of the same order, one obtains: y S 1 e/ , 31 10 −6 . However these bounds are relaxed by many orders of magnitude for quarks and leptons of second and third generations. Notice that the previous bound requires the presence of two different Yukawa interactions, thus, in the absence of two couplings, the contribution is instead suppressed by the mass of the lepton, eliminating the tension for the case of the electron due to its small mass.
At first order the effect of LQs on the Higgs physics is independent of these couplings, thus we will assume that they are small enough to avoid issues with flavor. In the limit of vanishing couplings the LQs are stable, by demanding their decay length to be smaller than min requires y, z 10 −6 (TeV/m LQ ) 2 (0.1mm/ min )(E/6.5TeV), where E is the energy of the LQ that has been taken of order 6.5 TeV, that is the full available energy at LHC13 in pair production. Thus we will assume that at least one of the couplings that allows the decay of the LQs is on the proper window to satisfy flavor bounds and simultaneously allow a decay length smaller than 0.1mm. In this manner, the LQ masses considered in scenarios 1 and 2 are allowed by the direct searches at colliders since we avoid the more stringent constraints arising from long-lived particles searches (see Sec. 4.3).
In Secs. 5 and 6 we will discuss some specific scenarios, showing configurations of couplings with fermions that are compatible with these conditions.

Scenario with a light LQ
In this section we show the predictions for several interesting quantities in the scenario 1 of Sec. 4.4.1. This scenario allows a light LQ, with a mass ∼ 400 − 600 GeV, decaying to dijets. Since its decay to bottom quarks must be suppressed, the simplest choice is to take zS 1 d,i3 = 0 in Eq. (2.18), with zS 1 d,12 driving the decays ofS 1 andR 2 . The up-type components are mixed by µ 2 , whereas the down-type LQs do not mix, since µ 1 = 0. S 1 interacts with the other scalars through quartic couplings or exchange of other particles. The simplest choice is to assume interactions of S 1 with fermions of the third generation only.
We have kept only points with the lightest up-type LQ being heavier than 400 GeV and the rest of the LQs being heavier than 800 GeV. For the points that we show the mass of lightest state takes values 400 − 837 GeV, with most of the points corresponding to masses lighter than 600 GeV, as can be seen in Figs. 8 and 10. We find that the heavy up-type state has a mass ∼ 1.9 − 2.5 TeV, the light down-type state in general is identified with S 1 , except in some cases for masses above 2 TeV, and it has mass ∼ 0.86 − 2.52 TeV, whereas the heavy down-type LQ has mass ∼ 1.9 − 3.5 TeV, and in general corresponds to R d 2 . The mixing angle in the up-sector is smaller than 0.26.

Double Higgs production
We show in this section the most interesting results from the numerical scan for the scenario with a light LQ using the physical notion from the LET limit where possible. In Fig. 8 on the left we show the ratio between the total cross section and the SM cross section σ hh /σ SM hh against the lightest LQ mass, m light LQ , in TeV. It is clear from the figure that, as the lightest LQ mass decreases, the possibility of enlarging the total cross section increases, reaching values of order 2σ SM hh for LQ masses m light LQ < 0.5 TeV. However, lighter LQ masses do not necessarily imply larger cross sections. Also notice that when m light LQ > 0.7 TeV the cross sections are less than a 20% larger than the SM cross section. We have colored the points in order of increasing value of |µ 2 | from green to blue. Notice that the low values of |µ 2 | (green points) imply smaller cross sections whereas the points with the largest cross sections are associated with the largest values of |µ 2 |. In Fig. 8   In Fig. 9 we show σ hh /σ SM hh vs. µ 2 in TeV. Notice that the increments in the total cross section can be accomplished for values of |µ 2 | > 3 TeV, roughly independently of the sign of µ 2 . As larger values of µ 2 are considered, the increment in the total cross sections can be sharp for some of the points, following a quartic dependence on µ 2 , as we saw as well for the LET case.
In Fig. 10 we show σ hh /σ SM hh vs. δκ g , coding the colored points in order of increasing value of m light LQ from green to blue, including the LET result depicted by the blue line for σ LET hh (δκ g ) in the cubic dominated case as given in Eq. (3.23) with the replacements of Eqs. (3.27) and (3.28). We can clearly see that the LET describes nicely the main trend of the figure, in particular for most of the points with m light LQ > 0.51 TeV, for which the LET works well and in which it has been proven the strong correlation between the double Higgs cross section and the modifications to the Higgs-gluon coupling. The deviations from this behavior can be seen for the points that spread vertically for a fixed value of δκ g , and are represented by the green points which form a sort of a triangular region above the curve. This difference from the LET behavior happens due to the lightest LQ mass decreasing, entering a regime in which the LET breaks down and is unable to capture the full picture. Note also that the vertical break at δκ g ≈ −0.12 corresponds to the constraint on single Higgs production as depicted in Fig. 11 and that for m light LQ > 0.7 TeV, σ SM hh σ hh 1.2σ SM hh which agrees with what we find using the LET as shown in Fig. 6 once the quartic contributions are accounted for, which tend to contribute oppositely to the cubic, lowering the cross section.
We can estimate at what luminosities the √ s = 14 TeV LHC would be providing the first hints of double Higgs production for the largest cross sections predicted in this scenario with a light LQ. For that purpose, we use the projections by ATLAS [84][85][86] and CMS [87] that consider different final states but under the best case scenarios seem to be able to constrain σ hh 1.5 σ SM hh at 95% Confidence Level and at a luminosity of L ∼ 3 ab −1 . Given that the largest cross sections we can accomplish within this scenario are σ hh ∼ 2.3 σ SM hh , a relation we expect not to change sizeably after taking higher order corrections, and that the final states are the same as in the SM case, we estimate that for a luminosity L ∼ 2 ab −1 , the first hints (or the exclusion at 95% Confidence Level) should be able to be seen at the LHC.

Higgs couplings
We study next the correlation between κ g and κ γ . We also consider the dependence of κ g on the cubic coupling of the light up-type state which dominates the LQs contributions within this scenario. We show in Fig. 11 δκ γ vs. δκ g and the corresponding 1σ and 2σ experimental contours [59]. A first observation is that κ g enforces the stronger constraints, whereas κ γ plays a weaker role in imposing restrictions. In order to evaluate the double Higgs cross section, we have allowed deviations up to the 2σ level in this plane. Being the blue points the ones with largest values for this cross section, it is apparent from the figure that this is achieved at the expense of reaching the 2σ contour in the region of lower values of δκ g , in accordance to Fig. 10. As it can be also seen from this figure, it is pointless to explore larger positive values of δκ g because they are located in a region where the corrections to the double Higgs cross section are small, as it is shown by the green points.
In Fig. 12 we plot δκ g vs. C (u) 22 , the cubic coupling of the light up-type state. As it is shown in this figure, the blue points corresponding to the larger values of the double Higgs cross section have the larger negative values of this coupling, effect that can be followed from Eq. (3.26). In agreement with Fig. 11, positive values of δκ g are associated to smaller cross sections and are originated from positive values of C (u) 22 . The horizontal break around δκ g ≈ −0.12 corresponds to the constraint on the single Higgs production as shown in Fig. 11.

T parameter
We consider now the value of the T -parameter for this scenario. As discussed in Sec. 4.1, T is driven by λ 1 and µ i , thus we show in Fig. 13 this dependence. For each plot, we have checked that the approximations of Eq. (4.3) reproduce quite well the curve wrapping the cloud of points from below, while the points with larger values of T correspond to regions of the parameter space where the approximation is not valid, due to large µ 2 and/or λ 1 . For |µ 2 | 6 TeV most of the points exceed the upper limit on T , while one could take larger values of λ 1 , at the price of lowering |µ 2 |. We do not find correlation between T Figure 11: Values of δκ γ and δκ g for the points of the scan, in colors the values of σ hh /σ SM hh . The blue and violet lines correspond to the 1σ and 2σ experimental contours from [59]. The red dot is the central value of the measurement.

Scenario with heavy LQs
We briefly discuss now the scenario 2 of Sec. 4.4.1, where up to a sign, all the LQs have the same baryon and lepton numbers. The only coupling with fermions is y 2 and the simplest situation is to consider interactions with the fermions of the third generation only. Since in this case only one coupling to fermions is present, many flavor bounds, as well as bounds from direct searches, are less stringent than in the most general case. It is possible to relax even more flavor bounds by taking small couplings with the fermions of the first and second generation. We have made a scan on the parameters of the model with the only restriction that the couplings are perturbative, |µ 1,2 | ≤ 2π TeV and λ 1,2,3 , λ 1 ∈ (0, 2π), and then selected points where all the LQ masses are larger than 800 GeV. The output of the scan has shown that for most of the points the corrections to the double cross section are smaller than 20%, and for a few points they can reach 30%. The T -parameter is similar to the previous scenario.
In Fig. 14 we show σ hh /σ SM hh vs. δκ g . As opposed to the scenario with a light LQ, we observe a stronger correlation among both quantities which can be followed from the LET since now all the LQs are heavy. Besides, compared to that scenario, we have lower values of σ hh for similar δκ g . As an example, for δκ g = −0.1, the effect on the cross section is of order 20 − 30%, that is similar or smaller than the lowest values obtained in the scenario with a light LQ. However, if we compare both scenarios for similar values of m light LQ , we see that slightly larger values of σ hh are obtained for the heavy LQ scenario due to the fact that both µ 1 and µ 2 are present in this case. The blue points in this figure correspond to larger values of the lightest LQ mass, m light LQ . We see that for these points the correction to the SM cross section is smaller than the one corresponding to the green points for which we allow lower values of m light LQ , as expected. As a consequence of the moderate correction to σ SM hh generated by the LQs, it is clear that this scenario will be more difficult to test in double Higgs production at the LHC.

B anomalies
Besides their interesting effects on the Higgs phenomenology, LQs can also play a role explaining a set of deviations related with flavor physics. Anomalies in the decays of B mesons have been recently measured at several experiments [17-19, 23, 24], the largest ones were observed in the ratios R K ( * ) and R D ( * ) , via neutral and charged currents, respectively. One of the most appealing possibilities to explain these deviations is the presence of LQ states at the TeV scale. There have been many proposals and studies, in particular Ref. [88] has summarized the effect of individual LQs on these observables. For the content of LQs of the present paper, only S 1 can accommodate R K ( * ) or R D ( * ) , while simultaneously satisfying the bounds from other flavor observables, by adding small couplings with the second generation. For a combined explanation of both anomalies, Refs. [89,90] have shown that S 1 can, at best, explain them at the level of 2σ. 9 Besides, given the proper couplings, it can also accommodate the anomalous magnetic moment of the muon.
In the scenario 1 of Sec. 4.4.1 (conserving B and L) the impact of S 1 on the Higgs phenomenology that we have studied is subleading. Therefore its mass and couplings with fermions can be taken to accommodate R K ( * ) or R D ( * ) , or both of them at 2σ level. To leading order in this scenario the effects on the Higgs phenomenology and on the B anomalies are decoupled, with the first ones being dominated byR 2 andS 1 , and the second ones by S 1 . Whereas current limits do not exclude an explanation of the B anomalies [94], one may wonder to which extent this possibility might be affected by single resonant LQ production searches at the LCH in the future [95][96][97]. Since the constraints arising from these searches would only impact on LQs coupled to leptons, they would not impose restrictions on the lightest LQ which in this scenario drives the increase of the double Higgs production and is mainly identified withS 1 , whose couplings to leptons we put to zero. On the other hand, this would not be the case for the heavier LQs, as S 1 . However, since their masses can be above the TeV, one could expect that couplings to leptons of second and third generation of order 0.1-1, as those required to explain the low-energy anomalies, could satisfy the present bounds. Definitely, a dedicated analysis, that is beyond the scope of our work, is required to carefully assess an explanation of the low-energy anomalies given the bounds from single resonant production, explanation that is already in tension with the present LQ content [89,90]. 10 Generically, effects on the Higgs physics are dominated by cubic and quartic couplings between the scalar bosons, while explanations of the B-anomalies and the anomalous magnetic moment of the muon involve couplings to fermions. Besides, sizable effects on double Higgs production require masses below 1 TeV, whereas explanations of the anomalies can be accomplished by LQs with masses above the TeV. For these reasons, the absence of effects on double Higgs production would not discard the LQ explanation of the B-anomalies and/or anomalous magnetic moment of the muon.

Conclusions
In the absence of direct signals of new physics at the LHC, precision measurements in the Higgs sector could provide a portal to physics beyond the SM. In this paper we have studied the Higgs phenomenology in the presence of scalar LQs at the TeV scale. Interest in these new states has been boosted in the last years due to the anomalies observed in the decays of B mesons in several experiments, as well as quite recently due to the new measurement of the anomalous magnetic moment of the muon at Fermilab, given that LQs are one of the preferred options to explain both phenomena that seem to point to BSM physics. Besides their possible impact on flavor physics, their colored nature implies that they may also modify Higgs production at hadron colliders such as the LHC, since this production is dominated by gluon fusion, in which loops of colored particles enter. Although single Higgs production has been measured at the LHC with a precision of order 10% in agreement with the SM, double Higgs production is naturally expected to be probed at much larger luminosities (of order 3 ab −1 ) and thus there is room for O(1) modifications with respect to the SM predictions. It is however a non-trivial task for the new physics responsible for potentially large modifications in double Higgs production to produce at the same time small deviations in single Higgs production. This is one of the main issues that we have considered in this paper.
We have studied a model with three scalar LQs, a weak doublet and two singlets, with hypercharge 1/6, -2/3 and 1/3, respectively. We have shown that this set of fields, alongside their cubic and quartic interactions with the Higgs, provides an opportunity for enhancing double Higgs production at LHC while simultaneously keeping under control single Higgs production. We have considered the most general renormalizable potential, analyzing under some simplifying assumptions the stability of the EW vacuum and estimating the size of the couplings where the perturbative expansions are under control. At the leading order in which we work, the contributions from the LQs to Higgs production are decoupled from the LQs interactions with SM fermions. We have determined the conditions for baryon and lepton number conservation from the LQs interactions with fermions and we have estimated the bounds on these interactions from flavor physics. Another important constraint that we have considered comes from the T -parameter, that receives contributions from the cubic and some of the quartic interactions that induce a splitting between up-and down-type LQs. These contributions have been computed and we have selected the regions of the parameter space that satisfy the bounds. It is worth mentioning that when considering the bounds from the current LHC measurements on single Higgs production and decay, we added the LQs contributions to Higgs decays that involve loop-level generated couplings to the weak sector such as to diphotons or Zγ, finding, in agreement with the literature, that these corrections are in general below the current bounds and are less constrained than the gluon fusion contributions to single Higgs production.
A full numerical calculation of the cross section of double Higgs production at LHC, at leading order, is performed in the model. Comparing our results with similar studies in the literature, both in presence of LQs and in the presence of supersymmetric particles, we find good agreement in the corresponding parameter space. Since the full calculation must be carried at the numerical level, is very expensive computational-wise and sometimes even the physics may be hard to extract, we made extensive use of the Higgs low energy theorems for heavy LQ masses. Considering a correction to the interference term between the new physics and the SM by a factor 2 due to the relatively low top quark mass compared to others scales of the process, the analysis using the LET allowed us to understand several interesting relations, as the dependence of the cross sections with some of the parameters of the potential and with some physical masses and couplings. It also helped us in the search for the desired regions of the parameter space with better efficiency. In this context, we were able to find a rather simple relation between σ hh and δκ g that is valid in the limit of dominance of cubic couplings and describes well the trend shown in the full calculation with deviations arising mostly from quartic couplings but with a clear discrepancy in the regions of small LQ masses in which the LET breaks down.
Our study focuses on two different scenarios. The first one consists of a light LQ with a mass of order 400 − 600 GeV and fermionic interactions such that it only decays to dijets, evading the strongest bounds from direct detection, while the other LQs considered have masses larger than 800 GeV satisfying all direct search bounds. Following the indications from LET and using the full calculation, we find regions of the parameter space where the cross section can be more than two times the SM one, satisfying at the same time all constraints, in particular the ones coming from Higgs single production. In these regions of large enhancement, the cross section is dominated by the contribution from the light state. Using the expected sensitivities for double Higgs production from both ATLAS and CMS at the high luminosity LHC, we estimate that hints of the largest enhancements to di-Higgs production in this scenario should start to become apparent at luminosities L ∼ 2 ab −1 . The second scenario considers all the LQs to have masses larger than 800 GeV. In this case the corrections to the double Higgs production cross section, compatibles with the bounds from single Higgs production, are in general not larger than 20% than that of the SM, although in some cases they may reach 30%.
Finally, we notice that the separation between the LQs interactions with the Higgs sector (and their consequences on Higgs production) and the LQs interactions with the fermionic sector (which would provide higher order corrections to Higgs productions and are expected to be small) implies that an explanation of the B anomalies, and/or the anomalous magnetic moment of the muon, together with an enhancement of double Higgs productions, can be accomplished in the model we study. The best scenario for the explanation of the B anomalies with scalar LQs contains a weak triplet and a singlet with hypercharges 1/3. Likewise, there is a possibility to address (g − 2) µ within a similar scenario. Although the lack of any relevant impact on the double Higgs production would not prevent to explain the B-anomalies and/or (g − 2) µ in terms of LQs, it may be interesting to perform a complete study for the possible explanations of these anomalies within our setup and we leave it for future work.

Acknowledgments
We thank Javier Mazzitelli for assistance on the numerical calculations as well as Carlos Wagner for useful discussions. L.D. thanks IFLP for its hospitality during the accomplishment of this work. This work has been partially supported by CONICET

A Higgs-LQs interactions in the physical basis
After the EWSB, the potentials V 3 and V 4 provide the Higgs-LQ interaction. These could be written as next where we use the φ = (v + h)/ √ 2 prescription and Once mass matrices are simultaneously diagonalized, the EWSB picture of the combination V 3 and V 4 can be rewritten in terms of the physical fields, χ u 1 , χ u 2 , χ d 1 and χ d 2 as follows The new matrices are defined from C (u) , C (d) , Q (u) and Q (d) , and the rotation matrices presented in Eqs. (2.14) and (2.15) The matrix elements of these matrices are the physical field coupling constants that rule the interaction between the LQ mass eigenstate fields and the Higgs boson. They could be explicitly expressed in terms of the gauge parameters and the mixing angles as follows It is worth noting that, since we assume all parameters are real, C (u) , C (d) , Q (u) and Q (d) turn out to be symmetrical.

B Perturbative expansion
In this Appendix we make some estimations of the radiative corrections and bounds on the couplings from perturbativity. Neglecting the mixing between the LQs, the cubic couplings contribute to double Higgs production but not to single Higgs production at LHC, therefore it is interesting to consider large cubic couplings. This couplings can induce large radiative corrections, and eventually break the perturbative expansion. Thus we study NDA estimates of the allowed size of the cubic couplings.
We consider the terms of V 3 in Eq. (2.4), that are corrected at one loop by the Feynman diagrams of Fig. 15(left). At this level the correction to µ-itself is given by , where we have considered the limit of vanishing external momenta. Taking the limit of m LQ m H and demanding the one loop correction to be smaller than the tree level coupling we obtain For degenerate masses it reduces to µ 1 4πm LQ , whereas in the case of a large splitting m heavy m light : µ 1 4πm heavy /(log m heavy /m light ) 1/2 . A similar result is obtained for µ 2 , changing S 1 byS 1 .
Large cubic couplings can also have an impact on the cross section of double Higgs production. In Fig. 15 we show diagrams at one and two loops with the largest number of cubic interactions. These integrals can be computed analytically in the limit of negligible Higgs mass and zero external momentum, leading to: I (2) /I (1) = (1/2)(µ 1 /4πmR 2 ) 2 , where the superindex labels the number of loops. We have also computed these integrals numerically for external momentum of order 500 GeV and different angles. For the regions of the parameter considered in the phenomenological analysis we find this ratio to be in general 1%, with a few points being of order 5-10%. A full two loop calculation of these corrections is beyond the scope of this work. We take cubic couplings below 1/3 of the bound given above.