COMETA-2024-06

TTP24-010

Virtual QCD corrections to $gg\to ZZ$: top-quark loops from a transverse-
momentum expansion

Giuseppe Degrassia***email<EMAIL_ADDRESS>Ramona
Gröberb†††email<EMAIL_ADDRESS>Marco Vittic,d‡‡‡email:
<EMAIL_ADDRESS>
(a) Dipartimento di Matematica e Fisica, Università di Roma Tre and
INFN, sezione di Roma Tre, I-00146 Rome, Italy
(b) Dipartimento di Fisica e Astronomia ’G. Galilei’, Università di Padova and
INFN, sezione di Padova, I-35131 Padova, Italy
(c) Institute for Theoretical Particle Physics, Karlsruhe Institute of
Technology (KIT), D-76131 Karlsruhe, Germany
(d) Institute for Astroparticle Physics, Karlsruhe Institute of Technology
(KIT), D-76344 Eggenstein-Leopoldshafen, Germany

We present the virtual corrections due to the top-quark loops for the process
$gg\to ZZ$ at next-to-leading order in QCD. The associated two-loop box
diagrams are computed using a small-transverse-momentum expansion. Our results
are then merged with those available in the complementary energy region,
obtained via a high-energy expansion, in order to provide an analytic result
that is valid in the whole phase space. The results presented allow for an
efficient modelling of the signal–background interference as well as the
irreducible background in off-shell Higgs production.

##  1 Introduction

The quest for the determination of the properties of the Higgs boson, after
its discovery at the LHC, has been very rewarding. While the main Higgs
production channels have been established and new challenging decay modes are
being studied (see e.g. refs. [1, 2]), a large variety of measurements have
entered the precision domain, requiring improved predictions for Higgs-related
processes within the Standard Model (SM) [3, 4]. A significant role in this
program is played by the production of a pair of $Z$ bosons, $pp\rightarrow
ZZ$, which has been one of the discovery channels of the Higgs and today is
important both as a probe of the electroweak (EW) symmetry breaking and for
precision Higgs physics [5, 6].

In the theoretical SM prediction of $pp\rightarrow ZZ$ production two partonic
sub-processes have to be considered. The first one is quark-antiquark
annihilation, $q\bar{q}\rightarrow ZZ$, which gives the largest contribution
to the hadronic cross section. The leading-order (LO) amplitude for this
channel is related to purely EW tree-level diagrams [7], and corrections
through next-to-next-to-leading order (NNLO) in QCD [8, 9, 10, 11, 12, 13] and
through next-to-leading order (NLO) in the EW theory [14, 15, 16] are
available. The second partonic contribution, which is the main topic of this
paper, comes from the gluon-initiated channel, $gg\rightarrow ZZ$. The gluon-
initiated channel at LO is associated to one-loop diagrams, which contribute
as an important $\mathcal{O}(\alpha_{s}^{2})$ correction, accounting for about
10% of the hadronic cross section at NNLO for $\sqrt{s}=13$ TeV [17].

The one-loop diagrams for $gg\rightarrow ZZ$ at LO have been computed for the
first time in refs. [18, 19]. They feature two topologies: the triangles (see
fig. 1(a)) are associated to Higgs production via the sub-process
$gg\rightarrow H\rightarrow ZZ$, while the box diagrams (fig. 1(b)) are
related to the process of non-resonant (a.k.a.  _continuum_) $ZZ$ production,
which constitutes an irreducible background in experimental Higgs searches.
Moreover, continuum production plays a relevant role in the indirect
determination of the Higgs total decay width, $\Gamma_{H}$ [20, 21]. Indeed,
in refs. [22, 23] it has been suggested that upper limits on $\Gamma_{H}$ can
be obtained from the investigation of the invariant-mass distribution of the
two $Z$ bosons system, $M_{\scriptscriptstyle ZZ}$, away from the region where
the Higgs is produced on shell111Off-shell Higgs production can also provide
an excellent probe of new physics, for instance of light-quark Yukawa
couplings [29], modified trilinear Higgs self-coupling [30] or Higgs portal
models [31, 32] [24]. Under the assumption that no beyond-SM physics spoils
the known correlation between on- and off-shell amplitudes [25, 26], the
latest experimental measurements could exclude, for the first time, deviations
in $\Gamma_{H}$ below $\mathcal{O}(100\%)$ from its SM prediction [27, 28].

In this situation, a good theoretical control over the destructive
interference between the Higgs-mediated and continuum amplitudes is crucial.
This motivates the inclusion of NLO QCD corrections to the gluon-initiated
channel. The NLO corrections to the Higgs-mediated diagrams are known exactly,
adapting results for the production of a single Higgs with virtuality
$M_{\scriptscriptstyle ZZ}$ [33, 34, 35]. Concerning the continuum term, exact
analytic results for the related two-loop box integrals have been obtained in
the case of loops of massless quarks [36, 37]. We remark that, at the level of
the inclusive $gg\rightarrow ZZ$ cross section, the contribution from light
quarks is the dominant one, and furthermore it constitutes more than 50% of
the $\mathcal{O}(\alpha_{s}^{2})$ corrections to $ZZ$ production [10, 17].

The top-quark contribution to continuum $gg\rightarrow ZZ$ starts to become
relevant for invariant masses in the range $M_{\scriptscriptstyle
ZZ}>2m_{\scriptscriptstyle Z}$ and it is expected to be substantial in the
region where $M_{\scriptscriptstyle ZZ}$ is larger than twice the mass of the
top, $m_{t}$. At present, the top-quark contribution is not known in full
analytic form, as the scale associated to the mass of the heavy quark
complicates the calculation of the two-loop box integrals. Several approximate
analytic evaluations of this contribution have been discussed in the
literature. The large-mass expansion (LME) has been used in refs. [38, 39, 40]
to obtain reliable predictions in the region $2\,m_{\scriptscriptstyle
Z}<M_{\scriptscriptstyle ZZ}\lesssim 2\,m_{t}$. In ref. [40] an improvement of
the LME results by means of conformal mapping and Padé approximants has also
been studied. A further improvement of the LME via an expansion around the top
threshold has been presented in ref. [41]. Analytic predictions that are
reliable for large invariant masses have been obtained using the so-called
High-Energy (HE) expansion [42]. Still, the latter approach is expected to
break down in the region $M_{\scriptscriptstyle ZZ}\lesssim 750$ GeV, and in
ref. [42] Padé approximants were used in order to improve the expansion.
Finally, exact results based on numerical approaches have been obtained in
refs. [43, 44], showing good agreement with ref. [42] in the ranges of small
or large invariant masses. Very recently, the numerical results of ref. [43]
have been used to compute the full NLO QCD corrections to $gg\to ZZ$ [45].

In this paper, we present the computation of the contribution from top-quark
loops to the virtual NLO QCD corrections to $gg\rightarrow ZZ$. In the
calculation we employ an expansion in the transverse momentum of the $Z$,
$p_{\scriptscriptstyle T}$, following refs. [46, 47]. Our first goal is to
provide an accurate approximation of the virtual QCD corrections that can be
valid in the invariant-mass region that so far has not been covered by any of
the analytic approaches discussed above, namely the region $350\lesssim
M_{\scriptscriptstyle ZZ}\lesssim 750$ GeV. Since the $p_{\scriptscriptstyle
T}$ expansion “contains” the LME result, the region of validity of our
approach is actually given by $2\,m_{\scriptscriptstyle
Z}<M_{\scriptscriptstyle ZZ}\lesssim 750$ GeV.

Furthermore, it has been previously shown that an expansion in the forward
kinematics [48, 49], like the $p_{\scriptscriptstyle T}$ expansion, and the HE
expansion can be combined in order to approximate the two-loop box amplitudes
with good accuracy over the whole phase space. In the second part of the
paper, then, we show that this can be done also for $gg\to ZZ$, and we discuss
the combination of our new results with those of ref. [42].

The paper is organized as follows. In section 2 we set the notation and we
comment on the application of the $p_{\scriptscriptstyle T}$ expansion to the
$gg\to ZZ$ case. The next section is devoted to a comparison between the known
LO amplitude and our approximation. In section 4 we then consider the
application of the $p_{\scriptscriptstyle T}$ expansion at NLO and and discuss
how our $p_{\scriptscriptstyle T}$-expanded results can be merged with those
derived via the HE expansion. Section 5 contains our NLO results, and we
present our conclusions in section 6. The paper is complemented by two
appendices. In appendix A we report the explicit expressions for the
projectors we employ in the calculation. We present also the relation between
our form factors and the ones used in ref. [42]. In appendix B we report the
exact results for the NLO triangle and the reducible double-triangle
contributions.

##  2 Definitions and the $p_{\scriptscriptstyle T}$-expansion method

###  2.1 Definitions

We consider the process $g_{a}^{\mu}(p_{1})g_{b}^{\nu}(p_{2})\rightarrow
Z^{\rho}(p_{3})Z^{\sigma}(p_{4})$. The amplitude can be defined as

$\mathcal{A}=\sqrt{2}m_{\scriptscriptstyle
Z}^{2}G_{F}~{}\frac{\alpha_{s}(\mu_{R})}{\pi}\delta_{ab}~{}\epsilon^{a}_{\mu}(p_{1})\epsilon^{b}_{\nu}(p_{2})\epsilon_{\rho}^{*}(p_{3})\epsilon_{\sigma}^{*}(p_{4})~{}\hat{\mathcal{A}}^{\mu\nu\rho\sigma}(p_{1},p_{2},p_{3}),$
(1)

where $G_{F}$ is the Fermi constant, $\alpha_{s}(\mu_{R})$ is the strong
coupling constant evaluated at a renormalisation scale $\mu_{R}$ and the
polarization vectors of the gluons and the $Z$ bosons are
$\epsilon^{a}_{\mu}(p_{1}),\epsilon^{b}_{\nu}(p_{2})$ and
$\epsilon_{\rho}(p_{3}),\epsilon_{\sigma}(p_{4})$, respectively. The Lorentz
structure of the amplitude is encoded in the tensor
$\hat{\mathcal{A}}^{\mu\nu\rho\sigma}(p_{1},p_{2},p_{3})$, whose most general
decomposition consists of 138 Lorentz structures. However, by imposing the
transversality of the external polarization vectors w.r.t. the relative four-
momentum

$\epsilon(p_{i})\cdot p_{i}=0\qquad i=1,\dots,4,$ (2)

and by fixing the gauge of the external gluons with

$\epsilon(p_{1})\cdot p_{2}=0\qquad\epsilon(p_{2})\cdot p_{1}=0,$ (3)

$\hat{\mathcal{A}}^{\mu\nu\rho\sigma}$ can be written as a linear combination
of 20 Lorentz structures [36, 50, 42]

$\hat{\mathcal{A}}^{\mu\nu\rho\sigma}(p_{1},p_{2},p_{3})=\sum_{i=1}^{20}S_{i}^{\mu\nu\rho\sigma}f_{i}(\hat{s},\hat{t},\hat{u},m_{t},m_{\scriptscriptstyle
Z}),$ (4)

where the scalar form factors $f_{i}$ depend, besides $m_{t}$ and
$m_{\scriptscriptstyle Z}$, on the partonic Mandelstam variables. Assuming all
momenta to be incoming, the latter are defined as

$\hat{s}=(p_{1}+p_{2})^{2},\quad\hat{t}=(p_{1}+p_{3})^{2},\quad\hat{u}=(p_{2}+p_{3})^{2}$
(5)

and the relation $\hat{s}+\hat{t}+\hat{u}=2\,m_{\scriptscriptstyle Z}^{2}$ is
satisfied. We checked that the Lorentz structures that give a nonzero
contribution to the amplitude are

$\displaystyle S_{1}^{\mu\nu\rho\sigma}$
$\displaystyle=g^{\mu\nu}g^{\rho\sigma}$ $\displaystyle
S_{2}^{\mu\nu\rho\sigma}$ $\displaystyle=g^{\mu\rho}g^{\nu\sigma}$
$\displaystyle S_{3}^{\mu\nu\rho\sigma}$
$\displaystyle=g^{\mu\sigma}g^{\nu\rho}$ $\displaystyle
S_{4}^{\mu\nu\rho\sigma}$ $\displaystyle=p_{1}^{\rho}p_{3}^{\nu}g^{\mu\sigma}$
$\displaystyle S_{5}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{2}^{\rho}p_{3}^{\nu}g^{\mu\sigma}$ $\displaystyle
S_{6}^{\mu\nu\rho\sigma}$ $\displaystyle=p_{1}^{\rho}p_{3}^{\mu}g^{\nu\sigma}$
$\displaystyle S_{7}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{2}^{\rho}p_{3}^{\mu}g^{\nu\sigma}$ $\displaystyle
S_{8}^{\mu\nu\rho\sigma}$ $\displaystyle=p_{3}^{\mu}p_{3}^{\nu}g^{\rho\sigma}$
$\displaystyle S_{9}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\rho}p_{1}^{\sigma}g^{\mu\nu}$ $\displaystyle
S_{10}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\rho}p_{2}^{\sigma}g^{\mu\nu}$ $\displaystyle
S_{11}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\sigma}p_{2}^{\rho}g^{\mu\nu}$ $\displaystyle
S_{12}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{2}^{\rho}p_{2}^{\sigma}g^{\mu\nu}$ (6) $\displaystyle
S_{13}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\sigma}p_{3}^{\nu}g^{\mu\rho}$ $\displaystyle
S_{14}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{2}^{\sigma}p_{3}^{\nu}g^{\mu\rho}$ $\displaystyle
S_{15}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\sigma}p_{3}^{\mu}g^{\nu\rho}$ $\displaystyle
S_{16}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{2}^{\sigma}p_{3}^{\mu}g^{\nu\rho}$ $\displaystyle
S_{17}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\rho}p_{1}^{\sigma}p_{3}^{\mu}p_{3}^{\nu}$
$\displaystyle S_{18}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\rho}p_{2}^{\sigma}p_{3}^{\mu}p_{3}^{\nu}$
$\displaystyle S_{19}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{1}^{\sigma}p_{2}^{\rho}p_{3}^{\mu}p_{3}^{\nu}$
$\displaystyle S_{20}^{\mu\nu\rho\sigma}$
$\displaystyle=p_{2}^{\rho}p_{2}^{\sigma}p_{3}^{\mu}p_{3}^{\nu}$

where we follow the numbering of ref. [42].

In order to simplify the evaluation of the cross section, in our work we
express the amplitude in terms of a set of orthonormal projectors

$\hat{\mathcal{A}}^{\mu\nu\rho\sigma}(p_{1},p_{2},p_{3})=\sum_{i=1}^{20}\mathcal{P}_{i}^{\mu\nu\rho\sigma}\mathcal{A}_{i}(\hat{s},\hat{t},\hat{u},m_{t},m_{\scriptscriptstyle
Z}),$ (7)

where the tensors $\mathcal{P}_{i}^{\mu\nu\rho\sigma}$ are constructed as
linear combinations of the Lorentz structures defined in eqs. (6), using a
Gram-Schmidt procedure for the orthogonalization. In appendix A we give the
explicit expressions for the projectors. Here we point out that, to
efficiently perform the $p_{\scriptscriptstyle T}$ expansion [47], we choose
the projectors to be either symmetric or anti-symmetric under the interchange
$\\{\mu\leftrightarrow\nu,p_{1}\leftrightarrow p_{2}\\}$. This choice also
allows to reduce the number of relevant form factors222In fact, we observed
that enforcing the additional Bose symmetry
$\\{\rho\leftrightarrow\sigma,p_{3}\leftrightarrow p_{4}\\}$ further reduces
the relevant form factors to 12. While this may be used for improving the
practical implementation of our results, in this paper we use the 16 form
factors. from 20 to 16. We present our results in terms of the
$\mathcal{A}_{i}$ form factors of eq. (7), while in appendix A we include the
relations to obtain the latter as a combination of the $f_{i}$ in eqs. (4).

(a)

(b)

Figure 1: Representative Feynman diagrams contributing to the $gg\rightarrow
ZZ$ amplitude at LO. Only the contribution from top-quark loops is shown.

We consider a perturbative expansion of the form factors in the strong
coupling

$\mathcal{A}_{i}=\mathcal{A}_{i}^{(0)}+\frac{\alpha_{s}}{\pi}\mathcal{A}_{i}^{(1)}+\mathcal{O}(\alpha_{s}^{2})$
(8)

where one- and two-loop diagrams contribute respectively to the LO
$(\mathcal{A}_{i}^{(0)})$ and NLO $(\mathcal{A}_{i}^{(1)})$. According to the
topology of the relevant Feynman diagrams (see fig. 1), we identify a triangle
and a box contribution to the LO form factors

$\mathcal{A}_{i}^{(0)}=\mathcal{A}_{i}^{(0,\triangle)}+\mathcal{A}_{i}^{(0,\square)}.$
(9)

The above classification is modified at NLO, where the two-loop triangle and
box topologies (see fig. 2) are supplemented with one-particle-reducible
double-triangle diagrams as in fig. 2(c). Therefore the NLO form factors are
defined as

$\mathcal{A}_{i}^{(1)}=\mathcal{A}_{i}^{(1,\triangle)}+\mathcal{A}_{i}^{(1,\square)}+\mathcal{A}_{i}^{(1,\bowtie)}.$
(10)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 2: Representative Feynman diagrams contributing to the $gg\rightarrow
ZZ$ amplitude at NLO. Loops of bottom quarks are included only in the double-
triangle diagrams (c).

The main result of this paper is the evaluation of the
$\mathcal{A}_{i}^{(1,\square)}$ using the $p_{\scriptscriptstyle T}$
expansion, which is described below. For completeness, we provide also the
results for $\mathcal{A}_{i}^{(1,\triangle)}$ and
$\mathcal{A}_{i}^{(1,\bowtie)}$. We notice that only for
$\mathcal{A}_{i}^{(1,\bowtie)}$ we include the contributions from loops of
bottom as well as top quarks, since the whole fermion generation must be
considered to remove the axial anomaly.

With the previous definitions, the partonic cross section at LO is expressed
as

$\hat{\sigma}^{(0)}(\hat{s})=\frac{m_{\scriptscriptstyle
Z}^{4}G_{F}^{2}}{512\pi\hat{s}^{2}}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{2}\int_{\hat{t}^{-}}^{\hat{t}^{+}}d\hat{t}\sum_{i}\left|\mathcal{A}_{i}^{(0)}\right|^{2}$
(11)

with $\hat{t}^{\pm}=1/2\left[-\hat{s}+2m_{\scriptscriptstyle
Z}^{2}\pm\sqrt{\hat{s}^{2}-\hat{s}~{}4m_{\scriptscriptstyle Z}^{2}}\right]$.

###  2.2  The $p_{\scriptscriptstyle T}$-expansion method

We now briefly recall the main points of the $p_{\scriptscriptstyle T}$
expansion. The details of the method are presented in refs. [46, 47]. The
Feynman amplitude is expanded in the transverse momentum at the integrand
level via the introduction of the vector $r^{\mu}=p_{1}^{\mu}+p_{3}^{\mu}$,
which satisfies the relations

$r^{2}=\hat{t}\qquad r\cdot p_{1}=t^{\prime}\qquad r\cdot p_{2}=-t^{\prime}$
(12)

and can be written as

$r^{\mu}=\frac{t^{\prime}}{s^{\prime}}(p_{2}-p_{1})^{\mu}+r_{\perp}^{\mu},$
(13)

where the space-like vector $r_{\perp}^{\mu}$ is such that

$r_{\perp}\cdot p_{1}=r_{\perp}\cdot p_{2}=0,\qquad
r_{\perp}^{2}=-p_{\scriptscriptstyle T}^{2},$ (14)

with

$p_{\scriptscriptstyle T}^{2}=\frac{\hat{t}\hat{u}-m_{\scriptscriptstyle
Z}^{4}}{\hat{s}}.$ (15)

In eqs. (12,13) we introduce the primed Mandelstam variables

$s^{\prime}=p_{1}\cdot p_{2}=\frac{\hat{s}}{2},\quad t^{\prime}=p_{1}\cdot
p_{3}=\frac{\hat{t}-m_{\scriptscriptstyle Z}^{2}}{2},\quad
u^{\prime}=p_{2}\cdot p_{3}=\frac{\hat{u}-m_{\scriptscriptstyle Z}^{2}}{2},$
(16)

such that

$s^{\prime}+t^{\prime}+u^{\prime}=0,\qquad p_{\scriptscriptstyle
T}^{2}=\frac{2t^{\prime}u^{\prime}}{s^{\prime}}-m_{\scriptscriptstyle Z}^{2}.$
(17)

The above relations allow to rewrite $t^{\prime}$ as

$t^{\prime}=-\frac{s^{\prime}}{2}\left\\{1\pm\sqrt{1-2\frac{p_{\scriptscriptstyle
T}^{2}+m_{\scriptscriptstyle Z}^{2}}{s^{\prime}}}\right\\}.$ (18)

We are interested in the expansion of the amplitude in the forward limit,
corresponding to $t^{\prime}\sim 0$, i.e the minus-sign case in eq. (18). As
discussed in refs. [46, 47], the forward expansion is sufficient to obtain the
correct result for the cross section, if the (anti-)symmetry of the form
factors with respect to the exchange $t^{\prime}\leftrightarrow u^{\prime}$ is
ensured. The region of validity of the expansion is given by the condition
$p_{\scriptscriptstyle T}^{2}/(4\,m_{t}^{2})\lesssim 1$.

The $p_{\scriptscriptstyle T}$ expansion of the relevant Feynman diagrams
returns a result in terms of the ratios of small vs large quantities, $x/y$,
where $x\in\\{p_{\scriptscriptstyle T}^{2},m_{\scriptscriptstyle Z}^{2}\\}$
and $y\in\\{\hat{s},m_{t}^{2}\\}$. We notice that, in the
$p_{\scriptscriptstyle T}$ expansion, the transverse momentum and
$m_{\scriptscriptstyle Z}$ are treated on the same footing, i.e. a term
$\mathcal{O}(p_{\scriptscriptstyle T}^{2}/m_{\scriptscriptstyle Z}^{2})$ is
assumed to be $\mathcal{O}(1)$, not to be $\mathcal{O}(p_{\scriptscriptstyle
T}^{2})$. Therefore the expansion can be also considered in terms of the
quantities $m_{\scriptscriptstyle Z}^{2}$ and of the combination
$p_{\scriptscriptstyle T}^{2}+m_{\scriptscriptstyle Z}^{2}$ that enters in
$t^{\prime}$ (see eq. (18)).

After the scalar form factors $\mathcal{A}_{i}$ are expanded in the small
parameters above, and after the scalar integrals are decomposed along a basis
of master integrals (MI) using Integration-by-Parts (IBP) identities, the form
factors can be written as the following series

$\mathcal{A}_{i}=\mathcal{N}(p_{\scriptscriptstyle
T}^{2},m_{\scriptscriptstyle
Z}^{2})\sum_{N=0}^{\infty}\sum_{i+j=N}c_{ij}(p_{\scriptscriptstyle
T}^{2})^{i}(m_{\scriptscriptstyle Z}^{2})^{j},$ (19)

where the $c_{ij}$ coefficients are linear combinations of the MIs resulting
from the IBP reduction, which in turn depend only on $\hat{s}$ and
$m_{t}^{2}$, while $\mathcal{N}(p_{\scriptscriptstyle
T}^{2},m_{\scriptscriptstyle Z}^{2})$ is an overall normalization factor which
may depend on $p_{\scriptscriptstyle T}^{2}$ and $m_{\scriptscriptstyle
Z}^{2}$.

##  3 Validation of the $p_{\scriptscriptstyle T}$ expansion at LO

(a)

(b)

Figure 3: Absolute value of the form factor $\mathcal{A}_{9}^{(0,\square)}$
for moderate (a) and high (b) partonic centre-of-mass energies as a function
of the transverse momentum. The exact result and the results obtained at
various orders in the $p_{\scriptscriptstyle T}$ expansion are shown.

In this section we assess the level of accuracy of the $p_{\scriptscriptstyle
T}$ expansion in reproducing the known exact result at LO. We used FeynArts
[51] to generate the relevant amplitudes, which we expanded in the limit of
small $p_{\scriptscriptstyle T}$ using private code. The latter relies on
several functions implemented in FeynCalc [52, 53]. The expanded amplitude has
been decomposed along a basis of MIs with LiteRed [54, 55]. The MIs can be
expressed in terms of the $B_{0}$ and $C_{0}$ Passarino-Veltman functions [56]
of argument $m_{t}^{2}$ and $\hat{s}$ or $-\hat{s}$. We also used FeynCalc to
recompute the exact LO result for the amplitude and found it in
agreement333The comparison has been done for each form factor, using the
conversion formulas listed in eqs. (A.21)-(A.38) in appendix A. with the one
provided in ref. [42]. Since at NLO only the box contributions will be
computed in an approximate way via the $p_{\scriptscriptstyle T}$ expansion,
in the following we focus on the results of the box form factors,
$\mathcal{A}_{i}^{(0,\square)}$. In particular, we consider the absolute value
of the form factor $\mathcal{A}_{9}^{(0,\square)}$ which provides the largest
contribution to the cross section. We notice that the
$\mathcal{P}_{9}^{\mu\nu\rho\sigma}$ projector, as well as others, exhibits a
$\mathcal{O}(p_{\scriptscriptstyle T}^{-2})$ term (see eq. (A.9)) that can
question the fact that the corresponding form factor,
$\mathcal{A}_{9}^{(0,\square)}$, must be regular in the limit
$p_{\scriptscriptstyle T}\to 0$. However, the
$\mathcal{O}(p_{\scriptscriptstyle T}^{-2})$ terms that appear in the
projectors are always multiplied by combination of the
$S_{i}^{\mu\nu\rho\sigma}\>(i=17,\,18,\,19,\,20)$ Lorentz structures that,
once the gauge condition eq. (3) is enforced, exhibit an
$r_{\perp}^{\mu}\,r_{\perp}^{\nu}$ dependence. This dependence ensures that
the form factors are regular in the limit $p_{\scriptscriptstyle T}\to 0$.

In fig. 3 we plot the exact result for the $\mathcal{A}_{9}^{(0,\square)}$
form factor, compared to various $p_{\scriptscriptstyle T}$-expanded results
for two values of the partonic centre-of-mass energy. The expanded-over-exact
ratio in the bottom parts of the figure shows that the convergence of the
expansion is rather fast and that already the
$\mathcal{O}(p_{\scriptscriptstyle T}^{2})$ result can reproduce the exact one
with an accuracy below the percent level. The spikes that are visible in the
expanded-over-exact ratio are due to the fact that the form factor is zero
near $p_{\scriptscriptstyle T}\sim 100$ GeV, and the ratio is not numerically
stable. From a comparison between figs. 3(a) and 3(b) we observe that,
independently on the value of $\sqrt{\hat{s}}$ the expansion is convergent
only for values $p_{\scriptscriptstyle T}\lesssim 300$ GeV. This is consistent
with the limit $p_{\scriptscriptstyle T}^{2}\lesssim 4\,m_{t}^{2}$ that is
assumed in the expansion.

$M_{\scriptscriptstyle ZZ}(\text{GeV})$ | $\hat{\sigma}^{(0)}(\text{fb})$ |  |  |
---|---|---|---|---
| $\mathcal{O}(p_{\scriptscriptstyle T}^{0})$ | $\mathcal{O}(p_{\scriptscriptstyle T}^{2})$ | $\mathcal{O}(p_{\scriptscriptstyle T}^{4})$ | Exact
250 | 0.12486 | 0.12418 | 0.12412 | 0.12412
345 | 0.26619 | 0.26460 | 0.26448 | 0.26447
475 | 0.44887 | 0.44387 | 0.44368 | 0.44367
659 | 0.43135 | 0.38051 | 0.38644 | 0.38509
808 | 0.44108 | 0.30622 | 0.34712 | 0.32024
1255 | 0.70815 | 1.0303 | 0.66139 | 0.1883
Table 1: The partonic LO cross section as a function of $M_{\scriptscriptstyle
ZZ}$. The exact result and the ones obtained at various orders in the
$p_{\scriptscriptstyle T}$ expansion are shown.

In table 1 we present the values of the LO partonic cross section for $gg\to
ZZ$ as a function of the invariant mass of the two $Z$ bosons system. The
exact result is compared to various orders in the $p_{\scriptscriptstyle T}$
expansion. The triangle contribution (see fig. 1(a)) is evaluated always
exactly. The value of the strong coupling is fixed at $\alpha_{s}=0.118$,
while we use as other parameters $G_{F}=1.1663786\cdot
10^{-5}~{}\text{GeV}^{-2},m_{H}=125.1~{}\text{GeV},m_{Z}=91.1876~{}\text{GeV},m_{t}=173.21~{}\text{GeV}$.
The table shows that for $M_{\scriptscriptstyle ZZ}\lesssim 500$ GeV already
the $\mathcal{O}(p_{\scriptscriptstyle T}^{2})$ term is in agreement with the
exact result at the level of less than 1 per mille. Increasing
$M_{\scriptscriptstyle ZZ}$ up to 650 GeV the agreement between the exact and
the $p_{\scriptscriptstyle T}$-expanded result starts to deteriorate reaching
a $\sim 1$ per cent difference. However, the inclusion of the
$\mathcal{O}(p_{\scriptscriptstyle T}^{4})$ term brings back the difference to
a few per mille level. For higher values of $M_{\scriptscriptstyle ZZ}$ the
most important contributions to the cross section come from regions in the
phase space where the $p_{\scriptscriptstyle T}$ expansion is not valid and
the agreement between the exact and the $p_{\scriptscriptstyle T}$-expanded
results is poor.

##  4 Virtual corrections at NLO

Having showed that the $p_{\scriptscriptstyle T}$ expansion can provide
accurate results for the LO contribution, we move now to discuss its
application to the NLO QCD corrections, see fig. 2 for representative Feynman
diagrams. In fact, we are going to use the $p_{\scriptscriptstyle T}$
expansion only for the evaluation of the two-loop box diagrams. For the two-
loop one-particle-irreducible triangles we have adapted the full analytic
results of ref. [34], while the reducible double triangles have been computed
exactly using FeynCalc and checked with the results of ref. [40]. With
reference to eq. (10), we list the results for the form factors
$\mathcal{A}_{i}^{(1,\triangle)}$ and $\mathcal{A}_{i}^{(1,\bowtie)}$ in
appendix B.

###  4.1 Box diagrams

The evaluation of the two-loop box diagrams was carried out as follows. We
generated the diagrams with FeynArts and performed all the manipulations of
the amplitude within Mathematica and FeynCalc. The amplitude was contracted
with the 18 projectors in appendix A and each form factor expanded up to
$\mathcal{O}(p_{\scriptscriptstyle T}^{4})$. This led to the identification of
a set of 9 families of scalar integrals. These families were analyzed using
LiteRed and reduced to a basis of 52 known MIs [34, 57, 33, 58, 59, 60], which
we found to be the same as those encoutered in $HH$ and $ZH$ production. Among
the 52 MI, fifty can be expressed in terms of (generalised) harmonic
polylogarithms while two are elliptic integrals. The evaluation of the
(generalised) harmonic polylogarithms was done using the code handyG [61],
while the elliptic integrals were evaluated using the routines of ref. [62].
The top quark mass was renormalised in the on-shell scheme and the infra-red
(IR) poles were subtracted as in refs. [63, 64].

As a first check of our computation, we compared our $p_{\scriptscriptstyle
T}$-expanded result with the LME one presented in ref. [42]. It should be
noticed that $p_{\scriptscriptstyle T}$-expanded result actually “contains”
the LME one. The LME differs from the expansion in $p_{\scriptscriptstyle T}$
by the fact that $\hat{s}$ is treated as a small parameter with respect to
$m_{t}^{2}$ , and not on the same footing as in the latter case. This implies
that when the $p_{\scriptscriptstyle T}$-expanded result is further expanded
in terms of the $\hat{s}/m_{t}^{2}$ ratio the LME result has to be recovered.
This way, we were able to reproduce, at the analytic level, the LME result of
ref. [42] up to the order of our calculation.

$\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}$ | Numerical | $\mathcal{O}(p_{\scriptscriptstyle T}^{4})$
---|---|---
$++++$ | 3.15549(8) + i 0.47235(8) | 3.16038 + i 0.46980
$+++-$ | 0.15950(7) + i 0.14052(8) | 0.16305 + i 0.13885
$+-+-$ | -0.38609(7) + i 0.10539(7) | -0.38617 + i 0.11085
$-+++$ | -0.46990(8) + i 0.40207(8) | -0.46956 + i 0.40506
$+++0$ | 1.1248(2) - i 0.0805(2) | 1.1256 - i 0.0811
$+-+0$ | -1.4803(2) + i 0.4940(2) | -1.4799 + i 0.4977
$++00$ | 17.2585(6) + i 29.5669(6) | 17.2602 + i 29.5643
$+-00$ | 10.2869(5) - i 1.0571(6) | 10.2840 - i 1.0640
Table 2: Comparison of the NLO helicity amplitudes between the exact numerical
values and the values given by the $p_{\scriptscriptstyle T}$ expansion at
$\mathcal{O}(p_{\scriptscriptstyle T}^{4})$ in the kinematic point provided in
ref. [43].

As a second check of our computation we compare, in table 2, the contribution
of the two-loop box diagrams at $\mathcal{O}(p_{\scriptscriptstyle T}^{4})$
with the results of ref. [43], which rely on the numerical evaluation of the
scalar integrals. The comparison is done at the level of the helicity
amplitudes444The helicities of the polarization vectors are defined as:
$\epsilon_{\mu}^{\lambda_{1}}(p_{1})\,\epsilon_{\nu}^{\lambda_{2}}(p_{2})\,\epsilon^{*\,\lambda_{3}}_{\rho}(p_{3})\,\epsilon^{*\,\lambda_{4}}_{\sigma}(p_{4})$.
for the fixed phase-space point presented in ref. [43], defined as
$\hat{s}=142/17~{}m_{t}^{2}$ and
$\hat{t}=-125/22~{}m_{t}^{2}\sim-5.7~{}m_{t}^{2}$. We remark that the IR
subtraction scheme employed in ref. [43] ($q_{T}$ subtraction) differs from
ours, so that to reproduce the form factors of ref. [43],
$\mathcal{A}_{i}^{(1)q_{T}}$, a shift in our form factors has to be introduced
as follows:

$\mathcal{A}_{i}^{(1)q_{T}}=\mathcal{A}_{i}^{(1)}+\frac{C_{A}}{4}\pi^{2}~{}.$
(20)

We notice that the point chosen in ref. [43] lies somewhat beyond the region
of validity of the $p_{\scriptscriptstyle T}$ expansion. Yet, table 2 shows
that the relative difference between the $\mathcal{O}(p_{\scriptscriptstyle
T}^{4})$ values and the numerical ones is, in general with the exception of
only one helicity amplitude, at the per mille level and often better, as in
the case of the amplitudes involving two longitudinally-polarized $Z$ bosons.
The latter amplitudes are the dominant ones, as consequence of the Goldstone
Boson Equivalence Theorem.

###  4.2 Merging the $p_{\scriptscriptstyle T}$ and High-Energy expansions

The results we have presented so far allow to efficiently approximate the two-
loop box integrals with full top-mass dependence over a specific region of the
phase space, which spans the $ZZ$ production threshold up to moderate partonic
centre-of-mass energies (for the LHC). In the high-energy region the
$p_{\scriptscriptstyle T}$ expansion is still accurate for any values of
$\hat{s}$ and $|\hat{t}|\lesssim\,4m_{t}^{2}$ or
$|\hat{u}|\lesssim\,4m_{t}^{2}$ but it cannot cover the complementary region
$|\hat{t}|,|\hat{u}|\gtrsim 4m_{t}^{2}$ that becomes allowed by the
kinematics.

In ref. [48] we showed that the $p_{\scriptscriptstyle T}$ expansion and the
HE expansion are accurate in complementary phase-space regions. The HE
expansion is an expansion in terms of the ratios $x/y$ where
$x\in\\{m_{t}^{2},m_{\scriptscriptstyle Z}^{2}\\}$ and
$y\in\\{\hat{s},\hat{t},\hat{u}\\}$ and it is valid in phase-space regions
where both $|\hat{t}|/(4\,m_{t}^{2})\gtrsim 1$ and
$|\hat{u}|/(4\,m_{t}^{2})\gtrsim 1$. We also showed that if each expansion is
extended beyond its border of validity via the use of Padé approximants the
two results can be merged into a single prediction that is accurate over the
complete phase space.

Results for the evaluation of the two-loop box contribution in the $gg\to ZZ$
process via the HE expansion were presented in ref. [42]. Applying the method
of ref. [48] we merge our results for the box contribution with those obtained
in ref. [42]. Below we recall the main steps of the merging method:

  1. 1.

We use the formulas in eqs. (A.21)-(A.38) to adapt the form factors of ref.
[42] to our decomposition of the amplitude.

  2. 2.

We rewrite the small quantities on which the $p_{\scriptscriptstyle T}$ and
the HE expansions are based in terms of a scaling parameter $x$.

  3. 3.

We construct Padé approximants for each form factor, defined as

$[m/n](x)=\frac{p_{0}+p_{1}x+\dots+p_{m}x^{m}}{1+q_{1}x+\dots q_{n}x^{n}},$
(21)

where $p_{i},q_{i}$ are expressed as linear combinations of the coefficients
of each expansion. We refer to these approximants as the
$p_{\scriptscriptstyle T}$-Padé and HE-Padé.

  4. 4.

In evaluating the amplitude, we use the HE-Padé approximant for phase-space
points such that $|\hat{t}|>4m_{t}^{2}$ and $|\hat{u}|>4m_{t}^{2}$ whereas for
all the remaining phase-space points we use the $p_{\scriptscriptstyle
T}$-Padé approximant.

With our three terms in the $p_{\scriptscriptstyle T}$ expansion we
constructed a [1/1] $p_{\scriptscriptstyle T}$-Padé. The several terms in the
HE expansion presented in ref. [42] allow to construct different [x/y] HE-
Padé. We construct both the [5/5] and [6/6] HE-Padé and compare the results
using the different orders. As a result we found that both the [5/5] and [6/6]
HE-Padé give very similar results.

##  5 Results

(a)

(b)

Figure 4: (a) Absolute value of a form factor for a fixed $\sqrt{\hat{s}}$
using the $p_{\scriptscriptstyle T}$ and the HE expansion and the best Padé
approximants. (b) Partonic cross section at LO using the exact result, the
$p_{\scriptscriptstyle T}$ expansion and the prediction using the merged
$p_{\scriptscriptstyle T}$ and HE Padé approximants.

Before applying the merging approach for the form factors at NLO, we verified
its reliability testing it against the exact LO contribution. In fig. 4(a) the
exact result (solid black line) for the form factor,
$\mathcal{A}_{9}^{(0,\square)}$ as a function of $-\hat{t}/(4\,m_{t}^{2})$ is
compared to different approximate evaluations. In the figure the fixed order
results $\mathcal{O}(p_{\scriptscriptstyle T}^{4})$ in the
$p_{\scriptscriptstyle T}$ expansion (solid blue) and
$\mathcal{O}(m_{t}^{32})$ in the HE expansion (solid orange) are shown against
the Padé-improved versions, [1/1] $p_{\scriptscriptstyle T}$-Padé (dashed
light blue) and $[6/6]$ HE-Padé (dashed yellow). The bottom panel shows the
ratio of the different approximations over the exact result. It is evident
that the region $|\hat{t}|\sim 4\,m_{t}^{2}$, where both the
$p_{\scriptscriptstyle T}$ and HE expansions (solid blue and orange) begin to
diverge, is well covered by the respective Padé approximants, so that one has
a way to accurately reproduce the exact result for any value of $|\hat{t}|$.
In fig. 4(b) we compare different predictions at the level of the LO partonic
cross section. One can see that the accuracy of the merged result (dashed
green) is at the level of per mille or below, an order of magnitude better
than what observed in ref. [48]. We notice that, this remarkable degree of
accuracy is not achieved for each individual form factor. However, the impact
of larger deviations is suppressed with the integration over $\hat{t}$.

Figure 5: NLO form factor for $\sqrt{\hat{s}}=1200\text{ GeV}$ in various
approximations: $p_{\scriptscriptstyle T}$ expansion (solid blue), [1/1] Padé
approximant based on the $p_{\scriptscriptstyle T}$ expansion (dashed, light
blue), high-energy expansion (solid, pink) and [6/6] HE Padé approximant
(dashed, rosa). The lower panel shows the relative difference $\Delta$ (see
text).

(a)

(b)

Figure 6: Relative difference between several phase-space points of ref. [43]
and different approximations for the helicity amplitude ++00. Points in the
shaded region are outside the formal limit of validity of the
$p_{\scriptscriptstyle T}$ expansion. (a) [1/1] $p_{T}$-Padé only; (b) merging
of the $p_{T}$-Padé and the HE-Padé.

We now present our results at NLO. First, we checked their the validity in two
ways:

  * •

We compared the $p_{\scriptscriptstyle T}$-Padé and HE-Padé in the vicinity of
$|\hat{t}|\sim 4m_{t}^{2}$, for several values of $\hat{s}$ in the range [800,
2000] GeV, ensuring that the two analytic approximations are indeed
complementary and that the matching point is chosen well. In fig. 5 the NLO
form factor $\mathcal{A}_{9}^{(1,\square)}$ is shown as a function of
$-\hat{t}/(4\,m_{t}^{2})$ for both fixed-order $p_{\scriptscriptstyle T}$ and
HE expansions and Padé-improved results. The lower panel of the figure shows
the relative difference, $\Delta$, between the $p_{\scriptscriptstyle T}$-Padé
and the HE-Padé, defined as

$\Delta=2\left|\frac{[1/1]_{p_{\scriptscriptstyle
T}}-[6/6]_{\text{HE}}}{[1/1]_{p_{\scriptscriptstyle
T}}+[6/6]_{\text{HE}}}\right|,$ (22)

where one can see that at $|\hat{t}|\sim 4m_{t}^{2}$ the two expansion are in
good agreement (on the level of 0.5 per mille).

  * •

We compared our merged results with the helicity amplitudes of ref. [43],
obtained from a numerical evaluation of the scalar integrals. The results for
several phase-space points were provided privately by the authors of ref.
[43]. In fig. 6 we plot the results for the amplitude $++00$ using different
colors depending on the level of agreement $\Delta$ between the numerical and
merged prediction. The shaded region in the plots corresponds to the high-
energy region, where $|\hat{t}|>4m_{t}^{2}$ and $|\hat{u}|>4m_{t}^{2}$. In
fig. 6 (a) we compare the numerical evaluation with our results using only the
[1/1] $p_{\scriptscriptstyle T}$-Padé. Whereas the agreement is generally
better than per mille in the region of validity of the $p_{\scriptscriptstyle
T}$ expansion, deviations larger than 5% are visible in the high-energy region
for the real part. When we consider the merged prediction, the agreement with
the numerical evaluation is improved, and the majority of the phase space
points have differences below the per mille level. It must be noted that, as
seen in table 2, this level of agreement is not shared among all the helicity
amplitudes. This is however a consequence of the relations used to convert our
results in terms of form factors into helicity amplitudes. Indeed, for some
phase-space points, we observed delicate cancellations in the combination of
the form factors. At the same time, for these particular phase-space points,
the corresponding contributions to the helicity amplitudes are tiny.
Therefore, we expect an overall accuarcy at the subpercent level for our
results.

In fig. 7 one can see our result for the finite part of the virtual
corrections defined as

$\Delta\sigma_{\text{virt}}=\int_{\hat{t}^{-}}^{\hat{t}^{+}}d\hat{t}\frac{1}{2}\frac{1}{16\pi\hat{s}^{2}}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)\mathcal{V}_{\text{fin}}(\hat{t})$
(23)

where the finite part of the virtual correction is

$\mathcal{V}_{\text{fin}}=\frac{G_{F}^{2}m_{\scriptscriptstyle
Z}^{4}}{16}\left(\frac{\alpha_{s}(\mu_{R})}{\pi}\right)^{2}\left\\{\sum_{i}\left|\mathcal{A}_{i}^{(0)}\right|^{2}\frac{C_{A}}{2}\left(\pi^{2}-\log^{2}\left(\frac{\mu_{R}^{2}}{\hat{s}}\right)\right)+2\sum_{i}\text{Re}\left[\mathcal{A}_{i}^{(0)}\left(\mathcal{A}_{i}^{(1)}\right)^{*}\right]\right\\}.$
(24)

We used $\mu_{R}=M_{\scriptscriptstyle ZZ}/2$ as renormalisation scale. The
pink line shows the merged result of a [1/1] $p_{\scriptscriptstyle
T}$-expanded Padé approximant and a [6/6] HE Padé approximant. For comparison
we show the large mass expanded results of ref. [42]. Within the range of
validity of the large mass expansion the results agree very well. We have
checked that using a [5/5] HE Padé approximant would change the result only
below the integration error.

Figure 7: The partonic virtual corrections as a function of
$M_{\scriptscriptstyle ZZ}$ merging the $p_{\scriptscriptstyle T}$ expansion
and the HE (pink solid line) compared to the virtual corrections obtained when
using the large mass expanded results for the NLO box and triangles (dashed
blue line).

##  6 Conclusion

The accurate description of the process $gg\to ZZ$ is needed in order to
improve the extraction of the total width of the Higgs boson, as well as a
precision test for the SM. In this paper we have presented the calculation of
the top-quark loops in the virtual QCD corrections at NLO. Our main result is
the computation of the box diagrams using the $p_{\scriptscriptstyle T}$
expansion, which gives accurate results in a phase-space region that so far
has not been covered by other analytic approximations.

As a by-product of our calculation, we have used the results obtained in ref.
[42] to complement our analytic calculation, showing that a combination of the
two approaches can provide an efficient approximation for the cross section,
at a level of accuracy that is more than adequate for phenomenological
applications. We emphasise that our analytic approach is very flexible and can
be included into Monte Carlo programmes, as was already done for the case of
Higgs pair production, $gg\to HH$ [67]. It allows to change easily the input
parameters, which is for instance necessary to evaluate the top-mass
renormalisation scheme uncertainty that is expected to be large in the top-
mediated contributions. Indeed, it is currently the largest uncertainty for
the $gg\to HH$ [65, 66] and $gg\to ZH$ [68] processes. In addition, the
flexibility of the approach allows for a straightforward application to
beyond-the-SM scenarios.

We recall that for a full description of $gg\to ZZ$ at NLO in QCD also the
contribution mediated by loops of light quarks needs to be considered [36,
37], as well as the one from real-emission diagrams [17]. We leave the
combination of these effects with our results to future work.

## Acknowledgments

We thank Luigi Bellafronte and Pier Paolo Giardino for participating in the
early stages of this project. We are greatful to the authors of ref. [43] for
providing us with the results for additional phase-space points and for useful
discussions. The work of RG is supported in part by the Italian MUR
Departments of Excellence grant 2023-2027 ”Quantum Frontiers” and the ICSC –
Centro Nazionale di Ricerca in High Performance Computing, Big Data and
Quantum Computing, funded by European Union – NextGenerationEU. The authors
acknowledge support from the COMETA COST Action CA22130. The diagrams in this
paper have been drawn using JaxoDraw [69].

##  A Projectors

We present the explicit expressions of the orthonormal projectors
$\mathcal{P}_{i}^{\mu\nu\rho\sigma}$ appearing in eq. (7). Following the
notation of section 2, the antisymmetric projectors under the exchange
$\\{\mu\leftrightarrow\nu$ , $p_{1}\leftrightarrow p_{2}\\}$ are

$\displaystyle\mathcal{P}_{1}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}-S_{20}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.1) $\displaystyle\mathcal{P}_{2}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}-\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.2) $\displaystyle\mathcal{P}_{3}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{S_{5}^{\mu\nu\rho\sigma}}{u^{\prime}}-\frac{S_{6}^{\mu\nu\rho\sigma}}{t^{\prime}}-\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}-S_{20}^{\mu\nu\rho\sigma}+\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}-\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.3) $\displaystyle\mathcal{P}_{4}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{1}{s^{\prime}}\left(S_{9}^{\mu\nu\rho\sigma}-S_{12}^{\mu\nu\rho\sigma}\right)+\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}-S_{20}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.4) $\displaystyle\mathcal{P}_{5}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{1}{s^{\prime}}\left(\frac{t^{\prime}}{u^{\prime}}S_{11}^{\mu\nu\rho\sigma}-\frac{u^{\prime}}{t^{\prime}}S_{10}^{\mu\nu\rho\sigma}\right)+\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}-\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.5) $\displaystyle\mathcal{P}_{6}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{S_{13}^{\mu\nu\rho\sigma}}{u^{\prime}}-\frac{S_{16}^{\mu\nu\rho\sigma}}{t^{\prime}}-\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}-S_{20}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}-\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.6) $\displaystyle\mathcal{P}_{7}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle
Z}^{2}}{\sqrt{2}p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{\left(m_{\scriptscriptstyle Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)}{2m_{\scriptscriptstyle
Z}^{2}}\left(\frac{S_{13}^{\mu\nu\rho\sigma}}{u^{\prime}}-\frac{S_{16}^{\mu\nu\rho\sigma}}{t^{\prime}}\right)+\frac{1}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}}\left(u^{\prime}S_{14}^{\mu\nu\rho\sigma}-t^{\prime}S_{15}^{\mu\nu\rho\sigma}\right)$
(A.7) $\displaystyle+\frac{1}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}-S_{20}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}-\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
$\displaystyle\mathcal{P}_{8}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle
Z}^{2}}{\sqrt{2}p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{1}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}}\left(u^{\prime}S_{4}^{\mu\nu\rho\sigma}-t^{\prime}S_{7}^{\mu\nu\rho\sigma}\right)+\frac{\left(m_{\scriptscriptstyle
Z}^{2}-p_{\scriptscriptstyle T}^{2}\right)}{2m_{\scriptscriptstyle
Z}^{2}}\left(\frac{S_{5}^{\mu\nu\rho\sigma}}{u^{\prime}}-\frac{S_{6}^{\mu\nu\rho\sigma}}{t^{\prime}}\right)$
(A.8) $\displaystyle+\frac{1}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}-S_{20}^{\mu\nu\rho\sigma}+\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}-\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}.$

The symmetric projectors are

$\displaystyle\mathcal{P}_{9}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}^{2}}{\sqrt{p_{\scriptscriptstyle
T}^{4}+m_{\scriptscriptstyle Z}^{4}}}\Biggl{[}\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.9) $\displaystyle\mathcal{P}_{10}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{\sqrt{p_{\scriptscriptstyle T}^{4}+m_{\scriptscriptstyle
Z}^{4}}}{2p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{\left(m_{\scriptscriptstyle Z}^{4}-p_{\scriptscriptstyle
T}^{4}\right)}{\left(m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4}\right)p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}\right)+\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.10) $\displaystyle\mathcal{P}_{11}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}^{2}}{\sqrt{p_{\scriptscriptstyle
T}^{4}+m_{\scriptscriptstyle
Z}^{4}}}\Biggl{[}\frac{1}{s^{\prime}}\left(S_{9}^{\mu\nu\rho\sigma}+S_{12}^{\mu\nu\rho\sigma}\right)+\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.11) $\displaystyle\mathcal{P}_{12}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{\sqrt{p_{\scriptscriptstyle T}^{4}+m_{\scriptscriptstyle
Z}^{4}}}{2p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{\left(m_{\scriptscriptstyle Z}^{4}-p_{\scriptscriptstyle
T}^{4}\right)}{\left(m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4}\right)s^{\prime}}\left(S_{9}^{\mu\nu\rho\sigma}+S_{12}^{\mu\nu\rho\sigma}\right)+\frac{1}{s^{\prime}}\left(\frac{u^{\prime}\
}{t^{\prime}}S_{10}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{11}^{\mu\nu\rho\sigma}\right)$
(A.12) $\displaystyle+\frac{\left(m_{\scriptscriptstyle
Z}^{4}-p_{\scriptscriptstyle T}^{4}\right)}{\left(m_{\scriptscriptstyle
Z}^{4}+p_{\scriptscriptstyle T}^{4}\right)p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}\right)+\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
$\displaystyle\mathcal{P}_{13}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle Z}}{\sqrt{2}p_{\scriptscriptstyle
T}}\Biggl{[}\frac{S_{13}^{\mu\nu\rho\sigma}}{u^{\prime}}+\frac{S_{16}^{\mu\nu\rho\sigma}}{t^{\prime}}-\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}\right)-\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
(A.13) $\displaystyle\mathcal{P}_{14}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{m_{\scriptscriptstyle
Z}^{2}}{\sqrt{2}p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{\left(m_{\scriptscriptstyle Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)}{2m_{\scriptscriptstyle
Z}^{2}}\left(\frac{S_{13}^{\mu\nu\rho\sigma}}{u^{\prime}}+\frac{S_{16}^{\mu\nu\rho\sigma}}{t^{\prime}}\right)+\frac{1}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}}\left(u^{\prime}S_{14}^{\mu\nu\rho\sigma}+t^{\prime}S_{15}^{\mu\nu\rho\sigma}\right)$
(A.14) $\displaystyle-\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}+\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
$\displaystyle\mathcal{P}_{15}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{\sqrt{2}m_{\scriptscriptstyle
Z}^{3}}{p_{\scriptscriptstyle T}(p_{\scriptscriptstyle
T}^{2}+m_{\scriptscriptstyle Z}^{2})}\Biggl{[}\frac{1}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}}\left(u^{\prime}S_{4}^{\mu\nu\rho\sigma}+t^{\prime}S_{7}^{\mu\nu\rho\sigma}\right)$
(A.15) $\displaystyle+\frac{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}{2m_{\scriptscriptstyle
Z}^{2}p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}+\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)\Biggr{]}$
$\displaystyle\mathcal{P}_{16}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}+m_{\scriptscriptstyle
Z}^{2}}{2\sqrt{2}p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{2m_{\scriptscriptstyle
Z}^{2}}{s^{\prime}p_{\scriptscriptstyle T}^{2}(p_{\scriptscriptstyle
T}^{2}+m_{\scriptscriptstyle
Z}^{2})}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}+\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\right)$
(A.16) $\displaystyle-\frac{2\left(p_{\scriptscriptstyle
T}^{2}-m_{\scriptscriptstyle Z}^{2}\right)}{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)^{2}s^{\prime}}\left(u^{\prime}S_{4}^{\mu\nu\rho\sigma}+t^{\prime}S_{7}^{\mu\nu\rho\sigma}\right)+\frac{S_{5}^{\mu\nu\rho\sigma}}{u^{\prime}}+\frac{S_{6}^{\mu\nu\rho\sigma}}{t^{\prime}}\Biggr{]}$
$\displaystyle\mathcal{P}_{17}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}+2m_{\scriptscriptstyle
Z}^{2}}{2m_{\scriptscriptstyle Z}^{2}}\Biggl{\\{}\frac{m_{\scriptscriptstyle
Z}^{2}}{2m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}}\left(S_{2}^{\mu\nu\rho\sigma}+S_{3}^{\mu\nu\rho\sigma}\right)-\frac{m_{\scriptscriptstyle
Z}^{2}}{(p_{\scriptscriptstyle T}^{2}+2m_{\scriptscriptstyle
Z}^{2})s^{\prime}p_{\scriptscriptstyle
T}^{2}}\left(u^{\prime}S_{4}^{\mu\nu\rho\sigma}+t^{\prime}S_{7}^{\mu\nu\rho\sigma}\right)$
(A.17) $\displaystyle-\frac{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)m_{\scriptscriptstyle
Z}^{2}}{p_{\scriptscriptstyle T}^{2}\left(2m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}\Biggl{[}\frac{S_{5}^{\mu\nu\rho\sigma}}{2u^{\prime}}+\frac{S_{6}^{\mu\nu\rho\sigma}}{2t^{\prime}}-\frac{S_{13}^{\mu\nu\rho\sigma}}{2u^{\prime}}-\frac{S_{16}^{\mu\nu\rho\sigma}}{2t^{\prime}}+\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\Bigl{(}S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}$
$\displaystyle+\frac{u^{\prime}}{t^{\prime}}S_{18}^{\mu\nu\rho\sigma}+\frac{t^{\prime}}{u^{\prime}}S_{19}^{\mu\nu\rho\sigma}\Bigr{)}\Biggr{]}+\frac{m_{\scriptscriptstyle
Z}^{2}}{p_{\scriptscriptstyle T}^{2}s^{\prime}\left(2m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}\left(u^{\prime}S_{14}^{\mu\nu\rho\sigma}+t^{\prime}S_{15}^{\mu\nu\rho\sigma}\right)\Biggr{\\}}$

$\displaystyle\mathcal{P}_{18}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{S_{8}^{\mu\nu\rho\sigma}}{p_{\scriptscriptstyle
T}^{2}}-\frac{\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}{2p_{\scriptscriptstyle
T}^{4}s^{\prime}}\left(S_{17}^{\mu\nu\rho\sigma}+S_{20}^{\mu\nu\rho\sigma}\right)+\frac{u^{\prime}\
\left(-m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4}-2p_{\scriptscriptstyle
T}^{2}(s^{\prime}-t^{\prime}+u^{\prime})\right)S_{18}^{\mu\nu\rho\sigma}}{2p_{\scriptscriptstyle
T}^{4}\ \left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}t^{\prime}}$ (A.18) $\displaystyle-\frac{t^{\prime}\
\left(m_{\scriptscriptstyle Z}^{4}-p_{\scriptscriptstyle
T}^{4}+2p_{\scriptscriptstyle
T}^{2}(s^{\prime}+t^{\prime}-u^{\prime})\right)S_{19}^{\mu\nu\rho\sigma}}{2p_{\scriptscriptstyle
T}^{4}\ \left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}\ u^{\prime}}.$

Finally, we include the last two projectors, which have null norm in $D=4$
dimensions and thus do not contribute to the amplitude.

$\displaystyle\mathcal{P}_{19}^{\mu\nu\rho\sigma}$
$\displaystyle=S_{1}^{\mu\nu\rho\sigma}-\frac{S_{2}^{\mu\nu\rho\sigma}}{2}-\frac{S_{3}^{\mu\nu\rho\sigma}}{2}+\frac{1}{2\
p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(u^{\prime}S_{4}^{\mu\nu\rho\sigma}+t^{\prime}S_{7}^{\mu\nu\rho\sigma}\right)+\frac{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}{4\ p_{\scriptscriptstyle
T}^{2}}\left(\frac{S_{5}^{\mu\nu\rho\sigma}}{u^{\prime}}+\frac{S_{6}^{\mu\nu\rho\sigma}}{t^{\prime}}\right)$
(A.19) $\displaystyle+\frac{S_{8}^{\mu\nu\rho\sigma}}{p_{\scriptscriptstyle
T}^{2}}-\frac{\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}{2p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(S_{9}^{\mu\nu\rho\sigma}+S_{12}^{\mu\nu\rho\sigma}\right)+\frac{u^{\prime}\
\left(-m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4}-2p_{\scriptscriptstyle
T}^{2}(s^{\prime}-t^{\prime}+u^{\prime})\right)S_{10}^{\mu\nu\rho\sigma}}{2p_{\scriptscriptstyle
T}^{2}\ \left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}t^{\prime}}$ $\displaystyle-\frac{t^{\prime}\
\left(m_{\scriptscriptstyle Z}^{4}-p_{\scriptscriptstyle
T}^{4}+2p_{\scriptscriptstyle
T}^{2}(s^{\prime}+t^{\prime}-u^{\prime})\right)S_{11}^{\mu\nu\rho\sigma}}{2p_{\scriptscriptstyle
T}^{2}\ \left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}\ u^{\prime}}-\frac{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}{4p_{\scriptscriptstyle T}^{2}\
}\left(\frac{S_{13}^{\mu\nu\rho\sigma}}{u^{\prime}}+\frac{S_{16}^{\mu\nu\rho\sigma}}{t^{\prime}}\right)$
$\displaystyle-\frac{1}{2p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(u^{\prime}S_{14}^{\mu\nu\rho\sigma}+t^{\prime}\
S_{15}^{\mu\nu\rho\sigma}\right)+\frac{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}-s^{\prime}+t^{\prime}-u^{\prime}\right)u^{\prime}S_{18}^{\mu\nu\rho\sigma}}{p_{\scriptscriptstyle
T}^{2}\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)\
s^{\prime}t^{\prime}}$ $\displaystyle+\frac{t^{\prime}\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}-s^{\prime}-t^{\prime}+u^{\prime}\right)\
S_{19}^{\mu\nu\rho\sigma}}{p_{\scriptscriptstyle
T}^{2}\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}\ u^{\prime}}$
$\displaystyle\mathcal{P}_{20}^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{1}{2m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}}\Biggl{[}S_{3}^{\mu\nu\rho\sigma}-S_{2}^{\mu\nu\rho\sigma}-\frac{1}{p_{\scriptscriptstyle
T}^{2}s^{\prime}}\left(u^{\prime}S_{4}^{\mu\nu\rho\sigma}-t^{\prime}S_{7}^{\mu\nu\rho\sigma}+u^{\prime}S_{14}^{\mu\nu\rho\sigma}-t^{\prime}S_{15}^{\mu\nu\rho\sigma}\right)$
(A.20) $\displaystyle-\frac{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}{2p_{\scriptscriptstyle
T}^{2}}\left(\frac{S_{5}^{\mu\nu\rho\sigma}}{u^{\prime}}-\frac{S_{6}^{\mu\nu\rho\sigma}}{t^{\prime}}+\frac{S_{13}^{\mu\nu\rho\sigma}}{u^{\prime}}-\frac{S_{16}^{\mu\nu\rho\sigma}}{t^{\prime}}\right)\Biggr{]}$

The form factors $\mathcal{A}_{i}$ associated to the above projectors can be
expressed in terms of the form factors $f_{i}$ defined in eq. (2.8) of ref.
[42] via the following relations

$\displaystyle\mathcal{A}_{1}$ $\displaystyle=0$ (A.21)
$\displaystyle\mathcal{A}_{2}$ $\displaystyle=\frac{p_{\scriptscriptstyle
T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}\frac{2(t^{\prime}-u^{\prime})}{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}\left(p_{\scriptscriptstyle
T}^{2}f_{8}-f_{1}\right)+s^{\prime}p_{\scriptscriptstyle
T}^{2}\left(\frac{u^{\prime}}{t^{\prime}}f_{19}-\frac{t^{\prime}}{u^{\prime}}f_{18}\right)$
(A.22)
$\displaystyle+2t^{\prime}f_{4}-2u^{\prime}f_{5}+2t^{\prime}f_{6}-2u^{\prime}f_{7}+\frac{s^{\prime}t^{\prime}}{u^{\prime}}f_{10}-\frac{s^{\prime}u^{\prime}}{t^{\prime}}f_{11}\Biggr{]}$
$\displaystyle\mathcal{A}_{3}$ $\displaystyle=\frac{p_{\scriptscriptstyle
T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}f_{3}-f_{2}+\frac{\left(m_{\scriptscriptstyle
Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}}{2}\left(\frac{f_{7}}{t^{\prime}}-\frac{f_{4}}{u^{\prime}}\right)+u^{\prime}f_{5}-t^{\prime}f_{6}\Biggr{]}$
(A.23) $\displaystyle\mathcal{A}_{4}$ $\displaystyle=0$ (A.24)
$\displaystyle\mathcal{A}_{5}$ $\displaystyle=\frac{p_{\scriptscriptstyle
T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}\frac{2(t^{\prime}-u^{\prime})}{\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}f_{1}+\frac{s^{\prime}u^{\prime}}{t^{\prime}}f_{11}-\frac{s^{\prime}t^{\prime}}{u^{\prime}}f_{10}\Biggr{]}$
(A.25) $\displaystyle\mathcal{A}_{6}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}f_{3}-f_{2}+t^{\prime}f_{4}-u^{\prime}f_{7}+\frac{\left(m_{\scriptscriptstyle
Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}}{2}\left(\frac{f_{6}}{u^{\prime}}-\frac{f_{5}}{u^{\prime}}\right)\Biggr{]}$
(A.26) $\displaystyle\mathcal{A}_{7}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}f_{3}-f_{2}+\frac{\left(3s^{\prime}m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}s^{\prime}-2\left(u^{\prime}m_{\scriptscriptstyle
Z}^{2}+t^{\prime}\left(m_{\scriptscriptstyle
Z}^{2}+u^{\prime}\right)\right)\right)p_{\scriptscriptstyle
T}^{2}}{4m_{\scriptscriptstyle
Z}^{2}}\left(\frac{f_{5}}{t^{\prime}}-\frac{f_{6}}{u^{\prime}}\right)\Biggr{]}$
(A.27) $\displaystyle\mathcal{A}_{8}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}f_{3}-f_{2}\frac{\left(3s^{\prime}m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle
T}^{2}s^{\prime}-2\left(u^{\prime}m_{\scriptscriptstyle
Z}^{2}+t^{\prime}\left(m_{\scriptscriptstyle
Z}^{2}+u^{\prime}\right)\right)\right)p_{\scriptscriptstyle
T}^{2}}{4m_{\scriptscriptstyle
Z}^{2}}\left(\frac{f_{4}}{u^{\prime}}-\frac{f_{7}}{t^{\prime}}\right)\Biggr{]}$
(A.28) $\displaystyle\mathcal{A}_{9}$
$\displaystyle=\frac{1}{m_{\scriptscriptstyle
Z}^{4}}\Biggl{[}-2p_{\scriptscriptstyle
T}^{4}\left(t^{\prime}f_{4}+u^{\prime}f_{5}+t^{\prime}f_{6}+u^{\prime}f_{7}\right)-\left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)p_{\scriptscriptstyle
T}^{2}\left(f_{2}+f_{3}\right)$ (A.29)
$\displaystyle-\frac{\left(s^{\prime}(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2})(m_{\scriptscriptstyle
Z}^{2}-s^{\prime})+2m_{\scriptscriptstyle Z}^{2}(t^{\prime 2}+u^{\prime
2})\right)}{s^{\prime}}\left(p_{\scriptscriptstyle
T}^{2}f_{8}-f_{1}\right)-s^{\prime}\left(m_{\scriptscriptstyle
Z}^{4}+p_{\scriptscriptstyle T}^{4}\right)\left(f_{9}-p_{\scriptscriptstyle
T}^{2}f_{20}\right)$
$\displaystyle+\frac{2\left(u^{\prime}m_{\scriptscriptstyle
Z}^{2}+t^{\prime}\left(m_{\scriptscriptstyle
Z}^{2}+u^{\prime}\right)\right)}{s^{\prime}}\left(p_{\scriptscriptstyle
T}^{2}\left(t^{\prime 2}f_{18}+u^{\prime 2}f_{19}\right)-\left(t^{\prime
2}f_{10}+u^{\prime 2}f_{11}\right)\right)\Biggr{]}$
$\displaystyle\mathcal{A}_{10}$ $\displaystyle=\frac{2p_{\scriptscriptstyle
T}^{2}}{(m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4})}\Biggl{[}s^{\prime}p_{\scriptscriptstyle
T}^{2}\left(p_{\scriptscriptstyle
T}^{2}\left(\frac{t^{\prime}}{u^{\prime}}f_{18}+\frac{u^{\prime}}{t^{\prime}}f_{19}\right)-\left(\frac{t^{\prime}}{u^{\prime}}f_{10}+\frac{u^{\prime}}{t^{\prime}}f_{11}\right)\right)-2p_{\scriptscriptstyle
T}^{2}\left(t^{\prime}f_{4}+u^{\prime}f_{5}+t^{\prime}f_{6}+u^{\prime}f_{7}\right)$
(A.30) $\displaystyle+\frac{\left(\left(s^{\prime 2}-t^{\prime 2}-u^{\prime
2}\right)m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\left(s^{\prime 2}+t^{\prime 2}+u^{\prime
2}\right)\right)}{\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}}\left(f_{8}p_{\scriptscriptstyle
T}^{2}-f_{1}\right)-\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)\left(f_{2}+f_{3}\right)\Biggr{]}$
$\displaystyle\mathcal{A}_{11}$ $\displaystyle=\frac{1}{m_{\scriptscriptstyle
Z}^{4}}\Biggl{[}s^{\prime}\left(m_{\scriptscriptstyle
Z}^{4}+p_{\scriptscriptstyle
T}^{4}\right)f_{9}+\frac{2\left(u^{\prime}m_{\scriptscriptstyle
Z}^{2}+t^{\prime}\left(m_{\scriptscriptstyle
Z}^{2}+u^{\prime}\right)\right)}{s^{\prime}}\left(t^{\prime 2}f_{10}+u^{\prime
2}f_{11}\right)$ (A.31) $\displaystyle-\frac{\left(2u^{\prime
2}m_{\scriptscriptstyle Z}^{2}+s^{\prime}\left(m_{\scriptscriptstyle
Z}^{2}+u^{\prime}\right)\left(s^{\prime}m_{\scriptscriptstyle
Z}^{2}+2t^{\prime 2}+p_{\scriptscriptstyle
T}^{2}s^{\prime}\right)\right)}{s^{\prime}}f_{1}\Biggr{]}$
$\displaystyle\mathcal{A}_{12}$ $\displaystyle=\frac{2p_{\scriptscriptstyle
T}^{2}}{\left(m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4}\right)}\Biggl{[}\frac{\left(\left(s^{\prime 2}-t^{\prime 2}-u^{\prime
2}\right)m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\left(s^{\prime 2}+t^{\prime 2}+u^{\prime
2}\right)\right)}{\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}}f_{1}+p_{\scriptscriptstyle
T}^{2}s^{\prime}\left(\frac{t^{\prime}}{u^{\prime}}f_{10}+\frac{u^{\prime}}{t^{\prime}}f_{11}\right)\Biggr{]}$
(A.32) $\displaystyle\mathcal{A}_{13}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}\frac{\left(m_{\scriptscriptstyle Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)s^{\prime}}{2}\left(\frac{f_{5}}{t^{\prime}}+\frac{f_{6}}{u^{\prime}}\right)-\left(f_{2}+f_{3}+t^{\prime}f_{4}+u^{\prime}f_{7}\right)\Biggr{]}$
(A.33) $\displaystyle\mathcal{A}_{14}$
$\displaystyle=-\frac{p_{\scriptscriptstyle T}^{2}}{m_{\scriptscriptstyle
Z}^{2}}\Biggl{[}f_{2}+f_{3}+p_{\scriptscriptstyle
T}^{2}s^{\prime}\left(\frac{f_{5}}{t^{\prime}}+\frac{f_{6}}{u^{\prime}}\right)\Biggr{]}$
(A.34) $\displaystyle\mathcal{A}_{15}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}^{2}}{2m_{\scriptscriptstyle
Z}^{4}}\Biggl{[}\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)\left(f_{2}+f_{3}+t^{\prime}f_{4}+u^{\prime}f_{7}\right)-\left(m_{\scriptscriptstyle
Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)\left(u^{\prime}f_{5}+t^{\prime}f_{6}\right)\Biggr{]}$ (A.35)
$\displaystyle\mathcal{A}_{16}$ $\displaystyle=\frac{2p_{\scriptscriptstyle
T}^{2}}{m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}}\Biggl{[}f_{2}+f_{3}+\frac{2p_{\scriptscriptstyle
T}^{2}}{m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}}\left(u^{\prime}f_{5}+t^{\prime}f_{6}\right)\Biggr{]}$ (A.36)
$\displaystyle\mathcal{A}_{17}$ $\displaystyle=\frac{2m_{\scriptscriptstyle
Z}^{2}}{2m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}}\Biggl{[}\frac{\left(s^{\prime 2}\left(t^{\prime 2}+u^{\prime
2}\right)-\left(t^{\prime 2}-u^{\prime 2}\right)^{2}\right)}{s^{\prime
3}\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}f_{1}+f_{2}+f_{3}\Biggr{]}$ (A.37)
$\displaystyle\mathcal{A}_{18}$
$\displaystyle=\frac{\left(m_{\scriptscriptstyle Z}^{2}\left(2s^{\prime
2}-t^{\prime 2}-u^{\prime
2}\right)+s^{\prime}(s^{\prime}(p_{\scriptscriptstyle
T}^{2}-s^{\prime})+t^{\prime 2}+u^{\prime 2})\right)}{m_{\scriptscriptstyle
Z}^{2}s^{\prime}\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)}\left(p_{\scriptscriptstyle T}^{2}f_{8}-f_{1}\right)$ (A.38)

##  B Analytical Results

#### Results for Double-Triangle Diagrams

With reference to eq. (10), we present here the exact results for the double-
triangle contributions to the $\mathcal{A}^{(1,\bowtie)}_{i}$ with
$i=1,\dots,18$. We keep the dependence of the final result on the mass of the
bottom quark, $m_{b}$. We find

$\displaystyle\mathcal{A}^{(1,\bowtie)}_{1}$ $\displaystyle=0$ (B.1)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{2}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}\left(\left(m_{\scriptscriptstyle
Z}^{2}+2t^{\prime}\right)\ \Delta(t^{\prime})^{2}-(t^{\prime}\leftrightarrow
u^{\prime})\right)}{32\sqrt{2}m_{\scriptscriptstyle Z}\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}$ (B.2)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{3}$
$\displaystyle=-\frac{p_{\scriptscriptstyle
T}^{3}\left(u^{\prime}\Delta(t^{\prime})^{2}-(t^{\prime}\leftrightarrow
u^{\prime})\right)}{16\ \sqrt{2}m_{\scriptscriptstyle
Z}\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)^{2}}$
(B.3) $\displaystyle\mathcal{A}^{(1,\bowtie)}_{4}$ $\displaystyle=0$ (B.4)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{5}$
$\displaystyle=-\frac{p_{\scriptscriptstyle T}\ \left(t^{\prime}\
\Delta(t^{\prime})^{2}-(t^{\prime}\leftrightarrow
u^{\prime})\right)}{16\sqrt{2}m_{\scriptscriptstyle Z}\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)}$ (B.5)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{6}$
$\displaystyle=-\frac{p_{\scriptscriptstyle
T}\left(\left(m_{\scriptscriptstyle Z}^{2}+2t^{\prime}\right)u^{\prime 2}\
\Delta(t^{\prime})^{2}-(t^{\prime}\leftrightarrow
u^{\prime})\right)}{16\sqrt{2}m_{\scriptscriptstyle Z}\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)^{2}s^{\prime}}$ (B.6)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{7}$
$\displaystyle=\frac{\left(s^{\prime}\ p_{\scriptscriptstyle T}^{2}+t^{\prime
2}\right)u^{\prime 2}\ \Delta(t^{\prime})^{2}-(t^{\prime}\leftrightarrow
u^{\prime})}{16\sqrt{2}\ \left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)^{2}s^{\prime 2}}$ (B.7)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{8}$
$\displaystyle=\frac{u^{\prime}\left(t^{\prime}(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2})-2s^{\prime}p_{\scriptscriptstyle
T}^{2}\right)\Delta(t^{\prime})^{2}-(t^{\prime}\leftrightarrow
u^{\prime})}{32\sqrt{2}\ \left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)^{2}\ s^{\prime}}$ (B.8)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{9}$
$\displaystyle=\frac{\left(m_{\scriptscriptstyle Z}^{2}-p_{\scriptscriptstyle
T}^{2}\right)\left(\left(m_{\scriptscriptstyle Z}^{2}+2t^{\prime}\right)\
\Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})\right)}{64m_{\scriptscriptstyle Z}^{2}\
\sqrt{m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle T}^{4}}}$ (B.9)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{10}$
$\displaystyle=-\frac{p_{\scriptscriptstyle T}^{2}\
\left(\left(m_{\scriptscriptstyle Z}^{2}+2t^{\prime}\right)\
\Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})\right)}{32\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)\ \sqrt{m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle
T}^{4}}}$ (B.10) $\displaystyle\mathcal{A}^{(1,\bowtie)}_{11}$
$\displaystyle=\frac{u^{\prime}\left(m_{\scriptscriptstyle Z}^{2}+2\
t^{\prime}\right)\left(2s^{\prime}\ m_{\scriptscriptstyle
Z}^{2}+\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)(t^{\prime}+2\
u^{\prime})\right)\Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})}{32m_{\scriptscriptstyle Z}^{2}\ \left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)\sqrt{m_{\scriptscriptstyle
Z}^{4}+p_{\scriptscriptstyle T}^{4}}\ s^{\prime}}$ (B.11)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{12}$
$\displaystyle=-\frac{\left(s^{\prime}m_{\scriptscriptstyle
Z}^{4}+4p_{\scriptscriptstyle T}^{2}s^{\prime}\ m_{\scriptscriptstyle
Z}^{2}+3p_{\scriptscriptstyle T}^{4}s^{\prime}+4p_{\scriptscriptstyle T}^{2}\
\left(t^{\prime 2}+u^{\prime 2}\right)\right)\
\left(\Delta(t^{\prime})^{2}+\Delta(u^{\prime})^{2}\right)}{64\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)\sqrt{m_{\scriptscriptstyle Z}^{4}+p_{\scriptscriptstyle T}^{4}}\
s^{\prime}}$ (B.12) $\displaystyle\mathcal{A}^{(1,\bowtie)}_{13}$
$\displaystyle=\frac{p_{\scriptscriptstyle T}\left(\left(m_{\scriptscriptstyle
Z}^{2}+2t^{\prime}\right)\ u^{\prime
2}\Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})\right)}{16\sqrt{2}m_{\scriptscriptstyle Z}\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)^{2}s^{\prime}}$ (B.13)

$\displaystyle\mathcal{A}^{(1,\bowtie)}_{14}$
$\displaystyle=-\frac{\left(s^{\prime}\ p_{\scriptscriptstyle T}^{2}+t^{\prime
2}\right)u^{\prime 2}\ \Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})}{64\sqrt{2}\ \left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)^{2}s^{\prime 2}}$ (B.14)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{15}$
$\displaystyle=-\frac{p_{\scriptscriptstyle T}\
\left(\left(m_{\scriptscriptstyle Z}^{2}+2t^{\prime}\right)u^{\prime 2}\
\Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})\right)}{16\sqrt{2}m_{\scriptscriptstyle Z}\
\left(m_{\scriptscriptstyle Z}^{2}+p_{\scriptscriptstyle
T}^{2}\right)^{2}s^{\prime}}$ (B.15)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{16}$
$\displaystyle=\frac{\left(s^{\prime}\ p_{\scriptscriptstyle T}^{2}+t^{\prime
2}\right)u^{\prime 2}\ \Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow
u^{\prime})}{16\sqrt{2}\ \left(m_{\scriptscriptstyle
Z}^{2}+p_{\scriptscriptstyle T}^{2}\right)^{2}s^{\prime 2}}$ (B.16)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{17}$ $\displaystyle=0$ (B.17)
$\displaystyle\mathcal{A}^{(1,\bowtie)}_{18}$
$\displaystyle=\frac{1}{64}\left(\left(\frac{s^{\prime}\ p_{\scriptscriptstyle
T}^{2}}{t^{\prime 2}}+1\right)\
\Delta(t^{\prime})^{2}+(t^{\prime}\leftrightarrow u^{\prime})\right)$ (B.18)

where

$\Delta(x)=F(x,m_{t})-F(x,m_{b})$ (B.19)

and

$\displaystyle F(x,M)$ $\displaystyle=\frac{m_{\scriptscriptstyle
Z}^{2}}{x}\left[B_{0}(2x+m_{\scriptscriptstyle
Z}^{2},M^{2},M^{2})-B_{0}(m_{\scriptscriptstyle Z}^{2},M^{2},M^{2})\right]$
(B.20) $\displaystyle+4M^{2}C_{0}(0,m_{\scriptscriptstyle
Z}^{2},2x+m_{\scriptscriptstyle Z}^{2},M^{2},M^{2},M^{2})+2$

#### Results for Triangle Diagrams at NLO

With reference to eq. (10), we present here the exact results for the two-loop
one-particle-irreducible triangle contributions to the NLO form factors. We
obtain

$\displaystyle\mathcal{A}^{(1,\triangle)}_{1}$ $\displaystyle=0$ (B.21)
$\displaystyle\mathcal{A}^{(1,\triangle)}_{2}$
$\displaystyle=\frac{p_{\scriptscriptstyle
T}^{2}(\hat{t}-\hat{u})}{m_{\scriptscriptstyle Z}^{2}(p_{\scriptscriptstyle
T}^{2}+m_{\scriptscriptstyle
Z}^{2})}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.22) $\displaystyle\mathcal{A}^{(1,\triangle)}_{3}$ $\displaystyle=0$ (B.23)
$\displaystyle\mathcal{A}^{(1,\triangle)}_{4}$ $\displaystyle=0$ (B.24)
$\displaystyle\mathcal{A}^{(1,\triangle)}_{5}$
$\displaystyle=-\frac{p_{\scriptscriptstyle
T}^{2}(\hat{t}-\hat{u})}{m_{\scriptscriptstyle Z}^{2}(p_{\scriptscriptstyle
T}^{2}+m_{\scriptscriptstyle
Z}^{2})}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.25) $\displaystyle\mathcal{A}^{(1,\triangle)}_{6}$ $\displaystyle=0$ (B.26)
$\displaystyle\mathcal{A}^{(1,\triangle)}_{7}$ $\displaystyle=0$ (B.27)
$\displaystyle\mathcal{A}^{(1,\triangle)}_{8}$ $\displaystyle=0$ (B.28)
$\displaystyle\mathcal{A}^{(1,\triangle)}_{9}$
$\displaystyle=\frac{(\hat{s}(p_{\scriptscriptstyle
T}^{2}-m_{\scriptscriptstyle Z}^{2})+2m_{\scriptscriptstyle
Z}^{2}(p_{\scriptscriptstyle T}^{2}+m_{\scriptscriptstyle
Z}^{2})}{2m_{\scriptscriptstyle
Z}^{4}}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.29) $\displaystyle\mathcal{A}^{(1,\triangle)}_{10}$
$\displaystyle=\frac{2p_{\scriptscriptstyle T}^{2}(p_{\scriptscriptstyle
T}^{2}\hat{s}+m_{\scriptscriptstyle Z}^{4}-p_{\scriptscriptstyle
T}^{4})}{(p_{\scriptscriptstyle T}^{2}+m_{\scriptscriptstyle
Z}^{2})(p_{\scriptscriptstyle T}^{4}+m_{\scriptscriptstyle
Z}^{4})}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.30) $\displaystyle\mathcal{A}^{(1,\triangle)}_{11}$
$\displaystyle=-\frac{(\hat{s}(p_{\scriptscriptstyle
T}^{2}-m_{\scriptscriptstyle Z}^{2})+2m_{\scriptscriptstyle
Z}^{2}(p_{\scriptscriptstyle T}^{2}+m_{\scriptscriptstyle
Z}^{2})}{2m_{\scriptscriptstyle
Z}^{4}}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.31) $\displaystyle\mathcal{A}^{(1,\triangle)}_{12}$
$\displaystyle=-\frac{2p_{\scriptscriptstyle T}^{2}(p_{\scriptscriptstyle
T}^{2}\hat{s}+m_{\scriptscriptstyle Z}^{4}-p_{\scriptscriptstyle
T}^{4})}{(p_{\scriptscriptstyle T}^{2}+m_{\scriptscriptstyle
Z}^{2})(p_{\scriptscriptstyle T}^{4}+m_{\scriptscriptstyle
Z}^{4})}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.32) $\displaystyle\mathcal{A}^{(1,\triangle)}_{13}$ $\displaystyle=0$
(B.33) $\displaystyle\mathcal{A}^{(1,\triangle)}_{14}$ $\displaystyle=0$
(B.34) $\displaystyle\mathcal{A}^{(1,\triangle)}_{15}$ $\displaystyle=0$
(B.35) $\displaystyle\mathcal{A}^{(1,\triangle)}_{16}$ $\displaystyle=0$
(B.36) $\displaystyle\mathcal{A}^{(1,\triangle)}_{17}$
$\displaystyle=-\frac{2m_{\scriptscriptstyle Z}^{2}}{p_{\scriptscriptstyle
T}^{2}+2m_{\scriptscriptstyle
Z}^{2}}\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.37) $\displaystyle\mathcal{A}^{(1,\triangle)}_{18}$
$\displaystyle=\frac{\hat{s}}{\hat{s}-m_{H}^{2}}\left(C_{F}\mathcal{F}_{1/2}^{(2l)}+C_{A}\mathcal{G}_{1/2}^{(2l,CA)}\right)$
(B.38)

where the functions $\mathcal{F}_{1/2}^{(2l)}$ and
$\mathcal{G}_{1/2}^{(2l,CA)}$ are defined in eqs. (2.11) and (3.8) in ref.
[34].

## References

  * [1] ATLAS Collaboration, “A detailed map of Higgs boson interactions by the ATLAS experiment ten years after the discovery,” Nature 607 no. 7917, (2022) 52–59, arXiv:2207.00092 [hep-ex]. [Erratum: Nature 612, E24 (2022)].
  * [2] CMS Collaboration, A. Tumasyan et al., “A portrait of the Higgs boson by the CMS experiment ten years after the discovery,” Nature 607 no. 7917, (2022) 60–68, arXiv:2207.00043 [hep-ex].
  * [3] G. Heinrich, “Collider Physics at the Precision Frontier,” Phys. Rept. 922 (2021) 1–69, arXiv:2009.00516 [hep-ph].
  * [4] A. Huss, J. Huston, S. Jones, and M. Pellen, “Les Houches 2021—physics at TeV colliders: report on the standard model precision wishlist,” J. Phys. G 50 no. 4, (2023) 043001, arXiv:2207.02122 [hep-ph].
  * [5] ATLAS Collaboration, G. Aad et al., “Evidence of pair production of longitudinally polarised vector bosons and study of CP properties in $ZZ\to 4\ell$ events with the ATLAS detector at $\sqrt{s}=13$ TeV,” arXiv:2310.04350 [hep-ex].
  * [6] ATLAS Collaboration, G. Aad et al., “Measurement of $ZZ$ production cross-sections in the four-lepton final state in $pp$ collisions at $\sqrt{s}=13.6$ TeV with the ATLAS experiment,” arXiv:2311.09715 [hep-ex].
  * [7] R. W. Brown and K. O. Mikaelian, “${W}^{+}{W}^{-}$ and ${Z}^{0}{Z}^{0}$ pair production in ${e}^{+}{e}^{-}$, $\mathrm{pp}$, and $\overline{p}p$ colliding beams,” Phys. Rev. D 19 (Feb, 1979) 922–934. https://link.aps.org/doi/10.1103/PhysRevD.19.922.
  * [8] J. Ohnemus and J. F. Owens, “Order-${\mathrm{\alpha}}_{\mathit{s}}$ calculation of hadronic zz production,” Phys. Rev. D 43 (Jun, 1991) 3626–3639. https://link.aps.org/doi/10.1103/PhysRevD.43.3626.
  * [9] B. Mele, P. Nason, and G. Ridolfi, “Qcd radiative corrections to z boson pair production in hadronic collisions,” Nuclear Physics B 357 no. 2, (1991) 409–438. https://www.sciencedirect.com/science/article/pii/055032139190475D.
  * [10] F. Cascioli, T. Gehrmann, M. Grazzini, S. Kallweit, P. Maierhöfer, A. von Manteuffel, S. Pozzorini, D. Rathlev, L. Tancredi, and E. Weihs, “ZZ production at hadron colliders in NNLO QCD,” Phys. Lett. B 735 (2014) 311–313, arXiv:1405.2219 [hep-ph].
  * [11] T. Gehrmann, A. von Manteuffel, L. Tancredi, and E. Weihs, “The two-loop master integrals for $q\overline{q}\to VV$,” JHEP 06 (2014) 032, arXiv:1404.4853 [hep-ph].
  * [12] F. Caola, J. M. Henn, K. Melnikov, A. V. Smirnov, and V. A. Smirnov, “Two-loop helicity amplitudes for the production of two off-shell electroweak bosons in quark-antiquark collisions,” JHEP 11 (2014) 041, arXiv:1408.6409 [hep-ph].
  * [13] T. Gehrmann, A. von Manteuffel, and L. Tancredi, “The two-loop helicity amplitudes for $q\overline{q}^{\prime}\to{V}_{1}{V}_{2}\to 4$ leptons,” JHEP 09 (2015) 128, arXiv:1503.04812 [hep-ph].
  * [14] E. Accomando, A. Denner, and A. Kaiser, “Logarithmic electroweak corrections to gauge-boson pair production at the LHC,” Nucl. Phys. B 706 (2005) 325–371, arXiv:hep-ph/0409247.
  * [15] A. Bierweiler, T. Kasprzik, and J. H. Kühn, “Vector-boson pair production at the LHC to $\mathcal{O}(\alpha^{3})$ accuracy,” JHEP 12 (2013) 071, arXiv:1305.5402 [hep-ph].
  * [16] J. Baglio, L. D. Ninh, and M. M. Weber, “Massive gauge boson pair production at the LHC: a next-to-leading order story,” Phys. Rev. D 88 (2013) 113005, arXiv:1307.4331 [hep-ph]. [Erratum: Phys.Rev.D 94, 099902 (2016)].
  * [17] M. Grazzini, S. Kallweit, M. Wiesemann, and J. Y. Yook, “$ZZ$ production at the LHC: NLO QCD corrections to the loop-induced gluon fusion channel,” JHEP 03 (2019) 070, arXiv:1811.09593 [hep-ph].
  * [18] D. A. Dicus, C. Kao, and W. W. Repko, “Gluon production of gauge bosons,” Phys. Rev. D 36 (Sep, 1987) 1570–1572. https://link.aps.org/doi/10.1103/PhysRevD.36.1570.
  * [19] E. Glover and J. Van der Bij, “Z-boson pair production via gluon fusion,” Nuclear Physics B 321 no. 3, (1989) 561–590. https://www.sciencedirect.com/science/article/pii/0550321389902629.
  * [20] ATLAS Collaboration, M. Aaboud et al., “Constraints on off-shell Higgs boson production and the Higgs boson total width in $ZZ\to 4\ell$ and $ZZ\to 2\ell 2\nu$ final states with the ATLAS detector,” Phys. Lett. B 786 (2018) 223–244, arXiv:1808.01191 [hep-ex].
  * [21] CMS Collaboration, A. M. Sirunyan et al., “Measurements of the Higgs boson width and anomalous $HVV$ couplings from on-shell and off-shell production in the four-lepton final state,” Phys. Rev. D 99 no. 11, (2019) 112003, arXiv:1901.00174 [hep-ex].
  * [22] F. Caola and K. Melnikov, “Constraining the Higgs boson width with ZZ production at the LHC,” Phys. Rev. D 88 (2013) 054024, arXiv:1307.4935 [hep-ph].
  * [23] J. M. Campbell, R. K. Ellis, and C. Williams, “Bounding the Higgs Width at the LHC Using Full Analytic Results for $gg->e^{-}e^{+}\mu^{-}\mu^{+}$,” JHEP 04 (2014) 060, arXiv:1311.3589 [hep-ph].
  * [24] N. Kauer and G. Passarino, “Inadequacy of zero-width approximation for a light Higgs boson signal,” JHEP 08 (2012) 116, arXiv:1206.4803 [hep-ph].
  * [25] C. Englert and M. Spannowsky, “Limitations and Opportunities of Off-Shell Coupling Measurements,” Phys. Rev. D 90 (2014) 053003, arXiv:1405.0285 [hep-ph].
  * [26] C. Englert, Y. Soreq, and M. Spannowsky, “Off-Shell Higgs Coupling Measurements in BSM scenarios,” JHEP 05 (2015) 145, arXiv:1410.5440 [hep-ph].
  * [27] CMS Collaboration, A. Tumasyan et al., “Measurement of the Higgs boson width and evidence of its off-shell contributions to ZZ production,” Nature Phys. 18 no. 11, (2022) 1329–1334, arXiv:2202.06923 [hep-ex].
  * [28] ATLAS Collaboration, G. Aad et al., “Evidence of off-shell Higgs boson production from $ZZ$ leptonic decay channels and constraints on its total width with the ATLAS detector,” arXiv:2304.01532 [hep-ex].
  * [29] E. Balzani, R. Gröber, and M. Vitti, “Light-quark Yukawa couplings from off-shell Higgs production,” JHEP 10 (2023) 027, arXiv:2304.09772 [hep-ph].
  * [30] U. Haisch and G. Koole, “Off-shell Higgs production at the LHC as a probe of the trilinear Higgs coupling,” JHEP 02 (2022) 030, arXiv:2111.12589 [hep-ph].
  * [31] U. Haisch and G. Koole, “Probing Higgs portals with matrix-element based kinematic discriminants in $ZZ\to 4\ell$ production,” JHEP 04 (2022) 166, arXiv:2201.09711 [hep-ph].
  * [32] U. Haisch, M. Ruhdorfer, K. Schmid, and A. Weiler, “Quantum collider probes of the fermionic Higgs portal,” arXiv:2311.03995 [hep-ph].
  * [33] C. Anastasiou, S. Beerli, S. Bucherer, A. Daleo, and Z. Kunszt, “Two-loop amplitudes and master integrals for the production of a Higgs boson via a massive quark and a scalar-quark loop,” JHEP 01 (2007) 082, arXiv:hep-ph/0611236.
  * [34] U. Aglietti, R. Bonciani, G. Degrassi, and A. Vicini, “Analytic Results for Virtual QCD Corrections to Higgs Production and Decay,” JHEP 01 (2007) 021, arXiv:hep-ph/0611266.
  * [35] R. Harlander and P. Kant, “Higgs production and decay: Analytic results at next-to-leading order QCD,” JHEP 12 (2005) 015, arXiv:hep-ph/0509189.
  * [36] A. von Manteuffel and L. Tancredi, “The two-loop helicity amplitudes for $gg\to V_{1}V_{2}\to 4~{}\mathrm{leptons}$,” JHEP 06 (2015) 197, arXiv:1503.08835 [hep-ph].
  * [37] F. Caola, K. Melnikov, R. Röntsch, and L. Tancredi, “QCD corrections to ZZ production in gluon fusion at the LHC,” Phys. Rev. D 92 no. 9, (2015) 094028, arXiv:1509.06734 [hep-ph].
  * [38] K. Melnikov and M. Dowling, “Production of two Z-bosons in gluon fusion in the heavy top quark approximation,” Phys. Lett. B 744 (2015) 43–47, arXiv:1503.01274 [hep-ph].
  * [39] F. Caola, M. Dowling, K. Melnikov, R. Röntsch, and L. Tancredi, “QCD corrections to vector boson pair production in gluon fusion including interference effects with off-shell Higgs at the LHC,” JHEP 07 (2016) 087, arXiv:1605.04610 [hep-ph].
  * [40] J. M. Campbell, R. K. Ellis, M. Czakon, and S. Kirchner, “Two loop correction to interference in $gg\to ZZ$,” JHEP 08 (2016) 011, arXiv:1605.01380 [hep-ph].
  * [41] R. Gröber, A. Maier, and T. Rauh, “Top quark mass effects in $gg\to ZZ$ at two loops and off-shell Higgs boson interference,” Phys. Rev. D 100 no. 11, (2019) 114013, arXiv:1908.04061 [hep-ph].
  * [42] J. Davies, G. Mishima, M. Steinhauser, and D. Wellmann, “$gg\to ZZ$: analytic two-loop results for the low- and high-energy regions,” JHEP 04 (2020) 024, arXiv:2002.05558 [hep-ph].
  * [43] B. Agarwal, S. P. Jones, and A. von Manteuffel, “Two-loop helicity amplitudes for $gg\to ZZ$ with full top-quark mass effects,” JHEP 05 (2021) 256, arXiv:2011.15113 [hep-ph].
  * [44] C. Brønnum-Hansen and C.-Y. Wang, “Top quark contribution to two-loop helicity amplitudes for $Z$ boson pair production in gluon fusion,” JHEP 05 (2021) 244, arXiv:2101.12095 [hep-ph].
  * [45] B. Agarwal, S. Jones, M. Kerner, and A. von Manteuffel, “Complete NLO QCD Corrections to $ZZ$ Production in Gluon Fusion,” arXiv:2404.05684 [hep-ph].
  * [46] R. Bonciani, G. Degrassi, P. P. Giardino, and R. Gröber, “Analytical Method for Next-to-Leading-Order QCD Corrections to Double-Higgs Production,” Phys. Rev. Lett. 121 no. 16, (2018) 162003, arXiv:1806.11564 [hep-ph].
  * [47] L. Alasfar, G. Degrassi, P. P. Giardino, R. Gröber, and M. Vitti, “Virtual corrections to $gg\to ZH$ via a transverse momentum expansion,” JHEP 05 (2021) 168, arXiv:2103.06225 [hep-ph].
  * [48] L. Bellafronte, G. Degrassi, P. P. Giardino, R. Gröber, and M. Vitti, “Gluon fusion production at NLO: merging the transverse momentum and the high-energy expansions,” JHEP 07 (2022) 069, arXiv:2202.12157 [hep-ph].
  * [49] J. Davies, G. Mishima, K. Schönwald, and M. Steinhauser, “Analytic approximations of $2\to 2$ processes with massive internal particles,” JHEP 06 (2023) 063, arXiv:2302.01356 [hep-ph].
  * [50] F. Caola, J. M. Henn, K. Melnikov, A. V. Smirnov, and V. A. Smirnov, “Two-loop helicity amplitudes for the production of two off-shell electroweak bosons in gluon fusion,” JHEP 06 (2015) 129, arXiv:1503.08759 [hep-ph].
  * [51] T. Hahn, “Generating Feynman diagrams and amplitudes with FeynArts 3,” Comput. Phys. Commun. 140 (2001) 418–431, arXiv:hep-ph/0012260.
  * [52] V. Shtabovenko, R. Mertig, and F. Orellana, “FeynCalc 9.3: New features and improvements,” Comput. Phys. Commun. 256 (2020) 107478, arXiv:2001.04407 [hep-ph].
  * [53] R. Mertig, M. Bohm, and A. Denner, “FEYN CALC: Computer algebraic calculation of Feynman amplitudes,” Comput. Phys. Commun. 64 (1991) 345–359.
  * [54] R. N. Lee, “LiteRed 1.4: a powerful tool for reduction of multiloop integrals,” J. Phys. Conf. Ser. 523 (2014) 012059, arXiv:1310.1145 [hep-ph].
  * [55] R. N. Lee, “Presenting LiteRed: a tool for the Loop InTEgrals REDuction,” arXiv:1212.2685 [hep-ph].
  * [56] G. Passarino and M. J. G. Veltman, “One Loop corrections for $e^{+}e^{-}$ annihilation into $\mu^{+}\mu^{-}$ in the Weinberg Model,” Nucl. Phys. B160 (1979) 151.
  * [57] R. Bonciani, P. Mastrolia, and E. Remiddi, “Vertex diagrams for the QED form-factors at the two loop level,” Nucl. Phys. B 661 (2003) 289–343, arXiv:hep-ph/0301170. [Erratum: Nucl.Phys.B 702, 359–363 (2004)].
  * [58] M. Becchetti and R. Bonciani, “Two-Loop Master Integrals for the Planar QCD Massive Corrections to Di-photon and Di-jet Hadro-production,” JHEP 01 (2018) 048, arXiv:1712.02537 [hep-ph].
  * [59] S. Caron-Huot and J. M. Henn, “Iterative structure of finite loop integrals,” JHEP 06 (2014) 114, arXiv:1404.2922 [hep-th].
  * [60] A. von Manteuffel and L. Tancredi, “A non-planar two-loop three-point function beyond multiple polylogarithms,” JHEP 06 (2017) 127, arXiv:1701.05905 [hep-ph].
  * [61] L. Naterop, A. Signer, and Y. Ulrich, “handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran,” Comput. Phys. Commun. 253 (2020) 107165, arXiv:1909.01656 [hep-ph].
  * [62] R. Bonciani, G. Degrassi, P. P. Giardino, and R. Gröber, “A Numerical Routine for the Crossed Vertex Diagram with a Massive-Particle Loop,” Comput. Phys. Commun. 241 (2019) 122–131, arXiv:1812.02698 [hep-ph].
  * [63] G. Degrassi, P. P. Giardino, and R. Gröber, “On the two-loop virtual QCD corrections to Higgs boson pair production in the Standard Model,” Eur. Phys. J. C 76 no. 7, (2016) 411, arXiv:1603.00385 [hep-ph].
  * [64] R. Gröber, A. Maier, and T. Rauh, “Reconstruction of top-quark mass effects in Higgs pair production and other gluon-fusion processes,” JHEP 03 (2018) 020, arXiv:1709.07799 [hep-ph].
  * [65] J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, M. Spira, and J. Streicher, “Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme,” Eur. Phys. J. C 79 no. 6, (2019) 459, arXiv:1811.05692 [hep-ph].
  * [66] J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca, and M. Spira, “$gg\to HH$ : Combined uncertainties,” Phys. Rev. D 103 no. 5, (2021) 056002, arXiv:2008.11626 [hep-ph].
  * [67] E. Bagnaschi, G. Degrassi, and R. Gröber, “Higgs boson pair production at NLO in the POWHEG approach and the top quark mass uncertainties,” Eur. Phys. J. C 83 no. 11, (2023) 1054, arXiv:2309.10525 [hep-ph].
  * [68] G. Degrassi, R. Gröber, M. Vitti, and X. Zhao, “On the NLO QCD corrections to gluon-initiated ZH production,” JHEP 08 (2022) 009, arXiv:2205.02769 [hep-ph].
  * [69] D. Binosi and L. Theussl, “JaxoDraw: A Graphical user interface for drawing Feynman diagrams,” Comput. Phys. Commun. 161 (2004) 76–86, arXiv:hep-ph/0309015.