Hadroproduction of four top quarks in the POWHEG BOX

We present a new Monte Carlo event generator for the hadronic production of four top quarks in the POWHEG BOX framework. Besides the dominant next-to-leading order QCD corrections at $\mathcal{O}(\alpha_s^5)$ we also include all subleading electroweak productions channels at leading-order accuracy. We validate our theoretical predictions by comparing to parton-shower matched predictions obtained within the MC@NLO framework for stable top quarks. Furthermore, we investigate in detail the various sources of theoretical uncertainties. Finally, we investigate a single lepton plus jets signature to study for the first time the impact of the electroweak production modes as well as spin-correlation effects at the fiducial level.

We present a new Monte Carlo event generator for the hadronic production of four top quarks in the Powheg Box framework. Besides the dominant next-to-leading order QCD corrections at O(α 5 s ) we also include all subleading electroweak productions channels at leading-order accuracy. We validate our theoretical predictions by comparing to parton-shower matched predictions obtained within the MC@NLO framework for stable top quarks. Furthermore, we investigate in detail the various sources of theoretical uncertainties. Finally, we investigate a single lepton plus jets signature to study for the first time the impact of the electroweak production modes as well as spin-correlation effects at the fiducial level.

I. INTRODUCTION
Since the discovery of the top quark in 1995 by the CDF [1, 2] and DØ [3] experiments at the Tevatron the landscape of top-quark phenomenology has changed dramatically. With the advent of the Large Hadron Collider (LHC) top-quark pairs are so abundantly produced that the top-quark physics program has entered the precision era. Due to its high center-of-mass energy, the LHC offers the unique possibility to produce directly final states involving multiple heavy particles. Thus, in recent years the LHC experiments were able to measure the associated production of top-quark pairs with gauge bosons (W, Z, γ) [4-7] culminating in the discovery of the ttH production process [8,9].
However, there is one rare production process that is of particular interest but for which a discovery has not been claimed so far: the production of four top quarks, pp → tttt, which can only be produced in very high energetic collisions, due to the large threshold of √ s ∼ 700 GeV. Searches for this process have been conducted by both the ATLAS [10, 11] and the CMS [12-16] collaborations already for quite some time. In the recent ATLAS measurement [17] of the cross section an observed (expected) significance with respect to the background-only hypothesis of 4.7 (2.6) standard deviations has been reported.
In recent years, a lot of attention has been devoted to the four top-quark final state as its cross section can be significantly modified by possible physics beyond the Standard Model (BSM) [18][19][20][21][22][23][24][25][26][27]. For instance, in supersymmetric theories the signal is enhanced by cascade decays of gluino-pair production [28][29][30][31]. Due to the large Yukawa coupling y t of the top quark, the production of the pp → tttt process is also of high importance to study modifications of the Higgs sector. For example, the production of heavy Higgs bosons in association with top-quark pairs in two-Higgs-doublet models [32][33][34] can have a big impact on this production rate. Similar effects have been also found in certain top-philic Dark Matter models [35,36]. Furthermore, in some composite-Higgs models the top quark is not a fundamental particle and its compositeness can be studied at lower energies via the impact of higher-dimensional operators of an effective field theory and to which the pp → tttt process is particularly sensitive [37][38][39].
Given the strong sensitivity of the four top-quark production process to modifications of the Standard Model (SM) dynamics, a measurement of this process can also provide stringent constraints on four-fermion operators if interpreted within an effective field theory approach [39][40][41][42][43] or within simplified models [44,45]. However, the former should be taken with care in the presence of onshell intermediate resonances [46]. Finally, the four topquark final state can also be utilized to determine CP properties of the SM Higgs boson [47,48].
In summary, even though the production of four topquark events is one of the rarest and most energetic signatures at the LHC, it represents an extremely versatile process for testing the consistency of the SM and put further constraints on BSM physics. Nonetheless, top quarks are short-living unstable particles that decay predominantly into a W boson and a bottom quark giving rise to a W + W − W + W − bbbb final state. After accounting for the decays of the W bosons the tttt final state yields signatures of unprecedented complexity with at least 12 final state particles, which typically comprise four b jets and additional light jets or multiple leptons. Therefore, in order to harness the full potential of the pp → tttt production process, reliable theoretical predictions for the SM process are necessary.
The dominant next-to-leading order (NLO) QCD corrections at O(α 5 s ) have been calculated in Ref. [49] for the first time and later revised in Ref. [50]. The complete NLO corrections including electroweak (EW) corrections as well as all subleading contributions to the tttt production process at perturbative orders from O(α 5 ) to O(α 5 s ) have been computed for the first time in Ref. [51]. Furthermore, a comparison of parton shower matched computations in the MC@NLO [52,53] framework as provided by MG5 aMC@NLO [54,55] and Sherpa [56,57] has been presented in Ref. [58].
In this paper we present a state-of-the-art Monte Carlo event generator implemented in the Powheg Box V2 framework [59,60]. We consider next-to-leading QCD corrections for the pp → tttt production at O(α 5 s ), while also including all subleading EW production modes at leading-order (LO) accuracy. Furthermore, we model top-quark decays at LO while retaining spin-correlation effects. Having multiple event generators at hand allows to study in more detail various sources of theoretical uncertainties intrinsic to the different approaches of matching fixed-order NLO QCD calculations to parton showers. Our generator is publicly available on the Powheg Box website 1 .
The paper is organized as follows. In section II we provide a brief overview of the technical aspects of this work together with a short review of the structure of higher-order corrections for the four top-quark production at hadron colliders. Next, we discuss in section III the numerical setup for the theoretical predictions shown in section IV. At last, we will give our conclusions in section V.

II. TECHNICAL ASPECTS
In this section we introduce our new event generator and elaborate on the technical aspects of its implementation in the Powheg Box framework. Before doing so, we review the general structure and importance of the various higher-order corrections at the one-loop level.
A. Anatomy of higher-order corrections to pp → tttt We start by reviewing the findings of Ref. [51] to motivate our choice to also include a subset of the most important subleading contributions. At tree-level the pp → tttt process receives contributions at various orders of the strong (α s ) and electromagnetic (α) coupling. Therefore, the perturbative expansion of the leading-order cross section takes the form where Σ LO 1 , Σ LO 3 and Σ LO 5 correspond to squared amplitude contributions and are therefore strictly positive. In Fig. 1 we depict a few sample Feynman diagrams for these contributions. On the contrary, the terms Σ LO 2 and Σ LO 4 1 https://powhegbox.mib.infn.it originate mostly from the interference of Feynman diagrams from different perturbative orders, for the gg and qq initial states, and thus are not necessarily positive. Furthermore, they receive contributions from photon initiated processes.
An atypical feature of the pp → tttt process is the presence of both the qq and the gg initiated production channels in the EW contributions up to O(α 2 s α 2 ). The subleading terms at O(α s α 3 ) and O(α 4 ) contribute only to qq channels. The size of these EW contributions has been studied in Ref. [51], where it has been found that the Σ LO 2 and Σ LO 3 contributions can be as large as −30% and +40% respectively relative to the leading term, Σ LO 1 . The large contributions stem from the tt → tt scattering, for which a representative Feynman diagram is shown in the middle panel of Fig. 1. The exchange of massive bosons, such as the Higgs or Z boson, between non-relativistic top quarks gives rise to Sommerfeld enhancements, which have been thoroughly studied for top-quark pair production in Refs. [61,62]. Contrarily, the terms Σ LO 4 and Σ LO 5 only give rise to corrections below 1%. Even though, some of these, formally subleading, contributions can be very large the total correction amounts to roughly positive 5% − 15% depending on the chosen renormalization scale as there are strong cancellations among them. At the next-to-leading order the situation becomes rather complex and the perturbative expansion of the cross section ranges from O(α 5 ) to O(α 5 s ): As all LO production channels are non-vanishing the classification into pure QCD and pure EW corrections breaks down. The complexity of the computation of NLO corrections in this case increases tremendously because QCD and QED infrared singularities have to be subtracted simultaneously. Only the O(α 5 s ) contribution can be considered a pure QCD correction and represents the bulk of the NLO corrections. The sum of all subleading NLO corrections is below 5% at the inclusive level. However, partial contributions can be sizable as they compensate for the large scale dependence at the leading order. At the differential level NLO corrections to EW channels can become sizable nonetheless, as demonstrated in Ref. [51], especially in the thresold region. In this region, however, we expect the results based on a fixed-order calculation to be less reliable due to the appearance of large logarithms that ultimately need to be resummed. At future colliders operating at higher energies the subleading EW corrections contribute at the same level to inclusive cross sections, while their relative contributions for differential observables can be much more important if these are sensitive to the threshold region.
In summary, besides the NLO QCD corrections at O(α 5 s ), the dominant contribution for the production of four top quarks are the subleading production channels at LO that amount to roughly a +10% correction followed by the remaining NLO corrections. The theoretical uncertainties due to the renormalization and factorization scale dependence are also dominated by the leading NLO QCD corrections. In this work we consider NLO QCD corrections and the subleading EW production channels but not the higher-order corrections to those channels. Furthermore, we neglect photon initiated processes as their contributions are suppressed by the photon parton distribution function. In view of the complexity and size of NLO EW corrections we deem them as dispensable for now. Their inclusion will be beneficial at a later stage as they are expected to reduce the theoretical uncertainties of the subleading contributions here included at LO accuracy.

B. Implementation and validation
We proceed with the implementation details and cross checks that we performed to validate our calculation.
The Powheg Box framework already provides all process independent ingredients for computing NLO QCD corrections. The necessary tree-level and one-loop matrix elements at O(α 4 s ) and O(α 5 s ) are taken from OpenLoops2 [63][64][65], which also provides the necessary spin-and color-correlated born matrix elements to construct the infrared subtraction terms. We are using the OpenLoops2 interface to Powheg Box introduced in Ref. [66]. The matrix elements for the subleading leading-order electroweak contributions for the gg → tttt and qq → tttt processes are instead extracted from MG5 aMC@NLO. Technically, they are included as part of the virtual corrections and, therefore, do not affect the generation of radiation in the Powheg Box framework. Neither are they considered in the generation of the color assignment of the underlying born event. This is expected to have a very limited impact, because the EW contributions do not introduce any new basis elements in the colorflow decomposition of the tree-level matrix elements. Thus, our approach is an approximation that treats the subleading EW contributions as an inclusive correction applied to the leading QCD matrix elements. We have confirmed the validity of this approximation by an explicit leading-order parton-shower matched computation using full tree-level matrix elements but assigned colorflows either according to the QCD or the full matrix elements. We found very good agreement between these two options with differences smaller than 2% at the differential level.
The decays of the top-quarks are included using the algorithm presented in Ref. [67], which retains spin correlations at leading-order accuracy. On the technical side, the implementation is a straightforward adaptation of the algorithm used in Ref. [68], which also allows to include off-shell virtualities for resonant particles. For more details on the decay implementation, we refer the reader to the aforementioned references. All necessary QCD decay-chain matrix elements are also taken from MG5 aMC@NLO. Again, we have checked explicitly at leading-order that the structure of the EW contributions, does not modify the spin-correlations effects in the top decay.
We have validated our implementation by performing several cross checks. For instance, we compared all tree-level and one-loop matrix elements against MG5 aMC@NLO at a few phase space points. At fixed-order we reproduced the inclusive cross sections at NLO QCD as given in Ref. [49] and Ref. [51]. The integrated leading-order cross section including all electroweak production channels, as given by Eq. (1), has been checked against Helac-Phegas [69][70][71]. And finally, the algorithm for the decay of the four top final state has been validated at the differential level against MG5 aMC@NLO in conjunction with MadSpin [72]. To be precise, we compared leptonic observables in the fully leptonic decay channel at LO accuracy, ignoring the EW tree-level contributions for a moment. We find excellent agreement in all cases once spin correlations have been taking into account.

III. COMPUTATIONAL SETUP
We consider the production of four top-quarks at the LHC with a center-of-mass energy of √ s = 13 TeV. For this study we fix the SM parameters to the following values (3) The electromagnetic coupling is derived from the input parameters in the G µ -scheme [73] and given by We calculate the top-quark width at NLO accuracy from all the other input parameters by computing the threebody decay widths Γ(t → ff b) into any light fermionpair ff and a massive b quark. To this end, we employ a numerical routine of the MCFM implementation of Ref. [74]. If not explicitly mentioned otherwise all presented results have been obtained for the NNPDF3.1 [75] 2 parton distribution function (PDF) as provided through the LHAPDF interface [76]. We adopt the dynamical scale choice of Ref. [50] and choose for the renormalization and factorization scale where In order to estimate the theoretical uncertainty due to our particular choice of renormalization and factorization scales, we vary them independently in the range of (1, 1), (1, 2), (2, 1), (2, 2) , (7) and take the envelope as uncertainty estimate. Let us note, that we vary these scales only in the calculation of the hard matrix elements, thus neither the generation of the hardest emission nor the consecutive parton shower evolution are directly affected by these variations. The matching in the Powheg Box framework depends on the two damping parameters h damp and h bornzero that split the real matrix elements into a finite and an infrared singular contribution. While the singular piece is used to resum soft and collinear QCD splittings the finite part, containing only hard emissions, is treated at fixed-order. For more details on these parameters we 2 The LHAPDF ID numbers for the PDF sets are 315200 at LO and 303400 at NLO refer the reader to Refs. [60,[77][78][79]. Based on experience drawn from previous work [68,79,80] these parameters are chosen as and h damp is evaluated on the underlying Born kinematics. We study the impact of our choice by considering the envelope of the following independent variations of these parameters As part of the validation of our computation, we perform a comparison of our results with those obtained with the MG5 aMC@NLO framework that employs the MC@NLO matching to parton showers. We use the same input parameters and renormalization and factorization scales as discussed above. Furthermore, the matching in the MC@NLO scheme depends crucially on the choice of the initial shower scale µ Q . Here we keep the MG5 aMC@NLO default choice of We study the dependence on this scale by varying it by a factor of 2 up and down.
Finally, for all of the following theoretical predictions we use Pythia8 [81,82] (v.8.306) to perform the shower evolution. However, effects from matrix element corrections to the decays, hadronization and multiple interactions are not addressed in this work. The showered events are analyzed using the Rivet [83,84] framework.

IV. PHENOMENOLOGICAL RESULTS
In this section we present our theoretical predictions. We start by investigating the different sources of theoretical uncertainties as well as the impact of subleading EW contributions for the four top-quark production at the inclusive level, i.e. for stable top quarks. Afterwards, we present a sample study for decayed tops in the singlelepton decay channel with a particular focus on the impact of spin correlations contributions at the differential level.

A. Total cross sections
We start our discussion with a detailed investigation of the uncertainty budget of inclusive cross sections. In Table I we show the integrated cross sections at LO accuracy for the leading QCD contribution in the first column and the same also including subleading EW channels in the second column. In the fourth column predictions where NLO QCD corrections on top of the full LO contributions are taken into account are shown. Besides our default choice (NNPDF3.1), we also report integrated cross sections for the MMHT 2014 [85] and CT18 [86] PDF sets. In addition, we also provide the corresponding theoretical uncertainties from scale variations, denoted with δ scale and internal PDF uncertainties denoted with δ PDF . The size of the subleading EW contributions ranges from +5% to +10% depending on the PDF set employed, where for MMHT they have the smallest impact and for CT18 the largest. At leading order the scale uncertainties are large and of the order of 70%, which easily accounts for the differences between predictions based on different PDF sets. Including NLO QCD corrections reduces the scale uncertainties by more than a factor of 2 to 22%. The Kfactor strongly depends on the PDF set employed at the leading order. For instance, if LO PDF sets are used we find +33% corrections in the case of the NNPDF PDF set and only between +4% − 6% corrections for MMHT and CT18. However, when NLO PDF sets are used for both, LO and NLO predictions, we find a rather large K-factor equal to 1.54 stable with respect to the choice of PDF set. Note that the latter K-factor enters the calculation of the Sudakov form factor of the hardest emission. In view of its large size the necessity of including higher order corrections in fixed order as well as parton-shower matched predictions becomes apparent. The PDF uncertainties are much smaller than the scale uncertainties and have been estimated to be ±2% up to ±5% depending on the chosen PDF set. Let us note here, that we follow Ref. [87] and rescale the PDF uncertainties of the CT18 PDF set, which are provided at the 90% confidence level (CL), by a factor of 1/1.645 in order to make them comparable with the other PDF sets that provide uncertainties only at the 68% CL. Finally, we also observe that once NLO QCD corrections are taken into account differences for the central prediction for each PDF are at the 1% level, while there are sizable differences at LO.

B. Inclusive differential distributions
Let us now turn to the discussion of differential distributions at the fully inclusive level, i.e. for stable top quarks. We do not impose any selection cuts on the final state top quarks. Jets are defined via the anti-k T jet algorithm [88] with R = 0.4 as provided by FastJet [89,90] and subject to the following cuts Furthermore, we order top quarks, irrespective of them being particle or anti-particle, according to their transverse momenta in decreasing order. Differential cross section distributions are shown in the following as plots containing three panels. The upper panel shows theoretical predictions at NLO accuracy, the middle panel depicts scale uncertainties computed from an envelope of independent variations of renormalization and factorization scales, and the bottom panel illustrates matching uncertainties estimated by varying the various Powheg Box specific damping parameters or the initial shower scale µ Q in the case of the MC@NLO matching scheme in MG5 aMC@NLO. All theoretical predictions including NLO QCD corrections are labeled with QCD, while we label predictions that include EW Born contributions as QCD+EW. Furthermore, we do not show theoretical uncertainties for the QCD-only predictions obtained with the Powheg Box as we find that they are very similar to those when EW contributions are taken into account. Note that even though the subleading EW channels have separately a sizable scale dependence due to α n s (µ R ) with n ≤ 3, the sum of all EW production modes has a reduced dependence and is smaller than the corresponding scale uncertainty of the dominant NLO QCD corrections at O(α 5 s ). For the transverse momentum of the hardest top quark shown on the left of Fig. 2, we find very good agreement between MG5 aMC@NLO and the QCD only prediction of Powheg Box. Only in the threshold region, where the transverse momenta of all four top quarks become small, differences at the level of 15% are visible. However, this phase space region is also subject to sizable matching uncertainties of the order of 10% for all predictions. In addition, we observe that the EW LO processes give rise to sizable corrections in the threshold region of up to nearly +40%. However, for transverse momenta larger than 300 GeV the EW contributions modify the spectrum by less than 3%. For all theoretical predictions the scale uncertainties due to missing higherorder corrections are the dominant source of uncertainty over the whole range of the distribution. In the case of MG5 aMC@NLO they range from ±20% in the threshold region to ±30% in the tail of the distribution. In Powheg Box, the uncertainties are slightly smaller and their pattern is inverted ranging from ±28% at the beginning of the plotted spectrum to ±15% at the end. We checked explicitly that the increased scale dependence in the threshold region is not due to the inclusion of the EW contributions at leading order accuracy. Therefore, the differences should be rather associated to the different matching frameworks.
Let us now turn to the transverse momentum of the combined third and fourth hardest top quark depicted in the right panel of Fig. 2. Here we find a remarkable agreement between MG5 aMC@NLO and Powheg Box over the whole plotted range if only QCD contributions are taken into account. However, also here we observe sizable corrections of +12% from the EW contributions for small values of the transverse momentum while they are at most of the order of +3% above p T ≈ 250 GeV. As in the previous case, the uncertainties due to missing higher-order corrections dominate over those due to matching for most of the plotted spectrum. For the  Powheg Box prediction the scale uncertainties start out very symmetric with ±20% but grow more and more asymmetric towards the end of the plotted range with estimated uncertainties of +5% and −20%. A similar trend is not visible for MG5 aMC@NLO predictions, as uncertainties are slightly asymmetric from the start with −25% and +16% for small transverse momenta and −30% and +12% uncertainties in the tail of the distribution. For Powheg Box matching uncertainties never exceed ±3%, while for MG5 aMC@NLO they become as large as 10% in the tail of the distribution.
Next we discuss the invariant mass distribution of the four top quarks shown in the left panel of Fig. 3. We observe that the EW contributions are sizable not only at the production threshold M tttt ∼ 4m t but also in the high-energy tail of the distribution. To be specific, we find sizable corrections of the order of +20% in the threshold region that then decrease down to +5% at intermediate value of the invariant mass of around 1.2 TeV. However, the leading order EW production channels give rise to +10% corrections in the very high-energetic tail of the distribution. The matching uncertainties are for all theoretical predictions small and estimated to be below 5%, where in the case of Powheg Box they are the smallest. However, the uncertainties originating from the choice of renormalization and factorization scales are between 20% − 30% over the whole spectrum.
We now turn to the H T distribution, as defined in Eq. (6), depicted in the right panel of Fig. 3. Independently of the employed matching scheme and the inclusion of EW contributions, we find excellent agreement between all predictions above H T ≈ 1.2 TeV with differences below 3%. Only in the first three bins we encounter sizable differences between the various approaches. For instance, the MG5 aMC@NLO prediction is up to 20% smaller at the beginning of spectrum if compared to the corresponding Powheg Box prediction that takes into account the same perturbative corrections. We observe that EW contributions are sizable here as well and yield a +23% enhancement in the low tail. Matching uncertainties are larger than in the previous observables and reach up to ±10% in the case of Powheg Box, and up to ±25% in the case of MG5 aMC@NLO. On the other hand, scale uncertainties are at most ±25% for Powheg Box predictions and in the case of MG5 aMC@NLO they reach ±30%.
Finally, we discuss the transverse momentum and the pseudorapidity distribution of the leading jet as shown in Fig. 4. Contrary to all previous observables, these two distributions are only predicted at leading order accuracy. Nonetheless, they are useful for exploring differences between the POWHEG and the MC@NLO matching schemes. Let us start with the transverse momentum distribution. Here we find large differences between MG5 aMC@NLO and Powheg Box up to 35% at transverse momenta of the order of 150 GeV. Furthermore, MG5 aMC@NLO predicts a softer spectrum in the tail of the distribution. The corresponding theoretical uncertainties are large as well and their bands include predictions from the other event generators, respectively, throughout the majority of the spectrum. In all cases scale uncertainties start around ±25% and grow up to ±55% at the end of the plotted range. In addition, the spectrum exhibits a sizable dependence on the parton-shower matching related parameters. In the case of Powheg Box matching uncertainties are estimated to be of the order of ±6% at the beginning and increase up to ±25% at the end of the distribution. On the other hand, for MG5 aMC@NLO the initial shower scale dependence is rather different, while the beginning and end of the spectrum have modest uncertainties of the order of ±10% they grow as large as 95% in the intermediate range. Thus, the choice of initial shower scale µ Q has a severe impact on the shape of the distribution.
At last, let us turn to the pseudorapidity distribution as depicted on the right of Fig. 4. Also for this observable, we find sizable difference between the various predictions. The impact of the EW production modes is modest and yields rather constant positive corrections at the level of 5% − 7% over the whole spectrum. The largest dif- ferences are found when different matching schemes are employed. The MG5 aMC@NLO predictions is considerably larger by nearly 25% in the central rapidity region as compared to the Powheg Box prediction that also includes only pure QCD corrections. In all cases, the scale uncertainty is the dominant contribution to the theoretical uncertainty and amounts to a constant ±20% over whole plotted range. Similar differences in the modeling of the leading jet between MG5 aMC@NLO and the Powheg Box have been already observed for the pp → ttbb and pp → ttW ± processes, as discussed in Refs. [68,80,91,92].

C. Single lepton plus jets signature
In the following we study a single lepton plus jets signature in order to investigate the impact of spin-correlated top-quark decays and the impact of the leading order EW contributions at the fiducial level. The signature is characterized by the presence of exactly one charged lepton , with = e, µ, at least 4 b jets and at least 6 light jets. The lepton has to fulfill p T ( ) > 15 GeV and |y( )| < 2.5. Jets are formed using the anti-k T jet algorithm with R = 0.4 and a jet is labeled a b jet if at least one of its constituents is a heavy b quark. Light as well as b jets have to pass the p T (j) > 25 GeV and |y(j)| < 2.5 cuts. The definition of the fiducial phase space volume is inspired by Ref. [17].
We show in the following only theoretical predictions obtained with our Powheg Box implementation. We consider three predictions: one prediction that includes both spin correlations in the decay of the top quark as well as the subleading EW channels, and two predictions with either the first or the second improvement switched off. If spin correlations are omitted the decays of top quarks and W bosons are generated via independent 1 → 2 decays. We do not discuss matching uncertainties here anymore as we have seen in the previous section that theoretical uncertainties are dominated by missing higher-order corrections. Moreover, matching uncertainties are expected to be very similar between the various predictions as they are all based on Powheg Box.
For the integrated fiducial cross section we obtain for the three approaches the following results: We observe that EW contributions and spin-correlated decays have opposite effects on the fiducial cross section.
For spin-correlated decays the electroweak production modes increase the cross section by 5% with respect to the QCD only predictions. On the other hand, if spin correlations are ignored then cross sections decreases by 4%. Even though all predictions for integrated cross sections are compatible within the scale uncertainty of roughly ±20%, we will see in the following that this is not necessarily the case at the differential level.
We next inspect the inclusive cross section as a function of the number of light and b jets as shown in Fig. 5. For both observables we recognize essentially the same pattern, independently of the number of jets, with respect to spin correlations and EW production modes as discussed before. However, theoretical uncertainties grow with increasing number of jets. For instance, uncertainties for the cross section as a function of the number of light jets are of the order of ±22% for at least 6 jets and they increase up to ±36% for events with at least 14 jets. Furthermore, nearly 15% of all pp → tttt events are accompanied by at least 9 light jets and only 1% of all events in this signature are associated with at least 12 light jets. On the contrary, the corresponding distribution as a function of the number of b jets falls off more steeply. This is expected, as the dominant source of additional b jets, besides the 4 b jets originating from the topquark decays, are g → bb splittings in the parton-shower evolution. Additionally, initial state b → gb splittings occur as part of the real radiation contribution of the NLO QCD corrections as well as during the parton shower evolution. Nonetheless, these contributions are heavily suppressed by the b-quark parton distribution function. For instance, only 5% of all events have one additional b jet. The theoretical uncertainties as estimated from independent scale variations start again from ±22% and increase up to ±30% for at least 7 b jets. Note that for more realistic estimates of perturbative uncertainties of both cross sections as function of light and b jets in bins beyond the first one shower scale variations should be considered.
Next we discuss hadronic observables such as the transverse momentum of the hardest b jet as shown in the left panel of Fig. 6. First of all we notice that the spectrum is extremely hard. Between the peak of the distribution for transverse momenta around 150 GeV and the tail at 1 TeV the cross section does not drop even 3 orders of magnitude. Similar features have been observed for the pp → ttbb process in Refs. [93,94] for final states with 4 b jets. The scale uncertainties are rather constant over the whole plotted range and are estimated to be between ±25% and ±20%, where the tail of the distribution exhibits smaller uncertainties. Furthermore, we find a significant contribution of up to +10% due to the inclusion of the EW production channels for transverse momenta below 200 GeV.
We now turn to the H T observable as depicted on the right of Fig. 6. At the fiducial level we define H T via where we do not distinguish between light and b jets. For this observable we find an even larger impact due to the leading order EW production channels of up to +25% below 1 TeV. Above, they only contribute a rather constant +3% amount for the rest of the plotted range. Spin correlation effects have also only a mild impact on the tail of the distribution. They soften the spectrum by roughly 5%. Scale uncertainties are estimated to be ±22% at the beginning of the spectrum which then increases up to ±35% at the end of the plotted range. Now we discuss two angular observables that are used for the discrimination of the signal from the background in the experimental analysis of Ref. [17]. To this end, we show in Fig. 7 the minimal distance ∆R among all pairs of b jets, ∆R min bb , as well as among all b jets and the lepton, ∆R min b . For both angular distributions we find strong enhancements in the tails of the distributions due to the EW production of the pp → tttt process. To be specific, we find corrections up to +27% for ∆R min bb and up to +11% for ∆R min b . The scale uncertainties are below ±26% in the whole plotted range and spin correlations only have a very mild impact that lead to minor shape differences.
At last we turn to leptonic observables. We show the transverse momentum of the charged lepton on the left and the missing transverse momentum in the right panel of Fig. 8. For both observable we notice that spin correlations have to be taken into as they have a tremendous impact on the shape of the differential distribution. The transverse momentum distribution of the charged lepton is much harder if correlations are omitted and overshoots the tail of the distribution by nearly 60%. These deviations are not covered by the estimated theoretical uncertainties which are of the order of 20% − 25% over the whole plotted range.
Also for the missing transverse momentum we find large shape differences. However, in this case the distribution becomes softer as compared to the case when spin-correlated top-quark decays are taken into account. While there are significantly more events with low values of p miss T if uncorrelated decays are considered, for p miss T 250 GeV the spectrum is nearly constant −20% smaller than predictions that take spin correlations into account. Furthermore, the scale uncertainties are independent of the treatment of decays and EW production modes of the order of ±20%.

V. CONCLUSIONS
In this article we presented an implementation of the production of four top quarks at hadron colliders in the Powheg Box framework. Besides taking into account the leading NLO QCD corrections at O(α 5 s ) we also include formally subleading EW production channels at LO accuracy. Furthermore, our implementation allows to decay top quarks at LO accuracy retaining spin-correlation effects. This feature made it possible to study the effects  of spin-correlations in top-quark decays in this process for the first time.
We first investigated the impact of leading NLO QCD corrections and subleading EW channels on the total cross sections. We find that QCD corrections contribute at up to slightly over +50% and EW modes at the +5% to +10% level. The inclusion of leading NLO corrections leads to a reduction of scale uncertainties from up to 70% down to at most 22%. Thus the inclusion of NLO QCD corrections as well as the subleading EW channels is essential for reliable predictions of four top production.
We investigated modeling differences for the inclusive pp → tttt production process with stable top-quarks at the 13 TeV LHC by comparing MG5 aMC@NLO and  Powheg Box. We find very good overall agreement between the two frameworks for observables at NLO accuracy, with only minor differences due to the shower evolution in the threshold region. We do also find notable deviations in the hardest light jet spectra, predicted at LO accuracy, which coincides with findings in other production modes of associated top pair production. We also investigated the impact of the EW production channels at the differential level and estimated matching as well as scale uncertainties. We observe that the impact of these subleading channels is generally below 10% but can exceed that near the production threshold. Nonetheless, their inclusion represents a systematic improvement over pure NLO QCD predictions.
Furthermore, we also investigated a single lepton plus jets signature as it is currently employed for pp → tttt cross section measurements and addressed for the first time the size of the subleading EW contributions at the fiducial level. In this case they can reach up to nearly 40% in distributions used for signal/background discriminants in experimental analyses. We also studied the impact of spin-correlated top-quark decays and found that they are essential to obtain a reliable description of leptonic observables. We want to stress that spin-correlation effects are indispensable also for multi-lepton signatures as their impact is a general feature for leptonic observables and signature independent. Finally, we also estimated theoretical uncertainties due to missing higher-order corrections.
The pp → tttt process still awaits its full discovery at the LHC. Once, it is observed it will be instrumental in constraining possible new physics. We are hopeful, that this new tool will help to study more accurately the SM dynamics of the four top quark production process.