Contents lists available at ScienceDirect

Physics Letters B

www.elsevier.com/locate/physletb

Gluon fusion contribution to VHj production at hadron colliders

Pankaj Agrawala, Ambresh Shivajib *

a Institute of Physics, Sachivalaya Marg, Bhubaneswar, Odisha 751005, India

b Regional Centre for Accelerator-based Particle Physics, Harish-Chandra Research Institute, Chhatnag Road, Jhusi, Allahabad 211019, India

CrossMark

A R T I C L E I N F 0

Article history:

Received 30 September 2014

Received in revised form 2 December 2014

Accepted 8 December 2014

Available online 12 December 2014

Editor: G.F. Giudice

Keywords: LHC

Gluon fusion

One-loop

A B S T R A C T

We study the associated production of an electroweak vector boson and the Higgs boson with a jet via gluon-gluon fusion. At the leading order, these processes occur at one-loop level. The amplitudes of these one-loop processes are gauge invariant and finite. Therefore, their contributions towards the corresponding hadronic cross sections and kinematic distributions can be calculated separately. We present results for the Large Hadron Collider and its discussed upgrades. We find that the gluon-gluon one-loop process gives dominant contribution to the y Hj production. We observe a destructive interference effect in the gg ^ ZHj amplitude. We also find that in the high transverse momentum and central rapidity region, the ZHj production cross section via gluon-gluon fusion becomes comparable to the cross section contributions coming from quark-quark and quark-gluon channels.

© 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/3.0/). Funded by SCOAP3.

With the discovery of the Higgs boson [1,2] at the Large Hadron Collider (LHC), the standard model and its symmetry-breaking mechanism have been validated. The CMS and ATLAS Collaborations are measuring its properties.The Higgs boson was discovered through gg ^ H production mechanism. Since then signals of other production mechanisms have also been examined. In particular, the associated production of the Higgs boson with an elec-troweak boson has been explored [3,4]. To ensure that the discovered Higgs boson is indeed the standard model Higgs boson, there is a need to detect as many of its signatures as possible. Furthermore, since there is no signal that points to new physics beyond the standard model, there is a need to look for standard model processes that do not have large but accessible cross sections, and can be enhanced/modified by new physics effects.

The LHC and its proposed upgrades provide us an opportunity to explore two types of processes in more detail: (a) the processes which, in the standard model, begin at the one-loop level (one-loop being the leading order (LO)); (b) gluon-gluon scattering processes. As the centre-of-mass energy increases, the gluon-gluon luminosity increases, making many more processes observable. Study of such relatively rare processes is complementary to new physics searches at high energy colliders. Thus, the LHC provides a unique opportunity for testing many of the standard model predictions which was not possible at earlier high energy colliders.

* Corresponding author.

E-mail addresses: agrawal@iopb.res.in (P. Agrawal), ambreshkshivaji@hri.res.in (A. Shivaji).

For example, many one-loop gluon-gluon fusion processes are/will be accessible only at the LHC [5-10]. These gluon fusion one-loop processes can be studied both in the signal and background categories.

In this Letter, we are interested in the gluon-gluon contribution to the pp ^ VHj + X hadronic processes, i.e.,

g + g ^ VHj, (1)

where V is an electroweak vector boson and 'j' stands for a light-quark, or gluon initiated jet. The amplitude of the process gg ^ W± Hj is trivially zero due to the electromagnetic charge conservation. Therefore, V would refer to a photon or a Z boson. This process occurs at the one-loop via triangle, box, and pentagon diagrams. Since quarks carry both the electroweak and color charges, the leading order contribution comes from quark-loop diagrams. The prototype quark-loop diagrams are shown in Fig. 1 [11]. (We have not displayed some of the triangle and bubble diagrams which do not contribute due to the vanishing color factor.) For the Y Hj case, only pentagon class of diagrams contribute. The box diagrams with qqH coupling do not contribute because of the Furry's theorem. Due to the same reason, the leading order gg ^ y H amplitude vanishes. There are a total 24 pentagon diagrams. However, due to charge conjugation symmetry only 12 of those are independent. The full amplitude is proportional to the symmetric color factor dabc.

In the ZHj case, box and triangle diagrams also contribute. Once again there are 12 independent pentagon (PEN) diagrams. There are two types of box diagrams depending on the nature of

http://dx.doi.org/10.1016/j.physletb.2014.12.021

0370-2693/© 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/). Funded by SCOAP3.

Fig. 1. Classes of diagrams contributing to the gg ^ VHj processes. The complete set of diagrams can be generated by permuting the external legs. The leading order gg ^ YHj amplitude receives contribution only from the pentagon diagrams.

the Higgs coupling. The number of independent box diagrams with ZZH coupling (BX2) is 3. Because of the nonsymmetric color factors, in the case of PEN and BX2 amplitudes both the vector and axial-vector pieces of the Z boson coupling with quarks contribute. In the case of box diagrams with qqH coupling (BX1), there are 9 (3 x 3) independent diagrams due to self gluon coupling. There are also 3 independent triangle (TR) diagrams. In the case of BX1 and TR diagrams, contribution of the vector part of the qqZ coupling vanishes due to the Furry's theorem. The vector and axial-vector parts of the full amplitude are proportional to the symmetric color factor dabc and antisymmetric color factor fabc respectively. Except the top quark, all other quarks in the loop are treated massless. Therefore, in the diagrams involving qqH coupling, we take only the top quark contribution. Since the axial-vector part of the amplitude is proportional to the isospin quantum number f', only third generation quarks (which have large mass difference) give contribution. The triangle and box diagrams involving ZZH coupling are anomalous, therefore both the top and bottom quark contributions are taken. We ignore the contribution of the top quark loop of the BX2 diagrams in the case of vector ttZ coupling due to its large mass and no compensating factor from the couplings.

The amplitudes of the gluon fusion processes under consideration have following structure:

M(gg ^ YHj) = i— mVEN, (2)

dabc f abc

M(gg ^ ZHj) = i —Mv + — Mav, (3)

Mv = mVEN - mVX2, (4)

Mav = mPVN + mBX1 - Mbx2 - mAV, (5)

where the subscripts 'V' and 'AV' stand for the vector and axial-vector parts of the amplitude. Since the vector and axial-vector parts of the amplitude have symmetric and antisymmetric color factors, they do not interfere at the amplitude-squared level. Note that the leading order gluon fusion process gg ^ VHj contributes to the corresponding hadronic process at next-to-next-to-leading order (NNLO) in strong coupling parameter as. It also contributes to the QCD next-to-leading order (NLO) correction of gg ^ VH [12-14].

To calculate this amplitude, we follow a semi-numerical approach. As a first step, the traces of gamma-matrices in the prototype diagrams are calculated using FORM [15] in n dimensions.

This demands an implementation of a suitable n-dimensional prescription for y 5 in FORM. The gamma-matrices traces involving y 5 can also be calculated in four dimensions. This will lead to spurious/anomalous terms in amplitudes suffering from chiral anomaly. It is known that such anomalous terms contribute to the quark mass independent rational piece of the amplitude [16,17]. Since the standard model is free from anomaly, the spurious terms are expected to get canceled between the up-type and down-type quark contributions of a given quark generation. Since all the gamma-matrices in the trace are contracted either with the momenta or polarizations, the trace calculation does not lead to any explicit n dependence. The amplitude at this stage can be cast in terms of tensor integrals. It has rank-four, five-point functions as the most complicated tensor integrals. In standard notations, a rank-four, five-point tensor integral can be written as

E^vpa =

dni l^iv IP la

(2n)n dod1d2d3d4'

where l is the loop momentum, n is the space-time dimension and di 's are loop momentum dependent scalar propagators. The tensor integrals are known to be expressible in terms of scalar integrals [18]. This is the most important part of any multileg one-loop computation. We reduce pentagon tensor integrals into box tensor and scalar integrals using the fact that in four dimensions the loop momentum (l) can be expressed as a linear combination of four independent external momenta which are available in a five-point (2 ^ 3) process. The same fact is used to write a pentagon scalar integral in terms of five box scalar integrals [17,19]. Since the pentagon tensor integrals are ultraviolet (UV) finite, their reduction into box tensor and scalar integrals can be carried out in four dimensions consistently [20,21]. We follow the methods of Oldenborgh and Vermaseren (OV) [22] to reduce the box and lower-point tensor integrals. This part of the reduction is performed in n = 4 — 2e dimensions to regulate the bad ultraviolet behaviour of tensor integrals. In OV method, the organization of tensor integrals in terms of generalized Kronecker deltas promises a greater numerical stability of the reduction procedure. The complete one-loop tensor reduction library, OVReduce, is developed (in FORTRAN) to carry out the numerical reduction of one-loop tensor integrals required in 2 ^ 3 processes [7]. The singular structure of any one-loop amplitude is, thus, encoded in various scalar integrals of box, triangle, bubble and tadpole types. We use OneLOop library [23] to compute all the required scalar integrals. It uses dimensional regularization to regulate both the UV and infrared (IR) singularities. We require only one prototype pentagon and BX2 amplitudes to generate the vector part of the full amplitude. On the other hand, to generate the full axial-vector part of the amplitude, we require three prototype pentagon amplitudes, three BX1 amplitudes, one BX2 amplitude and one TR amplitudes. This is to ensure that permutation of momenta and polarizations is consistent with the permutation of external legs and their couplings with the massive quarks.

We perform many checks to ensure the correctness of the VHj amplitudes. The amplitudes are expected to be both ultraviolet and infrared (due to massless quarks) finite. This is an important check on the amplitudes. The IR finiteness holds for each quark loop diagram [24]. The UV finiteness is expected only in gauge invariant subamplitudes, however, there are amplitudes which are individually UV finite from naive power counting. For example, each pentagon amplitude and each axial-vector part of BX1, BX2 and TR amplitudes are UV finite. As an ultimate check, we have checked the gauge invariance of amplitudes with respect to all the gauge currents. We do it numerically by replacing their polarization vectors with their respective 4-momenta, e^(k) ^ k^. The vector part of the amplitude is gauge invariant with respect

Table 1

A comparison of parton channels contributing to pp -— yHj hadronic cross section. The numbers in bracket in 100 TeV column are with a minimum p'T cut of 50 GeV.

VS (TeV) 8 13 14 33 100

Y Hj, LO r , -, ffGG № YHj,LO r , -, aYQ+QG tab] 72 26 268 81 320 96 2029 491 13 435 (12 734) 2667 (1312)

to the Z boson. However, due to the explicit breaking of the chi-ral symmetry in presence of massive quarks, the axial-vector part of the amplitude is not gauge invariant with respect to the Z boson. We have checked explicitly that the quantum anomaly is canceled between the top and bottom quark contributions. The three-body phase space is generated using the phase space generator RAMBO [25]. Because of a very large and complicated expression of the amplitudes, we calculate the (helicity/polarized) amplitudes before squaring them. The parton level phase space integration and convolution of partonic cross section with the gluon distribution functions are done using the Monte Carlo integration method based on VEGAS [26]. To minimize the computation time, we perform the evaluation in a parallel environment with the help of the AMCI package [27] based on the PVM [28] parallel networking software. The computation of multileg, one-loop amplitudes often suffers from the issue of numerical instability due to the small Gram determinants for certain phase space points. Since the number of such phase space points is not very large, we systematically ignore their contribution by implementing a set of Ward identities. More details on this can be found in Refs. [29,30].

We now present the numerical results for the VHj processes. We have taken mt = 173 GeV and MH = 126 GeV. Following kinematic cuts are applied to obtain the results discussed below,

pJT > 30 GeV, pT > 20 GeV,

|yj|< 4.5, |yY|< 2.5, ARyj > 0.4. (7)

We use cteq6l1 parton distribution functions [31] and choose [R = [if = MH as the common central scale for renormalization and factorization. In Table 1, we compare the LO gluon-gluon (GG) contribution to pp — y Hj with the contribution from the leading order quark-quark (QQ) and quark-gluon (QG) initiated process at various collider centre-of-mass energies. The contribution from QQ+QG channel involving YYH and ZyH one-loop couplings is expected to be negligible because of the large electroweak coupling suppression. The leading order QQ+QG channel contribution is calculated using MG5 [32]. Due to small bottom quark mass

and low bottom quark flux in proton, the QQ+QG contribution is smaller than the GG contribution. However, the hadronic cross section in the gluon-gluon fusion channel itself is very small (~0.3 fb at 14 TeV) due to the top quark propagator suppression in the pentagon amplitudes. With this cross section, in run II of the LHC, we would expect about 100 such events. We also see that the ratio of GG channel and QQ+QG channel cross sections increases as the centre-of-mass energy increases. It is expected due to the faster increase in the GG luminosity with centre-of-mass energy. With proper background suppression strategies, one may be able to observe this process at HL-LHC and HE-LHC. In Fig. 2, we display the transverse momentum ( pT ) and rapidity ( y) distributions of the final state particles in the y Hj case. Note that the pT distributions for both the photon and the jet are well defined in the collinear region (pT ^ 0). This is so because the amplitude for gg ^ y Hj receives contribution only from massive quark loop diagrams.1

We have seen that unlike in the case of y Hj, the ZHj production via gluon-gluon fusion also involves box and triangle diagrams with both massive and massless quarks in the loop. At 14 TeV, the ZHj cross section is about 73 (64) fb with cteq6l1 (cteq6m) parton distributions and fixed scale choice ^ = MH. We observe that almost all the contribution to the cross section comes from the axial-vector or color antisymmetric part of the amplitude. At 14 TeV LHC, the vector and axial-vector contributions are, aV = 0.10 fb and aAV = 73.36 fb respectively. There are many gauge invariant subamplitudes in the gg ^ ZHj amplitude. The nature of interference among them can be considered as an important prediction of the standard model. We have two sets of gauge invariant subamplitudes - one involving ttH coupling and the other involving ZZH coupling. Like in many associated Higgs production processes, we find that there is a very strong destructive interference between the two sets of amplitudes [33,34]. To illustrate this feature, in Fig. 3, we have given pT distributions of the Higgs boson considering only one gauge invariant set of amplitudes at a time. This is very similar to the destructive interference effect seen in gg ^ ZH amplitude; it is also shown in the right panel of Fig. 3. To quantify the interference effect in the ZHj production, we have calculated the contributions of different gauge invariant sets towards the hadronic cross section. At 14 TeV, the cross sections are, attH = 52.56 fb, ctzzh = 211.78 fb, and CttH+ZZH = 73.56 fb. In presence of new physics in the Higgs sector, such a strong interference effect may lead to a very dif-

1 A collinear singular configuration involves three massless lines at a vertex.

gg>zHj

g g >z h

ttH+ZZH ttH ZZH

14 TeV LHC cteq6M h = mh

prh [gev]

pth [gev]

Fig.3. pT distributions of the Higgs boson displaying the strong interference effect between the gauge invariant subamplitudes involving ttH and ZZH couplings in gg ^ ZHj and gg ^ ZH processes.

gg>ZHj

Fig. 4. Variation of ZHj cross section with the anomalous parameters CzzH and CttH at VS = 14 TeV. The central line corresponds to CttH = 1, while the vertical spread about this line is due to the variation of CttH between 0.80 and 1.20, CttH = 1.20 being at the lower end.

ferent predictions for both the cross section and kinematic distributions. For example, if we introduce deviations (due to new physics) in the ttH and ZZH couplings through the scale factors CttH and Czzh respectively, the ZHj cross section at 14 TeV can be parametrized as

_ pp— ZHj GG

14 TeV

= 52.56Ct2tH + 211.78C|ZH - 190.78CttH C ZZH.

The standard model prediction corresponds to C ttH = CZZH = 1. In Fig. 4, we give the variation of this cross section by changing the anomalous parameters by 20% around their standard model values. Note that the effects of these parameters on the cross section is opposite in nature within the range they are varied. Also, because of a larger contribution from ZZH coupling diagrams which is due to the inclusion of light quark diagrams, the net cross section is more sensitive to variation in CZZH.

In Fig. 5, the pT and rapidity distributions for the jet is shown. Since the final state gluon can originate as radiation from the initial state gluons (as in BX1 and TR diagrams), the majority of events come from low pT region of the jet. We have calculated the gg ^ ZHj hadronic cross section at 8, 13, 14, 33 and 100 TeV collider centre-of-mass energies. In going from 8 TeV to 14 TeV, the cross section increases by about a factor of 5 due to an increase

in the gluon flux at a higher energy. In Table 2, we have compared the gluon-gluon (GG) contribution to pp — ZHj with the LO and NLO quark-quark (QQ) and quark-gluon (QG) predictions at various collider centre-of-mass energies.The NLO correction to this process has been computed in recent past [35,36]. We have obtained the NLO results using MG5_aMC@NLO [37] with cteq6m parton distributions. As expected, the relative importance of the GG contribution increases with the increasing collider energy. To illustrate this, we have displayed two ratios (= ^QH^/^GG^'1^) and R2 (= (aQQ+QGL° - aQQ+Qc)/aGGHj'LO) in Table 2. The ratio R2 suggests that even at 8 TeV, the GG contribution becomes comparable to the NLO enhancement of QQ+QG contribution. Nevertheless, the net contribution from the QQ+QG channel remains dominant at all energies. In the 14 TeV column, the numbers are presented with the percentage scale uncertainties. For that we have varied the scale by a factor of 2 about the central value ¡ = MH. A large scale uncertainty in the GG predictions is typical to many other gluon-gluon fusion processes which can be reduced by including the higher order corrections. One example of such a process is gg — H which has been widely studied in the literature [38]. The GG numbers for 100 TeV should be interpreted carefully. In the last row of the table, we also report gg — ZH hadronic cross sections. Note that, the gg — ZHj cross section exceeds gg — ZH cross section at VS = 100 TeV for a minimum p]T cut of 30 GeV. This is a reflection of the fact that at such high energies, the cross section is also influenced by the large-log terms like ln( p2/S). Due to this, the fixed-order perturbative cross section prediction cannot be reliable below certain pT and one needs to resum these large logarithms. Therefore, in the table, we have also reported the cross section at 100 TeV with p'T > 50 GeV.

We have also compared the normalized pT and rapidity distributions of the Z boson at LO(GG) and NLO(QQ+QG). It is shown in Fig. 6. The difference between the two contributions is consistent with the fact that at higher pT, the region of parton distributions with smaller momentum fraction x is easily accessible and gluon distributions are more important in this region. This comparison suggests that at high pT and low rapidity the GG contribution is relatively more important along with the dominant QQ+QG channel contribution. The pT and rapidity distributions of the Higgs also show a similar characteristic difference, however, the same is not true for the distributions of the jet.

In conclusion, we have studied the gluon-gluon fusion contribution to the pp — VHj process, where V e[y , Z}. It is a NNLO contribution in as. We find that the hadronic cross section for

gg>zHj

gg>zHj

y.O —1- -1- -i-1-1-1--- -1—1—1—1—1— ■ i i i i i i i i i i i 18

1.8 ■ 1 16

1.6 ■ 14

1.4 1.2 14 TeV LHC ■ 12

1.0 0.8 ■ cteq6M . H = MH e T3 10 8

0.6 ■ 6

0.4 ■ 4

0.2 0.0 ■ 2

3 50 100 150 200 250 300 0

pT[GeV]

2 3 4 5

Fig. 5. Transverse momentum and rapidity distributions of the jet in gg — ZHj.

Table 2

A comparison of partonic channels contributing to pp -— ZHj hadronic cross section. We use cteq6l1 (cteq6m) for LO (NLO) results and = = MH. The numbers in bracket in 100 TeV column are with a minimum p'T cut of 50 GeV.

VS (TeV)

13 14 33 100

63.93 226.80 266.90 3.55 73.56+525 co+11.7% 256.52-9.9% 302.60+3.4% 3.49 535.30 955.00 1076.00 1.78 4270.71 (2788.90) 4143.34 (2715.82) 4414.00 (3087.21) 0.97 (0.97)

0.63 0.63 0.23 0.06 (0.13)

82.38 97.70+32.1% 97.70-22.7% 569.90 3764.63

ZHj,LO .r. . aQQ+QG [fb]

ZHj, NLO aQQ+QG

ZHjjLO n aQQ+QG A1 = ZHjjLO

ZHj, NLO_ ZHjjLO QQ+QG aQQ+QG

14.96 95.12 114.20 6.36

1.27 24.60

pp>ZHj

pp>ZHj

£ ■a e

0.35 0.30 0.25 0.20

■g 0.15

50 100 150 200 250 300 350 400 450 500

pTz[GeV]

0.10 0.05 0.00

LO(GG) - :

■ 14 TeV LHC NLO(QQ+QG) ---------- :

: n = MH H

■ .-J l""L-' :

: J" i-T 1 Li '

: rJ H Vi ^

Fig. 6. A comparison of the LO (GG) and NLO (QQ+QG) contributions to the transverse momentum and rapidity distributions (normalized by cross sections) of the Z boson.

Y Hj is only about 0.3 fb at 14 TeV LHC, however, it is much larger than the leading order quark-quark and quark-gluon contribution which is only about 0.01 fb at 14 TeV. This raises the possibility of observing the production of the Higgs boson in association with a photon with sufficient integrated luminosity at the LHC. In the ZHj case, apart from the GG fusion contribution, we have used MG5_aMC@NLO to estimate the NLO corrections to the QQ+QG contribution. We find that at 13 or 14 TeV centre of mass energy, the gluon-gluon contribution is comparable to the NLO corrections. At higher centre-of-mass energies, the GG contri-

bution becomes even more important and can become comparable to the QQ+QG contributions. In the ZHj case, we also observe a destructive interference between the two gauge invariant sets of amplitudes involving ttH and ZZH couplings. Any modification to these couplings due to new physics effects can spoil the interference effect and lead to a very different prediction. A suitable pT and rapidity cuts can be employed to enhance the signal from this channel while comparing it with quark-quark and quark-gluon channels. For realistic predictions the parton shower effects should also be included which is beyond the scope of the present work.

ZHj, LO

ZHj LO GG

ZH, LO

Like many other gluon-gluon fusion LO processes, the processes considered here too suffer from a large scale uncertainty. By computing higher order QCD corrections, the scale uncertainty can be reduced. However, this will be quite challenging at the NLO as it requires evaluation of one-loop six-point functions and two-loop five-point functions.

Acknowledgements

We would like to acknowledge the use of Cluster Computing facility available at the Institute of Physics, Bhubaneswar. A.S. would like to acknowledge the hospitality of the institute where the present work was initiated. AS would also like to thank Manoj Mandal, Marco Zaro and Rikkert Frederix for their technical help with the MG5_aMC@NLO package. The work of A.S. is partially supported by funding available from the Department of Atomic Energy, Government of India, for the Regional Centre for Accelerator-based Particle Physics, Harish-Chandra Research Institute.

References

[1] G. Aad, et al., ATLAS Collaboration, Phys. Lett. B 716 (2012) 1, arXiv:1207.7214 [hep-ex].

[2] S. Chatrchyan, et al., CMS Collaboration, Phys. Lett. B 716 (2012) 30, arXiv: 1207.7235 [hep-ex].

[3] See, e.g S. Chatrchyan, et al., CMS Collaboration, Phys. Rev. D 89 (2014) 012003.

[4] See, e.g., P. Agrawal, S.D. Ellis, Phys. Lett. B 229 (1989) 145;

P. Agrawal, D. Bowser-Chao, K. Cheung, Phys. Rev. D 51 (1995) 6114.

[5] P. Agrawal, A. Shivaji, Phys. Rev. D 86 ( 2012 ) 073013, arXiv:1207.2927 [hep-ph].

[6] D. de Florian, Z. Kunszt, Phys. Lett. B 460 (1999) 184, arXiv:hep-ph/9905283.

[7] P. Agrawal, G. Ladinsky, Phys. Rev. D 63 (2001) 117504, arXiv:hep-ph/0011346.

[8] J.M. Campbell, R.K. Ellis, E. Furlan, R. Rontsch, arXiv:1409.1897 [hep-ph].

[9] T. Melia, K. Melnikov, R. Rontsch, M. Schulze, G. Zanderighi, J. High Energy Phys. 1208 (2012) 115, arXiv:1205.6987 [hep-ph].

[10] F. Campanario, Q. Li, M. Rauch, M. Spira, J. High Energy Phys. 1306 (2013) 069, arXiv:1211.5429 [hep-ph].

[11] D. Binosi, J. Collins, C. Kaufhold, L. Theussl, Comput. Phys. Commun. 180 (2009) 1709, arXiv:0811.4113 [hep-ph].

[12] D.A. Dicus, C. Kao, Phys. Rev. D 38 (1988) 1008;

D.A. Dicus, C. Kao, Phys. Rev. D 42 (1990) 2412 (Erratum).

13] B.A. Kniehl, Phys. Rev. D 42 (1990 ) 2253.

14] L. Altenkamp, S. Dittmaier, R.V. Harlander, H. Rzehak, T.J.E. Zirke, J. High Energy Phys. 1302 (2013) 078, arXiv:1211.5015 [hep-ph].

15] J.A.M. Vermaseren, arXiv:math-ph/0010025.

16] A. Shivaji, V. Ravindran, P. Agrawal, J. High Energy Phys. 1202 (2012) 057, arXiv:1111.6479 [hep-ph].

17] A.K. Shivaji, arXiv:1305.4926 [hep-ph].

18] G. Passarino, M.J.G. Veltman, Nucl. Phys. B 160 (1979) 151.

19] W.L. van Neerven, J.A.M. Vermaseren, Phys. Lett. B 137 (1984) 241.

20] Z. Bern, L.J. Dixon, D.A. Kosower, Nucl. Phys. B 412 (1994) 751, arXiv:hep-ph/ 9306240.

21] T. Binoth, J.P. Guillet, G. Heinrich, E. Pilon, C. Schubert, J. High Energy Phys. 0510 (2005) 015, arXiv:hep-ph/0504267.

22] G.J. van Oldenborgh, J.A.M. Vermaseren, Z. Phys. C 46 (1990) 425.

23] A. van Hameren, Comput. Phys. Commun. 182 (2011) 2427, arXiv:1007.4716 [hep-ph].

24] A. Shivaji, arXiv:1008.4792 [hep-ph].

25] R. Kleiss, W.J. Stirling, S.D. Ellis, Comput. Phys. Commun. 40 (1986) 359.

[26] G.P. Lepage, VEGAS: an adaptive multi-dimensional integration program, Cornell preprint CLNS 80-447, March 1980.

27] S. Veseli, Comput. Phys. Commun. 108 (1998) 9.

28] A. Geist, et al., PVM: Parallel Virtual Machine, MIT Press, Cambridge, MA, 1994.

29] F. Campanario, J. High Energy Phys. 1110 (2011) 070, arXiv:1105.0920 [hep-ph].

30] P. Agrawal, A. Shivaji, J. High Energy Phys. 1301 (2013) 071, arXiv:1208.2593 [hep-ph].

31] P.M. Nadolsky, H.L. Lai, Q.H. Cao, J. Huston, J. Pumplin, D. Stump, W.K. Tung, C.-P. Yuan, Phys. Rev. D 78 (2008) 013004, arXiv:0802.0007 [hep-ph].

32] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, T. Stelzer, J. High Energy Phys. 1106 (2011) 128, arXiv:1106.0522 [hep-ph].

33] P. Agrawal, S. Mitra, A. Shivaji, J. High Energy Phys. 1312 (2013) 077, arXiv:1211.4362 [hep-ph].

34] K. Nishiwaki, S. Niyogi, A. Shivaji, J. High Energy Phys. 1404 (2014) 011, arXiv:1309.6907 [hep-ph].

[35] S. Ji-Juan, G. Lei, M. Wen-Gan, Z. Ren-You, L. Liu-Sheng, H. Liang, J. High Energy Phys. 1203 (2012) 059, arXiv:1202.6422 [hep-ph].

36] G. Luisoni, P. Nason, C. Oleari, F. Tramontano, J. High Energy Phys. 1310 (2013) 083, arXiv:1306.2542 [hep-ph].

[37] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H.-S. Shao, T. Stelzer, et al., J. High Energy Phys. 1407 (2014) 079, arXiv:1405.0301 [hep-ph].

[38] S. Dittmaier, et al., LHC Higgs Cross Section Working Group Collaboration, arXiv:1101.0593 [hep-ph], and relevant references therein.