Search for top-philic heavy resonances in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

A search for the associated production of a heavy resonance with a top-quark or a top-antitop-quark pair, and decaying into a $t\bar{t}$ pair is presented. The search uses the data recorded by the ATLAS detector in $pp$ collisions at $\sqrt{s}= 13$ TeV at the Large Hadron Collider during the years 2015-2018, corresponding to an integrated luminosity of 139 fb$^{-1}$. Events containing exactly one electron or muon are selected. The two hadronically decaying top quarks from the resonance decay are reconstructed using jets clustered with a large radius parameter of $R=1$. The invariant mass spectrum of the two top quark candidates is used to search for a resonance signal in the range of 1.0 TeV to 3.2 TeV. The presence of a signal is examined using an approach with minimal model dependence followed by a model-dependent interpretation. No significant excess is observed over the background expectation. Upper limits on the production cross section times branching ratio at 95% confidence level are provided for a heavy $Z^\prime$ boson based on a simplified model, for $Z^\prime$ mass between 1.0 TeV and 3.0 TeV. The observed (expected) limits range from 21 (14) fb to 119 (86) fb depending on the choice of model parameters.


Introduction
The discovery of a new particle consistent with the Standard Model (SM) Higgs boson by the ATLAS [1] and CMS [2] collaborations was a major milestone in high-energy physics.However, the underlying nature of electroweak symmetry breaking remains unknown.Naturalness arguments suggest that the large quantum corrections to the Higgs boson mass are cancelled by some new mechanism to prevent an excessive level of fine-tuning.Such a mechanism has been proposed in several theories of Beyond the SM (BSM) physics.The large Yukawa coupling of the top quark to the Higgs boson motivates various top-quark-based resonance searches.In many BSM theories, such as composite Higgs scenarios [3-6], new 'top-philic' vector resonances are predicted.These resonances couple more strongly to the top quark than to the light quarks, such that all the other couplings except for the one with top quarks can be neglected.Typical t t resonance searches target new resonances produced through q q annihilation, assuming sizeable couplings to light quarks [7-10].Top-philic resonances, on the other hand, require different production modes, such as the production of the heavy resonance in association with a top-quark or a t t pair, resulting in three or four top quarks in the final state.Figure 1 shows tree-level diagrams for the production of a top-philic resonance Z .
This analysis searches for a new top-philic heavy resonance above a mass of 1.0 TeV.The ATLAS and CMS collaborations have previously published measurements of SM four-top-quark (t tt t) production [11][12][13][14][15]. with the first observation of this process reported by both experiments recently Fig. 1 Examples of tree-level Feynman diagrams for Z production in association with (a) t t, (b) t j (where j refers to any light quark), and (c) t W .The Z generation modes are derived from top quark final states produced via (a) strong, (b) electroweak, and (c) mixed QCD and electroweak interactions [16,17].Three-top-quark production, with a much smaller cross section of O(1 fb) predicted by the SM [18][19][20], has not been measured.In resonant BSM models, the two top quarks from the resonance decay, referred to as the 'resonance top quarks', are expected to be highly boosted.The other top quarks are referred to as the 'spectator top quarks' and are expected to have lower momenta.While previous BSM searches explored the t tt t final state for several target models [21,22], this analysis uniquely attempts to reconstruct the resonance explicitly, enabling a search of a new physics contribution with minimal model dependence.
In addition, a simplified model [23] considering a coloursinglet vector particle Z is used to generate simulated samples for further model-dependent interpretations, however without relying on how the resonance couples to other particles in a specific model.The interaction Lagrangian reads: L = tγ μ (c L P L + c R P R ) t Z μ = c t tγ μ (cos θ P L + sin θ P R ) t Z μ , where γ μ are the Dirac matrices and P L/R = (1 ∓ γ 5 )/2 are the chirality projection operators with γ 5 = iγ 0 γ 1 γ 2 γ 3 .The coupling of the vector singlet to the top quarks is defined as c t = c 2 L + c 2 R , with components that couple only to lefthanded and right-handed top-quarks c L and c R , respectively.The tangent of the chirality angle is defined as tan θ = c R /c L .For the mass range which is explored in this search with a Z mass (m Z ) much larger compared to the top-quark mass (m t ), the Z decay width can be approximated by /m Z ≈ c 2 t /(8π).To minimise model dependence, only the tree-level production of a heavy top-philic resonance is considered.At the LHC a large contribution comes from t t Z (Fig. 1a) which is independent of θ .Further contributions arise from t j Z (Fig. 1b) and t W Z (Fig. 1c) production modes, which are negligible for θ = π/2, and are largest for θ = 0 with similar magnitude as t t Z .
This search uses the data recorded by the ATLAS detector in ppcollisions at √ s = 13 TeV between 2015 and 2018, corresponding to 139 fb −1 .In the SM, the top quark is expected to decay into a W boson and a b-quark with a branching ratio of approximately 100% with subsequent hadronic or leptonic decay of the W boson.This analysis targets final states for which both resonance top quarks decay hadronically and one of the spectator top quarks decays leptonically.The singlelepton channel is preferred over the fully hadronic channel due to the significantly lower background arising from multijet processes, and the ability to trigger on events with at least one lepton.The top quarks from the resonance decay are expected to be highly boosted, and therefore their hadronic decays are reconstructed using jets with a large radius parameter and identified by requiring the jets to have a large mass and momentum.This final state with one lepton 1 suffers from a large background that is mostly composed of t t production in association with additional jets (t t+jets), especially when the associated jets contain b-hadrons.A data-driven technique assisted by simulations is used to overcome the challenge of modelling the t t+jets background.The presence of a signal is tested in signal-enriched regions using an approach with minimal model dependence followed by a model-dependent interpretation.

The ATLAS detector
The ATLAS detector [24] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry with a nearly 4π coverage in solid angle. 2 It consists of an inner tracking detector (ID) surrounded by a 1 In the rest of this article, 'lepton' refers exclusively to an electron or a muon.Tau leptons are not explicitly considered although their decay products can be accepted by the electron, muon and jet selection criteria. 2ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the zaxis 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 superconducting solenoid, electromagnetic (EM) and hadron calorimeters, and a muon spectrometer (MS).
The ID covers the pseudorapidity range |η| < 2.5.The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track, the first hit normally being in the insertable B-layer (IBL) installed before Run 2 [25,26].It is surrounded by a silicon microstrip detector and a straw tube transitionradiation tracking detector.The calorimeter system covers the pseudorapidity range |η| < 4.9.It consists of lead/liquidargon (LAr) sampling calorimeters which provide EM energy measurements with high granularity.A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|η| < 1.7).The endcap and forward regions are instrumented with LAr calorimeters for EM and hadronic energy measurements up to |η| = 4.9.The MS surrounds the calorimeters and is based on three large air-core toroidal superconducting magnets with eight coils each.The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector.The MS includes a system of precision tracking chambers and fast detectors for triggering.The ATLAS trigger and data acquisition system [27] consists of a two-level trigger system in order to select events.The firstlevel trigger is implemented in hardware and uses a subset of the detector information to accept events at a maximum rate of nearly 100 kHz.The second-level is a software-based trigger that reduces the accepted event rate to 1 kHz, on average, depending on the data-taking conditions [27].An extensive software suite [28] is used for real and simulated data reconstruction and analysis, for operation and in the trigger and data acquisition systems of the experiment.

Object reconstruction and event selection
This analysis uses a set of data events collected by the ATLAS detector between 2015 and 2018 at √ s = 13 TeV.Only events for which all detector subsystems were operational are considered.The data set corresponds to an integrated luminosity of 139 fb −1 [29,30].
This analysis is based on events where the detector readout is triggered by the presence of at least one electron or one muon, referred to as single-lepton triggers.Events are selected if the leptons satisfy either low transverse momentum p T thresholds, with identification and isolation requirements, or a looser identification criterion with no isolation requirement and higher p T thresholds.The lowest p T requirement used in the single-lepton triggers varies from 20 to 26 GeV depending on the data-taking period and the lepton flavour [31,32].
Events are required to have at least one primary vertex reconstructed from at least two ID tracks with p T > 500 MeV.For events with several primary vertices, the one with the largest sum of the squared transverse momenta of the associated tracks is taken [33].
Electron candidates are reconstructed from an isolated electromagnetic calorimeter energy deposit which is matched to a track in the ID [34].The pseudorapidity of the calorimeter energy cluster, η cluster , must satisfy |η cluster | < 2.47, excluding the transition region between the barrel and the end-cap EM calorimeter (|η cluster | / ∈ [1.37, 1.52]).Muon candidates are reconstructed by combining tracks in the ID with tracks in the MS [35] and are required to have |η| < 2.5.Both electron and muon candidates are required to have a p T above 28 GeV.Further selections on the longitudinal and transverse impact parameters are imposed.The transverse impact parameter divided by its estimated uncertainty, |d 0 |/σ (d 0 ), has to be less than five (three) for electron (muon) candidates.The longitudinal impact parameter z 0 is required to be |z 0 sin(θ )| < 0.5 mm for both electron and muon candidates.Electron candidates must pass a 'Tight' likelihood-based identification working point [34] employing calorimeter and tracking information that provides separation between electrons and jets.They are also required to be isolated according to the 'Tight' selection, a criteria based on the properties of the topological clusters in the calorimeter and of the ID tracks around the reconstructed electron [34].Muon candidates must satisfy the 'Medium' cut-based identification working point [36].They are also required to be isolated using the 'TightTrackOnly' selection, a criteria based on the properties of the ID tracks around the reconstructed muon [36].
Jets are reconstructed from topological clusters of calorimeter cells and tracking information using the anti-k t algorithm [37,38] with a radius parameter of R = 0.4, processed using a particle-flow algorithm [39].They are referred to as 'small-R jets'.Jets are required to have p T > 25 GeV and |η| < 2.5 and are calibrated as described in Ref. [40].The Jet Vertex Tagger (JVT) discriminant [41] is used to identify jets originating from the hard-scatter interaction through the use of tracking and vertexing information.Jets with p T < 60 GeV and |η| < 2.4 are required to satisfy the requirement JVT > 0.50, corresponding to a selection efficiency for hard-scatter jets of about 96% [42].The DL1r classification algorithm based on recurrent neural networks [43,44] is used to identify jets containing b-hadrons (bjets).This algorithm identifies b-jets against backgrounds of light-flavour and charm-quark-initiated jets using information about the impact parameters of tracks associated with the jet and the topological properties of the displaced vertices reconstructed within the jet.This analysis uses a b-jet efficiency working point of 77% as measured for jets with p T > 20 GeV and |η| < 2.5 in simulated t t events [45].The tagging algorithm gives an expected rejection factor (defined as the inverse of the efficiency) of about 192 against light-flavour jets, and about 5.6 against jets originating from charm quarks [44,46].
To resolve the potential ambiguities of a single detector response being assigned to two objects by the reconstruction algorithm, a sequential overlap removal procedure is applied.First, any electron found to share a track with another electron with higher p T is removed.Second, electrons sharing their track with a muon candidate are removed.Third, the closest jet in rapidity y and in φ using R y = ( y) 2 + ( φ) 2 = 0.2 of an electron is removed. 3Fourth, the electrons within R y = 0.4 of a remaining jet are removed.Fifth, the jets with less than three associated tracks that are within R y = 0.2 of a muon are removed.Finally, remaining muons are removed if their track is within R y = min(0.4,0.04+10 GeV/ p Tμ ) of a remaining jet.
After the overlap removal, the selected and calibrated small-R jets are used as inputs for jet reclustering [47] using the anti-k t algorithm with a radius parameter of R = 1.These reclustered jets are referred to as 'large-R jets', and are used as proxies of the hadronically decaying top quarks in this analysis.The calibration corrections and uncertainties for the reclustered large-R jets are inherited from the small-R jets [47].A trimming procedure [48] is applied to reclustered large-R jets which removes all the associated small-R jets that have a p T below 5% of the p T of the reclustered jet to suppress gluon radiation and mitigate pile-up effects.The large-R jets are required to have p T > 300 Gev, |η| < 2, mass m > 100 Gev, and at least two constituent jets.
At preselection level, the events are required to have exactly one lepton with p T > 28 Gev that matches the lepton that fired the trigger, and at least two large-R jets to capture the two resonance top quarks that decay hadronically.Additional jets are expected from the decay of the spectator top quarks for t t Z signal events, as well as from the associated production in the t j Z and t W Z production modes.The number of additional jets is counted using small-R jets with R > 1 to any of the selected large-R jets in an event.Events are further required to have at least two additional small-R jets and at least two b-tagged small-R jets.The btagged jets can be either inside or outside of the large-R jets and therefore can also count towards additional jets.

Event simulation
Monte Carlo (MC) samples of simulated events were produced to model the signal and background processes.The samples are normalised using the best theory predictions available.The generated event samples are processed through the full ATLAS detector [49] based on Geant4 [50].Only the MC samples describing SM t tt t production and systematic variations of the t t simulation are processed through a faster simulation making use of parameterised showers in the calorimeters [51].Additional simulated ppcollisions generated using Pythia8.186[52] using the A3 set of tuned parameters (tune) [53] and the MSTW2008LO parton distribution function (PDF) set [54] were overlaid to model the effects of multiple interactions in the same and nearby bunch crossings (pile-up).The distribution of the number of additional ppinteractions in the MC samples is re-weighted to match the one observed in data.All simulated samples were processed through the same reconstruction algorithms and analysis chain as the data.For all samples of simulated events, except those generated using Sherpa [55], the EvtGen 1.2.0 program [56] was used to describe the decays of bottom and charm hadrons.The NNPDF3.0NLO [57] PDF set was used in all matrix element (ME) calculations if not stated otherwise.
The signal samples are simulated using the MadGraph5_aMC@NLO 2.8.1 [58] generator at leadingorder (LO) in the five-flavour scheme with the NNPDF3.1LO[57] PDF set.The model implementation is based on the simplified top-philic resonance of Ref. [59] where a colour singlet vector particle is considered to exclusively couple to top and anti-top quarks.Six mass points are generated, spaced according to the experimental resolution.The Z masses are: 1.0 TeV, 1.25 TeV, 1.5 TeV, 2.0 TeV, 2.5 TeV, and 3.0 TeV.Interference effects between s-and t-channel processes for t t Z production are simulated.The interference with SM four-top-quark production was found to be negligible.Only the resonant Z production mode is included in the simulation of the t W Z and t j Z channels.The signal samples are generated with baseline choices of the chirality parameter θ = π/4 and the coupling between Z and top quarks of c t = 3.Other choices of these parameters are realised using MadGraph matrix element reweighting [60].The width of the resonance is computed automatically using MadGraph5_aMC@NLO and corresponds to a relative width of about 4% for c t = 1 and increases to 64% for c t = 4.
The production of t t events is modelled using the HVQ program [61,62] in the Powheg Box 2 [61,[63][64][65] generator at next-to-leading order (NLO) in QCD.The h damp parameter, which controls the transverse momentum p T of the first additional emission beyond the Born configuration, was set to 1.5 m t [66].
The t t+jets MC events are classified according to the flavour of the particle jets.The particle jets are reconstructed from the simulated stable particles using the anti-k t algorithm with a radius parameter R = 0.4, and are required to have p T > 15 GeV and |η| < 2.5.Events are labeled as t t+≥1b if at least one particle jet is matched within R < 0.3 to any b-hadron with p T > 5 GeV.In the remaining events, if at least one particle jet is matched within R < 0.3 to any chadron with p T > 5 GeV, the events are labeled as t t+≥1c.Only hadrons not associated with b-and c-quarks from topquark and W -boson decays are considered.All other events are labeled as t t+light.Events categorised as t t+≥1b and t t+≥1c are collectively referred to as t t +heavy flavour (HF) events.
Samples of single-top quark production backgrounds corresponding to the W t associated production, s-channel and tchannel production, were modelled using the Powheg Box 2 generator at NLO in QCD.Overlaps between the t t and W t final states are removed using the 'diagram removal' scheme [67].
A sample of t t + H events is modeled using the Powheg Box generator at NLO.The production of t t + V (with V = W, Z , including non-resonant Z /γ * contributions and t t + W W ) events are modeled using the Mad-Graph5_aMC@NLO 2.3.3 generator at NLO.Samples of V +jets (V = W, Z ) events are generated with the Sherpa 2.2.1 [68] generator using NLO-accurate MEs for up to two partons and LO-accurate MEs for up to four partons.Samples of diboson final states (V V ) were also simulated using the Sherpa 2.2.1 generator.The production of SM t tt t events is modeled using the MadGraph5_aMC@NLO 2.6.2 generator at NLO with the NNPDF3.1NLO[57] PDF.The rare process of single-top-quark associated Z boson production, t Z, is modeled using the MadGraph5_aMC@NLO generator at LO.
All events generated using Powheg Box or MadGraph5_aMC@NLO are interfaced with Pythia8.230[69] using the A14 tune [70] and the NNPDF2.3LO[57] PDF set.Events generated using Sherpa employ the set of tuned parameters developed by the Sherpa authors.
To assess the uncertainty due to the choice of generator, the t t and single-top-quark samples produced with the nominal generator set-ups are compared with alternative samples generated with MadGraph5_aMC@NLO (interfaced with Pythia8) for the calculation of the hard-scattering processes.These alternative samples were generated using the same PDF in the ME as in the nominal samples.Additional t t and single-top-quark samples were produced by replacing Pythia8 with Herwig 7.04 [71,72] for parton showering and hadronisation, using the H7UE tune [72] and the MMHT2014LO [73] PDF set.These samples are used to evaluate uncertainties due to the choice of parton shower and hadronisation model.

Analysis strategy
This analysis targets events where both resonance top quarks decay hadronically.For a heavy resonance signal with m Z ≥ 1.0 TeV, the two resonance top quarks are expected to be highly boosted.Large-R jets are used as proxies for the two highly boosted top quarks.The invariant mass distribution of the two large-R jets with the highest p T , labelled as m JJ , is scanned for an excess over the background prediction.
The m JJ distributions for the top-philic Z signal with different m Z and c t values are shown in Fig. 2. A peak near the generated value of m Z is typically observed for masses m Z ≤ 2 TeV and coupling strengths c t ≤ 1.The peak at the resonance mass in general is wider for higher Z masses and for higher c t values due to the increase of the Z natural width.All distributions are skewed towards the lower mass end for multiple reasons.One factor is the contamination from large-R jets that do not capture or only capture part of the decay products from a hadronically-decaying resonance top quark.Moreover, the parton luminosity effect and QCD radiation from the highly boosted top quarks also contribute to the lower mass side.The difference between the two chirality hypotheses, θ = 0 and θ = π/2, is small, with a slightly better mass resolution in case of θ = 0.This is due to the maximal contribution from t j Z and t W Z production modes, which have smaller ambiguity in reconstructing the Z resonance using the two large-R jets given the lower expected jet and b-jet multiplicities.
Despite the above effects that degrade the resolution for reconstructing the Z mass, the m JJ distributions for the signal with small m Z and c t are still distinct from those for background, as shown in Fig. 2. For the signal hypotheses with larger m Z and c t , the distributions become more background-like given the increasing resonance width and the low-mass shoulder.According to MC simulations, about 90% of background after the preselection consists of t t+jets events.Other background events mainly come from the production of t t Z, t t W , t t H, single top-quark and W/Z bosons in association with jets.The SM t tt t, diboson and other rare processes contribute less than 1% of all background events.For background events, a smoothly falling m JJ distribution is expected once m JJ passes the turn-on threshold due to the minimum requirements on the large-R jet mass and p T .
The search is performed in the range of the m JJ distribution between 1.0 TeV and 3.2 TeV.The lower bound was chosen to minimise the impact due to the turn-on threshold on the background modelling, whereas the upper bound was motivated by the vanishing amount of expected SM background events at the high end of the m JJ spectrum.The product of the efficiency and acceptance of the signal in the t t Z production mode after the preselection varies from 3.1% to 7.2% for the different signal hypotheses, estimated using the simulated samples described in Sect. 4.
In the single-lepton final state, Z signal events are expected to have a high multiplicity of jets in addition to the two large-R jets (N add.The predictions are given by MC simulations.The contributions from the SM t tt t, diboson, and V H processes are combined in the category 'Other'.Four signal samples are presented with m Z = 1.5 and 3.0 TeV and θ = 0 and π/2.All distributions are normalised to unit area discriminating power to separate these contributions.These two variables are used to categorise events into regions with different signal-to-background ratios.As illustrated in Fig. 4, nine regions are defined, denoted by 'Na', 'Mb', where Na and Mb represent N add.-jets and N b-jets , respectively, and both range from 2 to ≥4.The region with the lowest signal contamination, (2a, 2b), is referred to as the 'source region' and used to derive the data-driven background estimate detailed in Sect.6.The highest N add.-jets and N b-jets regions are chosen according to the expected signature of the t t Z events and help to maintain a statistically significant background estimation.All regions with at least three b-tagged jets have relatively higher expected signal contributions compared to other regions.They are used for signal extraction and are referred to as the 'signal regions'.In the (3a, 2b) and (≥4a, 2b) regions, the small signal contribution provides little gain in sensitivity, and the expected signal contamination is not negligible for the data-driven background estimate.Therefore they are not used as either signal nor source region.However, they are used as validation regions to verify the modelling of the background.Scaling the top-philic simplified signal model to the expected exclusions achieved by this search, the largest signal contamination is below 2% in the source region and below 5% in the validation regions, depending on the model parameter choices.Figure 3 illustrates the background composition in bins of N add.-jets and N b-jets expected from the MC predictions.The fractions of non-t t background are similar in most bins.The most dramatic change is in the fraction Fig. 4 An illustration of the analysis regions defined using N add.-jets and N b-jets .The source, validation and signal regions are highlighted in blue, gray and red, respectively.The names of the regions used in the rest of the document are also shown of t t+≥1b events.From MC simulations, the fraction of t t+≥1b events with respect to the total sum of background events is expected to increase from 8% in regions with two b-tagged jets to 68% in the most signal-like regions with at least four b-tagged jets.

Background estimation
The search for the resonance in the m JJ distribution heavily relies on the modelling of the continuum background.However, the MC prediction is expected to be unreliable and susceptible to large modelling systematic uncertainties due to various factors discussed in the following.The top-philic Z signal in the one-lepton channel features a large jet multiplicity.For the main background, t t+jets, events with a similar final state contain multiple partons beyond Born level.This means that in the nominal t t+jets sample generated using NLO ME, most of the jets in the final state are from the parton shower with limited precision [74].These jets form the large-R jets used to construct m JJ .Additional mismodelling and uncertainties come from the prediction of the boosted large-R jets, which are subject to the often badly modelled collinear radiation.Finally, the t t +HF background that populates the signal regions is underestimated by the current MC predictions [75,76].
A data-driven technique is used to provide a reliable background prediction.It is based on the similarity of the background m JJ shape in all analysis regions observed in MC simulations.Therefore, the shape of the background is obtained from data in the source region where there is negligible signal contamination.The background shape is then extrapolated Fig. 5 The distribution of data in the source region, along with the fitted function to data.The bands show the three independent variations obtained by diagonalising the covariance matrix of the fit.The third eigen variation is not visible due to its smallness to the signal regions.The extrapolation factors are derived from the MC simulations, accounting for the reduced event rate in higher N add.-jets and N b-jets bins and the minor differences in the shape of m JJ between the regions.The final background model is obtained using a profile likelihood fit via further adjustments according to all systematic variations considered in this analysis.
Firstly, the shape of the background is obtained by fitting the data in the source region, (2a, 2b), using the following parameterization [77,78]: where x = m JJ / √ s, with √ s being the pp collision centerof-mass energy, and p 1 to p 3 are the fit parameters controlling the shape of the function.The fit result is presented in Fig. 5, showing that the fitted function describes the data well within the uncertainties due to the limited number of data events.Alternative parameterizations considering additional parameters or alternative functional forms have been tested, but were found to provide no improvement in the performance of the background estimate.
The fitted function from the previous step is scaled by extrapolation factors for each signal region, in order to obtain the background prediction in that region.The extrapolation factors are derived from MC simulation in bins of m JJ , as the ratio of the expected number of events in the region of interest to that in the source region.To minimize the impact due to the limited number of MC events, the function in Eq. 1 is used, fitting the m JJ distribution predicted by MC simulation in both the region of interest and the source region.The extrapolation factors are then obtained by taking the ratio of these fitted functions.The shape of the background in the various analysis regions differs by up to ±15% with respect to that in the source region, according to the prediction of the extrapolation factors.The extrapolation is performed with all backgrounds combined.The shape of the m JJ distribution from the different background processes are compared using MC simulations.For the t t+jets, single top-quark and t t Z/W/H processes, constituting about 95% of the total background, the shapes are consistent within the uncertainties due to the limited number of MC events.Furthermore, the composition and the shape of the small backgrounds, i.e. those other than the t t+jets backgrounds, are found to be similar across all analysis regions.As an additional check, alternative extrapolation factors were derived with t t+jets MC simulation only, to account for the effect due to variation in the non-t t backgrounds.The resulting difference was found to be smaller than the uncertainty due to the limited number of MC events, and is therefore neglected.
The extrapolation factors are based on the MC predictions in each of the analysis regions.Systematic uncertainties affecting the predicted distributions (m JJ , N add.-jets and N b-jets ) are therefore propagated to the background estimate via the extrapolation factors.These include modelling uncertainties in all background MC samples and experimental uncertainties.For each systematic variation, the extrapolation factors are rederived using the same procedure as for the nominal background prediction.To account for the deficit of t t + HF events in MC simulation, additional normalisation uncertainties are assigned to t t+≥1c and t t+≥1b events.Furthermore, dedicated uncertainties in the data-driven background estimate and the extrapolation are considered.The details on the systematic uncertainties included in this analysis are described in Sect.7. All systematic uncertainties are incorporated in a profile likelihood fit to data in the signal regions, as discussed in Sect.8.

Systematic uncertainties
Different sources of systematic uncertainty affect the search presented here, including those related to the luminosity, the identification and reconstruction of the physics objects, the MC simulation of the signal and background processes, and the method used to estimate the background.In the following, a brief description of the sources of systematic uncertainty is provided.A particular emphasis is put on the ones related to the t t background prediction, which will be shown to have the largest impact on the sensitivity of the search.The systematic variations can affect the normalisation of the signal and background templates estimated in the different regions as well as the shape of the m JJ distributions.All systematic uncertainties on the background prediction enter the analysis exclusively via extrapolation factors, except for the uncertainties related to the functional fits to data and MC described in Sect.7.2, and the signal bias uncertainty described in Sect.7.5.The luminosity uncertainty only applies to the signal.

Experimental uncertainties
The uncertainty in the combined 2015-2018 integrated luminosity is 1.7 % [29], obtained using the LUCID-2 detector [30] for the primary luminosity measurements.This uncertainty affects the signal prediction in the model-dependent interpretation.
Other experimental uncertainties arise from corrections and calibrations applied to MC simulations.These uncertainties affect the background estimate via the extrapolation factors, as well as the MC prediction of the Z signals.An uncertainty is considered for the reweighting factors that correct the pile-up profile in MC simulations to match those in data.Uncertainties on the modelling of leptons arise from their momentum and energy scale calibration and resolution, as well as the trigger, reconstruction, identification, isolation efficiencies [34,35].Uncertainties on the modelling of jets mainly come from their energy scale (JES) and resolution (JER), containing effects from jet flavour composition, single-particle response, and pile-up [79,80].An uncertainty is assigned for the efficiencies of the JVT requirement on jets [42].Uncertainties are considered for the calibration of the b-tagging efficiencies, including the efficiencies of tagging bjets as well as the rates of mis-tagging c-jets and light-flavour jets [45,46,81].

Uncertainties on functional fit and extrapolation
Dedicated uncertainties are assigned for the data-driven background estimate and the MC extrapolation from the source region to other regions.The uncertainties from each functional fit are obtained from the covariance error matrix of the fit.Three statistically independent variations are extracted from an eigen-decomposition of the covariance matrix.This leads to three uncertainties associated to the functional fit in the source regions to data and MC, respectively.In addition, three uncertainties are also associated to the fit of the MC in each of the six signal regions which leads to 24 components of uncertainty in total.The resulting uncertainties due to the functional fit to data and MC in the source region are both correlated across the signal regions.

Theoretical uncertainties for the t t background
An uncertainty of 50% in the normalisation of the t t+≥1b events as well as the t t+≥1c events is applied [75,76], and an uncertainty of 10% is considered for the t t+light events.The uncertainties due to the choice of generator and parton shower model used to simulate the inclusive t t sample are evaluated by comparing the nominal t t sample with alternative t t samples, detailed in Sect. 4. The uncertainties associated with the choice of NLO generator are estimated by comparing the predictions of Powheg Box and MadGraph5_aMC@NLO.The effect of the choice of parton shower and hadronisation model is estimated from comparing the prediction of Powheg + Pythia 8 to Powheg + Herwig 7.04.These two uncertainties are split into four components, each affecting either only the shape of the m JJ distributions or the acceptance and migration for the 3b regions and the ≥4b regions separately.This treatment is motivated by the different compositions of t t+light, t t+≥1c and t t+≥1b events in the two b-jet multiplicity bins and the large effect on the overall acceptance and migration of the events across the regions due to different parton shower and hadronisation models.Uncertainties due to missing higherorder QCD corrections are estimated by separately varying the renormalisation and the factorisation scales by factors of 2.0 and 0.5 in the nominal t t sample and taking the envelope.Additionally, uncertainties in the amounts of initialand final-state radiation (ISR and FSR) from the parton shower (PS) are assessed by respectively varying the corresponding parameter of the A14 PS tune4 and by varying the FSR renormalisation scale by factors of 2.0 and 0.625.The uncertainty related to the parton distribution function (PDF) is evaluated by using the PDF4LHC systematic variations [82].

Modelling uncertainties for non-t t backgrounds
Non-t t background processes represent a minor fraction of the total background.Therefore, the relevant uncertainties have a small impact on the result of the analysis.Similar to the t t background, uncertainties due to missing higher-order QCD corrections, the amounts of initial-and final-state radiation, and the parton distribution function are also considered for most of the non-t t backgrounds.
An uncertainty of 30% in the total cross section of the three single-top-quark production modes is included.This conservative number is motivated by the uncertainties due to the modelling of several associated jets and heavy flavour jets in the production.Uncertainties associated with the parton shower and hadronisation model and the generator choice are evaluated by comparing the nominal Powheg + Pythia 8 sample for each process with alternative samples produced with Powheg + Herwig 7 and aMC@NLO + Pythia 8.An additional uncertainty on the interference between t W and t t production at NLO is evaluated by comparing the nominal t W sample produced using the diagram removal scheme with an alternative sample produced with the same generator but using the 'diagram subtraction' scheme [67].
Modelling uncertainties in the t t W , t t Z, and t t H processes are evaluated in a similar way.Uncertainties of 60%, 15%, and 20% are applied to the t t W , t t Z, and t t H cross sections, respectively [83][84][85].
An uncertainty of 60% is assumed for the V +jets production cross section.It is estimated by adding a 24% uncertainty in quadrature for each additional jet based on a comparison among different algorithms for merging LO matrix elements and parton showers [86].An uncertainty of 20% is assumed for the SM t tt t production cross section [87].A conservative cross section uncertainty of 50% is applied to all other small background processes.

Signal bias uncertainty
A dedicated uncertainty is considered in the model-dependent interpretation of the results, accounting for the bias on the extraction of the signal yields caused by the background model.The details on the signal extraction and the modeldependent interpretation are described in Sect.8. To evaluate the size of the bias, a large number of pseudo-datasets are sampled using the expected background directly from MC simulations.The model-dependent signal extraction (as described in Sect.8.3) is performed for each of these pseudodatasets, taking into account all systematic uncertainties described previously.The resulting distribution of the signal strengths is fitted with a Gaussian function.The deviation of the fitted central value from zero is taken as the size of the bias.This procedure is repeated for all signal mass points.The three points with the largest bias are used to compute a second order polynomial function.The eventual size of the bias for each mass point is evaluated from the fitted polynomial function.To obtain dedicated uncertainties for all choices of model parameters, this procedure is repeated for all distinct c t and chirality parameter θ values explored in this search.The uncertainty is determined using the MC signal template scaled by the corresponding value of bias, and is implemented as an uncertainty on the background expectation.Depending on the analysis region and signal point, this translates to uncertainties up to 14% in the individual signal regions, with the majority of values being a few percent in size.

Statistical analysis
The search for a resonance signal is conducted using the binned m JJ distributions in all signal regions.An approach with minimal model dependence is adopted, followed by a model-dependent interpretation.A profile likelihood fit is used to obtain the final background model and make further statistical inference regarding the presence of a signal.The search with minimal model dependence is performed using BumpHunter [88], a hypothesis-testing tool that searches for local data excesses compared to the expected background.To obtain the expected background, the profile likelihood fit is performed with a background-only hypothesis.In the model-dependent interpretation, the Z signal samples described in Sect. 4 are used to interpret the data.Profile likelihood fits are performed with the signal-plus-background hypothesis, testing the compatibility between data and the models with different m Z , c t and chirality parameter θ .

Profile likelihood fit
The statistical analysis is based on a binned profile likelihood function.The bin width is chosen to be 100 GeV, with two large bins of 500 GeV and 700 GeV at the high mass end of the m JJ spectrum to avoid empty bins.Each bin in the signal regions is represented by a Poisson probability term for the observed data, with the total expected yield given by the background model described in Sect.6 and the Z signal samples.The likelihood function L(μ, θ θ θ) is constructed as the product of the Poisson probability terms over all bins.This function depends on the signal-strength parameter μ, a multiplicative factor to the assumed pre-fit signal cross section, and θ θ θ , a set of nuisance parameters encoding the effect of systematic uncertainties.The nuisance parameters are implemented in the likelihood function as Gaussian or log-normal constraints.The total number of expected events in a given bin therefore depends on μ and θ θ θ.The fit is performed by finding the values of μ and θ θ θ that maximise the likelihood function.In the case of a fit with the background-only hypothesis μ is fixed to 0. The nuisance parameters θ θ θ allow variations of the expectations for signal and background according to the corresponding systematic uncertainties, whilst penalising the likelihood function for any deviation according to the specified constraints.The fit also reduces the impact of systematic uncertainties on the search sensitivity by exploiting the highly populated and background-dominated regions.The result of the profile likelihood fit is further used for the signal extraction.

Search with minimal model dependence
In the search using BumpHunter, the data and the expected background in each of the six fitted regions are compared in sliding windows of variable sizes.The smallest window is required to contain two bins with the binning shown in Figs.7 and 8 in Sect.9, corresponding to a width that is slightly smaller than the expected resolution extracted from the Z MC samples.The Poisson probability is evaluated for all windows.For each region, the window with the smallest Poisson probability is chosen to be the most prominent window.The corresponding probability p min is used to construct the BumpHunter test statistic t defined as where d and b represent the number of observed data and the expected background events in the window, respectively.A deviation of data from background is only considered as evidence against the background-only hypothesis if data exceed the expectation.To obtain the local p-value and the significance of the most interesting bump, 10 5 pseudo-experiments are sampled from the expected background, and the t value from data is compared to the distribution of the t values from the pseudo-experiments.The global p-value and significance are then computed for the most prominent window of each region, taking into account the trial factors that incorporate the look-elsewhere effect [89].

Model-dependent interpretation
The model-dependent interpretation is based on the signal samples described in Sect. 4. The signal strength is extracted by performing the profile likelihood fit with the signal-plus-background hypothesis.The fit and the subsequent statistical analysis are performed for each mass point.The test statistics are defined based on the likelihood ratio λ(μ) = L(μ, θ θ θ μ )/L( μ, θ θ θ), where μ and θ θ θ are the values of the parameters that maximise the likelihood function, and θ θ θ μ are the values of the nuisance parameters that maximise the likelihood function for a given value of μ.To evaluate the compatibility of the observed data with the background-only hypothesis, the test statistic q 0 , defined as is used.The resulting p-value represents the compatibility p 0 .A small p 0 indicates that the background-only model does not describe the data well, and the significance of the corresponding signal is further examined.In the absence of any significant excess above the background expectation, upper limits on the signal production cross section at 95% confidence level are derived.The test statistic q μ , defined as is used with the CL S method [90,91].Both p 0 and the upper limits are computed using the asymptotic approximation [92].All statistical analyses are performed using the RooStats framework [93][94][95].

Results
The results of the statistical analysis are presented in this section.Prior to the analysis using real data, the profile likelihood fit model and the statistical methods were tested against pseudo-datasets constructed using simulated events.Both the background modelling and the signal extraction were validated, using background-only and signal-plus-background pseudo-datasets.In particular, alternative t t background predictions were considered when building the pseudo-datasets.They were used to stress-test the fit, especially the incorporated systematic uncertainties related to t t background modelling.These alternative predictions include an enhanced composition of the t t + HF backgrounds, as well as a t t sample generated with Sherpa2.2.10,5 which is not used in the definition of any uncertainty.These stress-tests demonstrated the capability of the fit to constrain the relevant systematic uncertainties, and to model the data correctly in case of mismodelling in the extrapolation factors from MC predictions.The search result with minimal model-dependence is derived based on a simultaneous fit to all signal regions with a background-only hypothesis.The background estimate is examined in validation regions.This is done by propagating to the validation regions the post-fit nuisance parameters obtained from the fit to the signal regions.The resulting agreement between data and the estimated background in the validation regions is shown in Fig. 6.The largest discrepancies are evaluated and found to be insignificant after taking into account the look elsewhere effect.Figure 7 shows the post-fit m JJ distribution with the background-only hypothesis compared to the observed data in the two most sensitive signal regions, along with the results from BumpHunter.The rest of the signal regions are shown in Fig. 8.All fitted distributions show reasonable agreement with data within uncertainties.The goodness of fit was evaluated using a likelihood-ratio test in which the likelihood of the nominal fit is compared to that of a saturated model [102].The result indicates a good description of the data by the background model.The post-fit nuisance parameter values are all within 1σ of the prior uncertainty.No significant excess was found in the BumpHunter search.The most significant deviation between data and the background expectation was observed in the (2a, ≥4b) region at 1.2 TeV.The window with the largest deviation has a local significance of 1.04σ ; the corresponding global p-value is 0.15.The expected distribution with the presence of a top-philic Z signal from MC simulation is illustrated in the top panels of Figs. 7 and 8.The top-philic Z signal considered has m Z = 1.5 TeV with c t = 1 and θ = π/2, and is normalised to have an arbitrary large cross section of 51 fb, obtained using the cross section predicted by the simplified model [59] with these parameters scaled by a factor of 100.
The model-dependent signal extraction also reveals no significant excess in data.For all model parameter choices, the fitted signal strength μ is compatible with zero.The smallest local p 0 value is 3.8%, obtained with the 1.25 TeV mass point for c t = 1 and θ = 0.In addition, the fitted values of the nuisance parameter are compared to those obtained in the fit with the background-only hypothesis and found to be consistent.The observed and expected upper limits on the signal production cross section at 95% confidence level are shown for specific signal scenarios in Fig. 9.The observed (expected) limits range from 21 (14) fb to 119 (86) fb depending on the choice of model parameters.The strongest observed exclusion is observed at m Z = 2.5 TeV with c t = 1 and θ = π/2.At the higher mass range cross section limits are stronger for the smaller couplings c t due to the narrower signal width.While for the cross section limits there is no strong dependence on the chirality parameter θ , one can observe a stronger model constraint for θ = 0 because of the additional contributions from t j Z and t W Z production to the expected signal cross section.For c t = 4 top-philic Z mass of 1.0 TeV is excluded for both chirality parameter choices (θ = 0 and π/2).
Table 1 lists the impacts of various groups of uncertainties relative to the total uncertainty on the fitted signal strength for two examples of the top-philic Z signal: (m Z = 1.5 TeV, c t = 1, θ = 0) and (m Z = 3.0 TeV, c t = 4, θ = 0).The results are dominated by systematic uncertainties.The most important source of systematic uncertainty is the modelling of the t t+jets background.The uncertainties due to the migration and acceptance effect of the ME and PS generator choices in the ≥4b regions play a dominant role.They contribute 26% and 31% (36% and 48%) of the total uncertainty, respectively, for the signal with m Z = 1.5 TeV and c t = 1 (m Z = 3.0 TeV and c t = 4).This reflects the large difference between the predictions of the different generators in the phase space explored by this search.The uncertainties in the normalisation of t t+≥1b also contribute 18% of the total uncertainty for both signal hypotheses shown in Table 1.The uncertainties from the jet energy scale and resolution are also an important source, given the large jet multiplicity in the signal regions.This is followed by the uncertainties from the data-driven background estimate.In each plot two coupling-strength scenarios are shown, for c t = 1 (in gray) and c t = 4 (in black).The light (dark) blue curves illustrate the signal cross sections based on the model described in Ref. [59] for c t = 1 (c t = 4).The green bands surrounding the expected limits for c t = 1 and c t = 4 correspond to the 68% confidence intervals.All limit lines are obtained by interpolating linearly between the different mass hypotheses Table 1 The contribution from different systematic uncertainties relative to the total uncertainty on the fitted signal strength for two signal scenarios with θ = 0, grouped into categories.For each category, the fit is repeated with the corresponding group of nuisance parameters fixed to their best-fit values.The contribution from each category is then evaluated by subtracting in quadrature the uncertainty on the signal strength obtained in this fit from that of the full fit with all uncertainties.The percentage is calculated relative to the total uncertainty from the full fit.The contribution from the statistical uncertainty is also shown.The total systematic uncertainty is different from the sum in quadrature of the different groups due to the correlations among the nuisance parameters in the fit

Uncertainty categories
Relative contribution to the total uncertainty [%]

Conclusion
A search for a top-philic heavy resonance produced in association with a top quark or a top-quark pair and decaying to t t is presented.The search is performed using 139 fb −1 of pp collision data at √ s = 13 TeV collected by the ATLAS detector at the LHC.Large-R jets are used as proxies of the top quarks from the heavy resonance decay.The invariant mass spectra of the two large-R jets with the largest p T are used to test for the presence of a resonance signal in the range of 1.0 TeV to 3.2 TeV.Events in the single-lepton final state are selected and categorised according to the number of additional jets and b-tagged jets.The search is conducted in regions with at least three b-tagged jets, with and without assuming a specific Z signal model.No significant excess was observed above the expected background.The upper limits on the Z production cross section at 95% confidence level are computed for signals with six values of m Z between 1.0 TeV and 3.0 TeV based on a simplified model.The observed (expected) limits range from 21 (14) fb to 119 (86)  Trust, United Kingdom.The crucial computing support from all WLCG partners is acknowledged gratefully, in particular from CERN, the ATLAS Tier-1 facilities at TRIUMF (Canada), NDGF (Denmark, Norway, Sweden), CC-IN2P3 (France), KIT/GridKA (Germany), INFN-CNAF (Italy), NL-T1 (The Netherlands), PIC (Spain), ASGC (Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers.Major contributors of computing resources are listed in Ref. [103].

Data Availability Statement
This manuscript has no associated data or the data will not be deposited.[Authors' comment: All ATLAS scientific output is published in journals, and preliminary results are made available in Conference Notes.All are openly available, without restriction on use by external parties beyond copyright law and the standard conditions agreed by CERN.Data associated with journal publications are also made available: tables and data from plots (e.g.cross section values, likelihood profiles, selection efficiencies, cross section limits, …) are stored in appropriate repositories such as HEPDATA (http:// hepdata.cedar.ac.uk/).ATLAS also strives to make additional material related to the paper available that allows a reinterpretation of the data in the context of new theoretical models.For example, an extended encapsulation of the analysis is often provided for measurements in the framework of RIVET (http://rivet.hepforge.org/)."This information is taken from the ATLAS Data Access Policy, which is a public document that can be downloaded from http://opendata.
-jets ) and the b-tagged jets (N b-jets ).The distributions of N add.-jets and N b-jets are shown for the signals and the background in Fig. 3, demonstrating clear

Fig. 2
Fig. 2 The m JJ distributions from the top-philic Z signals compared to the total background after the event preselection.The signal distributions are shown for two values of the chirality angle, (a) θ = 0 and (b) θ = π/2.In each case, four signal samples with different mass resolutions are presented for the possible combinations of m Z = 1.5 TeV,