Di-vector boson production in association with a Higgs boson at hadron colliders

We consider the production of a Higgs boson in association with two electroweak vector bosons at hadron colliders. In particular, we examine $\gamma\gamma H$, $\gamma ZH$, $ZZH$, and $W^{+}W^{-}H$ production at the LHC (14 TeV), HE-LHC (27 TeV), and FCC-hh (100 TeV) colliders. Our main focus is to estimate the gluon-gluon ($gg$) channel ($gg \to VV^\prime H$) contributions to $pp \to VV^\prime H~(V,V^\prime=\gamma,Z,W)$ and compare them with corresponding contributions arising from the quark-quark ($qq$) channel. Technically, the leading order $gg$ channel contribution to $pp \to VV^\prime H$ cross section is an NNLO correction in $\alpha_s$. In the processes under consideration, we find that in the $gg$ channel, $W^{+}W^{-}H$ has the largest cross section. However, relative contribution of the $gg$ channel is more important for the $pp \to ZZH$ production. At the FCC-hh, $gg \to ZZH$ contribution is comparable with the NLO QCD correction to $qq \to ZZH$. We also compute the cross sections when $W$ and $Z$-bosons are polarized. In the production of $W^{+}W^{-}H$ and $ZZH$, we find that the $gg$ channel contributes more significantly when the vector bosons are longitudinally polarized. By examining such events, one can increase the fraction of the $gg$ channel contribution to these processes. Further, we have studied beyond-the-standard-model effects in the $\kappa$-framework. We find that the $gg$ channel processes $ZZH$ and $WWH$ have very mild dependence on $\kappa_\lambda$, but strong dependence on $\kappa_t$ and $\kappa_V$. The $qq$ channel processes mainly depend on $\kappa_V$. Dependence of the $gg$ channel contribution on $\kappa_V$ is stronger than that of the $qq$ channel contribution. Therefore focusing on events with longitudinally polarized $W$ and $Z$-bosons, one can find stronger dependence on $\kappa_V$ that can help us measure this parameter.


Introduction
After the discovery of a Higgs-like resonance, with a mass of 125 GeV, at the Large Hadron Collider (LHC) in 2012, various properties of this new particle have been studied. The spin and parity measurements have established it as a 0 + state at 99.9% CL against alternative scenarios [1]. Couplings of this new particle with the fermions and gauge bosons predicted in the standard model are getting constrained as more and more data are being analyzed by the LHC experiments [2][3][4]. To this end, the vector-boson fusion production of the Higgs boson, associated production of V H(V = Z, W ), and Higgs boson's decay into vector bosons set limits on the HV V couplings [5,6]. The gluon-gluon (gg) channel production of the Higgs boson helps in constraining the Htt coupling [6]. In addition, the evidence for the associated production of Higgs boson with a top-quark pair [7,8] will provide the direct measurement of Htt coupling. We still need to measure the trilinear and quartic Higgs self-couplings in order to know the form of the Higgs potential which will in turn reveal the exact symmetry breaking mechanism. The Higgs self-couplings can be probed directly in multi-Higgs production processes [9][10][11]. Recently, indirect methods of probing them at hadron and lepton colliders have also been proposed [12][13][14]. Similarly, the quartic couplings involving Higgs and vector bosons HHV V are also not constrained independently. This coupling can be probed in the vector-boson fusion production of a Higgs boson pair [15,16]. In order to find the signals of new physics, it is important that we improve our theoretical predictions for the processes involving Higgs boson at current and future colliders.
Loop-induced decay and scattering processes can play an important role in searching for new physics.
In the presence of new physics (new particles and/or interactions), the rates for such processes can differ significantly from their standard model predictions. In this regard, many gg channel scattering processes in 2 → 2 and 2 → 3 category have been studied [11,. In the present work, we are interested in loop-induced gg channel contribution to V V H (γγH, γZH, ZZH, and W + W − H) production. In QCD perturbation theory, the leading order gg channel contribution to pp → V V H is an NNLO contribution at the cross section level. Due to many electroweak couplings involved and loop-induced nature of gg → V V H processes, their cross sections are expected to be small. However, they can be important at high energy hadron colliders like 100 TeV pp collider such as proposed hadronic Future Circular Collider (FCC-hh) facility at CERN [42] and Super Proton-Proton Collider (SPPC) facility in China [43]. At such energy scale, the gluon flux inside the proton becomes very large. In fact, for γγH, the gg channel gives the dominant contribution.
Unlike the quark-quark contributions, which are mainly sensitive to HV V couplings, the gluongluon contribution allows access to Htt, HHH, and HHV V couplings as well. Note that the processes under consideration are background to pp → HH when one of the Higgs bosons decays into γγ/γZ/ZZ * or W W * final states. The process pp → ZZH is also a background to pp → HHH when two of the three Higgs bosons decay into bb final states. In this work, we present a detailed study of gg → γγH and γZH for the first time in the SM. The gg channel contribution to ZZH and W W H in the SM have been studied in the past [26,44,45]. We have presented the ZZH and W W H calculations in detail and have proposed methods to enhance the relative contribution of gluon-gluon channel over quark-quark channel. Since loop-induced processes are sensitive to new physics, we also study the effect of new physics in all V V H processes using a common BSM framework -the κ-framework. Going beyond the κ-framework, we have treated the HHV V coupling independently and emphasized its effect in ZZH and W W H processes. BSM study in a more sophisticated framework is desirable but it is beyond the scope of the present work.
Experimentally, W and Z-boson polarizations have been measured at hadronic colliders [46][47][48]. We also compute the cross sections for the processes when these bosons are polarized. For each process, the different production channels contribute predominantly to specific polarization configurations. This can help in enhancing the contribution of the gg channel, as compared to the qq channel. The gg channel have sometimes stronger dependence on the kappa parameters, in particular on κ V . Therefore, an event sample with larger gg channel contribution can be helpful.
The paper is organized as follows. In Sec. 2, we discuss the Feynman diagrams which contribute to gg → V V H amplitudes. The model independent framework to study new physics effects is outlined in Sec. 3. In Sec. 4, we provide details on the calculation techniques and various checks that we have performed in order to ensure the correctness of our calculation. In Sec. 5, we present numerical results in SM and BSM scenarios for all the V V H processes. Finally, we summarize our results and conclude in Sec. 6.

Gluon fusion Contribution to V V H
The gg channel contribution to pp → V V H is due to a loop-induced scattering process mediated by a quark-loop. The classes of diagrams contributing to gg → V V H processes are shown in Fig. 1 1 . For convenience, the diagrams contributing to gg → W W H process are shown separately in Fig. 2. The gg → γγH process receives contribution only from the pentagon diagrams, while, γZH receives contribution from both pentagon and box class of diagrams. In case of gg → ZZH, W W H processes, triangle class of diagrams also contribute. We have taken all quarks but the top-quark as massless. Therefore, the top-quark contribution is relevant in diagrams where Higgs boson is directly attached to the quark loop. In the diagrams where Higgs boson does not directly couple to the quark loop, light quarks can also contribute. The complete set of diagrams for each process can be obtained by permuting external legs. These permutations imply that there are 24 diagrams in pentagon topology, 6 diagrams in each box topology and 2 diagrams in each triangle topology. The diagrams in which only one type of quark flavor runs in the loop, are not independent. Due to Furry's theorem only half of them are independent [50]. This observation leads to a significant simplification in the overall calculation. This simplification, however, is not applicable to the W W H case, where flavor changing interaction is involved in the quark loop. For example, see (a) and (b) in Fig. 2.
Thus, there are 12 independent pentagon diagrams ( Fig. 1(a)) due to top-quark loop contributing to gg → γγH process. Similarly, the gg → γZH process receives contribution from 12 independent pentagon diagrams ( Fig. 1(a)) due to top-quark loop and 3 independent box diagrams ( Fig. 1(b)) for each quark flavor. In principle, 5 light quarks (u, d, c, s, b) and 1 heavy quark (t) contribute. The box class of diagrams arise due to ZZH coupling and has effective box topology of gg → γZ * amplitude. Furry's theorem, in this case, implies that the axial vector coupling of Z boson with quark does not contribute to gg → γZH amplitude.
Like the gg → γZH process, the gg → ZZH amplitude receives contribution from 12 independent pentagon diagrams with top-quark in the loop ( Fig. 1(a)). However, there are 6 independent box diagrams with effective box topology of gg → ZZ * amplitude for each quark flavor which covers the possibilities of H coupling with any of the two external Z bosons ( Fig. 1(b)). Further, a new box type contribution arises which has effective box topology of gg → HH * amplitude ( Fig. 1(c)). Once again there are 3 such independent diagrams with only top-quark in the loop. In addition to that, there are 4 independent triangle diagrams with top-quark in the loop and having effective triangle topology of gg → H * amplitude ( Fig. 1 (d), (e), (f)). In gg → ZZH amplitude, the Furry's theorem implies that the vector and axial vector coupling of Z boson with quarks can contribute at quadratic level only. Process gg → γγH receives contribution only from (a) type diagrams, while gg → γZH gets contribution from both (a) and (b) type diagrams. In the case of ZZH, all the diagrams contribute; the diagrams (b) and (f) cover the situation in which H is attached to a Z boson. Among all V V H amplitudes, the structure of gg → W W H amplitude is the most complex. Due to the involvement of flavor changing interactions in Fig. 2 (a) and (b), the Furry's theorem is not applicable to these diagrams. Therefore, 24 independent pentagon diagrams contribute to gg → W W H process for each generation of quarks. However, since we neglect Higgs coupling with light quarks including the b quark, there are only 12 non-zero independent pentagon diagrams. In Fig. 2 (b), all the three quark generations contribute. Taking into account the possibility of Higgs boson coupling with any of the two external W bosons, there are total 12 independent box diagrams of type (b) for each generation. In diagrams (a) and (b), the axial vector coupling of W with quarks contributes at quadratic as well as at linear level. Like in the gg → ZZH process, there are 3 independent box diagrams of type (c). Due to ZW W coupling, a new box contribution of type (d) having effective box topology of gg → HZ * amplitude appears. Furry's theorem for diagram (d) implies that the vector coupling of Z with quarks does not contribute to the amplitude. The same explains the absence of similar box diagram due to γW W coupling. Further, there are 4 independent triangle diagrams with top-quark loop (Fig. 2 (e), (f) (g)) as in case of the gg → ZZH process. A new type of 3 independent triangle diagrams for each quark flavor with effective triangle topology of gg → Z * amplitude appears, once again due to ZW W coupling ( Fig. 2 (h), (i)). These triangle diagrams are anomalous and they can receive contribution only from the third generation quarks as the bottom and top-quarks have very different masses. This is indeed the case for (h) type diagrams. However, we find that (i) type diagrams do not contribute. This is explained in the appendix A.

BSM Parametrization
Measuring the couplings of the Higgs boson with fermions, gauge bosons and with itself is an important aspect of finding the signatures of new physics at colliders. With the help of the data collected so far at the LHC, we now know couplings of the Higgs boson with the top quark with an accuracy of 10-20% and with vector bosons with an accuracy of 10% at 1σ [51]. The Higgs self-couplings, on the other hand, are practically unconstrained [52].
To study the new physics effects in V V H processes, we take the simplest approach of modifying the SM like couplings only, also known as the kappa framework for the parametrization of new physics [53,54]. In this framework, no new Lorentz structures and no new interaction vertices appear. The LHC experiments have interpreted the data using this framework so far. The couplings of our interest are Htt, HV V , HHH and HHV V . Out of these couplings, gg → γγH is sensitive to only Htt coupling. The HV V coupling affects all other processes. The couplings HHH and HHV V affect only gg → V V H, V = Z, W processes.
The modification in these couplings due to new physics is implemented through scale factor κ i for various couplings of the Higgs boson in the SM. In kappa framework, there are three such scale factors namely κ t for Higgs coupling with top-quark, κ V for Higgs coupling with vector bosons (κ HZZ = κ HW W = κ V ) 2 and κ λ for Higgs coupling with itself. Since in the SM both HV V and HHV V couplings are related, the scaling of HHV V coupling is also parametrized by κ V . In a more generic BSM framework, the HHV V coupling, in principle, can be independent of HV V coupling.
In the presence of BSM effects, the amplitudes for the gg channel processes depend on κ t , κ V , and κ λ as follows.
In the above, the amplitude M SM i is related to one of the diagram classes displayed in Fig. 1  does not affect the gauge invariance of the amplitudes with respect to the gluons as it will become clear in the next section. The standard model prediction can be obtained by setting κ t = κ V = κ λ = 1. Thus, except in gg → γγH, we can expect nontrivial interference effects on total and differential cross sections for gg → V V H processes due to new physics in κ-framework.

Calculation and Checks
The calculation of quark-loop diagrams is carried out using a semi-automated in-house package OVReduce [55] which allows the calculation of any one-loop amplitude with maximum five propagators in the loop. The main steps involved in our calculation are: quark-loop trace evaluation, one-loop tensor reduction to master integrals and evaluation of master integrals. Trace calculation and simplification of the amplitude are done using symbolic manipulation software FORM [56]. Tensor reduction of one-loop amplitudes into one-loop master integrals is done numerically following the method of Oldenborgh-Vermaseren [57]. Further, the one-loop master integrals are also calculated numerically using the OneLOop package [58]. More details on this can be found in [23]. We perform the calculation in 4 − 2 space-time dimensions to regulate ultraviolet (UV) and infrared (IR) singularities of one-loop master integrals. Since the couplings of Z and W bosons with quarks involve γ 5 , the trace calculation needs special care. We have used 4-dimensional properties of γ 5 in the calculation. This works because the SM is anomaly free. We have chosen Unitary gauge for the calculation of the amplitudes.
The amplitude calculation for each process can be efficiently organized using prototype amplitudes for each class of diagrams. For example, amplitudes for all the 12 independent pentagon diagrams in gg → γγH process can be obtained using only one prototype pentagon amplitude. Similarly, prototype amplitudes can be identified for each topology contributing to each process. The full amplitude for each process is a function of external momenta and polarization vectors/helicities. Due to huge expressions of the amplitudes, we calculate helicity amplitudes and the squaring of the amplitude for each process is done numerically. The number of helicity amplitudes for gg → γγH, γZH, ZZH, W W H processes are 16, 24, 36, and 36, respectively.
There are a number of checks that we have performed in order to ensure the correctness of the amplitudes. We have checked that the amplitudes are separately UV and IR finite. In 4 − 2 dimensions, these divergences appear as poles in 1/ (for UV and IR) and 1/ 2 (for IR only). Each pentagon diagram is UV finite. This we expect from the naive power counting. The individual box diagram is not UV finite, however, the full box amplitude, in each class, is UV finite. The UV finiteness of triangle amplitudes holds for each diagram. One-loop diagrams with all massive internal lines are IR finite, as expected. Thus, IR finiteness check is relevant to the diagrams with massless quarks in the loop. This includes box class of diagrams of Fig. 1(b) in gg → γZH and ZZH. In gg → W W H case, potentially IR divergent diagrams include Fig. 2(a), (b), (h) and (i). Unlike UV, the IR finiteness holds for each diagram [23].
We have also checked the gauge invariance of the amplitudes with respect to the external gluons. For that we numerically replace the gluon polarization vector µ (p) by its four momentum p µ and expect a gauge invariant amplitude to vanish. We find that the gauge invariance check holds for each class of diagrams. This is expected because different box and triangle topologies for each process arise due to the existence of various electroweak couplings. This is a very strong check on our calculation which is organized using only prototype amplitudes. However, this check cannot verify relative signs between different classes of diagrams. In order to verify such relative signs, one needs to perform gauge invariance check in electroweak theory which is a non-trivial task 3 . We rather rely on cross-checking the calculation using different methods and tools. We have compared our matrix element for each process with those calculated using MadLoop [60] and have found an excellent agreement. Being process specific, our code is efficient and provides greater flexibility when producing phenomenological result. Numerical predictions for cross section and kinematic distributions are obtained using Monte Carlo techniques for phase space integration. We use AMCI [61] package for Monte Carlo phase space integration which is based on VEGAS [62] algorithm and allows parallelization of phase space point generation and matrix-element computation using PVM software [63].

Numerical Results
The cross section and kinematic distributions for pp → V V H processes in SM and in BSM constitute the main results of this section. The numerical results are produced using following basic selection cuts unless stated otherwise, The results for the gg channel processes are calculated using CT14LO [64] parton distribution function (PDF) and partonic center-of-mass energy ( √ŝ ) is chosen as common scale for renormalization (µ R ) and factorization (µ F ). The results are obtained for three different choices of collider energies: √ s = 14, 27, and 100 TeV. From phenomenological point of view we will focus on p T (H) and We compare the gg channel contribution to pp → V V H with contribution arising from the qq channels. The qq channel contribution at LO and NLO (QCD) is calculated using MadGraph5 aMC@NLO [60] in five flavor scheme for all but WWH production. The qq channel contribution to WWH production is instead calculated in four flavor scheme 4 . The LO qq channel contributions are pure electroweak processes and they do not depend on α s . For LO and NLO (QCD) results, we use CTEQ14LO and CT14NLO PDFs, respectively [64]. The scale choice is same as in the gg channel calculation. In both gg and qq channel calculations, the scale uncertainties are estimated by varying µ R and µ F independently by a factor of two. We quote only minimum and maximum uncertainties thus obtained.
To quantify the relative importance of the gg channel contribution in processes dominated by the qq channel, we define following ratio, This ratio compares the leading order gg channel contribution with NLO QCD correction in the qq channel. Recall that technically gg channels contribute at NNLO. Similarly, at differential level we define another ratio, where, X denotes a kinematic variable.
As mentioned in section 3, the BSM effects are parametrized in terms of scale factors κ t , κ V and 4 For WWH production, currently MadGraph5 aMC@NLO cannot produce NLO correction to the bb channel. κ λ . In order to compare their relative importance, we vary them independently by 10% about their SM values. Further, we comment on the effect of κ λ and κ HHV V (the scale factor for the HHV V coupling 5 ) which are least constrained at present, in ZZH and W W H processes. The cross section for this process is dominated by the gg channel. In the qq channel, only bottomquark initiated subprocess contribute to γγH production. However, this cross section is quite small, owing to small bottom Yukawa coupling. In Tab. 1, we compare the gg and qq channel contributions to the hadronic cross section at 14, 27 and 100 TeV colliders. The results are with minimum 50 GeV transverse momentum of photons. We find that the gg channel contribution increases 40 times as the collider center-of-mass energy goes from 14 TeV to 100 TeV. Due to a small cross-section, this process cannot be observed at the HL-LHC; FCC-pp will be more suitable. The gg channel contribution becomes important at higher center-of-mass energy collider, as in this case smaller partonic momentum fractions (x) are accessible, where gluon flux is significantly large. The scale uncertainties on the cross sections for the gg channel are in the range of 20-30%. It is clear from the table that the qq channel contribution is negligible compared to the gg channel contribution. It is merely 1% of the gg channel contribution even after including the NLO-QCD corrections.

Predictions for the
In Fig. 3, we have plotted p T distributions for hardest photon, next-to-hardest photon, and Higgs in the left figure, and diphoton invariant mass distribution (M (γγ)) in the right figure for the 100 TeV collider (FCC-hh). The p T distributions for them peak around 150 GeV, 90 GeV, and 70 GeV, respectively. We find that the tail of p T (H) is softer than that of photons. The M (γγ) distribution shows an interesting feature -it has two peaks. The right peak occurs at around 350 GeV, exhibiting the tt threshold effect in the distribution. To verify that the second peak is indeed due to tt threshold effect, we changed in our code the value of m t from 173 GeV to 200 GeV and the second peak was found to get shifted to 400 GeV.
As mentioned before, this process is a background to double Higgs production process when one of the Higgs bosons decays into a photon pair. To manage the background one usually looks at 'γγbb' final state, instead of 'bbbb', as the signature of the double Higgs boson production. At a 100 TeV collider, while the cross section for the gg → γγH production, with the cuts in Eq. 5, is about 220 ab, the cross section for gg → HH → γγH, with the same set of cuts, is about 2600 ab. From the right GeV] s = 100 TeV ×200 Figure 3: Kinematic distributions for gg → γγH process in the SM at 100 TeV. In the p T distribution plot, γ 1 and γ 2 refer to the hardest and second hardest photons in p T , respectively. In the right plot, we show M (γγ) distribution for gg → γγH. In addition, the total cross section for the gg → HH → γγH process has been shown at 125 GeV. "×200" implies that the height of purple vertical line should be multiplied by a factor of 200 in order to get the correct cross section for the gg → HH → γγH process.
panel of Fig. 3, it can be seen that the cross section for γγH production in the bin from 120 GeV to 140 GeV is about 3 ab. On the other hand, all the cross section for gg → HH → γγH is concentrated in a very narrow width around the mass of Higgs, 125 GeV 6 . As a result, gg → γγH is an insignificant background to the process gg → HH → γγH.
Regarding anomalous coupling contributions, we note that as only pentagon diagrams contribute to the process gg → γγH, its cross section scales as κ 2 t . So a 10% change in κ t will change the cross section and distributions by about 20%. For the qq channel process, the cross section is too small. It depends on κ b , which we do not change from the standard model value.

Predictions for the pp → γZH process
Unlike γγH case, the γZH production receives dominant contribution from the qq channel. With p γ T > 50 GeV, the gg channel contributions to γZH production at 14, 27, and 100 TeV colliders are 4 ab, 16 ab, and 168 ab, respectively. The corresponding values for the LO qq channel contribution are 689 ab, 1733 ab, and 7498 ab, respectively. From Tab. 2, it can be seen that R 1 , which is the ratio of the gg channel contribution to NLO correction in the qq channel, is as small as 0.06 for 100 TeV collider, and even smaller for HE-LHC (27 TeV) and LHC (14 TeV). The scale uncertainties for the gg channel are around 20% while those for the qq channel at NLO are in the range of 2 − 3%. A larger scale dependence in the gg channel contribution can be attributed to the presence of higher power of α s factor in the gg amplitudes.  In Tab. 3, the effect of various p γ T cuts in gg and qq channels has been shown. As the cut on p γ T increases, the qq channel cross section decreases faster than the gg channel. In going from 50 GeV to 200 GeV cut, the cross section of the gg channel decreases roughly by a factor of 6, while that of the qq channel decreases by a factor of 9. Thus relative contribution from the gg channel can be enhanced with the help of harder p γ T cut. We find that the p T (H) cuts have opposite effect i.e. the gg channel is favored at low p T (H). 168  7498  10430  100  95  2812  4072  150  47  1366  2069  200  28  765  1190   Table 3: Effect of p γ T cut on the cross section of pp → γZH production at the 100 TeV collider (FCC-hh).
In Fig. 4, we have displayed p T distributions for the final state particles on the left, and γZ pair invariant mass distribution on the right for the 100 TeV collider. The p T distributions peak around 100 GeV while the M (γZ) distribution peaks around 200 GeV. Like the case of gg → γγH process as a background to gg → HH → γγH, the gg → γZH process is also an insignificant background to gg → HH → γZH. This is because at a 100 TeV collider, with the cuts in Eq. 5, the cross section for gg → HH → γZH is about 2000 ab, while the cross section for gg → γZH process is about 170 ab. Moreover, all the cross section for the gg → HH → γZH process congregates around the mass of the decaying Higgs boson, 125 GeV 7 , while, as can be seen from the right panel of the Fig. 4, the cross section for the gg → γZH process in the bin from 120 GeV to 140 GeV is about 3 ab. However, the qq channel for γZH production may act as an important background for the gg → HH → γZH process.
In Fig. 5, we show p T (H) distributions for different classes of diagrams -pentagon, box, and sum of their individual contributions, their interference, and total at the 100 TeV collider. The contribution of the box diagrams is more than the pentagon diagrams mainly because of the light quark contributions.
The interference effect between the pentagon and box diagrams has kinematic dependence. We find that in the region of our kinematic interest, it is always destructive and, near the peak, its effect is close to -30%. GeV] s = 100 TeV without top in Fig. 1(b) with top in Fig. 1 Fig. 1(b) without top in Fig. 1(b) Figure 5: Left: The contribution of pentagon (blue) and box (green) diagrams, as well as their squared sum, interference, and total contribution to p T (H) distributions for the gg → γZH process at the 100 TeV FCC-hh collider. Right: The effect of excluding top-quark contribution from the diagrams in Fig. 1(b) to the full amplitude.
Since the ggγZ * type box amplitude does not depend on the axial-vector coupling of the off-shell longitudinal Z-boson with the quarks, the top-quark contribution is not very significant at the level of total cross section. This is shown in the right panel of the Fig. 5. We can see that in the tail where top quark is effectively light, the cross section increases by about 20%. We have noted that the relative importance of gluon fusion channel can be enhanced by applying higher p T (γ) cuts. To distinguish the gg channel contribution from the dominant qq channel, one can use the polarized cross sections and distributions. In Fig. 6, we have displayed the LO cross sections for various helicity states of the final state particles, γ and Z boson. The figure also shows the contribution of various polarization states of the initial state particles. We cannot measure the polarization of the initial state particles that are in a bound state, proton. However, experimentally, one can measure the Z-boson polarization [46][47][48]. The polarization of photon has been measured by the LHCb collaboration in b-baryon's decay [65][66][67][68][69][70][71]. At a 100 TeV collider, the contribution of the gg channel process to the production of γZH is only 2.2%. However, if we look at those final states where photon and Z-boson have the same transverse polarization, then this ratio increases to 10 − 11%. (The qq channel makes largest contribution when the Z boson is longitudinally polarized.) This is a non-trivial contribution, and can be measured, if enough integrated luminosity is available. In Fig.7, we have plotted the Higgs boson and Z-boson p T distributions. By making appropriate cuts on the small and large p T of these particles, we can further enhance the gg channel contribution.
Turning to the effect of anomalous couplings, we find that the gg channel shows very small dependence on the κ t , as it is present only in pentagon diagrams whose contribution is small (see Fig. 5). However, it strongly depends on κ V , as the box-diagrams contribution is much more than the pentagondiagram contribution. We find that the change in cross section for κ t = 1.1(0.9) is 5.4% (-1.2%). On the other hand, for κ V = 1.1(0.9) the cross section changes by 18% (15%). We do not show the effect of anomalous couplings on the distribution. It can be understood qualitatively from Eq. 2 and Fig. 5 in the gg channel. The qq channel is sensitive to κ V only. The amplitude has overall linear dependence on κ V due to which the effect of anomalous coupling k V is flat for both total and differential cross sections.

Predictions for pp → ZZH
The cross sections for ZZH production via various channels have been tabulated in Table 4 along with the corresponding scale uncertainties. The total cross section for gg → ZZH is significantly larger than that of gg → γZH. This increase is mainly due to the contribution from axial-vector coupling of Z with quarks. The gg channel contributions to ZZH production at 14, 27, and 100 TeV colliders are 124 ab, 579 ab, and 7408 ab, respectively. The corresponding values of the LO qq channel contributions are 2184, 5997, and 36830 ab, respectively. The ratio, R 1 , is found to be 0.25, 0.4, and 1.05, respectively. Thus at 100 TeV, the gg channel contribution is as important as the QCD NLO correction. As has already been discussed, this increase in ratio R 1 with collider energy is due to the large gluon flux.  In the gg channel, the scale uncertainties of the total cross sections are in the range of 20-30% which is similar to the scale uncertainties observed for γγH and γZH. We find that the uncertainty due to the renormalization scale variation is more than that due to the factorization scale variation. While the change in the renormalization scale mainly changes α s , the change in the factorization scale changes the parton distribution function. The uncertainty for the renormalization scale variation is nearly same at all the collider energies. This happens as the contribution to the total cross section comes from nearly same region of partonic center of mass energy of the process and in every bin of this region, α s changes by nearly same factor for the change in the renormalization scale. However, uncertainty for the factorization scale variation is different for different colliders. This happens as for different collider energies, different x regions contribute to the process and for different x regions change in parton distribution function with the factorization scale is different, where x is partonic momentum fraction. We have also observed that with an increase in the factorization scale, for 14 and 27 TeV colliders, the cross-section decreases; however for 100 TeV collider the cross-section increases.  In the tree level qq channel, there is no QCD vertex. So here change in the renormalization scale does not affect the cross section. But, the change in the factorization scale can affect the cross section, and uncertainty increases with collider energy. However, when NLO QCD correction is considered, change in either of renormalization and factorization scales changes the cross section. The uncertainty in the cross section due to the renormalization scale variation is small as NLO QCD correction is much smaller than the LO results. The overall uncertainty in this case is smaller than the LO case, which is expected for higher order calculation. Ratio with top in Fig. 1(b) without top in Fig. 1(b) Figure 9: Left : SM contribution of pentagon (blue), box (green), triangle (gray) diagrams, as well as their squared sum (black), interference (orange) and total (red) contribution to p T (H) distributions in gg → ZZH at 100 TeV collider (FCC-hh). Right: The effect of excluding top-quark contribution from Fig. 1(b) to full amplitude. In Fig. 8, we have plotted p T distributions for leading p T (Z 1 ), next-to-leading p T (Z 2 ), and Higgs boson in the left figure, and Z-pair invariant mass distribution in the right figure for the 100 TeV collider. The p T distributions peak around 100 GeV, 60 GeV, and 80 GeV, respectively. The M (ZZ) distribution peaks around Z pair threshold.
Interference of various diagrams plays a major role in gg → ZZH production. In Fig. 9, we have shown the p T (H) distributions for penta, box, triangle, sum of their individual contributions, interference, and total at the 100 TeV collider (FCC-hh). As can be seen, the box diagrams give the largest contribution, then comes the triangle contribution and the penta contributes the least. Like in γZH case, the large box contribution is due to the light quarks in the loop. Further, because of large destructive interference, the total contribution is smaller by about a factor of five than the box contribution. We have found that the top-quark contribution in ggZZ * -type box diagram is quite significant despite the propagators suppression. This is due to the coupling of off-shell longitudinal Z boson (ef-fectively the Goldstone boson) with top-quark and it is proportional to m t 8 . We show the effect of excluding the top-quark contribution in ggZZ * -type box diagram ( Fig.1(b)) on p T (H) distribution in the right panel of Fig. 9. As we expect, excluding top-quark contribution in ggZZ * -type box diagram leads to non-unitary behavior in the full amplitude.
In the left figure of Fig. 10, we see that the shape of p T distribution for Higgs boson in the gg and qq channel processes is nearly same at 100 TeV collider (FCC-hh). The relative importance of the gg channel over the qq channel is visible in the tail. In the right plot, we give p T (H) distribution combining gg and qq (NLO) contributions as the best prediction from our calculations. In the bottom panel of the plot, R 2 signifies the ratio of differential cross section from the gg channel to that from NLO qq channel process. The dashed line shows the ratio of corresponding total cross sections, which is 0.17. At the tail of the distribution, we see the gg channel contribution becomes further important, but there differential cross section itself is quite small. Once again we find that if we categorize events based on the helicity states of the two Z bosons, the relative importance of the gg channel contribution over the qq channel contribution can be increased. From Fig. 11, we see that in the gg channel the longitudinal Z bosons contribute the most, while in the qq channel their transverse helicity states give dominant contribution. The relative cross section of the gg channel with respect to the qq channel is about 20%. However, if we restrict ourselves to the case when both Z-bosons are longitudinally polarized, then this ratio almost doubles. Since the cross section for these polarized states for the gg channel is about 2000 ab, there will be enough events to observe this process at a 100 TeV machine. At the distribution level, from the Fig. 12, we observe that if we restrict ourselves to the contributions from the longitudinal Z bosons with p T (H) beyond 150 GeV, the relative contribution of the gg channel increases significantly. Experimentally, one may look at the signature l + l − l + l − bb. This signature is obtained when Z → l + l − (l = e/µ) and for H →bb. Taking into account the branching ratios, and b-tagging efficiency, one may expect about 75 events at the FCC-hh collider (with 30 ab −1 integrated luminosity) from gg channel and about 210 events from qq channel. This is when both Z bosons are longitudinally polarized. This number will go down when detection and kinematic-cut efficiency factors are included. However if in future, one could use hadronic decay modes of a Z boson to measure its polarization, then the number of events would increase.
As can be seen from Eq. 3, the gg channel depends on κ t , κ V , and κ λ . We vary these κ's by 10% from their SM values. The gg channel strongly depends on both κ t and κ V . In the gg channel, ±10% change in κ t causes 68% and -18% variations in the cross section, respectively. And ±10% change in κ V causes 45% and -28% changes in the cross section, respectively. Similar variation in κ λ does not lead to much variation in the total cross section. Since this coupling is not yet well constrained, we will discuss it in detail in subsection 5.5. In Fig. 13, we display the effect of κ t and κ V on p T (H) distribution. We show the absolute distribution in the top panel, while in the bottom panel we show the ratio of distribution with anomalous coupling to that with the SM coupling. We can see that in the presence of anomalous κ t and κ V , the shape of the distribution remains more or less same. However, due to non-trivial interference effects, the modifications in presence of anomalous couplings are not same in all the bins. We see that for κ t = 1.1 the cross section in the bins near tail of the distribution increases by a factor of 2. On the other hand for κ V = 1.1, the maximum change in the cross section is around 1.5. Thus tail of the distributions are more sensitive to modifications in couplings due to high scale new physics. The qq channel depends mainly on κ V . However, as we have considered bottom quark contribution also, the qq channel depends on κ λ as well. In the qq channel, κ V comes as an overall factor both for LO and NLO amplitude, and so the effect of 10% change in κ V causes around 20% change in the cross section, both at total and differential levels. We find a very mild dependence on κ λ .

Predictions for pp → W W H
The cross section for this process is the largest among all the V V H processes considered in this paper. In Tab. 5, we report the cross section predictions for W W H process at different collider center-of-mass energies. The gg channel contributions to W W H production at 14, 27, and 100 TeV colliders are 290 ab, 1344 ab, and 17403 ab, respectively. These numbers are roughly 2.3 times higher than ZZH cross sections. As regards scale uncertainties, the gg → W W H cross sections follow the same pattern as observed in gg → ZZH. The corresponding values of the LO qq channel cross sections are 8658, 23040, and 128000 ab, respectively 9 . The ratio, R 1 , is found to be 0.15, 0.19, and 0.43, respectively. Unlike ZZH production, the contribution of the gg channel is relatively smaller.  In the left figure of Fig. 14, we can see that the p T distribution of W + and W − overlap with each other, which is expected in the case of the gg channel. The p T (H) distribution peaks around 100 GeV, and its fall in the tail is slower than that of p T (W ± ) distributions. In the right of Fig. 14, the distribution for invariant mass of W + and W − has been shown, which peaks around 200 GeV. Like gg → ZZH production case, in gg → W W H production also, interference of various diagrams plays a major role. On the left of Fig. 15, we have shown p T (H) distributions for individual topologies as well as for their interference at a 100 TeV collider. The box contribution is the largest in all the bins while the pentagon contribution is the lowest beyond p T > 100 GeV. The total contribution is much smaller than the box contribution because of strong destructive interference effect which is shown by orange line in the figure. GeV] s = 100 TeV without 3rd gen. in Fig. 2(b) with 3rd gen. in Fig. 2 Ratio with 3rd gen. in Fig. 2(b) without 3rd gen. in Fig. 2 Due to the presence of top quark propagators in ggW W * type box diagram, one may naively think of a suppressed contribution from the third generation quarks at low p T (H). In Fig. 15, we show the effect of excluding the third generation quark contribution from the ggW W * type box diagram, on the p T (H) distribution. Like in gg → ZZH, the third generation quark contribution in ggW W * type box diagram is necessary for the unitarization of the full amplitude.
In the left plot of Fig. 16, the normalized p T distributions for Higgs boson in the gg and qq channel processes have been shown for 100 TeV collider (FCC-hh). The p T (H) distribution in the gg channel peaks slightly on the harder side making the channel more relevant in higher p T (H) bins. To quantify it better we also plot the the ratio of distributions due to qq (NLO) + gg (LO) and qq (NLO). At differential level the ratio varies between 1.05 and 1.18 compared to its value (1.1) for the total cross section. Once again, we find that the gg channel contribution is more relevant at higher p T where its contribution reaches 18%. Similar to the case of ZZH, for this process also, the cross section in the gg channel is dominated by longitudinally polarized W -bosons (Fig. 17). The relative contribution of this channel is about 13%, with respect to the qq channel. However, when both W -bosons are longitudinally polarized, then this ratio increases to 32%. There will also be enough events at a 100 TeV collider to observe the gg channel contribution. The relative contribution of the gg channel over the qq channel can be further increased by requiring the p T (W ) to be beyond a certain value between 50 and 100 GeV, see Fig. 18. Here also one may consider leptonic decay channel for W bosons, as that will help in the measurement of its polarization. We consider l + ν l l −ν l bb final state as the signature. Here, as before l = e/µ. In the literature, various techniques, including Neural Network methods have been discussed to measure the W boson momentum [72]. Taking into account the branching ratios and the b-tagging efficiency, one may expect about 1750 events from gg channel and 5900 events from the qq channel at the FCC-hh collider with 30 ab −1 integrated luminosity. The number of these events would change depending on the detector and kinematic-cut efficiency factors. Next, we focus on the effect of anomalous couplings on the total and differential cross sections. The gg channel depends on κ t , κ λ , and κ V (see Eq. 4). We find that the channel is mostly sensitive to κ V and κ t . For κ V = 1.1(0.9) the cross section changes by about 38%(-26%). While, for κ t = 1.1(0.9) the cross section changes by about 54% (-3%). The dependence on κ λ is found to be relatively small. In Fig. 19, we show the effect of κ t and κ V on the p T (H) distribution for the gg channel. We do not show the distribution for anomalous κ λ as its effect on cross section is very small for 10% variation. We see that the shape remains more or less same in presence of anomalous couplings. We see that in the bins around 400 GeV, this ratio is around 1.5 for κ t = 1.1 and κ V = 1.1. For κ t = 0.9, the ratio remains close to 1 throughout all the bins and for κ V = 0.9, it is in the range 0.7-0.8. Similar to the case for qq → ZZH, the qq → W W H cross section is also proportional to κ 2 V at LO and NLO(QCD). So here as well, a 10% change in κ V gives around 20% obvious change in cross section, both at the total and differential levels.

Remarks on anomalous HHH and HHV V couplings
We have seen that the gluon fusion ZZH and W W H processes are most relevant for BSM physics due to their large cross sections. We found that their cross sections do not change much for a 10% variation in κ λ . However, we know that this coupling is presently unconstrained by the experimental data. According to the future projections for HL-LHC, only values κ λ −2 and κ λ 8 can be ruled out [9]. In this range the cross section for ZZH and W W H processes in the gg channel varies significantly. In fact, it can change maximum by a factor of 3. This is shown in the left panel of Fig. 20. Notice that the W W H process is more affected by anomalous HHH coupling than ZZH process. Although in SM model, the HHV V (V = Z, W ) coupling is correlated to the HV V coupling, in presence of new physics this correlation may not exist. Keeping this possibility in mind, we have varied the HHV V coupling independently 10 and we find that the cross section changes very strongly. This is shown in the right panel of Fig. 20. We can see that the effect of the HHV V coupling is relatively   larger on gg → ZZH than on gg → W W H. Close to SM values, the difference is negligible. According to a recent search for Higgs boson pair production via vector boson fusion carried out by the ATLAS collaboration using 126 f b −1 data collected at 13 TeV LHC, the allowed values of κ HHV V lie in the range (-0.56, 2.89) at 95% confidence level [16].
The quantity plotted in Fig. 20 is known as signal strength (µ) which has been utilized by experimentalists as observable for data analyses. The signal strength for each process can be parametrized as where κ i = κ λ , κ HHV V . In table 6, we have provided the values of c i 1 and c i 2 for ZZH and W W H processes for the 14 TeV LHC and a 100 TeV pp collider. We note that c k λ 2 is smaller by an order of magnitude than c k λ 1 , suggesting a strong interference effect mentioned before. Therefore, c κ λ 2 is relevant mostly for large values of κ i . On the other hand, c κ HHV V 2 is of the same order as c κ HHV V 1 . Since c i 1 is negative, the cross section increase observed in the figures for κ i < 1 is quite significant, which causes the (negative) lower bound on κ to be tighter than the (positive) upper one. At a 100 TeV pp collider, while the other c i s remain more or less same as that in 14 TeV collider, c κ HHV V 2 increases by around a factor of two, implying the possibility of a far more stringent bound on the HHV V couplings.
Since the gg fusion channel contribution to ZZH and W W H processes cannot be fully separated from the corresponding contributions from the qq channel, the above result should be interpreted care-fully. A realistic estimate of the BSM effects discussed above must include the contributions from qq channel. Since qq channel contributions are insensitive to κ λ and κ HHV V , they can be seen as one of the major backgrounds to the gluon fusion processes. As we have pointed out, the measurement of the polarization of the W/Z boson can help in reducing this background. A systematic signal-background analysis is beyond the scope of the present work. For the benefit of the reader, in Fig 21, we present the ratio σ/σ SM for pp → ZZH, W W H which includes both qq and gg channel contributions as functions of κ λ and κ HHV V . In obtaining these results only standard cuts mentioned in the previous sections have been applied. We can see that at the 14 TeV, the ratio of BSM and SM cross sections due to qq + gg channels is significantly smaller than that due to gg channel alone. Moreover, the ZZH process turns out to be more affected by κ λ and κ HHV V than the W W H. To be more precise, we find that by changing κ λ in the range (−2, 8), the cross section for ZZH process changes in the range 7 − 4% at the 14 TeV. The corresponding change at the 100 TeV falls in the range of 20 − 8%. On the other hand, when changing κ HHV V in the range (−0.56, 2.89), the maximum cross section change in ZZH process is found to be ∼ 8% and ∼ 46% at the 14 TeV and the 100 TeV, respectively. Again, we may mention that the polarization measurements would increase the fraction of gg channel events, thus increasing the dependence on κ HHV V .

Conclusions
In this paper, we have considered production of V V H (γγH, γZH, ZZH, and W W H) at protonproton colliders. We investigated the sensitivity of these processes to various couplings of the Higgs boson, in particular to HHH and HHV V couplings which are practically unconstrained. Our main focus was the gg channel contribution, which occurs at NNLO in α s . The scale uncertainties on the total cross sections are found to be of the order of 20%. A number of checks like UV and IR finiteness and gauge invariance of the amplitudes with respect to the gluons have been performed to ensure the correctness of the calculation. At a 100 TeV collider, the cross sections for these processes via the gg channel range from 0.2 fb to 17 fb, gg → W W H being the dominant channel among all. We have seen the gg → γγH and gg → γZH processes are insignificant background to gg → HH → γγH and gg → HH → γZH, respectively.
We have also compared the gg channel contribution with the fixed order NLO QCD correction to pp → V V H in order to emphasize their relative importance. For γγH production, the gg channel can be said to be the only production channel, as the bb channel process contribution is negligibly small. At a 100 TeV collider, the gg → γZH channel contribution is around 6% of the NLO QCD correction in the qq channel. The γZH production shows one interesting feature: with an increase in the p T cut on photon, the qq channel contribution decreases faster than the gg channel contribution. At this collider, the contribution of the gg channel to ZZH production is as important as the fixed order QCD NLO correction to the qq channel. On the other hand, the gg → W W H channel cross section is around half the fixed order NLO QCD correction to the qq channel. We have observed strong destructive interference effects among various classes of diagrams in gg → γZH, ZZH, W W H. Besides total cross sections at the LHC, HE-LHC, and FCC-hh, we have obtained relevant kinematic distributions at FCC-hh in the gg channel. We find that the p T (H) spectrum from the gg channel is harder than that from the qq channel for ZZH and WWH productions. We have also shown that by selecting events based on the polarization of final state vector bosons, the relative contribution of the gg channel over the qq channel can be enhanced.
In addition to the SM results, effect of anomalous couplings (κ t , κ V , and κ λ ) for Htt, HV V , HHV V , and HHH vertices have been studied in the kappa framework. We find that the new physics effects are quite important in gg → ZZH, W W H processes due to non-trivial interference effects in these processes. A 10% change in κ t on the higher side can enhance the gg → ZZH and W W H cross sections by 68% and 54%, respectively. Similar change in κ V enhances these cross sections by about 40%. Unlike in qq channels, the kinematic distributions in gg channels display non-trivial changes in presence of new physics. The dependence of the gg channel on the κ V is stronger than that of the qq channel. By considering events with longitudinally polarized vector bosons for the processes pp → ZZH, W W H, we can enhance the fraction of the gg channel contribution. This event sample will have even stronger dependence on κ V . Since the HHH and HHV V couplings are not well constrained, we have also considered larger independent variations in κ λ and κ HHV V . We find that the effect of κ HHV V on the cross section is much stronger than that of κ λ . Therefore the process pp → ZZH, W W H with longitudinally polarized Z and W bosons can help in determining the HHV V coupling.
Using momentum conservation p 12 = −p 35 − p 4 and transversality conditions e 3 .p 3 = e 4 .p 4 = 0, we get Using on-shell condition p 2 4 = M 2 W , we arrive at Following similar steps, it can be shown that contraction of current J 2 with p 12 leads to,