New Physics in Double Higgs Production at Future $e^+ e^-$ Colliders

We study the effects of new physics in double Higgs production at future $e^+ e^-$ colliders. In the Standard Model the chiral limit ($m_e=0$) plays an important role for this process, being responsible for the smallness of the tree-level diagrams with respect to the 1-loop contributions. In our work, we consider the possibility of an enhancement due to the contribution of Standard Model dimension-six effective operators. We show that there are only two relevant operators for this process that are not yet (strongly) constrained by other data. We perform a sensitivity study on the operator coefficients for several benchmark values of energy and integrated luminosity related to the proposed linear colliders such as CLIC, ILC and FCC-ee and we derive expected 95% CL limits for each benchmark scenario.


Introduction
The discovery of the Higgs boson during the Run I of the LHC [1,2] has put in place the main building block that was missing for the experimental validation of the Standard Model (SM). Since then, a great effort has been made by the experimental collaborations in the attempt to define the properties of this new particle, namely its mass, spin, parity and coupling to itself and the other particles of the SM, mainly through global fits in the so-called kappa framework [3,4].
These analyses are crucial to pin down the Higgs boson properties and to understand the nature of electroweak symmetry breaking (EWSB). These studies will play a fundamental role especially during the high-luminosity run of the LHC (HL-LHC) as well as for the future hadron and lepton colliders. Any deviations from SM predictions would unravel the presence of new physics.
As the second run of the LHC is coming to an end, no clear signs of new physics have been found yet. This fact points to a scenario in which new physics is most probably out of the reach of the LHC and in this case the best way to search for it is through indirect effects via precision measurements.
Precision studies of the properties of the Higgs boson and the nature of the electroweak symmetry breaking strongly motivate the construction of a lepton collider which benefits from a cleaner environment with respect to hadron colliders. There have been several proposals for a future electron-positron collider, such as the Compact Linear Collider (CLIC) [5], the International Linear Collider (ILC) [6], the Circular Electron Positron Collider (CPEC) [7] and the Future Circular Collider with e + e − (FCC-ee) at CERN, previously known as TLEP [8].

JHEP05(2019)020
The main production mechanism of the Higgs boson at e + e − colliders is the bremsstrahlung process e + e − → hZ (Higgsstrahlung). At a center-of-mass energies of 240-250 GeV, close to the maximum of the Higgsstrahlung cross section, this process will allow to determine Higgs couplings to gauge bosons with unprecedented precision. In addition there are also weak boson fusion production processes e + e − → W * W * /Z * Z * → hνν/he + e − which provide an increasingly powerful handle at higher center-of-mass energies. Finally, also the process e + e − → tth benefits from high energies and represent an important measurement to directly constrain the top Yukawa coupling. A comprehensive sensitivity study about the effect of new physics, parametrized by higher-dimensional operators, affecting these production mechanisms for the different proposed e + e − machines have been performed in [9][10][11].
Some of the Higgs boson couplings can also be tested in higher order processes involving for instance Higgs pair production. In this case, the Higgs self coupling and the couplings to gauge bosons can be measured in the so-called double higgstrahlung (e + e − → hhZ) and vector boson fusion (e + e − → e + e − (νν)hh) processes [12][13][14]. The top Yukawa can be measured in double Higgs production in association with top quarks (e + e − → hhtt). These processes are tree-level dominated processes but compared to the previous ones they are characterized by higher orders in the coupling constants.
On the other hand, the process e + e − → hh, where only two Higgs bosons are actually produced in the final state, is completely dominated by the contribution of one-loop diagrams and therefore one can test higher order effects in a clean way because they are not masked by tree diagram contributions. For instance, it can be useful to discriminate between the Higgs sector of the Standard Model from the more complicated scalar sectors belonging to possible extensions, e.g. two Higgs doublet model [15,16].
At hadron colliders, double Higgs production via gluon fusion at LHC has been exhaustively studied as a probe of physics beyond the SM [14,[17][18][19][20]. The sensitivity to new physics is enhanced due to a cancellation between triangle and box contributions in the gluon fusion process in the SM [21,22]. It is well known that Higgs pair production at hadron colliders is sensitive to new physics effects parametrized by higherdimensional operators [23,24]. On the other hand, an enhancement in the cross section can also arise from the presence of an hidden sector, as studied in [25]. Double Higgs production has also been studied as probe for Higgs anomalous couplings at future electron-proton colliders [26].
The SM cross section for double Higgs production at the LHC is not very large (approximately 37 fb at 14 TeV at NNLO) and the background can be challenging even for HL-LHC. Therefore, the cleaner environment of an electron-positron collider could be very helpful to find deviations from the SM or to improve bounds on new physics.
In this work we will proceed in that direction and focus on the process e + e − → hh at future lepton colliders as a probe of new physics which we take to be parametrized by the presence of dimension-six effective operators of the SM effective field theory (SMEFT). This paper is organized as follows: in section II we revise the SM computation for the process e + e − → hh. In section III we highlight the relevant SMEFT contributions that we consider in our study. In section IV we discuss different benchmark scenarios for future e + e − colliders, we present the analysis strategy and we report the 95% CL bounds on the operator coefficients for each benchmark scenario. In section V we conclude.

SM double Higgs production at e + e − colliders
The process e + e − → hh is an interesting one from the theoretical point of view because SM tree level diagrams (see figure 1) give a negligible contribution to the cross section since they are proportional to m e /υ, where m e is the electron mass and υ = 246 GeV is the Higgs vacuum expectation value (VEV). This fact has been recognized long ago and as a consequence the cross section is quite small both in the SM and MSSM extensions [15,27,28]. Non-negligible contributions to e + e − → hh can therefore only come from one-loop diagrams. In the SM, all one-loop diagrams involving theēeh vertex must give zero contributions in the chiral limit m e = 0, to all orders in perturbation theory. Furthermore, because of CP invariance, the diagrams containing intermediate γ and Z boson which give rise to two Higgs bosons, also vanish (see figure 2 (b)). Additional contributions from triangle diagrams involving the quartic W + W − hh/ZZhh and triple hhh couplings are also related to the renormalization of theēeh vertex (when one Higgs is taken to its vev) and hence negligible (see figure 2 (a) and (c)).
Therefore, the only contribution to Higgs pair production in the SM comes from W and Z box diagrams of figure 3. Notice that, contrary to double Higgs production in gluon fusion gg → hh, there is no such feature as the cancellation between triangle and box diagrams because the triangle ones are subleading and vanish in the m e = 0 limit. Moreover, the dependence of the SM cross section on the triple Higgs coupling λ is also negligible because  it enters only in diagrams that vanish in the m e = 0 limit (see figure 2 (c)) where the triangle loop is related to the renormalization of theēeh coupling. The energy dependence of the leading order SM cross section for e + e − → hh is shown in figure 4. The cross section acquires its maximum value of approximately 0.015 fb at around √ s = 500 GeV.

EFT contributions to e + e − → hh
Double Higgs production at e + e − colliders in the SM has been shown to have a tiny cross section of the order of fraction of femtobarns (see figure 4) as discussed in the previous section. With large luminosities expected at future e + e − colliders, order one hundred events might eventually be collected in the course of a few years. We acknowledge that this low number of events may force realistic analysis to face experimental challenges related to background discrimination and b-tagging efficiencies. The discussion of these issues is beyond the scope of this paper. On the other hand, cross sections can be enhanced by contributions coming from physics beyond the SM and in this paper we want to entertain this possibility. In particular JHEP05(2019)020 we will consider effects of new physics parametrized by the presence of higher dimensional operators in the SMEFT framework. The general SMEFT lagrangian can be written as where Λ is the mass scale of new physics, c (n) i are dimensionless coefficients and n is the dimension of the gauge invariant operators O (n) i built up with SM fields. It allows for a systematic study of deviations from the SM while respecting established symmetry principles.
In this work we focus on the contributions of dimension-six operators of the SMEFT because they give the leading contributions in the systematic expansion E/Λ, where E is the typical energy of the process (the unique dimension-five operator does not contribute to the process e + e − → hh). In this work we use the parametrization of [29]. In principle, all dimension-six operators that are relevant for the electron and Higgs sector should be considered. However, several of these operators are already constrained from other observables and therefore will not be taken into account in this study. In particular, dimension-six operators that modify theēeZ, eνW , hZZ and hW W vertices are already (strongly) constrained by electroweak precision data and LHC Higgs measurements [30][31][32][33][34][35][36] and we will safely ignore their effects. We are then left with two classes of effective operators that can give sizable contributions: operators that induce an effectiveēehh coupling and operators that generate an effectiveēett coupling. The first class enters at tree-level while the second class operators only contribute at one-loop.
There is a unique operator belonging to the first class On the other hand there are seven four-fermion operators belonging to the second class, however six of them give zero contribution because of their chirality structure and in the end we are left with just one four-fermion operator In the equations above c eϕ and c et are dimensionless coefficients, Λ is the scale of new physics, l = (ν e), q = (t b), ϕ is the Higgs doublet and ij is the total antisymmetric tensor of rank 2.
The operator in eq. (3.2) has been written with the constant piece υ 2 /2 subtracted to the invariant ϕ † ϕ term in order to formally maintain the tree level relation m e = y e υ/ √ 2 also in the effective theory. This mass relation is however altered by the potentially sizable loop correction to the electron mass coming from the top-quark loop induced by the effective operator in eq. (3.3). The contribution of this effective operator to the electron self energy in dimensional regularization is given by t

JHEP05(2019)020
where 1/¯ = 1/ − γ + log 4π. Thus the inverse electron propagator reads In MS the Yukawa counterterm is chosen to be such that the physical electron mass is given by 1 From the theoretical point of view this mass correction may introduce a fine tuning problem and in order to avoid it one must require that |δm e | m e . In this case we have that By inverting eq. (3.7) it is possible to express the relation between the Yukawa coupling and the c et coefficient as follows Therefore, thanks to this relation, tree level diagrams of figure 1 proportional to y e are not negligible anymore if c et = 0. Notice from eq. (3.9) that, contrary to the SM case, the limit of vanishing electron mass does not imply a vanishing Yukawa coupling. The scale µ entering in eq. (3.9) will be set equal to 2m h in the computation of e + e − → hh. The operator in eq. (3.2) introduces a tree level coupling of the electron to the Higgs given by After taking into account all contributions to theēeh vertex, it is possible to consider the Higgs decay to electrons [37] to extract a direct bound on our effective operator coefficients. We have where the explicit form of f (m 2 h , m 2 t ) is given in appendix A. Recent LHC measurements provide a coefficient κ LHC e 600. This upper bound is expected to be improved at future lepton colliders by roughly two orders of magnitude (κ e = 10), as shown in [37]. The operator in eq. (3.2), besides modifying theēeh vertex, induces also an effectiveēehh coupling given by which is not present in the SM. This operator contributes at tree level to e + e − → hh, as shown in figure 5. On the other hand, the operator in eq. (3.3) contributes to e + e − → hh through the counterterm related to the redefinition of the Yukawa coupling of eq. (3.9) and it also enters directly at one loop, as shown by the diagrams of figure 6. Notice that the operator in eq. (3.2) plays also the role of the counterterm needed to absorb the divergence produced by the one-loop insertion of the operator in eq. (3.3) and its coefficient c eϕ has to be formally taken as function of the renormalization scale µ. For the explicit derivation of the counterterm see appendix B. Therefore, in the process we are studying the coefficients c eϕ and c et are both formally evaluated at the scale µ = 2m h .
In our computation we consider just the leading contributions of the operator of eq. (3.2) which arise at tree level while the contributions of the operator of eq. (3.3) comes at one loop. The total cross section turns out to be a pure quadratic function of the coefficients c eϕ , c et , namely the only sizable new physics contributions are of order c 2 eϕ , c 2 et and c eϕ c et because linear terms coming from the interference between SM diagrams, which are helicity conserving, and new physics diagrams, which are helicity flipping, turns out to be proportional to m e /υ and therefore negligible. Helicity selection rules and non-interference effects in the context of dimension-six operators have been studied in [38]. Notice that the EFT expansion is under control because possible interference terms expected from dimension-eight operators which would give comparable contribution in term of the 1/Λ 4 expansion are proportional to m e /υ as well.

JHEP05(2019)020 4 Analysis and results
We compute the e + e − → hh cross section σ = σ( ceϕ Λ 2 , cet Λ 2 ) as function of the effective couplings as discussed in the previous section. In order to perform this calculation, we first implemented the effective lagrangian in FeynRules [39] and generate the corresponding FeynArts model output. We used FeynArts 3.10 [40] and FormCalc 8.4 [41] to compute the tree and one-loop amplitudes relevant for the process in the chiral limit (m e = 0). Finally, we use LoopTools [41] to compute numerically the cross section as a function of the center of mass energy and effective couplings. We further checked the cross section computation by means of the development version of NLOCT [42] and MadGraph5 aMC@NLO [43].
In order to extract the expected 95% CL limits on the effective operators couplings we assume the measured cross section to coincide with the SM predictions and we construct the following χ 2 function where σ SM = σ(0, 0). The total cross section uncertainty δσ that enters in the χ 2 computation is given by the combination of the expected experimental δσ exp and theoretical uncertainties δσ th . In our analysis we assume the theoretical uncertainty to be negligible such that the total uncertainty coincides with the expected experimental one, namely δσ = δσ exp , which is given by the sum in quadrature of statistic δσ stat and systematic uncertainties δσ sys The statistical uncertainty is taken to be δσ stat = σ SM /L, where L is the integrated luminosity. The systematic uncertainty has been parametrized by δσ sys = α σ SM , in analogy to the study performed in [10], where α is a dimensionless coefficient that represents the magnitude of the systematic error in relation to the SM cross section. We take a conservative value and we fix α = 0.1, which corresponds to a 10% error. However, the impact of the systematic uncertainty will be marginal since our total uncertainty turns out to be statistics dominated due to the expected smallness of the SM cross sections. We consider different benchmark values of the center of mass energy and luminosity that have been proposed for the future e + e − machines (see table 1) and for each configuration we determine 95% CL limits on the operator coefficients. Values of the coefficients for which χ 2 > 3.84 are considered excluded.
To perform a more realistic investigation we have to consider a set of possible final states that are assumed to be measured at future e + e − colliders in order to reconstruct the Higgs particle through its decay channels. Once a set of final states (ff ) and the corresponding branching ratio BR(h → ff ) have been identified, then we need to properly rescale the cross section and uncertainty that enter in the chi-squared function of eq. (4.1) by a factor k = BR(h → f 1f1 ) × BR(h → f 2f2 ). For instance, if we assume that each Higgs particle is going to be reconstructed only through its decay to bb then k ∼ 0.35. 3.0 < 0.002 (< 0.002) < 0.012 (< 0.015) Table 1. Table of the different benchmark scenarios (B) considered in our analysis. Each benchmark consists of a specific value of the center of mass energy ( √ s) and luminosity (L) that has been proposed for the future e + e − colliders. The last two columns represent the 95 % CL intervals for each operator coefficient taken individually in the analysis with k = 1 (k = 0.35).

JHEP05(2019)020
The results for k = 1 (k = 0.35) in which each operator is considered individually are reported in the last two columns of table 1. The table shows that all benchmark configurations considered in our study provide the same order of magnitude bound for the coefficient c eϕ /Λ 2 , which is |c eϕ /Λ 2 | 3 × 10 −3 TeV −2 for k = 1. This behaviour is expected since the contribution to the total cross section of the operator in eq. (3.2) is almost insensitive to the energy in the process. Assuming an order one coefficient for c eϕ implies a quite strong bound on the new physics scale of the order Λ 18 TeV. It is important to notice that this scale changes if other assumptions about the UV theory responsible for this kind of NP are made. For instance, in the context of minimal flavor violation models one expects c eϕ to be naturally suppressed by the small electron Yukawa, namely c eϕ =c eϕ y e withc eϕ ∼ O (1). In this case the bound on the new physics scale reduces to Λ 30 GeV, which is clearly outside the validity of EFT.
Assuming flavor universality in the dimension-six lepton Yukawa operators, namley c eϕ = c µϕ = c τ ϕ , it is possible to use recent LHC Higgs measurements to derive a bound on c eϕ in this more constrained scenario [36]. In this case we have which is weaker than what we obtained from our analysis. On the other hand, the bound on the coefficient c et /Λ 2 turns out to be weaker than the bound on c eϕ /Λ 2 . This is expected since the c et /Λ 2 contribution enters at one-loop compared to c eϕ /Λ 2 which enters at tree level. Moreover, the bound depends on the benchmark configuration considered, because the contribution to the total cross section of the operator in eq. analysis, however one has to keep in mind that the fine tuning bound is based on theoretical considerations while our bound is based on experimental measurements. Moreover, the actual fine tuning could be milder thanks to cancellations induced by additional operators that we are not considering in our study.
The results for k = 1 in which both effective operator coefficients are taken into account are shown in figure 7. For each benchmark configuration, the exclusion region is represented by an ellipse. Points that lie outside the ellipse are considered excluded at 95% CL. By inspection of figure 7, we can infer that the best sensitivity is given by benchmark scenario number 5 which is characterized by the highest, among the considered configurations, center of mass energy of 3000 GeV. In figure 7 we show also the expected exclusion region (light pink strip) obtained from Higgs decay measurements at future lepton colliders in eq. (3.11). This contribution is complementary and allows to improve the bounds we obtain just by considering e + e − → hh.

JHEP05(2019)020
The results for k = 0.35 are not presented since they differ from the results in figure 7 by ∼30% and the corresponding ellipses do not present significant modifications.

Conclusions
Double Higgs production at future e + e − colliders offers the possibility to explore the sensitivity to dimension-6 operators involving electrons that have not been constrained yet. The small SM cross section and the clean environment make this process an ideal laboratory for these studies. In particular, two operators are relevant for this process and are characterized by dimensionless Wilson coefficients c eϕ and c et . By including their contributions to the double Higgs cross section we derived 95% bounds based on several benchmarks for these future colliders under certain assumptions of final decay channels to be reconstructed and the errors. We found that the bounds on c eϕ typically probe scales of O(10 TeV) while the c et operator is less constrained since it enters only at one-loop level. More stringent limits on c et /Λ 2 of O(10 −3 ) TeV −2 can be obtained by studying top quark pair production at future e + e − colliders, as shown in [44], and using that result would give marginalized bounds on c eϕ of the same order than the ones obtained in table 1. In conclusion, searches for e + e − → hh should also be pursued in addition to the more traditional double Higgs production in double higgstrahlung and vector boson fusion in order to explore these possible new couplings. which agrees with [45] and [46].