Probing the Top-Higgs Sector with Composite Higgs Models at Present and Future Hadron Colliders

We study the production of $t{\bar t}h$ and $t{\bar t}hh$ at hadron colliders, in the minimal Composite Higgs Models, based on the coset $SO(5)/SO(4)$. We explore the fermionic representations ${\bf 5}$ and ${\bf 14}$. A detailed phenomenological analysis is performed, covering the energy range of the LHC and its High Luminosity upgrade, as well as that of a future 100 TeV hadron collider. Both resonant and non-resonant production are considered, stressing the interplay and complementary interest of these channels with each other and double Higgs production. We provide sets of representative points with detailed experimental outcomes in terms of modification of the cross sections as well as resonance masses and branching ratios. For non-resonant production, we gauge the relative importance of Yukawa, Higgs trilinear, and contact $t\bar{t}hh$ vertices to these processes, and consider the prospect for distinguishing the fermion representations from each other and from the Standard Model. In the production of top partners, we find that the three-body decay channel $W^+ W^- t$ becomes significant in certain regions of parameter space having a degenerate spectrum, and is further enhanced with energy. This motivates both higher energy machines as well as the need to go beyond the current analysis performed for the searches for these resonances.


Introduction
The discovery of the Higgs boson [1,2] , has aroused the interest in measuring with high precision the Higgs properties. This has led to greater emphasis on a strong joint effort between theory and experiments to explore the Higgs sector. The Large Hadron collider, LHC, at CERN, which is also a Top factory, offers an unique playground to perform the searches in the Top-Higgs sector both now and with its upgrade into the High-Luminosity LHC.
Besides, the importance of knowing in depth the Higgs sector, led to studying a Higgs factory as the next machine in Particle Physics as well as the importance of a very high energy hadron collider. While all current measurements within their present precision are consistent with a Standard Model (SM) Higgs boson, there is room for deviations from the SM that could hold the key to a deeper understanding of the phenomenon of electroweak symmetry breaking (EWSB). Furthermore, as the luminosity increases and experiments are progressively upgraded, a number of Higgs-related processes are becoming more and more accessible, in the next years to come. An important open question is whether the observed 125 GeV scalar is a composite bound state of more fundamental constituents, or whether it is elementary down to distances much shorter than currently explored. The experimental program will be essential in illuminating this fundamental question.
The idea of Higgs compositeness received an important boost with the construction of rather complete models that can be consistent with all current measurements [3]. Such constructions incorporated a geometric solution to the hierarchy problem [4], a dynamical mechanism for EWSB, and an appealing understanding of the flavor structure observed in the SM. In these scenarios, the Higgs boson is understood as a pseudo-Nambu Goldstone boson (pNGB), somewhat analogous to the pions in QCD, but also with important differences. While much model building and extensions have been proposed in the literature, in this work we focus on the minimal setup, based on the symmetry breaking pattern SO(5) → SO (4). These are referred to as "Minimal Composite Higgs Models" (MCHM).
It is known that there is considerable model-dependence associated with the fermion sector of the MCHM, which is of relevance to our work. In particular, the top sector is expected to play a crucial role in these models, given that the top quark couples most strongly to the Higgs boson. We focus on the production of one or two Higgs bosons in association with a top/anti-top pair. The tth cross section is being actively measured [5,6], in many different channels, with values compatible with the SM prediction, within approximately 20% uncertainty. An experimental search for the tthh process has been just performed, for the first time at the LHC [7]. Such a process is of particular interest in the present class of models, due to the generic prediction of charge 2/3 vectorlike "top partners" that can decay in the th channel, thus leading to the previous final state. These resonances are already constrained by previous searches in this channel [8] as well as in combination with tZ and bW decay channels [9,10].
Since the top partner resonances are under active search with already defined bounds above 1 TeV, we bring here attention to the fact that frequently a large fraction of the tthh cross section is composed of non-resonant production, which then becomes an important discriminator for the MCHM models. This component is also more directly related to the pNGB nature of the Higgs boson, unlike the vectorlike fermionic resonances. This work highlights the complementary role between non-resonant tthh production and the tth and hh channels, specially in terms of measuring the trilinear Higgs couplings. We gauge the relative importance of Yukawa, trilinear Higgs and new contributions to the non-resonant tthh channel, accessing quantitatively the correlation between tth and tthh cross sections.
To illustrate the interest and importance of these two processes for the Higgs sector and within the Composite Higgs context, this work presents a phenomenological analysis at the parton level, with two specific realizations of the MCHM, in the framework of the present and future hadron colliders.
The goal is to establish phenomenological differences between a model with minimal embedding of fermions (the MCHM 5 ) and the simplest one that can allow for an increased top Yukawa coupling (the MCHM 14 ). We also provide sets of representative points with a large coverage of the parameter space of those models. Beyond these two specific MCHM realizations used as a showcase to scan and thus explore experimentally the MCHM, we link these two cases to the effective field theory (EFT)appli cable in the limit of decoupled resonances. To do so and complete this study, we also present our results in terms of modifications to the SM couplings which are more directly comparable to experiment.
Therefore this overall study can be used to guide searches for new physics, bridging the BSM phenomenological and experimental goals. Indeed, starting from the current results being achieved at the LHC (and promptly evolving), this study provides some defined directions and predictions to look for, over the next decades, at the 14 TeV High Luminosity phase of the LHC and at a future High energy Hadron collider in the range of 100 TeV to possibly 150 TeV, as in the proposed FCC-hh at CERN and SppC in China [11][12][13].
The features of the composite Higgs scenarios directly relevant for this work are reviewed in Section 2. Section 3 sets the overall phenomenological analysis strategy. It includes the definition of the parameter space for the two MCHM scales considered in this work, the simple physical observables relevant here, the strategy and tools to extract the representative points in the parameter space, the implementation of both models into the event generator and the operator analysis with its effects on the tth and tthh processes. The MCHM low scale scenario is described in Section 4 with the corresponding sets of representatives points in both models, their detailed experimental outcomes in terms of modification of the cross-sections, of the non-resonant components and new resonances. Similar outcomes are shown in Section 5 for the MCHM at high scale that will be only reachable with the future high energy hadron colliders. Section 6 describes the EFT perspectives, stressing for instance the correlation between some selected EFT parameters in the MCHM 5 and the MCHM 14 both at low and high scales. Section 7 concludes by stressing how this detailed phenomenological analysis of these two MCHM scenarios, serves as an useful showcase for the exploration of the BSM world by experiments, with well-defined but still broad scope guidelines. It highlights in this way, what can be already achieved at the presently running LHC and the forthcoming HL-LHC era. It emphasizes the unique importance of a future high energy hadron collider to explore in details the Top-Higgs sector.

The Higgs as a Pseudo-Nambu Goldstone Boson
We will present only the features of composite Higgs scenarios that are directly relevant to our study. In particular, we will not describe the spin-1 sector, which is fixed by the pattern of symmetry breaking, in our case SO(5) → SO (4). We refer the reader to the complete review [14] where the general construction and results for the MCHM bosonic sector are presented. For our purposes, it is sufficient to know that the breaking of the electroweak (EW) symmetry can be parametrized in unitary gauge through the orthogonal matrix where h is the Higgs boson, h 0 = h is its expectation value, and f is the scale at which the breaking SO(5) → SO(4) occurs. In Eq. (2.1), 1 3×3 is the 3 × 3 identity matrix and 0 is a 3-dimensional null vector. We also recall that in this model one finds We will often use the variable to express our results. It characterizes the deviations from a SM Higgs due to compositeness. The current bound consistent with Higgs data is f 800 GeV, or ξ 0.1, [15][16][17][18][19][20][21]. Our main focus is on the fermion sector and its interplay with the Higgs boson. This depends on the composite resonances associated with the top quark. Up to EWSB effects, these resonances fall into representations of the unbroken SO (4). One expects to find that sets of these SO (4) representations build up representations of SO (5) with the non-degeneracy arising from the spontaneous breaking of SO (5) (characterized by the scale f ). For concreteness, we will focus on two possibilities: • Resonances falling into a 5 of SO (5), which split into SO(4) multiplets as 5 = 4+1.
The first case allows for the minimal number of extra fermionic degrees of freedom, while imposing the custodial protection of the Zb L b L coupling [22], which is important when considering EW precision measurements [23][24][25][26][27][28][29]. The second possibility is non-minimal and has been considered in [16,19,[30][31][32][33][34][35]. It has been pointed out that it allows for an enhancement of certain Higgs couplings w.r.t. to the SM [19,32], in contrast to the minimal case (and other possibilities) which always leads to suppressions. One of our goals in this work is to contrast these two scenarios from the point of view of searches for the tth and tthh processes at the 14 TeV High Luminosity phase of the LHC, and the High Energy pp collider projected to run at 100 TeV, and even up to 150 TeV. We describe the relevant features of the previous fermion embeddings in the following sections.

The Fermion Sector of the MCHM 5
Considering only the top and its partners, the model is comprised by the "elementary" fields q L = (t L , b L ) and t R together with a set of "composite" fermionic resonances. The elementary fields have the SU (3) C × SU (2) L × U (1) Y quantum numbers of the SM left-handed (LH) top-bottom SU (2) L doublet and of the SM right-handed (RH) top SU (2) L singlet, respectively. The composite resonances fall into SO (5) representations that split into SO (4) representations due to the spontaneous breaking at the scale f . They contain a number of vector-like SU (2) L × U (1) Y representations, that depend on the fermion embedding. The smallest representation compatible with custodial symmetry is the 5 of SO (5), whose decomposition under SO(4) is given by a fourplet, Ψ 4 , and a singlet, Ψ 1 . We denote the corresponding states by Ψ 4 ∼ (X 5/3 , X 2/3 , T, B) ,

4)
The subindex denotes the electric charge. While not explicitly indicated, the T andT states have charge 2/3 and B has charge −1/3. The states (T, B) transform as a SU (2) L doublet and have hypercharge Y = 1/6, i.e. they have the same SM quantum numbers as the elementary field q L , while an exotic SU (2) L doublet with Y = 7/6 is composed of the (X 5/3 , X 2/3 ) states. Finally,T has the SM quantum numbers of t R . The elementary sector is simply described by where D stands for the SU (3) C × SU (2) L × U (1) Y covariant derivative. We write the composite sector directly in terms of the SO(4) multiplets: Here, the covariant derivative contains only the gluon and hypercharge fields, that is, The remaining electroweak interactions are inside the d µ and e µ symbols, which are defined in terms of the Maurer-Cartan form iU −1 (∂ µ + ig a A a µ T a )U = d µ,â Tâ + e µ,a T a ≡ d µ + e µ , (2.7) where T a are the generators of the unbroken SO (4) and Tâ are the broken generators. The gauge fields A a µ belong to the algebra of SO(4) ∼ = SU (2) L × SU (2) R , in which the electroweak fields are embedded as follows: the W fields gauge SU (2) L while B µ gauges T 3 R and the remaining generators are ungauged. The hypercharge is then given by Y = 2/3 + T 3 R . 1 This covariant derivative allows for the non-linear realization of the full SO (5) symmetry in the kinetic terms, even though the Lagrangian (2.6) exhibits explicitly only the SO(4) symmetry [36,37]. However, note that the SM covariant derivatives break the SO(5) global symmetry explicitly. In App. A we give the SO(5) generators, including those of the gauged SU (2) L ×U (1) Y subgroup. The e µ symbol term contains corrections to the electroweak interactions of the resonances, due to compositeness. These are detailed in App. C. Apart from the "kinetic" terms, we include separate mass terms for Ψ 4 and Ψ 1 . The difference M 4 − M 1 arises from the spontaneous breaking of SO (5), which we do not describe here. For our purposes, M 4 and M 1 can be treated as independent phenomenological parameters.
As mentioned above, some of the fermionic resonances have the same SM quantum numbers as the elementary fields, leading to the possibility of mixing between them. In order to write the elementary-composite mixing terms, it is convenient to embed the elementary states using SO (5) notation as follows while Ψ 4 and Ψ 1 are similarly written in 5-plet notation as (see App. B for further details) (2.9) In terms of these definitions, the mass mixing Lagrangian takes the form thus implementing the idea of partial compositeness [38]. Finally, we include in our Lagrangian additional Higgs interactions involving the d µ symbol, which are allowed by the symmetries at the lowest order of derivatives, and are expected to arise from integrating out heavy resonances not included in our low energy theory (see [14] for further details). These are given by (2.11) where P L,R = (1 ∓ γ 5 )/2 are the left and right projectors and c L and c R are couplings, expected to be order one. 2 These terms contain extra Higgs interactions 3 , detailed in App. C. We should emphasize, however, that these additional couplings do not affect the top Yukawa at tree level, as long as c L R are taken to be real. This can be seen by noting that the operator in Eq. (2.11) is antisymmetric in Ψ 4 , Ψ 1 , such that terms with the same mass eigenstate fermion will cancel out between the operator and its complex conjugate. This holds in fact for similar operators built out of any fermion representation, see [32], which will be relevant for the 14 representation described below. For this reason, these modifications will only be important in characterizing the extended fermionic sector of the models. We find that the shape of distributions is largely unaffected by these terms, baring the possibility of a tuned cancellation between the vertices of Eq. (2.11) and those of Eq. (2.10), as we will see in section 3.
The complete Lagrangian (in the top sector) is The charge 2/3 mass matrix in the {t L , T L , X 2/3,L ,T L } vs {t R , T R , X 2/3,R ,T R } basis is then given by Diagonalization of this matrix leads to the physical fermion eigenstates, which are in general admixtures of the original elementary and composite states. The lightest one is identified with the observed top quark. Our numerical analysis follows from this mass matrix, as described in subsequent sections. The remaining resonances have masses 14)

The Fermion Sector of the MCHM 14
In the second scenario, instead of assuming that the composite states span a 5 of SO (5), we assume that they span a 14 of SO (5). This multiplet decomposes under SO(4) as a fourplet and a singlet, as in Eq. (2.4), plus a nonet Ψ 9 . We denote the corresponding states by Under SU (2) L , this nonet breaks into three triplets, U with Y = 5/3, V with Y = 2/3 and F with Y = −1/3. Adding to Eq. (2.6) the nonet Ψ 9 , we obtain (the precise structure of Ψ 9 is given in App. B): To write the mixing between the elementary and composite sectors, it is convenient to formally embed all the elementary and composite states into "14" representations of SO (5), in analogy to what was done for the 5 case. We denote the elementary embeddings by Q 14 L and T 14 R and continue using the notation Ψ 9 , Ψ 4 and Ψ 1 for the composite embeddings. All of these become 5 × 5 traceless symmetric matrices, whose precise form is given in App. B. The mixing Lagrangian is then written as We also include extra d µ symbol interactions allowed by the symmetries (2.19) where c 4 , c 9 and c T 9 are order one couplings and i, j are SO(4) indices. Here, for simplicity we take the strong sector to be parity symmetric. We also expect the two derivatives term with c T 9 to be subdominant in most channels, since it is suppressed by an extra power of the cutoff Λ 4πf . The explicit form of these vertices as well as those arising from the e µ symbol in the kinetic term are reported in App. C. The complete Lagrangian (in the top sector) is which leads to the charge 2/3 mass matrix in the {t L , T L , X 2/3,L ,T L , U 2/3,L , V 2/3,L , F 2/3,L } vs {t R , T R , X 2/3,R ,T R , U 2/3,R , V 2/3,R , F 2/3,R } basis: where we defined The remaining states have masses As in the 5-plet case, the previous mass matrices form the fundamental input to our phenomenological analysis.

Partial Compositeness and Higgs Couplings
The above models incorporate the partial compositeness paradigm of [38], via linear mixing of the elementary fields q L and t R with composite operators transforming as singlets, 4-plets or nonets of the SO(4) symmetry as described by Eqs. (2.10) and (2.18). In addition to giving rise to the top mass, the same operators are responsible for the top-Higgs Yukawa coupling, which is of central importance to this work. 4 The mechanism is illustrated diagrammatically in Fig. 1, where H is the Higgs doublet [see Eq. (B.2)]. Each green box represents an insertion of the corresponding operator in Eqs. (2.10) or (2.18), to leading order in H/f . For example, the mixing of the singlet Ψ 1 with T R can happen at 0-th order in H, while the Ψ 1 -Q L mixing requires an insertion of the Higgs field, which transforms as a 4 of SO(4): 4 Q L ⊗ 4 H ⊃ 1. Similarly, the mixing of the 4-plet Ψ 4 with Q L can happen at 0-th order in H, but requires an The mixing with a nonet resonance is qualitatively different, requiring one H insertion for the Q L -Ψ 9 mixing and two insertions for the T R -Ψ 9 mixing: 4 Q L ⊗ 4 H ⊗ 9 Ψ 9 ⊃ 1 and 4 H ⊗ 4 H ⊗ 9 Ψ 9 ⊃ 1, respectively. As a result, the leading order coupling thus induced is the non-SM like, non-renormalizable operator The top mass is obtained by replacing H by its vev, leading to a factor of three in the ratio of the top mass to the top Yukawa (associated to the tth operator) when it is induced by the "cubic" interaction than for the "linear" interaction. The presence of the various channels simultaneously can then lead to an enhancement of the top Yukawa coupling w.r.t. the SM. One should also notice that when M 1 = M 4 , the linear coupling cancels out, 5 and the leading order is cubic [19].
Although we will not use it in our numerical analysis, useful approximate expressions for the top mass in the MCHM 5 and in the MCHM 14 are given by where we defined the mass ratios r 1 = M 1 /M 4 and r 9 = M 9 /M 4 and took y Li ≡ y L and y Ri ≡ y R for i = 1, 4, 9 (see Section 3.1). The "wavefunction renormalization" factors have the form Z L,R = 1 + y 2 L,R f 2 /M 2 4,1 + O(ξ). The expression (2.25) displays explicitly the behavior described above.
The same effect has an impact on the coupling of the Higgs to two gluons (ggh), normalized to the SM top contribution, which is given by: where M is the fermion mass matrix, assuming all states are much heavier than the Higgs boson (for our purposes, light state contributions can be neglected).
For the MCHM 5 , using M = M 5 2/3 in Eq. (2.13), this gives For the MCHM 14 , one gets , (2.29) arises from the charge 2/3 sector, M 14 2/3 in Eq. (2.21), and arises from the charge −1/3 sector, M 14 −1/3 in Eq. (2.22). For the latter, we explicitly removed the zero-mode (the physical bottom quark) and neglected its contribution to the ggh coupling. When in all other cases, reflecting the underlying cubic versus linear coupling of the top quark to the Higgs field. In all cases, c b g ∼ ξ 1. Global constraints on the gluon fusion process from Higgs measurements, which allow for about 20% deviations from unity in c g at the 95% C.L. [40], will then also impose constraints on the allowed deviations in the top Yukawa coupling from the SM limit.

Higgs Decays
While the amplitudes of the processes tth and tthh are determined by the Lagrangians written in the previous sections, it is necessary to specify how the light families of the SM are treated in the context of the composite Higgs models in order to take into account the possible modifications in Higgs decays. Such deviations are expected to be small, since the observed 125 GeV resonance is known to exhibit SM-like properties. The dominant Higgs decay channels are then as in the SM: h → bb, W + W − , gg, τ + τ − , cc, ZZ. We neglect the decays into γγ, γZ and µ + µ − , which have branching fractions ranging from 0.23%, 0.15% down to 0.01%, as well as even rarer decay channels.
Given the importance of the bb channel and the fact that b L is embedded into SO (5) together with t L , we must also specify how b R fits into the models. There are several ways to proceed. Rather than exploring the various possibilities, we choose to supplement the RH bottom with composite resonances that span a 10 of SO (5) for both the MCHM 5 and the MCHM 14 . Such models, called MCHM 5,5,10 and MCHM 14,14,10 , were introduced in [16]. For the remaining fermions, we choose to replicate the scheme employed for the third family. We also assume that the lepton sector follows the same scheme as the quark sector. Furthermore, we assume that the mixing angles between the elementary and composite states associated with the light families are small, as in "anarchy" models of flavor (see e.g. [41,42]).
Under the previous assumptions, it is found that the couplings of the composite Higgs to the vector bosons and to light f f pairs are controlled by two model-dependent functions that depend only on ξ. In the MCHM 5 the coupling of the Higgs to a pair of gluons depends also only on ξ, but in the MCHM 14 it depends on additional microscopic parameters, as shown in Subsection 2.3.
The partial widths are then simply obtained by rescaling the SM ones. For the MCHM 5 , one finds [16] For the MCHM 14 , the bottom channel is controlled by F 1 instead of F 2 , and the ggh coupling is controlled by c 14 g of Eq. (2.28) instead of F 1 (ξ). The total Higgs width in the MCHM models under consideration can then be written as

Parameter Space
The top sector of the models described in Sections 2.1 and 2.2 is controlled by the (vector-like) mass parameters, M i , and by several dimensionless couplings, y Li , y Ri , as well as c L and c R for the 5 and c 4 , c 9 and c T 9 for the 14. These parameters are a priori complex. However, not all phases are physical. In order to identify the number of physical phases we can proceed as follows. We can start by absorbing the phase of each M i (thus making it real and positive) by redefining the phases of Ψ i,L or Ψ i,R . This leaves one free phase in each such pair, say in Ψ i,R , that can be adjusted to absorb the phase of the corresponding y Li (thus making all of the y Li real and positive). Finally, we can absorb the phase of one of the y Ri into T 5 R or T 14 R . We conclude that there are three (five) physical phase(s) in the MCHM 5 (MCHM 14 ). Alternatively, we can choose all the y Li and y Ri to be real and positive, putting the physical phases in M 1 , c L,R for the MCHM 5 , and in M 1 , M 4 , c 4 , c 9 , c T 9 for the MCHM 14 . In this work, for simplicity, we will assume that all parameters are real, which amounts to imposing CP conservation in the strong sector 6 . This leaves three physical signs in the case of the MCHM 5 and five signs in the MCHM 14 . We will choose these signs to be sign(M 1 ), sign(c L ) and sign(c R ) in the first case and sign(M 1 ), sign(M 4 ), sign(c 4 ), sign(c 9 ) and sign(c T 9 ) in the second. Finally, in order to simplify our analysis, we will disregard the derivative couplings in eqs. 2.11 and 2.19 until the end of Sec. 3, which leaves us with M i , y Li and y Ri . The effect of the neglected operators will be considered separately in Sec. 3.8.
While the above parameters respect the SO(4) symmetry, in general they violate the SO(5) symmetry. The SO(5) symmetric limit corresponds to M 1 = M 4 and y L1 = y L4 , y R1 = y R4 for the MCHM 5 , and M 1 = M 4 = M 9 and y L1 = y L4 = y L9 , y R1 = y R4 = y R9 for the MCHM 14 . It turns out that deviations from the SO(5) symmetric limit in the dimensionless couplings corresponds to a "hard" breaking of the symmetry, in the sense that the Higgs effective potential is finite when the SO(5) symmetry relations are satisfied, but becomes UV sensitive when not 7 . Deviations from the SO(5) symmetric limit in the M i , on the contrary, correspond to a "soft" breaking, and the Higgs potential remains IR dominated in that case. This motivates us to focus on the case where for all i, as a way to reduce the number of independent parameters. Small deviations from this limit mean that one can expect additional UV dependent contributions to the Higgs potential. Such contributions can affect the region of parameter space that leads to EWSB and to a Higgs mass m h = 125 GeV. Thus, we take the point of view that imposing that the Higgs mass be reproduced by strictly adhering to the case of a calculable Higgs potential is overly restrictive in the context of our collider study. For this reason we will simply fix the mass of the Higgs, and leave the study of points that reproduce the Higgs mass in the strictly calculable limit to future work.
In conclusion, we are left with the following set of parameters: • MCHM 5 : f , |M 1 |, |M 4 |, sign(M 1 ), y L and y R .
One of these parameters can be further fixed by requiring that the top mass be reproduced. We choose to fix y R in this way. Our procedure is to require that 2), we take m t = 150 GeV, the running top mass at the scale of the resonances, which will be typically around 2 − 3 TeV 8 . On the other hand, in the tth and tthh production the relevant scales are of the order of a couple hundred GeV. We therefore distinguish between the high-scale running top mass (relevant for the diagonalization of the mass matrix), and a low scale running top mass, relevant to the physical processes of interest. We take for the latter the pole top mass of m t = 173 GeV, which also enters in kinematical quantities. To first approximation, this takes into account the running between the two scales.
Our strategy to extract the physical quantities is straightforward: given values for the M i and for y L , we find y R from Eq. (3.2). We then diagonalize numerically the fermion mass matrix to obtain the spectrum, and the unitary transformations U L and U R such that with all the physical masses real and positive. This is done with Mathematica [44]. We also treat the charge −1/3 sector in the MCHM 14 numerically. 9 The physical spectrum and the rotation matrices are the main input to the rest of the numerical analysis.
For example, we can obtain other quantities of interest, such as the Yukawa matrix in the mass eigenbasis is the Yukawa matrix in the gauge eigenbasis. The most relevant quantity will be the (1, 1) entry in Y mass 2/3 , which corresponds to the top Yukawa coupling. It includes exactly all treelevel effects arising from the Higgs compositeness (the dependence through √ ξ = sin h 0 f ), as well as the mixing with the vector-like resonances.
The non-linear dependence on h in composite Higgs models leads to interactions between top pairs and a number of Higgs bosons. The tthh vertex, whose Feynman rule is given by i times the (1, 1) entry of 1 2 d 2 M 2/3 /dh 2 0 , after rotating to the mass eigenbasis, also enters in our analysis.

MCHM Scales, low versus high
The parameter space of MCHM 5 and MCHM 14 is explored here in two steps. The first step considers the parameter space relevant for the reach of the LHC machine. This step includes two running operation stages of the LHC, i.e from now until 2023-2024 with about 400 f b −1 total integrated luminosity at 13 and 14 TeV, and from 2026 to about 2038, with 10 times more luminosity and 14 TeV CM energy (may be slightly more) with the HL-LHC. The region of the parameter space corresponding to the overall LHC machine operation (from now until the end of the HL-LHC) is labelled as the "Low Scale MCHM", as the dimensionful parameters will take values of a few TeV. This will be the focus for the remainder of this section and section 4.
The second step of this analysis is extended to the "High Scale MCHM". This relates to the future hadron colliders in project, expected to run at CM energies around 100 TeV or even higher [11][12][13]. For the MCHM high scale regime, the starting hypotheses are either: 1. No new physics is discovered at the HL-LHC, i.e. within the possible reach in mass and/or precision of this collider, or 2. Some evidence (3σ effect) is found such as new high mass resonance(s) or a deviation from the SM for µ(tth), i.e., the top Yukawa, or 3. A deviation on the µ(tth) is present at 5σ and µ(tthh) is observed, but with sizeable uncertainty. Thus a higher energy pp collider would allow higher precision measurements and looking for further effects.
This is the subject of section 5. But it is worth stressing already here, and when looking to the results in section 4, that the points generated in the MCHM low scale parameter space, are generated both at 14 and at 100 TeV. This is simply because of the reasons listed above.
Here below are specified the ranges considered in each case and the reasoning behind them.

Definition of the ranges for the Low Scale parameters
For the MCHM 5 , we consider the following ranges for the parameters: For the MCHM 14 , we use: In order to remain in a perturbative regime, justifying the present tree-level analysis, we will take y L < 3. For the same reason, we also check that y R , as determined by the top mass, is below 4. The distribution of points within those ranges was not uniform, due to computing constraints, but we strive to cover most of the parameter space.

Definition of the ranges for the High Scale parameters
The range to be covered by the parameters for the High scale case is defined, for each considered MCHM scenario, by the possible reach of a high energy hadron collider order 100 TeV in CM and at least 20 ab −1 total luminosity, taking also into account previous studies [45]. Moreover this parameter space range must be linked continuously to the one defined for the Low Scale, which will already be mainly tackled by the HL-LHC; indeed, some showcase scenarios in the Low scale will remain of interest in the High scale as pointed out in the next two sections.
Given the above, for the MCHM 5 we consider:

Simple Physical Observables
The diagonalization of the Q = 2/3 mass matrix leads to several physical quantities of interest. There is a rich spectrum of vectorlike states. Up to small EW symmetry breaking effects, these vector-like masses are approximately given by: • MCHM 14 : We use the mass of the lightest state, which we call T (1) , as a proxy for the scale of the new physics. There are important lower bounds from direct searches for such resonances by both the ATLAS [9] and CMS [10] collaborations. These searches consider vectorlike top partner resonances decaying exclusively in the bW , tZ and th channels. These bounds depend mildly on the decay branching fractions of the heavy state and are roughly around 1.3 TeV. However in the models under consideration here, one finds sometimes additional decay channels, as we will show in section 4, so these bounds must be taken with care. We also consider the deviations from the SM of the top quark decay width. However, the previous direct bounds imply that such deviations are well within the experimental uncertainties, and therefore do not impose additional bounds on the models.
As mentioned earlier, a quantity of direct interest in our study is the top Yukawa coupling. Figure [19]. Specifically, for the smaller range, only Region III allows an enhancement, while, for the larger scan range, the top Yukawa can also be enhanced in Region I. We also show red contours of constant mass of the lightest top partner resonance M T (1) . In the dark overlayed regions, the lightest top partner is approximately excluded by direct searches, assuming it decays only in the bW , tZ and th channels [9,10], while the green overlay shows the region in which µ(ggh) differs from unity by 20% or more, which is in tension with current Higgs coupling measurements [40]. In the white area, the top mass cannot be reproduced by values of y R within our considered perturbative range.

Strategy to select example points and benchmark points
We select two classes of points to study the physics, in each of the two MCHM scenarios, for both the MCHM at low and at high scale. These are respectively labeled the "example points" and the "benchmark points".
The selected example points are chosen based on striking features or on how accessible they are in the near future (at the start of the HL-LHC) or towards the end of the HL-LHC or in the much longer term with a 100 TeV pp collider. The striking features are defined with the present results from the LHC experiments on the Top-Higgs sector or the prospect studies achieved for the HL-HE/LHC scenarios or the FCC-hh project. However they do not necessarily represent what is the typical behaviour of the parameter space we explored. They are picked according to our criteria of what is interesting and thus carry a bias. While that is very useful to see what kind of phenomenology can be produced by the model, and give insight on what is happening or could happen in the future, it is not the most comprehensive and extensive way of looking into the possibilities of the models.
On the other side, looking individually on all kinematic distributions for hundreds of points in parameter space is just unfeasible. An interesting compromise can be reached using the approach proposed in [46]: one can use a statistical test to group points of the parameter space into "clusters" based on the similarity between the kinematic distributions of the final or intermediary states produced by them. This is called the "clustering strategy".
Following this strategy, one can then use the same test statistic to choose one point from each cluster as a typical representative of that behaviour, a benchmark point. Analysis designed to search for those benchmark points will be guaranteed to cover all possible phenomena in the region of parameter space considered. We outline the main steps of this algorithm below and refer the reader to [46] for further details and discussion.
The first task is, given two different points in parameter space and one or more kinematic distributions generated by these points (the p T of the Higgs or top quark, invariant masses, etc.), to decide how similar the distributions are. We organize each of these sets of distributions in samples S a , where a identifies the point in parameter space, so that running over the bins of the sample is the same as running over all bins of all distributions included in the analysis. In order to decide on similarity between the samples we will use the following log-likelihood ratio: where n (i,a) is the number of event counts in the i-th bin of sample S a , n (i,b) is the same for sample S b and N bins is the number of bins in the sample. T S ab is zero for identical samples and increasingly more negative for increasingly different distributions, so if T S ab > T S cd it means S a and S b are more mutually similar than S c and S d are. Now, starting from a number of samples N sample , we follow the steps below to organize them into clusters: 1. We start with a number of clusters that is equal to the number of samples, thus with N cluster = N sample , each cluster containing exactly one sample.
2. We obtain the cluster-to-cluster similarity between two clusters, defined as the minimum T S ab between members of those clusters: T S min = min ab ({T S ab }) (a runs over all samples in the first cluster and b does the same for the second cluster).
3. We calculate T S min between all possible pair of clusters, and find the two clusters with the highest T S min . Merge these two clusters into one. The number of clusters N cluster diminishes by one.
4. We repeat step 3 until the desired N cluster is obtained.
For each step in the clustering algorithm we can also choose one special sample within each cluster that is the best representative of its behaviour. We do that by calculating T S min a = min b ({T S ab }) for each sample S a in the cluster, where b runs over every element in the cluster except a. The sample with the highest T S min a is the benchmark sample and the equivalent point in parameter space will be called a benchmark point.
The appropriate final number of clusters is a compromise between a very fine grained view, with a unwieldy number of very homogeneous clusters (in the limit we go back to N cluster = N sample ), and the opposite extreme, with just a few clusters that contain a huge number of samples that are very heterogeneous in behaviour (thus with a benchmark sample that will not be a very good representative of the whole group). What can be done is to run this algorithm all the way down to N cluster = 2 clusters, keeping record of each step. That way one can examine each of the different realizations and decide on the ideal number. The same can be said about which kinematic distributions to include in the samples used for the clustering: we find the most interesting ones by experimentation.

Implementation of both MCHM models into the Event Generator
We have implemented both models in FeynRules (v2.3) [47] and produced an associated UFO file for each model, that can be interfaced with MadGraph 5 (v2.6.2) [48]. The numerical input from the diagonalization of the mass matrices is then fed via a customwritten Python script into the param card.dat for processing within MG5. We simulate the tth and tthh processes in MG5. We have also checked that the deviations from the SM in the top quark properties are negligible, since the new physics is rather heavy.
It is important to stress here that the generation framework we fully developed at the parton and LO level can be connected directly to detector simulations such as DELPHES (fast simulation) or detailed experiment full simulations such as the ones of ATLAS or CMS.

The tth Process
We start by considering the tth process in the MCHM scenarios. This is related in a very simple way to the same process in the SM. The Feynman diagrams (at tree level) are identical in all the models, involving top/anti-top pair production, with a Higgs boson radiated from the top lines. Radiation of the Higgs boson from initial state qq lines can be neglected due to the small Yukawa couplings (and PDF suppressions for the heavier flavors). As a result the amplitude is simply proportional to the top Yukawa coupling. The cross section in the MCHM scenarios can then be simply expressed in terms of the SM cross section as All the modifications due to Higgs compositeness, or mixing with vector-like fermions, enter only through the top Yukawa coupling. Therefore, as in other BSM cases, a modification in the total rate w.r.t. the SM is expected, but not in kinematic distributions.
Besides as we will see in some example cases this deviation can be very small, and still compatible with MCHM. This means that a high precision (order 1% or less) measurement of µ(tth) might be required (see sections 4 and 5)

The tthh Process
Next, we consider the tthh process. The additional radiated Higgs boson allows for a richer dependence on the new physics than in tth. There are two qualitatively different contributions: 1. Resonant processes, in which vectorlike charge 2/3 resonances are produced, and subsequently decay in the th channel. The resonances can appear either in pairs (QCD pair production) or singly (as an intermediate state in a fermion line, involving a flavor-changing Yukawa interaction). The process is, however, largely dominated by QCD pair production.
2. Non-resonant processes, in which only the diagrams without intermediate top partners are included, see Fig. 3. Figure 3. Representative diagrams for the non-resonant tthh process, illustrating the three distinct physical subprocesses: the Yukawa vertex, the Higgs trilinear self-coupling and the "double Higgs" Yukawa vertex arising in composite Higgs scenarios.
The presence of resonant processes can lead to important enhancements in the tthh cross section w.r.t. the SM, depending on their mass. The non-resonant process carries information that is distinct from the resonant part, as discussed below. It is therefore useful to define a "non-resonant cross section" as obtained from this subset of diagrams, which we label as "NR-tthh". One can similarly define a resonant cross section in terms of the diagrams involving QCD vector-like pair production. We find that, to an excellent approximation, the total tthh cross section is given by the sum of these two cross sections. The dependence of these two subprocesses with kinematic variables (such as the th invariant mass) is different, so in principle they can be separated experimentally.
We show in Fig. 4 the invariant mass distribution of th for a point in our scan of the MCHM 5 , where we display both the resonant and non-resonant contributions to the tthh cross section, as well as the corresponding SM process for comparison. We notice that the NR-tthh follows the SM cross section but displays a suppression and the relative importance of the resonant process w.r.t. the non-resonant one increases with larger CM energies. The cross section for both processes also increases significantly with the CM energy increase from 14 to 100 TeV. Likewise (µ(tthh) increases by a factor of 4 for the same increase in CM energy, while µ(NR-tthh) does not change with energy.

The Non-Resonant tthh Process
The p T distributions for the NR-tthh process show that the typical scales involved are around 100 GeV. As explained before, the top Yukawa coupling should be evaluated at around that scale for the Yukawa vertices appearing in the NR-tthh cross section. For simplicity, we use the top quark pole mass, which also appears in kinematic quantities.
Depending on the vertices they contain, we may divide the diagrams in the MCHM scenarios into three categories, as illustrated in Fig. 3: 1. Diagrams involving only the tth vertex.

Diagrams involving both the tth vertex and the trilinear Higgs self-interaction:
3. Diagrams involving the tthh vertex ("double Higgs" Yukawa vertex).
While the first two categories are composed of the same diagrams which appear in the corresponding SM process, the third category is particularly interesting, since the contact tthh vertex has no counterpart in the SM, and is a direct consequence of the non-linear realization of the Higgs sector in composite models [50]. For this reason, it would be extremely interesting if one could get experimental evidence on this coupling. In order to get a sense for the relative importance of the different physical subprocesses, we have simulated the NR-tthh cross section turning off, in turn, the double Yukawa coupling and the trilinear coupling. Results of this study are shown in subsection 4.2.3.

Operators Analysis
We move on to studying the effects of the high dimension operators in eqs. 2.11 (controlled by the couplings c L and c R ) and 2.19 (controlled by c 4 , c 9 and c T 9 ) on the two processes under consideration.

Effect on the tth Process
As mentioned in section 2.1 the derivative operators of eqs. 2.11 and 2.19 do not modify the tth vertex at tree level, appearing only in vertexes involving the Higgs and one or more of the new fermionic resonances. Since the leading diagrams in this process do not involve these resonances the changes in their decay widths do not affect it. Therefore, the cross section of the tth process is not modified by the new operators at tree level.

Effect on the tthh Process
The tthh process cross section is affected in two ways. First, the leading diagrams in this process contain resonances as intermediate states. The change in their decay widths will therefore have an effect on the total cross section. The other source of modification arises from the Yukawa type vertex between the Higgs and two different flavour states which enters directly into the main diagrams.
In order to have a better comprehension of the implications of the operators in the tthh total cross section, we took some example points in the parameter space and did a scan in the region of the (c L , c R ) plane compatible with perturbativity (c L , c R ∈ [−3, 3]). In Fig. 5 we show the scan for one of the points. There, we see that there is a modification of the tthh cross section by a factor that lies in the range [0. 3,2]. This modification is dominated by the change in the branching ratio of the decay channel T (1) → th, which is compatible with the fact that in most of the region T (1) is narrow. We can also see that the derivative operators can be tuned to strongly suppress the branching ratio of T (1) into th, setting it to essentially zero in a small region (in red on Fig. 5). It is also worth noting that the tthh cross section does not go to zero in that tuned region, because even if we suppress the resonant production, there is still the NR-tthh contribution, thus the absence of a red region in the left plot of Fig. 5.  of (c L , c R ), as well as the SM distributions for comparison. One notices that the different combinations of (c L , c R ) produce similar kinematic distributions, except for the one with (c L = 0.5, c R = 0.5) which is close to the fine tuned combination which suppresses T (1) → th. An identical situation is seen in other kinematic distributions and other points of the MHCM 5 parameter space, with the fine tuned combination of (c L , c R ) that suppresses the width being different for every MCHM 5 point.
The conclusion is that we have two qualitatively different situations: 1. For most values of (c L , c R ) the shapes of the distributions do not change and there is a change in the total cross section of a factor between 0.5 to 6 for 14 TeV and up to 8 for 100 TeV. This is equivalent to a k-factor and there is no advantage in including two extra parameters to that effect, since we are working at LO and a rescaling of cross sections is needed for comparison with experiment anyway.
2. For a small region around a particular value of (c L , c R ), which is different for every point in the MHCM 5 parameter space, the branching ratio BR[T (1) → th] becomes very small (sometimes even zero), and the distributions rapidly change into those of the non-resonant tthh production. It would be interesting to explore what happens to other decay channels near these points, but that is beyond the scope of this paper. We intend to return to this point in a future work focusing on the search for top-partners.  For the reasons given above, we will take c L = c R = 0, and ignore the effects of the derivative operators for the rest of our analysis. We have verified that the situation is qualitatively the same in the MCHM 14 .

MCHM at Low Scale
We now proceed to the scan of the parameter space for the MCHM. As mentioned before, we start by looking at the "Low Scale" region, which was characterized in section 3.2.1. We repeat the ranges here for convenience. For the MCHM 5 , we scanned over: For the MCHM 14 , we used: The results of the analysis of the scan over the parameters are presented in two ways, namely, the plots of a number of selected observables (Sec. 4.1), completed by the selection of a number of example-points including some of their relevant physics characteristics (Sec. 4.2) and a broader survey of the parameter space and its benchmark points, made in Sec. 4.3.

Scanning Over Parameter Space
We show in this section a number of selected observables that highlight the behaviour of the MCHM in the scanned regions of the parameter space. For completing the results of this scanning over parameters analysis, these plots include the selected example-points, discussed in details in the next subsection, with their labeling as defined in Tables 1 and 2 (Sec. 4.2). We will conveniently use the M 1 −M 4 space displayed in Fig. 2 to define regions in the parameter space for both scenarios under consideration. For the MCHM 5 we will define the following two regions 10 : and for the MCHM 14 we will define the following four regions: Each region is populated with about 200 points chosen at random, with points violating our conditions (see Sec. 3.1) being disregarded. Each point in the MHCM 5 is given by a choice of (f, |M 1 |, |M 4 |, sign(M 1 ), y L ) and a point in MHCM 14 by a choice of (f, |M 1 |, |M 4 |, |M 9 |, sign(M 1 ), sign(M 4 ), y L ). Each point is then passed to our implementation of the model in MadGraph, which gives us cross sections and distributions.
As shown in Eq. (3.7), the tth cross section is simply related to the SM one by a rescaling of the top Yukawa coupling. The deviations in the top Yukawa coupling from the SM limit have two distinct origins: • Deviations due to the compositeness nature of the Higgs boson, which arise from the dependence on the Higgs through trigonometric functions. This depends only on ξ, but is model-dependent and can in principle be used to distinguish the MCHM 5 from the MCHM 14 .
• Deviations arising from the mixing of the top quark with the new Q = 2/3 resonances. This effect depends on all the microscopic parameters of the model in a complicated manner through the diagonalization of the mass matrix. However, due to the fact that the resonances must be much heavier than the top quark, the deviations arising from the mixing are typically subdominant to the ones arising from Higgs compositeness.
One concludes that the tth cross section in the MCHM scenarios is largely controlled by a single parameter, which we can take to be the scale of global symmetry breaking, f . This is illustrated in Fig. 7. The signal strength µ(tth) of the production cross section is shown for √ s = 14 TeV but it does not depend on the CM energy at tree level, so the same results apply to √ s = 100 TeV. The CMS and ATLAS collaborations have presented results on the experimental measurements of µ(tth). Their reported best fits are: µ(tth) = 1.14 +0. 31 −0.27 for the combined 13 TeV result at an integrated luminosity of 35.9 fb −1 given by CMS [6] and µ(tth) = 1.32 +0. 28 −0.26 for the combined 13 TeV result at an integrated luminosity of up to 79.8 fb −1 given by ATLAS [5]. In Fig. 7 we show the 1σ and 2σ limits from CMS (as ATLAS has not reported on their 2σ limits).
In Fig. 7, one can see that, in Region I of the MCHM 5 , the points follow two distinct behaviours. The lower curve has most of its points (blue dots) corresponding to rather low M T (1) masses (below 1.5 or 1.6 TeV) and µ(tth) below 0.8, thus, in larger tension with the observed value while still within 2σ of it. Points with higher M T (1) (brown dots) are spread around a second curve with µ(tth) larger than 0.8 and with M T (1) greater than 2 TeV.
In contrast, in Region II of the MCHM 5 , there is only one smooth curve with points with a relatively small dispersion and equally distributed over the whole scanned range in M T (1) . The selected example points in the MCHM 5 are indicated in both Regions (P 1 , P 3 and P 4 in Region I and P 2 and P 5 in Region II). The important result is that in the MCHM 5 there is always a deficit in the tth production cross section as compared to the SM. This is, indeed, a main feature of the MCHM 5 .
For the MCHM 14 , three different cases are identified concerning the evolution of this variable versus f and M T (1) . Region I has some similarity with the corresponding Region I of the MCHM 5 . The main difference is that, in the MCHM 14 , µ(tth) can reach much smaller values (down to 0.2 if M T (1) is smaller than ∼ 1.6 TeV (blue dots), and down to 0.4 even for higher M T (1) masses). Thus, a fair fraction of all these scanned points have more than 2σ tension with the observed data. This case is represented by the example point Q 5 (see Table 2 in Sec. 4.2).
Regions II and IV of the MCHM 14 are very similar to each other, and also to Region II of the MCHM 5 , and are thus included in the same plot of Fig. 7. The main difference with the MCHM 5 case lies in a larger dispersion of the points and again the larger range in µ(tth) they cover (down to 0.2). Two example points, Q 6 and Q 7 , have been selected and they are shown in this Figure (see Table 2 in Sec. 4

.2)
The last MCHM 14 scenario for this observable refers to Region III, which deserves special attention. Fig. 2 shows an increase of tth production cross section as compared to the SM in a fraction of Region III. This is a main feature of the MCHM 14 as compared to the MCHM 5 or the MCHM 10 . This feature is clearly visible in Fig. 7, where Region III of MCHM 14 is the only one containing points with µ(tth) > 1.
To further explore Region III, a special scan with 100 additional points was performed, extending M 9 down to 1.3 TeV. All of them are gathered in the lower right plot in Fig. 7. There are two curves. One curve gathers a major part of the low M T (1) cases (blue dots). Some correspond to a dramatic deficit in tth production cross section getting near to zero. Most of the points corresponding to µ(tth) larger or close to 1 correspond to larger M T (1) masses i.e. masses larger than 1.8 -2 TeV (brown points). However some rare blue points (lower M T (1) masses) can also correspond to µ(tth) greater than 1. The highest µ(tth), above 1, are at relatively small f value, as expected. The selected Q 1 , Q 2 , Q 3 , Q 4 example points correspond to those cases but with a µ(tth) still within 1σ of the current LHC results.
The observable µ(tth) is thus a basic and key observable, not only to indicate that there is some BSM effect, but also to reject the MCHM 5 while keeping the MCHM 14 as still possible, if an enhancement w.r.t the SM is confirmed. If a deficit is instead observed, both MCHM scenarios will be possible, but the distinction between them is tricky and will depend on detailed phenomenology. More details are presented in Subsection 4.2 with the selected example points.
Turning now to the tthh process, we show in Fig. 8 how the signal strength µ(tthh) depends on the mass of the lightest Q = 2/3 resonance, for both MCHM scenarios and for different CM energies. There is a larger dispersion in the points of the MCHM 14 . However it must be noted that all the points in the four MCHM 14 regions (about 1000 scanned points) are included in a single plot whereas only 400 (2 regions) are gathered in the MCHM 5 case. These plots show the expected result that for lighter resonances the resonant production can result in a significant enhancement of the total tthh cross section. This effect becomes more prominent for larger CM energies.
14 TeV 100 TeV 150 TeV μ(tthh) 14 14 TeV 100 TeV 150 TeV These effects are highlighted by the example points. For instance, P 1 and P 3 (in the MCHM 5 case) are showing a large enhancement in µ(tthh) when increasing the CM energy whereas P 2 , P 4 and P 5 are not showing such an effect.
To complete the results shown in Fig. 8, Fig. 9 presents the ratio between the nonresonant contribution and the total cross section as a function of M T (1) merging the points of all the corresponding regions for each of the MCHM cases. The trends are quite similar between each MCHM scenario with a larger dispersion of the points in the MCHM 14 (with again the caveat of 1000 scanned points for MCHM 14 , versus 400 points for MCHM 5 ).  MCHM 14 14 TeV 100 TeV 150 TeV For heavier masses, the resonant production decreases and the total cross section can be dominated by the non-resonant contribution. Thus one expects deviations from the SM even when the resonances are rather heavy (say 3 TeV).
It is important to stress here that, except for cases with resonances close to the current direct search limit, we see that the non-resonant cross section accounts for a significant fraction of the total cross section. It is therefore of interest to search for deviations from the SM in this quantity, in addition to the dedicated resonant searches.
Finally, we show in Fig. 10 the NR-tthh cross section (normalized to the SM tthh cross section) as a function of the normalized tth cross section. There is a clear correlation, which reflects the fact that both are mainly controlled by the top Yukawa coupling, as explained above. One can note that the dispersion about the general trend is larger for MCHM 14 than for MCHM 5 .

Selection of some example points for each MCHM scenario
In order to illustrate the physics of the MCHM scenarios we chose a number of example-points. The selection criteria in the M 1 -M 4 plane (Fig. 2), takes into account the present experimental results including the LHC measurement of the tth production MCHM 5 14 TeV 100 TeV 150 TeV process [5,6], its prospects at the start of the HL-LHC, after the end of Run 3 with 300 fb −1 [51], and the exclusion limits on the pair production of heavy vector-like top partners currently obtained by ATLAS and CMS [8][9][10]. The plots resulting from the parameter scans in section 4.1 add an important input to this selection. The MCHM parameters characterizing each of the example points and the main observable quantities are summarized in Table 1 and Table 2, where we list the signal strengths for relevant energies, the spectra of vector-like fermionic resonances and the branching ratios of the lightest top partner. 5   Table 1 lists five selected points belonging either to Region I or Region II of the MCHM 5 , with the f scale ranging from about 900 GeV ("strong" compositeness) up to 2.45 TeV ("loose" compositeness) and with different values of y L . Note that similar scenarios can be found in either region. The chosen points are thus selected independently in one or the other case as reflected in Table 1.

Selected example points and their main features for the MCHM
The first point in Table 1 (P 1 ) shows a non-resonant contribution accounting for almost half the total cross section with a strong compositeness. The deficit in µ(tth) is relatively large and a bit borderline in regards to the estimated 1σ uncertainty on this measurement by the end of Run 3 (with at least 300 fb −1 ) [52,53], and to the rather low masses of the two lightest heavy top partner and the charge 5/3 resonance. This P 1 -scenario will be fully scanned (also including its overall resonances spectrum) at the HL-LHC where a 1σ uncertainty of 4.3% is expected on µ(tth) [54]. On the contrary, the relatively slight increase of the tthh production cross section with respect to the SM at 14 TeV, might not be reachable even at high luminosity and may even be visible only when the discrepancy with the SM further increases at higher CM energies. This makes this eventual scenario interesting to look at, even if possibly quickly disregarded at a certain stage of the HL-LHC run. Contrary to point P 1 , point P 2 shows a high NR-tthh and a clear deficit in µ(tthh); both effects are also visible for points P 4 and P 5 . However, P 2 , as it happens with P 1 , presents a deficit in µ(tth) because both have a rather low f value (strong compositeness). Points P 4 and P 5 , instead, have a high f value which translates into a µ(tthh) very close to 1.
Point P 3 has still a rather low f value with, as striking features, the strong increase in µ(tthh) at the expense of the low NR-tthh contribution (see dominant decay of the lightest resonance into th) and µ(tth) getting close to 1. All the expected resonances, in this case, have relatively low mass well reachable at the HL-LHC (and even may be before, i.e. by end of the forthcoming Run 3). The HL-LHC increased luminosity will allow to measure the branching decay especially into th, predicted to be dominant with respect to tZ. Besides, the full HL-LHC dataset could indicate a possible excess in µ(tthh).
Points P 4 and P 5 are rather similar in terms of all the measurable quantities listed in this Table. The NR-tthh contribution is 100% for both cases and will remain dominant even at higher energy accelerators. Moreover, while µ(tth) stays very close to 1 (due to a high f value, especially for the point P 4 ), the deficit in µ(tthh) could already be evidenced with the HL-LHC. Therefore, even if the points P 4 and P 5 look more like scenarios for the higher CM hadron colliders, a first breakthrough on such scenarios, especially for P 5 could be achieved by the end of the HL-LHC. Finally, note that the lightest resonance in both cases has a low branching ratio into W b (especially P 4 ), whereas a more important 3-body decay. This enhanced W + W − t channel occurs when T (1) comes from the fourplet and is thus almost degenerate with X 5/3 (both are controlled by M 4 ) as can be seen in the table. This effect will be specially important at higher energies, as we will discuss in detail in Sec. 5.
It is interesting to note that, as expected in these models (see Subsections 2.2 and 3.3), in all cases, the resonances show a mass degeneracy between two or three of them (MCHM 5 ) or even more (MCHM 14 ), in many cases due to EWSB. The separation in mass between those states can be of a few tens of GeV down to even a few hundreds of MeV.
The different scenarios described as example points for the MCHM 5 , present interesting features that allow distinguishing them from each other. They represent a variety of cases, covering different locations of the MCHM 5 parameter space; thus, they are interesting for exploring this Minimal Composite Higgs Model.

Selected example points and their main features for the MCHM 14
The MCHM 14 parameter space involves four different cases in the M 1 -M 4 plane (Fig. 2). A special attention is given to the Region III, as it contains the area with µ(tth) larger than 1. In this region, the M 9 value ranges from 1.3 up to 4 TeV, whereas in the other ones the lowest M 9 value is 2 TeV. The first four points (Q 1 to Q 4 ) in Table 2 are the selected example points for this region. The selection of the MCHM 14 points in the Region III requests, in addition to the criteria listed at the beginning of Subsection 4.2, µ(tth) to be larger than 1, this is the main difference between the two MCHM scenarios and also an important observable for the exclusion of the MCHM 5 . Note that all these points correspond to a relatively low f value, i.e., high compositeness.
The two first points are rather similar; they both correspond to low f and M 9 values. Q 1 is close to the 1σ experimental limits w.r.t the µ(tth) value. The three lightest associated resonances are already almost within the limits published by ATLAS and CMS. However, this point is a good example of a high increase in µ(tthh) (like also observed in some examples of the MCHM 5 ) and it includes 50% of non-resonant contribution at 14 TeV. If one disregards the resonant contributions, it primarily differs from an equivalent MCHM 5 scenario by the increase in µ(tth), as compared to the SM.
The model parameters of point Q 2 have values very close to the ones of Q 1 but a slightly higher f value (lower compositeness); it thus translates into a smaller µ(tth) value. Moreover, µ(tthh) decreases compared to Q 1 , although it still stays relatively high.
Points Q 3 and Q 4 , have both a relatively low f value but higher M 9 values with about 2 TeV (for Q 3 ) and 2.7 TeV (for Q 4 ). Both have µ(tth) greater than 1 but well within the current experimental limits. Their selection is also based on the request for a high non-resonant contribution.
For completeness, the remaining three regions of the M 1 -M 4 parameter space were considered. In each of them, a representative point is selected as summarized in Table 2: The point Q 5 in Region I, Q 6 in Region II and Q 7 in Region IV.
In the overall covered space these 3 regions provide quite similar cases. Moreover, as shown in Fig. 2, a subregion of Region IV is excluded because of the constraint on the ggh coupling (0.8 < ggh/SM < 1.2). The three example points have different f values (around 1 TeV for Q 5 and around 2 TeV for Q 6 and Q 7 ). They all have a relatively large M 9 (3 TeV). The example points in these three regions show a µ(tth) value smaller than one and no strong increase in µ(tthh). The selected points Q 5 and Q 7 have large NR-tthh contributions. Point Q 6 has only 50% NR-tthh contribution. All the NR-tthh relative contributions decrease sharply at higher CM energies, as more phase space becomes available for the production of resonances.
The use of these preliminary observables shows that it will be difficult to disentangle between both MCHM scenarios if a deficit in µ(tth) is measured. A much more detailed analysis will be required. In some cases, the HL-LHC will perhaps provide a first indication, but a potential discovery will likely need higher energy together with higher luminosity.

NR-tthh contributions in the MCHM 5 and the MCHM 14
In order to clarify the different contributions to the non resonant tthh production, we simulated these contributions separately and summarized the results in Tables 3 and 4. The ratios in those tables are obtained by turning off one or more couplings in the model in order to disregard particular classes of diagrams, which are indicated in the table.
The σ hh /σ NR ratios show that the effects of the double Higgs Yukawa coupling are typically at the couple to few percent level in the MCHM 5 and the MCHM 14 and hardly show any variation with CM energy increase. We also find, by examining the σ Yuk /σ NR ratios, that the effect of the trilinear Higgs self-interaction can be around 15% in both MCHM 5 and MCHM 14 . For comparison, the effect of the trilinear Higgs self-interaction in the SM tthh cross section is about 20%, with a very mild CM energy dependence. Thus, it is largely the top Yukawa that governs the NR-tthh (just as in the SM), which, to a first approximation, then scales as (y t /y SM t ) 4 . This correlation explains the behavior  Table 3. Study of NR-tthh for the MCHM 5 points in Table 1. The cross sections σ hh and σ Yuk are obtained by disregarding the classes of diagrams on the last column and σ NR is the total NR-tthh. The LO SM tthh production is indicated by σ SM and σ SM Yuk means we disregarded the SM trilinear Higgs coupling. The top Yukawa couplings are indicated by y t and y SM t in the MCHM and SM respectively.  Table 4. Study of NR-tthh for the MCHM 14 points in Table 2. The cross sections σ hh and σ Yuk are obtained by disregarding the classes of diagrams shown on Table 3 and σ NR is the total NR-tthh. The LO SM tthh production is indicated by σ SM and σ SM Yuk means we disregarded the SM trilinear Higgs coupling. The top Yukawa couplings are indicated by y t and y SM t in the MCHM and SM respectively.
in Fig. 4 and is evident in the last three lines of Tables 3 and 4, where the enhancement or suppression in the cross section directly arises from the change in the top Yukawa.
The low contributions of the trilinear Higgs and the double Higgs Yukawa couplings present challenges. It makes the measurement of these two couplings harder, and both are important to characterize the compositeness of the Higgs (as opposed to the spectra of fermionic resonances) and are even harder to probe on the resonant production. Further-more the shape of all kinematic distributions of the NR-tthh final states will be almost identical in shape to the SM one, changing only on the magnitude of the integrated cross section. The top Yukawa can be more directly accessed in the tth channel, but it will be necessary to develop combined analyses between the tth and the NR-tthh to isolate the non-resonant contribution and extract information about these couplings.
Let's stress here a key-role of tthh production process in the study of the self Higgs coupling and the important interplay between the measurement of the κ λ , the tri-linear Higgs coupling normalized to the SM value, through the hh process (both gluon fusion and Vector Boson Fusion) and of the tthh process. Even if the trilinear Higgs contribution to the tthh process represents only about 15% (slightly lower than in the SM case) of the total cross-section in the MCHM considered scenarios, it is not as difficult to access as the double Higgs Yukawa contribution at the percent level, and indeed it is worth it. The measurement of the κ λ parameter of the Higgs sector via several processes is becoming of increasing importance as it remains the experimentally least constrained Higgs parameter. This is due to the fact that the "traditional" way to access it at LHC, via the hh production, is currently challenging, because of the still relatively low crosssections and signature efficiency. In the years to come and over the whole HL-LHC era this will be indeed an essential experimental goal. In order to increase the experimental sensitivity reach, ATLAS and CMS experiments are already now searching, in addition to the hh production through the gluon fusion process [55,56], for the production process through the Vector Boson Fusion, VBF(hh) [57] and [58]. Recently the need to look for the complementary contribution of tthh on this specific topic is outlined both in view of the HL-LHC and even more of the FCC-hh at 100 TeV [59,60]. It is worth noting that the SM cross-section of the VBF(hh) and of tthh processes are of the same order, i.e. at the fb level at 14 TeV. Besides, while both ATLAS and CMS are conducting searches for the VBF (hh), a search for the tthh production has been carried on, for the first time, at CMS [7], stressing the growing interest on the experimental side.
To conclude, from the phenomenological viewpoint the tthh channel does not include destructive interference among diagrams unlike in the hh case (see e.g. [61]), and if λ > λ SM , it provides the leading channel where to observe an excess over the SM expectation. From the experimental viewpoint, the two tops in addition to the Higgs pair strengthen the signal efficiency with respect to the hh and hhjj signatures, making it accessible already now at LHC and furthermore at HL-LHC. Besides, the large increase with energy of its cross-section in the SM and even more in the MCHM case makes it an essential channel to study at FCC-hh for high precision and BSM measurements of the Higgs sector.
As a final remark, we wish to point out that the reader should be careful when looking for the separation between the resonant and the non-resonant production in Fig. 4. The separation looks clear for √ s = 14 TeV and difficult for √ s = 100 TeV, giving the wrong impression that there is little hope for exploring the NR-tthh at future accelerators. That only happens because the same point in parameter space was used for both plots, and that is a point with low lying fermionic resonances which are produced abundantly at higher energies and dominate over everything else. We explored also the hypothetical situation where no resonance was found at the end of the HL-LHC run, forcing us into regions of the parameter space where the top partners are heavier and the non-resonant production is more pronounced, as we shall see in Section 5.

Cluster Analysis applied to the MCHM at low scale
In the next two subsections we apply the clustering idea to the MCHM 5 and the MCHM 14 , using the parton level kinematic distributions of the tthh process to do the cluster analysis described in Section 3.4.

Clustering of the MCHM 5
In the case of the MCHM 5 we start with the points generated by the scan of Sec. 4.1: 400 points divided between Regions I and II. We then apply the following "cuts" to remove points that are already constrained at a 3σ level 11 0.33 ≤ µ(tth) ≤ 2.07, (4.1) We also check if the points in the scan are excluded or not by the experimental measurements of the ttZ and ttW production cross sections. The latest measurement of the ttZ production, performed by the CMS collaboration, corresponding to an integrated luminosity of 77.5 fb −1 and at 1σ level is σ(ttZ) = 0.95 ± 0.08 pb [62]. The ttZ signal strength was calculated dividing the measured cross section by the SM prediction σ SM = 0.84 ± 0.10 pb, obtaining µ(ttZ) = 1.13 ± 0.16. The latest measurement of the signal strength of the ttW production, performed by the ATLAS collaboration, corresponding to an integrated luminosity of 36.1 fb −1 is µ(ttW ) = 1.44 ± 0.32 at 1σ level [63]. All the scanned points survived these constraints at 3σ level 12 .
That leaves us with 348 points at the start of the clustering algorithm. Each one of these points was implemented in MadGraph 5 and events were generated for the production of tthh at a fixed luminosity of 3000 fb −1 . At parton level we have only three different particles present: the top, the anti-top and two Higgs bosons (we focused on the most energetic one). Using MadAnalysis we obtained histograms for the following kinematic distributions:  11 Neither CMS nor ATLAS report the 3σ error directly, so the best we can do here is to assume a Gaussian error and estimate the 3σ threshold simply by multiplying their 1σ intervals by 3, using the ATLAS/CMS combined value when available and the smallest one otherwise. 12 Once again 3σ is roughly estimated as three times the 1σ intervals.
We then constructed samples from these distributions using single distributions or combinations containing two or three distributions. Since the main contribution to the total tthh cross section is coming from the decays of top partners, all these kinematic variables are strongly correlated and that means that most combinations led to very similar clusters.
Initially one is led to consider only the invariant mass, as this is the one variable that makes the resonant structure evident, but we find that using only M [t, h 1 ] for the clustering puts too much weight on the exact position of the peaks produced by the resonances, instead of more general behaviours like the two peak structure that shows up in the p T distribution of the Higgs. The result would be to have many benchmark points in the region with lighter and narrower resonances, and just one or two in the rest of the parameter space. Including at least one of the angular variables brings extra physical information and takes away that emphasis, resulting in benchmark samples that are more evenly spread. There is very little difference in regards to what angular distribution we choose, but the results with θ[t] resulted in samples being more evenly distributed among clusters.
We also experimented with the number of clusters and found out that with a small number of clusters (N cluster < 10) we obtain one highly populated cluster that contains all the samples characterized by heavier (M T (1) 1.5 TeV) and broader resonances. This cluster results from the merger of two large clusters when we go from 11 to 10 clusters, but these two are already generated early on the clustering process, which means that stopping with N cluster much bigger than 11 only changes the low population clusters, that are already very homogeneous, so there is not much gain in increasing N cluster .
We finally decided to stick with N cluster = 11 and on using M [t, h 1 ] and θ[t] for the cluster analysis of the MCHM 5 , the resulting clusters are shown on Figs. 11 and 12, where we also included p T [h 1 ] distributions to show their typical behaviour. We first note that most of the clusters are very homogeneous in distributions, in the sense that all curves in each plot are very similar. The benchmark points (black lines in the plots) will then be very good representatives of the behaviour of each cluster. This is true even for the p T [h 1 ] distributions, which were not used as a criteria for clustering, and also for all the distributions we have checked (listed above). The only striking exception is cluster 3 (in Fig. 11), which contains all samples with heavy fermionic resonances, with M T (1) 1.65 TeV, or no resonances at all decaying to tthh. That is a consequence of the comparatively low count of events in the resonant region, when compared to the non-resonant part of the distribution, which is very similar to them all. This could probably be fixed by doing a dedicated scan for points in that region and a separate clusterization, producing more clusters and benchmark points. We decided against it because we already covered that region with example points P 2 , P 4 and P 5 of table 1, which were in fact grouped in cluster 3 and are shown as green curves in Fig. 11. Figures 11 and 12 also allow us to survey the main features of the whole parameter space of the model. In particular one can see that points in Region I (colored in red) got separated from those of Region II (in blue). Region I is concentrated in clusters 3 to 7, while Region II dominates the rest of the clusters. One can also verify that the fermionic resonances produced in Region I are generally wider than resonances in Region   Table 1 are shown in green (P 3 is in cluster 8 and P 1 is in cluster 11).
II, this can be more clearly seen in clusters 4, 5 and 7 where resonances from both regions overlap. In order to understand this, it is useful to look at an approximate expression for the top mass, given in Eq. (2.25). This mass is generated in the MCHM 5 by mixing of the elementary fermion fields with the fourplet and singlet resonances. Looking at the first two diagrams in Fig. 1, as well as the Lagrangian of Eq. (2.10), we see this mixing leads to an insertion of y L times y R , times a mass insertion of M 1 or M 4 for chirality.
These diagrams must interfere destructively, since the pNGB Higgs vacuum misalignment is generated by SO(5) breaking and hence must vanish in the SO(5) symmetric M 4 = M 1 limit. This leads to a dependence m t ∼ y L y R |M 4 − M 1 |, as shown explicitly in Eq. (2.25). In Region I, there is a cancellation between same sign M 1 and M 4 , such that larger values of y L,R are typically needed to generate the top mass 13 . This enhanced mixing also leads, upon mass diagonalization, to a greater value for the t T (1) h vertex, and hence, to wider resonances.
Another interesting general feature is the presence in many cases of more than one peak in the M [h 1 , t] distribution. While cluster 9 represents well the usual simplifying assumption used in top partner searches, namely the presence of only one resonance decaying to the th, tZ or bW channels, many of the other clusters contain a sizeable presence of more complicated peak structures, coming specially from Region II. Cluster 10 is the perfect example of this, as its benchmark point has a double peak structure with the second resonance giving a stronger contribution than the lightest one. Most exclusion limits for top partners are obtained through analyses optimized for the situation in cluster 9, and it would be interesting to see how those limits change if more resonances are considered, specially if they overlap significantly.
In Fig. 13 we show how the benchmark points are placed in the parameter space, together with the example points of Table 1 and the rest of the points in the scan not excluded by constraints. We can see that the benchmark points, complemented by the example ones, are well distributed in the parameter space. We finally list the benchmark points and their main features in Table 5, where we can verify many of the features visible in the distribution plots of Figs. 11 and 12. Points in Region I (C 4 , C 5 , C 6 and C 7 ) have on average higher couplings y L and y R and wider T (1) than those in Region II, although the extreme cases in each region can be similar. The point C 7 , which contains a narrow T (1) for Region I standards, is quite similar to C 10 , in which T (1) is exceptionally wide for Region II.
Another striking feature of Table 5 is that for all of the benchmark points (but C 3 ) there is at least a 10% branching ratio in 3-body decays. Similarly striking and linked to the previous observation, we note that for all points but this same C 3 point, the 2-body decay into bW is very small, namely between 0% and 2.4%. That is relevant as most top partner searches were done under the assumption that the three 2-body channels (th, tZ and bW ) comprise the full width. We will explore the phenomenology of these non-standard branching ratios in a future work. 14 We applied the same clustering method to the MCHM 14 , using again Eqs. 4.1 and 4.2 as constraints. The ideal clustering, following the criteria of homogeneity within clusters while keeping the number of clusters small, was obtained including M [t, h 1 ], p T [t] and ∆R[h 1 , t] in the samples and stopping at 12 clusters. The distribution of the benchmark points can be seen in Figure 14 and their main properties are listed in Table 6. We omit the plots of all distributions and clusters as the general features are very similar to  Table 1, benchmark points of Table 5 and the points consistent with constraints in Eqs. 4.1 and 4.2.  the MCHM 5 , with the points without light resonances decaying to tthh being grouped into less homogeneous clusters (specially clusters 4 and 9), but all the plots are available online [64]. There is still a tendency towards narrower resonances when M 1 and M 4 have opposite signs (Regions II and IV), as expected due to the 1−r 1 term in Eq. 2.25, but now there is a ξ-proportional correction that is also sensitive to the sign of M 4 (as M 9 > 0) and makes this tendency weaker (and hard to notice if one looks only to the benchmark points).  Figure 14. Low scale scan of the MCHM 14 parameter space, including example points of Table 2, benchmark points of Table 6 and the points consistent with constraints in Eqs. 4.1 and 4.2. Figure 14 shows that, between the example points of section 4.2.2 and the benchmark points obtained here we are covering the parameter space quite well. One can see that Table 6 contains many points that have a sizeable 3-body decay in all regions and, like in the case of the MCHM 5 , a fair amount of the scanned points have overlapping resonances (which is the case for the benchmark points of clusters 3, 8 and 11), both of which are relevant concerns for top partner searches and constraints.

Clustering of the MCHM
We would finally like to stress the presence of the benchmark points D 1 and D 9 , both with a marked increase in σ(tth), showing that such a situation is not uncommon in the MCHM 14 and is a really interesting possibility to evidence this realization of the MCHM.

MCHM at High Scale
We continue the analysis of the two studied MCHM scenarios, extending the dimensionful parameters to higher scales that can be accessible only in the context of the high CM energy pp colliders in project [11][12][13][65][66][67]. A higher CM energy around 100 TeV is, indeed, requested to confront scenarios with very high mass resonances, with a possible deviation of the tth cross section from the SM value at the percent level (or even less),  (GeV)  -1382  2000  3889  3236  2836  3473  1500  1467  2000  3230  2965  1329  f(GeV)  881  1275  1012  1288  1550  1912  863  931  1071  1648  1155  1244   with the requested precision to study the tthh process (e.g. the various NR components) and to measure the branching ratios of various decays of the produced resonances, even if, for instance, the lightest one is detected at the HL-LHC.

Scanning Over Parameter Space
The numerical strategy we use here is the same as the one we followed for the Low Scale analysis and described in Subsection 3.1 but with the following extended range of parameters in order to perform a high scale scan. For the MCHM 5 we consider: We divide the parameter space of the MCHM 5 and the MCHM 14 in the same regions we used for the Low Scale analysis (Sec. 4.1). In each region for each model we scanned 200 random points. We present the main results of the scans in Fig. 15 and Fig. 17, where we joined all the possible regions for each model.   (1) , with color coded values of f . As explained in Subsection 3.7, at low scales the main contribution to this cross section comes from the QCD vector-like pair production. However, the later process undergoes a quick drop for heavy resonances, so it is expected that beyond certain value of M T (1) the QCD resonance pair production becomes a subdominant part of the total tthh cross section. In the plots of Fig. 17, we see that this value is about M T (1) ≈ 4 TeV as the cross section becomes basically independent from M T (1) beyond that point. In the left plot, corresponding to the MCHM 5 , above this value there is no point with a µ(tthh) bigger than one, meaning that the main contributor is now the non-resonant part of the process, which, as we saw in Fig. 10, is directly related to µ(tth) and always suppressed w.r.t the SM. Since the later is, in turn, mainly controlled by the compositeness scale f (see Fig.15), the points with smaller values of f have a more strongly suppressed cross section and as it increases, the cross section approaches the SM limit.
We see that, in the right plot of Fig. 17 (corresponding to the MCHM 14 ) there are still points with considerable enhancements with respect to the SM beyond M T (1) ≈ 4 TeV. There are two kinds of points with this behaviour: some points present enhancements in the tth Yukawa coupling which reflect in enhancements in the NR-tthh and consequently in an increase of order 10% in tthh. The second case is more subtle and is responsible for the biggest cross sections in the region above M T (1) ≈ 4 TeV. Those are fined tuned points, in which the parameters lead to strong coupling interaction between the top, a vector-like resonance and the Higgs. This vertex enters in diagrams like the one in Fig. 16, in which a single intermediate top partner is produced (not to be confused with weak single production), increasing the cross section as well. It is important to mention that the second type of points usually has a low compositeness scale (f < 2 TeV) and can be constrained by the increased precision on the measurement of the top Yukawa attainable by the HL-LHC (see Fig 15).

Cluster analysis applied to the MCHM at High scale
Following the analysis of Secs. 4.2 and 4.3 we look now at specific points in the parameter space in order to examine the phenomenological features of the MHCM at high scale. As we show below, the increased mass range in the scan demanded a slightly different strategy in regards to the clustering technique of Section 3. 4. In what follows, we divided the points of the scan in smaller sets before clustering. This approach, together with the fact the phenomenological behaviour at higher scales is more uniform, provides enough points to showcase all the interesting features of the model. For that reason we have not chosen additional example points for the High Scale scan, and we proceed directly to the results of the clustering algorithm. 5 We start with the 400 points of the scan in the previous section (200 in each region) and discard all points violating the following constraint on κ t , the top Yukawa coupling normalized w.r.t. to the SM: 0.9 < κ t < 1.1 (5.1) This constraint is based on the projected precision on the top Yukawa measurement in the High Luminosity phase of the LHC. In [51] the 1σ uncertainty on κ t is projected to be 3.4% for an accumulated luminosity of 3000 fb −1 . We constrain our points to the region that is 3 times that uncertainty around the SM value κ t = 1, as a rough estimate of the 3σ region. We do this in order to keep points that have a smaller chance of being constrained by measurements at the time when the 100 TeV pp collider starts its operation. We can directly translate eq. 5.1 into constraints in µ(tth), and we show these limits as red dashed lines in Fig. 15, where one can see that this allows for fairly small values of µ(tth) ( 0.8).

Clustering of the MCHM
By the end of the HL-LHC phase the points close to that limit will be of course at the edge of the allowed region, and could be even excluded, specially if the central value of µ(tth) turns out to be above 1 or the precision is better than expected. But on the other hand, it is still a possibility that the opposite happens (central value of µ(tth) < 1 for instance) so we understand it is too early to disregard those points. The remaining points are clustered following the method described in Section 5.1. Different combinations of kinematic distributions are considered but, in all of them, the clustering procedure was giving too much weight to the position of peaks on points with lighter resonances, creating many low population clusters for those and leaving the ones with heavier resonances combined in a unique cluster. This is similar to what happened in the clusterization of the low Scale MCHM 5 (see cluster 3 in Fig, 12), but the problem is exacerbated by the fact that, in this case, we are working with a bigger range of parameters and, consequently, with a wider variety of resonance masses. As we are interested in selecting points distributed in the whole parameter space, we decide to first group the points in slices, guided by the mass of the lightest resonance, and then to apply the clustering algorithm to each slice individually. The following five slices, are used: In this way we ensure we are clustering the points in a more homogeneous way. The last slice groups all the points with resonances heavier than 6 TeV since, as we see in Fig. 17, it is expected that all of them possess similar kinematic characteristics because the contribution of the resonances for them are negligible. The selected clustering is realized by using the M [t, h 1 ], p T [t] and θ [t] distributions and stopping at 2 clusters in each slice, for a total of 10 clusters. The plots showing the clustering in full detail are available online [64]. The benchmark points corresponding to each cluster are listed in Table 7.  The first general feature to be noticed in Table 7 is the importance of the NRtthh contribution to the total cross-section. Only in point E 2 it is not the dominant contribution, although it is still about half of the cross section. Even in the points with resonances in the 2 ∼ 3 TeV range the NR contribution is high (0.7 in E 1 , ≥ 0.9 in E 3 and E 4 ). As the resonances get heavier, the NR component becomes almost all of the cross section. The increasing suppression of the T (1) pair production also causes µ(tthh) to be close to or below 1 for the points with heavier ressonances (E 3 to E 10 ), as expected for the NR contribution in the MCHM 5 .
The second striking feature is the number of points with a sizeable or 3-body decay (all points with the exception of E 6 and E 9 ). The effect happens when both M T (1) and M X 5/3 are largely controlled by M 4 , and are thus almost degenerate (whereas in E 6 and E 9 , M T (1) is closer to M 1 ). In these cases, there is also a marked suppression of BR(T (1) →W + b). This was already present in the low scale scan but becomes more prominent here, and one can check in the Table that the 3-body decay becomes more important as M T (1) increases. Considering the points with this property (which are the majority of the benchmark points) we see that the 3-body branching ratio increases from ∼ 10% at M T (1) ≈ 1.3 TeV (see Table 5) to quickly becoming the dominant decay channel around M T (1) ≈ 3 TeV, being already 25% in point E 2 (with M T (1) = 2.1 TeV). This makes it extremely important to take three body decays into consideration if one wishes to push the direct search exclusion of top partners beyond the present ballpark of 1.3 TeV. Additionally, the fact that this decay happens through the exotically charged X 5/3 , also highlights the importance of examining more complete models of the fermionic sector, as opposed to simplified models that only include a single electroweak doublet or singlet top partner.
On the other hand, points E 6 and E 9 exhibit the behaviour typically well captured by simplified models. In these, M T (1) ≈ M 1 is well separated from the rest of the spectrum, which is dominated by masses close to a much higher M 4 . The branching ratios follow then the expected pattern: the W + b channel branching ratio is twice that of the th and tZ channels.
Other features already present in the low scale scan are also present at higher scale, specially the presence of double peaked structures in many points. Among the benchmark points, we highlight points E 5 and E 10 where we find two top partners close in mass. In fact in these two points, the exotic X 5/3 and the bottom partner B (1) are also close by, allowing the 3-body decay to go through either of those channels. This causes the width of T (1) to increase a lot: (Γ/M ) T (1) = 24% for E 5 and even up to (Γ/M ) T (1) = 87% for E 10 . Points E 7 and E 8 also exhibit similar behaviour, but there is a bigger mass gap and smaller (Γ/M ) T (1) . Figure 19 shows the scanned points as well as the benchmark points in the M 1 -M 4 plane with color coded f . There, we see that the benchmark points obtained with the clustering technique are covering both regions of the model. Besides verifying the homogeneity of each cluster we also searched for points that have behaviours that deviate significantly from the benchmarks, and found none. This indicates that the behaviour of the model is quite uniform and well represented by the benchmark points, even in the region with a small number of them, as for instance in Region II (negative M 1 ). The same is true for regions with high M 1 , M 4 , and specially high f , as we are moving towards the decoupling limit of the model and the points become more similar to each other and to the SM (which explains why only one point, E 3 , was chosen to represent the high f region). For these reasons, it would be superfluous to choose and include additional example points.

Clustering of the MCHM 14
In the case of the MCHM 14 at high scale we started with 800 points evenly distributed in the four regions of the M 1 -M 4 space as described in Section 5.1. Then, we selected all points allowed by the constraint in Eq. 5.1 and proceed with the clustering technique. As in the MCHM 5 , we divide the points in the five slices of Eq. 5.2. We stopped at 2  Table 7 represented by triangles and the points satisfying the constraint in Eq. 5.1. The compositeness scale f is color coded.
or 3 clusters in each slice, depending on how homogeneous were the obtained clusters, and the corresponding benchmark points are listed in Table 8. The plots showing the clustering in full detail are available online [64]. Many features are similar to the MCHM 5 : the NR-tthh once again becomes dominant once M T (1) goes above 3 TeV, and both µ(tth) and µ(tthh) go below 1 when that happens. In principle, the MCHM 14 could allow for an increased top Yukawa and thus µ(tth) > 1 and µ(tthh) > 1 even in the non-resonant regime. Only one point with that behaviour was selected by the algorithm (F 11 ), but all scales for that point are so high that it is almost SM-like.
In regard to the spectrum and decay of the lightest top partner, many points mimic the behaviour of the MCHM 5 , but the situation is richer here: • The points F 1 , F 3 , F 7 , F 9 and F 10 reproduce the dominant behaviour in the MCHM 5 , with the masses of T (1) and X 2/3 essentially set by M 4 , large branching ratios in 3-body decays and suppressed W b decays. All the important remarks made in that case apply.
• The points F 5 and F 8 have the split spectrum similar to E 6 and E 9 . The top partner has a mass set by M 1 and the rest of the spectrum is much heavier. Decay channels follow the pattern assumed in simplified models.
• The points F 2 , F 4 and F 6 are qualitatively new. The mass scale is largely set by M 9 and we have three top partners that are degenerate not only among themselves,  but also coincide with the masses of the lightest exotically charged fermion and the bottom partner. The width of T (1) does not increase as in the degenerate cases of the MCHM 5 , and there is a suppression of the th channel. A complete survey of all the decay channels of T (2) and T (3) would be needed to disentangle these states (and determine which are contributing to tthh and in what degree), but that is beyond the scope of this work. We intend to explore this point in a future study.
Moroever, unlike the points F 4 and F 6 , the point F 2 , if still valid after HL-LHC, will be evidenced at a 100 TeV machine thanks to the high precision measurement of µ(tth) and, more importantly, to the large deviation from the SM in µ(tthh).
• The point F 11 is really SM-like from several experimental viewpoints; even with the high precision reachable on µ(tth) and µ(tthh) measurements at a 100 TeV collider, it will not be possible to disentangle them from SM. Besides, the very high mass and highly degenerate in M 4 and M 9 resonance spectrum, only shows a highly dominant 3 body decay of T (1) , as particular features. Let's note that points F 5 and F 8 have the same highly degenerate in M 4 and M 9 resonance spectrum but, as mentioned above, with a pattern of decay channels similar to simplified models. The F 11 point would deserve some specific phenomenological and experimental study, that goes beyond the scope of this paper. Figure 19 shows the benchmark points in the M 1 -M 4 plane with color coded values of f . The benchmark points obtained with the clustering cover this space homogeneously and then we decided also not to add example points.  Table 8 represented by triangles and the points satisfying the constraint in Eq. 5.1. The compositeness scale f is color coded.

Effective Field Theory Perspective
In the analysis reported in the previous sections, we always consider the effect of the fermionic resonances in the cross sections, however, for the points in parameter space with larger masses, decoupling occurs and the processes are then largely non-resonant. In this case, it can be useful to present expressions for the modifications of the SM couplings, in the context of an effective field theory, in which the heavy degrees of freedom are out of experimental reach and can be safely integrated out 14 , in order to facilitate a comparison with other non-resonant studies in the literature.
The aim of this section is to more generally show the interplay between the considered MCHM scenarios in their overall scanned parameter spaces, and the EFT framework. To do so we present our results in terms of modifications to the SM couplings, which are more directly comparable to the increasing number of experimental results on these parameters. It is straightforward to express these results in terms of a non-redundant operator basis such as the strongly-interacting light Higgs (SILH) basis [68]. The Higgs effective Lagrangian we consider here, after EWSB and neglecting the light fermion interactions, is given below: where κ λ , κ t , c 2 , c g and c 2g are coefficients that encode the modifications to the SM. We do not present here the expression for c 2g as this operator will not contribute to the tth and tthh processes at tree level. The coupling κ λ is the same in both the MCHM 5 and the MCHM 14 : 15 For the MCHM 5 , the remaining couplings are given by: Here we defined θ L and θ R by tan θ L = y L f /M 4 and tan θ R = y R f /M 1 . For the MCHM 14 , the corresponding couplings are given by: with c t g and c b g defined in equations 2.29 and 2.30, respectively. In Figs. 20 and 21, we have plotted selected EFT parameters against each other, for the points scanned in Sections. 4 and 5. The colors identify different regions in the parameter space of each model for Fig. 20, which covers the low scale region of the 15 Some of the expressions below are reported also in [19], and presented here again for completeness. the MCHM 14 . For instance, it can be seen that the aforementioned enhancement in κ t is only present for Region III in the 14, while the relation κ λ = c g is only satisfied for the MCHM 5 ; therefore, a measurement of κ t > 1 or κ λ = c g would argue strongly against the MCHM 5 . Evidence of non-zero c 2 is also particularly interesting since this vertex is absent in the SM. As expected, for smaller values of the compositeness scale f , the deviations from the SM are of course larger, and the separation between the different models is easier.

Conclusions and Outlook
The Minimal Composite Higgs Models MCHM 5 and MCHM 14 are studied here as an important show-case for the exploration of the beyond standard Model world, in the top-Higgs sector at high energy hadron colliders. The full generation and simulation framework is developed which allows a phenomenological analysis over the whole parameter space of these two models. Besides, the developed generation and simulation software defines a preliminary experimental framework; they are ready for use by more detailed data analysis at the LHC and HL-LHC experiments and by advanced detector simulations developed for the 100 TeV machines in project.
The focus is on the analysis of the tth and tthh production processes, covering in the later both the production of heavy fermions and the non-resonant contributions. In this paper we have gauged the relative contribution of these two cases to tthh and shown that the non-resonant is sizeable, in fact dominant once the fermions are heavier than 4 TeV (even in the case with 100 TeV of center of mass energy), and gives access not only to the trilinear Higgs self-coupling, as it does in the SM case, but also to the double Yukawa coupling (tthh vertex) which is introduced by the MCHM. The relative contributions of the top-Yukawa, the trilinear Higgs and the double Yukawa couplings were also analysed, showing that the trilinear Higgs coupling contributes approximately with 15% of the non-resonant cross section (see tables 3 and 4).
A systematic exploration of the parameter space of both models was conducted. Using a clustering algorithm, it was possible to find a small number of benchmark points that showcase the phenomenology in a comprehensive way. These points were then complemented by exceptional points where needed, and the most important phenomenological data is summarized in tables 1, 5 and 7 for the MCHM 5 , and tables 2, 6 and 8 for the MCHM 14 . This exploration was done in two steps, first covering physics at "Low Scale", where we focus on the part of the parameter space that is within reach of the HL-LHC. On a second step we extend the analysis to the "High Scale" region, with mass parameters reaching up to 30 TeV, which will be of interest for planned future colliders such as the 100 TeV FCC-hh and the SppC.
The first observations from the phenomenological and preliminary experimental analysis are that a deviation from the SM in the tth production is also an essential measurement for MCHM. An increase will reject the MCHM 5 scenario and greatly refine the areas of the parameter space where MCHM 14 would be valid. A deficit instead, would make MCHM 5 and MCHM 14 both possible. The measurement of this observable is expected to be achieved with 3.4% accuracy at the HL-LHC and thus with a very high accuracy (at least at the percent level) at the future high energy hadron collider in project at 100 TeV.
Regarding the search for fermionic resonances that are present in the partial compositeness scenario we explored, we find that in most of the benchmark points the 3-body decay channel of the lightest top-partner starts to become increasingly important as its mass grows. For all the benchmark points in this category, the 3-body decay is already 10% of the branching ratio when the mass is around 1.3 TeV, increasing to around 25% at 2 TeV and becomes the main decay channel when the top-partner reaches 3 TeV. The same points also have a marked suppression of decays to W + b. As expected (see e.g. [39]), we find that it is not so common to have a "split spectrum" with a low lying toppartner separated in mass from other fermionic resonances since a full SO(4) multiplet is controlled by a single mass parameter. This results in complicated interplay between the many states present, and is part of the reason for the increase in three 3-body decay channel. This makes it extremely important to take 3-body decays into consideration in future top-partner searches, and we intend to explore the decay patterns of all the fermionic resonances, and the resulting search strategies, in a upcoming work.
A comparison with the EFT, valid in energy regimes where the masses of the fermionic resonances are not reached, is also provided (see figures 20 and 21). This is specially interesting in the situation we are right now, and most probably will be at least until the end of the HL-LHC phase, in which no new states have been discovered but great improvements on the measurements of the coupling constants are being made. If deviations from the SM in more than one of these couplings are found, the specific combination of deviations can be used to differentiate not only between the MHCM 5 and the MCHM 14 , but also between different regions of those models, as these regions generate different correlations between the effective couplings.
For the process tthh, both MCHM scenarios can deviate significantly from the SM expectation, either as a deficit or an increase. For this reason also this channel will play an important role in searches of composite Higgs models. The issue is the relatively low SM cross section of about 1 fb at tree level at 14 TeV, whereas tth is roughly a factor 500 higher. Thus, the aim at HL-LHC will be to evidence this process and to get a first indication of a strong deviation from the SM. For really exploring MCHM, higher energy together with higher luminosity as foreseen in the pp colliders in project towards the second half of this century (100 TeV and exceeding 20 ab −1 total integrated luminosity) are not only a plus but indeed even a necessity.

A Representations of SO(5)
We use the following 5 × 5 matrix representation of the generators T B of SO (5): where we used the normalized basis v 1 and v T 0 = (0, 0, 0, 0, 1) .

B Embeddings of SO(4) into SO(5)
The four NGBs resulting from the spontaneous breaking of SO(5) → SO(4) = SU (2) L × SU (2) R can be parametrized as where Tâ are the (four) broken generators in SO(5)/SO (4) given in Eq (A.1) and f is the scale of spontaneous symmetry breaking. The hâ transform as a 4-plet of SO (4), and can be arranged in a doublet of SU (2) L as: One can assume that EWSB proceeds through a non-vanishing vev h 0 = h 4 , with the vev's of the other components vanishing. In unitary gauge, h 1 = h 2 = h 3 = 0 and h 4 = h 0 + h, where h is the physical Higgs boson. Using the explicit form of the broken generators results in the matrix given in Eq. (2.1). It is also easy to embed the various fermion SO(4) multiplets used in the main text into (incomplete) representations of SO(5) that simplify the writing of the Lagrangian. For example, • For the MCHM 5 : The elementary fermions are written as , and the composite fermions are written as These result in Eqs. (2.8) and (2.9).