Stealth dark matter confinement transition and gravitational waves

We use non-perturbative lattice calculations to investigate the finite-temperature confinement transition of stealth dark matter, focusing on the regime in which this early-universe transition is first order and would generate a stochastic background of gravitational waves. Stealth dark matter extends the standard model with a new strongly coupled SU(4) gauge sector with four massive fermions in the fundamental representation, producing a stable spin-0 'dark baryon' as a viable composite dark matter candidate. Future searches for stochastic gravitational waves will provide a new way to discover or constrain stealth dark matter, in addition to previously investigated direct-detection and collider experiments. As a first step to enabling this phenomenology, we determine how heavy the dark fermions need to be in order to produce a first-order stealth dark matter confinement transition.

We use non-perturbative lattice calculations to investigate the finite-temperature confinement transition of stealth dark matter, focusing on the regime in which this early-universe transition is first order and would generate a stochastic background of gravitational waves. Stealth dark matter extends the standard model with a new strongly coupled SU(4) gauge sector with four massive fermions in the fundamental representation, producing a stable spin-0 'dark baryon' as a viable composite dark matter candidate. Future searches for stochastic gravitational waves will provide a new way to discover or constrain stealth dark matter, in addition to previously investigated directdetection and collider experiments. As a first step to enabling this phenomenology, we determine how heavy the dark fermions need to be in order to produce a first-order stealth dark matter confinement transition.

I. INTRODUCTION AND OVERVIEW
The confining gauge-fermion theory of quantum chromodynamics (QCD) produces the massive stable protons and nuclei of the visible universe, making it compelling to hypothesize that new strong dynamics could also underlie the dark sector. Stealth dark matter [1,2] is a particularly attractive model of composite dark matter, based on a new strongly interacting SU(N D ) gauge sector with even N D ≥ 4, which is coupled to four massive fermions in the fundamental representation. As detailed in Ref. [1], the four 'dark fermions' transform in non-trivial vector-like representations of the electroweak group, in order to generate the correct cosmological dark matter abundance while also satisfying all experimental constraints. Although these 'dark' fermions are electrically charged and couple to the standard model (SM) Higgs boson, following the dark-sector confinement transition they give rise to a composite dark matter candidate in the form of the lightest spin-0 SU(N D ) 'dark baryon', which is a singlet under the entire SM gauge group. This dark matter candidate is automatically stable on cosmological time scales due to the conservation of dark baryon number, and it acquires mass both from confinement and from the masses of its fermion constituents. * david.schaich@liverpool.ac.uk Experimental constraints on the stealth dark matter model come from both direct-detection searches and collider experiments, with direct-detection cross sections arising from non-perturbative form factors of the dark baryon. For example, direct detection through Higgs boson exchange depends on the dark baryon's scalar form factor, as well as on the relative sizes of the vectorlike and electroweak-breaking fermion mass terms that appear in the model's lagrangian [1]. Existing directdetection searches, combined with lattice calculations of that scalar form factor, require that the vector-like contributions to the dark fermions' masses dominate over the electroweak-breaking contributions [1,3]. Those lattice calculations considered the minimal case N D = 4, which is also the case we will consider in this work. This choice minimizes the computational costs of our lattice calculations, while still being large enough for large-N scaling relations to recast results to larger N D ≥ 6 with reasonable reliability. (See Ref. [4] for a thorough review of the large-N framework.) Direct detection can also proceed through photon exchange, and the symmetries of the model strongly suppress this cross section by forbidding the leading magnetic moment and charge radius contributions to it. The contribution from the dark baryon's electromagnetic polarizability is unavoidable, and provides a lower bound on direct-detection signals for the entire class of dark matter models featuring neutral dark baryons with charged constituents (reviewed in Ref. [5]). Lattice calculations of that polarizability [2], again for the case N D = 4, obtain the constraint M DM 0.2 TeV from existing directdetection searches. 1 The steep dependence of the cross section on the dark baryon mass, σ ∝ 1/M 6 DM , causes the predicted signal to fall below the irreducible neutrino background for M DM 0.7 TeV [2].
Stronger constraints on stealth dark matter currently come from collider searches for vector (V ) and pseudoscalar (P ) 'dark mesons', some of which are electrically charged. If M P /M V < 0.5 so that V → P P decay is possible, the dark vector meson becomes a broad resonance and masses as light as M P 0.13 TeV and M V 0.3 TeV remain viable [7]. Lattice calculations of the meson and baryon spectrum can translate these bounds into constraints on M DM > M V . In this work we will focus on the heavy-mass regime, M P /M V > 0.5, where V → P P decays are kinematically forbidden. The dominant decay process is then V → + − , which could be observed in searches for Z → + − . This produces the constraint M V 2 TeV reported by Ref. [7], assuming this process is dominated by a single dark vector meson. In the heavy-mass regime, we can approximate M DM N D 2 M V to turn this into a lower bound on the dark baryon mass.
It is difficult to set an upper bound on the mass of the dark baryon, though some very rough estimates can be made by requiring that the stealth dark matter model produces the observed cosmological dark matter abundance. Specifically, Ref. [1] estimates that a predominantly symmetric thermal abundance of stealth dark matter would match cosmology for M DM of order tens to hundreds of TeV, while M DM smaller than a few TeV would require a predominantly asymmetric abundance. There is therefore a significant allowed range of stealth dark matter masses up to hundreds of TeV, which will be very challenging for direct detection or collider experiments to constrain.
This makes the possibility of using gravitational waves to constrain or discover stealth dark matter particularly exciting. There is increasing interest in probing dark sectors by searching for a stochastic background of gravitational waves that would be produced by a firstorder phase transition in the early universe [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. Such searches are an important component of the science programs for future space-based facilities including the LISA observatory [23,24], DECIGO [25] and AEDGE [26]. This approach has the advantage of involving only gravity, the force that provides the existing astrophysical and cosmological evidence for dark matter. In the context of strongly coupled composite models such as stealth dark matter, the transition of interest is the confinement transition through which the state of the system changes from a high-temperature deconfined plasma of 'dark gluons' and dark fermions to stable SM-singlet dark baryons. If this confinement transition was first order, its properties including the nucleation temperature and latent heat govern the stochastic spectrum of the gravitational waves it produced, making reliable knowledge of these properties a crucial ingredient to extract constraints from future observations [8,23,24].
In this paper we use non-perturbative lattice calculations to investigate the finite-temperature confinement transition of SU(4) stealth dark matter. We focus on the first goal of determining the region of parameter space for which the confinement transition of this gaugefermion theory is first order, in contrast to the continuous crossover of QCD. Achieving this first goal is a necessary step to enable more detailed future lattice investigations of the resulting gravitational waves. Some preliminary results from this work previously appeared in Ref. [27].
The 'Columbia plot' [28] shown in Fig. 1 illustrates what we can expect based on symmetry arguments and continuum effective models [29][30][31]. Although we specialize this version of the plot to the two pairs of degenerate fermions that stealth dark matter involves, for generic SU(N ) gauge theories with N ≥ 3 and N f 2N fundamental fermions, first-order transitions are expected in two regimes: where the fermions are sufficiently heavy or sufficiently light. These expectations have been supported by lattice calculations, though at present controlled continuum extrapolations have been achieved by lattice analyses of only two points on the Columbia plot. One of these is the (2+1)-flavor physical point of SU (3) QCD-see the recent review Ref. [32] and references therein. The other is the SU(N ) pure-gauge system that corresponds to the infinite-mass limit in the upper-right corner of the plot, for 3 ≤ N ≤ 10 [4,[33][34][35][36].
In this work we will focus on the SU(4) heavy-mass first-order transition region connected to the pure-gauge limit. Compared to the light-mass region, this both reduces computational costs and limits the reach of collider constraints on M DM , which become more powerful as the ratio M DM /M P grows towards the M P → 0 chiral limit. Although stealth dark matter requires at least a small mass splitting between the two pairs of degenerate fermions, in order to guarantee that all 'dark mesons' are unstable and do not disrupt Big Bang nucleosynthesis [1], for simplicity we will consider in this work only the limit of four degenerate flavors, corresponding to the "N f = 4" diagonal line in Fig. 1. In Section V we will discuss prospects for future investigations of the more general non-degenerate situation.
The first goal mentioned above is now a matter of determining how heavy these N f = 4 dark fermions need to be in order to produce a first-order stealth dark matter confinement transition. This investigation is the first lattice study of the heavy-mass region of the Columbia plot for SU(4) gauge theory with dynamical fermions. Even for the case of SU(3) this region has received relatively little attention compared to the QCD physical point and chiral limit. See Ref. [37] (and references therein) for a recent SU(3) investigation, which finds that very large masses are needed to produce a first-order transition. These masses are parameterized by the ratio M P /T c 10, where T c is the equilibrium critical temperature and the extrapolation to the continuum limit is not yet under control. The need for very large masses for a first-order SU(3) transition may be related to the known weakness of the first-order confinement transition in the SU(3) pure-gauge limit [4,37]. Since this pure-gauge confinement transition strengthens significantly with increasing N ≥ 4 [4,34,35], stealth dark matter may exhibit qualitatively different behavior, motivating our dedicated lattice calculations.
We begin in the next section by explaining the strategy of our lattice calculations, including our nHYP-improved unrooted-staggered lattice action, the SU(4) ensembles we have generated using it, and the observables we focus on to analyze the confinement transition. In Section III we test our methods by considering the pure-gauge limit, which provides a less-expensive means to assess the discretization artifacts of our lattice action. We also exploit our prior knowledge that the pure-gauge SU(4) transition is strongly first-order, which allows us to view our pure-gauge results as a guide to the signals we should expect for a first-order confinement transition with dynamical fermions. In Section IV we add those N f = 4 degenerate dynamical fermions, and supplement our finitetemperature analyses with zero-temperature meson spectroscopy calculations. These ingredients allow us to determine the ratio of dark pion and dark vector meson masses, M P /M V > 0.9, required for the stealth dark matter confinement transition to be first order.
We discuss our conclusions in Section V, and look ahead to our follow-up work that will investigate this first-order transition in more detail, in order to predict more detailed features of the gravitational waves it would produce. Key parameters that need to be computed or estimated to predict the gravitational-wave spectrum are the latent heat (or vacuum energy fraction), the phase transition duration, and the bubble wall velocity [38,39]. Only the first of these is straightforward to determine through lattice calculations, and this will be the next focus of our investigations. Even without a careful continuum-extrapolated analysis of the latent heat, our results reported in this paper will allow future searches for stochastic gravitational waves (resulting in either detections or exclusions) to set novel constraints on stealth dark matter and similar models. For example, the gravitational-wave spectrum also depends on the transition temperature T * , which may differ from the equilibrium critical temperature T c used to set the scale of our lattice calculations, due to possible supercooling. If we can assume T * T c or estimate how they differ, then our results for the mass dependence of the stealth dark matter transition will translate information on T * from gravitational-wave searches into predictions for both the approximate mass scale of the dark baryons as well as the minimum masses of the dark mesons being searched for at colliders.

A. Context and lattice action
As usual [40], our SU(4) lattice calculations employ a hypercubic grid of L 3 × N t sites defining a discrete euclidean space-time. We impose thermal boundary conditions (periodic for bosons, antiperiodic for fermions) in the temporal direction, while all fields are subject to periodic boundary conditions in the spatial directions. The lattice spacing 'a' between neighboring lattice sites is set through the input bare gauge coupling β F ∝ 1/g 2 0 , which we discuss in more detail below. The temperature in 'lattice units' is the inverse temporal extent of the lattice, T = 1/(a·N t ), and in the finite-temperature context we are interested in N t < L.
For a fixed lattice volume L 3 ×N t we proceed by varying the bare coupling β F to scan the temperature. Below we discuss the observables we monitor as functions of the coupling, which reveal the critical β (c) F corresponding to T c , and provide information about the order of the transition. Setting the lattice scale by taking T c = 1/(a c · N t ) to be a fixed physical temperature means that the lattice spacing at the transition decreases as N t increases, identifying the a → 0 continuum limit with the limit N t → ∞. If N t is too small, the large lattice spacing may result in significant systematic errors from discretization artifacts.
At the same time, the aspect ratio L/N t must be suf-ficiently large to ensure that systematic errors from the finite spatial volume are also under control. This motivates keeping N t as small as discretization artifacts allow. The large lattice spacings at small-N t thermal transitions correspond to strong bare gauge couplings, and studies spanning many years [41,42] have observed that such strong couplings for can produce a bulk (zerotemperature) transition into a lattice phase with no continuum limit. For SU(4) Yang-Mills theory with a lattice action that includes both fundamental and adjoint plaquette terms, with respective couplings β F and β A , this bulk transition is first order for sufficiently large β A > 0, with a cross-over persisting when β A = 0. 2 With N t 4, thermal transitions for β A = 0 effectively merge with this bulk crossover, resulting in unmanageable discretization artifacts.
In an attempt to ameliorate this problem, we follow Ref. [43] and use a negative adjoint coupling β A = −β F /4 in the fundamental-adjoint gauge action. At tree level for SU(N ) gauge theory, requiring β A > −β F /2. This tree-level relation is not accurate at the critical β (c) F of the thermal transitions with N t ≤ 12, which can be seen by contrasting our pure-gauge results in the next section against past studies of SU(4) lattice gauge theory using β A = 0 [34-36, 44, 45]. Continuum extrapolations would therefore be required to quantitatively compare our puregauge results (e.g., for the latent heat) with that earlier work. The same is true for comparisons with Ref. [33], which avoids strong-coupling bulk transitions by modifying the lattice action to restrict the fundamental plaquette to a single Z 4 vacuum.
Unlike those prior pure-gauge studies, we also carry out calculations with four dynamical fermions in the fundamental representation of SU(4). As discussed in Section I, for simplicity we consider only four degenerate flavors, which allows us to use an unrooted staggeredfermion lattice action. To reduce discretization artifacts for the relatively large fermion masses a·m that we will consider, we also improve the fermion action by incorporating smearing. Again following Ref. [43], we use a single nHYP smearing step [46,47] with parameters (0.5, 0.5, 0.4).

B. Strategy
The considerations above lead us to the following strategy for the ensembles of gauge configurations we generate to map out the finite-temperature SU(4) phase diagram.
• We need to consider several fermion masses a·m in order to determine the regime in which the stealth dark matter confinement transition is first order. 3 Our smallest fermion mass a·m = 0.05 is chosen to overlap the mass range considered in Refs. [1,2]. We also carry out pure-gauge calculations corresponding to the a·m → ∞ quenched limit. In total we consider a·m = {0.05, 0.1, 0.2, 0.4, ∞}.
• For each of those five a·m, we want at least three N t in order to enable N t → ∞ continuum extrapolations. In total we consider N t = {4, 6, 8, 12}, but we will see in the next section that N t = 4 may suffer from large discretization artifacts despite our improved lattice action. We will therefore use N t ≥ 6 to carry out continuum extrapolations, which remain work in progress. While these continuum extrapolations will be important for our subsequent studies of (e.g.) the latent heat, they are not crucial for our present task of determining the SU(4) phase diagram. In this work we will focus on N t = 8, the largest temporal extent for which a large amount of data is available, using the other N t primarily to assess discretization artifacts.
• For each {a·m, N t } we want at least three aspect ratios L/N t ≥ 2 in order to enable extrapolations to the thermodynamic limit of infinite spatial volume. In our present work these multiple spatial volumes are most useful for distinguishing between firstorder transitions and continuous crossovers, for instance from the L dependence of relevant susceptibilities or kurtoses. More careful infinite-volume extrapolations will again feature in our upcoming detailed studies of transition properties. So far we have considered aspect ratios L/N t = {2, 3, 4, 6, 8}.
• Finally, for each {a·m, N t , L/N t }, we scan in temperature by varying the input bare fundamental coupling β F . We begin at a high value of β F deep in the deconfined phase and systematically lower the temperature through the transition and into the confined phase, starting each lower-temperature calculation from a thermalized gauge-field configuration generated at slightly higher β F . Once we are deep in the confined phase we reverse this process and also scan from low to high temperature in order to check for possible hysteresis. Following these initial coarse scans with relatively large ∆β F = 1-2 between subsequent calculations, we carry out one or two rounds of refined scans around the transition region with smaller 0.02 ≤ ∆β F ≤ 0.2.
For both pure-gauge and dynamical calculations we use the hybrid Monte Carlo (HMC) algorithm [48], employing QHMC/FUEL [49] on top of the USQCD SciDAC software stack, 4 which provides efficient performance for arbitrary SU(N ) gauge groups. We use a second-order Omelyan integrator [50] with multiple time scales [51] and (for a·m < ∞) an additional heavy pseudofermion field [52], fixing a trajectory length of τ traj = 1 molecular dynamics time unit (MDTU) and tuning molecular dynamics step sizes to target roughly 60%-80% acceptance rates [50]. We monitor the 'Creutz equality' [53] e −∆H = 1 to ensure that our HMC parameter choices are appropriate. We also accumulate a similar number of MDTU for both pure-gauge and dynamical calculations. While larger volumes and higher statistics could be obtained with more efficient algorithms in the pure-gauge case, our goal here is to use this known first-order transition to illuminate the signal quality we may expect from the algorithms and statistics available to us in the more expensive dynamical case.
In total, with a · m = {0.05, 0.1, 0.2, 0.4, ∞}, N t = {4, 6, 8, 12} and L/N t = {2, 3, 4, 6, 8} we have generated 1,369 finite-temperature HMC Markov chains (or 'streams'), each with at least 2,000 MDTU and up to 50,000 MDTU. We use the same HMC parameters for both high-and low-start streams, which allows us to combine 1,154 of these streams into 577 joint ensembles with approximately doubled statistics. In addition, we generated 12 zero-temperature ensembles with lattice volume 24 3 × 48, at the critical coupling β (c) F and at β (c) F ± 0.2 for each a · m < ∞. We use these zero-temperature ensembles to compute the meson spectrum and relate a · m to the ratio of dark pion and dark vector meson masses, M P /M V . This provides a convenient parameterization of the fermion masses that can easily be compared to previous quenched lattice studies of stealth dark matter [1,2], which used valence Wilson fermions with 0.55 M P /M V 0.77. The variation in the number of MDTU per finitetemperature stream is driven by auto-correlations that increase significantly around the transition (even if the 'transition' is a continuous crossover), requiring longer HMC streams in this region. For each stream we set a thermalization cut by hand based on human inspection of time-series plots, and use the 'autocorr' module in emcee [54] to estimate auto-correlation times τ for selected non-topological observables discussed below. We then divide our measurements into bins for jackknife analyses, with bin sizes larger than τ and at least 100 MDTU, collecting sufficient data to ensure that at least ten such statistically independent bins are available. The maximum auto-correlation time we observe, 4 usqcd-software.github.io τ ≈ 5400 MDTU, produces 16 jackknife bins, eight from each of the high-and low-start streams.

C. Observables
The key observable signalling the confinement transition is the Polyakov loop (P L), the gauge-invariant trace of the product of gauge links wrapping around the temporal extent of the lattice. In the pure-gauge SU(N ) theory, the Polyakov loop is an order parameter of the (temporal) Z N center symmetry, which breaks spontaneously in the high-temperature deconfined phase where the magnitude |P L| → N as β F → ∞ and the argument is restricted to lie near any one of the N degenerate vacua oriented at e iφ = e 2πik/N with k = 0, · · · , N −1. Dynamical fermions in the fundamental representation explicitly break this center symmetry, picking out the positive real axis (φ = 0) as the preferred vacuum. In order to apply identical analyses to both the pure-gauge and dynamical theories, we focus on the magnitude |P L| as the most useful observable.
We improve the signal for the Polyakov loop by computing it after smoothing the lattice gauge fields by applying the Wilson flow, a continuous transformation that systematically removes short-distance lattice cutoff effects [55,56]. This Wilson-flowed Polyakov loop P L W is a modern variant of the RG-blocked Polyakov loop investigated in older works [57,58], and has previously been used in Refs. [59][60][61][62][63]. The removal of short-distance fluctuations significantly enhances the signal without affecting the physics of the transition, producing much clearer contrasts between confined systems with small |P L W | 1 and deconfined systems with large |P L W | ∼ N . We restrict the 'flow time' t by requiring c ≡ √ 8t/N t ≤ 0.5 or equivalently t ≤ N 2 t /32. Since 4 ≤ N t ≤ 12, this maximal c = 0.5 still corresponds to modest flow times 0.5 ≤ t ≤ 4.5, respectively. In this paper we will therefore only show results obtained with c = 0.5. Behind the scenes we also monitor c = 0.2, 0.3 and 0.4 to check that our focus on c = 0.5 doesn't introduce systematic errors. In particular, |P L W | with c = 0.5 is the main observable whose auto-correlation time we monitor to set jackknife bin sizes. 5 In addition to the expectation value |P L W | itself, we also compute the susceptibility 5 We also monitor the auto-correlation time of the chiral condensate ψψ , but the relatively large masses we consider strongly break chiral symmetry and leave ψψ of little use for analyzing the confinement transition. and kurtosis (equivalent to the Binder cumulant) for the (volume-averaged) Wilson-flowed Polyakov loop magnitude O = |P L W |. This susceptibility exhibits a peak at the confinement transition, with the order of the transition reflected by the L-dependence of the peak height and of the kurtosis [64]. We will similarly use the plaquette susceptibility χ = L 3 N t 2 − 2 to identify the zero-temperature bulk phase transition. Because the plaquette is much less noisy than the Polyakov loop, there is no need to improve its signal with the Wilson flow.
Another quantity sensitive to the confinement transition is the spatial/temporal anisotropy of the Wilsonflowed energy density t 2 E(t) [60][61][62]. Following Ref. [62] we consider the ratio introduced by Ref. [65] for tuning anisotropic lattice spacings. Here the 'space-space' E ss (t) is computed from 'clover' terms built out of four plaquettes oriented in the purely spatial planes x-y, x-z and y-z, while the clover terms contributing to the 'space-time' E sτ (t) are oriented in the x-τ , y-τ and z-τ planes. We will again focus on values of the flow time t corresponding to c = 0.5. In the low-temperature confined phase, the system is isotropic and R E ≈ 1, while the breaking of temporal (but not spatial) center symmetry in the hightemperature deconfined phase produces R E > 1.
Finally, we also monitor the 'deconfinement fraction' discussed in Refs. [44,66], which measures the proportion of Polyakov loop measurements whose arguments fall within a certain (tunable) angle θ < π/4 around any of the Z 4 vacua. As above, we consider the Wilson-flowed arg(P L W ) at flow times corresponding to c = 0.5. With N in of N tot measurements suitably aligned along the Z 4 axes, we define the deconfinement fraction Using these observables, we will now reproduce the well-studied first-order confinement transition in puregauge SU(4) Yang-Mills theory, and use that experience to investigate the mass dependence of the stealth dark matter confinement transition with N f = 4 degenerate dynamical fermions.

III. PURE-GAUGE LIMIT
Over the years there have been several lattice investigations of the SU(N ) Yang-Mills confinement transition with N > 3, primarily exploring the approach to the large-N limit. See Refs. [33-36, 44, 45] for work with a focus on N = 4 (building on much earlier studies [67][68][69][70]) and Ref. [4] for a broader review. We revisit this calculation with two main goals, in addition to confirming that our code and algorithms are working correctly. First, we will use the computationally inexpensive puregauge limit to check the discretization artifacts of our improved fundamental-adjoint gauge action, and assess which N t will be safe to use in dynamical calculations without complications from the bulk transition discussed above. Second, our prior knowledge that the pure-gauge SU(4) transition is strongly first-order allows us to observe the quality of signals we should expect for a firstorder transition with dynamical fermions, which will be useful to distinguish this case from a continuous crossover in Section IV.

A. Discretization artifacts
In Fig. 2 we show how the critical coupling β (c) F of the pure-gauge thermal confinement transition depends on the temporal extent of the lattice N t , to clarify the more abstract discussions in Section II above. With fixed aspect ratio L/N t = 2 for N t = 4, 6, 8 and 12, the transition is clear in both the Wilson-flowed Polyakov loop magnitude |P L W | and the Wilson-flowed E(t) anisotropy, illustrating the behaviour described in Section II. As N t decreases, the fixed critical temperature T c = 1/(a c ·N t ) implies a larger lattice spacing, which in turn corresponds to the stronger bare coupling (smaller β (c) F ) shown in Fig. 2. This larger lattice spacing results in larger discretization artifacts, which only become unmanageable if the coupling becomes sufficiently strong to cause a zerotemperature bulk transition into a lattice phase with no continuum limit. This zero-temperature transition occurs around the same β F ≈ 13 for all N t , and is signalled by a peak in the plaquette susceptibility χ , as opposed to the peak in the Wilson-flowed Polyakov loop susceptibility χ |P L W | that is one signal of the confinement transition. In Fig. 3 we compare these two susceptibilities on the same set of axes for for lattice volumes 16 3 × 4, 24 3 × 6 and 32 3 × 8, each with aspect ratio L/N t = 4. Because the height of the peak in χ |P L W | is orders of magnitude larger than that in χ , we plot the relative  F on the temporal extent of the lattice Nt, comparing lattice volumes 8 3 ×4, 12 3 ×6, 16 3 ×8 and 24 3 ×12 with aspect ratio α ≡ L/Nt = 2. As Nt decreases, confinement occurs at stronger couplings (smaller β (c) F ), as shown by both the Wilson-flowed Polyakov loop magnitude |P LW | (top) and the Wilson-flowed E(t) anisotropy (bottom). We plot separate results for the high-and low-start streams, with lines connecting points to guide the eye, to show the absence of hysteresis. The small differences between the high-and low-start Nt = 12 results for 16.2 ≤ βF ≤ 17 are due to these streams freezing in different topological sectors, with topological charges Q = 0 and Q = −1, respectively. susceptibilities obtained by normalizing each data set by the maximum height of its respective peak.
In Fig. 3 we can see that the N t = 4 confinement transition at β (c) F ≈ 13.6 is dangerously close to the bulk transition at β F ≈ 13.2. We will therefore need to be wary of including N t = 4 in N t → ∞ continuum extrapolations, which was also the case for older studies using β A = 0 [35,44]. So although we can expect reduced discretization artifacts thanks to our improved fundamental-adjoint gauge action with negative β A = −β F /4, this improvement appears insufficient to allow us to rely on smaller, cheaper lattice volumes. Al- ready for N t = 6 we can see a much healthier separation between the two transitions in Fig. 3, which improves as N t increases thanks to the N t -dependence of the thermal confinement transition in contrast to the N t -independence of the bulk transition. For our ongoing studies of the latent heat and other properties of the stealth dark matter confinement transition, we therefore plan to carry out continuum extrapolations using N t = 6, 8 and 12. These continuum extrapolations are not crucial for our present task of determining the dynamical SU(4) phase diagram, so for the remainder of this work we will focus on N t = 8 as the largest temporal extent for which we have already accumulated a great deal of numerical data.

B. Order of the transition
The final goal of our small-scale pure-gauge calculations is to confirm our prior knowledge that the SU(4) confinement transition seen above is indeed strongly first order rather than continuous. We do this employing the same HMC algorithm, lattice volumes and statistics that we will use in the dynamical case, in order to illuminate the quality of signals we may expect to see for a firstorder transition with heavy dynamical fermions.
Already in Fig. 2 we can see that the Wilson-flowed Polyakov loop magnitude and the Wilson-flowed E(t) anisotropy do not show any sign of hysteresis for aspect ratio L/N t = 2. This remains true for larger aspect ratios as well. While hysteresis in the thermodynamic limit L → ∞ can be expected for a strongly first-order transi- tion, its absence for these lattice volumes does not imply a continuous transition in the infinite-volume continuum theory of interest. Indeed, from other observables we do have evidence confirming the known first-order nature of the pure-gauge SU(4) confinement transition. In particular, Fig. 4 shows the histogram of Wilson-flowed Polyakov loop magnitude |P L W | measurements on 24 3 ×8 lattices at β F = 15.0 near the confinement transition. The histogram features two clearly separated peaks, with approximately the same height, which is characteristic of the confined/deconfined phase coexistence at a first-order transition. This doublepeaked structure is clear confirmation that our calculations suffice to reproduce the known first-order SU(4) confinement transition.
A familiar means of determining the order of a con- FIG. 6. The deconfinement fraction f from Eq. (5) with uncertainties obtained as described in the text, for puregauge SU(4) lattice ensembles with Nt = 8 and aspect ratios L/Nt = 2, 3 and 4. The more rapid change from the f → 1 deconfined limit to the f → 0 confined limit with increasing L is consistent with a discontinuous first-order transition in the L → ∞ thermodynamic limit.
finement transition is to investigate how the maximum height χ max of the (Wilson-flowed) Polyakov loop susceptibility peak scales with the spatial lattice volume L 3 . A first-order transition is characterized by direct volume scaling χ max ∝ L 3 , in contrast to both the critical scaling χ max ∝ L 3b of a second-order transition with critical exponent b = 1 and the L-independence of a continuous crossover [71][72][73]. In Fig. 5 we present the |P L W | susceptibility peaks for our pure-gauge N t = 8 ensembles with aspect ratios L/N t = 2, 3 and 4, which are consistent with the expected first-order volume scaling. However, with the lattice volumes and statistics available to us it is difficult to quantitatively verify the volume scaling that would confirm a first-order transition. In addition to the large uncertainties around the transition, 6 the peak will occur at slightly different critical β (c) F for each different L, and the values of β F we have sampled may not exactly match these critical couplings. The situation is similar for the Wilson-flowed Polyakov loop kurtosis [Eq. (3)], which suffers from even larger uncertainties. Robustly determining these peak locations and heights is usually done through multi-ensemble reweighting [64,75], which we have not yet attempted. Instead, we will rely on our other evidence for a first-order transition, and take Fig. 5 as an indication of the behavior we should expect to see for a first-order confinement transition in stealth dark matter with dynamical fermions.
To the same end, in Fig. 6 we show the L dependence of the deconfinement fraction f for the same N t = 8 ensembles with aspect ratios L/N t = 2, 3 and 4. In Eq. (5) we normalized the deconfinement fraction so that f → 1 in the deconfined phase and f → 0 in the confined phase. These limits are clearly seen in Fig. 6, up to some residual fluctuations around zero in the β F < β (c) F confined regime. The key feature consistent with the first-order nature of the pure-gauge SU(4) confinement transition is that the change between these two limits becomes more rapid as L increases, eventually becoming discontinuous in the L → ∞ thermodynamic limit. This is another feature of a first-order transition that we will monitor in the case of the stealth dark matter confinement transition, to which we now turn.

IV. DYNAMICAL N f = 4 MASS DEPENDENCE
We now consider the more challenging task of studying stealth dark matter by coupling SU(4) lattice gauge theory to N f = 4 degenerate dynamical fermions. Compared to pure-gauge SU(N ) theories, much less work has been done to investigate finite-temperature dynamics with N > 3 and dynamical fermions. Ref. [62] investigates N f = 2 for 3 ≤ N ≤ 5 to explore the approach to the large-N limit, while Ref. [76] also considers N f = 2 for SU(4), as a limit of a theory with multiple fermion representations motivated by a composite Higgs model with partial compositeness.
Compared to composite Higgs studies in which some of the fermions must be massless and others are generically light in order to produce near-conformal dynamics, our task is simplified by considering relatively heavy fermions corresponding to the upper-right corner of the 'Columbia plot' in Fig. 1. As described in Section II, we consider a·m = {0.05, 0.1, 0.2, 0.4}, with the smallest a·m = 0.05 chosen to overlap with the masses considered by previous lattice studies of stealth dark matter [1,2]. The largest a·m = 0.4 turns out to be the only one for which we observe a first-order confinement transition. After presenting our results for the mass dependence of the transition, we will convert these values of a·m into ratio of dark pion and dark vector meson masses, M P /M V , for more direct comparison with Refs. [1,2].

A. Nt = 8 transition results
As for the pure-gauge limit in Section III, we begin by briefly considering the critical coupling β (c) F of the thermal confinement transition of stealth dark matter. Since the dependence on the temporal extent of the lattice N t is similar in both cases, in Fig. 7 we focus on the bare fermion mass a·m dependence of β (c) F for N t = 8, including the a·m → ∞ pure-gauge limit. As expected, lighter dynamical fermions more effectively screen the gauge interactions, requiring stronger bare couplings (smaller β F ) to produce the transition. Figure 7 shows this for both the Wilson-flowed Polyakov loop magnitude |P L W | and the Wilson-flowed E(t) anisotropy. From both these results and the corresponding |P L W | susceptibility peaks discussed below we can easily read off β  [77], it indicates that the fermions are not so heavy as to be effectively quenched.
Since we observed no hysteresis for these quantities in the pure-gauge case in Fig. 2, it is not surprising that none of our dynamical N f = 4 streams exhibit any hysteresis, either. For this reason we have simplified Fig. 7 by including only high-start results. An initial sign of a first-order transition for a · m = 0.4 comes from Fig. 8, which shows a double-peaked structure consistent with confined/deconfined phase coexistence at a first-order transition. Compared to the pure-gauge histogram in Fig. 4, the valley between the two peaks is much less dramatic in this dynamical case, and we see no two-peak structure for any of our a·m ≤ 0.2 ensembles. This suggests that a · m = 0.2 is sufficiently small to move the system out of the heavy-mass first-order region that appears to contain a·m = 0.4.
In Figs. 9 and 10 we more comprehensively compare our four dynamical masses a · m = {0.05, 0.1, 0.2, 0.4}, considering the same L dependence of the Wilson-flowed Polyakov loop susceptibility χ |P L W | and deconfinement fraction f as shown for the pure-gauge theory in Figs. 5 and 6, respectively. We again focus on N t = 8 with aspect ratios L/N t = 2, 3 and 4, generating a higher density of ensembles around the transition for each case, except a·m = 0.05 which is clearly a smooth crossover.
For the susceptibility χ |P L W | in Fig. 9, the height of the peaks increases by an order of magnitude as the mass increases from a·m = 0.05 to 0.4, though that last case still remains significantly below the scale of the pure-gauge peaks in Fig. 5 (again indicating that the fermions are not so heavy as to be effectively quenched). In combination with the fixed width of the horizontal axes, the increasing range of the vertical axes produces narrower-looking peaks as a·m increases. The key feature is the L dependence of the maximum peak heights, which as discussed in Section III is difficult to determine given the increasing uncertainties around the transition and the non-zero ∆β F = 0.02 separating ensembles in the transition region. Within these limitations, only the a·m = 0.4 results could be consistent with the volume scaling χ max ∝ L 3 of a first-order transition, with behavior qualitatively similar to that of the pure-gauge theory in Fig. 5. The smaller masses a·m ≤ 0.2 all produce susceptibility peaks without clear L dependence, which is the behavior expected for a continuous crossover. Figure 10 supports the same conclusion, most notably in the very slow decrease of the 32 3 × 8 deconfinement fraction on the β F < β (c) F confined side of the critical coupling indicated by the |P L W | susceptibility peaks for a · m ≤ 0.2. Empirically, we also observe f ≈ 1 for the 32 3 ×8 ensembles that produce the largest susceptibilities for each a·m ≤ 0.2, with only a·m = 0.4 and the puregauge theory producing (respectively) f = 0.906(41) and f = 0.795 (14) significantly different from unity. Again, the a·m = 0.4 results are the only ones qualitatively consistent with the pure-gauge behavior in Fig. 6. While the development of a discontinuity in the L → ∞ thermodynamic limit is not obvious in this case, the clear contrast with the a · m ≤ 0.2 results still suggests a change to a first-order transition for a·m = 0.4.

B. Zero-temperature spectroscopy
Our final task in this work is to parameterize the a · m discussed above in a convenient form for comparison with previous lattice studies of stealth dark matter [1,2]. We do this by computing the ratio of dark pion and dark vector meson masses, M P /M V , which requires 'zero-temperature' lattice calculations with N t > L. We carry out these zero-temperature calculations at the β (c) F of the N t = 8 transitions discussed above, for each bare fermion mass a·m = {0.05, 0.1, 0.2, 0.4}. As N t increases and the corresponding lattice spacing a c 1/(T c ·N t ) decreases, we will need to consider correspondingly smaller bare masses a · m in order to take the a → 0 continuum limit along a 'line of constant physics' with fixed M P /M V . In this work we restrict ourselves to determining the M P /M V corresponding to the N t = 8 transitions.
To determine M P and M V we carry out correlated fits of the corresponding two-point staggered correlation functions, over appropriate fit ranges [t min , t max ]. We do not include any excited states in our fits, instead considering relatively large t min to reduce any possible excitedstate contamination. For a · m = 0.05 and 0.1, we fix t max = N t /2 = 24 and combine results for all t min in the range 10 < t min < 18. For the larger masses a·m = 0.2 and 0.4, the exponential decay of the correlation functions C(t) ∼ e −M t at large times t can cause the signal in the vector channel to be overwhelmed by statistical noise for t < N t /2. This requires that we set a smaller t max = 16, which in turn demands a smaller range of 6 < t min < 12.
Our results for a·M P and a·M V are compiled in Table I, where for reference we also include results for the scale √ 8t 0 introduced in Ref. [78] and defined through the Wilson flow discussed in Section II C. Following   Refs. [76,79,80], we define this scale through the condition t 2 E(t) t=t0 = 0.4, where the energy density E(t) is evaluated after flow time t using the standard clover construction mentioned in Section II C. This choice incorporates the leading-order scaling t 2 E(t) ∼ N to generalize the canonical SU(3) value of 0.3 to our SU(4) theory. For convenience we also record the ratio of the pseudoscalar meson mass to the N t = 8 critical temperature, M P /T c = a·M P N t . The results shown in Table I do not include systematic uncertainties related to the choice of fit ranges and possible excited-state contamination or finite-volume effects. Based on our expectation that the overall uncertainty in the M P /M V ratio of interest will be dominated by its dependence on the coupling β F , we simply set that overall uncertainty by varying β We can also see that larger β F (smaller lattice spacings) produce larger M P /M V , confirming that smaller a · m will be needed to stay on a line of constant physics when taking the N t → ∞ continuum limit in future work.
Previous lattice studies of stealth dark matter [1,2] considered the mass range 0.55 M P /M V 0.77, using valence Wilson fermions on quenched gauge field configurations. For the N t = 8 transition, our spectrum results for a · m = 0.1 lie just above this range, which was our motivation for investigating the a·m = 0.05 case with M P /M V = 0.65 (3). In the bigger picture, we see that the M P /M V > 0.9 required for stealth dark matter to produce a first-order transition in the early universe is significantly larger than the masses previously considered. This may have non-trivial implications for the phenomenology of the theory, which we will discuss below and could be explored in future research.

V. CONCLUSIONS AND NEXT STEPS
We have presented non-perturbative lattice investigations of the finite-temperature confinement transition of SU(4) stealth dark matter, motivated by the possibility that this early-universe phase transition could have produced a stochastic background of gravitational waves that may be constrained or discovered by future searches. A first-order transition is required to produce such a stochastic background of gravitational waves, so we have focused on determining the region of parameter space for which the stealth dark matter confinement transition is first order, considering relatively heavy dynamical fermions corresponding to the upper-right corner of the Columbia plot (Fig. 1). The infinite-mass limit reduces to pure-gauge SU(4) Yang-Mills theory, which is known to exhibit a strongly first-order confinement transition [4,34,35]. We analyzed both the puregauge theory and a range of dynamical-fermion masses 0.05 ≤ a·m ≤ 0.4, finding that heavy masses corresponding to a dark meson mass ratio M P /M V > 0.9 are required to produce a first-order stealth dark matter confinement transition.
Focusing on finite-temperature transitions for temporal lattice extent N t = 8, we identified three signals of a first-order transition for which our a · m = 0.4 results exhibit the same qualitative behavior as we observe for the known first-order transition in the puregauge limit, in contrast to our other calculations with a · m ≤ 0.2. First, Figs. 8 and 4 show double-peaked structures in the histogram of Wilson-flowed Polyakov loop magnitude |P L W | measurements, indicating confined/deconfined phase coexistence. Second, the a·m = 0.4 case is the only one for which the |P L W | susceptibility peaks in Fig. 9 grow with the spatial lattice volume L 3 , similar to the pure-gauge peak in Fig. 5 and as required to be consistent with first-order volume scaling χ max ∝ L 3 Finally, the a · m = 0.4 deconfinement fraction results in Fig. 10 are the only set that resemble the pure-gauge case in Fig. 6 and could be consistent with a discontinuity developing in the L → ∞ thermodynamic limit as required for a first-order transition.
Concluding that heavy bare fermion masses a·m > 0.2 are required in order to obtain a first-order N t = 8 confinement transition, we carried out zero-temperature dark meson spectroscopy calculations to translate this into the constraint M P /M V > 0.9 for the dimensionless dark meson mass ratio. We therefore predict that stealth dark matter will produce a stochastic gravitational wave background only for dark fermion masses significantly heavier than those considered by previous lattice studies of stealth dark matter [1,2], which corresponded to 0.55 M P /M V 0.77. Even in this heavy-mass regime the dynamical fermions play a significant role, as shown by the mass dependence of the critical coupling in Fig. 7 and the height of the |P L W | susceptibility peaks in Fig. 9 compared to Fig. 5. However, such dark fermion masses much larger than the confinement scale, as implied by these large M P /M V > 0.9, may result in stable dark glueballs that contribute to the relic density [5], potentially requiring reconsideration of the phenomenology and constraints reported by Refs. [1,2].
Of course, as discussed in Section I, we are considering stealth dark matter in the N f = 4 limit where all four dark fermions have the same mass. While only a small splitting between two pairs of degenerate fermions is required by Big Bang nucleosynthesis, such a splitting could in principle be quite large, without running afoul of other constraints. For such an N f = 2 + 2 theory, the lighter pair of fermions should produce a smaller meson mass ratio M P /M V in the first-order transition region, which needs to be kept in mind when applying collider constraints. In the future it may be interesting to carry out dedicated finite-temperature lattice calculations exploring this more general N f = 2 + 2 setup. While this could be done by taking the square root of the staggered-fermion lattice action used in this work, switching to domain-wall fermions should also be considered.
Turning back to the N f = 4 case, we can compare our results for M P /T c in Table I with the SU(3) endpoint value M P /T c 10 reported by Ref. [37] (for N f = 2 and N f = 2 + 1). Based on our conclusion that the heavymass line of N t = 8 first-order transitions turns into a continuous crossover between 0.2 < a·m < 0.4, we predict an SU(4) endpoint value between 7 M P /T c 10. This is not significantly different than the SU(3) value, though only rough comparisons are possible given that different lattice actions are used and continuum extrapolations have not yet been completed in either case. In particular, Ref. [37] reports significant changes upon moving from N t = 4 to N t = 6, both of which are smaller than the N t = 8 we consider here.
With M P /M V > 0.5, the strongest constraint on the stealth dark matter model is M V 2 TeV coming from Z → + − searches [7]. By using the result a·M V ≈ 1.3 for a·m = 0.4 in Table I, we can relate to translate this constraint into an estimate for the minimum stealth dark matter critical temperature required to produce gravitational waves, Recalling from Section I that the dark baryon may have a mass of hundreds of TeV, we can consider a rough upper bound for M V M DM /2 also in the range of hundreds of TeV, which would imply a critical temperature of tens of TeV. If supercooling effects are mild enough that the transition temperature T * is not too much lower than this equilibrium critical temperature, then the peak frequency of the gravitational wave spectrum would likely correspond to a range of frequencies well suited to be probed by the LISA observatory [8,23,24] and the proposed future Einstein Telescope [81]. While the discovery of such stochastic gravitational waves from the early universe would of course be very exciting, even constraints on their spectrum would place novel new bounds on the viable parameter space of stealth dark matter, likely going beyond what may be possible at colliders and directdetection experiments.
Looking beyond the predictions discussed above, now that we have located a first-order stealth dark matter confinement transition, the next stage of our work will be to study it in more detail in order to more robustly predict the spectrum of gravitational waves it would produce. The key parameters we will investigate are the latent heat, the phase transition duration, and the bubble wall velocity. To this end, we have begun nonperturbative lattice analyses of the latent heat, which will reuse some of the ensembles we have presented here, in addition to more calculations with larger L and N t in order to extrapolate to the L → ∞ thermodynamic limit and the N t → ∞ continuum limit. As N t increases and the lattice spacing at the confinement transition decreases, we will need to work with smaller a·m to stay on a line of constant physics with approximately fixed M P /M V , which will also add to the numerical costs of these larger-volume calculations. It will also be challenging to establish robust non-perturbative constraints on the phase transition duration and bubble wall velocity for the first-order stealth dark matter transition, but even so our lattice calculations should be able to provide new insight into those quantities. In parallel, it will also be valuable to explore alternative approaches to analyzing such first-order phase transitions, such as density-ofstates techniques [74].