Gluon-induced Higgs-strahlung at next-to-leading order QCD

Gluon-induced contributions to the associated production of a Higgs and a Z-boson are calculated with NLO accuracy in QCD. They constitute a significant contribution to the cross section for this process. The perturbative correction factor (K-factor) is calculated in the limit of infinite top-quark and vanishing bottom-quark masses. The qualitative similarity of the results to the well-known ones for the gluon-fusion process $gg\to H$ allows to conclude that rescaling the LO prediction by this K-factor leads to a reliable NLO result and realistic error estimate due to missing higher-order perturbative effects. We consider the total inclusive cross section as well as a scenario with a boosted Higgs boson, where the Higgs boson's transverse momentum is restricted to values ptH>200GeV. In both cases, we find large correction factors $K\approx 2$ in most of the parameter space.


Introduction
With the recent observation of a new particle at the LHC [1,2] and the related evidence at the Tevatron [3], efforts to determine its identity are of highest priority. Among the most important observables are the total and differential cross sections. First measurements of these quantities indicate that the new particle is indeed the long-sought Higgs boson of the Standard Model (SM). In order to definitely confirm or exclude this hypothesis, accurate measurements and corresponding precision calculations of the cross section in the various production modes are required.
The current theoretical knowledge of the SM cross sections is in general quite impressive and documented in Refs. [4,5]. Subject of the current paper is a particular contribution to the so-called Higgs-strahlung process pp → HV (V = W, Z). While it has been a major search mode for Higgs bosons at the Tevatron, it used to be considered of minor importance at the LHC due to its small cross section and large background. However, it belongs to the channels that were analysed by the ATLAS and the CMS experiments already with the first data. The signal-to-background ratio for Higgs-strahlung can be significantly enhanced when cutting on events where the Higgs boson is produced at large p T,H [6].
The leading-order (LO) cross section for this process can be written as a convolution of the cross section for the Drell-Yan process pp → V * with the decay rate for V * → HV , where V * denotes an off-shell gauge boson of momentum k: σ HV,DY (pp → HV ) = dk 2 This relation holds exactly through next-to-leading order (NLO) QCD, i.e. O(α s ), and approximately through next-to-next-to-leading order (NNLO). The QCD effects of Eq. (1) are therefore strongly dominated by the Drell-Yan corrections to the cross section σ DY ; they are known through NNLO QCD for the total HW/HZ cross sections [7][8][9], and for HW production also differentially [10]. Typical Feynman diagrams for the Drell-Yan type contribution are shown in Fig. 1 (a-f). They contribute to the cross section at order g 4 α n s (n = 0, 1, 2) and increase it by about 30% with respect to LO. Here and in what follows, α s = g 2 s /(4π), with g s the strong and g the weak coupling constant. Apart from the Drell-Yan-like QCD corrections at NNLO, there are top-loop-induced contributions such as the ones shown in Fig. 1 (g-j). Their interference with the LO and the real-emission NLO amplitude is of order λ t g 3 α 2 s , with λ t the top Yukawa coupling, and their numerical impact is at the percent level [11].
In contrast to the NLO QCD and dominant NNLO QCD corrections, electroweak (EW) corrections do not respect a factorization into Drell-Yan-like production and decay, since irreducible (box) corrections to qq ( ) → HV already contribute at NLO. The NLO EW corrections have been evaluated in Ref. [12] for the total HV cross sections, where they amount to −(5−10)%, and in Ref. [13] for differential distributions as part of the HAWK Monte Carlo program, which fully includes all decays and off-shell effects of the weak boson V = W, Z. In distributions the EW corrections can grow to −(10−20)%. As suggested in Ref. [14], NLO EW and Drell-Yan-like NNLO QCD corrections can be conveniently combined in factorized form, where the EW corrections modify the QCD prediction by a relative correction factor that is rather insensitive to the parton luminosities.
Recently, QCD corrections to the H → bb decay have been considered as well [15]. These final-state corrections should be carefully taken into account in the Higgs reconstruction.
In this paper we focus on another type of contribution which is specific to HZ production, namely gluon fusion, mediated by top-and bottom-quark loops. Typical diagrams of this channel are shown in Fig. 2. Owing to the initial-state gluons, it cannot interfere with the LO amplitude and therefore contributes to the cross section at order λ 2 t g 2 α 2 s . For M H = 125 GeV, at leading, i.e., one-loop order it amounts to about 4% (6%) of the total Higgs-strahlung cross section at the LHC with 8 TeV (14 TeV) [9]. Since it has no lower-order correspondence, it is separately gauge invariant and IR and UV finite. The two initial-state gluons lead to a rather strong renormalization and factorization scale dependence of about 30%, thus increasing the theoretical uncertainty of the HZ relative to the HW process, where the gg channel does not exist (at this order). Experience from the gluon-fusion process gg → H shows, however, that the LO scale uncertainty drastically underestimates the actual size of the higher-order corrections. Owing to the similarity of the gg → H and the gg → HZ processes in their QCD structure (same initial states and colour structure, both loop-induced), we expect a similar phenomenon in the latter.
The goal of the present paper is to improve on the theory uncertainty of the gg → HZ process by calculating its NLO QCD corrections. Note that they are of order α 3 s and thus formally contribute to the N 3 LO corrections of the Higgs-strahlung process. Technically the described NLO calculation involves massive, multi-scale two-loop diagrams that are beyond present calculational techniques, so that we are forced to employ asymptotic expansions in the limit of a large top-quark mass. We note that the same strategy was already successfully applied to the calculation of NLO corrections to the related process of Higgs pair production via gluon fusion, gg → HH [16].
The paper is organized as follows: In Section 2 we briefly outline the problem, before describing the details of our calculation in Section 3. Our numerical results are discussed in Section 4, and our conclusions given in Section 5.
2 Outline of the problem 2.1 Leading order At LO and in covariant R ξ gauge, the Feynman diagrams contributing to the gluon-induced Higgs-strahlung process can be divided into three types, shown in Fig. 2: (a) Box diagrams for gg → HZ: Only massive quarks run in the loop due to the proportionality to the respective Yukawa coupling. Note, however, that these graphs tend to zero also in the heavy-quark limit.
(b) Triangle diagrams for gg → Z * → HZ: Owing to Furry's theorem, all contributions from vector couplings compensate each other, so that only the axial-vector coupling of the Z boson needs to be taken into account. Since the axial-vector coupling is proportional to the third component of the weak isospin of the quark (± 1 2 ), the contribution of a single quark generation vanishes in the equal-mass case. Assuming massless quarks in the first two fermion generations, this leaves a non-vanishing contribution only from the third generation. The amplitudes tend to zero in the heavy-quark limit.
It is interesting to note that only the longitudinal part of the Z-boson propagator contributes, while all contributions of the transverse part vanish. This consequence of the Landau-Yang theorem [17,18] can be used at NLO to facilitate the calculation significantly, as will be described below.
(c) Triangle diagrams for gg → G 0 → HZ: Only the massive-quark loops contribute here, where G 0 is the would-be Goldstone boson partner to the Z boson. The graphs are both proportional to the respective Yukawa coupling and to the third component of the weak isospin of the quark and tend to a constant in the heavy-quark limit.
While the box diagrams (a) are gauge-parameter independent in the R ξ gauge, both the vertices (b) and (c) depend on the gauge-parameter of the Z boson. The sum of (b) and (c) for each quark generation is, of course, gauge-parameter independent.
The full result for the LO amplitudes for the process gg → HZ can be found in Ref. [19]; the hadronic cross section can be easily obtained using the program vh@nnlo [9,20]. We have rederived the LO cross section with the full dependence on the top-and bottom-quark masses as a basic ingredient of our NLO calculation.

Next-to-leading order
The Feynman diagrams for the NLO QCD corrections to the gg → HZ process are obtained from the LO gluon-fusion diagrams shown in Fig. 2 (a-c) by attaching virtual and real gluons and quarks to internal and external quark and gluon lines in all possible ways: Figure 1: Representative diagrams to hadronic HZ production of Drell-Yan type up to NNLO (a-f) and non-Drell-Yan-like NNLO graphs with Higgs radiation off top-quark loops; both types of corrections (up to NNLO) are not considered in this publication.
(m) Figure 2: Representative diagrams to hadronic HZ production via quark-loop-induced gluon fusion. It is understood that crossed diagrams have to be taken into account as well.
• For the real corrections, this results in triangle diagrams with two massive external momenta, Fig. 2 (g,h), box diagrams with one or three massive external momenta, Fig. 2 (i) and (d,e), and pentagon diagrams with two massive external momenta, Fig. 2 (f) (counting off-shell gluons as massive lines). Note that Fig. 2 (e) is a crossed version of Fig. 1 (j), for example; as pointed out above, it can interfere at O(λ t g 3 α 2 s ) with the NLO real-emission Drell-Yan type amplitude, which has already been taken into account in Ref. [11]. In the present paper, we work at O(λ 2 t g 2 α 3 s ) and need to evaluate the square of such terms.
• For the virtual corrections, we encounter two-loop vertex and box diagrams, Fig. 2 (j,l,m), as well as one-particle-reducible diagrams with two one-loop triangle insertions, Fig. 2 (k).

Details of the calculation and effective-field-theory approach
While the majority of the integrals could be calculated using well-known techniques, a general result for the massive double-box integrals shown in Fig. 2 (j) is beyond current technology. However, motivated by observations made in Refs. [21][22][23][24][25][26][27], for example, we follow a strategy that has been successfully applied to higher-order corrections to Higgs production via gluon fusion. Instead of calculating the Feynman integrals in full generality, we determine the perturbative correction factor in the limit of infinite top-quark and vanishing bottom-quark masses (referred to as "effective theory" in what follows). For the gluon-fusion process, both inclusive and differential, it turns out that this factor is rather insensitive to the top-quark mass effects [21][22][23][24][25][26][27].
Using asymptotic expansion of Feynman diagrams [28,29], the heavy-top limit can be interchanged with the loop integration, which simplifies the calculation enormously.

(i) LO amplitude
At LO, the diagrams with top-quark loops reduce to vacuum diagrams (integrals with vanishing external momenta) in the effective theory. Because already at LO the loop integrals are UV divergent, some care is needed in the calculation of the Dirac traces that involve the matrix γ 5 . Both at LO and NLO, we consistently use the 't Hooft-Veltman scheme [30,31], where γ 5 anticommutes with the first four, but commutes with all other Dirac matrices. In practice, we insert γ 5 = − i 4! µνρσ γ µ γ ν γ ρ γ σ ( 0123 = +1), keep the -tensor outside of the D-dimensional integration, and project onto four dimensions only after all divergent terms have cancelled among each other.
The result for the LO amplitude M 0 and its polarization-and colour-averaged square is with M H the Higgs mass and M Z the Z mass, α and α s the electromagnetic and the strong coupling constants, s 2 w = 1 − c 2 w the sine of the weak mixing angle, and λ(a, b, c) = , with a, b denoting colour indices andŝ = (p 1 + p 2 ) 2 the usual Mandelstam variable. The polarization vectors ε i (p i ) (i = 1, 2) correspond to the respective incoming gluons, and ε * Z (p Z ) to the outgoing Z boson. The shorthand (ε 1 , ε 2 , p 1 , p 2 ) stands for the contraction of the -tensor with the 4-vectors in the argument.

(ii) NLO virtual corrections
The large-mass expansion [28,29], formulated in terms of the "method of regions", states that the loop integrand is to be Taylor-expanded in all relevant regions of loop momenta, and the final result of each diagram is the sum over all those expansions. If a single internal mass m is considered large, a loop momentum q µ i can either be large, In this way, the virtual NLO (two-loop) diagrams with top-quark loops reduce to either massive twoloop vacuum integrals or products of massless one-loop triangles with massive one-loop vacuum integrals. As a result of the large-mass expansion, all one-particle irreducible twoloop graphs involving top-quark loops vanish, except for the ones with corrections to the ggG 0 vertex. Note that the latter as well as the axial vector part of the ggZ vertex each receive an anomalous counterterm [32] to restore chiral symmetry in the massless-quark limit in the 't Hooft-Veltman scheme [30,31] for γ 5 , in the following denoted as δZ P 5 and δZ A 5 , respectively. Since we set m q = 0 for q = t both in the propagators and in the Yukawa couplings, the set of diagrams to be evaluated with internal massless quarks does not contain any two-loop box diagrams. The only genuine two-loop integrals with non-vanishing external momentum are three-point functions as the one shown in Fig. 2 (l). As their LO counterparts at one loop, their contribution to the cross section vanishes when working in Landau gauge, where the Z propagator D µν (k) is proportional to the polarization sum λ for a physical vector particle Z * of mass k 2 , The ggZ * subamplitude therefore corresponds to the decay of a massive into two massless vector particles which is forbidden due to the Landau-Yang theorem [17,18].
In summary, the one-particle-irreducible two-loop diagrams comprise only non-vanishing contributions from top-quark loops in the ggG 0 vertex and massless quark loops in the ggZ * vertex. Since the contribution of the latter vanishes in Landau gauge for the Z boson, the genuine two-loop calculation is particularly simple in that gauge. It is, however, instructive to inspect the situation in general R ξ gauge for the Z boson as well, a task that is pursued in the Appendix.
The reducible two-loop graphs involve the product of the one-loop induced gg * H and gg * Z vertices and can be calculated with conventional one-loop techniques. These graphs by themselves form a gauge-invariant, UV-finite and IR-finite subset of diagrams.
In practice, we have performed two completely independent calculations of all loop contributions. Version 1 follows basically the strategy described in Ref. [16] for the related process of scalar-pseudoscalar Higgs-boson pair production. Here the Feynman diagrams are generated with the program FeynArts [33], and the large-mass expansion of the diagrams involving top-quark loops is performed by inhouse Mathematica routines. The calculation is carried out in Landau gauge, so that no two-loop diagrams with massless-quark loops contribute. In the second calculation, the diagrams are generated by QGRAF [34] and expanded using EXP/Q2E [35,36]. The massive two-loop vacuum integrals resulting from top-quark loops are calculated using MATAD [37]. This calculation is carried out both in Landau and in unitarity gauge. In the latter, the Goldstone bosons are absent, but massless-quark loops contribute to the ggZ * vertex; these graphs are calculated with the program MINT [38]. While the results of the two different calculations in Landau gauge are in mutual agreement term by term, for the calculation in unitarity gauge only the full result agrees (after the renomalization procedure).
Including all required counterterms, the virtual contribution is given by 1 where M 0 ( ) is the LO amplitude in D = 4 − 2 dimensions, M 1 ( ) the amplitude of the virtual corrections at NLO, and dPS 2 denotes the 2-particle phase-space element. We renormalize the strong coupling in the MS scheme assuming n f = 6 flavours, and the top-quark mass in the on-shell scheme. Note that the renormalization factor is gauge dependent; in Landau gauge, it is where C F = 4 3 , T R = 1 2 , and C A = 3 are the QCD colour factors. In unitarity gauge, the mass counterterm δZ m is absent, but δZ A in Eq. (5) needs to be added according to the dipole subtraction method [39], where n l = 5 denotes the number of light flavours.
In order to be consistent with the currently available PDF sets, we express our final result in terms of α (5) s , the strong coupling with five active flavours using the matching relation [40] α (5) so that the logarithmic dependence on m t vanishes. This procedure is equivalent to decoupling the top quark in the α s renormalization by subtraction of the top-quark loop in the gluon self energy at zero-momentum transfer, instead of using the MS prescription. In the remainder of this paper, we set α s ≡ α s . Inserting the QCD colour factors, the final result for the virtual contribution can be written as where (cf. Eq. (3)) and σ (virt,red) denotes the contribution from the type of reducible diagrams shown in Fig. 2 (k). It is not proportional to σ LO ; as a function of the partonic Mandelstam variablesŝ = (p 1 + p 2 ) 2 ,t = (p 1 − p Z ) 2 , andû = (p 1 − p H ) 2 , it reads (iii) NLO real corrections The real corrections are induced by the partonic channels gg → HZg, gq → HZq, gq → HZq, and qq → HZg, where in the channels involving external quarks only the squares of the diagrams with closed quark loops are taken into account. At first sight, the most complicated one-loop diagrams are pentagon graphs with a top-quark loop. However, in the large-top-mass limit, these graphs vanish. The algebraically most complicated diagrams are the box graphs with external Z * ggg fields; they are the only ones that receive contributions from the vector-coupling of the Z boson, while all other diagrams (summed in pairs of opposite charge flow in the loop) depend only on the Z-boson axial-vector coupling.
The actual calculation of the diagrams can be performed using standard one-loop calculational techniques and has been carried out in three completely independent ways. The first approach builds on graphs from FeynArts 1.0 [33] and reduces or expands the full amplitudes with inhouse Mathematica routines, which produce output in the form of Fortran code. The occurring one-loop tensor and scalar integrals are numerically evaluated with the (not yet public) library Collier that is based on the results of Refs. [41][42][43][44][45][46]. The second calculation is based on the program packages FeynArts 3.2 [47] and Form-Calc/LoopTools [48,49], as far as the calculation of one-loop graphs is concerned that do not involve top-quark loops. The large-mass expansion of the top-quark loops here again is carried out using inhouse routines (independent from the ones of version 1).
The third approach again uses the QGRAF/EXP/Q2E/MATAD-setup for the generation and expansion of the diagrams and the evaluation of the massive vacuum diagrams. The calculation of the massless triangle and box diagrams is performed by an extended version of the FORM [50] routine previously used in Ref. [11], which implements algebraic Passarino-Veltman reduction [42] and analytic results of the scalar integrals given in Ref. [51]. Here we explicitly verified the gauge invariance of the result with respect to the Z propagator as well as to the external gluons. For the latter, we assumed a general axial gauge, where the polarization sum reads with an arbitrary light-like vector n which drops out in the squared amplitude.
All real-emission channels contain IR singularities in their integration over phase space. More precisely, gg → HZg becomes IR singular if the emitted gluon becomes soft or collinear to one of the incoming gluons. The other channels involve only collinear singularities. The separation of the IR singularities in the phase-space integration is achieved using the standard dipole subtraction method [39], where an auxiliary cross section is subtracted from the full real-emission part and added back after an analytical integration (in D dimensions) over the one-particle emission phase that contains the IR singularity (cf. Eq. (7)).

Input values
We use the following input parameters: In the effective-theory approximation, we set m t → ∞ and m b = 0, as discussed above. For the electromagnetic coupling constant α we employ the G µ scheme, where the coupling constant is defined as As the default PDF sets, we use MSTW2008 with n = 0 (n = 1) at LO (NLO). Both the PDFs and the α s evolution are implemented with the help of the LHAPDF library [53]. Our default choice for the renormalization and the factorization scales µ R and µ F is the invariant mass of the HZ system:

Leading-order considerations
In order to get an idea about the quality of the effective theory, we show some studies at LO before presenting our NLO results. Figure 3 shows the partonic cross section both for the exact top-mass dependence and in the effective theory. The exact result exhibits a kink at the top-quark pair threshold √ŝ = 2m t = 344 GeV which clearly cannot be reproduced by the effective-theory approach. For larger values of √ŝ , we do not expect an expansion in 1/m t to converge. In fact, higher-order terms in this expansion would most likely worsen the prediction in the region of larger √ŝ .   The situation becomes more problematic in the boosted regime which we study by imposing a lower cut on the Higgs' transverse momentum, requiring p T,H > 200 GeV, see Fig. 4 (b). In this case, the minimal value for √ŝ is already above the top-quark threshold when M H = 100 GeV. Consequently, the direct application of the effective-theory approximation is off by almost a factor of five to ten, which is clearly unacceptable.
A direct evaluation of the NLO contribution in the effective theory is therefore not possible. However, in Refs. [21][22][23][24][25][26][27] it was shown for the process gg → H at NLO and NNLO that the perturbative correction factor, defined at NLO in Eq. (2), depends only very weakly on the top-quark mass. To some degree, this holds even far outside the convergence region of the heavy-top expansion, as long as only the leading term in 1/m t is taken into account. Motivated by this observation, we move on to NLO and present our results in the next

Correction factor
As outlined above, we evaluate the NLO hadronic cross section by rescaling the full LO result by the perturbative K-factor calculated in the effective theory: Since we are aiming at a NLO quantity, it actually might be more appropriate to evaluate the formally LO cross sections in Eq. (17) with NLO PDFs. We checked that the effect of this is much smaller than the uncertainty due to variations of the renormalization and factorization scale, which is why we stick to LO PDFs in σ LO .  Figure 6: Individual contributions to the NLO hadronic K-factor at √ s = 14 TeV as discribed in the main text. Figure 5 shows the gluon-induced cross section obtained in this way for √ s = 8 TeV and √ s = 14 TeV hadronic centre-of-mass energy, together with the corresponding perturbative correction factor K. Part (a) of Fig. 5 shows the total inclusive cross section, while in part (b) the boosted scenario with p T,H > 200 GeV is shown. In both cases, we observe a K-factor of the order of two, almost independent of M H , with a slight increase towards lower centre-of-mass energies. This behaviour is very similar to the one observed for gluon-induced single-Higgs [21][22][23][24][25][26][27] and Higgs pair production [16]. The correction even slightly exceeds the well-known correction factor for gg → H.
A breakdown into individual contributions to K is shown in Fig. 6 for the total inclusive cross section at √ s = 14 TeV: • K LO -change of PDF sets from LO to NLO • ∆K ggvirt -virtual corrections including integrated dipole terms according to Ref. [39] • ∆K ggreal -correspondingly regularized real corrections  • ∆K qg , ∆K qq -contributions from qg and qq initial states The sum of all these terms results in K tot , the total K-factor.

Residual scale uncertainty
As described in the introduction, the LO scale dependence for this purely gluon-induced process is quite large. NLO corrections typically decrease this uncertainty. Let us recall the situation in the gluon-fusion process gg → H, however: For the LO result, the usually adopted scale variation by a factor of two around the central scale leads to a gross underestimation of the size of the higher-order effects. At NLO, the scale uncertainty is not significantly smaller than at LO, but it does provide a good estimate of the NNLO effects. Consistently, inclusion of the NNLO corrections leads to a significant reduction of the scale uncertainty.
Expecting a similar behaviour for the gg → HZ process, it is not surprising to see the result shown in Fig. 7: Both for the inclusive and the boosted scenario the scale dependence decreases from more than 100% at LO to 60% at NLO when the renormalization and factorization scales are varied simultaneously by a factor of six around their central value µ 0 , see Eq. (16). As for the process gg → H, the behaviour in µ/µ 0 is strictly monotonous, and the LO and NLO curves do not intersect. Therefore, a preferred value for µ F and µ R cannot be deduced from these plots. The radiative corrections increase with µ/µ 0 , so there is a slight tendency towards choosing smaller values of µ. Nevertheless, in our numerical analysis we stick to the "natural" value µ 0 as the central choice.
Note that, also similar to what is observed in gg → H, variation by a factor of two would not lead to any overlap between the LO and the NLO predictions.
The similarity between the processes gg → H and gg → HZ suggests that the NLO error estimate due to scale variation is quite reliable for the process gg → HZ. In order to take into account the fact that the effective theory is expected to work not quite as well in gg → HZ as in gg → H, we determine this uncertainty by varying µ within a factor of three rather than two around the central value µ 0 . The numerical results are listed in Table 1.

Total inclusive cross section
In this section we provide the most up-to-date numbers for the total inclusive cross section for the Higgs-strahlung process at the LHC with 8 and 14 TeV, including • NNLO Drell-Yan terms σ HV,DY of order g 4 α n s (n = 0, 1, 2) [7-9]; • electroweak corrections which are applied as an overall factor to the Drell-Yan terms [12,14]; • top-loop-induced corrections of order O(λ t g 3 α 2 s ) [11]; • gluon-induced terms of order λ 2 t g 2 α n s (n = 2, 3); n = 3 corresponds to the newly calculated terms of this paper.
For the non-gluon-fusion part of the cross section, the scale variation is obtained by using the MSTW2008NNLO PDF set and varying µ F and µ R independently within the interval (µ R , µ F )/M H ∈ [1/3, 3] × [1 /3, 3], which results in a cross section interval [σ The central value of the total cross section is then obtained as where σ gg are the boundaries of the scale uncertainty interval of the gluon-induced component which can be obtained from Table 1. Accordingly, the scale uncertainty is calculated as   Table 2: Total inclusive cross section for the processes pp → HW and pp → HZ. The latter includes the newly calculated NLO gluon-induced terms. The evaluation of the scale and PDF+α s uncertainties is described in the main text.
The influence of the newly evaluated NLO gluon-induced terms on the overall PDF+α s uncertainty is rather small, since this contribution comprises only about 5% of the total cross section, and will be neglected. Therefore, we base the estimate of the PDF+α s uncertainty solely on what is currently contained in vh@nnlo (i.e., LO gluon-induced terms are taken into account). Following the PDF4LHC recommendations [54] by using the NNLO PDF sets from MSTW2008 [52], CT10 [55], and NNPDF23 [56], we obtain a cross section interval [σ PDF+αs ] and calculate the resulting uncertainty as Our results are shown in Table 2. We find that the NLO gluon-induced terms calculated in this paper increase the central values of the HZ cross section by about 4% (7%) at 8 TeV (14 TeV). Since the K-factor for these terms is of the order of two, and their scale uncertainty decreases by roughly the same factor when going from LO to NLO (see Table 1), the overall scale uncertainty on the total inclusive cross section remains almost unaffected by the inclusion of the new terms. For completeness, we also include updated numbers for HW production in Table 2, even though they are not affected by the NLO gluon-fusion terms calculated in this paper.
As a side remark, we note that the PDF+α s uncertainties of Table 2 are significantly smaller than in Ref. [11]. This is due to the use of only NNLO PDF sets in this newer version, while the previous numbers were based on a rescaling of the NNLO MSTW2008 uncertainty by the NLO PDF error. For HW production, also the scale uncertainty is slightly smaller in Table 2 than in Ref. [11]. This is because these previous numbers were obtained by linearly adding uncertainties of the "top-induced" terms to the rest, while here we vary the scale in both contributions simultaneously. For the HZ cross section, this is overcompensated by the uncertainty of the NLO gluon-fusion component, see above.

Conclusions
The gluon-induced corrections to the Higgs-strahlung process have been calculated through NLO, i.e. O(α 3 s ). The perturbative correction factor is obtained in the limit m t → ∞ and m b = 0, and then used to rescale the full LO cross section in order to obtain a prediction for the NLO result. The behaviour of the perturbative series and the residual scale variation are found to be comparable to the gluon-fusion processes of single-Higgs and Higgs pair production. The success of the effective-field-theory approach in describing corrections to single-Higgs productions, where exact calculations for higher-order corrections are available, gives confidence in a reliable prediction of the theoretical uncertainty due to higher-order effects for gluon-induced HZ production.
Numerical results were provided for current and future LHC energies, both for the fully inclusive cross section and for boosted Higgs kinematics. We use these results in order to provide the most up-to-date predictions for the hadronic Higgs-strahlung process at relevant collider energies. The large corrections on the gluon-induced terms, combined with their large scale uncertainty, increases the overall uncertainty on the total HZ cross section by about 1%.
In the near future, the results will be included in the publically available numerical program vh@nnlo.
Inserting the latter into Eq. (21) and exploiting the fact that the transversal part of D ρσ ξ (k) (first term in Eq. (22)) does not contribute, we obtain where we have used k = p Z + p H , p Z · ε * Z = 0, and k 2 =ŝ. The terms in square brackets can be simplified using the well-known form of the ABJ anomaly.
To this end, we recall that the ABJ anomaly relation for the axial current j µ f,5 = ψ f γ µ γ 5 ψ f and the pseudo-scalar operator p f,5 = ψ f γ 5 ψ f for a quark q reads ∂ µ j µ q,5 = 2im q p q,5 − α s 4π F a µνF a,µν , whereF a,µν = 1 2 µνρσ F a ρσ is the dual of the gluonic field-strength tensor F a,µν . Eq. (24) is an operator relation, valid for bare quantities, expressing chiral symmetry. The Adler-Bardeen theorem [58] states that it is correct to all orders in regularization schemes that respect chiral symmetry. In regularizations that are not chirally symmetric, Eq. (24) has to be restored by extra counterterms from evanescent operators. For the 't Hooft-Veltman γ 5 scheme [30,31], which is employed in our calculation, these are the counterterms δZ A 5 and δZ P 5 calculated in Ref. [32] and already used in Section 3. The relevant momentum-space Feynman rules for the composite operators j µ q,5 , p q,5 , and FF are given by The dotted line in the Feynman rules indicate that the momentum p flows into the vertex. The operator FF induce Feynman rules with three and four gluon lines as well, but those will not contribute in the following. The fermionic Feynman rules are related to the couplings of Z/G 0 to the quark: 2s w c w γ µ γ 5 + vector part, G 0q q : − eI 3 w,q where I 3 w,q is the third component of the weak isospin of q.
Now we can make contact with the Green functions G g a g b Z and G g a g b G 0 introduced above. Since we consider only QCD corrections, in the relevant graphs the couplings of the external Z/G 0 lines to an internal quark line represent the only electroweak coupling in the ggZ/G 0 vertex functions. These electroweak couplings can be interpreted as the insertions of the operators j µ q,5 and p q,5 , because the vector part of the Z coupling does not contribute, as explained in Section 2. Thus, we obtain FF FF FF FF FF where the dependence on the gauge parameter ξ cancels, as it should be. For the LO matrix element, the Green function G g a g b (FF ) just has to be replaced by the Feynman rule for the FF operator with two external gluon legs, yielding in agreement with Eq. (3). The NLO QCD corrections to G g a g b (FF ) are induced by the diagrams shown in Fig. 8. The actual calculation of these one-loop diagrams is very simple.
where C 0 is the one-loop scalar 3-point integral in the convention of Refs. [44][45][46] (see e.g. Ref. [60] for the explicit result). This is the unrenormalized result for the virtual correction which still has to be renormalized. As stated above, the anomaly equation (24) is valid for bare quantities, so that the FF term receives renormalization contributions from α s and the gluon field. Since in our case only the two-gluon contribution of FF is relevant, we get the simple factorizing contribution to the NLO amplitude, Combining the renormalized virtual amplitude M 1PI,virt 1 + M 1PI,ct 1 with the contribution δ CS , see Eq. (7), of the Catani-Seymour I-operator of the subtraction function yields the 1PI part of the correction to the cross section, as given in Eq. (9) by first term on the r.h.s.