Testing the dark SU(N) Yang-Mills theory Confined Landscape: From the Lattice to Gravitational Waves

We pave the way for future gravitational-wave detection experiments, such as the Big Bang Observer and DECIGO, to constrain dark sectors made of SU(N) Yang-Mills confined theories. We go beyond the state-of-the-art by combining first principle lattice results and effective field theory approaches to infer essential information about the non-perturbative dark deconfinement phase transition driving the generation of gravitational-waves in the early universe, such as the order, duration and energy budget of the phase transition which are essential in establishing the strength of the resulting gravitational-wave signal.


I. INTRODUCTION
The null search results for dark matter (DM) via direct detection and colliders suggest that it is likely that DM resides in a hidden sector which couples weakly to the Standard Model (SM) [1][2][3][4][5][6][7]. Yet, very little is known about the dark side of the Universe and it is therefore highly desirable to be able to test the immense landscape of available dark/hidden sectors. Here we concentrate on the well-motivated scenario that the dark side features composite sectors made by non-abelian Yang-Mills theories which are mainly gravitationally coupled. These theories are physically motivated because the dynamics of the dark sector very naturally mimics the SM QCD featuring strong interactions. Furthermore, these theories are well-behaved at short distance denoted as asymptotically freedom [8,9], meaning that the theories are, per se, ultraviolet complete before coupling to gravity, and they do not introduce new types of hierarchies beyond the SM one. The latter means that the theories are stable against quantum corrections. This dynamics has been widely implemented to DM models [10][11][12][13][14][15][16][17][18].
These type of theories are unfortunately inaccessible to current colliders or direct searches, limiting our ability to test them and therefore pin down the model underlying the dark sector. Here we propose to investigate the dynamics of this hidden sector via the detection of gravitational waves (GWs) in a model-independent fashion. We make as few assumptions as possible on the specific extensions of the SM, and assume the minimal interaction of gravity between the SM and new strongly-coupled sectors. Intriguingly, GWs provide a unique window for detecting the dark deconfinement phase transition 1 . * Email: huang@cp3.sdu.dk; ORCID: 0000-0001-7939-3246 † Email: m.reichert@sussex.ac.uk; ORCID: 0000-0003-0736-5726 ‡ Email: sannino@cp3.sdu.dk; ORCID: 0000-0003-2361-5326 § Email: wang@cp3.sdu.dk; ORCID: 0000-0002-5602-6897 1 Confinement occurs at sufficiently low temperatures when gluons To set the stage, we assume that the dark landscape is constituted by n copies of SU (N ) Yang-Mills confined theories for a given confinement scale. To determine the relevant physical information in strongly-coupled theories we often require lattice simulations and/or effective field theory approaches. We adopt state-of-the-art results of lattice simulations [19,20] combined with welldefined effective approaches [21][22][23][24] to precisely pin down the nonperturbative physics involved in the (dark) deconfinement phase transition as functions of the temperature and number of dark colours. For the effective description we marry the Polyakov loop action [21][22][23][24] with lattice simulations. The N = 3 case was extensively investigated in the literature [25][26][27] while here we go beyond the state-of-the-art by incorporating the lattice results [20] at the effective action level for N = 4, 5, 6, and 8. These cases are phenomenologically motivated as they arise in a variety of Grand Unified Theories and composite models. Our work thus covers a wide range of theories. We carefully analyse the dark phase transition for the first few numbers of colours and then generalise it to arbitrarily large numbers. This allows us to acquire an unprecedented eagle view on the dynamics involved in phase transitions of dark composite models by generalising the results to arbitrary numbers of colours. In our work, we go beyond the state-of-the-art by connecting different research fields from (astro-)particle physics over first principle numerical simulations to GW astronomy. It underscores the necessity of orchestrated plans and efforts to unravel the enigma on the nature of DM.
We investigate the GW generation triggered by the dark confinement phase transition discovering, for our generic setup, that: (i) The strength parameter α, reform composite states known as glue balls. At high temperature the theory deconfines the gluons by melting the composite states. Therefore we indicate by Dark confinement-deconfinement to the expected phase transition as function of the temperature occurring during the evolution of the universe and taking place in the hidden sector. lated to the energy budget of the phase transition, takes values around α ≈ 1/3, while the parameter β, that measures the inverse duration of the phase transition, assumes values of the order of 10 4 − 10 5 in units of the Hubble time. (ii) The GW signal emerging from sound waves dominates over the bubble collision and turbulence due to the impact of the friction term [28,29] related to the bubble-wall velocity. (iii) The strength of the induced GW signal is nearly independent of the number of colours for N ≥ 6. That is because the strength depends on the jump in the entropy across the deconfinement phase transition per degree of freedom rather than on the overall jump in entropy, which is inevitably proportional to N 2 . The strength of the GW signal culminates at the case of SU (6) and then gradually decreases with an increasing number of dark colours. Also the peak frequency increases with the number of colours.
The bubble profile and the nucleation rate can be also directly computed in the thin-wall approximation, which allows us to have an independent check of our results from the effective Polyakov loop model. The analysis procedure is neatly summarised by the flow chart in Fig. 1.
In the thin-wall approximation, the nucleation rate is directly obtained from the latent heat and the surface tension, which have been computed with lattice simulations [19,20]. The results from both methods are in qualitative agreement, reinforcing the validity and consistency of our work.
This work constitutes a stepping stone towards embarking in a careful analysis of dark sectors featuring both dark gluons and quarks 2 . In this case, the relevant phase transitions include the dark deconfinement and the dark chiral phase transition. We can take into account these transitions by extending the current work to properly marrying lattice data with the appropriate effective actions introduced first in [53,54].

A. Polyakov Loop
In this work, we consider SU (N ) Yang-Mills theory at finite temperature T . The dynamics is purely gluonic and no fermions are involved. Following 't Hooft [55,56], in any SU (N ) gauge theory, a global Z N symmetry, called the central symmetry, naturally emerges from the associated local gauge symmetry. It is possible to construct a number of gauge invariant operators charged under this global Z N symmetry. Among them, the most notable one is the Polyakov loop, where is the thermal Wilson line, P denotes the path ordering, g is the SU (N ) gauge coupling, and A 0 is the vector potential in the time direction. The symbols x and τ denote the three spatial dimensions and the Euclidean time, respectively. The Polyakov loop can be transformed under the Z N symmetry The phase φ shows the discrete symmetry Z N . From (3), it is clear that is real when N = 2 and otherwise is complex. An important feature of the Polyakov loop is that its expectation value vanishes below the critical temperature T c , i.e. T <Tc = 0, while it possesses a finite expectation value above the critical temperature, i.e.
T >Tc > 0 . In fact, at very high temperature, the allowed vacua exhibit a N -fold degeneracy and we have where 0 is defined to be real and 0 → 1 as T → ∞. Thus, the Polyakov loop is a suitable order parameter in the finite temperature phase transition of the SU (N ) gauge theory.

B. Effective Potential of the Polyakov Loop Model
The Polyakov Loop Model (PLM) was proposed by Pisarski in [21,22] as an effective field theory to describe the confinement-deconfinement phase transition of the SU (N ) gauge theory. The Polyakov loop (1) plays the role of an order parameter. The simplest effective potential preserving the Z N symmetry is given by Gravitational-wave spectrum Figure 1. Schematic summary of our work flow from lattice simulations to gravitational waves. The lattice results for the latent heat and the surface tension are taken from [19], while the ones for the energy-momentum tensor and the pressure are taken from [20].
We have chosen the coefficients b 3 and b 4 to be temperature independent following the treatment in [25,27], which studied the SU (3) case, and also neglected higher orders in | | in (5a). Note that there is no a 4 term in the parameterize of b 2 (T ) in (5b) in [25,27] while we find it can improve the chi-square fitting discussed below. The a 2 term in (5b) has the physics meaning of the "fuzzy bag" term in the "fuzzy bag" model 3 proposed in [57] as a generalization of the famous MIT bag model [58]. On the other hand the a 4 term actually captures the low temperature information and is equivalent to the P[ ] contribution 4 in the model proposed in [24]. Note that the above PLM potential (5a) is the minimal case since we have only considered the Polyakov loop with charge one. For higher charge cases, say charge two cases, the effective potential will be similar to a multi-scalar fields Higgs portal model (see e.g. [23]). However, in the special case where the higher charge Polyakov loop is heavy and can be integrated out, the low energy effective field theory shares a similar form as the current setting of the PLM potential in (5).
With the set-up of the PLM effective potential using (5), we study the SU (N ) Yang-Mills theory with N = 3, 4, 5, 6, and 8. By choosing these numbers of colours, we can take the advantage of the existing lattice data [20]. In the following, we explicitly list the PLM potential corresponding to the number of colours. Extra terms 3 In the Fuzzy Bag model the pressure as a function of temperature p(T ) is written as p(T ) = fpertT 4 −B fuzzy T 2 −B MIT where fpert denotes the perturbative contributions, B fuzzy is the "fuzzy bag" term and B MIT is the term associated with the usual MIT bag model. 4 In [24], it was proposed that the total effective potential are added to some of the cases such that the potential is bounded from below (the SU (4) case) or the fit to the data is decent -chi-squared per degree of freedom is around or below one (the SU (6), SU (8) cases).
For the SU (3) and SU (5) cases, the PLM potential is exactly given by the formula (5a) with N = 3, 5. For the SU (3) case, there is also an alternative logarithmic parameterization, see e.g. [27,59], given by The coefficients inside the logarithm are determined by the Haar measure for which the explicit form for SU (N ) with N > 3 is unknown. Thus, we do not have a logarithmic parameterization for N > 3. For the SU (4) case, the PLM potential is more subtle since the b 3 term is given by 4 + * 4 and thus of the same order as the b 4 term. As consequence, their effects are indistinguishable for real values of and we have to introduce an | | 6 term to properly parameterize the lattice results [19,20]. Thus, the PLM potential for SU (4) is given by where b 2 (T ) is given by (5b). For the SU (6) and SU (8) cases, the PLM potentials are parametrized in the same way as  where b 2 (T ) is again given by (5b). We emphasise that we could include higher-order terms such as | | 8 in the potentials for SU (3) and SU (4) ((5a) and (8)) but they would not improve the fit on the Lattice data and the respective coefficient b 8 would be strongly suppressed.
C. Fitting the PLM potential to lattice data With the explicit PLM effective potential for different colours, we are now able to determine the parameters b i by fitting the potential to the lattice results in [20]. The thermodynamical observables measured on the lattice are the pressure p, the energy density e and the trace of the energy-momentum tensor θ 5 and the entropy density s. The lattice simulations compute the difference between the finite temperature expectation value and the zero temperature one. The energy density e and the entropy density s can be written as linear combinations of the pressure p and the trace of energy-momentum tensor θ Thus we only use the lattice data of θ and p from [20] to determine the coefficients of a i and b i in the above PLM potential setting. We only have access to the statistical 5 Our θ is defined the same as ∆ in paper [20].
uncertainties and therefore we inflated them by a factor of two to mimic the effect of the systematic uncertainties.
During the chi-square (χ 2 ) analysis, we impose the Stefan-Boltzmann (SB) limit: | | → 1 for T → ∞ and p/T 4 | T →∞ → 1.21 · (N 2 − 1) · π 2 /45 [20], which provides two constraints on the parameters of the polynomial parameterizations but only one constraint for the logarithmic case. The above guarantees that the pressure approaches the ideal gas law at infinite temperature. Additionally, the parameters of a i and b i need to fulfil the constraint that T c , which is a priori only a parameter in (5), is indeed the critical temperature, i.e., the temperature at which two minima are degenerate.
We employ the Python package emcee [60], which is based on Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler, to find favourable regions of the parameter space. In Fig. 2, with the help of the analysis tool for MCMC samples, GetDist [61], we display the best-fit regions of the SU (3) log case where a 0 is fixed by the SB limit with a 0 = 4.26. In Fig. 3, we demonstrate how well the best-fit point, with a reduced χ 2 = 0.70, can fit both, p and θ. We present the best-fit values of the potential parameters for the colours N = 3, 4, 5, 6, and 8 in Tab. 1.

III. FIRST-ORDER PHASE TRANSITION AND GRAVITATIONAL WAVES
In this section, we discuss the order of the confinementdeconfinement phase transition and the resulting GW signal. We start with a brief review of the bubble nucleation process and the computation of the GW parameters α (strength parameter), β (inverse duration time), and v w (bubble-wall velocity). Then, we discuss the analytic results obtained from the thin-wall approximation and compare with the PLM fitting results. Remarkably, the analytic results from the thin-wall approximations show interesting patterns that are consistent with those of the fitting to the lattice results. For reviews on GWs from first-order phase transitions see, e.g., [62][63][64][65][66][67].

SU(3) log
Best-fit Data Figure 3. Comparison between the data (red) and best-fit curves (blue) for the pressure (left panel) and trace of the energy momentum tensor (right panel). Note that we normalize both p and θ with respect to the SB limit as done in [20].

A. Bubble nucleation
In this section, we briefly review the generic picture of bubble nucleation processes where some subtleties related to our models are emphasized.
The conventional picture of a first-order phase transition is that, as the universe cools down, a second minimum with a non-zero vacuum expectation value (broken phase) develops at a critical temperature. This triggers the tunnelling from the false vacuum (unbroken phase) to the stable vacuum (broken phase) below the critical temperature. In our model, this picture is reversedin a sense, as the universe cools down, the tunnelling occurs from the broken phase (deconfinement phase) to the unbroken phase (confinement phase). The underlying reason behind this reversed phenomenon is that the discrete symmetry Z N is broken in the deconfinement phase at high temperature while it is preserved at the confinement phase at low temperature.
The tunnelling rate due to thermal fluctuations per unit volume as a function of the temperature from the metastable vacuum to the stable one is suppressed by the three-dimensional Euclidean action S 3 (T ) [68][69][70][71] and we have The three-dimensional Euclidean action reads where ρ is a scalar field with the effective potential V eff . The scalar field ρ has mass dimension one, [ρ] = 1, in contrast to the Polyakov loop , which is dimensionless. Furthermore, V (ρ, T ) has mass dimension four. After rewriting the scalar field as ρ = T and converting the radius into a dimensionless quantity r = r T , the action becomes which has the same form as (12).
Keep in mind that r in the bubble-profile solution is not the physical bubble radius but the product of bubble radius and the temperature. The bubble profile (instanton solution) is obtained by solving the equation of motion of the action in (13) with the associated boundary conditions To attain the solutions, we used the method of overshooting/undershooting and employ the Python package CosmoTransitions [72]. A sample plot is shown in Fig. 4 where we also indicate the thickness of the bubble wall, which we will use later. We substitute the bubble profile (r , T ) into the three-dimensional Euclidean action (13) and, after integrating over r , S 3 depends only on T .

B. Inverse duration time of the phase transition
An important parameter for the computation of the GW signal from a first-order phase transition is the inverse duration time β. For sufficiently fast phase transitions, the decay rate can be approximated by where t * is the characteristic time scale for the production of GWs. The inverse duration time then follows as  The dimensionless version is obtained by dividing with the Hubble parameter H where we used that dT /dt = −H(T )T . Note that in the above analysis, we have implicitly assumed that the temperature in the strongly coupled hidden sector (denoted as T d ) and the temperature in the visible sector (T v ) are the same (i.e. ξ ≡ T d /T v = 1). In general, these two temperatures can be different. In this case, the inverse duration is given bỹ with Hubble parameter H(T v , T d ) given Here, g * ,d and g * ,v are the effective number of relativistic degrees of freedom in the hidden and visible sector, respectively. The phase transition temperature T * is often taken as the nucleation temperature T n , which is defined as the temperature at which the rate of bubble nucleation per Hubble volume and time is approximately one, i.e. Γ/H 4 ∼ O(1). A more accurate definition is to use the percolation temperature T p , which is defined as the temperature at which the probability to have the false vacuum is about 0.7. For very fast phase transitions, as in our case, the nucleation and percolation temperature are almost identical T p T n . However, even a small change in the temperature leads to an exponential change on the vacuum decay rate Γ, see (16), and consequently we use the percolation temperature throughout this work. We write the false-vacuum probability as [73,74] with the weight function [75] The percolation temperature is defined by I(T p ) = 0.34, corresponding to P (T p ) = 0.7 [76]. Using T * = T p in (18) yields the dimensionless inverse duration time.

C. Strength Parameter α
Many analysis have used the MIT bag model to obtain the strength parameter α of the phase transition. As already mentioned in Sec. II, the bag model is not sufficient to precisely describe the confinementdeconfinement phase transition, and the Fuzzy Bag model [57] is required. In the bag model, the bag constant is used to parameterize the strength of the phase transition The bag constant parameterizes the jump in both the pressure and energy density across the phase boundary For work that defines the strength parameter beyond the bag model, see, e.g., [77,78]. Here, we define the strength parameter α from the trace of the energy-momentum tensor θ α = 1 3 where ∆X = X (+) − X (−) for X = (θ, e, p) and (+) denotes the meta-stable phase (outside of the bubble) while (−) denotes the stable phase (inside of the bubble). The enthalpy density w ± is defined by which encodes the information of the number of relativistic degrees of freedom (d.o.f). It is intuitive to use the trace of the energy momentum tensor θ to quantify the strength of the phase transition α. In the limiting case when θ = 0, the system possesses conformal symmetry and there is a smooth second-order phase transition occurring. θ is a quantity to measure the deviation from the conformal symmetry and thus also measures the deviation from the second-order phase transition. The larger θ is, the further away from the conformal symmetry and second-order phase transition and thus the stronger the first-order phase transition is.
In the case of the confinement-deconfinement phase transition, α can be directly computed from the lattice results of ∆e and ∆p of [20], see Sec.II C. In our language of the PLM potential, we set the pressure and energy in the symmetry-broken phase to zero and measure energy and pressure relative to this phase, e − ∼ p − ∼ 0. Thus, α can be rewritten in terms of the V where we have used as well as (26). Furthermore, at the percolation temperature (which is close to T c ), we always have e + p + [20], leading to α ≈ 1/3. Note that our definition of α only depends on the degrees of freedom in the hidden sector. In other works [79][80][81], two different α have been introduced where one of them is denoted by α d , identical to the one defined in (27), and the other is α tot = ∆θ/3w + tot , in which w + tot is the total enthalpy including the visible and dark relativistic degrees of freedom. The parameter α d is then used to compute the wall velocity and efficiency factors, while α tot is used in the GW formula for the peak amplitude. To avoid the confusion, we only define a single α but take into account the dilution effect on the GW signals due to the presence of other degrees of freedom, see Sec. III G for more details.

D. Bubble-wall velocity
The bubble-wall velocity v w is another important parameter, which determines the strength of the GW signal. The bubble-wall velocity requires a detailed analysis of the forces that act on the bubble wall. The forces can be divided into two parts. The first force arises from the difference of the vacuum potential (pressure) between the confinement and deconfinement phases. This force accelerates the wall and causes the bubble to expand. The second force is the friction on the wall, which can be further divided into two kinds as discussed below [28,29,82,83]. For more recent work which calculates the bubble wall velocity beyond the leading-log approximation see e.g. [84].
Direct Mass Change The first kind of friction is due to the direct mass change of a particle when passing through the interface between the two phases (first proposed in [28]). The mass change results in a momentum change along the bubble moving direction, leading to a friction force on the bubble wall where F 1 , A, p f denotes respectively the friction force, the surface area and the pressure on the wall associated with the friction force. ∆m 2 a represent the mass square difference between the stable phase and metastable phase for the particle species a. f a (p, in) is the distribution function for the incoming particles i.e. in the deconfinement phase. In the phase transition from deconfinement to confinement, the gluons will confine to glueballs and become massive. A detailed estimate of this friction force relies on the estimate of the glueball mass of different numbers of colours. More importantly, (29), rigorously speaking, is derived from 1 → 1 process (one incoming particle and an outgoing one) whereas the formation of glueballs from gluons is more complicatedprocesses such as 2 → 1 and 3 → 1 may take place. In this case, a generalization of (29) would be required. A more detailed study on the glueball formation is beyond the scope of this work and will be pursued in the future. Nevertheless, there certainly exists friction in light of the direct mass change from gluons to glueballs.
Particle Splitting The second kind of friction (first discussed in [29]) is through the particle splitting (transition radiation process) where an incoming particle changes its momentum (along the bubble wall direction) through emitting another particle that exerts a friction force on the bubble wall. It was shown in [29] that in a large class of transition radiation processes such as S → V T S, F → V T F, V → V T V , (where S, V, F, T denote respectively scalar, vector, fermion and transverse modes), the friction is given by where γ is the Lorentz factor since the friction scales with the incoming particle density and g is the coupling of the involved interaction. ∆m represents the mass change of the particle at the interface, which implies that the friction resulting from the particle splitting process will always be associated with the above-mentioned friction of the direct mass change at the interface. In a weakly coupled theory, this second kind of friction is sub-leading compared with the previous one. However, in our strongly coupled system, this second friction can be equally important.
Wall Velocity In summary, we can write the total pressure on the bubble wall as (see also [85]) where ∆V denotes the pressure due to the difference of the vacuum potential between the confinement and deconfinement phases, which accelerates the wall. Assuming ∆V > p f1 , we can obtain the equilibrium γ (denoted as γ eq below) when the net force on the bubble wall becomes zero and the wall velocity (also γ) ceases to grow Using (32), we can obtain the terminal wall velocity. Relation to Energy Budget In the end, it is important to consider the fraction E wall /E V where E wall corresponds to the wall energy at the terminal velocity and E V is the total vacuum energy. This fraction describes how much of the total vacuum energy goes into accelerating the bubble wall. This part will eventually contribute to the GW signal from the bubble collisions. The remaining part of the energy 1 − E wall EV goes into the surrounding plasma and contributes to the generation of GWs via sound waves and turbulence, where typically the sound wave contribution dominates. In the case of the deconfinement phase transition, both friction terms p f1 and p f2 are non-perturbative due to the strong gauge coupling. Thus the main part of the energy will be stored in the plasma surrounding the bubble wall and, in consequence, we can focus on the GW production from sound waves. Due to the non-perturbative nature of the friction terms, it is highly challenging to determine them quantitatively. Instead, we treat the terminal bubble wall velocity as an input parameter and investigate the impact of different values.

E. Thin-Wall Approximation
The advantage of the thin-wall approximation is that we can calculate analytically the decay rate of the false vacuum in terms of the latent heat and the surface tension. The latter are provided from lattice results as a function of the number of colours N [19,20]. The thinwall formula for the Euclidean action is shown in [71,86] and we briefly review it below. The three-dimensional Euclidean action is written as where p de and p co denote respectively the pressure in the deconfinement and confinement phase, σ is the surface tension of the nucleation bubble, and r c is the critical radius of the nucleation bubble defined by On the other hand, the difference of the pressure between the deconfinement and confinement phase is also linked to the latent heat L via Finally, by using (34) and (35), the three-dimensional Euclidean action (33) can be written as a function of latent heat L and surface tension σ where the latent heat from the lattice results [19] is The lattice error on the 1/N 2 coefficient bares the largest uncertainty and will eventually contribute the most to the uncertainty of the GW parameters, as we will see later. The surface tension on the other hand can be either proportional to N or N 2 due to indecisive lattice results. Intuitively, one may expect that the strength of the phase transition increases with N . The strength however depends on both L and σ as shown in (36), where L/N 2 is related to the latent heat per d.o.f -(37) becomes independent of N for N 1 -and σ ∝ N or N 2 . As a result, the strength of the GW only grows with N if σ ∝ N 2 . As we shall see in the next section, Sec. III F, the PLM fitting prefers the σ ∝ N case. Nonetheless, we discuss both scaling behaviours of the surface tension in the following.
σ proportional to N : In this case, the lattice fitting function of the surface tension is [19] By implementing (38) and (37) to the Euclidean action (36), we obtain This function has an interesting behaviour: for fixed temperature factor T 3 c (Tc−T ) 2 , S 3 has a maximum at N ∼ 11. In the large-N limit, the Euclidean action behaves as S 3 ∼ 1/N , which implies that the effective PLM potential scales as V eff ( ) ∼ N 2 , see (13) and (14).
From (39) together with (22), we determine the percolation temperature. As a rule of thumb, the phase transition occurs around S 3 /T ∼ 150 for T c in the GeV range. For other T c , this criterion changes with logarithmically with T c . We observe that the difference between percolation temperature T p and critical temperature T c denoted as δT = |T p − T c | starts to increase from N = 3 until it reaches a maximum at N = 11 and then gradually decreases. As mentioned above, one might naively expect that the strength of the phase transition increases with N . Thus δT , which relates to the strength of the phase transition, should also increase with N . However, this pattern only corresponds to σ ∝ N 2 case. It should be noted that δT at N = 11 is around 20 times bigger than δT at N = 3. A bigger value of the temperature difference δT implies a longer duration and a stronger 20 40 60 80 100 10 3 10 4 10 5 first-order phase transition 6 and a stronger GW signal. Thus, we expect an increasing GW signal from N = 3 to N = 8 using the aforementioned method of the PLM effective potential if surface tension σ is proportional to N .
We can go one step further to derive the dimensionless inverse durationβ. Using (39), (22) and (17), we compute values of the dimensionless inverse durationβ which are shown in Fig. 5 with solid blue line and dashed green line. Theβ shares exactly the inverse pattern as the more intuitive parameter δT discussed above i.e.β first decreases to around N = 11 and then increases with N . Since the gravitational wave peak amplitude is inversely proportional toβ, it is consistent with the above discussion using δT .
σ proportional to N 2 : In this case, the lattice fitting function of surface tension is [19] while the latent heat is still following (37). By substituting (40) and (37) into the Euclidean action (36), we obtain This function has the following behaviour: For fixed temperature factor T 3 c (Tc−T ) 2 , S 3 keeps increasing with the number of colours N . In the large-N limit, the Euclidean action behaves as S 3 ∼ N 2 , which implies V eff ( ) ∼ 1/N 4 as can be seen from (13) and (14). 6 Similar features sometimes are shown in the case of supercooling with a strong first-order phase transition and a longer duration [87][88][89][90][91][92][93].
The pattern in this scenario is different from the previous case of σ ∝ N . By using (41) and (22), we find that δT = |T p − T c | monotonically increases with N . Thus, the larger N the bigger δT , resulting in a stronger firstorder phase transition and GW signals. Nonetheless, for both cases, σ ∼ N and σ ∼ N 2 , the thin-wall approximation gives consistent results for small N , i.e., N 7. The ambiguity of the scaling behaviour of the surface tension can only be resolved in a strict sense by more accurate lattice results at large N . However, as we will show in the next section, the PLM fitting procedure seems so be only consistent with σ ∼ N .
The dimensionless inverse duration timeβ is again calculated by (41), (22) and (17), and the results are summarized in Fig. 5 with dashed red and orange lines. Thẽ β shares exactly the inverse pattern as the more intuitive parameter δT discussed above i.e.β keeps decreasing with N . Since the gravitational wave peak amplitude is inversely proportional toβ, it is consistent with the above discussion using δT .

F. Thin-Wall Approximation vs Fitting of PLM Potential
In this section, we compare the results of the dimensionless inverse durationβ between the two methods: the thin-wall approximation and the PLM potential fitting. We also comment on the wall thickness, which depends on the number of dark colours, and relates to the validity of the thin-wall approximation.
The comparison ofβ is displayed in Fig. 6. The result of the fitted PLM potential, marked by green triangles, shows the following pattern: apart from N = 3,β first decreases with N and then increases after N = 6. Interestingly, this pattern qualitatively agrees with that of the thin-wall approximation with σ ∝ N , although there the turning point is located around N ∼ 11. For the thin-wall approximation, we include error bands due to the lattice error displayed in (37), (38), and (40). The main error stems from the N 2 coefficient in the latent heat (37). Note that we do not display the statistical uncertainties ofβ associated with preferred regions of the χ 2 fits in Fig. 6 since they are small compared to the dot size (of the order of 10%). However, the systematic uncertainties of the lattice data discussed in [20] have not been included in our fitting procedure. They may give rise to larger uncertainties onβ. In the later computation of GW signals, we try to include those uncertainties by enhancing the statistical error by a generous factor of five.
On the right panel of Fig. 6, for points of N > 8 without available lattice data, we assume that the energy and the pressure normalised to the SB limit become independent on N in the large-N limit. These assumptions are supported by lattice data for the pressure and energy [20]. This entails that p ∼ e ∼ N 2 in the large-N limit and thus the effective PLM potential scales as In this case, the potential for N > 8 can be obtained by a simple rescaling that of N = 8, i.e., V eff (N ) = N 2 V eff (N = 8)/8 2 . As discussed in the previous section, the scaling of the potential with N 2 corresponds to the scenario of σ ∝ N in the thin-wall approximation.
We observe thatβ from rescaled PLM potentials has a power-law behaviour as a function of N (linear function in the log-log plot in Fig.6). Intriguingly, the N = 6 data point, which is obtained by using the direct lattice results rather than through rescaling, is in good agreement with the rescaling results. This seems to indicate that the information encoded in the lattice results for N = 6 and N = 8 favours the scenario of σ ∝ N rather than σ ∝ N 2 . Note that the blue curve in Fig. 6 also becomes a linear function in the log-log plot for large N but with a slightly smaller slope compared to that of the PLM fitting potentials.
The fact that the PLM fitting favours σ ∝ N over σ ∝ N 2 has a direct impact on the peak amplitude of the GW signal as discussed in the next two sections. It implies that the GW peak amplitude first increases (corresponding to a decreasing peak frequency) from N = 4 to N = 6 which has the lowest frequency and then gradually decreases (while the peak frequency increases) with increasing N . On the other hand, for the case σ ∝ N 2 , which is not favoured by the PLM fitting potential, the signal becomes monotonically stronger with larger N (≥ 4).
Before discussing the GW spectrum, we comment on the wall thickness. The wall thickness ∆r and the bubble radius r * can be directly computed from the instanton solution at the percolation temperature, see Fig. 4. We choose the wall thickness definition that the two wall boundaries are located 10% away from the broken and unbroken Polyakov loop vacuum expectation values 7 . We show the ratio of the wall thickness to the bubble radius in Fig. 7 where the values for cases of N > 8 are obtained via rescaling of the SU (8) potential as discussed above. The wall is relatively thick for cases of N = 3, 4, 5 and becomes thinner for N ≥ 6. It approaches a constant in the large-N limit. This is also consistent with what we have observed in Fig. 6; the results of PLM fitting potential are in agreement with or close to those of the thin-wall approximation for N = 6, 8 while noticeable deviations between two methods are present for N = 3, 4, 5.

G. Gravitational-wave spectrum
We briefly review the computation of the GW spectrum from the parameters α,β, and v w . In general, there are three contibutions to the GW spectrum: collisions of bubble walls [94][95][96][97][98][99][100][101][102][103], sound waves in the plasma after bubble collision [104][105][106][107][108] and magnetohydrodynamic turbulence in the plasma [109][110][111][112][113][114][115][116]. As discussed in Sec. III D, in the case of the deconfinement phase transition, the contributions from sound waves are dominating and thus we focus on this contribution. Following [65,66,117], the GW spectrum from sound waves is given by with the peak frequency f peak 1.9 · 10 −5 Hz g * 100 and the peak amplitude Here, h = H/(100km/s/Mpc) is the dimensionless Hubble parameter and g * is the effective number of relativistic degrees of freedom including the the SM degrees of freedom g * ,SM = 106.75 and the dark sector ones g * ,SU (N ) = 2(N 2 − 1) · n, where n is the number of copies of dark sectors. The factor Ω 2 SU (N ) accounts for the dilution of the GWs by the visible matter which does not participate in the phase transition. The factor reads Ω SU (N ) = ρ rad,SU (N ) ρ rad,tot = g * ,SU (N ) g * ,SU (N ) + g * ,SM , where we again assumed that both sectors have the same temperature. In other works [79][80][81], two different strength parameters α tot and α d were introduced as discussed in Sec. III C. In this case, the peak amplitude can be expressed in terms of these two quantities without involving the dilution factor: Notice that the efficiency factor and the wall velocity depend on α d only.
In the last sections, we have detailed the computation of the parameters α andβ, and argued that we use the wall velocity as a free input parameter. The last important ingredient is the efficiency factor κ, which describes the fraction of energy that is used to produce GWs. The efficiency factor is made up of the efficiency factor κ v [101] and an additional suppression factor due to the length of the sound-wave period [85,90,118]. In total, the efficiency factor κ sw is given by Note that we measure τ sw in units of the Hubble time and thus it is dimensionless. We first discuss the contribution from κ v where we use the results from [101]. This efficiency factor depends on the wall velocity and the strength parameter. While it increases for larger α, it typically has a maximum when the wall velocity assumes the Chapman-Jouguet detonation velocity v J , which is given by For the deconfinement phase transition where α ≈ 1/3, the detonation velocity takes the value v J = 0.866. Due to the complicated dependence of the efficiency factor on the bubble-wall velocity, we simply display it for the wall velocities used here. In the next section, we test the impact of the wall velocity on the GW spectrum employing the values v w = (1, v J , 0.2). For v w = 1, κ v is given by which implies κ v ≈ 0.3 for α ≈ 1/3. At the Chapman-Jouguet detonation velocity v J , the efficiency factor reads and for α ≈ 1/3 we have κ v ≈ 0.45. As expected, this value is larger that for v w = 1. For smaller wall velocities v w < c s , where c s = 1/ √ 3 is the speed of sound, the efficiency factor decreases rapidly and the generation of GW from sound waves is suppressed [119]. For example, for v w = 0.2, we have which implies κ v ≈ 0.19 for α ≈ 1/3. The second contribution to the efficiency factor κ sw stems from a suppression due to the length of the soundwave period τ sw , see (47). In [85,90], the length of the sound-wave period was given by The recent work [118] has analysed the length of the sound-wave period in an expanding universe and there the suppression was given by  They depend on the root-mean-square fluid velocity [85,107], which is given bȳ In our case,β 1 and thus (52) and (53) lead to almost identical suppression factors. For the subsequent analysis in the next section, we use (53).
An important quantity that determines the detectability of a GW signal at a given detector is the signal-tonoise-ratio (SNR) [120,121] given by Here, h 2 Ω GW is the GW spectrum given by (42) the sensitivity curve of the detector, and T the observation time, for which we assume T = 3 years. We compute the SNR of the GW signals for the next generation of GW observatories which are LISA [30][31][32], BBO [33][34][35][36][37], DECIGO [37][38][39][40], ET [41][42][43][44], and CE [45,46] 8 . The sensitivity curves of these detectors are nicely summarised and provided in [123]. It is in general a difficult question, from which SRN onwards a GW signal will be detectable, a typical estimate being SNR > 1 − 10. This issue is also linked to how well the astrophysical foreground such as gravitational radiation from inspiralling compact binaries is understood and can be subtracted from the signal, see, e.g., [124,125]. Here, we make the optimistic assumption that a signal with SNR > 1 is detectable.

IV. RESULTS
Here we discuss our main results on testing the dark confinement landscape using the next generation of GW observatories which are LISA, BBO, DECIGO, ET, and CE. We focus on the results obtained through the fitting of the effective PLM potential. The results from the thin-wall approximation are in qualitative agreement with those of the PLM potential fitting in particular in the case of the surface tension proportional to N as discussed above 9 . The values of the dimensionless inverse durationβ for the different values of N are displayed in on the wall velocity is only mild. The strength parameter takes values α ≈ 1/3, see Sec. III C. From the PLM fitting, we obtain the statistical uncertainty on α and β, which we inflate by a generous extra factor of five to account for hidden systematic errors. We display this uncertainty with a band on the GW spectrum. For the GW experiments, we display in all figures the power-law integrated sensitivity curves, see e.g. [123,126]. The powerlaw integrated sensitivity curves can differ from standard sensitivity curves by several orders of magnitude. For the detectability of the GWs, we therefore refer strictly to the SNR, which is displayed in Fig. 12. As input parameters for our computation, we have the wall velocity v w , the confinement temperature T c , the number of dark colours N , and the number of dark SU (N ) copies n. Let us start by discussing Fig.8 where we show how different bubble wall velocities affect the GW spectrum using SU (6) with T c = 1 GeV as a testbed example. There are two competing effects at play here. The first is that the efficiency factor is maximal at the Chapman-Jouguet detonation velocity v J , see (50). The second is that the amplitude itself is proportional to the wall velocity, see (44). This means that at fixedβ and apart from the v J case, higher wall velocities tend to provide higher peak amplitudes and lower peak frequencies of GWs.
In Fig. 9, we compare the SU (6) GW spectrum with and without the suppression factor given in (53). The suppression factor leads in our case typically to a suppression of 10 3 − 10 4 or so. The suppression is significant for weak phase transitions and small for strong phase transitions. Excitingly the GW signal with suppression may still be detectable by BBO and DECIGO. Should the suppression of the GW signal due to the length of the sound-wave period be smaller than expected, then even LISA may be able to detect a signal from a deconfinement phase transition.
The dependence of the GW spectrum on the number of dark colours is shown in Fig. 10 for the values of N = 3, 4, 5, 6, 8. All spectra are plotted with the bubblewall velocity set to the Chapman-Jouguet detonation velocity and with T c = 1 GeV. To make the figure concise, we do not show the error bands of the GW spectrum but only display the GW spectrum using the central values instead. From the plot, we learn that the peak amplitude of the induced GW signal is nearly independent of the number of colours for N ≥ 6. That is due to the fact the strength depends on the jump in the entropy across the deconfinement phase transition per d.o.f. rather than the overall jump in entropy, which is inevitably proportional to N 2 . This argument in principle also applies to small numbers of colors N = 3, 4, 5, which is reflected in the overall mild dependence ofβ on N seen in Fig. 6. However, for small N , the GW signal is more strongly diluted by the d.o.f. of the SM, see (45). The strength of the GW signal first increases (corresponding to decreasing peak frequency) starting from N = 4 until reaching its maximal amplitude for N = 6 (lowest frequency) and then the GW amplitude slowly decreases (with increasing frequency) when increasing N . This is in agreement with our expectations presented in Sec. III F stemming from the dependence of the inverse duration timeβ with respect to the dark colours N shown in Fig. 6.
In Fig. 11, we present how different confinement scales (including T c = 1 GeV, 1 TeV, and 1 PeV) affect the GW spectrum. As expected, a higher confinement scale leads to a higher GW peak frequency. On the other hand, the shape of the GW spectrum is independent of the confinement scale and also the peak amplitude depends only mildly on the confinement scale. Interestingly, BBO and DECIGO will test confinement phase transitions in the GeV range.
In Fig. 12, we show the SNR of the phase transition of one dark SU (6) sector for different GW detectors as a function of the confinement temperature T c . Due to the minor dependence for N ≥ 6, as displayed in Fig. 10, it is expected that cases of larger N will fea- ture similar SNRs shown here. Assuming that the signal is detectable for SNR > 1, we find that BBO and DE-CIGO will test theories with a confinement scale within 1 GeV T c 100 GeV. Other GW detectors such as LISA, CE, and ET manage to achieve an SNR of O(10 −3 ) in the GeV/TeV range. This analysis includes the suppression factor due to the short sound-wave period (53).
The GW experiments will test a lager part of the landscape if the suppression factor is smaller than expected.
In Fig. 13, we show how future GW observatories will constrain the dark deconfinement landscape. We make the assumptions that there are n non-interacting copies of SU (N ) gauge theories and all of them undergo the phase transition at the same scale T c . The GW signals of these phase transitions add up linearly. However, due to BBO DECIGO 10 1 10 2 1 2 3 4 5 T c in GeV n Figure 13. We display the exclusion curves of n dark SU (N ) phase transitions from the experiments BBO, DECIGO as a function of the confinement temperature Tc. We assume an observation time of three years and that the signal is detectable for a signal-to-noise ratio SNR > 1.
the dilution of the GW signal over the non-participating degrees of freedom, the GW signal of each sector is suppressed by a factor of approximately 1/n 2 . Summing up the GW signal of all sectors leads to a total suppression of roughly 1/n. In other words, adding more independent sectors with the same scale of phase transition weakens the experimental constraints. For each dark copies n and T c , we compute the SNR with respect to the future GW detectors BBO and DECIGO, see (55). We assume that the signal is detectable for SNR > 1, and thus those theories will be tested in the future. Excitingly, BBO and DE-CIGO will cover the range of 1 GeV T c 100 GeV of this landscape and can maximally test four dark SU (N ) copies. The results in Fig. 13 again apply to scenarios of N ≥ 6. For N < 6, the GW signal is slightly suppressed, see Fig. 10 -the resulting SNR is smaller but the qualitative features of Fig. 13 still hold.

V. CONCLUSIONS AND OUTLOOK
In this work, we explored the landscape of the strongly coupled dark sectors composed of n-copy SU (N ) Yang-Mills confined theories coupled mainly gravitationally to our world. We employed state-of-the-art lattice results combined with effective field theory (PLM) approaches to investigate the GW signal arising from the dark deconfinement-confinement phase transitions in the early universe. As a comparison, we have also applied the analytic thin-wall approximation, which yields consistent results with those of the PLM approach. Our procedure is summarized in Fig. 1.
We discovered that the strength of the GW signal only depends mildly on the number of colours N for N ≥ 6. We find that the strength parameter of the phase transition is α ≈ 1/3, while the inverse duration time is β = 10 4 − 10 5 in units of the Hubble time. Because of the fact that n copies of a theory with the same confining scale need to share the same universe energy budget, we find that the next generation of gravitational waves observatories are sensitive only to a very small number of copies with confining scales from fraction to hundreds of GeVs as shown in Fig. 12.
We consider our work a natural stepping stone towards the inclusion of matter fields in different representations of the gauge group. In particular, it is interesting to consider the interplay of the dark confinement phase transition and the one stemming from dark chiral symmetry breaking using the methodology of [53,54] and the lattice results from [127]. Another avenue would be to extend the analysis beyond gauge-fermion theories to complete asymptotically free ones that feature composite dynamics including elementary scalars. These theories are well behaved at high energies while still featuring compositeness at low energies. Last but not the least it would be intriguing to investigate for these theories the GW imprint that could come from a dark symmetry, broken at arbitrary high temperatures as shown to exist in [128] following Weinberg's seminal work for UV incomplete theories [129].
Note added: shortly after this paper was submitted, we became aware of the nice complementary work [130]. They investigated the deconfinement phase transition using instead the Matrix Model rather than Polyakov Loop Model yielding compatible results with ours.