Finite-temperature phase structure of SU(4) gauge theory with multiple fermion representations

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.

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 [1,2]. Our previous work [3,4] 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 [5,6] 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 [7][8][9][10], 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 [11], but it is possible that the theory explored in that work-the SU(3) theory with adjoint fermions-is infrared conformal rather than confining [12]. 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. 1 and 13) below.
The paper is organized as follows: In Section II, we discuss the theory, its symmetries, and the expected behavior of its confinement and chiral transitions. We discuss our lattice methodology in Section III. We then present our numerical findings. In Section IV, 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 V. We conclude in Section VI. A preliminary account of this work can be found in [13]. For lattice work on a composite Higgs model with a different gauge group, see [14].
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.

A. 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 where Ut( 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 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 [15]. 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.

B. 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. [5]), or fundamental first, or a simultaneous transition. As discussed below in Section V, our results indicate that the third possibility is what occurs for our system.
One of us [16] has carried out a Pisarski-Wilczek stability analysis [17] 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 singlerepresentation theories coupled together. The prediction of this calculation [16], carried out to first order in = 4 − d, is a first-order phase transition. Fig. 1 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 [18]. 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 [19,20] 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. [21]), so we expect the fundamentalonly theory to behave similarly. Fig. 1 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.

A. Action
For this study, we employ the same lattice action as in Ref. [3]. For the fermions, we use a clover-improved Wilson action built from fat gauge links constructed by normalizedhypercubic (nHYP) smearing [22,23]. We construct the action for the sextet fermions by promoting the smeared links to the sextet representation [15]. We set both clover coefficients equal to unity, c SW = 1, a choice motivated by results from Ref. [24]. The action for the gauge sector is the usual plaquette action with gauge coupling β augmented by an nHYPdislocation suppression (NDS) term [25], constructed from the nHYP-smeared links. As in our zero-temperature study, we fix the NDS parameter γ such that β/γ = 125. Altogether, Columbia plot illustrating expectations for the order of the finite-temperature phase transition. The axes are the masses of the two fermion species in the theory, with m 4 on the x-axis and m 6 on the y-axis. The upper right corner is the pure-gauge limit; the lower left corner is the double chiral limit; the upper-left corner is the fundamental-only chiral limit; the lower-right corner is the sextet-only chiral limit. Green fields indicate regions of parameter space where the theory is predicted to exhibit a first order transition. Blue lines indicate regions of parameter space where the theory is predicted to exhibit a second-order transition.
the simulation parameter space is three dimensional: β, and two hopping parameters κ 4 and κ 6 .

B. Spectroscopy
Fermion masses m 4 and m 6 for the two representations are defined through the axial Ward identity (AWI), where 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 O 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 [26][27][28]. These are called "P+A correlators" in the literature.

C. 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 I. To render this exploration tractable, we found it necessary to automate much of our data generation and analysis; see [29] 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 IV in the Appendix summarizes these zero-temperature data sets. We use the lattice scale to derive the physical temperature along the phase boundary.

D. 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 [30,31]: 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 [32][33][34]. Previous applications of anisotropy have used the quantity 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 (cf. Ref. [35]). 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 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 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. [36].

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
We begin with the gauge theory coupled only to sextet fermions. Figure 2 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 [15] 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. 2 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 3 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. 3, we ran zero-temperature simulations at the same bare couplings in order to determine the lattice scale (see Table IV). As we describe in Appendix A, 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 t a) −1 . As can be seen in Fig. 3, 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. 4. 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 [37], while the crossover in QCD at physical quark masses occurs near 150 MeV [21].
While the transition seen in Fig. 2 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 [15], 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. 2, 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. The other limiting case of our model contains only fundamental fermions and no sextets. Figure 5 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. 5 is typical for this theory: everywhere we have investigated, we observe only a single unified chiral and confinement crossover. Figure 6 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. 6). See Table IV for the required zero-temperature data. We plot T c versus the fermion mass m 4 in Fig. 7. 1 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.  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. To the left is the phase diagram for N t = 6 lattices, while to the right is the same region of bare parameter space for N t = 8 lattices. Blue dots indicate confined and chirally broken ensembles. Yellow stars indicate deconfined and chirally restored ensembles with m r > 0 for both species; red Xs are in regions where m 4 < 0 or m 6 < 0 or both. In the right figure, the transition region from N t = 6 is overlaid in gray, demonstrating that the transition moves as N t is varied. The circled ensembles have matching zero-temperature ensembles available. Ensembles enclosed by diamonds are where the phase changed when volume was changed (see text). Figures 9 and 10 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 IV). Comparing the left and right panels in Figs. 9 and 10, 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. 8, 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 A. Tracing along the transition bands in Figs. 9 and 10, 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. 11 and 12. 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. [3], 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. and N t = 8. Lattice spacings used to determine the temperature are computed using the fit to t 1 /a 2 with t 1 = 1/(780 MeV) 2 as discussed in Appendix A. Axes are matched between N t = 6 and N t = 8. Black lines with error bands indicate the temperature on ensembles on either end of the transition bands in Fig. 10. The transition temperature lies in the span indicated by the double-headed arrows.

VI. CONCLUSIONS
We summarize our results for the phase structure of our theory in two Columbia plots in Fig. 13. The axes in Fig. 13(a) are plotted in lattice units, while the axes in Fig. 13(b) are plotted in physical units, as defined in Appendix A. Each plot has two kinds of symbols. Open symbols show quark masses from matching zero-temperature simulations: am 4 , am 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 am 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 A. If there were no lattice artifacts, the open and closed symbols for each (N t , β) set would coincide. The solid lines in Fig. 13(b), 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. 13 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 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 [16].
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 [15] makes this clear.  Fig. 13(a), the quark masses are in lattice units. In Fig. 13(b), quark masses are in MeV [t 1 ≡ 1/(780 MeV) 2 ], as defined in Appendix A. Each color and symbol is associated with a different β and N t . Closed symbols are finite-temperature quark masses, while hollow symbols are zero-temperature quark masses from ensembles on the transition boundary The lattice spacings for the zero temperature quark masses are computed directly from t 1 /a 2 on that ensemble. The lattice spacings for the finite-temperature quark masses are computed from the model described in Appendix A.
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 [38,39]. These signals may be accessible to near-future gravitational wave detectors such as LISA [39].   . (A4) 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) [42]. Our data set, however, samples the transition at relatively large lattice spacings a √ t 0 . In this regime, t 0 develops nonlinear latticespacing 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. 14 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 [3]. 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/ √ 2t 1 (the dashed line in Fig. 14), 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 t a) −1

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 where and α(β), γ(β), r 2 0 (β), 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 2 0 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 2 0 /γ) correction]. To obtain a concrete realization of the abstract model (A5), we approximate α(β), γ(β), and r 2 0 (β) as linear functions, α(β) = α 0 + α 1 β 8 γ(β) = γ 0 + γ 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, 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. (A5)-(A8), we obtain the fit parameters in Table II 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 15 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 2 0 (β), and C(β) are constant, providing a simplified fourparameter realization of Eq. (A5). Fitting to 12 ensembles at fixed β = 7.75 yields the model parameters in the second column of Table III with χ 2 /12 = 0.78 and Q = 1 − P = 0.67. By comparison, the model parameters in the third column of Table III are predictions obtained by plugging the best-fit parameters in Table II into Eqs. (A7) and (A8) for β = 7.75. Even though the full-dataset fit includes more than three times as many ensembles, the parameters