Matrix-Element Corrections to Parton Shower Simulations for Higgs Hadroproduction

We implement matrix-element corrections to HERWIG parton shower simulations for Standard Model Higgs boson production at hadron colliders. We study the Higgs transverse momentum distribution and compare different versions of HERWIG and resummed calculations. The HERWIG results exhibit a remarkable improvement as many more events are generated at large transverse momentum after the inclusion of matrix-element corrections.


Introduction
The Standard Model (SM) of electroweak interactions predicts the existence of the Higgs boson, which is responsible for the mechanism of mass generation. However, such a particle has not been experimentally discovered yet. Searches for the Higgs boson are one of the main goals of the current experiments at the Tevatron accelerator and, ultimately, at future ones, like the Large Hadron Collider (LHC).
In order to perform such searches, precise QCD calculations are mandatory. Other mechanisms may turn out to be useful at hadron colliders [1], yet Higgs production via parton fusion is numerically the most important one, especially at the LHC. Here, the leading-order processes in the strong coupling constant are qq → H and gg → H, with the gluon-gluon fusion mechanism overwhelming the quark-antiquark annihilation channel.
In order to investigate the phenomenology of the Higgs boson, fixed-order calculations may be reliable as long as one only considers inclusive observables, such as total cross sections; for less inclusive quantities, one needs to account for multi-parton radiation in order to perform trustworthy phenomenological analyses [2]. Standard Monte Carlo (MC) algorithms [3,4] describe parton radiation in the soft or collinear approximation, but can have regions of the phase space, so-called 'dead zones', where no radiation is allowed.
In the dead zone, one can however rely on the fixed-order result, as in this region the radiation is neither soft nor collinear enhanced. Several methods have recently been suggested in order to match parton showers and fixed-order matrix elements [5][6][7][8]. In this paper we follow the strategy which has been already used to implement matrix-element corrections to the HERWIG event generator [3] for several processes: e + e − annihilation into quark pairs [9], Deep Inelastic Scattering (DIS) [10], top quark decay [11] and vector boson hadroproduction [12]. The dead zone is here populated by using the exact next-to-leading order (NLO) tree-level matrix element result and the shower in the already-populated region is corrected using the exact amplitude any time an emission is capable of being the hardest so far.
In Section 2 we review the HERWIG parton shower algorithm for the initial-state radiation. In Sections 3 and 4 we discuss the implementation of hard and soft matrixelement corrections, respectively. In Section 5 we present results on the Higgs transverse momentum distributions using different versions of HERWIG and also partonic calculations from the literature. In Section 6 we summarize the main results of our work.

The HERWIG parton shower algorithm
HERWIG simulates the initial-state radiation in hadron collisions according to a 'backward evolution', in which the scale is reduced away from the hard vertex and traces the hard-scattering partons back into the incoming hadrons [13]. The algorithm relies on the universality of the elementary branching probability in the soft or collinear approximation. The probability of the emission of an additional soft/collinear parton from a parton i is given by: The HERWIG ordering variable is Q 2 i = E 2 ξ i , where E is the energy of the parton that splits and ξ i = p·p i EE i , with p and p i being the four-momenta of the splitting and of the emitted parton, respectively; z i is the energy fraction of the outgoing space-like parton with respect to the incoming one; P ab (z) is the Altarelli-Parisi splitting function for a parton a evolving in b. In the massless approximation, ξ i = 1 − cos θ, where θ is the emission angle to the incoming-hadron direction. For soft emission (E i ≪ E), ordering according to Q 2 i corresponds to angular ordering. In (1) f a (x i , Q 2 i ) is the parton distribution function for the partons of type a in the initial-state hadron, x i being the parton energy fraction. The quantity is the Sudakov form factor, expressing the probability of evolution from Q 2 2 to Q 2 1 with no resolvable branching. Unitarity dictates that the Sudakov form factor sums up allorder virtual and unresolved real contributions. In (1), Q imax is the maximum value of Q, fixed by the hard process, and Q c is the value at which the backward evolution is terminated, corresponding, in the case of HERWIG, to a cutoff on the transverse momentum of the showering partons. However, since Q c is smaller than the minimum scale at which the parton distribution functions are evaluated, an additional cutoff on the evolution variable Q i has to be set. If the backward evolution has not resulted in a valence quark, an additional non-perturbative parton emission is generated to evolve back to a valence quark. Such a valence quark has a Gaussian distribution with respect to the non-perturbative intrinsic transverse momentum in the hadron, with a width q T int , that is an adjustable parameter and whose default value is zero.
As the variables Q 2 i and z i are frame-dependent, one needs to specify the frame where the evolution occurs. One can show that, as a result of the Q 2 -ordering, the maximum Qvalues of two colour-connected partons i and j having momenta p i and p j are related via Q imax Q jmax = p i · p j , which is Lorentz-invariant. For Higgs boson production in hadron collisions, symmetric limits are set in HERWIG: Q imax = Q jmax = √ p i · p j . Furthermore, the energy of the parton which initiates the cascade is set to E = Q max = √ p i · p j . Such conditions define the HERWIG frame. It follows that ordering according to Q 2 implies that, in the showering frame, ξ < z 2 .
The region ξ > z 2 is therefore a 'dead zone' for the shower evolution. In such a zone the physical radiation is not soft or collinear enhanced, but not completely absent as it happens in the standard algorithm. It is indeed the purpose of this paper to improve the HERWIG parton shower algorithm and populate the dead zone by including matrixelement corrections.

Hard matrix-element corrections
The Born processes leading to Higgs production via parton fusion at hadron colliders are gg → H (gluon-gluon fusion), of O(α 2 S α), and qq → H (quark-antiquark annihilation), of O(α). In the SM, Higgs production via gluon-gluon fusion is mediated by a quark loop. Hereafter, we shall consider only the top quark contribution in the loop, which is largely dominant in the SM, with finite top mass.
The matrix elements squared for the corrections to gg → H can be found in [14,15], where top mass effects are fully included, and in [16], where the authors have considered the infinite top mass approximation in the loop. HERWIG uses the formulae of Ref. [14]. The actual expressions for such amplitudes are rather involved once the top mass is fully taken into account; therefore, we do not report them here for the sake of brevity. The real NLO corrections to qq → H are instead rather straightforward: the formulae we use can be read from Eq. (3.62) of [17], with appropriate Yukawa couplings and crossing.
In order to implement matrix-element corrections to Higgs production, we follow the prescription of Ref. [18], where a method to include higher-order corrections from real radiation to a generator of the lowest-order process has been proposed.
Referring, e.g., to the correction to gluon-gluon fusion given by we define the partonic Mandelstam variables of process (3) As in [18], we consider a generic Higgs decay H → ℓ 1 ℓ 2 and compute the differential cross sections dσ 3 for the 2 → 3 process gg → gℓ 1 ℓ 2 and dσ 2 for the 2 → 2 one gg → ℓ 1 ℓ 2 via an intermediate Higgs boson.
As done in [12], if we assume that the Higgs squared momentum and rapidity in the Born process are conserved in the transition from the 2 → 2 to the 2 → 3 process, we find that a factorization formula holds. We obtain: where M is the matrix element for the process of Eq. For the gg → H processes we evaluate the parton distribution functions and the strong coupling constant at m 2 H , which is the hard scale of the gg, qq → H processes. For the gg → Hg case, we set such scales to the Higgs transverse mass, m 2 T = q 2 T + m 2 H , q T being the H transverse momentum, which accounts for the additional parton emission via matrix-element corrections. Such choices are the same as the ones done for W/Z production [12].
To include matrix-element corrections, one will have to implement the weight function given by Eq. (3). In order to get the phase space which is populated by HERWIG and the dead zone, where no radiation is allowed, we can repeat all the steps which have been employed for Drell-Yan interactions and report here just the final results obtained in [12]. The total phase-space limits, in terms of the variablesŝ andt, read: where s is the hadronic centre-of-mass energy squared. The valueŝ = m 2 H corresponds to soft radiation, and the linest = 0 andt = m 2 H −ŝ to collinear emission.
As shown in Ref. [12], the HERWIG phase space, which corresponds to ξ < z 2 in terms of the showering variables, is given by: witht min = m 2 H −ŝ −t max . In Fig. 1 it is plotted the physical phase space, along with the region which the showering algorithm populates, for m H = 115 GeV, which is the HERWIG default value, and √ s = 300 GeV. The soft and collinear regions are covered by the shower and one has an overlapping region, where radiation may come from either parton. As expected, Fig. 1 exhibits the presence of a dead zone, where the standard HERWIG algorithm allows no radiation.
Similar arguments are also valid for the other corrections to the gg → H Born process. Nevertheless, a few comments are in order. The processes qg → qH and gq → qH are not soft divergent, but only collinear. This implies that the lower limit for, e.g., gq → qH ist min = m 2 H −ŝ and there is no overlapping region. Likewise, for qg → qH, the upper limit ist max = 0. The process qq → gH, when taken as a correction to gg → H, cannot be interpreted in the parton-shower language as it is neither soft nor collinear singular. However, we included this process as well via matrixelement corrections: the dead zone will then be the full physical phase space given by Eqs. (6)- (7).
In order to implement hard matrix-element corrections, we populate the dead zone using the probability distribution given by the exact matrix element, i.e., Eq. (5), where M will have to be the correct amplitude squared of the hard-scattering process. Moreover, in order to maximize the efficiency of the event generation, the fraction of events generated according to the different subprocesses gg, qg, gq or qq is proportional to the corresponding cross sections in the so-called 'H + jets' process, where the hard process is always one of the corrections to gg → H.
Hard corrections to qq → H processes are similar to what discussed for the gluongluon fusion channel. In particular, the phase-space configuration for qq → gH is as in Fig. 1; the one for qg → qH or gq → gH is the same as for the analogous corrections to gluon-gluon fusion.

Soft matrix-element corrections
The implementation of soft-matrix-element corrections can be performed using the general method of Ref. [5]. Any time an emission in the HERWIG phase space is capable of being the hardest so far, we use the exact matrix element instead of the standard HERWIG algorithm. The hardness of a radiation is measured in terms of the transverse momentum of the emitted parton with respect to the emitting one.
The inclusion of the soft correction is performed by multiplying the parton shower branching probability by a factor which is the ratio of the HERWIG to the matrixelement distribution. It reads: In Eq. (10), J(ŝ,t; z, ξ) is the Jacobian factor for the transformation (z, ξ) → (ŝ,t). As the kinematics for Higgs production are the same as for vector boson production, we can just use for the Jacobian factor the result reported in Ref. [12].
Before closing this section, we would like to point out that matrix-element corrections to Higgs production are differently implemented in the PYTHIA event generator [4]. In fact, in PYTHIA the parton shower approximation is used in all the physical phase space and the exact matrix element corrects only the first emission, rather than the hardestso-far one. Furthermore, the approximation of a top quark of infinite mass is used in PYTHIA to define the ratios of the gg → Hg and qg → Hq to gg → H cross sections, while the latter does use the complete expressions 2 . A comparison of phenomenological results obtained running HERWIG and PYTHIA is given in [2,19].

Transverse momentum distribution
We would now like to present results for the Higgs transverse momentum distribution and investigate the impact of matrix-element corrections. In particular, we wish to compare HERWIG results without and with such corrections as well as the improved HERWIG version with the resummed calculation of Ref. [20] and the so-called 'Monte Carlo at next-to-leading order' (MC@NLO) implementation of [21]. We shall consider Higgs production at the Tevatron and LHC and we shall always assume that the intrinsic transverse momentum is q T,int = 0. Unless otherwise stated, we shall use the default parton distribution functions of HERWIG, but we can anticipate that the relative effect of matrix-element corrections is basically the same, regardless of the chosen set of parton distribution functions.
In Fig. 2 we consider Higgs production at the Tevatron, i.e., pp collisions at √ s = 2 TeV, which is the centre-of-mass energy of the current Run 2. In Fig. 3 we plot instead the q T distribution at the LHC, i.e. pp collisions at √ s = 14 GeV. We consider the HERWIG prediction with (solid histogram) and without (dotted histogram) matrixelement corrections. Beyond q T ≃ m H /2 the matrix-element corrected version allows for many more events. In fact, one can prove that, within the standard algorithm, q T is constrained to be q T < m H : the events at large q T are therefore generated via the same exact amplitude. At small q T , the prediction which includes hard and soft corrections displays a suppression. By default, after matrix-element corrections, the total normalization is still the same and equal to the LO result. It is therefore reasonable that the enhancement at large q T implies a reduction of the number of events which are generated at small transverse momentum. This result was already found for W/Z production [12] and is visible especially at the LHC.
In Fig. 4 we plot the improved HERWIG spectrum (solid line) for the LHC, along    with the result obtained running the H + jets process (dotted line). In order for such a comparison to be reliable, we have turned the qq → H hard process off, as 'H + jets' does not currently implement the corrections to quark-antiquark annihilation. We use the cutoff q T,min = 30 GeV for the 'H + jets' generation. Fig. 4 shows that at small q T the two predictions are fairly different, but at large transverse momentum they agree. This is a reasonable result, since, after matrix-element corrections, large-q T events of both processes are generated via the exact amplitude. This is a check that the implementation of the hard correction is reliable.
Next we compare the new HERWIG version with the resummed partonic calculation of Ref. [20], where the authors resummed soft higher-order contributions to the gg → H process. In fact, the differential distribution dσ/dq T for the production of a Higgs boson of transverse momentum q T presents terms ∼ 1/q 2 T α n S ln m (m 2 H /q 2 T ) which get arbitrarily large at small q T , i.e., for q T ≪ m H . The leading logarithms (LL) correspond to m = n + 1, the next-to-leading logarithms (NLL) to m = n, the next-to-next-toleading logarithms (NNLL) to m = n − 1.
The authors of Ref. [20] have resummed such enhanced logarithms to NNLL level at small q T and matched them to the NNLO result at large q T in order to obtain a reliable result over the full q T range. However, for the sake of comparison with HERWIG, which includes leading logarithms and only some subleading terms (see, e.g., the discussion in Ref. [22] on the comparison of HERWIG and resummation for W/Z production), we use the results of [20] in the NLL approximation, matched to the NLO prediction 3 .
In order for such a comparison to be trustworthy, we have to do the same assumptions as [20]. We use a Higgs mass value m H = 125 GeV and, once we apply matrix-element corrections, in Eq. (5) we evaluate the strong coupling constant and the parton distribution functions at a scale given by the Higgs mass m 2 H . As in [20], we also use the approximation of a top quark with infinite mass in the loop, which is a user-defined option of the HERWIG event generator, and choose the MRST2001 leading-order parton distribution functions [23]. We finally turn the Born qq-initiated processes off.
The result of the comparison is shown in Fig. 5. The normalization and the small-q T behaviour are clearly different. In fact, the total cross section is LO in HERWIG and NLO in our approximation of the work in [20]. The discrepancy at small transverse momenta is instead due to the different logarithmic accuracy which the two considered approaches implement. However, the two curves agree at large q T , where the NLO calculation dominates. This is another independent check of the reliability of the implementation of matrix-element corrections.
Finally, we would like to compare the results of standard HERWIG after matrixelement corrections with the MC@NLO event generator (version 2.2) [21], which implements numerically the method of Ref. [7] to simulate Higgs boson production at hadron colliders. As discussed in [7], the MC@NLO approach implements both real and virtual corrections to the hard-scattering process, in such a way that predicted observables, like total cross sections, are correct to NLO accuracy. Moreover, the MC@NLO showers turn soft matrix-element corrections off.
Version 2.2 of the MC@NLO [21] currently includes only the corrections to Higgs production in the gluon-gluon fusion channel, hence we shall again have to turn the quark-annihilation processes off for the sake of a reliable comparison. In both matrixelement corrected HERWIG and MC@NLO generators we set factorization and renormalization scales for the NLO processes like gg → Hg equal to the Higgs transverse mass.
In Fig. 6 we show the result of this comparison, which exhibits similar features to the one with the resummed calculation. The normalization and the small-q T behaviour of the two curves are different, but nevertheless the large transverse momentum predictions are in good agreement.

Conclusions
We have implemented matrix-element corrections to HERWIG parton shower simulations for direct Higgs production at hadron colliders. We have considered corrections to gg → H and qq → H and have used the exact tree-level NLO matrix element to populate the HERWIG dead zone of the physical phase space and for every hardest-so-far emission in the already-populated region.
We have considered the Higgs transverse momentum q T distribution and have compared HERWIG predictions with and without matrix-element corrections. We have found a remarkable effect of such corrections at large q T , as many more events are generated here via the exact amplitude. We have compared the matrix-element corrected result with the q T prediction yielded by the HERWIG 'H + jets' process and have found agreement at large transverse momentum.
We have compared HERWIG provided with hard and soft corrections with a recent NLL+NLO soft-resummation calculation and, after making consistent choices for the Higgs masses and the scales entering the calculation, we have found very good agreement in the large-q T range.
Moreover, the HERWIG implementation with NLO real corrections has feared rather well against the MC@NLO result, as obtained by using both real and virtual NLO corrections to the hard partonic process. Besides obvious differences in the total normalization (which is LO in standard HERWIG and NLO in the MC@NLO approach) and at small q T , the large-q T spectra agree well, which is another consistency check of the reliable inclusion of matrix-element corrections.
Between the described implementation and the one available within the MC@NLO option, we believe that HERWIG is presently a reliable event generator for Higgs production at hadron colliders both at small and large transverse momentum and that the currently-available matrix-element corrections will play an important role to perform any analysis on Higgs searches at present and future colliders. In particular, the option described here may be the most convenient choice for transverse momentum values q T > ∼ m H . Finally, it should be noted that the HERWIG implementation presented in this paper includes both a finite value for the quark mass in the loop of the relevant gluon-gluon fusion subprocesses and the corrections to the quark-antiquark annihilation channels, thus lending itself to a generalisation of the algorithm to the case of the Minimal Supersymmetric Standard Model, work which is currently in progress.