$ZH$ production in gluon fusion: two-loop amplitudes with full top quark mass dependence

We present results for the two-loop helicity amplitudes entering the NLO QCD corrections to the production of a Higgs boson in association with a $Z$-boson in gluon fusion. The two-loop integrals, involving massive top quarks, are calculated numerically. Results for the interference of the finite part of the two-loop amplitudes with the Born amplitude are shown as a function of the two kinematic invariants on which the amplitudes depend.


Introduction
The production of a Higgs boson in association with a vector boson, also called Higgs-Strahlung, is an important process at the LHC. For example, this process was used to discover the decay of Higgs bosons to b-quark pairs [1,2], and it is very well suited to constrain anomalous couplings of the Higgs boson in both the Yukawa and the gauge boson sector. The production of a Higgs boson in association with a leptonically decaying Z-boson facilitates triggering independently of the Higgs boson decay. This makes this channel especially attractive in combination with challenging Higgs decays, like invisible or hadronic decays, in particular H → bb. The regime of boosted Higgs bosons is particularly interesting, as the large-p T,H region is sensitive to new physics [3][4][5][6][7][8][9]. The loop-induced gluon-initiated contributions, calculated at LO in Refs. [10,11], are finite and enter at order α 2 s , i.e. formally at NNLO considering the pp → V H process. However, due to the dominance of the gluon PDFs at the LHC they are sizeable, contributing about 6% of the total NNLO cross section, and the contribution can be twice as large in the boosted Higgs boson regime p H T 150 GeV [12,13]. As only the LO results for gg → ZH are known and implemented in current Monte Carlo generators, the large scale uncertainties pertaining to this channel constitute a significant part of the uncertainties in experimental analyses of the V H process [1,2,14,15]. Furthermore, the gluon-initiated subprocess is very sensitive to modified Yukawa couplings and/or non-SM particles running in the loop, For these reasons the NLO corrections to this process are very important. However, the NLO corrections contain two-loop integrals involving m t , m H and m Z , and such integrals are currently unknown analytically. Therefore the QCD corrections have been calculated in various approximations so far. In Ref. [16] they have been calculated in the m t → ∞ limit, resulting in a K-factor used to reweight the full one-loop result. In addition, top quark mass effects at NLO QCD have been considered in the framework of a 1/m t -expansion in Ref. [17], including Padé approximants constructed from expansion terms up to 1/m 8 t . However, the 1/m t -expansion becomes invalid for invariant masses of the HZ system larger than 2m t , which is the interesting region with regards to new physics searches. Soft gluon resummation for the gg → ZH process has been calculated in Ref. [18]. In Ref. [19], a data-driven strategy to extract the gluon-initiated component (or, more precisely, the non-Drell-Yan component) for ZH production has been suggested, based on the comparison of the ZH to the W H cross section and the corresponding invariant mass distributions of the V H system. Considering the full process pp → ZH, inclusive NNLO QCD corrections have been available for quite some time [20] and are implemented in the program VH@NNLO [20][21][22][23]. NLO electroweak corrections have been calculated in Refs. [24,25], combined NLO QCD+EW corrections are also available [26][27][28]. Differential QCD corrections to the process pp → ZH have been calculated up to NNLO, including H → bb decays at NLO [29,30] and at NNLO [31]. In Ref. [32], fully differential NNLO results for V H observables including the decays of the vector boson into leptons and the Higgs boson into b-quarks with off-shell propagators of the vectorand Higgs bosons have been calculated. The combination of fixed-order QCD computations with parton showers has been studied at NLO+PS in association with up to one jet [27,[33][34][35] and at NNLOPS [36,37]. In addition, the NNLO corrections have been combined with NNLL resummation in the 0-jettiness variable and matched to a parton shower within the Geneva Monte Carlo framework [38]. Threshold corrections up to N 3 LL for the Drell-Yan-type part of the inclusive cross section have been calculated in Ref. [39], soft-gluon resummation of both threshold logarithms and logarithms which are important at low transverse momentum of the V H pair have been considered up to N 3 LL in Ref. [40] and have been found to be very close to the fixed order NNLO result. The process bb → ZH in the five-flavour scheme, but with a non-vanishing bottom-quark Yukawa coupling, has been calculated in the soft-virtual approximation at NNLO QCD in Ref. [41], the polarised qq → ZH two-loop amplitudes have been calculated in Ref. [42]. In this paper, we calculate the two-loop virtual corrections to the process gg → ZH, including massive top quarks in the loops. We focus on the description of the calculation and display numerical results for the two-loop amplitudes. A phenomenological study based on these results is postponed to a subsequent publication. The structure of this work is as follows. In Section 2 we describe details of the calculation such as the integral families, our projection operators, the reduction to master integrals and the evaluation of the master integrals. In Section 3 we show results for the finite part of the virtual amplitude, for some benchmark points as well as in terms of two-dimensional grids, before we conclude in Section 4.

Details of the calculation 2.1 Amplitude definition
We consider the process with p 2 1 = p 2 2 = 0, p 2 3 = m 2 Z and p 2 4 = m 2 H . The Mandelstam invariants are defined by The diagrams contributing to the process (2.1) at leading order are shown in Fig. 1. In our calculation we use the 't Hooft-Feynman gauge, therefore the diagram involving the exchange of a would-be Goldstone boson G 0 needs to be taken into account. We treat all quarks except the top-quark as massless, therefore only top-quark loops contribute in those diagrams where the Higgs boson couples directly to the fermion loop. The effect of a finite bottom quark mass on the total LO cross section is at the per mille level [17]. In the triangle diagrams, only the axial vector part of the Z-boson coupling contributes due to Furry's theorem. In addition, the massless quark contributions cancel in each isospin doublet, such that only the third generation of quarks give a non-zero contribution for this class of diagram. The leading order amplitude can be expressed in terms of seven form factors [11] containing one-loop three-and four-point functions. Some of the form factors can be related by crossing p 1 ↔ p 2 such that only four form factors remain to be calculated. However we choose not to express our amplitude in terms of these form factors, as will be explained in Section 2.1.1.
In the m t → ∞ limit the amplitude simplifies considerably and is given by [16] where is the totally antisymmetric Levi-Civita symbol with (a, b, c, d) ≡ µνρσ a µ b ν c ρ d σ and ε 1 , ε 2 are the polarisation vectors of the incoming gluons, carrying colour indices a and b, while ε * Z is the polarisation vector of the outgoing Z-boson. The spin-and colour-averaged Born matrix element then can be written as where λ denotes the Källén function, λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc. We calculate using conventional dimensional regularisation (CDR) with D = 4 − 2 .
For the treatment of γ 5 within dimensional regularisation we use the 't Hooft-Veltman scheme [43,44] in the variant of Refs. [45,46], i.e. we use and add finite renormalisation terms to restore chiral symmetry in the massless quark limit, see Section 2.2.3. The contraction of two -symbols leads to linear combinations of metric tensors which are treated as D-dimensional.

Tensor structures and projection to a basis of linear polarisations
We define the tensor amplitude A µ 1 µ 2 µ 3 by extracting the polarisation vectors from the amplitude of process (2.1), where the ε µ i λ i denote the polarisation vectors. The Lorentz tensor structures appearing in the amplitude A µ 1 µ 2 µ 3 were discussed in Ref. [11]. However, we do not use form factors related to these Lorentz structures here, but rather use projections based on the momentum-basis representations of the linear polarisation vectors of external particles as suggested in Ref. [47]. We also use the fact that only one axial current is involved in the QCD corrections to this amplitude, therefore all relevant Lorentz structures contain only a single Levi-Civita symbol. In addition, conditions such as transversality and Bose symmetry regarding the two external gluons further constrain the possible Lorentz structures. Following the procedure of Ref. [47], we define the following normalised linear polarisation vectors, where the frame that has been used is shown in Fig. 2, The normalisation factors N i for i ∈ {x, y, T, l}, associated with each of these spacelike polarisation vectors in Eq. (2.7), can be determined from their (negative) norm squares ε 2 i , i.e. N i = 1/ −ε 2 i , which we choose to include only at the very end of the calculation. The physical meaning of these vectors is apparent in the center-of-mass frame of p 1 and p 2 , where the beam axis determined by { p 1 , p 2 } is taken as the z-axis and the plane determined by { p 1 , p 3 } defines the x-O-z plane, as illustrated in Fig. 2. Then the polarisation vector ε x is orthogonal to the beam axis and lies within the x-O-z plane, while ε y is perpendicular to this plane. The vector ε T lies within the x-O-z plane but points to a direction orthogonal to p 3 , and ε l , the longitudinal polarisation of the Z-boson, points along the direction of p 3 in the center-of-mass frame. Our projectors are given by products of three linear polarisation vectors, which are further re-written solely in terms of external momenta and the Levi-Civita symbol, all with D-dimensional Lorentz indices. However, only the six projectors with an odd number of Levi-Civita symbol, i.e. an odd number of ε y , are relevant for the process (2.1). This is because there are only three linearly independent external momenta in A µ 1 µ 2 µ 3 and all Lorentz structures appearing in the amplitude contain a Levi-Civita symbol. Concretely, the six projectors we use for A µ 1 µ 2 µ 3 are given by The projections made with P µ 1 µ 2 µ 3 4 and P µ 1 µ 2 µ 3 5 are related to those with P µ 1 µ 2 µ 3 2 and P µ 1 µ 2 µ 3 3 , respectively, through crossing of p 1 and p 2 . Therefore we only need to compute four projections in total. Note that all open Lorentz indices in Eq. (2.9) are understood to be D-dimensional, and pairs of Levi-Civita symbols in P µ 1 µ 2 µ 3 6 need to be contracted before being used for the projection of the amplitude (see Ref. [47] for a more detailed discussion). From the point of view of a form factor decomposition, the set of linear-polarisation projectors in Eq. (2.9) represent precisely the Lorentz tensor decomposition basis in use. This is a consequence of the orthogonality of the projectors (although the linear completeness is in general only ensured in the 4-dimensional limit [41,47]). This projection defines six quantities (2.10) The physical interpretation of the quantities A n as linearly polarised amplitudes offers us a convenient short-cut to transform to a helicity basis defined w.r.t the same reference frame. The usual helicity amplitudes for the process (2.1) can be constructed from the linear ones, using the relations and the longitudinal polarisation of the Z-boson is given by ε µ 3 l in Eq. (2.7).

Kinematics
The kinematic invariants in terms of the scattering angle θ between the beam axis and p 3 in the centre-of-mass frame are given by (2.13) The limit β ZH → 1 corresponds to s m 2 + , m 2 − . In the forward scattering region θ → 0, if in addition β ZH → 1 is fulfilled, we have |t| |u|. Analogously, for θ → π and β ZH → 1, the ratio t/u is very small. The virtual amplitude has a threshold at s = 4m 2 t , therefore it is also useful to define at the top quark pair production threshold s = 4m 2 t and β t → 1 for s → ∞. For 2 → 2 kinematics, the transverse momentum p T of the Higgs-or Z-boson (with Therefore, in the high energy limit, p 2 T → s (1 − cos 2 θ).

Diagram Generation, reduction and calculation of the master integrals
We generate the 1-and 2-loop diagrams using QGRAF [48]. The projectors described in Section 2.1.1 are applied and the Feynman rules are inserted using FORM [49][50][51]. For the reduction we have defined eight integral families F i , five planar (i = 1, . . . , 5) and three non-planar (i = 6, . . . , 8) families. We also use five additional families obtained from the original families by exchanging p 1 and p 2 . Each family contains nine propagators which allows all irreducible scalar products in the numerator to be expressed in terms of inverse propagators. Concretely, the occurring integrals have the form 15) where the N j denote genuine propagators of the form 1/(k 2 −m 2 ) with exponents r i ≥ 1 and s i ≥ 0. The families are listed in Tables 1 and 2. Integral families which differ from the listed ones by exchanging p 1 and p 2 are not shown. The integrals appearing in the amplitude are reduced to a minimal set of master integrals as described in the following.

Choice of master integrals and reduction procedure
The numerical evaluation of finite integrals is typically much simpler than that of their divergent counterparts, we therefore follow the strategy of Ref. [52] to obtain a quasifinite basis of master integrals. This choice of master integrals is not unique and the size of the coefficients of the integrals after reduction depends on this choice. In particular, it is known that it can be advantageous to choose a basis where the denominators in the reduction tables factorizes into factors containing only the space-time dimension D and factors depending on the kinematic invariants only. While algorithms to find such a basis have been presented in Refs. [53,54], we obtained such a basis following a different approach, namely, by iterating over different combinations of master integral candidates and analysing the resulting reduction tables, restricted to a small subset of the full IBP system and neglecting subsectors. This allows us to define additional criteria for the selection of the preferred masters, such as the size of the appearing denominator factors or the order in = (4 − D)/2 at which an integral starts contributing to the amplitude. With our choice of master integrals, the D-dependence of the denominator factors of the reduction rules factors for all integrals, except for some one-particle reducible integrals. Furthermore, all seven-propagator integrals only start contributing to the amplitude at order −1 and the size of the amplitude coefficients reduces by about a factor of 5 compared to a default choice of finite masters with minimal propagator powers. In Table 1: Planar integral families used for the reduction.
In order to perform the reduction to our chosen set of master integrals, we utilize the Kira package [55,56] in combination with the rational function interpolation library FireFly [57,58]. A crucial benefit of finite-field interpolation techniques is that one circumvents large intermediate expressions even during the calculation of IBP tables, which leads to a significant runtime and memory reduction in general. As a consequence, we interpolate the required reduction tables that relate the integrals of the amplitudes that contribute up to r = 7 and s = 4 to the desired set of master integrals that requires IBP relations up to r = 9 and s = 2 in D dimensions. To simplify the reduction of the integrals occurring in the amplitude, we scale the Z and H mass w.r.t. m t and approximate m 2 Z /m 2 t = 23/83 and m 2 H /m 2 t = 12/23, corresponding to m t = 173.21 GeV, m Z = 91.18 GeV, and m H = 125.1 GeV. Note that although we retain the dependence on m t , the ratio of the Z-boson mass and Higgs boson mass to the top quark mass is fixed in our calculation. By setting m t = 1 we remove an additional scale that can be restored by dimensional analysis, which leaves us with a 3-parameter problem. As we use quasi-finite master integrals in six dimensions, the set of dimensional recurrence relations (DRRs), which connects integrals in six and four dimensions, has to be calculated as well. Hence, we utilize LiteRed [59] and Reduze [60] in order to obtain these relations. The DRRs are subsequently related to our set of master integrals using Kira with FireFly in the same setup as described above. We note that the calculation of the relations between the DRRs and our basis of master integrals was the most demanding step in the whole calculation. It took about four days of wall-clock time running on a machine with two Intel ® Xeon ® Silver 4116 and hyper-threading. The required memory never exceeded about 100 GB of memory. Afterwards the DRRs and the reduction tables obtained by reducing the integrals of the amplitudes are combined to a custom system of equations that fills roughly 9 GB of disk space. This system is again solved by employing Kira with FireFly in order to interpolate the final set of replacement rules. Finally, the latter set of replacements is inserted into the amplitudes with the help of the ff insert executable of FireFly. It is worth mentioning that our calculation, which was split into several substeps, could have been performed in a single run by using the Kira option iterative reduction: masterwise with a custom system of equations that also includes the amplitudes. However, as one needs to hold the whole system of all integral families in memory at once in this case, we observed faster runtime and lower memory consumption by splitting the calculation as described above and running all steps on different machines in parallel. Additionally, splitting the reduction into several steps is convenient when studying different bases of master integrals and their impact of the total file size of the resulting amplitudes.

Evaluation of the two-loop amplitude
To evaluate the master integrals, we first apply sector decomposition as implemented in the program pySecDec [61,62]. For integrals which diverge in the limit of 4 spacetime dimensions ( → 0), sector decomposition resolves singularities in the regulator leaving only finite integrals over the Feynman parameters which can then be integrated numerically. In the physical region, singularities can appear on the real axis of the Feynman parameters and a causal prescription to avoid the singularities is required. We deform the integration contour of the Feynman parameters into the complex plane as described in Refs. [63][64][65][66].
In the present calculation we have selected a basis of quasi-finite integrals. This means that poles in can appear only in the prefactor of our integrals after Feynman parametrisation and, thus, sector decomposition is not required to resolve singularities in the regulator. Nevertheless, we observe that applying a sector decomposition greatly simplifies the numerical evaluation of the finite integrals. This observation can be partly understood by noting that sector decomposition also removes integrable singularities appearing at the boundaries where one or more Feynman parameters vanish. We therefore process all integrals with pySecDec before numerical integration. The numerical integration itself is performed using the quasi-Monte Carlo (QMC) algorithm [62,67]. For all integrals we apply a Korobov periodising transform with weight 3. We observe that the use of the rank-1 shifted lattice rules greatly reduces the number of samples required to obtain the integrals to sufficient precision for the computation of our amplitude, as compared to straightforward Monte Carlo sampling. In line with our previous work on the processes pp → HH [68,69], pp → HJ [70] and gg → γγ [71], we extend the pySecDec program such that it can produce a code capable of evaluating the entire amplitude, rather than computing each integral separately. The advantage of this structure is that the number of sampling points used for each sector of each integral can be dynamically set according to its contribution to the total uncertainty on the amplitude. We utilise a variant of the procedure described in Ref. [69] to minimise the total time taken to obtain a given relative accuracy on the amplitude.

Renormalisation
We may expand each of the form factors, A i=1,...,n , in the bare strong coupling a 0 = α 0 /(4π) according to We renormalise the strong coupling in the MS scheme with the heavy quark loop in the gluon self-energy subtracted at zero momentum. The heavy quark mass is renormalised in the OS scheme. The UV renormalised amplitude is obtained via the replacement Here, n g = 2 is the number of external gluon legs, a s = α s (µ 2 R )/(4π) S −1 (µ 2 R /µ 2 0 ) with S = (4π) e − γ E , and the renormalisation constants are expanded according to We handle the γ 5 matrices which appear in this calculation using the Larin scheme [46]. According to this scheme an additional finite renormalisation is required, it is denoted in the equations above by δZ 5 . At the level of individual form factors, the procedure outlined corresponds to the following relations

Definition of the finite part of the virtual two-loop amplitude
In order to obtain IR finite amplitudes we use the subtraction scheme described in Ref. [72] We present results at renormalisation scale µ 2 R = s and after changing to the helicity basis, we define for the square of the amplitude. The electroweak coupling α is set to 1.

Checks of the calculation
We have verified that our results have the expected universal pole structure and crossing symmetries. We also compared our exact result with the expansion in the large top quark mass limit presented in [17]. The virtual contributions V n are expanded up to order 1/m 2n t and reweighted with the exact Born term B. The one-particle-reducible double triangle contributions V red are included with full top quark mass dependence.
In Fig. 3 the ratio of expanded to full virtual contribution V n /V is shown for expansion orders 0 ≤ n ≤ 4 for a fixed scattering angle cos(θ) ≈ 0.052. For energies close to the production threshold at β t = −1 (s = (m Z + m H ) 2 ) the expanded result approximates the exact calculation well with a ratio V 4 /V ≈ 0.9989, while the agreement worsens closer to the top quark pair threshold at β t = 0 (s = 4m 2 t ), where the large m t -expansion is expected to break down. In addition, we have compared our results to a recent calculation in the high energy limit [73] and found agreement to the extent expected from previous comparisons for the process pp → HH performed in Ref. [74].  [17] for fixed scattering angle cos(θ) ≈ 0.052. The range of the horizontal axis −1 ≤ β t ≤ 0 corresponds to the energy range (m Z + m H ) 2 ≤ s ≤ 4m 2 t between production threshold and top quark pair threshold.

Numerical results for the two-loop amplitudes
For the presentation of our results, we have evaluated a total of 460 phase-space points at 2-loop. We request per mille precision for each of the linearly polarised amplitudes, this is obtained for most phase-space points with between 45 minutes and 24 hours of run time using 2 x Nvidia Tesla V100 Graphics processing units (GPUs). In Fig. 4, we show results for the unpolarised modulus of the Born amplitude, as well as for the finite part of the virtual two-loop amplitude interfered with the Born amplitude. We see that the unpolarised Born result shows a rather flat dependence on the scattering angle around and below the top quark pair threshold, due to the dominating s-wave contribution in this region, and starts to curve as energy further increases because partial waves with higher angular momentum play an increasingly important role. We also observe that in the two-loop case, the top quark pair production threshold region is much more peaked than at leading order, due to log β t -terms appearing for the first time at two-loop order. In Fig. 5, the ratio of the two-loop amplitude to the Born-amplitude is shown separately. In the following we will show results for five independent helicity amplitudes in the β t -cos θ-plane. The remaining amplitudes can be obtained from overall helicity flips as well as crossing between the two initial gluons. Bose-and CP-symmetry (no CKMmatrix is involved in our calculation) imply [75,76] A λ 1 ,λ 2 ,λz (t, u) = (−1) λz A λ 2 ,λ 1 ,λz (u, t) , To understand the qualitative behaviour, we can expand the amplitude in terms of Wigner d-functions, where s i and s f denote the initial and final state total spins, respectively, and J denotes the total angular momentum of the system. For gg → ZH, the initial state has total spin s i = 0 for equal helicities of the initial state gluons, (λ 1 , λ 2 ) = (+, +) or (−, −), while for the case (λ 1 , λ 2 ) = (+, −) or (−, +), the initial state has total spin s i = 2. Therefore the amplitude A ++0 is dominated by the partial wave d 0 00 (θ) and provides the largest contribution to the total squared amplitude. In particular in the low energy region, where the ZH system has relatively small kinetic energy, the s-wave contribution should dominate, reflected in the homogeneity of A ++0 in cos θ. As the center-of-mass energy increases, the contributions from partial waves with higher angular momenta also start to play a role, leading to a non-flat behaviour in cos θ. From Fig. 6, we further observe that the helicity amplitudes with the polarisations Z ± are suppressed compared to those with a longitudinally polarised Z-boson. The amplitudes A ++± are antisymmetric under exchange of t and u, so antisymmetric in cos θ, and therefore have no s-wave contribution. The d-wave contributions d 2 0,±1 (θ) are proportional to ± cos θ sin θ, which vanish at cos θ = 0, ±1, a behaviour that can be observed in A ++− . However, kinematics encoded in the coefficients of the partial waves also plays a major role, such that the shapes cannot be explained by partial waves alone. Note that A ++− is about five orders of magnitude smaller than A ++0 , and A +++ is also suppressed, therefore the amplitudes A ++± give a very minor contribution in the sum of all polarisation configurations. For different helicities of the initial state gluons, (λ 1 , λ 2 ) = (+, −) or (−, +), the initial state has total spin s i = 2, such that the partial wave contributions start from J = 2. Therefore the amplitude A +−0 is much smaller than A ++0 . This is also reflected in the fact that the amplitudes are basically zero except at very high energies. The amplitude A +−0 has no contribution from d J 2,0 (θ) with even J as it is antisymmetric in cos θ [76], therefore its leading partial wave is given by d 3 2,0 (θ) ∼ cos θ sin 2 θ, vanishing at cos θ = ±1 and cos θ = 0. The amplitudes A +−± seem to have their leading partial waves given by d 2 2,±1 (θ) ∼ sin θ(1 ± cos θ). Consequently A +−+ is highly suppressed in the backward direction, while A +−− is highly suppressed in the forward direction. In Table 3 we list numerical results of the amplitude to facilitate comparisons with our results.

Conclusions and Outlook
We have numerically calculated the two-loop amplitudes for the production of a Higgsand a Z-boson in gluon fusion with massive top quark loops. The results for the finite part of the two-loop amplitude interfered with the Born amplitude are plotted as a function of the scattering angle and the centre-of-mass energy, for the total unpolarised amplitude as well as for individual helicity amplitudes. The projection of the amplitudes to scalar quantities has been carried out with projectors onto linear polarisation states, from which helicity amplitudes are constructed [47]. The reduction to master integrals has been performed with the program Kira [55,56] in combination with the rational function interpolation library FireFly [57,58], using in addition LiteRed [59] and Reduze [60] to obtain dimensional recurrence relations. The master integrals have been calculated using pySecDec [61,62]. The integration is sufficiently stable and accurate, also in the near-threshold and forward scattering regions, to perform phenomenological studies based on these results, after including the real radiation contributions. We postpone such a phenomenological analysis to a subsequent publication. Our method for the first time has been applied to a process with three different mass scales, m t , m H and m Z . We expect that it can be applied successfully to other two-loop amplitudes involving several mass scales in the future.