Finite-temperature phase structure of SU(4) gauge theory with multiple fermion representations
Venkitesh Ayyar Thomas DeGrand Daniel C. Hackett William I. Jay Ethan T. Neil Yigal Shamir Benjamin Svetitsky
























Department of Physics, University of Colorado, Boulder, Colorado 80309, USA


Department of Physics, University of Colorado, Boulder, Colorado 80309, USA


Department of Physics, University of Colorado, Boulder, Colorado 80309, USA


Department of Physics, University of Colorado, Boulder, Colorado 80309, USA

ethan.neil@colorado.edu
Department of Physics, University of Colorado, Boulder, Colorado 80309, USA
RIKEN-BNL Research Center, Brookhaven National Laboratory, 
 Upton, New York 11973, USA


Raymond and Beverly Sackler School of Physics and Astronomy,
Tel Aviv University, 69978 Tel Aviv, Israel


Raymond and Beverly Sackler School of Physics and Astronomy,
Tel Aviv University, 69978 Tel Aviv, Israel

 


We investigate the phase structure of SU(4) gauge theory with the gauge field simultaneously coupled to two flavors of fermion in the fundamental representation and two flavors of fermion in the two-index antisymmetric representation.
We find that the theory has only two phases, a low-temperature phase with both species of fermion confined and chirally broken, and a high-temperature phase with both species of fermion deconfined and chirally restored.
The single phase transition in the theory appears to be first order, in agreement with theoretical predictions.







Finite-temperature phase structure of SU(4) gauge theory with multiple fermion representations
    Benjamin Svetitsky
    December 30, 2023
==============================================================================================













§ INTRODUCTION

	
	We have been performing numerical simulations of SU(4) gauge theory coupled to
two flavors of  Dirac fermion in the fundamental (quartet) representation and two flavors
of Dirac fermion in the two-index antisymmetric (sextet) representation, which is a real representation.
These studies have been motivated by the use of a related theory—with three flavors of 
Dirac fundamentals and five Majorana sextets—as a model for a composite Higgs boson alongside a partially composite top quark <cit.>. 
Our previous work <cit.> was concerned with the mesonic and baryonic properties of this system.
Here we describe our studies of its finite temperature behavior.

The presence of multiple fermion representations (a multirep theory) opens the possibility of dynamical scale separation between the confinement and chiral transitions for each representation.
Scale separation may arise from a variety of different mechanisms.
One possibility is separation between the chiral symmetry breaking scales associated with the two representations, as in the “tumbling" scenario <cit.> in which chiral symmetry is broken spontaneously at a hierarchy of different mass scales. 
Quenched studies on very small lattices in the early 80's indicated the existence of 
separated chiral transitions for different fermion representations <cit.>,
but it is not known whether these results persist in the presence of dynamical fermions.
It also may be possible for the scales of the chiral and confinement transitions to be 
different.
Previous work with dynamical fermions has pointed to this 
possibility <cit.>, but it is possible that the theory explored in that work—the SU(3) theory with adjoint fermions—is infrared conformal rather than confining <cit.>.
A final possibility is that the confinement transitions of different representations can be separated: if center symmetry breaks in several stages, there may exist phases where some representations of charge are deconfined while others remained confined. 

Our numerical data lead us to the conclusion that 
there is only a single  finite-temperature phase transition in this theory, 
with the characteristics of chiral restoration and deconfinement for both fermion species.
We find this to hold in 
limit theories where one or the other fermion species is decoupled, as well as in the full theory coupled to both species simultaneously.
The sextet-only theory has an order parameter which characterizes its phase,
and so it presumably possesses a real phase transition line. We could not determine its order.
The fundamental-only theory has no order parameter for confinement and appears to exhibit 
crossover behavior.
For the full theory, we find only first-order transitions in the region we explore.
We summarize our knowledge in the “Columbia plots” (Figs. <ref> and <ref>) below.

The paper is organized as follows:
In Section <ref>, we discuss the theory, its symmetries, and the expected 
behavior of its confinement and chiral transitions.
We discuss our lattice methodology in Section <ref>.
We then present our numerical findings.  In Section <ref>, we 
examine the phase structure in two limiting cases of the theory, keeping
only fundamental fermions, or only sextet ones.
We then examine the phase structure
of the full theory with both species of fermion in Section <ref>.
We conclude in Section <ref>. A preliminary account of this work can be found in <cit.>. For lattice work on a composite Higgs model with a different gauge group, see <cit.>.

In what follows, a quantity labeled m_4, P_4, etc., corresponds to that quantity as measured for the fundamental fermions, while m_6, P_6, etc., correspond to the sextet fermions.




§ CONTINUUM THEORY

	
	


 §.§ Polyakov Loops and Multiple Representations



The usual diagnostic for confinement is the Polyakov loop, which can be constructed from gauge links in any representation of the group.
For our study, the relevant Polyakov loops are defined as

    P_4(x⃗)    = Ω_x⃗
    
    	P_6(x⃗)    = 12[
    		(Ω_x⃗)^2 -  (Ω_x⃗)^2
    	]

where

    Ω_x⃗ = ∏_t=1^N_t U_t̂(x⃗, t).

In SU(3), it is possible to write any higher-representation Polyakov loop in terms of the fundamental Polyakov loop and its complex conjugate.
The behavior of higher-representation Polyakov loops is thus completely determined if the fundamental Polyakov loop is known.
In SU(4), P_4 and P_6 are a sufficient set.

The expectation value of a Polyakov loop in a representation R measures the free-energy F_R of a static charge in that representation via

    |P_R| ∼exp[-F_R/T].

This Polyakov loop is an order parameter if the fermion action preserves enough center symmetry (invariance under global Z_N gauge transformations) to protect it in the unbroken phase.
In an SU(4) gauge theory, sextet fermions break Z_4 to Z_2.
The residual Z_2 symmetry is not enough to protect the sextet Polyakov loop, but the fundamental Polyakov loop remains an order parameter when only sextet fermions are present <cit.>.
Adding fundamental fermions breaks the Z_4 center symmetry completely, so in the full theory neither Polyakov loop is an order parameter.

Even when Polyakov loops are not order parameters, they are observed to jump from a small value to a large value as the temperature rises.
This jump signals a qualitative change in the physics of color charge screening.
It could be that this occurs at different temperatures for fermions in different representations.





 §.§ Order of Chiral Phase Transitions



Each representation of fermion could have its own distinct chiral transition.
There are three possible orders in which the chiral transitions may be encountered
when cooling from T=∞ to T=0: sextet first (as predicted by
arguments along the lines of the ones in Ref. <cit.>), 
or fundamental first, or a simultaneous transition.
As discussed below in Section <ref>, our results indicate that
the third possibility is what occurs for our system.

One of us <cit.> has carried out a Pisarski–Wilczek stability analysis <cit.> for this single phase transition, based on the pattern of chiral symmetry breaking of the two fermion species.
This method consists of examining the critical behavior of a three-dimensional effective theory of the chiral condensates of the theory—a linear sigma model.
If the effective theory has any stable fixed points, it may undergo a second-order phase transition governed by one of the fixed points, but the transition may also be first-order.
In the absence of stable fixed points, the analysis predicts a first-order transition.
This procedure generalizes straightforwardly to theories where there are multiple 
representations of fermion, as long as the chiral symmetry breaking pattern can be written as a direct product of the patterns of the single-representation sectors (up to additional U(1) factors due to non-anomalous axial symmetries).
The effective theory for the multirep theory is then simply the single-representation theories coupled together.
The prediction of this calculation <cit.>, carried out to first order in ϵ=4-d, is a first-order phase transition.

Fig. <ref> is a rough sketch of a “Columbia plot” summarizing the theoretical predictions for the nature of the finite-temperature transition in the various fermion-mass regimes (with some inputs from our results discussed below).
This sketch is made in analogy with the QCD Columbia plot, where the order of the phase transition encountered is plotted as a function of m_u=m_d and m_s. Here, we plot the order of the phase In the pure-gauge limit, the transition is first order
<cit.>.
The stability analysis predicts that the transition in the double chiral limit m_4=m_6=0 will also be first order.
First-order transitions are generically robust against small perturbations, so these transitions presumably extend into the regions around m_4=m_6=∞ and m_4=m_6=0.
High-order Pisarski–Wilczek calculations <cit.> indicate that the transitions in the massless limits of the fundamental-only and sextet-only theories can be second-order
[(m_4,m_6)=(0,∞) or (∞,0)].
The fundamental-only limit, SU(4) with two flavors of fundamental fermion, 
is similar to QCD with m_u=m_d=0 and m_s=∞.  QCD in this limit is
believed to exhibit a second order phase transition with O(4) critical exponents
(compare the discussion in Ref. <cit.>),
 so we expect the fundamental-only theory to behave similarly.
Fig. <ref> shows a second-order phase transition in either of these limits.
In the single-species limits, nonzero fermion mass will convert a second-order  chiral
transition into a  crossover.
Adding heavy fermions of the other species can leave the second-order transition undisturbed or convert it to first order, as shown, but the phase transition cannot disappear as long as one species is exactly massless.
If either single-species transition were first order,
 there would be a first order region in the corresponding corner.
In the pure sextet theory there will be a confinement transition for all values of m_6 because of the residual Z_2 center symmetry; we indicate a second-order transition there, though it could be first-order.


There are no analytical predictions for the intermediate region,
 where m_4 and m_6 are neither light nor heavy; 
we thus have no predictions for whether the first order regions connect, or whether there is an intermediate continuous crossover region.



	


§ LATTICE METHODS

	
	

 §.§ Action



For this study, we employ the same lattice action as in Ref. <cit.>.
For the fermions, we use a clover-improved Wilson action built from fat gauge links constructed by normalized-hypercubic (nHYP) smearing <cit.>.
We construct the action for the sextet fermions by promoting the smeared links to the sextet representation <cit.>.
We set both clover coefficients equal to unity, c_SW = 1, a choice motivated by results from Ref. <cit.>.
The action for the gauge sector is the usual plaquette action with gauge coupling β augmented by an nHYP-dislocation suppression (NDS) term <cit.>, constructed from the nHYP-smeared links.
As in our zero-temperature study, we fix the NDS parameter γ such that β/γ=125.
Altogether, the simulation parameter space is three dimensional: β, and two hopping parameters  κ_4 and κ_6.	





 §.§ Spectroscopy





Fermion masses m_4 and m_6 for the two representations
 are defined through the axial Ward identity (AWI),

    ∂_μ⟨ 0 | A_μ a^(r)(x) 𝒪_r(0) | 0 ⟩
    = 2 m_r ⟨ 0 | P_a^(r)(x) 𝒪_r(0) | 0 ⟩,

where a is an isospin index.
We use the local unimproved axial current A_μ a^(r) and pseudoscalar
 density P_a^(r) in each representation r.
For 𝒪_r we take a pseudoscalar source.
The axial Ward identity is a statement of current conservation and is thus local and relatively insensitive to finite-volume effects as long as we stay in a confined phase.

We extract meson screening masses in the scalar, pseudoscalar, vector, and pseudovector channels from two-point correlation functions extending in a spatial lattice direction.
This is a standard technique in a finite temperature simulation.
We construct propagators with composite boundary conditions to double the effective 
length of the lattice <cit.>.
These are called “P+A correlators” in the literature.



 §.§ Data sets


In order to search for the various possible phase transitions, we explored a wide region of the three-dimensional bare parameter space.
This required an unusually large and heterogeneous data set, summarized in Table <ref>.
To render this exploration tractable, we found it necessary to automate much of our data generation and analysis; see <cit.> for a description of our methods. 

For the full theory we focused predominantly on β = 7.4 and β=7.75,
 mostly on 12^3 × 6 and 16^3 × 8 volumes but with some additional data
 on 18^3 × 6 and 24^3 × 8 to check for finite-volume effects.
For the single representation theories, we ran only on a 12^3 × 6 volume.
We also made use of zero-temperature data to determine the lattice scale,
the fermion masses, and the pseudoscalar-to-vector mass ratio for some bare couplings near the transition.
Table <ref> in the Appendix summarizes these zero-temperature data sets.
We use the lattice scale to derive the physical temperature along the phase boundary.






 §.§ Phase diagnostics


With two species of fermion, there are (in principle) four distinct transitions to be considered, namely, the confinement and chiral transitions for each representation.
We need independent observables for each of these.




  §.§.§ Confinement transition



Polyakov loops in the  fundamental and sextet  representations are used
 to tell whether each fermion representation is confined or deconfined.
We also employ two additional quantities based on the Wilson flow <cit.>: 
the flowed anisotropy and Polyakov loops at long flow time.


The flowed observable t^2 E(t), where E is the energy density at flow time t, is commonly used to determine the scale for zero-temperature lattices.
E is typically defined by summing over all orientations of clover terms or of plaquettes.
When measured on finite temperature lattices, spatial–temporal anisotropy 
in t^2 E(t) can be employed to determine the phase <cit.>. 
Previous applications of anisotropy  have used the quantity

    Δ(t) = t^2 ⟨ E_ss(t) - E_st(t) ⟩,

where E_ss(t) and E_st(t) represent the contribution from space–space and space–time clovers.
In this work we look at a related quantity, the flowed anisotropy R_E, defined as

    R_E(t) = E_ss(t)/E_st(t) = 1+
    Δ(t) /⟨ t^2 E_st(t) ⟩

(cf. Ref. <cit.>).
In the low-temperature confined and chirally broken phase, the gauge fields are roughly isotropic and the observable R_E(t) ≈ 1 for any reasonable flow time t.
In the high-temperature deconfined and chirally restored phase, hypercubic symmetry is broken
 strongly: the temporal center symmetry is broken, while the spatial center 
symmetries are still (almost) preserved.
In such anisotropic phases, R_E(t) departs from unity even at small flow times.
In this paper, we always measure R_E(t) at flow time t/a^2 = 1.

The behavior of Polyakov loops at long flow times provides a sharp diagnostic of confinement.
For a lattice with temporal extent N_t, the flow time
 ratio c_t = √( 8 t ) / (N_ta) is a rough measure of the
 Wilson-flow smearing in the temporal direction. 
Defining “long flow time” as c_t > 1, we find that flowed Polyakov loops 
exhibit nearly binary behavior at long flow times, depending on the phase.
On deconfined configurations, volume-averaged Polyakov loops rapidly reach their maximal values max |P_R| = d(R), where d(R) is the dimension of the representation R.
On confined lattices, volume-averaged Polyakov loops wander or move only
 very slowly towards their maximal values.

All of our phase diagnostics (unflowed Polyakov loops, flowed
anisotropy, Polyakov loops at long flow time) agree everywhere in our data set, within our resolution in coupling space.
The flowed anisotropy and flowed Polyakov loops may be used to determine the phase of an ensemble without comparing it with nearby ensembles or picking some arbitrary threshold value, as is required when using unflowed Polyakov loops.
Such ensemble-local observables are better suited for automation.

For further discussion of these flow-based diagnostics, see Ref. <cit.>.




  §.§.§ Chiral transition



Because we are using Wilson fermions, we 
use an indirect probe to determine
 whether chiral symmetry is broken in each sector: parity doubling in the meson sector.
In the chirally restored phase we expect parity partners to be degenerate.
We thus examine the mass splittings between the scalar and pseudoscalar states and between the 
vector and axial vector states, for each species of fermion.




§ PHASE STRUCTURE OF LIMITING CASES

	

	

 §.§ Sextet-Only Theory: κ_4=0

    	
    	

  §.§.§ Phase structure




We begin with the gauge theory coupled only to sextet fermions.
Figure <ref> shows the behavior of our various observables 
 along a typical slice through bare parameter space.
The top panel shows the behavior of the fundamental and sextet Polyakov loops.
Dynamical sextets break the center symmetry from Z_4 → Z_2; the fundamental loop
 P_4 is an order parameter for the spontaneous breaking of the
 residual Z_2 center symmetry <cit.>
 and thus for confinement of static charges in the fundamental representation.
Although it is not an order parameter, we also examine the sextet Polyakov loop.
In the top panel, we see that all three  confinement diagnostics jump simultaneously and vary only smoothly elsewhere: there is only a single confinement transition in this theory.
The middle panel shows the  mass splittings of the sextet mesons.
The parity partners become degenerate simultaneously, indicating the restoration of chiral symmetry.
Within our resolution in κ_6, 
the confinement transition coincides with the chiral transition.



The behavior shown on the slice in Fig. <ref> is typical for this theory: everywhere we have looked in parameter space, we see a single, unified confinement and chiral transition, as in QCD.
Figure <ref> summarizes our findings for the phase diagram for the sextet-only theory in the β–κ_6 plane.

For the four points marked by circles in Fig. <ref>, we ran zero-temperature simulations at the same bare couplings in order to determine the lattice scale (see Table <ref>).
As we describe in Appendix  <ref>, we set the lattice scale in each zero-temperature ensemble through calculation of the flow scale t_1/a^2.
Choosing the fiducial value 1/√(t_1)≡780 MeV gives a physical value to the lattice spacing a in each ensemble, and hence to the temperature T=(N_ta)^-1.
As can be seen in Fig. <ref>, one of the ensembles is a blue point on the confined side of the transition, while the other three are orange points on the deconfined side.
These provide lower and upper bounds, respectively, on T_c at the corresponding m_6 values.
We plot these temperatures in Fig. <ref>.
The transition temperature curve must pass below the upper bounds (downward arrows) and above the lower bound (upward arrow).
We can compare the transition temperatures seen here to those in more familiar theories: the transition in the pure SU(3) gauge theory occurs near 280 MeV <cit.>, while the crossover in QCD at physical quark masses occurs near 150 MeV <cit.>.



While the transition seen in Fig. <ref> may look continuous, we can make no claim concerning the order of the transition for any value of m_6 since we have only a single volume.  As the fermion mass m_6→∞ we expect to obtain SU(4) pure gauge theory. In this limit there is a first order transition, which should persist for large values of m_6.



  §.§.§ Effects of the NDS action


Our collaboration has previously examined the phase diagram of the sextet-only theory <cit.>, but without the NDS term in the gauge action.
In this previous study, the sextet-only theory was found to have a bulk transition.
The plaquette showed a large discontinuity at a κ_6 value below the thermal transition.
This large discontinuity is absent in our data. 
In Fig. <ref>, the plaquette shows structure occurring simultaneously with the response of the Polyakov loops, but no additional structure.

To search for a bulk transition, we ran a grid of 4^4 ensembles over the same region of bare 
parameter space covered by our 12^3×6 data. We found that the finite-temperature transition shifts substantially when changing N_t, thus confirming that it is not a softened bulk transition.
The NDS term appears to have completely banished the bulk transition, at least from the region of bare parameter space that we have explored.

    	
	

 §.§ Fundamental-Only Theory: κ_6=0

		
		


The other limiting case of our model contains only fundamental fermions and no sextets.
Figure <ref> shows the behavior of our  observables along a typical slice
 through bare parameter space, varying κ_4 over the transition at fixed β=9.2.
The top panel shows our confinement diagnostics, the fundamental Polyakov loop and the flowed anisotropy.
The two quantities change simultaneously and only once, varying smoothly otherwise: there is only a single crossover in each observable.
The middle panel shows our chiral diagnostics, the mass splittings of parity-partner mesons.
The splittings smoothly go to zero beyond κ_4 = 0.127, indicating chiral restoration. 
Comparing the top and middle panels, we see that the chiral and confinement crossovers overlap.
In the bottom panel, we see that the quark masses and plaquette vary smoothly, as expected for a crossover.



The behavior observed on the slice in Fig. <ref> is typical for this theory: everywhere we have investigated, we observe only a single unified chiral and confinement crossover.
Figure <ref> summarizes our findings for the β–κ_4 phase diagram for the fundamental-only theory.
As we did for the sextet theory, we have determined the physical temperature at three points along the transition, this time choosing three points inside the crossover region (points enclosed by circles in Fig. <ref>).
See Table <ref> for the required zero-temperature data.
We plot T_c versus the fermion mass m_4 in Fig. <ref>.
[Note that our rough estimate of T_c(m_4) does not correspond to the peak of any susceptibility in the crossover region.]

Here also, we expect that the first-order deconfinement transition of the pure gauge theory will reappear as the fermion mass rises towards infinity. Our simulations have not yet encountered this transition for masses as large as m_4≈ 0.32 / √(t_1) = 250 MeV.
 





§ PHASE STRUCTURE OF THE FULL THEORY

   	
   	

Figure <ref> depicts a slice through bare parameter space in the full theory with both fermion species dynamical, varying κ_6 while holding β=7.4 and κ_4=0.1285 fixed.
The top panel shows the behavior of our confinement diagnostics.
The Polyakov loops for both representations and the flowed anisotropy all jump simultaneously.
This means that the two species confine simultaneously.
The middle panel shows the behavior of the chiral diagnostics for both 
representations, the mass splittings of parity partner states.
The parity partners of both representations are split significantly at small κ_6, but simultaneously become nearly degenerate as κ_6 is increased.
This indicates that chiral symmetry restoration occurs simultaneously for the two representations.
Comparing the top and middle panels, we see that the combined confinement transition and the combined chiral transition coincide.
Away from the single jump, all quantities vary smoothly.
This behavior is typical throughout the region of bare parameter space that we have investigated: all phase diagnostics jump simultaneously and only once.
Thus, we find only two phases: a low-temperature phase where all fermions are confined and chirally broken and a high-temperature phase where all fermions are deconfined and chirally symmetric.






Figures <ref> and <ref> show phase diagrams for the theory at β=7.75 and β=7.4 for N_t=6 and N_t=8.
In these plots, confined and chirally broken regions of parameter space are highlighted in blue while deconfined and chirally restored regions of the parameter space are highlighted in orange.
The transition thus lies somewhere in the white band in each phase diagram.
Points enclosed by diamonds are where the phase diagnostics change when varying N_s from 12 to 18 when N_t=6 and from 16 to 24 when N_t=8.
The absence of many such points indicates that the location of the transition is insensitive to finite-volume effects.
Points enclosed by circles have zero-temperature data available at parameters near the transition (see Table <ref>).
Comparing the left and right panels in Figs. <ref> and <ref>, we see that the transition moves substantially in bare parameter space as we vary N_t.
This behavior is consistent with a thermal transition, and inconsistent with a bulk transition.

As shown in the bottom panel of Fig. <ref>, the plaquette and quark masses show a discontinuous jump at the transition, providing strong evidence that the observed transition is first order.
In support of this finding, we have also observed several tunneling events in the process of equilibrating new ensembles near the transition.
We have observed this to occur after more than 1000 trajectories, much longer than the typical equilibration time for this volume.

We are interested in determining whether the transition temperature T_c is comparable to its value in QCD, how strongly T_c depends on fermionic effects, and how significant are lattice spacing artifacts in T_c.
In the full theory, T_c is a function of the fermion masses and 1/N_t.
We do not have sufficient data to constrain the location of the transition with any of m_4, m_6, or a held fixed.
Instead, we examine the behavior of T_c as we interpolate along the transition at fixed β and 1/N_t. 
For simplicity, we use κ_4 to parameterize each transition curve, along which all of m_4, m_6, and T_c vary.
We have estimated the lattice spacing and thus T_c using the fit to t_1/a^2 described in Appendix <ref>.
Tracing along the transition bands in Figs. <ref> and <ref>, wherever there are ensembles matched in κ_4 on the edges of the transition band, we use our t_1/a^2 fit to estimate upper and lower bounds for T_c.
The resulting bounds on T_c as a function of κ_4 are the horizontal black dashes with error bands in Figs. <ref> and <ref>.
As indicated by the spanning arrows, the transition temperature must lie between these bounds.
Note that there is an uncontrolled systematic error for these bounds: the phases of some ensembles near the transition edge may be misdiagnosed if they have not been equilibrated long enough to tunnel to the correct phase.

As in the sextet-only and fundamental-only theories, T_c is comparable to its value in QCD (150 MeV) and SU(3) pure gauge theory (280 MeV).
For both N_t=6 curves, T_c may not remain constant as we vary κ_4 to interpolate along the transition.
We cannot exclude that the dependence may be a lattice artifact, but as one might expect from fermionic influence on the transition, T_c appears to depend more strongly on κ_4 at β=7.75 than at β=7.4.
Comparing with the κ_c curves of Fig. 19 of Ref. <cit.>, we see that at β=7.4 the transition curve is roughly parallel to κ_c and thus traces lines of approximately constant quark mass for whichever species is lighter; this slow variation in the masses is consistent with the observed slow variation in T_c. Meanwhile, at β=7.75 the transition curve moves further away from κ_c when κ_4 ≈κ_6, leading to heavier fermions; when the fermions are heavier, the system becomes more pure-gauge-like and T_c increases.








§ CONCLUSIONS

	
    	We summarize our results for the phase structure of our theory in two Columbia plots  in Fig. <ref>.
The axes in Fig. <ref> are plotted in lattice units, while the axes in Fig. <ref> are plotted in physical units, as defined in Appendix <ref>.
Each plot has two kinds of symbols.
Open symbols show quark masses from matching zero-temperature simulations: a m_4, a m_6, and t_1/a^2 are all taken from zero-temperature ensembles run at bare parameters matched to finite-temperature ensembles near the transition.
Closed symbols show the quark masses from finite-temperature simulations: am_4 and a m_6 are measured on finite-temperature ensembles along the deconfined side of the transition curve, and t_1/a^2 is obtained from the fit described in Appendix <ref>.
If there were no lattice artifacts, the open and closed symbols for each (N_t, β) set would coincide.
The solid lines in Fig. <ref>, which come from the interpolation, lie close to the open symbols, which were measured directly and did not require a fitting function.
The curves in Fig. <ref> where m_4∞ and m_6∞ indicate the masses at which we have examined the phase structure of the full theory in detail and found only a single first-order thermal transition.
The points where m_4=∞ indicate where we have examined the sextet-only theory and found only a single transition.
The points where m_6=∞ indicate where we have examined the fundamental-only theory and found only a single crossover or transition.







Our numerical investigation of the full multirep theory finds only a single, first-order thermal transition.
This non-observation of separated chiral phase transitions is in direct contradiction to predictions of the Most Attractive Channel hypothesis, according to which the sextet fermions should condense before the fundamentals as the temperature is lowered.
While our exploration of the three-dimensional parameter space is by no means exhaustive, and separated phase transitions might exist for some values of fermion mass, we have examined the theory at masses ranging from 50 to 400 MeV and ruled them out in this domain.
We similarly find that the fundamental-only and sextet-only theories appear to be QCD-like, with a combined chiral and confinement transition.

In the multirep theory we find a strongly first-order phase transition.
This is consistent with a one-loop Pisarski-Wilczek scaling analysis
appropriate to the limit where m_4 = m_6 = 0 <cit.>.

Increasing the quark masses in the single-species theories, we expect eventually to run into the first-order region of the pure-gauge transition.
Our exploration of these theories stops short of the quark masses at which this occurs.
Similarly, because we have only explored relatively light quark masses, we are unable to determine whether the first-order regions surrounding the double-chiral limit and the pure-gauge limit are connected.


We observe no bulk transitions in the full  theory or in either of its single-species limits.
We attribute this to our use of the NDS action.
Comparison to the previous study of the sextet-only theory <cit.> makes this clear.

The finite-temperature properties of this model may have important implications for cosmology
if it or something like it happens to be realized in nature.
First order transitions in the early Universe give rise to gravitational 
waves with distinctive properties <cit.>.
These signals may be accessible to near-future gravitational wave detectors such as LISA <cit.>.














Research was supported by U.S. Department of Energy Grant
Number  under grant DE-SC0010005 (Colorado)
and  by the Israel Science Foundation under grants 449/13 and 491/17.
Brookhaven National Laboratory is supported
by the U. S. Department of Energy under contract DE-SC0012704.
Some computations were performed on the University of Colorado cluster.
This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award number CNS-0821794) and the University of Colorado Boulder. The Janus supercomputer is a joint effort of the University of Colorado Boulder, the University of Colorado Denver and the National Center for Atmospheric Research.
Additional computations were done on facilities of the USQCD Collaboration at Fermilab,
which are funded by the Office of Science of the U. S. Department of Energy.
The computer code is based on the publicly available package of the
 MILC collaboration <cit.>. 








§ SCALE SETTING AND TRANSITION TEMPERATURE

	
	


 §.§ Setting the scale using Wilson flow


In order to estimate the transition temperature, we need the lattice spacing.
To this end, we generated zero-temperature ensembles at points in bare parameter space along the transition surface and measured the Wilson flow scale t_1 <cit.>.

We can define a general flow scale t_* via

    t^2 E(t)|_t=t_* = C

where E(t) is a discretization of the gluonic action density G_μν G^μν at flow time t, and C is a constant defining the scale.
In QCD, for t_0 one sets C = 0.3 and for t_1 one sets C = 2/3; the physical scale enters by fixing
√(t_0)≃ 0.142 fm <cit.>.
To translate quantities to physical values for comparison with QCD, we need some way of matching our definition of t_0 in N_c=4 with the definition in N_c=3.

Large-N_c scaling arguments provide a translation procedure <cit.>.
The observable t^2 E(t) can be used to define the renormalized coupling at the scale t <cit.>,

    g^2_wf(t)≡128 π^2/3 (N_c^2-1)t^2 E(t)  .

The usual large-N_c scaling argument holds the 't Hooft coupling λ_0 = g_0^2 N_c constant in N_c.
Thus, g^2 ∼ 1/N_c at leading order, regardless of renormalization scheme.
We can immediately read off from Eq. (<ref>) that, for g^2 ∼ 1/N_c to hold, t^2 E(t)∼ N_c at leading order.
Comparing Eq. (<ref>), we conclude that a reasonable choice is C∝ N_c,
and so for SU(4) we define the scales t_0 and t_1 via

    t_0^2 E(t_0)   = 0.4 
    t_1^2 E(t_1)   = 8/9.


Assuming that our N_c scaling has produced an equivalent quantity in SU(4), we may take the value t_0 = (0.142 fm)^2 from SU(3) <cit.>.
Our data set, however, samples the transition at relatively large lattice spacings a ≳√(t_0).
In this regime, t_0 develops nonlinear lattice-spacing artifacts; in practice, this is the failure of t^2 E(t) to reach the linear regime before it exceeds 0.4, and so we are sampling the “knee” in the typical t^2 E trajectory.
Because t_1 is larger, it can be measured on larger lattice spacings than t_0.
In order to use t_1 in lieu of t_0, we need a conversion factor.

In Fig. <ref> we examine the behavior of the ratio of √(t_1/t_0) as a function of the lattice spacing over our entire zero-temperature data set <cit.>.
The quantities √(t_1) and √(t_0) are fixed lengths, so their ratio should be some constant independent of lattice spacing.
We see, however, that at large lattice spacing the ratio is not constant.
Making the cut a < 1/√(2 t_1) (the dashed line in Fig. <ref>), we find √(t_1/t_0)≃ 1.77, and so √(t_1)≃ 0.252 fm or equivalently 1/√(t_1)≃ 783 MeV.
For simplicity we take 1/√(t_1)≡ 780 MeV for our fiducial value.
Determining the physical temperature of any given ensemble merely amounts to measuring the dimensionless quantity t_1/a^2 in order to calculate T=(N_ta)^-1=(N_t√(t_1))^-1×√(t_1/a^2).






 §.§ Fitting the lattice spacing


We have only a few zero-temperature ensembles in the region of bare parameters relevant to the transition.
In order to perform a more detailed temperature analysis, we modeled t_1/a^2 as a function of the bare parameters.
We do not have any theoretical expectations for the form of t_1/a^2, so this modeling is entirely empirical.

To include in our fit, we have 39 zero-temperature ensembles where t_1/a^2 > 1 to avoid discretization effects; of these, 12 ensembles are at β=7.75.
Lattice spacings computed using the Wilson flow are easily determined to very high statistical precision.
Without knowing the true form of t_1/a^2 as a function of β, κ_4, and κ_6, and with our limited dataset, we are unable to produce a model that can fit t_1/a^2 with convincingly small χ^2.
We are thus led to inflate the errors on t_1/a^2 to 2.5% of the value.
Such a model predicts the value of t_1/a^2 to within 2.5% for any (β,κ_4,κ_6) in the region of interest.

The form of our model is motivated by several observations about the behavior of t_1/a^2 as a function of the bare parameters (β, κ_4, κ_6).
At fixed β, we observe that (1) as κ_4 or κ_6 is increased, t_1/a^2 increases monotonically; (2) t_1/a^2 varies smoothly as a function of κ_4/κ_6, in such a way that there are smooth curves of constant t_1/a^2; (3) these curves of constant t_1/a^2 are roughly elliptical in shape; and (4) as both κ's go to zero and the fermions decouple, t_1/a^2 settles to a constant.
This motivates the functional form

    t_1/a^2 = exp[ r^2 - r_0^2(β)/γ(β)]
    	+ C(β) ,

where

    r^2 ≡ (8κ_4)^2 + α(β) (8κ_6)^2

and α(β), γ(β), r_0^2(β), and C(β) are functions of β only and thus constant at fixed β.
An equally reasonable functional form would be a power law rather than an exponential; however, fits to such functional forms produce unreasonably large powers and do not model the data as well.
Note that we write the model in terms of 8κ_4, 8κ_6, and β/8 (below) which are all O(1), so that the size of the fit parameters may be compared easily.
Physically, r_0^2 quantifies the value of r^2 where fermionic effects are frozen out; α projects the elliptical curves of constant t_1/a^2 in the κ_4–κ_6 plane to circles; γ quantifies how quickly t_1/a^2 increases as fermionic effects become strong; and C is the value of t_1/a^2 in the pure-gauge limit where κ_4=κ_6=0 [up to a small exp(-r_0^2/γ) correction].

To obtain a concrete realization of the abstract model (<ref>), we approximate α(β), γ(β), and r_0^2(β) as linear functions,

    α(β)    = α_0 + α_1 (β/8) 
    γ(β)    = γ_0 + γ_1 (β/8) 
    
    		r_0^2    = R_0 + R_1 (β/8)
 
where α_0, α_1, γ_0, γ_1, R_0, and R_1 are fit parameters.
In the pure-gauge limit where κ_4=κ_6=0, as β→ 0, we expect a →∞ and thus t_1/a^2 → 0.  We thus demand that C(0)=0, and model C(β) as a power,

    C(β) = C_0 (β/8)^C_1,

where C_0 and C_1 are fit parameters.
Fitting the 39 ensembles of our zero-temperature t_1/a^2 dataset to the model defined by Eqs. (<ref>)–(<ref>), we obtain the fit parameters in Table <ref> with χ^2/31=0.87 and Q=1-P=0.67.
The resulting model predicts the value of t_1/a^2  within 5% for all 39 ensembles included in the fit, and within 2.5% for 28 of these ensembles.
Figure <ref> shows the predictions of the model versus data at β=7.75.








To cross-check our model, we compare with fits of subsets of the dataset to simpler models.
At fixed β, all of α(β), γ(β), r_0^2(β), and C(β) are constant, providing a simplified four-parameter realization of Eq. (<ref>).
Fitting to 12 ensembles at fixed β=7.75 yields the model parameters in the second column of Table <ref> with χ^2/12=0.78 and Q=1-P=0.67.
By comparison, the model parameters in the third column of Table <ref> are predictions obtained by plugging the best-fit parameters in Table <ref> into Eqs. (<ref>) and (<ref>) for β=7.75.
Even though the full-dataset fit includes more than three times as many ensembles, the parameters agree closely; if this did not hold, it would suggest that the model is overfitting the data.
For the 12 ensembles at β=7.75, the predictions of these two fits agree within 1.5%.



	





§ SUMMARY OF ENSEMBLES NEAR THE TRANSITION

	
	



	



