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

TL;DR

This paper investigates how new physics, via dimension-six operators, could enhance double Higgs production at future electron-positron colliders, identifying key operators and estimating experimental sensitivities.

Contribution

It identifies two relevant effective operators affecting double Higgs production and evaluates their potential detectability at upcoming colliders.

Findings

Two key operators are not yet strongly constrained by existing data.

Sensitivity estimates for operator coefficients at various future collider scenarios.

Expected 95% confidence level limits for different energy and luminosity benchmarks.

Abstract

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.

Tables

Table 1. Table 1 : Table of the different benchmark scenarios considered in our analysis. Each benchmark consists of a specific value of the center of mass energy ( s ) 𝑠 (\sqrt{s}) and luminosity ( L ) 𝐿 (L) that has been proposed for the future e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} colliders. The last two columns represent the 95 % CL intervals for each operator coefficient taken individually in the analysis with k = 1 𝑘 1 k=1 ( k = 0.35 ) 𝑘 0.35 (k=0.35) .

Benchmark	Experiment	$\sqrt{s}$ (GeV)	$L$ (ab-1)	$\left|c_{e\phi }/{\Lambda }^2\right|\left({\mathrm{TeV}}^{-2}\right)$	$\left|c_{et}/{\Lambda }^2\right|\left({\mathrm{TeV}}^{-2}\right)$
1	FCC-ee	350	2.6
2	CLIC	380	0.5
3	ILC	500	4
4	CLIC	1500	1.5
5	CLIC	3000	3.0

Equations

$${\mathcal{L}}_{\mathrm{SMEFT}}={\mathcal{L}}_{\mathrm{SM}}+\sum _i\frac{c_i^{(n)}}{{\Lambda }^{n-4}}{\mathcal{O}}_i^{(n)}+\dots$$

$$\frac{c_{e\phi }}{{\Lambda }^2}({\phi }^{†}\phi -\frac{{\upsilon }^2}{2}){\bar{l}}_L\phi e_R+h.c.$$

$$\frac{c_{et}}{{\Lambda }^2}{ϵ}_{ij}{\bar{l}}_L^i{e}_R{\bar{q}}_L^j{t}_R+h.c.$$

$$\begin{gathered}\fmfgraph*(30,30)\fmfset{thin}{0.02pt}\fmfset{arrow_{l}en}{2mm}\fmfleft{i}\fmfright{o}\fmf{fermion,label=$t$,l.side=right,tension=0.7}{v,v}\fmf{fermion}{i,v,o}\fmffreeze\fmfv{label={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}$c_{et}$},l.angle=-90,l.dist=5thick,d.fill=full,fore=(0,,0,,150),d.shape=circle,d.size=7}{v}\fmflabel{$e^-$}{i}\fmflabel{$e^+$}{o}\end{gathered}\qquad\quad=\quad-i\Sigma_{e}\quad=\quad-i\frac{6}{\left(4\pi\right)^{2}}\frac{c_{et}}{\Lambda^{2}}m_{t}^{3}\left(1+\frac{1}{\bar{\epsilon}}+\log\frac{\mu^{2}}{m_{t}^{2}}\right)$$

$$/p-y_e\frac{\upsilon }{\sqrt{2}}-\delta y_e\frac{\upsilon }{\sqrt{2}}-{\Sigma }_e$$

$$\delta y_e=-\frac{6}{(4\pi {)}^2}\frac{\sqrt{2}}{\upsilon }\frac{c_{et}}{{\Lambda }^2}{m}_t^3\frac{1}{\overline{ϵ}}$$

$$m_e=y_e\frac{\upsilon }{\sqrt{2}}+\frac{6}{(4\pi {)}^2}\frac{c_{et}}{{\Lambda }^2}{m}_t^3\left(1+\log \frac{{\mu }^2}{m_t^2}\right)$$

$$\left|\frac{c_{et}}{{\Lambda }^2}\right|\lesssim \frac{8{\pi }^2}{3}\frac{m_e}{m_t^3}\simeq 2\times 10^{-3}{\mathrm{TeV}}^{-2}$$

$$y_e(\mu )=\frac{\sqrt{2}}{\upsilon }{m}_e-\frac{6}{(4\pi {)}^2}\frac{\sqrt{2}}{\upsilon }\frac{c_{et}}{{\Lambda }^2}{m}_t^3\left(1+\log \frac{{\mu }^2}{m_t^2}\right)$$

$$g_{\bar{e}eh}=\frac{c_{e\phi }{\upsilon }^2}{{\Lambda }^2\sqrt{2}}$$

$$\left|-\frac{m_e}{\upsilon }+\frac{c_{e\phi }(\mu ){\upsilon }^2}{{\Lambda }^2\sqrt{2}}-\frac{3}{(4\pi {)}^2}\frac{y_t}{\sqrt{2}}\frac{c_{et}}{{\Lambda }^2}(4m_t^2-m_h^2)\left[f(m_h^2,m_t^2)+\log \frac{{\mu }^2}{m_t^2}\right]\right|\lesssim 600\frac{m_e}{\upsilon }$$

$$g_{\bar{e}ehh}=\frac{3c_{e\phi }\upsilon }{2{\Lambda }^2\sqrt{2}}$$

$${\chi }^2={\chi }^2(\frac{c_{e\phi }}{{\Lambda }^2},\frac{c_{et}}{{\Lambda }^2})=\frac{{\left[\sigma (\frac{c_{e\phi }}{{\Lambda }^2},\frac{c_{et}}{{\Lambda }^2})-{\sigma }_{\mathrm{SM}}\right]}^2}{\delta {\sigma }^2}$$

$$\delta \sigma =\sqrt{\delta {\sigma }_{\mathrm{stat}}^2+\delta {\sigma }_{\mathrm{sys}}^2}=\sqrt{\frac{{\sigma }_{\mathrm{SM}}}{L}+{\alpha }^2{\sigma }_{\mathrm{SM}}^2}.$$

$$-\frac{m_e}{\upsilon }\to -\frac{m_e}{\upsilon }+\frac{c_{e\phi }(\mu ){\upsilon }^2}{{\Lambda }^2\sqrt{2}}-\frac{3}{(4\pi {)}^2}\frac{y_t}{\sqrt{2}}\frac{c_{et}}{{\Lambda }^2}(4m_t^2-q^2)\left[f(q^2,m_t^2)+\log \frac{{\mu }^2}{m_t^2}\right]$$

$$f(q^2,m_t^2)=2+\sqrt{1-\frac{4m_t^2}{q^2}}\log \frac{2m_t^2-q^2+\sqrt{q^2(q^2-4m_t^2)}}{2m_t^2}$$

$${\mathcal{M}}^{(a)div}+{\mathcal{M}}^{(b)div}$$

$${\mathcal{M}}^{(c)ct}+{\mathcal{M}}^{(d)ct}+{\mathcal{M}}^{(e)ct}$$

$$\delta c_{e\phi }=\frac{6}{(4\pi {)}^2}{c}_{et}{y}_t\left(y_t^2-\lambda \right)\frac{1}{\overline{ϵ}}$$

$$\delta y_e=-\frac{3}{(4\pi {)}^2}{c}_{et}{\upsilon }^2{y}_t^3\frac{1}{\overline{ϵ}}$$

$$\mu \frac{\partial c_{e\phi }}{\partial \mu }=\frac{12}{(4\pi {)}^2}{c}_{et}{y}_t\left(y_t^2-\lambda \right)$$

Full text

I Introduction

The discovery of the Higgs boson during the Run I of the LHC Chatrchyan:2012ufa 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 (Khachatryan:2016vau, ; Aaboud:2018xdt, ).

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) CLIC , the International Linear Collider (ILC) ILC , the Circular Electron Positron Collider (CPEC) CPEC and the Future Circular Collider with $e^{+}e^{-}$ (FCC-ee) at CERN, previously known as TLEP fcc-ee .

The main production mechanism of the Higgs boson at $e^{+}e^{-}$ colliders is the bremsstrahlung process $e^{+}e^{-}\to 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^{-}\to W^{*}W^{*}/Z^{*}Z^{*}\to h\nu\bar{\nu}/he^{+}e^{-}$ which provide an increasingly powerful handle at higher center-of-mass energies. Finally, also the process $e^{+}e^{-}\to 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 Ellis:2017kfi ; Durieux:2017rsg ; DiVita:2017vrr .

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^{-}\rightarrow hhZ$) and vector boson fusion ($e^{+}e^{-}\rightarrow e^{+}e^{-}(\nu\bar{\nu})hh$) processes Boudjema:1995cb ; Haba:2013xla ; Panico:2015cst . The top Yukawa can be measured in double Higgs production in association with top quarks ($e^{+}e^{-}\rightarrow hht\bar{t}$). 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^{-}\rightarrow 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 Djouadi:1996hp ; LopezVal:2009qy .

At hadron colliders, double Higgs production via gluon fusion at LHC has been exhaustively studied as a probe of physics beyond the SM hh-pp . The sensitivity to new physics is enhanced due to a cancellation between triangle and box contributions in the gluon fusion process in the SM cancellation . It is well known that Higgs pair production at hadron colliders is sensitive to new physics effects parametrized by higher-dimensional operators Goertz:2014qta ; effective . On the other hand, an enhancement in the cross section can also arise from the presence of an hidden sector, as studied in Oliveira:2010uv . Double Higgs production has also been studied as probe for Higgs anomalous couplings at future electron-proton colliders Kumar:2015kca .

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^{-}\rightarrow 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^{-}\rightarrow 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.

II SM double Higgs production at $e^{+}e^{-}$ colliders

The process $e^{+}e^{-}\to hh$ is an interesting one from the theoretical point of view because SM tree level diagrams (see Fig.  1) give a negligible contribution to the cross section since they are proportional to $m_{e}/\upsilon$, where $m_{e}$ is the electron mass and $\upsilon=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 Gaemers:1984vw ; Djouadi:1996hp ; LopezVillarejo:2008xw .

Non-negligible contributions to $e^{+}e^{-}\to hh$ can therefore only come from one-loop diagrams. In the SM, all one-loop diagrams involving the $\bar{e}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 $\gamma$ and $Z$ boson which give rise to two Higgs bosons, also vanish (see Fig. 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 $\bar{e}eh$ vertex (when one Higgs is taken to its vev) and hence negligible (see Fig. 2 (a) and (c)).

Therefore, the only contribution to Higgs pair production in the SM comes from W and Z box diagrams of Fig. 3. Notice that, contrary to double Higgs production in gluon fusion $gg\to 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 $\lambda$ is also negligible because it enters only in diagrams that vanish in the $m_{e}=0$ limit (see Fig. 2 (c)) where the triangle loop is related to the renormalization of the $\bar{e}eh$ coupling.

The energy dependence of the leading order SM cross section for $e^{+}e^{-}\to hh$ is shown in Fig. 4. The cross section acquires its maximum value of approximately $0.015$ fb at around $\sqrt{s}=500$ GeV.

III EFT contributions to $e^{+}e^{-}\to 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 Fig. 4) as discussed in the previous section. However, with large luminosities expected at future $e^{+}e^{-}$ colliders, a few hundred events might eventually be collected in the course of a few years, allowing for the experimental study of this final state. 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 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

[TABLE]

where $\Lambda$ is the mass scale of new physics, $c_{i}^{(n)}$ are dimensionless coefficients and $n$ is the dimension of the gauge invariant operators ${\cal O}_{i}^{(n)}$ 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/\Lambda$, where $E$ is the typical energy of the process (the unique dimension-five operator does not contribute to the process $e^{+}e^{-}\to hh$). In this work we use the parametrization of Grzadkowski:2010es . 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 $\bar{e}eZ$, $e\nu W$, $hZZ$ and $hWW$ vertices are already (strongly) constrained by electroweak precision data and LHC Higgs measurements Han:2004az ; deBlas:2016ojx ; deBlas:2017wmn ; Ellis:2014jta ; Englert:2015hrx ; Butter:2016cvz ; Ellis:2018gqa 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 $\bar{e}ehh$ coupling and operators that generate an effective $\bar{e}e\bar{t}t$ 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

[TABLE]

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

[TABLE]

In the equations above $c_{e\varphi}$ and $c_{et}$ are dimensionless coefficients, $\Lambda$ is the scale of new physics, $l=(\nu\,\,e)$, $q=(t\,\,b)$, $\varphi$ is the Higgs doublet and $\epsilon_{ij}$ is the total antisymmetric tensor of rank 2.

The operator in Eq. (2) has been written with the constant piece $\upsilon^{2}/2$ subtracted to the invariant $\varphi^{\dagger}\varphi$ term in order to formally maintain the tree level relation $m_{e}=y_{e}\upsilon/\sqrt{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). The contribution of this effective operator to the electron self energy in dimensional regularization is given by {fmffile}selfelectron

[TABLE]

where $1/\bar{\epsilon}=1/\epsilon-\gamma+\log 4\pi$. Thus the inverse electron propagator reads

[TABLE]

In $\overline{\rm MS}$ the Yukawa counterterm is chosen to be

[TABLE]

such that the physical electron mass is given by111If one had used an on-shell scheme, then the electron mass definition would have remained unchanged while the Yukawa counterterm would have been modified including also the finite and $\mu$-dependent piece. In the end the two schemes give the same result, as it should be.

[TABLE]

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 $|\delta m_{e}|\lesssim m_{e}$. In this case we have that

[TABLE]

By inverting Eq. (7) it is possible to express the relation between the Yukawa coupling and the $c_{et}$ coefficient as follows

[TABLE]

Therefore, thanks to this relation, tree level diagrams of Fig. 1 proportional to $y_{e}$ are not negligible anymore if $c_{et}\neq 0$. Notice from eq. (9) that, contrary to the SM case, the limit of vanishing electron mass does not imply a vanishing Yukawa coupling. The scale $\mu$ entering in Eq. (9) will be set equal to $2m_{h}$ in the computation of $e^{+}e^{-}\to hh$.

The operator in Eq. (2) introduces a tree level coupling of the electron to the Higgs given by

[TABLE]

After considering all contributions to the $\bar{e}eh$ vertex, it is possible to show that the recent upper bound on the electron Yukawa coupling obtained from Higgs decay Altmannshofer:2015qra $y_{e}<600\,y_{e}^{\rm SM}$ implies that

[TABLE]

where $f(m_{h}^{2},m_{t}^{2})$ is given in Appendix A. The operator in Eq. (2), besides modifying the $\bar{e}eh$ vertex, induces also an effective $\bar{e}ehh$ coupling given by

[TABLE]

which is not present in the SM. This operator contributes at tree level to $e^{+}e^{-}\to hh$, as shown in Fig. 5.

On the other hand, the operator in Eq. (3) contributes to $e^{+}e^{-}\to hh$ through the counterterm related to the redefinition of the Yukawa coupling of eq. (9) and it also enters directly at one loop, as shown by the diagrams of Fig. 6.

Notice that the operator in Eq. (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) and its coefficient $c_{e\varphi}$ has to be formally taken as function of the renormalization scale $\mu$. For the explicit derivation of the counterterm see Appendix B. Therefore, in the process we are studying the coefficients $c_{e\varphi}$ and $c_{et}$ are both formally evaluated at the scale $\mu=2m_{h}$.

In our computation we consider just the leading contributions of the operator of Eq. (2) which arise at tree level while the contributions of the operator of Eq. (3) comes at one loop. The total cross section turns out to be a pure quadratic function of the coefficients $c_{e\varphi}$, $c_{et}$, namely the only sizable new physics contributions are of order $c_{e\varphi}^{2}$, $c_{et}^{2}$ and $c_{e\varphi}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}/\upsilon$ and therefore negligible. Helicity selection rules and non-interference effects in the context of dimension-six operators have been studied in Azatov:2016sqh . 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/\Lambda^{4}$ expansion are proportional to $m_{e}/\upsilon$ as well.

IV Analysis and results

We compute the $e^{+}e^{-}\to hh$ cross section $\sigma=\sigma(\frac{c_{e\varphi}}{\Lambda^{2}},\frac{c_{et}}{\Lambda^{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 Alloul:2013bka and generate the corresponding FeynArts model output. We used FeynArts 3.10 FeynArts and FormCalc 8.4  FormCalc to compute the tree and one-loop amplitudes relevant for the process in the chiral limit ($m_{e}=0$). Finally, we use LoopTools Hahn:1998yk 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 Degrande:2014vpa and MadGraph5_aMC@NLO Alwall:2014hca .

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 $\chi^{2}$ function

[TABLE]

where $\sigma_{{\rm SM}}=\sigma(0,0)$. The total cross section uncertainty $\delta\sigma$ that enters in the $\chi^{2}$ computation is given by the combination of the expected experimental $\delta\sigma_{\rm exp}$ and theoretical uncertainties $\delta\sigma_{\rm th}$. In our analysis we assume the theoretical uncertainty to be negligible such that the total uncertainty coincides with the expected experimental one, namely $\delta\sigma=\delta\sigma_{\rm exp}$, which is given by the sum in quadrature of statistic $\delta\sigma_{\rm stat}$ and systematic uncertainties $\delta\sigma_{\rm sys}$

[TABLE]

The statistical uncertainty is taken to be $\delta\sigma_{\rm stat}=\sqrt{\sigma_{\rm SM}/L}$, where $L$ is the integrated luminosity. The systematic uncertainty has been parametrized by $\delta\sigma_{\rm sys}=\alpha\,\sigma_{\rm SM}$, in analogy to the study performed in Durieux:2017rsg , where $\alpha$ 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 $\alpha=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 $\chi^{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 $(f\bar{f})$ and the corresponding branching ratio ${\rm BR}(h\to f\bar{f})$ have been identified, then we need to properly rescale the cross section and uncertainty that enter in the chi-squared function of Eq. (13) by a factor $k={\rm BR}(h\to f_{1}\bar{f}_{1})\times{\rm BR}(h\to f_{2}\bar{f}_{2})$. For instance, if we assume that each Higgs particle is going to be reconstructed only through its decay to $b\bar{b}$ then $k\sim 0.35$.

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\varphi}/\Lambda^{2}$, which is $|c_{e\varphi}/\Lambda^{2}|\lesssim 3\times 10^{-3}$ TeV*-2* for $k=1$. This behaviour is expected since the contribution to the total cross section of the operator in Eq. (2) is almost insensitive to the energy in the process. Assuming an order one coefficient for $c_{e\varphi}$ this implies a quite strong bound on the new physics scale of the order $\Lambda\gtrsim 18$ TeV.

On the other hand, the bound on the coefficient $c_{et}/\Lambda^{2}$ turns out to be weaker than the bound on $c_{e\varphi}/\Lambda^{2}$. This is expected since the $c_{et}/\Lambda^{2}$ contribution enters at one-loop compared to $c_{e\varphi}/\Lambda^{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. (3) turns out to be quite sensitive of the energy of the process. Assuming an order one coefficient for $c_{et}$ and $k=1$, the weakest bound $|c_{et}/\Lambda^{2}|\lesssim 0.15$ TeV*-2* is obtained from one of the benchmark configurations with lowest center of mass energy and luminosity and can be translated into $\Lambda\gtrsim 2.5$ TeV, while the strongest bound $|c_{et}/\Lambda^{2}|\lesssim 0.01$ TeV*-2* is obtained from the benchmark configuration with highest center of mass energy and can be translated into $\Lambda\gtrsim 10$ TeV. The case $k=0.35$ shows modifications of the bounds of the order of 25-50 % with respect to the $k=1$ case. The fine tuning bound on $c_{et}/\Lambda^{2}$ in Eq. (8) is one order of magnitude stronger than the best expected bound coming from our 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 Fig. 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 Fig.  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 Fig 7 we do not show the bound of Eq. (11) coming from Higgs decay since it is weaker than the expected bounds we obtained from $e^{+}e^{-}\to hh$ and the exclusion region would lie outside the range of the plot.

The results for $k=0.35$ are not presented since they differ from the results in Fig. 7 by $\sim$30% and the corresponding ellipses do not present significant modifications.

V 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\varphi}$ 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\varphi}$ typically probe scales of ${\cal O}(10\;\mbox{TeV})$ while the $c_{et}$ operator is less constrained since it enters only at one-loop level (of course, more stringent limits on $c_{et}/\Lambda^{2}$ of ${\cal O}(10^{-3})$ TeV*-2* can be obtained by studying top quark pair production at future $e^{+}e^{-}$ colliders, as shown in Durieux:2018tev ). In conclusion, searches for $e^{+}e^{-}\rightarrow 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.

Acknowledgments

AT would like to thank ICTP-SAIFR and IFT-UNESP for hospitality. The work of AV was supported by a CAPES MSc fellowship. RR is partially supported through a CNPq grant 307925/2017-0 and a Fapesp 2016/01343-7 grant. AV and CD are supported by the Fund for Scientific Research F.N.R.S. through the F.6001.19 convention. We thank E. Pontón for a valuable discussion about the correction to the electron mass induced by higher dimensional operators and the renormalization procedure. AV would like to thank F. Maltoni for discussion.

Appendix A $\bar{e}eh$ coupling modification

The effective operators in Eq. (2) and (3) modify the $\bar{e}eh$ coupling with respect to the SM case as follows

[TABLE]

where $q$ is the Higgs momentum and

[TABLE]

Eq. (15) has been obtained by taking into account the tree level contribution to the $\bar{e}eh$ vertex coming from the $c_{e\varphi}$ operator, the redefinition of the Yukawa coupling in eq. (9), the top-loop diagram induced by the $c_{et}$ operator and the proper counterterms.

Appendix B Divergent $e^{+}e^{-}\to hh$ diagrams and counterterms

The one-loop diagrams (b) and (c) of Fig. 6 are proportional to $c_{et}$ and UV divergent. The computation in dimensional regularization of the divergent part of these diagrams gives

[TABLE]

where $y_{t}$ is the top Yukawa, $\lambda$ the Higgs self-coupling, $u_{e}$ and $v_{e}$ are the electron and positron Dirac spinors and $1/\bar{\epsilon}=1/\epsilon-\gamma+\log 4\pi$. Let us now consider the counterterm diagrams proportional to $\delta c_{e\varphi}$ and $\delta y_{e}$ needed to cancel this divergence, as shown in diagrams (c), (d) and (e) of Fig. 8. We have that

[TABLE]

By comparing eq. (17) and (18) we obtain the explicit form of the counterterms in $\overline{\rm MS}$

[TABLE]

[TABLE]

where we have used $m_{t}=y_{t}\upsilon/\sqrt{2}$ and $m_{h}^{2}=2\lambda\upsilon^{2}$. Notice that the counterterm in Eq. (20) coincides with the one derived in Eq. (6), as it should be. From the explicit form of the counterterm $\delta c_{e\varphi}$ we can read off the contribution of $c_{et}$ to the RG equation of $c_{e\varphi}$

[TABLE]

which agrees with Jenkins:2013zja and Jenkins:2013wua .