aainstitutetext: Institut für Theoretische Teilchenphysik und Kosmologie, RWTH
Aachen University,
52056 Aachen, Germanybbinstitutetext: Max-Planck-Institut für Physik,
Boltzmannstraße 8, 85748 Garching, Germany.ccinstitutetext: The Henryk
Niewodniczański Institute of Nuclear Physics,
ul. Radzikowskiego 152, 31-342 Krakow, Poland

# Quark Mass Effects in Higgs Production

Michał Czakon a Felix Eschment b Marco Niggetiedt c Rene Poncelet a Tom
Schellenberger<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>tom.schellenberger@rwth-
aachen.de

###### Abstract

We examine the effect of finite top- and bottom-quark masses on the Higgs
production cross section in the gluon-gluon fusion channel. We employ both
$\overline{\text{MS}}$ and on-shell renormalisation for the quark masses and
provide a thorough comparison. Furthermore, we explore alternative treatments
of quark masses, in particular in the four-flavour scheme, and investigate
their impact on the cross section. Our work also presents novel predictions
for differential cross sections in the Higgs rapidity. The results lead to a
significant reduction of scale uncertainties, and our analysis enables us to
offer well-grounded recommendations for future research in this area.

††preprint: P3H-24-048, TTK-24-27, IFJPAN-IV-2024-9, MPP-2024-137

##  1 Introduction

The groundbreaking discovery of the Higgs boson in 2012 ATLAS:2012yve ;
CMS:2012qbp  led to an ongoing effort to scrutinise whether this newfound
particle behaves as predicted by the Standard Model (SM). At the Large Hadron
Collider (LHC), the Higgs boson is primarily produced in gluon-gluon fusion,
where two gluons produce the Higgs via a quark-triangle loop. The gluon-gluon-
fusion cross section is therefore one of the central observables in Higgs-
physics phenomenology. Currently, its most precise experimental determination
CMS:2022dwd  is

$\sigma_{ggH}=47.1\pm 3.8\ \mathrm{pb}$ (1)

at a center of mass energy of $13\ \mathrm{TeV}$, and theoretical predictions
with uncertainties negligible in comparison are desirable.

Since the top quark is by far the heaviest quark in the SM, it also has the
largest impact on the cross section, hence Higgs couplings to lighter quarks
are often neglected. Moreover, the top quark’s large mass enables the use of
the heavy-top limit (HTL) approximation, where the top quark is treated as
infinitely massive. This method is valid up to power corrections of order
$m_{H}^{2}/4m_{t}^{2}$, and significantly streamlines computations by reducing
the loop order by one and removing the top-quark mass as a scale. These
simplifications have allowed for calculations at next-to-leading Dawson:1990zj
; Djouadi:1991tka , next-to-next-to-leading Ravindran:2002dc ;
Harlander:2002wh ; Anastasiou:2002yz , and next-to-next-to-next-to-leading
Anastasiou:2015vya ; Mistlberger:2018etf ; Anastasiou:2016cez  order in
$\alpha_{s}$ (NLO, NNLO, $\text{N}^{3}$LO), yielding state-of-the-art
predictions with a theoretical uncertainty of approximately $6\%$. Studies
have shown that the gluon-gluon-fusion channel receives particularly large
perturbative corrections. In fact, the cross section is more than doubled
going from LO to NLO, and NNLO corrections remain substantial, contributing
roughly 20% to the total cross section. Because of the large size of the
corrections, accounting for the previously omitted finite top-quark mass is
strongly advisable also at higher perturbative orders. In Higgs effective
field theory (HEFT), finite-top-quark-mass effects are included by rescaling
the HTL cross section to recover the LO result with finite top-quark mass:

$\sigma^{\text{N}^{n}\text{LO}}_{\text{HEFT}}=\frac{\sigma^{\text{LO}}}{\sigma^{\text{LO}}_{\text{HTL}}}\sigma_{\text{HTL}}^{\text{N}^{n}\text{LO}}\approx
1.065\times\sigma^{\text{N}^{n}\text{LO}}_{\text{HTL}}.$ (2)

This approach assumes that the ratio $\sigma/\sigma_{\text{HTL}}$ remains
relatively constant across perturbative orders, providing a significant
improvement over HTL by capturing much of the finite-top-quark-mass effects
for fully inclusive observables. Nevertheless, with the theoretical precision
of the N3LO predictions, this simple rescaling method is insufficient, and a
more comprehensive treatment of the top-quark mass is indispensable.

Besides to the top quark, the Higgs boson can also couple to other quark
flavours in the loop. Gluon-gluon-fusion amplitudes in which the Higgs couples
to light-quark flavours are suppressed by powers of $m_{q}^{2}/m_{H}^{2}$. The
leading contribution of these amplitudes to the cross section arises via
interference with top-quark induced gluon-gluon fusion. Despite the strong
power suppression, the interference receives large logarithmic enhancements of
the form $\log^{2}(m_{q}^{2}/m_{H}^{2})$, rendering the total effect of
lighter quark flavours on the cross section quite sizeable.

The effect of finite top and bottom masses on the total Higgs production cross
section has been studied up to NNLO in refs. Czakon:2021yub ; Czakon:2023kqm ;
Niggetiedt:2024nmp , thereby reducing the previous theory uncertainty of 6%
Anastasiou:2016cez  by roughly a third. In the present study we aim to expand
on this work and provide a more in-depth analysis. In particular, we want to
address the questions left open in ref. Czakon:2023kqm  and investigate the
impact of different mass-renormalisation and flavour schemes (FS). The
treatment of heavy quarks is not unambiguous. A theory with $n_{l}$ light
quarks taken as massless is called $n_{l}$-flavour scheme, $n_{l}$FS. In ref.
Czakon:2023kqm , the calculation was performed in the 5FS, ergo the bottom
quark was treated as massless. For a non-vanishing contribution, the bottom
quark must still couple to the Higgs. Hence, it is considered massive inside
closed quark loops that couple to the Higgs. Although this method ensures that
the resulting amplitude is gauge invariant, the retrospective introduction of
mass requires stronger justification. We want to address this issue and
explore the impact on the cross section of a consistent bottom-quark treatment
as a massive particle, i.e. the impact of working in the 4FS.

At NLO, the predictions for the top-bottom-interference contribution differ by
approximately 30% between the OS and $\overline{\text{MS}}$ schemes, raising
concerns about the reliability of cross-section predictions at this order.
Furthermore, in the OS scheme, the NNLO contribution to top-bottom
interference is almost of the same size as the NLO correction but with
opposite sign Czakon:2023kqm . Alternating corrections of similar size
typically hint at poor perturbative convergence. Combined with the fact that
the scale uncertainties increase going from NLO to NNLO, an alternative
computational framework is desirable. An alternative mass-renormalisation
scheme could mitigate the poor convergence and yield more reliable
predictions.

Finite-quark-mass effects can also be leveraged to determine the Yukawa
couplings to the respective quarks. The sensitivity of the measurements can be
significantly enhanced by additionally considering differential cross
sections, e.g. the Higgs-$p_{T}$ Bishara:2016jga ; Bonner:2016sdg  or the
Higgs-rapidity Soreq:2016rae  distribution. This is particularly useful for
lighter quarks like the charm quark. CMS studies CMS:2018gwt ; CMS:2023gjz
have found an upper bound for the charm Yukawa coupling of about five times
its SM value, by fitting NLO predictions to the measured distributions.
Meanwhile, NNLO predictions for $p_{T}$ distributions with finite top and
bottom mass have become available Grazzini:2013mca ; Melnikov:2016emg ;
Lindert:2017pky ; Caola:2018zye ; Bonciani:2022jmb  and could potentially
yield improved bounds. Rapidity distributions on the other hand are even more
challenging, as they require a contribution from the three-loop virtual
corrections absent in the $p_{T}$ distribution. Our study aims to bridge this
gap, potentially enhancing the sensitivity for future analyses of this kind.

This paper is organised as follows: In section 2 we explain our computational
framework. We put a special emphasis on the modifications required for
different mass-renormalisation and flavour schemes. In section 3 we present
the results for the fully-inclusive cross sections as well as differential
distributions for the Higgs-$p_{T}$ and Higgs-rapidity. In section 4 we give
our conclusions and recommendations for future phenomenological studies.
Finally, appendix A contains a comparison of our results to existing Higgs+jet
analyses.

##  2 Computational Details

The general workflow mirrors that of ref. Czakon:2023kqm . In essence, we
calculate the gluon-gluon-fusion amplitudes for both real and virtual
corrections, i.e. the amplitudes for three-loop double-virtual, two-loop real-
virtual and one-loop double-real corrections, and then perform the phase-space
integration using Monte Carlo techniques. Soft and collinear divergences which
cancel between the different-multiplicity phase spaces as well as those
absorbed in the initial state renormalisation are handled with Stripper – the
C++ implementation of the sector-improved residue subtraction scheme
Czakon:2010td ; Czakon:2014oma ; Czakon:2019tmo .

Amplitudes for double-real radiation have been calculated in ref.
Budge:2020oyl  and we use the implementation available in MCFM
Campbell:2019dru ; Campbell:1999ah ; Campbell:2011bn , which in turn uses
QCDLoop Ellis:2007qk ; Carrazza:2016gav  for the evaluation of scalar
integrals. The three-loop double-virtual corrections are known in terms of
large-mass and/or high-energy expansions to very high precision for both the
single-quark Czakon:2020vql  as well as the mixed-flavour contributions
Niggetiedt:2023uyk .

The two-loop real-virtual corrections can be categorised into mixed- and
single-flavour diagrams. The former consist of diagrams with two different
quark flavours (see fig. 3) and they can be written as a product of simple
one-loop integrals which we solve analytically. The single-flavour diagrams on
the other hand involve genuine two-loop integrals, which we solve numerically.
To do this, the amplitude is first projected to form factors and the required
set of master integrals is determined by numerically solving the system of
integration-by-parts identities. All appearing loop integrals are reduced to a
predefined basis of master integrals with the public software
Kira$\oplus$FireFly Maierhofer:2017gsa ; Maierhofer:2018gpa ; Klappert:2020nbg
. The reduction to master integrals is computationally very expensive, mostly
because of the tremendous complexity in intermediate expressions, which is why
we employ finite field methods for the interpolation of multivariate rational
expressions as implemented in FireFly Klappert:2019emp ; Klappert:2020aqs .
Even still, a full reconstruction in terms of all three dimensionless
kinematic variables needed for the $gg\longrightarrow gH$,
$q\bar{q}\longrightarrow gH$, $qg\longrightarrow qH$ and
$\bar{q}g\longrightarrow\bar{q}H$ processes, which for example can be chosen
as

$z=1-\frac{m_{H}^{2}}{\hat{s}},\quad\lambda=\frac{\hat{t}}{\hat{t}+\hat{u}},\quad\text{and}\quad
x=\frac{m_{q}^{2}}{\hat{s}},$ (3)

with $\hat{s},\hat{t}$ and $\hat{u}$ the usual partonic Mandelstam variables,
remains challenging. In order to minimise the computational effort, a
particular basis of master integrals involving dimension shifts for the most
complicated integrals is chosen. Furthermore, we use a fixed ratio of the
quark and Higgs mass square, $m_{q}^{2}/m_{H}^{2}$, to eliminate one of the
scales in favour of a rational number, which further simplifies the rational-
function coefficients in front of the master integrals. In fact, we compute
the master integrals and the complete amplitudes for five different bottom-
quark masses and four different top-quark masses, listed in tab. 1. After the
reduction to master integrals, we derive a set of coupled differential
equations. The integrals for the boundary condition at $x\rightarrow\infty$
can be solved with the help of the large-mass expansion. We then transport the
solution to physical values of $x$, by solving the differential equations
numerically using the Bulirsch-Stoer algorithm (see fig. 1).

Figure 1: Illustration of the strategy for numerically solving the
differential equations of the master integrals. Integrals on the boundary at
$x\rightarrow\infty$ are computed with the large-mass expansion. The
differential equation is then solved in $x$ for multiple $z$ and a fixed value
of $\lambda$ close to the symmetry axis. The red plane indicates the physical
plane, with a fixed quark mass.

This is done multiple times for many values of $z$ and a fixed value of
$\lambda$ close to the symmetry axis of $\lambda=1/2$. Subsequently, we solve
the differential equations in $\lambda$ to map out the entire physical plane,
marked red in fig. 1, such that we obtain the amplitude as a dense two
dimensional grid. During the Monte Carlo integration we interpolate between
the grid points using cubic splines. The accuracy of the splines was verified
by removing a mass grid and recalculating the resulting cross sections. The
discrepancies between the recalculated cross sections and the original ones
were below the Monte Carlo uncertainties, confirming the robustness of our
interpolation method. We note that we only calculated the interference of all
bottom-quark amplitudes with the OS-renormalised, as well as with the
$\overline{\mathrm{MS}}$-renormalised top-quark amplitudes at the same
renormalisation scale, to construct grids for cubic-spline interpolation.

To further improve the quality of the spline fits, we subtract all infrared
divergences from the amplitudes and add them back only after performing the
interpolation. The subtraction of infrared divergences can be simplified by
first subtracting a rescaled cross section of the HTL,

$\begin{split}&\braket{M^{(0)}_{gg\rightarrow Hg}}{M^{(1)}_{gg\rightarrow
Hg}}|_{\text{regulated}}\equiv\braket{M^{(0)}_{gg\rightarrow
Hg}}{M^{(1)}_{gg\rightarrow Hg}}\\\
&\quad-\left[r\braket{M^{(0)}_{gg\rightarrow
Hg,\text{HTL}}}{M^{(1)}_{gg\rightarrow
Hg,\text{HTL}}}+\frac{8\pi\alpha_{s}}{\hat{t}}\braket{P_{gg}^{(0)}\\!\left(\frac{\hat{s}}{\hat{s}+\hat{u}}\right)}\braket{M^{(0)}_{gg\rightarrow
H}}{M^{(1)}_{gg\rightarrow H}-M^{(1)}_{gg\rightarrow H,\text{HTL}}}\right]\\\
&\braket{P_{gg}^{(0)}\\!\left(z\right)}=2C_{A}\left(\frac{z}{1-z}+\frac{1-z}{z}+z(1-z)\right).\end{split}$
(4)

The rescaling factor is chosen in such a way as to match the LO cross section
with finite quark mass

$r=\frac{\sigma^{\text{LO}}_{i}}{\sigma^{\text{LO}}_{\text{HTL}}},$ (5)

where $\sigma_{i}$ refers to the cross section of the interference
contribution or the cross section with finite top-quark mass. This way, the
first subtraction term in eq. (4) already captures most of the infrared
singularities. The only infrared divergence which is not yet removed stems
from a soft or collinear splitting of the two-loop form factors. After
subtracting these left-over divergences with the corresponding splitting
functions, the regulated squared matrix element is free of infrared
singularities apart from logarithmic ones. Note that the splitting function
also removes the soft divergences. The regulated squared amplitude for
$qg\longrightarrow qH$ is defined analogously with the above splitting
function replaced by

$-\braket{P_{qq}^{(0)}(z)}=-T_{F}\left(1-2z(1-z)\right).$ (6)

The amplitude for $q\bar{q}\longrightarrow gH$, on the other hand, has no
additional singularities and no subtraction of splitting functions is needed.

###  2.1 $\overline{\text{MS}}$ scheme

Compared to the calculation in the OS scheme, for the $\overline{\text{MS}}$
scheme we need to introduce two modifications: first, the running of the mass
has to be considered, and second, the renormalisation of the amplitudes has to
be altered. The latter does not require us to redo the entire calculation from
scratch, in fact the $\overline{\text{MS}}$-renormalised amplitudes can be
easily derived from the OS-renomalised amplitudes. To do this we first use the
renormalisation-scheme invariance of the bare quark mass,

$Z_{m}^{\overline{\text{MS}}}\overline{m}=Z_{m}^{\text{OS}}m^{\text{OS}},$ (7)

to derive the relation between the OS- and the
$\overline{\text{MS}}$-renormalised masses:

$m^{\text{OS}}=\overline{m}\left(1+c_{1}\frac{\alpha_{s}}{\pi}+c_{2}\left(\frac{\alpha_{s}}{\pi}\right)^{2}+\mathcal{O}(\alpha^{3})\right).$
(8)

The mass-renormalisation constants were first derived in refs. Gray:1990yh ;
TARRACH1981384 . The renormalisation procedure must also be consistent with
the used FS, that means in the 5FS, the bottom quark has to be massless inside
closed bottom loops even when renormalising the bottom mass. Similarly, we
have to decouple the top quark, and in the 4FS also the bottom quark, from the
running of $\alpha_{s}$:

$\begin{split}&\alpha_{s}^{(n_{l})}(\mu)=\zeta\left(\alpha_{s}^{(n_{l}+1)}(\mu),\ln\frac{m^{2}}{\mu^{2}}\right)\alpha_{s}^{(n_{l}+1)}(\mu),\\\
&\zeta(\alpha_{s},L)=1+\frac{\alpha_{s}}{4\pi}\frac{4}{3}T_{F}L+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\left[-\frac{14}{3}+T_{F}\left(\frac{20}{3}C_{A}+4C_{F}\right)L+\frac{16}{9}T_{F}^{2}L^{2}\right].\end{split}$
(9)

Here $C_{F}=4/3,\ C_{A}=3,\ T_{F}=1/2$ and $n_{l}$ is the number of light
quarks. We then use Eq. (8) recursively to eliminate any residual logarithmic
dependence on the OS mass and finally obtain

$\begin{split}c_{1}&=C_{F}\left(1-\frac{3}{4}\ln\\!\frac{\overline{m}^{2}}{\mu^{2}}\right),\\\
c_{2}&=-\frac{1}{384}C_{F}\bigg{(}C_{A}\left(-132\log^{2}\\!\frac{\overline{m}^{2}}{\mu^{2}}+740\log\\!\frac{\overline{m}^{2}}{\mu^{2}}+144\zeta(3)-1111+\pi^{2}(32-96\log(2))\right)\\\
&\hskip
28.45274pt+3C_{F}\bigg{(}-36\log^{2}\\!\frac{\overline{m}^{2}}{\mu^{2}}-36\log\\!\frac{\overline{m}^{2}}{\mu^{2}}-96\zeta(3)+71+8\pi^{2}(8\log(2)-5)\bigg{)}\\\
&\hskip
28.45274pt+4n_{l}T_{F}\left(12\log^{2}\\!\frac{\overline{m}^{2}}{\mu^{2}}-52\log\\!\frac{\overline{m}^{2}}{\mu^{2}}+8\pi^{2}+71\right)\bigg{)}\\\
&\quad-\sum_{i}\frac{1}{96}C_{F}T_{F}\bigg{(}12\log^{2}\\!\frac{\overline{m}^{2}}{\mu^{2}}-24\log\\!\frac{\overline{m}^{2}}{\mu^{2}}\log\\!\frac{\overline{m}_{i}^{2}}{\mu^{2}}-52\log\\!\frac{\overline{m}^{2}}{\mu^{2}}+32\log\\!\frac{\overline{m}_{i}^{2}}{\mu^{2}}\\\
&\hskip
28.45274pt+48(x_{i}-1)^{2}\left(x_{i}^{2}+x_{i}+1\right)\text{Li}_{2}(1-x_{i})-48\left(x_{i}^{4}+x_{i}^{3}+x_{i}+1\right)\left(\text{Li}_{2}(-x_{i})+\log\\!(x_{i})\log\\!(x_{i}+1)\right)\\\
&\hskip
28.45274pt+48x_{i}^{4}\log^{2}(x_{i})+48x_{i}^{2}\log(x_{i})-16\pi^{2}(x_{i}^{3}+x_{i})+72x_{i}^{2}+71\bigg{)},\end{split}$
(10)

where we defined $x_{i}=\frac{\overline{m}_{i}}{\overline{m}}$. The sum runs
over all heavy-quark flavours which are decoupled from the running of
$\alpha_{s}$. This includes $\overline{m}$ itself in case the corresponding
quark is decoupled. Subsequently, the amplitudes can simply be expanded to
find the relation between the two schemes:

$\begin{split}&\mathcal{M}^{\overline{\text{MS}}}=\mathcal{M}^{\text{OS}}+\delta\mathcal{M},\\\
&\delta\mathcal{M}^{(1)}=\overline{m}c_{1}\frac{\alpha_{s}}{\pi}\frac{\mathrm{d}\mathcal{M}^{\text{OS},(0)}}{\mathrm{d}m}\bigg{|}_{m=\overline{m}},\\\
&\delta\mathcal{M}^{(2)}=\overline{m}\left[c_{1}\frac{\alpha_{s}}{\pi}\frac{\mathrm{d}\mathcal{M}^{\text{OS},(1)}}{\mathrm{d}m}\bigg{|}_{m=\overline{m}}+c_{2}\left(\frac{\alpha_{s}}{\pi}\right)^{2}\frac{\mathrm{d}\mathcal{M}^{\text{OS},(0)}}{\mathrm{d}m}\bigg{|}_{m=\overline{m}}\right]+\frac{1}{2}\left(\overline{m}c_{1}\frac{\alpha_{s}}{\pi}\right)^{2}\frac{\mathrm{d}^{2}\mathcal{M}^{\text{OS},(0)}}{\mathrm{d}m^{2}}\bigg{|}_{m=\overline{m}},\end{split}$
(11)

where the superscript indicates the perturbative order in $\alpha_{s}$.

To compute the running of the quark masses, we use the public software package
CRunDec Schmidt:2012az  at four-loop accuracy.

Table 1: Mass ratios squared, $m_{b}^{2}/m_{H}^{2}$ and $m_{t}^{2}/m_{H}^{2}$, used for the generation of numerical values of the two-loop four-point amplitudes. We chose values corresponding to the OS masses, according to PDG recommendations Workman:2022ynf , and $\overline{\text{MS}}$-masses at scales relevant for 7-point variation w.r.t. to a central scale of $m_{H}/2$. For the latter, we used $\overline{m}_{b}(\overline{m}_{b})=4.18\ \text{GeV}$, as recommended by the PDG and $\overline{m}_{t}(\overline{m}_{t})=162.7\ \mathrm{GeV}$. The dynamic scale (12) used for our differential cross-section predictions can exceed $m_{H}$, and reaches values of up to $1$ TeV. At this high scale, the bottom-quark mass is $\overline{m}_{b}(1\ \mathrm{TeV})=2.41~{}\mathrm{GeV}$. To avoid extrapolations outside our numerical grids for the amplitudes, we therefore included another mass value, $\overline{m}_{b}^{\text{min}}$, below this minimum. Grid mass | Value [GeV] | Approximate ratio $m_{b,t}^{2}/m_{H}^{2}$ | Relative error [‰]
---|---|---|---
$m_{b}$ | 4.78 | $\frac{1}{684}$ | 0.2
$\overline{m}_{b}(m_{H})$ | 2.789 | $\frac{1}{2011}$ | 0.9
$\overline{m}_{b}(m_{H}/2)$ | 2.961 | $\frac{1}{1782}$ | 0.2
$\overline{m}_{b}(m_{H}/4)$ | 3.170 | $\frac{1}{1557}$ | 1.0
$\overline{m}_{b}^{\text{min}}$ | 1.67 | $\frac{1}{5602}$ | 0.1
$m_{t}$ | $172.4$ | $\frac{23}{12}$ | 4
$\overline{m}_{t}(m_{H})$ | $166.1$ | $\frac{136}{77}$ | 0.0
$\overline{m}_{t}(m_{H}/2)$ | $176.2$ | $\frac{149}{75}$ | 0.0
$\overline{m}_{t}(m_{H}/4)$ | $188.2$ | $\frac{213}{94}$ | 0.0

###  2.2 4-flavour scheme

In the 5FS, the bottom quark is treated as massless except in closed fermion
loops that couple to the Higgs. This treatment is formally equivalent to
introducing two versions of the bottom quark: one with zero mass and another
with finite mass which is absent from the proton. Since the Yukawa coupling to
the Higgs is an arbitrary parameter in QCD, all infrared poles and gauge
dependence must cancel independently of this parameter. This simple argument
ensures that, even though we only consider a subset of all Feynman diagrams
allowed in the SM, the resulting cross section will be internally consistent.
Also, note that if the bottom-quark mass is discarded while maintaining a
finite Yukawa coupling, the cross section will vanish all together, as the
scalar nature of the coupling requires a helicity flip inside fermion loops.

Figure 2: Example diagram contributing to the three-loop double-virtual
corrections to Higgs production containing two fermion loops of different
flavours. A full list of diagrams can be found in fig. 1 of ref.
Czakon:2020vql . With the thick fermion line corresponding to a top quark and
the thin line corresponding to a bottom quark in the 4-flavour scheme, this
diagram contains quarks of two different masses. Figure 3: Diagrams with two
fermion loops of different flavours contributing to the real-virtual
corrections to Higgs production. The set is complete except for permutations
of external gluons and reversal of the thick fermion-loop-line direction.
Initial-state quarks are massless. Quark-gluon initial states can be obtained
from the quark-antiquark case by crossing. Relevance to the 4-flavour scheme
as in fig. 2. Figure 4: Diagrams corresponding to the real radiation of a
quark-antiquark pair in Higgs production. Permutation of external particles
and reversal of fermion-loop direction contribute additional diagrams. In the
4-flavour scheme, radiation of massive bottom quarks has to be taken into
account. Infrared divergences are regulated by the finite mass and cancel with
the contributions from figs. 2 and 3. The shaded region in the first diagram
might be shifted to the PDFs in the 5-flavour scheme, whereas bottom quarks
are absent from the initial state in the 4-flavour scheme.

Although the above approach is consistent in the framework of this
calculation, the inhomogeneous treatment of light quark masses depending on
the loop structure is rather ad hoc. It is therefore desirable to have an
alternative computational framework without these contentions. For that reason
we also perform the calculation in the 4FS, in which the bottom quark is
always treated as a massive particle, and hence excluded from the proton PDFs.
Then, amplitudes like those depicted in figs. 2 or 3 exhibit logarithmic mass
enhancements, which cancel once the radiation of a massive quark pair (see
fig. 4) is taken into account as well. These one-loop amplitudes are computed
with Recola Actis:2016mpe ; Denner:2017wsf . Note that we still exclude
diagrams in which external bottom quarks couple to the Higgs, corresponding to
the tree-level amplitude $gg\longrightarrow b\bar{b}H$, as these diagrams as
well as all the interference contributions are accounted for in corrections to
the bottom-bottom fusion production channel Dittmaier:2003ej  in the 4FS (see
fig. 5).

Figure 5: Example interference diagram for a contribution to the bottom-bottom fusion channel at NLO in the 4FS. The contribution is proportional to $Y_{b}Y_{t}$. Because the Higgs couples to an external bottom quark, this kind of diagrams is not considered in the present calculation as a correction to gluon-gluon fusion. Table 2: HEFT cross section in the 5-flavour scheme and for different bottom-quark masses in the 4-flavour scheme. In the last column the cross section and the scale variation are computed with the $\overline{\text{MS}}$-mass. The results are computed for LHC @ 13 TeV using the NNPDF31_nnlo_as_0118 PDF set. The central scale is chosen at $\mu_{R}=\mu_{F}=m_{H}/2$. The scale uncertainties are determined with seven-point variation. Results for the 5FS are computed with SusHi Harlander:2012pb ; Harlander:2016hcx . Order |  $\sigma_{\text{HEFT}}$ [pb]
---|---
$\sqrt{s}=13$ TeV
| 5FS | 4FS | 4FS | 4FS | 4FS
|  |  $m_{b}=0.01$ GeV |  $m_{b}=0.1$ GeV |  $m_{b}=4.78$ GeV |  $\overline{m}_{b}(\overline{m}_{b})=4.18$ GeV
$\mathcal{O}(\alpha_{s}^{2})$ | $+16.30$ | +16.27 | +16.27 | +16.27 | $16.27$
LO | $16.30^{+4.36}_{-3.10}$ | $16.27^{+4.63}_{-3.22}$ | $16.27^{+4.63}_{-3.22}$ | $16.27^{+4.63}_{-3.22}$ | $16.27^{+4.63}_{-3.22}$
$\mathcal{O}(\alpha_{s}^{3})$ | +21.14 | +20.08(3) | +20.08(3) | +20.08(3) | +20.08(3)
NLO | $37.44^{+8.42}_{-6.29}$ | $36.35(3)^{+8.57}_{-6.32}$ | $36.35(3)^{+8.57}_{-6.32}$ | $36.35(3)^{+8.57}_{-6.32}$ | $36.35(3)^{+8.57}_{-6.32}$
$\mathcal{O}(\alpha_{s}^{4})$ | +9.72 | +10.8(4) | +11.1(4) | +9.5(2) | $+9.6(2)$
NNLO | $47.16^{+4.21}_{-4.77}$ | $47.2(4)^{+5.4}_{-5.4}$ | $47.5(4)^{+5.4}_{-5.5}$ | $45.9(2)^{+4.3}_{-4.9}$ | $46.0(2)^{+4.4}_{-5.0}$

To validate that the cross section is free from any enhancements due to the
bottom quark mass, we computed the Higgs production cross section in top-
induced gluon-gluon fusion within the HEFT framework for various bottom-quark
masses111We note that the rescaling factor used to define HEFT does not depend
on the flavour scheme in use. As indicated in eq. (2), only the effect of a
finite top-quark mass is taken into account in the rescaling.. The results are
presented in tab. 2. We found that the 4FS results for extremely low bottom
quark masses are in very good agreement with the 5FS result, and show no sign
of divergence. The findings also highlight the sizeable impact of finite-
quark-mass effects, even if the quark does not directly couple to the Higgs,
as the cross section is shifted by more than 2% for realistic bottom-quark
masses (last two columns of tab. 2). This few percent effect seems to be in
rough agreement with the order of magnitude of this effect in Higgs+jet
production Pietrulewicz:2023dxt .

##  3 Results

In this section we present our results for the finite-quark-mass effects in
the gluon-gluon-fusion channel for various computational setups. Following the
Higgs Working Group recommendations LHCHiggsCrossSectionWorkingGroup:2016ypw ,
our reference computation uses the 5FS, an OS-renormalised top quark with a
pole mass of $m_{t}=173.06$ GeV222The Higgs-Working-Group recommends using an
OS-mass of $m_{t}=172.4$. The different value was chosen because of the
rationalisation procedure described in sec. 2. We have verified, that choosing
a different mass has a negligible impact on the final result. and an
$\overline{\text{MS}}$-renormalised bottom mass of
$\overline{m}_{b}(\overline{m}_{b})=4.18$ GeV. We then explore the impact of
alternative mass-renormalisation schemes for both the top and the bottom
quark. In addition, we perform the computation in the 4FS to examine if the
flavour schemes are compatible.

We use two scale choices: for fully-inclusive cross sections we use a fixed
scale with a central value of $\mu=m_{H}/2$, while differential cross sections
are computed with a dynamic scale and a central value of

$\mu=\frac{H_{T}}{2}\equiv\frac{1}{2}\left(\sqrt{m_{H}^{2}+p_{T}^{2}}+\sum_{i}p_{i,T}\right)$
(12)

where $p_{T},\ p_{i,T}$ are the transverse momenta of the Higgs and the $i$-th
final state parton, respectively. In the $\overline{\mathrm{MS}}$ scheme, the
mass of the bottom quark is also evaluated at $\mu$. The scale was chosen to
match previous studies Lindert:2017pky ; Bonciani:2022jmb ; Jones:2018hbb .

Table 3: Effect of the finite top-quark mass on the gluon-gluon fusion cross section for two different computational setups. The results are computed for LHC @ 13 TeV using the NNPDF31_nnlo_as_0118 PDF set. The central scale is chosen at $\mu_{R}=\mu_{F}=m_{H}/2$. The scale uncertainties are determined with seven-point variation. Order |  $(\sigma_{t}-\sigma_{\text{HEFT}})$ [pb]
---|---
$\sqrt{s}=13$ TeV
| 5FS | 4FS
|  $m_{t}=173.06$ GeV |  $m_{t}=173.06$ GeV
|  |  $\overline{m}_{b}(\overline{m}_{b})=4.18$ GeV
LO | - | -
$\mathcal{O}(\alpha_{s}^{3})$ | $-0.30$ | $-0.27$
NLO | $-0.30^{+0.10}_{-0.17}$ | $-0.27^{+0.09}_{-0.16}$
$\mathcal{O}(\alpha_{s}^{4})$ | $+0.14$ | $+0.12$
NNLO | $-0.16^{+0.13}_{-0.03}$ | $-0.15^{+0.10}_{-0.02}$

Tab. 3 illustrates the impact of a finite top-quark mass on the cross section,
i.e. the difference between the result in the full theory restricted to Higgs
production through a closed top-quark loop, $\sigma_{t}$, and in HEFT. The
difference between the 5FS and 4FS (second and third column of tab. 3) only
amounts to $-0.01\ \mathrm{pb}$ at NNLO. We remind, however, that there is a
significant $-2.5\%$ shift between the 5FS and 4FS results in HEFT itself (see
tab. 2). The influence of the flavour scheme beyond HEFT is tiny, because the
effects of a finite top-quark mass are power-suppressed and at NNLO only
amount to roughly 3‰ of the total cross section.

In tab. 4 the difference of the cross section $\sigma_{t}$ in the two
renormalisation schemes considered for the top-quark mass is presented. We
find that it only amounts to $-0.01\ \mathrm{pb}$ at NNLO.

Table 4: Difference of cross sections for Higgs production through a closed top-quark loop with the top-quark mass defined either in the $\overline{\mathrm{MS}}$ or the OS scheme. The results are computed for LHC @ 13 TeV using the NNPDF31_nnlo_as_0118 PDF set. The central scale is chosen at $\mu_{R}=\mu_{F}=m_{H}/2$. The scale uncertainties are determined with seven-point variation. Order |  $(\sigma_{t}^{\overline{\mathrm{MS}}}-\sigma_{t}^{\mathrm{OS}})$ [pb]
---|---
$\sqrt{s}=13$ TeV
$\mathcal{O}(\alpha_{s}^{2})$ | $-0.04$
LO | $-0.04^{+0.12}_{-0.17}$
$\mathcal{O}(\alpha_{s}^{3})$ | $+0.02$
NLO | $-0.02^{+0.14}_{-0.30}$
$\mathcal{O}(\alpha_{s}^{4})$ | $+0.01$
NNLO | $-0.01^{+0.12}_{-0.24}$

Tab. 5 presents the results for the top-bottom interference contribution to
the cross section $\sigma_{t\times b}$. As already noticed in ref.
Czakon:2023kqm , the interference converges quite badly if the bottom-quark
mass is OS-renormalised, as indicated by the alternating perturbative
corrections of a similar order of magnitude (see third column of tab. 5).
Furthermore, scale uncertainties fail to correctly estimate the effects of
higher orders and actually increase from NLO to NNLO, reaching 15%. However,
results for $\overline{\text{MS}}$-renormalised bottom-quark masses (second
column) show a much better perturbative convergence and additionally, we
observe a decrease in scale uncertainties to below 8%. Adapting a mixed
renormalisation scheme, in which only the Yukawa coupling is renormalised in
the $\overline{\text{MS}}$ scheme, while the mass in the propagators is kept
in the OS scheme, also results in a better convergence and reduced scale
uncertainties as noted in ref. Czakon:2023kqm . Without the inclusion of
electroweak corrections this mixed approach is valid, but if these effects are
taken into account as well, the Higgs vacuum expectation value acquires an
additional finite renormalisation that propagates to the masses of gauge
bosons. Since the masses of the gauge bosons are defined in the on-shell
scheme by default, the mixed approach is inconsistent. Ignoring this issue, we
note that the main improvements to the perturbation series of our problem are
already captured by the mixed scheme. For the top-quark mass, the
renormalisation scheme seems to have very little impact on the final result.
Lastly, the NNLO result in the 4FS is 5.2% lower than in the 5FS and agrees
well within scale uncertainties. The magnitude of this decrease is of a
similar size as observed in HEFT (see tab. 2).

Table 5: Top-bottom interference contribution to the gluon-gluon fusion cross section for various computational setups. The results are computed for LHC @ 13 TeV using the NNPDF31_nnlo_as_0118 PDF set. The central scale is chosen at $\mu_{R}=\mu_{F}=m_{H}/2$. The scale uncertainties are determined with seven-point variation. Order |  $\sigma_{t\times b}$ [pb]
---|---
$\sqrt{s}=13$ TeV
| 5FS | 5FS | 5FS | 4FS
|  $m_{t}=173.06$ GeV |  $m_{t}=173.06$ GeV |  $m_{t}(m_{t})=162.7$ GeV |  $m_{t}=173.06$ GeV
|  $\overline{m}_{b}(\overline{m}_{b})=4.18$ GeV |  $m_{b}=4.78$ GeV |  $\overline{m}_{b}(\overline{m}_{b})=4.18$ GeV |  $\overline{m}_{b}(\overline{m}_{b})=4.18$ GeV
$\mathcal{O}(\alpha_{s}^{2})$ | $-1.11$ | $-1.98$ | $-1.12$ | $-1.15$
LO | $-1.11^{+0.28}_{-0.43}$ | $-1.98^{+0.38}_{-0.53}$ | $-1.12^{+0.28}_{-0.42}$ | $-1.15^{+0.29}_{-0.45}$
$\mathcal{O}(\alpha_{s}^{3})$ | $-0.65$ | $-0.44$ | $-0.64$ | $-0.66$
NLO | $-1.76^{+0.27}_{-0.28}$ | $-2.42^{+0.19}_{-0.12}$ | $-1.76^{+0.27}_{-0.28}$ | $-1.81^{+0.28}_{-0.30}$
$\mathcal{O}(\alpha_{s}^{4})$ | $+0.02$ | $+0.43$ | $-0.02$ | $-0.02$
NNLO | $-1.74(2)^{+0.13}_{-0.03}$ | $-1.99(2)^{+0.29}_{-0.15}$ | $-1.78(1)^{+0.15}_{-0.03}$ | $-1.83(2)^{+0.14}_{-0.03}$

Similar to the total cross section, the Higgs-rapidity distribution due to the
top-bottom interference (right panel of fig. 6) also shows that the
$\overline{\text{MS}}$ renormalisation generally yields smaller scale
uncertainties. The different schemes are compatible within the estimated error
bands. The Higgs-$p_{T}$ distribution (left panel of fig. 6) reveals that the
main improvement of the uncertainties comes from the low-$p_{T}$ region,
specifically the first bin below 10 GeV. For higher $p_{T}$, the scale
uncertainties are of very similar size, and even slightly smaller for the OS
scheme. At low $p_{T}$, the effect of finite bottom-quark masses is highly
relevant, reaching almost 8% of the total cross section, whereas the effect is
almost negligible (below 1%) above 50 GeV. The rapidity distribution shows
less-pronounced features, closely resembling a constant shift of about $-4\%$
across all rapidity bins. The $p_{T}$ distributions were compared to ref.
Caola:2018zye  for the OS-renormalised cross section (see fig. 8). We find
good agreement except for very low $p_{T}$, where the central value is
identical but our scale uncertainties are slightly smaller. The
$\overline{\text{MS}}$-renormalised $p_{T}$ distribution was also successfully
checked against ref. Bonciani:2022jmb  (see fig. 9).

Figure 6: Interference contribution to the Higgs-$p_{T}$ (left) and Higgs-
rapidity (right) distribution. Compared are the results for an
$\overline{\text{MS}}$\- and OS-renormalised bottom-quark mass. The lower
panel shows the ratio to the distribution including bottom- and top-quark
effects at NNLO. Coloured bands indicate the envelope of the seven-point scale
variation and the error bars represent the Monte Carlo uncertainties.

In figs. 7, we compare distributions in HEFT against the cross sections in
full QCD333Contributions without top quark couplings are excluded, since they
are negligible. We also neglect couplings to quarks lighter than the bottom
quark.. As is well-known, the $p_{T}$ distribution is very sensitive to the
inclusion of finite-quark-mass effects, as the large-$p_{T}$ tail has a
different scaling behaviour in HEFT. This can be understood from a simple
dimensional analysis. At very high $p_{T}$, the only relevant mass scale is
$p_{T}$ itself. $\mathrm{d}\sigma/\mathrm{d}p_{T}^{2}$ has mass dimension
$-4$, ergo in full QCD the tail must scale as

$\frac{\mathrm{d}\sigma}{\mathrm{d}p_{T}^{2}}\sim p_{T}^{-4}.$ (13)

But in HEFT the cross section must always be proportional to $1/v^{2}$, where
$v$ is the vacuum expectation value, because the Higgs-gluon coupling in HEFT
always gives a factor of $1/v$. Therefore the large-$p_{T}$ tail in HEFT
scales as

$\frac{\mathrm{d}\sigma}{\mathrm{d}p_{T}^{2}}\sim p_{T}^{-2}v^{-2}$ (14)

and is hence less suppressed. Scale-uncertainty bands at large $p_{T}$ show
some slight tension between NLO and NNLO, suggesting that the standard central
scale choice (12) might not be optimal for this process.

Figure 7: Higgs-$p_{T}$ (left) and Higgs-rapidity (right) distribution for the
gluon-gluon fusion channel. Compared are results obtained in HEFT with those
including finite-quark-mass effects in the square of the top-quark-loop
amplitude and in the interference of the top-quark-loop and bottom-quark-loop
amplitudes. The bottom-quark mass is renormalised in the
$\overline{\text{MS}}$ scheme, while the top-quark mass is renormalised in the
OS scheme.

For the Higgs-rapidity distribution in fig. 7 we also observe that the scale-
uncertainty bands of the different pertubative orders do not overlap. For both
full QCD and HEFT, the difference between NLO and NNLO amounts to a rapidity
independent upwards shift of about 20%. Similarly, finite quark masses also
only shift the distribution, this time by around $-4\%$, an effect which is
directly derived from the interference contribution in fig. 6.

##  4 Conclusions

In this work, we studied the impact of using variable flavour number and
different mass-renormalisation schemes on the Higgs production cross section
in the gluon-gluon fusion channel. We found that the two different flavour
schemes are in excellent agreement, finally justifying our treatment of the
bottom-quark mass in the 5FS. Our computation in two different renormalisation
schemes showed that renormalising the bottom-quark mass in the
$\overline{\text{MS}}$ scheme results in significantly smaller scale
uncertainties and displays a better perturbative convergence. The
renormalisation of the top-quark mass on the other hand does not impact the
cross section significantly. We therefore recommend using the 5FS, an
$\overline{\text{MS}}$-renormalised bottom mass and an OS-renormalised top
mass going forward. The cross section for this set-up reads444The
$\mathrm{N}^{3}\mathrm{LO}$ HEFT cross section was computed using SusHi
Harlander:2012pb ; Harlander:2016hcx .

$\boxed{\sigma_{ggH}=48.81(1)^{+0.65}_{-2.02}(\text{N}^{3}\text{LO
HEFT})-0.16^{+0.13}_{-0.03}(\text{NNLO }t)-1.74(2)^{+0.13}_{-0.03}(\text{NNLO
}t\times b)~{}\text{pb}.}$ (15)

We further provided new results for the differential cross section
distribution for the Higgs rapidity, which will be useful in upcoming LHC
analyses, especially for determining the Yukawa couplings of the lighter
quarks.

###### Acknowledgements.

We would like to thank Jonas Lindert for providing us with raw data from ref.
Caola:2018zye . This work was supported by the Deutsche Forschungsgemeinschaft
(DFG) under grant 396021762 - TRR 257: Particle Physics Phenomenology after
the Higgs Discovery, and grant 400140256 - GRK 2497: The Physics of the
Heaviest Particles at the LHC. The authors gratefully acknowledge the
computing time provided to them at the NHR Center NHR4CES at RWTH Aachen
University (project number p0020025). This is funded by the Federal Ministry
of Education and Research, and the state governments participating on the
basis of the resolutions of the GWK for national high performance computing at
universities (www.nhr-verein.de/unsere-partner).

##  Appendix A Comparison to Higgs+jet analyses

In this appendix, we compare our results for the Higgs-$p_{T}$ spectrum with
existing analyses of the Higgs+jet process.

In the OS scheme, we tested our results for the Higgs-$p_{T}$ spectrum against
the fixed-order results presented in ref. Caola:2018zye . In that reference,
the authors computed the two-loop real-virtual corrections under the
assumption of an infinitely heavy top quark and a nearly massless bottom
quark. Fig. 8 shows that these approximations perform well for $p_{T}<40$ GeV,
but result in an error of approximately 20% above this threshold. The nearly
massless bottom quark approximation is expected to be valid if $m_{b}\ll
m_{H},p_{T}$, suggesting that most of the error arises from the HTL
approximation. Additionally, we observe that the upwards scale variations in
ref. Caola:2018zye  are significantly larger in the first two bins ($p_{T}<14$
GeV). After reaching out to the authors, we identified that this discrepancy
is due to one specific scale variation ($\mu_{R}/\mu=1,\mu_{F}/\mu=2$), but we
have not found any errors on our end that would account for this difference.

Figure 8: Top-bottom interference contribution to the Higgs production cross
section. Displayed are our results for the Higgs-$p_{T}$ spectrum and the
results presented in ref. Caola:2018zye , obtained from private communication.
Our data was computed with $m_{t}=173.05$ GeV, $m_{b}=4.78$ GeV and
$m_{H}=125$ GeV, whereas the authors of ref. Caola:2018zye  used $m_{t}=173.2$
GeV, $m_{b}=4.75$ GeV and the same Higgs mass. For the sake of this comparison
we used the same PDF and scale choice as the authors of ref. Caola:2018zye ,
namely PDF4LHC15_nlo_30 and $\mu=\frac{1}{2}\sqrt{m_{H}^{2}+p^{2}_{T}}$.
Transparent bands indicate scale uncertainties, whereas the error bars in the
lower plot indicate the Monte Carlo uncertainties dominated by the
uncertainties of ref. Caola:2018zye .

We also compared our $\overline{\mathrm{MS}}$-renormalised results, i.e. with
both the top- and bottom-quark mass in the $\overline{\mathrm{MS}}$ scheme,
for the Higgs-$p_{T}$ distribution with those in ref. Bonciani:2022jmb . For
the real-virtual corrections, we only computed the amplitudes with a closed
top-quark loop for scales up to $\mu_{R}=125\ \mathrm{GeV}$. Above this scale,
we have to extrapolate outside of our numerical grids. As we can see from fig.
9, we find excellent agreement between the predictions up to transverse
momenta of $400\ \mathrm{GeV}$. Above this scale, we find that the predictions
diverge, which is caused by the poor precision of the extrapolation.

Figure 9: Result for the differential cross section of gluon-gluon fusion and
comparison to ref. Bonciani:2022jmb . The result is computed with finite top-
and bottom-quark mass, where we neglect the cross-section contribution that
does not contain any top-Higgs coupling at NNLO. For the comparison we use the
same computational set-up as Bonciani:2022jmb , i.e. mass values of
$\overline{m}_{t}=163.4$ GeV, $\overline{m}_{b}=4.18$ GeV, $m_{H}=125.25$ GeV,
while the running of the mass is computed at two-loop accuracy. We apply the
anti-$k_{T}$ jet algorithm with $R=0.4$ and use the central scale of
$\mu=\frac{1}{2}\left(\sqrt{m_{H}^{2}+p_{T}^{2}}+\sum_{i}p_{i,T}\right)$,
where the sum runs over all jet momenta. We apply a cut on the transverse
momentum of the jets, requiring that at least one of the jet transverse
momenta satisfies $p_{i,T}>20$ GeV. We use the NNPDF40_nlo_as_01180 PDF set.

## References

  * (1) ATLAS collaboration, _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 [1207.7214].
  * (2) CMS collaboration, _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30 [1207.7235].
  * (3) CMS collaboration, _A portrait of the Higgs boson by the CMS experiment ten years after the discovery._ , _Nature_ 607 (2022) 60 [2207.00043].
  * (4) S. Dawson, _Radiative corrections to Higgs boson production_ , _Nucl. Phys. B_ 359 (1991) 283.
  * (5) A. Djouadi, M. Spira and P.M. Zerwas, _Production of Higgs bosons in proton colliders: QCD corrections_ , _Phys. Lett. B_ 264 (1991) 440.
  * (6) V. Ravindran, J. Smith and W.L. Van Neerven, _Next-to-leading order QCD corrections to differential distributions of Higgs boson production in hadron hadron collisions_ , _Nucl. Phys. B_ 634 (2002) 247 [hep-ph/0201114].
  * (7) R.V. Harlander and W.B. Kilgore, _Next-to-next-to-leading order Higgs production at hadron colliders_ , _Phys. Rev. Lett._ 88 (2002) 201801 [hep-ph/0201206].
  * (8) C. Anastasiou and K. Melnikov, _Higgs boson production at hadron colliders in NNLO QCD_ , _Nucl. Phys. B_ 646 (2002) 220 [hep-ph/0207004].
  * (9) C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, _Higgs Boson Gluon-Fusion Production in QCD at Three Loops_ , _Phys. Rev. Lett._ 114 (2015) 212001 [1503.06056].
  * (10) B. Mistlberger, _Higgs boson production at hadron colliders at N 3LO in QCD_, _JHEP_ 05 (2018) 028 [1802.00833].
  * (11) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog et al., _High precision determination of the gluon fusion Higgs boson cross-section at the LHC_ , _JHEP_ 05 (2016) 058 [1602.00695].
  * (12) M. Czakon, R.V. Harlander, J. Klappert and M. Niggetiedt, _Exact Top-Quark Mass Dependence in Hadronic Higgs Production_ , _Phys. Rev. Lett._ 127 (2021) 162002 [2105.04436].
  * (13) M. Czakon, F. Eschment, M. Niggetiedt, R. Poncelet and T. Schellenberger, _Top-Bottom Interference Contribution to Fully-Inclusive Higgs Production_ , 2312.09896.
  * (14) M. Niggetiedt and M. Wiesemann, _Higgs-boson production in the full theory at NNLO+PS_ , 2407.01354.
  * (15) F. Bishara, U. Haisch, P.F. Monni and E. Re, _Constraining Light-Quark Yukawa Couplings from Higgs Distributions_ , _Phys. Rev. Lett._ 118 (2017) 121801 [1606.09253].
  * (16) G. Bonner and H.E. Logan, _Constraining the Higgs couplings to up and down quarks using production kinematics at the CERN Large Hadron Collider_ , 1608.04376.
  * (17) Y. Soreq, H.X. Zhu and J. Zupan, _Light quark Yukawa couplings from Higgs kinematics_ , _JHEP_ 12 (2016) 045 [1606.09621].
  * (18) CMS collaboration, _Measurement and interpretation of differential cross sections for Higgs boson production at $\sqrt{s}=$ 13 TeV_, _Phys. Lett. B_ 792 (2019) 369 [1812.06504].
  * (19) CMS collaboration, _Measurements of inclusive and differential cross sections for the Higgs boson production and decay to four-leptons in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, _JHEP_ 08 (2023) 040 [2305.07532].
  * (20) M. Grazzini and H. Sargsyan, _Heavy-quark mass effects in Higgs boson production at the LHC_ , _JHEP_ 09 (2013) 129 [1306.4581].
  * (21) K. Melnikov and A. Penin, _On the light quark mass effects in Higgs boson production in gluon fusion_ , _JHEP_ 05 (2016) 172 [1602.09020].
  * (22) J.M. Lindert, K. Melnikov, L. Tancredi and C. Wever, _Top-bottom interference effects in Higgs plus jet production at the LHC_ , _Phys. Rev. Lett._ 118 (2017) 252002 [1703.03886].
  * (23) F. Caola, J.M. Lindert, K. Melnikov, P.F. Monni, L. Tancredi and C. Wever, _Bottom-quark effects in Higgs production at intermediate transverse momentum_ , _JHEP_ 09 (2018) 035 [1804.07632].
  * (24) R. Bonciani, V. Del Duca, H. Frellesvig, M. Hidding, V. Hirschi, F. Moriello et al., _Next-to-leading-order QCD corrections to Higgs production in association with a jet_ , _Phys. Lett. B_ 843 (2023) 137995 [2206.10490].
  * (25) M. Czakon, _A novel subtraction scheme for double-real radiation at NNLO_ , _Phys. Lett. B_ 693 (2010) 259 [1005.0274].
  * (26) M. Czakon and D. Heymes, _Four-dimensional formulation of the sector-improved residue subtraction scheme_ , _Nucl. Phys. B_ 890 (2014) 152 [1408.2500].
  * (27) M. Czakon, A. van Hameren, A. Mitov and R. Poncelet, _Single-jet inclusive rates with exact color at $\mathcal{O}$ (${\alpha}_{s}^{4}$)_, _JHEP_ 10 (2019) 262 [1907.12911].
  * (28) L. Budge, J.M. Campbell, G. De Laurentis, R.K. Ellis and S. Seth, _The one-loop amplitudes for Higgs + 4 partons with full mass effects_ , _JHEP_ 05 (2020) 079 [2002.04018].
  * (29) J. Campbell and T. Neumann, _Precision Phenomenology with MCFM_ , _JHEP_ 12 (2019) 034 [1909.09117].
  * (30) J.M. Campbell and R.K. Ellis, _An Update on vector boson pair production at hadron colliders_ , _Phys. Rev. D_ 60 (1999) 113006 [hep-ph/9905386].
  * (31) J.M. Campbell, R.K. Ellis and C. Williams, _Vector boson pair production at the LHC_ , _JHEP_ 07 (2011) 018 [1105.0020].
  * (32) R.K. Ellis and G. Zanderighi, _Scalar one-loop integrals for QCD_ , _JHEP_ 02 (2008) 002 [0712.1851].
  * (33) S. Carrazza, R.K. Ellis and G. Zanderighi, _QCDLoop: a comprehensive framework for one-loop scalar integrals_ , _Comput. Phys. Commun._ 209 (2016) 134 [1605.03181].
  * (34) M. Czakon and M. Niggetiedt, _Exact quark-mass dependence of the Higgs-gluon form factor at three loops in QCD_ , _JHEP_ 05 (2020) 149 [2001.03008].
  * (35) M. Niggetiedt and J. Usovitsch, _The Higgs-gluon form factor at three loops in QCD with three mass scales_ , 2312.05297.
  * (36) P. Maierhöfer, J. Usovitsch and P. Uwer, _Kira—A Feynman integral reduction program_ , _Comput. Phys. Commun._ 230 (2018) 99 [1705.05610].
  * (37) P. Maierhöfer and J. Usovitsch, _Kira 1.2 Release Notes_ , 1812.01491.
  * (38) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024 [2008.06494].
  * (39) J. Klappert and F. Lange, _Reconstructing rational functions with FireFly_ , _Comput. Phys. Commun._ 247 (2020) 106951 [1904.00009].
  * (40) J. Klappert, S.Y. Klein and F. Lange, _Interpolation of dense and sparse rational functions and other improvements in FireFly_ , _Comput. Phys. Commun._ 264 (2021) 107968 [2004.01463].
  * (41) N. Gray, D.J. Broadhurst, W. Grafe and K. Schilcher, _Three Loop Relation of Quark (Modified) Ms and Pole Masses_ , _Z. Phys. C_ 48 (1990) 673.
  * (42) R. Tarrach, _The pole mass in perturbative qcd_ , _Nuclear Physics B_ 183 (1981) 384.
  * (43) B. Schmidt and M. Steinhauser, _CRunDec: a C++ package for running and decoupling of the strong coupling and quark masses_ , _Comput. Phys. Commun._ 183 (2012) 1845 [1201.6149].
  * (44) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (45) S. Actis, A. Denner, L. Hofer, J.-N. Lang, A. Scharf and S. Uccirati, _RECOLA: REcursive Computation of One-Loop Amplitudes_ , _Comput. Phys. Commun._ 214 (2017) 140 [1605.01090].
  * (46) A. Denner, J.-N. Lang and S. Uccirati, _Recola2: REcursive Computation of One-Loop Amplitudes 2_ , _Comput. Phys. Commun._ 224 (2018) 346 [1711.07388].
  * (47) S. Dittmaier, M. Krämer and M. Spira, _Higgs radiation off bottom quarks at the Tevatron and the CERN LHC_ , _Phys. Rev. D_ 70 (2004) 074010 [hep-ph/0309204].
  * (48) R.V. Harlander, S. Liebler and H. Mantler, _SusHi: A program for the calculation of Higgs production in gluon fusion and bottom-quark annihilation in the Standard Model and the MSSM_ , _Comput. Phys. Commun._ 184 (2013) 1605 [1212.3249].
  * (49) R.V. Harlander, S. Liebler and H. Mantler, _SusHi Bento: Beyond NNLO and the heavy-top limit_ , _Comput. Phys. Commun._ 212 (2017) 239 [1605.03190].
  * (50) P. Pietrulewicz and M. Stahlhofen, _Two-loop bottom mass effects on the Higgs transverse momentum spectrum in top-induced gluon fusion_ , _JHEP_ 05 (2023) 175 [2302.06623].
  * (51) LHC Higgs Cross Section Working Group collaboration, _Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector_ , 1610.07922.
  * (52) S.P. Jones, M. Kerner and G. Luisoni, _Next-to-Leading-Order QCD Corrections to Higgs Boson Plus Jet Production with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 120 (2018) 162001 [1802.00349].


11institutetext:  Karlsruhe Institute of Technology (KIT) 22institutetext:
Observatorio Pierre Auger, Av. San MartÃn Norte 304, 5613 MalargÃe, Argentina

# Scintillator Surface Detector simulations for AugerPrime

David Schmidt for the Pierre Auger Collaboration Authors list:
http://www.auger.org/archive/authors_2018_10.html11<EMAIL_ADDRESS>22
<EMAIL_ADDRESS>
###### Abstract

Knowledge of the mass composition of ultra-high-energy cosmic rays is
understood to be a salient component in answering the open questions in the
field. The AugerPrime upgrade of the Pierre Auger Observatory aims to enhance
its surface detector with the hardware necessary to reconstruct primary mass
for individual events. This involves placing a scintillation-based detector
with an active area of $3.8\text{\,}{\mathrm{m}}^{2}$ on top of each existing
water-Cherenkov detector in its surface detector array. Here, we present the
methods for simulating this Scintillator Surface Detector. These simulations
have and will continue to aid in the interpretation of measurements with
AugerPrime as well as the development and improvement of event reconstruction
algorithms including primary mass.

##  1 Introduction

The Pierre Auger Observatory (Auger) Auger , with its unprecedented exposure
of over
$60\,000\text{\,}{\mathrm{km}}^{2}\text{\,}\mathrm{sr}\text{\,}\mathrm{yr}$
acquired during more than decade of data collection, has led to progress in
the field of Ultra-High-Energy Cosmic Rays (UHECR) physics in a number of
ways. To name a couple, the suppression in the energy spectrum above
approximately $40\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$ has been confirmed
to high precision Cutoff1 ; Cutoff2 , and a large-scale dipole anisotropy
above $8\text{\,}\mathrm{E}\mathrm{e}\mathrm{V}$ has been clearly observed
Dipole . With these advances, however, Auger has also uncovered additional and
unexpected complexity. Perhaps most notable is the trend towards heavier
composition at the highest energies, which reopens questions regarding the
origins of the flux suppression. Additionally, the long standing mystery as to
the origins of UHECRs remains unsolved, and the additional complication of
mass now presents further challenges. For this reason, the observatory is
embarking on its next phase by equipping its Surface Detector (SD) array with
the hardware necessary to estimate the mass of UHECRs on an event-by-event
basis.

This upgrade, known as AugerPrime AugerPrime ; AugerPrimeICRC ;
AugerPrimeUHECR , capitalizes on the fact that the magnitude of the muonic
component of Extensive Air Showers produced by UHECRs scales with primary
mass. AugerPrime aims to disentangle the contributions of the electromagnetic
and muonic shower components in SD measurements, thereby obtaining a mass-
sensitive estimator. This is being accomplished through the installation of an
Scintillator Surface Detector (SSD) on top of each of the existing Water-
Cherenkov Detectors that make up the SD array. The relatively higher
sensitivity of the SSDs (WCDs) to particles of the electromagnetic (muonic)
air shower components allows for this disentanglement.

A prototype array of SSDs has been operating in the field since September of
2016, and deployment of production design SSDs for the complete array is
currently underway. In order to aid in understanding and interpreting
measurements of the upgraded SD, as well as develop reconstruction algorithms,
we have developed comprehensive SSD simulations. The methods used in these
simulations are described in this proceeding and build upon AugerPrimeOffline
.

##  2 Scintillator Surface Detector

The active area of an SSD is comprised of 48 polystyrene scintillator bars
measuring $1600\text{\,}\mathrm{mm}$ by $50\text{\,}\mathrm{mm}$ with a
thickness of $10\text{\,}\mathrm{mm}$. This makes for an active area of just
over $3.8\text{\,}\mathrm{m}^{2}$, which is split between two wings of 24 bars
each (see Fig. 1). $1\text{\,}\mathrm{mm}$ diameter Wave-Length Shifting (WLS)
fibers of $5.8\text{\,}\mathrm{m}$ length are routed through two bean shaped
holes inside each bar with a U-bend of $10\text{\,}\mathrm{cm}$ diameter at
the ends furthest from the central axis of the SSD such that both fiber ends
may be routed into a cookie situated between the two scintillator wings. There
are approximately $1.25\text{\,}\mathrm{m}$ of fiber between the central ends
of the scintillator bars and the cookie, and the lengths of the fibers are
uniform for all scintillator bars through which they are routed. The bars and
fibers are housed within an aluminum casing whose bottom consists of a
composite panel with $22\text{\,}\mathrm{mm}$ of Extruded Polystyrene (XPS)
situated between two $1\text{\,}\mathrm{mm}$ sheets of aluminum. The top of
the housing consists of a $1\text{\,}\mathrm{mm}$ sheet of aluminum. The space
within the SSD housing not occupied by the bars and fibers is mostly filled
with Expanded Polystyrene (EPS), resulting in a volume of air inside the
detector of less than $10\text{\,}\mathrm{L}$. For a detailed description of
the detector materials, geometry, and construction, see SSDICRC2017 .

Figure 1: Top-down schematic depiction of the arrangement of bars and fibers
within the SSD casing. The written dimensions describe the active surfaces of
the detectors, namely the scintillator bars.

##  3 Simulation Application

The SSD simulation application is housed within Auger’s official simulation
and reconstruction framework, $\overline{\text{Off}}$
$\underline{\text{line}}$ Offline , which was upgraded to accommodate for the
demands of the large-scale detector upgrade that is AugerPrime
AugerPrimeOffline . The SSD simulations are based in the field-standard
simulation software Geant4 Geant4  (version 4.10), but, where it is possible,
make use of measurements of SSDs performed with a centimeter-precision muon
telescope for enrichment and tuning. The combined use of Geant4 and the muon
telescope measurements has made highly detailed simulation of the detector
possible. The simulations account for the energy-loss processes of traversing
particles, the decays occurring from excited states of atoms within the
scintillator bars and WLS fibers, the attenuation of photons along the fibers,
the Photoelectrons generated at the photo-cathode of the Photomultiplier Tube
(PMT), the current at the base of the PMT, and the response of the AugerPrime
station electronics UUB . The following sections step through the construction
of the detector in Geant4 and each stage of the simulation.

Figure 2: Visualization of both the SSD and WCD volumes defined in the
Geant4-based simulation application with a close-up of one corner of the SSD

###  3.1 Detector Construction

SSDs are constructed within the same Geant4 world volume as their partner
WCDs. The volumes for both detectors may be observed in Fig. 2. Regarding the
SSD, volumes for the aluminum casing, XPS, EPS, and scintillator bars
themselves are defined according to the dimensions and material properties of
the production SSD design, described in Sec. 2. This is important to
accurately reproduce the shielding above and below the active area of the
scintillator bars. Precise definition of the position of the SSD relative to
the WCD is important for preserving signal correlations, as these depend on
the fraction of particles that intersect both detectors.

Figure 3: Depiction of the time evolution of the signal from single, through-
going muons in an SSD. Muon telescope measurements are shown in the plots on
the left and the parameterization of these measurements are shown in those on
the right. Middle: Density plot of the signal evolution over one scintillator
wing of the SSD. The bottoms of the histograms correspond to a particle
crossing at the non-central edges of the scintillator bars, where photons
corresponding to the two legs of the signal traverse approximately equal
distance and therefore produce a unimodal distribution. The tops of the
histograms correspond to the central edges of the scintillator bars, where the
photons corresponding to the two legs of signal traverse significantly
different distances and therefore result in a bimodal distribution. Slices in
the horizontal axis have been normalized to the same maximum. Top and bottom:
Signal shapes at the central and non-central ends of the scintillator bars.

###  3.2 Photoelectrons

For each particle that intersects or is generated within the active
scintillator volumes, the energy deposited within the bars is extracted from
Geant4 and converted into a number of photoelectrons at the photo-cathode of
the PMT. This conversion must take into account the attenuation of photons
within the bars and along the fibers as well as the excitations and decays
corresponding to scintillation and wavelength-shifting processes. The
conversion takes the form of

$\bar{N}(\textbf{x})=N_{\mathrm{ref}}\frac{E_{\mathrm{dep}}}{E_{\mathrm{ref}}}f_{\mathrm{att}}(\textbf{x})$
(1)

where $E_{\mathrm{dep}}$ is the energy deposit in the scintillator bars,
$E_{\mathrm{ref}}$ and $N_{\mathrm{ref}}$ are the reference energy and PE
number derived from simulations, and $f_{\mathrm{att}}(\textbf{x})$ describes
the degree of photon attenuation as a function of particle crossing position
x. $E_{\mathrm{ref}}$ is the peak in the distribution of energy deposited in
the SSD by an ensemble of vertical, through-going muons simulated for the
volumes described in Sec. 3.1. $N_{\mathrm{ref}}$ is a constant derived such
that average value of Eqn. 1 when integrated over all possible values of x
inside the scintillator bars is equal to the the average number of PEs
observed in real measurements of the SSDs with the aforementioned muon
telescope. For a more detailed description of how these parameters were
obtained, see AugerPrimeOffline . With this construction, the mean number of
PEs for a single, through-going muon is guaranteed to match real measurements.
The attenuation function $f_{\mathrm{att}}(\textbf{x})$ is written as

$f_{\mathrm{att}}(\textbf{x})=A(\textbf{x})L(\textbf{x})\ ,$ (2)

where $A(\textbf{x})$ describes the attenuation along the WLS fibers. This
attenuation is of the form

$A(\textbf{x})=e^{-\ell_{\mathrm{c}}(\textbf{x})/\lambda_{\mathrm{f}}}+e^{-\ell\mathrm{f}(\textbf{x})/\lambda_{\mathrm{f}}}\
,$ (3)

where the two exponential terms correspond to the two paths photons may travel
to reach the PMT due to the U-bend routing of the fibers. $\ell_{\mathrm{c}}$
($\ell_{\mathrm{f}}$) corresponds to the shorter (longer) travel path and is
calculated for a given crossing position x. $\lambda_{\mathrm{f}}$ is the
effective attenuation length of the fiber obtained through analysis of muon
telescope SSD measurements. The term $L(\textbf{x})$ corresponds to a
decreased yield at crossing positions very close the edges of the scintillator
bars. Thus, for any energy deposit and particle crossing position, an expected
number of PEs corresponding to each of the two paths photons may travel to
reach the PMT may be calculated. A Poisson randomization is then performed on
each of these expectations. For each PE coming out of these randomizations,
two draws, one from each of two exponential distributions are performed. This
model of two exponential decays, which may arbitrarily be attributed to the
excitation associated with scintillation in the scintillator bars and the
wave-length shifting process in the WLS fibers, was found to provide a
sufficient description of the data AugerPrimeOffline . The production time for
each PE is then obtained via

$t_{\mathrm{PE}}=t_{\mathrm{global}}+t_{\mathrm{kin}}+t_{\mathrm{bar}}+t_{\mathrm{fiber}}$
(4)

where $t_{\mathrm{global}}$ is the global time obtained from Geant4,
$t_{\mathrm{kin}}$ is the kinematic delay of the photon as it traverses its
respective path to the PMT, and $t_{\mathrm{bar}}$ and $t_{\mathrm{fiber}}$
are the times arising from the two aforementioned exponential decays. The
average measured and parameterized time distributions for crossing positions
along the complete active surface of one wing of the SSD, which compound Eqn.
1, Eqn. 4, and the convolution of the two exponential decays, may be observed
in Fig. 3.

Figure 4: Left: Average charge measured for single, through-going muons as a
function of the crossing position for one wing of an SSD. Measurements with
the muon telescope are shown by the square, red markers, whereas simulated
data is represented by the circular, black markers. Larger values on the
horizontal axis correspond to crossing positions further from the central edge
of the scintillator, as depicted by the semi-transparent schematic in the
backdrop. Right: Average charge as a function of zenith angle. Once more, muon
telescope measurements are given by the red markers, whereas simulations are
shown in black. Both data sets have been scaled such that they equal one at an
angle of $35\text{\,}\mathrm{\SIUnitSymbolDegree}$ to compare their evolution
with tracklength ($1/\mathrm{cos}\theta$).

###  3.3 PMT and electronics

With a trace of PE production times in hand, the PMT is simulated as follows.
For each PE, the average pulse shape at the PMT base corresponding to an
Single Photoelectron (SPE) is scaled by a random draw from the corresponding
charge distribution. Both this pulse shape and the charge distribution were
measured for SPE events for the PMT in question, as described in
AugerPrimeOffline . Performing this process for each PE results in a time
distribution of the electrical current at the base of the PMT, which is the
input to the upgraded SD station electronics UUB .

The effect of the electronics on the time trace is simulated by convolving the
base current distribution with the measured electronics transfer function.
This results in an analog voltage trace, which is then scaled, sampled, and
digitized. The scaling is performed such that the resulting peak in the charge
distribution for a simulated ensemble of through-going muons equals that which
has been measured by prototype SSDs deployed in an Engineering Array (EA) EA
at the observatory site near MalargÃe since 2016. The sampling and
digitization are performed in accordance with the $120\text{\,}\mathrm{MHz}$
12 bit design of the AugerPrime electronics.

##  4 Validation

Validation of the output of the simulations was performed on a per-particle
basis by comparing the position-dependent charge and time distributions of
simulated muons to those measured by the muon telescope. In Fig. 4, the
dependencies of charge on crossing position and zenith angle of the particles
are compared between measurements and simulations. In both cases, the
simulations provide a reasonable description of the data. A full validation
closely examining the shapes of the charge and time distributions is currently
underway as are additional comparisons with measurements of SSDs deployed in
the EA.

Figure 5: Simulated SSD measurement. Time traces for PE production, current at
the base of the PMT, post-electronics digitized signals, and calibrated to the
signal of one VMIP are shown.

##  5 Example Application

For a given simulated extensive air shower, the particles reaching ground near
each SD station are injected into a cylindrical volume that completely encases
both the SSD and WCD. Propagation of these incident particles into the shared
volume and subsequent interactions and generation of new particles is then
left to the Geant4 software, and signals are simulated as described in Sec. 3.
A depiction of time distributions of photoelectrons, current at the base of
the PMT, and digital signals after electronics processing are shown in Fig. 5
for an example station simulation. An example event, where SSD and WCD pairs
sample the lateral distribution of particles at ground at various locations,
based on Auger’s hexagonal grid, is depicted in Fig. 6.

Figure 6: Simulated signals of SSDs alongside partner WCDs for various samples
of the lateral distribution of the particles reaching ground in a simulated
EAS.

##  6 Conclusions

Detector simulations for the SSD of the AugerPrime upgrade have been developed
and the methods, making use of Geant4 and real detector measurements, have
been described here. Initial validation has shown general agreement between
the detector simulations and centimeter-precision muon telescope measurements
for the charge and time distributions produced by single, through-going muons.
Additional, more detailed, validation is currently being performed. The
combined SSD and WCD simulation application, housed in the
$\overline{\text{Off}}$ $\underline{\text{line}}$ detector simulation and
reconstruction software framework of Auger, has already and will continue to
aid in the interpretation of measurements performed with AugerPrime as well as
in the development and improvement of event reconstruction algorithms which
include primary mass.

## References

  * (1) The Pierre Auger Collaboration, Nucl. Instrum. Meth. A 798, 172-213 (2015)
  * (2) The Pierre Auger Collaboration, Phys. Rev. Lett. 101, 061101 (2008)
  * (3) The Pierre Auger Collaboration, Phys. Lett. B 685, 239-246 (2010)
  * (4) The Pierre Auger Collaboration, Science 357, 1266-1270 (2017)
  * (5) The Pierre Auger Collaboration, arXiv:1604.03637 (2016)
  * (6) D. Martello for the Pierre Auger Collaboration, Proc. 35th Int. Cosmic Ray Conf. PoS(ICRC2017) 486, (2017)
  * (7) A. Castellina for the Pierre Auger Collaboration, these Proceedings
  * (8) D. Schmidt for the Pierre Auger Collaboration, PoS ICRC2017, 353 (2017)
  * (9) R. Smida for the Pierre Auger Collaboration, PoS ICRC2017, 390 (2017)
  * (10) S. Argiro et al., Nucl. Instrum. Meth. A 580, 1485-1496 (2007)
  * (11) S. Agostinelli et al., Nucl. Instrum. Meth. A 506, 250-303 (2003)
  * (12) T. SuomijÃrvi for the Pierre Auger Collaboration, PoS ICRC2017, 450 (2017)
  * (13) Z. Zong for the Pierre Auger Collaboration, PoS ICRC2017, 449 (2017)

  *[Auger]: The Pierre Auger Observatory
  *[UHECR]: Ultra-High-Energy Cosmic Rays
  *[UHECRs]: Ultra-High-Energy Cosmic Rays
  *[SD]: Surface Detector
  *[SSD]: Scintillator Surface Detector
  *[SSDs]: Scintillator Surface Detector
  *[WCDs]: Water-Cherenkov Detector
  *[WLS]: Wave-Length Shifting
  *[XPS]: Extruded Polystyrene
  *[EPS]: Expanded Polystyrene
  *[PMT]: Photomultiplier Tube
  *[WCD]: Water-Cherenkov Detector
  *[PE]: Photoelectron
  *[PEs]: Photoelectron
  *[SPE]: Single Photoelectron
  *[EA]: Engineering Array
  *[VMIP]: Vertical Minimum Ionizing Particle
  *[EAS]: Extensive Air Shower