Production of $HHH$ and $HHV(V=\gamma,Z)$ at the hadron colliders

We consider the production of two Higgs bosons in association with a gauge boson or another Higgs boson at the hadron colliders. We compute the cross sections and distributions for the processes $ p p \to H H H $ and $H H Z$ within the standard model. In particular, we compute the gluon-gluon fusion one-loop contributions mediated via heavy quarks in the loop. It is the leading order contribution to $ p p \to H H H $ process. To the process $ p p \to H H Z $, it is next-to-next-to-leading-order (NNLO) contribution in QCD coupling. We also compare this contribution to the next-to-leading-order (NLO) QCD contribution to this process. The NNLO contribution can be similar to NLO contribution at the Large Hadron Collider (LHC), and significantly more at higher center-of-mass energy machines. We also study new physics effects in these processes by considering $ttH, HHH, HHHH, HZZ$, and $HHZZ$ interactions as anomalous. The anomalous couplings can enhance the cross sections significantly. The $gg \to HHH$ process is specially sensitive to anomalous trilinear Higgs boson self-coupling. For the $gg \to HHZ$ process, there is some modest dependence on anomalous $HZZ$ couplings.


Introduction
The CMS and ATLAS Collaborations have been collecting data at the Large Hadron Collider (LHC) for several years [1, 2]. Their major discovery has been of much-anticipated Higgs boson in 2012 [3,4]. There are many reasons to go beyond the standard model (SM). Since 2012, the search has been going on for any hint for physics beyond the standard model. There have been a number of anomalies that were suspected at different times [5][6][7], however none has stood test of times. There is no evidence of any signal for beyond the standard model scenarios. A number of popular scenarios involving supersymmetry and large extra dimensions are getting severely constrained [8][9][10][11]. Various processes are being analyzed for any hint of a new scenario [1,2]. In such a situation, exploration of rare processes and radiative corrections provide promising avenues to explore. We should note that some other experiments, e. g. LHCb, have reported some unexplained phenomena [12].
At a hadron collider, as centre-of-mass energy increases, so does gluon-gluon luminosity. Therefore, at the LHC and at future probable hadron colliders, the gluons initiated processes would play important role. In this paper, we are considering few such processes. The processes that we consider occur at oneloop level. The gluon-gluon initiated 2 → 3 one-loop processes have been considered in the literature. First full calculation of gg → γγg was presented in [13]. Many other authors have also computed the contribution of the gluon-gluon initiated processes on 'multi-bosons +jets' [14][15][16][17][18][19][20][21][22][23][24][25][26]. These processes are important ingredients for the NLO QCD predictions for gg → BB(B = γ, Z, W, H) processes [27][28][29][30][31]. These different calculations use different reduction techniques, different packages for computing scalar integrals, and overall different philosophy for the computation. We use our own tensor-reduction code, and have developed a comprehensive package for such calculations.
In this paper, we consider the processes gg → HHH, HHγ, and HHZ, and their contribution to hadron level processes -pp → HHH, HHγ, and HHZ. We presented some preliminary results on these processes along with other Higgs processes in 2 → 3 category in a conference proceeding [32]. Triple Higgs production via gluon fusion has been studied by many authors [33][34][35][36][37][38]. Plehn and Rauch [33] considered the possibility of measuring quartic Higgs boson coupling in the HHH production. In the reference [35], authors have considered the QCD correction to the HH and HHH production in a Higgs effective field theory approach. The authors in [36][37][38], have considered the possibility of observing the HHH production at a 100 TeV collider including new physics effects. In Ref. [25] standard model cross sections for a number of loop-induced gluon fusion processes including gg → HHH, HHZ are reported. While our work on pp → HHH has some overlap with these papers as discussed below, our detailed study of pp → HHZ process in SM and beyond is new and being presented here for the first time in the literature. We have also computed using different tools and looked at different aspects of the processes.
The exploration of the production of HHH is important, as it is one of the very few processes where quartic Higgs boson coupling is involved. This process may allow the direct measurement of this coupling. With the measurement of self-couplings of the Higgs boson, one can confirm the form of the Higgs potential. Unlike the HHH production, the process pp → HHZ gets contribution from the tree-level processes. One can also compute next-to-leading order (NLO) QCD corrections to this tree-level process [39]. In this paper, our focus is on gluon-gluon annihilation contribution, but we also compare it with the LO and NLO contributions. The gg → HHZ contribution can be thought of as next-to-next-to-leading order (NNLO) corrections [40]. The production of HHZ is important, as it involves HHH and HHZZ couplings. It is also a background to triple Higgs production process.
In search for new physics scenarios, the use of anomalous interactions may play an important role. This is a model-independent approach which can systematize the search. They may point towards a model that would be suitable to go beyond the standard model. In this paper, we consider possible modification of standard model interactions, inspired by dimension-six operators in effective field theory approach. In particular, we consider the modifications of ttH, HHH, HHHH, HZZ, and HHZZ interactions. We study the effects of these anomalous interactions on the production cross sections and on kinematic distributions in gg → HHH, HHZ processes. Note that in Refs. [36,37], the HHH process is studied considering only SM-like deviations in trilinear and quartic couplings. We have in addition considered derivative couplings which have different effects on some distributions as we demonstrate. In Ref. [38], all the CP-even dimension-six operators relevant to HHH process has been considered. In present work, our approach towards new physical effects in HHH process is more phenomenological and we have considered modifications to only those couplings which are present in the standard model at tree-level. For example, we do not consider ttHH, ttHHH, ggH and ggHH interactions. However, for ttH coupling we have considered both CP-even and CP-odd anomalous interactions. Similar approach is taken to study new physics effects in HHZ process.
The paper is organized as follows. In the next section 2, we discuss these various processes and gluon-gluon annihilation contributions to them. In the section 3, the anomalous contribution to various relevant vertices is discussed. In the section 4, we provide details on the method of calculation and numerous checks. Our numerical results for SM and BSM are presented in section 5. In the last section 6, we present our conclusions.

Processes
We consider the gluon-gluon annihilation contribution to the following processes: Let us first consider the process p p → H H γ. In the standard model, at tree-level, this process has vanishingly small cross section due to very small light-quark and Higgs boson coupling. At one-loop, there is no contribution to this process form gluon-gluon fusion channel. Using Furry's theorem, one can see that all diagrams contributing to this process add up to zero. The p p → H H H process occurs primarily via g g → H H H and leading contribution comes from one-loop diagrams. Because of the very small light-quark and Higgs boson couplings, the tree level process qq → HHH makes very small contribution. Therefore, we only consider the contribution of the g g → H H H process. This is a one-loop process. There are 24 pentagon, 18 box, and 8 triangle diagrams contributing to the process for each quark flavor in the loop. There are some more diagrams in the box and triangle categories at one-loop level, but their contributions are zero, as they do not conserve color charge (Tr(λ a )=0). In the loop, a light-quark would not contribute due to very small Yukawa coupling; so we keep only diagrams with a top-quark and a bottom quark. We do not need to numerically compute all the diagrams separately, as many diagrams are related to one another by charge conjugation symmetry, or crossing. We can compute all the diagrams using six prototype diagrams as shown in Fig. 1. By permuting the legs in the prototype diagrams, and using charge conjugation (Furry's theorem), all the amplitudes can be calculated. Out of 24 (=4!) pentagon diagrams, we had to numerically calculate only 12 diagrams as other 12 diagrams can be related to the previous diagrams by Furry's theorem. Out of these 12 pentagon diagrams, 6 (=3!) diagrams can be obtained from PENTA1 prototype diagram by permuting the Higgs bosons; other 6 diagrams can be obtained similarly from PENTA2 prototype diagram. Similarly, out of 18 box diagrams, we need to numerically compute only 9 diagrams, and from these the rest can be found out using Furry's theorem only. Out of these 9 box diagrams, 6 can be calculated from BOX1 prototype diagram by permuting Higgs bosons in six different ways. Other 3 box diagrams can be obtained by the permutation of external Higgs boson in BOX2 prototype diagram. Out of 8 triangle diagrams, we need to numerically compute only 4; TRIANGLE1 prototype diagram gives one of these and TRIANGLE2 prototype diagram gives the rest 3 diagrams by permuting external Higgs bosons. The p p → H H Z can occur at tree level through quark-antiquark annihilation. The cross section is large enough for it to be observable with high luminosity option. One can also compute NLO QCD corrections to this process easily. Our focus will be contribution of gluon-gluon fusion process: g g → H H Z which occurs at the one-loop level. Formally this contribution is of NNLO order. However, as we shall see, at the LHC, this contribution can be of the same order as the NLO corrections. Due to enhanced gluon-gluon luminosity at larger center-of-mass energies, this NNLO correction will even dominate over NLO correction at higher energy hadron colliders. There are 24 pentagon, 18 box, and 8 triangle diagrams contributing to the process g g → H H Z. Like in the case g g → H H H, here also the diagrams are divided into many prototype classes as shown in Fig. 2 . Here, as before, out of 24 pentagon, 18 box, and 8 triangle diagrams, we need to numerically compute only 12 pentagon, 9 box, and 4 triangle diagrams respectively. Out of these 12 pentagon diagrams, each of PENTA1 and PENTA2 prototype diagrams gives four diagrams, and each of PENTA3 and PENTA4 prototype diagrams gives two diagrams. In the case of box diagrams, BOX1, BOX2, BOX3, and BOX4 prototype diagrams give 4, 2, 2, and 1 diagrams respectively. There are three prototype triangle diagrams -TRIANGLE1, TRIANGLE2, and TRIANGLE3; these prototype diagrams give 1, 1, and 2 diagrams respectively. Note that due to Furry's theorem only the axial-vector part of the Z boson coupling with the quarks in the loop contributes. Feynman diagrams in Fig. 1 and Fig. 2 have been made using JaxoDraw [41]. Computation of one-loop diagrams in gg → HHH, HHZ processes are described in section 4.

Anomalous interactions of Higgs boson
There are a number of arguments for going beyond the standard model. In absence of any new resonance at the LHC, one way to consider the effects of beyond-the-standard-model scenarios is to consider the possible modification of standard model vertices. We are mainly interested in anomalous couplings of the Higgs boson which would affect the processes under consideration. These include, ttH, HZZ, HHZZ, HHH and HHHH interactions. Some of them, for example, ttH and HZZ are already constrained by the existing LHC data [42]. The trilinear Higgs self coupling is very weekly constrained by the data [43] and couplings HHZZ and HHHH are unconstrained at present. In the following, we consider most general interaction Lagrangian incorporating the BSM physics which would lead to deviations in the Higgs couplings of our interests.

AnomalousttH Vertex
In the standard model, the top quark couples with the Higgs boson via the Yukawa coupling. This leads to a scalar ttH coupling. The most general vertex forttH interaction can be parametrized as, Here v is electroweak symmetry breaking scale and it is approximately 247 GeV. In the standard model, κ t =κ t = 0. We use following bounds for κ t andκ t [44,45]: This vertex contributes to both HHH and HHZ production. As we shall see, the scaling of the scalar coupling, as parametrized by κ t can change the cross section significantly.

Anomalous HHH and HHHH Vertices
After the discovery of the Higgs boson, one of the important task is to determine the form of the Higgs potential. As discussed before, one of the characteristic feature of the Higgs potential in the standard model is specific form of Higgs boson self-couplings. The Higgs boson self-interactions can be expressed in terms of the anomalous couplings as In the SM, g 4H just scale the trilinear and quartic Higgs boson self-couplings, respectively. Both of these couplings occur in HHH production. In the HHZ production, only trilinear coupling occurs. As we shall see, both the processes are sensitive to modification of trilinear coupling. These couplings are poorly determined, as they occur in processes with small cross section. It may take a decade or more before a serious bound can be put on these couplings. For illustration, we take these parameters in the range between -1.0 to 1.0.

Anomalous HZZ and HHZZ vertices
These two vertices occur in the process gg → HHZ. While HZZ vertex occurs in other processes, like pp → HZ, the vertex HHZZ occurs mainly in processes involving double Higgs boson production, so is poorly constrained. The most general HZZ interaction that can be written has following form, In the above g is the coupling parameter of the SU (2) L group and c W = cosθ W , θ W being the weak angle. The SM HZZ interaction corresponds to g (0) HZZ = 0. This interaction has derivative couplings which lead to momentum dependence. This is unlike standard model.
The following are the various bounds on the anomalous coupling parameters [46] −0.10 ≤ g (0) There is also modification of HHZZ coupling. This modification only scales the standard model interactions: In the SM, g HHZZ = 0. In absence of any available bound for this coupling, we allow the parameter g (0) HHZZ to vary between -0.1 and 0.1.
The anomalous interactions mentioned above are well motivated within the framework of an effective field theory in which the new physics effects are parametrized in terms of higher dimensional operators. These operators are constructed from the SM fields and respect the symmetries of the SM. A complete list of independent dimension-six operators is now available [47][48][49]. The anomalous couplings introduced above are related to the Wilson coefficients of these operators [50][51][52][53][54].
The Feynman rules for the anomalous Higgs vertices are listed in the appendix A. As we shall see, in the allowed range of the parameter values, the contribution of the anomalous vertices can be important in our processes. Note that even in the presence of these anomalous couplings the amplitude for gg → HHγ process does not receive any contribution.

Calculation and Checks
As discussed in the section 2, we compute prototype pentagon, box, and triangle diagrams. Then by using crossing and Furry's theorem, we can compute rest of the diagrams. All the diagrams have a fermion loop, so to calculate the amplitude we need to compute trace of a string of gamma matrices. We calculate the traces of the prototypes diagrams using FORM [55]. The process gg → HHZ includes ttZ coupling, which has both vector and axial vector parts. The presence of γ 5 requires special care, due to the potential presence of anomalies. We have handled this situation in two different ways. Since the process is free from UV divergences, we can take trace in four dimensions. We have also calculated the trace in n dimensions using Larin's prescription for γ 5 [56,57]. Both methods, in the end give same results when the contributions from both the top and bottom quark loops are considered. When we include ttH anomalous pseudo-scalar coupling, the trace will include more γ 5 matrices for both HHH and HHZ production.
In the first step, we use FORM to take the trace and write the amplitude in terms of tensor and scalar integrals. We use an in-house package OVReduce to reduce tensor integrals to lower-point tensor integrals and scalar integrals in dimensional regularization . This package is based on the methods of Oldenborgh and Vermaseren [21,58]. After the reduction, all that we need to do is to compute various scalar integrals of box, triangle, and bubble types. To compute these scalar integrals, we use OneLOop library by Andreas van Hameren [59]. It uses dimensional regularization to regulate UV and IR divergences. For the pentagon scalar integrals, we use van Neerven-Vermaseren technique [21,60]. As the amplitudes are quite large and complicated, we first compute helicity amplitude for a process numerically before squaring it. The Monte Carlo integration over the entire phase space has been performed using VEGAS code [61] as implemented in AMCI (Advanced Monte Carlo Integration) package [62] . AMCI implements a parallel version of VEGAS algorithm which makes use of Parallel Virtual Machine (PVM) software [63].
While performing multi-leg one loop calculations, the issue of numerical instability comes in for certain phase space points. As the number of such points are not large, and contribution of these phase space points is not expected to be large, as is the practice, we exclude these phase space points by setting a suitable upper bound on the amplitude-squared 1 . This upper bound is chosen after finding out possible values the amplitude-squared can have by running the code. This upper bound is increased until we hit the unstable phase space point. The cross section remains stable and does not change, even when this upper bound is increased by several orders of magnitude. For the HHH process, we do not come across any unstable phase space point. For the gg → HHZ process, the number of unstable points are well below 0.002%.
We have performed many checks to verify the correctness of the amplitude for each process. These checks include verification of cancellation of UV & IR divergences and gauge invariance. Because of the presence of the Higgs boson in the final state, we have only top and bottom quarks in the loop. In both the processes, all pentagon, box, and triangle diagrams are therefore separately IR finite. Overall amplitude at any phase space point is UV finite. Each pentagon diagram is UV finite, which is again as expected as simple power counting reveals it. Individual box diagrams are also UV finite, since one will have at most two (three)-tensor box integrals in gg → HHH(HHZ). Each triangle diagram is also UV finite. In both these processes, we have verified that in the large top quark mass limit the amplitude becomes constant implying non-decoupling of top quark in m t → ∞ limit [64].
To check gauge invariance, we replace the polarization vector of a gluon with its momenta in the amplitude. In gg → HHH, the overall amplitude has been checked to be gauge invariant with respect to both the gluons. We find that each triangle diagram is individually gauge invariant, while each pentagon and box diagram is not. However, all pentagon diagrams together, and all box diagrams together are gauge invariant. Here interesting point is that the pentagon, box, and triangle diagrams are separately gauge invariant with respect to the gluons 2 . So it may be tempting to use only one class of diagrams to compute the cross section. However as will see below, it can lead to serious errors. Here we have done the calculation only in four dimension. The amplitude is found to be gauge invariant for both gluons in the presence of pseudo-scalar coupling even when we consider only one quark.
We have calculated the gg → HHZ amplitude treating γ 5 in four-dimension and in n-dimension using Larin's prescription [56]. In four-dimension, if we consider only one quark, each triangle diagram is gauge invariant with respect to one gluon but not for the other. All the triangle diagrams taken together are also not gauge invariant for the other gluon. However, if we consider both top and bottom quarks in the amplitude, each triangle diagram is gauge invariant for both gluons. This is related to the quantum anomaly associated with the axial vector coupling of Z boson with fermions. None of the pentagon or box diagram is individually gauge invariant for any gluon. However, referring to Fig. 2, all the pentagon diagrams taken together, or all the diagrams of BOX1 and BOX2 classes together, or all the diagrams of BOX3 and BOX4 classes together are gauge invariant for both gluons. In n-dimension, all the pentagon diagrams together, all the BOX1 and BOX2 diagrams together, all the BOX3 and BOX4 diagrams together, and each triangle diagram are gauge invariant with respect to the both gluons for top and bottom quarks separately. In any case, to remove the anomaly associated with the chiral current of Z boson the contributions from both the top and bottom quarks in the loop must be included.

Numerical results
In our computation of cross sections and distributions we have used: p H,Z T > 1 GeV, |y H,Z | < 5. The cut on the p T of the Higgs boson and the Z boson (of 1 GeV) is there to reduce the number of phase space points which introduce numerical instability. We have checked and it can be understood from the p T distributions presented below that the effect of removing the cut is negligible. The results for gluon fusion processes are obtained using cteq6l1 parton distribution functions [65], and using µ R = µ F = √ŝ (partonic center-of-mass energy ) as the renormalization and factorization scales. We have also included uncertainties in the results for the LHC by varying the renormalization and factorization scales by a factor of 2. The scale uncertainties are listed in the table 3 .

5.1
The process pp → HHH In Table 1, we present the cross section for this process at the LHC center-of-mass energy, and at other proposed hadronic colliders. The cross section at 13 TeV center-of-mass energy is 32.0 attobarn. So as of now only 2-3 such events may have been produced at the LHC. With even expected 3 ab −1 luminosity, there will be only about 100 events. At 100 TeV, thanks to a larger gluon flux, the cross section will be 100 times larger. Note that these cross sections suffer from large scale uncertainty (-22% to 31% at 13 TeV) which is typical to gluon fusion processes. These uncertainties are large due to the significant dependence of the strong coupling constant, α s (µ), on the renormalization scale.  In the Table 2, we have displayed the values of cross sections when only pentagon, or box, or triangle type diagrams are considered. These categories of diagrams are separately gauge invariant with respect to the gluons. So if we wish to estimate cross section by only keeping a category of diagrams, we will make an order of magnitude error. This is because there is large destructive interference among these categories of diagrams. If for simplicity, we include only triangle diagrams in the calculation, we will underestimate the cross section, while inclusion of only box or pentagon diagrams will overestimate the cross section. For 13 TeV centre-of-mass energy we see that the total cross section is about 32.0 attobarn, whereas penta, box, and triangle class contribute 94.4, 53.6, and 3.5 attobarn respectively.   This shows that there is strong destructive interference between the different categories of diagrams. At higher center-of-mass energies the interference effect becomes smaller. If one takes the Higgs effective field theory approach, there would be inclusion of triangle type diagrams only, and the cross section would be underestimated. In Fig. 3, we have plotted the contribution of various category of diagrams with respect to the p T of the leading Higgs boson at √ s = 13 TeV and 100 TeV. We see that it is the pentagon type diagrams, which give harder Higgs boson, than other categories of diagrams. Interference kills such events and the p T peaks between 130 and 160 GeV. In Fig. 4, we have plotted a number of physical quantities involving leading, next-to-leading, and next-to-next-to-leading Higgs bosons arranged according to their transverse momenta. As would be expected leading Higgs boson's p T is harder and peaks around 140-160 GeV. Softest Higgs boson p T is mostly around 50 GeV with a large tail. However all three types are produced mainly centrally. The leading and next-to-leading Higgs bosons are produced more backto-back than other Higgs bosons. Two of the softer Higgs bosons are produced closer to each other. The masses of the Higgs boson pair also show expected behavior. The mass of the two larger p T Higgs bosons has a peak about 375 GeV, while for the two softer Higgs boson, it is near the twice of the mass   of the Higgs boson. The invariant mass of the three Higgs bosons peak around 550 GeV. At higher center-of-mass energy, 100 TeV, the behavior of physical quantities is largely the same, so we have not given separate plots.
In Fig. 5, we have plotted the ratios of the cross section with various anomalous couplings and the standard model cross section. We have plotted for the range of parameters mentioned in section 3, except forκ t for which we have doubled the range. We see that the cross section is not sensitive to the pseudo-scalar ttH coupling and modification of quartic Higgs boson coupling. It is not surprising because of three Higgs boson coming out of same vertex, requiring the fourth Higgs boson to be far off-shell. However, it is sensitive to the scaling of the scalar ttH coupling. Cross section can significantly change with the change in this coupling -by a factor of 3-4. The cross section is also sensitive to the sign of this coupling. Interestingly, the cross section is sensitive to trilinear Higgs boson coupling too. Cross section can increase by an order of magnitude by the change of derivative coupling, g (1) 3H . This coupling is momentum dependent and will have larger effect at higher center-of-mass energy machine. The scaling of the coupling, g (0) 3H , can also change the cross section by about a factor of 3. If one could detect HHH events at a future collider, one can probe trilinear coupling easily. Looking at Fig. 6, we see that the scaling of the coupling changes the low p T (of the leading Higgs boson) events more; while derivative coupling tends to also change large p T events. Therefore, one can probe both couplings by focusing on low p T events in one case and higher p T events in another case. The p T distributions for the sub-leading Higgs boson and mass distributions of two or three Higgs bosons show similar sensitivity. Note that the cross-section is symmetric forκ t .
Let us now consider the possibility of observing the production of the three Higgs bosons, HHH. At the LHC, the cross section is of the order of 32 ab. It leads to too few events with even the highest possible integrated luminosity of 3 ab −1 to be observable in the sea of the background. However we have seen that this process is sensitive to trilinear Higgs boson self-coupling. Such anomalous couplings can enhance the cross section by a factor of 3 to 8. Even with this enhancement, once we include branching ratios, kinematic cuts, tagging, and other efficiencies, there will be too few events to be visible. Given that it is one of the few processes to help determine the quartic Higgs boson selfcoupling, it will be worthwhile to look for this process. At a 100 TeV machine, it could be possible. The cross section there is about 3100 ab. With large enough data, this process might be observable. The process pp → HHH will give rise to 'multileptons' and 'few leptons with jets' signatures. There will be irreducible background from ZZZ production, and an array of reducible backgrounds depending on the signature. The ZZZ production cross section is about 136 fb which can be controlled if we include branching ratios and construct Higgs boson masses. In the case of reducible backgrounds there will be background from a top-pair production with jets or vector bosons, multi vector boson production with jets, and multijets. To tame these backgrounds, one may require tagging of bottom jets and tau jets. Given the severity of the background, one may still need a multivariate analysis. In [36], authors have studies the the channel HHH → bbbbγγ in the standard model and its simple deformation with marginal success. More recently the authors of [37] have studied 'bbbbτ τ ' channel and authors in [38] have similarly studied 'bbl + l − + 4 jets' channel with modest success. Therefore, a modification of the interactions that will enhance the signal significantly and improved search strategies will be needed to determine Higgs boson self-couplings.

The process pp → HHZ
Unlike the pp → HHH process, this process can occur at the tree level. In this section, we will mainly focus the NNLO contribution to this process that occur via gg → HHZ process. We have estimated tree level value and one-loop QCD corrections, i. e., NLO contribution using MadGraph5 aMC@NLO [66]. We will see that NNLO contribution is comparable to NLO contribution at 13 TeV. NNLO contribution becomes even more important as the center-of-mass energy increases, and can become comparable to the tree level value.   Table 3: A comparison of different perturbative orders in QCD coupling contributing to pp → HHZ hadronic cross section at √ s = 8, 13, 33, and 100 TeV. We also calculate ratios R 1 , R 2 and R 3 which quantify the GG contribution with respect to the QQ(LO) and QQ(NLO) contributions.
In Table 3, we have given LO, NLO, and NNLO contribution to this process at different center-ofmass energies. We have used cteq6l1 parton distribution for LO and NNLO calculation and cteq6m for the NLO calculation [65], withŝ as renormalization/factorization scales for NNLO calculation and sum of transverse mass for LO and NLO calculation (in MadGraph5 aMC@NLO). Uncertainties are estimated by using other cteq6 parton distributions and changing the renormalization/factorization scales by a factor of 2. At the LHC, the leading order contribution is about 237 ab. NLO corrections are about 24%, and add 57.8 ab to the LO value. NNLO corrections are about 18%, and add 42.3 ab to the LO value. The cross section including LO, NLO, and NNLO values is about 336.8 ab, leading to about 1000 events at the maximal proposed integrated luminosity. As in the case of HHH production, the value of NNLO contribution becomes more significant as center-of-mass energy increases. At 100 TeV, NLO correction is only 13%, while NNLO contribution is 81% of the LO value due to an increase in the gluon flux. We see that NNLO contribution approaches the LO value. It is however not alarming, only more useful. NNLO contribution is from gluon-gluon annihilation, while LO result is from quarkantiquark scattering. At 13 TeV, the scale uncertainties for GG channel are in the range of -21% to 31%. The reason for this large uncertainty and its decrease with center of mass energy is same as that for gg → HHH process. Uncertainties in the QQ processes are smaller, as these are primarily electroweak processes. √ s (TeV) 8 Table 4: SM contribution of pentagon, box, and triangle diagrams to the total cross section in gg → HHZ at different collider center-of-mass energies, displaying a destructive interference effect.  Table 4 demonstrates interference effect in the process HHZ. Unlike HHH process, here triangle diagram contribution is larger and interference effect is more severe. If we keep only one class of diagrams which are gauge invariant with respect to the gluons, we will overestimate the cross section by an order of magnitude at 13 TeV. We see that the total cross section is 42.3 attobarn, whereas penta, box, and triangle class contribute 148.1, 434.7, and 475.6 attobarn respectively. This shows there is a strong destructive interference between the different class of diagrams. The destructive interference effect becomes stronger at higher √ s. In Fig. 7, we have plotted the contribution of the individual class of diagrams with respect to the p T of the leading (in p T ) Higgs boson. Here effects are different. Triangle and box diagrams have larger cross section and contribute more to higher p T events. The interference effects seem to cut off large p T contribution of the individual class of diagrams. The p T of the leading Higgs boson, after the interference, shifts to lower values and peaks around 120 GeV.
In Fig. 8, we have plotted a few physical quantities, as in the case of HHH production for 13 TeV LHC. The p T distribution of the leading Higgs bosons and Z boson are similar. Both peaks around 110 GeV. As would be expected, the leading Higgs boson's p T is harder than next-to-leading Higgs  Table 3). compare the GG(LO) and QQ (LO and NLO) contributions in kinematic distributions for p T (H 1 ) and p T (Z) at 13 TeV LHC. In Fig. 9 we see that contribution of GG channel is characteristically different from that of QQ channel. The GG channel gives rise to softer events in p T (H 1 ), while harder events in p T (Z). These features remain true at higher centre-of-mass energies. In Fig. 10, we give the p T (H 1 ) and p T (Z) distributions combining QQ(NLO) and GG(LO) channels at 13 and 100 TeV. At the level of total cross section, the GG contribution is about 14% with respect to the QQ(NLO) contribution at 13 TeV. However, in the distributions the GG contribution can reach more than 20 % in certain bins. Similarly, at 100 TeV collider, although the contribution of the GG process to the total cross section is about 72% of the QQ(NLO) contribution (see Table 3), in the case of p T (H 1 ) between 100 − 200 GeV, the enhancement due to GG channel could be more than 100%.
The process gg → HHZ has four different couplings -ttH, HZZ, HHH, and HHZZ. In Fig. 11, we probe sensitivity of the production cross section to these anomalous couplings. The cross section is mainly sensitive to ttH and HZZ couplings. The scalar Yukawa coupling, κ t , and the scaling coupling, g HZZ , lead to largest modification. The cross section could double with allowed parameter range of these parameters. The effects of other anomalous couplings are rather modest. Like triple Higgs boson production, here also the gg → HHZ is not sensitive toκ t in the allowed parameter range mentioned in section 3. We find that tri-linear Higgs boson anomalous self-couplings do not play any significant role in this NNLO process. The derivative HZZ couplings and HHZZ coupling also don't play any significant role in the allowed range mentioned in section 3. In Figs. 12, 13, 14, the effects of various anomalous couplings relevant to gg → HHZ on leading p T (H) distribution are shown. For that we take some benchmark values. We find that both κ t andκ t lead to harder tail compared to the SM prediction. The contribution to the cross section at higher p T is significantly large for higher κ t . This may be a better avenue to see the effect of such interactions. Similar features are found for parameters g values. There seems to be harder tail of this p T distribution for larger values of g (1) 3H , specially for the positive values. The distribution for g (0) 3H doesn't seem to show any special feature. The sensitivity of the NNLO process to anomalous interactions is similar at 100 TeV.
In this paper our focus has been the NNLO gg → HHZ process. That is why we have presented detailed results for this process. One may ask how sensitive is LO process to anomalous interactions. We have explored this sensitivity by using Madgraph. By including anomalous vertices, we find that LO quark-antiquark annihilation process is quite sensitive to anomalous derivative HZZ interaction. The cross section can increase by an order of magnitude. The increase is more at higher center-of-mass energy. This is unlike NNLO gg → HHZ process.
The process pp → HHZ is likely to be observable at the LHC. Including the NLO and NNLO corrections, the cross section is about 320 ab at 13 TeV. With the full integrated luminosity, one may expect around 1000 events. This process should be visible using various multilepton + jets signature. The main irreducible background of ZZZ has the cross section of 9.2 fb. But Z → bb, τ τ branching ratios are smaller by a factor of 2-3 as compared to the Higgs boson decay. Then by restricting the number of jets in the signature, one may be able to detect this process. To look for evidence beyond the standard model, one may look for enhancement in large p T (Z) events. Reducible backgrounds, as mentioned above, may be tamed by flavor-tagging of jets. At higher energy machines, this process would definitely be observable.

Conclusions
In this paper, we have considered the processes -pp → HHH, HHγ, and HHZ. Our focus was on the gluon-gluon fusion contribution to them. The one-loop amplitude for the process gg → HHγ vanishes exactly. The process pp → HHH is important as it involves both trilinear and quartic Higgs boson couplings. A measurement of this process along with di-Higgs production can help in determining the form of the Higgs potential. It may be seen only if there exists anomalous interactions. This process is specially sensitive to trilinear Higgs boson couplings. This process can be observed at large centerof-mass energy machines with high luminosity. It will be challenging. The process pp → HHZ may be observable at the LHC after accumulation of 3 ab −1 luminosity. The GG(LO) contribution to this process is actually a NNLO contribution in α s , and due to a large gluon flux it is 14% of the QQ(NLO) contribution to pp → HHZ at 13 TeV LHC. In certain kinematic windows GG(LO) contribution can be more than 20%. At a 100 TeV machine, gg → HHZ can be as important as qq → HHZ. This process is important, as it involves HHH and HHZZ couplings and is background to triple Higgs production. The effect of ttH and HZZ anomalous couplings are more significant in the distributions than in the total cross section. This process can definitely be observed at higher energy, such as 100 TeV, machines with enough luminosity.