#  TTP20-041, P3H-20-074
Virtual corrections to $gg\to ZH$ in the high-energy and large-$m_{t}$ limits

Joshua Daviesa, Go Mishimab, Matthias Steinhauserb
a Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH,
UK
b Department of Physics, Tohoku University, Sendai, 980-8578 Japan
c Institut für Theoretische Teilchenphysik
Karlsruhe Institute of Technology (KIT)
Wolfgang-Gaede Straße 1, 76128 Karlsruhe, Germany

###### Abstract

We compute the next-to-leading order virtual corrections to the partonic
cross-section of the process $gg\to ZH$, in the high-energy and large-$m_{t}$
limits. We use Padé approximants to increase the radius of convergence of the
high-energy expansion in $m_{t}^{2}/s$, $m_{t}^{2}/t$ and $m_{t}^{2}/u$ and
show that precise results can be obtained down to energies which are fairly
close to the top quark pair threshold. We present results both for the form
factors and the next-to-leading order virtual cross-section.

##  1 Introduction

At the CERN Large Hadron Collider (LHC), gluon fusion processes play an
important role due to the large gluon luminosities at high collision energies.
As a consequence one often observes that gluon fusion–induced processes
provide a numerically large contribution to the theory predictions of
production cross-sections. This is true even for processes for which the
leading-order (LO) contribution consists of one-loop diagrams. A prime example
of such a process is inclusive Higgs boson production, where the gluon-fusion
channel is about an order of magnitude larger than all other contributions.

In this paper we consider the associated production of a $Z$ and a Higgs
boson, $pp\to ZH$, often called “Higgs Strahlung”. This process was of
particular importance in the observation of the Higgs boson’s decay into
bottom quarks at ATLAS [1] and CMS [2]. At LO it is mediated by a tree-level
process in which a quark and anti-quark annihilate. For this channel,
corrections up to next-to-next-to-next-to-leading (N3LO) order are available
[3] (see also Ref. [4] and references therein) and are included in the program
vh@nnlo [5, 6]. Electroweak and QCD corrections are also included in the
program HAWK [7].

Associated $ZH$ production can also occur via the loop-induced gluon fusion
process. Although formally of next-to-next-to-leading order (NNLO) with
respect to $pp\to ZH$, it provides a sizeable contribution, in particular in
the boosted Higgs boson regime in which the transverse momentum of the Higgs
boson is large [8, 9]. Furthermore, the process $gg\to ZH$ provides sizeable
contributions to the uncertainties of $ZH$ production with a subsequent decay
of the Higgs boson in a pair of bottom quarks, e.g., Ref. [10]. Being loop-
induced, $gg\to ZH$ is sensitive to physics beyond the Standard Model. In Ref
[11] it was suggested that the gluon-initiated component of $pp\to ZH$ can be
extracted by comparing to $WH$ production, which only has a Drell-Yan–like
component. It is thus important to consider NLO QCD corrections to $gg\to ZH$,
requiring the computation of two-loop box-type Feynman diagrams with two
different final-state masses ($m_{Z}$ and $m_{H}$) and the massive top quark
propagating in the loops.

An exact LO (one-loop) calculation was performed in [12]. At NLO only
approximations in the large $m_{t}$ limit are known [13, 14]. In this work we
consider the two-loop NLO virtual corrections in the high-energy limit,
expanding the two-loop master integrals for $m_{t}^{2}\ll s,t$, where $s$ and
$t$ are the Mandelstam variables. Furthermore, we also provide analytic
results for the form factors in the large-$m_{t}$ limit.111In [14] only the
squared amplitude has been computed. A similar approach has been applied to
the related process of Higgs boson pair production, $gg\to HH$, where a
comparison to exact numerical calculations [15] could be performed and good
agreement was found, even for relatively small values for the Higgs boson
transverse momentum [16]. In Ref. [17] the high-energy expansion was
successfully applied to the top quark contribution of the two-loop diagrams
contributing to $gg\to ZZ$.

The remainder of the paper is organized as follows: In Section 2 we introduce
our notation, and give our definitions for the form factors and the virtual
finite cross-section. In Section 3 we consider the quality of our
approximations by comparing with the exact LO expressions. In Section 4, we
briefly discuss the two-loop form factors in the large-$m_{t}$ limit and in
Section 5 we discuss the form factors in the high-energy limit, and
investigate the behaviour of Padé approximants constructed using this
expansion. In Section 6 we study the NLO virtual finite cross-section and
apply our Padé scheme to extend the approxmation to a larger kinematic region.
Finally in Section 7 we conclude our findings. Auxiliary material can be found
in the Appendix. In Appendix A we present analytic results for the one-
particle reducible double-triangle contribution and in Appendix B we briefly
discuss out treatment of $\gamma_{5}$ and the application of projectors to
obtain the form factors. In Appendix C we present the relations between the
form factors and helicity amplitudes for the process $gg\to ZH$. The
expansions of the form factors are provided in an analytic, computer readable
form in the ancillary files of this paper [18].

##  2 Notation and technicalities

We consider the amplitude for the process $g(p_{1})g(p_{2})\to
Z(p_{3})H(p_{4})$ where all momenta are assumed to be incoming. This leads to
the following definitions for the Mandelstam variables,

$\displaystyle s$ $\displaystyle=$ $\displaystyle(p_{1}+p_{2})^{2}\,,$
$\displaystyle t$ $\displaystyle=$ $\displaystyle(p_{1}+p_{3})^{2}\,,$
$\displaystyle u$ $\displaystyle=$ $\displaystyle(p_{1}+p_{4})^{2}\,,$ (1)

with

$\displaystyle s+t+u$ $\displaystyle=$ $\displaystyle m_{Z}^{2}+m_{H}^{2}\,.$
(2)

Additionally, $p_{1}^{2}=p_{2}^{2}=0$, $p_{3}^{2}=m_{Z}^{2}$ and
$(p_{1}+p_{2}+p_{3})^{2}=m_{H}^{2}$. We also introduce the transverse momentum
($p_{T}$) and the scattering angle ($\theta$) of the final-state bosons, which
are related to the Mandelstam variables as follows,

$\displaystyle t$
$\displaystyle=-\frac{s}{2}\left(1-\beta\cos\theta\right)+\frac{m_{H}^{2}+m_{Z}^{2}}{2}\,,$
$\displaystyle u$
$\displaystyle=-\frac{s}{2}\left(1+\beta\cos\theta\right)+\frac{m_{H}^{2}+m_{Z}^{2}}{2}\,,$
$\displaystyle p_{T}^{2}$ $\displaystyle=\frac{ut-
m_{H}^{2}m_{Z}^{2}}{s}=\frac{s}{4}\beta^{2}\sin^{2}\theta\,,$ (3)

where

$\displaystyle\beta$
$\displaystyle=\sqrt{1-2\frac{m_{Z}^{2}+m_{H}^{2}}{s}+\frac{\left(m_{Z}^{2}-m_{H}^{2}\right)^{2}}{s^{2}}}\,.$
(4)

We denote the polarization vectors of the gluons and the $Z$ boson by
$\varepsilon_{\lambda_{1},\mu}(p_{1})$, $\varepsilon_{\lambda_{2},\nu}(p_{2})$
and $\varepsilon_{\lambda_{3},\rho}(p_{3})$, in terms of which the amplitude
can be written as

$\displaystyle\mathcal{M}_{\lambda_{1},\lambda_{2},\lambda_{3}}$
$\displaystyle=$ $\displaystyle
A^{\mu\nu\rho}\>\varepsilon_{\lambda_{1},\mu}(p_{1})\>\varepsilon_{\lambda_{2},\nu}(p_{2})\>\varepsilon_{\lambda_{3},\rho}(p_{3})\,.$
(5)

Due to charge-conjugation invariance, the vector coupling of the $Z$ boson to
the quarks in the loop does not contribute. The axial-vector coupling is
proportional to the third component of the isospin and thus mass-degenerate
quark doublets also do not contribute. Since we consider the lightest five
quarks to be massless, only contributions from the top-bottom doublet remain.
Furthermore, each individual term of $A^{\mu\nu\rho}$ is proportional to the
totally anti-symmetric $\varepsilon$ tensor from the axial-vector coupling.

Figure 1: LO and NLO Feynman diagrams contributiong to $gg\to ZH$. Curly, wavy
and dashed lines represent gluons, $Z$ bosons and scalar particles (Higgs or
Goldstone bosons), respectively. Solid lines stand for top or bottom quarks.
Note that the latter are only present in the triangle diagrams.

At LO and NLO both triangle- and box-type diagrams have to be considered.
Examples of these diagram classes are shown in Fig. 1. In the box-type
diagrams the Higgs boson couples directly to the quark loop, so diagrams
involving the bottom quark are suppressed by their Yukawa coupling with
respect to diagrams involving the top quark. This is not the case for the
triangle-type diagrams; here contributions from both the top and bottom quark
loops must be considered. At NLO there is also the contribution from the one-
particle reducible double-triangle contribution, shown as first diagram in the
second row of Fig. 1. Note that in the numerical results discussed in the main
part of this paper these contributions are excluded, however, we present exact
analytical results in Appendix A. The contribution from the reducible double-
triangle diagrams to the (differential) partonic cross- section is implemented
in the computer program which comes together with Ref. [14].

In total one can form 60 tensor structures from the indices $\mu,\nu,\rho$,
the independent momenta $p_{1},p_{2}$ and $p_{3}$ and an $\varepsilon$ tensor.
Details of the computation of these 60 structures and our treatment of
$\gamma_{5}$ in $d=4-2\epsilon$ dimensions are given in Appendix B. After
applying transversality conditions
($p_{i}\cdot\varepsilon_{\lambda_{i}}(p_{i})=0$), gauge invariance w.r.t. the
gluons ($p_{1\mu}A^{\mu\nu\rho}=p_{2\nu}A^{\mu\nu\rho}=0$) and Bose symmetry
($A^{\mu\nu\rho}(p_{1},p_{2},p_{3})=A^{\nu\mu\rho}(p_{2},p_{1},p_{3})$), 14
tensor structures remain which can be grouped such that one has to introduce 7
form factors. We follow the decomposition of Ref. [12] and write

$\displaystyle
A_{ab}^{\mu\nu\rho}(p_{1},p_{2},p_{3})=i\delta_{ab}\frac{\sqrt{2}G_{F}M_{Z}}{s}\frac{\alpha_{s}(\mu)}{\pi}\tilde{A}^{\mu\nu\rho}(p_{1},p_{2},p_{3})\,,$
$\displaystyle\tilde{A}^{\mu\nu\rho}(p_{1},p_{2},p_{3})=$ (6)
$\displaystyle\Bigg{\\{}\left(\frac{s}{2}\varepsilon^{\mu\nu\rho\alpha}p_{2\alpha}-p_{2}^{\mu}\varepsilon^{\nu\rho\alpha\beta}p_{1\alpha}p_{2\beta}\right)F_{1}(t,u)-\left(\frac{s}{2}\varepsilon^{\mu\nu\rho\alpha}p_{1\alpha}-p_{1}^{\nu}\varepsilon^{\mu\rho\alpha\beta}p_{1\alpha}p_{2\beta}\right)F_{1}(u,t)$
$\displaystyle\mbox{}+\left(p_{3}^{\mu}+\frac{m_{Z}^{2}-t}{s}p_{2}^{\mu}\right)\varepsilon^{\nu\rho\alpha\beta}p_{2\alpha}\Big{[}p_{1\beta}F_{2}(t,u)+p_{3\beta}F_{3}(t,u)\Big{]}$
$\displaystyle\mbox{}+\left(p_{3}^{\nu}+\frac{m_{Z}^{2}-u}{s}p_{1}^{\nu}\right)\varepsilon^{\mu\rho\alpha\beta}p_{1\alpha}\Big{[}p_{2\beta}F_{2}(u,t)+p_{3\beta}F_{3}(u,t)\Big{]}$
$\displaystyle\mbox{}+\left(\frac{s}{2}\varepsilon^{\mu\nu\rho\alpha}p_{3\alpha}-p_{2}^{\mu}\varepsilon^{\nu\rho\alpha\beta}p_{1\alpha}p_{3\beta}+p_{1}^{\nu}\varepsilon^{\mu\rho\alpha\beta}p_{2\alpha}p_{3\beta}+g^{\mu\nu}\varepsilon^{\rho\alpha\beta\gamma}p_{1\alpha}p_{2\beta}p_{3\gamma}\right)F_{4}(t,u)\Bigg{\\}},$

where $a,b$ are colour indices. Note that while the decomposition is the same,
the form factors in Ref. [12] have dimension $1/\mbox{GeV}^{2}$ whereas we
pull out an overall factor of $1/s$ such that our form factors are
dimensionless. Only $F_{1}$ receives contributions from the triangle-type
diagrams discussed above. We represent the expansion coefficients of the form
factors in the strong coupling constants with the following notation,

$\displaystyle F$ $\displaystyle=$ $\displaystyle
F^{(0)}+\frac{\alpha_{s}(\mu)}{\pi}F^{(1)}+\cdots\,.$ (7)

At this point it is convenient to define new form factors which are linear
combinations of those of Eq. (6),

$\displaystyle F_{12}(t,u)$ $\displaystyle=$ $\displaystyle
F_{1}(t,u)-\frac{t-m_{Z}^{2}}{s}F_{2}(t,u)\,,$ $\displaystyle F_{12}(u,t)$
$\displaystyle=$ $\displaystyle F_{1}(u,t)-\frac{u-m_{Z}^{2}}{s}F_{2}(u,t)\,,$
$\displaystyle F_{2}^{-}(t,u)$ $\displaystyle=$ $\displaystyle
F_{2}(t,u)-F_{2}(u,t)\,,$ $\displaystyle F_{2}^{+}(t,u)$ $\displaystyle=$
$\displaystyle F_{2}(t,u)+F_{2}(u,t)\,.$ (8)

It is easy to see that $F_{2}^{+}(t,u)$ drops out in the squared amplitude and
thus does not contribute to physical quantities. It is furthermore convenient
to introduce

$\displaystyle F_{12}^{-}(t,u)$ $\displaystyle=$ $\displaystyle
F_{12}(t,u)-F_{12}(u,t)\,,$ $\displaystyle F_{12}^{+}(t,u)$ $\displaystyle=$
$\displaystyle F_{12}(t,u)+F_{12}(u,t)\,,$ $\displaystyle F_{3}^{-}(t,u)$
$\displaystyle=$ $\displaystyle F_{3}(t,u)-F_{3}(u,t)\,,$ $\displaystyle
F_{3}^{+}(t,u)$ $\displaystyle=$ $\displaystyle F_{3}(t,u)+F_{3}(u,t)\,,$ (9)

leaving six physically relevant functions:

$\displaystyle F_{12}^{+}(t,u),\quad F_{12}^{-}(t,u),\quad
F_{2}^{-}(t,u),\quad F_{3}^{+}(t,u),\quad F_{3}^{-}(t,u),\quad F_{4}(t,u),$
(10)

where only $F_{12}^{+}(t,u)$ has contributions from triangle-type diagrams.
$F_{k}^{-}(t,u)$ (with $k=12,2,3$) and $F_{4}(t,u)$ are anti-symmetric w.r.t.
the exchange of the arguments $t$,$u$, and $F_{k}^{+}(t,u)$ (with $k=12,3$)
are symmetric. At leading order $F_{3}^{-}(t,u)=0$, however it has non-zero
contributions starting at NLO.

For the computation of the one- and two-loop Feynman diagrams (some examples
are shown in Fig. 1) in the high-energy limit, we proceed as follows: After
producing the amplitude we perform a Taylor expansion in the $Z$ and Higgs
boson masses (since $m_{Z}^{2},m_{H}^{2}\ll m_{t}^{2}$), leaving one- and two-
loop integrals which depend only on $s$, $t$ and $m_{t}$. Using integration-
by-parts reduction techniques with the programs FIRE 6 [19] and LiteRed [20],
these integrals can be reduced to the same basis of 161 two-loop master
integrals as in Refs. [21, 22], in which they were computed as an expansion in
the high-energy limit to order $m_{t}^{32}$. Inserting these expansions into
the amplitude yields its high-energy approximation. We use FORM 4.2 [23] for
most stages of the computation. The calculation is performed in the covariant
$R_{\xi}$ gauge and we allow for a general electroweak gauge parameter
$\xi_{Z}$ which appears in the propagators of the $Z$ boson and Goldstone
boson $\chi$. Thus only the triangle-type diagrams depend on $\xi_{Z}$, and
this dependence cancels upon summing the $Z$ and $\chi$ contributions.

We also investigate the large-$m_{t}$ limit of the form factors. This
expansion is straightforward and proceeds in analogy to [22, 17]. The programs
q2e and exp [24, 25] are used to produce an asymptotic expansion for $m_{t}\gg
q_{1},q_{2},q_{3}$, again performed using FORM, yielding products of massive
vacuum integrals and massless three-point integrals. Results for the $gg\to
ZH$ amplitude, expanded to order $1/m_{t}^{8}$, have been previously published
in [14]; here we provide one additional term in this expansion, to order
$1/m_{t}^{10}$, and provide results at the level of the form factors.

The renormalization of the ultra-violet (UV) divergences proceeds in a similar
way as for the processes $gg\to HH$ [22] and $gg\to ZZ$ [17]. In particular,
we work in the six-flavour theory and renormalize the top quark mass on shell
and the strong coupling $\alpha_{s}$ in the $\overline{\rm MS}$ scheme. In
addition, our treatment of $\gamma_{5}$ (see Appendix B for more details)
requires that we apply additional finite renormalization to the axial and
pseudo-scalar currents [26],

$\displaystyle Z_{5A}$ $\displaystyle=$ $\displaystyle
1-\frac{\alpha_{s}}{\pi}C_{F}+{\cal O}(\epsilon)\,,$ $\displaystyle Z_{5P}$
$\displaystyle=$ $\displaystyle 1-2\frac{\alpha_{s}}{\pi}C_{F}+{\cal
O}(\epsilon)\,.$

The subtraction of the infra-red poles proceeds according to Ref. [27], using
the conventions of Refs. [22, 17]. The subtraction has the form [27]

$\displaystyle F^{(1)}$ $\displaystyle=$ $\displaystyle F^{(1),\rm
IR}-K_{g}^{(1)}F^{(0)}\,\,=\,\,\tilde{F}^{(1)}+\beta_{0}\log\left(\frac{\mu^{2}}{-s-i\delta}\right)F^{(0)}\,,$
(11)

with $\beta_{0}=11C_{A}/12-Tn_{f}/3$ and $F^{(1),\rm IR}$ is UV renormalized
but still infra-red (IR) divergent. $F^{(1)}$ is finite. An explicit
expression for $K_{g}$ is given in Eq. (2.3) of Ref. [22]. After the second
equality sign in Eq. (11) we make the $\mu$ dependence explicit. Note that
below we only need $\tilde{F}^{(1)}$ to construct the NLO cross-section.

In analogy to the process $gg\to HH$ [28] we define the NLO virtual finite
cross-section for $gg\to ZH$ as

$\displaystyle\tilde{{\cal V}}_{\rm fin}$ $\displaystyle=$
$\displaystyle\frac{G_{F}^{2}m_{Z}^{2}}{16s^{2}}\left(\frac{\alpha_{s}}{\pi}\right)^{2}\sum_{\lambda_{1},\lambda_{2},\lambda_{3}}\Bigg{\\{}\left[\tilde{A}_{\rm
sub}^{\mu\nu\rho}\tilde{A}_{\rm
sub}^{\star,\mu^{\prime}\nu^{\prime}\rho^{\prime}}\right]^{(1)}+\frac{C_{A}}{2}\left(\pi^{2}-\log^{2}{\frac{\mu^{2}}{s}}\right)\left[\tilde{A}_{\rm
sub}^{\mu\nu\rho}\tilde{A}_{\rm
sub}^{\star,\mu^{\prime}\nu^{\prime}\rho^{\prime}}\right]^{(0)}\Bigg{\\}}$
(12)
$\displaystyle\mbox{}\qquad\times\varepsilon_{\lambda_{1},\mu}(p_{1})\>\varepsilon^{\star}_{\lambda_{1},\mu^{\prime}}(p_{1})\>\varepsilon_{\lambda_{2},\nu}(p_{2})\>\varepsilon^{\star}_{\lambda_{2},\nu^{\prime}}(p_{2})\>\varepsilon_{\lambda_{3},\rho}(p_{3})\>\varepsilon^{\star}_{\lambda_{3},\rho^{\prime}}(p_{3})\,,$

where the form factors entering the amplitude $\tilde{A}_{\rm
sub}^{\mu\nu\rho}$ are the IR-subtracted finite form factors $\tilde{F}^{(1)}$
of Eq. (11). The superscripts “(0)” and “(1)” after the square brackets in Eq.
(12) indicate that we take the coefficients of $(\alpha_{s}/\pi)^{0}$ and
$(\alpha_{s}/\pi)^{1}$, respectively, of the squared amplitude, in accordance
with Eq. (7). For the discussion in Section 6 it is convenient to introduce
the $\alpha_{s}$-independent quantity

$\displaystyle{\cal V}_{\rm fin}$ $\displaystyle=$
$\displaystyle\frac{\tilde{{\cal V}}_{\rm fin}}{\alpha_{s}^{2}}\,.$ (13)

##  3 Comparison at leading order

In this section we compare our high-energy expansion with the exact LO result.
We implement this by using FeynArts [29] to generate the amplitude and
FormCalc [30] to compute it, performing a Passarino-Veltman reduction to
three- and four-point one-loop scalar integrals. Schouten identities allow us
to write the result in terms of the tensor structures and form factors of Eqs.
(6) and (8). We use Package-X [31] to evaluate the Passarino-Veltman functions
with high precision and to produce an analytic expansion in the limit of a
large top quark mass. We have verified that our implementation of the exact LO
result reproduces the cross-sections provided in the literature [32, 13, 14].
Additionally we have compared the large-$m_{t}$ expansion derived from this
result with a direct expansion of the amplitude as described above.

(a) (b) (c)

Figure 2:  (a) LO squared amplitude and (b) ratio of expansions to the exact
result. Plot (c) is a zoomed-in version of (b). Note that for better
readability in (a) no quartic corrections are shown; the black curve in (a)
refers to the exact result and the blue curve (and associated uncertainty
band) is the result obtained from Padé approximation.

In Fig. 2 we show the squared amplitude at LO, as a function of $\sqrt{s}$,
for a fixed scattering angle $\theta=\pi/2$. The solid black lines correspond
to the exact result. The coloured lines correspond to different expansion
depths in $m_{Z}$ and $m_{H}$, as detailed in the plot legend. All of the
latter include high-energy expansion terms up to $m_{t}^{32}$. The red curve
demonstrates that after including quadratic terms both in $m_{H}$ and $m_{Z}$,
the deviation from the exact result is below 1% (for $\sqrt{s}\;\hbox
to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise
1.0pt\hbox{$>$}\;1000$ GeV). The $m_{H}^{2}$ corrections appear to be more
important than the $m_{Z}^{2}$ corrections. The inclusion of the quartic terms
(shown as green and pink curves, best visible in Fig. 2(c)) improves the
accuracy of the expansion, leading to an almost negligible deviation from the
exact expression. These terms are much harder to compute, however, so in the
NLO results we will restrict the approximation to the quadratic corrections
only. The solid blue curve and associated uncertainty band in Fig. 2(a) shows
the result of a procedure to improve the expansions based on Padé
approximants, which is discussed in more detail in Section 5. Here it is based
on the expansion to quadratic order in $m_{H}$ and $m_{Z}$, confirming that
our computation of the NLO amplitude only to this order is sufficient.

In Fig. 2 and in the following sections, we use the following parameter
values:

$\displaystyle
m_{Z}=91.1876~{}\mbox{GeV},\,m_{t}=172.9~{}\mbox{GeV},\,m_{H}=125.1~{}\mbox{GeV},\,G_{F}=1.16638\\!\times\\!10^{-5}\mbox{GeV}^{-2}.$
(14)

##  4 NLO form factors: large-$m_{t}$ limit

In this section we discuss the large-$m_{t}$ expansion of the form factors at
NLO. While the expansion of $\tilde{F}_{12}^{+}$ starts at $m_{t}^{0}$, we
find that the other form factors of Eq. (10) exhibit some cancellation in the
leading contributions. In particular, $\tilde{F}_{12}^{-}$,
$\tilde{F}_{2}^{-}$, $\tilde{F}_{3}^{+}$ and $\tilde{F}_{4}$ start only at
$1/m_{t}^{4}$, and $\tilde{F}_{3}^{-}$ starts at $1/m_{t}^{6}$.

We now present the leading terms of the large-$m_{t}$ expansion to establish
our notation. For brevity, we restrict ourselves to $\tilde{F}_{12}^{+}$.
Expressions for the expansion of all form factors to $1/m_{t}^{10}$ are
provided in the ancillary files of this paper [18]. Our result reads

$\displaystyle\tilde{F}_{12}^{+,(1)}$
$\displaystyle=-2C_{A}+\frac{s}{m_{t}^{2}}\,\left(\frac{11C_{A}}{72}-\frac{C_{F}}{4}\right)+\frac{1}{m_{t}^{4}}\left(C_{F}\,s\left[\frac{222m_{H}^{2}+162m_{Z}^{2}-263s}{2160}\right]\right.$
$\displaystyle+C_{A}\left(\frac{1}{129600}\left[-75\left\\{m_{H}^{4}+m_{Z}^{4}-\frac{6}{5}m_{H}^{2}m_{Z}^{2}\right\\}+4011m_{H}^{2}s+6713m_{Z}^{2}s\right.\right.$
$\displaystyle\left.\left.\vphantom{\frac{1}{1}}-2646s^{2}+60\left\\{m_{H}^{2}+m_{Z}^{2}-s-t\right\\}t\right]-s\,\log{\left(\frac{-s}{m_{t}^{2}}\right)}\left[\frac{9m_{H}^{2}+26m_{Z}^{2}-9s}{2160}\right]\right)\,,$
(15)

where $C_{A}=3$ and $C_{F}=4/3$ are the quadratic Casimir invariants of
$SU(3)$.

We have verified that after constructing $\tilde{{\cal V}}_{\rm fin}$ in Eq.
(12) we find agreement with the results of Ref. [14] up to order
$1/m_{t}^{8}$.

##  5 NLO form factors: high-energy limit

|
---|---
(a) | (b)
Figure 3:  Results for $\tilde{F}_{3}^{+}$ as a function of $\sqrt{s}$ for
fixed $\theta=\pi/2$ (a) and for fixed $p_{T}$ (b). In (a) dashed and solid
lines correspond to the imaginary and real parts, respectively. The red curves
in (a) represent the Padé results. In (b) only the real part is shown, and the
Padé results are shown as green dashed lines. The high-energy expansions up to
order $m_{t}^{30}$ and $m_{t}^{32}$ are shown in blue. The widely-separated
pair of curves correspond to $p_{T}=400$ GeV.

We now turn to the high-energy expansion of the NLO form factors. The analytic
expressions are large, so we refrain from showing them here but we provide
analytic expressions for all form factors in the ancillary files of this paper
[18]. For illustration we discuss in the following the results for
$\tilde{F}_{3}^{+}$.

In Fig. 3(a) we plot $\tilde{F}_{3}^{+}$ as a function of $\sqrt{s}$. We can
see that the high-energy expansions of both the real (blue solid lines) and
imaginary parts (green dashed lines) converge well for $\sqrt{s}\;\hbox
to0.0pt{\lower 3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise
1.0pt\hbox{$>$}\;800$ GeV, which is in analogy to $gg\to HH$ [22] and $gg\to
ZZ$ [17]. Note that the lighter-coloured curves include fewer $m_{t}$
expansion terms; the darkest lines show the expansion up to $m_{t}^{32}$.

As in previous publications on $gg\to HH$ and $gg\to ZZ$, we make use of Padé
approximants to improve the description of the high-energy expansion. The
methodology used follows that of Section 4 of Ref. [17] and so is not
described in detail, but is only summarized briefly here. The expansion is
used to produce 28 different Padé approximants, which are combined to produce
a central value and error estimate for the approximation. In Fig. 3(a) the
Padé results are shown as red bands. The width of the bands denote the
uncertainty estimates. For regions in which the high-energy expansion
converges, the Padé-based approximation reproduces the expansion. However, it
also produces reliable results for smaller values of $\sqrt{s}$, as can be
expected from the comparison with the LO result shown in Fig. 2(a).

In Fig. 3(b) we show $\tilde{F}_{3}^{+}$ for the fixed values of
$p_{T}=200,\ldots,800$ GeV. The blue (dashed and dotted) curves correspond to
the high-energy expansions and the green (dashed) curves to the Padé results.
For all Padé curves we also show the corresponding uncertainty band, which for
$p_{T}=200$ GeV is relatively large but for $p_{T}=250$ GeV the uncertainty
band is already quite small; it is completely negligible for higher values of
$p_{T}$. Note that the high-energy expansions are only shown for $p_{T}\geq
400$ GeV; for lower $p_{T}$ values the curves lie far outside of the plot
range. For $p_{T}\;\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar
0\relax\sim$}\hss}\raise 1.0pt\hbox{$>$}\;450$ GeV the expansions converge and
are very close to the Padé results. For $p_{T}=400$ GeV, while the expansions
initially agree with each other and the Padé close to $\sqrt{s}=800$ GeV, they
diverge for larger values of $\sqrt{s}$. We recall here that the high-energy
expansion is an expansion in $m_{t}^{2}/s$, $m_{t}^{2}/t$ and $m_{t}^{2}/u$.
For a fixed value of $p_{T}$, increasing $\sqrt{s}$ can lead to values of $t$
or $u$ which are not large enough for convergence.

In Section 6 we will apply the Padé procedure to the virtual finite cross-
section, in order to compare our results with a state-of-the-art numerical
evaluation at NLO [33].

##  6 NLO virtual finite cross-section

Figure 4: ${\cal V}_{\rm fin}$ as a function of $\sqrt{s}$ for eight values of
$p_{T}$, at the scale $\mu^{2}=s$. The uncertainty estimate of the Padé
procedure is displayed as light-coloured bands. For $p_{T}\geq 350$ GeV, we
also show two high-energy expansion curves including terms up to order
$m_{t}^{30}$ and $m_{t}^{32}$.

Our starting point is $\tilde{{\cal V}}_{\rm fin}$ in Eq. (12). For the LO
form factors we use the exact results and for the two-loop form factors we use
the high-energy expansion. This allows us to write ${\cal V}_{\rm fin}$ in Eq.
(13) as an expansion in $m_{t}$ up to order $m_{t}^{32}$.222Note the double-
triangle contribution, which is known analytically (see Appendix A and Ref.
[14]), is not included in our numerical results for ${\cal V}_{\rm fin}$. At
this point we can apply the Padé approximation procedure as described in
Section 5. In Fig. 4 we show ${\cal V}_{\rm fin}$ for several fixed values of
$p_{T}$ as a function of $\sqrt{s}$. For $p_{T}=400$ GeV and larger we also
show two high-energy expansion curves, which include terms up to $m_{t}^{30}$
and $m_{t}^{32}$. These curves agree well with each other and with the Padé
approximation which they produce. For lower values of $p_{T}$, these curves do
not agree with each other, and are not visible within the plot range.

For $p_{T}=150$ GeV Padé procedure produces stable results with an uncertainty
of about 10%. For $p_{T}=200$ GeV the uncertainties are notably smaller and
for higher $p_{T}$ values they are negligible. We advocate to use results
based on our approach for $p_{T}\;\hbox to0.0pt{\lower 3.5pt\hbox{$\mathchar
0\relax\sim$}\hss}\raise 1.0pt\hbox{$>$}\;200$ GeV. For $p_{T}\approx 150$ GeV
the Padé approach provides important results for cross-checks. For lower
values of $p_{T}$ other methods have to be used for the calculation of $gg\to
ZH$, see Ref. [33].

We have compared our results with those of Ref. [33], which are obtained
numerically, but without making any expansions. In the high-energy region, for
104 kinematic points with $p_{T}\geq 200$ GeV we agree to within 2-3%, and for
42 points with $150\leq p_{T}<200$ GeV we agree to within 10% and our values
are consistent to within the uncertainty of the Padé procedure. We also
construct ${\cal V}_{\rm fin}$ using the large-$m_{t}$ expanded NLO form
factors of Section 4 and find that for 120 points with $\sqrt{s}\leq 284$ GeV,
we agree to within 1%. For larger values of $\sqrt{s}$, as expected, the
large-$m_{t}$ expansion starts to diverge significantly from the numerical
results as one approaches the top quark threshold at $\sqrt{s}\approx 346$
GeV.

##  7 Conclusions

In this paper we consider two-loop NLO corrections to the Higgs Strahlung
process $gg\to ZH$. The corresponding Feynman integrals involve five
dimensionful parameters ($s$, $t$, $m_{t}$, $m_{Z}$ and $m_{H}$) which makes
an analytic calculation impossible with the current state-of-the-art
techniques. Numerical computations are also very challenging, however they
have recently been achieved in Ref. [33]. In Section 6 we have discussed the
cross-check of our results against these.

In our approach we use the hierarchy between the top quark, Higgs and $Z$
boson masses and perform an expansion for $m_{t}^{2}\gg m_{Z}^{2},m_{H}^{2}$,
effectively eliminating the dependence on $m_{Z}$ and $m_{H}$ from the
integrals. We show at one-loop order that the expansion converges quickly,
which allows us to truncate the expansion at NLO after the quadratic terms. As
far as the remaining scales are concerned, we concentrate on the high-energy
region where $s,t\gg m_{t}^{2}$. We expand our master integrals in this limit
such that we obtain results for the form factors including expansion terms up
to $m_{t}^{32}$. This allows us to construct, for each phase-space point
(e.g., a particular pair of $\sqrt{s}$, $p_{T}$ values) a set of Padé
approximants, which considerably extend the region of convergence of our
expansion. Our approach provides both central values and uncertainty
estimates.

We provide results for all form factors involved in the $gg\to ZH$ amplitude,
both at one- and two-loop order. Furthermore, we provide relations between the
form factors and helicity amplitudes in Appendix C. The main emphasis of the
paper is the IR subtracted virtual corrections to the partonic cross-section,
which can be combined with other (e.g. Drell-Yan–like) corrections to $gg\to
ZH$. Our method provides precise results for $p_{T}\;\hbox to0.0pt{\lower
3.5pt\hbox{$\mathchar 0\relax\sim$}\hss}\raise 1.0pt\hbox{$>$}\;200$ GeV with
almost negligible uncertainties. For lower values of $p_{T}$ our uncertainties
increase, in particular for smaller values of $\sqrt{s}$. The analytic results
for the form factors obtained in this paper can be obtained in electronic form
from Ref. [18].

Acknowledgements. We thank the authors of Ref. [33] for making their results
available for comparison, prior to publication. The work of J.D. was in part
supported by the Science and Technology Research Council (STFC) under the
Consolidated Grant ST/T00102X/1. The work of G.M. was in part supported by
JSPS KAKENHI (No. JP20J00328). This research was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762
— TRR 257 “Particle Physics Phenomenology after the Higgs Discovery”.

##  Appendix A Double-triangle contribution

In this appendix we present results for the one-particle reducible double-
triangle contribution, as shown in the first diagram in the second row of Fig.
1. We have computed the six form factors in Eq. (10). The results can be
expressed in terms of the functions (see also Ref. [34] where the
corresponding contributions for $gg\to HH$ were computed)

$\displaystyle{B_{0}(x,0,0)}$ $\displaystyle=$ $\displaystyle
2-\log\frac{-x}{\mu^{2}}\,,$ $\displaystyle{B_{0}(x,m^{2},m^{2})}$
$\displaystyle=$ $\displaystyle
2+\beta_{x}\log\frac{\beta_{x}-1}{\beta_{x}+1}-\log\frac{m^{2}}{\mu^{2}}\,,$
$\displaystyle{C_{0}(x,y,0,m^{2},m^{2},m^{2})}$ $\displaystyle=$
$\displaystyle\frac{1}{2(x-y)}\left(\log^{2}\frac{\beta_{x}+1}{\beta_{x}-1}-\log^{2}\frac{\beta_{y}+1}{\beta_{y}-1}\right)\,,$
(16)

with $\beta_{x}=\sqrt{1-4m^{2}/x}$ and $\beta_{y}=\sqrt{1-4m^{2}/y}$. Our
explicit calculation shows that $F_{3}^{+}(t,u)=F_{3}^{-}(t,u)=F_{4}(t,u)=0$.
We furthermore find that $F_{12}(t,u)$ is equal to $F_{2}(t,u)$ with

$\displaystyle F_{2}(t,u)=$
$\displaystyle\frac{8m_{t}^{2}m_{Z}^{2}s}{\left(m_{H}^{2}-u\right)^{2}\left(m_{Z}^{2}-u\right)^{2}}\bigg{\\{}2u\left[B_{0}(m_{H}^{2},m_{t}^{2},m_{t}^{2})-B_{0}(u,m_{t}^{2},m_{t}^{2})\right]$
$\displaystyle+\left(m_{H}^{2}-u\right)\left[\left(-m_{H}^{2}+4m_{t}^{2}+u\right)C_{0}(u,m_{H}^{2},0,m_{t}^{2},m_{t}^{2},m_{t}^{2})+2\right]\bigg{\\}}$
$\displaystyle\quad\times\bigg{[}B_{0}(u,m_{t}^{2},m_{t}^{2})-B_{0}(m_{Z}^{2},m_{t}^{2},m_{t}^{2})+B_{0}(m_{Z}^{2},0,0)-B_{0}(u,0,0)$
$\displaystyle\quad\quad+2m_{t}^{2}\frac{u-m_{Z}^{2}}{m_{Z}^{2}}C_{0}(u,m_{Z}^{2},0,m_{t}^{2},m_{t}^{2},m_{t}^{2})\bigg{]}\,,$
(17)

which can be used to construct $F_{12}^{+}(t,u)$, $F_{12}^{-}(t,u)$ and
$F_{2}^{-}(t,u)$ (cf. Section 2). We note that it is straightforward to expand
the results in Eq. (17) in the large- and small-$m_{t}$ limit.

The contribution of the double-triangle diagrams to the squared matrix element
have been computed in Ref. [14] and implemented in the corresponding computer
program, see Appendix of Ref. [14].

##  Appendix B Projectors and $\gamma_{5}$

Each LO and NLO diagram contributing to the $gg\to ZH$ amplitude contains a
quark with either an axial-vector coupling to a $Z$ boson, or a pseudo-scalar
coupling to a Goldstone boson. The $\gamma_{5}$ matrix present in these
couplings is not defined in the $d=4-2\epsilon$ space-time dimensions of
dimensional regularization. The couplings are re-written in terms of anti-
symmetric tensors, according to [26], as

$\displaystyle\gamma^{\mu}\gamma^{5}$ $\displaystyle=$
$\displaystyle\frac{i}{12}\,\varepsilon^{\mu\nu\rho\sigma}\left(\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}-\gamma_{\sigma}\gamma_{\rho}\gamma_{\nu}\right)\,,$
$\displaystyle\gamma^{5}$ $\displaystyle=$
$\displaystyle\frac{i}{96}\,\varepsilon^{\mu\nu\rho\sigma}\left(\gamma_{\mu}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}+\gamma_{\sigma}\gamma_{\rho}\gamma_{\nu}\gamma_{\mu}-\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}\gamma_{\mu}-\gamma_{\mu}\gamma_{\sigma}\gamma_{\rho}\gamma_{\nu}\right)\,,$
(18)

which allow one to compute the loop integrals in two ways. The first is to
ensure the $\varepsilon$ tensors are not contracted until one can safely work
in four dimensions, by solving tensor loop integrals, performing UV
renormalization and IR subtraction (as detailed in Eq. (11)) and then finally
contracting with the $\varepsilon$ tensors.

An alternative approach, which avoids the need to compute tensor integrals, is
to project the amplitude onto the 60 possible rank-three Lorentz structures
which can be formed from the three independent momenta and the $\varepsilon$
tensor. Since each term of the projectors onto these structures contains an
$\varepsilon$ tensor, the result of their contraction with the r.h.s. of Eq.
(18) can be treated correctly in $d$ dimensions. The 60 projectors act on the
amplitude as

$\displaystyle P_{i}^{\mu_{1}\mu_{2}\mu_{3}}\>A_{\mu_{1}\mu_{2}\mu_{3}}=F_{i}$
(19)

to produce 60 form factors $F_{i}$, which can be reduced to a minimal set (see
the discussion around Eq. (6)). Each projector $P_{i}^{\mu_{1}\mu_{2}\mu_{3}}$
can be written generically as

$\displaystyle
C_{i,1}\varepsilon^{\mu_{1}\mu_{2}\mu_{3}\nu_{1}}q_{1\,\nu_{1}}+C_{i,2}\varepsilon^{\mu_{1}\mu_{2}\mu_{3}\nu_{1}}q_{2\,\nu_{1}}+\cdots+C_{i,60}g^{\mu_{1}\mu_{2}}\varepsilon^{\mu_{3}\nu_{1}\nu_{2}\nu_{3}}q_{1\,\nu_{1}}q_{2\,\nu_{2}}q_{3\,\nu_{3}}\,,$
(20)

and we contract this general form with our amplitude. In the final result, we
specify each of the 60 sets of coefficients $\\{C_{i,j}\\}$ in order to obtain
the 60 form factors $F_{i}$.

##  Appendix C Helicity amplitudes

In this appendix we describe how one can obtain the helicity amplitudes for
the process $gg\to ZH$ from the tensor decomposition which we have introduced
in Section 2. In analogy to Eq. (5) we introduce

$\displaystyle\tilde{\mathcal{M}}_{\lambda_{1},\lambda_{2},\lambda_{3}}$
$\displaystyle=$
$\displaystyle\tilde{A}^{\mu\nu\rho}\varepsilon_{\lambda_{1},\mu}(p_{1})\varepsilon_{\lambda_{2},\nu}(p_{2})\varepsilon_{\lambda_{3},\rho}(p_{3})\,.$
(21)

We furthermore specify the (contravariant) external momenta and polarization
vectors as follows:

$\displaystyle p_{1}=\frac{\sqrt{s}}{2}\left(\begin{array}[]{c}1\\\ 0\\\ 0\\\
1\end{array}\right)\\!,\,\,p_{2}=\frac{\sqrt{s}}{2}\left(\begin{array}[]{c}1\\\
0\\\ 0\\\
-1\end{array}\right)\\!,\,\,p_{3}=\frac{\sqrt{s}}{2}\left(\begin{array}[]{c}-y_{1}\\\
0\\\ -y_{2}\sin\theta\\\
y_{2}\cos\theta\end{array}\right)\\!,\,\,p_{4}=\frac{\sqrt{s}}{2}\left(\begin{array}[]{c}-2+y_{1}\\\
0\\\ y_{2}\sin\theta\\\ -y_{2}\cos\theta\end{array}\right),$ (38)
$\displaystyle\varepsilon_{+}(p_{1})=[\varepsilon_{-}(p_{1})]^{\star}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
i\\\ 1\\\
0\end{array}\right)\\!,\,\,\qquad\varepsilon_{+}(p_{2})=[\varepsilon_{-}(p_{2})]^{\star}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
-i\\\ 1\\\ 0\end{array}\right)\\!,\,\,$ (47)
$\displaystyle\varepsilon_{+}(p_{3})=[\varepsilon_{-}(p_{3})]^{\star}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}0\\\
i\\\ \cos\theta\\\
\sin\theta\end{array}\right)\,,\quad\varepsilon_{0}(p_{3})=\frac{\sqrt{s}}{2m_{Z}}\left(\begin{array}[]{c}y_{2}\\\
0\\\ y_{1}\sin\theta\\\ -y_{1}\cos\theta\end{array}\right),$ (56)
$\displaystyle y_{1}=1+\frac{m_{Z}^{2}-m_{H}^{2}}{s}\,,\qquad
y_{2}=\beta=\sqrt{1-2\frac{m_{Z}^{2}+m_{H}^{2}}{s}+\frac{\left(m_{Z}^{2}-m_{H}^{2}\right)^{2}}{s^{2}}}\,,$
(57)

where $\varepsilon_{0}$ denotes the longitudinal components of polarization
vectors. Recall that all external momenta are defined as incoming, see Eq.
(2). We have chosen the convention for the polarisation vector of $p_{2}$,
following Ref. [35], such that
$\varepsilon_{+}(p_{2})\to\varepsilon_{-}(p_{1})$ in the center-of-momentum
frame. Furthermore, the polarization vectors satisfy

$\displaystyle\sum_{\lambda_{1}}\varepsilon_{\lambda_{1},\mu}\left(p_{1}\right)\varepsilon_{\lambda_{1},\mu^{\prime}}^{\star}\left(p_{1}\right)$
$\displaystyle=-g_{\mu\mu^{\prime}}+\frac{p_{1,\mu}p_{2,\mu^{\prime}}+p_{2,\mu}p_{1,\mu^{\prime}}}{p_{1}\cdot
p_{2}}\,,$
$\displaystyle\sum_{\lambda_{2}}\varepsilon_{\lambda_{2},\nu}\left(p_{2}\right)\varepsilon_{\lambda_{2},\nu^{\prime}}^{\star}\left(p_{2}\right)$
$\displaystyle=-g_{\nu\nu^{\prime}}+\frac{p_{2,\nu}p_{1,\nu^{\prime}}+p_{1,\nu}p_{2,\nu^{\prime}}}{p_{1}\cdot
p_{2}}\,,$
$\displaystyle\sum_{\lambda_{3}}\varepsilon_{\lambda_{3},\rho}\left(p_{3}\right)\varepsilon_{\lambda_{3},\rho^{\prime}}^{\star}\left(p_{3}\right)$
$\displaystyle=-g_{\rho\rho^{\prime}}+\frac{p_{3,\rho}p_{3,\rho^{\prime}}}{m_{Z}^{2}}\,,$
(58)

which means that we have fixed the gauge for the external gauge bosons. With
the above choice, some of the tensor structures in
$A_{ab}^{\mu\nu\rho}\left(p_{1},p_{2},p_{3}\right)$, which are proportional to
either $p_{1}^{\nu}$ or $p_{2}^{\mu}$, are irrelevant because

$\displaystyle\varepsilon_{\lambda_{1}}\left(p_{1}\right)\cdot
p_{2}=\varepsilon_{\lambda_{2}}\left(p_{2}\right)\cdot p_{1}=0\,.$ (59)

This reduces the number of Lorentz structures in Eq. (6) from 14 to 8. We note
that, as we will see in Eq. (63), the dependence on $F_{2}^{+}$ drops out in
the helicity amplitude.

In total there are $2\times 2\times 3=12$ helicity amplitudes. However, due to
various symmetries, only 4 of them are independent. With explicit calculation
we find

$\displaystyle\tilde{\mathcal{M}}_{-\lambda_{1},-\lambda_{2},-\lambda_{3}}=-\tilde{\mathcal{M}}_{\lambda_{1},\lambda_{2},\lambda_{3}}\,,$
$\displaystyle\tilde{\mathcal{M}}_{++-}=-\tilde{\mathcal{M}}_{+++}\Big{|}_{\theta\to\theta+\pi,y_{2}\to-
y_{2}}\,,$
$\displaystyle\tilde{\mathcal{M}}_{+--}=-\tilde{\mathcal{M}}_{+-+}\Big{|}_{\theta\to\theta+\pi,y_{2}\to-
y_{2}}\,,$ (60)

yielding the 4 independent helicity amplitudes. Note that the Mandelstam
variables are invariant under the simultaneous replacements
$\theta\to\theta+\pi,y_{2}\to-y_{2}$ and thus the form factors do not change.

For the evaluation of the $\varepsilon$ tensor, we use the convention that

$\displaystyle\varepsilon^{0123}=+1,$ (61)

and so some of the relevant contractions are as follows:

$\displaystyle\varepsilon^{\mu\rho\alpha\beta}\varepsilon_{\mu+}(p_{1})\varepsilon_{\rho+}(p_{3})p_{1\alpha}p_{2\beta}$
$\displaystyle=-i\frac{s}{4}(1-\cos\theta)\,,$
$\displaystyle\varepsilon^{\mu\rho\alpha\beta}\varepsilon_{\mu+}(p_{1})\varepsilon_{\rho+}(p_{3})p_{1\alpha}p_{3\beta}$
$\displaystyle=i\frac{s}{8}(1-\cos\theta)(y_{1}-y_{2})\,,$
$\displaystyle\varepsilon^{\rho\alpha\beta\gamma}\varepsilon_{\rho+}(p_{3})p_{1\alpha}p_{2\beta}p_{3\gamma}$
$\displaystyle=-i\frac{s\sqrt{s}}{4\sqrt{2}}y_{2}\sin\theta\,,$
$\displaystyle\varepsilon^{\mu\rho\alpha\beta}\varepsilon_{\mu+}(p_{1})\varepsilon_{\rho
0}(p_{3})p_{1\alpha}p_{2\beta}$
$\displaystyle=i\frac{s\sqrt{s}}{4\sqrt{2}m_{Z}}y_{1}\sin\theta\,,$
$\displaystyle\varepsilon^{\mu\rho\alpha\beta}\varepsilon_{\mu+}(p_{1})\varepsilon_{\rho
0}(p_{3})p_{1\alpha}p_{3\beta}$
$\displaystyle=-i\frac{m_{Z}\sqrt{s}}{2\sqrt{2}}\sin\theta\,,$
$\displaystyle\varepsilon^{\mu\nu\rho\alpha}\varepsilon_{\mu+}(p_{1})\varepsilon_{\nu+}(p_{2})\varepsilon_{\rho+}(p_{3})p_{1\alpha}$
$\displaystyle=i\frac{\sqrt{s}}{2\sqrt{2}}\sin\theta\,,$
$\displaystyle\varepsilon^{\mu\nu\rho\alpha}\varepsilon_{\mu+}(p_{1})\varepsilon_{\nu+}(p_{2})\varepsilon_{\rho
0}(p_{3})p_{1\alpha}$
$\displaystyle=-i\frac{s}{4m_{Z}}(y_{1}\cos\theta+y_{2})\,.$ (62)

In terms of the form factors, the 4 independent helicity amplitudes read

$\displaystyle\tilde{\mathcal{M}}_{+-0}$
$\displaystyle=-i\frac{s}{16m_{Z}}y_{2}\sin^{2}\theta(sy_{1}F_{2}^{-}-2m_{Z}^{2}F_{3}^{-})$
$\displaystyle\tilde{\mathcal{M}}_{+-+}$
$\displaystyle=i\frac{s\sqrt{s}}{16\sqrt{2}}y_{2}\sin\theta(1-\cos\theta)(2F_{2}^{-}-y_{1}F_{3}^{-}-y_{2}F_{3}^{+})$
$\displaystyle\tilde{\mathcal{M}}_{+++}$
$\displaystyle=i\frac{s\sqrt{s}}{16\sqrt{2}}\sin\theta\left[4F_{12}^{-}-(y_{1}-y_{2})(2F_{2}^{-}+4F_{4}+y_{2}F_{3}^{-}+y_{2}\cos\theta
F_{3}^{+})\right]$ $\displaystyle\tilde{\mathcal{M}}_{++0}$
$\displaystyle=-i\frac{s}{8m_{Z}}\left[(sy_{1}F_{12}^{-}-2m_{Z}^{2}F_{2}^{-}-4m_{Z}^{2}F_{4})\cos\theta-
sy_{2}F_{12}^{+}+m_{Z}^{2}y_{2}\sin^{2}\theta F_{3}^{+}\right]\,,$ (63)

where we have omitted the arguments of the form factors, $F(t,u)$.

These helicity amplitudes have some general properties which are valid at any
loop order. It is useful to introduce the partial wave decomposition of the
amplitude [35]

$\displaystyle\mathcal{M}_{\lambda_{a},\lambda_{b},\lambda_{c}}=\frac{1}{4\pi}\sum_{J}(2J+1)\bra{\lambda_{c}}S^{J}\ket{\lambda_{a}\lambda_{b}}d^{J}_{\lambda_{a}-\lambda_{b},\lambda_{c}}(\theta)\,,$
(64)

where $\lambda_{a},\lambda_{b}$ ($=\pm 1$ in this case) are the helicities of
initial state particles and $\lambda_{c}$ $(=\pm 1,0)$ is the helicity of the
$Z$ boson. $J$ is the total angular momentum of the system, $S^{J}$ is the
$J$-component of the S-matrix, and $d^{J}_{M,M^{\prime}}(\theta)$ is the
Wigner small-$d$ function (see, e.g. [36]). In the case of
$\mathcal{M}_{++0}$,

$\displaystyle
d^{J}_{\lambda_{a}-\lambda_{b},\lambda_{c}}(\theta)=d^{J}_{0,0}(\theta)=P_{J}(\cos\theta)\,,$
(65)

where $P_{J}(x)$ are the Legendre polynomials, which are even functions for
even $J$, and odd functions for odd $J$. Taking into account that
$F_{12}^{-},F_{2}^{-},F_{3}^{-}$ and $F_{4}$ are anti-symmetric in
$t\leftrightarrow u$ exchange and thus odd functions in $\cos\theta$, whereas
$F_{12}^{+}$ and $F_{3}^{+}$ are symmetric in $t\leftrightarrow u$ exchange
and thus even functions in $\cos\theta$, we find that $\mathcal{M}_{++0}$ is
an even function of $\cos\theta$. Using Eq. (64) and the property of the
Legendre polynomials mentioned above, we conclude that only $J$-even
components components contribute to $\mathcal{M}_{++0}$ and thus this
amplitude is a $J$-even channel. In a similar discussion, it is also
straightforward to show that $\mathcal{M}_{+-0}$, expanded in terms of
$d^{J}_{2,0}(\theta)=P_{J}^{2}(\cos\theta)$ where $P_{J}^{n}$ is the
associated Legendre polynomial, is a $J$-odd channel, and that
$\mathcal{M}_{+++}$, expanded in terms of
$d^{J}_{0,1}(\theta)=P_{J}^{1}(\cos\theta)$, is a $J$-even channel. On the
other hand, $\mathcal{M}_{+-+}$ does not have such a feature. In total, the
squared amplitude should be symmetric in $t\leftrightarrow u$ exchange due to
the fact that the two initial state gluons are indistinguishable, and this
symmetry is made apparent when summing all of the squared helicity amplitudes.

The contribution from the triangle diagrams is present only in
$\mathcal{M}_{++0}$ via $F_{12}^{+}$, and this can be understood in the
following way. Due to the Landau-Yang theorem, the on-shell mode of the
mediating $Z$-boson is forbidden; only the off-shell mode and the Goldstone
boson propagate. Since the off-shell mode and the Goldstone boson behave as
scalars under rotation, they do not appear in $\mathcal{M}_{+-0}$ ($J$-odd) or
$\mathcal{M}_{+-+}$ ($J$-indefinite). Furthermore, using the transversal
condition of the $Z$-boson, $\varepsilon_{\pm}(p_{3})\cdot p_{3}=0$ one can
also derive that $\varepsilon_{\pm}(p_{3})\cdot p_{4}=0$ and thus the
transverse components of the final state $Z$-boson do not couple to the
scalar-scalar-vector interaction (see, e.g., the explicit form of the Feynman
rule for the Goldstone-Higgs-$Z$ boson vertex). Because of this property, the
contribution from triangle diagrams is absent in $\mathcal{M}_{+++}$; only
$\mathcal{M}_{++0}$ (and $\mathcal{M}_{--0}$, due to Eq. (60)) contains such
contributions.

## References

  * [1] M. Aaboud et al. [ATLAS], Phys. Lett. B 786 (2018), 59-86 [arXiv:1808.08238 [hep-ex]].
  * [2] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 121 (2018) no.12, 121801 doi:10.1103/PhysRevLett.121.121801 [arXiv:1808.08242 [hep-ex]].
  * [3] M. C. Kumar, M. K. Mandal and V. Ravindran, JHEP 03 (2015), 037 [arXiv:1412.3357 [hep-ph]].
  * [4] G. Heinrich, [arXiv:2009.00516 [hep-ph]].
  * [5] O. Brein, R. V. Harlander and T. J. E. Zirke, Comput. Phys. Commun. 184 (2013), 998-1003 [arXiv:1210.5347 [hep-ph]].
  * [6] R. V. Harlander, J. Klappert, S. Liebler and L. Simon, JHEP 05 (2018), 089 [arXiv:1802.04817 [hep-ph]].
  * [7] A. Denner, S. Dittmaier, S. Kallweit and A. Mück, Comput. Phys. Commun. 195 (2015), 161-171 [arXiv:1412.5390 [hep-ph]].
  * [8] C. Englert, M. McCullough and M. Spannowsky, Phys. Rev. D 89 (2014) no.1, 013013 [arXiv:1310.4828 [hep-ph]].
  * [9] R. V. Harlander, S. Liebler and T. Zirke, JHEP 02 (2014), 023 [arXiv:1307.8122 [hep-ph]].
  * [10] G. Aad et al. [ATLAS], [arXiv:2007.02873 [hep-ex]].
  * [11] R. V. Harlander, J. Klappert, C. Pandini and A. Papaefstathiou, Eur. Phys. J. C 78 (2018) no.9, 760 [arXiv:1804.02299 [hep-ph]].
  * [12] B. A. Kniehl, Phys. Rev. D 42 (1990) 2253.
  * [13] L. Altenkamp, S. Dittmaier, R. V. Harlander, H. Rzehak and T. J. Zirke, JHEP 02 (2013), 078 doi:10.1007/JHEP02(2013)078 [arXiv:1211.5015 [hep-ph]].
  * [14] A. Hasselhuhn, T. Luthe and M. Steinhauser, JHEP 1701 (2017) 073 [arXiv:1611.05881 [hep-ph]].
  * [15] S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk and T. Zirke, JHEP 1610 (2016) 107 [arXiv:1608.04798 [hep-ph]].
  * [16] J. Davies, G. Heinrich, S. P. Jones, M. Kerner, G. Mishima, M. Steinhauser and D. Wellmann, JHEP 11 (2019), 024 [arXiv:1907.06408 [hep-ph]].
  * [17] J. Davies, G. Mishima, M. Steinhauser and D. Wellmann, JHEP 04 (2020), 024 [arXiv:2002.05558 [hep-ph]].
  * [18] `https://www.ttp.kit.edu/preprints/2020/ttp20-041/`.
  * [19] A. V. Smirnov and F. S. Chuharev, [arXiv:1901.07808 [hep-ph]].
  * [20] R. N. Lee, J. Phys. Conf. Ser. 523 (2014), 012059 [arXiv:1310.1145 [hep-ph]].
  * [21] J. Davies, G. Mishima, M. Steinhauser and D. Wellmann, JHEP 03 (2018), 048 [arXiv:1801.09696 [hep-ph]].
  * [22] J. Davies, G. Mishima, M. Steinhauser and D. Wellmann, JHEP 01 (2019), 176 [arXiv:1811.05489 [hep-ph]].
  * [23] B. Ruijl, T. Ueda and J. Vermaseren, [arXiv:1707.06453 [hep-ph]].
  * [24] R. Harlander, T. Seidensticker and M. Steinhauser, Phys. Lett. B 426 (1998) 125 [hep-ph/9712228].
  * [25] T. Seidensticker, hep-ph/9905298.
  * [26] S. Larin, Phys. Lett. B 303 (1993), 113-118 [arXiv:hep-ph/9302240 [hep-ph]].
  * [27] S. Catani, Phys. Lett. B 427 (1998) 161 [hep-ph/9802439].
  * [28] G. Heinrich, S. P. Jones, M. Kerner, G. Luisoni and E. Vryonidou, JHEP 08 (2017), 088 [arXiv:1703.09252 [hep-ph]].
  * [29] T. Hahn, Comput. Phys. Commun. 140 (2001), 418-431 [arXiv:hep-ph/0012260 [hep-ph]].
  * [30] T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. 118 (1999), 153-165 [arXiv:hep-ph/9807565 [hep-ph]].
  * [31] H. H. Patel, Comput. Phys. Commun. 218 (2017), 66-70 [arXiv:1612.00009 [hep-ph]].
  * [32] O. Brein, M. Ciccolini, S. Dittmaier, A. Djouadi, R. Harlander and M. Kramer, [arXiv:hep-ph/0402003 [hep-ph]].
  * [33] L. Chen, G. Heinrich, S.P. Jones M. Kerner, J. Klappert and J. Schlenk, 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.
  * [34] G. Degrassi, P. P. Giardino and R. Gröber, Eur. Phys. J. C 76 (2016) no.7, 411 doi:10.1140/epjc/s10052-016-4256-9 [arXiv:1603.00385 [hep-ph]].
  * [35] M. Jacob and G. C. Wick, Annals Phys. 7 (1959), 404-428
  * [36] J. J. Sakurai and J. Napolitano, “Modern Quantum Mechanics,” Cambridge University Press.