Multi-Representation Dynamics of SU(4) Composite Higgs Models: Chiral Limit and Spectral Reconstructions

We present a lattice study of the $SU(4)$ gauge theory with two Dirac fermions in the fundamental and two in the two-index antisymmetric representation, a model close to a theory of partial compositeness. Focus of this work are the methodologies behind the computation of the spectrum and the extrapolation of the chiral point for a theory with matter in multiple representations. While being still technical, this study provides important steps towards a non-perturbative understanding of the spectrum of theories of partial compositeness, which present a richer dynamics compared to single-representation theories. The multi-representation features are studied first in perturbation theory, and then non-perturbatively by adopting a dual outlook on lattice data through a joint analysis of time-momentum correlation functions and smeared spectral densities.


Introduction
The Standard Model (SM) of particle physics describes the strong and electro-weak (EW) interactions with remarkable accuracy, and no clear deviation from its predictions has been observed.There are however open problems pointing towards the idea that the SM effectively describes Nature only up to a cutoff energy scale Λ SM located at least in the TeV range.One of these issues, known as "Naturalness problem", lies in the Higgs sector, where quantum corrections are expected to push the mass of the Higgs boson towards Λ SM .The experimental value [2,3], however, lies well within the EW range.This value could only be understood within the SM by relying on fine-tuned cancellations, which are considered unsatisfactory from a theoretical perspective.The solutions that have been proposed to tackle this fine-tuning issue have been generating a vast and fertile literature, including ideas such as supersymmetry and technicolor.Another popular solution, and focus of this work, is the "composite Higgs" scenario [4] where the Higgs boson's mass is explained in terms of Goldstone dynamics.A gauge theory is postulated to describe a new strongly-interacting sector and its dynamics.Depending on the fermionic content of the theory, a global flavor symmetry is realised: if the fermions condensate, the spontaneous breaking of such symmetry generates Goldstone bosons, among which there would be the Higgs doublet, the lightest state of the new sector.
A number of models are compatible with the composite scenario.A set of minimal candidate gauge theories have been identified in Ref. [5].Such theories exhibit a doublet compatible with the SM Higgs boson in the low energy theory, together with other interesting phenomenological features such as asymptotic freedom and the presence of a custodial symmetry.Remarkably, these theories yield a composite state with the same SM quantum numbers as the heavy quarks.This crucial property could clarify the hierarchical structure of the quark masses: if the composite partner of the top quark has a large enough anomalous dimension, the mass hierarchy arises naturally.This idea goes under the name of partial compositeness [6] and it has been the subject of several lattice studies in recent years [7][8][9][10][11][12].The phenomenology of theories of partial compositeness is rich and can differ from Quantum Chromodynamics (QCD) in many aspects.It is common for these theories to involve fermions transforming in multiple representations of the gauge group.While QCD and its fermions in the fundamental representation have been extensively studied for several decades leading to tremendous improvement in the field, theories with multiple representations are still at an earlier stage, despite recent notable progresses [7][8][9][10][11].A promising model of partial compositeness has been proposed by Ferretti [1].Its UV completion is an SU (4) gauge theory featuring three Dirac fermions in the fundamental (Fund) and antifundamental representation, and five Majorana fermions in the two-index antisymmetric (2AS) representation of the gauge group.In this work, following the work presented in Refs.[7][8][9][10], we perform a lattice study of a simplified version of the Ferretti model, containing two Dirac fermions in the fundamental and antifundamental and two Dirac fermions in the AS representation of the gauge group SU (4).This model does not form a bound state compatible with the Higgs boson in the low-energy limit, but it is expected to maintain some of its non-perturbative features, and it therefore represents a solid starting point towards a better understanding of this class of theories.
Our investigations focus on the methodologies required to address theories with matter in multiple representations.Although lattice simulations provide a flexible framework which has been well established over many decades, theories of partial compositeness present challenging features, and a dedicated study is an important step in order to ultimately make contact with their phenomenology.We will show that a richer dynamics, a distinguishing mark of these theories, complicates tasks such as the computation of the spectrum or the extrapolation to the chiral limit.For the model under consideration, we have used numerical simulations and perturbative calculations in order to clarify such dynamics.In particular, we have generalised previous perturbative results to the case of multirepresentation theories in order to predict the critical mass of Wilson fermions.We have generated gauge configurations at different fermionic masses, enabling the analysis of the theory at the chiral point.For the first time, we performed a comparison of mesonic masses estimated from correlation functions and spectral densities: these two quantities, related by a Laplace transform, provide a dual outlook on non-perturbative data obtained on the lattice.The computation of spectral densities from lattice correlators has been the subject of several studies leading to interesting progress in the field [13][14][15][16][17][18][19] and it is here performed for the first time in the context of BSM physics, where the lack of phenomenological inputs and a more sophisticated dynamics call for new computational strategies.
The rest of this work is structured as follows: Sec. 2 examines the model introduced by Ferretti and its simplified version with less matter content.Sec. 3 contains perturbative results for the critical mass of Wilson fermion, with a focus on the multi-representation dynamics.The lattice setup is outlined in Sec. 4, and the observables targeted in our numerical simulations are described in Sec. 5.In Sec. 6 we present results for the extrapolations of our data to the chiral point.Sec.7.1 is dedicated to the extraction of spectral densities from lattice correlators.Sec.7.2 contains a discussion on features of the spectrum of multi-representation theories.Sec.7.3 concludes the discussion with a description of methodology and results of fits of spectral densities, which are compared with the study of lattice correlators in the time-momentum representation.

The Ferretti model
We briefly recall the features of the model introduced by Ferretti in Ref. [1] following the conventions of [10].Its UV completion is described by the gauge group G = SU (4).The gauge field is coupled to three Dirac fermions that we express in terms of Weyl doublets χ a m , χa m respectively in the fundamental and antifundamental representation of the gauge group, together with five Majorana fermions ψ I mn in the 2AS representation, which is real and dimension 6.The indices a = 1, 2, 3 and I = 1, . . . 5 are flavor indices, while m, n = 1, 2, 3, 4 denote the color.This matter content induces a global symmetry described by the group The charges of the fermions with respect to the flavor group are described in Table 1.
This symmetry breaking pattern is interesting for several reasons.Given that SO(5) ⊃ SU (2)×SU (2), the pattern is compatible with the requirement of custodial symmetry The unbroken group SU (3) c , related to the fundamental sector, is responsible for the strong interaction of QCD once it is gauged.The unbroken group related to the 2AS fermions, SO(5), contains the EW group SU (2) L × U (1) Y .In fact, SO(5) ⊃ SO(4) SU (2) L × SU (2) R .We then define an U (1) R as the subgroup of SU (2) R generated by the generator of isospin rotation T R : the correct hypercharges Y are then obtained by Y = T (3) R + X, X being the charge under U (1) X .The quotient SU (5)/SO (5) is therefore the relevant one for EW symmetry breaking: by writing its 14 Goldstone bosons in terms of SM charges, 14 we identify an SU (2) doublet 2 ±1/2 , H, that is compatible with the Higgs boson.
Turning to the composite partner for the top quark, this is introduced as a Dirac fermion Ψ [1] in the low energy theory that has charges (5, 3) 2/3 with respect to the unbroken subgroup H. States with these quantum numbers, relevant for partial compositeness, are obtained in this theory by color singlet combinations of fermions in different representations.Such baryonic content is therefore typical of theories with multi-representation matter.

Two-flavor Ferretti model
In this work we will focus on a simplified version of the Ferretti model.We will consider two Dirac fermions in the fundamental and two Dirac fermions in the 2AS representation of the gauge group SU (4), a model that has been already studied in [7][8][9][10]12].While retaining the multi-representation dynamics and some non-perturbative features, this choice changes the global symmetries of the theory.
It is important to understand the discrete symmetries of each sector in order to give the correct interpretation to the lattice data.Isospin, in particular, is useful in classifying scattering processes.The isospin group in the fundamental sector is the well known SU (2).Since the symmetry breaking pattern is the same as in massless two-flavor QCD, characterised by the quotient group and does not require to be discussed here.For completeness, we only mention that the three Goldstone bosons π 1 , π 2 , π 3 arising from these cosets can be labelled with eigenvalues of the azimuthal component of the isospin, ±1, 0: The multiplet (π + , π 0 , π − ) has eigenvalue −1 under the G-parity defined by combining charge conjugation C with an SU(2)-isospin rotation, τ i being SU (2) generators.In the 2AS sector, instead, the symmetry breaking pattern yields the cosets SU (4) SO( 4) .
The generators T i of SU (4) can be found in Appendix A. These cosets are characterised by 9 Goldstone bosons Π i , i = 1, . . .9, that we represent exponentially as where Ti , i = 1, . . .9 are the broken generators of SU (4).Under an isospin transformation, where X n are the unbroken generators of SU (4) generating SO(4), represented as 9 × 9 matrices.We give such a representation in Appendix F. The maximum set of commuting generators here is two, therefore we choose to diagonalise X 1 and X 6 .The Goldstone bosons in the isospin basis can be then labelled as Π a1,a6 , where a n are eigenvalues of the generators X n From these expressions it can be shown that the operation of charge conjugation acts on this multiplet as a transformation of SO(4).As a consequence, any G-parity is equivalent to an isospin rotation and does not provide selection rules for transition amplitudes.Implications of this feature will be discussed in Sec.7.2.

Perturbative Results
In this work we adopt a Wilson-type discretisation of the Dirac action with a clover term, the details of which will be given in Sec. 4. This choice for the action breaks chiral symmetry.The critical mass, i.e. the value of the bare mass of a fermion corresponding to a vanishing renormalised mass can be first estimated in perturbation theory, suggesting values for the numerical simulations, and providing first insights on the multi-representation dynamics.The perturbative expansion of the Wilson action generates, for each representation, the same vertices that appear in lattice QCD up to group theoretical factors, making it easy to generalise the existing result.To this end, we will refer to the calculation of the critical mass of Wilson fermions at two loops [20] in lattice QCD, the cactus resummation at one loop [21] and their generalisation to a generic representation of SU (N ) [22].These results will be extended to the case of multiple representations in the reminder of this section.

Multiple representations
In this section we analyse the effect of multiple representations in the computation of the fermionic self-energy.The motivation is to gain insights about the way a representation can affect the other.For simplicity we use, for this task, Wilson-type fermions, delegating the discussion of the clover term and cactus improvement to the next section where we estimate critical masses.For a given representation, we write the perturbative expansion of the one particle irreducible two-point function as The one-loop contribution Σ 1 , for instance, takes contributions from two diagrams, a tadpole and a sunset.It is useful to parametrise Σ 1 in terms of powers of the lattice spacing, Imposing the vanishing of the renormalised mass for the fermion yields an expression for the critical mass at one loop The same can be done order by order in powers of the coupling g 2 0 , which is related to the usual lattice coupling by The result in a given representation R can be written in terms of the contribution of each diagram [20] Σ (1) (R) = 2C 2 (R) c (1) , c (1) This result, plugged into Eq.( 13), provides the one-loop estimate of the critical mass.However, it contains no information about the multi-representation dynamics, for which we need at least the O(g 4 0 ) result.In Wilson lattice QCD, Σ (2) takes contribution from 26 diagram [20].If the theory has fermions in two different representation R and R , there are four additional diagrams contributing to the self-energy of R due to loops of R .These are shown in Fig. 1 and their contribution is the same as their single-representation counterparts [20] up to some group theory factors.The diagrams I 1 and Figure 1: Two-loop diagrams in which the representation R (solid line) can contribute to the self energy, and therefore the critical mass, of the fermion in the representation R (double solid line).These diagrams, apart from group theoretical invariants, depend on factors that have been computed in Ref. [20].
I 2 in Fig. 1 must be evaluated together in order to be infrared finite.Their value is where c (2) 1 = 0.00079263(8) [21] is representation independent, and n f is the number of fermions.Similarly, I 3 + I 4 gives the infrared-finite result where c (2) 2 = 0.000393556 (7).The two-loop part of the self energy for the representation R due to the presence of R is then We then add to this contribution the 26 single-representation diagrams Σ (2) with The total two-loop self energy is given by the sum of Eqs. ( 18) and ( 19) Unsurprisingly if R = R the extra term is equivalent to the addition of extra flavor content.
The multi-rep contributions alone are small compared to the rest of the terms, providing about 10 − 13% of the two-loop part.Although these results are only perturbative, they will find a nonperturbative counterpart in Sec. 6.

Critical mass
A perturbative prediction that compares to our numerical results can be obtained by including the clover term into the action, resulting into an additional interaction vertex.The prediction can be further improved by resumming an infinite series of a specific type of gauge invariant (cactus) diagrams [21].For this analysis, we consider one-loop results.
The critical mass m c (R) of a Wilson fermion can be computed from Eq. (15).Considering the clover term yields new contributions in powers of c sw [23], where the coefficients i can be read from Table 2 of Ref. [23].Concerning the Wilson term, the cactus resummation is performed by rescaling the self-energy Σ → Σ/ c0 , where c0 is a function of N and g 2 0 but not of the representation.c0 can be found by solving the following equation, where L α N are generalised Laguerre polynomials of degree N .The solution for β = 11 and N = 2, 3, 4 is shown in Fig. 2 as the intersection between the various curves and the vertical line.For SU (4) we find c0 = 0.731607 .
The resummation with the clover improvement term corresponds to also rescale g 2 0 → g 2 0 /c 0 and c sw → c sw c0 [24].The prediction for the critical mass is then The numerical values for SU (4) with two fundamental and two antisymmetric fermions are listed in Table 2.We find the one-loop result with cactus resummation to be the closest one to non-perturbative results. W.
W These values are to be compared with the non-perturbative results of Sec. 6.

Lattice Setup
The lattice setup has already been discussed in Ref. [10], and here we only recall the main ideas.We generate gauge configurations by using the HMC algorithm [25] with a multi-level second order integrator [26].The lattice action can be decomposed into gauge and fermionic part.The latter is further decomposed into contributions from each representation [27] For the gauge part we use the Wilson plaquette action The fermionic action for the representation R is We adopt a Wilson discretisation for the Dirac action D R in each representation R, plus a clover improvement term The Wilson term in position space is where and κ R is related to the bare mass m R 0 of the fermion in the representation The order a clover improvement term is where Q µν (x) being the clover combination of plaquettes around the point x [28].The improvement coefficient c R sw can be computed in powers of the coupling, A discussion on the O(g 2 0 ) coefficient can be found in [29].In this work, we set the coefficient to its tree level value, c R sw = 1, for both representations.This choice for the discretisation action breaks explicitly chiral symmetry, resulting in an additive term to the renormalisation of the fermion masses.Simulations at the chiral point are performed indirectly by extrapolating the value of the critical masses m R c , i.e. the value of the bare masses at which each fermion has a renormalised vanishing mass Direct simulations at the chiral point are in fact impossible due to the spectral properties of the Wilson operator.The lowest eigenvalue of the Wilson operator approaches zero in the chiral limit, meaning its inverse becomes increasingly ill-conditioned.When extrapolating to the critical point (see Eq. ( 39)) one needs to care that exceptional configurations do not occur in the gauge average.These are configurations especially close to the chiral point that can jeopardise the inversion of the fermionic operator.For this reason, we monitor the gauge distribution of the lowest eigenvalues of the Wilson-clover operator [30] as in Fig. 3, making sure that they remain sufficiently far from the origin.In order to parametrise the breaking of chiral symmetry, we compute the PCAC mass, which yields a definition of the quark mass through the axial Ward identity.Our starting point is the work done in [10], where the space of bare parameters for this model has been explored.The lightest ensemble of that work is present in this analysis under the name A0.
Details about the ensembles generated in this work are found in Appendix D. We have performed simulations for increasingly lighter masses, allowing the extrapolation of chiral points for both the fermionic representations.We generate gauge configurations at single lattice spacing, using β = 11, and at a single volume corresponding to a lattice of dimensions L 3 × T = 16 3 × 32.For the production of gauge configurations and the measurements of observables, we use the software Grid [31] and Hadrons [32].

Observables
In this section, we give a description of the observables that are targeted in this work.These include two point functions that will enable us to compute parameters for the chiral symmetry breaking, mesonic masses and smeared spectral densities.We compute correlation functions of the following interpolators which are defined for each representation In this work we only consider isospin-vector operators, which are used to build two-point functions The two-point functions encode information about finite volume matrix elements and energy levels.This can be seen by expanding the previous expression on a complete set of states: While we will assume, to slim the notation, a positive sign between the exponentials in the previous expression, it is understood that this may vary depending on the quantum numbers a, b.The energy of the ground state with given quantum numbers can be obtained by the asymptotic behaviour in the Euclidean time of the corresponding correlator.For example, the effective mass of a pseudoscalar meson can be computed from the pseudoscalar-pseudoscalar correlator Similarly the PCAC mass, defined through the axial Ward identity, is obtained from the pseudoscalar and axial correlators which has O(a) effects for our choice of the unimproved pseudo-axial operator.
In the chiral limit, the pseudoscalar mass aM R PP for each representation vanishes.In our ensembles we generated configurations describing mesons of antisymmetric fermions being generally heavier: this Figure 4: Probability densities for the lowest eigenvalue of the fermionic operator (left), the PCAC mass (centre) and the pseudoscalar mass (right) for both representations.These values are especially interesting since they are based on the ensemble S0, where we find compatible values for the lowest eigenvalues and, within 1σ, the PCAC masses.Nonetheless, a difference arises in the pseudoscalar masses.The distributions of the PCAC and the pseudoscalar masses are obtained from a resampled set of configurations.can be understood from Fig. 4 where we show the gauge distributions of the lowest eigenvalue of the fermionic operator, the fermionic masses and the pseudoscalar masses for both representations from the ensemble S0.In this ensemble, the lowest eigenvalues of the Wilson operator are the same, and the fermionic masses are also compatible within one σ.Nonetheless, an unambiguous gap appears in the masses of the mesons.The Gell-Mann-Oakes-Renner relation predicts this behaviour to be the result of different chiral condensates and pseudoscalar decay constants for the two representations.
Ground state energies can be estimated from lattice correlators according to Eq. (43).Other energy levels and matrix elements can be estimated by fitting sums of exponentials, as in Eq. ( 42).The extraction of excited states from these fits is in general hindered by the increasing number of degrees of freedom that are needed in order to perform the fit when more excited states are being targeted, with a given number of data points.This becomes particularly problematic when dealing with highly correlated data, which limits the effective number of degrees of freedom.Moreover, as the infinite volume limit is approached, the spectrum above the multi-particle threshold becomes continuous and resolving energy levels becomes exponentially harder.For these reasons, it is desirable to have many correlators in order to perform simultaneous fits.Alternatively, variational methods such as the generalised eigenvalue problem (GEVP) are a well established way to obtain finite volume spectra [33].
Other observables that allow the extraction of finite volume quantities are spectral densities.These are related to lattice correlators by a Laplace transform Spectral densities contain the same information as lattice correlators, with the difference that for the spectral densities the information is encoded in a function of the energy rather than Euclidean time.
A finite volume spectral density ρ L ab (E) can in fact be expanded as The distributional nature of the spectral density occurs due to the finite volume of the simulation and it is mostly lost in the continuum limit, where above the multi-particle threshold the spectral density becomes continuous, while it retains Dirac deltas in presence of single-particle or bound states.In order to treat this distributional character we smear the finite volume spectral density [13] with a Gaussian function, so that the smeared spectral density ρ L, R σ,ab (E) is a continuous function even at finite L. In this work, we will focus on the extraction of the finite volume energies and matrix elements from smeared spectral densities.This task does not require us to take the infinite volume limit, neither to remove the regularisation by extrapolating at zero smearing radius σ.By rewriting Eq. ( 46) for the smeared spectral density, it is clear that this function can be fitted against sums of Gaussians in order to obtain the energies and the overlaps, similarly to how Eq. ( 42) is commonly used to extract the same quantities by fitting sums of exponentials.This offers a dual picture with quite different features that will be described in the next sections.

Chiral Limit
At the chiral point, both pseudoscalar mesons become massless.The discretisation of the lattice action that we adopt breaks chiral symmetry explicitly.We therefore compute the PCAC masses from the axial Ward identity in order to quantify the breaking of chiral symmetry, and we scan the space of bare parameters until we are able to locate the bare masses for which the PCAC masses vanish.At the chiral point we also expect the Gell-Man-Oakes-Renner relation to be valid, with the pseudoscalar mass squared (M R PP ) 2 scaling as m R PCAC for the representation R. The correlation functions of fermions in the fundamental representation that we use to extract the masses are more affected by autocorrelation, which translates into usually larger statistical errors for the masses of the mesons related to that representation.The measurements of the PCAC masses are listed in Tables 6  and 7 of Appendix E. In order to perform a linear extrapolation for the vanishing of the PCAC mass in a given representation, the bare mass of the other must be kept fixed.Fig. 5 shows three of these extrapolations.In the left panel, we extrapolate the chiral point for the fundamental fermions by fixing am (49) Figure 5: Three chiral extrapolations, one on the left for fundamental fermions and two on the right for antisymmetric fermions.The bands correspond to linear fits.The heavier point in the plot on the left was also present in [10].
Comparing these values to the perturbative results in Table 2, we see that the prediction improved by resumming cactus diagrams is the closest to the non-perturbative values.We can also observe from Tables 6 and 7 that taking a fermion in the representation R towards the chiral point pushes the fermion in the representation R to be also lighter.Figure 6 shows the dependence of the pseudoscalar masses M R PP with respect to the quark masses.For both representations R the scaling is compatible with the one predicted by the chiral Lagrangians, (M R PP ) 2 ∼ m R PCAC , with the pseudoscalar mesons becoming massless in the chiral limit.The values of the masses are reported in Table 8.

Smeared Spectral Densities from Lattice Correlators
The non-perturbative calculation of spectral densities has been receiving increasing attention [13][14][15][16][17][18][19].In this work we will analyse spectral densities obtained from lattice correlators, for the first time to our knowledge, in the context of BSM and multi-representation theories.In order to facilitate the discussion, it is useful to recall the computational details of the calculation, aiming to a self-contained discussion.

The numerical procedure
Since we are interested in finite volume quantities we will omit, in this section, the dependence on the spatial volume L that we will assume to be finite.Computing spectral densities ρ R ab (E) from Euclidean correlators C R ab (t) involves an inverse Laplace transform, where E min can range between zero and the energy of the ground state, since ρ R ab vanishes in that interval.The inversion of Eq. ( 50) is numerically an ill-posed problem which needs to be regularised.Recalling from Sec. 5 that we are interested in the smeared version of ρ R ab (E), it is especially convenient to approach the inverse problem using the Backus-Gilbert type regularisation [34] introduced in [14], which yields spectral densities smeared with a chosen smearing function f (E), An important observation is that a fixed smearing kernel is crucial in order to perform fits and extrapolations of the results.The idea of the algorithm is to generate the target smearing kernel, the Gaussian ∆ σ (E) of Eq. (47), as a linear combination of the same exponentials appearing in the Laplace transform of Eq. ( 50), g τ e −τ aE + e (−T +τ )aE , where t = τ a, a being the lattice spacing.Once the coefficients g τ ≡ g τ (σ, E ) are known one can simply obtain an estimator for the smeared spectral density, On the lattice the correlators C R ab are available for a finite number of times, therefore the sum in Eqs. ( 52) and (53) must be truncated at the appropriate cutoff τ max .Since our lattices have temporal length T and periodic boundaries, we have aτ max = T /2.The reconstructed smearing kernel, will necessarily differ from the Gaussian ∆ σ (E) at finite τ max , inducing a systematic error on the final result.The computation of the coefficients g is achieved through the minimisation of the functional where λ is a trade-off input parameter that we will discuss later in this section.The functional A α [g], introduced in [14], measures the difference between the exact smearing kernel and the one we can reconstruct with the available data The parameter α < 2 enables the selection between a class of norms in order to measure the distance between the target and the exact function.Choosing larger values of α allows for the integrand to decay faster at high energies.The functional B[g] is needed to regularise the problem [34], making it numerically stable.We define B[g] to be dimensionless, namely where Cov is the covariance matrix of the correlator C(aτ ) estimated over N bins The algorithmic parameters can be gathered to simplify the notation The minimisation of W α [g] corresponds to solving the following linear problem, which has to be performed at each energy and smearing radius for which we want to estimate ρ R ab, σ (E ).The nature of the functionals appearing in these definitions is intimately related to the uncertainties on the estimator ρ R ab, σ (E ).The statistical error is in fact estimated by The systematic error, unavoidable at finite τ max , is estimated by monitoring the quantity as we now describe.Regions where d(g p ) is small are dominated by the statistical uncertainty; on the other hand, the reconstructed smearing function of Eq. ( 52) will be as close as possible to the exact one.For these reasons, it is not surprising that where d(g p ) is small, the results for ρ R ab, σ are stable2 in response to variations of the algorithmical, unphysical parameters p, within statistical error.If we define the coefficients g as we find, with the given normalisations of the functionals and with the quality of our data, that in the range 0.1 < k < 5 and α = 1 the outcome of the reconstruction is in a region of algorithmical stability: d(g p ) is small, and systematic fluctuations are well within the statistical uncertainty.When this is not realised, the systematic error is estimated as where g is defined through Eq. (63) at a different value of k that allows us to account for the systematic fluctuations.In this work, we find that a definition of g at k/10 provides a conservative estimate of ∆ sys (E ) for a reconstruction performed with g at the value k.Fig. 7 shows an example, for a specific energy E , of a stability regime where the spectral reconstruction for different values of λ (black points) does not fluctuate outside the statistical error (black bars).Three points on the plot are highlighted: the ones corresponding to the choices of k = 0.1, k = 2.5 and k = 1 in Eq. ( 63), showing that indeed no systematic component affects the uncertainty on the result.The value corresponding to k = 1 is associated to the value of λ at which one achieves the optimal balance in agreement with the prescription from Refs.[14,16].

Excited states in the antisymmetric sector
The excited states created by hadronic interpolators have a big impact on the extraction of the effective masses, since it can be hard to distinguish them from the ground state.In theories with multiple representations (or even just more flavors of a single representation) this effect is amplified when all the fermions are approaching the chiral limit.An interpolator can now create states containing particles from both representations.Since we simulate lighter fundamental fermions, we can expect this feature to be more visible in the 2AS sectors rather the in the fundamental one.Moreover, the 2AS sector does not have G-parity selection rules preventing certain states to mix, and we thus expect this channel to have a richer dynamics.In order to go into more details, it is useful to set the notation at least for the pseudoscalar interpolators and the correlation functions The importance of both identifying and controlling the excited states can be appreciated by looking at the smeared spectral density, for instance, in the pseudoscalar channel The magnitude of each matrix element n| ŌR P (0)|0 L determines the weight of each Gaussian ∆ σ located at the energy E n (L).If these energies are too close, or the matrix elements are too large, resolving different states can become laborious.With this motivation, in order to control the excited states we build correlation functions including two different types of fermionic field: local, pointlike operators and Gaussian-smeared ones 3 .Operator smearing allows working with interpolating operators that have a weaker overlap with excited states.Details concerning the measurements of the correlation functions and operator smearing can be found in Appendix C.
The states created by each operator can be identified, as we have seen, by means of its symmetries.In the fundamental sector, at zero angular momentum, a pseudoscalar meson can induce the following transitions from the vacuum 0|π|π , 0|π|πππ , 0|π|πΠΠ , . . .
that will enter in our analysis through Eq. ( 67).Since in our simulations PP this phenomenology is reminiscent, up to E 3π , of QCD, and one can expect computational aspects to be also similar.Conversely, due to the triviality of its G-parity, the 2AS sector has a multi-particle threshold located at E 2Π .In addition, since the other representation has lighter particles, states containing pseudoscalar mesons from the fundamental sector are not guaranteed to have energies far from the ground state.Possible overlaps with a pseudoscalar meson are in fact 0|Π|Π , 0|Π|Π Π , 0|Π|Πππ , 0|Π|Πππππ , . . .
The aforementioned features complicate the extraction of M (2AS) PP , as it can be understood from Fig. 8 where we show two different types of signal from a correlator built with local, unsmeared operators.The excited states contaminate the signal for the ground state resulting, in the left panel, in an effective mass that does not reach a clear plateau.The problem is also manifest in the energy picture, as it is shown in the right panel of the same figure, where the smeared spectral density ρ(2AS) PP, σ does not exhibit the expected Gaussian peak around the mass of the pseudoscalar meson, but it rather grows monotonically.Indeed, by decreasing the smearing radius σ of Eq. (47) one should be able to resolve such peak, but this cannot be realised with the current quality of the data.While the temporal length of the lattice poses an intrinsic obstacle to the thermalisation of the effective mass, the spectral reconstruction in principle allows obtaining smaller smearing radii also by increasing the number of configurations, since the systematic and the statistical error are related by Eq. (55).: Results from a two point function of pseudoscalar operators built with point-like antisymmetric fermionic fields.The correlator is estimated from the ensemble B1.The left panel exhibits the effective mass as a function of time.Due to the nature of the excited states in the 2AS sector, the mass does not reach a plateau in the available time.The dominance of the excited states can also be understood from the smeared spectral density in the right panel.The overlap between the interpolator and the excited states it creates is too large: the spectral density smeared according to Eq. ( 47) is dominated by contributions above the multi-particle threshold, preventing the identification of the ground state.
Having established that the excited states present a challenge in the 2AS sector, it is natural to look at correlation functions of smeared operators defined in Appendix C, which have suppressed overlap with the excited states.Fig. 9 shows the effective mass and the spectral reconstruction from the correlation function of such operators.On the left panel, the plateau in the effective mass shows an improvement compared to the corresponding result in Fig. 8.The effective mass is independent of time, within its statistical errors, for t/a > 10.The smeared spectral density, shown in the right panel, demonstrates again the suppression of the excited states, with contributions from higher-energy states becoming smaller.As a result, a single peak is clearly visible at E 2M , which cannot be resolved because the smearing radius is too large, σ M (2AS) PP .These energies can be nonetheless estimated by fitting the spectral density to a sum of Gaussians.This idea will be expanded in the next section.47).The peak includes contributions from mainly M Π and E ΠΠ .

Fits of spectral densities
In this section we describe fit strategies for spectral densities.A parallel discussion on fits of correlators will highlight the differences between the two methodologies and will lead to a quantitative comparison between predictions for the pseudoscalar masses obtained from the two approaches, which is presented in detail at the end of the section.The model functions used in the fits are c (k) (t) for correlators and f (k) σ (t) for the smeared spectral densities where σ is the smearing radius defined by Eq. ( 47).The integer k encodes how many states are included in our model function.E n , a n and b n are the fit parameters which relate to finite volume energies and matrix elements.These are estimated by minimising appropriate χ 2 functions where covariance matrices are estimated as in Eq. ( 58) both for correlators and spectral densities.(t).The points at which the spectral density is evaluated are chosen in order to minimise the condition number of its covariance matrix.Due to this freedom, we obtain a matrix for the spectral density that is better conditioned than the one for the correlator.
On the lattice, the temporal length T constrains the maximum number of data points and hence degrees of freedom for fitting a correlator C R PP (t).The effective number of degrees of freedom may be further constrained in the presence of correlated data.The effect of correlation between times t and t is taken into account by the covariance matrix appearing in χ 2 g (k) .A smeared spectral density ρ R PP, σ (E) can be in principle evaluated for any energy E, yet being it derived from the correlator itself, it seems unrealistic to think it could lead to a larger number of degrees of freedom to be used in a fit.The information, however, is mixed non-trivially: the spectral density at a single point depends on the correlator at each lattice time t.A consequence of this feature is manifested in the extraction of the ground state, where the problem of having uncontaminated signal at large times is shifted to having a small enough resolution σ in the energy.
With this motivation, it is interesting to study the number of degrees of freedom we can exploit in each correlated fit.While the correlators can only be evaluated at integers 0 ≤ t < T , we have freedom to choose the energies in a given interval at which we evaluate the spectral densities.Our criterion is to select those that minimise the condition number of Cov[ρ σ ] in order to maximise the information that is passed to the χ 2 .We observe that correlated data in time does not necessarily translate into correlated data in energy space, as shown in Fig. 10, where we compare covariance matrices for C (2AS) PP (t) and ρ (2AS) PP, σ (E) from the ensemble B3 with a smearing radius σ = 0.2/a.In both cases, the degrees of freedom are enough to fit at least two states.The matrix Cov −1 tt [C], in particular, is evaluated from correlators measured from t/a = 8 to t/a = 16 with an interval of two.Regardless of the thinning in time, adjacent points shows a very high correlation.If the covariance matrix of the correlator Cov[C] is ill-conditioned, it needs to be regularised by applying a cutoff on its smaller eigenvalues before χ 2 g k is evaluated, a problem that is not faced when fitting spectral densities, whose covariance matrix is easier to invert.Indeed, one has to invert the covariance of the correlators in order to compute the spectral density, but only in the combination defined by Eq. ( 55): the matrices obtained from the functionals A α [g] and B[g], if both ill conditioned, regularise each other for suitable values of λ.The choice of a cutoff for the covariance matrix Cov[C] is therefore absorbed into the choice of the parameter λ.
Fig. 11 shows two examples of correlated fits of the smeared spectral densities.On the left panel, the correlator used to extract ρ (2AS) PP, σ is built with smeared interpolating fields.The model function is Figure 11: Examples of fits of spectral densities, showing the breakdown of the contribution of each Gaussian.On the left-panel, a three Gaussian fit of a smeared spectral density extracted from a correlator that uses smeared, non-local fields.Due to this choice for the interpolator, the Gaussians in the plot contribute less as we go higher in the energy range.On the right-panel, we show a similar plot obtained from local interpolators: the two Gaussian fit is still able to isolate the ground state even if the effective mass does not plateau.Both plots correspond to the pseudoscalar 2AS channel.The left panel is obtained within the ensemble B3 and the spectral density has smearing radius σ = 0.24/a.The right panel is derived from the ensemble B2 and the smearing radius of the spectral density is σ from Eq. (70).The plot shows the explicit breakdown of each Gaussian, which correspond to the contribution of different energy levels: as each of them becomes decreasingly important, we clearly see the excited state suppression achieved by the choice of smeared operators in the correlation function.At low energies, the spectral density is almost entirely dominated by the first energy level, therefore the corresponding fit parameter is mainly constrained by energies near the origin.Cancellations between the three correlated contributions combine in an error on the fit result that is generally smaller than the one on the single Gaussians.We identify the first peak as the value of aM (2AS) PP and the second with aE ΠΠ .The right panel of Fig. 11 is instead obtained from local interpolators, which have a larger overlap with excited states.As shown in Fig. 8, the effective mass of Eq. (43) does not reach a clear plateau in this case, yet the fit of the spectral density is able to isolate the ground state.The error on the fit, however, is at best one order of magnitude larger than the corresponding one obtained in the left panel.The choice of smeared operators is therefore preferred.
We now turn to the comparison of the fit results.As described in Sec.7.2, the extraction of the ground state presents challenges in the antisymmetric sector, which therefore provides an interesting testbed to compare the two frameworks.We use for the comparison the ensembles B1-B4, where both representations are fairly light.We begin by discussing the estimate of aM (2AS) PP from lattice correlators.In some cases, the covariance matrix of the two point function had to be regularised by introducing a cutoff on its lower eigenvalues, in order to invert it in the χ 2 .The choice of this cutoff can translate into fluctuations in the estimate of the pseudoscalar mass outside the statistical error, that have been accounted for by adding a systematic component to the uncertainty.This problem does not appear while fitting spectral densities, because the covariance of the spectral density is better conditioned.We have also ensured that no contamination was present in the estimate for the pseudoscalar mass by comparing fits of one and more exponentials.The values of the pseudoscalar mass obtained with this approach can be found in Table 8.In this framework the smearing of the  62)).The blue band is the fit result of the smeared spectral density to a sum of Gaussians at the point E .Encouragingly, the fit is compatible with all points in the scan, showing that stability in the reconstruction translates into stability for the fits.
interpolators has been crucial, as it is clear from Fig. 8 where the effective mass of Eq. ( 43) does not plateau due to the short temporal extent of the lattice.The spectral density, on the other hand, does not rely on any large time behaviour.Its limit lies in the high energy range, where its error becomes large.In order to perform the reconstruction, the algorithmic inputs are chosen in the region of stability as described in Sec. 7. Fig. 12 updates the plot from Fig. 7 showing that stability for the reconstruction translates into stability for the fits: the blue band, which is the fit result for the smeared spectral density at the energy E , is compatible with all the values generated by different choices for the unphysical parameter λ.
The choice of the smearing radius is dictated by the quality of the data.In this work, we managed to obtain values ranging from σ = 0.18/a to 0.33/a.Since the separation between finite volume energies should be roughly 2π/L 0.4, these values can be considered acceptable.For each ensemble, we have performed the fit at a fixed value of σ, obtaining a prediction for the pseudoscalar mass aM (2AS) PP, σ .We then performed a scan over different smearing radii to check for systematic effects.As shown in the example of Fig. 13, the smearing radii adopted were found to be small enough to identify the ground state; in most cases, still, we have observed fluctuations for aM (2AS) PP, σ as σ varies.When they occurred, we added half the spread of these fluctuations as a systematic error.Fig. 13 also shows a comparison between fits that include two and three states: the results are in good agreement, signalling that no contamination from excited states is affecting the estimate of the mass.
The comparison between fits of spectral densities and correlators is shown in Fig. 14.The predictions are always compatible, and the errors lie in the same order of magnitude, but the uncertainty on the correlator is generally smaller, up to a factor of approximately two.
The numerical values used in this comparison are listed in Table 3.It should be noticed that the outcome of this comparison holds with the given amount of statistics and time extent of the lattice.These quantities, in fact, heavily influence the analysis both on the side of the correlators and the spectral density, yet the way in which the two methodologies are affected can be different.

Conclusions
The model studied in this work allows us to understand the dynamics of partial compositeness better, by extending the effort started in Ref. [10] at a reasonable computational cost.In order to make contact with the phenomenology of these theories from a lattice perspective, the systematics related to the computation of the spectrum and to the extrapolation to the chiral limit need to be under control.Our analysis has been developed in this context.We have explored the perturbative structure of the theory, extending previous computations of the critical mass [20][21][22][23][24] to the case of multiple fermionic representations.We have computed the self-energy of a fermion in a given representation, focusing on the dynamical effects due to the presence of more fermionic fields in a different representation of the gauge group.While leaving a clear footprint, the second representation has only a minor numerical impact on the value of the critical mass.This result, while perturbative, found a counterpart at the non-perturbative level in the outcome of our simulations, as it is shown in Tables 6-8.By approaching the chiral limit in a given representation, in fact, we found a weak dependence on the bare parameters of the other one.The extrapolation to the chiral limit was based on several ensembles that were generated for this study.This task was hindered by a strong autocorrelation affecting observables built from the fermion fields in the fundamental representation, the lightest in terms of pseudoscalar masses.A great portion of this work focused on the extraction of the spectrum, a crucial task in order to understand if these models are realistic from a phenomenological perspective.We have shown that the 2AS sector is characterised by complicated dynamics, due to the interplay between different representations, and the lack of selection rules dictated by discrete symmetries.We have found operator smearing to be essential, in this context, to provide a reliable analysis of the lattice data.In our analysis, we have taken advantage of recent progresses [13,[15][16][17][18][19] in the numerical solution of the inverse problem that allow precise reconstructions of smeared spectral densities.Our approach, introduced in Ref. [14], yields spectral densities smeared with a chosen function and with controlled systematics.Due to these features it was possible, for the first time, to explore in this work the extraction of finite volume energies from fits of spectral densities.Our results, compared to other established approaches, provide complementary insights and fully compatible results.
Encouraged by these conclusions, we plan to explore the possibility of extending this computational setup to the study of baryonic operators, as well as to other theories of composite Higgs.Our choices, both in terms of software [32,35] and computational framework, are general and easy to adapt for the study of other theories, whether they contain different number of colors, fermionic content or gauge groups.Moreover, our work clarified technical aspects of the partial compositeness dynamics: this important step will allow moving towards studies of more phenomenological relevance.Theories of partial compositeness provide a rich set of new particles, from pseudo Goldstone bosons to heavyquark partners, that are charged under the Standard Model and could be important for direct and indirect search of new physics.The extraction of spectral densities, validated in this context by our work, can be used not only to extract energy levels but also to directly compute inclusive crosssections [13,36], a possibility that we leave for future studies.The knowledge of inclusive processes can be important for the indirect search of new physics, with particles from the new sector leaving footprints in observables precisely measured both at present and future colliders.representation R as C 2 (R): Comparing Eq. (74) and Eq.(75) one obtains dim(R) being the dimension of the representation R. With these conventions, the values of these invariants for the fundamental and 2AS representation are listed in Table 4.The explicit form of the generators that were used in this work is as follows.For the fundamental representation, the N − 1 Cartan generators are defined as , −k, 0, 0, 0, . . ., 0, 0 whereas the remaining N 2 − N non-diagonal generators are defined as for 1 ≤ i < j ≤ N .Note that, with these conventions, for N = 2 the generators in the fundamental representation reduce to 1 2 σ a (σ a denoting a Pauli matrix), while for N = 3 the generators in the fundamental representation are 1  2 λ a (where λ a denotes a Gell-Mann matrix).Given a generic element u of the SU (N ) group in the fundamental representation, the same group element in the 2AS representation is a complex-valued matrix of size (N (N − 1)/2) × (N (N − 1)/2) whose entries (which can be labelled by pairs of indices (i, j), with 1 ≤ i < j ≤ N ) are defined as Interpreting Eq. (79) as a map R from the group elements in the fundamental representation to group elements in the two-index antisymmetric representation, the generators in the latter representation can then be obtained from the pushforward R of the generators in the fundamental representation.

B Autocorrelation Times
In an HMC algorithm, gauge configurations are generated through a Markov chain process.As a consequence, subsequent configurations can be correlated.Accounting for autocorrelation is essential in the estimation of observables through the Monte Carlo average.Consider a succession a i of measurements of an observable A. We estimate the integrated autocorrelation time τ int of the observable A as where t is the Monte Carlo time, W is the summation window, and Γ(t) is our approximation of the autocorrelation function of the series: a − and a + are the average over N − t measurements, respectively the first and the last ones.In order to choose a summation window, we check at which Monte Carlo time t the function Γ(t)/Γ(0) is smaller than twice its variance [37], which is estimated within the Madras Sokal approximation [38].Typical plots that we obtain for these quantities are shown in Fig. 15 for the pseudoscalar-pseudoscalar  82) for an increasing summation window W .In both panels, the correlation functions are computed for operators at different lattice times in order to monitor both the short and long distance behaviour.By knowing the correlation of an observable between trajectories in the HMC, we can establish how many of them have to be skipped performing the measurement.This is our strategy for the computation of the correlators on our ensembles.Alternatively, the naive estimate for the statistical error σ 0 of a given observable can be corrected to be σ 2 = 2τ int σ 2 0 .Accounting for autocorrelation has been essential in the following analysis, especially for those quantities computed in the proximity of the chiral limit, where the autocorrelation was more significant.

C Measurements of Correlation Functions
In this section we discuss the details behind the measurement of correlation functions C R ab (t) ≡ C R ab (t, 0) C R ab (t, 0) = where R denotes the representation of the fermionic fields.For each gauge configuration n = 1, . . .N cf g we average measurements performed every 4 source-times t i , resulting in the correlator C R ab, n (t).from smeared spectral densities for different smearing radii σ and different number of states n.These values appear in Fig. 13.

F Isospin Generators
The isospin group of the 2AS sector in the 2-flavor Ferretti model is SO(4).The Goldstone bosons arising from the SU (4)/SO(4) cosets transform in a 9-dimensional representation of the isospin, whose generators X n , n = 1, . . .6 are listed in this appendix.These are obtained according to the convention of [10] regarding the generators and the structure constants of SU (4).

Figure 2 :
Figure 2: The solution of Eq. (26) for N = 2, 3, 4 and β = 11, the value used in the numerical simulations, as will be described in Sec. 4.

Figure 3 :
Figure 3: Distribution of the lowest eigenvalue of the Wilson-clover operator for bare parameters am F 0 = −0.45,β = 11 and different masses for the 2AS fermions: am 2AS 0 = −0.59 on the left and am 2AS 0 = −0.60 on the right.Both distributions are far enough from the origin to ensure no exceptional configurations enter the ensemble.

=
−0.45.In the right panel, the two lines represent different extrapolations for the 2AS chiral point taken at different values of the bare fundamental mass, am (F) 0 = −0.45 and am (F) 0 = −0.47.The chiral point of antisymmetric fermions does not show a strong response to the shift in the bare mass of the fundamental fermions.This is in line with the perturbative prediction at two loops of Sec. 3, where we have found the critical mass to only mildly depend on the other representation.The chiral points for β = 11.0 =−0.47 = −0.630(4) .

Figure 6 :
Figure 6: Scaling predicted by the Gell-Man-Oakes-Renner relation, with the squared mass of the pseudoscalar Goldstone boson scaling linearly with the quark mass, here estimated from the PCAC relation.The extrapolation is compatible with the Goldstone bosons becoming massless at the chiral point.

Figure 7 :
Figure 7: Example of region of algorithmical stability at a given energy E .Different values of λ, which translate into different values of d(g p ) on the x-axis, produce predictions for the smeared spectral density ρ R ab, σ (E ) that are compatible within statistical error (black bars).In this case, R = 2AS, σ = 0.21/a and E M 2AS PP .The green and orange points correspond, according to Eq. (63), to values of k = 0.1 and k = 2.5 respectively.The red point, extended in the horizontal band, is obtained at k = 1 and it corresponds to the value of λ at which one achieves the optimal balance A α [g ]/A α [0] = B[g ].

Figure 8
Figure8: Results from a two point function of pseudoscalar operators built with point-like antisymmetric fermionic fields.The correlator is estimated from the ensemble B1.The left panel exhibits the effective mass as a function of time.Due to the nature of the excited states in the 2AS sector, the mass does not reach a plateau in the available time.The dominance of the excited states can also be understood from the smeared spectral density in the right panel.The overlap between the interpolator and the excited states it creates is too large: the spectral density smeared according to Eq. (47) is dominated by contributions above the multi-particle threshold, preventing the identification of the ground state.

.
The observed smeared spectral density is the result of two contributions coming mainly from the energy levels M

Figure 9 :
Figure 9: Results from Fig. 8, this time using Gaussian-smeared interpolators according to Appendix C.These operators are tuned to have smaller overlaps with the excited states.Consequently, the effective mass plot on the left reaches a plateau, providing an estimate for aM (2AS) PP .The right panel similarly shows how suppressed excited states allow for a clear peak to emerge in the spectral reconstruction smeared with σ = M (2AS) PP according to Eq. (47).The peak includes contributions from mainly M Π and E ΠΠ .

1 Figure 12 :
Figure12: The plot updates Fig.7, which shows that the reconstruction (red band) does not change outside the statistical error (black bars) for different choices of the unphysical parameter λ, in the given range of d(g p ) (cf.Eq. (62)).The blue band is the fit result of the smeared spectral density to a sum of Gaussians at the point E .Encouragingly, the fit is compatible with all points in the scan, showing that stability in the reconstruction translates into stability for the fits.

Figure 13 :
Figure 13: Fit results for aM (2AS) PP from the ensemble B3, obtained from two (green) and three (red) Gaussian fits of smeared spectral densities at different smearing radii σ.Fluctuation at different values of the smearing radius translate into a systematic component of the uncertainty.This is summed in quadrature to the statistical error in the gray, horizontal band, the estimate for the pseudoscalar mass aM (2AS) PP = 0.3550(31).

Figure 14 :
Figure 14: Graphical comparison between the two predictions for aM (2AS) PP for the ensembles B1-B4.

Figure 15 :
Figure 15: Left panel: autocorrelation function Γ(t)/Γ(0) compute for the two point function of two pseudoscalar mesons with fermions in the fundamental representation.Different colors represent different intervals in lattice time between the source and the sink.Right panel: integrated autocorrelation time as a function of the summation window, computed for the same correlator on the left side at several times.The two point function is obtained from the ensemble A0.

Table 1 :
Flavor charges of the fermions in the Ferretti model.Neglecting couplings with the SM fermions, G is an exact symmetry.Spontaneous symmetry breaking happens once the bilinears for both representations acquire a non-vanishing expectation value, leaving the unbroken subgroup H.The quotient group determining the low-energy dynamics is

Table 2 :
Critical masses for different lattice actions and improvements.W. is an abbreviation for Wilson fermion, W.c. for Wilson clover.When not specified, the results are obtained at 1 loop.

Table 3 :
Predictions for the pseudoscalar mass in the 2AS sector from different ensembles.The values are depicted in Fig. 14.

Table 4 :
Group-theoretical factors used in this work.

Table 5 :
Ensembles used to extrapolate the chiral limit of the SU (4) gauge theory with two fundamental and two two-index antisymmetric fermions.The coupling is β = 11, and the volume of the lattice is 16 3 × 32.

Table 9 :
Fit results for aM