ZH production in gluon fusion: two-loop amplitudes with full top quark mass dependence
Long Chen, Gudrun Heinrich, Stephen P. Jones, Matthias Kerner, Jonas Klappert, Johannes Schlenk
PPrepared for submission to JHEP
ZU-TH 45/20CERN-TH-2020-199IPPP/20/57P3H-20-076KA-TP-21-2020TTK-20-42PSI-PR-20-21
Z H production in gluon fusion: two-loop amplitudes with full top quark massdependence
Long Chen, a,f
Gudrun Heinrich, b Stephen P. Jones, c,d
Matthias Kerner, e JonasKlappert, f Johannes Schlenk ga Max Planck Institute for Physics, F¨ohringer Ring 6, 80805 M¨unchen, Germany b Institute for Theoretical Physics, Karlsruhe Institute of Technology, 76128 Karlsruhe, Ger-many c Theoretical Physics Department, CERN, Geneva, Switzerland d Institute for Particle Physics Phenomenology, Durham University, Durham DH1 3LE, UK e Physik-Institut, Universit¨at Z¨urich, Winterthurerstrasse 190, 8057 Z¨urich, Switzerland f Institute for Theoretical Particle Physics and Cosmology, RWTH Aachen University, 52056Aachen, Germany g Theory Group LTP, Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland
E-mail: [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] Abstract:
We present results for the two-loop helicity amplitudes entering the NLOQCD corrections to the production of a Higgs boson in association with a Z -bosonin gluon fusion. The two-loop integrals, involving massive top quarks, are calculatednumerically. Results for the interference of the finite part of the two-loop amplitudeswith the Born amplitude are shown as a function of the two kinematic invariants onwhich the amplitudes depend. Keywords:
LHC, QCD phenomenology, two-loop computations, Higgs boson, vectorbosons, top quark mass a r X i v : . [ h e p - ph ] D ec ontents 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 todiscover the decay of Higgs bosons to b -quark pairs [1, 2], and it is very well suitedto constrain anomalous couplings of the Higgs boson in both the Yukawa and thegauge boson sector. The production of a Higgs boson in association with a leptonicallydecaying Z -boson facilitates triggering independently of the Higgs boson decay. Thismakes this channel especially attractive in combination with challenging Higgs decays,like invisible or hadronic decays, in particular H → b ¯ b . The regime of boosted Higgsbosons is particularly interesting, as the large- p T,H region is sensitive to new physics [3–9].The loop-induced gluon-initiated contributions, calculated at LO in Refs. [10, 11], arefinite and enter at order α 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,– 1 –ontributing about 6% of the total NNLO cross section, and the contribution canbe twice as large in the boosted Higgs boson regime p HT (cid:38)
150 GeV [12, 13]. Asonly the LO results for gg → ZH are known and implemented in current MonteCarlo generators, the large scale uncertainties pertaining to this channel constitute asignificant 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 modifiedYukawa couplings and/or non–SM particles running in the loop, For these reasons theNLO corrections to this process are very important. However, the NLO correctionscontain two-loop integrals involving m t , m H and m Z , and such integrals are currentlyunknown analytically. Therefore the QCD corrections have been calculated in variousapproximations 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 quarkmass effects at NLO QCD have been considered in the framework of a 1 /m t -expansionin Ref. [17], including Pad´e approximants constructed from expansion terms up to1 /m t . However, the 1 /m t -expansion becomes invalid for invariant masses of the HZ system larger than 2 m t , which is the interesting region with regards to new physicssearches. Soft gluon resummation for the gg → ZH process has been calculated inRef. [18].In Ref. [19], a data-driven strategy to extract the gluon-initiated component (or, moreprecisely, the non-Drell-Yan component) for ZH production has been suggested, basedon the comparison of the ZH to the W H cross section and the corresponding invariantmass distributions of the
V H system.Considering the full process pp → ZH , inclusive NNLO QCD corrections have beenavailable for quite some time [20] and are implemented in the program VH@NNLO [20–23]. NLO electroweak corrections have been calculated in Refs. [24, 25], combined NLOQCD+EW corrections are also available [26–28].Differential QCD corrections to the process pp → ZH have been calculated up toNNLO, including H → b ¯ b decays at NLO [29, 30] and at NNLO [31]. In Ref. [32], fullydifferential NNLO results for V H observables including the decays of the vector bosoninto leptons and the Higgs boson into b -quarks with off-shell propagators of the vector-and Higgs bosons have been calculated.The combination of fixed-order QCD computations with parton showers has been stud-ied at NLO+PS in association with up to one jet [27, 33–35] and at NNLOPS [36, 37].In addition, the NNLO corrections have been combined with NNLL resummation inthe 0-jettiness variable and matched to a parton shower within the Geneva
MonteCarlo framework [38].Threshold corrections up to N LL for the Drell-Yan-type part of the inclusive crosssection have been calculated in Ref. [39], soft-gluon resummation of both threshold– 2 –ogarithms and logarithms which are important at low transverse momentum of the
V H pair have been considered up to N LL in Ref. [40] and have been found to bevery close to the fixed order NNLO result. The process b ¯ b → ZH in the five-flavourscheme, but with a non-vanishing bottom-quark Yukawa coupling, has been calculatedin the soft-virtual approximation at NNLO QCD in Ref. [41], the polarised q ¯ q → 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 calculationand display numerical results for the two-loop amplitudes. A phenomenological studybased on these results is postponed to a subsequent publication. The structure of thiswork is as follows. In Section 2 we describe details of the calculation such as the integralfamilies, our projection operators, the reduction to master integrals and the evaluationof the master integrals. In Section 3 we show results for the finite part of the virtualamplitude, for some benchmark points as well as in terms of two-dimensional grids,before we conclude in Section 4. We consider the process g ( p ) + g ( p ) → Z ( p ) + H ( p ) , (2.1)with p = p = 0, p = m Z and p = m H . The Mandelstam invariants are defined by s = ( p + p ) , t = ( p − p ) , u = ( p − p ) , (2.2)with s + t + u = m Z + m H . The diagrams contributing to the process (2.1) at lead-ing 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 needsto be taken into account. We treat all quarks except the top-quark as massless, there-fore only top-quark loops contribute in those diagrams where the Higgs boson couplesdirectly to the fermion loop. The effect of a finite bottom quark mass on the total LOcross section is at the per mille level [17].In the triangle diagrams, only the axial vector part of the Z -boson coupling contributesdue to Furry’s theorem. In addition, the massless quark contributions cancel in eachisospin doublet, such that only the third generation of quarks give a non-zero contri-bution for this class of diagram. – 3 – G Z Z ZH H
Figure 1 : Diagrams contributing to gg → ZH at leading order. Diagrams related bycrossings are not shown.The leading order amplitude can be expressed in terms of seven form factors [11] con-taining one-loop three- and four-point functions. Some of the form factors can berelated by crossing p ↔ p 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 willbe explained in Section 2.1.1.In the m t → ∞ limit the amplitude simplifies considerably and is given by [16] M = − α s α sin θ w cos θ w m Z δ ab (cid:15) ( ε , ε , p , p ) p · ε ∗ Z s , (2.3)where (cid:15) is the totally antisymmetric Levi-Civita symbol with (cid:15) ( a, b, c, d ) ≡ (cid:15) µνρσ a µ b ν c ρ d σ and ε , ε 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- andcolour-averaged Born matrix element then can be written as | M | = α s α
256 sin θ w cos θ w m Z λ ( s, m Z , m H ) , (2.4)where λ denotes the K¨all´en function, λ ( a, b, c ) = a + b + c − ab − ac − bc .We calculate using conventional dimensional regularisation (CDR) with D = 4 − (cid:15) .For the treatment of γ within dimensional regularisation we use the ’t Hooft–Veltmanscheme [43, 44] in the variant of Refs. [45, 46], i.e. we use γ = i (cid:15) µνρσ γ µ γ ν γ ρ γ σ ,
12 ( γ µ γ − γ γ µ ) = i (cid:15) µνρσ γ ν γ ρ γ σ (2.5)and add finite renormalisation terms to restore chiral symmetry in the massless quarklimit, see Section 2.2.3. The contraction of two (cid:15) -symbols leads to linear combinationsof metric tensors which are treated as D -dimensional.– 4 – .1.1 Tensor structures and projection to a basis of linear polarisations We define the tensor amplitude A µ µ µ by extracting the polarisation vectors from theamplitude of process (2.1), A = ε µ λ ( p ) ε µ λ ( p ) ( ε µ λ ( p )) (cid:63) A µ µ µ , (2.6)where the ε µ i λ i denote the polarisation vectors. The Lorentz tensor structures appearingin the amplitude A µ µ µ were discussed in Ref. [11]. However, we do not use formfactors related to these Lorentz structures here, but rather use projections based on themomentum-basis representations of the linear polarisation vectors of external particlesas suggested in Ref. [47]. We also use the fact that only one axial current is involvedin the QCD corrections to this amplitude, therefore all relevant Lorentz structurescontain only a single Levi-Civita symbol. In addition, conditions such as transversalityand Bose symmetry regarding the two external gluons further constrain the possibleLorentz structures.Following the procedure of Ref. [47], we define the following normalised linear polari-sation vectors, where the frame that has been used is shown in Fig. 2, ε µx = N x ( − s p µ − s p µ + s p µ ) ,ε µy = N y (cid:0) (cid:15) µµ µ µ p µ p µ p µ (cid:1) ,ε µT = N T (cid:0)(cid:0) − s ( s + s ) + 2 m Z s (cid:1) p µ + (cid:0) s ( s + s ) − m Z s (cid:1) p µ + s ( − s + s ) p µ ) ,ε µl = N l (cid:0) − m Z ( p µ + p µ ) + ( s + s ) p µ (cid:1) , (2.7)with s ≡ p · p = s , s ≡ p · p = s + t − m H , s ≡ p · p = m Z − t .The normalisation factors N i for i ∈ { x, y, T, l } , associated with each of these space-like polarisation vectors in Eq. (2.7), can be determined from their (negative) normsquares ε i , i.e. N i = 1 / (cid:112) − ε i , which we choose to include only at the very end of thecalculation.The physical meaning of these vectors is apparent in the center-of-mass frame of p and p , where the beam axis determined by { (cid:126)p , (cid:126)p } is taken as the z -axis and theplane determined by { (cid:126)p , (cid:126)p } defines the x -O- z plane, as illustrated in Fig. 2. Thenthe 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 (cid:126)p , and ε l , the longitudinal polarisationof the Z -boson, points along the direction of (cid:126)p in the center-of-mass frame.Our projectors are given by products of three linear polarisation vectors, ε µ i ε µ j ε µ k with i , j ∈ { x, y } and k ∈ { T, y, l } , (2.8)– 5 – igure 2 : The coordinate system in the center-of-mass frame of the two incominggluons.which are further re-written solely in terms of external momenta and the Levi-Civitasymbol, all with D -dimensional Lorentz indices. However, only the six projectors withan odd number of Levi-Civita symbol, i.e. an odd number of ε y , are relevant forthe process (2.1). This is because there are only three linearly independent externalmomenta in A µ µ µ and all Lorentz structures appearing in the amplitude contain aLevi-Civita symbol. Concretely, the six projectors we use for A µ µ µ are given by P µ µ µ = ε µ x ε µ x ε µ y , P µ µ µ = ε µ x ε µ y ε µ T , P µ µ µ = ε µ x ε µ y ε µ l , P µ µ µ = ε µ y ε µ x ε µ T , P µ µ µ = ε µ y ε µ x ε µ l , P µ µ µ = ε µ y ε µ y ε µ y . (2.9)The projections made with P µ µ µ and P µ µ µ are related to those with P µ µ µ and P µ µ µ , respectively, through crossing of p and p . Therefore we only need to computefour projections in total. Note that all open Lorentz indices in Eq. (2.9) are understoodto be D -dimensional, and pairs of Levi-Civita symbols in P µ µ µ need to be contracted before being used for the projection of the amplitude (see Ref. [47] for a more detaileddiscussion).From the point of view of a form factor decomposition, the set of linear-polarisationprojectors in Eq. (2.9) represent precisely the Lorentz tensor decomposition basis inuse. This is a consequence of the orthogonality of the projectors (although the linearcompleteness is in general only ensured in the 4-dimensional limit [41, 47]).This projection defines six quantities A n ≡ P µ µ µ n A µ µ µ ( n = 1 , . . . , . (2.10)The physical interpretation of the quantities A n as linearly polarised amplitudes offersus a convenient short-cut to transform to a helicity basis defined w.r.t the same reference– 6 –rame. The usual helicity amplitudes for the process (2.1) can be constructed from thelinear ones, using the relations ε µ ± ( p ) = √ (cid:0) ε µ x ± iε µ y (cid:1) ,ε µ ± ( p ) = √ (cid:0) ε µ x ∓ iε µ y (cid:1) ,ε µ ± ( p ) = √ (cid:0) ε µ T ± iε µ y (cid:1) (2.11)and the longitudinal polarisation of the Z -boson is given by ε µ l in Eq. (2.7). The kinematic invariants in terms of the scattering angle θ between the beam axis and (cid:126)p in the centre-of-mass frame are given by t = − (cid:18) s − m Z − m H + cos θ (cid:113) λ ( s, m Z , m H ) (cid:19) ,u = − (cid:18) s − m Z − m H − cos θ (cid:113) λ ( s, m Z , m H ) (cid:19) , (2.12)where λ ( s, m Z , m H ) = ( s − ( m Z + m H ) ) ( s − ( m Z − m H ) ) is the K¨all´en function.Defining m + = m Z + m H , m − = m Z − m H and β ZH = (cid:16) − m s (cid:17) (cid:16) − m − s (cid:17) , we canalso write t = 12 (cid:0) m Z + m H (cid:1) − s β ZH cos θ ) ,u = 12 (cid:0) m Z + m H (cid:1) − s − β ZH cos θ ) . (2.13)The limit β ZH → s (cid:29) m , m − . In the forward scattering region θ →
0, if in addition β ZH → | t | (cid:29) | u | . Analogously, for θ → π and β ZH →
1, the ratio t/u is very small.The virtual amplitude has a threshold at s = 4 m t , therefore it is also useful to define β t = s − m t s + 4 m t − m Z + m H ) , (2.14)such that − ≤ β t ≤ β t = − s = s min = ( m Z + m H ) , β t = 0 at the topquark pair production threshold s = 4 m t and β t → s → ∞ .For 2 → p T of the Higgs- or Z -boson (with (cid:126)p HT = − (cid:126)p ZT ), fulfils 0 ≤ p T ≤ √ λ/ (2 √ s ) and s p T = tu − m Z m H . Therefore, in the highenergy limit, p T → s (1 − cos θ ). – 7 – .2 Diagram Generation, reduction and calculation of the master integrals We generate the 1- and 2-loop diagrams using QGRAF [48]. The projectors describedin Section 2.1.1 are applied and the Feynman rules are inserted using FORM [49–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 additionalfamilies obtained from the original families by exchanging p and p . Each familycontains nine propagators which allows all irreducible scalar products in the numeratorto be expressed in terms of inverse propagators. Concretely, the occurring integralshave the form I = ( µ ) (cid:15) (cid:90) d D k iπ D (cid:90) d D k iπ D N r j . . . N r t j t N − s j t +1 . . . N − s n − t j n , (2.15)where the N j denote genuine propagators of the form 1 / ( k − m ) with exponents r i ≥ s i ≥ p and p are not shown. The integrals appearing in the amplitudeare reduced to a minimal set of master integrals as described in the following. The numerical evaluation of finite integrals is typically much simpler than that of theirdivergent counterparts, we therefore follow the strategy of Ref. [52] to obtain a quasi-finite basis of master integrals. This choice of master integrals is not unique and the sizeof 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 thereduction tables factorizes into factors containing only the space-time dimension D andfactors depending on the kinematic invariants only. While algorithms to find such abasis have been presented in Refs. [53, 54], we obtained such a basis following a differentapproach, namely, by iterating over different combinations of master integral candidatesand analysing the resulting reduction tables, restricted to a small subset of the full IBPsystem and neglecting subsectors. This allows us to define additional criteria for theselection of the preferred masters, such as the size of the appearing denominator factorsor the order in (cid:15) = (4 − D ) / D -dependence of the denominator factors of thereduction rules factors for all integrals, except for some one-particle reducible integrals.Furthermore, all seven-propagator integrals only start contributing to the amplitudeat order (cid:15) − and the size of the amplitude coefficients reduces by about a factor of 5compared to a default choice of finite masters with minimal propagator powers. In– 8 – F F k − m t k − m t k k − m t k − m t ( k − k ) − m t ( k − k ) ( k − k ) ( k + p ) ( k + p ) − m t ( k + p ) − m t ( k + p ) − m t ( k + p ) − m t ( k + p ) − m t ( k − p ) ( k − p ) − m t ( k − p ) − m t ( k − p ) − m t ( k − p ) − m t ( k − p ) − m t ( k − p − p ) − m t ( k − p − p ) − m t ( k − p − p ) − m t ( k + p + p ) ( k − p − p ) − m t ( k − p − p ) − m t ( k + p − p ) F F k − m t k k k − m t ( k − k ) − m t ( k − k ) − m t ( k + p ) − m t ( k + p ) ( k + p ) ( k + p ) − m t ( k − p ) − m t ( k − p ) ( k − p ) ( k − p ) − m t ( k − p − p ) − m t ( k − p − p ) ( k − p − p ) ( k − p − p ) − m t Table 1 : Planar integral families used for the reduction. F F F k k − m t k − m t k k − m t k − m t ( k + p ) ( k − k ) ( k − k ) ( k − p ) ( k + p ) − m t ( k − p − p − p ) − m t ( k + p ) − m t ( k + p ) − m t ( k − p − p − p ) − m t ( k − k ) − m t ( k − p ) − m t ( k − p ) − m t ( k − p − p ) − m t ( k − p ) − m t ( k − p ) − m t ( k − k + p ) − m t ( k − k + p ) ( k − k + p ) ( k − k + p + p ) ( k − p − p ) − m t ( k − p − p ) − m t Table 2 : Non-planar integral families used for the reduction.particular the size of the coefficient of the two-propagator (double-tadpole) integralreduced from about 150 MB to less than 5 MB.– 9 –n 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 onecircumvents large intermediate expressions even during the calculation of IBP tables,which leads to a significant runtime and memory reduction in general. As a conse-quence, we interpolate the required reduction tables that relate the integrals of theamplitudes that contribute up to r = 7 and s = 4 to the desired set of master inte-grals that requires IBP relations up to r = 9 and s = 2 in D dimensions. To simplifythe reduction of the integrals occurring in the amplitude, we scale the Z and H massw.r.t. m t and approximate m Z /m t = 23 /
83 and m H /m t = 12 /
23, corresponding to m t = 173 .
21 GeV, m Z = 91 .
18 GeV, and m H = 125 . m t , the ratio of the Z -boson mass and Higgs boson mass to thetop quark mass is fixed in our calculation. By setting m t = 1 we remove an additionalscale that can be restored by dimensional analysis, which leaves us with a 3–parameterproblem.As we use quasi-finite master integrals in six dimensions, the set of dimensional recur-rence relations (DRRs), which connects integrals in six and four dimensions, has to becalculated as well. Hence, we utilize LiteRed [59] and
Reduze [60] in order to obtainthese 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 calcula-tion of the relations between the DRRs and our basis of master integrals was the mostdemanding step in the whole calculation. It took about four days of wall-clock timerunning on a machine with two Intel ® Xeon ® Silver 4116 and hyper-threading. Therequired memory never exceeded about 100 GB of memory.Afterwards the DRRs and the reduction tables obtained by reducing the integrals ofthe amplitudes are combined to a custom system of equations that fills roughly 9 GBof disk space. This system is again solved by employing
Kira with
FireFly in orderto interpolate the final set of replacement rules. Finally, the latter set of replacementsis 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, couldhave 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 memoryat once in this case, we observed faster runtime and lower memory consumption bysplitting the calculation as described above and running all steps on different machinesin parallel. Additionally, splitting the reduction into several steps is convenient whenstudying different bases of master integrals and their impact of the total file size of theresulting amplitudes. – 10 – .2.2 Evaluation of the two-loop amplitude
To evaluate the master integrals, we first apply sector decomposition as implementedin the program py
SecDec [61, 62]. For integrals which diverge in the limit of 4 space-time dimensions ( (cid:15) → (cid:15) leaving only finite integrals over the Feynman parameters which can then be integratednumerically. In the physical region, singularities can appear on the real axis of theFeynman parameters and a causal prescription to avoid the singularities is required.We deform the integration contour of the Feynman parameters into the complex planeas described in Refs. [63–66].In the present calculation we have selected a basis of quasi-finite integrals. Thismeans that poles in (cid:15) can appear only in the prefactor of our integrals after Feynmanparametrisation and, thus, sector decomposition is not required to resolve singulari-ties in the regulator. Nevertheless, we observe that applying a sector decompositiongreatly simplifies the numerical evaluation of the finite integrals. This observation canbe partly understood by noting that sector decomposition also removes integrable sin-gularities appearing at the boundaries where one or more Feynman parameters vanish.We therefore process all integrals with py SecDec before numerical integration.The numerical integration itself is performed using the quasi-Monte Carlo (QMC) algo-rithm [62, 67]. For all integrals we apply a Korobov periodising transform with weight3. We observe that the use of the rank-1 shifted lattice rules greatly reduces the numberof samples required to obtain the integrals to sufficient precision for the computationof 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 py SecDec program such that it can produce acode capable of evaluating the entire amplitude, rather than computing each integralseparately. The advantage of this structure is that the number of sampling points usedfor each sector of each integral can be dynamically set according to its contribution tothe total uncertainty on the amplitude. We utilise a variant of the procedure describedin Ref. [69] to minimise the total time taken to obtain a given relative accuracy on theamplitude.
We may expand each of the form factors, A i =1 ,...,n , in the bare strong coupling a = α / (4 π ) according to A i = a A (0) i + a A (1) i + O ( a ) . (2.16)We renormalise the strong coupling in the MS scheme with the heavy quark loop in thegluon self-energy subtracted at zero momentum. The heavy quark mass is renormalised– 11 –n the OS scheme. The UV renormalised amplitude is obtained via the replacement A UV i = Z Z n g / A A i ( a → a s Z α , m → m Z m , y t → y t Z m )= a s A (0) i + a s ( n g δZ A + δZ α + δZ ) A (0) i + a s δZ m A mct , (0) i + a s A (1) + O ( a s )(2.17)Here, n g = 2 is the number of external gluon legs, a s = α s ( µ R ) / (4 π ) S − (cid:15) ( µ R /µ ) (cid:15) with S (cid:15) = (4 π ) (cid:15) e − (cid:15)γ E , and the renormalisation constants are expanded according to Z j = 1 + a s δZ j + O ( a s ) ( j = α, A, m, δZ α = − (cid:15) β + δZ hq α , β = 113 C A − T R n F , (2.18) δZ A = − δZ hq α = (cid:18) µ R m (cid:19) (cid:15) (cid:18) − (cid:15) T R (cid:19) , (2.19) δZ m = (cid:18) µ R m (cid:19) (cid:15) C F (cid:18) − (cid:15) − (cid:19) , (2.20) δZ = (cid:40) δZ ,ns = − C F , diagrams involving a Z-top vertex δZ ,ρ = − C F , diagrams involving a Goldstone-top vertex . (2.21)We handle the γ matrices which appear in this calculation using the Larin scheme [46].According to this scheme an additional finite renormalisation is required, it is denotedin the equations above by δZ .At the level of individual form factors, the procedure outlined corresponds to the fol-lowing relations A UV i = (cid:16) α s π (cid:17) A (0) , UV i + (cid:16) α s π (cid:17) A (1) , UV i + O ( α s ) , A (0) , UV i = S − (cid:15) (cid:18) µ R µ (cid:19) (cid:15) A (0) i , A (1) , UV i = S − (cid:15) (cid:18) µ R µ (cid:19) (cid:15) A (1) i + (cid:18) − β (cid:15) + δZ (cid:19) A (0) , UV i + δZ m A mct , (0) i . (2.22) In order to obtain IR finite amplitudes we use the subtraction scheme described inRef. [72] A (0) , fin i = A (0) , UV i , (2.23) A (1) , fin i = A (1) , UV i − I A (0) , UV i , (2.24)– 12 –ith I = I soft1 + I coll1 , (2.25) I soft1 = − e (cid:15)γ E Γ(1 − (cid:15) ) (cid:18) µ R s (cid:19) (cid:15) (cid:18) (cid:15) + iπ(cid:15) (cid:19) C A , (2.26) I coll1 = − β (cid:15) (cid:18) µ R s (cid:19) (cid:15) . (2.27)We present results at renormalisation scale µ R = s and after changing to the helicitybasis, we define B = 8 T R (2 · (cid:88) { λ i } A (0) , fin { λ i } A (cid:63) (0) , fin { λ i } , (2.28) V = 8 T R (2 · (cid:88) { λ i } (cid:16) A (0) , fin { λ i } A (cid:63) (1) , fin { λ i } + A (1) , fin { λ i } A (cid:63) (0) , fin { λ i } (cid:17) (2.29)for the square of the amplitude. The electroweak coupling α is set to 1. We have verified that our results have the expected universal pole structure and crossingsymmetries. We also compared our exact result with the expansion in the large topquark mass limit presented in [17]. The virtual contributions (cid:101) V n are expanded up toorder 1 /m nt and reweighted with the exact Born term B . The one-particle-reducibledouble triangle contributions V red are included with full top quark mass dependence. V n = BB n (cid:101) V n + V red (3.1)In Fig. 3 the ratio of expanded to full virtual contribution V n / V is shown for expansionorders 0 ≤ n ≤ θ ) ≈ . β t = − s = ( m Z + m H ) ) the expanded result approximatesthe exact calculation well with a ratio V / V ≈ . β t = 0 ( s = 4 m t ), where the large m t -expansionis expected to break down.In addition, we have compared our results to a recent calculation in the high energylimit [73] and found agreement to the extent expected from previous comparisons forthe process pp → HH performed in Ref. [74].– 13 – . − . − . − . − . . β t . . . . . . V n / V n = 0 n = 1 n = 2 n = 3 n = 4 Figure 3 : Comparison of exact virtual contributions V with the expansion in the largetop quark mass limit V n from [17] for fixed scattering angle cos( θ ) ≈ . − ≤ β t ≤ m Z + m H ) ≤ s ≤ m t between production threshold and top quark pair threshold. For the presentation of our results, we have evaluated a total of 460 phase-space pointsat 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 ofrun 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 thefinite part of the virtual two-loop amplitude interfered with the Born amplitude. Wesee that the unpolarised Born result shows a rather flat dependence on the scatteringangle around and below the top quark pair threshold, due to the dominating s-wavecontribution in this region, and starts to curve as energy further increases becausepartial waves with higher angular momentum play an increasingly important role. Wealso observe that in the two-loop case, the top quark pair production threshold region ismuch more peaked than at leading order, due to log β t -terms appearing for the first timeat two-loop order. In Fig. 5, the ratio of the two-loop amplitude to the Born-amplitudeis shown separately. – 14 – t − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) -1.00-0.75-0.50-0.250.000.250.500.751.00 B total β t − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) -1.00-0.75-0.50-0.250.000.250.500.751.00 V total Figure 4 : Dependence of the leading order (left) and virtual contribution (right) onthe parameters β t and cos θ , summing over all polarisation. β t − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) -1.00-0.75-0.50-0.250.000.250.500.751.00 V / B− total Figure 5 : Dependence of the ratio V / B on the parameters β t and cos θ , summing overall polarisations.In the following we will show results for five helicity amplitudes in the β t –cos θ –plane.The remaining amplitudes can be obtained from overall helicity flips as well as cross-ing between the two initial gluons. Bose-symmetry and the behaviour under paritytransformations imply A λ ,λ ,λ z ( t, u ) = ( − λ z A λ ,λ ,λ z ( u, t ) , (3.2) A λ ,λ ,λ z ( t, u ) = −A − λ , − λ , − λ z ( t, u ) . To understand the qualitative behaviour, we can expand the amplitude in terms of– 15 –igner d -functions, A { λ i } ( θ ) = ∞ (cid:88) J =0 (2 J + 1) A J { λ i } d Js i ,s f ( θ ) , (3.3)where s i and s f denote the initial and final state total spins, respectively, and J denotesthe total angular momentum of the system. For gg → ZH , the initial state has totalspin s i = 0 for equal helicities of the initial state gluons, ( λ , λ ) = (+ , +) or ( − , − ),while for the case ( λ , λ ) = (+ , − ) or ( − , +), the initial state has total spin s i = 2.Therefore the amplitude A ++0 is dominated by the partial wave d ( θ ) and providesthe largest contribution to the total squared amplitude. In particular in the low energyregion, where the ZH system has relatively small kinetic energy, the s-wave contributionshould dominate, reflected in the homogeneity of A ++0 in cos θ . As the center-of-massenergy increases, the contributions from partial waves with higher angular momentaalso start to play a role, leading to a non-flat behaviour in cos θ . Note that eq. (3.2)implies that A ++0 is symmetric under exchange of t ↔ u and therefore is symmetric incos θ . As A ++0 is composed of partial waves d J ( θ ), which are even in cos θ for even J ,no partial wave components with odd J can contribute to A ++0 .From Fig. 6, we further observe that the helicity amplitudes with the polarisations (cid:15) Z ± are suppressed compared to those with a longitudinally polarised Z -boson. Theamplitudes A ++ ± are antisymmetric under exchange of t and u , so antisymmetric incos θ . The contributions with J = 1, being proportional to sin θ , therefore do notoccur in A ++ ± . The d-wave contributions d , ± ( θ ) 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, suchthat the shapes cannot be explained by partial waves alone. Note that A ++ − is aboutfive orders of magnitude smaller than A ++0 , and A +++ is also suppressed, thereforethe amplitudes A ++ ± give a very minor contribution in the sum of all polarisationconfigurations.For different helicities of the initial state gluons, ( λ , λ ) = (+ , − ) or ( − , +), the initialstate has total spin s i = 2, such that the partial wave contributions start from J = 2.Therefore the amplitude A + − is much smaller than A ++0 . This is also reflected inthe fact that the amplitudes are basically zero except at very high energies. Theamplitude A + − has no contribution from d J , ( θ ) with even J as it is antisymmetricin cos θ , therefore its leading partial wave is given by d , ( θ ) ∼ cos θ sin θ , vanishingat cos θ = ± θ = 0. The amplitudes A + −± have their leading partial wavesgiven by d , ± ( θ ) ∼ sin θ (1 ± cos θ ). Consequently A + − + is highly suppressed in thebackward direction, while A + −− is highly suppressed in the forward direction.– 16 – t − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) -1.00-0.75-0.50-0.250.000.250.500.751.00 B ++0 β t − . − . − . − .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) -1.00-0.75-0.50-0.250.000.250.500.751.00 V ++0 β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . B . . . . . +++ β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . V +++ β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . B . . . . . . . . ++- β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . V . . . . . . . . . ++- Figure 6 : Dependence of the leading order (left) and virtual contribution (right) onthe parameters β t and cos( θ ) for the individual helicity amplitudes with s i = 0.In Table 3 we list numerical results of the amplitude to facilitate comparisons with ourresults. – 17 – t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . B +-0 β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . V +-0 β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . B . . . . . . . +-+ β t - . - . - . - .
25 0 .
00 0 .
25 0 .
50 0 .
75 1 . c o s ( θ ) - . - . - . - . . . . . . V +-+ Figure 7 : Dependence of the leading order (left) and virtual contribution (right) onthe parameters β t and cos( θ ) for the individual helicity amplitudes with s i = 2. s/m t t/m t B V
Table 3 : Numerical results for various phase space points at the scale µ R = s . Thenumber in parentheses gives the numerical uncertainty on the last digit of the virtualamplitude. – 18 – Conclusions and Outlook
We have numerically calculated the two-loop amplitudes for the production of a Higgs-and a Z -boson in gluon fusion with massive top quark loops. The results for the finitepart of the two-loop amplitude interfered with the Born amplitude are plotted as afunction of the scattering angle and the centre-of-mass energy, for the total unpolarisedamplitude as well as for individual helicity amplitudes.The projection of the amplitudes to scalar quantities has been carried out with projec-tors 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], usingin addition
LiteRed [59] and
Reduze [60] to obtain dimensional recurrence relations.The master integrals have been calculated using py
SecDec [61, 62]. The integrationis sufficiently stable and accurate, also in the near-threshold and forward scatteringregions, to perform phenomenological studies based on these results, after includingthe real radiation contributions. We postpone such a phenomenological analysis to asubsequent publication.Our method for the first time has been applied to a process with three different massscales, m t , m H and m Z . We expect that it can be applied successfully to other two-loopamplitudes involving several mass scales in the future. Acknowledgements
We would like to thank Joshua Davies, Go Mishima and Matthias Steinhauser for thecomparison of results prior to publication. We also would like to thank Tom Zirkefor useful conversation about the large mass expansion. This research was supportedby the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) undergrant 396021762 - TRR 257, by the COST Action CA16201 (‘Particleface’) of theEuropean Union and by the Swiss National Science Foundation (SNF) under grantnumber 200020-175595. MK acknowledges support by the Forschungskredit of theUniversity of Zurich, grant no. FK-19-102. The research of JS was supported bythe European Union through the ERC Advanced Grant MC@NNLO (340983). SJ issupported by a Royal Society University Research Fellowship (Grant URF/R1/201268).We gratefully acknowledge resources provided by the Max Planck Computing and DataFacility (MPCDF). Some of the calculations were performed with computing resourcesgranted by RWTH Aachen University under project rwth0541.– 19 – eferences [1]
ATLAS collaboration, M. Aaboud et al.,
Observation of H → b ¯ b decays and V H production with the ATLAS detector , Phys. Lett. B (2018) 59 [ ].[2]
CMS collaboration, A. M. Sirunyan et al.,
Observation of Higgs boson decay to bottomquarks , Phys. Rev. Lett. (2018) 121801 [ ].[3] K. Mimasu, V. Sanz and C. Williams,
Higher Order QCD predictions for AssociatedHiggs production with anomalous couplings to gauge bosons , JHEP (2016) 039[ ].[4] A. Greljo, G. Isidori, J. M. Lindert, D. Marzocca and H. Zhang, Electroweak Higgsproduction with HiggsPO at NLO QCD , Eur. Phys. J. C (2017) 838 [ ].[5] S. Alioli, W. Dekens, M. Girard and E. Mereghetti, NLO QCD corrections to SM-EFTdilepton and electroweak Higgs boson production, matched to parton shower inPOWHEG , JHEP (2018) 205 [ ].[6] S. Banerjee, C. Englert, R. S. Gupta and M. Spannowsky, Probing ElectroweakPrecision Physics via boosted Higgs-strahlung at the LHC , Phys. Rev. D (2018)095012 [ ].[7] S. Banerjee, R. S. Gupta, J. Y. Reiness and M. Spannowsky, Resolving the tensorstructure of the Higgs coupling to Z -bosons via Higgs-strahlung , Phys. Rev. D (2019) 115004 [ ].[8] J. Baglio, S. Dawson, S. Homiller, S. D. Lane and I. M. Lewis,
Validity of standardmodel EFT studies of VH and VV production at NLO , Phys. Rev. D (2020)115004 [ ].[9] K. Rao, S. D. Rindani and P. Sarmah,
Study of anomalous gauge-Higgs couplings usingZ boson polarization at LHC , .[10] D. A. Dicus and C. Kao, Higgs Boson - Z Production From Gluon Fusion , Phys. Rev.D (1988) 1008.[11] B. A. Kniehl, Associated Production of Higgs and Z Bosons From Gluon Fusion inHadron Collisions , Phys. Rev. D (1990) 2253.[12] R. V. Harlander, S. Liebler and T. Zirke, Higgs Strahlung at the Large Hadron Colliderin the 2-Higgs-Doublet Model , JHEP (2014) 023 [ ].[13] C. Englert, M. McCullough and M. Spannowsky, Gluon-initiated associated productionboosts Higgs physics , Phys. Rev. D (2014) 013013 [ ].[14] ATLAS collaboration, G. Aad et al.,
Measurement of the associated production of aHiggs boson decaying into b -quarks with a vector boson at high transverse momentumin pp collisions at √ s = 13 TeV with the ATLAS detector , . – 20 – CMS collaboration, A. M. Sirunyan et al.,
Search for the associated production of theHiggs boson and a vector boson in proton-proton collisions at √ s =
13 TeV via Higgsboson decays to τ leptons , JHEP (2019) 093 [ ].[16] L. Altenkamp, S. Dittmaier, R. V. Harlander, H. Rzehak and T. J. Zirke, Gluon-induced Higgs-strahlung at next-to-leading order QCD , JHEP (2013) 078[ ].[17] A. Hasselhuhn, T. Luthe and M. Steinhauser, On top quark mass effects to gg → ZH at NLO , JHEP (2017) 073 [ ].[18] R. V. Harlander, A. Kulesza, V. Theeuwes and T. Zirke, Soft gluon resummation forgluon-induced Higgs Strahlung , JHEP (2014) 082 [ ].[19] R. Harlander, J. Klappert, C. Pandini and A. Papaefstathiou, Exploiting the WH/ZHsymmetry in the search for New Physics , Eur. Phys. J. C (2018) 760 [ ].[20] O. Brein, A. Djouadi and R. Harlander, NNLO QCD corrections to the Higgs-strahlungprocesses at hadron colliders , Phys. Lett. B (2004) 149 [ hep-ph/0307206 ].[21] O. Brein, R. Harlander, M. Wiesemann and T. Zirke,
Top-Quark Mediated Effects inHadronic Higgs-Strahlung , Eur. Phys. J. C (2012) 1868 [ ].[22] O. Brein, R. V. Harlander and T. J. Zirke, vh@nnlo - Higgs Strahlung at hadroncolliders , Comput. Phys. Commun. (2013) 998 [ ].[23] R. V. Harlander, J. Klappert, S. Liebler and L. Simon, vh@nnlo-v2: New physics inHiggs Strahlung , JHEP (2018) 089 [ ].[24] M. Ciccolini, S. Dittmaier and M. Kr¨amer, Electroweak radiative corrections toassociated WH and ZH production at hadron colliders , Phys. Rev. D (2003) 073003[ hep-ph/0306234 ].[25] A. Denner, S. Dittmaier, S. Kallweit and A. M¨uck, Electroweak corrections toHiggs-strahlung off W/Z bosons at the Tevatron and the LHC with HAWK , JHEP (2012) 075 [ ].[26] A. Denner, S. Dittmaier, S. Kallweit and A. M¨uck, HAWK 2.0: A Monte Carloprogram for Higgs production in vector-boson fusion and Higgs strahlung at hadroncolliders , Comput. Phys. Commun. (2015) 161 [ ].[27] F. Granata, J. M. Lindert, C. Oleari and S. Pozzorini,
NLO QCD+EW predictions forHV and HV +jet production including parton-shower effects , JHEP (2017) 012[ ].[28] P. Obul, S. Dulat, T.-J. Hou, A. Tursun and N. Yalkun, Next-to-leading order QCDand electroweak corrections to Higgs-strahlung processes at the LHC , Chin. Phys. C (2018) 093105 [ ]. – 21 –
29] G. Ferrera, M. Grazzini and F. Tramontano,
Associated ZH production at hadroncolliders: the fully differential NNLO QCD calculation , Phys. Lett. B (2015) 51[ ].[30] J. M. Campbell, R. K. Ellis and C. Williams,
Associated production of a Higgs boson atNNLO , JHEP (2016) 179 [ ].[31] G. Ferrera, G. Somogyi and F. Tramontano, Associated production of a Higgs bosondecaying into bottom quarks at the LHC in full NNLO QCD , Phys. Lett. B (2018)346 [ ].[32] R. Gauld, A. Gehrmann-De Ridder, E. Glover, A. Huss and I. Majer,
Associatedproduction of a Higgs boson decaying into bottom quarks and a weak vector bosondecaying leptonically at NNLO in QCD , JHEP (2019) 002 [ ].[33] G. Luisoni, P. Nason, C. Oleari and F. Tramontano, HW ± /HZ + 0 and 1 jet at NLOwith the POWHEG BOX interfaced to GoSam and their merging within MiNLO , JHEP (2013) 083 [ ].[34] D. Goncalves, F. Krauss, S. Kuttimalai and P. Maierh¨ofer, Higgs-Strahlung: Mergingthe NLO Drell-Yan and Loop-Induced 0+1 jet Multiplicities , Phys. Rev. D (2015)073006 [ ].[35] B. Hespel, F. Maltoni and E. Vryonidou, Higgs and Z boson associated production viagluon fusion in the SM and the 2HDM , JHEP (2015) 065 [ ].[36] W. Astill, W. Bizo´n, E. Re and G. Zanderighi, NNLOPS accurate associated HZproduction with H → bb decay at NLO , JHEP (2018) 157 [ ].[37] W. Bizo´n, E. Re and G. Zanderighi, NNLOPS description of the H → bb decay withMiNLO , JHEP (2020) 006 [ ].[38] S. Alioli, A. Broggio, S. Kallweit, M. A. Lim and L. Rottoli, Higgsstrahlung atNNLL (cid:48) + NNLO matched to parton showers in GENEVA , Phys. Rev. D (2019)096016 [ ].[39] M. Kumar, M. Mandal and V. Ravindran,
Associated production of Higgs boson withvector boson at threshold N LO in QCD , JHEP (2015) 037 [ ].[40] S. Dawson, T. Han, W. Lai, A. Leibovich and I. Lewis, Resummation Effects inVector-Boson and Higgs Associated Production , Phys. Rev. D (2012) 074007[ ].[41] T. Ahmed, A. Ajjath, L. Chen, P. K. Dhani, P. Mukherjee and V. Ravindran, Polarised Amplitudes and Soft-Virtual Cross Sections for b ¯ b → ZH at NNLO in QCD , JHEP (2020) 030 [ ].[42] T. Ahmed, W. Bernreuther, L. Chen and M. Czakon, Polarized q ¯ q → Z + Higgsamplitudes at two loops in QCD: the interplay between vector and axial vector form – 22 – actors and a pitfall in applying a non-anticommuting γ , JHEP (2020) 159[ ].[43] G. ’t Hooft and M. Veltman, Regularization and Renormalization of Gauge Fields , Nucl. Phys. B (1972) 189.[44] P. Breitenlohner and D. Maison, Dimensional Renormalization and the ActionPrinciple , Commun. Math. Phys. (1977) 11.[45] S. A. Larin and J. A. M. Vermaseren, The α s corrections to the Bjorken sum rule forpolarized electroproduction and to the Gross-Llewellyn Smith sum rule , Phys. Lett.
B259 (1991) 345.[46] S. Larin,
The Renormalization of the axial anomaly in dimensional regularization , Phys. Lett. B (1993) 113 [ hep-ph/9302240 ].[47] L. Chen,
A prescription for projectors to compute helicity amplitudes in D dimensions , .[48] P. Nogueira, Automatic Feynman graph generation , J. Comput. Phys. (1993) 279.[49] J. Vermaseren,
New features of FORM , math-ph/0010025 .[50] J. Kuipers, T. Ueda, J. Vermaseren and J. Vollinga, FORM version 4.0 , Comput. Phys.Commun. (2013) 1453 [ ].[51] B. Ruijl, T. Ueda and J. Vermaseren,
FORM version 4.2 , .[52] A. von Manteuffel, E. Panzer and R. M. Schabinger, A quasi-finite basis for multi-loopFeynman integrals , JHEP (2015) 120 [ ].[53] A. Smirnov and V. Smirnov, How to choose master integrals , Nucl. Phys. B (2020)115213 [ ].[54] J. Usovitsch,
Factorization of denominators in integration-by-parts reductions , .[55] P. Maierh¨ofer, J. Usovitsch and P. Uwer, Kira—A Feynman integral reduction program , Comput. Phys. Commun. (2018) 99 [ ].[56] J. Klappert, F. Lange, P. Maierh¨ofer and J. Usovitsch,
Integral Reduction with Kira 2.0and Finite Field Methods , .[57] J. Klappert and F. Lange, Reconstructing rational functions with FireFly , Comput.Phys. Commun. (2020) 106951 [ ].[58] J. Klappert, S. Y. Klein and F. Lange,
Interpolation of Dense and Sparse RationalFunctions and other Improvements in
FireFly , .[59] R. N. Lee, LiteRed 1.4: a powerful tool for reduction of multiloop integrals , J. Phys.Conf. Ser. (2014) 012059 [ ]. – 23 –
60] A. von Manteuffel and C. Studerus,
Reduze 2 - Distributed Feynman IntegralReduction , .[61] S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner, J. Schlenk et al., pySecDec: atoolbox for the numerical evaluation of multi-scale integrals , Comput. Phys. Commun. (2018) 313 [ ].[62] S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner and J. Schlenk,
A GPUcompatible quasi-Monte Carlo integrator interfaced to pySecDec , Comput. Phys.Commun. (2019) 120 [ ].[63] D. E. Soper,
Techniques for QCD calculations by numerical integration , Phys. Rev. D (2000) 014009 [ hep-ph/9910292 ].[64] T. Binoth, J. Guillet, G. Heinrich, E. Pilon and C. Schubert, An Algebraic/numericalformalism for one-loop multi-leg amplitudes , JHEP (2005) 015 [ hep-ph/0504267 ].[65] Z. Nagy and D. E. Soper, Numerical integration of one-loop Feynman diagrams forN-photon amplitudes , Phys. Rev. D (2006) 093006 [ hep-ph/0610028 ].[66] S. Borowka, J. Carter and G. Heinrich, Numerical Evaluation of Multi-Loop Integralsfor Arbitrary Kinematics with SecDec 2.0 , Comput. Phys. Commun. (2013) 396[ ].[67] Z. Li, J. Wang, Q.-S. Yan and X. Zhao,
Efficient numerical evaluation of Feynmanintegrals , Chin. Phys. C (2016) 033103 [ ].[68] S. Borowka, N. Greiner, G. Heinrich, S. Jones, M. Kerner, J. Schlenk et al., HiggsBoson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-QuarkMass Dependence , Phys. Rev. Lett. (2016) 012001 [ ].[69] S. Borowka, N. Greiner, G. Heinrich, S. Jones, M. Kerner, J. Schlenk et al.,
Full topquark mass dependence in Higgs boson pair production at NLO , JHEP (2016) 107[ ].[70] S. Jones, M. Kerner and G. Luisoni, Next-to-Leading-Order QCD Corrections to HiggsBoson Plus Jet Production with Full Top-Quark Mass Dependence , Phys. Rev. Lett. (2018) 162001 [ ].[71] L. Chen, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk et al.,
Photon pairproduction in gluon fusion: Top quark effects at NLO with threshold matching , JHEP (2020) 115 [ ].[72] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, Universality oftransverse-momentum resummation and hard factors at the NNLO , Nucl. Phys. B (2014) 414 [ ].[73] J. Davies, G. Mishima and M. Steinhauser,
Virtual corrections to gg → ZH in thehigh-energy and large- m t limits , . – 24 –
74] J. Davies, G. Heinrich, S. P. Jones, M. Kerner, G. Mishima, M. Steinhauser et al.,
Double Higgs boson production at NLO: combining the exact numerical result andhigh-energy expansion , JHEP (2019) 024 [ ].[75] S. Dittmaier, Weyl-van der Waerden formalism for helicity amplitudes of massiveparticles , Phys. Rev. D (1998) 016007 [ hep-ph/9805445 ].[76] B. A. Kniehl and C. P. Palisoc, Associated production of Z and neutral Higgs bosons atthe CERN Large Hadron Collider , Phys. Rev. D (2012) 075027 [ ].].