Search for pair production of up-type vector-like quarks and for four-top-quark events in final states with multiple $b$-jets with the ATLAS detector

A search for pair production of up-type vector-like quarks ($T$) with a significant branching ratio into a top quark and either a Standard Model Higgs boson or a $Z$ boson is presented. The same analysis is also used to search for four-top-quark production in several new physics scenarios. The search is based on a dataset of $pp$ collisions at $\sqrt{s}=13$ TeV recorded in 2015 and 2016 with the ATLAS detector at the CERN Large Hadron Collider and corresponds to an integrated luminosity of 36.1 fb$^{-1}$. Data are analysed in the lepton+jets final state, characterised by an isolated electron or muon with high transverse momentum, large missing transverse momentum and multiple jets, as well as the jets+$E_{T}^{miss}$ final state, characterised by multiple jets and large missing transverse momentum. The search exploits the high multiplicity of jets identified as originating from $b$-quarks, and the presence of boosted, hadronically decaying top quarks and Higgs bosons reconstructed as large-radius jets, characteristic of signal events. No significant excess above the Standard Model expectation is observed, and 95% CL upper limits are set on the production cross sections for the different signal processes considered. These cross-section limits are used to derive lower limits on the mass of a vector-like $T$ quark under several branching ratio hypotheses assuming contributions from $T \rightarrow Wb$, $Zt$, $Ht$ decays. The 95% CL observed lower limits on the $T$ quark mass range between 0.99 TeV and 1.43 TeV for all possible values of the branching ratios into the three decay modes considered, significantly extending the reach beyond that of previous searches. Additionally, upper limits on anomalous four-top-quark production are set in the context of an effective field theory model, as well as in an universal extra dimensions model.


Introduction
The discovery of a new particle consistent with the Standard Model (SM) Higgs boson by the ATLAS [1] and CMS [2] experiments at the Large Hadron Collider (LHC) represents a milestone in high-energy physics. A comprehensive programme of measurements of the Higgs boson properties to unravel its nature is underway at the LHC, so far yielding results compatible with the SM predictions. This makes it more urgent than ever before to provide an explanation for why the electroweak mass scale (and the Higgs boson mass along with it) is so small compared to the Planck scale, a situation known as the hierarchy problem. Naturalness arguments [3] require that quadratic divergences that arise from radiative corrections to the Higgs boson mass are cancelled out by some new mechanism in order to avoid fine-tuning. To that effect, several explanations have been proposed in theories beyond the SM (BSM).
One such solution involves the existence of a new strongly interacting sector, in which the Higgs boson would be a pseudo-Nambu-Goldstone boson [4] of a spontaneously broken global symmetry. One particular realisation of this scenario, referred to as Composite Higgs [5,6], addresses many open questions in the SM, such as the stability of the Higgs boson mass against quantum corrections, and the hierarchy in the mass spectrum of the SM particles, which would be explained by partial compositeness. In this scenario, the top quark would be a mostly composite particle, while all other SM fermions would be mostly elementary. A key prediction is the existence of new fermionic resonances referred to as vector-like quarks, which are also common in many other BSM scenarios. Vector-like quarks are defined as colourtriplet spin-1/2 fermions whose left-and right-handed chiral components have the same transformation properties under the weak-isospin SU(2) gauge group [7,8]. Depending on the model, vector-like quarks are classified as SU(2) singlets, doublets or triplets of flavours T, B, X or Y , in which the first two have the same charge as the SM top and bottom quarks while the vector-like Y and X quarks have charge −4/3e and 5/3e. In addition, in these models, vector-like quarks are expected to couple preferentially to third-generation quarks [7,9] and can have flavour-changing neutral-current decays in addition to the charged-current decays characteristic of chiral quarks. As a result, an up-type T quark can decay not only into a W boson and a b-quark, but also into a Z or Higgs boson and a top quark (T → W b, Zt, and Ht). Similarly, a down-type B quark can decay into a Z or Higgs boson and a b quark, in addition to decaying into a W boson and a top quark (B → Wt, Z b and Hb). Vector-like Y quarks decay exclusively into W b and vector-like X quarks decay exclusively into Wt. To be consistent with the results from precision electroweak measurements a small mass-splitting between vector-like quarks belonging to the same SU(2) multiplet is required, but no requirement is placed on which member of the multiplet is heavier [10]. At the LHC, vector-like quarks with masses below ∼1 TeV would be predominantly produced in pairs via the strong interaction. For higher masses, single production, mediated by the electroweak interaction, may dominate depending on the coupling strength of the interaction between the vector-like quark and the SM quarks.
Another prediction of the Composite Higgs paradigm, as well as other BSM scenarios, such as Randall-Sundrum extra dimensions, is the existence of new heavy vector resonances, which would predominantly couple to the third-generation quarks and thus lead to enhanced four-top-quark production at high energies [11][12][13][14][15]. In particular, the class of models where such vector particles are strongly coupled to the right-handed top quark are much less constrained by precision electroweak measurements than in the case of couplings to the left-handed top quark [16]. In the limit of sufficiently heavy particles, these models can be described via an effective field theory (EFT) involving a four-fermion contact interaction [17]. The corresponding Lagrangian is where t R is the right-handed top quark spinor, γ µ are the Dirac matrices, C 4t is the coupling constant, and Λ is the energy scale above which the effects of direct production of new vector particles must be considered. Anomalous four-top-quark production also arises in Universal Extra Dimensions (UED) models, which involve new heavy particles. For instance, in an UED model with two extra dimensions that are compactified using the geometry of the real projective plane (2UED/RPP) [18], the momenta of particles are discretised along the directions of the extra dimensions. A tier of Kaluza-Klein (KK) towers is labelled by two integers, k and , referred to as "tier (k, )". Within a given tier, the squared masses of the particles are given at leading order by m 2 = k 2 /R 2 4 + 2 /R 2 5 , where πR 4 and πR 5 are the sizes of the two extra dimensions. The model is parameterised by R 4 and R 5 or, alternatively, by m KK = 1/R 4 and ξ = R 4 /R 5 . Four-top-quark production can arise from tier (1,1), where particles from this tier have to be pair produced because of symmetries of the model. Then they chain-decay into the lightest particle of this tier, the heavy photon A (1,1) , by emitting SM particles. The branching ratios of A (1,1) into SM particles are not predicted by the model, although the decay into tt is expected to be dominant [19]. This paper presents a search for TT production with at least one T quark decaying into Ht with H → bb, or into Zt with Z → νν, as well as for anomalous four-top-quark production within an EFT model and within the 2UED/RPP model (see Figure 1). Recent searches for TT production have been performed by the ATLAS [20, 21] and CMS [22,23] collaborations using up to 36.1 fb −1 of pp collisions at √ s = 13 TeV. The most restrictive 95% CL lower limits on the T quark mass obtained are 1.35 TeV and 1.16 TeV, corresponding to branching ratio assumptions of B(T → W b) = 1 and B(T → Zt) = 1, respectively. Previous searches for anomalous tttt production have been performed by the ATLAS Collaboration using the full Run-1 dataset [24,25], where 95% CL limits of |C 4t |/Λ 2 < 6.6 TeV −2 and m K K > 1.1 TeV were obtained in the case of the EFT and the 2UED/RPP models, respectively. A recent search by the CMS Collaboration [26] using 35.9 fb −1 of pp collisions at √ s = 13 TeV has set an upper limit of 41.7 fb on the SM tttt production cross section, about 4.5 times the SM prediction, thus placing some constraints on anomalous production with kinematics like in the SM. This search uses 36.1 fb −1 of data at √ s = 13 TeV recorded in 2015 and 2016 by the ATLAS Collaboration, and it closely follows the strategy developed in Run 1 [25], although it incorporates new ingredients, such as the identification of boosted objects, to substantially enhance sensitivity for heavy resonances. Data are analysed in the lepton+jets final state, characterised by an isolated electron or muon with high transverse momentum, large missing transverse momentum and multiple jets and, for the first time in searches for vector-like quarks, also in the jets+E miss T final state, characterised by multiple jets and large missing transverse momentum.

ATLAS detector
The ATLAS detector [27] at the LHC covers almost the entire solid angle around the collision point,1 and consists of an inner tracking detector surrounded by a thin superconducting solenoid producing a 2 T axial magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large toroid magnet assemblies. The inner detector consists of a high-granularity silicon pixel detector, including the insertable B-layer [28], installed in 2014, and a silicon microstrip tracker, together providing a precise reconstruction of tracks of charged particles in the pseudorapidity range |η| < 2.5, complemented by a transition radiation tracker providing tracking and electron identification information for |η| < 2.0. The calorimeter system covers the pseudorapidity range |η| < 4.9. Within the region |η| < 3.2, electromagnetic (EM) calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) electromagnetic calorimeters, with an additional thin LAr presampler covering |η| < 1.8, to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by a steel/scintillator-tile calorimeter, segmented into three barrel structures within |η| < 1.7, and two copper/LAr hadronic endcap calorimeters. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimised for electromagnetic and hadronic measurements, respectively. The muon spectrometer measures the trajectories of muons with |η| < 2.7 using multiple layers of high-precision tracking chambers located in a toroidal field of approximately 0.5 T and 1 T in the central and endcap regions of ATLAS, respectively. The muon spectrometer is also instrumented with separate trigger chambers covering |η| < 2.4. A two-level trigger system [29], consisting of a hardware-based Level-1 trigger followed by a software-based High-Level Trigger (HLT), is used to reduce the event rate to a maximum of around 1 kHz for offline storage. 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector. The x-axis points from the IP to the centre of the LHC ring, the y-axis points upward, and the z-axis coincides with the axis of the beam pipe. Cylindrical coordinates (r,φ) are used in the transverse plane, φ being the azimuthal angle around the beam pipe. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of ∆R ≡ (∆η) 2 + (∆φ) 2 .

Object reconstruction
Interaction vertices from the proton-proton collisions are reconstructed from at least two tracks with transverse momentum (p T ) larger than 400 MeV that are consistent with originating from the beam collision region in the x-y plane. If more than one primary vertex candidate is found, the candidate whose associated tracks form the largest sum of squared p T [30] is selected as the hard-scatter primary vertex.
is the measured uncertainty in d 0 , and by requiring that their longitudinal impact parameter, z 0 , satisfies |z 0 sin θ| < 0.5 mm. To further reduce the background from non-prompt leptons, photon conversions and hadrons, lepton candidates are also required to be isolated. A lepton isolation criterion is defined by calculating the quantity I R = p trk T , where the sum includes all tracks (excluding the lepton candidate itself) within the cone defined by ∆R < R cut about the direction of the lepton. The value of R cut is the smaller of r min and 10 GeV/p T , where r min is set to 0.2 (0.3) for electron (muon) candidates, and p T is the lepton p T . All lepton candidates must satisfy I R /p T < 0.06. Candidate jets are reconstructed with the anti-k t algorithm [34-36] with a radius parameter R = 0.4 (referred to as "small-R jets"), using topological clusters [37] built from energy deposits in the calorimeters calibrated to the electromagnetic scale. The reconstructed jets are then calibrated to the particle level by the application of a jet energy scale derived from simulation and in situ corrections based on √ s = 13 TeV data [38]. Calibrated jets are required to have p T > 25 GeV and |η| < 2.5. Quality criteria are imposed to reject events that contain any jets arising from non-collision sources or detector noise [39]. To reduce the contamination due to jets originating from pile-up interactions, an additional requirement on the Jet Vertex Tagger (JVT) [40] output is made for jets with p T < 60 GeV and |η| < 2.4.
Jets containing b-hadrons are identified (b-tagged) via an algorithm [41,42] that uses multivariate techniques to combine information about the impact parameters of displaced tracks and the topological properties of secondary and tertiary decay vertices reconstructed within the jet. For each jet, a value for the multivariate b-tagging discriminant is calculated. In this analysis, a jet is considered b-tagged if this value is above the threshold corresponding to an average 77% efficiency to tag a b-quark jet, with a light-jet2 rejection factor of ∼134 and a charm-jet rejection factor of ∼6.2, as determined for jets with p T > 20 GeV and |η| < 2.5 in simulated tt events.
Overlaps between candidate objects are removed sequentially. Firstly, electron candidates that lie within ∆R = 0.01 of a muon candidate are removed to suppress contributions from muon bremsstrahlung. Overlaps between electron and jet candidates are resolved next, and finally, overlaps between remaining jet candidates and muon candidates are removed. Clusters from identified electrons are not excluded during jet reconstruction. In order to avoid double-counting of electrons as jets, the closest jet whose axis is within ∆R = 0.2 of an electron is discarded. If the electron is within ∆R = 0.4 of the axis of any jet after this initial removal, the jet is retained and the electron is removed. The overlap removal procedure between the remaining jet candidates and muon candidates is designed to remove those muons that are likely to have arisen in the decay chain of hadrons and to retain the overlapping jet instead. Jets and muons may also appear in close proximity when the jet results from high-p T muon bremsstrahlung, and in such cases the jet should be removed and the muon retained. Such jets are characterised by having very few matching inner-detector tracks. Selected muons that satisfy ∆R(µ, jet) < 0.04 + 10 GeV/p µ T are rejected if the jet has at least three tracks originating from the primary vertex; otherwise the jet is removed and the muon is kept.
The candidate small-R jets surviving the overlap removal procedure discussed above are used as inputs for further jet reclustering [43] using the anti-k t algorithm with a radius parameter R = 1.0. In this way it is possible to evaluate the uncertainty in the mass of the large-R jets that arises from the uncertainties in the energy scale and resolution of its constituent small-R jets. In order to suppress contributions from pile-up and soft radiation, the reclustered large-R (RCLR) jets are trimmed [44] by removing all small-R (sub)jets within a RCLR jet that have p T below 5% of the p T of the reclustered jet. Due to the pile-up suppression and p T > 25 GeV requirements made on the small-R jets, the average fraction of small-R jets removed by the trimming requirement is less than 1%. The resulting RCLR jets are required to have |η| < 2.0 and are used to identify high-p T hadronically decaying top quark or Higgs boson candidates by making requirements on their transverse momentum, mass, and number of constituents. Hadronically decaying top quark candidates are reconstructed as RCLR jets with p T > 300 GeV, mass larger than 140 GeV, and at least two subjets. Higgs boson candidates are reconstructed as RCLR jets with p T > 200 GeV, a mass between 105 and 140 GeV, and a p T -dependent requirement on the number of subjets (exactly two for p T < 500 GeV, and one or two for p T > 500 GeV). In the following, these are referred to as "top-tagged" and "Higgs-tagged" jets, respectively, while the term "jet" without further qualifiers is used to refer to small-R jets.
The missing transverse momentum ì p miss T (with magnitude E miss T ) is defined as the negative vector sum of the p T of all selected and calibrated objects in the event, including a term to account for energy from soft particles in the event which are not associated with any of the selected objects. This soft term is calculated from inner-detector tracks matched to the selected primary vertex to make it more resilient to contamination from pile-up interactions [45,46].

Data sample and event preselection
This search is based on a dataset of pp collisions at √ s = 13 TeV with 25 ns bunch spacing collected by the ATLAS experiment in 2015 and 2016, corresponding to an integrated luminosity of 36.1 fb −1 . Only events recorded with a single-electron trigger, a single-muon trigger, or a E miss T trigger under stable beam conditions and for which all detector subsystems were operational are considered.
Single-lepton triggers with low p T threshold and lepton isolation requirements are combined in a logical OR with higher-threshold triggers without isolation requirements to give maximum efficiency. For muon triggers, the lowest p T threshold is 20 (26) GeV in 2015 (2016), while the higher p T threshold is 50 GeV in both years. For electrons, triggers with a p T threshold of 24 (26) GeV in 2015 (2016) and isolation requirements are used along with triggers with a 60 GeV threshold and no isolation requirement, and with a 120 (140)  Events satisfying the trigger selection are required to have at least one primary vertex candidate. They are then classified into the "1-lepton" or "0-lepton" channels depending on the multiplicity of selected leptons. Events in the 1-lepton channel are required to satisfy a single-lepton trigger and to have exactly one selected electron or muon that matches, with ∆R < 0.15, the lepton reconstructed by the trigger. In the following, 1-lepton events satisfying either the electron or muon selections are combined and treated as a single analysis channel. Events in the 0-lepton channel are required to satisfy the E miss T trigger and to have no selected leptons. In addition, events in the 1-lepton (0-lepton) channel are required to have ≥5 (≥6) small-R jets. In the following, all selected small-R jets are considered, including those used to build large-R jets. For both channels, backgrounds that do not include b-quark jets are suppressed by requiring at least two b-tagged jets.
Additional requirements are made to suppress the background from multijet production. In the case of the 1-lepton channel, requirements are made on E miss T as well as on the transverse mass of the lepton and and each of the four highest-p T jets. The latter requirement in the 0-lepton channel is very effective in suppressing multijet events, where the large E miss T results from the mismeasurement of a high-p T jet or the presence of neutrinos emitted close to a jet axis.
The above requirements are referred to as the "preselection" and are summarised in Table 1.

Signal and background modelling
Signal and most background processes were modelled using Monte Carlo (MC) simulations. In the simulation, the top quark and SM Higgs boson masses were set to 172.5 GeV and 125 GeV, respectively.
All simulated samples, except those produced with the S [47] event generator, utilised E G v1.2.0 [48] to model the decays of heavy-flavour hadrons. To model the effects of pile-up, events from minimum-bias interactions were generated using the P 8.186 [49] event generator and overlaid onto the simulated hard-scatter events according to the luminosity profile of the recorded data. The generated events were processed through a simulation [50] of the ATLAS detector geometry and response using G 4 [51]. A faster simulation, where the full G 4 simulation of the calorimeter response is replaced by a detailed parameterisation of the shower shapes [52], was adopted for some of the samples used to estimate systematic uncertainties. Simulated events are processed through the same reconstruction software as the data, and corrections are applied so that the object identification efficiencies, energy scales and energy resolutions match those determined from data control samples.

Signal modelling
Samples of simulated TT events were generated with the leading-order (LO) generator4 P 2.2 [8,53] using the NNPDF2.3 LO [54] parton distribution function (PDF) set and passed to P 8.186 [49] for parton showering and fragmentation. The A14 [55] set of optimised parameters for the underlying event (UE) description using the NNPDF2.3 LO PDF set, referred to as the "A14 UE tune", was used. The samples were generated assuming singlet couplings and for heavy-quark masses between 350 GeV and 1.5 TeV in steps of 50 GeV. Additional samples were produced at three mass points (700 GeV, 950 GeV and 1.2 TeV) assuming doublet couplings in order to confirm that, at fixed branching fraction, kinematic differences arising from the different chirality of singlet and doublet couplings have negligible impact on this search. The vector-like quarks were forced to decay with a branching ratio of 1/3 into each of the three modes (W, Z, H). These samples were reweighted using generator-level information to allow results to be interpreted for arbitrary sets of branching ratios that are consistent with the three decay modes summing to unity. The generated samples were normalised to the theoretical cross sections computed using T ++ v2.0 [56] at next-to-next-to-leading order (NNLO) in quantum chromodynamics (QCD), including resummation of next-to-next-to-leading logarithmic (NNLL) soft gluon terms [57][58][59][60][61], and using the MSTW 2008 NNLO [62,63] set of PDFs. The predicted pair-production cross section at √ s = 13 TeV ranges from 24 pb for a vector-like quark mass of 350 GeV to 2.0 fb for a mass of 1.5 TeV, with an uncertainty that increases from 8% to 18% over this mass range. The theoretical uncertainties result from variations of the factorisation and renormalisation scales, as well as from uncertainties in the PDF and α S . The latter two represent the largest contribution to the overall theoretical uncertainty in the cross section and were calculated using the PDF4LHC prescription [64]  Samples of simulated four-top-quark events produced via an EFT and within the 2UED/RPP model were generated at LO with the M 5_aMC@NLO [67] generator (referred to in the following as MG5_aMC; the versions used are 2.2.3 and 1.5.14 for EFT and 2UED/RPP, respectively) and the NNPDF2.3 LO PDF set, interfaced to P 8 (the versions used are 8.205 and 8.186 for EFT and 2UED/RPP, respectively) and the A14 UE tune. The EFT tttt sample was normalised assuming |C 4t |/Λ 2 = 4π TeV −2 , where C 4t denotes the coupling constant and Λ the energy scale of new physics, which yields a cross section of 928 fb computed using MG5_aMC. In the case of the 2UED/RPP model, samples were generated for four different values of m KK (from 1 TeV to 1.8 TeV in steps of 200 GeV) and the B [68] generator was used to decay the pair-produced excitations from tier (1,1) generated by M 5. The corresponding predicted cross section ranges from 343 fb for m KK = 1 TeV to 1.1 fb for m KK = 1.8 TeV.

Background modelling
After the event preselection, the main background is tt production, often in association with jets, denoted by tt+jets in the following. Small contributions arise from single-top-quark, W/Z+jets, multijet and diboson (WW, W Z, Z Z) production, as well as from the associated production of a vector boson V (V = W, Z) or a Higgs boson and a tt pair (ttV and ttH). All backgrounds are estimated using samples of simulated events and initially normalised to their theoretical cross sections, with the exception of the multijet background, which is estimated using data-driven methods. The background prediction is further improved during the statistical analysis by performing a likelihood fit to data using multiple signal-depleted search regions, as discussed in Section 6.
The nominal sample used to model the tt background was generated with the NLO generator P -B v2 [69][70][71][72] using the CT10 PDF set [65]. The P -B model parameter h damp , which controls matrix element to parton shower matching and effectively regulates the high-p T radiation, was set to the top quark mass, a setting that was found to describe the tt system's p T at √ s = 7 TeV [73]. The nominal tt sample was interfaced to P 6.428 [74] with the CTEQ6L PDF set and the Perugia 2012 (P2012) UE tune [75]. Alternative tt simulation samples used to derive systematic uncertainties are described in Section 7.3.
All tt samples were generated inclusively, but events are categorised depending on the flavour content of additional particle jets not originating from the decay of the tt system (see Ref. [76] for details). Events labelled as either tt+≥1b or tt+≥1c are generically referred in the following as tt+HF events, where HF stands for "heavy flavour". A finer categorisation of tt+≥1b events is considered for the purpose of applying further corrections and assigning systematic uncertainties associated with the modelling of heavy-flavour production in different topologies [76]. The remaining events are labelled as tt+lightjets events, including those with no additional jets. In previous studies, better agreement between data and prediction was observed, particularly for the top quark p T distribution, when comparing to NNLO calculations [77]. These small improvements to the modelling are incorporated by reweighting all tt samples to match their top quark p T distribution to that predicted at NNLO accuracy in QCD [78,79]. This correction is not applied to tt+≥1b events, which instead are reweighted to an NLO prediction in the four-flavour (4F) scheme of tt+≥1b including parton showering [80], based on S +O L [47,81] (referred to as S OL in the following) using the CT10 PDF set. This reweighting is performed separately for each of the tt+≥1b categories in such a way that their inter-normalisation and the shape of the relevant kinematic distributions are at NLO accuracy, while preserving the nominal tt+≥1b cross section in P -B +P . The corrections described in this paragraph are applied to the nominal as well as the alternative tt samples.
Samples of single-top-quark events corresponding to the t-channel production mechanism were generated with the P -B v1 [82] generator that uses the 4F scheme for the NLO matrix element calculations and the fixed 4F CT10f4 [65] PDF set. Samples corresponding to the Wtand s-channel production mechanisms were generated with P -B v2 using the CT10 PDF set. Overlaps between the tt and Wt final states are avoided by using the "diagram removal" scheme [83]. The parton shower, hadronisation and the underlying event are modelled using P 6.428 with the CTEQ6L1 PDF set in combination with the P2012 UE tune. The single-top-quark samples were normalised to the approximate NNLO theoretical cross sections [84][85][86].
Samples of W/Z+jets events were generated with the S 2.2 [47] generator. The matrix element was calculated for up to two partons at NLO and up to four partons at LO using C [87] and O L [81]. The matrix element calculation was merged with the S parton shower [88] using the ME+PS@NLO prescription [89]. The PDF set used for the matrix-element calculation is NNPDF3.0NNLO [90] with a dedicated parton shower tuning developed for S . Separate samples were generated for different W/Z+jets categories using filters for a b-jet (W/Z+ ≥1b+jets), a c-jet and no b-jet (W/Z+ ≥1c+jets), and with a veto on band c-jets (W/Z+light-jets), which were combined into the inclusive W/Z+jets samples. Both the W+jets and Z+jets samples were normalised to their respective inclusive NNLO theoretical cross sections in QCD calculated with FEWZ [91].
Samples of WW/W Z/Z Z+jets events were generated with S 2.1.1 using the CT10 PDF set and include processes containing up to four electroweak vertices. The matrix element includes zero additional partons at NLO and up to three partons at LO using the same procedure as for the W/Z+jets samples. The final states simulated require one of the bosons to decay leptonically and the other hadronically. All diboson samples were normalised to their NLO theoretical cross sections provided by S .
Samples of ttV and ttH events were generated with MG5_aMC 2.3.2, using NLO matrix elements and the NNPDF3.0NLO [90] PDF set. Showering was performed using P 8.210 and the A14 UE tune. The ttV samples were normalised to the NLO cross section computed with MG5_aMC. The ttH sample was normalised using the NLO cross section [92][93][94][95][96] and the Higgs boson decay branching ratios calculated using H [97].
The production of four-top-quark events in the SM was simulated by samples generated at LO using MG5_aMC 2.2.2 and the NNPDF2.3 LO PDF set, interfaced to P 8.186 in combination with the A14 UE tune. The sample was normalised to a cross section of 9.2 fb, computed at NLO [67].
The background from multijet production ("multijet background" in the following) in the 1-lepton channel contributes to the selected data sample via several production and misreconstruction mechanisms. In the electron channel, it consists of non-prompt electrons (from semileptonic bor c-hadron decays) as well as misidentified photons (e.g. from a conversion of a photon into an e + e − pair) or jets with a high fraction of their energy deposited in the EM calorimeter. In the muon channel, the multijet background is predominantly from non-prompt muons. The multijet background normalisation and shape are estimated directly from data by using the "matrix method" technique [98], which exploits differences in lepton identification and isolation properties between prompt leptons and leptons that are either non-prompt or result from the misidentification of photons or jets. Further details can be found in Ref. [25]. The main type of multijet background that contributes to the 0-lepton channel are events in which the energy of a high-p T jet is mismeasured, consequently leading to a large missing transverse momentum in the final state. Most of this background is suppressed by selecting events satisfying ∆φ 4j min > 0.4. The remaining multijet background in each search region is estimated from a control region defined with the same selection as the search region, but with the selection on ∆φ 4j min changed to ∆φ 4j min < 0.1. The normalisation of the multijet background is extrapolated from the control region to its corresponding search region by performing an exponential fit to the ∆φ  Figure 2: Comparison of the distribution of (a) the jet multiplicity, and (b) the b-tagged jet multiplicity, between the total background (shaded histogram) and several signal scenarios considered in this search. The selection used in (a) corresponds to events in the 1-lepton channel satisfying the preselection requirements, whereas the selection used in (b) corresponds to events in the 0-lepton channel satisfying the preselection requirements and ≥7 jets. The signals shown correspond to: TT production in the weak-isospin doublet and singlet scenarios, and in the B(T → Zt) = 1 case, assuming m T = 1 TeV; and tttt production within an EFT model.

Search strategy
The searches discussed in this paper primarily target TT production where at least one of the T quarks decays into a Higgs boson and a top quark resulting in the following processes: TT → HtHt, Ht Zt and HtW b.5 For the dominant H → bb decay mode, the final-state signatures in both the 1-lepton and 0-lepton searches are characterised by high jet and b-tagged jet multiplicities, which provide a powerful experimental handle to suppress the background. The presence of high-momentum Z bosons decaying into νν or W bosons decaying leptonically, either to an electron or muon that is not reconstructed or to a hadronically decaying τ-lepton that is identified as a jet, yields high E miss T , which is exploited by the 0-lepton search. Both searches have some sensitivity to TT → Zt Zt and ZtW b, with Z → bb. Possible contributions from pair production of the B or X quarks that would be included, along with the T quark, in a weak-isospin doublet are ignored. Such particles are expected to decay primarily through X, B → Wt [8], and thus not lead to high b-tagged jet multiplicity, which is the primary focus of these searches. High jet and b-tagged jet multiplicities are also characteristic of tttt events (both within the SM and in BSM scenarios); this search is sensitive to these events. The four-top-quark production scenarios considered here do not feature large E miss T , so only the 1-lepton search is used to probe them. No dedicated re-optimisation for tttt events was performed.
In Figure 2(a) the jet multiplicity distribution in the 1-lepton channel after preselection (described in Section 4) is compared between the total background and several signal scenarios, chosen to illustrate Top-tagged jet multiplicity Figure 3: Comparison of the distribution of (a) the Higgs-tagged jet multiplicity and (b) the top-tagged jet multiplicity, between the total background (shaded histogram) and several signal scenarios considered in this search. The selection used in (a) corresponds to events in the 1-lepton channel satisfying the preselection requirements and ≥6 jets, whereas the selection used in (b) corresponds to events in the 0-lepton channel satisfying the preselection requirements and ≥7 jets. The signals shown correspond to: TT production in the weak-isospin doublet and singlet scenarios, and in the B(T → Zt) = 1 case, assuming m T = 1 TeV; and tttt production within an EFT model. differences among various types of signals the search is sensitive to. A similar comparison for the b-tagged jet multiplicity distribution is shown in Figure 2(b) for events in the 0-lepton channel after preselection plus the requirement of ≥7 jets.
Compared to Run 1, the larger centre-of-mass energy in Run 2 provides sensitivity to higher-mass signals, which decay into boosted heavy SM particles (particularly Higgs bosons and top quarks). These potentially give rise to a high multiplicity of large-R jets that capture their decay products (see Section 3). While tt+jets events in the 1-lepton and 0-lepton channels are expected to typically contain one top-tagged jet, the signal events of interest are characterised by higher Higgs-tagged jet and top-tagged jet multiplicities, as illustrated in Figures 3(a) and 3(b). The small fraction (about 5%) of background events with ≥2 top-tagged jets or ≥1 Higgs-tagged jets results from the misidentification of at least one large-R jet where initial-or final-state radiation was responsible for a large fraction of the constituents.
In order to optimise the sensitivity of the searches, the selected events are categorised into different regions depending on the jet multiplicity (5 and ≥6 jets in the 1-lepton channel; 6 and ≥7 jets in the 0-lepton channel), b-tagged jet multiplicity (3 and ≥4 in the 1-lepton channel; 2, 3 and ≥4 in the 0-lepton channel) and Higgs-and top-tagged jet multiplicity (0, 1 and ≥2). In the following, channels with N t top-tagged jets, N H Higgs-tagged jets, n jets, and m b-tagged jets are denoted by "N t t, N H H, nj, mb". Whenever the top/Higgs-tagging requirement is made on the sum N t + N H ≡ N tH , the channel is denoted by "N tH tH, nj, mb". In addition, events in the 0-lepton channel are further categorised into two regions according to the value of m b T, min , the minimum transverse mass of E miss T and any of the three (or two, in events with exactly two b-tagged jets) leading b-tagged jets in the event: m b T, min < 160 GeV (referred to as "LM", standing for "low mass") and m b T, min > 160 GeV (referred to as "HM", standing for "high mass"). This kinematic variable is bounded from above by the top quark mass for semileptonic tt background events, while the signal can have higher values of m b T, min due to the presence of high-p T neutrinos from T → Zt, Z → νν or T → W b, W → ν decays. Although the requirements of a minimum top/Higgs-tagged jet multiplicity reduces the value of m b T, min because of the resulting stronger collimation of the top quark decay products, this variable still provides useful discrimination between signal and tt background, as shown in Figure 4. While the 1-lepton channel only considers regions with exactly 3 or ≥4 b-tagged jets, the 0-lepton channel also includes regions with exactly two b-jets and m b T, min > 160 GeV, to gain sensitivity to TT → Zt Zt decays with at least one Z → νν decay.
To further improve the separation between the TT signal and background, the distinct kinematic features of the signal are exploited. In particular, the large T quark mass results in leptons and jets with large energy in the final state and the effective mass (m eff ), defined as the scalar sum of the transverse momenta of the lepton, the selected jets and the missing transverse momentum, provides a powerful discriminating variable between signal and background. The m eff distribution peaks at approximately 2m T for signal events and at lower values for the tt+jets background. For the same reasons, the various tttt signals from BSM scenarios also populate high values of m eff . An additional selection requirement of m eff > 1 TeV is made in order to minimise the effect of possible mismodelling of the m eff distribution at low values originating from small backgrounds with large systematic uncertainties, such as multijet production. Such a requirement is applied for regions with N t + N H ≤ 1 in the 1-lepton channel, and for all regions in the 0-lepton channel. Since the TT signal is characterised by having at least one top/Higgs-tagged jet and large values of m eff , this minimum requirement on m eff does not decrease the signal efficiency. In Figure 5, the m eff distribution is compared between signal and background for events in signal-rich regions of the 1-lepton and 0-lepton channels. The kinematic requirements in these regions result in a significantly harder m eff spectrum for the background than in regions without top/Higgs-tagged jets, but this variable still shows good discrimination between signal and background. Thus, the m eff distribution is used as the Fraction of events / 250 GeV 1000 1500 2000 2500 3000 3500 Fraction of events / 250 GeV Figure 5: Comparison of the distribution of the effective mass (m eff ), between the total background (shaded histogram) and several signal scenarios considered in this search. The selection used in (a) corresponds to events in the (1t, 1H, ≥6j, ≥4b) region of the 1-lepton channel, whereas the selection used in (b) corresponds to events in the (≥2tH, ≥7j, 2b, HM) region of the 0-lepton channel. The signals shown correspond to: TT production in the weak-isospin doublet and singlet scenarios, and in the B(T → Zt) = 1 case, assuming m T = 1 TeV; and tttt production within an EFT model. The last bin in each distribution contains the overflow.
final discriminating variable in all regions considered in this search.
The regions with ≥6 jets (≥7 jets) are used to perform the search in the 1-lepton (0-lepton) channel (referred to as "search regions"), whereas the regions with exactly 5 jets (6 jets) are used to validate the background modelling in different regimes of event kinematics and heavy-flavour content (referred to as "validation regions"). A total of 12 search regions and 10 validation regions are considered in the 1-lepton channel, whereas 22 search regions and 16 validation regions are considered in the 0-lepton channel, defined in Tables 2 and 3 respectively. In each channel, there are fewer validation regions than signal regions since some validation regions are merged to ensure a minimum of about 10 expected events. The level of possible signal contamination in the validation regions that have high event yields, and are therefore the regions that are most useful to validate the background prediction, depends on the signal scenario considered but is typically well below 10% for a 1 TeV T quark.
The overall rate and composition of the tt+jets background strongly depends on the jet and b-tagged jet multiplicities, as illustrated in Figure 6. The tt+light-jets background is dominant in events with exactly two b-tagged jets, which typically correspond to the two b-quarks from the top quark decays. It also contributes significantly to events with exactly three b-tagged jets, in which typically a charm quark from the hadronic W boson decay is also b-tagged. Contributions from tt+≥1c and tt+≥1b become significant as the b-tagged jet multiplicity increases, with the tt+≥1b background being dominant for events with ≥4 b-tagged jets. The regions with different top/Higgs-tagged jet multiplicities probe different kinematic regimes, both soft (e.g. low-mass T quark) and hard (e.g. high-mass T quark or BSM tttt production). The search regions with the higher multiplicities of top-/Higgs-tagged jets and b-tagged jets in both the 1-lepton and 0-lepton channels, as well as the HM regions in the 0-lepton channel, have the largest  (1 TeV)  T  T  + light-jets  t  t 1c 1-lepton 0-lepton Figure 6: Comparison between the data and the background prediction for the yields in the search regions considered in the 1-lepton and 0-lepton channels, before the fit to data ("Pre-fit"). The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt ". The expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The bottom panel displays the ratio of data to the SM background ("Bkg") prediction. The hashed area represents the total uncertainty of the background, excluding the normalisation uncertainty of the tt+ ≥ 1b background, which is determined via a likelihood fit to data. signal-to-background ratio, and therefore drive the sensitivity of the search. The remaining search regions have significantly lower signal-to-background ratios, but are useful for checking and correcting the tt+jets background prediction and constraining the related systematic uncertainties (see Section 7) through a likelihood fit to data (see Section 8). A summary of the signal-to-background ratio in the different search regions is displayed in Figure 7 for the T quark signal with various decay configurations. A similar fitting strategy was followed in the Run-1 search in the 1-lepton channel [25].
A summary of the observed and expected yields before the fit to data in five of the most sensitive search regions in the 1-lepton and 0-lepton channels can be found in Tables 4 and 5, respectively. The search regions shown in Table 4 for the 1-lepton channel are a selection of some of the regions with the highest S/ √ B ratio (where S and B are the expected signal and background yields, respectively) across several signal benchmark scenarios considered (TT in the B(T → Ht) = 1, T doublet, and T singlet scenarios, in all cases assuming m T = 1 TeV, and tttt within an EFT and the 2UED/RPP models). Similarly, the search regions shown in Table 5 for the 0-lepton channel are a superset of the regions with the highest S/ √ B ratio for different TT signal benchmark scenarios (T doublet, T singlet and B(T → Zt) = 1, also assuming m T = 1 TeV).

0-lepton channel
Search regions (≥7 jets) Validation regions (6 jets)     Table 5: Predicted and observed yields in the 0-lepton channel in five of the most sensitive search regions (depending on the signal scenario) considered. The background prediction is shown before the fit to data. Also shown are the signal predictions for different benchmark scenarios considered. The quoted uncertainties are the sum in quadrature of statistical and systematic uncertainties in the yields, excluding the normalisation uncertainty of the tt+ ≥ 1b background, which is determined via a likelihood fit to data.

Systematic uncertainties
Several sources of systematic uncertainty are considered that affect the normalisation of signal and background and/or the shape of their m eff distributions. Each source of systematic uncertainty is considered to be uncorrelated with the other sources. Correlations for a given systematic uncertainty are maintained across processes and channels, unless explicitly stated otherwise.
The leading sources of systematic uncertainty vary depending on the analysis region considered. For example, the total systematic uncertainty of the background normalisation in the highest-sensitivity search region in the 1-lepton channel (≥0t, ≥2H, ≥6j, ≥4b) is 25%, with the largest contributions originating from uncertainties in tt+HF modelling and flavour tagging efficiencies (b, c, and light). The above uncertainty does not include the uncertainty in the tt+ ≥ 1b normalisation, which is allowed to vary freely in the fit to data. However, as discussed previously, the joint fit to data across the 34 search regions considered in total in the 1-lepton and 0-lepton channels allows the overall background uncertainty to be reduced significantly, e.g., in the case of the search region specified above, down to 10% (including the uncertainty in the tt+ ≥ 1b normalisation). Such a reduction results from the significant constraints that the data places on some systematic uncertainties, as well as the correlations among systematic uncertainties built into the likelihood model.
The following sections describe the systematic uncertainties considered in this analysis.

Luminosity
The uncertainty in the integrated luminosity is 2.1%, affecting the overall normalisation of all processes estimated from the simulation. It is derived, following a methodology similar to that detailed in Ref. [99], from a calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016.

Reconstructed objects
Uncertainties associated with leptons arise from the trigger, reconstruction, identification, and isolation efficiencies, as well as the lepton momentum scale and resolution. These are measured in data using Z → + − and J/ψ → + − events [31, 33]. The combined effect of all these uncertainties results in an overall normalisation uncertainty in signal and background of approximately 1%.
Uncertainties associated with jets arise from the jet energy scale and resolution, and the efficiency to pass the JVT requirement. The largest contribution results from the jet energy scale, whose uncertainty dependence on jet p T and η, jet flavour, and pile-up treatment is split into 21 uncorrelated components that are treated independently in the analysis [38].
The leading uncertainties associated with reconstructed objects in this analysis originate from the modelling of the b-, c-, and light-jet-tagging efficiencies in the simulation, which is corrected to match the efficiencies measured in data control samples [41]. Uncertainties in these corrections include a total of six independent sources affecting b-jets and four independent sources affecting c-jets. Each of these uncertainties has a different dependence on jet p T . Seventeen sources of uncertainty affecting light jets are considered, which depend on jet p T and η. The sources of systematic uncertainty listed above are taken as uncorrelated between b-jets, c-jets, and light-jets. An additional uncertainty is included due to the extrapolation of these corrections to jets with p T beyond the kinematic reach of the data calibration samples used (p T > 300 GeV for band c-jets, and p T > 750 GeV for light-jets); it is taken to be correlated among the three jet flavours. This uncertainty is evaluated in the simulation by comparing the tagging efficiencies while varying e.g. the fraction of tracks with shared hits in the silicon detectors or the fraction of fake tracks resulting from random combinations of hits, both of which typically increase at high p T due to growing track multiplicity and density of hits within the jet. Finally, an uncertainty related to the application of c-jet scale factors to τ-jets is considered, but has a negligible impact in this analysis. The combined effect of these uncertainties results in an uncertainty in the tt background normalisation ranging from 4% to 12% depending on the analysis region. The corresponding uncertainty range for signal is 2-12%, assuming TT production in the weak-isospin doublet scenario and m T = 1 TeV.

Background modelling
A number of sources of systematic uncertainty affecting the modelling of tt+jets are considered. An uncertainty of 6% is assigned to the inclusive tt production cross section [56], including contributions from varying the factorisation and renormalisation scales, and from uncertainties in the PDF, α S , and the top quark mass, all added in quadrature. Since several search regions have a sufficiently large number of events of tt+≥1b background, its normalisation is completely determined by the data during the fit procedure. In the case of the tt+≥1c normalisation, since the fit to the data is unable to precisely determine it and the analysis has very limited sensitivity to its uncertainty, a normalisation uncertainty of 50% is assumed.
Alternative tt samples were generated using P -B interfaced to H ++ 2.7.1 [100] and MG5_aMC 2.2.1 interfaced to H ++ 2.7.1 in order to estimate systematic uncertainties related to the modelling of this background. The effects of initial-and final-state radiation (ISR/FSR) are explored using two alternative P -B +P samples, one with h damp set to 2m t , the renormalisation and factorisation scales set to half the nominal value and using the P2012 radHi UE tune, giving more radiation (referred to as "radHi"), and one with the P2012 radLo UE tune, h damp = m t and the renormalisation and factorisation scales set to twice the nominal value, giving less radiation (referred to as "radLow") [101].
Uncertainties affecting the modelling of tt+≥1b production include shape uncertainties (including intercategory migration effects) associated with the NLO prediction from S OL, which is used for reweighting the nominal P -B +P 6 tt+≥1b prediction. These uncertainties include different scale variations, a different shower-recoil model scheme, and two alternative PDF sets (see Ref. [102] for details), and are significantly smaller than those estimated by comparing different event generators. An uncertainty due to the choice of generator is assessed by comparing the tt+≥1b predictions obtained after reweighting P -B +P 6 to the NLO calculation from S OL and to an equivalent NLO calculation from MG5_aMC+P 8, which differs in the procedure used to match the NLO matrix element calculation and the parton shower (see Section 1.6.8 of Ref. [103]). The uncertainty from the parton shower and hadronisation model is taken from the difference between the MG5_aMC calculation showered with either P 8 or H ++. Additional uncertainties are assessed for the contributions to the tt+≥1b background originating from multiple parton interactions or final-state radiation from top quark decay products, which are not part of the NLO prediction. The latter are assessed via the alternative "radHi" and "radLow" samples, as discussed below. The nominal NLO corrections, as well as their variations used to propagate the theoretical uncertainties in the NLO prediction, are adjusted so that the particle-level cross section of the tt+≥1b background (i.e. prior to reconstruction-level selection requirements) is fixed to the nominal prediction, i.e. effectively only migrations across categories and distortions to the shape of the kinematic distributions are considered.
In the following, uncertainties affecting all tt+jets processes are discussed. Uncertainties associated with the modelling of ISR/FSR are obtained from the comparison of the P -B +P 6 "radHi" and "radLow" samples (see Section 5.2) with the nominal P -B +P 6 sample. An uncertainty associated with the choice of NLO generator is derived by comparing two tt samples, one generated with P -B +H ++ and another generated with MG5_aMC+H ++, and propagating the resulting fractional difference to the nominal P -B +P 6 prediction. An uncertainty due to the choice of parton shower and hadronisation model is derived by comparing events produced by P -B interfaced to P 6 or H ++. Finally, the uncertainty in the modelling of the top quark's p T , affecting only the tt+light-jets and tt+≥1c processes, is evaluated by taking the full difference between applying and not applying the reweighting to match the NNLO prediction. The above uncertainties are taken as uncorrelated between the tt+light-jets, tt+≥1c and tt+≥1b processes. In the case of tt+≥1b, in all instances the various HF categories and the corresponding partonic kinematics for the alternative MC samples are reweighted to match the NLO prediction of S OL so that only effects other than distortions to the inter-normalisation of the various tt+≥1b topologies and their parton-level kinematics are propagated. In the case of tt+light-jets and tt+≥1c the full effect of these uncertainties is propagated. Similarly to the treatment of the NLO corrections and uncertainties associated with tt+≥1b discussed above, in the case of the additional uncertainties derived by comparing alternative tt samples, the overall normalisation of the tt+≥1c and tt+≥1b background at the particle level is fixed to the nominal prediction. In this way, only migrations across categories and distortions to the shape of the kinematic distributions are considered. In order to maintain the inclusive tt cross section, the tt+light-jets background is adjusted accordingly.
Uncertainties affecting the modelling of the single-top-quark background include a +5%/−4% uncertainty in the total cross section estimated as a weighted average of the theoretical uncertainties in t-, Wtand s-channel production [84][85][86]. Additional uncertainties associated with the modelling of ISR/FSR are assessed by comparing the nominal samples with alternative samples where generator parameters were varied (i.e. "radHi" and "radLow"). For the tand Wt-channel processes, an uncertainty due to the choice of parton shower and hadronisation model is derived by comparing events produced by P -B interfaced to P 6 or H ++. These uncertainties are treated as fully correlated among single-top production processes, but uncorrelated with the corresponding uncertainty in the tt+jets background. The sum in quadrature of the above uncertainties on the single top normalisation at the preselection level is 20% in the 1-lepton channel and 20%(25%) in LM(HM) regions of the 0-lepton channel, respectively. An additional systematic uncertainty on Wt-channel production concerning the separation between tt and Wt at NLO [104] is assessed by comparing the nominal sample, which uses the so-called "diagram subtraction" scheme, with an alternative sample using the "diagram removal" scheme. This uncertainty, which is taken to be single-sided, has a strong shape dependence and affects the Wt normalisation by about −50% in the 1-lepton channel and LM regions of the 0-lepton channel, and by about −75% in HM regions of the 0-lepton channel. Due to the small size of the simulated samples, and hence limited statistical precision, these uncertainties cannot be reliably estimated in each analysis region and so their estimates at the preselection level are used instead. They are treated as uncorrelated across regions with different top-tagged jet and Higgs-tagged jet multiplicities and between the 1-lepton and 0-lepton channels.
Uncertainties affecting the normalisation of the V+jets background are estimated for the sum of W+jets and Z+jets, and separately for V+light-jets, V+≥1c+jets, and V+≥1b+jets subprocesses. The total normalisation uncertainty of V+jets processes is estimated by comparing the data and total background prediction in the different analysis regions considered, but requiring exactly 0 b-tagged jets. Agreement between data and predicted background in these modified regions, which are dominated by V+light-jets, is found to be within approximately 30%. This bound is taken to be the normalisation uncertainty, correlated across all V+jets subprocesses. Since S 2.2 has been found to underestimate V+heavy-flavour by about a factor of 1.3 [105], additional 30% normalisation uncertainties are assumed for V+≥1c+jets and V+≥1b+jets subprocesses, considered uncorrelated between them. These uncertainties are treated as uncorrelated across regions with different top-/Higgs-tagged jet multiplicities and between the 1-lepton and 0-lepton channels.
Uncertainties in the diboson background normalisation include 5% from the NLO theory cross sections [106], as well as an additional 24% normalisation uncertainty added in quadrature for each additional inclusive jet-multiplicity bin, based on a comparison among different algorithms for merging LO matrix elements and parton showers [107]. Therefore, normalisation uncertainties of 5% ⊕ √ 3 × 24% = 42% and 5% ⊕ √ 4 × 24% = 48% are assigned for events with exactly 5 jets and ≥6 jets, respectively (this assumes that two jets come from the W/Z decay, as in WW/W Z → ν j j). Recent comparisons between data and S 2.1.1 for W Z(→ ν )+ ≥4 jets show agreement within the experimental uncertainty of approximately 40% [108], which further justifies the above uncertainty. This uncertainty is taken to be uncorrelated across regions with different top-/Higgs-tagged jet multiplicities and between the 1-lepton and 0-lepton channels Uncertainties in the ttV and ttH cross sections are 15% and +10%/−13%, respectively, from the uncertainties in their respective NLO theoretical cross sections [109][110][111]. Finally, an uncertainty of 30% is estimated for the NLO prediction of the SM tttt cross section [67]. Since no additional modelling uncertainties are taken into account for these backgrounds, and the 1-lepton and 0-lepton channels cover different kinematic phase spaces, the above uncertainties in the ttV, ttH, and SM tttt cross sections are taken to be uncorrelated between the two channels.
Uncertainties in the data-driven multijet background estimate receive contributions from the limited sample size in data, particularly at high jet and b-tag multiplicities, as well as from the uncertainty in the misidentified-lepton rate, measured in different control regions (e.g. selected with a requirement on either the maximum E miss T or m W T ). Based on the comparisons between data and total prediction in multijet-rich selections, the normalisation uncertainties assumed for this background are 50% (100%) for electrons with |η cluster | ≤ 1 (|η cluster | > 1), and 50% for muons, taken to be uncorrelated across regions with different top-/Higgs-tagged jet multiplicities and between events containing electrons and events containing muons. In the case of the 0-lepton channel, the normalisation uncertainty assigned to the multijet background is 100%. No explicit shape uncertainty is assigned since the large statistical uncertainties associated with the multijet background prediction, which are uncorrelated between bins in the final discriminant distribution, are assumed to effectively cover possible shape uncertainties.

Statistical analysis
For each search, the m eff distributions across all regions considered are jointly analysed to test for the presence of a signal predicted by the benchmark scenarios. The statistical analysis uses a binned likelihood function L(µ, θ) constructed as a product of Poisson probability terms over all bins considered in the search. This function depends on the signal-strength parameter µ, which multiplies the predicted production cross section for signal, and θ, a set of nuisance parameters that encode the effect of systematic uncertainties in the signal and background expectations. Therefore, the expected total number of events in a given bin depends on µ and θ. With the exception of the parameter that controls the normalisation of the tt+≥1b background, all other nuisance parameters are implemented in the likelihood function as Gaussian or log-normal constraints. The above-mentioned tt+≥1b normalisation factor is a free parameter of the fit.
For a given value of µ, the nuisance parameters θ allow variations of the expectations for signal and background according to the corresponding systematic uncertainties, and their fitted values result in the deviations from the nominal expectations that globally provide the best fit to the data. This procedure allows a reduction of the impact of systematic uncertainties on the search sensitivity by taking advantage of the highly populated background-dominated regions included in the likelihood fit. To verify the improved background prediction, fits under the background-only hypothesis are performed, and differences between the data and the post-fit background prediction are checked using kinematic variables other than the ones used in the fit. The m eff distributions in validation regions not used in the fit are also checked. Statistical uncertainties in each bin of the predicted m eff distributions due to the limited size of the simulated samples are taken into account by dedicated parameters in the fit.
The test statistic q µ is defined as the profile likelihood ratio: q µ = −2 ln(L(µ,θ µ )/L(μ,θ)), whereμ andθ are the values of the parameters that maximise the likelihood function (subject to the constraint 0 ≤μ ≤ µ), andθ µ are the values of the nuisance parameters that maximise the likelihood function for a given value of µ. The test statistic q µ is evaluated with the RooFit package [112,113]. A related statistic is used to determine the probability that the observed data are compatible with the background-only hypothesis (i.e. the discovery test) by setting µ = 0 in the profile likelihood ratio and leavingμ unconstrained: q 0 = −2 ln(L(0,θ 0 )/L(μ,θ)). The p-value (referred to as p 0 ) representing the probability of the data being compatible with the background-only hypothesis is estimated by integrating the distribution of q 0 from background-only pseudo-experiments, approximated using the asymptotic formulae given in Refs. [114,115], above the observed value of q 0 . Some model dependence exists in the estimation of the p 0 , as a given signal scenario needs to be assumed in the calculation of the denominator of q µ , even if the overall signal normalisation is left floating and fitted to data. The observed p 0 is checked for each explored signal scenario. Upper limits on the signal production cross section for each of the signal scenarios considered are derived by using q µ in the CL s method [116,117]. For a given signal scenario, values of the production cross section (parameterised by µ) yielding CL s < 0.05, where CL s is computed using the asymptotic approximation [114,115], are excluded at ≥ 95% CL.

Results
This section presents the results obtained from searches in the 1-lepton and 0-lepton channels, as well as their combination, following the statistical analysis discussed in Section 8.

Likelihood fits to data
A binned likelihood fit under the background-only hypothesis is performed on the m eff distributions in all search regions considered. In this section, the results of the simultaneous likelihood fit to the search regions in the 1-lepton and 0-lepton channels are discussed. This combined fit is used to obtain results on TT production. In this combination, all common systematic uncertainties are considered fully correlated between the 1-lepton and 0-lepton channels, with the exception of those affecting non-tt backgrounds. To  obtain the results in the individual channels, separate fits are performed. In general, good agreement is found among the fitted nuisance parameters in the individual and combined fits.
A comparison of the distribution of observed and expected yields in the search regions in the 1-lepton and 0-lepton channels after the combined fit is shown in Figure 8 (see Figure 6 for the results before the combined fit). The post-fit yields in five of the most sensitive search regions in the 1-lepton and 0-lepton channels can be found in Tables 6 and 7, respectively. For the same search regions, the corresponding m eff distributions, both before and after the fit to data, are shown in Figures 9-13. The binning used for the m eff distributions in the different search regions represents a compromise between preserving enough discrimination between the background and the different signal hypotheses considered, and keeping the statistical uncertainty on the background prediction per bin well below 30%. While some of the systematic uncertainties from individual sources described in Section 7 vary across the m eff spectrum, the total pre-fit uncertainty is largely independent of m eff . The large number of events in the signal-depleted regions, together with their different background compositions, and the assumptions of the fit model, constrain the combined effect of the sources of systematic uncertainty. As a result, an improved background prediction is obtained with significantly reduced uncertainty, not only in the signal-depleted channels, but also in the signal-rich channels such as (≥0t, ≥2H, ≥6j, ≥4b) in the 1-lepton channel. In the combined fit, the channels with three b-tagged jets are effectively used to constrain the leading uncertainties affecting the tt+light-jets background prediction, while the channels with ≥4 b-tagged jets are sensitive to the uncertainties affecting the tt+HF background prediction. In particular, one of the main corrections determined in the fit is a scale factor that multiplies the tt+≥1b normalisation by 0.90 ± 0.23 relative to the nominal prediction.6 In addition, the nuisance parameter controlling the tt+≥1c normalisation is adjusted to scale this background by a factor of 1.3 ± 0.4 relative to its nominal prediction. The fit results in better agreement between data and prediction in the channels with ≥3 b-tagged jets, where the tt+HF background dominates. Detailed studies were performed to verify the stability of the fit against variations in the treatment of the systematic uncertainties affecting the tt+HF background (e.g. by decorrelating normalisation and shape uncertainties between different tt+≥1b categories, or by scaling the tt+≥1b and tt+≥1c backgrounds by a common factor), finding in all instances a robust post-fit background prediction. Furthermore, the impact on the background-only fit of injecting a TT signal (with m T = 1 TeV) in the doublet configuration was confirmed to be negligible. Although there is no single nuisance parameter directly responsible for the normalisation of tt+light-jets background, the yields for this contribution within each region are affected by systematic uncertainties in the tt modelling and the jet flavour tagging, and thus are changed after the fit.
A comparison of the distribution of observed and expected yields in all validation regions considered, before and after the combined fit in the search regions, is shown in Figure 14. Agreement between data and prediction in normalisation and shape of the m eff distribution for these regions, which are not used in the fit, is generally improved after the fit, giving confidence in the overall procedure. To increase the background yields and strengthen the validation of the fit strategy, comparisons between data and background prediction, before and after the fit, are performed for more-inclusive event selections. As an example, the distributions of two kinematic variables used to define the search strategy can be found in Figures 15 and 16. They display respectively the Higgs-tagged jet multiplicity in the 1-lepton channel, after requiring at least 6 jets and 3 b-jets, and the distribution of the m b T, min variable in the 0-lepton channel for events containing at least 7 jets and 2 b-jets, together with at least one top/Higgs-tagged jet. Although these variables are not directly used in the fit, a good description of the data by the post-fit background prediction is observed, which further validates the fitting procedure. The result of the background-only fit to data is used for the background prediction in the computation of the limits presented in the following subsections.  Table 6: Predicted and observed yields in the 1-lepton channel in five of the most sensitive search regions (depending on the signal scenario) considered. The multijet background is considered negligible in these regions and thus not shown. The background prediction is shown after the combined fit to data in the 0-lepton and 1-lepton channels under the background-only hypothesis. The quoted uncertainties are the sum in quadrature of statistical and systematic uncertainties in the yields, computed taking into account correlations among nuisance parameters and among processes.  Table 7: Predicted and observed yields in the 0-lepton channel in five of the most sensitive search regions (depending on the signal scenario) considered. The background prediction is shown after the combined fit to data in the 0-lepton and 1-lepton channels under the background-only hypothesis. The quoted uncertainties are the sum in quadrature of statistical and systematic uncertainties in the yields, computed taking into account correlations among nuisance parameters and among processes.
(d) Figure 9: Comparison between the data and prediction for the m eff distribution in some of the most sensitive search regions in the 1-lepton channel, before and after performing the combined fit to data in the 0-lepton and 1-lepton channels ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. Shown are the (≥2t, 0-1H, ≥6j, 3b) region (a) pre-fit and (b) post-fit, and the (1t, 0H, ≥6j, ≥4b) region (c) pre-fit and (d) post-fit. In the pre-fit figures the expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.
(d) Figure 10: Comparison between the data and prediction for the m eff distribution in some of the most sensitive search regions in the 1-lepton channel, before and after performing the combined fit to data in the 0-lepton and 1-lepton channels ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. Shown are the (1t, 1H, ≥6j, ≥4b) region (a) pre-fit and (b) post-fit, and the (≥2t, 0-1H, ≥6j, ≥4b) region (c) pre-fit and (d) post-fit. In the pre-fit figures the expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The blue triangles indicate points that are outside the vertical range of the figure. The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.
(d) Figure 11: Comparison between the data and prediction for the m eff distribution in some of the most sensitive search regions, before and after performing the combined fit to data in the 0-lepton and 1-lepton channels ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. Shown are the (≥2H, ≥6j, ≥4b) region in the 1-lepton channel (a) pre-fit and (b) post-fit, and the (≥2tH, ≥7j, 2b, HM) region in the 0-lepton channel (c) pre-fit and (d) post-fit. In the pre-fit figures the expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The small contributions from ttV, ttH, single top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The blue triangles indicate points that are outside the vertical range of the figure. The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.
(d) Figure 12: Comparison between the data and prediction for the m eff distribution in some of the most sensitive search regions in the 0-lepton channel, before and after performing the combined fit to data in the 0-lepton and 1-lepton channels ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. Shown are the (1t, 1H, ≥7j, 3b, HM) region (a) pre-fit and (b) post-fit, and the (≥2t, 0-1H, ≥7j, 3b, HM) region (c) pre-fit and (d) post-fit. In the pre-fit figures the expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.
(d) Figure 13: Comparison between the data and prediction for the m eff distribution in some of the most sensitive search regions in the 0-lepton channel, before and after performing the combined fit to data in the 0-lepton and 1-lepton channels ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. Shown are the (1t, 0H, ≥7j, ≥4b, HM) region (a) pre-fit and (b) post-fit, and the (≥2tH, ≥7j, ≥4b) region (c) pre-fit and (d) post-fit. In the pre-fit figures the expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.
1-lepton 0-lepton (b) Figure 14: Comparison between the data and background prediction for the yields in each of the validation regions considered in the 1-lepton and 0-lepton channels (a) before the fit ("Pre-fit") and (b) after the fit ("Post-fit"). The fit is performed on the data in 1-lepton and 0-lepton channels under the background-only hypothesis considering only the search regions. In the pre-fit figure the expected TT signal (solid red) corresponding to m T = 1 TeV in the T doublet scenario is also shown, added on top of the background prediction. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The bottom panels display the ratios of data to the total background prediction ("Bkg"). The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.
Higgs-tagged jets multiplicity   Figure 15: Comparison between the data and prediction for the Higgs-tagged jet multiplicity in the 1-lepton channel after preselection plus the requirement of ≥6 jets and ≥3 b-tagged jets, (a) before and (b) after performing the combined fit of the m eff spectrum to data in the 0-lepton and 1-lepton channels search regions ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.  T, min ) in the (≥1tH, ≥7j, ≥2b) region of the 0-lepton channel (a) before and (b) after performing the combined fit of the m eff spectrum to data in the 0-lepton and 1-lepton channels search regions ("Pre-fit" and "Post-fit", respectively) under the background-only hypothesis. The small contributions from ttV, ttH, single-top, W/Z+jets, diboson, and multijet backgrounds are combined into a single background source referred to as "Non-tt". The last bin in all figures contains the overflow. The bottom panels display the ratios of data to the total background prediction ("Bkg"). The hashed area represents the total uncertainty of the background. In the case of the pre-fit background uncertainty, the normalisation uncertainty of the tt+ ≥ 1b background is not included.

Limits on vector-like quark pair production
No significant excess above the SM expectation is found in any of the search regions. Upper limits at 95% CL on the TT production cross section are set in several benchmark scenarios as a function of the T quark mass m T and are compared to the theoretical prediction from Top++. The resulting lower limits on m T correspond to the central value of the theoretical cross section. The scenarios considered involve different assumptions about the decay branching ratios. The search in the 1-lepton (0-lepton) channel is particularly sensitive to the benchmark scenario of B(T → Ht) = 1 (B(T → Zt) = 1). In contrast, both the 1-lepton and the 0-lepton searches have comparable sensitivity to the weak-isospin doublet and singlet scenarios, and thus their combination represents an improvement of 60-70 GeV on the expected T quark mass exclusion over the most sensitive individual search. The limits corresponding to the weak-isospin doublet and singlet scenarios obtained for the combination of the 1-lepton and 0-lepton searches are shown in Figure 17. A summary of the observed and expected lower limits on the T quark mass in the different benchmark scenarios for the individual 1-lepton and 0-lepton searches, as well as their combination, is given in Table 8. As can be seen, the observed mass limits for the 1-lepton search are above the expected limits in all benchmark scenarios. Detailed studies on the statistical model found no sources of systematic bias and showed that the results are consistent with downward statistical fluctuations in data in some of the highest m eff bins in three search regions: (1t, 1H, ≥6j, ≥4b), (≥2t, 0-1H, ≥6j, 3b), and (≥0t, ≥2H, ≥6j, ≥4b). Several other regions with similar event kinematics and background composition to these three search regions show good agreement between data and expectations. In particular, additional regions with larger event yields were constructed to test this agreement by merging signal regions in certain categories, but retaining similar multiplicities of b-tagged jets or boosted objects as the original signal regions.   ) σ 1 ± Theory (NNLO prediction 95% CL observed limit 95% CL expected limit σ 1 ± 95% CL expected limit σ 2 ± 95% CL expected limit ) σ 1 ± Theory (NNLO prediction 95% CL observed limit 95% CL expected limit σ 1 ± 95% CL expected limit σ 2 ± 95% CL expected limit Zt) + B(T → Ht) = 1. To probe this branching ratio plane, the signal samples are reweighted by the ratio of the desired branching ratio to the original branching ratio in P , and the complete analysis is repeated. Owing to the complementarity of the 1-lepton and 0-lepton searches in probing the branching ratio plane, their combination represents a significant improvement over the individual results, as illustrated in Figure 18. In this case, the observed lower limits on the T quark mass range between 0.99 TeV and 1.43 TeV depending on the values of the branching ratios into the three decay modes. In particular, a vector-like T quark with mass below 0.99 TeV is excluded for any values of the branching ratios into the three decay modes. The corresponding range of expected lower limits is between 0.91 TeV and 1.34 TeV. Figure 19 presents the corresponding observed and expected T quark mass limits in the plane of B(T → Ht) versus B(T → W b), obtained by linear interpolation of the calculated CL s versus m T .

Limits on four-top-quark production
The 1-lepton search is used to set limits on BSM four-top-quark production by considering different signal benchmark scenarios (see Section 5.1 for details). In the case of tttt production via an EFT model with a four-top-quark contact interaction, the observed (expected) 95% CL upper limit on the production cross section is 16 fb (31 +12 −9 fb). The upper limit on the production cross section can be translated into an observed (expected) limit on the free parameter of the model |C 4t |/Λ 2 < 1.6 TeV −2 (2.3 ± 0.4 TeV −2 ). In the context of the 2UED/RPP model, the observed and expected upper limits on the production cross section times branching ratio are shown in Figure 20 as a function of m KK for the symmetric case Exp. 1-lepton limit Exp. 0-lepton limit Exp. combination limit Obs. combination limit Exp. 1-lepton limit Exp. 0-lepton limit Exp. combination limit Obs. combination limit Figure 18: Observed (red filled area) and expected (red dashed line) 95% CL exclusion in the plane of B(T → W b) versus B(T → Ht), for different values of the vector-like T quark mass for the combination of the 1-lepton and 0lepton searches. In the figure, the branching ratio is denoted "BR". The background estimate used in the computation of the limits is the result obtained from the background-only fit to data. Also shown are the expected exclusions by the individual searches, which can be compared to that obtained through their combination. The grey (light shaded) area corresponds to the unphysical region where the sum of branching ratios exceeds unity, or is smaller than zero. The default branching ratio values from the P event generator for the weak-isospin singlet and doublet cases are shown as plain circle and star symbols, respectively. (ξ = R 4 /R 5 = 1), assuming production by tier (1,1) alone. The comparison to the LO theoretical cross section translates into an observed (expected) 95% CL limit on m K K of 1.8 TeV (1.7 TeV). Theory (LO prediction) 95% CL observed limit 95% CL expected limit σ 1 ± 95% CL expected limit σ 2 ± 95% CL expected limit   (1,1) in the symmetric case (ξ = R 4 /R 5 = 1). The background estimate used in the computation of the limits is the result obtained from the background-only fit to data. The surrounding shaded bands correspond to ±1 and ±2 standard deviations around the expected limit. The thin red line shows the theoretical prediction, computed at LO in QCD, for the production cross section of four-top-quark events by tier (1,1) assuming B(A (1,1) → tt) = 1, where the heavy photon A (1,1) is the lightest particle of this tier.

Conclusion
A search for pair production of up-type vector-like quarks (T) with significant branching ratio into a top quark and either a Standard Model Higgs boson or a Z boson is presented. The same analysis is also used to search for four-top-quark production, in several new physics scenarios. The search is based on pp collisions at Additionally, upper limits on the four-top-quark production cross section are set in several new physics scenarios. In the case of tttt production from a contact interaction in an EFT model, the observed (expected) 95% CL upper limit on the production cross section is 16 fb (31 +12 −9 fb). In the context of a 2UED/RPP model, 95% CL observed (expected) lower limits on m K K of 1.8 TeV (1.7 TeV) are derived.
Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF; BSF, GIF and Minerva, Israel; BRF, Norway; CERCA Programme Generalitat de Catalunya, Generalitat Valenciana, Spain; the Royal Society and Leverhulme Trust, United Kingdom.