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

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.


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 → t tt 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 → t tH(H → t t) process can be used to constrain the Higgs boson Yukawa coupling [1,2], complementary to loop sensitivity through electroweak corrections in t t production [3][4][5] and tree-level sensitivity in the pp → t tH and/or pp → tH process [6][7][8][9][10][11][12].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 [13][14][15][16], heavy Kaluza-Klein gluons and quarks from extra dimensions [17], scalar gluons occurring in supersymmetric models and produced in pairs [18,19], singly produced top-philic resonances [20,21] or a heavy scalar or pseudoscalar boson from extended electroweak symmetry breaking sectors produced in association with a top-quark pair [22][23][24][25].In addition, the pp → t tt t process provides a direct way to constrain the SM Effective Field Theory Wilson coefficients sensitive to the quartic couplings between top-quarks [26][27][28].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 tt t cross section and cause modifications in various differential cross-section distributions.Given such rich phenomenology of the SM pp → t tt 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 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 → W b, its final states are classified via the W gauge boson decays into W → ℓν ℓ and W → q q ′ .Consequently, the SM pp → t tt t → W + W − W + W − b bb 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 (BR) of 40%.It is followed by the fully hadronic and opposite-sign dilepton modes with BR ≈ 20% each, the same-sign dilepton and the trilepton modes with BR ≈ 10% each and finally by the fully leptonic mode with the smallest contribution of the order of 1.2% [29].The last branching ratio is further reduced to 0.2% when ℓ ± = e ± , µ ± are only considered and the τ lepton is not taken into account.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 [30,31].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σ that is needed to claim observation.Independent of the experiment and taking into account current uncertainties of the measurements, the obtained t tt 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 tW ± production, the normalisation of the pp → t tZ background process and the modelling of the pp → t tt t signal process.Therefore, an accurate theoretical description of t tt 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 → t tt t process with stable top quarks have been calculated for the first time a while ago [32] and a few years later recomputed in Refs.[33,34].Besides NLO QCD corrections, a further step towards more precise predictions for this process has been achieved by evaluating the complete-NLO corrections [35] that additionally include NLO EW higherorder effects and subleading contributions at perturbative orders from O(α 5 s ) to O(α 5 ).However, it has been shown that there are cancellations between different terms in α s and α, so in general the higherorder effects for the pp → t tt t process are dominated by α s corrections to the leading QCD process at O(α 4 s ).Finally, very recently the calculation for pp → t tt t at the next-to-leading logarithmic accuracy has been carried out [36] including also constant O(α s ) non-logarithmic contributions that do not vanish at the absolute production threshold defined as M 2 /ŝ → 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. [37] by matching theoretical predictions for the pp → t tt t process with the dipole-style p T -ordered parton shower from the Pythia8 Monte Carlo (MC) program [38,39] using the Powheg framework [40,41].In this study, top-quark decays in the ℓ+jets channel have been modelled at LO. Consequently, spin correlations are also described at the LO level in the perturbative expansion.
For the first time, not only the impact of NLO QCD corrections and t 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 [33], that employs the MC@NLO matching to parton shower programs [42,43], 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 collinear 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 → t tt 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 → 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 → t tt t process can be decomposed into factorizable and non-factorizable contributions.The latter corrections imply a cross-talk between t tt t production and t/ t decays as well as between the t and t decays.In the limit when Γ/m → 0, these corrections must vanish.The exact way in which they vanish has already been described ample times in literature, see e.g.[44][45][46][47][48][49][50][51][52], and goes under the name of the narrow-width-approximation (NWA).The aim of this paper is to provide theoretical predictions for the pp → t tt t process that include NLO QCD corrections to t tt t production and four top-quark decays in the NWA.In our study, we consider the 4ℓ top-quark decay channel comprising the following decay chain pp , where ℓ ± = e ± , µ ± .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 infraredsafe observables.In addition, such theoretical predictions allow us to qualitatively increase the level of precision and accuracy for the SM pp → t tt 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 √ 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 crosssection 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.

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 q).For the fully leptonic top-quark decay channel that we are focusing on, it can be described as follows: gg/q q → t tt t 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 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.
In addition, we calculate next-to-leading-order (NLO) QCD corrections of the order of O(α 5 s α 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: where in the case of the gg/q q/qq 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 (Γ) of an unstable particle is much smaller than its mass (m), allowing us to safely take the limit Γ/m → 0 in the Breit-Wigner propagator Indeed, this is the case in our process where Γ t /m t ≈ 0.008 and Γ W /m W ≈ 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 σ NLO full and the latter as σ NLO LO dec .Additionally, we are going to assess the effects of expanding the top-quark width in our calculations, referred to as σ NLO 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 [54].In this framework the calculation of scattering amplitudes is based on the Dyson-Schwinger recursive algorithm [55][56][57][58].The Helac-NLO program includes Helac-1loop [59] for the evaluation of the numerators of the loop integrals and the rational terms, CutTools [60], which implements the OPP method for the numerical evaluation of one-loop amplitudes, based on a decomposition at the integrand level [61][62][63], and OneLoop [64] 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 q subprocess, on the other hand, we use the so-called scale test [65].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 [66] to extract the soft and collinear infrared singularities and to combine them with the virtual corrections.Specifically, the formulation presented in Ref. [67] for massive quarks as implemented in the Helac-Dipoles MC program [68] 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 [69] that uses the momentum mapping and the splitting functions derived in the context of an improved parton shower formulation by Nagy and Soper [70].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 [71].For both subtraction schemes a phase-space restriction on the contribution of the subtraction terms, called α max , is implemented [71][72][73][74].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 initialstate 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 [75], including Parni [76] 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 [77].
We have conducted multiple cross-checks to ensure the accuracy of our results.To begin with, we checked the coefficients of the poles in ϵ, i.e. 1/ϵ and 1/ϵ 2 , of the virtual amplitudes for gg and q q subprocesses with MadGraph − aMC@NLO [33] for a few random phase-space points for stable t tt t production and for the t → W + b → ℓ + ν ℓ b decay.In the next step, we repeated this procedure using the Recola MC program [78].This time, however, also the finite parts along with the coefficients of the poles in ϵ 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 [79][80][81], 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 → t tjj [82] and pp → t tγγ [83].For the real corrections, we have checked for several phase-space points that the poles in ϵ that appear in the 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 α 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 [84], which are subsequently converted into the Root Ntuple Files [85,86].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 [87].

Computational setup
We consider the process at NLO in QCD, where ℓ ± stands for ℓ ± = e ± , µ ± .If instead ℓ = e ± , µ ± , τ ± , then all the results presented in this paper should be multiplied by the following factor 81/16 ≈ 5.0625 to obtain the correct normalisation.To perform this calculation we use the NWA.For simplicity, we refer to this process as pp → t tt t in the 4ℓ top-quark decay channel.We do not consider τ 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 √ s = 13.6 TeV.The numerical values we use for the SM input parameters are given below: 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 µ R = m W .We utilise the Γ NLO W 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 Γ t are determined using the formulas provided in Ref.
[88], where α s (µ R = m t ) is used: The electromagnetic coupling α is evaluated in the G µ -scheme as follows: In accordance with the guidance for Standard Model processes provided by the PDF4LHC working group [89], we utilise three individual, modern PDF sets: MSHT20 [90] (our default PDF set), NNPDF3.1 [91] and CT18 [92].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 α s (m Z ) = 0.130 for our LO calculations [93].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 [94] involving five active flavours.
The limitations of fixed-order calculations lead to an inherent dependence on the renormalisation (µ R ) and factorisation scale (µ F ).The theoretical uncertainty of the cross section, associated with neglected higher-order terms in the perturbative expansion, can be estimated by varying µ R and µ F in α s and PDFs, up and down by a factor 2 around the central value of the scale of the process (µ 0 ) in the range Traditionally, the following additional condition is required We are looking for the minimum and maximum of the resulting cross sections, taking into account the following pairs ) , (0.5, 1) , (1, 2) , (1, 1), (1, 0.5), (2, 2), (0.5, 0.5) .
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 µ R scale in the calculation of Γ 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 which has already been employed in our previous studies for the pp → t tt t process with stable top quarks [32].On the other hand, for the dynamical scale setting we utilise with E T defined according to where the transverse masses m T (t) and m T ( t ) are given by This particular scale choice has already been explored in other studies for this process, see e.g.Refs.[32,35,37].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 |η| < 5 with the help of the infrared-safe anti-k T jet algorithm [95] 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 [30,31]: 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 → t tt t process including top-quark decays in the NWA can be given by, see e.g.[49-52, 96, 97] where dσ tt tt denotes the cross section for on-shell pp → t tt t production and dΓ t (dΓt) is the top-quark decay rate for Furthermore, the symbol × represents spin correlations.This formula holds to all orders in α s .However, all the terms on the right-hand side of Eq. (4.1) can be expanded in terms of α s to the required order in perturbative expansion.Thus, for example at the NLO level we have where dσ (0) , dΓ (0) and dσ (1) , dΓ (1) denote LO and NLO contributions to the production and decay processes, respectively.In addition, in Eq. ( 4.1) we have the top-quark width, Γ t , which enters the NWA calculation as a parameter through the denominator.There are different ways to treat the 1/Γ t factor.For example, it is possible to set the value of Γ t to a numerical value corresponding to the perturbative order of the full calculation, Γ NLO t .Alternatively, one can formally expand Γ t in powers of α s and keep only the terms consistent with the order of the performed calculation.Thus, for example at NLO in QCD we would have t .
(4.3)By incorporating Eq. ( 4.2) into Eq.( 4.1) we obtain the following expression for dσ LO at the LO level whereas, at the NLO QCD level we would have instead dσ NLO full given by The latter formula describes the full NLO QCD calculation with higher order QCD corrections incorporated in both t tt t production and top-quark decays.It is worth pointing out that the various terms in Eq. (4.5) 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 Γ NLO t in Eq. (4.5), which is a function of the strong coupling constant, spoils the rigorous expansion of the cross section in powers of α s .If also Eq. (4.3) is incorporated into Eq.(4.1) we would rather have the following expression that represents the corresponding expanded NLO QCD result, labelled as dσ NLO exp .We note that in Eq. (4.6) the dσ LO part is calculated with the NLO PDF set, thus, with the two-loop running of α 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σ NLO  LO dec .The contribution to this result comes from the first line of Eq. (4.5) only, while Γ LO t is used in the calculation.Thus, we can write for dσ NLO 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 σ NLO exp (our default findings at NLO QCD) and σ NLO LO dec at the integrated and differential cross-section level.Furthermore, we show the results for the σ NLO 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.

Integrated fiducial cross sections
We begin the presentation of our results for the pp → t tt t process in the 4ℓ top-quark decay channel with a discussion of the integrated fiducial cross section for the LHC Run III energy of √ 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, µ R = µ F = µ 0 = 2m t , and the default MSHT20 PDF set, are given by: We use the expanded version of the NWA at NLO in QCD, which for the sake of brevity is denoted simply by σ NLO .On the other hand, when the phase-space dependent scale choice is employed, , we obtain: 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 q/qq channels collectively make up about 12%.The NLO QCD corrections for this process are rather moderate, in the range of 9% − 12%, 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.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: µ 0 = 2m t and µ 0 = E T /4.The K factor, defined as K = σ NLO /σ 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 TeV. Results are presented for MSHT20, NNPDF3.1 and CT18 PDF sets.They are evaluated using µ 0 = 2m t and µ 0 = E T /4.Also given are theoretical uncertainties coming from the scale variation and from PDFs.MC integration errors are displayed in parentheses.In the last column the K-factor is shown.they are of the order of 6%.We observe substantial variations in the values of the LO cross section for the three LO PDF sets that are provided for α s (m Z ) = 0.130.These differences, ranging from 5% to 15%, result in a variety of the K factor: K = (1.25 − 1.30) for NNPDF3.1 and K = (1.03− 1.04) for CT18.As expected, using NLO PDF sets instead for the LO calculations would correspond to lower values of σ LO and, consequently, higher K factors.Similar effects would be observed when using LO PDF sets with α s (m Z ) = 0.118.In order to substantiate these statements in Table 2 we present LO integrated fiducial cross sections as obtained with the following PDF sets provided by the NNPDF collaboration: the LO NNPDF3.1 PDF set with α s (m Z ) = 0.130 and α s (m Z ) = 0.118 as well as the NLO NNPDF3.1 PDF set.For consistency and readability, the NLO cross section is also displayed.Again, the results are given for two scale settings µ 0 = 2m t and µ 0 = E T /4.We can observe that the LO results are simply not adequate to describe the pp → t tt t process in the 4ℓ channel.They are sensitive not only to scale variation, but also to the choice of PDF set and the value of α s (m Z ) used.This large dependence is reflected in the obtained K-factor, which for the process at hand ranges from K = 1.25 to K = 1.68 within a given PDF set and from K = 1.03 to K = 1.68 once the three PDF sets, which are recommended by the PDF4LHC working group, are employed.On the other hand, at the NLO QCD level and independent of the scale setting, the differences among σ NLO for the three PDF sets (MSHT20, NNPDF3.1 and CT18) are of the order of 1% only.Generally, as explained by the various PDF collaborations there are large differences between the LO and NLO PDF sets, both for central values and uncertainties.In particular, the LO uncertainties are significant due to the poor Also given are NLO predictions and theoretical uncertainties coming from the scale variation.MC integration errors are displayed in parentheses.In the last column the K-factor is also shown.quality of the fit.This causes the shift in quark PDFs that can be substantial, while the gluon PDF at small x is completely different between LO and NLO, see e.g.Ref. [91].This is due to the fact that the singular small-x behaviour of the quark to gluon splittings only starts at NLO.In addition, the PDF global analyses involve data from deep-inelastic scattering (DIS) experiments, fixed-target Drell-Yan (DY) data, and collider measurements from the Tevatron and LHC.However, gluon initiated DIS and DY processes, e.g.electroweak boson production processes from ATLAS, CMS and LHCb that are very important for the determination of the gluon PDF, vanish at LO.The graphical representation of the NLO results for the three PDF sets: MSHT20, CT18 and NNPDF3.1 for the pp → t tt t process in the 4ℓ top-quark decay channel is provided in Figure 2. Since the results at the integrated fiducial cross-section level are similar for µ 0 = 2m t and µ 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: µ 0 = 2m t and µ 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 µ R = µ F = ξµ 0 by a factor of ξ in the range ξ ∈ {0.125, ..., 8}.For each case of µ 0 also shown is the variation of µ R with fixed µ F and the variation of µ F with fixed µ 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 µ R (purple line).Changes in the factorisation scale µ F (green line) have a relatively small impact on the integrated fiducial cross section.Lastly, we remark that scale choices with ξ < 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 µ R and µ 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 µ R − µ F plane.In addition to the previous three special cases shown in Figure 3: µ R = µ F = ξµ 0 , µ R = µ 0 , µ F = ξµ 0 and µ R = ξµ 0 , µ F = µ 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 σ (N)LO (µ R, F = ξµ 0 ) to the one calculated using the central value of the scale σ (N)LO (µ R, F = µ 0 ).It should be noted that the region with µ R < µ 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 ξ, as we move toward lower values of µ R and µ F , we can see a significant increase, exceeding a factor of 6, in the σ LO (µ R, F = ξµ 0 )/σ LO (µ R, F = µ 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   where µ 0 = 2m t and µ 0 = E T /4.Both scales are varied simultaneously by a factor of ξ in the following range ξ ∈ {0.125, . . ., 8}.For each case of µ 0 also shown is the variation of µ R with fixed µ F and the variation of µ F with fixed µ R .and a lower limit of 0.46 for the σ NLO (µ R, F = ξµ 0 )/σ NLO (µ R, F = µ 0 ) ratio for both scale settings.Thus, when the whole range of ξ 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 µ R .
We can conclude this part by saying that scanning all potential (µ R , µ F ) pairs in the range 0.5 ≤ µ R /µ 0 ≤ 2 and 0.5 ≤ µ F /µ 0 ≤ 2, without applying the additional requirement of 0.5 ≤ µ R /µ F ≤ 2, would not change the size of the NLO scale uncertainties.Indeed, for the two additional extreme cases (µ R /µ 0 , µ F /µ 0 ) = {(2, 0.5), (0.5, 2)}, 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 → t tt t process in the 4ℓ channel would remain the same even if these two cases were included.
In the next step, we compare our NLO QCD predictions at the integrated fiducial cross-section level, labelled as σ NLO exp , to σ NLO LO dec and σ NLO full .As described earlier, σ NLO full corresponds to the full NWA calculation that takes into account QCD corrections in t tt t production and top-quark decays, while  Γ NLO t is treated as a fixed parameter.On the other hand, the result for σ NLO LO dec utilises NLO QCD corrections in the production stage only.In this case, top-quark decays are described with the LO accuracy using Γ LO t .The NLO QCD results for the three scenarios are summarised in Table 3.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 σ NLO full and σ NLO exp respectively.Secondly, the expansion of Γ NLO t in terms of α s results in the increased size of the scale dependence.Furthermore, we can compare the σ NLO LO dec result with the σ NLO 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 σ NLO LO dec with the full σ NLO 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 α s expansion of Γ t through NLO in QCD with our input parameters and for Decay treatment All results are given for our default setup with the MSHT20 PDF set and for two different scale choices: µ 0 = 2m t and µ 0 = E T /4.Also given are theoretical uncertainties coming from the scale variation.MC integration errors are displayed in parentheses.In the last column the ratio to the σ NLO exp result is provided.
α s (µ R = m t ) = 0.1076705 has the following form 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 → t tt t process in the 4ℓ 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/Γ t term not expanded has the smallest theoretical uncertainties, which are in the range of 15% − 16%.

Differential fiducial cross-section distributions
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 → t tt t process in the 4ℓ 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.
It should be mentioned here that measuring the differential cross-section distributions shown below will be rather difficult even at the High-Luminosity phase of the LHC.The reason behind this is the expected number of events for this process.Indeed, considering only ℓ ± = e ± , µ ± , the expected number of events is N = (15 − 20) per experiment, assuming an integrated luminosity of (3 − 4) ab −1 .On the other hand, when ℓ ± = e ± , µ ± , τ ± we can estimate N = (157 − 209) events for both ATLAS and CMS.We note that the numbers given here correspond to our definition of the fiducial phase-space region.These numbers will change when more inclusive event selection cuts are applied.Nevertheless, proper modelling of differential cross-section distributions is important for example for the many preparatory studies and various analyses that are currently being conducted for this process.In addition, the analysis at the differential cross-section level for the SM pp → t tt tt + X process is a necessary step towards a correct interpretation of possible signals of new physics that may arise in the 4ℓ channel.
We present our findings only for the dynamical scale choice µ R = µ F = µ 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 µ R = µ F = ξµ 0 , where µ 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 → t t + X process, where X = j, γ, Z, W ± , H, jj, γγ, b b, see e.g.Refs.[82,83,[98][99][100][101][102][103][104][105].Similar effects with a fixed scale setting have been examined for the pp → t t process at NLO in QCD in e.g.Ref. [88].
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 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 where µ 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 K-factor, on the other hand, changes from K = 0.98 to 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.TeV.The transverse momenta of the four b-jets, that are ordered according to their p T , are displayed for the dynamical scale setting µ R = µ F = µ 0 = E T /4 and the (N)LO MSHT20 PDF set.The blue (red) curve corresponds to the LO (NLO) result.Also shown the corresponding uncertainty bands resulting from scale variations.The lower panels present the differential 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.
Next, we shift our focus to certain dimensionless observables, highlighting the most noteworthy among those we have examined.In Figure 7  A distinct feature can be noticed for the ∆R b 1 ℓ 1 differential crosssection distribution.In particular, beyond the typical back-to-back configuration, a notable additional peak near ∆R b 1 ℓ 1 ≈ 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 TeV.The transverse momenta of the charged leptons, that are ordered according to their p T , are displayed for the dynamical scale setting µ R = µ F = µ 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 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.
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 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 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 µ R = µ F = µ 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 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.
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 veto T , with p veto T = (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  T cut large differences in the differential factor 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 veto T = 25 GeV we obtain negative values of the NLO cross section for p T ≥ 470 GeV for both observables.This pathological behaviour of fixed-order calculations can be cured by using shower effects or analytic resummation.It should be noted, however, that the jet veto might substantially increase the size of the NLO uncertainty bands.For example, for p veto T = 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. [106][107][108][109][110][111][112] would have to be used.
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, ℓ 1 , ∆R b 1 b 2 and ∆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 µ 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, ℓ 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, ℓ 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 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.
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σ    Results are given for µ R = µ F = µ 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σ NLO full and dσ NLO LO dec to dσ NLO 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.
three approaches together with their corresponding theoretical uncertainties as obtained from the scale dependence.The two middle panels display the ratios of dσ NLO full /dX and dσ NLO LO dec /dX to the dσ NLO 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  Results are given for µ R = µ F = µ 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σ NLO full and dσ NLO LO dec to dσ NLO 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.
their corresponding NLO results for µ 0 = E T /4.For the three observables studied M b 1 b 2 , p T, ℓ 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σ NLO exp /dp T, b 1 b 2 and dσ NLO 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 Figure 12 we additionally display differential cross-section distributions as a function of the difference in azimuthal angle between the two hardest (charged) leptons, ∆ϕ ℓ 1 ℓ 2 , and the two hardest b-jets, ∆ϕ b 1 b 2 .The ∆ϕ ℓ 1 ℓ 2 observable is sensitive to signals of numerous beyond the SM scenarios, where among others new heavy states (Y ) might be produced in the pp → t tY process with Y → t t decays.Indeed, angular distributions of charged leptons reflect spin correlations and can be employed to probe the CP numbers of Y .For X = ∆ϕ ℓ 1 ℓ 2 and X = ∆ϕ b 1 b 2 it is interesting to see the difference between dσ NLO exp /dX and dσ NLO LO dec /dX.In both cases we can observe differences in the range of 8%−11%.Furthermore, the differences between the full results and the expanded ones are up to 15%, still well within the NLO uncertainties that are of the order of 25%.Again, the smallest NLO uncertainties can be seen for the dσ NLO full /dX result.
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 corrections 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/Γ t factor in dσ NLO full /dX significantly reduces the size of the NLO uncertainties for all the observables that we have examined.Thus, the best NLO QCD predictions for the pp → t tt t process in the 4ℓ channel can be obtained by utilising the dσ NLO full /dX results.

Summary and Outlook
In this paper we have calculated NLO QCD corrections to pp → t tt t in the 4ℓ top-quark decay channel for the LHC Run III energy of √ 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 → t tt 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: µ 0 = 2m t and µ 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 Kfactor are driven by the underlining LO cross sections, and specifically by the value of α 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/Γ 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, µ 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 crosssection distributions for which this is the case: In the next step, we have compared the following NLO predictions: dσ NLO exp /dX, dσ NLO full /dX and dσ NLO 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/Γ t term in dσ NLO full /dX, significantly reduces the size of the NLO uncertainties for all the observables that we have scrutinised.Consequently, the best NLO QCD predictions for the process at hand can be obtained by employing σ NLO full and dσ NLO full /dX.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 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 fixedorder NLO QCD predictions.We plan to conduct such a comparison in the near future.
In addition, it would be advantageous to examine the so-called Matrix Element Corrections (MEC) to the Pythia8 parton shower program [113] that can upgrade the accuracy of top-quarks decays to approximately NLO QCD in the shower evolution.Quite recently, a study has been carried out taking into account the MEC for the simplest case of pp → t t production in di-lepton top-quark decay channel [114].Another study has just been performed for the pp → t tW ± process for the two same-sign leptons and jets (2SSℓ) as well as three-lepton (3ℓ) final states [115].Unfortunately, such analyses are not yet available for the pp → t tt t + X process.Although NLO plus parton shower predictions in the Powheg Box and Mc@Nlo framework for pp → t tt t + X production in the ℓ + jets channel are available in the literature [37], the MEC to the top-quark decays are not included there.As this is a much more complicated process with two t t pairs and at least 12 final state particles, it would be interesting to test the validity of such approximate methods in such a complex environment.Another very important issue to look at when using the MEC is spin correlations.Having fixed-order predictions that describe top-quark decays and spin correlations with NLO QCD accuracy, we are in an excellent position to test both aspects.Again, we plan to perform such a comparison in the near future.
The obtained NLO QCD results for the pp → t tt t + X process are especially important for the High-Luminosity phase of the LHC (HL-LHC).The upgrade is currently in progress and ATLAS and CMS experiments are expected to start taking data from the beginning of 2029.The HL-LHC will allow us to study pp → t tt t production in the 4ℓ channel in greater detail.We point out that our conclusions regarding the magnitude of NLO QCD corrections, the significance of higher-order effects in top-quark decays or the size of NLO uncertainties will remain the same if the center-of-mass energy of the LHC is increased from the current value to 14 TeV.Finally, we note that, despite its relatively small cross section, a good theoretical control over the pp → t tt t process in the 4ℓ channel is phenomenologically very relevant.Indeed, the NLO QCD uncertainties that we have estimated for this channel are at the level of 15% − 20%.For comparison, the projected uncertainties for the HL-LHC for the 4ℓ channel with ℓ ± = e ± , µ ± , combining ATLAS and CMS, are in the range of δ = 1/ √ N = 18% − 16%.On the other hand, if ℓ ± = e ± , µ ± , τ ± the later uncertainties decrease to 8% − 7%.We note that these uncertainties are estimated solely on the basis of an integrated luminosity planned for the HL-LHC that amounts to (3 − 4) ab −1 , see e.g.Ref. [116].Thus, at the HL-LHC, among other things, we will gain sensitivity to higher-order effects in top-quark decays in the pp → t tt t + X process.However, even now, pending new results from HL-LHC, the NLO QCD results obtained in this paper are highly relevant for various preparatory studies.Proper modelling of differential cross-section distributions, for example, is an invaluable contribution to the many Machine Learning studies that are currently being conducted for this process.We also believe that the NLO QCD analysis for the SM pp → t tt t + X (background) process at the LHC is a necessary step towards a correct interpretation of possible signals of new physics that may arise in the 4ℓ channel.
Finally, although the focus of this paper has been on providing NLO QCD results for the pp → t tt t process in the 4ℓ channel in the NWA, and the calculation of higher-order corrections to the full 2 → 12 process is a huge task at the moment, we can still estimate the impact of non-resonant and off-shell effects at the integrated and differential fiducial cross-section level with the help of the LO predictions.In fact, at the integrated cross-section level the impact of these omitted contributions is of the order of 1% only, and therefore well within the NWA uncertainty defined by O(Γ t /m t ) ≈ 0.9%, see e.g.Ref. [44].On the other hand, at the differential cross-section level, neglecting triple-, double-, single-and non-resonant top-quark contributions as well as Breit-Wigner top-quark propagators has a significant impact.In general, full off-shell effects are noticeable for dimensionful observables in the following phase-space regions: in high p T tails and in the vicinity of kinematical endpoints.Since the first type of observables is beyond the reach of the LHC, only the second type should be thoroughly investigated.Indeed, we have checked that differential cross-sectional distributions as a function of the minimum invariant mass of the b-jet and the positively charged lepton, M (ℓ + b), the reconstructed invariant mass of the top quark, M (t) as well as the stransverse mass of the top quark and W gauge boson, M T 2 (t) and M T 2 (W ), respectively, are among those observables that are very sensitive to the full off-shell effects in the pp → t tt t process.We refer the reader to Ref. [97] for detailed definitions of these observables.In such cases, it would be beneficial to merge the NLO QCD predictions as calculated in the NWA with the LO full off-shell results.To avoid double counting, the quadruple-resonant top-quark contributions must be properly subtracted from the full off-shell result.Such studies are left for the future.

Figure 1 :
Figure 1: Representative Feynman diagrams for the process pp → t tt t → W + W − W + W − b b b b → ℓ + ν ℓ ℓ − νℓ ℓ + ν ℓ ℓ − νℓ b b b b at the order of O(α 4 s α 8).For simplicity, the decays of top quarks are not shown in the diagrams.All graphs were generated with FeynGame[53].
b b b b, where q = u, d, s, c, b and ℓ ± stands for e ± , µ ± .The leading order (LO) contributions considered in this paper are of the order of O(α 4 s α 8 ) and they consist of Feynman diagrams containing four intermediate top quarks and W bosons along with their corresponding decay products.

Figure 2 :
Figure 2: Integrated fiducial cross sections at NLO in QCD for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.Scale and PDF uncertainties are depicted for the dynamical scale setting µ R = µ F = µ 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 :
Figure 3: Scale dependence of the integrated fiducial cross section at LO and NLO for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.The LO and the NLO MSHT20 PDF sets are employed.Renormalisation and factorisation scales are set to the common value µ R = µ F = µ 0where µ 0 = 2m t and µ 0 = E T /4.Both scales are varied simultaneously by a factor of ξ in the following range ξ ∈ {0.125, . . ., 8}.For each case of µ 0 also shown is the variation of µ R with fixed µ F and the variation of µ F with fixed µ R .

Figure 4 :
Figure 4: Contour plots for the scale dependence in the (µ R , µ F ) plane at LO (upper plots) and NLO (lower plots) for our default PDF set MSHT20 and two scale choices µ 0 = 2m t and µ 0 = E T /4.Results are presented for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.

Figure 5 :
Figure5: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ 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 µ R = µ F = µ 0 = E T /4 and the (N)LO MSHT20 PDF set.The blue (red) curve corresponds to the LO (NLO) result.Also shown the corresponding uncertainty bands resulting from scale variations.The lower panels present the differential 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.
we plot the angular separation between the first two hardest b-jets (∆R b 1 b 2 ) and the first and the third hardest b-jet (∆R b 1 b 3 ).Also displayed is the angular separation between the two softest b-jets (∆R b 3 b 4 ) as well as the hardest b-jet and the hardest lepton (∆R b 1 l 1 ).In the case of ∆R b 1 b 2 and ∆R b 1 b 3 the b-jets are mostly produced in back-to-back configurations.For the ∆R b 3 b 4 observable, on the other hand, a significant number of events can also be found for small values of ∆R b 3 b 4 .

Figure 6 :
Figure 6: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6TeV.The transverse momenta of the charged leptons, that are ordered according to their p T , are displayed for the dynamical scale setting µ R = µ F = µ 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 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.
, we show differential crosssection 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

Figure 7 :
Figure 7: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.The ∆R bb separations between b-jets are displayed for the dynamical scale setting µ R = µ F = µ 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 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 :
Figure 8: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ 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 µ R = µ F = µ 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 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.
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 veto T = (25, 50, 100, 200, 300, 500) GeV.For various values of the p veto

Figure 9 :
Figure 9: Normalised differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ 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 µ R = µ F = µ 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 veto T cut of p veto T = (25, 50, 100, 200, 300, 500) GeV.The lower panels present the ratio to the normalised LO cross section.

Figure 10 :
Figure 10: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.The transverse momentum of the hardest b-jet and charged lepton are presented together with the ∆R bb separation between the hardest and softest b-jets forµ R = µ F = µ 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.
NLO exp /dX, dσ NLO full /dX and dσ NLO LO dec /dX as a function of the observable X, where X = M b 1 b 2 , p T, ℓ 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 -23 -

Figure 11 :
Figure11: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.The upper panels show the absolute NLO QCD predictions for dσ NLO exp /dX (blue curve), dσ NLO full /dX (red curve) and dσ NLO LO dec /dX (green curve), whereX = M b 1 b 2 , p T, ℓ 1 , p T, b 1 and p T, b 1 b 2 .Results are given for µ R = µ F = µ 0 = E T /4and 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σ NLO full and dσ NLO LO dec to dσ NLO 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.

Figure 12 :
Figure12: Differential cross-section distributions for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.The upper panels show the absolute NLO QCD predictions for dσ NLO exp /dX (blue curve), dσ NLO full /dX (red curve) and dσ NLO LO dec /dX (green curve), where X = ∆ϕ ℓ 1 ℓ 2 and ∆ϕ b 1 b 2 .Results are given for µ R = µ F = µ 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σ NLO full and dσ NLO LO dec to dσ NLO 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.

Table 1 :
Integrated fiducial cross sections at LO and NLO in QCD for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6

Table 2 :
Integrated fiducial cross sections at LO and NLO in QCD for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV.The LO and NLO NNPDF3.1 PDF sets are employed for LO predictions with different values of α s

Table 3 :
Integrated fiducial cross sections at NLO in QCD for the pp → t tt t process in the 4ℓ channel at the LHC with √ s = 13.6 TeV for three different scenarios: σ NLO full , σ NLO LO dec and σ NLO exp .