Discovery potential of Higgs boson pair production through 4$\ell$+$E\!\!/$ final states at a 100 TeV collider

We explore the discovery potential of Higgs pair production at a 100 TeV collider via full leptonic mode. The same mode can be explored at the LHC when Higgs pair production is enhanced by new physics. We examine two types of fully leptonic final states and propose a partial reconstruction method. The reconstruction method can reconstruct some kinematic observables. It is found that the $m_{T2}$ variable determined by this reconstruction method and the reconstructed visible Higgs mass are important and crucial to discriminate the signal and background events. It is also noticed that a new variable, denoted as $\Delta m$ which is defined as the mass difference of two possible combinations, is very useful as a discriminant. We also investigate the interplay between the direct measurements of $t\bar{t} h$ couplings and other related couplings and trilinear Higgs coupling at hadron colliders and electron-positron colliders.


I. INTRODUCTION
The discovery of Higgs boson at the LHC has motivated the high energy community to think of the next generation p p colliders. A 100 TeV collider can offer us a huge potential to probe various new physics [1]. For example, new vector boson W and Z can be discovered up to [25][26][27][28][29][30][31][32][33][34][35] TeV [2]. A heavy Higgs bosons of the two Higgs doublet model can be probed up to 20 TeV or so via single associate production [3]. In the simplified model, the superpartners of top quark and the gluino can be probed up to 5 TeV and 10 TeV, respectively [4], which can make decisive evidences on the fate of the electroweak supersymmetry models. Dark matter candidate can be probed up to 10 TeV or higher [5][6][7][8]. A 100 TeV collider could also perform high precision measurements on Higgs properties [9], top quark properties, EW physics, and so on.
Among various new physics candidates, the shape of Higgs potential plays a very special role. As we know, the shape of Higgs potential is determined by Higgs fields and Higgs self-couplings, new Higgs fields and self-couplings exist in most of extensions of the standard model (SM). Therefore, it is well-known [10][11][12][13][14][15] that to probe Higgs self-couplings at colliders can offer us a way to understand the nature of Higgs bosons, to reconstruct the shape of Higgs potentials, to understand the mystery of electroweak symmetry breaking. These couplings could play crucial parts in the EW baryogenesis scenarios [16,17]. For example, they are crucial to determine whether CP violation is strong enough to produce a large enough matter anti-matter asymmetry, which is needed in terms of the one of three Sakharov criteria. In the two Higgs doublet model, it is possible to introduce meaningful complex Higgs self couplings which can induce a large enough CP violation which is needed for the EW baryogensis scenarios. These couplings are also important to determine whether the strong first-order phase transition could occur for a realistic EW baryogenesis scenarios. These couplings can affect the gravitational wave radiation in the process of bubble collisions [18][19][20], while the gravitational wave can induced a B mode which is detectable from the cosmological microwave background [21][22][23]. So it is well-motived to explore the shape of Higgs potential in our world.
Compared to the standard model, new physics could modify either effective trilinear (or cubic) or quartic couplings or both, either in 10 − 20% via loop corrections [24] or more than 100% − 300% via tree-level corrections (say adding a dimension-6 operator [25][26][27] or many higher dimensional operators [28]). The Lorentz structure of triple Higgs boson can even be modified in the Higgs-Gravity model [29,30], which could lead to energetic Higgs bosons in the final states [31]. Using the discovered Higgs boson as a probe, to measure the self-couplings of Higgs boson could help us to further understand the nature of the Higgs bosons and to extract the information of the shape of Higgs potentials which encodes the electroweak symmetry breaking. Therefore to measure trilinear couplings and quartic couplings [32][33][34][35] will be of great importance and could be one of the prime targets for both the LHC high luminosity runs and future collider projects.
The study on the di-Higgs boson final states in the SM and new physics models at hadronic colliders has been being a hot topic recently, various production processes and final states have been explored in literatures, bbγγ [31,36], W W bb [37], W W γγ [38], bbτ τ [36], and rare decay final states 3 2j [39] and others [40]. Both ATLAS and CMS collaborations of the LHC had performed realistic simulation and analysis on di-Higgs boson final states [41,42].
In this work, we extend our study in [39] to the pure leptonic mode, i.e. pp → hh → 4 + / E in a 100 TeV collider. To our best knowledge, this mode has not been carefully studied in literatures due to its tiny production rate in the SM at the collision energy of the LHC.
But for some new physics models, the production rate of di-Higgs can be enhanced by a factor from 10 to 100, then this mode could be accessible even at the LHC. For a 100 TeV collision, the production rate of this mode in the SM itself is large enough and is accessible.
Meanwhile, since it is pure leptonic final states, this mode can be searched by experimental groups relatively easy. Therefore, it is meaningful to perform a careful analysis on this mode either for the LHC runs or for a future 100 TeV collider project.
In order to determine the Higgs self-couplings at future hadron colliders, a precision measurement of top Yukawa coupling at the LHC and future colliders is very crucial. The top quark Yukawa coupling plays a remarkable role to help us to probe the properties of Higgs boson. It is the strongest Yukawa coupling, which almost saturates the perturbation bounds; it can affect the vacuum stability [43] much seriously than any other Yukawa couplings in the SM; it determines the multi-Higgs production at hadron colliders and affects the decay of Higgs boson to gluon pair, di-photon and Zγ final states much larger than the other fermions in the SM; it affects the Higgs self-coupling measurements at hadron colliders, both trilinear and quartic coupling measurements.
Therefore, in this work, to examine how top quark Yukawa coupling can affect the mea-surement of trilinear Higgs coupling, we take the following effective Lagrangian where the term Y t = √ 2m t /v is the Yukawa couplings of top quark in the SM, and both a and b are dimensionless parameters. The parameter b is related with the CP violation.
In the standard model, a = 1 and b = 0. In the two Higgs doublet model (2HDM) with no CP violation, a = ctgβ and b = 0. If there is CP violation in the 2HDM, the CP even and CP odd neutral scalars could mix which leads to a non-vanishing b. Early efforts to probe this coupling at hadron colliders could be found in [44]. The study of measurement of these couplings at linear colliders can be found [45]. A recent study at the LHC and future hadronic colliders how to measure these two free parameters can be found in [46,47], where a different but equivalent parametrisation was used. Theoretical calculation of loop corrections from tth can be found at [48]. A recent study on the CP properties of the 2HDM could be found in [49]. A systematic analysis on the constraints from the Higgs precision measurement to either a or b, interested readers can refer [50][51][52] for such a detailed study.
The term λ SM = m 2 h /2v 2 ≈ 0.13, while λ 3 is a free dimensionless parameter. In various new physics models, this coupling can vary in a large range. For example, in the framework of an effective operator, the strong first order electroweak phase transition has been explored in [27], where λ 3 can be in the range [5/3, 3]. In the model with a singlet + SM, the trilinear coupling could be larger than the value of the SM by more than 20% to 200% [53], and this deviation is dependent upon the mass of the singlet.
There are mainly two-fold aims for this work: 1) we explore the sensitivity of the pure leptonic mode pp → hh → 4 + / E at a 100 TeV collider, 2) we examine the complementarity of the direct measurement oftth and the direct measurement of λ 3 in the future colliders.
The work is organised as follows. In section II, we study the cross section of the process gg → hh. In section III, we analyse the sensitivity of two types of the same sign leptonic final states gg → hh → 4 + / E in a 100 TeV collider. In section IV, we examine the issue how tth measurement can affect the determination of the trilinear Higgs couplings. In section V, we examine the complimentarily to determine related Higgs couplings at hadronic colliders and electron-positron colliders. We end this work with a few discussions on the detector issues to probe pure leptonic modes in a 100 TeV collider. We provide an appendix to describe the quasi-Monte Carlo method implemented in our code "wat" which has been used in the work to evaluate cross section and to generate unweighted signal events.

LIDERS
We implemented the effective Lagrangian described in Eq. (1) as a new model file by modifying the one-loop SM model file in MadGraph5/aMC@NLO [54]. The parameters a, b, and λ 3 and the corresponding tree level vertices are added by following the UFO protocol [55]. According to the OPP method [56], we add two new R2 terms at one-loop level which are related with the process gg → hh by following the information given in the two-Higgs doublet model [57], shown as below: where a 1 , a 2 are the color indice, and µ 1 , µ 2 are the Lorentz indice, and g s is the QCD coupling constant.
Then we interface the loop matrix element produced by MadGraph5/aMC@NLO [54] to our integration and event generation code to obtain the leading order cross section and unweighted events of signal.
At the leading order, cross sections in hadronic colliders can be parametrised as the function of theoretical free parametersa, b, λ 3 in the following form It is noticed that this cross section is sensitive to the signs of a and λ 3 , respectively, but is insensitive to the sign of b. There are two types of diagrams contributing to the process gg → hh: the box diagrams and the triangle diagrams. The terms independent of λ 3 come from the squared amplitudes of box diagrams. The terms proportional to λ 2 3 are from the squared amplitudes of triangle diagrams, and the terms proportional to λ 3 is from the interference between the box and triangle diagrams.
We use the numerical approach to fit the coefficients G 1 − G 7 for the LHC at 14 TeV, 33 TeV and a 100 TeV collider. We generate more than 100 points in (a,b,λ 3 ) space for each collision energies. The coefficients of cross sections at the 14 TeV, 33 TeV and 100 TeV are presented in Table IX. Due to the enhancement of gluon flux at a 100 TeV, it is noticed that all coefficients are 40 times larger than G 14 4 and G 14 5 , this enhancement factor is smaller than that of box diagrams due to the s-channel suppression for energetic gluon fluxes; and the interference coefficients G 100 6 and G 100 7 are also 40 times larger than G 14 6 and G 14 7 .
Another interesting observation is that the coefficient G 3 is 7 times larger than the coefficient G 1 . Typically, when b is much smaller than 1, the contribution of G 3 term can not be large. But if b can be of order one, then the contribution of G 3 can be sizeable. In the works [52,58], more operators have been taken into account. We have compared our results presented there and found agreement.

III. FULL LEPTONIC MODES OF SIGNAL AND BACKGROUND EVENTS
There are quite a few advantages of the full leptonic mode of gg → hh → 4 + / E for a 100 TeV collision. On the first hand, it is relatively efficient to be searched by experiments.
Targeted objects in the final state are leptons and missing energy. They can be reconstructed efficiently by subdetector systems of LHC detectors. The lepton reconstruction efficiency with P t > 5 GeV is more than 90% and particle identification can be made at detector level. On the second hand, the relatively clean signal and the signal is robust against the contamination of pileups and underlying events, since primary collision vertices of the signal events can be reconstructed.
As explained in Sec. II, we interface loop matrix element from Madgraph5/aMC@NLO [54] with our code based on QMC to perform phase-space integration and event generation. The generated events are further showered by PYTHIA6 [59] and then used to perform physical analysis. We have not taken into account the detector simulation in this work.
We use Madgraph5 to generate background events in a collision energy For both signal and backgrounds, higher order corrections are taken into account by normalizing the total cross section to their (N)NLO results, which is indicated by K factor(K = Table III. We adopt the NNLO result for signal in Ref. [60], and K factor for processes hZ, ttZ, ZW + W − ,tth, hW + W − are obtained by MadGraph5/aMC@NLO [54] under on-shell approximation. The K factor for tttt and ttW + W − at 100 TeV collider is unknown, and we adopt a value K = 1.3, which is the K factor for tttt at 14 TeV LHC. The K factor for the process pp → W + W − W + W − is also unknown, and we use the K factor for ZW + W − since both of them have the same initial states at hadron colliders.
The backgrounds processes can be roughly categorized into the following three types: • 1) The single Z processes include pp → Zh, pp → ZW + W − , and pp → Ztt. The last process has a cross section around 35 times larger than the former two processes due to its QCD nature, while Zh and ZW + W − have similar cross sections.
• 2) The top pair processes include pp → tth, pp → tttt, and pp → ttW + W − . We notice that htt is the dominant background in this category.
• 3) The pure electroweak processes include pp → hW + W − and pp → W + W − W + W − . Each of their cross section is smaller than that of the signal, but the sum of them is comparable to that of the signal after taking into account the K factors.
The cross sections of these processes are listed in the Table III. It is worthy of remarking that single-Higgs associated processes in the categories above are the main background events, and it turns out that tth is the dominant background for the full leptonic mode, which can greatly affect the significance.
In order to select the most relevant events, we introduce the following preselection cuts:  • We demand the future detector can have a better coverage of η( ) as |η( )| < 4 for identified leptons. Due to the good space resolution power and fine granularity in ECAL detector, we require that the minimal angular separation between two leptons is ∆R min ( , ) ≥ 0.15. In Fig. 2, we demonstrate the distribution of the maximal η max ( ) and the minimal angular separation of any a pair of two leptons ∆R min ( ).
When these two cuts on leptons are applied, more than 90% signal events can be accepted. We will discuss how the coverage of η can affect our results in the discussion section.
• Consider the fact that tt processes have quite large contributions to the background events, we introduce the "b-jet veto" to reject events with a tagged b-jet with P t (b) > 40 GeV and |η(b)| < 5. We assume that the b tagging efficiency is 60% and simply time a factor 0.16 to this type of background processes.
• For the decay mode M3, in oder to suppress the background events from meson decays in the final states and Z boson decay, we demand the invariant mass of both lepton pairs should fall into two windows either 15 GeV ≤ m + − ≤ 80 GeV or m + − ≥ 100 GeV.
From the results given in Table III, it is noticed that the background from pp → ZZ is huge which is about 2000 times larger in magnitude. Furthermore, due to the off shell Z and Preselection Cuts low energy hadron veto and Z mass veto γ * , the cut of invariant mass of two pairs of four leptons can only suppress the background down by 100 order at most. Our previous experiences in studying the 3 + 2j [39] reveals that background events from Z, γ * -exchange processes are difficult to suppress. Therefore, the first two modes are challenging in the SM and will be neglected in this analysis. When there is a significant enhancement to the Higgs pair production by new physics, these two modes should also be considered.
In this work, we will focus on the third and forth modes (labelled as "M3" and "M4", respectively) and we will perform a detailed analysis on these two modes. It is noticed that compared with the third mode although the forth mode has a smaller production rate, it enjoys less background contributions from the SM.

M3 AND M4 CASES
In this section, we explore the kinematic features of the signal and background events. It is impossible to fully reconstruct all the final particles due to 4 neutrinos in the final state.
The physical observables can be divided into two types. The first type is defined as the global and topological event shape observables for each event, and the second type is defined from the partial reconstruction method introduced later. The first type includes the following observables listed below.
• o1) The missing transverse momentum spectrum. Since there are 4 neutrinos in the final states, we expect that there should be a large transverse momentum. For the signal, as shown in Fig. 3(a), the distribution peaks near 60-80 GeV.
• o2)) The invariant mass of four leptons. This quantity is expected to capture the mass of mother particles. For the signal, as shown in Fig. 3(b), the distribution peaks near 100-150 GeV. Since four neutrinos can take away half of the energy of Higgs pair, so this quantity is expected to be close the mass of one Higgs boson.
• o3) The transverse mass of each event is constructed from the sum of 4-momentum of leptons (denoted as P 4 , which has components (E 4 , P x 4 , P y 4 , P z 4 ) ) and the missing transverse momentum ( / P T which has components ( / E T , / P x , / P y , 0)) where To construct this observable, we boost the 4-momentum of leptons such that the P z 4 = 0, i.e. (Ẽ 4 , P x 4 , P y 4 , 0). Then we construct the observable as • o5) The number of jets in each event with P t (j) > 40 GeV. As demonstrated in Fig. 3(e), we noticed that when demand n j ≤ 2, around 90% signal events can be selected out, though with a considerable background events from tt processes.
Below we explain how to construct the second type of observables.
Obviously, the most important information about the the signal events is the mass of Higgs boson. So it is crucial to extract this useful observable. Due to the fact that four neutrinos can not be fully reconstructed, instead we can only determine the visible masses of each Higgs boson from the identified leptons.
To determine the visible mass of each Higgs boson mass, we encounter a minor combinatorics issue: there are two possible combinations in each event. To determine which one is correct, we follow the minimal mass method introduced in [39] to determine the right combination, which can yield a correctness up to 94% here as we have checked this by using the parton level data sample. The method evaluates the sum of two visible masses of Higgs bosons for each combination, and picks out the smaller one as the correct combination. In contrast, we examine a method by using the angular separation of leptons, which can only yield a correctness up to 85% at most.
After having found the visible mass of Higgs boson in the signal, by using the standard M T 2 method we can split the missing transverse momentum into two parts and exploit the kinematic feature of pair production to reconstruct the transverse mass of Higgs boson.
Then we can construct the second type of observables as listed below.
• o6) The observable ∆m, which is defined as the mass difference between two mass sums of reconstructed Higgs bosons in two possible combinations. The distribution is shown in Fig. 3(f), which shows a large shape difference between the signal and background processes in shape.
• o7)) The first reconstructed partial mass of Higgs boson with two leptons of the same flavour in M3 case, which is labelled as h 1 . In the M4 case, the one with the hardest lepton is labelled as h 1 . The distribution of this quantity is shown in Fig. 3(g) and It is noticed that all of these reconstructed observables are crucial and important for both M3 and M4 cases. There exist strong correlations among these observables for signal events and they are related to the mass of Higgs boson while for the background they are not necessarily related to the mass of Higgs boson. Such a fact can be utilised to separate signal and background, which is the spirit of multivariate analysis methods.
A. The M3 case In the M3 case, the dominant background processes is the single Z associated processes, as clearly demonstrated in the Table V in the cut-based analysis. From Fig. 3(g), it is noticed that the mass cut can clearly affect the single Z processes.
In the cut-based analysis, we choose the sequential 4 cuts: 1) a cut on the number of jets n j ≤ 2, which is supposed to suppress background processes associated with a top pair; 2) a cut on the variable m T 2 , as we demand m T 2 < 110 GeV, which can greatly suppress the background processes like three body and four body productions; 3) a cut on the reconstructed visible mass of Higgs boson is imposed as m h 1,2 ≤ 60 GeV, which can suppress background events from the Zh and htt processes; 4) a cut on the mass difference ∆m > 50 GeV. By using these cuts, we can achieve the S/B and significance 0.25 and 4.2, respectively. According to the results based on the cut-based method, it is noticed that the processes associated with top pair are the main background for the "M3" mode. In order to suppress this type of background, we introduce a few top veto observables in our BDT analysis. For example, the visible invariant mass of the whole event, as demonstrated in Fig. 3(d), can help to separate background events with and without top quark pair. Furthermore, the number of jets in each events can help to suppress the background events with top quark pair, as demonstrated in Fig. 3(e). It is noticed that these two observables are correlated with each other for the background processes.
In Table V, both a cut-based analysis and a BDT analysis are presented to compare. As demonstrated in Fig. 4(b), after taking into account these observables to reject top quarks final states, we can further suppress background and gain in the S/B and significance up to 0.38 and 6.1, respectively, as represented in the  Fig. 5(h), where the shapes of background and signal look similar.
Therefore, we impose a cut on m(h vis 1 ), instead of both as in the M3 mode case. A cut-based analysis and a BDT analysis are presented in Table VI. For the third cut, in this case, as we emphases that we only impose a cut on the visible mass m(h vis 1 ) < 60 GeV. As demonstrated in Fig. 6(b), after taking into account these observables to reject top quarks final states, we can further suppress the main background and gain in the S/B and significance up to 1.6 and 6.8, respectively, as represented in the Table VI. After using the same methods to veto top pair associated background in the MVA method, a better result is yield, the S/B and significance can reach 1.9 and 9.2, respectively. If an integrated luminosity 20-30 fb −1 is assumed, then a significance 30.0 is expected.
When comparing the results given in Table V and Table VI, we notice that the S/B of the M4 case is much better than that of the M3 case, due to the lack of huge background processes as in the M3 case. Furthermore, the analysis for the M4 case is relatively simpler than the M3 case due to the background from single Z processes can be efficiently suppressed.

V. INTERPLAY BETWEEN gg → hh AND pp → tth
It is noticed that events from the tth final state are the dominant background for the M4 case, which is also true for the pp → hh → 3 + 2j + / E T mode explored in [39] at the HCs. Therefore, the measurement on the top Yukawa coupling can significantly affect the detection of Higgs pair production. Below we examine how the measurement of tth couplings can affect the determination of λ 3 . The correlation between the measurement of tth and the measurement of Higgs pair production can be investigated by the cross section of gg → hh given in Eq. (4). This issue has not been addressed in literatures. In Fig. 7, we demonstrate the correlation between the determination of tth at the LHC and a future 100 TeV collider, where the bounds are estimated from the measurement of 3 2j + / E. These bounds have not optimised for each value of a and λ 3 , we simply use the information of cross sections in this projection. As demonstrated in [39], when the discrimination of signal and background is optimised, we can expect better bounds.
For the LHC with a luminosity 3 fb −1 , we assume that tth couplings can be determined up to 20% (10% is estimated by using the boosted techniques for tth and h → bb [62,63]). We deliberately take a larger value for this since there only statistics are taken into account, which is denoted by two solid lines in Fig. 7(a) as upper and lower bounds from tth measurements.
For a 100 TeV collider with a luminosity 3 fb −1 , we assume that a 5% precision can be achieved, which is denoted by two solid lines in Fig. 7(b) as upper and lower bounds from tth measurements. According to the study of [64], where by using the production ratio tth/ttZ measurement and boosted Higgs and top quark, it is argued that this coupling can be measured up to a precision 1% or so when only statistics is considered. In reality, more background processes and detector effects must be considered for each tt decay final states, so we assume a precision up to 5% as a relatively conservative and loose estimation.
Comparing Fig. 7(a) and Fig. 7(b), it is noticed that a 100 TeV collider can greatly shrink the uncertainty in determining the λ 3 and a parameters. Due to a 4 dependence of the cross section σ(pp → hh), a 5% uncertainty of δa can induce an uncertainty of λ 3 up to 20% or so. If the coupling a can be determined to a precision 1%, that will be undoubtedly crucial to pinpoint the value of λ 3 down to 5%. It is worthy of mentioning that the two-fold ambiguity in a with the same cross section can be removed by using the method to check the differential distribution of leptons in the final state, as demonstrated in [39], which is robust against the contamination of underlying events and pileup effects. In the effective Lagrangian given in Eq. (1), the parameter b is related with the strength of CP violation. We explore the determination of b parameter from Higgs pair production, which is given in Fig. 8(a). From Fig. 8(b), it is easy to read out that the potential of a 100 TeV collider to probe the parameter space expanded by b and λ 3 is obviously better than that of the LHC.
Higgs pair production at a 100 TeV collider can bound b down to 0.4. In contrast,the LHC runs can only constrain this parameter to 1.2 or so at most. When we quote these number, we have not taken into account the uncertainty in a. Obviously, the direct measurement from tth could impose a better constraint to the value of b, either at the LHC or at a 100 TeV collider.
The correlation of a and b in Higgs pair is provided in Fig. 9(a). We plot the dependence of cross section of tth upon a and b in the same plots. From Fig. 9(b), it is noticed that the LHC can be sensitive to the region where a < 0 and b is sizeable, which corresponds a large production rate of gg → hh due to a constructive interference.
To consider the bounds from the tth measurement in Fig. 9(a) and Fig. 9(b), we have parametrised the cross section of tth at the LHC 14 TeV and at a 100 TeV collider as where the values of t 1 and t 2 for 14 TeV and 100 TeV are computed by using fit and are provided in Table ( For a 100 TeV collider, it is noticed that there is a 3-fold ambiguity when we combine the measurements of tth and hh. To remove this 3-fold ambiguity, the differential distribution of tth final state should be carefully analysed, as demonstrated in [46] where quite a few differential observables are proposed. It is interesting to observe that both single Higgs production and Higgs pair production can indirectly help to determine the interaction between top quark and Higgs boson at a 100 TeV collider, and it is pointed out that the Higgs pair production can also help to disentangle and resolve the nature of ultraviolet contributions to Higgs couplings to two gluons [65]. It is also interesting to study the correlations of the cross section σ(pp → tth) and σ(pp → hh) when the uncertainty from parton distribution function will be considered. We use the dataset of CT10 [66] to examine the uncertainty of PDF. A recent analysis on CT14 NNLO PDF can be found in [67].
We consider the correlations between the processes gg → hh and pp → tth in three cases. The first case is (a, b, λ 3 ) = (1, 0, 1), the second case is (a, b, λ 3 ) = (1, 0, −1), and the third case (a, b, λ 3 ) = (0, 1, 1). For these three cases, the central values of cross sections of gg → hh and pp → tth are so different that the uncertainty from PDF can not lead to an overlap among them. In order to examine the correlation caused by the uncertainty of PDF, we normalise the cross sections with the central values in each case and plot the 90% confidence level contours in Fig. 10.
In the first case, it is noticed that the uncertainty of PDF at the LHC 14 TeV can at most lead to an uncertainty 5% in the cross section of pp → ttH and 4% or so in that of gg → hh.
The uncertainty of gg → hh is almost the same at a 100 TeV, while that of pp → tth is shrunk to 3%.
At the LHC 14 TeV, when comparing the second case with the first case, we notice that the uncertainty of pp → ttH is enhanced due to the cross section becomes smaller. The correlation angle of the second case is different from that of the first case due to form factors proportional to b 4 and b 2 are different from those proportional to a 4 , a 3 , and a 2 . When the collision energy is 100 TeV, the correlation starts to become stronger.
In the third case, the uncertainty of pp → ttH is smaller due to normalisation when compared with the first case. Meanwhile, the correlation becomes weaker due to the interference and form factors changing, the area of the third case is fatter than that of the first case, in both the 14 TeV and 100 TeV collisions.

VI. COMPLEMENTARITY BETWEEN HADRON AND ELECTRON-POSITRON COLLIDERS
Below we briefly comments on the sensitivity to probe the Higgs trilinear coupling λ 3 from other production processes at hadronic colliders simply based on the cross sections. A recent study on the interferences between the contribution of trilinear coupling and those independent of trilinear coupling can be found in [68], where it was found that the interference is almost insensitive to the change of collision energy of hadronic colliders.
We first examine the vector boson associated production processes. To examine how the couplings between Higgs boson and the weak bosons can affect the determination of λ 3 , we parameterise the couplings of the interactions V V h and V V hh as d 1 and d 2 , respectively.
In the standard model, both d 1 and d 2 are equal to one. According to the current analysis on the single Higgs boson production and decay final states, it is known that the d 1 is close to 1. But d 2 has not yet measured.
With this parameterisation given in Eq. (6) and Eq. (1), the cross section of W/Zhh at hadron colliders can be parametrised as where form factors V i are evaluated numericaly and are given in Table VIII. We notice that  TeV and 100 TeV collider are tabulated. R 33 (R 100 ) is defined as signs of V 2 and V 5 are different from those of others and there exist interference cancellations between V 1 (V 4 ) and V 2 (V 5 ). This process should be sensitive to λ 3 due to the large fraction of contribution from terms proportional to λ 3 , but the production rate is 3-order smaller than the process gg → hh. This process has recently been studied in [69]. It is sensitive to a specific window of λ 3 and is complimentary to other channels.
Here we examine the cross section of the process tthh, which can be parametrised similar to Eq. (4), but with a different set of free parameters T 1 − T 7 .
The values of these free parameters are evaluated numerically and are tabulated in Table   IX. It is noticed that all these parameters are positive and the interference can occur when a and λ 3 opposite signs. The interference terms proportional to T 4 and T 5 have largest enhancement factors than that of other terms. The feasibility for hadron colliders of this process has been studied in [70,71]. The cross section of vector boson fusion processes at hadron colliders can be parametrised as follows where f 1 − f 6 are given in Table X.
One interesting observation is that the interference among f 1 − f 3 and f 4 − f 5 is severe.
Although each term is increasing with the increase of collision energy and is large, the total cross section is smaller than that of the process pp → tthh at a 100 TeV collider due to the strong cancellation induced by the interference. There is one tradeoff from this interference: the term proportional to λ 2 , suppress the reducible backgrounds. It is also noticed that the interference between W W jj and ZZjj is suppressed by colour and in-flow quark flux and can be safely neglected. Our results are consistent with the results presented in [68].
As pointed in [72], this process can be used to probe the quartic coupling g V V hh vector boson and Higgs boson and the LHC can constrain this quartic coupling to a precision 50%.
A 100 TeV collider can help to fix the quartic coupling g V V hh .
It is also interesting to mention the process of single Higgs production associated with a single top. As demonstrated in [51,[73][74][75][76][77][78], single top and Higgs processes can be useful to probe the relative sign between a and d 1 . Here we can parameterise the cross section as As demonstrated in [51,75], LHC can determine the sign of sign(a d 1 ). The total cross section is mainly determined by the term proportional to a d 1 . As shown in Table XI, in the SM, there exists a strong cancellation between the term a d 1 and the rest. When this term switches its sign, the the cross section will be enhanced by one order or so. Therefore, this process is sensitive the sign of a d 1 . The single top and Higgs boson production could also be a signal of the flavour changed decay of top quark t → h + u and t → h + c from the processes pp → tt via the flavour changing neutral current processes, as shown in [79].
It is also interesting to examine the process of Higgs pair production associated with a single top, the cross section in the SM has been explored in [78]. Here we can parameterise the cross section as where we only keep those sizeable terms in this formula and drop those tiny terms. The values of J i , at the LHC 14 TeV, a 33 TeV collider, and a 100 TeV collider, are tabulated in Table XII. It is noticed that there exists a strong cancellation among terms, like (J 1 − J 3 , . If there is a new physics which can break this strong cancellation, it is expected that this process can a sensitive probe. Below we simply examine the complementarity of e + e − machines to probe trilinear coupling λ 3 . It is noticed that the vertex tth can be probed by measuring the production of e + e − → tth and its cross section can be parametrised as the following formula σ(e + e − → tth) = e 1 (s)a 2 + e 2 (s)b 2 + e 3 (s)ad 1 + e 4 (s)d 2 1 , where e i (s) , i = 1, 2, 3, 4 are form factors, and their dependence upon the collision energy s can explore numerically. In Fig. 11, we provide the numerical results. The analytic form of these form factors at tree level can be found in [80,81].
To extract the tth couplings at the ILC, the reference [82] had performed a detailed MC study by considering the + 6j and 8 j final states. The main background ttW/Z/(gg) and tW b are considered. It is observed that b jet veto is crucial to suppress the tW b background.
The sensitivity to a is around 10% at the ILC with a integrated luminosity 500 fb −1 , which is similar to the sensitivity of the LHC high luminosity runs.
It should be point out that at the hadron colliders, the cross section given in Eq. (5) should have the similar structure as given in Eq. (12). But due to the large gluon flux, the interference terms proportional to a and the terms independent of a and b are relatively tiny and can be safely neglected.
In contrast, the cross section of the vector boson fusion processes e + e − → ν eνe hh (mainly WW fusion processes) increases with the increase of collision energy √ s and can be parametrised as where form factors w 1 (s) − w 6 (s) are demonstrated in Fig. 11(b). It is noticed that w 2 (s) has a different sign from w 1 (s) and w 3 (s), and w 4 (s) has a different from w 5 (s). Similar to the VBF processes at hadron colliders, the large cancellations among w 1 (s) − w 3 (s) and w 4 (s) − w 5 (s) lead to a fact that the total cross section is mainly dependent upon the term proportional to λ 2 3 , which makes the process e + e − → ν eνe hh sensitive to the measurement of trilinear coupling λ 3 . A Monte Carlo study on the sensitivity of λ 3 via two jjbbbb and ννbbbb final states can be found in [83]. Both the processes e + e − → ν eνe hh and the processes e + e − → Zhh have no dependence upon the coupling of tth. The cross section of e + e − → Zhh can be expressed as The form factors are shown in Fig. 11(c), where it is noticed that z 2 and z 4 are positive at the very beginning but become negative (the sign is switched in the plot) when √ s > 600 GeV.
In the large √ s region, a strong cancellation among z 1 , z 2 , and z 3 occurs. A similar exact cancellation happens between z 4 and z 5 when √ s is large. Therefore the main contribution at high energy region is determined by z 1 − z 3 terms, which is much larger than the terms proportional to λ 3 .
Due to the huge irreducible background from the process of sequential emission of two Higgs bosons and detector effects, even near 500 GeV the uncertainty of λ 3 is larger than 30% [84,85]. With higher collision energy (say √ s ≥ 700 GeV), the production rate proportional to λ 3 decreases rapidly and the uncertainty to extract λ 3 increases. Therefore, high energy collisions (say higher than 800 GeV) for the process e + e − → Zhh could not offer us a better way to extract trilinear Higgs coupling λ 3 .
The structure of the cross section of the processes e + e − → tthh are more complicated, because all terms in the Lagrangian given by Eq. (1) and Eq. (6) come into play a role.
The leading dominant terms where the form factors h 1 − h 6 are displayed in Fig. 11(d). Although there are other terms, like terms proportional to b 2 λ 2 3 , d 1 a 2 λ 3 , and so on, but numerically they are small when compared with these 6 form factors so we can simply neglect them. Both this process and e + e − → bbhh probe trilinear Higgs coupling, as demonstrated in [86,87].
It is interesting to notice that the cancellations observed in processes pp → W/Z hh, pp → W W/ZZjj → hhjj, e + e − → Zhh, and e + e − → ν eνe hh, are not accidental and are tightly related with the unitarity in a theory with spontaneous symmetry breaking [88,89].
In the SM and SUSY, the process e + e − → hh via loop processes [90,91], the maximum production rate is around 0.01 fb or so when the collision energy is around 350-450 GeV. In the collision with polarised beans, the production rate can be enhanced by a factor 6-7. For the final state with 4b, the main background should be e + e − → ZZ → 4b, which is around 36.8 fb when collision energy is set to be 400 GeV. We expect a certain degree of feasibility if the invariant mass of two b jet can be reconstructed. This process can not directly help to probe trilinear Higgs coupling λ 3 . Nevertheless, it might be useful probe to the coupling of V V hh.
In the 2HDM with CP violation, the vertices Z(h i ∂h j − h j ∂h i ) can be induced at tree level [81,92,93]. If two lightest Higgs bosons of the model are kinematically reachable at a future e + e − collider, this process should be carefully studied. If the e + e − colliders could run in the photon-photon collision mode, then the process γγ → hh [94][95][96] can be useful to directly probe trilinear coupling and the new physics inside the loops.

VII. DISCUSSIONS AND CONCLUSIONS
In this work, we examine the discovery potential of the full leptonic final state of the Higgs pair production at a 100 TeV collider. Based on our provious analysis gg → hh → 3l2j / E , we also explore the correlation between the experimental bounds on the Yukawa couplings a and b and the bounds on λ 3 at both the LHC and a 100 TeV collider.
Due to the small production of the full leptonic final state in the SM, it is challenging to discover the Higgs pair production via the two modes explored in this work. But if there is new physics which might enhance the Higgs pair production by one or two order in magnitude, like some bench mark points allowed in the 2HDM [97], this mode could be considered.
It is obviously that if we can combining all final states of Higgs pair production, like bbγγ, W W * γγ, bbτ τ , etc, we can achieve a better bound on λ 3 .
Here we would like to briefly address the issues of detector to this type of signal. We show the acceptance efficiencies to the signal events in The second issue is related with the missing energy resolution, which can be related with HCAL coverage. As we demonstrated, the missing energy reconstruction is quite crucial in order to obtain kinematical observables which can suppress the background events efficiently.
Nevertheless, in the 100 TeV collider, if the HCAL detector coverage is less than 6, the missing energy from forward jets can be larger than 100 GeV, as demonstrated in [1]. So if HCAL detector coverage can reach to η(j) < 8, which can be good not only for the vector boson fusion processes but also for the Higgs pair production processes.

ACKNOWLEDGMENTS
We would like to thank Fapeng Huang  In this appendix, we explain the QMC and a new reweighting algorithm based on QMC.
Considering a d-dimensional integral This method is called quasi-Monte Carlo method, and the point set is called quasi-Monte Carlo rule [98].
Comparing to Monte Carlo integration, which using independently distributed random points and can only achieve a convergence speed proportional to O(N −0.5 ) , QMC integration can achieve a convergence speed near O(N −1 ), if a proper QMC rule is adopted. Two families of QMC rules attract most of interest: one is consist of digital nets and digital sequences; the other is lattice rules. In this work, we adopt the rank-1 lattice rule (R1LR): x i = i z n , i = 0, . . . , n − 1, where z, known as the generating vector, is required that all its components should be integer and relatively prime to n. The braces around the vector in the Eq. (A3) means only the fractional part of each component is taken.
3. Generate a point x with probability density function g( x) dxg( x) , and a random number r with uniform distribution in [0, 1). If this point is accepted, otherwise this point is rejected. This step can be repeated many times to obtain desired number of points.
The efficiency of this algorithm, which is defined by the ratio of accepted number of points and total number of points in the last step, is determined by the sampling function g( x), and this function can also be used in MC integration via importance sampling. The most widely used algorithm for the sampling function in high energy physics is the VEGAS algorithm proposed by Lepage [102], which has a factorized structure as the following where g i (x i ) are step functions. This algorithm starts with g( x) = 1, then samples with several points based on g( x) and optimizes g( x) based on sampled results. The sampling and optimizing procedure can be iterated for several times, to further optimize g( x). Although this function allow us to deal with function with large peak, a general function is unfactorizable and the efficiency of this algorithm is limited by this. Another disadvantage of this algorithm is that the optimization is based on information gained via MC sampling, which suffers from fluctuation.
As we found that QMC could give a better estimation of integral than MC methods, based on the structure of R1LR, we proposed a new g( x). The sampling function is constructed based on QMC points, so the fluctuation is much smaller. In addition, it do not assume a factorized structure, and the efficiency of this algorithm could be improved to arbitarily near 1 if a larger point set is adopted.
The sampling function in this algorithm is described in the following: 1. Any vector v ij = x i − x j is the period of the lattice, and we can choose d linearly independent vectors from them. Since there are multiple choices, we choose the shortest ones and mark them as e i (i = 1, 2, . . . , d).
2. With these linearly independent vector and one origin x k , we can build an affine coordinate system, and mark it as S k .
4. then (y 1 , y 2 , . . . , y n ) is the random point satisfying requested distribution in the affine coordinate system S k , i.e. the corresponding coordinates in original coordinates system is where the pair of braces means that only the fractional part is taken, as before in Eq.
(A3). In Table XIV, we demonstrate the efficiency of this algorithm in comparison with VEGAS algorithm. The function used in the test is the differential cross section of gluon fusion into Higgs pair at 14TeV LHC, where the dimension of integration is 4. We adopt the R1LR algorithm and perform an integration with a lattice containing 1024 points and repeat the computation with 4 shifts. After that we test the reweighting efficiency. For VEGAS algorithm, we compute 1024 points in each iteration and total number of iteration is 4, the same number of points as in the R1LR algorithm. It is observed that the QMC-based algorithm is better than VEGAS, both in the speed of convergence and reweighting efficiency.