aainstitutetext: Max Planck Institute for Physics, Föhringer Ring 6, 80805
München, Germanybbinstitutetext: Institute for Theoretical Physics, Karlsruhe
Institute of Technology, 76128 Karlsruhe, Germanyccinstitutetext: Theoretical
Physics Department, CERN, Geneva, Switzerlandddinstitutetext: Institute for
Particle Physics Phenomenology, Durham University, Durham DH1 3LE,
UKeeinstitutetext: Physik-Institut, Universität Zürich, Winterthurerstrasse
190, 8057 Zürich, Switzerlandffinstitutetext: Institute for Theoretical
Particle Physics and Cosmology, RWTH Aachen University, 52056 Aachen,
Germanygginstitutetext: Theory Group LTP, Paul Scherrer Institut, CH-5232
Villigen PSI, Switzerland

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

Long Chen b Gudrun Heinrich c,d Stephen P. Jones e Matthias Kerner f Jonas
Klappert g Johannes Schlenk<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

###### Keywords:

LHC, QCD phenomenology, two-loop computations, Higgs boson, vector bosons, top
quark mass

††preprint:  . ZU-TH 45/20
. CERN-TH-2020-199
. IPPP/20/57
. P3H-20-076
. KA-TP-21-2020
. TTK-20-42
. PSI-PR-20-21

##  1 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
Aaboud:2018zhk ; Sirunyan:2018kst , 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\to b\bar{b}$. The regime of boosted Higgs bosons is particularly
interesting, as the large-$p_{T,H}$ region is sensitive to new physics
Mimasu:2015nqa ; Greljo:2017spw ; Alioli:2018ljm ; Banerjee:2018bio ;
Banerjee:2019pks ; Baglio:2020oqu ; Rao:2020hel .

The loop-induced gluon-initiated contributions, calculated at LO in Refs.
Dicus:1988yh ; Kniehl:1990iva , are finite and enter at order
$\alpha_{s}^{2}$, i.e. formally at NNLO considering the $pp\to VH$ 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_{T}^{H}\gtrsim 150$
GeV Harlander:2013mla ; Englert:2013vua . As only the LO results for $gg\to
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 $VH$ process
Aaboud:2018zhk ; Sirunyan:2018kst ; Aad:2020eiv ; Sirunyan:2018cpi .
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.
Altenkamp:2012sx  they have been calculated in the $m_{t}\to\infty$ 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. Hasselhuhn:2016rqt , including Padé
approximants constructed from expansion terms up to $1/m_{t}^{8}$. 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\to ZH$ process has
been calculated in Ref. Harlander:2014wda .

In Ref. Harlander:2018yns , 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 $WH$
cross section and the corresponding invariant mass distributions of the $VH$
system.

Considering the full process $pp\to ZH$, inclusive NNLO QCD corrections have
been available for quite some time Brein:2003wg  and are implemented in the
program VH@NNLO Brein:2003wg ; Brein:2011vx ; Brein:2012ne ; Harlander:2018yio
. NLO electroweak corrections have been calculated in Refs. Ciccolini:2003jy ;
Denner:2011id , combined NLO QCD+EW corrections are also available
Denner:2014cla ; Granata:2017iod ; Obul:2018psx .

Differential QCD corrections to the process $pp\to ZH$ have been calculated up
to NNLO, including $H\to b\bar{b}$ decays at NLO Ferrera:2014lca ;
Campbell:2016jau  and at NNLO Ferrera:2017zex . In Ref. Gauld:2019yng , fully
differential NNLO results for $VH$ observables including the decays of the
vector boson into 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
studied at NLO+PS in association with up to one jet Luisoni:2013cuh ;
Goncalves:2015mfa ; Hespel:2015zea ; Granata:2017iod  and at NNLOPS
Astill:2018ivh ; Bizon:2019tfo . 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 Alioli:2019qzz .

Threshold corrections up to N3LL for the Drell-Yan-type part of the inclusive
cross section have been calculated in Ref. Kumar:2014uwa , soft-gluon
resummation of both threshold logarithms and logarithms which are important at
low transverse momentum of the $VH$ pair have been considered up to N3LL in
Ref. Dawson:2012gs  and have been found to be very close to the fixed order
NNLO result. The process $b\bar{b}\to 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. Ahmed:2019udm , the polarised
$q\bar{q}\to ZH$ two-loop amplitudes have been calculated in Ref.
Ahmed:2020kme .

In this paper, we calculate the two-loop virtual corrections to the process
$gg\to 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.

##  2 Details of the calculation

###  2.1 Amplitude definition

We consider the process

$\displaystyle g(p_{1})+g(p_{2})\to Z(p_{3})+H(p_{4}),$ (1)

with $p_{1}^{2}=p_{2}^{2}=0$, $p_{3}^{2}=m_{Z}^{2}$ and $p_{4}^{2}=m_{H}^{2}$.
The Mandelstam invariants are defined by

$\displaystyle
s=\left(p_{1}+p_{2}\right)^{2}\;,\;t=\left(p_{2}-p_{3}\right)^{2}\;,\;u=\left(p_{1}-p_{3}\right)^{2}\;,$
(2)

with $s+t+u=m_{Z}^{2}+m_{H}^{2}$. The diagrams contributing to the process (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 Hasselhuhn:2016rqt .

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
Kniehl:1990iva  containing one-loop three- and four-point functions. Some of
the form factors can be related by crossing $p_{1}\leftrightarrow 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}\to\infty$ limit the amplitude simplifies considerably and is
given by Altenkamp:2012sx

$\displaystyle
M_{0}=-\frac{\alpha_{s}\alpha}{\sin^{2}\theta_{\mathrm{w}}\cos^{2}\theta_{\mathrm{w}}\,m_{Z}}\,\delta^{ab}\epsilon(\varepsilon_{1},\varepsilon_{2},p_{1},p_{2})\,\frac{p_{4}\cdot\varepsilon_{Z}^{*}}{s}\;,$
(3)

where $\epsilon$ is the totally antisymmetric Levi-Civita symbol with
$\epsilon(a,b,c,d)\equiv\epsilon^{\mu\nu\rho\sigma}a_{\mu}b_{\nu}c_{\rho}d_{\sigma}$
and $\varepsilon_{1},\varepsilon_{2}$ are the polarisation vectors of the
incoming gluons, carrying colour indices $a$ and $b$, while
$\varepsilon_{Z}^{*}$ is the polarisation vector of the outgoing $Z$-boson.
The spin- and colour-averaged Born matrix element then can be written as

$\displaystyle\overline{|M_{0}|^{2}}=\frac{\alpha_{s}^{2}\alpha^{2}}{256\sin^{4}\theta_{\mathrm{w}}\cos^{4}\theta_{\mathrm{w}}m_{Z}^{4}}\,\lambda(s,m_{Z}^{2},m_{H}^{2})\;,$
(4)

where $\lambda$ denotes the Källén function,
$\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$.

Figure 1: Diagrams contributing to $gg\to ZH$ at leading order. Diagrams
related by crossings are not shown.

We calculate using conventional dimensional regularisation (CDR) with
$D=4-2\epsilon$. For the treatment of $\gamma_{5}$ within dimensional
regularisation we use the ’t Hooft–Veltman scheme tHooft:1972tcz ;
Breitenlohner:1977hr  in the variant of Refs. Larin:1991tj ; Larin:1993tq ,
i.e. we use

$\displaystyle\gamma_{5}$
$\displaystyle=\frac{i}{4!}\epsilon_{\mu\nu\rho\sigma}\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}\;,$
$\displaystyle\frac{1}{2}(\gamma_{\mu}\gamma_{5}-\gamma_{5}\gamma_{\mu})$
$\displaystyle=\frac{i}{3!}\epsilon_{\mu\nu\rho\sigma}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}$
(5)

and add finite renormalisation terms to restore chiral symmetry in the
massless quark limit, see Section 2.2.3. The contraction of two
$\epsilon$-symbols leads to linear combinations of metric tensors which are
treated as $D$-dimensional.

####  2.1.1 Tensor structures and projection to a basis of linear
polarisations

We define the tensor amplitude ${\mathcal{A}}_{\mu_{1}\mu_{2}\mu_{3}}$ by
extracting the polarisation vectors from the amplitude of process (1),

$\displaystyle{\mathcal{A}}{}$
$\displaystyle=\varepsilon_{\lambda_{1}}^{\mu_{1}}(p_{1})\,\varepsilon_{\lambda_{2}}^{\mu_{2}}(p_{2})\,(\varepsilon_{\lambda_{3}}^{\mu_{3}}(p_{3}))^{\star}\,{\mathcal{A}}_{\mu_{1}\mu_{2}\mu_{3}}\;,$
(6)

where the $\varepsilon_{\lambda_{i}}^{\mu_{i}}$ denote the polarisation
vectors. The Lorentz tensor structures appearing in the amplitude
${\mathcal{A}}_{\mu_{1}\mu_{2}\mu_{3}}$ were discussed in Ref. Kniehl:1990iva
. 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.
Chen:2019wyb . 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. Chen:2019wyb , we define the following
normalised linear polarisation vectors, where the frame that has been used is
shown in Fig. 2,

$\displaystyle\varepsilon^{\mu}_{x}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{x}\,\left(-s_{23}p^{\mu}_{1}-s_{13}p^{\mu}_{2}+s_{12}p^{\mu}_{3}\right)\,,$
$\displaystyle\varepsilon_{y}^{\mu}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{y}\,\left(\epsilon_{\mu_{1}\,\mu_{2}\,\mu_{3}}^{\mu}\,p^{\mu_{1}}_{1}\,p^{\mu_{2}}_{2}\,p^{\mu_{3}}_{3}\right),$
$\displaystyle\varepsilon^{\mu}_{T}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{T}\,\left(\left(-s_{23}(s_{13}+s_{23})+2m_{Z}^{2}s_{12}\right)p^{\mu}_{1}+\left(s_{13}(s_{13}+s_{23})-2m_{Z}^{2}s_{12}\right)p^{\mu}_{2}\right.$
$\displaystyle\left.+s_{12}(-s_{13}+s_{23})\,p^{\mu}_{3}\,\right)\,,$
$\displaystyle\varepsilon^{\mu}_{l}$ $\displaystyle=$
$\displaystyle\mathcal{N}_{l}\,\left(-2m_{Z}^{2}\left(p^{\mu}_{1}+p^{\mu}_{2}\right)+\left(s_{13}+s_{23}\right)p^{\mu}_{3}\right)\,,$
(7)

with $s_{12}\equiv 2p_{1}\cdot p_{2}=s\,,\,s_{13}\equiv 2p_{1}\cdot
p_{3}=s+t-m_{H}^{2}\,,\,s_{23}\equiv 2p_{2}\cdot p_{3}=m_{Z}^{2}-t$. The
normalisation factors $\mathcal{N}_{i}$ for $i\in\\{x,y,T,l\\}$, associated
with each of these space-like polarisation vectors in Eq. (2.1.1), can be
determined from their (negative) norm squares $\varepsilon^{2}_{i}$, i.e.
$\mathcal{N}_{i}=1/\sqrt{-\varepsilon^{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
$\left\\{\vec{p}_{1},\vec{p}_{2}\right\\}$ is taken as the $z$-axis and the
plane determined by $\left\\{\vec{p}_{1},\vec{p}_{3}\right\\}$ defines the
$x$-O-$z$ plane, as illustrated in Fig. 2.

Figure 2: The coordinate system in the center-of-mass frame of the two
incoming gluons.

Then the polarisation vector $\varepsilon_{x}$ is orthogonal to the beam axis
and lies within the $x$-O-$z$ plane, while $\varepsilon_{y}$ is perpendicular
to this plane. The vector $\varepsilon_{T}$ lies within the $x$-O-$z$ plane
but points to a direction orthogonal to $\vec{p}_{3}$, and $\varepsilon_{l}$,
the longitudinal polarisation of the $Z$-boson, points along the direction of
$\vec{p}_{3}$ in the center-of-mass frame.

Our projectors are given by products of three linear polarisation vectors,

$\varepsilon^{\mu_{1}}_{i}\,\varepsilon^{\mu_{2}}_{j}\,\varepsilon^{\mu_{3}}_{k}\,\text{
with $i\,,j\in\\{x,y\\}$ and $k\in\\{T,y,l\\}$},$ (8)

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
$\varepsilon_{y}$, are relevant for the process (1). This is because there are
only three linearly independent external momenta in
${\mathcal{A}}_{\mu_{1}\mu_{2}\mu_{3}}$ and all Lorentz structures appearing
in the amplitude contain a Levi-Civita symbol. Concretely, the six projectors
we use for ${\mathcal{A}}_{\mu_{1}\mu_{2}\mu_{3}}$ are given by

$\displaystyle\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{1}$ $\displaystyle=$
$\displaystyle\varepsilon^{\mu_{1}}_{x}\,\varepsilon^{\mu_{2}}_{x}\,\varepsilon^{\mu_{3}}_{y}\,,\quad\quad\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{2}=\varepsilon^{\mu_{1}}_{x}\,\varepsilon^{\mu_{2}}_{y}\,\varepsilon^{\mu_{3}}_{T},$
$\displaystyle\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{3}$ $\displaystyle=$
$\displaystyle\varepsilon^{\mu_{1}}_{x}\,\varepsilon^{\mu_{2}}_{y}\,\varepsilon^{\mu_{3}}_{l}\,,\quad\quad\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{4}=\varepsilon^{\mu_{1}}_{y}\,\varepsilon^{\mu_{2}}_{x}\,\varepsilon^{\mu_{3}}_{T},$
$\displaystyle\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{5}$ $\displaystyle=$
$\displaystyle\varepsilon^{\mu_{1}}_{y}\,\varepsilon^{\mu_{2}}_{x}\,\varepsilon^{\mu_{3}}_{l}\,,\quad\quad\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{6}=\varepsilon^{\mu_{1}}_{y}\,\varepsilon^{\mu_{2}}_{y}\,\varepsilon^{\mu_{3}}_{y}\,.$
(9)

The projections made with $\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{4}$ and
$\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{5}$ are related to those with
$\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{2}$ and
$\mathcal{P}^{\mu_{1}\mu_{2}\mu_{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.1.1) are understood to be
$D$-dimensional, and pairs of Levi-Civita symbols in
$\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{6}$ need to be contracted before being
used for the projection of the amplitude (see Ref. Chen:2019wyb  for a more
detailed discussion).

From the point of view of a form factor decomposition, the set of linear-
polarisation projectors in Eq. (2.1.1) 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 Chen:2019wyb ; Ahmed:2019udm ).

This projection defines six quantities

${\mathcal{A}}_{n}\equiv\mathcal{P}^{\mu_{1}\mu_{2}\mu_{3}}_{n}\,{\mathcal{A}}_{\mu_{1}\mu_{2}\mu_{3}}\quad(n=1,\ldots,6).$
(10)

The physical interpretation of the quantities ${\mathcal{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 (1) can be constructed from the linear ones, using
the relations

$\displaystyle\varepsilon^{\mu_{1}}_{\pm}(p_{1})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\varepsilon_{x}^{\mu_{1}}\pm
i\varepsilon_{y}^{\mu_{1}}\right)\,,$
$\displaystyle\varepsilon^{\mu_{2}}_{\pm}(p_{2})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\varepsilon_{x}^{\mu_{2}}\mp
i\varepsilon_{y}^{\mu_{2}}\right)\,,$
$\displaystyle\varepsilon^{\mu_{3}}_{\pm}(p_{3})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\varepsilon_{T}^{\mu_{3}}\pm
i\varepsilon_{y}^{\mu_{3}}\right)\,$ (11)

and the longitudinal polarisation of the $Z$-boson is given by
$\varepsilon^{\mu_{3}}_{l}$ in Eq. (2.1.1).

####  2.1.2 Kinematics

The kinematic invariants in terms of the scattering angle $\theta$ between the
beam axis and $\vec{p}_{3}$ in the centre-of-mass frame are given by

$\displaystyle t$
$\displaystyle=-\frac{1}{2}\left(s-m_{Z}^{2}-m_{H}^{2}+\cos\theta\sqrt{\lambda(s,m_{Z}^{2},m_{H}^{2})}\right)\;,$
$\displaystyle u$
$\displaystyle=-\frac{1}{2}\left(s-m_{Z}^{2}-m_{H}^{2}-\cos\theta\sqrt{\lambda(s,m_{Z}^{2},m_{H}^{2})}\right)\;,$
(12)

where
$\lambda(s,m_{Z}^{2},m_{H}^{2})=\left(s-(m_{Z}+m_{H})^{2}\right)\left(s-(m_{Z}-m_{H})^{2}\right)$
is the Källén function. Defining $m_{+}=m_{Z}+m_{H},m_{-}=m_{Z}-m_{H}$ and
$\beta_{ZH}^{2}=\left(1-\frac{m_{+}^{2}}{s}\right)\left(1-\frac{m_{-}^{2}}{s}\right)$,
we can also write

$\displaystyle t$
$\displaystyle=\frac{1}{2}\left(m_{Z}^{2}+m_{H}^{2}\right)-\frac{s}{2}\left(1+\beta_{ZH}\cos\theta\right)\;,$
$\displaystyle u$
$\displaystyle=\frac{1}{2}\left(m_{Z}^{2}+m_{H}^{2}\right)-\frac{s}{2}\left(1-\beta_{ZH}\cos\theta\right)\;.$
(13)

The limit $\beta_{ZH}\to 1$ corresponds to $s\gg m_{+}^{2},m_{-}^{2}$. In the
forward scattering region $\theta\to 0$, if in addition $\beta_{ZH}\to 1$ is
fulfilled, we have $|t|\gg|u|$. Analogously, for $\theta\to\pi$ and
$\beta_{ZH}\to 1$, the ratio $t/u$ is very small.

The virtual amplitude has a threshold at $s=4m_{t}^{2}$, therefore it is also
useful to define

$\displaystyle\beta_{t}=\frac{s-4m_{t}^{2}}{s+4m_{t}^{2}-2(m_{Z}+m_{H})^{2}}\;,$
(14)

such that $-1\leq\beta_{t}\leq 1$ with $\beta_{t}=-1$ for
$s=s_{min}=(m_{Z}+m_{H})^{2}$, $\beta_{t}=0$ at the top quark pair production
threshold $s=4m_{t}^{2}$ and $\beta_{t}\to 1$ for $s\to\infty$.

For $2\to 2$ kinematics, the transverse momentum $p_{T}$ of the Higgs- or
$Z$-boson (with $\vec{p}_{T}^{\,H}=-\vec{p}_{T}^{\,Z}$), fulfils $0\leq
p_{T}\leq\sqrt{\lambda}/(2\sqrt{s})$ and $s\,p_{T}^{2}=tu-m_{Z}^{2}m_{H}^{2}$.
Therefore, in the high energy limit, $p_{T}^{2}\to s\,(1-\cos^{2}\theta)$.

###  2.2 Diagram Generation, reduction and calculation of the master integrals

We generate the 1- and 2-loop diagrams using QGRAF Nogueira:1991ex . The
projectors described in Section 2.1.1 are applied and the Feynman rules are
inserted using FORM Vermaseren:2000nd ; Kuipers:2012rf ; Ruijl:2017dtg . For
the reduction we have defined eight integral families ${\cal F}_{i}$, five
planar $(i=1,\ldots,5)$ and three non-planar $(i=6,\ldots,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

$I=(\mu_{0}^{2})^{2\epsilon}\int\frac{\mbox{d}^{D}k_{1}}{i\pi^{\frac{D}{2}}}\int\frac{\mbox{d}^{D}k_{2}}{i\pi^{\frac{D}{2}}}\,N_{j_{1}}^{r_{1}}\ldots
N_{j_{t}}^{r_{t}}N_{j_{t+1}}^{-s_{1}}\ldots N_{j_{n}}^{-s_{n-t}}\;,$ (15)

where the $N_{j}$ denote genuine propagators of the form $1/(k^{2}-m^{2})$
with exponents $r_{i}\geq 1$ and $s_{i}\geq 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.

${\cal F}_{1}$ | ${\cal F}_{2}$ | ${\cal F}_{3}$
---|---|---
$k_{1}^{2}-m_{t}^{2}$ | $k_{1}^{2}-m_{t}^{2}$ | $k_{1}^{2}$
$k_{2}^{2}-m_{t}^{2}$ | $k_{2}^{2}-m_{t}^{2}$ | $(k_{1}-k_{2})^{2}-m_{t}^{2}$
$(k_{1}-k_{2})^{2}$ | $(k_{1}-k_{2})^{2}$ | $(k_{1}+p_{1})^{2}$
$(k_{1}+p_{1})^{2}-m_{t}^{2}$ | $(k_{1}+p_{1})^{2}-m_{t}^{2}$ | $(k_{2}+p_{1})^{2}-m_{t}^{2}$
$(k_{2}+p_{1})^{2}-m_{t}^{2}$ | $(k_{2}+p_{1})^{2}-m_{t}^{2}$ | $(k_{1}-p_{2})^{2}$
$(k_{1}-p_{2})^{2}-m_{t}^{2}$ | $(k_{1}-p_{3})^{2}-m_{t}^{2}$ | $(k_{2}-p_{2})^{2}-m_{t}^{2}$
$(k_{2}-p_{2})^{2}-m_{t}^{2}$ | $(k_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{2}-p_{2}-p_{3})^{2}-m_{t}^{2}$
$(k_{1}-p_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{1}-p_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{1}+p_{1}+p_{3})^{2}$
$(k_{2}-p_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{2}-p_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{2}+p_{1}-p_{2})^{2}$

${\cal F}_{4}$ | ${\cal F}_{5}$
---|---
$k_{1}^{2}-m_{t}^{2}$ | $k_{1}^{2}$
$k_{2}^{2}$ | $k_{2}^{2}-m_{t}^{2}$
$(k_{1}-k_{2})^{2}-m_{t}^{2}$ | $(k_{1}-k_{2})^{2}-m_{t}^{2}$
$(k_{1}+p_{1})^{2}-m_{t}^{2}$ | $(k_{1}+p_{1})^{2}$
$(k_{2}+p_{1})^{2}$ | $(k_{2}+p_{1})^{2}-m_{t}^{2}$
$(k_{1}-p_{2})^{2}-m_{t}^{2}$ | $(k_{1}-p_{3})^{2}$
$(k_{2}-p_{2})^{2}$ | $(k_{2}-p_{3})^{2}-m_{t}^{2}$
$(k_{1}-p_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{1}-p_{2}-p_{3})^{2}$
$(k_{2}-p_{2}-p_{3})^{2}$ | $(k_{2}-p_{2}-p_{3})^{2}-m_{t}^{2}$

Table 1: Planar integral families used for the reduction. ${\cal F}_{6}$ | ${\cal F}_{7}$ | ${\cal F}_{8}$
---|---|---
$k_{1}^{2}$ | $k_{1}^{2}-m_{t}^{2}$ | $k_{1}^{2}-m_{t}^{2}$
$k_{2}^{2}$ | $k_{2}^{2}-m_{t}^{2}$ | $k_{2}^{2}-m_{t}^{2}$
$(k_{1}+p_{1})^{2}$ | $(k_{1}-k_{2})^{2}$ | $(k_{1}-k_{2})^{2}$
$(k_{1}-p_{2})^{2}$ | $(k_{1}+p_{3})^{2}-m_{t}^{2}$ | $(k_{1}-p_{1}-p_{2}-p_{3})^{2}-m_{t}^{2}$
$(k_{2}+p_{1})^{2}-m_{t}^{2}$ | $(k_{2}+p_{3})^{2}-m_{t}^{2}$ | $(k_{2}-p_{1}-p_{2}-p_{3})^{2}-m_{t}^{2}$
$(k_{1}-k_{2})^{2}-m_{t}^{2}$ | $(k_{1}-p_{2})^{2}-m_{t}^{2}$ | $(k_{1}-p_{2})^{2}-m_{t}^{2}$
$(k_{2}-p_{2}-p_{3})^{2}-m_{t}^{2}$ | $(k_{2}-p_{2})^{2}-m_{t}^{2}$ | $(k_{2}-p_{1})^{2}-m_{t}^{2}$
$(k_{1}-k_{2}+p_{3})^{2}-m_{t}^{2}$ | $(k_{1}-k_{2}+p_{1})^{2}$ | $(k_{1}-k_{2}+p_{1})^{2}$
$(k_{1}-k_{2}+p_{2}+p_{3})^{2}$ | $(k_{2}-p_{1}-p_{2})^{2}-m_{t}^{2}$ | $(k_{2}-p_{1}-p_{2})^{2}-m_{t}^{2}$
Table 2: Non-planar integral families used for the reduction.

####  2.2.1 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.
vonManteuffel:2014qoa  to obtain a quasi-finite 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. Smirnov:2020quc ;
Usovitsch:2020jrk , 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 $\epsilon=(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 $\epsilon^{-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 particular
the size of the coefficient of the two-propagator (double-tadpole) integral
reduced from about 150 MB to less than 5 MB.

In order to perform the reduction to our chosen set of master integrals, we
utilize the Kira package Maierhoefer:2017hyi ; Klappert:2020nbg  in
combination with the rational function interpolation library FireFly
Klappert:2019emp ; Klappert:2020aqs . 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_{Z}^{2}/m_{t}^{2}=23/83$ and $m_{H}^{2}/m_{t}^{2}=12/23$, corresponding to
$m_{t}=173.21$ GeV, $m_{Z}=91.18$ GeV, and $m_{H}=125.1$ GeV. Our value for
$m_{W}$, which enters in electroweak couplings, amounts to
$m_{W}=80.379\,$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
Lee:2013mka  and Reduze vonManteuffel:2012np  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.

####  2.2.2 Evaluation of the two-loop amplitude

To evaluate the master integrals, we first apply sector decomposition as
implemented in the program pySecDec Borowka:2017idc ; Borowka:2018goh . For
integrals which diverge in the limit of 4 space-time dimensions
($\epsilon\rightarrow 0$), sector decomposition resolves singularities in the
regulator $\epsilon$ 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. Soper:1999xk ; Binoth:2005ff ; Nagy:2006xy ; Borowka:2012yc
.

In the present calculation we have selected a basis of quasi-finite integrals.
This means that poles in $\epsilon$ 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 Li:2015foa ; Borowka:2018goh . 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\rightarrow HH$
Borowka:2016ehy ; Borowka:2016ypz , $pp\rightarrow HJ$ Jones:2018hbb  and
$gg\rightarrow\gamma\gamma$ Chen:2019fla , 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.
Borowka:2016ypz  to minimise the total time taken to obtain a given relative
accuracy on the amplitude.

####  2.2.3 Renormalisation

We may expand each of the form factors, $\mathcal{A}_{i=1,\ldots,n}$, in the
bare strong coupling $a_{0}=\alpha_{0}/(4\pi)$ according to

$\mathcal{A}_{i}=a_{0}\mathcal{A}^{(0)}_{i}+a_{0}^{2}\mathcal{A}^{(1)}_{i}+\mathcal{O}(a_{0}).$
(16)

We renormalise the strong coupling in the $\overline{\mathrm{MS}}$ scheme with
the heavy quark loop in the gluon self-energy subtracted at zero momentum. The
heavy quark mass is renormalised in the $\mathrm{OS}$ scheme. The UV
renormalised amplitude is obtained via the replacement

$\displaystyle\mathcal{A}^{\mathrm{UV}}_{i}$
$\displaystyle=Z_{5}Z_{A}^{n_{g}/2}\mathcal{A}_{i}(a_{0}\rightarrow
a_{s}\,S_{\epsilon}^{-1}(\mu_{R}^{2}/\mu_{0}^{2})^{\epsilon}\,Z_{\alpha},m_{0}\rightarrow
m\,Z_{m},y_{t0}\rightarrow y_{t}\,Z_{m})$ $\displaystyle=$ $\displaystyle
a_{s}\,S_{\epsilon}^{-1}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{\epsilon}\mathcal{A}_{i}^{(0)}+a_{s}^{2}\,S_{\epsilon}^{-1}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{\epsilon}(\frac{n_{g}}{2}\delta
Z_{A}+\delta Z_{\alpha}+\delta Z_{5})\mathcal{A}_{i}^{(0)}$
$\displaystyle+a_{s}^{2}\,S_{\epsilon}^{-1}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{\epsilon}\delta
Z_{m}\mathcal{A}_{i}^{\mathrm{mct},(0)}+a_{s}^{2}\,S_{\epsilon}^{-2}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{2\epsilon}\mathcal{A}^{(1)}+\mathcal{O}(a_{s}^{3})$
(17)

Here, $n_{g}=2$ is the number of external gluon legs,
$a_{s}=\alpha_{s}(\mu_{R}^{2})/(4\pi)$ where $\alpha_{s}(\mu_{R}^{2})$ is the
renormalised coupling,
$S_{\epsilon}=(4\pi)^{\epsilon}\,e^{-\epsilon\gamma_{E}}$, and the
renormalisation constants are expanded according to $Z_{j}=1+a_{s}\delta
Z_{j}+\mathcal{O}(a_{s}^{2})\;(j=\alpha,A,m,5)$, with

$\displaystyle\delta Z_{\alpha}$
$\displaystyle=-\frac{1}{\epsilon}\beta_{0}+\delta
Z_{\alpha}^{\mathrm{hq}},\qquad\beta_{0}=\frac{11}{3}C_{A}-\frac{4}{3}T_{R}\,n_{F},$
(18) $\displaystyle\delta Z_{A}$ $\displaystyle=-\delta
Z_{\alpha}^{\mathrm{hq}}=\left(\frac{\mu_{R}^{2}}{m^{2}}\right)^{\epsilon}\left(-\frac{4}{3\epsilon}T_{R}\right),$
(19) $\displaystyle\delta Z_{m}$
$\displaystyle=\left(\frac{\mu_{R}^{2}}{m^{2}}\right)^{\epsilon}C_{F}\left(-\frac{3}{\epsilon}-4\right),$
(20) $\displaystyle\delta Z_{5}$ $\displaystyle=\begin{cases}\delta
Z_{5,ns}&=-4C_{F},\quad\textrm{diagrams involving a Z-top vertex}\\\ \delta
Z_{5,\rho}&=-8C_{F},\quad\textrm{diagrams involving a Goldstone-top
vertex}.\end{cases}$ (21)

We handle the $\gamma_{5}$ matrices which appear in this calculation using the
Larin scheme Larin:1993tq . According to this scheme an additional finite
renormalisation is required, it is denoted in the equations above by $\delta
Z_{5}$.

At the level of individual form factors, the procedure outlined corresponds to
the following relations

$\displaystyle\mathcal{A}^{\mathrm{UV}}_{i}$
$\displaystyle=a_{s}\mathcal{A}_{i}^{(0),\mathrm{UV}}+a_{s}^{2}\mathcal{A}_{i}^{(1),\mathrm{UV}}+\mathcal{O}(a_{s}^{3}),$
$\displaystyle\mathcal{A}_{i}^{(0),\mathrm{UV}}$
$\displaystyle=S_{\epsilon}^{-1}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{\epsilon}\mathcal{A}_{i}^{(0)},$
$\displaystyle\mathcal{A}_{i}^{(1),\mathrm{UV}}$
$\displaystyle=S_{\epsilon}^{-2}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{2\epsilon}\mathcal{A}_{i}^{(1)}+\left(-\,\frac{\beta_{0}}{\epsilon}+\delta
Z_{5}\right)\,\mathcal{A}_{i}^{(0),\mathrm{UV}}+S_{\epsilon}^{-1}\left(\frac{\mu_{R}^{2}}{\mu_{0}^{2}}\right)^{\epsilon}\delta
Z_{m}\,\mathcal{A}_{i}^{\mathrm{mct},(0)}.$ (22)

####  2.2.4 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. Catani:2013tia

$\displaystyle\mathcal{A}_{i}^{(0),\mathrm{fin}}$
$\displaystyle=\mathcal{A}_{i}^{(0),\mathrm{UV}},$ (23)
$\displaystyle\mathcal{A}_{i}^{(1),\mathrm{fin}}$
$\displaystyle=\mathcal{A}_{i}^{(1),\mathrm{UV}}-I_{1}\mathcal{A}_{i}^{(0),\mathrm{UV}},$
(24)

with

$\displaystyle I_{1}$
$\displaystyle=I_{1}^{\mathrm{soft}}+I_{1}^{\mathrm{coll}},$ (25)
$\displaystyle I_{1}^{\mathrm{soft}}$
$\displaystyle=-\frac{e^{\epsilon\gamma_{E}}}{\Gamma(1-\epsilon)}\left(\frac{\mu_{R}^{2}}{s}\right)^{\epsilon}\left(\frac{1}{\epsilon^{2}}+\frac{i\pi}{\epsilon}\right)2C_{A},$
(26) $\displaystyle I_{1}^{\mathrm{coll}}$
$\displaystyle=-\frac{\beta_{0}}{\epsilon}\left(\frac{\mu_{R}^{2}}{s}\right)^{\epsilon}.$
(27)

We present results at renormalisation scale $\mu_{R}^{2}=s$ and after changing
to the helicity basis, we define

$\displaystyle\mathcal{B}$ $\displaystyle=\frac{8\,T_{R}^{2}}{(2\cdot
8)^{2}}\sum_{\mathrm{\\{\lambda_{i}\\}}}\mathcal{A}_{\\{\lambda_{i}\\}}^{(0),\mathrm{fin}}\mathcal{A}_{\\{\lambda_{i}\\}}^{\star(0),\mathrm{fin}},$
(28) $\displaystyle\mathcal{V}$ $\displaystyle=\frac{8\,T_{R}^{2}}{(2\cdot
8)^{2}}\sum_{\mathrm{{\\{\lambda_{i}\\}}}}\left(\mathcal{A}_{{\\{\lambda_{i}\\}}}^{(0),\mathrm{fin}}\mathcal{A}_{{\\{\lambda_{i}\\}}}^{\star(1),\mathrm{fin}}+\mathcal{A}_{{\\{\lambda_{i}\\}}}^{(1),\mathrm{fin}}\mathcal{A}_{{\\{\lambda_{i}\\}}}^{\star(0),\mathrm{fin}}\right)$
(29)

for the square of the amplitude. We set the overall factor of
$e^{2}=4\pi\alpha$ to one in our amplitude.

##  3 Results

###  3.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 Hasselhuhn:2016rqt . The
virtual contributions $\widetilde{\mathcal{V}}_{n}$ are expanded up to order
$1/m_{t}^{2n}$ and reweighted with the exact Born term $\mathcal{B}$. The one-
particle-reducible double triangle contributions $\mathcal{V}^{\mathrm{red}}$
are included with full top quark mass dependence.

$\mathcal{V}_{n}=\frac{\mathcal{B}}{\mathcal{B}_{n}}\widetilde{\mathcal{V}}_{n}+\mathcal{V}^{\mathrm{red}}$
(30)

In Fig. 3 the ratio of expanded to full virtual contribution
$\mathcal{V}_{n}/\mathcal{V}$ is shown for expansion orders $0\leq n\leq 4$
for a fixed scattering angle $\mathrm{cos}(\theta)\approx 0.052$. For energies
close to the production threshold at $\beta_{t}=-1$
($s=\left(m_{Z}+m_{H}\right)^{2}$) the expanded result approximates the exact
calculation well with a ratio $\mathcal{V}_{4}/\mathcal{V}\approx 0.9989$,
while the agreement worsens closer to the top quark pair threshold at
$\beta_{t}=0$ ($s=4m_{t}^{2}$), where the large $m_{t}$-expansion is expected
to break down.

Figure 3: Comparison of exact virtual contributions $\mathcal{V}$ with the
expansion in the large top quark mass limit $\mathcal{V}_{n}$ from
Hasselhuhn:2016rqt  for fixed scattering angle $\mathrm{cos}(\theta)\approx
0.052$. The range of the horizontal axis $-1\leq\beta_{t}\leq 0$ corresponds
to the energy range $\left(m_{Z}+m_{H}\right)^{2}\leq s\leq 4m_{t}^{2}$
between production threshold and top quark pair threshold.

In addition, we have compared our results to a recent calculation in the high
energy limit Davies:2020drs  and found agreement to the extent expected from
previous comparisons for the process $pp\to HH$ performed in Ref.
Davies:2019dfy .

###  3.2 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\beta_{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.

Figure 4: Dependence of the leading order (left) and virtual contribution
(right) on the parameters $\beta_{t}$ and $\cos\theta$, summing over all
polarisation.  Figure 5: Dependence of the ratio $\mathcal{V}/\mathcal{B}$ on
the parameters $\beta_{t}$ and $\cos\theta$, summing over all polarisations.

In the following we will show results for five helicity amplitudes in the
$\beta_{t}$–$\cos\theta$–plane. The remaining amplitudes can be obtained from
overall helicity flips as well as crossing between the two initial gluons.
Bose-symmetry and the behaviour under parity transformations imply

$\displaystyle{\mathcal{A}}_{\lambda_{1},\lambda_{2},\lambda_{z}}(t,u)$
$\displaystyle=(-1)^{\lambda_{z}}{\mathcal{A}}_{\lambda_{2},\lambda_{1},\lambda_{z}}(u,t)\;,$
(31) $\displaystyle{\mathcal{A}}_{\lambda_{1},\lambda_{2},\lambda_{z}}(t,u)$
$\displaystyle=-{\mathcal{A}}_{-\lambda_{1},-\lambda_{2},-\lambda_{z}}(t,u)\;.$

To understand the qualitative behaviour, we can expand the amplitude in terms
of Wigner $d$-functions,

$\displaystyle{\mathcal{A}}_{\\{\lambda_{i}\\}}(\theta)=\sum_{J=0}^{\infty}(2J+1){\mathcal{A}}^{J}_{\\{\lambda_{i}\\}}d^{J}_{s_{i},s_{f}}(\theta)\;,$
(32)

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\rightarrow ZH$, the initial state has total spin $s_{i}=0$ for equal
helicities of the initial state gluons, $(\lambda_{1},\lambda_{2})=(+,+)$ or
$(-,-)$, while for the case $(\lambda_{1},\lambda_{2})=(+,-)$ or $(-,+)$, the
initial state has total spin $s_{i}=2$. Therefore the amplitude ${\cal
A}_{++0}$ is dominated by the partial wave $d^{0}_{00}(\theta)$ 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
${\cal A}_{++0}$ in $\cos\theta$. 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\theta$. Note that eq.
(31) implies that ${\cal A}_{++0}$ is symmetric under exchange of
$t\leftrightarrow u$ and therefore is symmetric in $\cos\theta$. As ${\cal
A}_{++0}$ is composed of partial waves $d^{J}_{00}(\theta)$, which are even in
$\cos\theta$ for even $J$, no partial wave components with odd $J$ can
contribute to ${\cal A}_{++0}$.

From Fig. 6, we further observe that the helicity amplitudes with the
polarisations $\epsilon^{Z}_{\pm}$ are suppressed compared to those with a
longitudinally polarised $Z$-boson. The amplitudes ${\cal A}_{++\pm}$ are
antisymmetric under exchange of $t$ and $u$, so antisymmetric in $\cos\theta$.
The contributions with $J=1$, being proportional to $\sin\theta$, therefore do
not occur in ${\cal A}_{++\pm}$. The d-wave contributions $d_{0,\pm
1}^{2}(\theta)$ are proportional to $\pm\cos\theta\sin\theta$, which vanish at
$\cos\theta=0,\pm 1$, a behaviour that can be observed in ${\cal 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 ${\cal A}_{++-}$ is about five orders of magnitude smaller
than ${\cal A}_{++0}$, and ${\cal A}_{+++}$ is also suppressed, therefore the
amplitudes ${\cal A}_{++\pm}$ give a very minor contribution in the sum of all
polarisation configurations.

For different helicities of the initial state gluons,
$(\lambda_{1},\lambda_{2})=(+,-)$ or $(-,+)$, shown in Fig. 7, the initial
state has total spin $s_{i}=2$, such that the partial wave contributions start
from $J=2$. Therefore the amplitude ${\cal A}_{+-0}$ is much smaller than
${\cal A}_{++0}$. This is also reflected in the fact that the amplitudes are
basically zero except at very high energies. The amplitude ${\cal A}_{+-0}$
has no contribution from $d_{2,0}^{J}(\theta)$ with even $J$ as it is
antisymmetric in $\cos\theta$, therefore its leading partial wave is given by
$d_{2,0}^{3}(\theta)\sim\cos\theta\sin^{2}\theta$, vanishing at
$\cos\theta=\pm 1$ and $\cos\theta=0$. The amplitudes ${\cal A}_{+-\pm}$ have
their leading partial waves given by $d_{2,\pm
1}^{2}(\theta)\sim\sin\theta(1\pm\cos\theta)$. Consequently ${\cal A}_{+-+}$
is highly suppressed in the backward direction, while ${\cal A}_{+--}$ is
highly suppressed in the forward direction.

Figure 6: Dependence of the leading order (left) and virtual contribution
(right) on the parameters $\beta_{t}$ and $\cos(\theta)$ for the individual
helicity amplitudes with $s_{i}=0$.

Figure 7: Dependence of the leading order (left) and virtual contribution
(right) on the parameters $\beta_{t}$ and $\cos(\theta)$ for the individual
helicity amplitudes with $s_{i}=2$.

In Table 3 we list numerical results of the amplitude to facilitate
comparisons with our results.

$s/m_{t}^{2}$ | $t/m_{t}^{2}$ | $\mathcal{B}$ | $\mathcal{V}$
---|---|---|---
1.707133657190554 | -0.441203767016323 | 0.4412287 | 35.429092(6)
3.876056604162662 | -1.616287256345735 | 37.2496999 | 4339.045(1)
4.130574250302561 | -1.750372271104745 | 66.3224970 | 6912.361(3)
4.130574250302561 | -2.595461551488002 | 67.1908198 | 6981.09(2)
134.5142052093564 | -70.34125943305149 | 4.1920928 | -153.9(4)
134.5142052093564 | -105.1770655376327 | 14.7405104 | 527(4)
Table 3: Numerical results for various phase space points at the scale
$\mu_{R}^{2}=s$. The number in parentheses gives the numerical uncertainty on
the last digit of the virtual amplitude.

##  4 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 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 Chen:2019wyb . The reduction to master integrals
has been performed with the program Kira Maierhoefer:2017hyi ;
Klappert:2020nbg  in combination with the rational function interpolation
library FireFly Klappert:2019emp ; Klappert:2020aqs , using in addition
LiteRed Lee:2013mka  and Reduze vonManteuffel:2012np  to obtain dimensional
recurrence relations. The master integrals have been calculated using pySecDec
Borowka:2017idc ; Borowka:2018goh . 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.

## Acknowledgements

We would like to thank Joshua Davies, Go Mishima and Matthias Steinhauser for
the comparison of results prior to publication. We also would like to thank
Tom Zirke for conversation about the large mass expansion and Ramona Gröber
for useful discussions. This research was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762
- TRR 257, by the COST Action CA16201 (‘Particleface’) of the European Union
and by the Swiss National Science Foundation (SNF) under grant number
200020-175595. MK acknowledges support by the Forschungskredit of the
University of Zurich, grant no. FK-19-102. The research of JS was supported by
the European Union through the ERC Advanced Grant MC@NNLO (340983). SJ is
supported by a Royal Society University Research Fellowship (Grant
URF/R1/201268). We gratefully acknowledge resources provided by the Max Planck
Computing and Data Facility (MPCDF). Some of the calculations were performed
with computing resources granted by RWTH Aachen University under project
rwth0541.

## References

  * (1) ATLAS collaboration, M. Aaboud et al., _Observation of $H\rightarrow b\bar{b}$ decays and $VH$ production with the ATLAS detector_, _Phys. Lett. B_ 786 (2018) 59 [1808.08238].
  * (2) CMS collaboration, A. M. Sirunyan et al., _Observation of Higgs boson decay to bottom quarks_ , _Phys. Rev. Lett._ 121 (2018) 121801 [1808.08242].
  * (3) K. Mimasu, V. Sanz and C. Williams, _Higher Order QCD predictions for Associated Higgs production with anomalous couplings to gauge bosons_ , _JHEP_ 08 (2016) 039 [1512.02572].
  * (4) A. Greljo, G. Isidori, J. M. Lindert, D. Marzocca and H. Zhang, _Electroweak Higgs production with HiggsPO at NLO QCD_ , _Eur. Phys. J. C_ 77 (2017) 838 [1710.04143].
  * (5) S. Alioli, W. Dekens, M. Girard and E. Mereghetti, _NLO QCD corrections to SM-EFT dilepton and electroweak Higgs boson production, matched to parton shower in POWHEG_ , _JHEP_ 08 (2018) 205 [1804.07407].
  * (6) S. Banerjee, C. Englert, R. S. Gupta and M. Spannowsky, _Probing Electroweak Precision Physics via boosted Higgs-strahlung at the LHC_ , _Phys. Rev. D_ 98 (2018) 095012 [1807.01796].
  * (7) S. Banerjee, R. S. Gupta, J. Y. Reiness and M. Spannowsky, _Resolving the tensor structure of the Higgs coupling to $Z$-bosons via Higgs-strahlung_, _Phys. Rev. D_ 100 (2019) 115004 [1905.02728].
  * (8) J. Baglio, S. Dawson, S. Homiller, S. D. Lane and I. M. Lewis, _Validity of standard model EFT studies of VH and VV production at NLO_ , _Phys. Rev. D_ 101 (2020) 115004 [2003.07862].
  * (9) K. Rao, S. D. Rindani and P. Sarmah, _Study of anomalous gauge-Higgs couplings using Z boson polarization at LHC_ , _Nucl. Phys. B_ 964 (2021) 115317 [2009.00980].
  * (10) D. A. Dicus and C. Kao, _Higgs Boson - $Z^{0}$ Production From Gluon Fusion_, _Phys. Rev. D_ 38 (1988) 1008.
  * (11) B. A. Kniehl, _Associated Production of Higgs and Z Bosons From Gluon Fusion in Hadron Collisions_ , _Phys. Rev. D_ 42 (1990) 2253.
  * (12) R. V. Harlander, S. Liebler and T. Zirke, _Higgs Strahlung at the Large Hadron Collider in the 2-Higgs-Doublet Model_ , _JHEP_ 02 (2014) 023 [1307.8122].
  * (13) C. Englert, M. McCullough and M. Spannowsky, _Gluon-initiated associated production boosts Higgs physics_ , _Phys. Rev. D_ 89 (2014) 013013 [1310.4828].
  * (14) ATLAS collaboration, G. Aad et al., _Measurement of the associated production of a Higgs boson decaying into $b$-quarks with a vector boson at high transverse momentum in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, 2008.02508.
  * (15) CMS collaboration, A. M. Sirunyan et al., _Search for the associated production of the Higgs boson and a vector boson in proton-proton collisions at $\sqrt{s}=$ 13 TeV via Higgs boson decays to $\tau$ leptons_, _JHEP_ 06 (2019) 093 [1809.03590].
  * (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_ 02 (2013) 078 [1211.5015].
  * (17) A. Hasselhuhn, T. Luthe and M. Steinhauser, _On top quark mass effects to $gg\to ZH$ at NLO_, _JHEP_ 01 (2017) 073 [1611.05881].
  * (18) R. V. Harlander, A. Kulesza, V. Theeuwes and T. Zirke, _Soft gluon resummation for gluon-induced Higgs Strahlung_ , _JHEP_ 11 (2014) 082 [1410.0217].
  * (19) R. Harlander, J. Klappert, C. Pandini and A. Papaefstathiou, _Exploiting the WH/ZH symmetry in the search for New Physics_ , _Eur. Phys. J. C_ 78 (2018) 760 [1804.02299].
  * (20) O. Brein, A. Djouadi and R. Harlander, _NNLO QCD corrections to the Higgs-strahlung processes at hadron colliders_ , _Phys. Lett. B_ 579 (2004) 149 [hep-ph/0307206].
  * (21) O. Brein, R. Harlander, M. Wiesemann and T. Zirke, _Top-Quark Mediated Effects in Hadronic Higgs-Strahlung_ , _Eur. Phys. J. C_ 72 (2012) 1868 [1111.0761].
  * (22) O. Brein, R. V. Harlander and T. J. Zirke, _vh@nnlo - Higgs Strahlung at hadron colliders_ , _Comput. Phys. Commun._ 184 (2013) 998 [1210.5347].
  * (23) R. V. Harlander, J. Klappert, S. Liebler and L. Simon, _vh@nnlo-v2: New physics in Higgs Strahlung_ , _JHEP_ 05 (2018) 089 [1802.04817].
  * (24) M. Ciccolini, S. Dittmaier and M. Krämer, _Electroweak radiative corrections to associated WH and ZH production at hadron colliders_ , _Phys. Rev. D_ 68 (2003) 073003 [hep-ph/0306234].
  * (25) A. Denner, S. Dittmaier, S. Kallweit and A. Mück, _Electroweak corrections to Higgs-strahlung off W/Z bosons at the Tevatron and the LHC with HAWK_ , _JHEP_ 03 (2012) 075 [1112.5142].
  * (26) A. Denner, S. Dittmaier, S. Kallweit and A. Mück, _HAWK 2.0: A Monte Carlo program for Higgs production in vector-boson fusion and Higgs strahlung at hadron colliders_ , _Comput. Phys. Commun._ 195 (2015) 161 [1412.5390].
  * (27) F. Granata, J. M. Lindert, C. Oleari and S. Pozzorini, _NLO QCD+EW predictions for HV and HV +jet production including parton-shower effects_ , _JHEP_ 09 (2017) 012 [1706.03522].
  * (28) P. Obul, S. Dulat, T.-J. Hou, A. Tursun and N. Yalkun, _Next-to-leading order QCD and electroweak corrections to Higgs-strahlung processes at the LHC_ , _Chin. Phys. C_ 42 (2018) 093105 [1801.06851].
  * (29) G. Ferrera, M. Grazzini and F. Tramontano, _Associated ZH production at hadron colliders: the fully differential NNLO QCD calculation_ , _Phys. Lett. B_ 740 (2015) 51 [1407.4747].
  * (30) J. M. Campbell, R. K. Ellis and C. Williams, _Associated production of a Higgs boson at NNLO_ , _JHEP_ 06 (2016) 179 [1601.00658].
  * (31) G. Ferrera, G. Somogyi and F. Tramontano, _Associated production of a Higgs boson decaying into bottom quarks at the LHC in full NNLO QCD_ , _Phys. Lett. B_ 780 (2018) 346 [1705.10304].
  * (32) R. Gauld, A. Gehrmann-De Ridder, E. Glover, A. Huss and I. Majer, _Associated production of a Higgs boson decaying into bottom quarks and a weak vector boson decaying leptonically at NNLO in QCD_ , _JHEP_ 10 (2019) 002 [1907.05836].
  * (33) G. Luisoni, P. Nason, C. Oleari and F. Tramontano, _$HW^{\pm}$ /HZ + 0 and 1 jet at NLO with the POWHEG BOX interfaced to GoSam and their merging within MiNLO_, _JHEP_ 10 (2013) 083 [1306.2542].
  * (34) D. Goncalves, F. Krauss, S. Kuttimalai and P. Maierhöfer, _Higgs-Strahlung: Merging the NLO Drell-Yan and Loop-Induced 0+1 jet Multiplicities_ , _Phys. Rev. D_ 92 (2015) 073006 [1509.01597].
  * (35) B. Hespel, F. Maltoni and E. Vryonidou, _Higgs and Z boson associated production via gluon fusion in the SM and the 2HDM_ , _JHEP_ 06 (2015) 065 [1503.01656].
  * (36) W. Astill, W. Bizoń, E. Re and G. Zanderighi, _NNLOPS accurate associated HZ production with $H\to b\overline{b}$ decay at NLO_, _JHEP_ 11 (2018) 157 [1804.08141].
  * (37) W. Bizoń, E. Re and G. Zanderighi, _NNLOPS description of the $H\to b\overline{b}$ decay with MiNLO_, _JHEP_ 06 (2020) 006 [1912.09982].
  * (38) S. Alioli, A. Broggio, S. Kallweit, M. A. Lim and L. Rottoli, _Higgsstrahlung at NNLL ${}^{\prime}+$NNLO matched to parton showers in GENEVA_, _Phys. Rev. D_ 100 (2019) 096016 [1909.02026].
  * (39) M. Kumar, M. Mandal and V. Ravindran, _Associated production of Higgs boson with vector boson at threshold N 3LO in QCD_, _JHEP_ 03 (2015) 037 [1412.3357].
  * (40) S. Dawson, T. Han, W. Lai, A. Leibovich and I. Lewis, _Resummation Effects in Vector-Boson and Higgs Associated Production_ , _Phys. Rev. D_ 86 (2012) 074007 [1207.4207].
  * (41) T. Ahmed, A. Ajjath, L. Chen, P. K. Dhani, P. Mukherjee and V. Ravindran, _Polarised Amplitudes and Soft-Virtual Cross Sections for $b\bar{b}\rightarrow ZH$ at NNLO in QCD_, _JHEP_ 01 (2020) 030 [1910.06347].
  * (42) T. Ahmed, W. Bernreuther, L. Chen and M. Czakon, _Polarized $q\bar{q}\rightarrow Z+$Higgs amplitudes at two loops in QCD: the interplay between vector and axial vector form factors and a pitfall in applying a non-anticommuting $\gamma_{5}$_, _JHEP_ 07 (2020) 159 [2004.13753].
  * (43) G. ’t Hooft and M. Veltman, _Regularization and Renormalization of Gauge Fields_ , _Nucl. Phys. B_ 44 (1972) 189.
  * (44) P. Breitenlohner and D. Maison, _Dimensional Renormalization and the Action Principle_ , _Commun. Math. Phys._ 52 (1977) 11.
  * (45) S. A. Larin and J. A. M. Vermaseren, _The $\alpha_{s}^{3}$ corrections to the Bjorken sum rule for polarized 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_ 303 (1993) 113 [hep-ph/9302240].
  * (47) L. Chen, _A prescription for projectors to compute helicity amplitudes in D dimensions_ , 1904.00705.
  * (48) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (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._ 184 (2013) 1453 [1203.6543].
  * (51) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (52) A. von Manteuffel, E. Panzer and R. M. Schabinger, _A quasi-finite basis for multi-loop Feynman integrals_ , _JHEP_ 02 (2015) 120 [1411.7392].
  * (53) A. Smirnov and V. Smirnov, _How to choose master integrals_ , _Nucl. Phys. B_ 960 (2020) 115213 [2002.08042].
  * (54) J. Usovitsch, _Factorization of denominators in integration-by-parts reductions_ , 2002.08173.
  * (55) P. Maierhöfer, J. Usovitsch and P. Uwer, _Kira—A Feynman integral reduction program_ , _Comput. Phys. Commun._ 230 (2018) 99 [1705.05610].
  * (56) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral Reduction with Kira 2.0 and Finite Field Methods_ , 2008.06494.
  * (57) J. Klappert and F. Lange, _Reconstructing rational functions with FireFly_ , _Comput. Phys. Commun._ 247 (2020) 106951 [1904.00009].
  * (58) J. Klappert, S. Y. Klein and F. Lange, _Interpolation of Dense and Sparse Rational Functions and other Improvements in FireFly_, 2004.01463.
  * (59) R. N. Lee, _LiteRed 1.4: a powerful tool for reduction of multiloop integrals_ , _J. Phys. Conf. Ser._ 523 (2014) 012059 [1310.1145].
  * (60) A. von Manteuffel and C. Studerus, _Reduze 2 - Distributed Feynman Integral Reduction_ , 1201.4330.
  * (61) S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner, J. Schlenk et al., _pySecDec: a toolbox for the numerical evaluation of multi-scale integrals_ , _Comput. Phys. Commun._ 222 (2018) 313 [1703.09692].
  * (62) S. Borowka, G. Heinrich, S. Jahn, S. Jones, M. Kerner and J. Schlenk, _A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec_ , _Comput. Phys. Commun._ 240 (2019) 120 [1811.11720].
  * (63) D. E. Soper, _Techniques for QCD calculations by numerical integration_ , _Phys. Rev. D_ 62 (2000) 014009 [hep-ph/9910292].
  * (64) T. Binoth, J. Guillet, G. Heinrich, E. Pilon and C. Schubert, _An Algebraic/numerical formalism for one-loop multi-leg amplitudes_ , _JHEP_ 10 (2005) 015 [hep-ph/0504267].
  * (65) Z. Nagy and D. E. Soper, _Numerical integration of one-loop Feynman diagrams for N-photon amplitudes_ , _Phys. Rev. D_ 74 (2006) 093006 [hep-ph/0610028].
  * (66) S. Borowka, J. Carter and G. Heinrich, _Numerical Evaluation of Multi-Loop Integrals for Arbitrary Kinematics with SecDec 2.0_ , _Comput. Phys. Commun._ 184 (2013) 396 [1204.4152].
  * (67) Z. Li, J. Wang, Q.-S. Yan and X. Zhao, _Efficient numerical evaluation of Feynman integrals_ , _Chin. Phys. C_ 40 (2016) 033103 [1508.02512].
  * (68) S. Borowka, N. Greiner, G. Heinrich, S. Jones, M. Kerner, J. Schlenk et al., _Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 117 (2016) 012001 [1604.06447].
  * (69) S. Borowka, N. Greiner, G. Heinrich, S. Jones, M. Kerner, J. Schlenk et al., _Full top quark mass dependence in Higgs boson pair production at NLO_ , _JHEP_ 10 (2016) 107 [1608.04798].
  * (70) S. Jones, M. Kerner and G. Luisoni, _Next-to-Leading-Order QCD Corrections to Higgs Boson Plus Jet Production with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 120 (2018) 162001 [1802.00349].
  * (71) L. Chen, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk et al., _Photon pair production in gluon fusion: Top quark effects at NLO with threshold matching_ , _JHEP_ 04 (2020) 115 [1911.09314].
  * (72) S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, _Universality of transverse-momentum resummation and hard factors at the NNLO_ , _Nucl. Phys. B_ 881 (2014) 414 [1311.1654].
  * (73) J. Davies, G. Mishima and M. Steinhauser, _Virtual corrections to $gg\to ZH$ in the high-energy and large-$m_{t}$ limits_, _JHEP_ 03 (2021) 034 [2011.12314].
  * (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 and high-energy expansion_ , _JHEP_ 11 (2019) 024 [1907.06408].