# Tautological rings of Hilbert modular varieties

Simon Cooper

###### Abstract

We compute the tautological ring for a Hilbert modular variety at an
unramified prime. The method generalises that of van der Geer from the Siegel
case.

##  1 Introduction

In this note we shall compute the tautological ring of Hilbert modular
varieties in characteristic $p$.

For a Shimura variety $S_{K}$, its tautological ring $T^{\bullet}(S_{K})$ is
the subring of its Chow ring with rational coefficients generated by the Chern
classes of automorphic vector bundles. This ring was first defined and studied
by van der Geer [vdG99] in the Siegel case and later by Wedhorn-Ziegler
[Zie18] for Shimura varieties admitting an embedding into a Siegel Shimura
variety (i.e. Hodge-type). Given an integral canonical model of a Shimura
variety, conjectured to exist at hyperspecial level by Milne [Mil92] (and
proven to exist for Hodge- and abelian- type by Kisin [Kis10] and Vasiu
[Vas99]), one can also define the tautological ring for the reduction of a
Shimura variety to characteristic $p$. When suitable toroidal
compactifications $S_{K}^{tor}$ exist there is the folk-lore conjecture that

$T^{\bullet}(S_{K}^{tor})\cong
H^{2\bullet}(\mathbf{X}^{\vee}(\mathbb{C}),\mathbb{Q})$

where $\mathbf{X}^{\vee}$ is the compact dual. Van der Geer proved this in the
Siegel case in all characteristics and Wedhorn-Ziegler proved that it holds
for Shimura varieties of Hodge-type in characteristic $p$ (where suitable
toroidal compactifications of integral models exist by work of Madapusi-Pera
[MP19]). We are interested in the related question:

######  Question 1.1.

Is there a similar description of $T^{\bullet}(S_{K})$?

In the Siegel case the answer is yes. Indeed, van der Geer proved that
$T^{\bullet}(A_{g})\cong H^{2\bullet}({\mathbf{X}_{g-1}}^{\vee})$, where
$\mathbf{X}_{g-1}$ is the double Siegel half-space. In this note we compute
$T^{\bullet}(S_{K})$ for a Hilbert modular variety in characteristic $p$ (see
the beginning of section 5 for a precise statement).

######  Theorem 1.

Let $X_{K}$ be a Hilbert modular variety at an unramified prime.

$T^{\bullet}(X_{K})\cong H^{2\bullet}(\textbf{X}^{\vee})/(\text{top degree
part})$

In particular, the only extra relation is the vanishing of the top Chern class
of the Hodge vector bundle pulled back from the Siegel case. It is not clear
what a general description for Hodge-type Shimura varieties could be. For
example, $S_{K}$ could itself be compact.

The proof of Theorem 1 generalises the method of van der Geer from the Siegel
case, using properness of subvarieties to obtain positivity of certain
classes. Pullback gives a surjective map of graded rings

$H^{2\bullet}(X^{\vee})\cong T^{\bullet}(S_{K}^{tor})\longrightarrow
T^{\bullet}(S_{K})$

Our aim is to determine the kernel of this morphism of graded rings. Section 2
recalls the Hilbert modular datum and in section 3 we review Wedhorn-Ziegler’s
definition of $T^{\bullet}(S_{K})$ and computation of
$T^{\bullet}(S_{K}^{tor})$ for the Hilbert modular datum. In section 4 we
choose an explicit Siegel embedding and explain that the pullback of the Hodge
vector bundle splits into line bundles whose Chern classes generate
$T^{\bullet}(S_{K}^{tor})$. This automatically gives a new relation coming
from the vanishing of the top Chern class of the Hodge vector bundle in the
Siegel case [vdG99, 1.2]. The final section tackles the problem of showing
that there are no further relations. This reduces to an analysis of the
highest degree part, which is at most $d$-dimensional. There are $d$ elements
given by multiples of $d-1$ generators and we show that they are linearly
independent. This requires the properness of certain strata closures with a
description of their cycle classes and is the novel difficulty in the Hilbert
modular case: in the Siegel case the highest degree part is $1$-dimensional so
it sufficed to find a non-zero class.

######  Acknowledgements.

The author would like to thank W. Goldring for his support and guidance and D.
Petersen for helpful discussions. The author is supported by grant KAW
2018.0356.

##  2 Hilbert modular datum

We begin by recalling the Shimura datum which gives a Hilbert modular variety
and describing the cocharacter datum one obtains from this in characteristic
$p$.

Let $F$ be a totally real field of degree $d>1$ over $\mathbb{Q}$. Define
$\mathbf{G}$ as the algebraic group

${\mathbf{G}}$${\mathbb{G}_{m,\mathbb{Q}}}$${Res_{F/\mathbb{Q}}GL_{2}}$${Res_{F/\mathbb{Q}}\mathbb{G}_{m}}$$\scriptstyle{det}$

Define $\mathbf{X}$ as the $\mathbf{G}(\mathbb{R})$-conjugacy class of the
cocharacter $h\colon\mathbb{S}\rightarrow G_{\mathbb{R}}$ given on
$\mathbb{R}$-points $\mathbb{C}^{\times}\rightarrow GL_{2}(\mathbb{R})^{d}$ by

$x+iy\mapsto\left(\begin{pmatrix}x&y\\\ y&-x\end{pmatrix}\right)_{i\leq d}$

$(\mathbf{G},\mathbf{X})$ is a Shimura datum, henceforth referred to as the
Hilbert modular datum given by $F$. The associated Shimura variety is the
Hilbert modular variety and is of PEL-type. Indeed, Raynaud constructed a
model over $\mathbb{Z}$ as a moduli space of principally-polarised abelian
varieties with real multiplication by $\mathcal{O}_{F}$.

######  Definition 2.1.

Let $p$ be a rational prime unramified in $F$. Fix an algebraic closure $k$ of
$\mathbb{F}_{p}$ and a hyperspecial level $K=K^{p}K_{p}$. Denote $X_{K}$ the
special fibre of the integral model of the Hilbert modular variety at level
$K$ and $X_{K,k}$ its base change to $k$. Fix a smooth proper toroidal
compactification $X_{K}^{tor}$, which exists by [MP19].

Let $G$ be the special fibre of a smooth reductive $\mathbb{Z}_{p}$-model of
$\mathbf{G}$. We obtain a cocharacter datum $(G,\mu)$ as in [GK19, 1.2.2]. If
$p$ is inert then $G$ is the fibre product

${G}$${\mathbb{G}_{m,\mathbb{F}_{p}}}$${Res_{\mathbb{F}_{p^{d}}/\mathbb{F}_{p}}GL_{2}}$${Res_{\mathbb{F}_{p^{d}}/\mathbb{F}_{p}}\mathbb{G}_{m}}$$\scriptstyle{det}$

and if $p$ splits then $G$ is the fibre product

${G}$${\mathbb{G}_{m,\mathbb{F}_{p}}}$${(GL_{2,\mathbb{F}_{p}})^{d}}$${(\mathbb{G}_{m,\mathbb{F}_{p}})^{d}}$$\scriptstyle{det}$

Fix an embedding $\mathbb{F}_{p}\hookrightarrow k$ so we can identify
$\\{\tau\colon\mathbb{F}_{p}\hookrightarrow k\\}\cong\mathbb{Z}/d\mathbb{Z}$.
For any unramified prime, the group $G_{\mathbb{F}_{p^{d}}}$ can be described
as the $d$-tuples in $GL_{2}^{d}$ which have the same determinant in each
factor. Any field extension $\mathbb{F}_{p^{d}}\subset K$ contains the field
of definition $\kappa=\mathbb{F}_{p}$ of $\mu$ and

$\mu_{K}\colon\mathbb{G}_{m,K}\longrightarrow G_{K}$

is given by $t\mapsto(\left(\begin{pmatrix}t&0\\\
0&1\end{pmatrix}\right)_{i\leq d},t)$.

Let $P$ be the parabolic subgroup associated with $\mu$ given on $k$-valued
points $P(k)=\\{(\begin{pmatrix}a_{i}&b_{i}\\\
0&d_{i}\end{pmatrix},a_{i}d_{i})\\}\subset G(k)$ (this is called $P_{-}(\mu)$
in [Zie18], $\mathbf{P}^{-}$ in [GK19]). This has a Levi subgroup
$L=Cent(\mu)=\\{(\left(\begin{pmatrix}a_{i}&0\\\
0&d_{i}\end{pmatrix}\right),a_{i}d_{i})\\}$ which is a maximal torus $T$.

##  3 Tautological ring of a toroidal compactification

In this section we review the computation of $T^{\bullet}(X_{K}^{tor})$.

For a Hodge-type Shimura variety $S_{K}$ in characteristic $p$, Wedhorn-
Ziegler [Zie18] proved that $T^{\bullet}(S_{K}^{tor})\cong
A^{\bullet}(\mathop{\mathop{\text{$G$-{\tt Zip}}}\nolimits^{\mu}}\nolimits)$,
where $\mathop{\mathop{\text{$G$-{\tt Zip}}}\nolimits^{\mu}}\nolimits$ is the
stack introduced in [PWZ15] defined from the cocharacter datum $(G,\mu)$ and
classifying certain semi-linear algebraic data. In characteristic $p$ the
morphism $\sigma\colon
S_{K}\rightarrow\mathop{\left[P\backslash*\right]}\nolimits$ given by the
$P$-torsor on $S_{K}$ appearing in the $G$-zip $H^{1}_{dR}$ factors as:

${S_{K}^{tor}}$${\mathop{\mathop{\text{$G$-{\tt
Zip}}}\nolimits^{\mu}}\nolimits}$${\mathop{\left[P\backslash*\right]}\nolimits}$${S_{K}}$$\scriptstyle{\zeta^{tor}}$$\scriptstyle{\beta}$$\scriptstyle{j}$$\scriptstyle{\zeta}$

######  Definition 3.1.

Let $S_{K}$ be a Hodge-type Shimura variety in characteristic $p$. Its
tautological ring is defined $T^{\bullet}(S_{K})\coloneqq
Im(\sigma^{*})=Im(\zeta^{*}\circ\beta^{*})\subset A^{\bullet}(S_{K})$.
Similarly, for a toroidal compactification $T^{\bullet}(S_{K}^{tor})\coloneqq
Im(\zeta^{tor,*}\circ\beta^{*})$.

In general, $\zeta$ is smooth ([Zha18, 3.1.2]) while smoothness of
$\zeta^{tor}$ is expected to follow from [LS18] (and is known in the PEL
case). Wedhorn-Ziegler proved that $\beta^{*}$ is surjective and
$\zeta^{tor,*}$ is injective, giving $T^{\bullet}(S_{K}^{tor})\coloneqq
Im(\zeta^{tor,*}\circ\beta^{*})\cong
A^{\bullet}(\mathop{\mathop{\text{$G$-{\tt Zip}}}\nolimits^{\mu}}\nolimits)$.

Now let $(G,\mu)$ be the cochacter datum associated with $X_{K}$ introduced in
section 2. Then $\pi\colon\mathop{\mathop{\text{$G$-{\tt
ZipFlag}}}\nolimits^{\mu}}\nolimits\longrightarrow\mathop{\mathop{\text{$G$-{\tt
Zip}}}\nolimits^{\mu}}\nolimits$ is an isomorphism, where
$\mathop{\mathop{\text{$G$-{\tt ZipFlag}}}\nolimits^{\mu}}\nolimits$ is the
stack of zip flags defined in [GK19, 2.1]. Thus, on Chow rings

$A^{\bullet}(\mathop{\mathop{\text{$G$-{\tt
Zip}}}\nolimits^{\mu}}\nolimits_{k})\cong
A^{\bullet}(\mathop{\mathop{\text{$G$-{\tt
ZipFlag}}}\nolimits^{\mu}}\nolimits_{k})$

The Weyl group $W=W_{G,T}=\\{\pm 1\\}^{\mathbb{Z}/d\mathbb{Z}}$ acts on the
abelian group $X^{*}(T_{k})$ via its action on $T_{k}$. This extends to an
action on the symmetric algebra $S\coloneqq
Sym(X^{*}(T_{k})\otimes_{\mathbb{Z}}{\mathbb{Q}})$ and we denote
$I=S_{+}^{W}$. Then [Zie18, 8.2] gives

$A^{\bullet}(\mathop{\mathop{\text{$G$-{\tt
ZipFlag}}}\nolimits^{\mu}}\nolimits_{k})=S/IS\cong\mathbb{Q}[z_{1},\ldots,z_{d}]/(z_{1}^{2},\ldots,z_{d}^{2})$

These generators $z_{i}\in S$ are given by characters

$z_{i}\colon\begin{pmatrix}a_{j}&0\\\ 0&d_{j}\end{pmatrix}\mapsto d_{i}$

in $X^{*}(T_{k})$. Thus,

$T^{\bullet}(X_{K,k}^{tor})\cong\frac{\mathbb{Q}[z_{1},\ldots,z_{d}]}{(z_{1}^{2},\ldots,z_{d}^{2})}$

######  Remark 3.1.

This description is essentially independent of the prime. However, one has an
action of Frobenius $\sigma$ on $S$ and this depends on the prime. For
example, if $p$ is inert then $\sigma\cdot z_{i}=z_{i+1}$ for all $i$, whereas
if $p$ is split then the action is trivial. Use $\sigma$ to also refer to the
corresponding permutation of $\\{1,\ldots,d\\}$.

######  Remark 3.2.

It is possible to interpret the relation $z_{i}^{2}$ geometrically. In the
notation of section 5, there is a now-where vanishing section in
$H^{0}(\overline{Z_{w_{i}}},L_{i}\otimes L_{\sigma(i)}^{p})$ giving
$0=c_{1}(L_{i}\otimes
L_{\sigma(i)}^{p}|_{\overline{Z_{w_{i}}}})=(z_{i}+pz_{\sigma(i)})(z_{i}-pz_{\sigma(i)})$
for each $i$.

##  4 A Siegel embedding

######  Notation.

Denote $x_{i}\coloneqq j^{*}(z_{i})\in T^{\bullet}(X_{K,k})$. Denote $V(\rho)$
the vector bundle on $[G\backslash*]$ given by a representation $\rho$ of $G$.

In this section we prove that the relation $x_{1}\ldots x_{d}=0$ holds in
$T^{\bullet}(X_{K,k})$. Fixing a Siegel embedding gives a description of
$z_{i}$ in terms of the Hodge vector bundle. Namely, the pullback of the Hodge
vector bundle $\mathbb{E}$ splits on $X_{K,k}$ into line bundles $\omega_{i}$
with $x_{i}=c_{1}(\omega_{i})$. The relation follows because
$c_{d}(\mathbb{E})=0$ [vdG99, 1.2]. Note that we have made a judicious choice
of embedding to give precisely $x_{i}=c_{1}(\omega_{i})$; this is simply to
ease computation. We will work with the following explicit description of the
general symplectic group involved in the Siegel Shimura datum.

######  Example 4.1.

Consider the symplectic vector space $V=\mathbb{F}_{p}^{2d}$ with symplectic
form given by $J\coloneqq\begin{pmatrix}0&-I_{d}\\\ I_{d}&0\end{pmatrix}$.
Then $GSp(V)=\\{M=\begin{pmatrix}A&B\\\ C&D\end{pmatrix}\mid M^{T}JM=\lambda
J\text{ for some }\lambda\in k\\}$. Let $\mathcal{A}_{d,K^{\prime}}$ the
moduli of principally-polarised abelian varieties with level $K^{\prime}$
hyperspecial at $p$. Then as in section 3 there is a morphism
$\sigma_{d}\colon\mathcal{A}_{d,K^{\prime}}\rightarrow\mathop{\text{$GSp(V)$-{\tt
Zip}}}\nolimits\rightarrow\mathop{\left[P_{d}\backslash*\right]}\nolimits$.

The Siegel parabolic $P_{d}=\\{\begin{pmatrix}A&B\\\ 0&D\end{pmatrix}\\}$ has
a Levi quotient isomorphic to $GL_{d}\times\mathbb{G}_{m}$. The Hodge vector
bundle on $\mathcal{A}_{d,K^{\prime}}$ is
$\mathbb{E}\coloneqq\sigma_{d}^{*}V(\rho_{d})$ given by the representation
$\rho_{d}\colon P_{d}\twoheadrightarrow
GL_{d}\times\mathbb{G}_{m}\xrightarrow{p_{1}}GL_{d}$. Explicitly

$\begin{pmatrix}A&B\\\ 0&D\end{pmatrix}\mapsto D$

Now consider the standard embedding $\iota\colon G\hookrightarrow
GSp(V)=GSp_{2d}$. This can be described by ”embedded squares” over $k$:

$(\begin{pmatrix}a_{i}&b_{i}\\\
c_{i}&d_{i}\end{pmatrix},\lambda)\mapsto\begin{pmatrix}A&B\\\
C&D\end{pmatrix}$

where $A=diag(a_{1},\ldots,a_{d})$, $B=diag(b_{1},\ldots,b_{d})$ etc. The
parabolic $P\subset G$ maps into the Siegel parabolic $P_{d}\subset GSp(V)$.

We can explicitly study the pullback of the Hodge vector bundle to $X_{K,k}$.
The following diagram commutes, giving
$\iota^{*}\mathbb{E}=\sigma^{*}V(\rho)$.

${\mathcal{A}_{d,K^{\prime}}}$${\mathop{\left[P_{d}\backslash*\right]}\nolimits}$${X_{K,k}}$${\mathop{\left[P\backslash*\right]}\nolimits}$

Explicitly, $\rho\colon P\rightarrow P_{d}\rightarrow GL_{d}$ is given by

$\begin{pmatrix}a_{i}&b_{i}\\\ 0&d_{i}\end{pmatrix}\mapsto D$

where $D=diag(d_{1},\ldots d_{d})$. Consider the characters of $P$ given by

$\rho_{j}\colon\begin{pmatrix}a_{i}&b_{i}\\\ 0&d_{i}\end{pmatrix}\mapsto
d_{j}$

Then $\rho=\rho_{1}\oplus\ldots\oplus\rho_{d}$ and so $V(\rho)\simeq
V(\rho_{1})\oplus\ldots\oplus V(\rho_{d})$ as vector bundles on
$\mathop{\left[P\backslash*\right]}\nolimits$. Considering the commutative
diagram

${X^{*}(P)}$${X^{*}(T)}$${A^{\bullet}(\mathop{\left[P\backslash*\right]}\nolimits)}$${S}$${S/IS}$${T^{\bullet}(X_{K,k})}$

we see that $x_{i}=\sigma^{*}c_{1}(V(\rho_{i}))$.

######  Proposition 4.1.

$x_{1}\cdot\ldots\cdot x_{d}=0$ in $T^{\bullet}(X_{K,k})$.

###### Proof.

On $X_{K,k}$ we have $\iota^{*}\mathbb{E}=\sigma^{*}V(\rho)$ so

$x_{1}\ldots
x_{d}=c_{d}(\sigma^{*}V(\rho))=c_{d}(\iota^{*}\mathbb{E})=\iota^{*}c_{d}(\mathbb{E})=0$

with the last equality given in [vdG99, 1.2]. ∎

##  5 Positivity of classes: proof of Theorem 1

In this section we prove that there are no further relations.

######  Theorem 1.

Let $F$ be a totally real field, $p$ an unramified prime, $k$ an algebraic
closure of $\mathbb{F}_{p}$, $K$ a hyperspecial level and $X_{K,k}$ the
Hilbert modular variety at level $K$ in characteristic $p$.

$T^{\bullet}(X_{K,k})=\frac{\mathbb{Q}[x_{1},\ldots,x_{d}]}{(x_{1}^{2},\ldots,x_{d}^{2},x_{1}\cdot\ldots\cdot
x_{d})}$

The proof of Theorem 1 follows from:

######  Proposition 5.1.

Let $p$ be a prime unramified in $F$. The $d$ classes

$\gamma_{i}\coloneqq x_{1}\ldots\widehat{x_{i}}\ldots x_{d}$

are linearly independent in $A^{d-1}(X_{K,k})$.

The closures of Ekedahl-Oort strata
$\overline{X_{K,w}}=\zeta^{-1}(\overline{Z_{w}})$ give cycles in
$T^{\bullet}(X_{K})$ [Zie18, 6.1]. The codimension $1$ strata are those
indexed by $w_{i}=(\epsilon_{j})\in\\{\pm 1\\}^{d}$ with
$\epsilon_{j}=(-1)^{\delta_{ij}+1}$ (length $d-1$). There are sections $s_{i}$
of $L_{i}\otimes L_{\sigma(i)}^{-p}$, where
$L_{i}\coloneqq\beta^{*}V(\rho_{i})$ on $\mathop{\mathop{\text{$G$-{\tt
Zip}}}\nolimits^{\mu}}\nolimits$, such that $Z(s_{i})=\overline{Z_{w_{i}}}$.
To be precise, there is a separating system of partial Hasse invariants for
the ordinary locus $Z_{ord}$ as defined in [GK18, 3.4.2]. The divisors
$D_{i}\coloneqq Z(\zeta^{*}s_{i})=\overline{X_{K,w_{i}}}$ are the Goren-Oort
strata closures defined and shown to be smooth and proper in [GO00] and
[Gor01]. Indeed, one can see by the modular description of $X_{K}^{tor}$ that
$\overline{X_{K,w_{i}}}$ doesn’t intersect the toroidal boundary: a semi-
abelian variety with real multiplication is either abelian or ordinary.

###### Proof of Proposition 5.1.

Denote $\gamma_{i}^{tor}\coloneqq z_{1}\ldots\widehat{z_{i}}\ldots z_{d}$ in
$A^{d-1}(X_{K,k}^{tor})$ so that $\gamma_{i}=j^{*}\gamma_{i}^{tor}$. We obtain
the following commutative diagram.

${X_{K,k}^{tor}}$${D_{i}}$${\mathop{\mathop{\text{$G$-{\tt
Zip}}}\nolimits^{\mu}}\nolimits_{k}}$${X_{K,k}}$$\scriptstyle{\zeta^{tor}}$$\scriptstyle{f_{i}}$$\scriptstyle{j}$$\scriptstyle{\zeta}$

Assume that $\sum_{i=1}^{d}a_{i}\gamma_{i}=0$ in $A^{d-1}(X_{K,k})$. For each
$i$ the projection formula [Ful98, 2.6(b)] for the morphism $f_{i}\circ j$
applied to the class $\sum_{i=1}^{d}a_{i}\gamma_{i}^{tor}$ in
$A^{d-1}(X_{K,k}^{tor})$ gives

$(f_{i}\circ j)_{*}(f_{i}\circ
j)^{*}\sum_{i=1}^{d}a_{i}\gamma_{i}^{tor}=[D_{i}].\sum_{i=1}^{d}a_{i}\gamma_{i}^{tor}$

$(f_{i}\circ
j)^{*}\sum_{i=1}^{d}a_{i}\gamma_{i}^{tor}=f_{i}^{!}\sum_{i=1}^{d}a_{i}\gamma_{i}=0$
by assumption. Since
$[D_{i}]=\zeta^{tor,*}[\overline{Z_{i}}]=z_{i}-pz_{\sigma(i)}$ this implies
that for all $i\in\mathbb{Z}/d\mathbb{Z}$

$\displaystyle 0$ $\displaystyle=[D_{i}].\sum_{i=1}^{d}a_{i}\gamma_{i}^{tor}$
$\displaystyle=(z_{i}-pz_{\sigma(i)}).\sum_{i=1}^{d}a_{i}\gamma_{i}^{tor}$
$\displaystyle=(a_{i}-pa_{\sigma(i)})(z_{1}\ldots z_{d})$

Since $deg(z_{1}\ldots z_{d})=d>0$ we have that $a_{i}-pa_{\sigma(i)}=0$ for
all $i$. It remains to show that the only solution is $a_{i}=0$ for all $i$.
For $p$ inert, $\sigma(i)=i+1$ (modulo $d$) and this system of linear
equations is given by the circulant matrix

$A\coloneqq\begin{pmatrix}1&-p&0&\ldots&0\\\ 0&1&-p&&\\\
\vdots&&\ddots&\ddots&\\\ &&&&\\\ &&&1&-p\\\ -p&&&&1\end{pmatrix}$

A finite sequence of row operations gives an upper triangular matrix with
diagonal entries $1,\ldots,1,1\pm p^{n}$ for some $n\neq 0$. The determinant
of this matrix is non-zero, so the only solution to the system of linear
equations is

$a_{1}=\ldots=a_{d}=0$

A similar argument yields the other cases. Hence,
$\gamma_{1},\ldots,\gamma_{d}$ are linearly independent. ∎

The proof of Theorem 1 can now be completed by showing that there are no
relations in lower degrees (below $d$).

###### Proof of Theorem 1.

$T^{\bullet}(X_{K,k})$ is a quotient of
$\frac{\mathbb{Q}[x_{1},\ldots,x_{d}]}{(x_{1}^{2},\ldots
x_{d}^{2},x_{1}\cdot\ldots\cdot x_{d})}$. The degree of the generators $x_{i}$
is $1$ so the maximum degree is degree $d-1$. Suppose that there is a relation

$\sum\limits_{I\subset\\{1,\ldots,d\\},|I|=l}a_{I}z_{I}=0$

in $T^{l}(X_{K,k})$, indexed over multisets, with $l<d$. For any $I_{0}$ there
exists some $i\notin I_{0}$. Let $J=\\{1,\ldots,i-1,i+1,\ldots,d\\}\setminus
I_{0}$. Multiplying by $z_{J}$ gives $\sum_{I}a_{I}z_{I\sqcup J}=0$ in
$T^{d-1}(X_{K,k})$. This can be rewritten uniquely as $\sum_{I}a_{I}z_{I\sqcup
J}=\sum b_{j}\gamma_{j}$ by Proposition 5.1. If $i$ is not in $I\sqcup J$ then
either $I=I_{0}$ or there is a repeated index so $z_{I\sqcup J}=0$. Hence,
$a_{I_{0}}$ is the only potentially non-zero coefficient of $\gamma_{i}$ and
by linear independence in $T^{d-1}(X_{K,k})$ we have $a_{I_{0}}=0$. Since this
holds for each $I_{0}$ we have that the relation is trivial.

Thus, there are no relations in lower degree and
$T^{\bullet}(X_{K,k})\cong\frac{\mathbb{Q}[x_{1},\ldots,x_{d}]}{(x_{1}^{2},\ldots
x_{d}^{2},x_{1}\cdot\ldots\cdot x_{d})}$

∎

## References

  * [Ful98] William Fulton.  Intersection theory, volume 2 of Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge. A Series of Modern Surveys in Mathematics [Results in Mathematics and Related Areas. 3rd Series. A Series of Modern Surveys in Mathematics].  Springer-Verlag, Berlin, second edition, 1998.
  * [GK18] Wushi Goldring and Jean-Stefan Koskivirta.  Automorphic vector bundles with global sections on $G$-zipZ-schemes.  Compos. Math., 154(12):2586–2605, 2018.
  * [GK19] Wushi Goldring and Jean-Stefan Koskivirta.  Strata Hasse invariants, Hecke algebras and Galois representations.  Invent. Math., 217:887–984, 2019.
  * [GO00] E. Z. Goren and F. Oort.  Stratifications of Hilbert modular varieties.  J. Algebraic Geom., 9:111–154, 2000.
  * [Gor01] Eyal Z. Goren.  Hasse invariants for Hilbert modular varieties.  Israel J. Math., 122:157–174, 2001.
  * [Kis10] Mark Kisin.  Integral models for Shimura varieties of abelian type.  J. Amer. Math. Soc., 23(4):967–1012, 2010.
  * [LS18] Kai-Wen Lan and Benoît Stroh.  Compactifications of subschemes of integral models of Shimura varieties.  Forum Math. Sigma, 6:Paper No. e18, 105, 2018.
  * [Mil92] James S. Milne.  The points on a Shimura variety modulo a prime of good reduction.  In The zeta functions of Picard modular surfaces, pages 151–253. Univ. Montréal, Montreal, QC, 1992.
  * [MP19] Keerthi Madapusi Pera.  Toroidal compactifications of integral models of Shimura varieties of Hodge type.  Ann. Sci. Éc. Norm. Supér. (4), 52(2):393–514, 2019.
  * [PWZ15] Richard Pink, Torsten Wedhorn, and Paul Ziegler.  $F$-zips with additional structure.  Pacific J. Math., 274:183–236, 2015.
  * [Vas99] Adrian Vasiu.  Integral canonical models of Shimura varieties of preabelian type.  Asian J. Math., 3(2):401–518, 1999.
  * [vdG99] Gerard van der Geer.  Cycles on the moduli space of abelian varieties.  In Moduli of curves and abelian varieties, pages 65–89. Friedr. Vieweg, Braunschweig, 1999.
  * [Zha18] Chao Zhang.  Ekedahl-Oort strata for good reductions of Shimura varieties of Hodge type.  Canad. J. Math., 70(2):451–480, 2018.
  * [Zie18] Torsten Wedhorn; Paul Ziegler.  Tautological rings of Shimura varieties and cycle classes of Ekedahl-Oort strata.  preprint, 2018.  http://arxiv.org/abs/1811.04843.


CERN-TH-2022-079

On the NLO QCD Corrections to Gluon-Initiated $ZH$ Production

Giuseppe Degrassia,b***email<EMAIL_ADDRESS>Ramona
Gröberc†††email<EMAIL_ADDRESS>Marco Vittib‡‡‡email:
<EMAIL_ADDRESS>and Xiaoran Zhaob§§§email<EMAIL_ADDRESS>
(a) CERN, Theoretical Physics Department, 1211 Geneva 23, Switzerland
(b) Dipartimento di Matematica e Fisica, Università di Roma Tre and
INFN, sezione di Roma Tre, I-00146 Rome, Italy
(c) Dipartimento di Fisica e Astronomia ’G. Galilei’, Università di Padova and
INFN, sezione di Padova, I-35131 Padova, Italy

We compute the QCD corrections at next-to-leading order for the process
$gg\rightarrow ZH$, including both the virtual two-loop terms and real-
emission contributions. The two-loop box diagrams in the virtual corrections
are approximated analytically over the complete phase space, combining the
results of an expansion in the limit of small transverse momentum and an
expansion in the regime of high energy. We obtain both inclusive and
differential results for the cross section. We find that the NLO QCD
corrections are of the same size as the LO contribution up to $ZH$ invariant
masses close to 1 TeV, but they increase significantly when higher energies
are considered, due to a class of real-emission diagrams in which the $Z$
boson is radiated from an open quark line. Finally, we estimate the
uncertainty due to the renormalization scheme used for the top-quark mass both
on the total and differential cross section.

##  1 Introduction

In the last ten years following the discovery of the Higgs boson [1, 2] a
great experimental effort has been undertaken at the Large Hadron Collider
(LHC) in order to measure the properties of this particle with high precision.
The increased statistics expected from the next LHC runs and from the High-
Luminosity phase will allow for more accurate comparisons with the Standard
Model (SM) predictions, and refined theoretical calculations are needed to
extract as much information as possible from the forthcoming measurements111A
comprehensive review on precision calculations for Higgs physics can be found
in Ref.[3].

The process in which a Higgs is produced together with a weak vector boson,
known as $VH$ associated production ($V=W,Z$), is of great relevance at the
LHC not only as a probe of the couplings between the Higgs and the weak
bosons, but also because of its sensitivity to the $H\rightarrow b\bar{b}$
decay. Indeed, the large QCD backgrounds affecting searches of $H\rightarrow
b\bar{b}$ can be reduced more efficiently in $VH$ than in other production
modes [4, 5, 6]. Recently $VH$ production has been considered also to improve
the constraints on the charm Yukawa coupling [7].

In current $VH$ analyses the impact of theoretical uncertainties from missing
higher-order terms in perturbative calculations depends on the choice of the
final-state boson $V$: while in the $WH$ case these uncertainties are around
1%, and are comparable to the uncertainties from the PDF and $\alpha_{s}$
determination, the corresponding uncertainties for $ZH$ production are
slightly larger, close to 3% [8]. In view of the progress in experimental
precision it is hence essential to improve the theoretical control over
$pp\rightarrow ZH$.

There are two partonic channels that contribute to the $pp\rightarrow ZH$
hadronic cross section. The $q\bar{q}$-initiated (Drell-Yan-like) channel
gives the leading contribution, and higher-order corrections are well under
control as they are known to next-to-next-to-leading order (NNLO) in QCD [9,
10, 11] and to NLO in the EW interactions [12, 13]. The $gg$-initiated channel
has been computed at LO in Refs. [14, 15], and it contributes for the first
time as a NNLO QCD correction to the hadronic process. However, this sub-
leading channel is enhanced by the large gluon luminosity at the LHC, and it
provides about 6% of the total cross section for 13 TeV collisions.
Differential analyses have pointed to significant differences in the shapes of
distributions between the $q\bar{q}$ and the $gg$ channels [16, 17, 18, 19],
showing an increased relative importance of the $gg$ channel in the boosted
regime [20]. The $gg\rightarrow ZH$ process has been considered also as a
probe of new physics effects, with examples including anomalous couplings[19,
20, 21, 22] and new degrees of freedom [23, 20]. Finally, the gluon-initiated
channel has the largest impact on the theoretical uncertainties for $pp\to
ZH$, as only the LO contribution is included in the Monte Carlo programs used
for experimental studies. Since an equivalent gluon-initiated contribution is
not allowed in $WH$ production because of electric charge conservation,
$gg\rightarrow ZH$ is responsible for the larger uncertainties in $ZH$
production compared to $WH$. The reduction of the current theoretical
uncertainties requires the calculation of the NLO QCD corrections to
$gg\rightarrow ZH$, which are the main subject of this paper.

Since the main contribution to $gg\to ZH$ at LO comes from one-loop diagrams
involving loops of massive quarks, the most challenging part of the NLO
calculation is associated to two-loop multi-scale integrals in the virtual
corrections. An estimate of these terms has been obtained in Ref. [24] in the
$m_{t}\rightarrow\infty$ limit and included in the program vh@nnlo [25, 18],
while the effects of a large but finite top-quark mass have been considered in
Refs. [26, 27]. An approximation that retains the effects of the top-quark
mass and that is accurate over the complete phase space has been presented in
Ref. [28], whereas the very recent calculation of Ref. [29] relies on a
combination of the Padé-improved high-energy expansion presented in Ref. [30]
and a numerical evaluation using sector decomposition [31]. All of the
previous results show that the NLO corrections are comparable in magnitude to
the LO contribution222In fact, the NLO corrections can become substantially
larger than the LO contribution in transverse-momentum distributions [29].,
and this fact makes the implementation of the NLO terms in a Monte Carlo code
(possibly interfaced with parton-shower generators [32, 33]) a priority for a
reliable interpretation of the experimental results333Additionally, in Ref.
[34] it has been shown that the resummation of soft-gluon effects can
significantly reduce the scale uncertainties for $gg\rightarrow ZH$.. In this
paper we take a step towards this goal, as we provide a fast and flexible way
to compute $gg\rightarrow ZH$ at NLO in QCD including the effects of a finite
top-quark mass .

We obtain a very reliable approximation of the two-loop virtual corrections by
merging analytic results that are accurate in two complementary phase-space
regions, namely the results of the transverse-momentum expansion of Ref. [35]
and the high-energy (HE) expansion of Ref. [30]. In Ref. [36] it has been
shown that this approach can provide a fast evaluation of the virtual
corrections with an accuracy of 1% or below. In this paper, we present a
complete assessment of the NLO QCD corrections to $gg\rightarrow ZH$ by
including the real-emission contributions, which are related to one-loop
diagrams with an additional parton in the final state. We use our results to
quantify the effects of the gluon-initiated channel on the hadronic cross
section, both at the inclusive and differential level, and we compare our
findings with previous independent calculations. Additionally, we discuss the
theoretical uncertainty related to the choice of renormalization scheme for
the top-quark mass. We finally point out a so-far unnoticed feature of the NLO
QCD corrections, as we observe that the contribution from the class of $2\to
3$ diagrams in which the $Z$ boson is emitted from an open quark line becomes
numerically relevant for collisions at very high energies, namely for
invariant masses in the range $M_{ZH}>1$ TeV. To the best of our knowledge,
the impact of these diagrams on $ZH$ production over this energy regime has
not been discussed in the literature.

This paper is structured as follows: in the next section we set our notation
and we describe the ingredients for the calculation of the NLO corrections to
$gg\rightarrow ZH$; in section 3 we present our results for the inclusive and
differential cross section and we assess the top-mass scheme uncertainty. In
section 4 we state our conclusions.

##  2 The $gg\to ZH$ channel at NLO

The cross section for the subprocess $gg\to ZH+X$ in the hadronic reaction
$pp\to ZH+X$ at center-of-mass energy $\sqrt{s}$, can be written as

$M_{ZH}^{2}\frac{d\,\sigma}{d\,M_{ZH}^{2}}=\sum_{a,b}\int_{0}^{1}dx_{1}dx_{2}\,\,f_{a}(x_{1},\mu_{F}^{2})\,f_{b}(x_{2},\mu_{F}^{2})\int_{0}^{1}dz~{}\delta\left(z-\frac{\tau}{x_{1}x_{2}}\right)M_{ZH}^{2}\frac{d\,\hat{\sigma}_{ab}}{d\,M_{ZH}^{2}}\,,$
(1)

where $M_{ZH}^{2}$ is the invariant mass of the $Z$-Higgs system,
$\tau=M_{ZH}^{2}/s$, $\mu_{F}$ is the factorization scale,
$f_{a}(x,\mu_{F}^{2})$, the parton density of the colliding proton for the
parton of type $a,\,(a=g,q,\bar{q})$ and $\hat{\sigma}_{ab}$ is the partonic
cross section for the subprocess $ab\to ZH+X$ at the partonic center-of-mass
energy $\hat{s}=x_{1}x_{2}s$. The partonic cross section can be written in
terms of the Born (i.e. LO) partonic cross section $\hat{\sigma}^{(0)}$ as:

$M_{ZH}^{2}\frac{d\,\hat{\sigma}_{ab}}{d\,M_{ZH}^{2}}=\hat{\sigma}^{(0)}(z\hat{s})\,z\,G_{ab}(z)\,,$
(2)

where, up to NLO terms in QCD,

$G_{ab}(z)=G_{ab}^{(0)}(z)+\frac{\alpha_{s}(\mu_{R})}{\pi}\,G_{ab}^{(1)}(z)\,$
(3)

with $\alpha_{s}(\mu_{R})$ the strong coupling constant defined at the
renormalization scale $\mu_{R}$.

The LO contribution $G_{ab}^{(0)}(z)$ is given by the $gg\to ZH$ channel only,
i.e.

$G_{ab}^{(0)}(z)=\delta(1-z)\,\delta_{ag}\,\delta_{bg}\,.$ (4)

The amplitude for $g^{\mu}_{a}(p_{1})g^{\nu}_{b}(p_{2})\to
Z^{\rho}(p_{3})H(p_{4})$ can be written as

$\displaystyle{\cal A}=i\sqrt{2}\frac{m_{\scriptscriptstyle
Z}G_{\mu}\alpha_{s}(\mu_{R})}{\pi}\delta_{ab}\epsilon^{a}_{\mu}(p_{1})\epsilon^{b}_{\nu}(p_{2})\epsilon_{\rho}(p_{3})\hat{{\cal
A}}^{\mu\nu\rho}(p_{1},p_{2},p_{3}),$ (5) $\displaystyle\hat{{\cal
A}}^{\mu\nu\rho}(p_{1},p_{2},p_{3})=\sum_{i=1}^{6}\mathcal{P}_{i}^{\mu\nu\rho}(p_{1},p_{2},p_{3}){\cal
A}_{i}(\hat{s},\hat{t},\hat{u},m_{t},m_{\scriptscriptstyle
H},m_{\scriptscriptstyle Z}),$ (6)

where $G_{\mu}$ is the Fermi constant and
$\epsilon^{a}_{\mu}(p_{1})\epsilon^{b}_{\nu}(p_{2})\epsilon_{\rho}(p_{3})$ are
the polarization vectors of the gluons and the $Z$ boson, respectively. The
tensors $\mathcal{P}_{i}^{\mu\nu\rho}$ are a set of orthogonal projectors
whose expressions can be found in Ref. [35]. The corresponding form factors
${\cal A}_{i}(\hat{s},\hat{t},\hat{u},m_{t},m_{\scriptscriptstyle
H},m_{\scriptscriptstyle Z})$ are functions of the masses of the top
quark444We neglect the masses of all quarks except for the top. ($m_{t}$),
Higgs ($m_{\scriptscriptstyle H}$) and $Z$ ($m_{\scriptscriptstyle Z}$)
bosons, and of the partonic Mandelstam variables

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

where $\hat{s}+\hat{t}+\hat{u}=m_{\scriptscriptstyle
Z}^{2}+m_{\scriptscriptstyle H}^{2}$ and we took all the momenta to be
incoming.

The ${\cal A}_{i}$ form factors can be expanded up to NLO terms as

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

with

$\displaystyle\mathcal{A}_{i}^{(0)}$ $\displaystyle=$
$\displaystyle\mathcal{A}_{i}^{(0,\triangle)}+\mathcal{A}_{i}^{(0,\square)},$
(9) $\displaystyle\mathcal{A}_{i}^{(1)}$ $\displaystyle=$
$\displaystyle\mathcal{A}_{i}^{(1,\triangle)}+\mathcal{A}_{i}^{(1,\square)}+\mathcal{A}_{i}^{(1,\bowtie)}~{},$
(10)

where the LO amplitude can be written in terms of two contributions, namely
the one-loop triangle diagrams ($\mathcal{A}_{i}^{(0,\triangle)}$) and the
one-loop box diagrams ($\mathcal{A}_{i}^{(0,\square)}$), while at NLO besides
the genuinely two-loop triangles ($\mathcal{A}_{i}^{(1,\triangle)}$) and the
two-loop boxes ($\mathcal{A}_{i}^{(1,\square)}$) also the contribution from
reducible double-triangle diagrams ($\mathcal{A}_{i}^{(1,\bowtie)}$) is
present. In the latter we include also diagrams featuring loops of bottom
quarks connecting to the $Z$ boson, and we set the bottom mass to zero. In Eq.
(10) the $\mathcal{A}_{i}^{(1,\triangle)}$ and $\mathcal{A}_{i}^{(1,\square)}$
contributions are understood as regularized with respect to the ultraviolet
(UV) and infrared (IR) singularities via the introduction of a counterterm as
in Ref. [27].

The Born partonic cross section can be written as

$\hat{\sigma}^{(0)}(\hat{s})=\frac{m_{\scriptscriptstyle
Z}^{2}G_{\mu}^{2}\alpha_{s}(\mu_{R})^{2}}{32\,\hat{s}\,\pi^{2}}\int{\rm
d}\Phi\,\sum_{i}\left|{\cal A}_{i}^{(0)}\right|^{2},$ (11)

where $\rm{d}\Phi$ is the two-particle Lorentz-invariant phase space.

The NLO terms $G_{ab}^{(1)}$ in Eq. (3) include, besides the $gg$ channel,
also the $2\to 3$ processes $gq\rightarrow ZHq$, $\bar{q}g\rightarrow
ZH\bar{q}$ and $q\bar{q}\rightarrow ZHg$. The NLO contribution to the $gg$
channel involves the two-loop virtual corrections to $gg\rightarrow ZH$
discussed above and one-loop real corrections from $gg\to ZHg$. As well known,
the individual contributions are IR divergent while their sum is finite. In
this work, at the level of cross section, the IR singularities in all channels
were treated via the dipole subtraction method [37]. The outcome of this
procedure for the $gg$ channel can be summarized as follows

$\displaystyle G_{gg}^{(1)}(z)$ $\displaystyle=$
$\displaystyle\delta(1-z)\left[C_{A}\,\frac{~{}\pi^{2}}{3}\,+\beta_{0}\,\ln\left(\frac{\mu_{R}^{2}}{\mu_{F}^{2}}\right)\,+{\cal
C}_{\scriptscriptstyle{\rm NLO}}\right]$ (12) $\displaystyle+$ $\displaystyle
P_{gg}(z)\,\ln\left(\frac{\hat{s}}{\mu_{F}^{2}}\right)+C_{A}\,\frac{4}{z}\,(1-z+z^{2})^{2}\,{\cal
D}_{1}(z)+{\cal R}_{gg}\,,$

where $C_{A}=N_{c}$ ($N_{c}$ being the number of colors),
$\beta_{0}=11/6\,C_{A}-1/2\,N_{f}\,C_{F}\,T_{R}$ ($N_{f}$ being the number of
active flavors, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and $T_{R}=1/2$) is the one-
loop $\beta$-function of the strong coupling in the SM, $P_{gg}$ is the LO
Altarelli-Parisi splitting function

$P_{gg}(z)~{}=~{}2\,C_{A}\,\left[{\cal
D}_{0}(z)+\frac{1}{z}-2+z(1-z)\right]\,,$ (13)

and

${\cal D}_{i}(z)=\left[\frac{\ln^{i}(1-z)}{1-z}\right]_{+}\,$ (14)

where the plus distribution is used.

The first line of Eq. (12) displays the two-loop virtual contribution, with

${\cal C}_{\scriptscriptstyle{\rm NLO}}=\,\frac{\int{\rm
d}\Phi\,2\,\sum_{i}\text{Re}\left[{\cal
A}_{i}^{(0)}\,\left(\mathcal{A}_{i}^{(1)}\right)^{*}\right]}{\int{\rm
d}\Phi\,\sum_{i}\left|{\cal A}_{i}^{(0)}\right|^{2}}\,.$ (15)

In the second line of Eq. (12), the term ${\cal R}_{gg}$ contains the
integration over the three-particle phase space of the $gg\to ZHg$ squared
amplitude minus two dipole subtractions555All the squared matrix elements are
understood as averaged over the initial-state spin, helicity and color and
summed over the final ones. which cure the IR singularities when the final-
state gluon becomes soft or becomes collinear with either of the initial-state
gluons. According to our normalization (see Eqs.(2, 3)) ${\cal R}_{gg}$ is
obtained dividing the latter quantity by
$\alpha_{s}(\mu_{R})\sigma^{(0)}(z\hat{s})\,z/\pi$.

The other NLO contributions to $G_{ab}$, i.e. the $gq\rightarrow
ZHq,\,\bar{q}g\to ZH\bar{q}$ and $q\bar{q}\to ZHg$ channels, require a little
discussion. The LO contribution in the hadronic reaction $pp\to ZH+X$ is given
by the tree-level Drell-Yan-like process $q\bar{q}\to Z^{*}\to ZH$. The NLO
real contribution, ${\cal O}(\alpha_{s})$, to this channel includes the tree-
level process $q\bar{q}\to ZHg$ and the crossed channels. The interference
between these tree-level diagrams and their one-loop corrections is usually
considered as an NNLO correction, ${\cal O}(\alpha_{s}^{2})$, to the Drell-
Yan-like contribution and taken into account in the NNLO evaluation of $pp\to
ZH+X$.
The NLO real corrections to the $gg\to ZH$ process are formally an
$\text{N}^{3}$LO contribution, ${\cal O}(\alpha_{s}^{3})$, to the $pp\to ZH+X$
reaction. In this paper we identify them as the square of the one-loop
diagrams containing a closed fermion loop to which a Higgs or a $Z$ boson, or
both particles, are attached in the $gq\rightarrow ZHq,\,\bar{q}g\to
ZH\bar{q}$ and $q\bar{q}\to ZHg$ channels. It should be remarked that our
definition of the NLO contribution to $gg\to ZH$ process differs from the one
employed in Refs. [28, 29]. In those references, diagrams in which a $Z$ boson
is emitted from an open fermion line, which we denote as $Z$-_radiated
diagrams_ for brevity (see subsection 2.2), were not taken into account
because they were assigned to uncalculated $\text{N}^{3}$LO corrections to the
Drell-Yan-like contribution.

The contribution of the $qg\to ZHq$ channel, and similarly for $\bar{q}g\to
ZH\bar{q}$, can be written as:

$G_{qg}^{(1)}(z)=P_{gq}(z)\left[\ln(1-z)+\frac{1}{2}\ln\left(\frac{\hat{s}}{\mu_{F}^{2}}\right)\right]+\frac{1}{2}C_{F}z+{\cal
R}_{qg}\,,$ (16)

where

$P_{gq}(z)~{}=~{}C_{F}\,\frac{1+(1-z)^{2}}{z}~{},$ (17)

and ${\cal R}_{qg}$ is obtained from the integration over the three-particle
phase space of the squared amplitude of the process minus one dipole
subtraction, which cures the IR singularity when the final quark becomes
collinear to the initial quark, normalized to
$\alpha_{s}(\mu_{R})\sigma^{(0)}(z\hat{s})\,z/\pi$.

Finally, the $q\bar{q}\to ZHg$ channel is IR safe. Its contribution can be
written as

$G_{q\bar{q}}^{(1)}(z)={\cal R}_{q\bar{q}}\,,$ (18)

where ${\cal R}_{q\bar{q}}$ is given by the square of the one-loop diagrams
integrated over the phase space and normalized to
$\alpha_{s}(\mu_{R})\sigma^{(0)}(z\hat{s})\,z/\pi$.

In the following subsections we discuss the method used for the calculation
and the implementation of the various contributions to Eq. (3). We evaluated
our results using different renormalization schemes for the top-quark mass,
namely the on-shell (OS) and the modified minimal subtraction
($\overline{\mathrm{MS}}$) scheme. In particular, for evaluating the top mass
in the $\overline{\mathrm{MS}}$ scheme, we first convert the OS mass to
$m_{t}^{\overline{\mathrm{MS}}}(\mu_{t}=m_{t}^{\textrm{OS}})$ using the three-
loop relation[38], and then run it at three-loop order [39] numerically to the
indicated scale $\mu_{t}$.

###  2.1 Virtual Corrections

Concerning the two-loop form factors $\mathcal{A}_{i}^{(1)}$ in Eq. (15), the
contributions from $\mathcal{A}_{i}^{(1,\triangle)}$ and
$\mathcal{A}_{i}^{(1,\bowtie)}$, as well as the LO contribution, were
evaluated in exact top-mass dependence using the results available from Ref.
[35]. Instead, the two-loop box integrals associated to
$\mathcal{A}_{i}^{(1,\square)}$ were evaluated combining two analytic
approximations corresponding to different kinematical regimes, and here we
briefly recall the main features of this approach, described in detail in Ref.
[36].

The box integrals depend on five scales, namely $m_{Z}$, $m_{H}$, $m_{t}$ and
the kinematic variables $\hat{s}$ and $\hat{t}$, where the latter can be
traded for the transverse momentum of the final-state particles, $p_{T}$. In
the $p_{T}$ expansion of Refs. [40, 35] it is assumed that the scales
associated to $m_{Z},m_{H}$ and to $p_{T}$ are small compared to the scales
set by $\hat{s}$ and $m_{t}$. Under this assumption, the box integrals are
expanded in ratios of small over large scales, and the resulting simplified
integrals are written as linear combinations of 52 master integrals (MI) using
Integration-by-Parts (IBP) identities obtained with LiteRed [41, 42]. On the
other hand, in the high-energy expansion used in Ref. [30] the two-loop box
integrals are first expanded in terms of small $m_{Z}$ and $m_{H}$, then an
IBP reduction is performed on the expanded integrals; the resulting MIs are
further expanded in the limit $m_{t}^{2}\ll\hat{s},|\hat{t}|$ and expressed in
terms of harmonic polylogarithms.

In Ref. [36] it has been shown that, in the forward kinematic regime666Here
and in Ref. [36] it is assumed that the form factors are (anti)symmetric under
$\hat{t}\leftrightarrow\hat{u}$. As a consequence, the discussion above can be
adapted to the backward regime simply replacing $|\hat{t}|$ with $|\hat{u}|$.
defined by $|\hat{t}|\leq|\hat{u}|$, the results of the $p_{T}$ expansion are
accurate in the phase-space region $|\hat{t}|\lesssim 4m_{t}^{2}$, while the
HE expansion is accurate in the complementary region $|\hat{t}|\gtrsim
4m_{t}^{2}$. If the convergence of both expansions is improved using Padé
approximants, then the combination of the results allows to approximate the
exact cross section with an accuracy of 1% or below over the whole phase
space. This procedure has been implemented in a FORTRAN program, which uses
the Padé-improved HE expansion for phase-space points such that
$|\hat{t}|>4m_{t}^{2}$ and $|\hat{u}|>4m_{t}^{2}$, and the Padé-improved
$p_{T}$ expansion for all the remaining phase-space points. When we compute
the results using the $\overline{\mathrm{MS}}$ scheme for the top-mass
renormalization, the form factors $\mathcal{A}_{i}^{(1,\triangle)}$ and
$\mathcal{A}_{i}^{(1,\square)}$ need to be shifted by a quantity that is
defined in Eqs.(14,15) in Ref. [36]. This shift is applied to the form factors
before the construction of the relative Padé approximants.

We use handyG [43] for the evaluation of polylogarithms and the routine of
Ref. [44] for the evaluation of the two elliptic integrals occurring in the
$p_{T}$ expansion results. As a result, the average timing for evaluating one
phase-space point using the Padé-improved HE expansion is around 0.004 s,
while using the Padé-improved $p_{T}$ expansion the timing ranges from 0.02 s
to 0.09 s. For comparison, the evaluation of a phase-space point for the
virtual terms using the results of a small-mass expansion requires on average
2 s [28].

###  2.2 Real Corrections

(a)

(b)

(c)

(d)

Figure 1: Representative Feynman diagrams for the $gg\to ZHg$ process.

For the $gg\to ZHg$ process, we adopt Recola2 [45, 46] to compute the one-loop
matrix element, and we cross checked the result with MadGraph5_aMC@NLO [47].
We include all diagrams with massive and massless closed quark loops. Some
representative diagrams are shown in Fig. 1. We note that in Recola2, the
value of the top-quark mass cannot be changed after process initialization,
and hence the process needs to be reinitialized each time if a dynamical top
mass is adopted (which is the case if a dynamical scale is chosen in the
$\overline{\mathrm{MS}}$ top-mass renormalization scheme, see the next
section). As a result, the average time to compute one phase-space point
increases from 0.2 s to 1.0 s.

For $qg\to ZHq$, and $q\bar{q}\to ZHg$, the one-loop matrix elements are
computed by MadGraph5_aMC@NLO, where we implement a filter to exclude diagrams
without a closed fermion loop. In other words, we include two classes of
Feynman diagrams: in the first class, examples of which are shown in Figs.
2(a) 2(c), both the $Z$ boson and Higgs boson are attached directly or
indirectly (i.e. by connecting to an intermediate virtual boson, similarly to
Fig. 1(c)) to a closed quark loop, while in the second class (as shown in Fig.
2(b) 2(d)) the Higgs boson is attached to a closed quark loop, but the $Z$
boson is radiated from an open fermion line. We note that both types of
diagrams can interfere with tree-level diagrams, hence produce
$\mathcal{O}(\alpha_{s}^{2})$ contributions. Such contributions were studied
in detail777They belong to the classes $R_{I}$ and $R_{II}$ for the _top-
mediated_ terms considered in Ref. [11]. in Ref. [11] and they were considered
as part of the NNLO corrections to $pp\to ZH$. On the other hand, in this
paper we compute the square of those diagrams, corresponding to
$\mathcal{O}(\alpha_{s}^{3})$ contributions that we consider as NLO
corrections to $gg\to ZH$.

(a)

(b)

(c)

(d)

Figure 2: Representative Feynman diagrams for the $qg$ and $q\bar{q}$
channels. In (b) and (d) examples of $Z$-radiated diagrams (see Sec. 3) are
depicted.

##  3 Results

In this section, we present our numerical results for a center-of-mass energy
$\sqrt{s}=13\text{ TeV}$. We adopt the following input parameters:
$m_{t}^{\textrm{OS}}=172.5\,\mathrm{GeV},\;m_{W}=80.385\,\mathrm{GeV},\;m_{Z}=91.1876\,\mathrm{GeV},\;m_{H}=125\,\mathrm{GeV},G_{\mu}=1.1663787\times
10^{-5}\,{\mathrm{GeV}^{-2}}$. We adopt the
NNPDF31_nnlo_as_0118 [48] parton distribution functions in a five flavour
scheme.

###  3.1 Inclusive Cross Section

In Table 1, we show the total cross section at LO and NLO adopting different
top-quark-mass renormalization schemes, i.e. OS and $\overline{\mathrm{MS}}$
with different scale choices. We fix the central value of the renormalization
and factorization scales to be $\mu_{C}=M_{ZH}/2$. The scale uncertainty is
obtained from the envelope of a 7-point variation of the central scale
according to
$(\mu_{R}/\mu_{C},\mu_{F}/\mu_{C})=(1,1),(1,\frac{1}{2}),(1,2),(\frac{1}{2},\frac{1}{2}),(\frac{1}{2},1),(2,1),(2,2)$.

Top-mass scheme | LO [fb] | $\sigma_{LO}/\sigma_{LO}^{OS}$ | NLO [fb] | $\sigma_{NLO}/\sigma_{NLO}^{OS}$ | $K=\sigma_{NLO}/\sigma_{LO}$
---|---|---|---|---|---
On-Shell | $64.01_{-20.3\%}^{+27.2\%}$ | - | $118.6_{-14.1\%}^{+16.7\%}$ | - | 1.85
$\overline{\mathrm{MS}},\mu_{t}=M_{ZH}/4$ | $59.40_{-20.2\%}^{+27.1\%}$ | 0.928 | $113.3_{-14.5\%}^{+17.4\%}$ | 0.955 | 1.91
$\overline{\mathrm{MS}},\mu_{t}=m_{t}^{\overline{\mathrm{MS}}}(m_{t}^{\overline{\mathrm{MS}}})$ | $57.95_{-20.1\%}^{+26.9\%}$ | 0.905 | $111.7_{-14.6\%}^{+17.7\%}$ | 0.942 | 1.93
$\overline{\mathrm{MS}},\mu_{t}=M_{ZH}/2$ | $54.22_{-20.0\%}^{+26.8\%}$ | 0.847 | $107.9_{-15.0\%}^{+18.4\%}$ | 0.910 | 1.99
$\overline{\mathrm{MS}},\mu_{t}=M_{ZH}$ | $49.23_{-19.9\%}^{+26.6\%}$ | 0.769 | $103.3_{-15.6\%}^{+19.6\%}$ | 0.871 | 2.10
Table 1: Total cross section at LO and NLO with full top-quark mass dependence
using different top-quark-mass renormalization schemes. The central value of
the renormalization and factorization scales is fixed to be
$\mu_{R}=\mu_{F}=M_{ZH}/2$. Scale uncertainties are taken from a 7-point scale
variation.

We find that the NLO corrections are large for each choice of the top-mass
renormalization scheme, with an approximate $K$-factor,
$K=\sigma_{NLO}/\sigma_{LO}$, of around 2. Moreover, the relative size of the
scale uncertainties is essentially the same regardless of the top-mass
renormalization scheme. We note that going from LO to NLO the relative size of
the scale uncertainties is reduced by a factor of about $2/3$. The OS scheme
leads to the largest value of the total cross section both at LO and NLO,
while in the $\overline{\mathrm{MS}}$ scheme for $\mu_{t}=M_{ZH}$ the smallest
cross section value is obtained. At LO, the difference between these two
schemes amounts to about 23%, while it decreases to 13% at NLO.

We notice that our OS results are about 20% larger at LO and 14% larger at NLO
than those of Ref. [29] (see Table 1 therein). This discrepancy is mainly due
to the different choice for $\mu_{C}$, and it is related only in a minor way
to the additional diagrams included in our calculation and to the different
input parameters adopted. To verify this, we have computed our results
including the same diagrams and adopting the same input parameters as in Ref.
[28] (which is in accordance with Ref. [29]) and we have found an agreement at
the level of the Monte Carlo error. Furthermore, when we consider the relative
importance of the scale uncertainties, we observe very similar results to Ref.
[29].

###  3.2 Differential Distributions

(a)

(b)

Figure 3: Invariant-mass distribution at LO (magenta) and NLO (green) for the
OS scheme (a) and the $\overline{\mathrm{MS}}$ scheme with the choice
$\mu_{t}=M_{ZH}/2$ (b). The scale uncertainties are depicted as shaded bands.
The lower panels show the $K$-factor.

(a)

(b)

Figure 4: Invariant-mass distribution at LO (magenta) and NLO (green) for the
OS scheme (a) and the $\overline{\mathrm{MS}}$ scheme (b) for a wide $M_{ZH}$
range. The NLO results in which the $Z$-radiated diagrams are excluded are
shown in blue. The lower panels show the $K$-factor.

In Fig. 3, we plot the $M_{ZH}$ distribution in both the OS scheme (3(a)) and
the $\overline{\mathrm{MS}}$ scheme with $\mu_{t}=M_{ZH}/2$ (3(b)) in the
region $M_{ZH}\in[200,800]$ GeV. In both schemes, the $K$-factor is about 3 in
the $ZH$ threshold region, then it decreases as $M_{ZH}$ increases. In the
top-pair threshold region ($M_{ZH}\sim 2\,m_{t}$), the OS scheme gives a peak
with $K$-factor slightly above 2, while the $\overline{\mathrm{MS}}$ scheme
shows a small dip followed by a peak instead. Increasing $M_{ZH}$ to about 800
GeV, the $K$-factor in the OS scheme decreases to about 1.5, while it remains
about 2 in $\overline{\mathrm{MS}}$ scheme.

If $M_{ZH}$ is increased to very large values, as shown in Fig. 4, we observe
that the $K$-factor starts to increase rapidly. At $M_{ZH}=2.5\,\mathrm{TeV}$,
the value of the $K$-factor can reach $\sim 6$ in the OS scheme, and $\sim 10$
in the $\overline{\mathrm{MS}}$ scheme. Such behaviour is due to the inclusion
of diagrams where the $Z$ boson is radiated from an open quark line, as in
Figs. 2(b) and 2(d). For comparison, in Fig. 4 we show also the NLO cross
section when these contributions are excluded. Indeed, one finds that in the
latter case the $K$-factor remains rather flat at high $M_{ZH}$ in both
schemes.

(a)

(b)

Figure 5: (a) Various contributions to $ZH$ production at
$\mathcal{O}(\alpha_{s}^{3})$: the real-emission contributions from the $qg$
and $q\bar{q}$ channel that include all the diagrams considered in this paper
are shown as orange and red lines, respectively. The contributions to the same
channels including only the $Z$-radiated diagrams are shown as green ($qg$)
and yellow ($q\bar{q}$) lines. The Drell-Yan-like contribution at NNLO is
shown as a black line. (b) Comparison of the invariant-mass distribution for
$gg\to ZH$ at LO and NLO with the Drell-Yan-like contribution at NNLO. For
$gg\to ZH$ at NLO, the results with (blue) and without (orange) the
contribution from the $Z$-radiated diagrams are shown.

To further assess the contribution of these $Z$-radiated diagrams, in Fig.
5(a) we show various pieces of the $\mathcal{O}(\alpha_{s}^{3})$ corrections
to the $pp\rightarrow ZH+X$ cross section: we compare the isolated
contribution of the $Z$-radiated diagrams for the $qg$ and $q\bar{q}$ initial
states (green and yellow lines, respectively) to the contribution from each of
the two partonic channels including all the relevant diagrams (orange and red
lines for $qg$ and $q\bar{q}$, respectively). We can see that in the $qg$
channel when $M_{ZH}>1\,\mathrm{TeV}$, the dominant contribution comes from
the square of the $Z$-radiated diagrams. We ascribe this feature to the
contribution from logarithmic terms of EW origin, of the form
$\log^{2}(m_{Z}^{2}/M_{ZH}^{2})$, which become large when the typical scale of
the process ($M_{ZH}$) and the EW scale (represented by $m_{Z}$) are very
different, see e.g. Ref. [49]. On the other hand, although the $q\bar{q}$
channel includes diagrams where the $Z$ boson is radiated from the external
quark lines, its size remains negligible with respect to the total cross
section, though also there we observe that the $Z$-radiated diagrams are
dominating the respective initial state at high $M_{ZH}$. This suppression can
be mainly attributed to the reduced partonic luminosity with respect to the
$qg$ channel. For comparison, in Fig. 5(a) we also report the size of the
Drell-Yan type contribution at NNLO (black line), which we obtained using
vh@nnlo [25, 18] with MCFM[50, 51, 52]. In the lower panel of Fig. 5(a) we
plot the ratio of the $\mathcal{O}(\alpha_{s}^{3})$ corrections computed by us
with respect to the NNLO Drell-Yan contribution. We can see that despite being
$\mathcal{O}(\alpha_{s}^{3})$, the relative importance of the $Z$-radiated
contribution can reach 2% when $M_{ZH}\sim 2\,\mathrm{TeV}$.

In Fig. 5(b) we compare our results for $gg\to ZH$ at LO (green line) and NLO
(blue line) with the Drell-Yan type contribution (black line). In the upper
panel we show the size of the differential cross section for the various
channels, while in the lower panel the ratio of the gluon-fusion with respect
to the NNLO Drell-Yan contribution is displayed. We can see that the gluon-
fusion contribution peaks around the top-pair threshold, which increases its
relative size over the Drell-Yan contribution by about 25% at LO, and about
45% at NLO. The relative size of the gluon-fusion contribution decreases above
the top-pair threshold as $M_{ZH}$ increases, and at NLO becomes dominated by
the $Z$-radiated terms for very large values of $M_{ZH}$. In particular, at 2
TeV the latter constitute more than half of the gluon-fusion contribution.

###  3.3 Change of Renormalization Scheme

(a)

(b)

Figure 6: Invariant-mass distribution at LO (a) and NLO (b) for different
choices of the top-mass renormalization scheme. The lower panels show the
ratio of the $\overline{\mathrm{MS}}$ results using various scales over the OS
results.

To assess the impact of the top-quark-mass renormalization scheme at
differential level, we show in Fig. 6 the differential cross section at LO
(6(a)) and NLO (6(b)) in various top-mass schemes. The lower panel shows the
ratio of the $\overline{\mathrm{MS}}$ results with respect to the OS scheme.

Before discussing the results, we note that the value of the top mass can
affect the shape of the distributions in two ways: first, in the high-energy
regime, it directly controls the overall size of the LO and NLO contributions
via the proportionality of the amplitude to the top-Yukawa coupling888In the
low-energy regime $M_{ZH}\lesssim 2m_{t}$, where an expansion in the
large-$m_{t}$ limit is accurate, the amplitude goes to a constant, as observed
for the first time in Ref. [53].; second, in the low $M_{ZH}$ region the value
of $m_{t}$ shifts the position of the peak associated to the top-pair
threshold, $M_{ZH}=2m_{t}$. The former effect has been already observed by the
authors of Ref. [29], and our results are in agreement when the same scale
choices ($\mu_{t}=M_{ZH}$ and $\mu_{t}=m_{t}^{\overline{\mathrm{MS}}}$) and
the same invariant-mass range ($M_{ZH}>400$ GeV) are considered. On the other
hand, since we can use our analytical results to investigate the low $M_{ZH}$
region, in Fig. 6 we are also able to observe the effect of the peak shift on
the invariant-mass distribution. In particular, we can see that the ratio
(orange line) between $\mu_{t}=M_{ZH}$ and the OS scheme is about 2 in the bin
$M_{ZH}\in[300,325]\,\mathrm{GeV}$ at LO, and decreases to about 1.75 at NLO.
Instead, for $\mu_{t}=M_{ZH}/4$ (green line) this effect is rather small.
Indeed, since renormalization-group evolution predicts that
$m_{t}^{\overline{\mathrm{MS}}}(\mu_{t})$ decreases monotonically as $\mu_{t}$
increases, we expect larger deviations from the OS results when larger values
of $\mu_{t}$ are chosen, see also Ref. [54]. In the region
$M_{ZH}>350\,\mathrm{GeV}$, the OS scheme provides the largest cross section,
while the cross section for the choice $\mu_{t}=M_{ZH}$ is the smallest.
Across all regions, we can see that going from LO to NLO reduces the
difference among different top-mass schemes.

To quantify the size of top-mass renormalization uncertainties, we follow the
procedure adopted in Refs. [55, 54]. In this approach, after a binning of the
$M_{ZH}$ range is chosen, the OS scheme result is taken as the central value,
while for each bin the uncertainty is obtained using the minimal and maximal
values from a set of results in different renormalization schemes. In
particular, for the top mass we considered the set
$\left\\{m_{t}^{\text{OS}};~{}m_{t}^{\overline{\mathrm{MS}}}(m_{t}^{\overline{\mathrm{MS}}});~{}m_{t}^{\overline{\mathrm{MS}}}(M_{ZH}/4);~{}m_{t}^{\overline{\mathrm{MS}}}(M_{ZH}/2);~{}m_{t}^{\overline{\mathrm{MS}}}(M_{ZH})\right\\}$.
In Table 2, we show our results under different choices of the bin size. Due
to the way in which the overall uncertainty is constructed, the latter becomes
of course bigger the smaller the bin size is chosen. This is though mainly due
to the top-pair threshold region, where the location of the bin can lead to a
bigger or smaller top-mass renormalization uncertainty in dependence on how
well the structure at the top-pair threshold999For instance, in Ref. [36] a
peak dip structure for the virtual corrections was observed at the top-pair
threshold. is resolved by the binning, while for large invariant masses the OS
scheme is always the largest and the scheme using
$m_{t}^{\overline{\mathrm{MS}}}(M_{ZH})$ is the smallest. In the very-low
invariant-mass region (i.e. $M_{ZH}\lesssim 275$ GeV) the scheme dependence is
small, as this region can be well described by a large mass expansion. We
hence expect that if the description of the top-pair threshold region was
improved by a Coulomb resummation, the dependence on the binning of the
uncertainty would become smaller.

We observe that the uncertainty nearly halves when going from LO to NLO for
choosing a bin size $\geq 100\,\mathrm{GeV}$. This is similar to what has been
observed for $gg\to HH$ in Ref. [55]. When smaller bin sizes are adopted, the
uncertainty does not half but still shows a sizable reduction.

Bin Width [GeV] | LO | NLO
---|---|---
1 | $64.01_{-35.9\%}^{+15.6\%}$ | $118.6_{-27.0\%}^{+17.2\%}$
5 | $64.01_{-35.6\%}^{+15.3\%}$ | $118.6_{-24.9\%}^{+14.7\%}$
25 | $64.01_{-33.1\%}^{+14.0\%}$ | $118.6_{-20.8\%}^{+10.9\%}$
100 | $64.01_{-25.3\%}^{+2.0\%}$ | $118.6_{-13.7\%}^{+0.6\%}$
$\infty$ | $64.01_{-23.1\%}^{+0\%}$ | $118.6_{-12.9\%}^{+0\%}$
Table 2: Inclusive OS results for the LO and NLO $gg\to ZH$ cross sections and
relative top-mass-scheme uncertainties. The uncertainties are computed
according to Refs. [55, 54].

##  4 Conclusions

In this paper we have presented the evaluation of the QCD corrections at NLO
for the $gg$-initiated channel to $ZH$ associated production. The computation
of the virtual corrections are based on the approach of Ref. [36], while for
this work we have computed the real corrections and included them in a fast
and flexible code. For the inclusive cross section, we have found that the NLO
QCD corrections are of the same size as the LO contribution, and they increase
the cross section by a factor of about 2. Our findings are in accordance with
independent results in the literature [28, 29]. We have also studied the
invariant-mass distribution of the $gg\to ZH$ channel at NLO, observing that
the perturbative $K$-factor is not flat across the $M_{ZH}$ range,
specifically in the region of threshold $ZH$ production and in the very-high-
energy tail ($M_{ZH}>1$ TeV). Furthermore, we have shown that in this latter
region the contribution from real-emission diagrams in which the $Z$ boson is
emitted from an open quark line is the dominant one, and it causes the
$K$-factor to rise up to 10 in the $\overline{\mathrm{MS}}$ scheme. We expect
that in $WH$ production an analogous class of $W$-radiated diagrams will give
a similar contribution in the high-energy regime. While the focus of our
analysis is on LHC phenomenology, this feature could be of interest for
studies at future colliders.

The $gg$-initiated channel is responsible for the larger theoretical
uncertainties in the prediction for $ZH$ production, compared to its $WH$
counterpart. We have shown that the inclusion of the NLO corrections bring
only a mild reduction of the scale uncertainties, about a factor 2/3,
suggesting that more accurate calculations are necessary in order to describe
the $gg$-initiated channel at a level that is adequate for experimental
studies in the future. Finally, the implementation used for our results
allowed us to study for the first time the impact of the uncertainty due to
the renormalization scheme for the top quark mass over the whole invariant-
mass range: we have found that different choices for the top mass scheme can
lead to substantially different results, and we suggest that this uncertainty
should be included in refined theoretical predictions.

## Acknowledgements

We thank Luigi Bellafronte and Pier Paolo Giardino for useful discussions. The
work of M.V. and X.Z. is supported by the Italian Ministry of Research (MUR)
under grant PRIN 20172LNEEZ.

## References

  * [1] ATLAS Collaboration, G. Aad et al., “Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC,” Phys. Lett. B 716 (2012) 1–29, arXiv:1207.7214 [hep-ex].
  * [2] CMS Collaboration, S. Chatrchyan et al., “Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC,” Phys. Lett. B 716 (2012) 30–61, arXiv:1207.7235 [hep-ex].
  * [3] G. Heinrich, “Collider Physics at the Precision Frontier,” Phys. Rept. 922 (2021) 1–69, arXiv:2009.00516 [hep-ph].
  * [4] ATLAS Collaboration, M. Aaboud et al., “Observation of $H\rightarrow b\bar{b}$ decays and $VH$ production with the ATLAS detector,” Phys. Lett. B 786 (2018) 59–86, arXiv:1808.08238 [hep-ex].
  * [5] CMS Collaboration, A. M. Sirunyan et al., “Observation of Higgs boson decay to bottom quarks,” Phys. Rev. Lett. 121 no. 12, (2018) 121801, arXiv:1808.08242 [hep-ex].
  * [6] ATLAS Collaboration, G. Aad et al., “Measurements of $WH$ and $ZH$ production in the $H\rightarrow b\bar{b}$ decay channel in $pp$ collisions at 13 TeV with the ATLAS detector,” Eur. Phys. J. C 81 no. 2, (2021) 178, arXiv:2007.02873 [hep-ex].
  * [7] ATLAS Collaboration, G. Aad et al., “Direct constraint on the Higgs-charm coupling from a search for Higgs boson decays into charm quarks with the ATLAS detector,” arXiv:2201.11428 [hep-ex].
  * [8] LHC Higgs Cross Section Working Group Collaboration, D. de Florian et al., “Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector,” arXiv:1610.07922 [hep-ph].
  * [9] T. Han and S. Willenbrock, “QCD correction to the p p $\to$ W H and Z H total cross-sections,” Phys. Lett. B 273 (1991) 167–172.
  * [10] O. Brein, A. Djouadi, and R. Harlander, “NNLO QCD corrections to the Higgs-strahlung processes at hadron colliders,” Phys. Lett. B 579 (2004) 149–156, arXiv:hep-ph/0307206.
  * [11] O. Brein, R. Harlander, M. Wiesemann, and T. Zirke, “Top-Quark Mediated Effects in Hadronic Higgs-Strahlung,” Eur. Phys. J. C 72 (2012) 1868, arXiv:1111.0761 [hep-ph].
  * [12] M. L. Ciccolini, S. Dittmaier, and M. Kramer, “Electroweak radiative corrections to associated WH and ZH production at hadron colliders,” Phys. Rev. D 68 (2003) 073003, arXiv:hep-ph/0306234.
  * [13] A. Denner, S. Dittmaier, S. Kallweit, and A. Muck, “Electroweak corrections to Higgs-strahlung off W/Z bosons at the Tevatron and the LHC with HAWK,” JHEP 03 (2012) 075, arXiv:1112.5142 [hep-ph].
  * [14] B. A. Kniehl, “Associated Production of Higgs and Z Bosons From Gluon Fusion in Hadron Collisions,” Phys. Rev. D 42 (1990) 2253–2258.
  * [15] D. A. Dicus and C. Kao, “Higgs Boson - $Z^{0}$ Production From Gluon Fusion,” Phys. Rev. D 38 (1988) 1008. [Erratum: Phys.Rev.D 42, 2412 (1990)].
  * [16] G. Ferrera, M. Grazzini, and F. Tramontano, “Associated ZH production at hadron colliders: the fully differential NNLO QCD calculation,” Phys. Lett. B 740 (2015) 51–55, arXiv:1407.4747 [hep-ph].
  * [17] B. Hespel, F. Maltoni, and E. Vryonidou, “Higgs and Z boson associated production via gluon fusion in the SM and the 2HDM,” JHEP 06 (2015) 065, arXiv:1503.01656 [hep-ph].
  * [18] R. V. Harlander, J. Klappert, S. Liebler, and L. Simon, “vh@nnlo-v2: New physics in Higgs Strahlung,” JHEP 05 (2018) 089, arXiv:1802.04817 [hep-ph].
  * [19] W. Bizoń, F. Caola, K. Melnikov, and R. Röntsch, “Anomalous couplings in associated VH production with Higgs boson decay to massive b quarks at NNLO in QCD,” Phys. Rev. D 105 no. 1, (2022) 014023, arXiv:2106.06328 [hep-ph].
  * [20] C. Englert, M. McCullough, and M. Spannowsky, “Gluon-initiated associated production boosts Higgs physics,” Phys. Rev. D 89 no. 1, (2014) 013013, arXiv:1310.4828 [hep-ph].
  * [21] C. Englert, R. Rosenfeld, M. Spannowsky, and A. Tonero, “New physics and signal-background interference in associated $pp\to HZ$ production,” EPL 114 no. 3, (2016) 31001, arXiv:1603.05304 [hep-ph].
  * [22] B. Yan and C. P. Yuan, “Anomalous $Zb\bar{b}$ Couplings: From LEP to LHC,” Phys. Rev. Lett. 127 no. 5, (2021) 051801, arXiv:2101.06261 [hep-ph].
  * [23] B. A. Kniehl and C. P. Palisoc, “Associated production of Z and neutral Higgs bosons at the CERN Large Hadron Collider,” Phys. Rev. D 85 (2012) 075027, arXiv:1112.1575 [hep-ph].
  * [24] L. Altenkamp, S. Dittmaier, R. V. Harlander, H. Rzehak, and T. J. E. Zirke, “Gluon-induced Higgs-strahlung at next-to-leading order QCD,” JHEP 02 (2013) 078, arXiv:1211.5015 [hep-ph].
  * [25] O. Brein, R. V. Harlander, and T. J. E. Zirke, “vh@nnlo - Higgs Strahlung at hadron colliders,” Comput. Phys. Commun. 184 (2013) 998–1003, arXiv:1210.5347 [hep-ph].
  * [26] A. Hasselhuhn, T. Luthe, and M. Steinhauser, “On top quark mass effects to $gg\to ZH$ at NLO,” JHEP 01 (2017) 073, arXiv:1611.05881 [hep-ph].
  * [27] 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].
  * [28] G. Wang, X. Xu, Y. Xu, and L. L. Yang, “Next-to-leading order corrections for $gg\to ZH$ with top quark mass dependence,” arXiv:2107.08206 [hep-ph].
  * [29] L. Chen, J. Davies, G. Heinrich, S. P. Jones, M. Kerner, G. Mishima, J. Schlenk, and M. Steinhauser, “ZH production in gluon fusion at NLO in QCD,” arXiv:2204.05225 [hep-ph].
  * [30] J. Davies, G. Mishima, and M. Steinhauser, “Virtual corrections to $gg\to ZH$ in the high-energy and large-$m_{t}$ limits,” JHEP 03 (2021) 034, arXiv:2011.12314 [hep-ph].
  * [31] L. Chen, G. Heinrich, S. P. Jones, M. Kerner, J. Klappert, and J. Schlenk, “$ZH$ production in gluon fusion: two-loop amplitudes with full top quark mass dependence,” JHEP 03 (2021) 125, arXiv:2011.12325 [hep-ph].
  * [32] S. Alioli, A. Broggio, S. Kallweit, M. A. Lim, and L. Rottoli, “Higgsstrahlung at NNLL’$+$NNLO matched to parton showers in GENEVA,” Phys. Rev. D 100 no. 9, (2019) 096016, arXiv:1909.02026 [hep-ph].
  * [33] S. Zanoli, M. Chiesa, E. Re, M. Wiesemann, and G. Zanderighi, “Next-to-next-to-leading order event generation for $VH$ production with $H\to b\bar{b}$ decay,” arXiv:2112.04168 [hep-ph].
  * [34] R. V. Harlander, A. Kulesza, V. Theeuwes, and T. Zirke, “Soft gluon resummation for gluon-induced Higgs Strahlung,” JHEP 11 (2014) 082, arXiv:1410.0217 [hep-ph].
  * [35] 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].
  * [36] 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,” arXiv:2202.12157 [hep-ph].
  * [37] S. Catani and M. H. Seymour, “A General algorithm for calculating jet cross-sections in NLO QCD,” Nucl. Phys. B 485 (1997) 291–419, arXiv:hep-ph/9605323. [Erratum: Nucl.Phys.B 510, 503–504 (1998)].
  * [38] K. Melnikov and T. v. Ritbergen, “The Three loop relation between the MS-bar and the pole quark masses,” Phys. Lett. B 482 (2000) 99–108, arXiv:hep-ph/9912391.
  * [39] K. G. Chetyrkin, “Quark mass anomalous dimension to O (alpha-s**4),” Phys. Lett. B 404 (1997) 161–165, arXiv:hep-ph/9703278.
  * [40] 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].
  * [41] 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].
  * [42] R. N. Lee, “Presenting LiteRed: a tool for the Loop InTEgrals REDuction,” arXiv:1212.2685 [hep-ph].
  * [43] 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].
  * [44] 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].
  * [45] S. Actis, A. Denner, L. Hofer, A. Scharf, and S. Uccirati, “Recursive generation of one-loop amplitudes in the Standard Model,” JHEP 04 (2013) 037, arXiv:1211.6316 [hep-ph].
  * [46] A. Denner, J.-N. Lang, and S. Uccirati, “Recola2: REcursive Computation of One-Loop Amplitudes 2,” Comput. Phys. Commun. 224 (2018) 346–361, arXiv:1711.07388 [hep-ph].
  * [47] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * [48] NNPDF Collaboration, R. D. Ball et al., “Parton distributions from high-precision collider data,” Eur. Phys. J. C 77 no. 10, (2017) 663, arXiv:1706.00428 [hep-ph].
  * [49] M. Ciafaloni, P. Ciafaloni, and D. Comelli, “Bloch-Nordsieck violating electroweak corrections to inclusive TeV scale hard processes,” Phys. Rev. Lett. 84 (2000) 4810–4813, arXiv:hep-ph/0001142.
  * [50] J. M. Campbell, R. K. Ellis, and C. Williams, “Vector boson pair production at the LHC,” JHEP 07 (2011) 018, arXiv:1105.0020 [hep-ph].
  * [51] J. M. Campbell, R. K. Ellis, and W. T. Giele, “A Multi-Threaded Version of MCFM,” Eur. Phys. J. C 75 no. 6, (2015) 246, arXiv:1503.06182 [physics.comp-ph].
  * [52] R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu, F. Petriello, and C. Williams, “Color singlet production at NNLO in MCFM,” Eur. Phys. J. C 77 no. 1, (2017) 7, arXiv:1605.08011 [hep-ph].
  * [53] B. A. Kniehl, “On the Decay Mode $Z\to Hgg$,” Phys. Rev. D 42 (1990) 3100–3106.
  * [54] J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca, M. Spira, and J. Streicher, “Higgs-Pair Production via Gluon Fusion at Hadron Colliders: NLO QCD Corrections,” JHEP 04 (2020) 181, arXiv:2003.03227 [hep-ph].
  * [55] 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].