1 Introduction

The experimental observation of a particle compatible with the Standard-Model Higgs boson at the Large Hadron Collider (LHC) [1, 2] and the lack of evidence of any New Physics are putting very tight constraints on theories beyond the Standard Model. Nevertheless, for all its shortcomings, it remains very hard to imagine that the Standard Model be the correct description of Nature up to energies much higher than the TeV scale.

An unsatisfactory aspect of the Standard Model is the fact that, among its parameters, it features a large number of Yukawa couplings, which cannot be derived from first principles, and which give rise to broadly separated masses for the fermions. Also in the fermionic sector, it does not account for the experimental evidence of neutrino oscillations [3, 4], implying that these particles are not massless (although it can be easily extended to accommodate massive neutrinos just by adding a handful further parameters, at least if they are Dirac particles). Even more remarkably, the Standard Model fails spectacularly at predicting \(95\%\) of the observed energy budget of the Universe [5], because it does not provide any explanation for Dark Matter or Dark Energy. Other unsatisfactory aspects of the Standard Model include the absence of unification of the gauge interactions, the “strong-\({{\mathcal {C}}}{\mathcal {P}}\) problem” of quantum chromodynamics (QCD), and the fact that it does not include a proper quantum formulation of gravity. Finally, as is well known, one of the major theoretical puzzles in the Standard Model is the lightness of the Higgs boson: being the only fundamental scalar in the theory, its mass receives contributions (of opposite signs) from quantum fluctuations at all energies up to the Planck scale, but their sum turns out to be surprisingly (“unnaturally”) small in comparison to the latter scale; for a recent review, see Ref. [6].

At least for the last of these issues, i.e. the “naturalness problem”, supersymmetry provides a conceptually very elegant solution: the (nearly) perfect cancellation of the contributions to the Higgs boson mass from quantum fluctuations of different fields is a consequence of the (only softly broken) symmetry relating bosonic and fermionic species in the theory. From a formal point of view, it is also worth remarking that supersymmetry is the only type of symmetry combining spacetime and internal degrees of freedom in a non-trivial way [7], evading the Coleman–Mandula theorem [8], and its experimental observation in elementary particle physics would be a major scientific discovery. In practice, however, its simplest realization in a framework compatible with the particle content of the Standard Model (the minimal supersymmetric Standard Model, MSSM), in which supersymmetry is necessarily broken, is far less aesthetically appealing: in particular, the MSSM has more than a hundred fundamental parameters, which, like their analogues in the non-supersymmetric Standard Model, cannot be derived from first principles. Despite the lack of predictive power due to this large number of free parameters, the MSSM (like most other New Physics models) generically predicts the existence of a host of new particles, including, in particular, four further Higgs particles, in addition to the Standard-Model one. All experimental searches in this direction so far, however, have come away empty-handed, indicating that supersymmetry, if exists, probably lies at an energy scale out of the reach of current accelerators.

Another popular theoretical framework that could explain the small mass of the Higgs boson is the one in which this particle is not considered as elementary, but rather as a composite state of some new, strongly coupled, elementary degrees of freedom, so that its lightness could be interpreted in terms of a Nambu–Goldstone mechanics – much like the pion, the lightest physical state in the QCD spectrum, is (nearly) massless because it can be interpreted as the Nambu–Goldstone boson associated with the breaking of chiral symmetry. This idea, dating back to more than thirty years ago [9, 10], has been studied in a large number of works [11]: the simplest models realizing this scenario can be constrained by severe phenomenological tests [12] and have been falsified by now, but more refined implementations of this idea remain theoretically attractive and could still be viable candidates for New Physics beyond the Standard Model.

Partial-compositeness models, in which some additional fermionic fields from this new strongly coupled sector are linearly coupled to the top quark, are particularly appealing [13, 14]. In this respect, a systematic, group-theoretical classification of the four-dimensional fermionic gauge theories providing an ultraviolet (UV) completion of composite-Higgs models was presented in Ref. [15], imposing the requirements related to the existence of a custodial symmetry, and the presence of top-quark partners. The simplest UV-complete model of this type was then discussed in Ref. [16]: it is a theory based on local invariance under an \(\mathrm {SU}(4)\) “hypercolor” group, featuring five flavors of massless Majorana fermions in the two-index antisymmetric representation, and three flavors of Dirac fermions in the fundamental representation of the gauge group. In the infrared limit, the formation of a condensate for the Majorana fermions in the two-index antisymmetric representation induces dynamical chiral-symmetry breaking according to the pattern \(\mathrm {SU}(5) \rightarrow \mathrm {SO}(5)\), and a composite state, embodying the Standard-Model Higgs boson doublet, arises then from the \(\mathrm {SU}(5)/\mathrm {SO}(5)\) coset [17]. The Dirac fermions in the fundamental representation bind with the Majorana fermions to form hypercolor-singlet states, that are interpreted as partners of the top quark, whereas the other massive Standard-Model fermions acquire their masses from quadratic coupling to the Higgs. This theory does not violate current experimental bounds e.g. on the decays of the Z boson, and is a viable UV-complete model for New Physics.

Since the crucial phenomena of chiral-symmetry breaking and hypercolor confinement in the model proposed in Ref. [16] are intrinsically non-perturbative in nature, a theoretical study of this theory from first principles requires lattice calculations. For technical reasons (related to the computational cost of the fermionic-matter content of the theory), however, it is more convenient to study first a closely related theory, with two flavors of Dirac fermions in the two-index antisymmetric representation of \(\mathrm {SU}(4)\), and two flavors of Dirac fermions in the fundamental representation of the gauge group. With such matter contents, the theory will undergo a different symmetry-breaking pattern (in particular, one which can not accommodate a state with quantum numbers compatible with those of the Standard-Model Higgs boson); nevertheless, it remains an interesting theoretical laboratory, in which the main features of the actual model discussed in Ref. [16] can be studied, at least at a qualitative or semi-quantitative level.

With this motivation, in the present work we present a detailed numerical investigation of the \(\mathrm {SU}(4)\) lattice gauge theory with two flavors of Dirac fermions in the two-index antisymmetric representation and two flavors of Dirac fermions in the fundamental representation of the gauge group, which recently has also been studied in a series of works [18,19,20,21]. The structure of this article is the following: in Sect. 2, we review the main features of the Ferretti model, in Sect. 3, we analyze in detail the symmetries of the Dirac operator (both in the continuum and in various lattice discretizations) in the two-index antisymmetric representation, and their implications for the spectrum supported by random matrix theory expectations. Next, in Sect. 4 we discuss the features of a hybrid Monte Carlo algorithm working with fermions in different representations, and in Sect. 5 we present our results, both as algorithmic checks and as first exploratory steps into the theory under consideration. Section 6 deals with the generalization of this type of studies to non-minimal partial. The concluding Sect. 7 presents a summary of this work, while the Appendices AB, and C respectively include our conventions for notations, the detailed proofs of some identities discussed in Sect. 3. and technical details about our hybrid Monte Carlo algorithm.

2 Overview of the model

Let us briefly review the model described in Ref. [16] which we refer to as “Ferretti model”. The UV completion is a gauge theory with \(G_{\mathrm{HC}}=\mathrm {SU}(4)\) “hypercolor” gauge group, coupled to five Weyl fermions \(\psi ^I_{mn}\) in the two-index antisymmetric representation of the hypercolor group (i.e. the dimension \(\mathbf {6}\) representation, that, in the following, we also call “sextet” representation: for a summary of group and group-representation properties, see, for instance, Ref. [22, appendix]) and three Dirac fermions written in terms of pair of Weyl fermions \(\chi ^a_m\),\({\bar{\chi }}^{a'}_m\) in the fundamental representation of the hypercolor group. Hence, in the field definition the indices \(I,a,a'\) run over the flavor and read respectively \(I=1,\ldots ,5\), whereas \(a,a'=1,\ldots ,3\); on the other hand, \(m,n = 1,\ldots ,4\) denote hypercolor indices. The global internal symmetry of the theory is

$$\begin{aligned} G_{\mathrm{F}}=\mathrm {SU}(5) \times \mathrm {SU}(3) \times \mathrm {SU}(3)' \times \mathrm {U}(1)_X \times \mathrm {U}(1)' \, . \end{aligned}$$
(1)

The charges of the various fields are listed in Table 1.

Table 1 \(G_{\mathrm{HC}}\) is the hypercolor gauge group and \(G_{\mathrm{F}}\) the global symmetry group before symmetry breaking

The symmetry-breaking pattern of the model can be described as

$$\begin{aligned} G_{\mathrm{F}}/H_{\mathrm{F}}= & {} \left( \frac{\mathrm {SU}(5)}{\mathrm {SO}(5)}\right) \times \left( \frac{\mathrm {SU}(3)\times \mathrm {SU}(3)'}{\mathrm {SU}(3)_c} \right) \nonumber \\&\times \left( \frac{\mathrm {U}(1) \times \mathrm {U}(1)'}{\mathrm {U}(1)_X} \right) \, , \end{aligned}$$
(2)

and is realized by the bilinear fermionic condensates \(\langle \epsilon ^{mnpq}\psi ^I_{mn}\psi ^J_{pq}\rangle \propto \delta ^{IJ}\), and \(\langle {\bar{\chi }}^{a'}_m \chi ^a_m \rangle \propto \delta ^{a' a}\). The symmetry-breaking pattern \(G_{\mathrm{F}}/H_{\mathrm{F}}\) is compatible with a custodial symmetry, described by the \(G_{\mathrm{cus}}\) group, such that \(H_{\mathrm{F}} \supset G_{\mathrm{cus}} \supset G_{\mathrm{SM}}\), with \(G_{\mathrm{cus}}=\mathrm {SU}(3)_c\times \mathrm {SU}(2)_L \times \mathrm {SU}(2)_R \times \mathrm {U}(1)_X\) and \(G_{\mathrm{SM}} = \mathrm {SU}(3)_c \times \mathrm {SU}(2)_L \times \mathrm {U}(1)_Y\) is the Standard Model gauge group. More in detail, the electroweak gauge group \(\mathrm {SU}(2)_L \times U(1)_Y\) is embedded in the unbroken \(\mathrm {SO}(5)\), by considering the subgroup \(\mathrm {SO}(4) \simeq \mathrm {SU}(2)_L\times \mathrm {SU}(2)_R\), identifying \(\mathrm {U}(1)_R\) as the subgroup of \(\mathrm {SU}(2)_R\) generated by the third generator \(T^3_R\), and setting the hypercharge \(Y=T^3_R + X\). The 14 Nambu–Goldstone bosons in the \(\mathrm {SU}(5)/\mathrm {SO}(5)\) coset can be classified according to their SM \(\mathrm {SU}(2)_L \times U(1)_R\) charges as:

$$\begin{aligned} \mathbf{14} \rightarrow \mathbf{1}_0 + \mathbf{2}_{\pm 1/2} + \mathbf{3}_0 + \mathbf{3}_{\pm 1} = (\eta , H ,\phi _0,\phi _{\pm })\, , \end{aligned}$$
(3)

where the field H can be interpreted as the Higgs field. Indeed this field is a doublet under \(\mathrm {SU}(2)_L\) and can therefore be written as a two-component complex field \(H=(H_+,H_0)\). The spin-1 / 2 states appear as a triplet of the hypercolor theory, and are natural candidates to play the rôle of top-quark partners: in the effective field theory description of the low-energy dynamics, the latter are introduced as a Dirac fermion field \(\Psi \) transforming in the \((\mathbf{5},\mathbf{3})_{2/3}\) of \(H_{\mathrm{F}}\). Such a field can be realized within the Standard Model, by decomposing the \((\mathbf{5},\mathbf{3})_{2/3}\) multiplet as

$$\begin{aligned} (\mathbf{5},\mathbf{3})_{2/3} \rightarrow (\mathbf{3},\mathbf{2})_{7/6} + (\mathbf{3},\mathbf{2})_{1/6} + (\mathbf{3},\mathbf{1})_{2/3} \, , \end{aligned}$$
(4)

where the numbers on the right-hand side label the irreducible representations of \(G_{\mathrm{SM}}\). The Nambu–Goldstone bosons can be combined into a \(\Pi \) field

$$\begin{aligned} \Pi = H + H^{\dag } + \phi _0 + \phi _+ + \phi ^{\dag }_+ , \end{aligned}$$
(5)

from which one can define

$$\begin{aligned} \Sigma = \exp \left( \frac{i\Pi }{f} \right) , \end{aligned}$$
(6)

with \(\Pi \) a real symmetric matrix. The matrix \(\Sigma \) defined in Eq. (6), however, transforms non-linearly under a transformation \(g \in \mathrm {SO}(5)\), so it is convenient to consider the field \(U=\Sigma \Sigma ^T=\exp ( 2i \Pi / f)\), which transforms linearly: \(U \rightarrow g U g^T\).

The couplings to vector bosons are obtained from the chiral Lagrangian

$$\begin{aligned} {\mathcal {L}} \supset \frac{f^2}{16}{{\,\mathrm{Tr}\,}}\left( (D_{\mu }U)^{\dag }D^{\mu }U \right) \end{aligned}$$
(7)

where the derivative is promoted to the covariant derivative, i.e.

$$\begin{aligned} D_{\mu } U = \partial _{\mu } U - igW^a_{\mu }[T^a_L,U] - ig'B_{\mu }[T^3_R,U]\, . \end{aligned}$$
(8)

The mass term for fermions is

$$\begin{aligned} {\mathcal {L}} \supset M{\bar{\Psi }}\Psi + \lambda _q f \bar{{\hat{q}}}_L \Sigma \Psi _R + \lambda _t f \bar{{\hat{t}}}_R\Sigma ^* \Psi _L\, , \end{aligned}$$
(9)

where \({\hat{q}}_L\) and \({\hat{t}}_R\) are the spurionic embedding of the SM quarks in the \(\mathbf{5}\) and \(\bar{\mathbf{5}}\) representations of \(\mathrm {SU}(5)\), respectively.

An important feature of such a model is the vacuum misalignment, which is responsible for electro-weak symmetry breaking. In particular, the SM fermionic couplings are responsible for negative contributions to the Coleman–Weinberg potential, which are necessary to generate a non-vanishing vacuum expectation value for the \(H_0\) component. Following Ref. [16], we set \(H_0 = h/\sqrt{2}\), while all other fields are set to zero. Then, the coupling of the h field to the SM gauge bosons and fermions reads

$$\begin{aligned} {{\,\mathrm{Tr}\,}}\left[ U(h)W_{\mu }U(h)^{\dag }W_{\mu } \right]&= \frac{1}{2} [1+\cos (2h/f)]W^c_{\mu }W^c_{\mu }\, , \end{aligned}$$
(10)
$$\begin{aligned} \bar{{\hat{q}}}_L U(h) {\hat{t}}_R + \bar{{\hat{t}}}_R U(h)^* {\hat{q}}_L&= \frac{1}{\sqrt{2}}\sin (2h/f)({\bar{t}}_Lt_R + {\bar{t}}_R t_L)\, . \end{aligned}$$
(11)

The potential can thus be parametrized by the two low-energy constants \(\alpha \) and \(\beta \) as

$$\begin{aligned} V(g) \propto \alpha \cos \left( \frac{2h}{f}\right) - \beta \sin ^2\left( \frac{2h}{f}\right) \end{aligned}$$
(12)

and a suitable electro-weak-breaking minimum can be obtained at \(\cos ( 2 \langle h \rangle /f ) = -\alpha /(2\beta )\) for \(|\alpha /\beta | \lesssim 2\). These two constants can be computed as described in Ref. [23]. In particular, one has

$$\begin{aligned} 2\beta = - y^2 C_{\mathrm{top}} \, , \end{aligned}$$
(13)

which, in principle can be computed on the lattice, as well as all the other low-energy constants relevant for the infrared physics of the theory.

3 Symmetries of the Dirac operator for fermions in the sextet representation

In order to construct the two-index antisymmetric representation for a generic \(\mathrm {SU}(N)\) group, we introduce a set \(\left\{ e^{(a,b)} \right\} \) of \(N(N-1)/2\) real and antisymmetric matrices of size \(N \times N\), which we label by strictly increasing pairs of indices \(1 \le a < b \le N\). We sort the set of (ab) pairs starting from \(a=1\) and \(b=2\), and then increasing b and letting a run from 1 to \(b-1\), so that the sorted list of (ab) pairs reads (1, 2), (1, 3), (2, 3), (1, 4), (2, 4), (3, 4), ..., \((N-1,N)\). The elements of the \(e^{(a,b)}\) matrices are defined by

$$\begin{aligned} \left( e^{(a,b)} \right) _{p q} = \frac{1}{\sqrt{2}} ( \delta _{p,a}\delta _{q,b} - \delta _{p,b}\delta _{q,a}). \end{aligned}$$
(14)

Then, given a generic element u of the \(\mathrm {SU}(N)\) group in the fundamental representation, the corresponding group element in the two-index antisymmetric representation is a complex-valued matrix of size \((N(N-1)/2) \times (N(N-1)/2)\), whose entries are defined as

$$\begin{aligned} U_{(a,b) (c,d)} = {{\,\mathrm{Tr}\,}}\left( e^{(a,b)\, {\tiny {T}}}\, u\, e^{(c,d)}\, {u}^{{\tiny {T}}} \right) . \end{aligned}$$
(15)

It is then trivial to work out the explicit form of an arbitrary generator in the two-index antisymmetric representation, that we denote as \(T^a_{\mathrm{2AS}}\), for example, by defining an infinitesimal real parameter \(\epsilon \), taking u to be the group element infinitesimally close to the \(N \times N\) identity matrix \(u=\mathbb {1}+ i \epsilon t^a + O(\epsilon ^2)\), and extracting the components of \(T^a_{\mathrm{2AS}}\) as the coefficients of the terms linear in \(i \epsilon \) in the resulting expression for \(U - \mathbb {1}\) (where now \(\mathbb {1}\) denotes the \((N(N-1)/2) \times (N(N-1)/2)\) identity matrix).

For the purposes of this work, let us focus on the \(\mathrm {SU}(4)\) group, whose generators in both the fundamental and in the two-index antisymmetric representation are reported in Appendix A. Consider now the totally antisymmetric four-index tensor \(\epsilon _{abcd}\), with \(\epsilon _{1\,2\,3\,4}=1\). Interpreting its indices pairwise, it can be used to construct a \(6 \times 6\) matrix W, acting on the antisymmetric two-index representation of the \(\mathrm {SU}(4)\) generators, whose rows (and columns) are labelled by the sorted (ab) (and (cd)) pairs introduced above. The elements of W are defined as

$$\begin{aligned} W_{(a,b)\, (c,d)} = \epsilon _{abcd}. \end{aligned}$$
(16)

Remembering that, in our conventions, the indices from 1 to 6 of the antisymmetric two-index representation of \(\mathrm {SU}(4)\) are associated with the sorted pairs (1, 2), (1, 3), (2, 3), (1, 4), (2, 4), (3, 4), in that order, W takes the form

$$\begin{aligned} W = \left( \begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 &{} -1 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} -1 &{} 0 &{} 0 &{} 0 &{} 0 \\ 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0 \end{array} \right) . \end{aligned}$$
(17)

Note that W is real, symmetric, and unitary, hence it squares to the identity matrix. It is easy to check that all generators in the antisymmetric two-index representation of \(\mathrm {SU}(4)\) satisfy

$$\begin{aligned} W^{-1} T^a_{\mathrm{2AS}} W = -KT^a_{\mathrm{2AS}}, \end{aligned}$$
(18)

where K denotes the complex-conjugation operator, defined by \(K\alpha =\alpha ^*\) for every \(\alpha \in {\mathbb {C}}\).

Having set our notations for the generators of the \(\mathrm {SU}(4)\) algebra in their antisymmetric two-index representation and the \(\gamma \) matrices (for their explicit forms, see Appendix A), let us now introduce the Euclidean Dirac operator for a fermionic Dirac field of (real) bare mass m, transforming under the antisymmetric two-index color representation in a theory with \(\mathrm {SU}(4)\) gauge symmetry. In the continuum, the Euclidean Dirac operator reads:

$$\begin{aligned} D_{{\tiny \mathrm{cont}}}= \gamma _\mu D_\mu + m = \gamma _\mu \left( \partial _\mu +ig A_\mu ^a T^a_{\mathrm{2AS}} \right) + m. \end{aligned}$$
(19)

Note that the kinetic (\(\gamma _\mu D_\mu \)) part of \(D_{{\tiny \mathrm{cont}}}\) is anti-Hermitian, whereas the mass term m is Hermitian, so that, in general, \(D_{{\tiny \mathrm{cont}}}\) is neither Hermitian, nor anti-Hermitian. However, the anti-commutation relations \(\left\{ \gamma _5, \gamma _\mu \right\} =0\) imply that the \(\gamma _5 D_{{\tiny \mathrm{cont}}}\) operator is Hermitian:

$$\begin{aligned} \left( \gamma _5 D_{{\tiny \mathrm{cont}}}\right) ^\dagger= & {} D_{{\tiny \mathrm{cont}}}^\dagger \gamma _5^\dagger = \left( -\gamma _\mu D_\mu +m \right) \gamma _5 \nonumber \\= & {} \gamma _5 \left( \gamma _\mu D_\mu +m \right) = \gamma _5 D_{{\tiny \mathrm{cont}}}. \end{aligned}$$
(20)

Let us introduce the notion of “anti-unitary operator”: given a complex Hilbert space \({\mathcal {H}}\) with inner product \(\langle \dots , \dots \rangle \), an invertible mapping

$$\begin{aligned} {\mathcal {U}} : {\mathcal {H}} \rightarrow {\mathcal {H}}, \qquad \phi \rightarrow {\mathcal {U}}(\phi ) \end{aligned}$$
(21)

(where \(\phi \) denotes an arbitrary element of \({\mathcal {H}}\)) is said to be “anti-unitary” if it is antilinear

$$\begin{aligned} {\mathcal {U}}\left( \alpha \phi + \beta \rho \right) = \alpha ^\star {\mathcal {U}}(\phi ) +\beta ^\star {\mathcal {U}}(\rho ) \end{aligned}$$
(22)

and satisfies

$$\begin{aligned} \langle {\mathcal {U}}(\phi ), {\mathcal {U}}(\rho ) \rangle = \langle \phi , \rho \rangle ^\star \end{aligned}$$
(23)

for every \(\phi \) and \(\rho \) in \({\mathcal {H}}\) and for every a and b in \({\mathbb {C}}\). It is possible to prove that, given a unitary operator \({\mathcal {V}}\), the \({\mathcal {V}}K\) operator is anti-unitary, and that, conversely, every anti-unitary operator \({\mathcal {U}}\) can be written as

$$\begin{aligned} {\mathcal {U}} = {\mathcal {V}} K, \end{aligned}$$
(24)

where \({\mathcal {V}}\) is a unitary operator.

Let us introduce the charge conjugation \({\mathcal {C}}\) and define the operator A as

$$\begin{aligned} A = W {\mathcal {C}} \gamma _5 K. \end{aligned}$$
(25)

The combination \(W {\mathcal {C}} \gamma _5\) appearing in Eq. (25) is a unitary operator, so it follows from Eq. (24) that A is anti-unitary. Moreover, it is trivial to show that A squares to minus the identity, because

$$\begin{aligned} A^2= & {} W {\mathcal {C}} \gamma _5 K W {\mathcal {C}} \gamma _5 K = W {\mathcal {C}} \gamma _5 W^\star {\mathcal {C}}^\star \gamma _5^\star \nonumber \\= & {} W^2 {\mathcal {C}}^2 \gamma _5^2 = - \mathbb {1}, \end{aligned}$$
(26)

having used the facts that W (acting only on the color indices) commutes with \({\mathcal {C}}\) and \(\gamma _5\) (which act only on the spinor indices), that W, \({\mathcal {C}}\) and \(\gamma _5\) are real, that \({\mathcal {C}}\) commutes with \(\gamma _5\), and that W, \(\gamma _5\) and K square to the identity, whereas \({\mathcal {C}}\) squares to minus the identity.

From the aforementioned properties of W, \({\mathcal {C}}\), \(\gamma _5\), and A it also follows that

$$\begin{aligned}{}[A, \gamma _5 D_{{\tiny \mathrm{cont}}}]= & {} 0. \end{aligned}$$
(27)

A detailed proof of the above relation is provided in Appendix B.

Now, let us introduce the Dirac operator for the lattice discretization of the theory with fermions in the antisymmetric two-index representation, on a hypercubic spacetime lattice of spacing a. Its matrix elements in the Wilson formulationFootnote 1 are of the form

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

Thus, one also has:

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

Defining the four, unitary, “positive-shift” operators \(P_\mu \), that act trivially on all internal degrees of freedom and have real matrix elements between sites x and y given by

$$\begin{aligned} (P_\mu )_{x,y}=\delta _{x+a{\hat{\mu }},y}, \end{aligned}$$
(30)

(while their inverses have elements \((P_\mu )^{-1}_{x,y}=\delta _{x-a{\hat{\mu }},y}\)), and the local “positively-oriented, parallel-transporter” operators \(U_\mu \) (having matrix elements \(U_\mu (x) \delta _{x,y}\) between sites x and y), the Wilson Dirac operator can be written as

$$\begin{aligned} D= & {} \frac{1}{a} \left\{ \mathbb {1}- \kappa \sum _{\mu =1}^4 \left[ ( \mathbb {1}- \gamma _\mu ) (U_\mu P_\mu ) \right. \right. \nonumber \\&\left. \left. + ( \mathbb {1}+ \gamma _\mu ) (U_\mu P_\mu )^\dagger \right] \right\} . \end{aligned}$$
(31)

We now prove that the \(\gamma _5 D\) lattice operator commutes with A, exactly as its continuum counterpart \(\gamma _5 D_{{\tiny \mathrm{cont}}}\) does. In order to prove this statement, we first study the transformation properties of the \(U_\mu (x)\) link variables under complex conjugation. When D is the Wilson Dirac operator for fermions in the antisymmetric two-index representation, a generic link variable \(U_\mu (x)\) can be written as the exponential of i times a linear combination with real coefficients (that is convenient to write as \(agA_\mu ^a(x)\)) of the \(T^a_{\mathrm{2AS}}\) generators defined by Eq. (15) and explicitly reported in Appendix A:

$$\begin{aligned} U_\mu (x)= & {} \exp \left( i a g A_\mu ^a(x) T^a_{\mathrm{2AS}} \right) \nonumber \\= & {} \sum _{n=0}^\infty \frac{(iag)^n}{n!} \left( A_\mu ^a(x) T^a_{\mathrm{2AS}} \right) ^n. \end{aligned}$$
(32)

As a consequence:

$$\begin{aligned} K U_\mu (x)= & {} U_\mu ^\star (x) = \exp \left( -i a g A_\mu ^a(x) T^a_{\mathrm{2AS}} \right) \nonumber \\= & {} \sum _{n=0}^\infty \frac{(-iag)^n}{n!} \left( A_\mu ^a(x) T^{a \star }_{\mathrm{2AS}} \right) ^n. \end{aligned}$$
(33)

Using Eq. (18), the latter equation can be rewritten as

$$\begin{aligned} K U_\mu (x)= & {} \sum _{n=0}^\infty \frac{(iag)^n}{n!} \left( - A_\mu ^a(x) W^{-1} T^a_{\mathrm{2AS}} W \right) ^n \nonumber \\= & {} W^{-1} \left\{ \sum _{n=0}^\infty \frac{(iag)^n}{n!} \left( A_\mu ^a(x) T^a_{\mathrm{2AS}} \right) ^n \right\} W \nonumber \\= & {} W^{-1} U_\mu (x) W. \end{aligned}$$
(34)

From the transpose of the latter identity, using the fact that W is symmetric and equal to its inverse, it follows that \(U_\mu ^\dagger (x)=W^{-1} \left( U_\mu ^\dagger (x) \right) ^\star W\), so that \(K U_\mu ^\dagger (x) = W^{-1} U_\mu ^\dagger (x) W\). This is actually a trivial implication of Eq. (34), since \(U_\mu (x)\) is an element of a unitary group, so \(U_\mu ^\dagger (x)=U_\mu ^{-1}(x)\) is still a group element, in the same representation. As a consequence, the Wilson Dirac operator D is such that

$$\begin{aligned}{}[A, \gamma _5 D]= & {} 0. \end{aligned}$$
(35)

with \(A^2=-\mathbb {1}\): this is a property that the Wilson Dirac operator shares with the continuum Dirac operator. A detailed proof of Eq. (35) is provided in Appendix B.

Equation (35) implies that \(\gamma _5 D\) can always be rewritten as a matrix whose elements are real quaternions of the form

$$\begin{aligned} q_0 + i \vec {\sigma } \cdot \vec {q}, \end{aligned}$$
(36)

where \(q_0\) and the components of \(\vec {q}\) are real. As a consequence, the eigenvalues of \(\gamma _5 D\) are pairwise-degenerate.

A second, more interesting, consequence is that certain universal features of the spectrum of eigenvalues of \(\gamma _5 D\) can be described by the chiral Gaußian symplectic ensemble (chSE) in random matrix theory – see Ref. [24] for an excellent review. In particular, the unfolded density of spacings s between subsequent eigenvalues of \(\gamma _5 D\) is expected to follow the Wigner surmise

$$\begin{aligned}&P(s) = N_\beta s^\beta \exp ({-} c_\beta s^2),\nonumber \\&\text{ with }~N_\beta = 2 \frac{\Gamma ^{\beta +1}\left( \frac{\beta }{2}+1\right) }{\Gamma ^{\beta +2}\left( \frac{\beta +1}{2}\right) }, \quad \nonumber \\&c_\beta = \frac{\Gamma ^{2}\left( \frac{\beta }{2}+1\right) }{\Gamma ^{2}\left( \frac{\beta +1}{2}\right) }, \end{aligned}$$
(37)

for the Dyson index corresponding to the symplectic ensemble, \(\beta =4\). This is expected to hold for the unfolded density of spacings, in which the spacing between subsequent eigenvalues of \(\gamma _5 D\) in one gauge-field configuration is rescaled by the local spectral density (obtained from an average over all configurations).

Note that, for the continuum and Wilson Dirac operators for fundamental \(\mathrm {SU}(4)\) fermions, no global anti-unitary symmetry like the one encoded in Eq. (35) exists. As a consequence, the unfolded density of spacings between eigenvalues of the Wilson Dirac operator for fermions in the fundamental representation of the \(\mathrm {SU}(4)\) gauge group is expected to be described by the Wigner surmise for the chUE, i.e. by Eq. (37) with Dyson index \(\beta =2\).

In passing, we also note that in the staggered formulation of the lattice Dirac operator \(D_{{\tiny {\mathrm{st}}}}\) defined as

$$\begin{aligned} D_{{\tiny {\mathrm{st}}}}= m \mathbb {1}+ \frac{1}{2a} \sum _{\mu =1}^4 \eta _\mu \left[ (U_\mu P_\mu ) - (U_\mu P_\mu )^\dagger \right] , \end{aligned}$$
(38)

where \(\eta _\mu \) has elements between sites x and y defined as

$$\begin{aligned} (\eta _\mu )_{x,y} = \delta _{x,y} (-1)^{\sum _{\nu <\mu } x_{\nu }}, \end{aligned}$$
(39)

and where \(\gamma _5\) is replaced by \(\epsilon \), having elements

$$\begin{aligned} \epsilon _{x,y} = \delta _{x,y} (-1)^{\sum _{\mu =1}^4 x_{\mu }}, \end{aligned}$$
(40)

the analogue of \(\gamma _5 D_{{\tiny \mathrm{cont}}}\) is

$$\begin{aligned} \epsilon D_{{\tiny {\mathrm{st}}}}= m \epsilon + \frac{1}{2a} \sum _{\mu =1}^4 \epsilon \eta _\mu \left[ (U_\mu P_\mu ) - (U_\mu P_\mu )^\dagger \right] . \end{aligned}$$
(41)

Now, consider the antiunitary operator

$$\begin{aligned} B = WK, \end{aligned}$$
(42)

which squares to the identity:

$$\begin{aligned} B^2 = (WK)^2 = W^2 = \mathbb {1}. \end{aligned}$$
(43)

Analogously to the continuum and Wilson formulation, also in this case it is possible to show that

$$\begin{aligned}{}[B, \epsilon D_{{\tiny {\mathrm{st}}}}] =0. \end{aligned}$$
(44)

As a consequence of the above relation (whose demonstration is provided in Appendix B), the staggered Dirac operator \(D_{{\tiny {\mathrm{st}}}}\) is such that \(\epsilon D_{{\tiny {\mathrm{st}}}}\) commutes with the antiunitary operator B, which squares to \(\mathbb {1}\). This property implies that \(\epsilon D_{{\tiny {\mathrm{st}}}}\) can always be rewritten as a matrix whose elements are real, and that its universal spectral properties are described in terms of the chiral Gaußian orthogonal ensemble (chOE) of random matrix theory. In particular, the unfolded eigenvalue spacing distribution is expected to be approximated by the Wigner surmise defined in Eq. (37), but with \(\beta =1\), instead of 4 (as for the continuum and Wilson Dirac operators). This difference between the anti-unitary symmetries of the staggered and the continuum Dirac operators is, in fact, unsurprising, given that a similar situation also occurs for the \(\mathrm {SU}(2)\) gauge group [25, 26], and the convergence of the staggered-spectrum results to the correct continuum limit occurs in a subtle way [27]. The investigation of the restoration of the continuum symmetry in the staggered discretization of fermions in the sextet representation of the \(\mathrm {SU}(4)\) group for \(a \rightarrow 0\), however, would require a dedicated investigation and lies clearly beyond the scope of the present study.

4 Lattice-calculation setup

The simulations for this project were performed using a hybrid Monte Carlo (HMC) algorithm implemented with the GRID lattice QCD library [28]. As discussed above, given the exploratory nature of this work, we considered an approximation of the Ferretti model, reducing its matter content down to two fundamental and two sextet fermions. This prescription greatly simplifies the computational cost of the theory allowing to use a two-flavor pseudofermion action in the two representations. While this matter content does not yield the same symmetry breaking pattern as in the original model, this theory still represents an interesting theoretical framework with rich non-perturbative dynamics, analogous to the one proposed in Ref. [16]. Moreover, the simulation code we developed admits a rational hybrid Monte Carlo implementation that allows to simulate any number of dynamical flavors in a generic representation.

As in a standard HMC algorithm, the main steps are the following:

  1. 1.

    generation of pseudofermion fields from a heat-bath distribution;

  2. 2.

    dynamical evolution of the gauge field configuration according to a fictitious Hamiltonian with randomly chosen initial momenta for each link;

  3. 3.

    “accept-reject” step, to correct for possible errors in the integration of the equation of motion of the previous step.

While several sophisticated techniques can considerably improve the algorithmic performance (in particular for the inversion of the Dirac operator), for the purposes of this work we limited ourself to a conjugate gradient solver, without preconditioning. Simulations of the theory on a much larger scale would, of course, require a careful optimization of the setup, which is not discussed in this work.

4.1 HMC with fermions in multiple representations

Earlier works addressing the simulation of gauge theories with dynamical fermions in a generic representation include Ref.  [29] and subsequent publications by the same authors. The exploration of models with fermions in multiple representations, however, is still ongoing, mainly due to the variety of different models and to the computational effort their study requires. Among recent works, we would like to mention Ref.  [18], which presents a substantial set of numerical results (including continuum and chiral extrapolations for various physical quantities), which are obtained from an algorithm similar to the one presented in detail in this section.

Let us write the gauge link variable defined in a generic representation R as

$$\begin{aligned} U_{\mu }(x)= \exp \left\{ i \alpha ^a(x)_{\mu } T_R^a \right\} . \end{aligned}$$
(45)

In order to define the molecular-dynamics (MD) force for both gauge and fermions, let us define the variation of the link variable as

$$\begin{aligned} \delta (U_{\mu }(x))= & {} \delta (\alpha _{\mu }(x))U_{\mu }(x) \quad \quad \text {with} \quad \quad \nonumber \\ \delta (\alpha _{\mu }(x))= & {} i\delta (\alpha _{\mu }^a(x))T_R^a \end{aligned}$$
(46)

and the conjugate momentum associated with each fundamental link

$$\begin{aligned} \pi (x,\mu )=i\pi ^a(x,\mu )T_F^a . \end{aligned}$$
(47)

Note that the full dependence on the representation is encoded into the generators \(T_R\), meaning that the algebra weights (i.e. the gauge field components) are the same in any representation of the gauge group. Generalizing the same idea as in Ref.  [29], we consider the following Hamiltonian:

$$\begin{aligned} H = H_{\pi } + H_g + \sum _R^{N_{\mathrm{rep}}} H^R_{f}, \end{aligned}$$
(48)

where

  • \(H_{\pi }\) is the kinetic contribution from the conjugate momenta associated with links in the fundamental representation,

  • \(H_g\) is the pure gauge contribution, also based on gauge fields in the fundamental representation, while

  • \(H^R_{f}\) is the fermionic contribution, which can be in an arbitrary representation.

In the present case, the latter is chosen to be \(H^F_f + H^{\mathrm{2AS}}_f\). These terms are formally defined in the same way, except that in \(H^F_f\) the links and the pseudofermion fields are in the fundamental representation, while in \(H^{\mathrm{2AS}}_f\) the same links are “promoted” to the two-index antisymmetric representation by Eq. (15), and the pseudofermions are generated by a different heat-bath distribution.Footnote 2

More in detail, the terms appearing in Eq. (48) are:

$$\begin{aligned} H_{\pi }&= \frac{1}{2} \sum _{x,\mu } \left( \pi (x,\mu )),\pi (x,\mu ) \right) = \frac{1}{2}T_F \sum _{x,a,\mu } \pi ^a(x,\mu )^2 \, , \end{aligned}$$
(49)
$$\begin{aligned} H_g&= S_g = \frac{\beta }{N_c} \sum _x \sum _{\mu < \nu } {{\,\mathrm{Re}\,}}{{\,\mathrm{Tr}\,}}\left( 1 - {\mathcal {P}}_{\mu \nu }(x) \right) \, , \end{aligned}$$
(50)
$$\begin{aligned} H^R_f&= S^R_f = \sum _x \phi ^{\dag ;R}(x)[D^{\dag ;R}\, D^R]^{-1}\phi ^R(x) \, . \end{aligned}$$
(51)

We emphasize that the superscript R means that Eq. (51) holds for an arbitrary representation R (Sanity checks of the algorithm are showed in Fig. 1, while in Table 2 the plaquette observable is showed to approach correctly the quenched results from Ref. [31] in the infinite mass limit).

Fig. 1
figure 1

Plaquette (right) and \(\Delta H\) (left) as a function of the MD trajectories for \(\beta =10.0\) and \(am_4=am_6=-\,0.55\)

Table 2 Benchmark comparison of the value of the average plaquette in the infinitely-heavy-fermion limit to the quenched results for \(\mathrm {SU}(4)\) from Ref. [31]

For this project we consider the discretized Dirac operator D (dropping the superscript R) as the Wilson operator with the \({\mathcal {O}}(a)\) clover improvement and bare fermion mass m (in unit of lattice spacing):

$$\begin{aligned} D = D_{\mathrm{Wilson}} + D_{\mathrm{clover}}, \end{aligned}$$
(52)

where the matrix element of the Wilson operator has been already introduced in Eq. (28) and the improvement term reads

$$\begin{aligned}{}[D_{\mathrm{clover}}^R]_{xy}&= \frac{ia}{2}c_{\mathrm{sw}}(g_0^2)\kappa ^R \sum _{\mu ,\nu } {\tilde{F}}^R_{\mu \nu }(x)\sigma _{\mu \nu }\delta _{xy} \, . \end{aligned}$$
(53)

Let us express the fermion masses in terms of the hopping parameter

$$\begin{aligned} \kappa&=\frac{1}{2(am_0 + 4)} \quad \quad \text {and} \quad \quad am^R_q=am_0-am^R_{\mathrm{crit}}\nonumber \\&=\frac{1}{2}\left( \frac{1}{\kappa } - \frac{1}{\kappa _{\mathrm{crit}}^R} \right) \end{aligned}$$
(54)

and \(\sigma _{\mu \nu }=(i/2)[\gamma _{\mu },\gamma _{\nu }]\). We stress that the critical value of the bare mass (or, equivalently, of the hopping parameter) which corresponds to a vanishing renormalized mass depends on the representation.

The gauge part entering the fermionic \({\mathcal {O}}(a)\) improvement is given by

$$\begin{aligned} {\hat{F}}_{\mu \nu }(x)= \frac{1}{8}\left[ {\mathcal {Q}}_{\mu \nu }(x) {-} {\mathcal {Q}}_{\nu \mu }(x) \right] \quad \text {with} \quad {\mathcal {Q}}_{\mu \nu }(x)={\mathcal {Q}}^{\dag }_{\nu \mu }(x) \end{aligned}$$
(55)

where \({\mathcal {Q}}_{\mu \nu }\) is the clover combination of plaquettes around the point x, while the improvement coefficient \(c_{\mathrm{sw}}\) can be expanded perturbatively as

$$\begin{aligned} c_{\mathrm{sw}}(g_0^2) = 1 + c^{(1);R}_{\mathrm{sw}}g_0^2 + {\mathcal {O}}(g_0^4) \end{aligned}$$
(56)
Fig. 2
figure 2

Left panel: example of Monte Carlo history of the gauge force. Right panel: example of Monte Carlo history of fermionic forces in the fundamental and two-index antisymmetric representations are displayed on the right. Note that we are employing a multi-level integration scheme with a relative factor among fermionic and gauge forces of 8

In this work, \(c_{\mathrm{sw}}(g_0^2)\) is fixed to its tree-level value.Footnote 3 Denoting the molecular-dynamics integration time by \(\tau \), the equations of motion can be written as

$$\begin{aligned}&\frac{d}{d\tau }U_{\mu }(x) = \pi (x,\mu )U_\mu (x)\, \end{aligned}$$
(57)
$$\begin{aligned}&\frac{d}{d\tau }\pi (x,\mu ) = F(x,\mu ), \end{aligned}$$
(58)

where the dynamics of the gauge link is governed by the force \(F(x,\mu )\) which reads

$$\begin{aligned} F(x,\mu ) = F_g(x,\mu ) + \sum _R^{N_{\mathrm{rep}}}F^R_f(x,\mu ). \end{aligned}$$
(59)

The force terms entering the HMC Hamilton equations are implicitly defined through

$$\begin{aligned} \delta S_g&= -(\delta \alpha ,F_g) \end{aligned}$$
(60)
$$\begin{aligned} \delta S_f&= -(\delta \alpha ,F_f). \end{aligned}$$
(61)

The variation of the gauge action (which is defined in terms of fundamental link variables) reads

$$\begin{aligned} \delta S_g = \frac{\beta }{N} \sum _x \sum _{\mu ,\nu } \delta \alpha _{\mu }^a(x) {{\,\mathrm{Re}\,}}{{\,\mathrm{Tr}\,}}\left[ iT_F^a U_{\mu }(x)V^{\dagger }_{\mu }(x) \right] \end{aligned}$$
(62)

with \(V_{\mu }(x)\) the sum of the forward and backward staples around the link \(U_{\mu }(x)\). The fermionic force is more intricate to derive. Dropping the R superscript and the site index to avoid cumbersome notation, the fermionic action variation is

$$\begin{aligned} \delta S_f = \kappa \sum _x \phi ^{\dag } (D^{\dag } D)^{-1} \delta (D^{\dag } D) (D^{\dag } D) \phi ; \end{aligned}$$
(63)

defining the modified pseudofermion fields

$$\begin{aligned} X&=(D^{\dag }D)^{-1} \phi \end{aligned}$$
(64)
$$\begin{aligned} Y&=DX, \end{aligned}$$
(65)

Eq. (63) simplifies to

$$\begin{aligned} \delta S_f = \kappa \sum _x \left( X^{\dag }(\delta D)Y + Y^{\dag }(\delta D)X \right) . \end{aligned}$$
(66)

In the case of the Wilson action (i.e. \(D=D_{\mathrm{Wilson}}\)), from Eq. (66) we have

$$\begin{aligned} \delta S_f^{\mathrm{Wilson}}&= i \kappa \sum _x {{\,\mathrm{Tr}\,}}\sum _{\mu } \delta \alpha ^a_{\mu }(x)T_R^a \left[ U_{\mu }(x)Y(x+{\hat{\mu }})\right. \nonumber \\&\quad \quad X^{\dag }(x)(\mathbb {1}{+} \gamma _{\mu })-Y(x)X^{\dag }(x{+}{\hat{\mu }})U^{\dag }_{\mu }(x)(\mathbb {1}{-} \gamma _{\mu }) \nonumber \\&\quad - X(x)Y^{\dag }(x+{\hat{\mu }})U^{\dag }_{\mu }(x)(\mathbb {1}+ \gamma _{\mu })\nonumber \\&\quad \left. +U_{\mu }(x)X(x+{\hat{\mu }})Y^{\dag }(x)(\mathbb {1}- \gamma _{\mu }) \right] \nonumber \\&= i \kappa \sum _x \sum _{\mu } {{\,\mathrm{Tr}\,}}_{\mathrm{color}} \left\{ \delta \alpha ^a_{\mu }(x)T_R^a \left[ U_{\mu }(x)W_{\mu }(x) {-} h.c. \right] \right\} \end{aligned}$$
(67)

with

$$\begin{aligned} W_{\mu }(x)= & {} {{\,\mathrm{Tr}\,}}_{\mathrm{spin}} \left[ Y(x+{\hat{\mu }})X^{\dag }(x) (\mathbb {1}+ \gamma _{\mu }) \right. \nonumber \\&\left. + X(x+{\hat{\mu }})Y^{\dag }(x)(\mathbb {1}- \gamma _{\mu }) \right] . \end{aligned}$$
(68)

On the other hand, the variation of the clover term defined in Eq. (53) reads

$$\begin{aligned} \delta S_f^{\mathrm{clov}}= & {} - \frac{i}{16} c_{\mathrm{sw}}(g_0^2)\kappa \sum _x \sum _{\mu ,\nu }{{\,\mathrm{Tr}\,}}_{\mathrm{color}} \left[ i \delta \alpha ^a_{\mu }(x) T_R^a U_{\mu }(x)\right. \nonumber \\&\left. C_{\mu }(x) + i C^{\dag }_{\mu }(x)U^{\dag }_{\mu }(x)\delta \alpha ^a_{\mu }(x)T_R^a \right] . \end{aligned}$$
(69)

The derivation of the above formula is reported in Appendix C. All equations above hold for a generic representation R; the dependence on the representation only enters \(\delta \alpha _{\mu }(x)\) (the forces for a specific set of bare parameters are displayed in Fig. 2 and the CG solver iterations in Fig. 3). In this way the MD equations can be easily generalized to arbitrary matter content, including for fields in multiple representations.

Fig. 3
figure 3

Comparison between conjugate-gradient-solver iterations for the fundamental and the two-index antisymmetric representation at degenerate bare fermion masses \(am_4=am_6=-\,0.55\). As expected, the Dirac operator in the fundamental representation at a fixed value of the bare mass has smaller eigenvalues than the one in the sextet representation. The Dirac operator for the fundamental representation is then more ill-conditioned than its sextet counterpart and requires more solver iteration to reach the same residual

5 Observables

Having discussed our results for elementary algorithmic quantities that can be monitored in the lattice simulation (such as plaquette expectation values, Monte Carlo histories of forces involved in the HMC algorithm, etc.), in this section we present our results from the computation of Dirac spectra, as discussed in Sect. 3, and of basic phenomenological observables which can be extracted from two-point correlation functions of “meson-like” and “baryon-like” states. With this terminology inspired by hadron physics, we respectively indicate hypercolor-singlet states built from a fermion and an anti-fermion, and from fermions only. We discuss in detail how these correlation functions can be built on the lattice and we provide numerical results for the “meson-like” states, while we restrict ourselves to a theoretical treatment of the “baryon-like” correlators, whose study is postponed to future works.

In particular, we focus on quantities providing information on the critical line of the theory. For this purpose, the best-suited quantities are the fermion masses defined in terms of the partially conserved axial current (PCAC), the masses of the “pion-like” states, that are interpreted as the pseudo-Nambu–Goldstone bosons associated with the breaking of chiral symmetry, and the distribution of the smallest eigenvalue of the Dirac operator, which is expected to get smaller when one approaches the critical line. Monitoring these quantities allows one to map out the phase structure of this lattice theory with clover-Wilson fermions in different representations, which is a necessary step before embarking in exhaustive investigation of its phenomenology – a task that we leave for future work.

Detailed results of the present study are shown in the figures and in the tables included here.

Fig. 4
figure 4

Unfolded density of eigenvalue spacings obtained from spectra of the Wilson Dirac operator with a clover improvement term, for the fundamental (left-hand-side panel) and sextet (right-hand-side panel) representations of the \(\mathrm {SU}(4)\) gauge group. The results were obtained from quenched simulations on a lattice of size \((L/a)^4=4\) for \(\beta =10.0\) and \(am_4=am_6=-\,0.2\). The numerical results are consistent with the Wigner surmise according to the symmetries of the operator, i.e. the chUE curve for the fundamental representation, and the chSE for the two-index antisymmetric representation. For completeness, the plots also show the curves corresponding to the chOE, and the Poissonian distribution that would be expected, if, instead of the eigenvalues of an operator, one were considering the spacings between uniformly distributed random real numbers

Fig. 5
figure 5

Same as in Fig. 4, but for the unfolded density of eigenvalue spacings obtained from spectra of the staggered Dirac operator in the fundamental (left-hand-side panel) and sextet (right-hand-side panel) representations of the \(\mathrm {SU}(4)\) gauge group. The results were obtained from quenched simulations on a lattice of hypervolume \((L/a)^4=4\). Also in this case, the numerical data follow the analytical curves predicted by the Wigner surmise, according to the anti-unitary symmetries of the operator: the results for the fundamental representation are in very good agreement with the prediction for the chUE, and those for the sextet representation match the prediction for the chOE

5.1 Unfolded distributions of Dirac-spectrum spacings

The analytical motivation for the study of unfolded distributions of the spacings between subsequent eigenvalues of the Dirac operator is discussed in detail in Sect. 3.

In our computation we define the unfolded density of eigenvalue spacings as follows. First, we compute the spectrum of \(\gamma _5 D\) on a set of \(n_{{\tiny {\mathrm{conf}}}}\) configurations, then we sort all non-degenerate eigenvalues in increasing order, labeling each of them by a positive integer that represents the eigenvalue position in the list. Then, the spacing s between subsequent eigenvalues in each configuration c is defined to be proportional to the difference of their positions in the list:

$$\begin{aligned} s = \frac{n_{i+1}^{(c)}-n_{i}^{(c)}}{{\mathcal {N}}}, \end{aligned}$$
(70)

where the normalization factor \(1/{\mathcal {N}}\) is fixed by requiring the average value of s to be equal to one, and the unfolded density of spacings, also normalized to one, is obtained by dividing the real non-negative half-axis into intervals of width \(\delta s\), and counting how many values of s are found in a generic interval \([k\delta s, (k+1)\delta s]\), with \(k \in {\mathbb {N}}\).

Figure 4 shows our results for the unfolded density of eigenvalue spacings that we extracted from an ensemble of spectra of the Wilson Dirac operator with clover improvement term, that we use in this work, which shares the same global anti-unitary symmetries as the continuum Dirac operator. The results were obtained from quenched simulations on a lattice of hypervolume \(L^4=(4a)^4\) at \(\beta =10.0\), with the choice \(am_4=am_6=-\,0.2\).

The results for the fundamental representation (in the left panel of the figure) and for the two-index antisymmetric representation (right panel) are in complete agreement with the predictions from the Wigner surmise, Eq. (37), for the expected Dyson indices (\(\beta =2\) for the fundamental representation and \(\beta =4\) for the sextet representation). For completeness, we also show the analytical predictions for the chOE, as well as the exponential distribution that would correspond to the unfolded spacing obtained from uniformly distributed random real numbers.

Similarly, Fig. 5 shows the results that we obtained from the same type of analysis, but using the staggered Dirac operator. As discussed in Sect. 3, the global anti-unitary symmetries of this operator for fermions in the two-index antisymmetric representation are different from those of the continuum Dirac operator, and this is confirmed by our numerical results shown in the right-hand side of this plot, which follow the chOE.

5.2 Meson-like observables

For a generic \(\mathrm {SU}(N)\) gauge group and Lorentz structure \(\Gamma _A\), the fermion bilinear with flavor indices \(f_1\),\(f_2\) in has the form

$$\begin{aligned}&O_A(x)={\overline{\psi }}_{f_1}(x)\Gamma _A \psi _{f_2}(x)\, \quad \text {with} \quad \nonumber \\&\quad \Gamma _A \in \left\{ \mathbb {1},\gamma _5,\gamma _{\mu },\gamma _{\mu }\gamma _5,\sigma _{\mu \nu }\right\} , \quad \text {and} \quad f_1\ne f_2. \end{aligned}$$
(71)

where the fermion field \(\psi \) can be in any representation of the gauge group. The two-point function can be written as

$$\begin{aligned} \langle O_A(x){\overline{O}}_B(y)\rangle =\langle {\overline{\psi }}_{f_2}(x)\Gamma _A\psi _{f_1}(x){\overline{\psi }}_{f_1}(y)\Gamma _B\psi _{f_2}(y)\rangle .\nonumber \\ \end{aligned}$$
(72)

Using Wick’s contractions, the above equation can be rewritten as

$$\begin{aligned} \langle O_A(x){\overline{O}}_B(y)\rangle&= -\Gamma _A^{ij}\Gamma _B^{kl} \langle \psi _{f_1}^l(y){\overline{\psi }}_{f_1}^i(x) \rangle \langle \psi _{f_2}^j(x){\bar{\psi }}_{f_2}^k(y)\rangle \nonumber \\&= - {{\,\mathrm{Tr}\,}}[\Gamma _A S(x,y)_{f_2} \Gamma _B S(y,x)_{f_1}], \end{aligned}$$
(73)

where S denotes the fermion propagator in coordinate space. Its \(\gamma _5\)-Hermiticity \(S^{\dagger }(y,x)=\gamma _5 S(x,y) \gamma _5\) implies

$$\begin{aligned} \langle O_A(x){\overline{O}}_B(y)\rangle = - {{\,\mathrm{Tr}\,}}[\gamma _5\Gamma _A S(x,y)_{f_2} \Gamma _B \gamma _5 S^{\dagger }(x,y)_{f_1}].\nonumber \\ \end{aligned}$$
(74)

This structure holds for fermions in any representation. In fact for a generic representation R we have \(R\otimes {\overline{R}}=\mathbb {1}\oplus \dots \), i.e. it is always possible to identify a hypercolor-singlet made of a fermion–antifermion pair.

5.3 Baryon-like observables

Let us refer to fermionic fields in the fundamental representation as \(q^a_i(x)\), where \(a=1,\ldots , N\) is a hypercolor index while i is a Dirac index, and to fields in the two-index antisymmetric representation as \({\mathcal {Q}}^{a b}_j(x)\) with spin j and \(a,b=1,\ldots ,N\). In order to avoid cumbersome notation we map the two-index into a single one \((a,b) \rightarrow \alpha =1,\ldots ,N(N-1)/2\) as discussed in Sect. 3, i.e. by sorting the two-index pairs as (1, 2), (1, 3), (2, 3), (1, 4), (2, 4), (3, 4), ..., \((N-1,N)\).

It is a trivial consequence of group-representation theory that the minimum number of fermions in the fundamental representation of the \(\mathrm {SU}(N)\) gauge group to construct a hypercolor-singlet state is N. In the current context, this corresponds to “baryon-like” states formed by four (fundamental) fermions, with a qqqq structure.Footnote 4 Similarly, hypercolor-singlet states can also be built from three fermions in the two-index antisymmetric representation fermions \(\mathcal {QQQ}\). A further, “hybrid” type of color-singlet states can be built by combining fermions in both representations, as in \(qq{\mathcal {Q}}\). In the present work we restrict ourselves to the study of this three-fermion baryon, which, playing the rôle of the top-quark partner in the model under investigation, is particularly interesting. Such a state is often referred to as a “chimera baryon”. The simplest interpolating operator for this stateFootnote 5 can be written as

$$\begin{aligned} O_{N_{\pm }}(x) = \epsilon _{abcd} P_{\pm } \Gamma _A q_a(x)( q_b^T(x) \Gamma _B {\mathcal {Q}}_{c d}(x)) \end{aligned}$$
(75)
$$\begin{aligned} {\overline{O}}_{N_{\pm }}(x) = \epsilon _{abcd} ({\overline{q}}_a(x) \Gamma _B \overline{{\mathcal {Q}}}_{bc}^T(x)) {\overline{q}}_d(x) \Gamma _A P_{\pm } \end{aligned}$$
(76)

where \(P_{\pm }=(\mathbb {1}\pm \gamma _0)/2\) projects onto the desired isospin channel, and \((\Gamma _A,\Gamma _B)\) define the spin content of the baryon. For the channel with angular momentum and parity quantum numbers \(J^P=1/2^+\), common choices are \((\Gamma _A,\Gamma _B) \in \left\{ (\mathbb {1},{\mathcal {C}}\gamma _5),(\gamma _5,{\mathcal {C}}),(\mathbb {1},i\gamma _0{\mathcal {C}}\gamma _5) \right\} \), where \({\mathcal {C}}=\gamma _0\gamma _2\) denotes the charge-conjugation matrix. The two-point contraction for these three-fermion objects can be written as

$$\begin{aligned}&\langle O_{N_{\pm }}(x) {\overline{O}}_{N_{\pm }}(y)\rangle \nonumber \\&\quad = -\epsilon _{abcd} \, \epsilon _{a'b'c'd'} \, \bigg [{\overline{q}}_a'^i(y) \Gamma _B^{ij}\overline{{\mathcal {Q}}}_{b'c'}^j(y)\bigg ]{\overline{q}}_d'^k(y)\Gamma _A^{kl} \, P_{\pm }^{lm}\nonumber \\&\qquad \Gamma _A^{mn}q_{a}^n(x)(q_{b}^o(y)\Gamma _B^{op}{\mathcal {Q}}_{cd}^p(x)) \nonumber \\&\quad = \epsilon _{abcd} \, \epsilon _{a'b'c'd'} \, \Gamma _B^{ij}(\Gamma _A P_{\pm }\Gamma _A)^{kn}\Gamma _B^{op} {\mathcal {K}}_{bcb'c'}^{pj}\nonumber \\&\qquad [S_{ad'}^{nk}S_{ba'}^{oi} - S_{aa'}^{ni}S_{bd'}^{ok}] \nonumber \\&\quad = \epsilon _{abcd} \, \epsilon _{a'b'c'd'} \, \{ {{\,\mathrm{Tr}\,}}[(\Gamma _A P_{\pm } \Gamma _A) S_{ad'}] \, \nonumber \\&\qquad {{\,\mathrm{Tr}\,}}[ S_{ba'} (\Gamma _B {\mathcal {K}}_{bcb'c'} \Gamma _B^T)^T] + \nonumber \\&\qquad -{{\,\mathrm{Tr}\,}}[(\Gamma _AP_{\pm }\Gamma _A)S_{aa'}(\Gamma _B{\mathcal {K}}_{bcb'c'}\Gamma _B^T)^TS_{bd'}] \}, \end{aligned}$$
(77)

where \(S^{ab}_{ij}\) is the fermionic propagator in the fundamental representation and \({\mathcal {K}}^{abcd}_{ij}\) is the one in the two-index antisymmetric representation, for the hypercolor indices (ab) and (cd).

By exchanging color indices, Eq. (77) can be recast into the form

$$\begin{aligned}&\langle O_{N_{\pm }}(x) {\overline{O}}_{N_{\pm }}(y)\rangle \nonumber \\&\quad =\epsilon _{abcd} \, \epsilon _{a'b'c'd'} \, \left\{ {{\,\mathrm{Tr}\,}}[(\Gamma _A P_{\pm } \Gamma _A) S_{aa'}] \, \right. \nonumber \\&\qquad {{\,\mathrm{Tr}\,}}[S_{bb'}(\Gamma _B {\mathcal {K}}_{cdc'd'} \Gamma _B^T)^T ] \nonumber \\&\left. \qquad + {{\,\mathrm{Tr}\,}}[(\Gamma _AP_{\pm }\Gamma _A)S_{aa'}(\Gamma _B{\mathcal {K}}_{bcb'c'}\Gamma _B^T)^TS_{dd'}] \right\} . \end{aligned}$$
(78)

Equation (78) is formally identical to the one relevant for the nucleon in quantum chromodynamics, where all quark fields are in the fundamental representation of the \(\mathrm {SU}(3)\) gauge group. It is well known that two-point functions interpolating baryonic states are typically very noisy, compared to the ones for mesons: this is mostly due to the presence of an additional propagator with respect to the mesonic case. To extract a clear signal from these correlation functions, several techniques have been developed (see Ref. [32] and references therein). In the theory investigated in this work, the problem is expected to be even more severe, due to the presence of the propagators in the sextet representation, hence we postpone a systematic study of baryon spectroscopy to a future publication.

Fig. 6
figure 6

On the top panels we show an example of the pseudoscalar correlators obtained with degenerate bare fermion masses \(am_4=am_6=-\,0.55\), on the bottom left panels the effective masses of “pion-like” states, for both representations for the pseudoscalar correlator in unit of lattice spacing while on the bottom right we display the PCAC fermion mass

5.4 Extraction of effective masses

Once the correlators are computed we project to zero-momentum by summing on the space directions \({\mathbf {x}}\) as

$$\begin{aligned} C(t)=\sum _{{\mathbf {x}}} \langle {\overline{O}}({\mathbf {x}},t)O({\mathbf {x}},0)\rangle . \end{aligned}$$
(79)

The masses of pseudoscalar (“pion-like”) and vector (“\(\rho \)-like”) states are respectively extracted from the asymptotic behavior of the \(C_{PP}(t)\) and \(C_{V_{i}V_{i}}(t)\) correlators. For large Euclidean-time separation, the former behave like

$$\begin{aligned} C_{PP}(t) \propto \exp \left\{ {-}M_{PP}t \right\} {+} \text {contribution from excited states}.\nonumber \\ \end{aligned}$$
(80)

In addition, in a system of finite Euclidean-time extent \(L_t\), where fermionic fields obey anti-periodic boundary conditions in the Euclidean-time direction, the correlator above also receives contributions from the periodic copies of the operators, resulting in additional terms like \(\exp \left\{ -M_{PP} (L_t -t) \right\} \), etc. on the right-hand side of Eq. (80).

The mass of the “meson-like” states is thus obtained by fitting the decay of the correlators at sufficiently large t, including the effect of the first periodic copy of the operators. That is, we define

$$\begin{aligned} aM_{\mathrm{eff}}= {{\,\mathrm{arcCosh}\,}}\left[ \frac{C_{PP}(t+a) + C_{PP}(t-a)}{2C_{PP}(t)} \right] . \end{aligned}$$
(81)

The same analysis is applied to the correlator involving the i-th component of the vector current \(C_{V_iV_i}(t)\). In order to study the distance from the critical line of the theory we also consider the PCAC fermion mass defined through the non-anomalous axial Ward identity

Fig. 7
figure 7

Distribution of the minimum eigenvalue of \(\gamma _5 D\) in the fundamental (blue) and two-index antisymmetric (red). The critical line is approached from the top to the bottom corresponding to bare fermion masses \(am_4=am_6=[-0.50,-0.55,-0.58]\). We observe that for degenerate bare fermion masses the fundamental representation is lighter than the two-index antisymmetric one, i.e. \(\langle |\lambda _{\mathrm{min}}^{\mathrm{Fund}}|\rangle < \langle |\lambda _{\mathrm{min}}^{\mathrm{2AS}}|\rangle \). This result is consistent with the hierarchies found in PCAC quark masses and in “pion-like” effective masses displayed in Fig. 6. This observation is also compatible with the results reported in Ref. [18]

$$\begin{aligned} am_{PCAC}= \frac{\tilde{\partial }_0 C_{AP(t)}}{2C_{PP}(t)} \end{aligned}$$
(82)

with \(\tilde{\partial }_0=(\partial _0 + \partial ^*_0)/2\) the symmetric derivative in the time-direction. Note that the PCAC fermion mass approaches to the continuum limit linearly in the lattice spacing. \({\mathcal {O}}(a)\) effect would be removed by considering the improved axial correlator \(C^{\mathrm{I}}_{AP}(t) = C_{AP}(t) + c_A(g_0) \, \tilde{\partial }_0 C_{PP}(t)\), with the (currently) unknown coefficient \(c_A(g_0)\) which depend on both number of colors and the representation of fermions. The top panels of Fig. 6 illustrates the typical hyperbolic-cosine shape of the pseudoscalar correlator in both representations, while bottom panels in Fig. 6 show fits to plateau region for the extraction of the two correspondent effective masses. Similar plot are provided for the PCAC fermion mass on the bottom right of Fig. 6.

5.5 Spectral observables and scale setting

As discussed in Sect. 3, a very interesting observable to probe the chiral regime of the theory is provided by the study of the low lying spectra of the Dirac operator in both representations under investigation. In this section rather than the Dirac operator itself, we prefer to consider the hermitian operator \(\gamma _5 D\), since the latter is Hermitian and hence has a real spectrum. On finite lattice the smallest eigenvalues of the Dirac operator defines a spectral gap

$$\begin{aligned} |\lambda _{\text {min}}| = \text {min}\{ |\lambda | : \lambda \text { is an eigenvalue of } \gamma _5D \}. \end{aligned}$$
(83)

As a further control on the critical line of the theory we observe the scaling of \(|\lambda _{\text {min}}|\) with the bare mass. An example of showing the drift of the smallest eigenvalues is showed in Fig. 7. Our results for the plaquette and the smallest eigenvalues in two different volumes and for several bare parameters are summarized in Tables 3 and 4. We note here that at degenerate bare fermion masses the spectral gap is much larger in the two-index representation respect to the fundamental one. This picture is consistent with both the PCAC fermion masses and pion masses. The scale is set using the Wilson flow introduced in Ref. [33]. The reference scale \(t_0\) is implicitly defined via the relation (generalized to \(\mathrm {SU}(N)\) as in Ref. [34])

$$\begin{aligned} \left. t^2 \langle E(t) \rangle \right| _{t_0} = 0.1125\frac{(N^2-1)}{N}, \end{aligned}$$
(84)

where the action density \(E(t)=\frac{1}{4}G^a_{\mu \nu }(t)G^a_{\mu \nu }(t)\) is constructed from the plaquette, formed by gauge links at flow time t. The r.h.s. of Eq. 84 is chosen to be a dimensionless number according to perturbative expansion at small t, reducing to 0.3 for \(N=3\).Footnote 6 Note that the (Gaußian) smearing radius of the Wilson flow scales with the flow time as \(\sqrt{8t}\). Hence, in order to avoid over-smearing we imposed \(t\le L^2/32\), with L as the shortest direction in our lattice. An example of fit used to extract the value of \(t_0/a^2\) is displayed in Fig. 8. We observe that for values of \(\beta <10.0\), where we expect a bulk phase transition the scale cannot be set since the reference scale is reached too fast and within the initial transient regime. This is a further confirmation pointing to an unphysical phase fully dominated by cutoff effects. However assessing the nature of such a phase would requires further investigations on larger volumes and more values of the bare gauge coupling.

Table 3 Table run, volume \((8^3\times 16)a^4\), plaquette gauge action and fermionic Wilson-clover \(N_f=2+2\) action. Runs \(A17-A26\) use the same bare parameters as in Ref. [35], however a direct comparison cannot be done, since in this work we use a different gauge action with respect to Ref. [35]. Nevertheless, the tension between our results and the ones in Ref.  [35] seems to indicate a surprisingly relevant shift of the line of constant physics due to the smearing procedure
Table 4 Table run, volume \((16^3\times 32)a^4\), plaquette gauge action and fermionic Wilson-clover \(N_f=2+2\) action
Fig. 8
figure 8

The purple band shows \(t^2\langle E(t)\rangle \) with its uncertainty, while the vertical green band is the value of \(t_0\) implicitly defined by requiring \(t_0^2\langle E(t_0)\rangle =0.421875\). As expected after a first transient induced by the plaquette discretization of the energy, the plotted quantity enters a polynomial regime

5.6 Discussion

The results presented here deserve some comments.

First of all, our data confirm that the simulation code that we used, featuring a Wilson Dirac operator with a clover improvement term, is a robust tool to explore the phase structure of this theory. Beside reproducing well-known results in the quenched limit, it also passes all other required algorithmic and physics consistency checks, and turns out to be efficient and easy to generalize to arbitrary matter content.

Our investigation of the spectrum of the Dirac operator in the \(\mathrm {SU}(4)\) theory with matter in \(2+2\) different representations confirms the non-trivial implications of the global anti-unitary symmetries of sextet fermions, and proves that the spectral properties of the continuum operator are correctly reproduced in our lattice simulation. Moreover, the distribution of the lowest eigenvalue of the Hermitian \(\gamma _5 D\) operator, which is a useful probe to study the chiral limit of the theory, follows what is expected from general arguments (e.g. the absolute value of the lowest eigenvalue of \(\gamma _5 D\) for fundamental fermions is always smaller than for sextet fermions, etc.).

Similarly, the investigation of “meson-like” hypercolor-singlet states that is summarized in Tables 56, and 7 provides useful information about the non-perturbative dynamics of this theory, and, again, confirms that states built from fermions in the two-index antisymmetric representation are generally heavier than those from fundamental valence fermions. Also, the mass hierarchies between pseudoscalar and vector states follow a pattern similar to the one familiar from quantum chromodynamics, and are consistent with how our results for PCAC masses for the fermions scale.

The plaquette expectation values reported in Tables 3 and 4 appear to reveal the presence of a rather large strong-coupling phase, likely dominated by quite severe, unphysical discretization effects: an important piece of information for future studies of this model with this lattice discretization. We also note a significant shift of the lines (or “surfaces”) of constant physics with respect to the analysis reported in Ref. [35] and in subsequent works by that group; however, it should be emphasized that any possible discrepancy between the parameters in our work and theirs does not imply that these studies are inconsistent with each other, simply because they are based on different lattice discretizations, and, by virtue of universality, only continuum-extrapolated physical results should agree. For our scale setting in terms of the \(t_0\) parameter, see also Table 7.

6 Generalization to other partial-compositeness models

While the numerical study reported in this work is restricted to (a slightly simplified version of) the theory proposed in Ref. [16], it should be remarked that this is only one in a broad class of partial-compositeness models potentially relevant to describe the electroweak-symmetry breaking mechanism and physics at the TeV scale. Hence, it would be interesting to study also other strong-dynamics models, with low-energy symmetries compatible with those of the Standard Model, but based on other gauge groups and/or with a different matter content.

In fact, the simulation code that we used in this work is very versatile and the exploration of the phase structure and physical observables that was carried out here could be easily repeated for other models.

As we mentioned, the model originally proposed in Ref. [16] features five Weyl fermions in the sextet representation, but in the present study we considered a closely related theory, which instead has two Dirac fermions in the sextet representation. Beside being simpler to simulate, the motivation underlying this choice is that the model with two sextet Dirac fermions (and two fundamental ones), which is an excellent proxy for the original model, has also been studied in other recent works [18,19,20,21], and, as usual, testing the universality of physical results obtained with a different lattice regularization is an important requirement of a lattice calculation. However, as our code includes numerical rational hybrid Monte Carlo routines, it can be used to repeat the calculation for any number of fermion flavors, in arbitrary combinations of multiple representations. The generalization to larger values of the number of hypercolor charges, too, is already implemented in our code, and the computational-cost scaling with this parameter does not involve particular subtleties (see, e.g., Ref. [36, section 3]).

Furthermore, our code can be readily adapted to different types of gauge groups. In this respect, a novel and interesting proposal for a different strongly coupled New Physics model has been recently put forward in Ref. [37]. Like in the model that we considered here [16], the idea underlying the construction of this model is that the contributions to the Higgs boson mass from its Yukawa coupling to the top quark can be partially compensated for by the presence of sufficiently light top partners. However, in contrast to the proposal of Ref.  [16], the model discussed in Ref. [37] is characterized by local invariance under a symplectic, rather than a special unitary, group.

More specifically, the model described in Ref. [37] is based on the \(\mathrm {SU}(4)/\mathrm {Sp}(4)\) symmetry-breaking scheme [10, 38] and its ultraviolet completion is a vector gauge theory with local internal invariance under the \(\mathrm {Sp}(6)\) group. In addition to the gauge bosons, the field content of the theory includes ten fermions in the fundamental representation, and one in the adjoint representation of the gauge group. The choice of this internal symmetry and matter fields comes from the requirements of a global symmetry sufficiently large to include the gauge group of the Standard Model, the existence of a non-linearly realized symmetry that could protect the mass of the Higgs boson from arbitrarily large quantum corrections, and the presence of massless fermions compatible with the ’t Hooft anomaly-matching conditions. As discussed in Ref. [37], this model is expected to present a rich low-energy phenomenology, which could include top-quark partners, scalar particles, and color-charged fermions. These features make it an interesting target for non-perturbative lattice calculations – a research program that could be a natural generalization of the present work.

Table 5 This table reports the value of the masses for the pseudoscalar (\(M_{P}\)) and the vector (\(M_{V}\)) states, together with the PCAC fermion mass. We note, as a consistency check, that the vector particle masses evaluated from correlators of their three different components are compatible with each other
Table 6 Table 5, continued
Table 7 Same as in Table 5, but including also meson-like states constructed from fermions in the sextet representation. We observe that at fix bare fermion mass these states are heavier than the ones built from fundamental fermions. This observation is consistently supported by the value of the pseudoscalar-state masses, PCAC fermion masses, as well as the average smallest eigenvalue of the Dirac-Wilson operator. In the last column, we also report the value of the scale-setting parameter \(t_0/a^2\)

It is worth remarking that the lattice investigation of \(\mathrm {Sp}(2N)\) gauge theories with dynamical fermions has already begun [39], and extending this type of calculations to the model described in Ref. [37] should be feasible with a minor effort with the technology already developed for the current project.

7 Concluding remarks and future perspectives

In the present article, we reported our results of a non-perturbative lattice investigation of a non-Abelian \(\mathrm {SU}(4)\) gauge theory with two dynamical flavors of fundamental Dirac fermions, and two dynamical flavors of Dirac fermions in the two-index antisymmetric representation. As discussed in the introduction, the main motivation to study this model is its close proximity to the simplest UV-complete partial-compositeness model, that was introduced in Ref. [16], and that may provide a solution to some of the tantalizing conundrums of the present state of affairs in theoretical elementary particle physics: in particular, it features a composite Higgs boson and a partially composite top quark. While the model studied in the present work has slightly different matter content with respect to the one advocated in Ref. [16], it is expected to capture its main features at least semi-quantitatively, and to provide useful guidance for future studies.

We carried out our Monte Carlo calculations by adapting existing code to a setup with fermionic matter in multiple, arbitrary representations; moreover, this code already supports rational hybrid Monte Carlo routines, so that an extension to an arbitrary number of fermion flavors would be straightforward. At the technical level, our lattice discretization of the continuum theory is based on a Wilson Dirac operator with clover improvement term. Our setup is, thus, slightly differentFootnote 7 with respect to the one used in Ref. [35] and in later works by that group [18,19,20, 40].

As discussed in detail in Sect. 5.6, the results that we presented here provide a clear picture of the phase structure of this lattice theory, and confirm important properties related to its global symmetries, as well as its non-perturbative dynamics. While this could already provide a useful roadmap for further lattice investigation of this model, it should be pointed out that the results of the very recent paper [40] appear to rule out the viability of this model for a partial-compositeness scenario: they indicate that the renormalized overlap factors relevant for the mixing of “chimera” states with the top quark are too small, and disfavor its phenomenological relevance for New Physics. The possibility that this problem could be evaded through a four-fermion coupling enhanced at low energies by a large, negative anomalous dimension was also ruled out, in particular in view of the QCD-like, rather than conformal, behavior of the spectroscopy of this theory, which our present results also confirm.

As we pointed out in Sect. 6, however, an interesting alternative partial-compositeness model has been recently proposed in Ref. [37], and the simulation algorithm that we used in the present study is sufficiently versatile to use it for the study of this model, too. The lattice investigation of strongly coupled models for New Physics, (see Ref. [41] for a very recent review) remains an active research field.