††institutetext: Institute for Theoretical Particle Physics and Cosmology,
RWTH Aachen University,
D-52056 Aachen, Germany

# Four top final states with NLO accuracy in perturbative QCD: 4 lepton
channel

Nikolaos Dimitrakopoulos  and Malgorzata Worek<EMAIL_ADDRESS>aachen.de<EMAIL_ADDRESS>
###### Abstract

Triggered by the observation of four top-quark production at the LHC by the
ATLAS and CMS collaboration we report on the calculation of the next-to-
leading order QCD corrections to the Standard Model process $pp\to
t\bar{t}t\bar{t}$ in the $4\ell$ top-quark decay channel. We take into account
higher-order QCD effects in both the production and decays of the four top
quarks. The latter effects are treated in the narrow width approximation,
which preserves top-quark spin correlations throughout the calculation. We
present results for two selected renormalisation and factorisation scale
settings and three different PDF sets. Furthermore, we study the main
theoretical uncertainties that are associated with the neglected higher-order
terms in the perturbative expansion and with the parameterisation of the PDF
sets. The results at the integrated and differential fiducial cross-section
level are shown for the LHC Run III center-of-mass energy of $\sqrt{s}=13.6$
TeV. Our findings are particularly relevant for precise measurements of the
four top-quark fiducial cross sections and for the modelling of top-quark
decays at the LHC.

###### Keywords:

Higher-Order Perturbative Calculations, Specific QCD Phenomenology, Top Quark

††dedication: TTK-24-02, P3H-24-001

##  1 Introduction

Being the heaviest known elementary particle of the Standard Model (SM), the
top quark has a special relation with the SM Higgs boson due to the magnitude
of its Yukawa coupling $(Y_{t})$ which is of the order of one. The top quark
is also predicted to have large couplings to hypothetical new particles that
might appear in many models beyond the SM (BSM). In this respect, rare
processes involving the top quark are particularly important to be studied at
the Large Hadron Collider (LHC). Not only can they explore the parameters of
the SM in regions hitherto inaccessible, but also, if the BSM couplings to the
top quark are indeed significant, serve as golden candidates to reveal signs
of new physics. Production of four top-quarks in $pp$ collisions, that
requires a partonic center-of-mass energy of at least $4m_{t}$, is one of the
rarest and heaviest processes accessible at the LHC. Although the $pp\to
t\bar{t}t\bar{t}$ process occurs predominantly through the strong interaction,
non-negligible contributions arise also from electroweak (EW) processes
occurring through the exchange of the electroweak gauge and/or SM Higgs boson.
Due to EW contributions involving the Higgs boson exchange, the $pp\to
t\bar{t}H(H\to t\bar{t})$ process can be used to constrain the Higgs boson
Yukawa coupling Cao:2016wib ; Cao:2019ygh , complementary to loop sensitivity
through electroweak corrections in $t\bar{t}$ production Schmidt:1992et ;
Kuhn:2013zoa ; Martini:2021uey  and tree-level sensitivity in the $pp\to
t\bar{t}H$ and/or $pp\to tH$ process Artoisenet:2013puc ; Demartin:2014fia ;
Demartin:2015uha ; Demartin:2016axk ; Bahl:2020wee ; Bortolato:2020zcg ;
Hermann:2022vit . Moreover, in various new physics scenarios, the small SM
cross section for the production of four top-quarks can be significantly
enhanced. These BSM scenarios include among others models with a composite
right-handed top quark Lillie:2007hd ; Pomarol:2008bh ; Kumar:2009vs , heavy
Kaluza-Klein gluons and quarks from extra dimensions Cacciapaglia:2011kz ,
scalar gluons occurring in supersymmetric models and produced in pairs
Plehn:2008ae ; Calvet:2012rk , singly produced top-philic resonances
Kim:2016plm ; Darme:2021gtt  or a heavy scalar or pseudoscalar boson from
extended electroweak symmetry breaking sectors produced in association with a
top-quark pair BhupalDev:2014bir ; Craig:2015jba ; Gori:2016zto ;
Craig:2016ygr . In addition, the $pp\to t\bar{t}t\bar{t}$ process provides a
direct way to constrain the SM Effective Field Theory Wilson coefficients
sensitive to the quartic couplings between top-quarks Zhang:2017mls ;
Ethier:2021bye ; Aoude:2022deh . If the scale of new physics is too high to be
observed directly at the LHC, it can manifest itself as a modification of the
SM $t\bar{t}t\bar{t}$ cross section and cause modifications in various
differential cross-section distributions. Given such rich phenomenology of the
SM $pp\to t\bar{t}t\bar{t}$ process at the LHC, the measurements of four-top
integrated and differential cross sections must be accompanied by very precise
and accurate theoretical predictions for the process. Such theoretical
predictions should be a prerequisite for shedding light on possible signals of
new physics that may emerge in this rare channel.

With its short lifetime, however, the top quark decays before top-flavoured
hadrons or $t\bar{t}$-bound states can be formed. We therefore only have
access to the properties of the top quark through its decay products. Given
that the top quark decays almost $100\%$ of the time into $t\to Wb$, its final
states are classified via the $W$ gauge boson decays into $W\to\ell\nu_{\ell}$
and $W\to q\bar{q}^{\,\prime}$. Consequently, the SM $pp\to
t\bar{t}t\bar{t}\to W^{+}W^{-}W^{+}W^{-}b\bar{b}b\bar{b}$ process produces a
plethora of final states, all of which yield interesting and distinctive
signatures at the LHC. The dominant $4t$ decay mode is the mono-leptonic
channel with a branching ratio $({\cal BR})$ of $40\%$. It is followed by the
fully hadronic and opposite-sign dilepton modes with ${\cal BR}\approx 20\%$
each, the same-sign dilepton and the trilepton modes with ${\cal BR}\approx
10\%$ each and finally by the fully leptonic mode with the smallest
contribution of the order of $1\%$ Alvarez:2016nrz . Both ATLAS and CMS have
recently observed the production of four top-quarks with the LHC Run II energy
of $13$ TeV using the integrated luminosity of about $140$ fb-1 ATLAS:2023ajo
; CMS:2023ftu . Events with two same-sign, three, or four charged leptons
(electrons and muons) and additional light- and/or $b$-jets have been analysed
to achieve the required significance of more than $5\sigma$ that is needed to
claim observation. Independent of the experiment and taking into account
current uncertainties of the measurements, the obtained $t\bar{t}t\bar{t}$
cross sections are in agreement with the SM prediction. At present, however,
both ATLAS and CMS measurements are limited by statistical uncertainties. The
largest systematic uncertainties, on the other hand, come from the modelling
of the additional jet activity in $t\bar{t}W^{\pm}$ production, the
normalisation of the $pp\to t\bar{t}Z$ background process and the modelling of
the $pp\to t\bar{t}t\bar{t}$ signal process. Therefore, an accurate
theoretical description of $t\bar{t}t\bar{t}$ production at the LHC comprising
all the features of the process is much needed and more relevant than ever.

On the theory side, NLO QCD corrections for the $pp\to t\bar{t}t\bar{t}$
process with stable top quarks have been calculated for the first time a while
ago Bevilacqua:2012em  and a few years later recomputed in Refs.
Alwall:2014hca ; Maltoni:2015ena . Besides NLO QCD corrections, a further step
towards more precise predictions for this process has been achieved by
evaluating the complete-NLO corrections Frederix:2017wme  that additionally
include NLO EW higher-order effects and subleading contributions at
perturbative orders from ${\cal O}(\alpha_{s}^{5})$ to ${\cal O}(\alpha^{5})$.
However, it has been shown that there are cancellations between different
terms in $\alpha_{s}$ and $\alpha$, so in general the higher-order effects for
the $pp\to t\bar{t}t\bar{t}$ process are dominated by $\alpha_{s}$ corrections
to the leading QCD process at ${\cal O}(\alpha_{s}^{4})$. Finally, very
recently the calculation for $pp\to t\bar{t}t\bar{t}$ at the next-to-leading
logarithmic accuracy has been carried out vanBeekveld:2022hty  including also
constant ${\cal O}(\alpha_{s})$ non-logarithmic contributions that do not
vanish at the absolute production threshold defined as
$M^{2}/\hat{s}\rightarrow 1$, where $M=4t$. In all these studies only the
magnitude of the higher-order effects in the total production rate of four top
quarks has been investigated, while the reliable description of the fiducial
cross sections is still lacking. For realistic studies, the final state with
$12$ particles must be considered that consists of four $b$-jets, additional
light jets and/or charged leptons together with the (significant) missing
transverse momentum. Here, too, reliable theoretical predictions are needed.
The first attempt in this direction has already been accomplished in Ref.
Jezo:2021smh  by matching theoretical predictions for the $pp\to
t\bar{t}t\bar{t}$ process with the dipole-style $p_{T}$-ordered parton shower
from the Pythia8 Monte Carlo (MC) program Sjostrand:2006za ; Sjostrand:2014zea
using the Powheg framework Nason:2004rx ; Alioli:2010xd . In this study, top-
quark decays in the $\ell+jets$ channel were modelled at the LO level,
preserving spin correlation effects. For the first time, not only the impact
of NLO QCD corrections and $t\bar{t}$ spin correlation effects, but also the
importance of the subleading EW production modes, that have been included with
the LO accuracy only, has been investigated at the differential fiducial
cross-section level. In addition, a comparison with the results obtained
within the MG5-aMC@NLO framework Alwall:2014hca , that employs the MC@NLO
matching to parton shower programs Frixione:2002ik ; Frixione:2003ei , has
been performed. Basically, parton showers represent higher-order corrections
to the hard process. However, in this approach an approximation scheme is used
in which only colinear parton splitting or soft gluon emission contributions
are included. Thus, to cover the entire fiducial phase space it is important
to incorporate higher-order effects in top-quark decays already at the matrix
element level. This would represent a further step towards the complete
description of the jet radiation pattern for this process.

Unfortunately, NLO QCD predictions for $pp\to t\bar{t}t\bar{t}$ that are valid
at the level of observable particles such as leptons, $b$-jets and light-jets
originating either from the production of top quarks or in their decays would
require the calculation of higher order corrections to $2\to 12$ process,
which is a formidable task at present. Fortunately, the problem can be
simplified by examining only the (quadruple) resonance contributions for top-
quarks and $W$ gauge bosons while including spin degrees of freedom exactly.
Indeed, all QCD corrections for the $pp\to t\bar{t}t\bar{t}$ process can be
decomposed into factorizable and non-factorizable contributions. The latter
corrections imply a cross-talk between $t\bar{t}t\bar{t}$ production and
$t/\bar{t}$ decays as well as between the $t$ and $\bar{t}$ decays. In the
limit when $\Gamma/m\to 0$, these corrections must vanish. The exact way in
which they vanish has already been described ample times in literature, see
e.g. Fadin:1993kt ; Melnikov:1993np ; Melnikov:1995fx ; Denner:1997ia ;
Beenakker:1999ya ; Melnikov:2009dn ; Campbell:2012uf ; Bevilacqua:2019quz ;
Czakon:2020qbd , and goes under the name of the narrow-width-approximation
(NWA). The aim of this paper is to provide theoretical predictions for the
$pp\to t\bar{t}t\bar{t}$ process that include NLO QCD corrections to
$t\bar{t}t\bar{t}$ production and four top-quark decays in the NWA. In our
study, we consider the $4\ell$ top-quark decay channel comprising the
following decay chain $pp\to t\bar{t}t\bar{t}\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}$,
where $\ell^{\pm}=e^{\pm},\mu^{\pm}$. This is the first time that such a
complete study for this process is conducted at the NLO level in QCD. With a
full set of higher-order QCD effects included, we can investigate their effect
at the integrated fiducial cross-section level as well as for various
infrared-safe observables. In addition, such theoretical predictions allow us
to qualitatively increase the level of precision and accuracy for the SM
$pp\to t\bar{t}t\bar{t}$ process, making future comparison with LHC data more
feasible.

The paper is organised as follows. In Section 2 we describe our calculation
and the Helac-Nlo MC program that we use in our studies. Input parameters and
fiducial phase-space cuts for simulating the ATLAS and/or CMS detector
response are given in Section 3. We briefly review the NWA framework in
Section 4, where we also provide the explicit expressions for the NLO QCD
cross section. Numerical results for the integrated fiducial cross sections
for the LHC Run III energy of $\sqrt{s}=13.6$ TeV are presented in Section 5.
In particular, we show the results for two renormalisation and factorisation
scale choices and three different PDF sets. On the other hand, differential
cross-section distributions, which are generated only for the dynamic scale
setting, are shown in Section 6. The theoretical uncertainties of the
integrated and differential fiducial cross sections, that are associated with
neglected higher order terms in the perturbative expansion and with different
parametrisations of the parton distribution functions (PDFs), are also given
in Section 5 and Section 6. In addition, we study there the differences and
similarities between the expanded and unexpanded NLO QCD results and
investigate the magnitude and impact of higher-order QCD effects in top-quark
decays. Finally, in Section 7 our conclusions are given.

##  2 Description of the calculation

In high-energy $pp$ collisions at the LHC, the production of four top quarks
can take place, followed by their decay into various final states. This
process is initiated either by the interaction of two gluons ($gg$) or by a
pair of a quark and an anti-quark ($q\bar{q}$). For the fully leptonic top-
quark decay channel that we are focusing on, it can be described as follows:
$gg/q\bar{q}\to t\bar{t}t\bar{t}\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}$,
where $q=u,d,s,c,b$ and $\ell^{\pm}$ stands for $e^{\pm},\mu^{\pm}$. The
leading order (LO) contributions considered in this paper are of the order of
$\mathcal{O}(\alpha_{s}^{4}\alpha^{8})$ and they consist of Feynman diagrams
containing four intermediate top quarks and $W$ bosons along with their
corresponding decay products. This implies that for the production of four top
quarks, we are only taking into account the dominant QCD contributions, while
the electroweak coupling becomes relevant only during the decay phase. At LO
the $gg$ subprocess involves a total of 72 Feynman diagrams, while for each
$q\bar{q}$ subprocess, the corresponding number is 14. Representative Feynman
diagrams are depicted in Figure 1, where, for simplicity, the decays of top
quarks are not shown.

Figure 1: Representative Feynman diagrams for the process $pp\to
t\bar{t}t\bar{t}\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}$
at the order of $\mathcal{O}(\alpha_{s}^{4}\alpha_{s}^{8})$. For simplicity,
the decays of top quarks are not shown in the diagrams. All graphs were
generated with FeynGame Harlander:2020cyh .

In addition, we calculate next-to-leading-order (NLO) QCD corrections of the
order of $\mathcal{O}(\alpha_{s}^{5}\alpha^{8})$. These higher-order
corrections comprise both virtual and real emission contributions. The virtual
corrections are derived by calculating the interference between the sum of all
one-loop diagrams with the Born amplitude. For the real emission
contributions, an additional QCD parton is emitted, leading to the following
subprocesses:

  * •

$gg\to tt\bar{t}\bar{t}\,(g)\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\,(g)\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}\,(g)\,,$

  * •

$\bar{q}q/q\bar{q}\to tt\bar{t}\bar{t}\,(g)\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\,(g)\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}\,(g)\,,$

  * •

$gq/qg\to tt\bar{t}\bar{t}\,q\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\,q\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}\,q\,,$

  * •

$g\bar{q}/\bar{q}g\to tt\bar{t}\bar{t}\,\bar{q}\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}\,b\bar{b}\,\bar{q}\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}\,b\bar{b}\,\bar{q}\,,$

where in the case of the $gg/q\bar{q}/\bar{q}q$ initiated subprocesses $(g)$
indicates the fact that the additional gluon can be emitted during both the
production and decay stages. Throughout our study, we adopt the NWA for the
treatment of the top quark and the $W$ gauge boson. This approximation is
valid when a width $(\Gamma)$ of an unstable particle is much smaller than its
mass $(m)$, allowing us to safely take the limit $\Gamma/m\to 0$ in the Breit-
Wigner propagator

$\lim_{\Gamma/m\to
0}\dfrac{1}{(p^{2}-m^{2})^{2}+m^{2}\Gamma^{2}}=\dfrac{\pi}{m\Gamma}\delta(p^{2}-m^{2})\,.\\\
\vspace{0.2cm}$ (1)

Indeed, this is the case in our process where $\Gamma_{t}/m_{t}\approx 0.008$
and $\Gamma_{W}/m_{W}\approx 0.026$, justifying our treatment. The NWA ensures
that the considered particles are kept on-shell, which allows us to separately
analyse the production and decay stages of the top quarks and investigate the
impact of higher-order corrections at each of these stages. Hence, in this
study, we examine various scenarios involving NLO QCD corrections applied
either in both the production and decay stages or solely during the production
of four top quarks. We label the former approach as $\sigma^{\rm NLO}_{\rm
full}$ and the latter as $\sigma^{\rm NLO}_{\rm LO_{dec}}$. Additionally, we
are going to assess the effects of expanding the top-quark width in our
calculations, referred to as $\sigma^{\rm NLO}_{\rm exp}$. More details about
these three approaches are going to be provided in Section 4.

Our NLO QCD calculations are performed with the help of the Helac-NLO program
Bevilacqua:2011xh . In this framework the calculation of scattering amplitudes
is based on the Dyson-Schwinger recursive algorithm Draggiotis:1998gr ;
Kanaki:2000ey ; Draggiotis:2002hm ; Cafarella:2007pc . The Helac-NLO program
includes Helac-1loop vanHameren:2009dr  for the evaluation of the numerators
of the loop integrals and the rational terms, CutTools Ossola:2007ax , which
implements the OPP method for the numerical evaluation of one-loop amplitudes,
based on a decomposition at the integrand level Ossola:2006us ; Ossola:2008xq
; Draggiotis:2009yb , and OneLoop vanHameren:2010cp  for the evaluation of the
scalar integrals. The preservation of gauge symmetries by this approach is
explicitly checked by studying Ward identities up to the one-loop level. For
the $gg$ subprocess we perform this test for every phase-space point.
Quadruple precision is used to recompute events which fail the gauge-
invariance check. For the $q\bar{q}$ subprocess, on the other hand, we use the
so-called scale test Badger:2010nx . The latter test, that is based on
momentum rescaling, is also performed for each phase-space point. Also in this
case, for the failed points the amplitude is recomputed using higher
precision. In order to optimise the performance of the Helac-1loop system
reweighting techniques, helicity and colour sampling methods are used. We
employ the Catani-Seymour dipole subtraction formalism Catani:1996vz  to
extract the soft and collinear infrared singularities and to combine them with
the virtual corrections. Specifically, the formulation presented in Ref.
Catani:2002hc  for massive quarks as implemented in the Helac-Dipoles MC
program Czakon:2009ss  is employed. It has been further extended by us to work
for arbitrary helicity eigenstates and colour configurations of the external
partons. Furthermore, Helac-Dipoles comprises another subtraction scheme
Bevilacqua:2013iha  that uses the momentum mapping and the splitting functions
derived in the context of an improved parton shower formulation by Nagy and
Soper Nagy:2007ty . Also in this case random polarisation and colour sampling
of the external partons are employed to optimise the calculations. Compared to
standard dipole subtraction, this scheme features a significantly smaller
number of subtraction terms and facilitates the matching of NLO calculations
with parton showers including quantum interference Czakon:2015cla . For both
subtraction schemes a phase-space restriction on the contribution of the
subtraction terms, called $\alpha_{\rm max}$, is implemented Nagy:1998bb ;
Nagy:2003tz ; Bevilacqua:2009zn ; Czakon:2015cla . After combining virtual and
real corrections, singularities connected to collinear configurations in the
final state as well as soft divergences in the initial and final states cancel
for collinear-safe observables after applying a jet algorithm. Singularities
connected to collinear initial-state splittings are removed via factorisation
by parton distribution function redefinitions. The phase space integration is
performed with the help of the Monte Carlo generator Kaleu vanHameren:2010gg ,
including Parni vanHameren:2007pt  for the importance sampling. For the real
corrections Kaleu is equipped with additional, special channels that proved to
be important for phase-space optimisation in the subtracted real emission part
of the calculation Bevilacqua:2010qb .

We have conducted multiple cross-checks to ensure the accuracy of our results.
To begin with, we checked the coefficients of the poles in $\epsilon$, i.e.
$1/\epsilon$ and $1/\epsilon^{2}$, of the virtual amplitudes for $gg$ and
$q\bar{q}$ subprocesses with MadGraph-aMC@NLO Alwall:2014hca  for a few random
phase-space points for stable $t\bar{t}t\bar{t}$ production and for the $t\to
W^{+}b\to\ell^{+}\nu_{\ell}\,b$ decay. In the next step, we repeated this
procedure using the Recola MC program Actis:2016mpe . This time, however, also
the finite parts along with the coefficients of the poles in $\epsilon$ were
cross-checked without separating the production stage from the top-quark
decays. Because Recola is able to provide matrix elements in the so-called
double-pole approximation Denner:2000bj ; Accomando:2004de ; Denner:2016jyo ,
it was rather straightforward to prepare this program for the use in the NWA
by interfacing it to the Helac-NLO system. This approach has already been
tested before in the case of NLO QCD calculations for $pp\to t\bar{t}jj$
Bevilacqua:2022ozv  and $pp\to t\bar{t}\gamma\gamma$ Stremmer:2023kcd . For
the real corrections, we have checked for several phase-space points that the
poles in $\epsilon$ that appear in the ${\cal I}$-operator and those appearing
in the virtual part of our calculations cancel. We have confirmed that the
real emission part of the calculation is independent of any phase-space
restriction imposed by the unphysical $\alpha_{\rm max}$ parameter. This check
has been applied to both the Catani-Seymour and Nagy-Soper subtraction
schemes. Finally, having two independent subtraction schemes allowed us to
cross-check the correctness of the real emission part of the calculation in an
even more robust way.

All our LO and NLO results are stored in the modified Les Houches Files
Alwall:2006yp , which are subsequently converted into the Root Ntuple Files
Antcheva:2009zz ; Bern:2013zja . This approach offers a significant advantage
as it allows us to efficiently reweight these events. As a consequence, we can
produce differential distributions for new observables, perform their
rebinning or generate theoretical predictions for various scenarios, including
different PDF sets, other scale settings, and various (more exclusive) cuts on
the final state particles, all without the need to rerun the entire process
from scratch. All the above tasks are carried out using the HEPlot program
Bevilacqua:HEPlot .

##  3 Computational setup

We consider the process

$pp\to t\bar{t}t\bar{t}+X\to
W^{+}W^{-}W^{+}W^{-}b\bar{b}b\bar{b}+X\to\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,\,\ell^{+}\nu_{\ell}\,\ell^{-}\bar{\nu}_{\ell}\,b\bar{b}b\bar{b}+X\\\
\vspace{0.2cm}$ (2)

at NLO in QCD, where $\ell^{\pm}$ stands for $\ell^{\pm}=e^{\pm},\mu^{\pm}$.
To perform this calculation we use the NWA. For simplicity, we refer to this
process as $pp\to t\bar{t}t\bar{t}$ in the $4\ell$ top-quark decay channel. We
do not consider $\tau$ leptons as they are often analysed separately at the
LHC due to their rather rich and complex decay pattern. Our results are
obtained for LHC Run III center-of-mass energy of $\sqrt{s}=13.6$ TeV. The
numerical values we use for the SM input parameters are given below:

$\begin{array}[]{lll}G_{\mu}=1.1663787\cdot 10^{-5}~{}{\rm
GeV}^{-2}\,,&&m_{t}=172.5~{}{\rm GeV}\,,\vspace{0.2cm}\\\ m_{W}=80.379~{}{\rm
GeV}\,,&&\Gamma_{W}^{\rm NLO}=2.0972~{}{\rm GeV}\,,\vspace{0.2cm}\\\
m_{Z}=91.1876~{}{\rm GeV}\,,&&m_{b}=0~{}{\rm GeV}\,.\end{array}$ (3)

Apart from the top quark and the $W$ boson, all the other particles appearing
in our process are treated as massless. Also, we keep the Cabibbo-Kobayashi-
Maskawa mixing matrix diagonal for all our calculations. Since the leptonic
$W$-boson decays do not receive NLO QCD corrections, to account for some
higher-order effects in their decays, we use the NLO QCD value for the $W$
gauge boson width as calculated for $\mu_{R}=m_{W}$. We utilise the
$\Gamma_{W}^{\rm NLO}$ value everywhere, i.e. in the calculation of LO and NLO
matrix elements. Regarding the width of the top quark, we adopt its LO (NLO)
value in LO (NLO) calculations. The LO and NLO values of $\Gamma_{t}$ are
determined using the formulas provided in Ref. Denner:2012yc , where
$\alpha_{s}(\mu_{R}=m_{t})$ is used:

$\begin{array}[]{lll}\Gamma_{t}^{\rm LO}=1.4806842~{}{\rm
GeV}\,,&&\Gamma_{t}^{\rm NLO}=1.3535983~{}{\rm GeV}\,.\end{array}$ (4)

The electromagnetic coupling $\alpha$ is evaluated in the $G_{\mu}$-scheme as
follows:

$\alpha_{G_{\mu}}=\frac{\sqrt{2}}{\pi}\,G_{F}\,m_{W}^{2}\,\sin^{2}\theta_{W}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\sin^{2}\theta_{W}=1-\frac{m_{W}^{2}}{m_{Z}^{2}}\,.\\\
\vspace{0.2cm}$ (5)

In accordance with the guidance for Standard Model processes provided by the
PDF4LHC working group Butterworth:2015oua , we utilise three individual,
modern PDF sets: MSHT20 Bailey:2020ooq  (our default PDF set), NNPDF3.1
NNPDF:2017mvq  and CT18 Hou:2019efy . Since there is no LO PDF set available
for CT18, we have substituted it with CT14. In particular, in this case we
employ CT14llo with $\alpha_{s}(m_{Z})=0.130$ for our LO calculations
Dulat:2015mca . The running of the strong coupling is performed with a two-
loop (one-loop) accuracy at NLO (at LO) and provided by the LHAPDF library
Buckley:2014ana  involving five active flavours.

The limitations of fixed-order calculations lead to an inherent dependence on
the renormalisation $(\mu_{R})$ and factorisation scale $(\mu_{F})$. The
theoretical uncertainty of the cross section, associated with neglected
higher-order terms in the perturbative expansion, can be estimated by varying
$\mu_{R}$ and $\mu_{F}$ in $\alpha_{s}$ and PDFs, up and down by a factor $2$
around the central value of the scale of the process $(\mu_{0})$ in the range

$\frac{1}{2}\leq\frac{\mu_{R}}{\mu_{0}},\frac{\mu_{F}}{\mu_{0}}\leq 2\,.\\\
\vspace{0.2cm}$ (6)

Traditionally, the following additional condition is required

$\frac{1}{2}\leq\frac{\mu_{R}}{\mu_{F}}\leq 2\,.\\\ \vspace{0.2cm}$ (7)

We are looking for the minimum and maximum of the resulting cross sections,
taking into account the following pairs

$\left(\frac{\mu_{R}}{\mu_{0}}\,,\frac{\mu_{F}}{\mu_{0}}\right)=\Big{\\{}\left(2,1\right),\left(0.5,1\right),\left(1,2\right),(1,1),(1,0.5),(2,2),(0.5,0.5)\Big{\\}}\,.\\\
\vspace{0.2cm}$ (8)

At this point, it is important to note that we use the same numerical value
for the top-quark width across all scale choices during the scale variation.
We have observed though that the uncertainties stemming from varying the
$\mu_{R}$ scale in the calculation of $\Gamma_{t}$ are of the order of $4\%$
only. As we will discuss in Section 5, these uncertainties are similar in
magnitude to the internal PDF uncertainties. In our studies, we make use of
two scale settings. For the fixed scale choice we use

$\mu_{0}=2m_{t}\,,\\\ \vspace{0.2cm}$ (9)

which has already been employed in our previous studies for the $pp\to
t\bar{t}t\bar{t}$ process with stable top quarks Bevilacqua:2012em . On the
other hand, for the dynamical scale setting we utilise

$\mu_{0}=E_{T}/4\,,\\\ \vspace{0.2cm}$ (10)

with $E_{T}$ defined according to

$E_{T}=\sum_{i=1,2}m_{T}(t_{i})+\sum_{i=1,2}m_{T}(\,\bar{t}_{i}\,)\,,\\\
\vspace{0.2cm}$ (11)

where the transverse masses $m_{T}(t)$ and $m_{T}(\,\bar{t}\,)$ are given by

$m_{T}(t)=\sqrt{m_{t}^{2}+p_{T}^{2}(t)}\,,\\\
\vspace{0.2cm}\quad\quad\quad\quad\quad
m_{T}(\,\bar{t}\,)=\sqrt{m_{t}^{2}+p_{T}^{2}(\,\bar{t}\,)}\,.\\\
\vspace{0.2cm}$ (12)

This particular scale choice has already been explored in other studies for
this process, see e.g. Refs. Bevilacqua:2012em ; Jezo:2021smh ;
Frederix:2017wme . We note that in our approach, we construct the momenta of
the top quarks from their decay products without accounting for the potential
presence of an additional light jet, even if it arises in top-quark decays. We
have checked, however, that its inclusion in the construction of top quark
momenta does not have any significant impact on our results. Jets are
constructed from all final-state partons with $|\eta|<5$ with the help of the
infrared-safe anti-$k_{T}$ jet algorithm Cacciari:2008gp  with a clustering
parameter $R=0.4$. We require at least four $b$-jets as well as four charged
leptons. All final states have to fulfil the following selection criteria that
closely mimic the ATLAS and CMS detector response ATLAS:2023ajo ; CMS:2023ftu
:

$\begin{array}[]{lll}p_{T,\,\ell}>25~{}{\rm
GeV}\,,&\quad\quad\quad\quad\quad|y_{\ell}|<2.5\,,&\quad\quad\quad\quad\quad\Delta
R_{\ell\ell}>0.4\,,\\\\[5.69046pt] p_{T,\,b}>25~{}{\rm
GeV}\,,&\quad\quad\quad\quad\quad|y_{b}|<2.5\,,&\quad\quad\quad\quad\quad\Delta
R_{bb}>0.4\,.\end{array}$ (13)

We set no restriction on the kinematics of the extra light jet (if resolved)
and employ no cut on the missing transverse momentum.

##  4 Explicit expressions for the NWA cross section through NLO

In the following, we briefly review some of the features of the NWA that will
be exploited throughout our study. The fully factorised form of the
differential cross section for the $pp\to t\bar{t}t\bar{t}$ process including
top-quark decays in the NWA can be given by, see e.g. Melnikov:2009dn ;
Campbell:2012uf ; Bevilacqua:2019quz ; Behring:2019iiv ; Czakon:2020qbd ;
Bevilacqua:2022twl

$d\sigma=d\sigma_{t\bar{t}t\bar{t}}\times\frac{d\Gamma_{t}}{\Gamma_{t}}\times\frac{d\Gamma_{\bar{t}}}{\Gamma_{t}}\times\frac{d\Gamma_{t}}{\Gamma_{t}}\times\frac{d\Gamma_{\bar{t}}}{\Gamma_{t}}\,,\\\
\vspace{0.2cm}$ (14)

where $d\sigma_{tt\bar{t}\bar{t}}$ denotes the cross section for on-shell
$pp\to t\bar{t}t\bar{t}$ production and $d\Gamma_{t}$ ($d\Gamma_{\bar{t}}$) is
the top-quark decay rate for $t\to W^{+}b\to\ell^{+}\nu_{\ell}b$ ($\bar{t}\to
W^{-}\bar{b}\to\ell^{-}\bar{\nu}_{\ell}\bar{b}$). Furthermore, the symbol
$\times$ represents spin correlations. This formula holds to all orders in
$\alpha_{s}$. However, all the terms on the right-hand side of Eq. (14) can be
expanded in terms of $\alpha_{s}$ to the required order in perturbative
expansion. Thus, for example at the NLO level we have

$\begin{split}d\sigma^{\rm
NLO}&=d\sigma_{tt\bar{t}\bar{t}}^{(0)}+\alpha_{s}d\sigma_{tt\bar{t}\bar{t}}^{(1)}\,,\\\\[5.69046pt]
d\Gamma_{t}^{\rm
NLO}&=d\Gamma_{t}^{(0)}+\alpha_{s}d\Gamma_{t}^{(1)}\,,\end{split}$ (15)

where $d\sigma^{(0)},\,d\Gamma^{(0)}$ and $d\sigma^{(1)},\,d\Gamma^{(1)}$
denote LO and NLO contributions to the production and decay processes,
respectively. In addition, in Eq. (14) we have the top-quark width,
$\Gamma_{t}$, which enters the NWA calculation as a parameter through the
denominator. There are different ways to treat the $1/\Gamma_{t}$ factor. For
example, it is possible to set the value of $\Gamma_{t}$ to a numerical value
corresponding to the perturbative order of the full calculation,
$\Gamma_{t}^{\rm NLO}$. Alternatively, one can formally expand $\Gamma_{t}$ in
powers of $\alpha_{s}$ and keep only the terms consistent with the order of
the performed calculation. Thus, for example at NLO in QCD we would have

$\Gamma_{t}^{\rm NLO}=\Gamma_{t}^{(0)}+\alpha_{s}\Gamma_{t}^{(1)}\,.\\\
\vspace{0.2cm}$ (16)

By incorporating Eq. (15) into Eq. (14) we obtain the following expression for
$d\sigma^{\textrm{LO}}$ at the LO level

$d\sigma^{\textrm{LO}}=~{}d\sigma^{\rm
LO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
LO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm
LO}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
LO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm LO}}\,,\\\
\vspace{0.2cm}$ (17)

whereas, at the NLO QCD level we would have instead
$d\sigma^{\textrm{NLO}}_{\rm full}$ given by

$\begin{split}d\sigma^{\textrm{NLO}}_{\rm full}&~{}=d\sigma^{\rm
NLO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm NLO}}\\\\[5.69046pt]
&~{}+d\sigma^{\rm
LO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{1}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm NLO}}\\\\[5.69046pt]
&~{}+d\sigma^{\rm
LO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{1}_{\bar{t}}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm NLO}}\\\\[5.69046pt]
&~{}+d\sigma^{\rm
LO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{1}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm NLO}}\\\\[5.69046pt]
&~{}+d\sigma^{\rm
LO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
NLO}}\times\frac{d\Gamma^{1}_{\bar{t}}}{\Gamma_{t}^{\rm NLO}}\,.\end{split}$
(18)

The latter formula describes the full NLO QCD calculation with higher order
QCD corrections incorporated in both $t\bar{t}t\bar{t}$ production and top-
quark decays. It is worth pointing out that the various terms in Eq. (18) are
separately infrared finite and do not interfere with each other. There is no
cross-talk between the production stage and decays, as well as between
different decays of the top quarks. The presence of $\Gamma_{t}^{\rm NLO}$ in
Eq. (18), which is a function of the strong coupling constant, spoils the
rigorous expansion of the cross section in powers of $\alpha_{s}$. If also Eq.
(16) is incorporated into Eq. (14) we would rather have the following
expression

$d\sigma^{\rm NLO}_{\textrm{exp}}=d\sigma^{\rm
NLO}=d\sigma^{\textrm{NLO}}_{\rm
full}\times\left(\dfrac{\Gamma^{\textrm{NLO}}_{t}}{\Gamma^{\textrm{LO}}_{t}}\right)^{4}-d\sigma^{\textrm{LO}}\times\dfrac{4(\Gamma^{\textrm{NLO}}_{t}-\Gamma^{\textrm{LO}}_{t})}{\Gamma^{\textrm{LO}}_{t}}\,,\\\
\vspace{0.2cm}$ (19)

that represents the corresponding expanded NLO QCD result, labelled as
$d\sigma^{\rm NLO}_{\textrm{exp}}$. We note that in Eq. (19) the
$d\sigma^{\textrm{LO}}$ part is calculated with the NLO PDF set, thus, with
the two-loop running of $\alpha_{s}$. Moreover, having all the components of
the NLO QCD calculation in the NWA at our disposal, we can also study the case
where the NLO QCD corrections in top-quark decays are simply neglected. This
last case we label as $d\sigma^{\rm NLO}_{\rm LO_{dec}}$. The contribution to
this result comes from the first line of Eq. (18) only, while $\Gamma_{t}^{\rm
LO}$ is used in the calculation. Thus, we can write for
$d\sigma^{\textrm{NLO}}_{\rm LO_{dec}}$

$d\sigma^{\textrm{NLO}}_{\rm LO_{dec}}=~{}d\sigma^{\rm
NLO}_{t\bar{t}t\bar{t}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
LO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm
LO}}\times\frac{d\Gamma^{0}_{t}}{\Gamma_{t}^{\rm
LO}}\times\frac{d\Gamma^{0}_{\bar{t}}}{\Gamma_{t}^{\rm LO}}\,.\\\
\vspace{0.2cm}$ (20)

In this work we consider predictions in the NWA based on the various levels of
accuracy, as described above. Specifically, we are going to provide the
results for $\sigma^{\textrm{NLO}}_{\textrm{exp}}$ (our default findings at
NLO QCD) and $\sigma^{\textrm{NLO}}_{\rm LO_{dec}}$ at the integrated and
differential cross-section level. Furthermore, we show the results for the
$\sigma^{\rm NLO}_{\rm full}$ case. The goal is to analyse the similarities
and differences between the three approaches and to estimate the size of
higher-order corrections in top-quark decays. As mentioned earlier, in the
literature these latter effects are approximated only by parton shower
programs.

##  5 Integrated fiducial cross sections

We begin the presentation of our results for the $pp\to t\bar{t}t\bar{t}$
process in the $4\ell$ top-quark decay channel with a discussion of the
integrated fiducial cross section for the LHC Run III energy of
$\sqrt{s}=13.6$ TeV. When considering the input parameters and the phase-space
cuts as detailed in Section 3, our LO and NLO findings as calculated in the
NWA for the fixed scale setting, $\mu_{R}=\mu_{F}=\mu_{0}=2m_{t}$, and the
default MSHT20 PDF set, are given by:

$\displaystyle\begin{split}\sigma^{\textrm{LO}}(\textrm{MSHT20},\mu_{0}&=2m_{t})=4.3868(3)^{+73\%}_{-40\%}\textrm{
[scale] ab}\,,\\\
\sigma^{\textrm{NLO}}(\textrm{MSHT20},\mu_{0}&=2m_{t})=4.885(2)^{+13\%}_{-20\%}\textrm{
[scale] }^{+4\%}_{-3\%}\textrm{ [PDF] ab}\,.\end{split}$ (21)

We use the expanded version of the NWA at NLO in QCD, which for the sake of
brevity is denoted simply by $\sigma^{\rm NLO}$. On the other hand, when the
phase-space dependent scale choice is employed,
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$, we obtain:

$\displaystyle\begin{split}\sigma^{\textrm{LO}}(\textrm{MSHT20},\mu_{0}&=E_{T}/4)=4.7479(3)^{+74\%}_{-40\%}\textrm{
[scale] ab}\,,\\\
\sigma^{\textrm{NLO}}(\textrm{MSHT20},\mu_{0}&=E_{T}/4)=5.158(3)^{+12\%}_{-20\%}\textrm{
[scale] }^{+4\%}_{-3\%}\textrm{ [PDF] ab}\,.\end{split}$ (22)

At the LO level, the dominant contribution comes from the $gg$ channel,
accounting for approximately $88\%$ of the full integrated fiducial cross
section, while the remaining contributions that come from the
$q\bar{q}/\bar{q}q$ channels collectively make up about $12\%$. The NLO QCD
corrections for this process are rather moderate, in the range of $9\%-11\%$,
depending on the scale setting. They are much smaller than the LO scale
uncertainties, which are of the order of $70\%$ for both scale choices. The
incorporation of NLO calculations significantly reduces the latter
uncertainties down to $20\%$. We note here, that for both scale settings and
independently of the order in perturbative expansion these uncertainties are
asymmetric. Apart from the theoretical error arising from the scale
dependence, we must also account for another source of uncertainty linked to
the parameterisation of PDFs. We use the prescription from the MSHT20 group to
provide the $68\%$ confidence level (C.L.) PDF uncertainties. These internal
PDF uncertainties are of the order of $4\%$ for both scale settings. Thus,
they are well below the theoretical uncertainties due to the scale dependence.
The latter uncertainties remain the dominant source of the systematics.

PDF |  $\sigma^{\textrm{LO}}$ [ab] | $\delta_{scale}$ |  $\sigma^{\textrm{NLO}}$ [ab] | $\delta_{scale}$ | $\delta_{PDF}$ | $\mathcal{K}=\sigma^{\textrm{NLO}}/\sigma^{\textrm{LO}}$
---|---|---|---|---|---|---
$\mu_{R}=\mu_{F}=\mu_{0}=2m_{t}$
MSHT20 | 4.3868(3) |  |  $+3.2237$ (73%)
---
$-1.7332$ (40%)
4.885(2) |  |  $+0.618$ (13%)
---
$-0.998$ (20%)
|  $+0.210$ (4%)
---
$-0.155$ (3%)
1.11
NNPDF3.1 | 3.7389(2) |  |  $+2.6811$ (72%)
---
$-1.4545$ (39%)
4.837(2) |  |  $+0.625$ (13%)
---
$-0.999$ (21%)
|  $+0.105$ (2%)
---
$-0.105$ (2%)
1.29
CT18 | 4.6757(3) |  |  $+3.3754$ (72%)
---
$-1.8311$ (39%)
4.847(2) |  |  $+0.613$ (13%)
---
$-0.989$ (20%)
|  $+0.288$ (6%)
---
$-0.236$ (5%)
1.04
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$
MSHT20 | 4.7479(3) |  |  $+3.5156$ (74%)
---
$-1.8855$ (40%)
5.158(3) |  |  $+0.630$ (12%)
---
$-1.051$ (20%)
|  $+0.218$ (4%)
---
$-0.161$ (3%)
1.09
NNPDF3.1 | 4.0930(3) |  |  $+2.9792$ (73%)
---
$-1.6063$ (39%)
5.114(3) |  |  $+0.627$ (12%)
---
$-1.050$ (21%)
|  $+0.110$ (2%)
---
$-0.110$ (2%)
1.25
CT18 | 5.0003(3) |  |  $+3.6151$ (72%)
---
$-1.9623$ (39%)
5.115(3) |  |  $+0.627$ (12%)
---
$-1.041$ (20%)
|  $+0.297$ (6%)
---
$-0.245$ (5%)
1.02

Table 1: Integrated fiducial cross sections at LO and NLO in QCD for the
$pp\to t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with
$\sqrt{s}=$ 13.6 TeV. Results are presented for MSHT20, NNPDF3.1 and CT18 PDF
sets. They are evaluated using $\mu_{0}=2m_{t}$ and $\mu_{0}=E_{T}/4$. Also
given are theoretical uncertainties coming from the scale variation and from
PDFs. In the last column the ${\cal K}$-factor is shown. Figure 2: Integrated
fiducial cross sections at NLO in QCD for the $pp\to t\bar{t}t\bar{t}$ process
in the $4\ell$ channel at the LHC with $\sqrt{s}=$ 13.6 TeV. Scale and PDF
uncertainties are depicted for the dynamical scale setting
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$. Results are evaluated using the following
NLO PDF sets: MSHT20, CT18 and NNPDF3.1. The red vertical line represents the
central prediction for our default MSHT20 PDF set.

Figure 3: Scale dependence of the integrated fiducial cross section at LO and
NLO for the $pp\to t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC
with $\sqrt{s}=$ 13.6 TeV. The LO and the NLO MSHT20 PDF sets are employed.
Renormalisation and factorisation scales are set to the common value
$\mu_{R}=\mu_{F}=\mu_{0}$ where $\mu_{0}=2m_{t}$ and $\mu_{0}=E_{T}/4$. Both
scales are varied simultaneously by a factor of $\xi$ in the following range
$\xi\in\left\\{0.125,\dots,8\right\\}$. For each case of $\mu_{0}$ also shown
is the variation of $\mu_{R}$ with fixed $\mu_{F}$ and the variation of
$\mu_{F}$ with fixed $\mu_{R}$.

Figure 4: Contour plots for the scale dependence in the ($\mu_{R},\mu_{F}$)
plane at LO (upper plots) and NLO (lower plots) for our default PDF set MSHT20
and two scale choices $\mu_{0}=2m_{t}\;\textrm{and}\;\mu_{0}=E_{T}/4$. Results
are presented for the $pp\to t\bar{t}t\bar{t}$ process in the $4\ell$ channel
at the LHC with $\sqrt{s}=$ 13.6 TeV.

Within a given PDF set we can not assess the differences that enter into the
parameterisation of the various PDF sets. To investigate this issue, following
the recommendation of the PDF4LHC group for the usage of PDFs suitable for
applications at the LHC that are connected to SM processes, we additionally
present our results for NNPDF3.1 and CT18 PDF sets. The LO and NLO integrated
fiducial cross sections for the three PDF sets are shown in Table 1 for two
different scale choices: $\mu_{0}=2m_{t}$ and $\mu_{0}=E_{T}/4$. The
$\mathcal{K}$ factor, defined as ${\cal K}=\sigma^{\rm NLO}/\sigma^{\rm LO}$,
is also provided in the last column. The $68\%$ C.L. PDF uncertainties for
NNPDF3.1 are up to $2\%$ only, whereas for the CT18 PDF set they are of the
order of $6\%$. As expected, we observe significant variations in the values
of the LO cross section for the three LO PDF sets that are provided for
$\alpha_{s}(m_{Z})=0.130$. These differences result in a variety of the
$\mathcal{K}$ factor: ${\cal K}=(1.25-1.29)$ for NNPDF3.1 and ${\cal
K}=(1.02-1.04)$. Using NLO PDF sets instead for the LO calculations would
correspond to lower values of $\sigma^{\rm LO}$ and, consequently, higher
${\cal K}$ factors. On the other hand, at the NLO QCD level and independent of
the scale setting, the differences among $\sigma^{\rm NLO}$ for the three PDF
sets are of the order of $1\%$ only.

The graphical representation of the NLO results for the three PDF sets:
MSHT20, CT18 and NNPDF3.1 for the $pp\to t\bar{t}t\bar{t}$ process in the
$4\ell$ top-quark decay channel is provided in Figure 2. Since the results at
the integrated fiducial cross-section level are similar for $\mu_{0}=2m_{t}$
and $\mu_{0}=E_{T}/4$, for simplicity we display the results only for the
dynamic scale choice. For comparison, we also show the central prediction for
our default MSHT20 PDF set, which is represented by a red vertical line.
Figure 2 shows a clear representation of the theoretical-error budget for our
NLO QCD predictions for this process. We can observe that even though the
scale uncertainties are asymmetric for the three PDF sets, they remain
consistent in magnitude and are of the order of $20\%$. In each case, they
also dominate the size of the final theoretical error.

The graphical representation of the result for the scale dependence for both
scale settings: $\mu_{0}=2m_{t}$ and $\mu_{0}=E_{T}/4$ is shown in Figure 3.
The behaviour of LO and NLO cross sections for the default MSHT20 PDF sets is
shown while varying the renormalisation and factorisation scales
simultaneously $\mu_{R}=\mu_{F}=\xi\mu_{0}$ by a factor of $\xi$ in the range
$\xi\in\left\\{0.125,...,8\right\\}$. For each case of $\mu_{0}$ also shown is
the variation of $\mu_{R}$ with fixed $\mu_{F}$ and the variation of $\mu_{F}$
with fixed $\mu_{R}$. As already discussed, the LO scale dependence is large
illustrating the well known fact that the LO prediction can only provide a
rough estimate of the cross section for the process at hand. A significant
reduction in the scale uncertainty can be seen when NLO QCD corrections are
incorporated. We can further observe that the primary source of NLO scale
uncertainty is related to the variation of $\mu_{R}$ (purple line). Changes in
the factorisation scale $\mu_{F}$ (green line) have a relatively small impact
on the integrated fiducial cross section. Lastly, we remark that scale choices
with $\xi<0.25$ might yield negative values for the integrated NLO fiducial
cross section and thus should be avoided since they correspond to unphysical
results.

The dependence of the LO and NLO cross sections on $\mu_{R}$ and $\mu_{F}$,
which are varied independently around a central value of the scale, is
presented in Figure 4. Thus, this time we plot the distribution of the LO and
NLO cross sections in the $\mu_{R}-\mu_{F}$ plane. In addition to the previous
three special cases shown in Figure 3: $\mu_{R}=\mu_{F}=\xi\mu_{0}$,
$\mu_{R}=\mu_{0},\,\mu_{F}=\xi\mu_{0}$ and
$\mu_{R}=\xi\mu_{0},\,\mu_{F}=\mu_{0}$, also other cases are depicted here.
These contour plots provide complementary information to the previous scale
dependence plots. We display the ratio of the integrated fiducial cross
section $\sigma^{\rm(N)LO}(\mu_{R,\,F}=\xi\mu_{0})$ to the one calculated
using the central value of the scale $\sigma^{\rm(N)LO}(\mu_{R,\,F}=\mu_{0})$.
It should be noted that the region with $\mu_{R}<\mu_{0}/4$ is not included in
the NLO plots due to the negative values of the integrated cross section, as
already explained above. Analysing the entire plotted range of $\xi$, as we
move toward lower values of $\mu_{R}$ and $\mu_{F}$, we can see a significant
increase, exceeding a factor of $6$, in the $\sigma^{\rm
LO}(\mu_{R,\,F}=\xi\mu_{0})/\sigma^{\rm LO}(\mu_{R,\,F}=\mu_{0})$ ratio. On
the other hand, larger scale values yield smaller cross section ratios. At
NLO, the scale variations are less pronounced, reaching an upper limit of
$1.22-1.24$ and a lower limit of $0.46$ for the $\sigma^{\rm
NLO}(\mu_{R,\,F}=\xi\mu_{0})/\sigma^{\rm NLO}(\mu_{R,\,F}=\mu_{0})$ ratio for
both scale settings. Thus, when the whole range of $\xi$ is considered the
uncertainties are substantial, even at the NLO level. Again we can observe
that the scale variation is driven by the changes in $\mu_{R}$.

We can conclude this part by saying that scanning all potential
$(\mu_{R},\mu_{F})$ pairs in the range $0.5\leq\mu_{R}/\mu_{0}\leq 2$ and
$0.5\leq\mu_{F}/\mu_{0}\leq 2$, without applying the additional requirement of
$0.5\leq\mu_{R}/\mu_{F}\leq 2$, would not change the size of the NLO scale
uncertainties. Indeed, for the two additional extreme cases
$(\mu_{R}/\mu_{0},\mu_{F}/\mu_{0})=\left\\{(2,0.5),(0.5,2)\right\\}$, that are
usually omitted, the magnitude of the scale uncertainties is in the range of
$8\%-10\%$ only, independently of the scale setting. Thus, the size of the NLO
QCD theoretical error for the $pp\to t\bar{t}t\bar{t}$ process in the $4\ell$
channel would remain the same even if these two cases were included.

Decay treatment |  $\sigma_{i}^{\rm NLO}$ [ab] |  $+\delta_{scale}$ [ab] |  $-\delta_{scale}$ [ab] | $\sigma_{i}^{\rm NLO}/\sigma^{\rm NLO}_{\rm exp}-1$
---|---|---|---|---
$\mu_{R}=\mu_{F}=\mu_{0}=2m_{t}$
${\rm full}$ | 5.448(3) |  $+0.147$ (3%) |  $-0.848$ (16%) | $+11.5\%$
${\rm LO_{dec}}$ | 5.280(3) |  $+1.114$ (21%) |  $-1.219$ (23%) | $+8.1\%$
${\rm exp}$ | 4.885(2) |  $+0.618$ (13%) |  $-0.998$ (20%) | $-$
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$
${\rm full}$ | 5.718(3) | $+0.138\;(2\%)$ | $-0.875\;(15\%)$ | $+10.9\%$
${\rm LO_{dec}}$ | 5.628(3) | $+1.213\;(22\%)$ | $-1.310\;(23\%)$ | $+9.1\%$
${\rm exp}$ | 5.158(3) | $+0.630\;(12\%)$ | $-1.051\;(20\%)$ | $-$

Table 2: Integrated fiducial cross sections at NLO in QCD for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with $\sqrt{s}=$
13.6 TeV for three different scenarios: $\sigma^{\rm NLO}_{\rm full}$,
$\sigma^{\rm NLO}_{\rm LO_{dec}}$ and $\sigma^{\rm NLO}_{\rm exp}$. All
results are given for our default setup with the MSHT20 PDF set and for two
different scale choices: $\mu_{0}=2m_{t}$ and $\mu_{0}=E_{T}/4$. Also given
are theoretical uncertainties coming from the scale variation. In the last
column the ratio to the $\sigma^{\rm NLO}_{\rm exp}$ result is provided.

In the next step, we compare our NLO QCD predictions at the integrated
fiducial cross-section level, labelled as $\sigma^{\rm NLO}_{\rm exp}$, to
$\sigma^{\rm NLO}_{\rm LO_{dec}}$ and $\sigma^{\rm NLO}_{\rm full}$. As
described earlier, $\sigma^{\rm NLO}_{\rm full}$ corresponds to the full NWA
calculation that takes into account QCD corrections in $t\bar{t}t\bar{t}$
production and top-quark decays, while $\Gamma_{t}^{\rm NLO}$ is treated as a
fixed parameter. On the other hand, the result for $\sigma^{\rm NLO}_{\rm
LO_{dec}}$ utilises NLO QCD corrections in the production stage only. In this
case, top-quark decays are described with the LO accuracy using
$\Gamma_{t}^{\rm LO}$. The NLO QCD results for the three scenarios are
summarised in Table 2. A few comments are in order. Firstly, we can notice
that the difference between the result in the full NWA and the expanded one is
in the range of $11\%-12\%$, thus, well within the corresponding NLO
uncertainties, which are at the $15\%-16\%$ and $20\%$ level, for $\sigma^{\rm
NLO}_{\rm full}$ and $\sigma^{\rm NLO}_{\rm exp}$ respectively. Secondly, the
expansion of $\Gamma_{t}^{\rm NLO}$ in terms of $\alpha_{s}$ results in the
increased size of the scale dependence. Furthermore, we can compare the
$\sigma^{\rm NLO}_{\rm LO_{dec}}$ result with the $\sigma^{\rm NLO}_{\rm exp}$
one. The omission of NLO QCD corrections in top-quark decays increases not
only the value of the integrated fiducial cross section by $8\%-9\%$, but also
the size of the theoretical uncertainties due to the scale dependence from
$20\%$ to $23\%$. Had we compared $\sigma^{\rm NLO}_{\rm LO_{dec}}$ with the
full $\sigma^{\rm NLO}_{\rm full}$ result rather than with the expanded one,
the difference in the absolute value of the integrated fiducial cross section
would be much smaller. Indeed, since various values for the top-quark width
are used in both results, the difference is at the level of $2\%-3\%$ only
depending on the scale setting. We note that NLO QCD corrections to the top-
quark width are negative and of the order of $9\%$. The $\alpha_{s}$ expansion
of $\Gamma_{t}$ through NLO in QCD with our input parameters and for
$\alpha_{s}(\mu_{R}=m_{t})=0.1076705$ has the following form

$\Gamma_{t}=(1.4806842-1.1803223\,\alpha_{s})\,{\rm GeV}\,.$ (23)

On the other hand, the change in the size of the scale dependence would be
more substantial as we would observe the increase from $15\%-16\%$ to $23\%$.

We conclude this section by stating that the higher-order QCD effects in top-
quark decays for the $pp\to t\bar{t}t\bar{t}$ process in the $4\ell$ channel
are substantial, and amount to half of the theoretical uncertainties due to
scale variation. Moreover, their inclusion reduces the overall size of
theoretical uncertainties. Finally, as expected, the NLO QCD result calculated
with the $1/\Gamma_{t}$ term not expanded has the smallest theoretical
uncertainties, which are in the range of $15\%-16\%$.

##  6 Differential fiducial cross-section distributions

Figure 5: Differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with
$\sqrt{s}=13.6$ TeV. The transverse momenta of the four $b$-jets, that are
ordered according to their $p_{T}$, are displayed for the dynamical scale
setting $\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$ and the (N)LO MSHT20 PDF set. The
blue (red) curve corresponds to the LO (NLO) result. Also shown are the
corresponding uncertainty bands resulting from scale variations. The lower
panels present the differential ${\cal K}$-factor together with its
uncertainty band and the relative scale uncertainties of the LO cross section.
Monte Carlo integration errors are displayed in both panels.

Figure 6: Differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with
$\sqrt{s}=13.6$ TeV. The transverse momenta of the charged leptons, that are
ordered according to their $p_{T}$, are displayed for the dynamical scale
setting $\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$ and the (N)LO MSHT20 PDF set. The
blue (red) curve corresponds to the LO (NLO) result. Also shown are the
corresponding uncertainty bands resulting from scale variations. The lower
panels display the differential ${\cal K}$-factor together with its
uncertainty band and the relative scale uncertainties of the LO cross section.
Monte Carlo integration errors are displayed in both panels.

Figure 7: Differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with
$\sqrt{s}=13.6$ TeV. The $\Delta R_{bb}$ separations between $b$-jets are
displayed for the dynamical scale setting $\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$
and the (N)LO MSHT20 PDF set. The blue (red) curve corresponds to the LO (NLO)
result. Also shown are the corresponding uncertainty bands resulting from
scale variations. The lower panels display the differential ${\cal K}$-factor
together with its uncertainty band and the relative scale uncertainties of the
LO cross section. Monte Carlo integration errors are displayed in both panels.

Figure 8: Differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with $\sqrt{s}=$
13.6 TeV. The transverse momentum of the $b_{1}b_{2}$ system and the system
comprising all $b$-jets are displayed for the dynamical scale setting
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$ and the (N)LO MSHT20 PDF set. The blue (red)
curve corresponds to the LO (NLO) result. Also shown are the corresponding
uncertainty bands resulting from scale variations. The lower panels display
the differential ${\cal K}$-factor together with its uncertainty band and the
relative scale uncertainties of the LO cross section. Monte Carlo integration
errors are displayed in both panels.

Figure 9: Normalised differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with $\sqrt{s}=$
13.6 TeV. The transverse momentum of the $b_{1}b_{2}$ system and the system
comprising all $b$-jets are displayed for the dynamical scale setting
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$ and the (N)LO MSHT20 PDF set. The blue (red)
curve corresponds to the LO (NLO) result. Also shown are the NLO results with
the $p_{T}^{veto}$ cut of $p_{T}^{veto}=(25,50,100,200,300,500)$ GeV. The
lower panels present the ratio to the normalised LO cross section.

While interesting, the NLO QCD corrections to the integrated fiducial cross
section do not give us a complete picture of the impact of the higher-order
effects on the $pp\to t\bar{t}t\bar{t}$ process in the $4\ell$ channel. To
understand how NLO QCD corrections affect specific observables and/or phase-
space regions, it is necessary to analyse the results at the differential
fiducial cross-section level. To this end, we start with the comparison
between LO and NLO QCD computations for several observables. Additionally, we
can inspect the size of theoretical uncertainties in different parts of the
phase space. We present our findings only for the dynamical scale choice
$\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$ and for the default (N)LO MSHT20 PDF set. It
is a well known fact that the fixed scale choice is not adequate at the
differential cross-section level. For example, for some values of
$\mu_{R}=\mu_{F}=\xi\mu_{0}$, where $\mu_{0}=2m_{t}$, the NLO results in the
high-energy tails of dimensionful distributions become negative. Furthermore,
in similar phase-space regions the LO and NLO scale variation bands no longer
overlap. Finally, it can even happen that the size of the scale variation at
the NLO level actually exceeds the one of the LO prediction. All these effects
have already been observed and discussed in detail in our previous studies for
the $pp\to t\bar{t}+X$ process, where
$X=j,\gamma,Z,W^{\pm},H,jj,\gamma\gamma,b\bar{b}$, see e.g. Refs.
Bevilacqua:2015qha ; Bevilacqua:2016jfk ; Bevilacqua:2018woc ;
Bevilacqua:2019cvp ; Bevilacqua:2020pzy ; Bevilacqua:2021cit ;
Stremmer:2021bnk ; Bevilacqua:2022ozv ; Bevilacqua:2022nrm ; Stremmer:2023kcd
. Similar effects with a fixed scale setting have been examined for the $pp\to
t\bar{t}$ process at NLO in QCD in e.g. Ref. Denner:2012yc .

In Figure 5 \- Figure 8 the upper panels show the absolute LO (blue curve) and
NLO (red curve) predictions together with their corresponding scale
uncertainty bands. The lower panels display the differential ${\cal K}$-factor
together with their scale uncertainty bands as well as the relative scale
uncertainties of the LO cross sections. Specifically, the results presented in
the lower panels are normalised to the corresponding LO cross sections
according to

${\cal K}^{\rm NLO}(\xi\mu_{0})=\frac{d\sigma^{\rm
NLO}(\xi\mu_{0})dX}{d\sigma^{\rm LO}(\mu_{0})/dX}\,,\quad\quad\quad\quad{\rm
and}\quad\quad\quad\quad{\cal K}^{\rm LO}(\xi\mu_{0})=\frac{d\sigma^{\rm
LO}(\xi\mu_{0})dX}{d\sigma^{\rm LO}(\mu_{0})/dX}\,,$ (24)

where $\mu_{0}$ is the central value of the scale and $X$ represents the
plotted observable. The error bands are determined similarly to the case of
the integrated cross-section case, i.e. with the help of the 7-point scale
variation that is utilised bin-by-bin.

We start our discussion with the transverse momentum of the $b$ jets that are
presented in Figure 5. The four $b$-jets are ordered according to their
$p_{T}$, thus, $p_{T,\,b_{1}}$ $(p_{T,\,b_{4}})$ corresponds to the transverse
momentum of the hardest (softest) $b$-jet. We notice that the size of the LO
scale uncertainties is similar to that at the integrated cross-section level.
The size of the ${\cal K}$-factor, on the other hand, changes from ${\cal
K}=0.98$ to ${\cal K}=1.23$ depending on the $b$-jet and the phase-space
region. Thus, distortions in the shape of the $b$-jet spectra up to $25\%$ can
be observed after NLO QCD corrections are taken into account. We can further
notice a substantial reduction in the size of theoretical uncertainties when
moving from LO to NLO. The latter uncertainties are at the $20\%$ level.

A fairly similar picture emerges for the transverse momenta of the four
charged leptons that are depicted in Figure 6. In this case, we receive NLO
QCD corrections in the range of $(-7\%,+15\%)$, while NLO QCD uncertainties
are consistently at the $20\%$ level, thus, again almost $4$ times smaller
than at the LO level.

Next, we shift our focus to certain dimensionless observables, highlighting
the most noteworthy among those we have examined. In Figure 7 we plot the
angular separation between the first two hardest $b$-jets $(\Delta
R_{b_{1}b_{2}})$ and the first and the third hardest $b$-jet $(\Delta
R_{b_{1}b_{3}})$. Also displayed is the angular separation between the two
softest $b$-jets $(\Delta R_{b_{3}b_{4}})$ as well as the hardest $b$-jet and
the hardest lepton $(\Delta R_{b_{1}l_{1}})$. In the case of $\Delta
R_{b_{1}b_{2}}$ and $\Delta R_{b_{1}b_{3}}$ the $b$-jets are mostly produced
in back-to-back configurations. For the $\Delta R_{b_{3}b_{4}}$ observable, on
the other hand, a significant number of events can also be found for small
values of $\Delta R_{b_{3}b_{4}}$. A distinct feature can be noticed for the
$\Delta R_{b_{1}\ell_{1}}$ differential cross-section distribution. In
particular, beyond the typical back-to-back configuration, a notable
additional peak near $\Delta R_{b_{1}\ell_{1}}\approx 0.5$ is present. This
configuration corresponds to the decay of highly energetic top quarks,
resulting in small angular separations between their decay products. In all
four cases, the NLO QCD corrections are in the range of $(3\%-22\%)$. Also in
the case of dimensionless observables, NLO uncertainties are rather constant
and at the level of $20\%$. This should be compared to the corresponding LO
uncertainties that are up to even $83\%$.

The observables discussed so far share a common characteristic, namely their
NLO QCD corrections are rather stable and of the order of $20\%$. Moreover,
the corresponding NLO uncertainties are of a similar magnitude. Finally, the
NLO uncertainty bands consistently coincide with the corresponding LO ones.
However, there are certain types of observables for which NLO QCD corrections
can become very large, particularly, in some phase-space regions and the NLO
predictions might no longer fall within the LO uncertainty bands. For
illustrative purposes, in Figure 8, we show differential cross-section
distributions for the transverse momentum of the $b_{1}b_{2}$ and
$b_{1}b_{2}b_{3}b_{4}$ systems. At LO, in the rest frame of the (anti-)top
quark the maximum transverse momentum of the resulting $b$-jet is given by
$p_{T,\,b}=(m_{t}^{2}-m_{W}^{2})/2m_{t}\approx 67.5$ GeV. In order to produce
$b$-jets with $p_{T,\,b}>67.5$ GeV they need to be boosted via the $p_{T}$ of
the parent particle. Consequently, at LO the $b\bar{b}b\bar{b}$ system with
high $p_{T}$ is kinematically strongly disfavoured. At NLO, on the other hand,
the $t\bar{t}t\bar{t}$ system can acquire large transverse momentum by
recoiling against the extra light jet, resulting in a significant enhancement
of events at large values of $p_{T,\,b_{1}b_{2}b_{3}b_{4}}$. Indeed, we can
observe that the differential ${\cal K}$-factor starts to rapidly increase for
$p_{T,\,b_{1}b_{2}b_{3}b_{4}}\gtrsim 270$ GeV and NLO QCD predictions can be
as much as $2.4$ times larger than the LO ones in the tails of the
$p_{T,\,b_{1}b_{2}b_{3}b_{4}}$ distribution. We can also notice that the NLO
theoretical uncertainties in those phase-space regions substantially increase
and are of the order of $50\%$. Finally, the NLO results are not within the
corresponding LO uncertainties that are at the level of $76\%$. Similar but
less pronounced effects can also be observed for the $p_{T,\,b_{1}b_{2}}$
observable in the same plotted range. In this case, the NLO QCD corrections
increase maximally to $45\%$, while NLO uncertainties to $30\%$. To illustrate
the strong dependence of these two observables on extra jet radiation, we can
analyse them again, this time also applying a veto on the additional light- or
$b$-jet. The difference between the NLO cross section without and with this
jet veto is given by a LO calculation, since in the latter case the existence
of the five well-separated jets is demanded. While the four $b$-jets still
need to have $p_{T,\,b}>25$ GeV, the additional jet, if resolved and passing
all the cuts that are also required for the $b$-jet, should additionally have
$p_{T,\,j}<p_{T}^{veto}$, with $p_{T}^{veto}=(25,50,100,200,300,500)$ GeV,
where $j=g,u,d,c,s,b$. Jet vetos significantly affect the normalisation of
differential cross-section distributions, consequently, we study only their
normalised versions. In Figure 9 we show the normalised LO and NLO QCD
differential cross-section distributions for $p_{T,\,b_{1}b_{2}b_{3}b_{4}}$
and $p_{T,\,b_{1}b_{2}}$. Also depicted are the NLO QCD results with the jet
veto of $p_{T}^{veto}=(25,50,100,200,300,500)$ GeV. For various values of the
$p_{T}^{veto}$ cut large differences in the differential factor ${\cal K}$ can
be observed. The substantial higher-order QCD corrections that we could
previously see in the tails of the two distributions have now been reduced to
just a few per cent. For $p_{T}^{veto}=25$ GeV we obtain negative values of
the NLO cross section for $p_{T}\geq 470$ GeV for both observables. Thus, this
particular value is too low and should not be considered in experimental
analyses. It should be noted, however, that the jet veto might substantially
increase the size of the NLO uncertainty bands. For example, for
$p_{T}^{veto}=50$ GeV the NLO uncertainties increase by a factor of $3-7$
depending on the observable. The observed large scale dependence for the more
exclusive NLO results for $p_{T,\,b_{1}b_{2}}$ and
$p_{T,\,b_{1}b_{2}b_{3}b_{4}}$ gives cause for concern, as jet vetoes are
widely used in experimental analyses at the LHC. For example, vetoing the
additional jet activity is crucial for suppressing various SM background
processes or in various searches for new physics effects. Indeed, in many
analyses carried out by the ATLAS and CMS collaborations the most common jet
veto scheme comprises imposing a maximum transverse momentum cut on
anti-$k_{T}$ jets. Therefore, in such cases, in order to reduce theoretical
uncertainties and to improve the convergence of perturbative predictions,
approaches similar to those described in e.g. Ref. Berger:2010xi ;
Stewart:2010tn ; Becher:2012qa ; Banfi:2012yh ; Banfi:2012jm ; Michel:2018hui
; Campbell:2023cha  would have to be used.

Figure 10: Differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with $\sqrt{s}=$
13.6 TeV. The transverse momentum of the hardest $b$-jet and charged lepton
are presented together with the $\Delta R_{bb}$ separation between the hardest
and softest $b$-jets for $\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$. The upper panels
show the absolute NLO QCD predictions for the three PDF sets MSHT20 (solid red
curve), CT18 (dashed blue curve) and NNPDF3.1 (dashed green curve). The middle
panels display the ratio to the results with the MSHT20 PDF set as well as
their scale dependence. The lower panels present the relative size of internal
PDF uncertainties calculated separately for each PDF set. Monte Carlo
integration errors are displayed in upper panels only.

Figure 11: Differential cross-section distributions for the $pp\to
t\bar{t}t\bar{t}$ process in the $4\ell$ channel at the LHC with
$\sqrt{s}=13.6$ TeV. The upper panels show the absolute NLO QCD predictions
for $d\sigma^{\rm NLO}_{\rm exp}/dX$ (blue curve), $d\sigma^{\rm NLO}_{\rm
full}/dX$ (red curve) and $d\sigma^{\rm NLO}_{\rm LO_{dec}}/dX$ (green curve),
where $X=M_{b_{1}b_{2}},\,p_{T,\,\ell_{1}},\,p_{T,\,b_{1}}$ and
$p_{T,\,b_{1}b_{2}}$. Results are given for $\mu_{R}=\mu_{F}=\mu_{0}=E_{T}/4$
and the NLO MSHT20 PDF set. Also provided are theoretical uncertainties as
obtained from the scale dependence. The two middle panels display the ratios
of $d\sigma^{\rm NLO}_{\rm full}$ and $d\sigma^{\rm NLO}_{\rm LO_{dec}}$ to
$d\sigma^{\rm NLO}_{\rm exp}$ together with their relative NLO scale
uncertainties. The lowest panels compare the relative size of NLO scale
uncertainties for the three approaches, normalised to their corresponding NLO
results. Monte Carlo integration errors are displayed in all panels except the
last one.

In the next step, we turn our attention to the PDF uncertainties. We have
already shown that at the integrated fiducial cross-section level, such
uncertainties are negligible when compared to the theoretical uncertainties
coming from scale variation. In the following, we would like to check whether
this is still the case at the differential cross-section level. We again
present NLO QCD results for the same three PDF sets: MSHT20, CT18, and
NNPDF3.1. We focus on investigating the consistency of these predictions and
the size of the corresponding PDF uncertainties in different regions of the
fiducial phase space. In Figure 10 we present afresh the following observables
$p_{T,\,b_{1}}$, $p_{T,\,\ell_{1}}$, $\Delta R_{b_{1}b_{2}}$ and $\Delta
R_{b_{3}b_{4}}$. This time, however, the upper panels display the absolute
differential cross-section distributions at NLO in QCD for the three PDF sets:
MSHT20 (solid red curve), CT18 (dashed blue curve) and NNPDF3.1 (dashed green
curve). The middle panels depict the ratio to the NLO result with the MSHT20
PDF set as well as its scale dependence. Finally, the lower panels present the
relative size of the internal PDF uncertainties separately for each PDF set.
In this case, all curves are normalised to their corresponding central
predictions for $\mu_{0}=E_{T}/4$.

We can confirm the stability of our predictions with respect to the different
PDF sets employed. For example, for the dimensionful observables, the
differences between the NLO results are up to $5\%$, $7\%$ for the
$p_{T,\,\ell_{1}}$, $p_{T,\,b_{1}}$ distribution, respectively. The results
with the NNPDF3.1 PDF set show the most significant differences with respect
to the distributions obtained with the default MSHT20 PDF set. Nevertheless,
any deviations observed are well within the range of the scale uncertainties,
represented by the red uncertainty bands, that are in the range of
$16\%-23\%$. By examining the lower panels of the $p_{T,\,\ell_{1}}$ and
$p_{T,\,b_{1}}$ distributions, we can conclude that the internal PDF
uncertainties tend to increase in the high-$p_{T}$ regions, where they account
for a substantial fraction of the theoretical error. Indeed, for all three PDF
sets these uncertainties are of the order of $10\%$. In the case of
dimensionless observables the differences among the NLO QCD results for
MSHT20, CT18 and NNPDF3.1 are negligible, while the internal PDF uncertainties
are well below $10\%$. In summary, after examining various observables, we
conclude that different PDF sets yield consistent results for the central
value of the scale. However, the internal PDF uncertainties may be more
significant in the tails of some dimensionful observables. In particular, they
can reach half the value of the theoretical error due to scale variation.

In the end, we examine the previously defined three approaches $d\sigma^{\rm
NLO}_{\rm exp}/dX$, $d\sigma^{\rm NLO}_{\rm full}/dX$ and $d\sigma^{\rm
NLO}_{\rm LO_{dec}}/dX$ as a function of the observable $X$, where
$X=M_{b_{1}b_{2}},\,p_{T,\,\ell_{1}},\,p_{T,\,b_{1}}$ and
$p_{T,\,b_{1}b_{2}}$. These results are shown in Figure 11. The upper panels
show the absolute NLO QCD predictions for the three approaches together with
their corresponding theoretical uncertainties as obtained from the scale
dependence. The two middle panels display the ratios of $d\sigma^{\rm
NLO}_{\rm full}/dX$ and $d\sigma^{\rm NLO}_{\rm LO_{dec}}/dX$ to the
$d\sigma^{\rm NLO}_{\rm exp}/dX$ result, together with their relative
uncertainties. The lowest panels compare again the relative size of NLO scale
uncertainties for the three approaches, this time, however, the findings are
normalised to their corresponding NLO results for $\mu_{0}=E_{T}/4$. For the
three observables studied $M_{b_{1}b_{2}},\,p_{T,\,\ell_{1}}$ and
$p_{T,\,b_{1}}$, our observations and conclusions are very similar. The
differences between the full results and the expanded ones are of the order of
$10\%$, which should be compared to the corresponding NLO uncertainties that
are in the range of $15\%-25\%$. On the other hand, NLO QCD corrections in
top-quark decays are consistently at the level of $10\%$. Larger effects could
be observed for $p_{T,\,b_{1}b_{2}}$, especially in the tail of the
distribution. In this case the differences between $d\sigma^{\rm NLO}_{\rm
exp}/dp_{T,\,b_{1}b_{2}}$ and $d\sigma^{\rm NLO}_{\rm
full}/dp_{T,\,b_{1}b_{2}}$ increase to $20\%$. However, the NLO uncertainties
are also increasing and are now in the range of $25\%-35\%$. Nevertheless, the
size of higher-order QCD effects in top-quark decays remains similar to the
case of the other three observables. In short, the overall picture is rather
similar to that already observed at the integrated fiducial cross-section
level. Firstly, the magnitude of the NLO QCD uncertainties increases when
higher-order QCD corrections to top-quark decays are neglected. The latter
effects are consistently at the $10\%$ level. Secondly, the differences
between the expanded and unexpanded NLO results are within the obtained NLO
uncertainties. However, the lack of the expansion of the $1/\Gamma_{t}$ factor
in $d\sigma^{\rm NLO}_{\rm full}/dX$ significantly reduces the size of the NLO
uncertainties for all the observables that we have examined.

##  7 Summary

In this paper we have calculated NLO QCD corrections to $pp\to
t\bar{t}t\bar{t}$ in the $4\ell$ top-quark decay channel for the LHC Run III
energy of $\sqrt{s}=13.6$ TeV. We have taken into account higher-order QCD
effects in both the production and decays of the top quarks. The top-quark
decays have been treated in the NWA, which maintains the spin correlations
throughout the calculation. This is the first time that such a complete study
for this process has been conducted at the NLO level in QCD. Indeed, up to
now, top-quark decays in the $pp\to t\bar{t}t\bar{t}$ process have been
treated only within parton shower frameworks.

In the current study, the integrated fiducial cross sections and their
theoretical uncertainties have been evaluated for two different scale
settings: $\mu_{0}=2m_{t}$ and $\mu_{0}=E_{T}/4$ as well as for the three PDF
sets: MSHT20, CT18, and NNPDF3.1. For the default MSHT20 PDF set the NLO QCD
corrections are rather modest, of the order of $10\%$ only, regardless of the
scale setting. For the other PDF sets, however, NLO QCD corrections might
increase up to even $30\%$. The differences in the size of the ${\cal
K}$-factor are driven by the underlining LO cross sections, and specifically
by the value of $\alpha_{s}(m_{Z})$ used in the LO PDF sets. The inclusion of
higher-order QCD effects has a significant impact on the theoretical error
from the $7$-point scale variation, which has decreased substantially from
$74\%$ at LO to $20\%$ at NLO. Apart from the theoretical error due to the
scale dependence, for all three PDF sets we have also calculated the internal
PDF uncertainties. These uncertainties are very small, of the order of
$2\%-6\%$, depending on the PDF set used. In addition, the differences among
the NLO results as evaluated with MSHT20, CT18, and NNPDF3.1 are of the order
of $1\%$ only. Consequently, the PDF uncertainties are well below the
theoretical uncertainties due to the scale dependence. The latter remains the
dominant source of the theoretical error. Furthermore, we have shown that the
inclusion of NLO QCD corrections in top-quark decays is necessary, as these
higher-order effects are not negligible. Indeed, not only are they at the
level of $10\%$, but they also impact the size of the NLO theoretical error.
Their omission increases the overall magnitude of this error. Finally, we have
checked that the NLO QCD result calculated with the $1/\Gamma_{t}$ term not
expanded has the smallest theoretical uncertainties, which are at the level of
$15\%$ instead of $20\%$.

At the differential cross-section level, the size of the NLO QCD corrections
depends on the studied observable and the examined region of the fiducial
phase space. For the majority of the observables, and when employing the
dynamical scale setting, $\mu_{0}=E_{T}/4$, the NLO QCD corrections are rather
stable and moderate, of the order of $20\%$. Moreover, the corresponding NLO
uncertainties are of similar magnitude. On top of that, the NLO uncertainty
bands consistently coincide with the corresponding LO ones. However, there are
certain types of observables for which not only the NLO QCD corrections can be
very large, even up to $140\%$, but also the NLO uncertainty bands are outside
of the LO ones. Finally, they comprise large NLO uncertainties. These are
dimensionful observables that are very sensitive to additional jet radiation.
We have been able to identify the following differential cross-section
distributions for which this is the case: $p_{T,\,b_{1}b_{2}b_{3}b_{4}}$,
$p_{T,\,\ell_{1}\ell_{2}\ell_{3}\ell_{4}}$, $p_{T}^{miss}$,
$p_{T,\,b_{1}b_{2}}$ and $p_{T,\,\ell_{1}\ell_{2}}$.

In the next step, we have compared the following NLO predictions:
$d\sigma^{\rm NLO}_{\rm exp}/dX$, $d\sigma^{\rm NLO}_{\rm full}/dX$ and
$d\sigma^{\rm NLO}_{\rm LO_{dec}}/dX$ as a function of various observables
$X$. By analysing these predictions we have concluded that the overall pattern
is rather similar to that which has already been observed at the integrated
fiducial cross-section level. Namely, the magnitude of the NLO uncertainties
increases when NLO QCD corrections to top-quark decays are neglected. The
higher-order corrections in top-quark decays are rather stable and of the
order of $10\%$. In addition, the difference between the expanded and
unexpanded NLO results is within the obtained NLO scale uncertainties.
Nonetheless, the lack of the expansion of the $1/\Gamma_{t}$ term in
$d\sigma^{\rm NLO}_{\rm full}/dX$, significantly reduces the size of the NLO
uncertainties for all the observables that we have scrutinised.

In conclusion, we would like to mention that it would be beneficial to conduct
a comparison of the results obtained in this work with the ones from
literature where top-quark decays are generated with the help of parton
showers. Such a comparison could assess the extent to which parton shower
effects can imitate higher-order effects in top-quark decays. We would also be
able to verify the importance of NLO QCD corrections to angular observables
that are sensitive to top-quark spin correlations. Furthermore, such a
comparison could help to identify the phase-space regions that are indeed
affected by resummed dominant soft-collinear logarithmic corrections from
parton showers, which are absent from our fixed-order NLO QCD predictions. We
plan to conduct such a comparison in the near future.

###### Acknowledgements.

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. Support by a grant of the Bundesministerium für Bildung
und Forschung (BMBF) is additionally acknowledged. The authors gratefully
acknowledge the computing time provided to them at the NHR Center NHR4CES at
RWTH Aachen University (project number p0020216). 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.

## References

  * (1) Q.-H. Cao, S.-L. Chen and Y. Liu, _Probing Higgs Width and Top Quark Yukawa Coupling from $t\bar{t}H$ and $t\bar{t}t\bar{t}$ Productions_, _Phys. Rev. D_ 95 (2017) 053004 [1602.01934].
  * (2) Q.-H. Cao, S.-L. Chen, Y. Liu, R. Zhang and Y. Zhang, _Limiting top quark-Higgs boson interaction and Higgs-boson width from multitop productions_ , _Phys. Rev. D_ 99 (2019) 113003 [1901.04567].
  * (3) C.R. Schmidt and M.E. Peskin, _A Probe of CP violation in top quark pair production at hadron supercolliders_ , _Phys. Rev. Lett._ 69 (1992) 410.
  * (4) J.H. Kühn, A. Scharf and P. Uwer, _Weak Interactions in Top-Quark Pair Production at Hadron Colliders: An Update_ , _Phys. Rev. D_ 91 (2015) 014020 [1305.5773].
  * (5) T. Martini, R.-Q. Pan, M. Schulze and M. Xiao, _Probing the CP structure of the top quark Yukawa coupling: Loop sensitivity versus on-shell sensitivity_ , _Phys. Rev. D_ 104 (2021) 055045 [2104.04277].
  * (6) P. Artoisenet et al., _A framework for Higgs characterisation_ , _JHEP_ 11 (2013) 043 [1306.6464].
  * (7) F. Demartin, F. Maltoni, K. Mawatari, B. Page and M. Zaro, _Higgs characterisation at NLO in QCD: CP properties of the top-quark Yukawa interaction_ , _Eur. Phys. J. C_ 74 (2014) 3065 [1407.5089].
  * (8) F. Demartin, F. Maltoni, K. Mawatari and M. Zaro, _Higgs production in association with a single top quark at the LHC_ , _Eur. Phys. J. C_ 75 (2015) 267 [1504.00611].
  * (9) F. Demartin, B. Maier, F. Maltoni, K. Mawatari and M. Zaro, _tWH associated production at the LHC_ , _Eur. Phys. J. C_ 77 (2017) 34 [1607.05862].
  * (10) H. Bahl, P. Bechtle, S. Heinemeyer, J. Katzy, T. Klingl, K. Peters et al., _Indirect $\mathcal{CP}$ probes of the Higgs-top-quark interaction: current LHC constraints and future opportunities_, _JHEP_ 11 (2020) 127 [2007.08542].
  * (11) B. Bortolato, J.F. Kamenik, N. Košnik and A. Smolkovič, _Optimized probes of $CP$ -odd effects in the $t\bar{t}h$ process at hadron colliders_, _Nucl. Phys. B_ 964 (2021) 115328 [2006.13110].
  * (12) J. Hermann, D. Stremmer and M. Worek, _$\mathcal{CP}$ structure of the top-quark Yukawa interaction: NLO QCD corrections and off-shell effects_, _JHEP_ 09 (2022) 138 [2205.09983].
  * (13) B. Lillie, J. Shu and T.M.P. Tait, _Top Compositeness at the Tevatron and LHC_ , _JHEP_ 04 (2008) 087 [0712.3057].
  * (14) A. Pomarol and J. Serra, _Top Quark Compositeness: Feasibility and Implications_ , _Phys. Rev. D_ 78 (2008) 074026 [0806.3247].
  * (15) K. Kumar, T.M.P. Tait and R. Vega-Morales, _Manifestations of Top Compositeness at Colliders_ , _JHEP_ 05 (2009) 022 [0901.3808].
  * (16) G. Cacciapaglia, R. Chierici, A. Deandrea, L. Panizzi, S. Perries and S. Tosi, _Four tops on the real projective plane at LHC_ , _JHEP_ 10 (2011) 042 [1107.4616].
  * (17) T. Plehn and T.M.P. Tait, _Seeking Sgluons_ , _J. Phys. G_ 36 (2009) 075001 [0810.3919].
  * (18) S. Calvet, B. Fuks, P. Gris and L. Valery, _Searching for sgluons in multitop events at a center-of-mass energy of 8 TeV_ , _JHEP_ 04 (2013) 043 [1212.3360].
  * (19) J.H. Kim, K. Kong, S.J. Lee and G. Mohlabeng, _Probing TeV scale Top-Philic Resonances with Boosted Top-Tagging at the High Luminosity LHC_ , _Phys. Rev. D_ 94 (2016) 035023 [1604.07421].
  * (20) L. Darmé, B. Fuks and F. Maltoni, _Top-philic heavy resonances in four-top final states and their EFT interpretation_ , _JHEP_ 09 (2021) 143 [2104.09512].
  * (21) P.S. Bhupal Dev and A. Pilaftsis, _Maximally Symmetric Two Higgs Doublet Model with Natural Standard Model Alignment_ , _JHEP_ 12 (2014) 024 [1408.3405].
  * (22) N. Craig, F. D’Eramo, P. Draper, S. Thomas and H. Zhang, _The Hunt for the Rest of the Higgs Bosons_ , _JHEP_ 06 (2015) 137 [1504.04630].
  * (23) S. Gori, I.-W. Kim, N.R. Shah and K.M. Zurek, _Closing the Wedge: Search Strategies for Extended Higgs Sectors with Heavy Flavor Final States_ , _Phys. Rev. D_ 93 (2016) 075038 [1602.02782].
  * (24) N. Craig, J. Hajer, Y.-Y. Li, T. Liu and H. Zhang, _Heavy Higgs bosons at low $\tan\beta$: from the LHC to 100 TeV_, _JHEP_ 01 (2017) 018 [1605.08744].
  * (25) C. Zhang, _Constraining $qqtt$ operators from four-top production: a case for enhanced EFT sensitivity_, _Chin. Phys. C_ 42 (2018) 023104 [1708.05928].
  * (26) SMEFiT collaboration, _Combined SMEFT interpretation of Higgs, diboson, and top quark data from the LHC_ , _JHEP_ 11 (2021) 089 [2105.00006].
  * (27) R. Aoude, H. El Faham, F. Maltoni and E. Vryonidou, _Complete SMEFT predictions for four top quark production at hadron colliders_ , _JHEP_ 10 (2022) 163 [2208.04962].
  * (28) E. Alvarez, D.A. Faroughy, J.F. Kamenik, R. Morales and A. Szynkman, _Four tops for LHC_ , _Nucl. Phys. B_ 915 (2017) 19 [1611.05032].
  * (29) ATLAS collaboration, _Observation of four-top-quark production in the multilepton final state with the ATLAS detector_ , _Eur. Phys. J. C_ 83 (2023) 496 [2303.15061].
  * (30) CMS collaboration, _Observation of four top quark production in proton-proton collisions at $\sqrt{s}=13$ TeV_, _Phys. Lett. B_ 847 (2023) 138290 [2305.13439].
  * (31) G. Bevilacqua and M. Worek, _Constraining BSM Physics at the LHC: Four top final states with NLO accuracy in perturbative QCD_ , _JHEP_ 07 (2012) 111 [1206.3064].
  * (32) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _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 [1405.0301].
  * (33) F. Maltoni, D. Pagani and I. Tsinikos, _Associated production of a top-quark pair with vector bosons at NLO in QCD: impact on $t\bar{t}H$ searches at the LHC_, _JHEP_ 02 (2016) 113 [1507.05640].
  * (34) R. Frederix, D. Pagani and M. Zaro, _Large NLO corrections in $t\bar{t}W^{\pm}$ and $t\bar{t}t\bar{t}$ hadroproduction from supposedly subleading EW contributions_, _JHEP_ 02 (2018) 031 [1711.02116].
  * (35) M. van Beekveld, A. Kulesza and L.M. Valero, _Threshold Resummation for the Production of Four Top Quarks at the LHC_ , _Phys. Rev. Lett._ 131 (2023) 211901 [2212.03259].
  * (36) T. Ježo and M. Kraus, _Hadroproduction of four top quarks in the powheg box_ , _Phys. Rev. D_ 105 (2022) 114024 [2110.15159].
  * (37) T. Sjostrand, S. Mrenna and P.Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05 (2006) 026 [hep-ph/0603175].
  * (38) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (39) P. Nason, _A New method for combining NLO QCD with shower Monte Carlo algorithms_ , _JHEP_ 11 (2004) 040 [hep-ph/0409146].
  * (40) S. Alioli, P. Nason, C. Oleari and E. Re, _A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX_ , _JHEP_ 06 (2010) 043 [1002.2581].
  * (41) S. Frixione and B.R. Webber, _Matching NLO QCD computations and parton shower simulations_ , _JHEP_ 06 (2002) 029 [hep-ph/0204244].
  * (42) S. Frixione, P. Nason and B.R. Webber, _Matching NLO QCD and parton showers in heavy flavor production_ , _JHEP_ 08 (2003) 007 [hep-ph/0305252].
  * (43) V.S. Fadin, V.A. Khoze and A.D. Martin, _How suppressed are the radiative interference effects in heavy instable particle production?_ , _Phys. Lett. B_ 320 (1994) 141 [hep-ph/9309234].
  * (44) K. Melnikov and O.I. Yakovlev, _Top near threshold: All alpha-S corrections are trivial_ , _Phys. Lett. B_ 324 (1994) 217 [hep-ph/9302311].
  * (45) K. Melnikov and O.I. Yakovlev, _Final state interaction in the production of heavy unstable particles_ , _Nucl. Phys. B_ 471 (1996) 90 [hep-ph/9501358].
  * (46) A. Denner, S. Dittmaier and M. Roth, _Nonfactorizable photonic corrections to $e^{+}e^{-}\to W^{+}W^{-}\to 4$ fermions_, _Nucl. Phys. B_ 519 (1998) 39 [hep-ph/9710521].
  * (47) W. Beenakker, F.A. Berends and A.P. Chapovsky, _One loop QCD interconnection effects in pair production of top quarks_ , _Phys. Lett. B_ 454 (1999) 129 [hep-ph/9902304].
  * (48) K. Melnikov and M. Schulze, _NLO QCD corrections to top quark pair production and decay at hadron colliders_ , _JHEP_ 08 (2009) 049 [0907.3090].
  * (49) J.M. Campbell and R.K. Ellis, _Top-Quark Processes at NLO in Production and Decay_ , _J. Phys. G_ 42 (2015) 015005 [1204.1513].
  * (50) G. Bevilacqua, H.B. Hartanto, M. Kraus, T. Weber and M. Worek, _Off-shell vs on-shell modelling of top quarks in photon associated production_ , _JHEP_ 03 (2020) 154 [1912.09999].
  * (51) M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to leptonic observables in top-quark pair production and decay_ , _JHEP_ 05 (2021) 212 [2008.11133].
  * (52) R.V. Harlander, S.Y. Klein and M. Lipp, _FeynGame_ , _Comput. Phys. Commun._ 256 (2020) 107465 [2003.00896].
  * (53) G. Bevilacqua, M. Czakon, M.V. Garzelli, A. van Hameren, A. Kardos, C.G. Papadopoulos et al., _HELAC-NLO_ , _Comput. Phys. Commun._ 184 (2013) 986 [1110.1499].
  * (54) P. Draggiotis, R.H.P. Kleiss and C.G. Papadopoulos, _On the computation of multigluon amplitudes_ , _Phys. Lett. B_ 439 (1998) 157 [hep-ph/9807207].
  * (55) A. Kanaki and C.G. Papadopoulos, _HELAC: A Package to compute electroweak helicity amplitudes_ , _Comput. Phys. Commun._ 132 (2000) 306 [hep-ph/0002082].
  * (56) P.D. Draggiotis, R.H.P. Kleiss and C.G. Papadopoulos, _Multijet production in hadron collisions_ , _Eur. Phys. J. C_ 24 (2002) 447 [hep-ph/0202201].
  * (57) A. Cafarella, C.G. Papadopoulos and M. Worek, _Helac-Phegas: A Generator for all parton level processes_ , _Comput. Phys. Commun._ 180 (2009) 1941 [0710.2427].
  * (58) A. van Hameren, C.G. Papadopoulos and R. Pittau, _Automated one-loop calculations: A Proof of concept_ , _JHEP_ 09 (2009) 106 [0903.4665].
  * (59) G. Ossola, C.G. Papadopoulos and R. Pittau, _CutTools: A Program implementing the OPP reduction method to compute one-loop amplitudes_ , _JHEP_ 03 (2008) 042 [0711.3596].
  * (60) G. Ossola, C.G. Papadopoulos and R. Pittau, _Reducing full one-loop amplitudes to scalar integrals at the integrand level_ , _Nucl. Phys. B_ 763 (2007) 147 [hep-ph/0609007].
  * (61) G. Ossola, C.G. Papadopoulos and R. Pittau, _On the Rational Terms of the one-loop amplitudes_ , _JHEP_ 05 (2008) 004 [0802.1876].
  * (62) P. Draggiotis, M.V. Garzelli, C.G. Papadopoulos and R. Pittau, _Feynman Rules for the Rational Part of the QCD 1-loop amplitudes_ , _JHEP_ 04 (2009) 072 [0903.0356].
  * (63) A. van Hameren, _OneLOop: For the evaluation of one-loop scalar functions_ , _Comput. Phys. Commun._ 182 (2011) 2427 [1007.4716].
  * (64) S. Badger, B. Biedermann and P. Uwer, _NGluon: A Package to Calculate One-loop Multi-gluon Amplitudes_ , _Comput. Phys. Commun._ 182 (2011) 1674 [1011.2900].
  * (65) S. Catani and M.H. Seymour, _A General algorithm for calculating jet cross-sections in NLO QCD_ , _Nucl. Phys. B_ 485 (1997) 291 [hep-ph/9605323].
  * (66) S. Catani, S. Dittmaier, M.H. Seymour and Z. Trocsanyi, _The Dipole formalism for next-to-leading order QCD calculations with massive partons_ , _Nucl. Phys. B_ 627 (2002) 189 [hep-ph/0201036].
  * (67) M. Czakon, C.G. Papadopoulos and M. Worek, _Polarizing the Dipoles_ , _JHEP_ 08 (2009) 085 [0905.0883].
  * (68) G. Bevilacqua, M. Czakon, M. Kubocz and M. Worek, _Complete Nagy-Soper subtraction for next-to-leading order calculations in QCD_ , _JHEP_ 10 (2013) 204 [1308.5605].
  * (69) Z. Nagy and D.E. Soper, _Parton showers with quantum interference_ , _JHEP_ 09 (2007) 114 [0706.0017].
  * (70) M. Czakon, H.B. Hartanto, M. Kraus and M. Worek, _Matching the Nagy-Soper parton shower at next-to-leading order_ , _JHEP_ 06 (2015) 033 [1502.00925].
  * (71) Z. Nagy and Z. Trocsanyi, _Next-to-leading order calculation of four jet observables in electron positron annihilation_ , _Phys. Rev. D_ 59 (1999) 014020 [hep-ph/9806317].
  * (72) Z. Nagy, _Next-to-leading order calculation of three jet observables in hadron hadron collision_ , _Phys. Rev. D_ 68 (2003) 094002 [hep-ph/0307268].
  * (73) G. Bevilacqua, M. Czakon, C.G. Papadopoulos, R. Pittau and M. Worek, _Assault on the NLO Wishlist: $pp\to t\bar{t}b\bar{b}$_, _JHEP_ 09 (2009) 109 [0907.4723].
  * (74) A. van Hameren, _Kaleu: A General-Purpose Parton-Level Phase Space Generator_ , 1003.4953.
  * (75) A. van Hameren, _PARNI for importance sampling and density estimation_ , _Acta Phys. Polon. B_ 40 (2009) 259 [0710.2448].
  * (76) G. Bevilacqua, M. Czakon, A. van Hameren, C.G. Papadopoulos and M. Worek, _Complete off-shell effects in top quark pair hadroproduction with leptonic decay at next-to-leading order_ , _JHEP_ 02 (2011) 083 [1012.4230].
  * (77) 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].
  * (78) A. Denner, S. Dittmaier, M. Roth and D. Wackeroth, _Electroweak radiative corrections to $e^{+}e^{-}\to W^{+}W^{-}\to 4$ fermions in double pole approximation: The RACOONWW approach_, _Nucl. Phys. B_ 587 (2000) 67 [hep-ph/0006307].
  * (79) E. Accomando, A. Denner and A. Kaiser, _Logarithmic electroweak corrections to gauge-boson pair production at the LHC_ , _Nucl. Phys. B_ 706 (2005) 325 [hep-ph/0409247].
  * (80) A. Denner and M. Pellen, _NLO electroweak corrections to off-shell top-antitop production with leptonic decays at the LHC_ , _JHEP_ 08 (2016) 155 [1607.05571].
  * (81) G. Bevilacqua, M. Lupattelli, D. Stremmer and M. Worek, _Study of additional jet activity in top quark pair production and decay at the LHC_ , _Phys. Rev. D_ 107 (2023) 114027 [2212.04722].
  * (82) D. Stremmer and M. Worek, _Associated production of a top-quark pair with two isolated photons at the LHC through NLO in QCD_ , _JHEP_ 08 (2023) 179 [2306.16968].
  * (83) J. Alwall et al., _A Standard format for Les Houches event files_ , _Comput. Phys. Commun._ 176 (2007) 300 [hep-ph/0609017].
  * (84) I. Antcheva et al., _ROOT: A C++ framework for petabyte data storage, statistical analysis and visualization_ , _Comput. Phys. Commun._ 180 (2009) 2499 [1508.07749].
  * (85) Z. Bern, L.J. Dixon, F. Febres Cordero, S. Höche, H. Ita, D.A. Kosower et al., _Ntuples for NLO Events at Hadron Colliders_ , _Comput. Phys. Commun._ 185 (2014) 1443 [1310.7439].
  * (86) G. Bevilacqua, “unpublished.” 2019.
  * (87) A. Denner, S. Dittmaier, S. Kallweit and S. Pozzorini, _NLO QCD corrections to off-shell top-antitop production with leptonic decays at hadron colliders_ , _JHEP_ 10 (2012) 110 [1207.5018].
  * (88) J. Butterworth et al., _PDF4LHC recommendations for LHC Run II_ , _J. Phys. G_ 43 (2016) 023001 [1510.03865].
  * (89) S. Bailey, T. Cridge, L.A. Harland-Lang, A.D. Martin and R.S. Thorne, _Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs_ , _Eur. Phys. J. C_ 81 (2021) 341 [2012.04684].
  * (90) NNPDF collaboration, _Parton distributions from high-precision collider data_ , _Eur. Phys. J. C_ 77 (2017) 663 [1706.00428].
  * (91) T.-J. Hou et al., _New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC_ , _Phys. Rev. D_ 103 (2021) 014013 [1912.10053].
  * (92) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , _Phys. Rev. D_ 93 (2016) 033006 [1506.07443].
  * (93) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J. C_ 75 (2015) 132 [1412.7420].
  * (94) M. Cacciari, G.P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (95) A. Behring, M. Czakon, A. Mitov, A.S. Papanastasiou and R. Poncelet, _Higher order corrections to spin correlations in top quark pair production at the LHC_ , _Phys. Rev. Lett._ 123 (2019) 082001 [1901.05407].
  * (96) G. Bevilacqua, H.-Y. Bi, H.B. Hartanto, M. Kraus, M. Lupattelli and M. Worek, _$t\bar{t}b\bar{b}$ at the LHC: On the size of off-shell effects and prompt b-jet identification_, _Phys. Rev. D_ 107 (2023) 014028 [2202.11186].
  * (97) G. Bevilacqua, H.B. Hartanto, M. Kraus and M. Worek, _Top Quark Pair Production in Association with a Jet with Next-to-Leading-Order QCD Off-Shell Effects at the Large Hadron Collider_ , _Phys. Rev. Lett._ 116 (2016) 052003 [1509.09242].
  * (98) G. Bevilacqua, H.B. Hartanto, M. Kraus and M. Worek, _Off-shell Top Quarks with One Jet at the LHC: A comprehensive analysis at NLO QCD_ , _JHEP_ 11 (2016) 098 [1609.01659].
  * (99) G. Bevilacqua, H.B. Hartanto, M. Kraus, T. Weber and M. Worek, _Hard Photons in Hadroproduction of Top Quarks with Realistic Final States_ , _JHEP_ 10 (2018) 158 [1803.09916].
  * (100) G. Bevilacqua, H.B. Hartanto, M. Kraus, T. Weber and M. Worek, _Towards constraining Dark Matter at the LHC: Higher order QCD predictions for $t\bar{t}+Z(Z\to\nu_{\ell}\bar{\nu}_{\ell})$_, _JHEP_ 11 (2019) 001 [1907.09359].
  * (101) G. Bevilacqua, H.-Y. Bi, H.B. Hartanto, M. Kraus and M. Worek, _The simplest of them all: $t\bar{t}W^{\pm}$ at NLO accuracy in QCD_, _JHEP_ 08 (2020) 043 [2005.09427].
  * (102) G. Bevilacqua, H.-Y. Bi, H.B. Hartanto, M. Kraus, M. Lupattelli and M. Worek, _$t\bar{t}b\bar{b}$ at the LHC: on the size of corrections and b-jet definitions_, _JHEP_ 08 (2021) 008 [2105.08404].
  * (103) D. Stremmer and M. Worek, _Production and decay of the Higgs boson in association with top quarks_ , _JHEP_ 02 (2022) 196 [2111.01427].
  * (104) G. Bevilacqua, H.B. Hartanto, M. Kraus, J. Nasufi and M. Worek, _NLO QCD corrections to full off-shell production of $t\overline{t}Z$ including leptonic decays_, _JHEP_ 08 (2022) 060 [2203.15688].
  * (105) C.F. Berger, C. Marcantonini, I.W. Stewart, F.J. Tackmann and W.J. Waalewijn, _Higgs Production with a Central Jet Veto at NNLL+NNLO_ , _JHEP_ 04 (2011) 092 [1012.4480].
  * (106) I.W. Stewart, F.J. Tackmann and W.J. Waalewijn, _N-Jettiness: An Inclusive Event Shape to Veto Jets_ , _Phys. Rev. Lett._ 105 (2010) 092002 [1004.2489].
  * (107) T. Becher and M. Neubert, _Factorization and NNLL Resummation for Higgs Production with a Jet Veto_ , _JHEP_ 07 (2012) 108 [1205.3806].
  * (108) A. Banfi, G.P. Salam and G. Zanderighi, _NLL+NNLO predictions for jet-veto efficiencies in Higgs-boson and Drell-Yan production_ , _JHEP_ 06 (2012) 159 [1203.5773].
  * (109) A. Banfi, P.F. Monni, G.P. Salam and G. Zanderighi, _Higgs and Z-boson production with a jet veto_ , _Phys. Rev. Lett._ 109 (2012) 202001 [1206.4998].
  * (110) J.K.L. Michel, P. Pietrulewicz and F.J. Tackmann, _Jet Veto Resummation with Jet Rapidity Cuts_ , _JHEP_ 04 (2019) 142 [1810.12911].
  * (111) J.M. Campbell, R.K. Ellis, T. Neumann and S. Seth, _Jet-veto resummation at N 3LLp \+ NNLO in boson production processes_, _JHEP_ 04 (2023) 106 [2301.11768].