$W^+W^-H$ production through bottom quarks fusion at hadron colliders

With the standard model working well in describing the collider data, the focus is now on determining the standard model parameters as well as for any hint of deviation. In particular, the determination of the couplings of the Higgs boson with itself and with other particles of the model is important to better understand the electroweak symmetry breaking sector of the model. In this letter, we look at the process $pp \to WWH$, in particular through the fusion of bottom quarks. Due to the non-negligible coupling of the Higgs boson with the bottom quarks, there is a dependence on the $WWHH$ coupling in this process. This sub-process receives largest contribution when the $W$ bosons are longitudinally polarized. We compute one-loop QCD corrections to various final states with polarized $W$ bosons. We find that the corrections to the final state with the longitudinally polarized $W$ bosons are large. It is shown that the measurement of the polarization of the $W$ bosons can be used as a tool to probe the $WWHH$ coupling in this process. We also examine the effect of varying $WWHH$ coupling in the $\kappa$-framework.


Introduction
Standard Model (SM) has been very successful. It has been tested in a wide variety of low energy and high energy experiments [1,2]. Although there is no firmly established conflict between the data and the standard model predictions, the model is not yet fully validated. In particular, the Higgs sector of the model is not yet fully explored. The Higgs potential can still have many allowed shapes [3]. Self-couplings of the Higgs boson and its couplings with some of the standard model particles are still loosely bound. The more precise measurement of the couplings can also lead to hints to beyond the standard model scenarios.
In this letter, we are interested in the coupling of the Higgs boson with the W and Z bosons (Collectively referred to as V ) in particular, we are interested in the quartic V V HH couplings. In the standard model, the V V H and V V HH couplings are related.
The experimental verification of this relationship is important to put the standard model on a firm footing. There are scenarios beyond the SM, where these couplings are either not related or have different relationship [4]. The ATLAS collaboration has put a bound on this coupling at the Large Hadron Collider (LHC). Using the VBF mechanism of a pair of Higgs boson, and using 126 fb −1 of data at 13 TeV, there is a bound of −0.43 < κ V 2 H 2 < 2.56 at 95% confidence level [5]. Here κ V 2 H 2 is the scaling factor for the V V HH coupling. However, in this process bound on W W HH and ZZHH couplings cannot be separated. The process pp → HHV , where a pair of Higgs bosons are produced in association with a W or a Z boson, allows us to separately measure HHW W and HHZZ couplings. Gluon-gluon fusion would contribute to HHZ production. This mechanism is important at HE-LHC and FCC-hh. However, dependence on the scaling of HHV V coupling is weak. The expected bound from the W HH production at the HL-LHC is −10.6 < κ V 2 H 2 < 11 [6], which is quite loose. Instead of these processes, we consider the process pp → HW W at hadron colliders. This process can help us in measuring HHW W coupling, independent of HHZZ coupling. This processes can take place by both quark-quark and gluon-gluon scattering. At a 100 TeV collider, gluon-gluon scattering and bottom-bottom quark scattering give important contributions. These contributions depend on HHW W coupling. The gluon-gluon contribution is discussed in [7]. This contribution is significantly lower than the contribution of bottom-bottom scattering. But the contribution of bottom-bottom scattering is only about 15 − 20% of the light quarks scattering contribution at the 100 TeV center of mass energy (CME) and at the leading order (LO), light quarks contribution does not depend on W W HH coupling. The dependence on this quartic coupling, W W HH, can be enhanced if we measure the polarization of the final state W bosons. There is significant enhancement of the fraction of the bottom-bottom scattering events, when both W bosons in the final states are longitudinally polarized. The ATLAS and CMS collaborations have measured the W polarization at the LHC [8][9][10]. We compute the one-loop QCD corrections to various combinations of final state W bosons polarization. The longitudinally polarized W boson final states receive largest corrections, leading to even larger fraction of events with bottombottom scattering. We also scale the W W HH coupling and examine the effect of the NLO QCD corrections and the measurement of the polarization of W bosons. It appears that an analysis of W W H events, when both the W bosons are longitudinally polarized, can help in determining the W W HH coupling.
The paper is organized as follows. The second and third sections describes the process and the details of the calculations. In the fourth section, we present the numerical results, and the last section has the conclusions.

The Process
We are interested in quark-quark scattering for the production of W W H. To study W W HH coupling, we consider this process in five-flavour scheme. We study the process bb → W + W − H at hadron colliders. We take bottom quarks as massless but at the same time, we consider bb H Yukawa coupling which is proportional to the mass of the bottom quark. With this consideration, the diagrams with W W HH coupling would appear, with the Higgs boson coupling to the bottom quark. This coupling would not appear at the leading order (LO) for the other quarks in the initial state. This channel has been discussed only with tt H Yukawa couplings [11] but not with bb H Yukawa couplings.
At the LO there are 20 diagrams -9 s-channel and 11 t-channel. A representative set of diagrams are displayed in Fig. 1. Only one of the diagrams has W W HH coupling which is one of our main points of interest. We vary W W HH coupling in order to see its impact on the cross section for the different center of mass energies. There is no strong coupling dependency in the LO diagrams; they solely depend on electroweak couplings. Some of the t-channel diagrams depend on tt H Yukawa couplings and give large contributions to the LO cross section, due to the top-quark mass dependency of tt H Yukawa coupling.
To compute the one-loop QCD corrections to this process, we need to include one-loop diagrams and next-to-leading order (NLO) tree level diagrams. The one-loop diagrams can be categorized as pentagon, box, triangle as well as bubble diagrams. There are 3 pentagon diagrams, 14 box diagrams, 34 triangle diagrams, and 14 bubble diagrams. A few representative NLO diagrams are displayed in Fig. 2. There is only one one-loop diagram (triangle) which has W W HH coupling. Bubble diagrams are UV divergent and a few triangle diagrams are also UV divergent. To remove UV divergence from the amplitude, counterterm (CT) diagrams need to be added to the virtual amplitudes. There are 15 vertex CT diagrams and 14 self energy CT diagrams. A set of CT diagrams are shown in Fig. 3. Also, most of the virtual diagrams are infrared (IR) singular. In order to remove IR singularities from the virtual diagrams, one needs to include real emission diagrams. There are three such processes. These processes are a) bb → W + W − Hg, b) gb → W + W − Hb and c) bg → W + W − Hb. There are 54 Feynman diagrams for the each of these processes. We have shown a few diagrams for the first sub-process in Fig. 3. All these diagrams have been generated using a Mathematica package, FeynArts [12].   Figure 3. Few sample real emission Feynman diagrams for the sub-process bb → W + W − H g and vertex CT diagrams and self energy CT diagrams for the process bb → W + W − H.

Calculations and Checks
We have to perform 2 → 3 and 2 → 4 tree level and 2 → 3 one loop calculations. For the calculation we use helicity methods. As a starting point, we consider a few prototype diagrams in each case. With suitable crossing, and coupling choices, we can compute rest of the diagrams. We compute helicity amplitudes at the matrix element level for the prototype diagrams. These helicity amplitudes can be used to probe the physical observables dependent on the polarization of external particles. As mentioned before, the b-quarks are treated as massless quarks because of their small mass and we use massless spinors for b-quarks. The tree level helicity amplitudes can be written in terms of the spinor products pq and [pq] [13]. For one-loop amplitudes, we use an extra object -the vector current pγ µ q]. We take the functional form of the spinor products pq or [pq] from Ref. [14] and we extend their treatment to calculate the functional form of the vector current pγ µ q]. We have checked that the calculated pγ µ q] satisfies various spinor identities. We adopt four-dimensional-helicity (FDH) scheme [15,16] to compute the amplitudes. In this scheme all spinors, γ-matrices algebra are computed in 4-dimensions. We use package FORM [17] to implement the helicity formalism.
Using FORM, we write helicity amplitude in terms of spinorial objects, scalar products of momenta and polarizations. For the one-loop calculations, we also have tensor and scalar integrals. The one-loop scalar integrals are computed using the package OneLoop [18]. We use an in-house reduction code, OVReduce [19,20], to compute tensor integrals in dimensional regularization. Finally, the phase space integrals have been done with the advanced Monte-Carlo integration (AMCI) package [21]. In AMCI, the VEGAS [22] algorithm is implemented using parallel virtual machine (PVM) package [23].
Few checks have been performed to validate the amplitudes. The one-loop amplitudes have both ultraviolet (UV) and infrared (IR) singularities. UV singularities are removed by using counter-term (CT) diagrams, and the IR divergences are removed using Catani-Seymour (CS) dipole substraction methods. Cancellation of these diverges are powerful checks on the calculation. All UV singularities are removed by fermionic mass and wave function renormalization. There are no UV singularities coming from pentagon and box diagrams as there are no 4-point box tensors in those amplitudes. UV singularities are coming from the triangle as well as bubble diagrams. The appropriate vertex and self energy counterterms (CT) diagrams have been added in total amplitude which gives renormalized amplitude. A few sample CT diagrams are depicted in Fig. 4. We use the MS scheme for massless fermions and the on-shell subtraction scheme for massive fermions.
The next check is infrared (IR) singularity cancellation. We implement the Catani-Seymour dipole subtraction method [24] for the cancellation of IR singularities. Except bubble diagrams, all other virtual diagrams are IR singular. Collectively all IR singularities coming from virtual diagrams cancel with IR singularities coming from real emission diagrams.
Following the Catani-Seymour method, the NLO cross section can be written as Where dσ R , dσ V and dσ A are exclusive cross section, one-loop virtual correction and approximation term respectively. dσ A has the same pointwise singular behaviour as dσ R and hence behaves as a local counterterm for dσ R and then first integration can be performed safely in → 0 limits. The second term of the second integral will give dipole I term which will remove all the infrared singularities from virtual correction and add a finite contribution. The dipole I factor comes from analytical integration of dσ A in d-dimensions over one-parton phase space. It can be written as Where dσ B is born level cross section and the symbol ⊗ describes phase space convolution and sum over spin and color indices. The term dσ B ⊗ I is evaluated over the rest of mparton phase space and cancels all singularities from renormalized virtual amplitudes. As discuss before, we use the FDH scheme, so we take I term in the FDH scheme. The term I given in Ref. [24] is in conventional dimensional regularization (CDR) scheme and in any other regularization scheme (RS) it is given as [25] In the FDH scheme, γ RS Now with this I term, we have checked that the integration in the second term of Eq. 3.1 is IR safe. Also, there are other terms in the dipole subtraction method, called P and K terms which will add finite contributions to σ N LO . These terms come from the factorization of initial-state singularities into parton distribution functions. The color operator algebra, explicit form of V ij,k , I, P and K are given in Ref. [24]. There are three real emission sub-processes that can contribute to σ N LO . These processes are as these processes mimic the Born level process in soft and collinear regions. Due to large contributions, top resonance in the last two processes jeopardizes the perturbative calculation. The cross sections for these two processes are five to six times higher than the Born level cross section. One can't remove those top resonant diagrams as it will affect the gauge invariance and we have checked that the interference between resonant and non-resonant diagrams coming from the off-shell region is large which will again ruin the perturbative computations. There are several techniques to remove these on-shell contributions safely [26][27][28][29]. One can also restrict resonant top momenta out of the on-shell region and can have contribution only from the off-shell region. To implement the last technique with a standard jet veto, one needs a very large number of phase space points to get a stable cross section. The implementation of these techniques is beyond the scope of this paper. Instead of these techniques, we exclude the last two channels by assuming b-quark tagging with 100% efficiency [11,30].

Numerical Results
The sub-process bb → W + W − H gives a significant contribution to the main process We calculate the NLO QCD contribution to this process. In particular we focus on the corrections to cross sections and distributions for various polarization configurations of the final state particles. We also probe variation of cross sections with W W HH anomalous coupling. Some of the Feynman diagrams, tree-level diagrams, as well as one-loop diagrams are heavy vector bosons, Higgs boson and top quark mediated. We use complex-mass scheme (CMS) [31] throughout our calculation to handel the resonance instabilities coming from these massive unstable particles. We take Weinberg angle as The input SM parameters are [32]: There are several pieces in the one-loop calculation which contribute to total σ N LO . As we have discussed above, virtual amplitudes, CT amplitudes, dipole I, P and K terms, dipole subtracted real emission amplitudes contribute to the finite part. We find that there are significant contributions from all these pieces except dipole subtracted real emission amplitudes which gives an almost vanishing contribution.
We use CT14llo and CT14nlo PDF sets [33] for LO (σ LO ) and NLO (σ N LO ) cross sections calculation. We use these PDF sets through LHAPDF [34] libraries. As mentioned before we calculate the cross sections in three different CMEs corresponding to current and proposed future colliders. We choose renormalization (µ R ) and factorization (µ F ) scales dynamically as where p T,W , p T,H are the transverse momenta and M W , M H are the masses of W and Higgs bosons. We measure the scale uncertainties by varying both µ R and µ F independently by a factor of two around the µ 0 given in Eq. 4.1.

Results for the SM
We have listed the cross sections for different CMEs with their respective scale uncertainties in Table 1. As we see in Table 1  The cross section rapidly increases with CME as PDFs for b-quarks are small for lower energies. The relative enhancements RE = σ NLO QCD −σ LO σ LO due to NLO QCD correction are also presented in that table. The RE also increases with CME and it is 19.0%, 25.9% and 31.7% for 14, 27 and 100 TeV CMEs respectively. We have calculated scale uncertainty as the relative change in the cross sections for the different choices of scales within bound 0.5µ 0 ≤ µ R /µ F ≤ 2µ 0 . We see that the NLO uncertainties are a little bit higher than the LO. As there is no strong coupling (α s ) at the Born level, the LO uncertainties come largely from the factorization scale whereas at the NLO the uncertainties come from both, factorization as well as renormalization scales. To see the different scale uncertainties separately, we vary µ R and µ F independently. We see the renormalization scale uncertainty varies from ∼ −10% to ∼ 3% and the factorization scale varies from ∼ −14% to ∼ 18% at NLO depending on CMEs from 14 to 100 TeV.

CME(TeV)
σ As discussed before, we probe the contributions from different polarization configurations of the final state W bosons to the LO and NLO cross sections. The right-handed, left-handed and longitudinal polarization of a W boson are denoted as '+', '-', and '0'. The contributions of different nine polarization combinations of final state W bosons are given in Table 2 for 14, 27 and 100 TeV CMEs. We see that the large contributions are coming from the longitudinal polarization states and among them, the '00' combination gives the largest contribution to the total cross sections. Relative enhancement (RE) for the '00' combination increases with the CME and it becomes ∼ 75% at 100 TeV. In the R ξ gauge, the pseudo Goldstone bosons couple to massive fermions with a coupling proportional to the mass of the fermion. These pseudo Goldstone boson represents the longitudinal polarization state of a W boson. This leads to larger values of the cross section in longitudinal polarization combinations due to heavy fermion mediated diagrams. These longitudinal polarization modes are useful for background suppression to this process. The background may come from the processes with gauge bosons or gluons or photons couplings with light quarks. The negligible masses of the light quarks (u, d, s and c) lead to the suppression of backgrounds in polarization combinations that includes longitudinal polarization.
To find the relative contribution of the bottom-bottom scattering to the pp → W + W − H process, we compute the cross sections in other qq channels along with the bb channel. The results are presented in Table 3. The cross sections in qq channels (4FNS) have been calculated using MagGraph5_aMC5@NLO [35]. MagGraph5_aMC5@NLO cannot compute the one-loop QCD corrections to the bb channel due to the presence of the resonances in the diagrams. As we see in Table 3, the bb channel gives significant contributions to the full process  7,11]. If one adds gg channel, these numbers will be changed accordingly. As we see in Table 3, the corrections are pretty high in qq channels (4FNS). In those channels, MadGraph5_aMC@NLO includes all real emission diagrams and the results are complete but we impose jet veto on b-quarks with 100% efficiency for real emission diagrams to overcome certain difficulties discussed in Sec. 3. The proper inclusion of all real emission diagrams may increase the QCD correction significantly in bb channel.  We have plotted a few different kinematical distributions at the NLO level in Fig. 4 and Fig. 5. In Fig. 4, the upper-panel histograms are for the transverse momentum(p T ) of final state particles at 14 and 100 TeV CMEs. As expected W bosons p T distributions almost coincide with each other. The p T distributions of the Higgs boson is a bit harder. The differential cross sections are maximum around p T = 100 TeV for the Higgs boson and near p T = 80 TeV for the W bosons. In the lower panel of Fig. 4, we have plotted the histograms for the different invariant masses (M ij,ijk ) at 14 and 100 TeV CMEs. Invariant mass thresholds are around ∼ 210, ∼ 170, ∼ 290 TeV and distributions are peaked around 250, 230, 490 TeV for M HW , M W W and M HW W respectively. In Fig. 5, we have plotted differential cross sections with respect to rapidity (η) of final state particles and cosine angle (cos θ) between the two final state particles for 100 TeV CME. The distributions have maxima around η = 0, −0.4 and 0.4 for the Higgs boson, W + and W − boson respectively. From the cos θ plot in Fig. 5, it is clear that maximum contributions come when two final state particles are near to collinear region i.e, θ ∼ 0 or π. In Fig. 6, we have plotted the LO and the NLO distributions to show the effect of the one-loop QCD corrections. The distributions are for only 100 TeV CME. The behavior for the 14 TeV CME is similar. In the upper half of Fig. 6, p T distributions are plotted and in the lower half of Fig. 6, invariant masses have been plotted at 100 TeV CME. We see a increase for the smaller values of the kinematic variables in all the plotted distributions.   Figure 6. The LO and NLO differential cross section distribution with respect to transverse momentums (p T ) and invariant masses (M ij/ijk ) for 100 TeV CME.

Anomalous coupling effect
As we discussed in the introduction, V V HH coupling in SM is only loosely bound so far. We allow W W HH coupling to deviate from the SM value in the search for new physics in the context of κ-framework [36,37]. In the κ framework, only SM couplings deviate by a scale factor. κ is defined as the deviation from the SM coupling. It is a scale factor. Although W W H and W W HH couplings in the SM are related but in many Effective Field Theory frameworks, these couplings can vary independently [4]. As there is no QCD correction to W W HH-vertex, the anomalous coupling will not affect the renormalization. We have checked that the UV and IR poles cancel with the same CTs and dipole terms as in the SM. We denote deviation of V V HH coupling from the SM as κ V 2 H 2 and κ V 2 H 2 = 1 in the SM. In this framework, we vary κ V 2 H 2 from −2.0 to 2.0 and calculate the relative increment RI = σκ−σ SM σ SM in the total cross section, whereas the κ for other SM couplings are set to 1. We choose κ V 2 H 2 = −2.0, −1.0, 1.5, 2.0 and tabulate the results for the LO and NLO cross sections at 14 and 100 TeV CMEs in Table 4. It is clear from Table 4 that cross sections are lower than SM prediction when κ V 2 H 2 is positive and higher than the SM predictions when κ V 2 H 2 is negative. There is not a significant relative increment (−0.8 to + 6.7%) at 14 TeV. At 100 TeV, relative increment vary from −5.6% to +27.8% for the LO cross section and from −7.0% to +26.2% for the NLO cross section. There is also HHH coupling involved in this process. We also observe the HHH anomalous coupling effect on the total cross sections. We vary corresponding κ H 3 from 0.5 to 2.0. We see that there is no significant change in the LO as well as the NLO cross sections and relative increase are smaller than 1% for 14 and 100 TeV CMEs. We see something very interesting in Table 5. The cross sections for the two longitudinally polarized W bosons configuration have stronger dependence on the κ V 2 H 2 . For the NLO cross sections the dependence is almost twice as strong as in the total cross sections. This again demonstrates the importance of measuring the polarization of the W bosons. However this dependence is weaker as compared to the LO cross sections. The difference in this dependence underlines the importance of considering the NLO corrections. In Fig. 7, we have plotted the NLO differential cross section distributions for the Higgs boson and W + boson transverse momenta, and different invariant masses. The maxima of the differential cross sections are about at the same value as for the SM. As there is not that Interesting fact about the negative κ V 2 H 2 is that the distribution are harder. This difference in the shape can be used in putting a strong bound on the coupling. One could put a cut on p W T , or one of the plotted invariant masses to select events with a larger component of anomalous events.

Conclusion
In this letter, we have focused on the NLO QCD corrections to bb → W W H. This process has significant dependence on W W HH coupling. But, the contribution of this process to pp → W W H is only about 15 − 20% of that of light quark scattering. This is where the consideration of the polarization of the W bosons helps. When both the W bosons are longitudinally polarized, then this fraction can increase to 70 − 80%. It turns out that the NLO QCD corrections are also largest for this polarization configuration, making the dependence on the W W HH coupling even stronger. For example, at the 100 TeV CME, the NLO corrections are about 32%, but the corrections are about 75%, when both final state W bosons are longitudinally polarized. Our study suggests that the measurement of the polarization of the final state W/Z bosons can be a useful tool to measure the couplings of the vector bosons and Higgs boson. We have also examined the effect of the variation of κ V 2 H 2 . The variation in the cross section can be twice as large when we consider longitudinally polarized W bosons. In addition, we find that the invariant mass and the p W T distributions are considerably harder for the negative values of κ V 2 H 2 . This can also be useful to put a stronger bound on the coupling. However, to find the bound, one would need to do a detailed background analysis which we leave for the future.