Gravitational Waves from Dark Yang-Mills Sectors

Dark Yang-Mills sectors, which are ubiquitous in the string landscape, may be reheated above their critical temperature and subsequently go through a confining first-order phase transition that produces stochastic gravitational waves in the early universe. Taking into account constraints from lattice and from Yang-Mills (center and Weyl) symmetries, we use a phenomenological model to construct an effective potential of the semi quark-gluon plasma phase, from which we compute the gravitational wave signal produced during confinement for numerous gauge groups. The signal is maximized when the dark sector dominates the energy density of the universe at the time of the phase transition. In that case, we find that it is within reach of the next-to-next generation of experiments (BBO, DECIGO) for a range of dark confinement scales near the weak scale.


Introduction
The first direct observations of gravitational waves (GWs) [1][2][3] and the prospects for increasing experimental sensitivity in the next decades have put us on the precipice of a new era of multi-messenger astrophysics and cosmology. Gravitational wave experiments are poised not only to provide direct probes of energetic astrophysical phenomena, such as binary black hole mergers, but also provide a new window into the early universe via the measurement of a stochastic background of gravitational waves. Such a background may have a number of origins, including inflation, topological defects and cosmological firstorder phase transitions (PTs) [4,5]. In the visible sector, both electroweak (EW) gauge symmetry and the approximate chiral symmetry of QCD were spontaneously broken during phase transitions; these are, however, known not to be first-order 1 .
Gravitational waves may also shed crucial light onto dark sectors. The existence of a dark matter component in our universe indicates that unknown particles might be hiding from observation. In optimistic scenarios, experiments may be sensitive to dark sectors that couple to the visible sector via portals that are not significantly suppressed. However, Nature may not be so forgiving: such portals may simply not exist, in which case the only interactions between the dark and visible sector are gravitational, a possibility that is unfortunately consistent with all current data. In this context, gravitational wave searches might be necessary to determine the properties of dark sectors.
In this paper we study stochastic gravitational waves produced during the confinement transition in pure Yang-Mills dark sectors. These are, of course, some of the simplest non-Abelian gauged dark sectors that might exist, but they are also well-motivated in string theory. For instance, dark gauge sectors naturally arise in the ten-dimensional E 8 × E 8 heterotic string itself [8], orbifold compactifications thereof [9][10][11][12][13][14], its free fermionic realizations [15,16], and it smooth Calabi-Yau compactifications [17][18][19]; on G 2 compactifications of M-theory [20][21][22]; and on seven-branes in F-theory. In fact, the latter provides the greatest evidence for dark gauge sectors: the three largest concrete F-theory ensembles [23][24][25], which dwarf the rest of the currently known string landscape, all exhibits tens or hundreds of gauged dark sectors 2 . In the simplest cases many of these factors are pure super-Yang-Mills sectors (from so-called non-Higgsable clusters) that flow to Yang-Mills sectors below the SUSY breaking scale; in the most common scenarios, the gauge groups are low rank SU (N ) groups, G 2 , F 4 , and E 8 . Therefore in addition to being particularly simple extensions of the SM, dark Yang-Mills sectors are also well-motivated by ultraviolet considerations.
Pure Yang-Mills theories are expected to produce a first-order confining phase transition for almost all gauge groups of interest. 3 However, since there is not a first-principles description (other than lattice simulations) of the phase transition, many effective models have been considered in the literature. Among those are quasi-particle models [28][29][30][31], approaches based on the functional renormalization group [32][33][34][35], Polyakov loop models [36,37], as well as so-called matrix models [38][39][40][41][42][43][44][45][46]. The latter are particularly interesting as they can be applied to any gauge group G, the only input needed being the structure of the Lie algebra associated with G.
The outline of this work is as follows. We first discuss, in Section 2, the relevant symmetries for the construction of a matrix model of the confinement phase transition of pure Yang-Mills dark sectors along with results of lattice simulations that can be used to constrain it. We apply these considerations to construct the effective potential in concrete examples in Section 3. The familiar case of SU(N ) is treated in detail and contrasted with the exceptional groups G 2 and F 4 , also expected to confine in first-order phase transitions. Section 4 then estimates the gravitational wave signal emitted during the confining transitions, accounting for theoretical uncertainties and determining their potential for detection in future experiments. We show that these transitions are not long-lasting, so that the GW signal emitted during the PT is suppressed, being only accessible to next-to-next generation searches. We end on Section 5 with a summary and conclusions.

Symmetry constraints and lattice
The confining phase transition in pure Yang-Mills theory can be described by an effective model based on the relevant order parameter, the Polyakov loop. In this section, we 2 In addition to motivating dark gauge sectors, these F-theory ensembles also motivate studies of axionlike particles, see e.g., [26,27]. 3 A notable exception being SU (2). set the stage for the construction of a matrix model of (de)confinement in the absence of quarks, discussing the relevant symmetries as well as describing how lattice observables can be used to constrain the form of the effective potential.

Symmetries of the effective potential
The order parameter for the confinement phase transition is the expectation value of the Polyakov loop l, the normalized trace of the thermal Wilson line L, in the fundamental representation of the gauge group with g being the gauge coupling, β the inverse of the temperature, T a the generators of g in the fundamental and d f the dimension of the fundamental representation. Following the phenomenological approach of Refs. [38][39][40], we consider an effective potential V (L) for which the variables are the eigenvalues of the Wilson line L, referred to as a matrix model for confinement. This type of model can correctly describe the order of the phase transition for the gauge groups and, as we will see, allows for appropriate fits of thermodynamic observables studied on the lattice. For simplicity, we take the time component of the vector potential to be constant This component can always be diagonalized by a gauge transformation, so we take it to be an element of the Cartan subalgebra h of the Lie algebra g associated with the gauge group G. The Cartan subalgebra is defined as the maximal subalgebra of mutually commuting generators. If {H 1 , H 2 , ..., H r } is a basis for h (with r being the rank of G), a general element H ∈ h can be written as H = q i H i , with q 1 , ...q r being coordinates in the Cartan subalgebra.
Below the critical temperature 5 T c , the system is at a confined phase and the expectation value of the Polyakov loop vanishes identically 6 l = 0, while above T c this order parameter becomes non-zero l = 0. Therefore, the effective potential has to be such that the (de)confinement phase transition is accompanied by spontaneous breakdown of center symmetry, so it should be invariant under center transformations. For SU(N ), center transformations are of the form z k = exp(2πik/N ) (2.2) at the Lie group level, with k = 0, 1, ..., N − 1. The thermal Wilson line in the fundamental transforms as L → z k L, so that such transformations act on the elements of the Cartan subalgebra as H → H + k diag(1, 1, ..., −(N − 1))/N . Center symmetry is, however, absent for gauge groups with trivial centers, such as G 2 , F 4 and E 8 .
inequality L(A0) ≤ L( A0 ) explicitly demonstrated in the functional renormalization approach (e.g., in [32-35]). A less simplified model should be able to distinguish these two order parameters. 5 This temperature is of the same order as the confinement scale Λ at which the running gauge coupling diverges. For example, in the case of SU(N ) one has Tc ∼ 1.5Λ [47,48]. 6 The fundamental Polyakov loop is related to the free energy F qq/2 of a static quark-antiquark pair at infinite distance by l ∼ exp(−βF qq/2 ). In the confined state, F qq/2 → +∞, so l → 0.
The roots α of the Lie algebra g are linear functions on the Cartan subalgebra h, defined by the commutation relations E α being elements of g (E α ∈ h for non-zero roots) denoted root vectors. Note that in the first equality the roots are elements of the dual space h * , but they can be mapped one-to-one into elements of h, as done in the second equality, if one takes α(·) ≡ α, · , with ·, · denoting the Killing form. Any root can be written as a linear combination of the elements in a set ∆ = {α 1 , ..., α r } with integer coefficients that are all either non-negative or non-positive. The elements of ∆ are called the positive simple roots of g.
With roots α as elements of the Cartan subalgebra, one can consider the group of reflections w α i about the hyperplanes perpendicular to each simple root α i , known as the Weyl group W . It can be shown that a Weyl transformation maps roots into roots, so that a particular choice of positive simple roots ∆ can be mapped into any other choice ∆ = w α ∆ by such reflections. Therefore, as all choices of ∆ are equivalent, any function with domain in h, such as the effective potential we wish to construct, has to be invariant under the Weyl group.
Another important concept is that of the Weyl chamber, given by the set of H ∈ h such that α i , H ≥ 0, for any positive simple root α i . This defines explicitly ∆-dependent upper half planes in the Cartan subalgebra, so a Weyl transformation maps a Weyl chamber into another. It can be shown that the orbit under Weyl reflections of any point in the interior of a Weyl chamber ( α i , H > 0) has a number of elements equal to the order of the Weyl group. In other words, no Weyl transformation (other than the identity) maps a Weyl chamber to itself.
Writing 7 A 0 ≡ 2πH/βg in terms of H ∈ h, the thermal Wilson line becomes We adopt a basis with elements H i such that their matrix exponential exp(2πiH i ) is either the identity or an element z k of the center for any i = 1, ..., r. With this choice, the effective potential becomes periodic in the q-coordinates with unit period and we identify q i ∼ q i + 1 for each i, as elements connected by a center transformation should give the same value of the potential; its domain can then be restricted to the subset 0 ≤ q i < 1.
Consider, as an illustration, the case of SU (3). The weights in the fundamental can be written as (2α 1 + α 2 )/3, (−α 1 + α 2 )/3 and −(α 1 + 2α 2 )/3 with α 1 and α 2 being the two positive simple roots. Then, a general element of the Cartan subalgebra is written, in the fundamental representation, as Taking q i = α i (H), H 1 = diag(2, −1, −1)/3 and H 2 = diag(1, 1, −2)/3 yields the desired form in Eq. 2.4. Note that this choice of q-coordinates implies that the interval 0 ≤ q i < 1 is entirely contained within a single Weyl chamber, with its boundaries having at least one vanishing q i . In the following, the choice of coordinates q i = α i (H) will be made for all gauge groups. As we have argued above, the effective potential V (q) ≡ V (q 1 , q 2 , ..., q r ) describing the confinement phase transition must be Weyl group-invariant. Weyl transformations generate permutations of all roots, hence the potential has to be invariant under the corresponding permutations of its arguments. For example, any positive root of SU(N ) can be written as α = n i=m α i , with 1 ≤ m ≤ n ≤ r. As a consequence, the potential has to be invariant under permutations 8 of the set {± n i=m q i }, with 1 ≤ m ≤ n ≤ r.

Thermodynamics of the gluon plasma
Our goal is to construct an effective potential that describes the semi quark-gluon plasma (semi-QGP) in the absence of dynamical quarks, i.e., in pure Yang-Mills theories. The region of semi-QGP, which occurs in a range of temperatures from the critical temperature T c to approximately 4T c , is characterized by a sharp increase of pressure starting from approximately zero in the confined phase (in units of the Stephan-Boltzmann limit, i.e., p/p SB ≈ 0, with p SB /T 4 ≡ d A π 2 /45 and d A being the dimension of the adjoint representation) and asymptotically approaching the equation of state for an ideal gas at increasing temperatures.
We focus our attention on the region close to the critical temperature, looking for effective potentials that give the appropriate order for the phase transition and reproduce the behavior of thermodynamic quantities measured on the lattice.
At high temperatures T T c , the effective potential is given by the free energy of a gas of gluons in a constant background field A 0 and can be found perturbatively as [40] V pt at one-loop order. The sum runs over all the roots α of g and the function B 4 (x) is a shifted Bernoulli polynomial with · being the floor function. We take the effective potential in the semi-QGP region to be the sum of the perturbative contribution V pt in Eq. (2.6) and a nonperturbative contribution V npt , that respects the symmetries discussed in Section 2.1. Note also that on the interval 0 < x < 1 the function B 4 is polynomial; it is, however, not analytic at the origin (its third derivative involves the divergent sum ∞ n=1 n −1 ) nor at any integer value of x. Therefore, the perturbative part of the effective potential is polynomial in the interior of a Weyl chamber 9 , but has singular behavior on the hyperplanes perpendicular to the roots, i.e., at the boundaries of the Weyl chambers. To avoid introducing additional singularities, we assume that the nonperturbative part of the potential is also polynomial in the interior of the Weyl chambers. As we will see, this assumption leads to Bernoulli polynomials of all (even) orders as building blocks of V npt . In particular, the shifted Bernoulli polynomial of degree two, given by 1/2 will be used extensively.

Lattice observables
The nature of the confinement phase transition, either continuous or not, can be determined on the lattice from the behavior of the order parameter l at the transition temperature. We are interested in first-order PTs, as these can potentially produce significant stochastic gravitational wave signals [49]. This type of transition involves a discontinuous change in the Polyakov loop at T = T c . Lattice simulations of pure SU(N ) gluodynamics have determined that the confinement phase transition is indeed first-order for N ≥ 3 colors 10 [50][51][52][53]. Similar behavior was also found for gauge groups of the Sp(N ) type [54]. In addition, the phase transition for the exceptional group G 2 was shown to be discontinuous (see, e.g., [55][56][57]), with l ≈ 0 below T c , even in the absence of center symmetry.
For gauge groups with a large number of gluons, there is a large mismatch between the number of degrees of freedom above (gluons) and below (color singlet glueballs) the critical temperature, as the latter is essentially independent of the dimension of the group. As such, one can expect, as conjectured in [58], the confinement phase transition to also be of first order in the case of larger gauge groups, such as F 4 and E 8 , not yet studied on the lattice.
Lattice simulations also seem to indicate the temperature dependence of the nonperturbative part of the effective potential. This can most clearly be seen in the behavior of the interaction measure ∆, defined as where e(T ) is the energy density and p(T ) = −V (q min (T ); T ) the pressure of the gas of gluons with q min (T ) denoting the coordinates of the global minimum of the potential at a temperature T . Above the phase transition, in the interval 1.1T c T 4T c , the interaction measure is observed to be directly proportional to T 2 for all the gauge groups studied on the lattice 11 , as shown in Fig. 1i for the groups SU(N ) with N = 3, 4, 6 and G 2 . In this interval, the interaction measure divided by the square of the temperature is approximately 9 Note that we also need the restriction α, H < 1 for any root α. Henceforth, we use the term Weyl chamber to describe the region defined by the set of inequalities 0 ≤ α, H ≤ 1, with α being any positive root. 10 These calculations were performed only up to N = 8. They, however, show that the first-order transition gets stronger with increasing N ; one then expects that the transition continues to be of first order for arbitrary values of N ≥ 3. 11 The case of SU(2), which displays a second-order phase transition, appears to slightly deviate from this behavior [59]. (ii) Figure 1: The behavior of thermodynamic quantities on the lattice; the data for the interaction measure (i) is taken from [60] for SU (3), from [61,62] for SU(4) and SU(6) (the figure only shows data points from [61] for clarity) and for the exceptional group G 2 adapted from [57]. The data for the renormalized Polyakov loop (ii) of SU(N ) is reproduced from Ref. [63] for N = 3, and [64] for N = 4 and N = 5.
constant with its value per gluon being approximately the same for each group: all data points fall (within error) in the range 0.38 ∆/d A T 2 0.5. In the following, we refer to this region as the expected region, shaded in purple. This behavior indicates that the dominant contribution to V npt should be, according to Eq. (2.8), proportional to T 2 , at least for temperatures right above the PT. Following Ref. [40], we also allow for a temperature independent constant. In addition, we continue the expansion in even power of the temperature and include a term proportional to T −2 with coefficient independent of the coordinates q. As we will see, this extra term allows for a better fit of our model to lattice thermodynamics results.
For SU(N ) groups, measurements of the renormalized Polyakov loop are available from the lattice in the cases N = 3 [63], N = 4 and N = 5 [64]. These are shown in Fig. 1ii. The data points show a similar trend for the different number of colors shown. Thus, we again select an expected region for the value of the Polyakov loop, shown in purple, and make the assumption that the renormalized Polyakov loop approximately falls within this region for an arbitrary number of colors, as well as for other gauge groups. It is fitted by l(q min ) in our model. In addition, the latent heat for the SU(N ) transition was determined in Ref. [61] and can be used to further constrain our effective model for SU(N ), being given by the discontinuity in the interaction measure across the transition . For other gauge groups, we do not impose a value for the latent heat as a constraint, as these are not yet available from lattice studies.
The approximately universal behavior of the interaction measure and of the renormalized Polyakov loop described in this section will be used as a guide for our effective description of the gluon plasma close to the critical temperature. In the next section, we combine the symmetry considerations of Section 2.1 with these lattice results to construct an effective model for the confinement phase transition. We focus on the interval T c T 2T c , since, as explained in the following, this allows for our simplified model to adequately fit the necessary observables.

The effective potential
Given the symmetry and lattice constraints introduced in Section 2, in this Section we construct the effective potentials describing the semi quark-gluon plasma phase that characterizes Yang-Mills theories just above T c . Specifically, we impose center and Weyl group invariance as well as the expectations for thermodynamic quantities such as the interaction measure and for the Polyakov loop inspired by the apparent universality of lattice results discussed above. In Section 3.1, we treat the more familiar case of SU(N ) and contrast it with the exceptional cases of G 2 and F 4 in Section 3.2.

SU(N )
We start by generalizing the choice in Eq. (2.5) of coordinates in the Cartan subalgebra of SU(3) to an arbitrary number of colors N ≥ 3. The positive simple roots of the Lie algebra can be written as Thus, a general element H ∈ h can be written as with the coordinates again chosen as q i ≡ α i (H).
In terms of these coordinates, the perturbative part of the effective potential from Eq. (2.6) can be written as It is explicitly invariant under the Weyl group, as the summatory runs over all roots of g. We can also check invariance under center symmetry. It was shown in Section 2.1 that a center transformation z k acts on the weights of the fundamental as µ i → µ i + k/N for i = 1, ..., N − 1 and µ N → µ N − k(N − 1)/N . One can check, using Eq. (3.1), that this transformation shifts the roots α i by an integer. The function B 4 (x) in Eq. (2.7) has unit period and the change in its argument under a center transformation is an integer for each element in the sum, so V pt (q) is indeed invariant under center symmetry. This agrees with the fact that the adjoint representation has zero N -ality, so that terms constructed from the adjoint Polyakov loop (involving all the roots) should be left invariant by center transformations.
We now consider the nonperturbative contribution to the potential. Combining the symmetries from Section 2.1 with the lattice results discussed in Section 2.3, we assume that it is a Weyl group and center-invariant almost-everywhere polynomial function, with a dominant component proportional to T 2 . By almost-everywhere polynomial we mean a function that, like the Bernoulli polynomial in (2.7), is polynomial except at the boundaries of Weyl chambers. For simplicity, we assume a polynomial of degree four, this being the lowest degree necessary to describe a first-order phase transition as a thermal transition from a metastable vacuum to the true vacuum of the theory, separated by a barrier of finite height 12 . To account for Weyl symmetry, we consider terms of the form α∈W ·α where W ·α denotes the orbit of a rootα under Weyl transformations and P 1 , P 2 , ... are polynomials in the interior of a Weyl chamber of degree (less or equal to) four . For su(N ), this orbit is the set of all roots. For other algebras, roots might have different lengths and, as we will see, one has to include terms summing over distinct orbits.
The periodicity of the coordinates in the Cartan subalgebra constrains the form of the polynomials P . Consider, for example, a term containing one of the coordinates, q i . If one performs a Weyl reflection w α i that takes q i into −q i , followed by the transformation q i → q i + 1, the potential should be left invariant. Note that the resulting transformation, q i → 1 − q i , keeps the coordinate within its restricted domain 0 ≤ q i < 1. Thus, the polynomials P should have the property P (1 − q) = P (q). Bernoulli polynomials of degree n obey B n (1 − q) = (−1) n B n (q), so the ones with even degree form the appropriate basis for our construction. Therefore, the most general terms of the form (3.4) obeying these symmetries are where the multiplying factors are chosen for convenience and the sums run over all roots of su(N ). For clarity, we write these terms explicitly in the case of SU (3), with q 3 ≡ q 1 + q 2 . Note that, since V 1 is of degree two, a term proportional to V 2 1 is also allowed by the symmetries. Such a term is, however, a linear combination of V 2 and V 3 above. The nonperturbative part of the effective potential is then taken to be of the form where c i , d j are coefficients still to be set. As explained in Section 2.3, the temperature dependence of V npt (q) right above the critical temperature is mainly given by a component proportional to T 2 and we write it as a linear combination of terms in Eq. (3.5). As these are the only q-dependent terms, it encodes the dynamics of the phase transition. As adopted in Ref. [40], we allow for a temperature-independent constant 13 d 1 T 4 c , where the factor T 4 c makes the coefficient d 1 dimensionless. Close to the critical temperature, it is reasonable to allow for some physics to give increasing contributions to the effective potential for decreasing T , corresponding to the appearance of terms proportional to negative powers of T ; if such physics does not exist, lattice will fit the coefficients to zero. Thus, we continue the expansion in even powers of the temperature and add a term proportional to T −2 . As we will see, the inclusion of such a term allows for a correct description of the evolution of the Polyakov loop as a function of temperature, while simultaneously fitting other thermodynamic quantities. The model, however, fails to do so if one sets d 2 = 0.
The confined state, the center-symmetric state with vanishing Polyakov loop, has coordinates (q c ) i = 1/N for all i = 1, ..., N − 1, which amounts to having all the eigenvalues of the thermal Wilson line L equally spaced along the unit circle. As observed on the lattice, the confinement transition at T c does not take the system directly to the perturbative vacuum (q d ) i = 0 (i.e., l(T + c ) = 1). Thus, the discontinuous transition happens between a metastable state at q c and another state q t = q d inside the Weyl chamber. The coefficients in Eq. (3.7) are not all independent. First, it is necessary to impose that the phase transition happens at T = T c . In addition, we assume that the pressure of the glueball gas in the confined state vanishes 14 , V (q c ; T c ) = 0.
Before explicitly imposing the constraints discussed above, let us comment on a simplifying assumption, termed the uniform eigenvalue ansatz, i.e., the assumption that the eigenvalues for the thermal Wilson line L on the minimum of the effective potential are 13 Note, however, that the authors of that work do not allow for a V2 term. 14 It is certainly true that the pressure in the confined state is always much smaller than at much higher temperatures, but a nonvanishing value can be measured on the lattice (see e.g, [65]). Its value, however, is small enough that this assumption should not change our results appreciably. equally displaced along a section of the unit circle, for all temperatures. This amounts to taking q i (T ) = (1 − r(T ))/N for any i, with r(T ≤ T c ) = 0 at the confined state and r(T → ∞) = 1 at the perturbative vacuum. This ansatz reduces the problem to a onedimensional thermal transition between two vacua 15 and it will be used for large numbers of colors, N ≥ 8, allowing for an estimation of the thermal transition rate in those cases. Under this assumption, the terms defined in (3.5) become (3.8) Note that the uniform eigenvalue ansatz amounts to having the global minimum of the potential always located at the line that is equidistant from the faces at the boundary of the Weyl chamber. In the next section, we use this fact to generalize the uniform eigenvalue ansatz.
As explained in Section 2.3, Polyakov loop data constrains the form of the effective potential. In particular, we chose the state q t above the transition in such a way that it matches the value for l observed on the lattice (e.g., l(T + c ) ≈ 0.4 for SU(3) and l(T + c ) ≈ 0.5 for SU (4)). For larger numbers of colors, we fit the Polyakov loop in our model to match the mid-sectional curve on the expected region shown in Fig. 1ii.
Finally, we impose that the latent heat of the transition agrees with the values found in Ref. [61], where we made the approximation ∆(T − c ) ≡ e(T − c ) − 3p(T − c ) ≈ 0, as both the energy density e and the pressure p are negligible in the confined phase. When applied to the nonperturbative potential in Eq. (3.7), these constraints, along with the ones discussed previously, reduce the number of independent coefficients from six to two. The remaining coefficients are then found by fitting to lattice data, the results for N = 3, 4, 6 being shown in Fig. 2. The numerical values for all coefficients in Eq. (3.7) are shown in the Appendix. It is clear that the model can quantitatively describe lattice thermodynamics in the interval of interest, from T c up to approximately 2T c . For temperatures not in this range, our model gives wildly incorrect or even unphysical results (e.g., negative pressure). This is, of course, a result of trying to describe a strongly interacting system with a potential that can be nicely written down as a sum of a non-perturbative and a one-loop order perturbative term. Thus, we content ourselves with the less ambitious goal of trying to model the region close to the PT.

G 2 and F 4
Now we generalize the discussion from the previous section to the exceptional groups G 2 and F 4 . First, these groups have trivial centers and, also, the set of all roots for the lie algebras g 2 and f 4 are now divided into two sets, of long roots and of short roots, which do no mix under Weyl reflections. Both these facts combined decrease the amount of symmetry that can be imposed in the structure of the effective potential and, as a consequence, more terms are allowed in its construction.
Starting with G 2 , the positive simple roots can be written as a linear combination of the weights in the lowest-dimensional representation 16 (7) as α 1 = µ 1 − µ 2 and α 2 = −µ 1 (with α 1 , α 1 > α 2 , α 2 and µ 1 , µ 2 are weights in the fundamental), so that a general element of the Cartan subalgebra is given by, in the fundamental representation, again with q i ≡ α i (H). Note that a trivial center requires the matrix exponentials exp(2πiH i ) to be the identity, so the diagonal entries of the matrices H i have to be integer numbers. When written in terms of the positive simple roots, the sets of positive long and short roots are 17 , respectively, α L = W · α 1 = {α 1 , α 1 + 3α 2 , 2α 1 + 3α 2 } and α S = W · α 2 = {α 2 , α 1 + α 2 , α 1 + 2α 2 }. Thus, the possible terms of the form (3.5) can now have a sum running on either one of these sets of roots, i.e., are the building blocks for the effective potential. The nonperturbative polynomial contribution to the effective potential can then be written as (3.12) The boundary of a Weyl chamber of g 2 is defined by the vanishing of the Killing form with the two positive simple roots, which are of different lengths, i.e., a point on the boundary obeys α i , H = 0 for i = 1 or 2. Therefore, the root system lacks the symmetry necessary for an assumption similar to the uniform eigenvalue ansatz, adopted in the previous section. Therefore, the effective potential V (q 1 , q 2 ) is necessarily two-dimensional; even if we impose that, initially, the global minimum of the potential lies equidistant from each hyperplane at the boundary of the Weyl chamber, the subsequent dynamics violates such condition.
The confined state, as in the case of SU(N ), is seen on the lattice [66,67] to have a very small value of the traced Polyakov loop in the fundamental representation, l(T − c ) 1, which we take to vanish identically. Note that, as opposed to the case of SU(N ), this order parameter does not necessarily vanish below the critical temperature, as center symmetry is absent 18 . A priori, any element of the Cartan subalgebra h with vanishing Polyakov loop can be taken as the confined state. This set defines a line in the Cartan subalgebra of g 2 on which we allow the confined state to be located, shown in blue in Fig. 3 along with the interior of a Weyl chamber.
Having constructed the potential, as the sum of Eq. (2.6) and (3.12), we then proceed as done in the case of SU(N ) and impose the following constraints. First, at T = T c the global minimum of the potential jumps discontinuously, as the temperature is raised, from the confined state q c to a state with coordinates q t , both with (approximately) vanishing pressure at that temperature. As mentioned previously, q c is randomly chosen subject to the condition l(q c ) = 0 and, based on the behavior of lattice data for SU(N ), we choose (also drawing randomly) the state q t so that 0.38 l(q t ) 0.55 (this interval is taken from the expected region of Fig. 1ii). The region inside the Weyl chamber in Fig. 3 that obeys this bound is shown in red. Once both q c and q t are chosen, these conditions reduce the ten coefficients in Eq. (3.12) to four, which are then fitted by the lattice data (only available in the case of G 2 ) or expected lattice behavior. Specifically, for the observables not yet calculated on the lattice, we fit the model to the midsection of the expected regions in Figs. 1i and 1ii. We also selected the potentials that give values for the pressure that are as close as possible to zero in a temperature range δT ∼ 0.1T c right below the critical temperature. This is imposed in an attempt to extrapolate the model to temperatures slightly below T c , so that the gravitational wave signal can be reliably calculated. A similar construction can be made for the group F 4 . Its positive simple roots can be written as [71,72] again with q i ≡ α i (H) and the basis {H 1 , H 2 , H 3 , H 4 } having integer elements on the diagonal (the numbers after the ellipsis are determined by writing the additional weights µ as linear combinations of the positive simple roots 19 ). One can then again divide the roots into sets of long roots α L and short roots α S and construct the possible terms in the nonperturbative potential as in Eq. (3.11).  The boundaries of the Weyl chamber are now defined by the vanishing of the Killing form with respect to the four simple roots, two of which are long (α 1 and α 2 ) and the other two short (α 3 and α 4 ). The potential should be invariant under the Weyl group, which includes transformations that permute each pair (long or short) of positive simple roots. This allows for a simplifying assumption generalizing the uniform eigenvalue ansatz described in the case of SU(N ): we can take the minima of the potential to be always located at the plane that is equidistant from the Weyl chamber boundary hyperplane defined by the two long roots and also from the hyperplane defined by the short roots 20 . In other words, we can project the potential to the plane defined by q 1 = q 2 and q 3 = q 4 , reducing the dimensionality of the effective potential from four to two. We emphasize that this is not a necessary assumption. However, both in the case of the uniform eigenvalue ansatz for SU(N ) as well as for its generalized version in the case of F 4 , the model can accurately fit the (expected) behavior of thermodynamic quantities from lattice, so hopefully not much is lost by our assumption.
The resulting curves for the fits to the interaction measure and the Polyakov loops for both G 2 and F 4 are presented in Fig. 4. These plots show the resulting best-fit curves for many different choices of confined state q c and q t . As the figure shows, we were able to construct a number of effective potentials that reproduce the (expected) lattice behavior.

Stochastic gravitational wave signal
Equipped with effective potentials for the semi-QGP phase in pure Yang-Mills theories, in this Section we compute stochastic gravitational wave spectra produced during the associated confinement transitions and study their possible observation in planned experiments.
First-order phase transitions in the early universe are well-known sources of a stochastic gravitational wave background [73][74][75]. This type of transition proceeds via nucleation and subsequent expansion of bubbles of the true ground state of the theory on a background in the metastable vacuum. There are different mechanisms that can generate gravitational radiation during a first-order PT (for a recent detailed description, see the reviews [76,77]); first, gravitational waves are produced during the collision of the expanding bubbles and, subsequently, the energy released to the thermal plasma by the transition generate sound wave and magnetohydrodynamic turbulence contributions. The scalar field contribution from bubble collisions is subdominant in the case of a nonrunaway PT, in which the bubble wall reaches a finite terminal velocity due to friction exerted by the thermal plasma. In that case, the fraction of the latent heat that becomes kinetic energy of the scalar field is vanishingly small, which renders the contribution of bubble collisions to the stochastic GW signal negligible.
In the case of a confining PT in a dark gauge sector without matter, i.e. dark Yang-Mills, the order parameter is the Polyakov loop, which is constructed out of the temporal component of the non-Abelian vector potential (see Eq. (2.1)). Thus, the scalar field generating the bubbles of true vacuum should interact strongly with the thermal plasma surrounding them. Therefore, we assume that the confinement phase transition proceeds via nonrunaway bubbles and only account for the contribution of sound waves and turbulence to the stochastic background of gravitational waves. This is in agreement with the results of Ref. [78], which shows that transition splitting, radiation emitted from gauge bosons acquiring a mass when traveling across the bubble wall from symmetric to broken phase, generates enough friction to impede the runaway of the bubble. In a confining PT, gauge bosons go from a deconfined to a bound state when crossing the bubble wall, and it is thus reasonable to expected that enough friction is generated by the plasma and the phase transition should follow the nonrunaway case.

Parameters of the phase transition
Before calculating the spectrum of gravitational waves from the first-order PT, a number of parameters, particular to each physical model, have to be determined.
First, the strength of the PT is encoded in the parameter α, determined by the ratio of the change in the interaction measure ∆ across the phase transition to the total thermal energy density of the universe in the symmetric phase, given in terms of the enthalpy w ≡ e + p, as calculated at the nucleation temperature T n , where the signs refer to the symmetric (+) and broken (−) phases. This temperature is the one at which there is on average one bubble of the confined phase nucleated per Hubble volume, which implies with M Pl being the reduced Planck mass and g * is the number of relativistic degrees of freedom at T = T n . For small amounts of supercooling T n T c (which turns out to be the case for the transitions considered here), the nucleation temperature is also very close to the bubble percolation temperature T * at which the PT can successfully complete and GWs are produced. Henceforth, we take T c ≈ T n ≈ T * .
Another important parameter is the inverse duration of the PT, defined as with H * the Hubble parameter at T * and S 3 is the action for the O(3)-symmetric bounce solution for a thermal transition between the metastable and the true vacua. In the case of pure Yang-Mills described in Section 3, this action can be obtained by going one step further in the operator expansion, adding the gauge kinetic term at leading order, as where ρ and Ω are, respectively, the three-dimensional radial coordinate and solid angle and F µν ≡ F a µν T a . In addition to α and β, one should also determine the bubble wall velocity v w as well as the efficiency factors κ v (v w , α d ) and κ tb (v w , α d ) for conversion of latent heat into bulk and turbulent motion, respectively. A proper determination of the bubble wall velocity necessitates a treatment of the dynamics of the bubble expansion, with an appropriate modeling of the friction terms (see, e.g. [79] for a discussion of effects that contribute to this dynamics). This is, however, outside of the scope of this work and we assume a relativistic bubble wall velocity v w 1, expected to hold for values of α not much smaller than O(1). On the other hand, the efficiency factors are not additional parameters, as they depend exclusively on the bubble wall velocity and on with w d denoting the enthalpy in the dark sector only. In the limit v w 1, one has [80] κ Note that, if the dark sector dominates the energy density of the universe at the time of GW production, we get α ≈ α d and a single parameter suffices. Moreover, if the phase transition is fast enough, i.e., β/H * 1, which we show to be the case for the pure Yang-Mills confining transition in the following sections, bulk motion quickly becomes turbulent and one can take κ tb κ v [81].

Energy budget and glueball-dominated phase
We now address some cosmological considerations that maximize the GW signal. For fixed values of the parameters discussed in the previous section, the energy density in the form of gravitational waves ρ gw, * emitted during the PT in a given dark sector is directly proportional to the radiation energy density in that sector ρ d rad, * prior to confinement. If the universe is radiation-dominated at the time of the transition, the GW density parameter right after production obeys with ρ other rad, * being the radiation energy density in other sectors (which include, of course, the visible sector). With the ratio ρ gw, * /ρ d rad, * fixed, a maximal signal is obtained if the dark sector going though the confinement phase transition dominates the energy density of the universe at the time of the transition, i.e., for ρ d rad, * ρ other rad, * , which we assume from here on 21 .
As mentioned in Section 4.1, if most of the energy density is in the dark sector at the time of the PT, one has α ≈ α d . Now, lattice results show that the pressure of the gluon gas in the semi-QGP is negligible below T c and that the pressure is continuous across the confinement PT, so the parameter α is reduced to as the energy density in the confined state is small compared to its value above the critical temperature and T n ≈ T c . In general, if the energy density in other sectors cannot be neglected, we get α < α d ≈ 1/3 as w d (T + n ) ≤ w(T + n ) in Eq. (4.5). This decrease in the value of α further suppresses the GW signal when compared to the case in which the dark sector is dominant.
The density parameter of gravitational waves redshifted to today, Ω gw , depends on the detailed evolution of the Hubble parameter since the time of GW production. In particular, if a confining dark sector dominates the energy density at the time of the PT, one expects to have a period of matter domination 22 after confinement occurs, with most of the energy density of the universe in the form of dark glueballs, ρ gb, * ≈ ρ d rad, * . If that happens before big bang nucleosynthesis (BBN), the glueballs ultimately have to decay (mostly) to radiation in the visible sector before the onset of BBN, as a persistent early matter domination phase would spoil its predictions. For simplicity, we assume that glueballs decay directly to visible sector radiation at some later time 23 , when the scale factor is a τ . During their lifetime, the 21 For cases in which the dark sector is cold with respect to the visible sector or when many sectors contribute to the total energy density, see [82,83] and [84], respectively. 22 That is not exactly true, as glueball 3 → 2 self-interactions, while still active, make them redshift slightly faster than matter; the correction factor is, however, a slowly varying logarithm in the scale factor, i.e., ∝ log(a) [85][86][87]. 23 We assume an instantaneous decay of the glueballs, as that is sufficient to estimate the order of magnitude of the entropy exchanged between the sectors. For a more careful treatment, see e.g. [88].
energy density in glueballs increases as ∝ a relative to the GW energy density, so that at the time of decay ρ gw,τ ρ gb,τ = a * a τ ρ gw, * ρ gb, * ≈ a * a τ ρ gw, * ρ d rad, * . (4.9) with a * being the scale factor at the time of bubble percolation. Then, the energy density in glueballs is transferred to visible sector radiation, so that right after their decay ρ v rad,τ = ρ gb,τ , and the density parameter in GWs becomes where we used Eq. (4.7) in the limit ρ d rad, * ρ other rad, * . Thus, a longer period of early matter domination means a stronger suppression of the GW signal 24 , compared to the case in which the dark glueballs decay to visible sector radiation almost immediately after the PT.
We estimate the maximum amplitude of the GW spectrum by assuming that the glueballs decay quickly to the visible sector, in such a way that the factor a * /a τ in Eq. (4.10) is approximately one. This amounts to a situation in which the lifetime of glueballs with respect to SM decays is much shorter than the age of the universe at BBN, T BBN ∼ O(min). Given the requirement of gauge symmetry, the lowest dimension operator connecting the dark and visible sectors is of dimension six [90], of the form with H the SM Higgs doublet and M being the mass scale of the degrees of freedom connecting visible and dark sectors 25 . For a confinement scale 26 of Λ ∼ 100 GeV, the lifetime of glueballs is smaller than 1s for M 10 8 GeV (see, e.g., Fig. 2 of [90]). Therefore, we assume the presence of the higher-dimensional operator in Eq. (4.11) with Λ ∼ 100 GeV << M 10 8 GeV, so that dark glueballs decay quickly enough and our description of the dark sector as pure Yang-Mills is justified. Given the assumptions above, of instantaneous confinement transition and glueball decay with negligible lifetime, we can relate the the temperature of the visible sector plasma right after glueballs decay (which also coincides with bubble percolation), T v * , to the confinement scale in the dark sector Λ d T c . The energy density originally in dark radiation is ultimately transformed into energy in the visible sector plasma, so that with g d * and g v * the number of relativistic degrees of freedom in the dark and visible sectors, respectively, at percolation. For finite glueball lifetimes, the temperature T v gb in the visible sector right after the decay is given by Eq. (4.12) multiplied by the factor (a * /a τ ) 3/4 . 24 For more on the effect of matter domination on the gravitational wave signal, see [89]. 25 Such an operator can be generated by integrating out either scalar and fermionic mediators with masses ∼ M that couple to the SM Higgs and are charged under the dark gauge group, see [90]. 26 This value is taken here for the dark confinement scale since it is the one that maximizes the projected reach of the future GW searches, i.e., BBO and DECIGO.

Energy density in gravitational waves
As discussed at the beginning of this section, gravitational waves in a nonrunaway PT are produced both by sound waves and turbulence, the former giving a larger contribution. The total energy density produced is the sum of the terms [76,77,91] with R * = (8π) 1/3 v w /β being the mean bubble separation at percolation, c s ∼ 1/ √ 20 [92] the speed of sound in the plasma at T c ,Ω gw ∼ 10 −2 a numerical factor obtained from simulations and K the fraction of kinetic energy in the plasma, given by Note that the expression for the sound wave contribution in Eq. (4.13) already takes into account the suppression factor for short-lasting PTs, recently discussed, e.g., in [93]. In addition, the spectral shape functions have the form with h * = a * H * /a 0 the inverse Hubble time at percolation redshifted to today and the the peaks are at the frequencies (redshifted to today) with z p 10 obtained in numerical simulations. The remaining factor in Eqs. (4.13) and (4.14), F gw , accounts for the redshift of the amplitude of the GW density parameters from the time of emission to today, being therefore sensitive to assumptions about the intermediate cosmic evolution of the universe. As discussed in Section 4.2, we neglect the lifetime of the glueballs produced in the confining transition, so that the universe follows the standard cosmic evolution after glueball decay and with Ω γ,0 the energy density in photons today and g v s, * , g v s,0 the effective numbers of entropic degrees of freedom in the visible sector at after glueball decay and today, respectively, and g v 0 is the number of relativistic degrees of freedom today.
Note that the expressions in Eqs. (4.13) and (4.14) are valid for sources that are not long-lasting, i.e., for β/H * O(1). To calculate the sound wave contribution, we use the web-based tool PTPlot introduced in [77].

Results
Before using the CosmoTransitions package, the action in Eq. (4.4) has to be written in terms of canonically normalized fields ϕ as For example, for the gauge groups discussed in Section 3, these fields can be written written as where the variable r introduced in the discussion above Eq. (3.8). The canonically normalized fields for SU(N ) without the uniform eigenvalue ansatz can also be easily obtained, but are not shown explicitly here. The action S 3 as a function of temperature can then be calculated for each value of the dark coupling constant α s (T c ) = g 2 /4π at the critical temperature, which we assume to be in the interval α s (T c ) ∈ (0.2, 0.4) for all gauge groups 27 . Examples are shown in Fig. 5 for SU(N ). Fig. 5i shows how the action S 3 changes with the choice of α s (T c ) and 5ii shows S 3 for N = 3 (in blue), 4 (red), 6 (green) and 8, 10, 16, 32 (all falling on the gray band) with α s (T c ) = 0.3. For N < 8, we find the action S 3 with the full potential, without the uniform eigenvalue assumption, while for N ≥ 8 that we use the simplified one-dimensional potential, using Eq. (3.8). From these results for the action S 3 as a function of temperature, one can determine the values of the inverse duration β from Eq. (4.3) for the confining transition in each case, the result being shown in Table 2. Note that for all numbers of colors ≥ 8, the curves for the action S 3 (T ) as a function of temperature are approximately identical, so that the value of β should not vary significantly for large values of N . This is to be expected; the Ndependent terms in Eqs. (3.8) and (3.9) that cannot be absorbed by the free coefficients in Eq. (3.7) are all O(1/N 2 ). As a consequence, each term in Eq. (3.8) is (up to normalization) approximately independent of N for large N .
For G 2 and F 4 , the value of the parameter β depends on the choice of states q c and q t , defined in Section 3. For the models shown in Fig. 4, the distribution of values of β is shown in Figs. 6i and 6ii for the values of the dark coupling constant α(T c ) = 0.1, 0.3 and 0.5. These distributions seem to be independent of the value of the coupling for G 2 and   β/H * SU(3) (5.1 ± 0.6) × 10 4 SU(4) (2.9 ± 0.6) × 10 4 SU(6) (7.7 ± 2.2) × 10 4 SU(8) (4.0 ± 0.8) × 10 4 Table 1: Values of the parameter β/H * for SU(N ). Errors are estimated by calculating β with different coefficients c i and d i which also give good fits to lattice data.
to have a slight dependence on the coupling for F 4 , with peaks at β/H * ∼ 10 4 in all cases except F 4 with α(T c ) = 0.1, which peaks at β/H * ∼ 10 5 . Similar values of β/H * were also observed in other effective models, e.g., describing the chiral phase transition in confining dark sectors with matter [96,97]. This gives a GW signal orders of magnitude smaller than estimated with more optimistic choices for the duration of the PT, e.g., in [49]. (See [98,99] for related works that also discuss gravitational waves from confining PTs.) Once the parameter β is determined, the gravitational wave signal can be calculated from Eqs. (4.13) and (4.14). The resulting spectra for T v * = 100 GeV and g v * = 100 (the approximate number of relativistic degrees of freedom in the standard cosmic evolution at T v * ) are shown in Fig. 7, along with projected experimental sensitivities for next-generation (LISA), taken from PTPlot, and next-to-next generation (BBO and DECIGO) gravitational wave searches, adapted from Refs. [100][101][102][103][104]. The shaded strips represent the uncertainty estimated by varying the best-fit coefficients of the model in Eqs. (3.7) and (3.12) as well as different choices of q c and q t in the case of G 2 and F 4 . Additional uncertainty comes from varying the dark coupling constant in the interval α s (T c ) ∈ (0.2, 0.4). The signal is within range of next-to-next generation searches if T v * ∼ 100 GeV, although many orders of magnitude out-of-reach of LISA. For different values of T v * , the maximum amplitude of the signal changes only slightly (due to a small change in T n ); from Eq. (4.13), the energy density of GWs can be seen to depend on the combination H * R * ∼ H * /β, which is independent of T v * in our model. On the other hand, the peak frequency varies linearly with T v * (Eq. (4.18)), so that any significant deviation from T v * ∼ 100 GeV pushes the signal out of the observable range.

Summary and conclusions
In this work we studied stochastic gravitational wave backgrounds produced by confining phase transitions in dark Yang-Mills sectors. This requires constraining effective potentials by symmetry and lattice considerations, constructing them for concrete simple Lie groups, and then computing the gravitational wave signal. We describe each in turn.
In Section 2 we set the stage for constructing an effective matrix model for the semi quark-gluon plasma in pure Yang-Mills theories. To do so, we discuss the necessary symmetries (center symmetry and Weyl group invariance) as well as lattice observables that constrain the effective potential on thermal Wilson line eigenvalues that take values in the Cartan subalgebra and serve as order parameters for the phase transition.  (6), G 2 and F 4 . This is of the same order as the dark confinement scale T c ∼ Λ, although the exact relationship depends on the gauge group and requires lattice calculations (see footnote 5). The uncertainty comes mostly from the different choices of good fits for SU(N ) and of states q c and q t in the case of G 2 and F 4 . The value of the dark coupling constant at the critical temperature is taken to be in the range 0.2 ≤ α s (T c ) ≤ 0.4, which also contributes to the uncertainty.
In Section 3, these constraints were implemented, yielding (for each group) an effective potential that models the behavior of a strongly coupled gas of gluons close to the confinement phase transition. Such behavior is determined by currently available lattice data for SU(N ) gauge groups with small numbers of colors (N = 3, 4, 6) as well as for G 2 . The universality observed in this data was used to extend the matrix model to describe the exceptional gauge groups F 4 , as well as SU(N ) with larger numbers of colors. Assuming universal thermodynamic behavior, we showed that a simple effective model can appropriately describe observable quantities such as the interaction measure and the renormalized Polyakov loop for all gauge groups considered.
Equipped with the effective potentials for the confinement transitions, we computed the stochastic gravitational wave background in Section 4. This requires the determination of the action of the bounce solution in a thermal transition between a confined state and a (partially) deconfined one, which in turn allowed for an estimation of the gravitational wave signal. For all gauge groups considered, the GW signal is only accessible to futuristic experimental searches such as BBO and DECIGO, being many orders of magnitude below the projected reach of LISA. This happens because the PT is not long-lasting, having an inverse duration parameter β/H * ∼ 10 4 or larger, suppressing the GW energy density emitted by sound waves in the plasma. In addition, this signal is only visible when the glueballs resulting in a dark sector decay to visible sector radiation at a temperature T v * ∼ 100 GeV. For temperatures not of this order of magnitude, the spectrum's peak frequency takes the signal out of the range of observation of both BBO and DECIGO. Interestingly, T v * ∼ 100 GeV occurs when the dark confinement scale is near the weak scale. Though this signal is relatively weak, its interest derives from the fact that the only surefire model-independent way to detect dark sectors is gravitationally. As much as we might wish for stronger portals, they simply may not exist, a stubborn fact that is unfortunately consistent with all current evidence for dark sectors. However, in spite of these sobering facts, the importance of dark sectors simply demands a deeper understanding of gravitational probes, even when potential signals are decades away.    Table 3: Coefficients in Eq. (3.12) for G 2 and F 4 best-fit curves.