Signatures of vector-like top partners decaying into new neutral scalar or pseudoscalar bosons

We explore the phenomenology of models containing one Vector-Like Quark (VLQ), $t'$, which can decay into the Standard Model (SM) top quark, $t$, and a new spin-0 neutral boson, $S$, the latter being either a scalar or pseudoscalar state. We parametrise the underlying interactions in terms of a simplified model which enables us to capture possible Beyond the SM (BSM) scenarios. We discuss in particular three such scenarios: one where the SM state is supplemented by an additional scalar, one which builds upon a 2-Higgs Doublet Model (2HDM) framework and another which realises a Composite Higgs Model (CHM) through partial compositeness. Such exotic decays of the $t'$ can be competitive with decays into SM particles, leading to new possible discovery channels at the Large Hadron Collider (LHC). Assuming $t'$ pair production via strong interactions, we design signal regions optimised for one $t'\rightarrow S t$ transition (while being inclusive on the other \bar{t'} decay, and vice versa), followed by the decay of $S$ into the two very clean experimental signatures $S\rightarrow \gamma \gamma$ and $S\rightarrow Z(\rightarrow \ell^+\ell^-)\gamma$. We perform a dedicated signal-to-background analysis in both channels, by using Monte Carlo (MC) event simulations modelling the dynamics from the proton-proton to the detector level. Under the assumption of BR$(t' \rightarrow S t) = 100\%$, we are therefore able to realistically quantify the sensitivity of the LHC to both the $t'$ and $S$ masses, assuming both current and foreseen luminosities. This approach paves the way for the LHC experiments to surpass current VLQ search strategies based solely on $t'$ decays into SM bosons ($W^\pm, Z$, $h$).


Introduction
During Run II at the LHC, the ATLAS and CMS experiments have collected almost 150 fb −1 and 180 fb −1 of data, respectively, at a centre-of-mass (CM) energy of 13 TeV. These data are now being analysed by the collaborations and, so far, no significant deviations from the SM have been recorded. This has significantly restricted the parameter space of the most common scenarios attempting to solve the hierarchy problem of the SM, such as supersymmetry and compositeness. Yet, it is important to find a viable solution to this flaw of the SM. This is inevitably connected to studying both top quark and Higgs boson dynamics, as the hierarchy problem of the SM originates from their mutual interactions. A pragmatic approach is to investigate BSM scenarios in which either of or both the top and Higgs sectors of the SM are enlarged through the presence of companions to the SM states (t and h), by which we mean additional spin-1/2 and spin-0 states, respectively, with the same electromagnetic (EM) charge but different mass (naturally heavier) and possibly different quantum numbers as well. Some guidance in exploring the various BSM possibilities in this respect is afforded by experimental measurements of observables where both the top quark and the SM-like Higgs boson enter. On the one hand, a sequential fourth family of chiral SM quarks is strongly constrained indirectly from Higgs boson measurements due to their non-decoupling properties [1], while VLQs (which transform as triplets under colour but whose left-and right-handed components have identical electroweak (EW) quantum numbers) can evade these bounds easily. On the other hand, the possibility of the existence of additional Higgs bosons has not been excluded by experimental data and may well be theoretically motivated by the fact that neither the matter nor the gauge sectors are minimal. Moreover, the Higgs sector is extended in any supersymmetric model or in the 2HDM.
Similarly, any model in which a Higgs boson arises as a pseudo-Nambu-Goldstone Boson (pNGB), other than the minimal model based on the symmetry breaking pattern SO(5)/SO(4), will include additional light (pseudo)scalars that might well have eluded direct searches due to their reduced couplings to the EW bosons and top quark.
Hence, it is of some relevance to assess the viability at the LHC of BSM models with both top quark partners (of VLQ nature) and companion scalar or pseudoscalar particles (both charged and neutral). In fact, it is particularly intriguing to investigate the possibility of isolating experimental signatures where the two particle species interact with each other, namely, when the t decays into a new (pseudo)scalar.
It is thus important to ask how the presence of exotic decay channels of VLQs can affect the current bounds and whether these might actually be promising discovery channels on their own. This question has been asked in similar contexts in various preceding works [14][15][16][17][18][19][20][21][22][23], each concentrating on a specific BSM construction. Here, in contrast, we follow the approach of [24], which adopts a set of simplified scenarios based on effective Lagrangians (motivated by compositeness).
In our paper, we build upon this last work, by adopting a simplified scenario which contains, above and beyond the SM particle spectrum, a top-like VLQ, t , as well as an additional scalar (or pseudoscalar) particle, S, in turn leading to the new decay channel t → S t. As for the decay modes of S, we will concentrate on two of the experimentally cleanest channels accessible at the LHC, namely, S → γ γ and S → Z γ, with the Z boson decaying in turn into electrons or muons. We will show in section 2 that there exist well motivated phenomenological scenarios where these can indeed be decay modes with significant BRs, for the case of both fundamental and composite Higgs states. In section 3 we estimate LHC constraints using published ATLAS and CMS searches in γ γ and Z γ final states while in section 4 we will describe our MC simulations, based on the pair production process p p → t t , followed by the decay chains t → S(→ γ γ) t or t → S(→ Z γ) t, with thet treated inclusively (and vice versa). Section 5 is then dedicated to interpreting the ensuing MC results in three theoretical scenarios embedding a t alongside additional (pseudo)scalar states focusing on cases with BR(t → S t) = 100%, while in section 6 we conclude.

The simplified model
The purpose of this section is to present the relevant details about the class of models whose phenomenology we aim to study. We begin with a general description of a simplified model that captures all relevant features. This is the model used for the analysis in section 4. We then justify the use of this simplified model by introducing three more specific models that can all be described with the same generic Lagrangian by a mapping of the fields and the couplings, provided that the processes considered in this paper are studied.
As discussed in the introduction, we are interested in exotic decays of a top partner t (of mass m t ) into the ordinary top quark t and a scalar (or pseudoscalar) generically denoted by S (of mass m S ) in the simplified model. We can thus augment the SM Lagrangian L SM by the following interaction Lagrangian with operators up to dimension five involving these two additional fields, Here κ S L and κ S R are the Yukawa couplings of the S to the t and t . In the second line, f sums over all SM fermions (including the top t) and κ f is the dimensionless reduced Yukawa coupling. In the last line V µν denotes the field strengths of the U(1) Y , SU(2) L and SU(3) C gauge bosons B µ , W µ , G µ in the gauge eigenbasis, g V is the associated gauge coupling (g , g, g s respectively) and V µν = (1/2) µνργ V ργ is the dual field strength tensor. The coefficientsκ V and κ V are couplings associated with dimension-five operators and are typically generated by loops of heavy particles or via anomalies. The couplings λ V for any gauge boson V are only generated if S is charged under some of the SM gauge groups and gets a vacuum expectation value (VEV) or if it mixes with such states, e.g., the Higgs boson. Since SU(3) C and U(1) EM are unbroken for the strong and EM interactions, λ V = 0 for the respective gauge bosons. We choose to normalise all terms with only one dimensionful parameter, the VEV v = 246 GeV.
In practice, we consider an S state of either scalar or pseudoscalar nature, but not a mixture. We therefore do not consider CP-violation in this paper. This means that either κ V orκ f are zero, in the scalar case, or κ V , λ V and κ f are zero, in the pseudoscalar case.
The total widths of t and S are kept as free parameters in the simulation as an indication that other interactions and other states might be present. These interactions are not explicitly required to describe the process p p → t t → SStt apart from their contribution to the total widths. Here we only report the analytic expression for the partial width of the exotic t decay, specifically.
where x t ≡ m t /m t and x S ≡ m S /m t . This formula is valid for decays into both scalar and pseudoscalar S. This defines the simplified model that will be used in the rest of this paper. Let us now briefly discuss three specific examples of models that motivate the use of the above simplified model and the mapping between the former and the latter. The results in this paper, given in terms of the simplified model above, can then easily be reinterpreted in terms of each model, if needed. In a forthcoming paper, we will specify these models in more detail and will discuss their specific phenomenology.

Example 1: adding a VLQ and a scalar to the SM
In order to illustrate how a particular model can be related to the phenomenological simplified model (eq. (2.1)), we will first present a simple model of top-quark partial compositeness (PC) in some detail. The model consists of the SM extended by a top partner VLQ and a scalar singlet. In this model the top quark acquires its mass via the mixing with the top partner. This model is not intended as a complete, realistic model, but provides an example of a model with an additional scalar S that is neutral under the SM gauge group. We will only be concerned with the couplings between the top quarks and S, leaving the coupling inducing the decay of the S to SM states as in eq. (2.1).
We denote the gauge eigenstates in the top sector by t L , t R and T . The notation t L/R is to prevent confusion with the mass eigenstates that are to be denoted by t and t . The Lagrangian for this model before EW symmetry breaking (EWSB) can be written as where the SM Higgs doublet is denoted by H with H = iσ 2 H * . The SM Yukawa coupling for the top quark is here denoted by y and Q L is the left-handed quark doublet of the third generation. The couplings λ a,b S are real if S is a scalar and purely imaginary if S is a pseudoscalar. The mass m 2 is a non-diagonal entry in the mass matrix of eq. (2.5). The remaining couplings are dimensionless. After EWSB, we have a mass matrix where we defined m t = yv/ √ 2 and m 1 = λ 1 v/ √ 2. The mass matrix can be diagonalised by bi-orthogonal rotations by the angles θ L,R , separately for left-and right-handed fermions, as follows (where s X ≡ sin θ X and c X ≡ cos θ X ) where {t, t } are the mass eigenstates and the mixing angles are given by The mass eigenvalues m t and m t are found by computing the eigenvalues. This model can be mapped to the simplified model Lagrangian in eq. (2.1) by performing the rotation in eq. (2.6) inside eq. (2.4). Focusing on the mixing terms yields while for the coupling to the top we have There is also a diagonal term involving the t , which is proportional to λ a S c L c R + λ b S c L s R . It is not included in the simplified model, but instead generates a contribution to the effective coefficients κ V and κ V from loop diagrams.
Let us also briefly discuss the decays of the t and S in this model. The t has both the standard and non-standard decay channels discussed above, where the width of the t → S t channel is given by eq. (2.2) with the couplings defined in eq. (2.8). The scalar can, in general, decay into the final states gg, γγ, Zγ, ZZ, W W and tt. We always assume m S < m t , which forbids the decay S → t t . Apart from the tt channel, all the other decays are generated by loops of the t and t .
We may now examine the decay of the t and S depending on the coupling of T L with T R and t R . The t → S t decay is induced by the λ a S and λ b S couplings. If we are interested in a large BR(t → St), we may achieve that easily in a wide region of parameter space by considering suitable values of these couplings. For example, when the T L couples to t R (i.e., S is sufficiently large, as the partial width is proportional to (λ a S /λ t ) 2 . However, this will also increase the s-channel production of S through gg fusion,  therefore, this scenario is heavily constrained by the gg → S → γγ resonance search data from the LHC. In figure 1, we show the BRs of t for a specific benchmark point where the t → S t channel has a BR of almost 100%. As for the S decay, the S → gg channel dominates if the tt decay is not kinematically allowed, m S < 2m t . The total decay width is governed by Γ S→gg , and hence the branching ratio in the γγ channel is approximately Despite the small BR, the S → γγ decay is a clean and well motivated channel. For instance, in the search for a VLQ decaying into a Higgs boson and a top, t → ht, the h → γγ decay channel (which has a BR of 0.23%) is still sensitive [25]. We also note that there is no dedicated di-jet search, t → St → ggt, although it has been recently proposed in ref. [26]. The current bounds estimated by a recast of R-parity violating (RPV) supersymmetry searches [27] are not competitive. Other loop induced channels are more suppressed than S → γγ. For example, the partial widths of S → Z γ and S → ZZ, modulo negligible m Z corrections, are 2 tan 2 θ W Γ S→γγ and tan 4 θ W Γ S→γγ , respectively. For m S > 2m t , the tree-level S → tt channel usually dominates over the loop induced decays. However, in a region of parameter space, the tt decay can be tuned down by suitable values of the off-diagonal entries in the mass matrix in eq. (2.5). We find that, when sin θ L sin θ R (or equivalently m 1 m 2 ), the effective Stt coupling, depending on the λ a S and λ b S couplings, is not sufficiently large to compete with the loop induced decays of S. The six tops final state via t → S t → ttt has been discussed in ref. [22] with both a recast from current searches and a dedicated analysis.

Example 2: adding a VLQ to the 2HDM
The 2HDM (see [28] for a review) is widely used as a minimal model for an extended Higgs sector that goes beyond additional singlet scalars. With additional vector-like top partners (see [29][30][31] for previous work), the 2HDM may be seen as the low-energy manifestation of a composite Higgs scenario, such as in [32]. Specifically, we here consider a vector-like top partner T with charge +2/3 in the singlet representation of the SM EW group. We further consider Yukawa couplings of the SM quarks of Type-II, i.e., that the up-and down-type quarks couple to different doublets.
The Higgs sector of the 2HDM has an additional neutral scalar H, a pseudoscalar A and a charged H ± state. This enables us to obtain simple formulae where either H or A can play the role of S in the simplified model Lagrangian in eq. (2.1). The details of the model and the involved parameters as well as the mapping onto the simplified model Lagrangian of eq. (2.1) are discussed in appendix A.1. Let us here only discuss the mixing of the physical top quark t and top partner t .
The physical mass of the heavy top, m t , is different from the mass M of the vector-like T due to t-T mixing. The mass matrix can be diagonalised in the same way as in eq. (2.6) to obtain the physical states (t L,R , t L,R ) in terms of the gauge eigenstates ( t L,R , T L,R ). The mixing angles θ L and θ R are not independent parameters and we can derive similar relations to eq. (2.6) (see eq. (A.7)), in terms of the Yukawa couplings y t and ξ T that couple the left-handed quark doublet Q L to the right-handed SM top t R and the vector-like T R , respectively (see eq. (A.4) and eq. (A.14)). The two mixing angles in this case satisfy [31] tan while the mass of the t is related to the Lagrangian parameters and the physical top quark mass via . (2.12) The t -t interaction can thus be described by three independent physical parameters: two quark masses, m t and m t , and a mixing angle, s L = sin θ L .
In the 2HDM with a VLQ, the scalar S is an additional Higgs boson. The dimension-five operators in eq. (2.1) are then generated through loops and in general S can be produced through gg → S. It can then decay in all the bosonic channels that we consider in this paper and, in addition, in fermionic ones. (The BRs in this model are discussed in section 5.) These channels give rise to constraints from all the usual collider observables. In addition, the scalar sector of this model is subject to the same unitarity, perturbativity and vacuum stability constraints as the usual 2HDM [28,33]. The Yukawa coupling y t is constrained from unitarity to be less than 4π, while ξ T is a derived quantity. Since the new top partner will contribute to gauge boson self energies, the mixing angle θ L can be constrained from EW Precision Tests (EWPTs) such as the S and T parameters. Based on ref. [31], such bounds require the mixing angle θ L to be in the range (−0.15, +0.15). However, the constraints coming from BR(b → sγ) are the most relevant ones, as the mixing angle is restricted to be in the range (−0.1, +0.1) for large m t , i.e., around 1 TeV.

Example 3: realisation in partial compositeness
Lastly, we present a Composite Higgs Model (CHM), which motivates the analysis in this paper by having a top partner with enhanced exotic decay mode and a pseudoscalar with dominant Z γ decay. The model is closely related to one of the earliest non-minimal models of composite Higgs with fermionic partial compositeness [34], based on the coset space SU(4)/Sp(4), where Sp is the symplectic group. The usual Higgs field H is a bi-doublet of SU(2) L × SU(2) R , which together with a singlet S (usually denoted by η in the CHM literature) forms the five dimensional anti-symmetric irreducible representation of Sp(4), This scenario has the further appeal of belonging to a class of models that can be obtained from an underlying gauge theory with fermionic matter [35,36] and the additional features arising from this fact have been studied in, e.g., [37]. Here, however, we want to focus on the bare bones of the model, namely the above-mentioned coset structure with the addition of one fermionic partner Ψ. (We only consider partial compositeness in the top sector).
The fermionic sector also consists of a bi-doublet and a singlet in the 5 of Sp(4). We will see that, as already anticipated in [24] (see also [14]), the possible decay patterns of the fermionic partners are richer than what is usually considered in current searches and, in particular, the lightest top-partner has an enhanced decay into the exotic channel t → S t.
To summarise, in addition to the SM fields the model has an additional pseudoscalar S, three top partners T, T , T (all of electric charge +2/3), a bottom partner B (charge −1/3) and an additional coloured fermion X of charge +5/3. Like in the previous example models, all of these fermions are vector-like Dirac spinors, to be thought of as in the gauge eigenbasis, i.e., before their mass matrices are diagonalised. The difference here is that there are more than one new fermion.
The mixing with the third family quarks of the SM depends on how they are embedded in a representation of SU(4). We choose this embedding such that the custodial symmetry of [38] is preserved, see appendix A.2 for details. In addition, the choice of having an elementary t R distinguishes this model from similar ones studied in [39], where the t R was taken to be fully composite. The elementary t R seems more appealing, since chiral fermions are notoriously difficult to obtain from underlying strongly coupled theories. We do not address the origin of the bottom quark mass in this work, which would add additional model dependence that is not relevant for the experimental signatures of interest. See appendix A.2 for more details on the construction of the model and the singular value decomposition of the mass matrix.
We end up with four top quark mass eigenstates, which we denote, in increasing mass order, by t, t , t and t . Here t is the known SM top quark of mass m t = 173 GeV. We diagonalise the mass matrix numerically, but a perturbative expansion for the masses gives some insight into the mass spectrum. We find (see appendix A.2) where M is the mass parameter of the Ψ, y L and y R are the respective couplings of the Q L and t R to the Ψ and pNGBs while f is the "pion decay constant" of the strongly coupled theory. We also defined M = M 2 + y 2 L f 2 . The mass of the bottom partner (mostly aligned with B) turns out to be of the same order as that of the heaviest top partner m t , while X has mass equal to M ≡ m t since it does not mix with anything.
Substituting the mass eigenstates (see appendix A.2) into the Lagrangian and considering the coupling that mixes the two lightest eigenstates t and t with the pNGBs, we see that no mixing with the Higgs field h arises, while the S couples, up to terms of order (2.16) allowing us to match the models with the parameters of the phenomenological Lagrangian eq. (2.1) From the analysis of the spectrum and of the couplings, we see that we can concentrate on a model with two mass degenerate VLQs t and X, with ∼ 100% branching ratios X → W + t and t → S t. The decay modes of t to SM vector bosons are highly suppressed, t being a singlet of SU(2) L × SU(2) R . For this model, it is thus crucial to understand whether the BSM decay t → S t can compete with the SM decay X → W + t whose signatures have been looked for at the LHC [6] providing bounds to the model parameter M > 1.2 TeV. We address this question in this work. Just above the t mass scale there is a further top partner, t , with more diverse and model dependent decay modes, so it is likely to be less relevant to experimental searches. The last top partner t and the B are heavy and can be ignored altogether. The coupling of the S to gauge bosons can be motivated by the analysis of the underlying gauge theory [35,36] and is given at leading order by the Lagrangian where the "Abelian" field strength tensors are defined as V µν = ∂ µ V ν − ∂ ν V µ , thus omitting the "non-Abelian" part, which would contribute to interactions with three and four gauge bosons that we ignore here. A is a model dependent dimensionless anomaly coefficient: For instance, in the model analysed in [24] A is given by the dimension of the representation of the hyper-fermions. Note that there are no couplings of type SSV since the S does not acquire a VEV. Also, there is no anomalous coupling SF µν F µν to the EM field, thus the decay S → γ γ is highly suppressed and for m S 2m W the decay S → Z γ has near 100% branching ratio. Once again, we can match the current model with the remaining couplings of the phenomenological Lagrangian in eq. (2.1): The mass of S is expected to be small m S m h and thus in the region where the decay into Z γ is motivated. In this particular model, it is given by m h /(2 cos θ) plus corrections proportional to explicit underlying fermions masses, which are disfavoured by fine tuning arguments. For t R symmetric for example, m η tends to vanish and should get its mass completely from underlying fermion masses. Other representations and other models give different expressions, but all agree on the approximate estimate that m S is light due to its pNGB nature.
As far as direct S production goes, we observe that, choosing the spurion embeddings as above, no diagonal coupling of type St i t i (t i = t, t , t , t ) is directly generated [34]. This means that the gluon fusion process is not present and the direct production proceeds mainly via EW vector bosons. Diagonal fermionic couplings for the top and for lighter fermions can be induced by further enlarging the model but we ignore them and consider the fermiophobic case. The coupling of S to fermions is nevertheless generated via loop of gauge bosons and might be relevant for low m S [40,41].

LHC constraints from γγ and Zγ resonance searches
To perform a phenomenological analysis of the γγ and Z γ final states it is necessary to estimate the allowed regions in the masses of the VLQ and (pseudo)scalar. This is done in this section by recasting one ATLAS and one CMS search at 13 TeV and providing the ensuing limits in the m t vs m S plane.
The searches used for the recast are briefly described in the following.
• An ATLAS "Search for new phenomena in high-mass diphoton final states" [42], used to set constraints for the γγ final state. This search looks for resonances with spin 0 or 2 decaying into two photons. For the spin 0 resonances (of interest for our analysis) the explored diphoton invariant mass region ranges from 200 GeV to 2700 GeV. The search cuts on the transverse energy of the leading and subleading identified photons, E T > 40 GeV and E T > 30 GeV, respectively, and requires E T to be larger than a fraction of the diphoton invariant mass, E T > 0.4m γγ GeV (leading photon) and E T > 0.3m γγ GeV (subleading photon).
• A CMS "Search for standard model production of four top quarks with same-sign and multilepton final states" [43], used to set constraints for the Z γ final state. This search looks for final states with two (same-charge) or three leptons, and different numbers of jets and b-jets, depending on the signal region. No cuts are imposed on photons in the final state. The most relevant cuts are applied to the jet and b-jet multiplicity and differ depending on the signal region.   [42] and CMS search [43], respectively. The solid black lines represents the bounds on the two masses obtained by comparing the upper limits with the pair production cross section of t at NLO+NNLL computed through Hathor [52] under the assumption of 100% BRs for both t and S in the respective channels and in the narrow width approximation (NWA).
The recast simulations are done using MadGraph5_aMC@NLO [44] with a dedicated UFO [45] model file corresponding to the simplified Lagrangian in eq. (2.1). Events are generated at leading order and interfaced with Pythia 8.2 [46] and Delphes 3 [47] for showering and fast detector simulation. As Parton Distribution Functions (PDFs), the NNPDF 3.1 at NLO set [48] has been chosen, obtained through the LHAPDF 6 library [49] using PDF ID 303400. The recast and validation of the searches is then performed through MadAnalysis 5 [50,51].
Simulations have been performed in a grid of t and S masses: m t has been varied in the range 400 GeV to 1000 GeV in steps of 100 GeV, while m S starts from a minimum value of 200 GeV and increases in steps of 100 GeV until reaching the kinematical limit m t − m S − m t = 0. A point in the small mass gap region m t − m S − m t = 10 GeV has been included as well.
The results are shown in figure 2 as upper limits on the cross section (in pb). The observed bound on the t and S masses, represented as a solid black contour, has been obtained by comparing the upper bounds on the cross section with the cross section for pair production of t obtained at NLO+NNLL through Hathor [52], under the assumption of 100% BR for t → S t and for S → γ γ (figure 2 left panel) or Z γ (figure 2 right panel) in the narrow width approximation (NWA). The range of validity of the NWA in terms of the ratio between the total width and mass of t is discussed in appendix B. In the γ γ channels the allowed region for m t is above ∼ 600 GeV almost independently of m S . In the Z γ channel the bounds are slightly more sensitive to the mass gap between the VLQ and the (pseudo)scalar, barring statistical fluctuations: the bound on m t is however between ∼ 700 GeV and ∼ 800 GeV for all the allowed m S .
The bounds obtained are typically weak compared to dedicated VLQ searches. We stress, however, that the bounds provided in this section are simply meant to give an idea about the optimal sensitivity of current searches for the final states considered above. In realistic scenarios the BRs of t and S into such final states will be likely smaller than 100%, which trivially implies that the bounds will get weaker. In this case, other channels might be more sensitive depending on the BRs of the t (and the recasting of different searches more sensitive to other final states has been performed, e.g. in [26], after the appearance of this analysis). Indeed, only a combination of bounds from different final states would give a full picture for any given benchmark point (defined in terms of masses and BRs of t and S). The way bounds are provided in figure 2, however, represents one of the elements of this picture. As a practical example, if a benchmark is considered in which the BRs of t → S t or S → γ γ or Z γ are smaller than 100%, the observed upper limits on the cross section represented by the grid of numbers in figure 2 can be directly compared with the σ× BRs of a given benchmark to determine the corresponding bound.
In the next section we propose a dedicated analysis to look for the signatures we are interested in leading to a much better sensitivity than the ones presented in figure 2.

Analysis
In its full generality, a top partner t may decay into the usual three SM channels W + b, Z t, h t or additional exotic channels. In this paper we are focusing our attention on the case of pair production p p → t t and subsequent decay into the BSM channels t → S t, where S is a neutral (pseudo)scalar decaying into SM EW diboson pairs. We have chosen the decays S → γ γ and S → Z γ as our target signal, since they are experimentally very clean bosonic decay channels. In the case of the Z γ channel we only consider further leptonic decays of the Z.
The analyses are optimised to look for only one pair of photons or Zγ final states originating from the same S. When limits from these analyses are reinterpreted in specific models, the BRs of the S can significantly affect the limits therein. In order to reinterpret the results in the models described in section 2, we need to evaluate the efficiencies of the signal region cuts while taking into consideration all possible decays of S. We assume t decays at 100% rate as t → S t. For S, we consider all the possible bosonic decay channels necessary to ensure gauge invariance in the CHM, 1 In this section we briefly define the objects used in the analyses (with a longer discussion for reproducibility in appendix C), then describe the tools and processes for the simulation of events to model signal and background (section 4.2), and finally we present event selections Figure 3: Pair production of t with decay of the t into (anti-)top and S in both branches. S is then decayed in one branch into γ γ or Z γ, depending on the signal pursued, and inclusively in the other branch.
to extract the signal in the two considered signal regions (SR): the γ γ SR in section 4.3 and the Z γ SR in section 4.4.

Object definition
In the following the definition and selection of objects at reconstructed level are briefly outlined. A more detailed account can be found in appendix C. The default ATLAS Delphes card [47] is used, with minor modifications and calorimeter objects that fall in the calorimeter transition region 1.37 < |η| < 1.53 are excluded. Isolation and overlap removal is done in the Delphes card for most of the objects. The basic objects used are photons (γ), leptons ( ), jets (j) and b-jets (j b ). Photons are required to have a p T > 30 GeV and |η| < 2.37. Leptons in this paper are understood to mean electrons or muons only, and not τ -leptons. Leptons must fulfil p T > 25 GeV and |η| < 2.47. Jets are reconstructed by using the FastJet [53] package and Delphes with the anti-k t algorithm [54] using R = 0.4. Jets are required to pass p T > 25 GeV and |η| < 2.47. In Delphes, a b-jets is a jet which contains a truth b-quark.
The compound objects used are Z bosons, missing transverse energy (E miss T ) and the scalar transverse energy (H T ). Z bosons are identified as two opposite-sign same-flavour [47], where i runs over the energy deposits in the calorimeter. H T is the scalar sum of the p T of all reconstructed basic objects used in the analysis (jets, muons, electrons and photons).

Simulations
All simulations in this study have been performed using the following framework: Mad-Graph 5_aMC@NLO [44] was used to generate events at leading order accuracy. Pythia 8.2 [46] and Delphes 3 [47] have been used for showering and fast detector simulation, respectively. For the signal simulations, the parton distribution function (PDF) NNPDF 3.1 at NLO set [48] set has been chosen, obtained through the LHPDF 6 library [49] using PDF ID 303400. For the background simulations instead the MadGraph default NNPDF 2.3 LO with PDF ID 230000 has been used. Hathor [52], with NNLO MSTW2008 [55] PDFs.
The numerical values of the pair production cross-sections, which only depend on m t , are shown in figure 4. They were computed through Hathor [52], with NNLO MSTW2008 [55] PDFs.
The background of the γ γ SR is dominated by pp → γ γ + jets mediated by QCD interactions. The backgrounds γ γ+t+jets and γ γ+tt were found to be negligible and hence are not considered for the diphoton analysis. Events from the pp → γ γ + jets process are generated with up to three jets, including jets initiated by b-quarks, in the matrix element. The final jets after showering and jet clustering are matched to the original partons with the MLM method [56] as implemented in Pythia. In the simulation of the initial state bquarks are explicitly considered as part of the incoming protons. This accounts for processes with an odd number of b-jets in the final state, such as those initiated by gb → γ γ + uūb. To ensure enough statistics in the high mass tail the events are generated in slices of the diphoton invariant mass M bkg γγ with ∼ 1 M events per slice, where M bkg γγ refers to the invariant mass of the generated (not reconstructed) photons. Table 1 lists the slices along with the fiducial cross section for each slice. The invariant mass of the two photons for all slices is shown in figure 5. If there are more than two photons in the event, the pair with invariant mass closer to 160 GeV is shown in this figure. The high-mass slices have small tails towards lower masses, which occurs when one or both of the hard photons is lost in the reconstruction and the selected photons originate from e.g. the hadronisation process. The contribution from these mis-reconstructions is typically small and can be mitigated 74.0 ± 0.6 4.45 ± 0.03 Table 1: Fiducial cross section for each mass slice of the two major background processes. For the γ γ + jets background the slices refer to M bkg γγ while for the Z γ + jets background the slices refer to M bkg Zγ at the generator level. The sums of the fiducial cross sections over all slices for each process are also listed together with their estimated value.
further with ∆R cuts on the photons. The small peak at 160 GeV is due to the selection requirement that the invariant mass of the photons is close to 160 GeV. The total fiducial cross section in the M bkg γγ > 50 GeV region is calculated by generating 25K events in the allowed range using the same setup as in the full event generation, resulting in 74.0 pb, in good agreement with the sum of the fiducial cross sections for the individual slices.
The dominant background in the S → Z γ final state is pp → Z γ+jets, with Z → + − . Events from this process are generated using the same setup as for the γ γ +jets background, with up to two hard jets in the matrix elements. For the same reason as for γ γ + jets the event generation for the Zγ + jets background is performed in slices of the invariant mass of the generator-level Z and γ, M bkg Zγ , with ∼ 2M events each, listed in table 1 together with their fiducial cross section. The latter at M bkg Zγ > 50 GeV is estimated to be 4.451 pb by generating 25K events in the allowed kinematic range, which, again, is in good agreement with the sum of the fiducial cross sections of the slices. SM top-quark pair production associated to a photon and to a Z and a photon can also give relevant contributions to the background. We generated 150K events of the process tt + Zγ and let the top decay inclusively and the Z leptonically via MadSpin. For tt + γ we generated 300K events and required the top quarks to decay leptonically to either electrons or muons. We use the LO cross sections 0.315 fb for decayed tt + Z + γ and 94 fb for decayed tt + γ events. The invariant mass of the Z γ system, for each of the mass slices of Z γ + jets, together with tt + γ and tt + Z + γ, is shown in figure 6. In that figure, at least one Z boson, one photon and one b-jet, according to the definitions in section 4.1, are required. If there are more than one Z and/or γ candidate we choose the system with invariant mass closer to 160 GeV to present in this specific plot.
In both final states, non-prompt backgrounds are also possible. These are expected to be reduced significantly since we use tight identification requirements for leptons and photons. Furthermore, in analyses with similar final states, the backgrounds with one or more jets mis-identified as photons was found to be significantly smaller than those with prompt photons [57]. Thus, we do not consider non-prompt background sources in either For the signal simulation and definition, we generated the process pp → t t with t → S t and S decaying into EW bosons, eq. (4.1). We define our signal samples as any possible decay combination, (S → X)(S → Y ) where X, Y ∈ {γ γ, Z γ, W W, ZZ}. Both the Z and W decay inclusively in our signal definition.
The UFO model for signal simulations is the same one used for recasting LHC bounds, corresponding to the simplified Lagrangian of eq. (2.1). Decays of interest are thus turned on or off by setting the corresponding couplings. In the following analysis, couplings are set such that the widths for the top partner t and scalar S are 0.1% of their mass, to allow the use of the NWA. A quantitative determination of this parameter, performed in appendix B, is essential to determine the range of validity of signal simulations in experimental analyses and also for the subsequent reinterpretation of results in terms of theoretical models.
For the simulations, we use κ R S = 0, keeping only the κ L S coupling. This is an important assumption, as fixing a different chirality of the top coupling can lead to observable differences. Indeed, it is known that the dominant chirality of the couplings of a VLQ interacting with the SM top quark can be probed by looking at the transverse momentum of the decay products of the W boson emerging from the top quark [58,59]. Differently from the SM case, however, here the kinematics of the decay products of t is not only affected by its mass, but also by the S mass.
Similarly we turn off the scalar S couplings, κ W = κ B = λ W = λ Z = 0, when we assume a pseudoscalar nature of the S state. The scalar or pseudoscalar nature of S can also in principle affect the kinematical distributions of its decay products. We have therefore performed simulations imposing specific decay channels, to check, at reconstruction level but without including detector effects, how large differences can be between the above scenarios in differential distributions. We found that there is no observable difference in our predictions with respect to a scalar S in terms of kinematical distributions. In view of this indistinguishability, in the 2HDM+VLQ case, we will assume the S state to represent alternatively a CP-even and CP-odd neutral Higgs states entering the t decay. 1

S → γγ signal region
In this section, the diphoton final state is presented. From an experimental point of view, the diphoton final state gives a very clean signature in the detector, which makes it attractive to study. We considered t masses m t = 600 to 1800 GeV in steps of 200 GeV, every kinematically allowed S mass is investigated, via the discrete values of m S = 100 GeV, 200 GeV, 400 GeV, and then in steps of 200 GeV up to the highest kinematically available mass, m S = m t − 200 GeV. The wide selection of S and t masses enables the possibility to study both threshold effects and highly boosted decay products.
To select the signal we demand the presence of 2 photons and 1 b-jet defined according to section 4.1. If more than one pair of photons is present we choose the pair whose invariant mass is closer to m S and define these photons as "best" photon candidates, γ 1 , γ 2 . Unless otherwise specified, a pair of photons is assumed to be the "best" pair. The invariant mass of the system with the two "best" photon candidates is required to be within 20 GeV from the nominal S mass, |M γγ − m S | < 20 GeV.
In order to further enhance the signal discrimination with respect to the background for low m S values we use the fact that the S is produced in a boosted regime. The top partners t andt will be produced nearly at rest and the pair will be back-to-back. The large difference in mass between t and S will make S boosted and thus also the photon pair from S will be collimated. In figure 7 we show the ∆R γγ distributions for different m S and for m t = 800 GeV fixed. We take advantage of this characteristic signal profile and require ∆R γγ < 2.3 from m S = 100 GeV to m S = 200 GeV.
The selection cuts are summarised in table 2. Note that, due to limitations in statistics, the cuts are sub-optimal. The discrimination between signal and background could be improved significantly by tightening the cuts in a real experimental analysis.
In table 3 we show the efficiencies (number of events left after the cut divided by the number before the cut) of the selection cuts numbered in table 2 for different m S values. In the upper part of the table, the signal process is defined with both S decaying into diphotons, i.e., ttS(→ γ γ)S(→ γ γ) in the final state. This is the process we use to optimise the selection cuts. We display only the m t = 1 TeV case in the table. In the lower part of the table, the efficiencies for the background sample are displayed. It can be noticed that the last two cuts are the most efficient ones in removing the background and keeping signal events.    Table 4: Selection cuts applied to the Z γ signal region. Details on the cuts are given in the text. Refer to section 4.1 for the definition of the objects.
The final efficiencies for the signal decay channel S(→ γ γ)S(→ γ γ) are discussed in section 4.5. The efficiencies for the other signal decay channels with at least one branch decaying into γ γ are presented in appendix D.

S → Z γ signal region
In the S → Z γ final state we require at least one Z boson candidate reconstructed according to the definitions in section 4.1. In addition to the Z candidate we require the presence of at least one isolated photon. The system of one isolated photon and one Z candidate whose invariant mass is closest to the nominal S mass is called the "best S candidate". To efficiently distinguish the signal from the background we exploit the high multiplicity of objects and high total energy of a typical signal event. We require H T +E miss T > 0.3m t , where H T is the scalar sum of the p T of all reconstructed basic objects and E miss T is the missing transverse energy of the event as described in section 4.1. We finally require the invariant mass of the S candidate to be within 15 GeV of the nominal S mass, i.e., |M Zγ − m S | < 15GeV. A summary of these selection cuts is presented in table 4, with some information on the object definitions for convenience.
The distributions of M Zγ before cut 5 and H T + E miss T before cut 4 and 5 are shown in figure 8, for the masses m S = 160 GeV and m t = 1400 GeV. There is a great discriminating power in the H T + E miss T observable due to the large multiplicity and energy of a typical signal event. We note that the used cut is not optimised to suppress the background due to lack of MC statistics. A realistic experimental analysis could harden this cut to further reduce the background and use data-driven methods to estimate it without relying too much on MC estimates.
For illustrative purposes, in table 5, we show the efficiencies of the selection cuts numbered in table 4 for different m S values. We display only the case m t = 1400 GeV in the table. In the upper subtable, the signal process is defined with both S decaying into Z γ, S(→ Z γ)S(→ Z γ) in the final state. This is the process we use to optimise the selection cuts. In the lower subtable, the efficiencies for the background sample are displayed. Except the mass-window cut for the S candidates, all cuts depend on m t .

Efficiencies
The signal efficiencies for the two different signal regions are the last piece of information necessary for reconstructing the number of signal events. In figure 9 we provide, as il-   lustrative examples, the efficiencies for the (γ γ)(γ γ) channel in the γ γ SR and for the (Z γ)(Z γ) channel in the Z γ SR, for which the selections have been optimised. Further efficiency plots for different channels are provided in appendix D. All efficiencies have been computed considering signal samples of 10 4 MC events, corresponding to a statistical uncertainty of the order of 10% which can affect the evaluation of efficiencies especially when they are small. The whole set of efficiencies, combined with the BRs chosen in section 4, allows one to compute the expected total number of events via eq. (5.2) in the following section, where the results of the study are discussed.
In the next section we will show how to estimate the number of events for both signal and backgrounds for different model assumptions and devise a simple statistical framework for model interpretation.

Results
In this section we discuss the discovery potential of LHC for the models introduced previously. Essentially, we propose a counting experiment comparing the number of expected background events with the number of signal events. The expected number of background events in one of the signal regions SR ∈ {γ γ, Z γ}, B SR , is given by with L the integrated luminosity, and σ Bγ γ = 74.0 pb and σ B Z γ = 4.58 pb our best estimate of the total background cross section for the γ γ and Z γ signal regions, respectively, and B SR the efficiency after all cuts in the corresponding SR. The number of background events can be extracted for arbitrary values of m S and m t by interpolating the data presented in tables 6-7.
It should be noted that we only present the estimates for the irreducible background. This turns out to be negligible in the high mass region and its values are presented only to show this fact and for completeness. Fake rates are also expected to be negligible in the high-mass region [60].
The number of expected signal events for each SR is given by where Y,X SR is the final efficiency in appropriate signal region SR for the signal sample with decay (S → X)(S → Y ) with X, Y ∈ {γ γ, Z γ, W W, ZZ}. (In these expressions we assume the validity of the NWA and assume 100% BR t → S t andt → St.)    In appendix D we tabulate the above efficiencies, allowing one to estimate the signal in any of the theoretical models discussed here by simply computing the corresponding BR. The discovery potential for a more generic model can be also estimated using the numbers provided as long as the efficiency times BR of any extra decay channel is known to be small.
Having computed the number of signal (S) and background (B) events, we estimate the significance by employing the formula [61][62][63] that is obtained by using the "Asimov" data-set into the profile likelihood ratio. The explicit expression above, containing the uncertainty σ b on the background, is found in ref. [64]. We consider an overall σ b = 10%B systematic uncertainty on B. This number is most likely a conservative estimate and it is estimated by comparing the systematic uncertainties of ATLAS and CMS analyses with similar final states, especially high-mass Z γ searches [65,66] and high mass γ γ searches [67][68][69].

Model interpretation
Recall that the main focus is the study of models where the top partner has 100% BSM BR t → S t and S decays into EW gauge bosons. Even within this limited framework, we still need to discuss the relative strengths of the various S decay channels, controlled by the couplings in eq. (2.1). We start by considering the optimal reaches for the two SR considered in this analysis, corresponding to scenarios where S decays fully either into γ γ or Z γ. Such scenarios are likely non-physical, but they allow to determine the maximum potential of the selections. The LHC reaches for this simplified scenario are presented in figure 10 for two different LHC luminosities, corresponding to the final luminosity at the end of Run II and the nominal final luminosity of Run III. It can be noticed that the sensitivity of the search diminishes for increasing m t due to the reduction of production cross section, but it improves with increasing m S because of the reduction of the background yields (see table 7).
We now move on to more theoretically motivated scenarios. We first consider the benchmark motivated by partial compositeness, where only the anomaly induced pseudoscalar couplingsκ B andκ W are non-zero.
In this case, the structure of the anomaly coefficients [24] in all explicit realizations gives κ B +κ W = 0, thus suppressing the S → γ γ decay. This leads to a 100% BR(S → Z γ) below the W W threshold and still an acceptably large value above it, as displayed in figure 11 (left). The LHC reaches for this scenario are presented in figure 11 (right) for two different LHC luminosities, corresponding to the final luminosity at the end of Run II and the nominal final luminosity of Run III. Here, we consider only the Z γ SR because of the negligible sensitivity of the γ γ SR. Different effects are present in the reach of figure 11. For m S 2m W the sensitivity is optimal due to a 100% decay rate of both S into Z γ (S → Z γ, S → Z γ) and a high efficiency ( figure 9 (right)). Above threshold the S → V V, S → V V (V = W, Z) decay channels kick in with ≈ 64% rate and negligibible efficiency, while the S → Z γ, S → Z γ rate reduces to ≈ 4%. The mixed decay S → V V, S → Z γ takes 16% of the branching ratio and have an efficiency approximately constant and near 40% compared to the pure Z γ case ( figure 21). This depletion in the signal explains the kink of sensitivity lost near the m S ≈ 2m W threshold. In both regions the sensitivity improves with increasing values of m S due to a rapid decrease of the background, as noticed in figure 10.
The interpretation for the composite Higgs model described in section 2.3 is straightforward. The S is photophobic and we can read the bounds directly from figure 11. It is encouraging to see that even for not optimised cuts this channel could be competitive with the search for the +5/3 charged partner [11]. Some more details for this model are given in appendix A.
For the 2HDM+VLQ case, the interpretation is somewhat more complicated because of the more numerous parameters, richer particle spectrum and, hence, the decay patterns of t and S. A scan has been performed by varying the 2HDM input parameters described in appendix A.1 to obtain benchmark points characterised by the highest BRs of t and S into the final states considered in this analysis in order to maximise the sensitivity. Such points are simply representative of the 2HDM spectrum, as we ignored the fact that HiggsBounds excludes the majority of them. In fact, the scope of this selection is to illustrate the potential of the model independent analysis developed in this paper rather than to constrain specific theoretical models. We first restricted the scan by enforcing an almost exclusive decay of the t into the CP-even scalar H by setting the masses of the CP- The decay of H into γ γ is around 0.3% below the hh threshold, while its decay into Z γ is ∼0.05%. The generically dominant decay of H is into gg, which is on average around 70%, followed by W W (∼20%) and ZZ (5% to 10%), while the BRs into bb and cc are 1% or less. Above the 2m h ≈ 250 GeV threshold, H → hh dominates and all other BRs drop significantly, until the (on-shell) tt channel opens and becomes dominant. Then, we do a second scan with the role of H and A interchanged (approximately 80,000 points) and compute the BRs as described above. Here, there cannot be W W and ZZ decays of the A state, so that gg and bb decays share the majority oft the decay rate (about 90% of it, with the remainder saturated by τ + τ − and Zh, which we then neglect in the MC generation) till the (off-shell) tt * channel opens (coincidentally enough, around m t + m b + m W ≈ 260 GeV), with the γ γ and Z γ rates being generally lower than in the previous case. Given the low BR into Z γ, in addition to the subleading BR of the Z into leptons, no significant sensitivity is expected in the Z(→ + − ))γ final state and, therefore, we will focus only on the γ γ SR for the 2HDM+VLQ in the case of both a light H and A. This has been done only in the region of parameter space where high sensitivity is obtained, i.e., for m t less than 1 TeV, and are in average around 20%. The efficiencies for the {γ γ, hh} channel have also been calculated above the H → 2h threshold, in the region of high sensitivity, and are found to be around 30%. Given the illustrative nature of this example, we assumed the efficiencies for the {γ γ, tt * } in the case of a light A channel to be flat and around 30%.
The results for the 2HDM+VLQ are shown in figure 13. For the case of a light H state, some discovery reach has been found for m t around 600 GeV and exclusion is possible up to m t around 700 GeV, almost independently of m S below the (on-shell) tt threshold. For the case of a light A state, the reach in m t for both discovery and exclusion is somewhat deeper than in the previous case, by some 50 GeV. In contrast, the one in m S is very similar, as it again collapses at approximately 2m t . 2

Conclusions
While the case for VLQs, especially those of top flavour, has already been well established from the theoretical side, the experimental pursuit of their signatures at the LHC has been somewhat limited, as ATLAS and CMS analyses have primarily been carried out under the assumption that such new states of matter decay into SM particles only, i.e., via t → W + b, Z t and h t. This approach clearly enables one to make the most in terms of optimising the signal-to-background ratio in an analysis, chiefly because one can attempt reconstructing the measured W + , Z and h masses. However, if one considers VLQ models with additional particles this is overly restrictive since the VLQ may decay via exotic channels involving scalars or pseudoscalars. While the kinematic handles available to enhance these exotic channels may be apparently limited in comparison (as the exotic scalar or pseudoscalar states may have not been discovered already and/or their mass not measured), the size of the associated BRs could be large enough so as to nonetheless enable sensitivity to these channels. Furthermore, if the companion Higgs states are heavier than the W + , Z and h objects of the SM, the signal would anyhow be present in a region of space where the background contamination is minimised. Based on this reasoning, in this paper, we have set out to assess the scope of the LHC to test t decays into neutral (pseudo)scalar states, whose nature could be either fundamental or composite. As an example of spin-0 fundamental states, we have assumed here a Higgs sector comprised of the SM state supplemented by a scalar boson as well as a 2HDM (Type-II) containing both a scalar and pseudoscalar state (which we have taken light one at a time). As an example of spin-0 composite states, we have looked at a CHM where an additional pseudoscalar state emerges as a pNGB of the underlying new strong dynamics. In fact, we have also shown how all such models can conveniently be parametrised in the form of a simplified model onto which they can be mapped.
Of the various possible decay modes of this additional neutral (pseudo)scalar bosons, which we have collectively labelled as S, we have considered here two of the cleanest probes possible at the LHC, i.e., S → γ γ and Z γ (with the Z decaying into electron/muon pairs). In doing so, we have performed a dedicated signal-to-background analysis exploiting parton level event generation, QCD shower and hadronisation effects as well as detector emulation aimed at establishing the sensitivity of the LHC experiments to such decays, where the S state emerges from a companion top decay, t → S t, following t t production (with thē t decay treated inclusively). In the case of both S signatures, we have not attempted any reconstruction of the SM top quark entering the t decay chain although, on a trialand-error basis, we have assumed knowledge of the S mass, to be able to exploit both the cleanliness of the two S decay channels and the ability of a standard LHC detector in sampling γ γ and Z(→ + − )γ invariant masses with high resolution. Indeed, this approach also enables us to compare on a more equal footing the scope of t → S t signatures with that of t → W + b, Z t and h t ones, where a mass reconstruction is normally imposed on the W ± , Z and h decay products.
As a result of this approach, we have found that the t → S t signatures give a level of sensitivity not dissimilar from that obtained through studies of t → W + b, Z t and h t. For specific regions of the parameter space of VLQ models with exotic Higgs states, which have survived all available constraints from both direct and indirect t and S searches (including those obtained by ourselves from recasting experimental studies for other sectors), we have found the following exclusion and discovery reaches. For a simplified model maximising both the t and S BRs, m t can be probed in both the γ γ and Z γ channels up to approximately 2 TeV for S masses well into the TeV region. In the CHM scenario considered, coverage is not dissimilar for the γ γ case but for the Z γ the t reach is limited to 1.6 TeV. Finally, in the 2HDM+VLQ, it is possible to exclude m t up to around 700(750) GeV and discover m t up to around 600(650) GeV almost independently of m S when S represents the CPeven(odd) H(A) state and below the (on-shell) tt threshold for the decay of S. This is limited to the γ γ case, though, as Z γ gives no sensitivity at both Run II and III.
Hence, in connection to all of the above, we can confidently conclude to have surpassed the state-of-the-art in VLQ searches in two respects: firstly, by testing the scope of non-SM decays of the t state and, secondly, by deploying a selection procedure which is model independent yet enables one to interpret its results in a variety of theoretical scenarios. Furthermore, it should be noted that, while restricting ourselves to the case of γγ and Zγ signatures of the (pseudo)scalar states emerging from the described VLQ decays, there is no reason why our procedure cannot be applied to other S decays. Indeed, it can also be further improved (e.g., by reconstructing top-quark decays).
In summary, we believe that there is significant margin for improving the sensitivity of the LHC to models with a heavy top partner, through the exploitation of its decay channels into exotic (i.e., non-SM-like) neutral (pseudo)scalar states, which are ubiquitous in BSM constructs containing such a new fermion. In fact, over sizeable regions of the parameter space of the realistic VLQ models considered here, we have found that sensitivity to both the t and S mass can extend well into the TeV region, thereby being competitive with the currently studied SM channels. While in this paper we have limited ourselves to illustrating this through a few benchmarks examples, in a forthcoming paper, we shall quantify the regions of parameter space of our models where such a phenomenology can be realised, including tensioning the scope of standard and exotic t decays against each other.

A Details of the models
In this appendix, additional details are given of the models; the 2HDM+VLQ model in appendix A.1 and the composite Higgs model in appendix A.2.

A.1 The 2HDM with an additional VLQ
The scalar potential of the model includes two identical scalar doublets (Φ 1 , Φ 2 ) and a discrete symmetry Φ i → (−1) i Φ i (i = 1, 2), which is only violated softly by dimension-two terms [28], We take all parameters in the above potential to be real (although m 2 12 and λ 5 could in principle be complex). The two complex scalar doublets may be rotated into a basis where only one doublet acquires a VEV, the Higgs basis, where G 0 and G ± are the would-be Goldstone bosons and H ± are a pair of charged Higgs bosons. A is the CP odd pseudoscalar, which does not mix with the other neutral states. The Goldstone bosons are aligned with the VEV in Higgs flavor space, while the A is orthogonal. The physical CP even scalars h and H are mixtures of ϕ 0 1,2 and the scalar mixing is parametrized as where tan β = v 1 /v 2 is the angle used to rotate Φ 1,2 to the Higgs basis fields H 1,2 , α is the additional mixing angle needed to diagonalize the mass matrix of the CP-even scalars, and s β−α = sin(β − α), c β−α = cos(β − α). The most general renormalisable interaction and mass terms involving the VLQ can be described by the following Lagrangian (where we only include the third generation SM quarks), where H i ≡ iσ 2 H * i (i = 1, 2), Q L is the SM quark doublet and M is a bare mass term for the VLQ, which is unrelated to the Higgs mechanism of EWSB. Note that often the Yukawa couplings of the 2HDM are written in terms of the fields Φ 1 , Φ 2 . In eq. (A.4) we use the Higgs basis fields, so the Yukawa couplings y T , ξ T must be defined accordingly. In a Type II-model, as we are considering in this paper, the up-type quarks only couple to the doublet Φ 2 , while down-type quarks only couple to Φ 1 . Additional mixing terms of the form T L t R can always be rotated away and reabsorbed into the definitions of the Yukawa couplings. In the weak eigenstate basis ( t, T ), where t is the SM top quark, the top quark and VLQ mass matrix is where y t is the Yukawa coupling of the top quark. It is clear from the above mass matrix that the physical mass of the heavy top, m t , is different from M due to the t-T mixing.
The mass matrix M can be diagonalised by a bi-unitary transformation in the same way as in section 2.1 to obtain the physical states (t L,R , t L,R ) in terms of the gauge eigenstates ( t L,R , T L,R ), The mixing angles θ L and θ R are not independent parameters. From the bi-unitary transformations we can derive the relations and by using the traces and determinants we end up with the relations and a relationship between θ L and θ R and the Yukawa couplings, The t -t interaction can thus be described by three independent physical parameters: two quark masses m t , m t and a mixing angle s L = sin θ L . After rotating the weak eigenstates ( t L , T L ) into the mass eigenstates, the Yukawa Lagrangian takes the following form [31]: where U L,R are the matrices appearing in eq. (A.6). The neutral Higgs couplings to top (t) and top partner (t ) pairs are in the notation of eq. (2.1) given by (with S = H or A) The couplings here are normalised to y t / √ 2, which is what the H SM tt coupling would be in the case of no mixing between the t and T and additionally, in the alignment limit of the 2HDM, sin(β − α) → 1 where the lightest neutral scalar h is the SM-like Higgs boson. Note that in eq. (2.1) the terms with diagonal couplings St t of the top partner to the scalars are not included, since they are not phenomenologically relevant in this paper. We include them in eq. (A.15) for completeness, however. Note also that the combination (c β−α − s β−α cot β) that occurs in eq. (A.15) is proportional to the 2HDM Type II Yukawa coupling of the heavier Higgs boson H.
In our analysis we have used a modified version of the public code 2HDMC [70] with a VLQ added according to the description above. We have scanned over the parameter space of the model, which is constrained by Higgs data from the LHC that can be evaluated using the public code HiggsBounds [71]. In addition, 2HDMC can evaluate oblique parameters and theoretical constraints on unitarity, perturbativity and positivity of the potential. However, since our aim here is rather to demonstrate the use of the method developed in this paper, we have not made a comprehensive scan to satisfy these bounds, but instead we have considered parameter points that provide large BRs of t → St and S → γγ for S = H or A. We have therefore chosen to make the Higgs boson that does not play the role of S as well as the charged Higgs boson heavy. We perform random scans over the parameters and generate 10 5 points for each of the scenarios with S = H, A. We then keep those points where the product BR(t → S t) × BR(S → γγ) > 10 −3 . The scalar S is taken in the range 180 GeV < m S < 350 GeV, while for S = H the other heavy scalar is taken in the range 600 GeV < m A < 1000 GeV. For S = A, instead we choose m H = 1 TeV. The charged Higgs mass is always m H ± = 1 TeV. The remaining Higgs sector parameters are in the ranges 0.99 < |s β−α | < 1, 0.1 < tan β < 1 and we take m 2 12 = m 2 A sin β cos β. Finally, the VLQ couplings are taken in the ranges 500 GeV < m t < 1500 GeV, −0.15 < s L < 0.15 and 10 < y T < 15.

A.2 The composite Higgs model
As mentioned in the main text, the SM Higgs H field in this model is a bi-doublet of SU(2) L ×SU(2) R , which together with a singlet S forms the five dimensional anti-symmetric irrep of Sp (4), The fermionic sector also consists of a bi-doublet and a singlet in the 5 of Sp (4), The new fermions mix with the third family quarks of the SM. The mixing is obtained by choosing to embed both the left-handed Q L = ( t L , b L ) T and the right-handed t R as spurions into the 6 of SU(4). The non-zero components of Q L fit into the bi-doublet of the SU(2) L × SU(2) R subgroup, while t R is in the singlet of the 6 → 5 + 1 decomposition of SU(4) → Sp(4). The choice for Q L is essentially dictated by the need to preserve the custodial symmetry of [38]. The construction of the interaction Lagrangian from the general formalism has been addressed in many papers and will not be reviewed here. Suffice it to say that we combine the five pNGBs into a 4 × 4 matrix Π and exponentiate it to obtain Π , transforming as: Σ → gΣh −1 , for g ∈ SU(4), h ∈ Sp(4), (A. 18) and use it to "dress" the fermionic field Ψ, written as a 4 × 4 anti-symmetric matrix. In this notation, the Lagrangian becomes where we indicated the dressing explicitly. (Note that Q L → gQ L g T and Ψ → hΨh T .) We allow only the Higgs field to acquire a VEV, and we denote the mixing angle by GeV. Generically f > 800 GeV from EWPT, although one can envisage mechanisms that would allow to lower that bound [72].
Computing eq. (A. 19) to all orders in θ and retaining only terms linear in h and S, h being the canonically normalised physical Higgs with VEV shifted to zero, we can write the part of eq. (A.19) concerning top partners as (see also [24]) where the mass and Yukawa matrices are given by The singular value decomposition of M is unwieldy, but can be performed numerically or perturbatively to order θ ≈ v/f . For the four top quark mass eigenstates t, t , t , t , the perturbative expressions for the masses are 22) The mass of the bottom partner (mostly aligned with B) turns out to be of the same order as that of the heaviest top partner m t , while X has mass equal to M ≡ m t since it does not mix with anything. For the top quarks, the conversion from gauge to mass eigenbasis reads, to O (v/f ), This spectrum justifies that choice of simplified model in the text where we neglect all the top partners other than the lightest one.
Regarding decays of the pseudoscalar in this model, in figure 14 we show the partial widths of S as a function of its mass, including the dominant loop induced fermionic channel S →b b relevant below the Z γ threshold. We use f /(A cos θ) = 500 GeV but all curves rescale by (500 GeV A cos θ/f ) 2 . We see that for all interesting regions of parameters the width is always very narrow, but still prompt.
The most promising parameter region for this class of models is m S 160 GeV, where the S decays dominantly to Z γ. This region is motivated from the model building perspective since it is expected m S < m h . From the experimental point of view it offers   Figure 15: Values of the κ S L,R coupling corresponding to fixed Γ t /m t ratios (0.1%, 1% and 10%) in the {m t , m S } plane. The blue contour corresponds to the kinematic limit m t − m t − m S = 0. The maximum value for the coupling to be in the perturbative region has been limited to 4π.

B Range of validity of the narrow-width approximation
In the processes under consideration both t and S are assumed to be in the narrow-width approximation (NWA), in order to factorise the production of the top partner from its decay chain. Such assumption, however, implies that the coupling t tS cannot exceed specific values which depend on the masses of t and S according to the relation in eq. (2.2). Considering as a simplifying and extreme assumption that the only available decay channel for t is into the SM top and S and that one chirality of the couplings is dominant with respect to the other, such that either κ S R κ S L or vice versa, the values of the coupling corresponding to different Γ t /m t ratios is shown in figure 15. For a specific {m t , m S } configuration, values of the coupling larger than those in the contours of figure 15 would produce a larger width. The determination of the validity of the NWA approximation To assess how the width of t affects the determination of the cross-section, the full 2 → 4 process pp → ttSS has been evaluated by imposing the presence of at least one t propagator in the topologies, in order to obtain the signal under the assumption of negligible Stt coupling. With such process, the off-shellness effects and contribution of topologies such as those in figure 16 are fully taken into account. Still under the assumption that t can only decay to S t and therefore that the only way to increase the total width of t for a given {m t , m S } configuration is by increasing κ S L,R , the ratio between the cross-sections of the full process and of the pair-production process in the NWA is shown in figure 17.
The effect of a large width is already noticeable when the Γ t /m t ratio reaches 1%, when the interference between the resonant channels and all the other contributions is negative and of the order of few percents in a region where m S + m t is around 80% of m t . If the Γ t /m t ratio is below 1% the relative ratio between cross-sections is dominated by the statistical fluctuations of the simulation. For this reason, the numerical results in the following sections assume Γ t /m t to be of order 0.1%.  Figure 17: Relative ratio of the cross-sections for the full process pp → ttSS (σ 2→4 ) and for the pair production process pp → t t → (S t)(St) where the t production and decay are factorised in the NWA approximation (σ Pair ). The ratio is shown for different values of the Γ t /m t ratios (0.1%, 1% and 10%), and the couplings κ S L,R are not allowed to exceed the perturbative limit 4π.

C Object definition
In the following, more details for the definition and selection of objects at reconstructed level are presented, elaborating on the brief description in section 4.1, in order to facilitate reproducibility and as a guide for possible future searches at colliders.
For all objects, the default ATLAS Delphes card [47] is used, with minor modifications in a few cases, as explained below. Objects that partially fall in the calorimeter transition region 1.37 < |η| < 1.53 are excluded, if they are reconstructed in the calorimeter, where η is the pseudorapidity. Relative angular distances in the detector are typically expressed as ∆R in the η-φ plane where φ is the azimuthal angle around the beampipe. A particle's transverse momentum p T is the momentum component in the plane transverse to the beam axis.
Isolation and overlap removal are needed to distinguish the objects from each other in the detector simulation, 3 which is done in the Delphes card, unless otherwise specified. This is achieved by creating the containers for the objects in mind: jets, photons, electrons and muons. In Delphes all objects passing their respective efficiency cut are first reconstructed as the respective object and as a jet. The object will then be put into the jet container and the container corresponding to the reconstructed object. By passing an isolation criterion the object is removed from the jet container and only kept in the container corresponding to the correct reconstruction. The criterion is met when an isolation variable I is within a certain constraint. The variable is defined by summing the p T of all objects, not including the candidate, within a cone of ∆R around the candidate and dividing by the candidate p T . That is, where the sum runs over all the objects i around the candidate within the ∆R cone. The objects used in the analysis are defined below.
Photons, γ, are reconstructed by considering energy deposits in the electromagnetic calorimeter (ECAL) and no tracks in the inner detector. Objects successfully reconstructed as photons are required to have a p T > 30 GeV and |η| < 2.37. Photons in the transition region are not taken into account. Overlap removals are done in the modified Delphes card as described above, where the photon candidate is identified and put in the correct container by passing the photon efficiency cut corresponding to the ATLAS tight quality efficiency cuts [73]. Isolation of the photon is done after the simulation and it is considered isolated when the isolation variable I < 0.008, where I is defined as described above.
Leptons, , are in the following understood to mean electrons or muons only, and not τ -leptons. Electrons are reconstructed by looking at both energy deposit in the ECAL and having a track in the inner tracking system. For the following, simulation in Delphes reconstruction of the electron is done by combining the reconstruction efficiency of the two subsystems and parametrise it as a function of energy and pseudorapidity. Muons pass the calorimeters and are reconstructed by combining the information from the inner tracker and the muon spectrometer. In Delphes, the user specifies the efficiency of the muons such that a muon is only reconstructed with a certain probability [47]. Leptons are required to pass an isolation criterion for which I < 0.12 within the cone ∆R < 0.2 for electrons and ∆R < 0.3 for muons. Furthermore, leptons are required to have p T > 25 GeV and be in the region of |η| < 2.47, excluding the transition region in the case of electrons. Further overlap removals of leptons are done in Delphes where the lepton candidate is identified and put into the correct container by passing the given lepton efficiency. For electrons, the efficiencies correspond to the ATLAS tight quality efficiency cut [74]. For muons, the default Delphes values are used.
Z bosons, Z, are identified as two leptons with same flavour and opposite signs, whose invariant mass fall within the window |M + − − m Z | < 10 GeV where M + − is the invariant mass of the reconstructed leptons.
Jets, j, are reconstructed by using the FastJet [53] package together with Delphes. Here the anti-k t algorithm [54] with a R parameter of R = 0.4 is in use for jet reconstruction. Jets are required to pass p T > 25 GeV and |η| < 2.47, excluding the transition region.
B-jets, j b , are jets which originate from the hadronisation of a b-quark. In Delphes this means a jet which contains a truth b-quark. The efficiency and misidentification rate is parametrised in Delphes based on estimates from ATLAS [47,75].
Missing transverse energy, E miss T , is computed in Delphes by taking the negative scalar sum of the transverse component of the momenta of all calorimeter towers (i.e., energy deposits in the calorimeter), E miss T = − i p T (i) [47]. The scalar transverse energy, H T , is computed by taking the scalar sum of the p T of all reconstructed basic objects used in the analysis, in this case: jets, muons, electrons and photons. All these objects which enter the H T definition are required to pass the stated analysis p T and η cuts.

D Signal efficiencies
In this appendix we present the signal efficiencies for each channel and mass point considered in the analysis, except those already shown in figure 9. Figures 18, 19, 20 and 21 show, respectively, the efficiencies for the γ γ SR and for channels where at least one of the two S decays into γ γ, γ γ SR and at least one S decaying to Z γ, Z γ SR and at least one S decaying to γγ and Z γ SR and at least one S decaying to Z γ.        Figure 19: Efficiencies for the γ γ SR and for channels where at least one of the two S decays into Z γ.     Figure 21: Efficiencies for the Z γ SR and for channels where at least one of the two S decays into Z γ.