Spectroscopy of SU(4) composite Higgs theory with two distinct fermion representations

We have simulated the SU(4) lattice gauge theory coupled to dynamical fermions in the fundamental and two-index antisymmetric (sextet) representations simultaneously. Such theories arise naturally in the context of composite Higgs models that include a partially composite top quark. We describe the low-lying meson spectrum of the theory and fit the pseudoscalar masses and decay constants to chiral perturbation theory. We infer as well the mass and decay constant of the Goldstone boson corresponding to the non-anomalous U(1) symmetry of the model. Our results are broadly consistent with large-Nc scaling and vector-meson dominance.


I. INTRODUCTION
Gauge theories coupled simultaneously to more than one fermion representation-"multirep" theories-open a new dimension in the study of gauge dynamics. Apart from the influence of each fermion species on the gauge field and vice versa, phase transitions and symmetry breaking in each species can affect the others dramatically. Of course, QCD already contains light quarks, strange quarks, and heavy quarks, and the influence of each species on the others is an old and continuing object of QCD calculations. The difference is that QCD's quarks are all equivalent, in that a tuning of the masses can change one into another. Fermions in inequivalent representations, on the other hand, enter the dynamics with different strengths irrespective of their masses.
As usual, symmetries offer the clearest perspective on the physics of inequivalent fermions. Each species has its maximal flavor symmetry, while no symmetries mix the different species. If all the fermions are made massless, the chiral symmetries of the species remain distinct. One symmetry could break spontaneously while others do not. This is a generalization of the old issue of scale separation, which was originally seen as a possible separation of a chiral scale from the confinement scale of the gauge theory [1][2][3]. It is possible that inequivalent representations, simultaneously coupled to the gauge field, define independent chiral scales. This might find expression in the finite-temperature physics of the theory, in the form of distinct phase transitions for each fermion species as well as for the confinement physics of the gauge field. Alternatively, one phase transition, possibly governed by the largest quadratic Casimir of the fermion representations, might trigger all the others to occur at the same scale.
We present here the first results of our work on the SU(4) gauge theory with N f = 2 Dirac fermions in each of two distinct representations, the fundamental 4 and two-index antisymmetric 6 (a real representation). We have chosen this model because it is close to a model proposed by Ferretti for a hypercolor theory that yields a composite Higgs boson [4,5] and a partially composite top quark [6]. Ferretti's model [7] contains 5 Majorana fermions in the sextet representation and 3 Diracs fermions in the fundamental; simulating this fermion content requires the costly rational hybrid Monte Carlo (RHMC) algorithm, and so, instead, we study the theory with 4 Majoranas (equivalent to 2 Dirac fermions) in the sextet and 2 Diracs in the fundamental. In Ferretti's model, the massless sextet Majorana fermions Ψ condense to break their chiral symmetry according to SU(5) → SO(5), whereupon the Standard Model's Higgs multiplet appears as Nambu-Goldstone (NG) bosons; our symmetry breaking scheme is 1 SU(4) → SO (4). The fundamental fermions ψ (4) in Ferretti's model are brought in so that the theory will possess fermionic baryons constructed as ψ (4) ψ (4) Ψ "chimera" bound states, to be used as top partners; they condense (again, in the chiral limit) according to SU(3) L ×SU(3) R → SU(3) V . In our model, the corresponding symmetrybreaking scheme is SU(2) L × SU(2) R → SU(2) V . We believe that our model contains all the qualitative physics of Ferretti's model while offering a laboratory for developing quantitative techniques.
Multirep theories of physical significance are not easy to come by. Apart from the phenomenological requirements, Ferretti's choice of model is constrained [8] by the simple fact that higher-representation matter fields push gauge theories into the conformal window unless the number of fermions is quite small. It is essential that the gauge theory of hypercolor exhibit confinement and the concomitant breaking of global symmetries.
In this work we present results from the mesonic sector of the theory, leaving baryonic observables for another paper. We have already explored the mesonic and baryonic spectrum of the SU(4) gauge theories with only fundamental [9] or only sextet fermions [10]. 2 Those results fit nicely into the body of work on QCD and its generalizations to larger values of N c . The analysis there, similar to QCD studies, related the gauge coupling β and hopping parameter κ to a physical scale r 1 (the Sommer scale) and the quark mass m q , and used the latter as an abscissa for plotting particle masses and decay constants. Here, of course, the space of bare couplings consists of the gauge coupling β and two hopping parameters κ 4 and κ 6 for the two fermion species. We translate these into the scale parameter t 0 , derived from the Yang-Mills gradient flow, and the two quark masses m 4 and m 6 .
Our main tool for understanding the meson spectrum is a recent generalization of chiral perturbation theory (χPT) to the low-energy sector of a two-representation theory [12]. This form of χPT provides formulas for masses, decay constants, and chiral condensates at next-to-leading order, with m 4 and m 6 as independent variables. These formulas contain an important qualitatively new piece of physics compared to QCD-communication between the different species. They describe, for instance, the dependence of the masses of the NG bosons of all the broken chiral symmetries on both fermion masses.
Another new feature of the two-representation theory is the existence of a non-anomalous singlet axial current, and a corresponding singlet NG boson that must be included in the low-energy chiral theory. This particle is denoted ζ in Ref. [12] and is of significant phenomenological interest for composite Higgs models [8,13,14]. In this work we do not probe this singlet pseudoscalar state directly. Nevertheless we extract information about it indirectly, via its virtual contributions to the properties of the flavored NG bosons associated with chiral symmetry breaking of the individual representations. In particular, its decay constant in the chiral limit is a parameter in the chiral Lagrangian and thus appears as a fit parameter, allowing us to infer its mass using the leading-order formula.
Besides the pseudoscalar channel, we calculate masses and matrix elements of the lightest vector bosons. The vector is the lightest narrow resonance in QCD, and its properties are closely related to those of the pseudoscalars within the framework of vector meson dominance (VMD). We explore the evidence for VMD in our theory and its consequences for the decay width of the vector. This is of particular phenomenological interest, since in composite Higgs models, the vector resonance is often one of the first signatures expected in collider searches.
The paper is organized as follows. In Sec. II we describe the lattice theory, the observables we use, and ensembles we generated. In Sec. III we describe our application of χPT, including the discretization effects of Wilson fermions, and our scale setting method which is based on t 0 . In Sec. IV we present our results for the pseudoscalar spectrum and decay constants, including the flavor singlet ζ. We present the vector particles in Sec. V and use VMD to estimate decay widths. In Sec. VI we discuss our results from the point of view of large-N c predictions, and present our overall conclusions.
The tables containing the various measured quantities have been collected together in Appendix A. In Appendix B we explain technical aspects of our analysis of lattice data. In Appendix C we review the definition of the U(1) axial current and of the mass parameter in Wilson χPT. Finally, Appendix D contains a calculation of perturbative Z-factors for the nHYP lattice action with dislocation-suppression.

A. Symmetries
The chiral symmetry of the fundamental fermions and its expected breaking are the same as in two-flavor QCD. The specifics of chiral symmetry breaking for the sextet representation are less well-known, so we will discuss them briefly; a more detailed explanation is given in [10,12].
The sextet representation of SU(4) is a real representation. Our model has two Dirac fermions charged under this representation, ψ (6) i ,ψ (6) i , i = 1, 2, which are equivalent to four Majorana fermions Ψ I , I = 1, . . . , 4. The global symmetry of the continuum theory is thus also SU(4). Using the language of Majorana fermions, the bilinear condensate Ψ I Ψ J is symmetric in its Majorana-flavor indices. Hence, after spontaneous symmetry breaking one expects the unbroken symmetry to be SO(4) [15]. One consequence of the enlarged symmetry is thatψ (6) ψ (6) mesons and ψ (6) ψ (6) diquarks (both gauge-singlet objects) are members of a degenerate multiplet of the unbroken group.
As usual, the chiral symmetries of the theory are explicitly broken by the Wilson term in the lattice action. The lattice theory thus has the same flavor symmetry as expected in the continuum theory after spontaneous symmetry breaking: SU(2) V × U(1) B for the fundamental representation and SO(4) for the sextet. Our use of Wilson fermions thus assumes that the spontaneous breaking of chiral symmetries is as would be forced by a bilinear condensate, and all measured correlation functions reflect this.
A special feature of the two-representation theory is the existence of a conserved U(1) axial current. While the individual U(1) currents J (4) 5µ and J (6) 5µ are anomalous, one can form a linear combination J 5µ of these currents that decouples from FF . Condensation of either fermion species then spontaneously breaks the non-anomalous axial symmetry, giving rise to a singlet NG boson that we denote ζ. We review the normalization of the U(1) current in Appendix C 1.

B. Lattice action and parameters
Our lattice action contains gauge-field terms and two fermion actions, one for each representation: Each fermion action is a Wilson-clover action built of gauge links constructed by nHYP smearing [16,17]. In S F the smeared links are promoted to the sextet representation [10]. There are two hopping parameters, κ 4 and κ 6 . We set both clover coefficients equal to unity, c SW = 1, a choice known to work well with nHYP smearing in QCD [18] and with fermions in higher representations [19].
The gauge-field action takes the form The first term is the usual plaquette action, while the second is an nHYP dislocationsuppression (NDS) term [20], constructed from the nHYP-smeared links. The NDS term is designed to avoid singularities in the nHYP smearing. For the present study, we hold the ratio γ/β fixed at 1/125 and use β as a free bare parameter. Concurrent with the work described here, we are also studying the finite-temperature phase structure of the theory [21,22]. Comparison of the sextet-only limit of this theory to earlier published results [10] shows that the use of the NDS action removes the previouslyobserved bulk transition from the interesting region of parameter space (see also Ref. [23]). In the multirep theory, we see no evidence for a bulk transition anywhere near the range of bare parameters at which we run, indicating that all of our ensembles correspond to the confined continuum phase with broken chiral symmetry.
We extract masses and decay constants in the usual way from two-point correlation functions. We denote pseudoscalar masses and decay constants in the representation r by M P r and F P r , respectively. The corresponding quantities in the vector channel are denoted by M V r and F V r .
We define the fermion masses m 4 and m 6 by imposing the axial Ward identity (AWI), where x = 0, and a is an isospin index. We use the local unimproved axial current A (r) µa and pseudoscalar density P (r) a in each representation r. For the determination of the AWI mass, we do not renormalize these currents because the mass itself is not a physical observables; based on our perturbative renormalization of these currents described in Appendix D (used for calculation of decay constants), the effect of including the renormalization would be small anyway, amounting to a few-percent shift of the masses. For O r we take a pseudoscalar source. When the distinction between representations is irrelevant, we will refer to the fermion mass defined by Eq. (2.3) as m AWI . Further information about our conventions and methods for spectroscopy is given in Appendix B.

C. Scale setting
We set the scale in our simulations using the flow scale, t 0 , introduced by Lüscher [24]. The flow scale is defined by the implicit equation where E(t) = 1 4 G a µν G a µν (t) is constructed from the clover form of the field strength G a µν at flow time t. Here C is a dimensionless number, conventionally [24] taken to be 0.3 in QCD. With this choice, √ t 0 corresponds to a length scale of 0.14 fm (i.e., an energy scale of 1.4 GeV) in QCD simulations [25,26].
For an arbitrary gauge theory, any value for C is a priori as good as any other. However, for comparison to existing studies with different gauge groups, it is useful to let C vary with N c . Arguments from large-N c QCD, supported by lattice data [26,27], suggest that t 0 ∼ N c at leading order. For the SU(4) simulations of this work we therefore use Lattice calculations give masses as dimensionless numbers M a and gradient-flow scales as t 0 /a 2 . Dimensionless products likeM ≡ M √ t 0 eliminate the lattice spacing a, and our tables and figures will display such quantities. To aid the intuition, one can mentally convert M √ t 0 to M/(1.4 GeV). We return to the subject of scale-setting and connect it to χPT in Section III below.

D. Ensembles
The ensembles used in this study are listed in Tables II-IV in Appendix A. They fall into three groups. The short runs with the smallest lattices, of size V = n 3 s × n t = 16 3 × 18, were used for general orientation in the three-dimensional coupling space (β, κ 4 , κ 6 ). The most important observables for this step were √ t 0 , the scale defined by gradient flow (see Tables V-VII), and the masses M P r of the pseudoscalars constructed respectively from fermions in the r = 4 and 6 representations (see Tables VIII-X).
The goal of this orientation was to find couplings that give t 0 /a 2 = O(1) along with pseudoscalar masses that are reasonably light, for subsequent comparison to χPT. It turned out that these short runs yielded results that are in themselves usable for the chiral fits to be presented below, and hence we include them in our analysis.
As can be seen in the tables, some ensembles differ in small changes to their κ r values. Our orientation runs found that t 0 /a 2 and aM P are often sensitive to these small changes.
We demanded that our ensembles satisfy the criterion M P r L > 4 for both representations, where L = n s a is the spatial size of the lattice. This is the familiar rule of thumb from QCD, based on the fact that leading-order finite-volume corrections are proportional to e −MπL ; a more detailed study of finite-volume effects in our data is given in Appendix B 4. We considered cutting data above a maximum value of t 0 /a 2 beyond which finite-volume effects severely contaminate determination of the flow scale; such a cut was found to be unnecessary following the cuts on M P L. We did eliminate ensembles with t 0 /a 2 < 0.94 because in these cases the flow did not enter a linear regime. These correspond to a large lattice spacing-in QCD language, 1/a < 1.3 GeV.
Having found interesting regions for study, we continued with high-statistics runs on lattices with V = 16 3 ×32. Finally, we have four extended runs on lattices with V = 24 3 ×48. These runs were done at large t 0 /a 2 and smallM P , so that the constraint M P L > 4 demanded an increase in L/a.
The pseudoscalar masses for all the ensembles are given in Tables VIII-X. To show our coverage of M P values, we map them in the (M P 4 , M P 6 ) plane in Figs. 1 and 2. The first shows the pseudoscalar masses obtained for 0.94 < √ t 0 /a < 1.41, which translates to a cutoff of 1.3 GeV < 1/a < 2 GeV in QCD language (most are in the neighborhood of √ t 0 /a = 1.05, or 1/a = 1.45 GeV). The second plot represents ensembles in the range 1.41 < √ t 0 /a < 1.64, or 2 GeV < 1/a < 2.3 GeV.

III. CHIRAL PERTURBATION THEORY
The standard framework for analyzing the light pseudoscalar sector is χPT. The generalization of χPT to a theory with fermions in two different representations was developed in Ref. [12], and the next-to-leading-order (NLO) results of this work provide the basis for our fits for the pseudoscalar masses and decay constants. We will also need Wilson chiral perturbation theory (WχPT), the extension of chiral perturbation theory to include the discretization errors of Wilson fermions [28][29][30][31][32][33].

A. Using a yardstick
We need a yardstick with which to measure dimensionful quantities as the fermion masses are varied. In this paper, we use √ t 0 for the characteristic length scale of every ensemble. To measure an observable in units of t 0 simply means to multiply it by the power of t 0 that renders it dimensionless. Since t 0 itself admits a chiral expansion [34], the resulting dimensionless quantity admits a chiral expansion whenever the original dimensionful observable does.
To see how this works, consider a gauge theory with mass-degenerate fermions of mass m, all in the same representation. In continuum χPT, the NLO expression for the decay constant is The remaining parameters in Eq. (3.1) are µ, the renormalization scale, and L, which is a (dimensionless) linear combination of the NLO low-energy constants (LECs), whose value depends on the choice of µ. The coefficient c of the logarithmic term is a calculable number that depends only on the fermion representation and on the number of flavors [35]. The NLO result for t 0 is where t 0,ch is the value of t 0 in the chiral limit, andk 1 is a new LEC. Notice that this expression depends analytically on the fermion mass m. As was shown in Ref. [34], logarithmic corrections to t 0 occur for the first time at the next-to-next-to-leading order (N 2 LO). Combining Eqs. (3.1) and (3.3) we obtain the NLO result for the dimensionless product Here we have chosen the renormalization scale to be µ = t −1/2 0,ch . LECs are independent of the fermion mass, and to preserve this feature we rescale them with t 0,ch , for example defining F = F √ t 0,ch . Equation (3.4) can then be written aŝ The expansion parameter is nowm, which is the fermion mass m measured in units of t 0,ch . Equation (3.5) is inconvenient becausem is not known for a given ensemble until t 0,ch is known. Finding t 0,ch (in units of t 0 of the given ensemble) requires a complicated fitting procedure that we wish to avoid. Instead, we opt for rescaling all observables of a given ensemble, including the fermion mass, with t 0 of the same ensemble. Introducingm ≡ m √ t 0 we now use Eq. (3.3) to relate the rescaled masses, which allows us rewrite Eq. (3.5) aŝ The transition fromm tom left no trace, because the difference is a higher-order correction.
More generally, at NLO the transition fromm tom can always be absorbed into a redefinition of the LECs. (A case where the redefinition is non-trivial is the NLO result for the pseudoscalar mass.) An appealing feature of Eq. (3.7) is that it looks the same as Eq. (3.1). In particular, the coefficient of the logarithmic term is unchanged. The only minor change is that the coefficient of the NLO analytic term is now L +k 1 /2 instead of L. (At N 2 LO things would become technically more complicated, because N 2 LO logarithmic corrections for t 0 would have to be incorporated as well.) It can be checked that this nice feature generalizes to an arbitrary fermion content. In the NLO fit formulae that we present below, all the logarithmic terms will thus have the same coefficients as in the usual continuum NLO expressions [12,35].

B. Wilson chiral perturbation theory
The extension of continuum chiral perturbation theory to include the discretization errors of Wilson fermions goes under the name of Wilson chiral perturbation theory, or WχPT. In the light pseudoscalar sector, WχPT allows us to extrapolate both to the chiral limit, m → 0, and to the continuum limit, a → 0. WχPT comes in two variants, depending on the choice of a power counting scheme. In this paper we follow the "generic small mass," or GSM, power counting, defined by where p is an external momentum, m is the fermion mass, and a is the lattice spacing, all measured in terms of a typical hadronic scale. The alternative power counting scheme, known as the "large cutoff effects," or LCE, power counting, is defined by The GSM scheme is appropriate when the fermion mass is not too small, and O(a 2 ) effects may be considered as subleading corrections. (We must, of course, remain within the chiral regime, meaning thatm = m √ t 0 is small.) In particular, our determination of the critical surface κ c r , where the mass of fermions in representation r vanishes, is done via extrapolation from the GSM regime. As a result, we do not probe the possible existence of an Aoki phase. For more details, see Appendices B 3 and C 2.
The fermion mass appearing in the LO lagrangian of WχPT is the so-called shifted mass, defined by m shifted = m ctm + aW 0 /B, (3.10) where m ctm is the fermion mass of continuum χPT, and W 0 is a new LEC from WχPT. The difference between the shifted and continuum masses vanishes in the continuum limit. For this lattice study, we need to know how the shifted mass m shifted compares to the fermion mass m AWI measured in our simulations via the axial Ward identity Eq. (2.3). As was shown in Ref. [36], m shifted = m AWI , up to corrections that are higher order in either of the above power counting schemes. In view of the important role that this result plays in our analysis, we briefly summarize the derivation of Ref. [36] in Appendix C 2. For our chiral fits we thus The last ingredient we need for our fits is the lattice spacing. Since we are measuring all dimensionful quantities in units of t 0 , it is natural to adopt a mass-dependent prescription, and to measure also the lattice spacing in units of t 0 . We thus introducê The Wilson discretization effects of any hatted (dimensionless) observable will be accounted for by an expansion inâ. In QCD, it is common to choose a mass-independent scale-setting prescription, whereby the lattice spacing is a function of the bare coupling β, but is independent of the bare fermion masses (see for example Refs. [37,38]). In brief, for every constant-β plane, this procedure requires finding the point where certain dimensionless quantities (such as M π /F π and M K /F K ) attain their real-world values. The value in lattice units of some dimensionful observable at the reference point is then used to determine the lattice spacing in physical units.
Here we have opted for mass-dependent scale setting because of several important differences. First, the BSM context does not provide us with any experimental results that could be used to define a reference point. This problem might be circumvented by invoking the chiral limit as a reference point on each constant-β plane. This, however, has the undesirable feature that the scale setting procedure would necessarily involve an extrapolation.
Second, in our model, as in many other models that have been studied in the BSM context, we observe a rapid change of t 0 /a 2 with the fermion mass, especially when the latter becomes light. Moreover, this phenomenon is quite general, and is seen for virtually any quantity that might be used to set the scale; its proper interpretation is thus that the lattice spacing itself is changing rapidly. The underlying reason is that, in comparison with QCD, BSM theories tend to have a large number of fermionic degrees of freedom, which have a strong screening effect on the bare coupling. When we consider the dependence of a hatted quantity, such asM P , on the hatted mass parameter,m, we expect to see some deviations from the continuum values, but such scaling violations should be small when the bare coupling is small enough. By contrast, as explained above, the lattice spacingâ itself can vary rapidly with the fermion mass(es). By using the mass-dependent scale-setting prescription of Eq. (3.12) we can incorporate this effect into our analysis. As we will see, the remaining scaling violations in the hatted quantities are small and amenable to WχPT.

C. Summary of χPT formulae
Our central fits below will include terms through NLO in the GSM power counting. These formulae depend exclusively on the dimensionless quantities we have introduced in the previous subsections. The NLO expressions for the pseudoscalar masses of the two representations are while the expressions for the decay constants are The one-loop chiral logarithms enter as where the dimensionless mass-squared of the singlet NG boson is defined bŷ This corresponds to the LO result of Ref. [12], rescaled by t 0 . Further technical details related to the ζ and our conventions for the conserved axial current appear in Appendix C 1.
The most important parameters in the expressions above are the LO LECS of the continuum two-representation theory (rescaled by . The general form of the analytic NLO continuum terms was discussed in [12]. Because we do not have enough independent quantities to distinguish the individual NLO LECs, we instead consider L M rs and L F rs as the parameters for the fit. Finally, the various W parameters account for the NLO analytic terms of WχPT in the GSM power-counting scheme. Overall, these formulae contain 21 undetermined parameters, which we will fit below using 172 correlated points of data: four data points for each of our 43 ensembles.
We have not presented NLO fit formulae for the mass and decay constant of the singlet NG boson ζ. We do not make use of these formulae in this work because we have not calculated fermion-disconnected diagrams, which is technically challenging, and so we do not have direct access to the singlet sector. Nevertheless, through their dependence on ∆ ζ , virtual ζ loops contribute to the masses and decay constants of the other NG bosons at NLO. In the next section we will explore what can be learned about the singlet sector from this effect.
Another interesting quantity is the chiral condensate in each representation. At lowest order in χPT (equivalently, in the corresponding chiral limit,m r → 0), the fermion condensate per flavor is given byΣ Instead of measuring the condensates directly-a formidable task with Wilson fermionswe will make use of Eq. (3.19) to extract their values in the (double) chiral limit from our analysis of the pseudoscalar masses and decay constants.

A. Masses and decay constants
We begin with the pseudoscalar mesons, which become NG bosons in the chiral limit. For a first look, we plot in Fig. 3 the squared massesM 2 P r . The sextet massM P 6 is plotted against the AWI massm 6 of the sextet fermion, ignoring the dependence on the fundamental fermion massm 4 , and likewise forM P 4 , plotted againstm 4 . As expected from leading-order chiral perturbation theory, the overall behavior of each squared mass is approximately linear. One supposes that the scatter around the straight lines is due to the hidden dependence on the other fermion mass, as well as corrections from NLO and from lattice artifacts. We will examine this hypothesis shortly.
The pseudoscalar decay constants are defined by at zero spatial momentum, which is the convention that gives F π 130 MeV in QCD. We calculate F P r with the procedure described in Appendix B 2, renormalizing according to the analogue of Eq. (B5). We plot the (rescaled) decay constantsF P r in Fig. 4. The data show a steady rise withm r . The same qualitative behavior is seen in QCD, where the pion decay constant is an increasing function of the quark mass. We have presented in Sec. III C the predictions of χPT in NLO for pseudoscalar observables. We conduct a joint fit of the four observablesM 2 P r andF P r to the NLO formulae of Eqs. (3.13)-(3.16). On each ensemble, we use single-elimination jackknife to construct the 6 × 6 correlation matrix among pseudoscalar masses, decay constants, and AWI masses of the fermions. These correlation matrices enter into the χ 2 that is minimized for the fit. We do not include correlations with the flow scale t 0 , which has negligible error compared to the other quantities we extract.
The full NLO fit to 21 parameters and 172 − 21 = 151 degrees of freedom gives χ 2 /DOF = 0.48. Table I  are also significant. From an empirical perspective, these four NLO Wilson terms form a necessary minimal set of artifact terms for modeling the data.
Figures 5 and 6 illustrate the sizes of the Wilson artifacts (red) emerging from this fit. In these figures, the "corrected" data (dark blue) result from subtracting the lattice artifacts from the data (light blue), allowing us to extrapolate to the continuum limit,â → 0. The corrected data follow fairly well the tree-level formula for the pseudoscalar masses and the continuum NLO result for the decay constants, respectively, both indicated by green bands. (The bands represent 1σ in the fit parameters.) In order to display a smooth curve for the continuum NLO result for the decay constants, we have included only the samerepresentation terms when drawing the green band (indicated by "continuum NLO SREP" in the figure). The remaining scatter and deviation in the subtracted data (dark blue) is evidence of coupling between the representations. Table I demonstrates that all five leading-order LECs are well-determined by the NLO fit. We note that the singlet decay constantF ζ is larger thanF 4 and similar in size toF 6 . Because measurement of chiral logarithms is known to be a difficult task in QCD studies, we return to the question of the stability of this result below. different from zero at the 2σ level, The converse influence of the fundamentals upon the sextets follows from exchanging (4 6). The ratios L M 64 /L M 66 , L F 64 /L F 66 , andW M 64 /W M 66 are all consistent with zero. Despite the large uncertainties, this suggests that the sextets influence the fundamentals significantly, while the converse is not true. The same qualitative conclusion is also evident, for instance, in the NLO continuum behavior of the decay constants. Figure 6 shows that subtracted data (dark blue) are, to good approximation, a smooth function ofm 6 only. In contrast, the corresponding fundamental result in Fig. 5 (also in dark blue) exhibits a conspicuous jaggedness, indicating important dependence on the sextet fermion mass.

B. Stability of the NLO fit
In this subsection we explore the stability of the NLO fit. First, since we are using priors to ensure convergence of the non-linear fitting procedure, it is important to verify that our results were not biased by them. To this end, we have redone the fit using the results of the first fit as initial guess, while multiplying the width of all priors by 10. Figure 7 shows the results of both fits for the 5 LO LECs in the two lines at the bottom. The results are indistinguishable, indicating that the LO LECs were not influenced by the priors. (The same is also true for the NLO LECs.)  The range of lattice spacings we explore is 0.4 â 2 1.1 . Measurement of the LECs also provides information about the convergence of the chiral expansion. With our convention for the decay constant, the expansion parameters of continuum χPT are ξ r ≡ 2B rmr /8π 2F 2 r . With the central-fit values for the LECs, these fermion masses correspond to the following ranges for the expansion parameters, We see that the the maximum of the sextet expansion parameter ξ 6 is smaller by a factor of 2.5 than the fundamental expansion parameter ξ 4 . The main reason is thatF 6 is significantly larger thanF 4 , as might be expected based on the relative dimension of the two representations (see Sec. VI A). It is quite plausible that ξ 6 is sufficiently small that the expansion in m 6 converges well over our entire ensemble set. The same may not be true for ξ 4 , whose value can be as large as 0.5. In the next three lines of Fig. 7 we study the influence on the LO LECs of dropping ensembles at the high end of them 4 range:m 4 > 0.09, > 0.07, and finally > 0.05. We see that truncating our data set has only a modest effect on theF r andB r parameters. On the other hand, since we only obtainF ζ through NLO logarithms, it is not surprising that the increase in the error bar ofF ζ is much more pronounced. Indeed, when we restrict tô m 4 < 0.05,F ζ is only 2σ away from zero. The next two lines in Fig. 7 investigate the possible influence of finite-volume effects on our central analysis; further discussion appears in Appendix B 4. The minimum cutoff on M P L in data used in the central fit is varied from its initial value of 4.0 in our main analysis to 4.5 and 5.0, excluding more data that may be expected to have the largest finite-volume contamination. Finally, in the top line we repeat our fit with all V = 16 3 × 18 ensembles excluded from the analysis, in order to test for systematic effects in our correlator analysis due to the smaller time extent. No significant change to our results is seen in any case.
The main systematic uncertainty about this non-QCD system is the neglect of N 2 LO corrections. We do not really know how high can we go in ξ 4 and ξ 6 if we want these corrections to remain below a certain level. While our stability tests give us some insight, we do not have enough data for a quantitative study of N 2 LO. Nevertheless, we take the smallness of ξ 6 and our stability tests onm 4 as evidence that the data are in the regime where NLO ChPT applies, even if we do not have enough information to quantify the corresponding systematic error.

C. The singlet Goldstone boson ζ
As explained in Sec. III C, the chiral fit in the fundamental and sextet sectors allows us to probe the ζ meson as well. We examine its mass in the chiral-sextet limit,m 6 → 0. Figure 8 showsM 2 ζ , constructed using Eq. (3.18) and the parameters of the central fit, in the continuum (â → 0) limit, as a function of the massm 4 of the fundamental fermions. The figure shows that the singlet boson is consistently lighter than the pseudoscalar of the fundamental sector in this limit.
We can make a conservative prediction regarding the ζ mass as follows. As we have just explained, we do not know how largem 4 can be while keeping the N 2 LO corrections below, say, 10% or 20%. Lowering the maximal value ofm 4 raises the uncertainty inF ζ , as seen in Fig. 7. Still, even if we lower the maximal value ofm 4 so as to, say, double the uncertainty ofF ζ , we would still find that M 2 ζ < M 2 4 at the 1σ level. The chiral-sextet limit is interesting for composite-Higgs models. In many models, including those proposed by Ferretti and Karateev [8], the symmetries of the Standard Model are embedded into the unbroken global symmetries, so that neither the fundamental nor the sextet fermions are required to be strictly massless. Nonetheless, successful models are likely to have very light sextet fermions, because a large sextet mass would prevent the Higgs field from condensing even after the generation of a potential from the coupling of the Higgs to Standard Model fields.

A. Masses and decay constants
We now turn to our results for vector masses and decay constants. Vector-meson decay constants appear in the literature with a variety of conventions. We define F V r to have units of energy, literature on precision electroweak observables, for example Ref. [39]. Figures 9 and 10 show results forM V r andF V r , respectively, each plotted against the fermion massm r in the same representation. As before, we measure all quantities in units of t 0 . The data for these plots are listed in Tables XI-XIII. Both quantities shows a modest, plausibly linear rise against the fermion mass, albeit with a large spread. Figure 11 shows the ratio of the pseudoscalar and vector masses, M P r /M V r , again plotted against the fermion massm r in the same representation. This ratio is greater than or equal to a half for all but the smallest masses. Because the decay V → P P is p-wave, the vector where k min = 2π/L is the minimum nonzero momentum. This condition is satisfied for both representations on all of our ensembles, so the vectors are indeed stable.
We modelM V r andF V r as linear functions of the fermion mass in the same representation and of the lattice spacing, for example, For this analysis, we restrict ourselves to the 30 ensembles for which we were able to measure the vector decay constants (see Tables XI-XIII

B. Vector meson dominance and the KSRF relations
The pseudoscalar and vector decay constants are related through the hypothesis of vector meson dominance (VMD). Kawarabayashi, Suzuki, Riazuddin, and Fayyazuddin (KSRF) showed long ago [40][41][42] that VMD leads to the prediction independent of representation. Figure 14 shows the ratio F V r /F P r in each representation, after subtracting lattice artifacts. The KSRF prediction is qualitatively successful. (In QCD, the experimental value is roughly 1.66.) Another result of KSRF is that the on-shell coupling constant g VPP mediating the decay of a vector into two pseudoscalars is given by We plot this ratio in Fig. 15. As already noted, in our ensembles the vector meson is stable. Nevertheless, we may use the KSRF result as a phenomenological estimate for the behavior close to the chiral limit. Using the tree-level formula for the V → P P decay width in the limit where M P r M V r , we can estimate the the mass-to-width ratio for each vector resonance, We want to put our results in context with comparisons to QCD. Large-N c counting (for N c colors) is a way to do that. Any quantity Q is expected to scale across N c as where p is some characteristic exponent determined by large-N c considerations, and the Q i are a set of expansion coefficients. Before we get to the (limited) comparisons between different N c 's we can make, our theory gives us a unique opportunity to compare the expansion coefficients for different representations. More specifically, if we neglect all the subleading corrections, our data allow us to compare the leading expansion coefficient Q 0 between the fundamental and two-index antisymmetric representations, for various obesrvables.
We start with meson masses, which are predicted to be independent of N c [p = 0 in Eq. (6.1)]. Figs. 3 and 9 reveal near-independence of representation of the pseudoscalar and vector masses when plotted against the corresponding fermion mass. This is further supported by the near equality ofB 4 andB 6 in Table I. We conclude that Q 0 is roughly independent of representation for the pseudoscalar and vector meson masses. Decay constants scale as √ N c for single-index fermions and as N c for two-index fermions. As for the leading expansion coefficient, a possible guess might be that the product N p c Q 0 follows the leading large-N c behavior of (dim r) 1/2 . This would imply that there exists a constant c such that Q 0 ≈ c for mesons made of fundamental-representation fermions, while Q 0 ≈ c/ √ 2 for mesons made out of fermions in the two-index antisymmetric representation. For N c = 4 we thus expect F 6 /F 4 ≈ N c /2 = √ 2. The ratio F 6 /F 4 is plotted in Fig. 17 for the pseudoscalar and vector mesons, showing good (perhaps too good) agreement with this expectation. Another consequence is that the ratio F V r /F P r is expected to be roughly independent of representation r, in agreement with Fig. 14 and the KSRF relation.
In this context, we can also compare the fermion condensates in the two representations. We can use the results of the chiral fit to calculate these (per flavor), each in its corresponding chiral limit. Using Eq. (3.19) we findΣ 6 /Σ 4 = 2.4(3) for the ratio of condensates in the (double) chiral limit, reflecting our empirical findings thatB 6 /B 4 1 andF 6 /F 4 √ 2.
Finally, we compare our results for the ratio M V r /F P r to the SU(3) case. We expect that this ratio for either SU(4) representation will be smaller than its value in QCD, which is 770 MeV/130 MeV = 5.9. This is borne out by Fig. 15. At a more quantitative level, large-N c scaling predicts this value to be 5.9 × 3/4 ≈ 5.1 for the fundamental representation, in good agreement with the data. If we regard the fundamental representation of SU(3) as the two-index antisymmetric one, we obtain a large-N c prediction of 5.9 × 3/4 ≈ 4.4 for the sextet. While the average value of this ratio is smaller for the sextet, one can say that this prediction is consistent with the general trend of our data. If we add to this the KSRF relation, we find correspondingly smaller values for the width to mass ratio of our vector mesons compared to the physical ρ meson.

B. Conclusions
In this paper we have described the low-lying mesonic spectrum of SU(4) gauge theory coupled to dynamical fermions in the fundamental and sextet representations. These multirep simulations are the first of their kind. Our choice of this theory was inspired by its close similarity to a composite-Higgs model first studied by Ferretti [7].
Our analysis focused on the masses of the pseudoscalar and vector states and their associated decay constants. Using the extension of χPT that accounts for the discretization errors of Wilson fermions, we carried out a global fit of the pseudoscalar masses and decay constants of the two representations, to NLO in the GSM power-counting scheme. We found significant lattice artifacts, which we were able to subtract, obtaining predictions for continuum-limit values. Our chiral fit provides mild evidence for coupling between the two fermion representations, a novel feature of multirep theories.
Through both the mass terms and the Wilson terms, our lattice setup incorporates the expected symmetry breaking patterns: SU(2) L × SU(2) R → SU(2) V in the fundamental sector, and SU(4) → SO(4) in the sextet sector. We did not carry out a dedicated study of alternative symmetry breaking patterns. Still, the success of the chiral fits provides some confirmation that the above symmetry breaking patterns are the right ones.
The main theoretical uncertainty of our chiral fits concerns the size of N 2 LO effects. Thanks to a large decay constant, the chiral expansion converges quickly in the sextet sector, supporting the hypothesis that N 2 LO effects are small in this sector. In the fundamental sector the chiral expansion converges more slowly. Hence, keeping N 2 LO effects below a certain comfortable level might require the exclusion of ensembles wherem 4 is on the high side. More quantitative statements cannot be made given our data. KSRF identify this quantity with the coupling g VPP . In QCD, this ratio is roughly 5.9.
The correlation functions that we calculated probe directly the pseudoscalar states made purely of fundamental or purely of sextet fermions. This is reflected in the stability of the LO parametersF r andB r if we constrain the maximal value ofm 4 to successively smaller values. The last LO parameter, the decay constantF ζ of the axion-like singlet NG boson, is not well-determined because we have not calculated propagators in the ζ channel. The ζ meson does contribute through virtual loops to the correlation functions we have studied. Accordingly, our fits depend onF ζ , but only through NLO logarithmic terms.F ζ is more sensitive to the upper limit onm 4 ; as a result, so is our prediction for the mass of the ζ boson. Nonetheless, we have argued that the ζ is lighter than the fundamental-sector NG bosons, M ζ < M 4 , in the sextet-chiral limitm 6 → 0, a limit which is interesting for the phenomenology of Ferretti's model. In a full composite-Higgs model, however, the masses of all pseudoscalar states can receive important corrections from the couplings to Standard Model fields.
In modeling our results for the vector mesons, we found that the ratio of pseudoscalar to vector decay constants agrees well with the KSRF result based on vector meson dominance. As discussed in Sec. VI A, comparing the KSRF prediction for the decay rate of the vector meson in the chiral limit to the QCD case shows reasonable agreement with large-N c counting.
Although our estimates for Γ V /M V depend on the well-motivated but non-rigorous assumption of vector meson dominance, the resulting narrowness is almost certainly generic. In large-N c , the widths of mesons made of fundamental-representation fermions scales as 1/ √ N c and thus they become narrower as N c increases. Insofar as large-N c provides the cleanest explanation for the narrowness and existence of mesons in QCD, the vector mesons should become narrower in theories with more colors. We proposed that in multirep theories, the generalization of √ N c is (dim r) 1/2 , a hypothesis supported by our data. This result is good news for phenomenologists looking to constrain models like the Ferretti model, since narrower states typically provide clearer signals in collider data.
As we have mentioned, we are also exploring the phase diagram of this multirep theory [21]. We have been looking for-and not finding-scale separation between the representations in the confinement and chiral transitions. We are also studying the baryon spectrum, a particularly interesting sector of the theory given its connection to top-quark physics and partial compositeness in the Ferretti model.
Other interesting avenues for the future work in this model (or multirep composite Higgs theories more generally) include quantities related to the Higgs potential. The contribution of the Standard Model's gauge fields to the Higgs potential, Π LR , is conceptually identical to the physics of electromagnetic splittings between pions in QCD and has been the subject of a recent pilot study on the lattice [43]. The top-quark contribution to the Higgs potential is considerably more challenging [44,45].

ACKNOWLEDGMENTS
Computations for this work were carried out with resources provided by the USQCD Collaboration, which is funded by the Office of Science of the U.S. Department of Energy; with the Summit supercomputer, a joint effort of the University of Colorado Boulder and Colorado State University, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University; by the Cy-Tera Project, which is co-funded by the   1.093 (9) 0.0422 (7) 0.0203 (10)  22 1.135(9) 0.0279 (11) 0.0251 (12)  23 1.128 (24) 0.0243 (7) 0.0326 (7)  24 1.132 (12) 0.0345 (8) 0.0323 (14)  25 1.100(10) 0.0331(5) 0.0325 (5)  26 1.111 (9) 0.0228(6) 0.0381 (8)  27 1.174(10) 0.0556 (7) 0.0220 (9)  28 1.095 (12) 0.0282 (7) 0.0427 (7)    In calculating correlation functions, we use a smeared source operator on the t = 0 time slice while using both point and smeared operators at the sink, projected onto zero spatial momentum. Smearing is done after fixing to the Coulomb gauge, always with smearing radius r 0 = 6a. On the large lattices we use antiperiodic boundary conditions in time for the fermion propagators. For the V = 16 3 × 18 ensembles, on the other hand, we superimpose propagators computed with periodic and anti-periodic boundary conditions, which effectively doubles the temporal size of the lattice-a technique sometimes called the "P+A trick" (see Ref. [46] and references therein).
After restricting each correlator to a range [t min , t max ] we find acceptable fits with singleexponential models, that is, without inclusion of excited states. In each representation, we extract the fermion mass m r from the axial Ward identity (2.3) via joint fits to the AP and P P correlators with a point sink.
We use the publicly available Python packages lsqfit [47] and gvar [48] for nonlinear fitting and classical error propagation. When computing ratios of quantities derived from different fits, we use single-elimination jackknife to propagate errors including correlations.
For each correlator, our fitting procedure is as follows. First, we vary the initial and final times [t min , t max ] used in the fits, amounting to a grid search over possible range fits. The best fit is chosen automatically using a criterion from the QCD literature with a preference for small χ 2 /dof, large fit ranges, and well-determined fit parameters [49]. We maximize where p is the unconstrained p-value, N dof denotes the number of degrees of freedom in the fit, and σ pn denotes the statistical error in the nth fit parameter. Although this criterion is ultimately arbitrary, it coincides with intuition about which fits ought to be considered good and removes subjective bias. We confirm that masses emerging from this procedure are consistent with expectations from effective mass plots; a representative comparison is shown in Fig. 18. We have also experimented with a "two-state" double-exponential ansatz, observing no significant changes in the ground-state masses within the combined statistical and systematic errors estimated using this procedure.
For an estimate of the systematic uncertainty associated with our fit-choice procedure, we compute the spread in the model parameters emerging from all nominally good fits satisfying Q ≥ 0.1. We then combine the statistical and systematic uncertainties conservatively using The systematic error assigned by this procedure is often comparable to the statistical error, and is occasionally significantly larger. The error estimates for the fermion masses m r in Tables V-XIII include this fit-range systematic.

Decay constants and operator renormalization
The lattice operators appearing in the correlation functions are subject to finite renormalization in order to obtain the continuum-normalized operators that are required for determination of decay constants. We carry out this procedure in lattice perturbation theory including tadpole improvement; the procedure is described in Appendix D. The explicit relationship between lattice and continuum operators is given by Eq. (D28).
Simultaneous fits to the smeared-source, point-sink (s, p) and smeared-source, smearedsink (s, s) correlation functions allow us to extract the mass, decay constant, and smeared amplitude. For example, the V V correlators in representation r give us the vector decay constant F V r defined in Eq. (5.1). The fit functions are giving the vector mass M V r and the point and smeared amplitudes A p V r and A s V r , respectively. In order to obtain decay constants with continuum normalizations, we apply the renormalization factors of Eq. (D28). The result is

Fermion mass determination and κ c r
To determine the critical values κ c r , which enter into the field normalization for decay constants defined in Appendix B 2, we perform a global fit to the AWI fermion masses in units of the flow scale t 0 as given in Tables V-VII. We use the model function and similarly for √ t 0 m 6 (with a separate set of coefficients). We find that these terms, which are a subset of all possible combinations of the bare parameters {β, κ 4 , κ 6 } through quadratic order, are sufficient to provide reasonable fit quality.
Since we are interested in the regions where m r → 0, we use only those ensembles that have √ t 0 m r < 0.08, a value determined empirically by inspecting our data for deviations from the simple analytic behavior of Eq. (B6). Our fits give χ 2 /dof of 16/24 and 21/24 for fitting √ t 0 m 4 and √ t 0 m 6 , respectively. The resulting κ c curves at two β values are shown in Fig. 19.
As noted above, because κ c is determined by extrapolation, we do not probe the existence of a possible Aoki phase.

Study of finite-volume corrections
All of our ensembles satisfy the criterion M P r L > 4 for both pseudoscalar meson masses. Using Eq. (3.18) with the results of our central fit, we further verify that M ζ L > 4 in all cases. Although this cut is known to provide a useful threshold for the suppression of finite volume effects in QCD, since we are studying a new system a more cautious treatment is worthwhile. Here we present three different analyses and arguments which lead us to estimate that finite-volume effects are no more than a few percent for our data, and are thus not resolved within the uncertainty of our results. One particular source of finite-volume effects can be the freezing of the evolution of topological charge [50,51]; we measure the topological charge Q on our ensembles using the Wilson flow to smear the gauge fields out to t/a 2 = 5.0 and find an acceptable distribution of Q in all cases.
First, to obtain a theoretical estimate of the expected size of finite-volume effects, we consider the size of leading-order finite-volume correction to tadpole diagrams in chiral perturbation theory [9,33,52]. The dimensionless figure of merit for this effect is 2I 1 (M P r , L)/F 2 P r , where This quantity gauges the effect of mesons interacting with their finite-volume image points. In Fig. 20 we plot 2I 1 (M P r , L)/F 2 P r for all of our ensembles, with each representation r plotted against the correspondingm r . We therefore expect that finite volume corrections do not exceed a few percent in the ensembles of this study. This formula assumes the applicability of chiral perturbation theory, which requires that F P r L 1; over all of our ensembles we find that F P 4 L 1.3 and F P 6 L 1.9.
Second, we recalculate the observables M P r , F P r , and t 0 on ensembles with several spatial volumes at two sets of bare parameters, (β, κ 4 , κ 6 ) = (7.75, 0.128, 0.131) and (7.75, 0.129, 0.1308). These bare parameters match the production V = 16 3 × 32 ensembles 36 and 37 as listed in Table III. Four ensembles hold N t = 32 fixed and vary the spatial volume as N s = 12, 14, 16, 18; the fifth ensemble at each point has N s = 24 and N t = 48.
Results of this test are shown in Figs. 21-23 below. For both sets of bare parameters, all observables down to the smallest N s = 12 are seen to be within ±5% of the central value obtained on N s = 24, and within 2-3% for N s = 16 which is the smallest spatial volume used in our central analysis.
Finally, we have included explicit variation of the finite-volume cut on pseudoscalar meson masses (i.e. minimum cut on M Pr L) in the stability analysis of our central chiral fit, as presented in Sec. IV B and Fig. 7. We also consider the effects of the finite temporal direction by cutting the N t = 18 ensembles out of the analysis. All of the fit results are seen to be stable at the one-sigma level as we vary the finite-volume cut. We conclude that finite-volume effects are not significant in our results at the level of precision we obtain. The conserved axial current is given by As usual, the normalization of U(1) currents is arbitrary. We normalize the individual axial currents J (r) 5µ such that all right-handed fields have unit charge. The ratio q 4 /q 6 is then fixed by the group traces of the two representations. For the normalization of J 5µ we adopt the same prescription as in Ref. [12]. The resulting charges are These charges were used in the χPT formulae of Sec. III C. Tracing Eqs. (3.13)-(3.16) back to the general NLO expressions of Ref. [12], one can check that they only depend on the ratios q r /F ζ . These ratios are independent of the choice of normalization for the axial current, because a rescaling of J 5µ by a factor λ implies a rescaling by the same factor of both the charges q 4 and q 6 and of the singlet's decay constant F ζ . All other LECs in Eqs. (3.13)-(3.16) are independent of this rescaling (for more details, see Ref. [12]). Of course, once the normalization of the charges is fixed according to Eq. (C2), the normalization of F ζ is fixed as well.

AWI mass and Wilson chiral perturbation theory
In this appendix we review the proof of Ref. [36] that the mass defined by imposing the axial Ward identity, m AWI , is equal to the shifted mass m shifted , which is the mass parameter occurring in the LO lagrangian of WχPT, up to higher-order corrections. For simplicity, we will consider the GSM power counting used throughout this paper.
A nice feature of the GSM scheme is that the LO lagrangian takes the same form as in the continuum. The reason is that the only new operator at O(a) is the Pauli term aψσ µν F µν ψ, which has the same chiral transformation properties as the fermion mass term. The Pauli term enters the Symanzik action with a coefficient that depends linearly on the parameter c SW of the clover term in the Wilson action. After the transition to the chiral effective theory, the non-derivative terms in the LO lagrangian take the form where the mass term has been mapped to the first term on the right-hand side, and the Pauli term to the second. χ is the usual spurion of continuum χPT, andÂ is a new spurion with the same chiral transformation properties as χ. The "expectation values" of the spurions are [53]  where B and W 0 are low-energy constants (LECs). Substituting back into Eq. (C3) gives where the shifted mass is defined by Eq. (3.10). For brevity, in the rest of this appendix we denote the shifted mass as m.
As explained in Sec. II B, in our numerical simulations we define the fermion mass m AWI by imposing the axial Ward identity, Eq. (2.3). Before we can use WχPT, we need to know the relation between m AWI and the shifted mass m of Eq. (3.10), in terms of which the chiral expansion is done.
This relation was analyzed carefully in Ref. [36] for the case of two-flavor QCD. The pseudoscalar density that was considered there, and which is also used in our work, is the usual local density, 3 P a loc =ψγ 5 T a ψ .
For the axial current, Ref. [36] considered where ∂ µ stands for a lattice derivative; the local axial current is For the purpose of this discussion, we may consider c A in Eq. (C7) as a free parameter. In our numerical simulations we use the naive axial current, which corresponds to c A = 0. To first order in the pion field, the lattice operators are mapped to the effective theory according to where π a is the effective pion field. The precise form of the O(a) corrections may be found in Ref. [36]. In both equations, they depend linearly on the clover parameter c SW . In addition, the O(a) correction in Eq. (C10) depends linearly on c A . Plugging this into Eq. (2.3) and using the LO pion mass, given by M 2 = 2Bm, we find To the order we are working, in general one expects also an O(m 2 ) correction. This correction vanishes, however, because the continuum theory satisfies m AWI = m identically. While the derivation of Ref. [36] was given for a complex representation, a similar argument applies to real (or pseudoreal) representations. Equation (C11) is robust in that changing the clover coefficient c SW or changing the parameter c A in Eq. (C7) will change the O(am) corrections, but will not affect the simultaneous vanishing of m AWI and the shifted mass m. This feature is disrupted only by O(a 2 ) effects, which is as it should be. Indeed, as shown in Ref. [28], depending on the sign of a particular O(a 2 ) LEC, in the region where m ∼ a 2 one either encounters the Aoki phase, or a first-order discontinuity line at which the pion mass reaches a non-zero minimum.
As a corollary of Eq. (C11), we may use the value of m AWI , taken from the numerical simulations, for the shifted mass m. At the order we are working, the differences between the two are absorbed into a redefinition of NLO parameters of the chiral expansion.

Appendix D: Gluon propagator and perturbative calculations for nHYP links
To perform one loop perturbation theory for the NDS action [20], we need the gluon propagator. This appendix describes its construction and gives perturbative results for current renormalization factors.
Normalized hypercubic links V x,µ are constructed from the dynamical gauge field U x,µ via three successive smearing steps [16,17]. Each step uses a weighted sum over staples, which is then reunitarized. Explicitly, The reunitarization operator P is defined as where The best way to understand the smearing is to go in reverse order. The staple sum extends into a different direction at each smearing step, such that each fat link V x,µ depends on a particular thin link U y,ν , if and only if there exists a hypercube to which both V x,µ and U y,ν belong.
The dislocation-suppressing action adds a new term to the pure-gauge action S g where the new term is In practice we fix the α i 's and take γ 1 = γ 2 = γ 3 = γ = zβ where z is held constant. The weak-coupling expansion to be sketched below gives the bare coupling g 2 0 as 1 To construct the gluon propagator we need to compute the gauge action in quadratic order. This is pretty standard; a good reference is Ref. [54]. We take the lattice action, S = a 4 x L(ψ, ψ, U ), and replace the link field by an expansion in terms of gauge fields where A µ (x) = a (λ a /2)A a µ gives the decomposition into color components. There is an identical expansion for the fat link, which we write as V µ (x) = exp[igaB µ (x)].
The action has an expansion in powers of A. In terms of the integral over the fourdimensional Brillouin zone, and the vector potential, the free gauge boson action is For the gauge boson, D ab µν = δ ab D µν . Just as in a continuum theory, the gauge boson action cannot be inverted to give the propagator without fixing the gauge. A conventional choice for a gauge fixing term is [introducingk µ = 2/a sin(ak µ /2)] Then the gauge boson propagator is found by solving the field equation We simply do this numerically, inverting the four by four matrix for each k value.
To perform the perturbative expansion of the NDS action, we look at each term in turn. Consider Q x,ρ;ξ = Ω † x,ρ;ξ Ω x,ρ;ξ .
Multiplying this out, we find that Q is a sum of loops of perimeter four and perimeter six, labeled P and E respectively. The planar E loops, E µν , are 1 × 2 loops extending in the ±ν direction from site x. Hence, Expanding each term out in terms of A's, we discover that Q −1 = 1 − S Q where now the expressions are quadratic functions of the gauge fields. Slightly abusing notation, we write P and E for the objects made of A µ 's, so that Nearly identical results obtain for the other two Q's, with several small exceptions. First, the analogs of the plaquettes are built of one thin link A µ (x) while the other three links are fattened. Second, the perimeter-6 links are built entirely of fat links. Finally, in addition to the 1 × 2 E µν plaquettes, there are "chair" plaquettes C µνρ which are "folded" about the µ axis. They extend in four directions (±ν, ±ρ).
Our simulations are done with z = 1/125. At that value of z, the perturbative properties of the NDS gauge action are almost indistinguishable from those of a pure Wilson action. Here are three examples.
First, the plaquette has an expansion Tr U plaq /N c = 1 − g 2 pC F where C F is the quadratic Casimir for fundamentals and p is a constant. For the Wilson action, p = 0.5. This expression is often replaced by − ln 1 N c Tr U plaq = g 2 pC F .
This defines a coupling g 2 in the so-called α V scheme, because the potential is written as The scale of the coupling is set by the Lepage-Mackenzie [55] prescription, log q * = d 4 q log qI(q) d 4 q I(q) . (D27) Results for p and q * for several values of z are given in Table XIV. With the gluon propagator in hand we can immediately compute the static Coulomb potential. With our conventions, the continuum potential is V (r) = 1/(4πr), and so plotting the rescaled lattice potential 4πrV (r) immediately exposes the lattice artifacts of a particular action. We show results for this quantity in Fig. 24 for z = 1/125, plotted together with the Wilson action result.
Last, we have the renormalization factors for currents. Calculations of matrix elements require a conversion to continuum regularization. We adopt the old tadpole-improved procedure of Lepage and Mackenzie [55]. In this scheme a continuum-regulated fermionic bilinear quantityQ with engineering dimension D [we have in mind finding the M S (modified minimal subtraction) value at scale µ] is related to the lattice value bȳ with of α V is carried out by numerical integration of the two-loop β-function, where the required coefficients for a theory with multiple fermion representations are [58] β 0 = 1 3(4π) 11C 2 (G) − 2 r N r C 2 (r) (D31) Here N r denotes the number of Dirac flavors in representation r, while T (r) and C 2 (r) are the standard trace and Casimir invariant for each representation. For our SU(4) multirep theory, we obtain β 0 = 53/(24π) and β 1 = 1531/(192π 2 ).