An Exploratory study of Higgs-boson pair production

Higgs-boson pair production is well known being capable to probe the trilinear self-coupling of the Higgs boson, which is one of the important ingredients of the Higgs sector itself. Pair production then depends on the top-quark Yukawa coupling $g_t^{S,P}$, Higgs trilinear coupling $\lambda_{3H}$, and a possible dim-5 contact-type $ttHH$ coupling $g_{tt}^{S,P}$, which may appear in some higher representations of the Higgs sector. We take into account the possibility that the top-Yukawa and the $ttHH$ couplings involved can be CP violating. We calculate the cross sections and the interference terms as coefficients of the square or the 4th power of each coupling $(g_t^{S,P}, \lambda_{3H}, g_{tt}^{S,P})$ at various stages of cuts, such that the desired cross section under various cuts can be obtained by simply inputing the couplings. We employ the $H H \to \gamma\gamma b \bar b$ decay mode of the Higgs-boson pair to investigate the possibility of disentangle the triangle diagram from the box digram so as to have a clean probe of the trilinear coupling at the LHC. We found that the angular separation between the $b$ and $\bar b$ and that between the two photons is useful. We obtain the sensitivity reach of each pair of couplings at the 14 TeV LHC and the future 100 TeV pp machine. Finally, we also comment on using the $b\bar b \tau^+ \tau^-$ decay mode in Appendix.


I. INTRODUCTION
A boson was discovered at the Large Hadron Collider (LHC) [1,2]. After almost all the Run I data were analyzed, the measured properties of the new particle are best described by the standard-model (SM) Higgs boson [3,4], which was proposed in 1960s [5]. The most constrained is the gauge-Higgs coupling C v ≡ g HW W = 0.94 +0. 11 −0.12 , which is very close to the SM value [6]. On the other hand, the relevant top-and bottom-Yukawa couplings are not determined as precisely as C v by the current data. Nevertheless, they are within 30 − 40% of the SM values [6].
Until now there is no information at all about the self-couplings of the Higgs boson, which emerges from the inner dynamics of the Higgs sector. For example, the trilinear couplings from the SM, two-Higgs doublet models (2HDM), and MSSM are very different from one another. Thus, investigations of the trilinear coupling will shed lights on the dynamics of the Higgs sector. One of the best probes is Higgs-boson-pair production at the LHC. There have been a large number of works of Higgs-pair production in the SM [7][8][9][10], in modelindependent fashion [11][12][13][14], and in special models beyond the SM [15] and in SUSY [16].
In the SM, Higgs-pair production receives contributions from two entangled sources, the triangle and box diagrams. The triangle diagram involves the Higgs self-trilinear coupling and the top-Yukawa coupling while the box diagram involves only the top-Yukawa coupling.
In order to probe the effects of the Higgs trilinear coupling, we have to disentangle the triangle diagram from the box diagram. We anticipate that the triangle diagram, which contains an s-channel Higgs propagator, does not increase as much as the box diagram as the center-of-mass energy √ŝ ≡ M HH increases. Therefore, the box diagram tends to give more energetic Higgs-boson pairs than the triangle diagram. Thus, the opening angle in the decay products of each Higgs boson can be used to isolate the triangle-diagram contribution.
Indeed, we found that the angular separation ∆R γγ and ∆R bb between the decay products of the Higgs-boson pair are very useful variables to disentangle the two sources.
Here we also entertain the possibility of a dimension-5 operator ttHH, which can arise from a number of extended Higgs models, including composite Higgs models or some general 2HDM's. For example, in a general 2HDM we can have a diagram with a (t L t R ϕ) vertex and a (ϕHH) vertex connected by the heavy ϕ. When the heavy ϕ is integrated out, we are left with the contact diagramt L t R HH. The anomalous ttHH coupling can contribute to Higgs-pair production via a triangle diagram. This triangle diagram is similar to the triangle diagram with the trilinear Higgs self coupling, except that it does not have the s-channel Higgs propagator. We shall show that the new contact diagram will give terms that can be combined with the terms of the triangle diagram, as in Eq. (4). We note that the kinematic behavior of the triangle diagram induced by the dim-5 ttHH contact interaction is different from that induced by the trilinear Higgs self coupling, because of the absence of the Higgs propagator in the contact diagram.
In this work, we adopt the effective Lagrangian approach, taking the liberty that the involved Higgs-boson couplings can be varied freely within reasonable ranges. The relevant couplings considered in this work are (i) the top-quark Yukawa coupling, (ii) the trilinear Higgs self-coupling, and (iii) the contact-type ttHH coupling. In the top-quark Yukawa and contact-type ttHH interactions, we take into account the possibility of the simultaneous presence of the scalar and pseudoscalar couplings which can signal CP violation. The rationale behind the CP-odd part is that the current data, other than the EDM constraints, cannot restrict the CP-odd part. The EDM constraints, however, depend on a number of assumptions and may therefore be weaken because of cancellation among various CPviolating sources [17]. On the other hand, ttHH is a dimension-5 operator, which may originate from a genuine dim-6 operator, e.g., (Q L Φt R )(Φ † Φ) after electroweak symmetry breaking, Φ = (0, (v + H)/ √ 2) T . This operator is thus suppressed only by two powers of higher scale, such that it can give a significant contribution at the LHC energies.
Our strategy is first to find a useful expression for Higgs-boson pair production cross sections in terms of these couplings, see Eq. (7). In this work, specifically, we employ the bbγγ decay mode of the Higgs-boson pair and look into the angular separation between the b andb and that between the photons. It is shown that one can map out the possible regions of Higgs couplings assuming certain values of measured cross sections, though it is channel dependent. Thanks to the largest branching 3 ratio of the Higgs boson into bb, the angular separation between the bottom-quark pair is an useful tool in for most of the proposed channel at the LHC.
In summary, the current work marks a number of improvements over previous published works as listed in the following: 1. We have included the CP-odd part in the top-Yukawa coupling. The CP-even and CP-odd parts are constrained by an elliptical-like equation by the current Higgs-boson data, as shown in Fig. 4 of Ref. [4]. Note that the effects of the CP-odd part of the top-Yukawa coupling at the LHC, ILC, and photon colliders were studied in Ref. [18], though we study the effects in depth for the LHC here.
Furthermore, we also include the CP-odd part of this contact coupling.
3. We have calculated an easy-to-use expression to obtain the cross sections as a function of the involved couplings at each center-of-mass energy. We also obtain similar expressions in various kinematic regions such that one can easily obtain the cross sections under the proposed experimental cuts for arbitrary values of the Higgs couplings.
4. With assumed uncertainties in the measurements of cross sections, we can map out the sensitivity regions of parameter space that can be probed at the LHC.
The organization of the work is as follows. In the next section, we describe the formalism for our exploratory approach and present an expression for the Higgs-boson pair production cross section in terms of various combinations of the Higgs couplings under consideration.
In Section III, we examine the behavior of each term of the cross section versus energies.
In Section IV, employing the HH → (γγ)(bb) decay mode, we illustrate how to extract the information on the Higgs couplings by exploiting the angular separations between the Higgs decay products. There we also discuss the prospect for the 100 TeV pp machine. We conclude in Section V and offer a few comments with regard to our findings. In Appendix, we compare the SM cross sections at 14 TeV with those at 100 TeV for the process pp → H → γγbb and give some comments on the τ + τ − bb decay mode.

II. FORMALISM
Higgs-boson pair production via gluon fusion goes through a triangle diagram with a Higgs-boson propagator and also through a box diagram with colored particles running in it. The relevant couplings involved are top-Yukawa and the Higgs trilinear self coupling. We further explore the possibility of a dim-5 anomalous ttHH contact coupling [14]. They are given in this Lagrangian: In the SM, λ 3H = g S t = 1 and g P t = 0 and g S,P tt = 0. The differential cross section for the process g(p 1 )g(p 2 ) → H(p 3 )H(p 4 ) in the SM was obtained in Ref. [7] as where functions F S = F , F SS = F , and G SS = G with F , and G given in Appendix A.1 of Ref. [7].
Here we extend the result to including the CP-odd top-Yukawa and the anomalous ttHH couplings: More explicitly in terms of each combination of couplings and ignoring the proportionality constant at the beginning of the equation, the above equation becomes where F P = F A , F SP = F , and G SP = G with F A , F and G given in Appendix A.2 of Ref. [7] while F P P = F , and G P P = G with F and G in Appendix A.3 of Ref. [7]. In the heavy quark limit, one may have [7] F S = + 2 leading to large cancellation between the triangle and box diagrams.
The production cross section normalized to the corresponding SM cross section, with or without cuts, can be parameterized as follows: +λ 3H e 1 (s)g S t g S tt + f 1 (s)g P t g P tt + g S tt e 2 (s)(g S t ) 2 + f 2 (s)(g P t ) 2 + e 3 (s)(g S tt ) 2 + f 3 (s)g S t g P t g P tt + f 4 (s)(g P tt ) 2 where the numerical coefficients c 1,2,3 (s), d 1,2,3,4 (s), e 1,2,3 (s), and f 1,2,3,4 (s) depend on s and 6 experimental selection cuts. Upon our normalization, the ratio should be equal to 1 when g S t = λ 3H = 1 and g P t = g S,P tt = 0 or c 1 (s) + c 2 (s) + c 3 (s) = 1. The coefficients c 1 (s) and c 3 (s) are for the SM contributions from the triangle and box diagrams, respectively, and the coefficient c 2 (s) for the interference between them.
Once we have the coefficients c i , d i , e i , and f i 's, the cross sections can be easily obtained for any combinations of couplings. Our first task is to obtain the dependence of the coefficients on the collider energy √ s, Higgs decay channels, experimental cuts, etc.

III. BEHAVIOR OF THE CROSS SECTIONS
We first examine the behavior of each piece of cross sections versus energies. We show the coefficients c i , d i , e i , f i 's at √ s = 8, 14, 33, 100 TeV in Table I  tremendous increases in f 4 . It is easy to see that the contact diagram is dim-5 and obviously grows with energy. At certain high enough energy, it may upset unitarity.
We can examine the validity of the anomalous ttHH contact coupling by projecting out the leading partial-wave coefficient for the scattering tt → HH. At high energy, the The leading partial-wave coefficient is given by 1 We note our results are in well accord with those in literature when the comparison is possible. The values of c 1,2,3 and e 1,2,3 at √ s = 100 TeV, for example, are in good agreement with those presented in Ref. [14].   TeV .
Therefore, the anomalous ttHH contact term can be safely applied at the LHC for g S tt 3−5 as most of the collisions occur at √ŝ a few TeV.
To some extent we have understood the behavior of the triangle, box, and contact diagrams with the center-of-mass energy, which is kinematically equal to the invariant mass M HH of the Higgs-boson pair. One can then uses M HH to enhance or reduce the relative contributions of triangle or box diagrams. The higher the M HH the relatively larger proportion comes from the box and contact diagrams. Since M HH correlates with the boost energy of each Higgs boson, a more energetic Higgs boson will decay into a pair of particles, which have a smaller angular separation between them than a less energetic Higgs boson. Therefore, the angular separation ∆R ij between the decay products i, j is another useful kinematic variable to separate the contributions among the triangle, box, and contact diagrams. 8

IV. NUMERICAL ANALYSIS
The Lagrangian in Eq. (1) consists of five parameters: the scalar and pseudoscalar parts of the top-Yukawa coupling g S,P t , the scalar and pseudoscalar parts of the anomalous contact coupling g S,P tt , and the Higgs trilinear coupling λ 3H . In order to facilitate the presentation and understanding of the physics, we study a few scenarios: 1. CPC1-the top-Yukawa coupling involves only the scalar part and the scale in the anomalous contact coupling is very large -only g S t and λ 3H are relevant. The relevant coefficients are c 1 , c 2 , and c 3 .
2. CPC2-the top-Yukawa and the anomalous contact couplings involve only the scalar part -g S t , g S tt , and λ 3H are relevant. The relevant coefficients are c 1 , c 2 , c 3 , e 1 , e 2 , and e 3 .
3. CPV1-the top-Yukawa coupling involves both the scalar and pseudoscalar parts - 4. CPV2-the contact ttHH coupling involves both the scalar and pseudoscalar partsg S tt , g P tt , and λ 3H are relevant while the top-Yukawa coupling is kept at fixed values. In this case, all the coefficients become relevant. In one of the simplest cases with g S t = 1 and g P t = 0, for example, the relevant coefficients are c 1 , c 2 , c 3 , e 1 , e 2 , e 3 , and f 4 .
Note that the above scenarios have been studied using different approaches and separately in literature: CPC1 in Ref. [12], CPC2 in Ref. [13], and CPV1 in Ref. [18]. The scenario CPV2 is new with the CP-odd component of the ttHH, in which the CP-odd component g P tt of the ttHH coupling appears either with g P t or in square. We used the CTEQ6L1 [19] with both the renormalization and factorization scales µ = M H for the parton distribution functions. Since we focus on the ratios relative to the SM predictions, we anticipate the uncertainties due to scale dependence, choice of parton distribution functions, and experimental acceptance are reduced to a minimal level 2 . For the branching ratios of the Higgs boson we employ the values for the SM Higgs boson listed in the LHC Higgs Cross Section Working Group [20]. Here we ignore the slight variation in the diphoton branching ratio due to the change in the top-Yukawa coupling. 3 This is because we do not want to interfere with the more important goal of the work -interference effects between the triangle and box diagrams and the sensitivity to the trilinear Higgs coupling.
For Higgs-pair production we used a modified MADGRAPH implementation [10,23], which allows us to vary the top-Yukawa, Higgs trilinear, and the contact ttHH couplings.
In this work, we make a working assumption that the bbγγ background can be estimated with a reasonable accuracy and can be extracted from the experimental data. Once the background is subtracted from the data, we are left with the signal events. In order to do so, we impose a set of basic cuts for the acceptance of the b quarks and photons, and use B-tagging according to the ATLAS template in Delphes v.3 [26], as well as an invariant mass window on the bb pair and the photon pair around the Higgs boson mass. The basic cuts are Note that in the appendix we shall discuss the feasibility of using the τ + τ − bb mode. With this set of basic cuts we continue to study the signal events in various kinematical regions separated by ∆R γγ and ∆R bb .
A. CPC1: g S t and λ 3H This is the simplest scenario to investigate the variation of the triangle and box diagrams with respect to changes in the Yukawa and trilinear couplings. The corresponding coefficients in Eq. (7) are c 1 , c 2 , c 3 . We let the Higgs boson pair decay into Note that if the process is studied without detector simulation, the distributions for other decay channels, like bbτ τ or γγτ τ would be the same. Nevertheless, the resolutions for b, τ , and γ are quite different in detectors, and so are the backgrounds considered for each decay channel. In the following, we focus on the γγbb channel, which has been considered in a number of works.
We use MADGRAPH v.5 [24] with parton showering by Pythia v.6 [25], detector simulations using Delphes v.3 [26], and the analysis tools by MadAnalysis 5 [27]. We have verified that the coefficients c 1,2,3 using the default, ATLAS, and CMS templates inside Delphes v.3 are within 10% of one another. From now on we employ the ATLAS template in detector simulations. We show the coefficients in Table II. We will come back to this table a bit later.
We Next, we look at the distributions versus ∆R γγ and ∆R bb in Fig. 3. The lines are the same as in Fig. 2. As we have explained in the previous section, ∆R γγ and ∆R bb between the decay products of each Higgs boson are useful variables to separate the triangle and box contributions. The angular distribution ∆R between the two decay products of each Higgs boson correlates with the energy of the Higgs boson, which in turns correlates with the invariant of the Higgs-boson pair. The higher the invariant mass, the more energetic the Higgs boson will be, and the smaller the angular separation between the decay products will be. Therefore, the triangle diagram has wider separation than the box diagram.
It is clear that the distributions of ∆R γγ and ∆R bb have similar behavior within uncertainties. The box diagram and also the SM, which is dominated by the box contribution, have a peak at ∆R γγ or ∆R bb less than 2.0, while the triangle diagram prefers to have the majority at larger ∆R γγ or ∆R bb , say between 2 and 3. We therefore come up with (i) ∆R γγ > 2 (< 2), (ii) ∆R bb > 2 (< 2), and (iii) ∆R γγ > 2 and ∆R bb > 2 (both < 2) to enrich the sample of triangle (box) contribution. In the following, we use ∆R to denote either ∆R γγ or ∆R bb , unless stated distinctively.
We can now look at Table II, where the coefficients for the ratio of the cross sections σ(gg → HH)/σ SM (gg → HH) as in Eq. (7) are shown. In the CPC1 case, the relevant coefficients are c 1 , c 2 , and c 3 in which c 1 is induced by the triangle diagram, c 3 by the box diagram, and c 2 by the interference between them. The rows labeled "Basic Cuts" are the ratio of cross sections under the set of cuts in Eq. (8). In the same Table, we also show the coefficients obtained after applying the angular-separation cuts of ∆R > 2 or < 2 and both ∆R γγ,bb > 2 or < 2. It is clear that ∆R γγ > 2 (< 2) enriches the triangle-diagram (box-diagram) contribution. Similar is true for ∆R bb > 2 (< 2). Further enhancement of triangle diagram can be obtained with both ∆R γγ > 2 and ∆R bb > 2, and vice versa for box diagram.
In the following, we investigate the sensitivity in the parameter space (g S t , λ 3H ) that one can reach at the 14 TeV with 3000 fb −1 luminosity by using the measurements of cross sections in various kinematical regions. Since we have found that the triangle and box contributions can be distinguished using the ∆R cuts, we make use of the measured cross sections in the kinematical regions separated by these cuts.
There are two issues that we have to considered when we take the measured cross sections in the kinematical regions. First, the SM backgrounds for the decay channel that we consider, and second, the Next-to-Leading-Order (NLO) corrections [28]. It was shown in Ref. [28] that the NLO and NNLO corrections can be as large as 100% with uncertainty of order 10-20%. The SM backgrounds, on the other hand, can be estimated with uncertainties less than the NNLO corrections. We therefore adopt an approach that the signal cross sections (after background subtraction) are measured with uncertainties of order 25-50%. About the signal cross sections, the first and the second columns of Table III show the SM cross sections for the process pp → HH → γγbb with detector simulations under various cuts at the 14 TeV LHC 4 . We have taken account of the SM NLO cross section σ SM (pp → HH) 34 fb, the Higgs branching fractions, and both the photon and b-quark reconstruction efficiencies with angular-separation cuts of ∆R bb ,γγ . In various kinematical 4 Before applying the basic cuts, we find the cross sections are 8.92 × 10 −2 and 2.41 in fb for  and SM-14(τ + τ − bb), respectively, which agree very well with those in Ref. [8]. While our cross sections after applying the basic cuts are smaller than those presented in Ref. [8] by a factor of ∼ 3 for γγbb and a factor of ∼ 30 for τ + τ − bb. This is basically because we have implemented full detector simulation to reconstruct b quarks, photons, and τ leptons from Higgs and partly due to different experimental cuts applied and different b-and τ -tagging efficiencies taken. One may need to optimize the cuts to increase the signal to background ratio but it is beyond the scope of this paper and we will pursue this issue later in our future publication. Incidentally the SM-100(γγbb) cross section is 3.73 fb before applying the basic cuts.
regions depending on the angular cuts, the cross sections range from ∼ 0.001 fb to ∼ 0.01 fb. With an integrated luminosity of 3000 fb −1 , we expect of order 30 signal events when the cross section is 0.01 fb. An estimate of the statistical error is given by the square root of the number of events √ N , which is then roughly 20% of the total number. Taking into account the uncertainty of order 10 − 20% from NLO and NNLO corrections, in this work, we use a total uncertainty of 25 − 50% in the signal cross section in the estimation of sensitivity of the couplings. Our approach is more or less valid except for the case in which both the ∆R bb > 2 and ∆R γγ > 2 cuts are imposed simultaneously. It would be challenging to measure this size of cross section only in the HH → γγbb mode and one may need to combine the measurements in different Higgs-decay channels. Or, one may rely on the future colliders such as a 100 TeV pp machine with larger cross sections and/or higher luminosities.
In Fig. 4, we show the contour lines of σ(gg → HH)/σ SM (gg → HH) = 0.5 and 1.5 with the Higgs boson pair decaying into (γγ)(bb). In each panel, we assume three measurements of the ratios corresponding to basic cuts (orange lines), ∆R > 2 (dashed black lines), and ∆R < 2 (solid black lines): here ∆R represents ∆R γγ (upper left), ∆R bb (upper right), or ∆R γγ,bb (lower). Therefore, for example, if the bacis-cuts cross section ratio is measured to be consistent with the SM prediction within 50% error, any points in the two bands bounded by the two pairs of orange lines are allowed. In each band, a rather wide range of g S t and λ 3H is allowed although they are correlated. Suppose we only make one measurement of the cross section without or with a cut on ∆R, we would not be able to pin down useful values for g S t and λ 3H . However, since the shape of the three bands are not exactly the same, we can make use of three simultaneous measurements in order to obtain more useful information for the couplings g S t and λ 3H . In the upper-left panel of Fig. 4, we suppose that one can make three measurements of cross sections: with basic cuts, ∆R γγ > 2, and ∆R γγ < 2. We assume that the measurements agree with the SM predictions within 25% or 50% uncertainty. The region of parameter space in (g S t , λ 3H ) bounded by the three measurements is shown by the lighter purple region for 50% uncertainty and darker purple region for 25% uncertainty. Similarly the upper-right panel is for the regions with the ∆R bb cut. In the lower panel, we show the regions with the combined cuts of ∆R γγ and ∆R bb : both larger than or smaller than 2. The implications from the measurements are very significant. First, all panels show that g S t is significantly away from zero if one can simultaneously measure the cross sections (no matter with 25% or 50% uncertainties) with basic cuts, ∆R γγ > 2, and ∆R γγ < 2; and similarly for ∆R bb and using both distributions. Second, as shown in the lower panel, the value for λ 3H is statistically distinct from zero if one measures the cross sections with a 25% uncertainty. This is achieved by using both ∆R γγ and ∆R bb > 2 or < 2. Furthermore, from the lower panel in Fig. 4 we can see that with 25% level uncertainty the values of λ 3H sensitivity regions are 0.3 |λ 3H | 2.6.
We can repeat the exercise with the measured cross sections being multiples of the SM predictions. We show the corresponding 25% and 50% regions in Fig. 5 for σ/σ SM = 0.5, 1, 2, 5, 10. Only with both ∆R γγ and ∆R bb , one can really tell if λ 3H is significantly distinct from zero. The sensitivity regions for λ 3H for various σ/σ SM are indicated by the darker color areas.
If the top-Yukawa coupling can be constrained more effectively by Higgs production, by ttH production, or by single top with Higgs production in the future measurements, say g S t = 1 ± 0.1 (10% uncertainty), it can help pinning down the acceptable range of λ 3H . However, even in this case, we emphasize the importance of simultaneous independent measurements, as illustrated in the following argument. In the limit of g S t = 1, the ratio of the cross sections is given by Suppose σ(gg → HH) is measured to be the same as σ SM (gg → HH) and then, using the relation c 1 (s)+c 2 (s)+c 3 (s) = 1, one may find the two solutions for λ 3H : 1 or −c 2 (s)/c 1 (s)−1.
For example, one may have λ 3H = 1 or 4 at most if only the basic-cuts ratio is measured, see Table II. Therefore, one cannot determine λ 3H uniquely with only one measurement even when the measurement is very precise and the exact value of g S t is known. It is unlikely to resolve this two-fold ambiguity at the LHC even we assume the three measurements of the ratios, as shown in Fig. 4. Also, the situation remains the same at the 100 TeV pp machine in which we have λ 3H = 1 or 5 in the bacis-cuts case when σ(gg → HH) = σ SM (gg → HH), see Table IV and Fig. 24. If a future e + e − linear collider and/or the 100 TeV pp machine are operating in the era of the high-luminosity LHC, combined efforts are desirable to determine the value of λ 3H uniquely [29]. B. CPC2: g S t , λ 3H , and g S tt This is the scenario that involves all scalar-type couplings in the triangle, box, and contact diagrams. The corresponding coefficients in Eq. (7) are c 1 , c 2 , c 3 , e 1 , e 2 , e 3 . Results at the detector level using the ATLAS template in Delphes v.3 are shown in Table II.
We first examine the cross section versus one input parameter at a time, shown in Fig. 6, while keeping the two parameters at their corresponding SM values. In the upper-left panel for σ/σ SM versus λ 3H , the lowest point occurs at λ 3H ≈ 2.5 when the interference term strongly cancels the triangle and box diagrams. Then the ratio increases from the lowest point on either side of λ 3H ≈ 2.5. Negative λ 3H s give constructive interference while positive λ 3H s give destructive interference. One may observe similar behavior when g S tt is varied as shown in the lower panel. Taking Since e 1 (s) > 0 and e 2 (s) < 0, we see that the contact diagram interferes constructively with the triangle diagram but destructively with the box diagram. The dominance of the box diagram leads to the totally destructive interference when g S tt > 0, resulting in the minimum at g S tt ≈ 0.5.
We show contours for the ratio σ/σ SM = 0.5, 1.5 in the plane of (λ 3H , g S t ) (upper-left), (λ 3H , g S tt ) (upper-right), and (g S t , g S tt ) (lower) in Fig. 7. The dashed lines denoted by −50% is for σ/σ SM = 0.5 and those by +50% for σ/σ SM = 1.5. In the upper-left panel in the plane of (λ 3H , g S t ), we show contours for g S tt = 0, 1. The g S tt = 0 is the same as the SM so that the contours are exactly the same as in Fig. 4, while with g S tt = 1, the contact diagram contributes significantly to the cross section, so that the contours shift more to negative (positive) λ 3H for positive (negative) g S t . In the upper-right panel, where we show the contours in the plane of (λ 3H , g S tt ), the g S tt negatively correlates with λ 3H because e 1 (s) > 0. In the lower panel, where we fix λ 3H = 1, somewhat nontrivial correlation between g S t and g S tt exists.
We use the same tools as in the CPC1 case to investigate the decay channel HH → γγbb We show the angular distributions ∆R γγ and ∆R bb between the two decay products of each Higgs boson in Fig. 9. The lines are the same as in Fig. 8. Similar to the CPC1 case, the higher the invariant mass, the more energetic the Higgs boson, and the smaller the angular separation between the decay products will be. Therefore, triangle diagram (red lines) has the widest separation, then followed by the box diagram (skyblue lines), and finally the contact diagram (green lines) has the smallest angular separation. We come up with the similar cuts as in the CPC1 case: ∆R larger or smaller than 2 to discriminate the triangle, box, and contact diagrams. We show in Table II the coefficients c 1 , c 2 , c 3 , e 1 , e 2 , e 3 such that the ratio of cross sections to the SM predictions can be given by Eq. (7).
Similar to what we have done for CPC1, we can make use of three simultaneous measurements of cross sections with basic cuts, ∆R > 2, and ∆R < 2. We show the region of parameter space that we can obtain using ∆R γγ (upper panels), ∆R bb (middle panels), and ∆R γγ and ∆R bb (lower panels) in the plane of (λ 3H , g S t ) in Fig. 10. Those on the left are for g S tt = 1 while those on the right are for g S tt = −1. Similarly, we show the parameter space in the plane of (λ 3H , g S t ) in Fig. 11 and in the plane of (g S t , g S tt ) in Fig. 12. for example, Higgs boson production cross section via gluon fusion and ttH production, both the real and imaginary parts of the coupling come in the form |g S t | 2 + |g P t | 2 , therefore one cannot tell the phase in the coupling 5 . The relevant coefficients for this CPV1 scenario are c 1 , c 2 , c 3 , d 1 , d 2 , d 3 , d 4 . They are shown in Table II at the detector simulation level (ATLAS).
We first show the variations of cross sections versus λ 3H with some fixed values of g S t and g P t in Fig. 13 6 . Also, the contours for the ratio σ/σ SM = 1 in the plane of (λ 3H , g S t ) (upperleft), (λ 3H , g P t ) (upper-right), and (g S t , g P t ) (lower) for a few values of the third parameter are shown in Fig. 14.
Similar to previous two scenarios, we use the same tools to analyze the decay channel HH → γγbb with parton showering and detector simulations. We show the invariant mass M γγbb and p Tγγ in Fig. 15, and the angular distributions ∆R γγ and ∆R bb between the two decay products of each Higgs boson in Fig. 16.
The terms by the triangle diagram (proportional to c 1 and d 1 in red and orange lines, respectively) give the widest separation among all the terms. The terms by the box diagram (proportional to c 3 and d 4 in skyblue and blue lines, respectively) give smaller angular separation. The full set of diagrams at the SM values (darkblue lines) and at g S t = g P t = 1/ √ 2 (grey lines) give similar results as the box diagram.
Similar to the CPC1 and CPC2 cases, we use the cuts ∆R larger or smaller than 2 to discriminate the triangle and box diagrams. We show in Table II the coefficients c 1 , c 2 , c 3 , d 1 , d 2 , d 3 , d 4 , which are relevant ones in the CPV1 scenario, such that the ratio of cross sections to the SM predictions can be obtained by Eq. (7). We show the region of parameter space that we can obtain using ∆R γγ , ∆R bb , and ∆R γγ and ∆R bb in the plane of (g S t , g P t ) in Fig. 17, in the plane of (λ 3H , g P t ) in Fig. 18, and in the plane of (λ 3H , g S t ) in Fig. 19. D. CPV2: g S tt , g P tt , and λ 3H Here we study another CP-violating scenario with the CP-even and CP-odd components of the ttHH coupling with the top-Yukawa couplings g S t and g P t at fixed values. Note that the CP-odd coupling g P tt only appears in the product with g P t or by itself squared. In this case, all the coefficients are relevant and they are shown in Table II at the detector simulation level (ATLAS).
Similar to previous scenarios we used HH → γγbb with parton showering and detector simulations. We use the cuts ∆R larger or smaller than 2 to discriminate the triangle and box diagrams. Using the coefficients presented in Table II, the ratio of cross sections to the SM predictions can be obtained by Eq. (7). We show the region of parameter space that we can obtain using ∆R γγ , ∆R bb , and ∆R γγ,bb in the plane of (g S tt , g P tt ) for fixed λ 3H = 1, g S t = 1, and g P t = 0 in Fig. 20; and similarly in the plane of (λ 3H , g P tt ) for g S t = 1, g P t = g S tt = 0 in Fig. 21; and finally in the plane of (λ 3H , g S tt ) for g S t = 1, g P t = 0, and g P tt = 0.5 in Fig. 22;

E. 100 TeV Prospect
All the results represented for the 14 TeV run were obtained by manipulating the coefficients represented in Table II. We represent the coefficients c 1,2,3 , d 1,2,3,4 , e 1,2,3 , f 1,2,3,4 for the 100 TeV pp machine in Table IV. Just for illustrations, we show the distributions of the invariant mass M γγbb and angular separation ∆R γγ for the CPC1 case at the 100 TeV machine in Fig. 23. We found that the behavior of the distributions at 100 TeV is very similar to those at 14 TeV. Therefore, the kinematic regions of interests separated by ∆R can be taken to be the same as 14 TeV. We can make simultaneous measurements of cross sections at 100 TeV pp machine to isolate the Higgs trilinear coupling. We show the sensitivity regions of parameter space in the CPC1 case at the 100 TeV pp machine in Fig. 24.
The regions are very similar to those in 14 TeV, though not exactly the same. Sensitivity reach for each coupling in other cases can be obtained by similar methods with the assumed luminosity.

V. CONCLUSIONS
In this work, we have studied the behavior of Higgs-boson pair production via gluon fusion at the 14 TeV LHC and the 100 TeV pp machine. We have performed an exploratory study with heavy degrees of freedom being integrated out and resulting in possible modifications of the top-Yukawa coupling, Higgs trilinear coupling, and a new contact ttHH coupling, as well as the potential CP-odd component in the Yukawa and contact couplings. We have identified useful variables -the angular separation between the decay products of the Higgs bosonto discriminate among the contributions from the triangle, box, and contact diagrams. We have successfully demonstrated that with three simultaneous measurements of the Higgspair production cross sections, defined by the kinematic cuts, one can statistically show a nonzero value for the Higgs trilinear coupling λ 3H if we can measure the cross sections with less than 25% uncertainty. This is the key result of this work.
We also offer the following comments with regards to our findings.
1. The triangle diagram, which contains an s-channel Higgs propagator, does not increase as much as the box diagram or the contact diagram with the center-of-mass energy √ŝ . This explains why the opening angle (∆R γγ or ∆R bb ) in the decay products of each Higgs boson is a useful variable to separate between the triangle and the box diagram. Thus, it helps to isolate the Higgs trilinear coupling λ 3H .
2. The contact diagram contains a dim-5 operator ttHH, which actually breaks the unitarity at about √ŝ ∼ 17.6/g S tt TeV. This implies that it could become dominant at high invariant mass.
3. Suppose we take a measurement of cross sections, we can map out the possible region of parameter space. Since in different kinematic regions the regions of parameter space are mapped out differently, such that simultaneous measurements can map out the intersected regions. With measurement uncertainties less than 25% one can statistically show a nonzero value for the Higgs trilinear coupling, and also obtain the sensitivity regions of λ 3H : 0.3 |λ 3H | 2.6. for σ/σ SM = 1. 4. We found that the behavior of the distributions of the invariant mass M γγbb and angular separation ∆R γγ or ∆R bb at 14 TeV are very similar to those at 100 TeV. We can then use the same method as in 14 TeV to isolate the Higgs trilinear coupling.
5. It is difficult, if not impossible, to determine the Higgs trilinear coupling uniquely at the LHC and 100 TeV pp machine even in the simplest case assuming very high luminosity and precise independent input for the top-Yukawa coupling. We suggest to combine the LHC results with information which can be obtained at a future e + e − linear collider.
6. If the couplings deviate from their SM values, the Higgs-boson pair production cross section can easily increase by an order of magnitude. For example, in the CPC2 case, σ/σ SM > 10 for λ 3H > 9 or < −4 when g S t = 1 and g S tt = 0, g S t > 1.7 or < −1.3 when λ 3H = 1 and g S tt = 0, and g S tt > 2.6 or < −1.4 when λ 3H = g S t = 1: see Fig. 6. The cross section larger than the SM prediction may reveal the new physics hidden behind the SM and we can have better prospect to measure the Higgs self coupling at the LHC.  Table III, we show the SM cross sections for pp → HH → γγbb at the 14 TeV LHC with and without angular-separation cuts. Note that the cross section before applying any cuts is about 0.09 fb and it becomes 0.005 fb after applying the basic cuts. In the region of ∆R γγ > 2(< 2), the cross section is 0.0013 fb (0.0038 fb) where it is dominated by the triangle (box) diagram. The ratio is about 1 : 2.8. We also show the cross sections for the 100 TeV pp machine, and the corresponding ratio is about 1 : 3.7. It shows the fact that the triangle diagram is more suppressed because of the s-channel Higgs propagator at higher energy. In the regions of ∆R bb larger and smaller than 2, the ratios are 1 : 5.7 and 1 : 8.2 at the 14 TeV LHC and the 100 TeV pp machine, respectively.
As we have promised, we are going to comment on the HH → τ + τ − bb decay mode. This mode has the obvious advantage of a larger branching ratio than the γγbb mode, but the identification efficiency and momentum measurements of τ leptons are much weaker than photons. In Table III, we show the SM cross sections for pp → HH → τ + τ − bb at the 14 TeV LHC with and without angular-separation cuts. Taking into account the branching TABLE II. 14 TeV LHC: The coeffcients for the ratio of the cross sections σ(gg → HH)/σ SM (gg → HH) as in Eq. (7) with and without the angular-separation cuts of ∆R γγ > or < 2; ∆R bb > or < 2; and ∆R γγ and ∆R bb both > 2 or < 2. The relevant coefficients for the CPC1 scenario are c 1 , c 2 , c 3 ; those for the CPC2 scenario are c 1 , c 2 , c 3 , e 1 , e 2 , e 3 ; and those for the CPV1 scenario are c 1 , c 2 , c 3 , d 1 , d 2 , d 3 , d 4 . All the coefficients are involved in the CPV2 case. Results are at the detector level using the ATLAS template in Delphes v.3.

√
s : 14 TeV   The SM cross sections for the process pp → HH → γγbb with various angularseparation cuts on ∆R bb ,γγ at the 14 TeV LHC (second column) and at the 100 TeV pp machine (third column). The last column shows them for the process pp → HH → τ + τ − bb at the 14 TeV LHC. We have taken account of the SM NLO cross section σ SM (pp → HH) 34 fb, the Higgs branching fractions, and both the photon and b-quark reconstruction efficiencies. The p T dependence of b-tagging efficiency is considered and 0.5 is taken for the τ -tagging efficiency. Also considered is the mis-tagging probability of P j→τ = 0.01. Results are at the detector level using the ATLAS template in Delphes v.3.

Cuts
SM -     The 25% and 50% sensitivity regions bounded by three measurements of cross sections with basic cuts, ∆R γγ > 2, and ∆R γγ < 2 (the upper-left panel); with basic cuts, ∆R bb > 2, and ∆R bb < 2 (the upper-right panel); and basic cuts, ∆R γγ , ∆R bb > 2, and ∆R γγ , ∆R bb < 2 (the lower panel). We assume that the measurements agree with the SM values with uncertainties of 25% and 50%, respectively.  The 25% and 50% sensitivity regions in (λ 3H , g S t ) bounded by three measurements of cross sections with basic cuts, ∆R γγ > 2, and ∆R γγ < 2 (the upper panels); with basic cuts, ∆R bb > 2, and ∆R bb < 2 (the middle panels); and basic cuts, ∆R γγ , ∆R bb > 2, and ∆R γγ , ∆R bb < 2 (the lower panels). The left panels are for g S tt = 1 while those on the right are for g S tt = −1. We assume that the measurements agree with the SM values with uncertainties of 25% and 50%, respectively. The 25% and 50% sensitivity regions in the plane (g S tt , g P tt ) (with fixed λ 3H = 1, g S t = 1, and g P t = 0) bounded by three measurements of cross sections with basic cuts, ∆R γγ > 2, and ∆R γγ < 2 (the upper-left panel); with basic cuts, ∆R bb > 2, and ∆R bb < 2 (the upper-right panel); and basic cuts, ∆R γγ , ∆R bb > 2, and ∆R γγ , ∆R bb < 2 (the lower panel). We assume that the measurements agree with the SM values with uncertainties of 25% and 50%, respectively.