Search for $t\bar{t}$ resonances in fully hadronic final states in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

This paper presents a search for new heavy particles decaying into a pair of top quarks using 139 fb$^{-1}$ of proton$-$proton collision data recorded at a centre-of-mass energy of $\sqrt{s}=13$ TeV with the ATLAS detector at the Large Hadron Collider. The search is performed using events consistent with pair production of high-transverse-momentum top quarks and their subsequent decays into the fully hadronic final states. The analysis is optimized for resonances decaying into a $t\bar{t}$ pair with mass above 1.4 TeV, exploiting a dedicated multivariate technique with jet substructure to identify hadronically decaying top quarks using large-radius jets and evaluating the background expectation from data. No significant deviation from the background prediction is observed. Limits are set on the production cross-section times branching fraction for the new $Z'$ boson in a topcolor-assisted-technicolor model. The $Z'$ boson masses below 3.9 and 4.7 TeV are excluded at 95% confidence level for the decay widths of 1% and 3%, respectively.


Introduction
Discovery of new phenomena beyond the Standard Model (SM) is of great importance for high-energy particle physics. Such discovery has the potential to shed light on unexplained observations in nature, e.g. the large difference between the scales of electroweak interactions at O(100) GeV and gravity at O (10 19 ) GeV, the Planck scale. The top quark, the heaviest elementary particle in the SM, could provide a window to this new physics through its large coupling to the scalar sector with a Higgs field. Resonant production of top and anti-top quarks (tt), if observed, strongly indicates the presence of new particles, such as those predicted by topcolor-assisted-technicolor (TC2) [1][2][3], the two-Higgs-doublet model (2HDM) [4] and Randall-Sundrum (RS) models of warped extra dimensions [5,6]. These new particles could appear in the mass spectrum of the tt system (m tt ) as a localized deviation from the SM prediction. This paper presents a search for tt resonances in the TeV mass range with subsequent decay into a fully hadronic final state (tt → W + bW −b with W → qq ), performed using 139 fb −1 of proton-proton (pp) collision data recorded in 2015-2018 at √ s = 13 TeV with the ATLAS detector at the Large Hadron Collider (LHC).
The fully hadronic final state benefits from the largest top-quark decay branching fraction. However, it poses challenges in reconstructing the tt system using hadronic jets, and in discriminating a new physics signal from SM production of multijet processes which have a high production cross-section. Exploiting dedicated top-quark identification techniques, both the ATLAS [7] and CMS [8] experiments demonstrated that the search in the fully hadronic final states can have comparable sensitivity at m tt above ∼ 1 TeV to the analysis in a lepton+jets final state (where one of the W-bosons decays into an electron or muon and a neutrino). This paper extends this further by adopting an advanced top-quark identification method based on a deep neural network. In this high-mass range, the top-quark decay products become close enough for the hadronically decaying top quark to be reconstructed as a single large-radius jet with a characteristic internal substructure. The search focuses on resonances with the intrinsic decay width comparable to or smaller than the detector resolution, by looking for a localized excess over the smoothly falling mass spectra of the reconstructed SM tt candidates. The background spectrum is derived from data by fitting a smoothly falling function to the m tt distributions.
Searches for heavy particles decaying into a tt pair have been performed with the ATLAS and CMS experiments using pp collisions at √ s = 7 TeV [9][10][11][12][13], 8 TeV [14-17] and 13 TeV [7,8,18,19]. Both experiments performed the searches in the fully hadronic and lepton+jets final states using 13 TeV data collected in 2015-2016, while CMS included the dileptonic final states as well. The spin-1 colour-singlet boson in a topcolor-assisted-technicolor model [1,2], used in the previous tt resonance search in the fully hadronic final state in ATLAS [7], continues to be used as a benchmark in the results presented in this paper. This leptophobic Z boson (denoted by Z TC2 ), referred to as Model IV in Ref.
[20], is mainly produced by qq annihilation and decays into first-and third-generation quarks. The model parameters are chosen to maximize the branching fraction of the Z TC2 → tt decay, which reaches 33%, and the intrinsic decay width of the Z TC2 boson divided by its mass m is set to Γ/m = 1% or 3%.1 Among several benchmark signal models, the ATLAS searches excluded masses below 3.1 and 3.0 TeV in the fully hadronic and lepton+jets final states, respectively, for the new Z TC2 boson with Γ/m = 1% [7,18]. The CMS search excluded the Z TC2 boson with Γ/m = 1% up to 3.8 TeV using the combination of all three final states including the dileptonic final state [8].

ATLAS detector
The ATLAS experiment uses a multipurpose, forward-backward symmetric detector2 with nearly 4π solid angle coverage, as described in Refs. [21][22][23]. It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer.
The ID consists of a silicon pixel tracker, a silicon microstrip tracker (SCT) and a transition radiation tracker, all immersed in a 2 T axial magnetic field, and provides charged-particle tracking in the range |η| < 2.5. The EM calorimeter is a lead/liquid-argon (LAr) sampling calorimeter, divided into a barrel section covering the range |η| < 1.475 and two endcap sections covering 1.375 < |η| < 3.2. In the region |η| < 1.8, an additional thin LAr presampler layer is used to correct for energy losses in the material upstream of the calorimeters. The hadronic calorimetry is provided by a steel/scintillator tile sampling calorimeter in the central region (|η| < 1.7) and by a copper/LAr calorimeter in the endcap regions (1.5 < |η| < 3.2). The forward region (3.1 < |η| < 4.9) is instrumented with copper/LAr and tungsten/LAr calorimeter modules optimized for electromagnetic and hadronic measurements, respectively. Surrounding the calorimeters is a muon spectrometer that consists of three air-core superconducting toroidal magnets and tracking chambers, providing precision tracking for muons with |η| < 2.7 and trigger capability for |η| < 2.4.
A two-level trigger system is used to select events for offline analysis [24]. Events are first selected by the level-1 trigger implemented with custom electronics, which uses a subset of the detector information to reduce the event rate to approximately 100 kHz. This is followed by a software-based trigger that reduces the accepted event rate to an average of 1 kHz by refining the level-1 trigger selection.

Data and simulation samples
The search is performed using 139 fb −1 of data recorded from pp collisions at √ s = 13 TeV. The analysis uses data collected during stable beam conditions with the relevant detectors operational. Events are required to have at least one pp interaction vertex with two or more tracks with transverse momentum (p T ) greater than 500 MeV. If more than one vertex is found in an event, the one with the largest p 2 T of associated tracks is chosen as the primary vertex. Simulated signal and background event samples are used to optimize the event selection, to validate the background estimation technique and to perform hypothesis testing of the benchmark Z TC2 signal model. The main backgrounds after applying the event selection criteria (Section 4) are expected to be composed of events from SM tt and multijet production processes. The total number of background events, including events from other minor backgrounds, is estimated directly from data using functional fits to the m tt spectra, as detailed in Section 5. However, simulated samples of SM tt and multijet events are used to establish the background estimation technique.
For SM tt production, the next-to-leading-order (NLO) Monte Carlo (MC) generator P -B v2 [25][26][27] was used with the NNPDF3.0 NLO [28] parton distribution function (PDF) set in the matrix element 2 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of ∆R ≡ (∆η) 2 + (∆φ) 2 .
calculations. The tt production cross-section is scaled to a next-to-next-to-leading-order (NNLO) calculation in QCD including resummation of next-to-next-to-leading logarithmic soft gluon terms with Top++2.0 [29][30][31][32][33][34][35]. Parton showering, hadronization and the underlying event were simulated using P v8.230 [36] with the leading-order (LO) NNPDF2.3 [37] PDF set and the A14 set of tuned parameters [38]. The h damp parameter, which controls the transverse momentum of the first additional parton emission beyond the Born level, was set equal to the top-quark mass [39]. The top-quark kinematics in tt events were corrected to account for electroweak higher-order effects [40]. The generated events were weighted by this correction factor as a function of the flavour and centre-of-mass energy of the initial partons, and of the decay angle of the top quarks in the centre-of-mass frame of the initial partons. The value of the correction factor decreases with increasing m tt from 0.98 at m tt = 0.4 TeV to 0.87 at m tt = 3.5 TeV. Multijet processes were simulated with P v8.186 [36] using the NNPDF2.3 LO PDF set and the A14 set of tuned parameters for the underlying event.
Simulated signal samples of a spin-1 Z TC2 boson decaying into a tt pair were generated using P v8.165 with the NNPDF2.3 LO PDF set and the A14 set of tuned parameters. The production cross-section is scaled to a NLO prediction by multiplying by a factor 1.3 [41] where g (x; µ, σ) represents a Gaussian function with parameters µ and σ, h (x; α CB , n CB , µ CB , σ CB ) represents a Crystal Ball function with parameters α CB , n CB , µ CB and σ CB and 0 < f < 1 is a fractional coefficient. The parameters µ and σ represent the mean value and the width of the Gaussian function, and the α CB and n CB are the threshold and exponent parameters of the CB function, respectively. Thus there are seven shape parameters and an additional normalization parameter. The parameters for the interpolated signal templates are obtained by using either a linear or polynomial interpolation between the values of parameters estimated using a CB + Gauss fit to signal MC m tt distributions. The interpolated signal m tt templates are used at masses of 1875, 2125, 2375, 2625, 2875, 3250, 3500, 3750, 4250, 4500 and 4750 GeV. The signal MC samples are also used to evaluate the acceptance and selection efficiencies for the signals considered in the search.
The E G v1.2.0 program [42] was used in all simulated samples to model the properties of bottom and charm hadron decays. All simulated samples include the effects of multiple pp interactions in the same and neighbouring bunch crossings (pile-up) and were processed through the ATLAS detector simulation [43] based on G 4 [44]. Pile-up effects were emulated by overlaying simulated minimumbias events generated with P v8.186 using the MSTW2008LO PDF set [45] and the A2 set of tuned parameters [46]. The number of overlaid minimum-bias events was adjusted to match the observed data. 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 control samples of data.

Event reconstruction and selection
At high transverse momentum with p T above approximately 350 GeV, the decay products of a hadronically decaying top quark can be reconstructed as a single large-radius (R = 1) jet. Such large-R jets are characterized by the presence of multiple cores associated with the 'subjets' from a b-hadron and a W → qq decay. One of the main backgrounds in this search, multijet production, is dominated by jets from light-flavour quarks and gluons. Typically, such jets have single cores, but would have enhanced contributions of multiple cores associated with large-angle emission when a conventional requirement of large jet mass is applied in the selection of top-quark candidate jets. Therefore, to discriminate genuine top-quark jets from multijet background with a mass around the top-quark mass, an improved analysis technique that takes advantage of more-detailed jet substructure information is desired. An additional challenge in high-mass tt resonance searches comes from the fact that the subjet containing a b-hadron (b-jet) in top-quark decays has, in close vicinity, hadronic activity from W → qq decays. Therefore, the identification of a b-jet inside a large-R jet will need to be optimized to reduce contributions from non-b-hadron decays. The analysis presented here addresses these challenges by adopting advanced top-quark tagging and b-jet identification techniques, as detailed below.

Object reconstruction
Large-R jets are built from three-dimensional topological clusters of energy deposits in the calorimeter, calibrated to the hadronic energy scale with the local cluster weighting (LCW) [47] procedure, using the anti-k t algorithm [48] with a radius parameter R = 1.0. The non-compensating calorimeter response and the energy loss in dead material or due to out-of-cluster leakage of deposited energy are accounted for in the LCW procedure. The reconstructed jets are 'trimmed' [49] to reduce contributions from pile-up and soft interactions. This is performed by reclustering the jet constituents into subjets using the k t algorithm [50-52] with a radius parameter R = 0.2 and discarding subjets with p T less than 5% of the p T of the parent jet [53]. The four-momenta of large-R jets are finally reconstructed from the momentum vectors of the remaining subjets and corrected using simulation [54,55]. The mass of the large-R jet, m J , is calculated by combining the calorimeter energy measurement with the track information from the ID to mitigate the effect on the mass resolution from the limited angular granularity of the calorimeter, as in Ref. [56]. The large-R jets considered in the analysis are selected by requiring p T > 200 GeV, |η| < 2.0 and m J > 50 GeV.
The identification of hadronically decaying top quarks that are reconstructed using large-R jets is performed using a multivariate classification algorithm employed in a deep neutral network [57]. In the kinematic regime of interest in this search, a single large-R jet captures the top-quark decay products, resulting in a characteristic multi-core structure within the jet, in contrast to a typical single-core structure associated with jets in multijet background processes. In order to exploit this characteristic behaviour for the top-quark identification, a multivariate top-tagging classifier has been developed [57]. The tagger uses multiple jet-level discriminants as inputs, e.g. calibrated jet p T and mass, information about the dispersion of the jet constituents such as N-subjettiness [58,59], splitting scales [60] and energy correlation functions [61,62].
The tagger used in this analysis is optimized for top-quark-initiated jets that satisfy the 'contained' criteria defined using the simulation as follows. First, particle-level large-R jets are built from all stable particles (with cτ > 10 mm) at the generator level using a radius parameter R = 1.0 and trimmed in the same way as for reconstructed jets in data. A trimmed particle-level jet is required to match a generator-level top quark within ∆R < 0.75, have a mass above 140 GeV and at least one ghost-associated [63] b-hadron; hereafter the selected particle-level jet is referred to as a truth-contained jet. A reconstructed detector-level large-R jet is then required to satisfy the same ∆R requirement with respect to the matched particle-level jet to ensure that the jet contains the top-quark decay products. The top-tagger calibrated for such large-R jets turns out to be less sensitive to generator differences in the determination of top-quark decay products falling inside a jet, compared with the definition used in Ref. [57]. In this analysis the 'top-tagged' large-R jets are selected using requirements on the classifier corresponding to an efficiency of 80% over the range of generator-level top-quark p T relevant for the considered Z TC2 signal. With the requirement of 80% efficiency, the rejection factor for light-flavour quark and gluon jets is approximately 30 (12) at a top-quark p T of 500 (3000) GeV.
Jets built from charged-particle tracks reconstructed in the ID (track-jets) are used to identify jets containing b-hadrons. The tracks are first selected by requiring them to be associated with the primary vertex and to contain a minimum number of hits in the pixel and SCT detectors, and then requiring them to have p T > 500 MeV and |η| < 2.5. In the topology of highly boosted top quarks [57], the b-hadron decay product is surrounded by other hadronic activity from W → qq decays and this additional contribution from nearby particles degrades the identification of charged-particle tracks from the b-hadron decay. To overcome this, the track-jets are built from the selected tracks using the anti-k t algorithm with R varying as a function of the jet p T . This 'variable-radius' (VR) track-jet [64] has an effective jet radius, R eff , proportional to the inverse of the p T of the jet in the jet-finding procedure: where the ρ-parameter that controls the effective radius is set to ρ = 30 GeV. There are two additional parameters, R min and R max , used to set the minimum and maximum bounds on the jet radius, and these are set to 0.02 and 0.4, respectively [65]. The values of these parameters are determined by examining the efficiency of identifying two b-jets within a large-radius jet associated with a high-p T Higgs boson decaying into a b-quark pair [65].
The VR track-jets are composed of at least two constituent tracks and are required to have p T > 10 GeV and |η| < 2.5. The VR track-jets containing b-hadrons are identified using the 'DL1' algorithm [66]. This algorithm is based on a multivariate classification technique with an artificial deep neural network to combine information from the impact parameters of displaced tracks, reconstructed muons in jets, and topological properties of secondary and tertiary decay vertices reconstructed within the jet. The b-jets are selected in the analysis using the requirement corresponding to an efficiency of 77% for identifying b-jets in simulated SM tt events. This requirement has corresponding rejection factors of 5 and 128 for jets containing c-hadrons and light-flavour jets, respectively. Efficiencies to identify b-jets, c-jets, and light-flavour jets are corrected in the simulation to account for deviations from the efficiencies observed in data [66]. Muons are reconstructed by matching tracks reconstructed in the ID and the muon spectrometer. The muon candidates are required to have p T > 25 GeV and |η| < 2.5, and satisfy the 'medium' quality requirements. The selected muons are calibrated in situ using Z → µµ and J/ψ → µµ decays [70].
Similarly to large-R jets, small-R jets are also built from three-dimensional topological clusters of energy deposits in the calorimeter, but reconstructed with a radius parameter of R = 0.4. Electron and muon candidate tracks are required to be associated with the primary vertex using criteria based on the longitudinal and transverse impact parameters. To avoid the misidentification of jets as electrons and electrons from heavy-flavour hadron decays, the closest small-R jet within a cone of size ∆R y = (∆y) 2 + (∆φ) 2 = 0.2 around a reconstructed electron is removed.4 If an electron is then found within ∆R y = 0.4 of a small-R jet, the electron is removed. If a muon is found within ∆R y = 0.04 + 10 GeV/p µ T of a small-R jet (where p µ T is the muon transverse momentum), the muon is removed if the jet contains at least three tracks, otherwise the jet is removed.

Event selection and categorization
The analysis uses events selected by triggers that require at least one large-R jet with p T > 360-460 GeV, depending on the data-taking period. The large-R jet with the highest p T in the event (referred to as the leading jet) is required to have p T > 500 GeV to ensure a nearly 100% trigger efficiency. Events are further required to contain at least one more large-R jet with p T > 350 GeV to enhance the presence of two jets that can each fully contain the top-quark decay products. Events containing leptons (electrons or muons) are removed to ensure that there is no overlap with events selected by other lepton+jets or dilepton analyses. In addition, events containing small-R jets consistent with originating from detector noise or calorimeter energy deposits associated with non-collision processes are removed. The invariant mass, m JJ , of the two highest-p T large-R jets is required to be larger than 1.4 TeV to avoid a kinematic bias caused by the jet p T requirements. The two leading jets are required to have a separation in azimuthal angle (|∆φ|) larger than 1.6 to ensure a back-to-back topology of top-quark jets. In addition, the rapidity distance between the two leading jets, |∆y JJ |, has to be less than 1.8. This requirement rejects multijet background events with large |∆y JJ | values, which are dominated by processes with t-channel gluon exchanges. Events selected after applying these cuts ('preselection') are further refined by the signal region (SR) selections described below.
In the analysis, the two leading large-R jets in the preselected events are required to be top-tagged. The selected events are further categorized into two orthogonal SRs by the number of b-tagged VR track-jets (n b ) associated with the top-tagged large-R jets by having a separation ∆R < 1.0. Events satisfying n b = 1, i.e. exactly one of the two top-tagged jets is associated with a b-tagged jet, are categorized into 'SR1b'. Events satisfying n b = 2, i.e. each of the two top-tagged jets is associated with a b-tagged jet, are categorized into 'SR2b'. The normalized distributions of the reconstructed m tt (m reco tt ) in the 1b and 2b signal events are shown in Figure 1 for different masses of the Z TC2 boson. The acceptance times efficiency as a function of the invariant mass of a top-quark pair at the generator level, m gen tt , is shown separately for SR1b and SR2b in Figure 2. The acceptance is measured as the fraction of events with two leading truth-contained large-R jets, both satisfying the kinematical requirements on the p T , η, ∆φ and ∆y described above, but not containing any generator-level leptons that satisfy certain kinematic selections.5 The acceptance × efficiency is obtained with respect to the full analysis selections including top-and b-tagging requirements on the two leading large-R jets. The acceptance increases with increasing m gen tt , largely due to the truth-contained requirements. The acceptance × efficiencies have different m gen tt dependence for the two signal regions, mainly caused by the different b-tagging requirements. 4 The rapidity is defined as y = 1 2 ln E+p z E−p z where E is the energy and p z is the longitudinal component of the momentum along the beam direction. 5 The generator-level leptons considered in the acceptance calculation are electrons or muons with p T > 25 GeV and |η| < 2.5. . The acceptance is measured as the fraction of events with two leading truth-contained large-R jets, both satisfying the kinematic requirements, but not containing generator-level electrons or muons, as described in Section 4.2. The acceptance × efficiency is calculated with respect to the full analysis selections including top-and b-tagging requirements on the two leading large-R jets. The m gen tt is calculated from the momenta of top and anti-top quarks at the generator level before final-state radiation. The branching fractions of the tt into all possible final states are included in the acceptance calculation.

Background estimation
The main backgrounds after applying the selection criteria described in Section 4 are expected to arise from SM production of tt pairs and multijet events. The background m reco tt distribution in the signal regions is estimated directly from data by performing a fit with a smoothly falling spectrum. The appropriate functional form of the spectrum is determined using combinations of data and simulated events, as described below.

Background modelling using data and simulation
In order to ensure the correct functional modelling of the background, the functional form is determined using expected background spectra in the signal regions. These are obtained by summing the expected distributions of the SM tt and multijet events.
For the tt background, the simulation samples described in Section 3 are used. For the multijet background, in addition to the multijet simulation samples in Section 3, dedicated control regions (CRs) in data are used to extract the expected m reco tt distributions at low masses. The distributions extracted from the data are combined with those from the simulation sample, such that the combined multijet sample provides more events than expected for the data at high masses. Table 1 shows the CRs defined in data according to whether the leading and subleading large-R jets pass or fail the top-and/or b-tagging requirements. The SRs defined in Section 4.2 correspond to the regions SR1b tbt ¡ b , SR1b t ¡ btb and SR2b tbtb . The first (last) two subscripts correspond to the leading (subleading) large-R jet, and they indicate that the jet is tagged (k) or not ( k) where k = t (b) corresponding to top-tagging (b-tagging). The SR1b signal region consists of the sum of the two regions, SR1b tbt ¡ b and SR1b t ¡ btb . The 'template region' TR t ¡ bt ¡ b with two top-tags and zero b-tags is dominated by multijet events, with a purity of about 97%, and hence it is used to extract the shape of the multijet distribution. The multijet background contributions expected in the SRi (i = 1b tbt ¡ b , 1b t ¡ btb , 2b tbtb ), N exp (SRi), is obtained by multiplying the observed events in the TR, N(TR t ¡ bt ¡ b ), by scale factors R for either one or both leading top-tagged jets to be b-tagged: where and Here the N without subscript represents the observed number of events in the respective region.
The scale factors from TR t ¡ bt ¡ b to SR2b tbtb have two values, R 1 and R 2 , corresponding to the two possible paths: TR t The two resulting N exp (SR2b tbtb ) values agree within 10% with the average of the two, and the average is used as the final prediction for the region SR2b tbtb . The κ t 1 b 2 (κ t 2 b 1 ) factor accounts for the correlation between the top-tagging of the leading (subleading) large-R jet and the b-tagging of the subleading (leading) large-R jet and is obtained as follows: Similarly, the κ t 1 b 1 (κ t 2 b 2 ) factor accounts for the correlation between the top-tagging of the leading (subleading) large-R jet and the b-tagging of the leading (subleading) large-R jet, as follows: These calculations in Eqs. The resulting multijet m reco tt distributions for the SRs are combined with the simulated sample of P multijet events to reduce statistical fluctuation at high mass. This is carried out by using the simulated P sample instead of the data-derived sample in the m reco tt region where the former has higher statistical power than the latter. It turns out that the simulated sample has different m reco tt shape than the data-derived sample and the difference is quantified as a linear fit to the ratio of data-derived to simulated event yields in m reco tt . Therefore, the m reco tt shapes of the simulated sample are corrected in the SRs by applying the fit values in bins of the m reco tt distributions. In order to avoid a dependence on the choice of m reco tt value where the data-derived and simulated multijet samples are combined, the switching point is varied in steps of  GeV between 2260 and 3030 GeV when the background modelling uncertainty is considered (Section 6.2). Table 1: Event categorization used to model the multijet background from data according to whether the leading and subleading large-R jets are top-tagged or b-tagged. If the large-R jet is top-tagged, it is denoted by t, and otherwise by £ t, as indicated in the left column or in the bottom row. Similarly, if the large-R jet is b-tagged, it is denoted by b, and otherwise by ¡ b. The percentages in parentheses show the expected fractions of SM tt events obtained using the tt and multijet simulation samples. Non-tt or non-multijet background events are negligible. The signal regions, SR1b and SR2b, are coloured in red, the template region (TR) in grey and the rest of the control regions A-I in light blue.
Leading large-R jet

Determination of background parameterization
The expected background m reco tt distributions obtained in Section 5.1 are used, along with the simulated SM tt samples, to determine the functional forms of the fits to the SR data. To do this, a set of 1000 background distributions (referred to as S bkg hereafter) is created for each SR by bin-wise variation of the nominal distribution assuming Poisson fluctuations. A typical bin width is chosen using the resolution of the reconstructed m reco tt distribution and is 60 GeV at m reco tt = 1.4 TeV, increasing to 100 (130) GeV at m reco tt = 4 (6) TeV. For each bin of the background m reco tt distribution, the square-root of the bin content is assigned as the uncertainty.
The following parameterized form [71] is used for the fit function: where x = m reco tt / √ s, p 0 is a normalization factor and p 1 through p 4 are free parameters controlling the shape of the m reco tt distributions. The function can be used as a two-or three-parameter function by setting p 3 = p 4 = 0 or p 4 = 0, respectively. The number of shape parameters is then determined by performing fits with individual functions to the S bkg set of background distributions and evaluating the goodness of the individual fits. The evaluation of the goodness of the fit is performed using χ 2 and BumpHunter [72] test statistics. The BumpHunter is a hypothesis-testing tool to look for local excesses or deficits in data relative to the background (obtained from fits in this analysis). After confirming that all the two-, three-and four-parameter functions are qualified, a Wilks' test [73] is repeated over the S bkg set of background distributions to determine how many parameters are needed to describe them. Finally, the functional form with three shape parameters (p 4 = 0) is selected for both SR1b and SR2b. The fit with the three-shape-parameter function is performed to the data to estimate background in each of the SRs.

Systematic uncertainties
For uncertainties in the signal prediction, experimental uncertainties associated with the reconstruction and calibration techniques are considered in the analysis. In addition, an uncertainty associated with the interpolated signal templates and the statistical uncertainty of the simulated samples are taken into account as detailed below. Since the background is directly estimated from data using a functional fit, only the uncertainties associated with the fit method and the statistical uncertainty of the data are taken into account for the background modelling. Each source of experimental systematic uncertainty in the signal is treated as fully correlated across signal regions. The background uncertainties are considered to be uncorrelated as they are estimated from statistically independent samples.

Signal modelling uncertainties
The experimental systematic uncertainties considered for the signal prediction are dominated by uncertainties associated with the large-R jet energy scales (JES), jet mass scale (JMS), jet energy resolutions (JER), jet mass resolution (JMR), and the top-tagging and b-tagging efficiencies.
The uncertainty in the scale of large-R jet p T and mass is evaluated by comparing the ratio of calorimeterbased to track-based measurements in multijet data and simulation [56,74]. The uncertainty in the jet p T resolution [75] is obtained by smearing the MC jet p T using a Gaussian function with a relative width of 0.02 and calculating the resulting change from the unsmeared distribution. This is obtained as a one-sided systematic uncertainty and then symmetrized in the rest of the analysis. The jet mass resolution [56] uncertainty, derived specifically in the context of top-quark-initiated jets, is obtained by smearing the MC jet mass using a Gaussian function with a width corresponding to a 20% relative uncertainty on the jet mass resolution.
The top-tagging efficiency for hadronically decaying top quarks is measured using both data and simulation samples enriched in tt events with one-lepton final states [57]. The difference between the data and simulation efficiencies is taken as a correction factor to the simulated signal yield. The uncertainty in the correction factor, estimated to be 10-15% (per jet) depending on the top-quark jet p T , is propagated through the signal yields in the SRs. An additional top-tagging systematic uncertainty is considered at high p T , beyond a top-quark jet p T of approximately 1 TeV, where the correction factor is not measured due to the limited number of one-lepton tt events in data. This uncertainty is evaluated as the variation of the efficiencies from simulated tt samples with different conditions applied to the G 4 calorimeter shower model and the detector material, and is 1-6% for the top-quark jet p T between 1 and 3 TeV. This additional uncertainty is added in quadrature to the correction factor uncertainty at the highest available p T bin to obtain the top-tagging systematic uncertainty beyond that p T bin. The components of jet p T scale and top-tagging uncertainties associated with the same sources of systematic uncertainties are varied together in the statistical analysis procedure.
The uncertainty in the b-tagging efficiency for b-quark-induced jets is accounted for in a way similar to that for the top-tagging uncertainty. The b-tagging efficiencies in data and simulated events are compared using tt events with final states containing two leptons, and the correction factor is derived as the difference between data and simulation efficiencies [76]. The b-tagging uncertainties in the correction factor are assessed in various kinematic regions, separately for b-jets, c-jets, and light-flavour jets. An additional b-tagging uncertainty at high p T beyond the measured p T range is also taken into account by an extrapolation technique. The uncertainties are then decomposed into various independent components for each of the flavour categories.
The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [77], obtained using the LUCID-2 detector [78] for the primary luminosity measurements. The pile-up modelling uncertainty for simulated events as well as the uncertainties associated with lepton reconstruction and identification are negligible compared with other uncertainties. The uncertainty associated with the interpolated signal is taken into account by comparing expected limits obtained using the simulated and interpolated signal samples at the same masses. The uncertainty is estimated to be at most 5%, and therefore a 5% uncertainty is assigned to the signal event yield, along with the statistical uncertainty due to the limited number of simulated events.
A summary of the post-fit impacts on the 2 and 4 TeV Z TC2 signal yields is given in Table 2. They are expressed as percentage deviations from the nominal yields due to the large-R jet energy scale and resolution, top-tagging and b-tagging uncertainties. The impacts on the signal yields are amplified from the per-jet p T -scale and top-tagging uncertainties due to the presence of two top-quark jets in the events.

Background modelling uncertainties
The uncertainty in the background estimate is primarily associated with the choice of functional form and the fit range used to estimate the background. With a given fit function with a certain number of shape parameters and fit range, there is an ambiguity in the fit result due to intrinsic limitations associated with the chosen fit. Because of this, when a fit is performed to the background-only distribution, the fit creates a systematic difference between the estimated and real backgrounds, producing a 'spurious' signal. This systematic difference, referred to as the spurious-signal uncertainty, is evaluated in the statistical analysis as a bias in the signal estimate obtained from a signal-plus-background fit to the m reco tt distribution constructed under the background-only hypothesis. The signal model used in the signal-plus-background fit is extracted from the simulated Z TC2 samples by performing fits to the reconstructed m reco tt spectra with the sum of a Gaussian function and a Crystal Ball function. The background model is the three-shape-parameter function  The spurious-signal uncertainty is obtained from a signal-plus-background fit performed over the S bkg set of background distributions used to determine the fit functions in Section 5.2. In the fit the switching point for the data-derived and simulated multijet samples is varied in each S bkg set of background distributions (Section 5.1). The spurious signal uncertainty (taken from the average result of the fits) is estimated as a function of the generated signal mass and quantified relative to the statistical uncertainty of the background prediction within ±20% around the signal peak. It decreases monotonically from 80% to 2% between 1.75 and 5 TeV for the 1b signal region. For the 2b signal region, it is approximately constant within 30-40% up to 4 TeV, then increases significantly and reaches 200% at 5 TeV. The increase of the spurious-signal uncertainty at masses above 4 TeV in the 2b signal region is associated with the fact that the background distribution falls more rapidly than that in the 1b signal region, causing the function to be less constrained in the high-mass region and thus producing an increased amount of spurious signal. An additional uncertainty associated with the statistical uncertainty of the expected background in the signal regions is also considered by propagating the impact of the uncertainties in the values of the fit parameters in Eq. (8) to the m reco tt spectra.6

Statistical analysis
The statistical interpretation of the data consists of two steps: a model-independent test for the presence of local deviations from smooth m reco tt spectra in data using BumpHunter, and hypothesis testing of data with the benchmark signal model using the profile likelihood ratio.
The BumpHunter test quantifies the significance of local excesses or deficits in the m reco tt distributions, taking into account the look-elsewhere effect [79, 80] associated with the scanned mass ranges.
The hypothesis testing with the benchmark signal is performed using a binned maximum-likelihood fit to the m reco tt distributions that is based on the expected signal and background yields. The likelihood model is defined as: where P pois (n i |λ i ) is the Poisson probability to observe n i events when λ i events are expected in bin i of the m reco tt distribution, and N (θ) is a series of Gaussian or log-normal distributions for the nuisance parameters, θ, corresponding to the systematic uncertainties related to the signal and background yields in each bin. The λ i is expressed as λ i = µs i (θ) + b i (θ) with µ being the signal strength, defined as a signal cross-section in units of the theoretical prediction, determined by the fit, and s i (θ) and b i (θ) being the expected numbers of signal and background events, respectively. Nuisance parameters are allowed to float in the fit and thus vary the normalization and shape of the signal m reco tt distribution as well as the shape of the background m reco tt distribution. The information about µ is extracted from a likelihood fit to data under the signal-plus-background hypothesis, using a test statistic based on the profile likelihood ratio. The distributions of the test statistic under the signal-plus-background and background-only hypotheses are obtained using the asymptotic formulae [81]. The systematic uncertainties with the largest post-fit impact on µ at m = 2 TeV and m = 4 TeV are the fit parameter uncertainty and the spurious-signal uncertainty for the background and the JES uncertainties for the signal. The level of agreement between the observed data and the background prediction is assessed by computing the local p 0 -value, defined as the probability to observe an excess at least as large as the one observed in data, under the background-only hypothesis. The global p 0 -value is computed by considering the look-elsewhere effect due to multiple testing on the signal mass points. Expected and observed upper limits are set at 95% confidence level (CL) on the cross-section times branching fraction (σ · B) of new particles decaying into a tt pair using the CL s prescription [82]. Cross-checks with the test statistic sampled using pseudo-experiments are performed to test the accuracy of the asymptotic formula in the high-mass region beyond 2 TeV. The σ · B limits from the asymptotic approximation are found to be stronger than those from the pseudo-experiments by at most 20% at masses above 4 TeV. The impact on the mass limit from this approximation is estimated to be below 100 GeV for the Z TC2 signal considered in this analysis.

Results
The observed m reco tt distributions in the two SRs with fits using the three-shape-parameter function are shown in Figure 3. The observed number of data events is 26 964 (8160) in SR1b (SR2b). The BumpHunter tests for the compatibility of the data and the background prediction show that the fit describes the data well for both SRs. The interval with the most significant deviation is 5.44-5.69 TeV for SR1b and 5.44-5.82 TeV for SR2b with the corresponding global p-values of 0.45 and 0.56, respectively. The parameter values of the fit functions determined from the fits to the data are provided in HEPData [83]. The parameter values are consistent with those obtained from fits to the S bkg set of expected background distributions, used in Section 5.2.
With the Z TC2 signal used in this analysis, the minimum local p 0 -value is found to be 0.06 (1.6σ) at Z TC2 mass of 1.88 TeV in the mass range between 1.75 and 5 TeV. In the absence of a significant excess above the background prediction, 95% CL upper limits on σ · B are calculated at each mass value of the Z TC2 signal model. The expected and observed upper limits on the σ · B of Z TC2 → tt are presented in Figure 4. The results from the two SRs are statistically combined to obtain these limits. From the comparison with the σ · B at NLO for the Z TC2 with Γ/m = 1% and 3%, the Z TC2 masses up to 3.9 and 4.7 TeV, respectively, are excluded at 95% CL. For the Z TC2 with Γ/m = 1.2% and the LO σ · B multiplied by 1.3 (scaled to NLO prediction), masses up to 4.1 TeV are excluded at 95% CL. The upper limits on σ · B are provided only up to 5 TeV for the Z TC2 signal mass because of the large spurious-signal uncertainty exceeding 200% at masses beyond 5 TeV, making the limit calculation unreliable at masses larger than ∼ 5.2 TeV. The expected sensitivity of the present analysis is limited by the statistical uncertainty of the background prediction over the full mass range, except at high mass beyond 4.5 TeV where the systematic uncertainty due to the spurious signal dominates the statistical uncertainty.
Observed 95% CL upper limit Expected 95% CL upper limit σ 1 ± Expected 95% CL upper limit σ 2 ± Expected 95% CL upper limit TeV, 139 fb s Figure 4: Observed and expected upper limits on the cross-section times branching fraction of the Z TC2 → tt as a function of the Z TC2 mass. The NLO theory cross-sections times branching fraction for the Z TC2 with Γ/m = 1% and 3% are shown by the red dotted and blue dashed lines, respectively. Shown by the red solid line is a NLO prediction, obtained by multiplying the LO theory cross-section times branching fraction by a factor 1.3 [41], for the Z TC2 with Γ/m = 1.2%.

Conclusion
This paper presents a search for new massive particles decaying into tt in the fully hadronic final state using 139 fb −1 of pp collision data recorded at √ s = 13 TeV with the ATLAS detector at the LHC. The search focuses on the mass range above 1.4 TeV. It uses an improved top-quark tagging based on a multivariate classification algorithm with a deep neural network and b-hadron identification with variable-radius track-jets for highly boosted top quarks. The background is estimated from a fit to data. No significant deviation from the background expectation is observed over the search range. Upper limits are set on the production cross-section times branching fraction for the Z TC2 boson in the topcolor-assisted-technicolor model, resulting in the exclusion of Z TC2 masses up to 3.9 and 4.7 TeV for decay widths of 1% and 3%, respectively. Compared to the previous analysis with the 2015-2016 data set, the improved analysis techniques presented here provide a 65% improvement in the expected cross-section limit at 4 TeV when using the same data set.         [79] L. [83] The Durham High-Energy Physics Database, https://www.hepdata.net.