Multi-representation dynamics of SU(4) composite Higgs models: chiral limit and spectral reconstructions

Abstract

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 proposed by G. Ferretti. 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.

1 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 $\varLambda _{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 $\varLambda _{SM}.$ The experimental value [1, 2], 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 [3] 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. [4]. 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 [5] and it has been the subject of several lattice studies in recent years [6,7,8,9,10,11]. 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 [6,7,8,9,10]. A promising model of partial compositeness has been proposed by Ferretti [12]. 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. [6,7,8,9], 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 multi-representation 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: Sect. 2 examines the model introduced by Ferretti and its simplified version with less matter content. Section 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 Sect. 4, and the observables targeted in our numerical simulations are described in Sect. 5. In Sect. 6 we present results for the extrapolations of our data to the chiral point. Section 7.1 is dedicated to the extraction of spectral densities from lattice correlators. Section 7.2 contains a discussion on features of the spectrum of multi-representation theories. Section 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.

2 Models

2.1 The Ferretti model

We briefly recall the features of the model introduced by Ferretti in Ref. [12] following the conventions of Ref. [9]. Its UV completion is described by the gauge group ${\mathscr {G}}=SU(4).$ The gauge field is coupled to three Dirac fermions that we express in terms of Weyl doublets $\chi ^a_m,$ ${\bar{\chi }}^a_m$ respectively in the fundamental and antifundamental representation of the gauge group, together with five Majorana fermions $\psi ^I_{mn}$ in the 2AS representation, which is real and dimension 6. The indices $a=1,2,3$ and $I=1,\dots 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 G

$$\begin{aligned} G = SU(5) \times SU(3) \times SU(3)^{\prime }\times U(1)_X \times U(1)^{\prime }. \end{aligned}$$

(1)

The charges of the fermions with respect to the flavor group are described in Table 1.

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

$$\begin{aligned} \frac{G}{H} = \frac{SU(5)\times SU(3) \times SU(3)' \times U(1)_X \times U(1)' }{SO(5) \times SU(3)_c \times U(1)_X }. \end{aligned}$$

(2)

This symmetry breaking pattern is interesting for several reasons. Given that $SO(5) \supset SU(2)\times SU(2),$ the pattern is compatible with the requirement of custodial symmetry $H \supset G_{\text {cust}} \supset G_{\text {SM}},$ with $G_{\text {cust}} = SU(3)_c \times SU(2)_L \times SU(2)_R \times U(1)_X$ and $G_{\text {SM}}$ being the SM gauge group $SU(3)_c \times SU(2)_L \times U(1)_Y.$ 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 \times U(1)_Y.$ In fact, $SO(5) \supset SO(4) \simeq SU(2)_L \times 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^{(3)}$: the correct hypercharges Y are then obtained by $Y=T_R^{(3)}+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,

$$\begin{aligned} \varvec{14} \rightarrow {\varvec{1}}_0 + {\varvec{2}}_{\pm 1/2} + {\varvec{3}}_0 \pm {\varvec{3}}_{\pm 1} \equiv (\eta ,h,\phi _0,\phi _\pm ), \end{aligned}$$

(3)

we identify an SU(2) doublet ${\varvec{2}}_{\pm 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 $\varPsi$ [12] in the low energy theory that has charges $(\varvec{5},\varvec{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.

2.2 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 [6,7,8,9, 11]. 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). The symmetry breaking pattern is the same as in massless two-flavor QCD, characterised by the quotient group

$$\begin{aligned} \frac{SU(2)_{\textrm{L}} \times SU(2)_{\textrm{R}}}{SU(2)_{\textrm{V}}}, \end{aligned}$$

(4)

and does not require to be discussed here. For completeness, we only mention that the three Goldstone bosons $\pi _1, \pi _2, \pi _3$ arising from these cosets can be labelled with eigenvalues of the azimuthal component of the isospin, $\pm 1, 0$:

$$\begin{aligned}{} & {} \pi _{+} = \frac{\pi _1 + i \pi _2}{\sqrt{2}},\nonumber \\{} & {} \pi _{-} = \frac{- \pi _1 + i \pi _2}{\sqrt{2}}, \nonumber \\{} & {} \pi _{0} = -\pi _3. \end{aligned}$$

(5)

The multiplet $\left( \pi _+, \pi _0,\pi _-\right)$ has eigenvalue $-1$ under the G-parity defined by combining charge conjugation ${\mathscr {C}}$ with an ${\textrm{SU}}(2)$-isospin rotation,

$$\begin{aligned} G = \exp \left( i \pi \tau _2 \right) \; {\mathscr {C}}, \end{aligned}$$

(6)

$\tau _i$ being SU(2) generators. In the 2AS sector, instead, the symmetry breaking pattern yields the cosets

$$\begin{aligned} \frac{SU(4)}{SO(4)}. \end{aligned}$$

(7)

The generators $T_i$ of SU(4) can be found in Appendix A. These cosets are characterised by 9 Goldstone bosons $\varPi _i,$ $i=1,\dots 9,$ that we represent exponentially as

$$\begin{aligned} U = \exp (i {\hat{T}}_i \varPi _i), \end{aligned}$$

(8)

where ${\hat{T}}_i,$ $i=1,\dots 9$ are the broken generators of SU(4). Under an isospin transformation,

$$\begin{aligned} U \rightarrow h U h^\dagger ,\quad h=\exp (i \omega _n X_n),\ n=1,\dots 6, \end{aligned}$$

(9)

where $X_n$ are the unbroken generators of SU(4) generating SO(4),  represented as $9\times 9$ matrices. We give such a representation in Appendix D. 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 $\varPi _{a_1,a_6},$ where $a_n$ are eigenvalues of the generators $X_n$

$$\begin{aligned}{} & {} \varPi _{-1,0}= -i\varPi _1+\varPi _2, \nonumber \\{} & {} \varPi _{1,0}= i\varPi _1+\varPi _2,\nonumber \\{} & {} \varPi _{-\frac{1}{2},-\frac{1}{2}}= \frac{-\varPi _3+i\varPi _4+i\varPi _6+\varPi _7}{2},\nonumber \\{} & {} \varPi _{\frac{1}{2},-\frac{1}{2}}= \frac{\varPi _3+i\varPi _4-i\varPi _6+\varPi _7}{2},\nonumber \\{} & {} \varPi _{-\frac{1}{2},\frac{1}{2}}= \frac{\varPi _3-i\varPi _4+i\varPi _6+\varPi _7}{2},\nonumber \\{} & {} \varPi _{\frac{1}{2},\frac{1}{2}}= \frac{-\varPi _3-i\varPi _4-i\varPi _6+\varPi _7}{2},\nonumber \\{} & {} \varPi _{0,-1} = -\frac{1}{\sqrt{2}} \varPi _5 -i\sqrt{\frac{3}{2}}\varPi _8+\varPi _9, \nonumber \\{} & {} \varPi _{0,1} = -\frac{1}{\sqrt{2}} \varPi _5 +i\sqrt{\frac{3}{2}}\varPi _8+\varPi _9, \nonumber \\{} & {} \varPi _{0,0}=\sqrt{2}\varPi _5 + \varPi _9. \end{aligned}$$

(10)

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 Sect. 7.2.

3 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 Sect. 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.

3.1 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 R,  we write the perturbative expansion of the one particle irreducible two-point function as

$$\begin{aligned} \varSigma _R\left( p,g_0,m^R_0\right)= & {} g_0^2 \varSigma _R^{(1)}\left( p,m^R_0\right) \nonumber \\{} & {} + g_0^4 \varSigma _R^{(2)}\left( p,m^R_0\right) + \mathscr {O}(g_0^6). \end{aligned}$$

(11)

In evaluating the Feynman diagrams contributing to the elements of the perturbative expansions, we set the fermion masses to their tree level values, $m^{R}_0=0.$ The one-loop contribution $\varSigma _R^1,$ for instance, comes from two diagrams, a tadpole and a sunset. It is useful to parametrise $\varSigma _R^1$ in terms of powers of the lattice spacing,

$$\begin{aligned} \varSigma _R^{(1)}(p,0) = \frac{\varSigma _{R, \, {\textrm{a}}}}{a}+ i {p}\!\!\!/ \varSigma _{R,\, {\textrm{b}}}. \end{aligned}$$

(12)

Imposing the vanishing of the renormalised mass for the fermion yields an expression for the critical mass at one loop

$$\begin{aligned} m_{{\textrm{c}}, \,R}^{(1)} = \frac{g_0^2 \varSigma _R^{{(1)}}(0,0)}{a} = \frac{g_0^2 \varSigma _{R,\, {\textrm{a}}}}{{\textrm{a}}}. \end{aligned}$$

(13)

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

$$\begin{aligned} \beta = \frac{2N}{g_0^2}. \end{aligned}$$

(14)

The result for $\varSigma _R^{(1)}$ can be written in terms of the contribution of each diagram [20]

$$\begin{aligned} \varSigma ^{(1)}_R(0,0) = 2C_2(R) \left[ c_1^{(1)} + c_2^{(1)} \right] , \end{aligned}$$

(15)

where $C_2(R)$ is the quadratic Casimir defined in Eq. (A.3), and

$$\begin{aligned} c_1^{(1)} + c_2^{(1)} = -0.162857058711(2). \end{aligned}$$

(16)

In the following we omit, in the expansion terms of the self energy, the dependence from p and $m_0$ which are set to zero. The result of Eq. (15), 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 $\mathscr {O}(g_0^4)$ result. In Wilson lattice QCD, $\varSigma ^{(2)}$ takes contribution from 26 diagram [20].

Fig. 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]

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 $I_2$ in Fig. 1 must be evaluated together in order to be infrared finite. Their value is

$$\begin{aligned} I^{R}_1+I^{R}_2 = n_f(R') 4 C_2(R) T_{R'} \, c_1^{(2)}, \end{aligned}$$

(17)

where $c^{(2)}_1=0.00079263(8)$ [21] is representation independent, and $n_f(R)$ is the number of fermions in the representation R. Similarly, $I_3+I_4$ gives the infrared-finite result

$$\begin{aligned} I_3^R + I_4^R = n_f(R') 4 C_2(R) T_{R'} \; c^{(2)}_2, \end{aligned}$$

(18)

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

$$\begin{aligned} \varSigma _{R'\rightarrow R}^{(2),\; \text {multi-rep}} = n_f(R') 4 C_2(R) T_{R'}\; 0.00118619(9). \end{aligned}$$

(19)

We then add to this contribution the 26 single-representation diagrams $\varSigma ^{(2), \; \text {one-rep}}_{R\rightarrow R}$ [22]

$$\begin{aligned}&\varSigma ^{(2), \; \text {one-rep}}_{R\rightarrow R} = C_2(R) N \, k_1 + 2C_2(R) T_R n_f(R)\, k_2 \nonumber \\&\quad + C_2(R) C_2({\textrm{Fund}})\, k_3 + C_2(R)^2 \,k_4, \end{aligned}$$

(20)

with

$$\begin{aligned}{} & {} k_1 = -0.001940(6),\quad k_2 = 0.00237236(16),\nonumber \\{} & {} k_3=-0.081429(8),\quad k_4 = 0.01516325(12). \end{aligned}$$

(21)

The total two-loop self energy is given by the sum of Eqs. (19) and (20)

$$\begin{aligned} \varSigma ^{(2)}_R&= C_2(R) N \, k_1 + 2C_2(R)\left[ T_R n_f(R) + T_{R'} n_f(R') \right] \, k_2 \nonumber \\&\quad + C_2(R) C_2({\textrm{Fund}})\, k_3 + C_2(R)^2 \,k_4. \end{aligned}$$

(22)

Unsurprisingly if $R=R'$ the extra term is equivalent to the addition of extra flavor content.

We can now list results for $N=4,$ $n_f({\textrm{Fund}})=n_f({2{\textrm{AS}}})=2.$ One loop values,

$$\begin{aligned}{} & {} \varSigma ^{(1)}_{\textrm{Fund}} = -0.610713970166(8),\nonumber \\{} & {} \varSigma ^{(1)}_{2{\textrm{AS}}} = -0.814285293555(10), \end{aligned}$$

(23)

and two-loop values,

$$\begin{aligned}{} & {} \varSigma ^{(2)}_{\textrm{Fund}} = -0.220826(53),\nonumber \\{} & {} \varSigma ^{(2)}_{2{\textrm{AS}}} = -0.270743(71). \end{aligned}$$

(24)

We are interested, in particular, in the contributions from one representation to the other. These are

$$\begin{aligned}{} & {} \varSigma _{{2{\textrm{AS}}}\rightarrow {\textrm{Fund}}}^{(2), \; \text {multi-rep}} = 0.0213512(14),\nonumber \\{} & {} \varSigma _{{\textrm{Fund}}\rightarrow {2{\textrm{AS}}}}^{(2), \; \text {multi-rep}} = 0.0355854(24). \end{aligned}$$

(25)

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 non-perturbative counterpart in Sect. 6.

3.2 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^{(1)}_{{\textrm{c}}, \, {\textrm{R}}}$ of a Wilson fermion can be computed from Eq. (15). Considering the clover term yields new contributions in powers of $c^R_{\textrm{sw}}$ [23],

$$\begin{aligned} m_{{\textrm{c}},\, R} = \frac{g_0^2 C_2(R)}{16 \pi ^2} \left( \epsilon _0 + \epsilon _1 c^R_{\textrm{sw}} + \epsilon _2 (c^R_{\textrm{sw}})^2 \right) , \end{aligned}$$

(26)

where the coefficients $\epsilon _i$ can be read from Table 2 of Ref. [23]. Concerning the Wilson term, the cactus resummation is performed by rescaling the self-energy $\varSigma \rightarrow \varSigma /\tilde{c_0},$ where $\tilde{c_0}$ is a function of N and $g_0^2$ but not of the representation. ${\tilde{c}}_0$ can be found by solving the following equation,

Fig. 2

The solution of Eq. (27) for $N=2,3,4$ and $\beta =11,$ the value used in the numerical simulations, as will be described in Sect. 4

$$\begin{aligned}{} & {} u e^{-u(N-1)/(2N)} \left[ \frac{N-1}{N} L^1_{N-1}(u) + 2L^2_{N-2}(u) \right] \nonumber \\{} & {} \quad = \frac{g_0^2(N^2-1)}{4}, {\tilde{c}}_0 \equiv \frac{g_0^2}{4u}, \end{aligned}$$

(27)

where $L_N^\alpha$ are generalised Laguerre polynomials of degree N. The solution for $\beta =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

$$\begin{aligned} {\tilde{c}}_0 = 0.731607. \end{aligned}$$

(28)

The resummation with the clover improvement term corresponds to also rescale $g_0^2\rightarrow g_0^2/{\tilde{c}}_0$ and $c^R_{\textrm{sw}} \rightarrow c^R_{\textrm{sw}} {\tilde{c}}_0$ [24]. The prediction for the critical mass is then

$$\begin{aligned} m_{{\textrm{c}}, \, R}^{\mathrm {1-loop+cactus}}= \frac{g_0^2}{{\tilde{c}}_0} \frac{C_2(R)}{16 \pi ^2} \left( \epsilon _0 + \epsilon _1 c^R_{\textrm{sw}} {\tilde{c}}_0 + \epsilon _2 (c^R_{\textrm{sw}})^2 {\tilde{c}}_0^2 \right) .\nonumber \\ \end{aligned}$$

(29)

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.

Table 2 Critical masses for different lattice actions and improvements. When not specified, the results are obtained at 1 loop. These values are to be compared with the non-perturbative results of Sect. 6

4 Lattice setup

The lattice setup has already been discussed in Ref. [9], and here we only recall the main ideas. We generate gauge configurations by using the HMC algorithm [25] with a second order integrator [26]. The integration scheme has two levels, one updating the gauge force and one updating the fermionic force. The lattice action can be decomposed into gauge and fermionic part. The latter is further decomposed into contributions from each representation [27]

$$\begin{aligned} S = S_{\textrm{g}} + S_{\textrm{f}},\quad S_{\textrm{f}} = S^{({{\textrm{Fund}}})} + S^{(2{\textrm{AS}})}. \end{aligned}$$

(30)

For the gauge part we use the Wilson plaquette action

$$\begin{aligned} S_{\textrm{g}} = \frac{\beta }{N_c} \sum _{x} \sum _{\mu < \nu } \; \text {Re}\; \text {Tr} \left\{ 1 - {\mathscr {P}}_{\mu \nu }(x) \right\} . \end{aligned}$$

(31)

The fermionic action for the representation R is

$$\begin{aligned} S_{\textrm{f}}^{R} = \sum _x {\bar{\psi }}^{R}(x) \, D^R_{x,y} \, \psi (x)^{R}, \end{aligned}$$

(32)

where, for both representations, we collect the fermionic degrees of freedom into a doublet of Dirac fermions $\psi ^R.$ We adopt a Wilson discretisation for the Dirac action $D^R$ in each representation R,  plus a clover improvement term

$$\begin{aligned} D^R = D_{\textrm{Wilson}}^R + D_{\textrm{clover}}^R. \end{aligned}$$

(33)

The Wilson term in position space is

$$\begin{aligned} D^R_{x,y}&= \delta _{x,y} - \kappa ^R \sum _{\mu =1}^4 \left[ \left( \mathbbm {1}- \gamma _{\mu } \right) U^R_{\mu }(x) \delta _{x+a{\hat{\mu }},y} \right. \nonumber \\&\quad \left. + \left( \mathbbm {1}+ \gamma _{\mu } \right) U_{\mu }^R\,^\dagger (y) \delta _{x-a{\hat{\mu }},y} \right] , \end{aligned}$$

(34)

where

$$\begin{aligned} U^R_{\mu }(x) = \exp (i\,\omega ^a_{\mu }(x) \, T_a^R), \end{aligned}$$

(35)

and $\kappa ^R$ is related to the bare mass $m_0^R$ of the fermion in the representation R

$$\begin{aligned} \kappa ^R = \frac{1}{2(am_0^R+4)}. \end{aligned}$$

(36)

The order a clover improvement term is

$$\begin{aligned} \left( D^R_{\textrm{clover}} \right) _{x,y} = \frac{ia}{2} \, c^R_{\textrm{sw}}(g_0^2) \, \kappa ^R \, \sum _{\mu \nu } {\tilde{F}}^R_{\mu \nu }(x) \, \sigma _{\mu \nu } \delta _{x,y}, \end{aligned}$$

(37)

where $\sigma _{\mu \nu }=\frac{i}{2}[\gamma _{\mu }, \gamma _\nu ]$ and

$$\begin{aligned} {\tilde{F}}_{\mu \nu }(x) = \frac{1}{8} \left[ Q_{\mu \nu }(x) - Q_{\nu \mu }(x) \right] ,\quad Q_{\mu \nu }(x) = Q_{\nu \mu }^\dagger (x), \end{aligned}$$

(38)

$Q_{\mu \nu }(x)$ being the clover combination of plaquettes around the point x [28]. The improvement coefficient $c^R_{\textrm{sw}}$ can be computed in powers of the coupling,

$$\begin{aligned} c^R_{\textrm{sw}}(g_0^2) = 1 + c_{\textrm{sw}}^{R\,(1)} \, g_0^2 + \mathscr {O}(g_0^4). \end{aligned}$$

(39)

A discussion on the $O(g_0^2)$ coefficient can be found in [29]. In this work, we set the coefficient to its tree level value, $c^R_{\textrm{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_{\textrm{c}}^R,$ i.e. the value of the bare masses at which each fermion has a renormalised vanishing mass

$$\begin{aligned} am^R = am_0^R - am_{\textrm{c}}^R = \frac{1}{2}\, \left( \frac{1}{\kappa ^R} - \frac{1}{\kappa ^R_{\textrm{c}}} \right) . \end{aligned}$$

(40)

Direct simulations at the chiral point are in fact impossible in our setup due to the spectral properties of the Dirac operator. The lowest eigenvalue of the Dirac operator approaches zero in the chiral limit, meaning its inverse becomes increasingly ill-conditioned. When extrapolating to the critical point (see Eq. (40)) one needs to care that exceptional configurations do not occur in the gauge average [30]. 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 operators in both representations 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.

Fig. 3

Distribution of the lowest eigenvalue of the Wilson-clover operator for bare parameters $am_0^{({{\textrm{Fund}}})}=-0.45,$ $\beta =11$ and different masses for the 2AS fermions: $am_0^{(2{\textrm{AS}})}=-0.59$ on the top and $am_0^{(2{\textrm{AS}})}=-0.60$ on the bottom. Both distributions are far enough from the origin to ensure no exceptional configurations enter the ensemble

Our starting point is the work done in [9], 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 E. 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 $\beta =11,$ and at a single volume corresponding to a lattice of dimensions $(L/a)^3 \times T/a = 16^3 \times 32.$ For the production of gauge configurations and the measurements of observables, we use the software Grid [31] and Hadrons [32].

5 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,

$$\begin{aligned}{} & {} O_{\;\textrm{P}}^{R}(x) = {\bar{\psi }}_{f_1}^R(x)\, \gamma _5 \,\psi _{f_2}^R(x),\nonumber \\{} & {} O_{\;\textrm{A}}^{R}(x) = {\bar{\psi }}_{f_1}^R(x) \,\gamma _0 \gamma _5 \, \psi _{f_2}^R(x),\quad f_1 \ne f_2. \end{aligned}$$

(41)

In the previous equation we use the same notation of Sect. 4, i.e. $\psi ^R$ collects the fermionic degrees of freedom of the representation R into a flavor doublet of Dirac spinors. The flavor indices are denoted by $f_1, f_2.$ In this work we only consider isospin-vector operators, and we build the following two-point functions

$$\begin{aligned} C_{ab}^{R}(t) = \frac{1}{L^3} \sum _{\varvec{x}} \langle {O^{R}_a(\varvec{x},t) O^{R^\dagger }_b(0)}\rangle ,\quad a,b = {\textrm{P}},{\textrm{A}}. \end{aligned}$$

(42)

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. Considering for instance the pseudoscalar–pseudoscalar correlator, we obtain

$$\begin{aligned} C_{PP}^{R}(t)&= \left( e^{-tM^R_{PP}} + e^{(-T+t)M^R_{PP}} \right) \nonumber \\&\quad \times \frac{\langle 0|O^{R}_P(0)|M^R_{PP} \rangle _L \langle M^R_{PP}|O^{R^\dagger }_P(0)|0\rangle _L}{2M^R_{PP}} +\cdots , \end{aligned}$$

(43)

where the ellipsis denotes terms that do not correspond to the state $\left| M_{PP}^R\right\rangle _L$ and that are exponentially suppressed. The mass $M_{PP}^R$ can be then obtained by the asymptotic behaviour in the Euclidean time of the effective mass

$$\begin{aligned} aM^{R}_{{\textrm{PP}}}(t) = \cosh ^{-1}\left[ \frac{ C^{R}_{{\textrm{PP}}}(t+a) + C^{R}_{{\textrm{PP}}}(t-a)}{2C^{R}_{{\textrm{PP}}}(t)} \right] . \end{aligned}$$

(44)

Similarly the PCAC mass, defined through the axial Ward identity, is obtained from the pseudoscalar–pseudoscalar and axial correlators

$$\begin{aligned} am^{R}_{\textrm{PCAC}} = \frac{\nabla _t C^{R}_{\textrm{AP}}(t)}{2 C^{R}_{{\textrm{PP}}}},\quad \nabla _t f(t) = \frac{ f(t+a) - f(t-a) }{2}, \end{aligned}$$

(45)

which has $\mathscr {O}(a)$ effects for our choice of the unimproved axial operator.

Fig. 4

Probability densities for the lowest eigenvalue of the fermionic operator (upper panel), the PCAC mass (central panel) and the pseudoscalar mass (lower panel) 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\sigma ,$ 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

In the chiral limit, the pseudoscalar mass $aM_{{\textrm{PP}}}^{R}$ for each representation vanishes. In our ensembles we generated configurations describing mesons of antisymmetric fermions being generally heavier: this 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, differently from the case shown in Fig. 3, the lowest eigenvalues of the Wilson operators in the two representations are the same, and the fermionic masses are also compatible within one $\sigma .$ 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 the large time behaviour of Eq. (44). Other energy levels and matrix elements can be estimated by fitting sums of exponentials, as in Eq. (43). 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 information provided. Moreover, as the infinite volume limit is approached, the spectrum above the multi-particle threshold becomes denser 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

$$\begin{aligned} C^{R}_{PP}(t) = \int _0^{\infty } dE \, \rho ^{R}_{PP}(E) \, \left( e^{-tE}+ e^{(-T+t)E} \right) , \end{aligned}$$

(46)

where we neglect the thermal effects inside $\rho _{PP}^R(E).$¹ 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 $\rho ^L_{ab}(E)$ can in fact be expanded as

$$\begin{aligned}&\rho _{PP}^{L,R}(E) \nonumber \\&\quad = \sum _n \frac{\langle 0|O^{R}_P(0)|n_\rangle L \langle n|O^{R^\dagger }_P(0)|0_\rangle L}{2E_n(L)}\, \delta \left( E-E_n(L) \right) . \end{aligned}$$

(47)

In order to cope with the distributional nature of the spectral density, we smear it [13] with a Gaussian function, $\varDelta _{\sigma }(E) = \exp {(-E^2/2\sigma ^2)}/\sqrt{2\pi }\sigma$

$$\begin{aligned} \rho ^{L,R}_{\sigma , PP}(E)=\int _0^\infty dE^{\prime }\, \varDelta _{\sigma }(E-E') \, \rho _{PP}^{L,R}(E'), \end{aligned}$$

(48)

so that the smeared spectral density $\rho ^{L,\,R}_{\sigma , 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 $\sigma .$ By rewriting Eq. (47) for the smeared spectral density,

$$\begin{aligned}{} & {} \rho _{\sigma , PP}^{L,R}(E) \nonumber \\{} & {} \quad = \sum _n \frac{\langle 0|O^{R}_P(0)|n_\rangle L \langle n|O^{R^\dagger }_P(0)|0_\rangle L}{2E_n(L)}\, \varDelta _{\sigma } \left( E-E_n(L) \right) ,\nonumber \\ \end{aligned}$$

(49)

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. (43) 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.

6 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. Around the chiral point we also expect the Gell-Man–Oakes–Renner relation to be valid, with the pseudoscalar mass squared $(M^{R}_{{\textrm{PP}}})^2$ scaling as $m_{\textrm{PCAC}}^{R}$ 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.

Fig. 5

Three chiral extrapolations, one on the top for fundamental fermions and two on the bottom for antisymmetric fermions. The bands correspond to linear fits. The heavier point in the plot on the top was also present in [9]. In the top panel, the fit has $\chi ^2/d.o.f.=0.23.$ In the bottom panel, $\chi ^2/d.o.f.=2.22$ for $am_0^{({\textrm{Fund}})}=-0.45$ and $\chi ^2/d.o.f.=0.48$ for $am_0^{({\textrm{Fund}})}=-0.47$

The measurements of the PCAC masses are listed in Tables 6 and 7 of Appendix F. 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. Figure 5 shows three of these extrapolations. In the top panel, we extrapolate the chiral point for the fundamental fermions by fixing $am_0^{(2{\textrm{AS}})}=-0.45.$ In the bottom panel, the two lines represent different extrapolations for the 2AS chiral point taken at different values of the bare fundamental mass, $am_0^{({{\textrm{Fund}}})}=-0.45$ and $am_0^{({{\textrm{Fund}}})}=-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 Sect. 3, where we have found the critical mass to only mildly depend on the other representation. The critical masses for $\beta =11.0$ are

$$\begin{aligned}{} & {} am_{\textrm{c}}^{({{\textrm{Fund}}})}\bigg |_{am_0^{(2{\textrm{AS}})}=-0.45} = -0.486(3)\nonumber \\{} & {} am_{\textrm{c}}^{(2AS)}\bigg |_{am_0^{({\textrm{Fund}})}=-0.45} = -0.637(3)\nonumber \\{} & {} am_{\textrm{c}}^{(2{\textrm{AS}})}\bigg |_{am_0^{({\textrm{Fund}})}=-0.47} =-0.630(4). \end{aligned}$$

(50)

Fig. 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. In the top panel, the fit has $\chi ^2/d.o.f.=1.68,$ while in the bottom panel $\chi ^2/d.o.f.=1.50$

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_{{\textrm{PP}}}^{R}$ with respect to the quark masses. For both representations R the scaling is compatible with the one predicted by the chiral Lagrangians, $(M^{R}_{{\textrm{PP}}})^2 \sim m^{R}_{\textrm{PCAC}},$ with the pseudoscalar mesons becoming massless in the chiral limit. The pseudoscalar masses, reported in Table 8, are obtained by fitting the effective masses of the pseudoscalar–pseudoscalar correlators.

7 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.

7.1 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 the spectral density $\rho _{PP}^{R}(E)$ from the Euclidean correlator $C_{PP}^{R}(t)$ involves an inverse Laplace transform,

$$\begin{aligned} C_{PP}^{R}(t) = \int _{E_{\textrm{min}}}^\infty \, dE\, \left( e^{-tE}+e^{(-T+t)E} \right) \, \rho _{PP}^{R}(E), \end{aligned}$$

(51)

where $E_{\textrm{min}}$ can range between zero and the energy of the ground state, since $\rho _{PP}^{R}$ vanishes in that interval. The inversion of Eq. (51) is numerically an ill-posed problem which needs to be regularised. Recalling from Sect. 5 that we are interested in the smeared version of $\rho _{PP}^{R}(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), 

$$\begin{aligned} \rho _{PP}^{R} [f] = \int _0^\infty \, dE \, f(E) \, \rho _{PP}^{R}(E). \end{aligned}$$

(52)

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 $\varDelta _{\sigma }(E)$ of Eq. (48), as a linear combination of the same exponentials appearing in the Laplace transform of Eq. (51),

$$\begin{aligned} \varDelta _{\sigma }(E-E^{\prime }) = \sum _{\tau =1}^{\infty } \, g_\tau \, \left( e^{-\tau a E}+e^{(-T+a\tau )E} \right) , \end{aligned}$$

(53)

where $t=\tau a,$ a being the lattice spacing. Once the coefficients $g_{\tau } \equiv g_\tau (\sigma ,E^{\prime })$ are known one can simply obtain an estimator for the smeared spectral density,

$$\begin{aligned} \rho _{PP, \, \sigma }^{R}(E^{\prime }) = \sum _{\tau =1}^{\infty } \, g_\tau \, C_{PP}^{R}(a \tau ). \end{aligned}$$

(54)

On the lattice the correlators $C_{PP}^{R}$ are available for a finite number of times, therefore the sum in Eqs. (53) and (54) must be truncated at the appropriate cutoff $\tau _{\textrm{max}}.$ Since our lattices have temporal length T and periodic boundaries, we have $a\tau _{\textrm{max}}=T/2.$ The reconstructed smearing kernel,

$$\begin{aligned} f(E,\varvec{g}) = \sum _{\tau =1}^{\tau _{\textrm{max}}} g_\tau \, \left( e^{-\tau a E}+e^{(-T+a\tau )E} \right) , \end{aligned}$$

(55)

will necessarily differ from the Gaussian $\varDelta _{\sigma }(E)$ at finite $\tau _{\textrm{max}},$ inducing a systematic error on the final result. The computation of the coefficients $\varvec{g}$ is achieved through the minimisation of the functional $W_{\alpha }[\varvec{g}]$

$$\begin{aligned} W_{\alpha }[\varvec{g}] = \frac{A_{\alpha }[\varvec{g}]}{A_{\alpha }[0]} + \lambda \,B[\varvec{g}], \end{aligned}$$

(56)

where $\lambda$ is a trade-off input parameter that we will discuss later in this section. The functional $A_{\alpha }[\varvec{g}],$ introduced in [14], measures the difference between the exact smearing kernel and the one we can reconstruct with the available data

$$\begin{aligned} A_{\alpha }[\varvec{g}] = \int _{E_{\textrm{min}}}^{\infty } dE \,e^{\alpha aE} \left| f(E,\varvec{g}) - \varDelta _{\sigma }(E-E^{\prime }) \right| ^2. \end{aligned}$$

(57)

The parameter $\alpha <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 $\alpha$ allows for the integrand to decay faster at high energies. The functional $B[\varvec{g}]$ is needed to regularise the problem [34], making it numerically stable. We define $B[\varvec{g}]$ to be dimensionless, namely

$$\begin{aligned} B[\varvec{g}] = \frac{E^2}{C_{PP}^{R}(a)^2}\; \sum _{\tau ,\tau ^{\prime }=1}^{\tau _{\textrm{max}}} g_\tau \, \text {Cov}_{\tau \tau ^{\prime }} \, g_{\tau ^{\prime }}, \end{aligned}$$

(58)

where $\text {Cov}$ is the covariance matrix of the correlator $C(a \tau )$ estimated over N bins (see Appendix C),

$$\begin{aligned}{} & {} C(a\tau )= \frac{1}{N}\sum _{n=0}^{N-1} C_n(a\tau ),\nonumber \\{} & {} \text {Cov}_{\tau \tau ^{\prime }}[C] =\frac{1}{N-1} \sum _{n=0}^{N-1} \left[ C_n(a\tau ) - C(a\tau ) \right] \,\nonumber \\{} & {} \qquad \qquad \qquad \quad \times \left[ C_n(a\tau ^{\prime }) - C(a\tau ^{\prime }) \right] . \end{aligned}$$

(59)

The algorithmic parameters can be gathered to simplify the notation

$$\begin{aligned} \varvec{p} = (\alpha ,\lambda , E_{\textrm{min}}, \tau _{\textrm{max}}). \end{aligned}$$

(60)

The minimisation of $W_{\alpha }[\varvec{g}]$ corresponds to solving the following linear problem,

$$\begin{aligned} \frac{\delta W_{\alpha }[\varvec{g}]}{\delta g_\tau }\bigg |_{\varvec{g} = \varvec{g}^{\varvec{p}}} = 0, \end{aligned}$$

(61)

which has to be performed at each energy and smearing radius for which we want to estimate $\rho _{PP, \, \sigma }^{R}(E^{\prime }).$ The nature of the functionals appearing in these definitions is intimately related to the uncertainties on the estimator $\rho _{PP, \, \sigma }^{R}(E^{\prime }).$ The statistical error is in fact estimated by

$$\begin{aligned} \varDelta _{\textrm{stat}}(E^{\prime },\varvec{g}^{\varvec{p}}) = \frac{C_{ab}^{R}(1)}{E^{\prime }} \sqrt{ B[\varvec{g}^{\varvec{p}}]}. \end{aligned}$$

(62)

The systematic error, unavoidable at finite $\tau _{max},$ is estimated by monitoring the quantity

$$\begin{aligned} d(\varvec{g}^{\varvec{p}}) = \sqrt{ \frac{A_0[\varvec{g}^{\varvec{p}}]}{A_0[\varvec{0}]} }, \end{aligned}$$

(63)

as we now describe. Regions where $d(\varvec{g}^{\varvec{p}})$ is small are dominated by the statistical uncertainty; on the other hand, the reconstructed smearing function of Eq. (53) will be as close as possible to the exact one. For these reasons, it is not surprising that where $d(\varvec{g}^{\varvec{p}})$ is small, the results for $\rho _{PP, \, \sigma }^{R}$ are stable² in response to variations of the algorithmical, unphysical parameters $\varvec{p},$ within statistical error. If we define the coefficients $\varvec{g}^{\star }$ as

$$\begin{aligned} A_{\alpha }[\varvec{g}^{\star }] / A_{\alpha }[\varvec{0}] = k \; B[\varvec{g}^{\star }], \end{aligned}$$

(64)

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 $\alpha =1$ the outcome of the reconstruction is in a region of algorithmical stability: $d(\varvec{g}^{\varvec{p}})$ is small, and systematic fluctuations are well within the statistical uncertainty. When this is not realised, the systematic error is estimated as

$$\begin{aligned} \varDelta _{\textrm{sys}}(E^{\prime }) = |\rho _{PP, \, \sigma }^{R}(E^{\prime }, \varvec{g}^\star )- \rho _{PP, \, \sigma }^{R}(E^{\prime }, \varvec{g}^{\star \star })|, \end{aligned}$$

(65)

where $\varvec{g^{\star \star }}$ is defined through Eq. (64) 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^{\star \star }$ at k/10 provides a conservative estimate of $\varDelta _{\textrm{sys}}(E^{\prime })$ for a reconstruction performed with $g^{\star }$ at the value k. Figure 7 shows an example, for a specific energy $E_{\star },$ of a stability regime where the spectral reconstruction for different values of $\lambda$ (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. (64), showing that indeed no systematic component affects the uncertainty on the result. The value corresponding to $k=1$ is associated to the value of $\lambda$ at which one achieves the optimal balance $A_{\alpha }[\varvec{g}^{\star }] / A_{\alpha }[\varvec{0}] =\; B[\varvec{g}^{\star }],$ in agreement with the prescription from Refs. [14, 16].

Fig. 7

Example of region of algorithmical stability at a given energy $E_{\star }.$ Different values of $\lambda ,$ which translate into different values of $d(\varvec{g}^{\varvec{p}})$ on the x-axis, produce predictions for the smeared spectral density $\rho _{PP, \, \sigma }^{R}(E_{\star })$ that are compatible within statistical error (black bars). In this case, $R={2{\textrm{AS}}},$ $\sigma =0.21/a$ and $E_{\star } \simeq M_{{\textrm{PP}}}^{(2{\textrm{AS}})}.$ The green and orange points correspond, according to Eq. (64), 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 $\lambda$ at which one achieves the optimal balance $A_{\alpha }[\varvec{g}^{\star }] / A_{\alpha }[\varvec{0}] =\; B[\varvec{g}^{\star }]$

7.2 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 for the pseudoscalar interpolators of Eq. (41)

$$\begin{aligned}{} & {} {\hat{\pi }}(x) = {\bar{\psi }}_{f_1}^{({{\textrm{Fund}}})}(x) \gamma _5\, \, \psi _{f_2}^{({{\textrm{Fund}}})}(x), \nonumber \\{} & {} {\hat{\varPi }}(x) = {\bar{\psi }}_{f_1}^{(2{\textrm{AS}})}(x) \gamma _5 \, \psi _{f_2}^{(2{\textrm{AS}})}(x),\quad f_1 \ne f_2, \end{aligned}$$

(66)

and the correlation functions

$$\begin{aligned}{} & {} C_{{\textrm{PP}}}^{({{\textrm{Fund}}})}(t) = \frac{1}{L^3} \sum _{\varvec{x}} \langle {\hat{\pi }}(\varvec{x},t) {\hat{\pi }}^\dagger (0),\rangle \nonumber \\{} & {} C_{{\textrm{PP}}}^{(2{\textrm{AS}})}(t) = \frac{1}{L^3} \sum _{\varvec{x}} \langle {\hat{\varPi }}(\varvec{x},t) {\hat{\varPi }}^\dagger (0).\rangle \end{aligned}$$

(67)

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

$$\begin{aligned} \rho _{\sigma , PP}^{L,R}(E)= & {} \sum _n \frac{\langle 0|O^{R}_{\textrm{P}}(0)|n_\rangle L \langle n|O^{R^\dagger }_{\textrm{P}}(0)|0_\rangle L}{2E_n(L)}\,\nonumber \\{} & {} \times \varDelta _{\sigma } \left( E-E_n(L) \right) . \end{aligned}$$

(68)

The magnitude of each matrix element $\langle n|\bar{O}^{R}_{\textrm{P}}(0)|0_\rangle L$ determines the weight of each Gaussian $\varDelta _{\sigma }$ 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, point-like operators and Gaussian-smeared ones.³ 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

$$\begin{aligned} \langle 0|{\hat{\pi }}|\pi ,\rangle \quad \langle 0|{\hat{\pi }}|\pi \pi \pi ,\rangle \quad \langle 0|{\hat{\pi }}|\pi \varPi \varPi ,\rangle \dots \end{aligned}$$

(69)

that will enter in our analysis through Eq. (68). Since in our simulations $M_{{\textrm{PP}}}^{(2{\textrm{AS}})} > M_{{\textrm{PP}}}^{({{\textrm{Fund}}})}$ this phenomenology is reminiscent, up to $E_{3\pi },$ 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\varPi }.$ 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

$$\begin{aligned}{} & {} \langle 0|{\hat{\varPi }}|\varPi ,\rangle \quad \langle 0|{\hat{\varPi }}|\varPi \,\varPi ,\rangle \quad \langle 0|{\hat{\varPi }}|\varPi \pi \pi ,\rangle \nonumber \\{} & {} \quad \langle 0|{\hat{\varPi }}|\varPi \pi \pi \pi \pi ,\rangle \dots \end{aligned}$$

(70)

The aforementioned features complicate the extraction of $M_{{\textrm{PP}}}^{(2{\textrm{AS}})},$ 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 ${\hat{\rho }}_{{\textrm{PP}}, \, \sigma }^{(2{\textrm{AS}})}$ 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 $\sigma$ of Eq. (48) 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. (56).

Fig. 8

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. (48) 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. Figure 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 \simeq 2 M_{{\textrm{PP}}}^{(2{\textrm{AS}})}.$ The observed smeared spectral density is the result of two contributions coming mainly from the energy levels $M_{{\textrm{PP}}}^{(2{\textrm{AS}})}$ and $E_{\varPi \varPi }\simeq 2 M_{{\textrm{PP}}}^{(2{\textrm{AS}})},$ which cannot be resolved because the smearing radius is too large, $\sigma \simeq M_{{\textrm{PP}}}^{(2{\textrm{AS}})}.$ 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.

Fig. 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_{{\textrm{PP}}}^{(2{\textrm{AS}})}.$ The right panel similarly shows how suppressed excited states allow for a clear peak to emerge in the spectral reconstruction smeared with $\sigma = M_{{\textrm{PP}}}^{(2{\textrm{AS}})}$ according to Eq. (48). The peak includes contributions from mainly $M_{\varPi }$ and $E_{\varPi \varPi }$

7.3 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 $g^{(k)}(t)$ for correlators and $f_ \sigma ^{(k)}(t)$ for the smeared spectral densities

$$\begin{aligned}{} & {} g^{(k)}(t) = \sum _{n=1}^k a_n \left( e^{-t E_n} + e^{(-T+t)E_n} \right) ,\nonumber \\{} & {} f^{(k)}_{\sigma }(E) = \sum _{n=1}^k b_n e^{-(E-E_k)^2/\, 2\sigma ^2}, \end{aligned}$$

(71)

where $\sigma$ is the smearing radius defined by Eq. (48). 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 $\chi ^2$ functions

$$\begin{aligned} \chi _{g^{(k)}}^2= & {} \sum _{t, r} \left( g^{(k)}(t) - C(t) \right) \text {Cov}^{-1}_{tr}[C] \left( g^{(k)}(r) - C(r) \right) ,\nonumber \\ \end{aligned}$$

(72)

$$\begin{aligned} \chi _{f_{\sigma }^{(k)}}^2= & {} \sum _{E,E'} \left( f_{\sigma }^{(k)}(E) - \rho _{\sigma }(E) \right) \nonumber \\{} & {} \times \text {Cov}^{-1}_{EE'}[\rho _{\sigma }] \left( f_{\sigma }^{(k)}(E') - \rho _{\sigma }(E') \right) , \end{aligned}$$

(73)

where covariance matrices are estimated as in Eq. (59) both for correlators and spectral densities.

On the lattice, the temporal length T constrains the maximum number of data points and hence degrees of freedom for fitting a correlator $C_{{\textrm{PP}}}^{R}(t).$ The effect of correlation between times t and $t'$ is taken into account by the covariance matrix appearing in $\chi _{g^{(k)}}^2.$ The smeared spectral density $\rho ^{R}_{{\textrm{PP}},\, \sigma }(E)$ depends on the correlator at all lattice times t,  and it can be in principle evaluated for any energy E. Not all points in energy, however, provide independent information. Indeed, the information contained in $\rho _{PP, \sigma }^R(E)$ and $\rho _{PP, \sigma }^R(E^\prime )$ for $\vert E-E^\prime \vert \ll \sigma$ is essentially the same, from both the physical and the statistical viewpoint. 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\le 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 $\text {Cov}[\rho _{\sigma }]$ in order to maximise the information that is passed to the $\chi ^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_{{\textrm{PP}}}^{(2AS)}(t)$ and $\rho _{{\textrm{PP}}, \, \sigma }^{(2{\textrm{AS}})}(E)$ from the ensemble B3 with a smearing radius $\sigma =0.2/a.$

In both cases, the numbers of data points are enough to fit at least two states. The matrix $\text {Cov}^{-1}_{t t^{\prime }}[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 $\text {Cov}[C]$ is ill-conditioned, it needs to be regularised by applying a cutoff on its smaller eigenvalues before $\chi ^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. (56): the matrices obtained from the functionals $A_{\alpha }[\varvec{g}]$ and $B[\varvec{g}],$ if both ill conditioned, regularise each other for suitable values of $\lambda .$ The choice of a cutoff for the covariance matrix $\text {Cov}[C]$ is therefore absorbed into the choice of the parameter $\lambda .$

Fig. 10

Covariance matrices for the lattice correlator $C_{{\textrm{PP}}}^{(2AS)}(t)$ at five time slices (top) and the smeared spectral density $\rho _{{\textrm{PP}}, \, \sigma }^{(2{\textrm{AS}})}(E)$ evaluated at six energies (bottom) from $C_{{\textrm{PP}}}^{(2{\textrm{AS}})}(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

Figure 11 shows two examples of correlated fits of the smeared spectral densities. On the top panel, the correlator used to extract $\rho _{{\textrm{PP}}, \, \sigma }^{(2{\textrm{AS}})}$ is built with smeared interpolating fields. The model function is $f_{\sigma }^{(3)}$ from Eq. (71). 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_{{\textrm{PP}}}^{(2{\textrm{AS}})}$ and the second with $a E_{\varPi \varPi }.$ The bottom 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. (44) 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 top panel. The choice of smeared operators is therefore preferred.

Fig. 11

Examples of fits of spectral densities, showing the breakdown of the contribution of each Gaussian. On the top 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 bottom 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 top panel is obtained within the ensemble B3 and the spectral density has smearing radius $\sigma =0.24/a.$ The bottom panel is derived from the ensemble B2 and the smearing radius of the spectral density is $\sigma =0.3/a$

We now turn to the comparison of the fit results. As described in Sect. 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 $a M_{{\textrm{PP}}}^{(2{\textrm{AS}})}$ 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 $\chi ^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 interpolators has been crucial, as it is clear from Fig. 8 where the effective mass of Eq. (44) 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 Sect. 7. Figure 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_{\star },$ is compatible with all the values generated by different choices for the unphysical parameter $\lambda .$

Fig. 12

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 $\lambda ,$ in the given range of $d(\varvec{g}^{\varvec{p}})$ (cf. Equation (63)). The blue band is the fit result of the smeared spectral density to a sum of Gaussians at the point $E_{\star }.$ Encouragingly, the fit is compatible with all points in the scan, showing that stability in the reconstruction translates into stability for the fits

The choice of the smearing radius is dictated by the quality of the data. In this work, we managed to obtain values ranging from $\sigma =0.18/a$ to 0.33/a. Since the separation between finite volume energies should be roughly $2\pi /L \simeq 0.4,$ these values can be considered acceptable. For each ensemble, we have performed the fit at a fixed value of $\sigma ,$ obtaining a prediction for the pseudoscalar mass $aM_{\mathrm {PP,\, \sigma }}^{(2{\textrm{AS}})}.$ 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_{\mathrm {PP,\, \sigma }}^{(2{\textrm{AS}})}$ as $\sigma$ varies. When they occurred, we added half the spread of these fluctuations as a systematic error. Figure 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.

Fig. 13

Fit results for $a M_{{\textrm{PP}}}^{(2{\textrm{AS}})}$ from the ensemble B3, obtained from two (green) and three (red) Gaussian fits of smeared spectral densities at different smearing radii $\sigma .$ 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 $a M_{{\textrm{PP}}}^{(2{\textrm{AS}})} = 0.3550(31)$

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.

Fig. 14

Graphical comparison between the two predictions for $a M_{{\textrm{PP}}}^{(2{\textrm{AS}})}$ for the ensembles B1–B4

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.

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

8 Conclusions

The model studied in this work allows us to understand the dynamics of partial compositeness better, by extending the effort started in Ref. [9] 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, 7 and 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 [14, 16] 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 [31, 32] 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 heavy-quark 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 cross-sections [13, 35], 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.

Appendices

Appendix A: Group-theoretical conventions

We denote the generators of SU(N),  the special unitary group of degree N,  as $T^a,$ with $1 \le a \le N^2-1.$ They are traceless Hermitian matrices satisfying the algebra

$$\begin{aligned} \left[ T^a,T^b\right] = i \sum _{c=1}^{N^2-1} f^{abc} T^c, \end{aligned}$$

(A.1)

where $f^{abc}$ are the structure constants of the su(N) Lie algebra, which are antisymmetric under permutations of every pair of indices, and are the same in all representations. The generators in an irreducible representation R,  denoted as $T^a_R,$ satisfy

$$\begin{aligned} {{\,\textrm{Tr}\,}}\left( T_R^a T_R^b \right) = \lambda _R \delta ^{ab}, \end{aligned}$$

(A.2)

where $\lambda _R$ is the Dynkin index of the representation R. In our normalisations, the Dynkin index of the fundamental representation is $\lambda _{{\tiny {\hbox {(Fund)}}}}=\frac{1}{2}.$ Another interesting group-theoretical invariant is the sum of the squares of the generators, which is the quadratic Casimir operator, and is proportional to the identity matrix in every representation; we denote the corresponding eigenvalue in the representation R as $C_2(R)$:

$$\begin{aligned} \sum _{a=1}^{N^2-1} T_R^a T_R^a = C_2(R) \mathbbm {1}. \end{aligned}$$

(A.3)

Comparing Eqs. (A.2) and (A.3) one obtains

$$\begin{aligned} \lambda _R = \frac{\text {dim}(R)}{N^2-1} \, C_2(R), \end{aligned}$$

(A.4)

$\text {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.

Table 4 Group-theoretical factors used in this work

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

$$\begin{aligned} T_{\tiny {\hbox {(Fund)}}}^k = \frac{1}{\sqrt{2k(k+1)}} \,{\textrm{diag}} ( \underbrace{1,\ldots ,1}_{k\,{\text {terms}}},-k,\underbrace{0,\ldots ,0}_{N-k-1\, {{\text {terms}}}}), \end{aligned}$$

(A.5)

for $1 \le k \le N-1,$ whereas the remaining $N^2-N$ non-diagonal generators are defined as

$$\begin{aligned} \left( T_{\tiny {\hbox {(Fund)}}}^{(i,j;1)}\right) _{a,b} = \frac{1}{2} \left( \delta _{a,i}\delta _{b,j} + \delta _{a,j}\delta _{b,i}\right) ,\end{aligned}$$

(A.6)

$$\begin{aligned} \left( T_{\tiny {\hbox {(Fund)}}}^{(i,j;2)}\right) _{a,b} = \frac{1}{2i} \left( \delta _{a,i}\delta _{b,j} - \delta _{a,j}\delta _{b,i}\right) , \end{aligned}$$

(A.7)

for $1 \le i < j \le N.$ Note that, with these conventions, for $N=2$ the generators in the fundamental representation reduce to $\frac{1}{2}\sigma ^a$ ($\sigma ^a$ denoting a Pauli matrix), while for $N=3$ the generators in the fundamental representation are $\frac{1}{2}\lambda ^a$ (where $\lambda ^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) \times (N(N-1)/2)$ whose entries (which can be labelled by pairs of indices (i, j),  with $1 \le i < j \le N$) are defined as

$$\begin{aligned} U_{(i,j) (k,l)} = -2 {{\,\textrm{Tr}\,}}\left( T_{\tiny {\hbox {(Fund)}}}^{(i,j;2)\, {\tiny {T}}}\, u\, T_{\tiny {\hbox {(Fund)}}}^{(k,l;2)}\, u^{\tiny {T}} \right) . \end{aligned}$$

(A.8)

Interpreting Eq. (A.8) as a map ${\mathscr {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 ${\mathscr {R}}_{\star }$ of the generators in the fundamental representation. Equivalently, if one considers group elements of the form $u=\mathbbm {1}_N + i \epsilon T_{\tiny {\hbox {(Fund)}}}^a + O(\epsilon ^2),$ the corresponding generator in the 2AS representation can be obtained as

$$\begin{aligned} \lim _{\epsilon \rightarrow 0} (-i/\epsilon )\left( U - \mathbbm {1}_{N(N-1)/2}\right) . \end{aligned}$$

(A.9)

Explicitly, with these conventions the SU(4) generators in the fundamental representation take the form:

$$\begin{aligned}{} & {} T^1_{\tiny {\hbox {(Fund)}}} = \frac{1}{2}\left( \begin{array}{cccc} 0 &{} 1 &{} 0 &{} 0 \\ 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^2_{\tiny {\hbox {(Fund)}}} = \frac{1}{2}\left( \begin{array}{cccc} 0 &{} -i &{} 0 &{} 0 \\ i &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\nonumber \\{} & {} T^3_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{ccccc} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^4_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\nonumber \\{} & {} T^5_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} -i &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ i &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^6_{\tiny {\hbox {(Fund)}}} = \frac{1}{2}\left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\nonumber \\{} & {} T^7_{\tiny {\hbox {(Fund)}}} = \frac{1}{2}\left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -i &{} 0 \\ 0 &{} i &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^8_{\tiny {\hbox {(Fund)}}} = \frac{1}{2\sqrt{3}} \left( \begin{array}{cccc} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} -2 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \end{array} \right) , \nonumber \\{} & {} T^9_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 1 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^{10}_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} -i \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ i &{} 0 &{} 0 &{} 0 \end{array} \right) , \nonumber \\{} & {} T^{11}_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 \end{array} \right) ,\quad T^{12}_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} -i \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} i &{} 0 &{} 0 \end{array} \right) , \nonumber \\{} & {} T^{13}_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 1 &{} 0 \end{array} \right) ,\quad T^{14}_{\tiny {\hbox {(Fund)}}} = \frac{1}{2} \left( \begin{array}{cccc} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} -i \\ 0 &{} 0 &{} i &{} 0 \end{array} \right) , \nonumber \\{} & {} T^{15}_{\tiny {\hbox {(Fund)}}} = \frac{1}{2\sqrt{6}} \left( \begin{array}{cccc} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} -3 \end{array} \right) , \end{aligned}$$

(A.10)

while the SU(4) generators in the two-index antisymmetric representation read:

$$\begin{aligned} T^1_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^2_{\tiny {\hbox {(2AS)}}} = \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -i &{} 0 &{} 0 &{} 0 \\ 0 &{} i &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} -i &{} 0 \\ 0 &{} 0 &{} 0 &{} i &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\end{aligned}$$

(A.11)

$$\begin{aligned} T^3_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^4_{\tiny {\hbox {(2AS)}}} = \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \end{array} \right) , \nonumber \\ T^5_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} i &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ -i &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -i \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} i &{} 0 &{} 0 \end{array} \right) ,\quad T^6_{\tiny {\hbox {(2AS)}}} = \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \end{array} \right) , \nonumber \\ T^7_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} -i &{} 0 &{} 0 &{} 0 &{} 0 \\ i &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -i \\ 0 &{} 0 &{} 0 &{} 0 &{} i &{} 0 \end{array} \right) ,\quad T^8_{\tiny {\hbox {(2AS)}}} = \frac{1}{\sqrt{12}}\left( \begin{array}{cccccc} 2 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -2 \end{array} \right) , \nonumber \\ T^9_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^{10}_{\tiny {\hbox {(2AS)}}} = \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} i &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} i \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ -i &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -i &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) , \nonumber \\ T^{11}_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 \\ 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^{12}_{\tiny {\hbox {(2AS)}}} = \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} -i &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} i \\ i &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -i &{} 0 &{} 0 &{} 0 \end{array} \right) , \nonumber \\ T^{13}_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) ,\quad T^{14}_{\tiny {\hbox {(2AS)}}} = \frac{1}{2}\left( \begin{array}{cccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} -i &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} -i &{} 0 \\ 0 &{} i &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} i &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) , \nonumber \\ T^{15}_{\tiny {\hbox {(2AS)}}}= & {} \frac{1}{\sqrt{6}}\left( \begin{array}{cccccc} 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 \end{array} \right) . \end{aligned}$$

(A.12)

Appendix 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 $\tau _{int}$ of the observable A as

$$\begin{aligned} \tau _{int} = \frac{1}{2}+ \sum _{t=1}^{W} \frac{\varGamma (t)}{\varGamma (0)}, \end{aligned}$$

(B.1)

where t is the Monte Carlo time, W is the summation window, and $\varGamma (t)$ is our approximation of the autocorrelation function of the series:

$$\begin{aligned}{} & {} \varGamma (t) \simeq \frac{1}{N_{cfg}-t} \sum _{n=1}^{N_{cfg}-t} \left( a_n - \langle a_-\right) \rangle \left( a_{n+t} - \langle a_+\right) \rangle , \nonumber \\{} & {} \quad 0 \le t < N_{cfg}. \end{aligned}$$

(B.2)

$\langle a_{-} \rangle$ and $\langle a_{+}\rangle$ are the averages over $N_{cfg}-t$ measurements, respectively the first and the last ones (see Section 4 of Ref. [36] for more details). In order to choose a summation window W of Eq. (B.1), we check at which Monte Carlo time t the function $\varGamma (t)/\varGamma (0)$ is smaller than twice its variance, which is estimated within the Madras Sokal approximation [36, 37].

Fig. 15

Top panel: Autocorrelation function $\varGamma (t)/\varGamma (0)$ computed 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. Bottom panel: Integrated autocorrelation time as a function of the summation window, computed for the same correlator on the top at several times. The two point function is obtained from the ensemble A0

Typical plots that we obtain for these quantities are shown in Fig. 15 for the pseudoscalar–pseudoscalar correlator of fundamental fermions at different lattice times. In the top panel, our estimation of the autocorrelation function $\varGamma (t)/\varGamma (0)$ is plotted as a function of Monte Carlo time, showing the typical decaying behaviour. The bottom panel shows the integrated autocorrelation time $\tau _{int}$ of Eq. (B.1) for an increasing summation window W. The error is computed in the Madras Sokal approximation,

$$\begin{aligned} \langle \delta \tau _{int}^2\simeq \rangle \frac{4W+2}{N_{cfg}} \, \tau _{int}^2. \end{aligned}$$

(B.3)

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 $\sigma _0$ of a given observable can be corrected to be $\sigma ^2=2\tau _{int} \sigma _0^2.$

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.

Appendix C: Measurements of correlation functions

In this section we discuss the details behind the measurement of correlation functions $C_{ab}^{R}(t) \equiv C_{ab}^{R}(t,0)$

$$\begin{aligned} C_{ab}^{R}(t,0) = \frac{1}{L^3} \sum _{\varvec{x}} \langle O^{R}_a(\varvec{x},t) O^{R^\dagger }_b(0,t_i),\rangle \quad a,b = P,A, \nonumber \\ \end{aligned}$$

(C.16)

where R denotes the representation of the fermionic fields. For each gauge configuration $n=1,\dots N_{cfg}$ we average measurements performed every 4 source-times $t_i,$ resulting in the correlator $C_{ab, \, n}^{R}(t),$

$$\begin{aligned} C_{ab, \, n}^{R} (t) = \frac{8}{T}\sum _{t_i=0,4,\dots }^{T/2} C_{ab,\, n}^{R}(t,t_i). \end{aligned}$$

(C.17)

The gauge average $C_{PP}^{R}(t)$ and its variance $\delta C_{ab}^{R}(t)^2$ are then estimated from a new set of N correlators, $C_{ab, \, n^{\prime }}^{R},$ resampled with a bootstrap procedure from $C_{ab, \, n}^{R}.$ In this step, trajectories are discarded in order to account for the autocorrelation that is monitored across observables and ensembles according to Appendix B.

$$\begin{aligned} C_{ab}^{R}(t)= & {} \frac{1}{N^{\prime }} \sum _{n^{\prime }=0}^{N^{\prime }-1} C_{ab, \, n^{\prime }}^{R}, \nonumber \\ \delta C_{ab}^{R}(t)^2= & {} \dfrac{1}{N^{\prime }} \sum _{n^{\prime }=0}^{N^{\prime }-1} \left[ C_{ab, \, n^{\prime }}^{R} - C_{ab}^{R}(t) \right] ^2 \end{aligned}$$

(C.18)

In order to suppress the excited states created by the hadronic operator in Eq. (C.16), we compute the two point functions by using different types of smearing on the fermionic fields. The procedure is not gauge invariant and requires working at fixed gauge. In particular, we adopt local $(\psi ^{R}(x))$ and smeared $({\tilde{\psi }}^{R}(x))$ fields

$$\begin{aligned}{} & {} \psi ^{R}(x)_{\alpha }^c = \int dy\, \delta (x-y) \delta _{\alpha \alpha ^{\prime }} \delta _{cc^{\prime }} \, \psi ^{R}(y)_{\alpha ^{\prime }}^{c^{\prime }};\nonumber \\{} & {} {\tilde{\psi }}_g^{R}(x)_{\alpha }^c = \int dy\, \frac{e^{-(x-y)^2/2g^2}}{\sqrt{2\pi }g} \delta _{\alpha \alpha ^{\prime }} \delta _{cc^{\prime }} \psi ^R(y)_{\alpha ^{\prime }}^{c^{\prime }}, \end{aligned}$$

(C.19)

where $c, c'$ are color indices, $\alpha , \alpha ^{\prime }$ are Dirac indices and a sum is intended over $\alpha ^{\prime }, c^{\prime }.$ The parameter g can is tuned according to Sect. 7.2 in order to suppress the excited states in the spectral reconstruction. By combining local and smeared operators, we obtain three types of correlators: local–local, smeared–smeared, local–smeared. In the latter, the operators at the source and the sink are different, and this can produce negative contributions to the spectral density. As demonstrated in the Appendix of [14], such terms do not jeopardise the extraction of the spectral density from the correlators.

Appendix D: 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,\dots 6$ are listed in this appendix. These are obtained according to the convention of [9] regarding the generators and the structure constants of SU(4).

$$\begin{aligned} X_1= & {} -\frac{i}{2} \begin{pmatrix}0 &{} -2 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 2 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{pmatrix}, \end{aligned}$$

(D.20)

$$\begin{aligned} X_2= & {} -\frac{i}{2} \begin{pmatrix}0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 &{} -\sqrt{3} &{} 0 &{} 0 &{} 0 &{} 0 \\ -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} \sqrt{3} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{pmatrix}, \end{aligned}$$

(D.21)

$$\begin{aligned} X_3= & {} -\frac{i}{2}\begin{pmatrix}0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} -\sqrt{3} &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \sqrt{3} &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{pmatrix}, \end{aligned}$$

(D.22)

$$\begin{aligned} X_4= & {} -\frac{i}{2}\begin{pmatrix}0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1/\sqrt{3} &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 &{} -1/\sqrt{3} &{} 0 &{} 0 &{} 0 &{} -2\sqrt{2/3} \\ -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 2\sqrt{2/3} &{} 0 &{} 0 &{} 0 \end{pmatrix}, \end{aligned}$$

(D.23)

$$\begin{aligned} X_5= & {} -\frac{i}{2}\begin{pmatrix}0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1/\sqrt{3} &{} 0 &{} 0 \\ -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} -1/\sqrt{3} &{} 0 &{} 0 &{} 0 &{} -2\sqrt{2/3} \\ 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 2\sqrt{2/3} &{} 0 &{} 0 \end{pmatrix}, \end{aligned}$$

(D.24)

$$\begin{aligned} X_6= & {} -\frac{i}{2} \begin{pmatrix}0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} -2/\sqrt{3} &{} 0 \\ 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 2/\sqrt{3} &{} 0 &{} 0 &{} 0 &{} -2\sqrt{2/3} \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 2\sqrt{2/3} &{} 0 \end{pmatrix}\end{aligned}$$

(D.25)

Appendix E: Ensembles

See Table 5.

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 $\beta =11,$ and the volume of the lattice is $16^3 \times 32$

Appendix F: Masses

See Tables 6, 7, 8 and 9.

Table 6 PCAC masses for the fundamental representation used in the chiral extrapolation. The bare mass of the antisymmetric fermions is fixed at $am_0^{(2AS)}=-0.45.$ They correspond to the ensembles A0–A4

Table 7 PCAC masses for the antisymmetric representation used in the chiral extrapolation. The values on the left are obtained with $am_0^{({\textrm{Fund}})}=-0.45$ (ensembles B0–B4), the ones on the right with $am_0^{({\textrm{Fund}})}=-0.47$ (ensembles C0–C4)

Table 8 Masses of the pseudoscalar mesons. On the left, we vary the bare fundamental mass while keeping $am_0^{(2{\textrm{AS}})}=-0.45$ (ensembles A0–A4). Conversely, on the right the fundamental bare mass is fixed at $am_0^{({\textrm{Fund}})}=-0.45$ (ensembles B0–B4)

Table 9 Fit results for $a M_{PP}^{(2AS)}$ from smeared spectral densities for different smearing radii $\sigma$ and different number of states n. These values appear in Fig. 13