Dark confinement-deconfinement phase transition: a roadmap from Polyakov loop models to gravitational waves

We explore the confinement-deconfinement phase transition (PT) of the first order (FO) arising in SU(N) pure Yang-Mills theory, based on Polyakov loop models (PLMs), in light of the induced gravitational wave (GW) spectra. We demonstrate that the PLMs with the Haar measure term, involving models successful in QCD with N = 3, are potentially incompatible with the large N scaling for the thermodynamic quantities and the latent heat at around the criticality of the FOPT reported from the lattice simulations. We then propose a couple of models of polynomial form, which we call the 4-6 PLM (with four- and six-point interactions among the basic PL fields which have center charge 1) and 4-8 PLM (with four- and eight-point interactions), and discuss how such models can naturally arise in the presence of a heavy PL with charge 2. We show that those models give the consistent thermodynamic and large N properties at around the criticality. The predicted GW spectra are shown to have high enough sensitivity to be probed in the future prospected interferometers such as BBO and DECIGO.


I. INTRODUCTION
Pure Yang-Mills (PYM) sector surviving at low energy is of particular interest in the new physics domain, in the context of dark sector of the universe, including dark matter (see, e.g., a review [1] and also [2][3][4]) and also in the string theory [5].However, if such a PYM sector is sufficiently secluded from the visible sector, then what is its detectable signal?
Due to the asymptotic freedom of non-Abelian gauge theory, the PYM theory in the ultraviolet scale lies in the deconfinement phase described by "free" gluons, whereas flowing down to the infrared scale, the interactions between gluons steadily get strong and eventually confine themselves to the colourless glueballs 1 .Therefore, there is a confinement-deconfinement transition at the critical temperature T c .It is well-known that for the SU (N ) PYM theory with N > 2, this phase transition (PT) is first order (FO).Then, one expects the corresponding FO cosmic/thermal PT in the early universe, during which gravitational waves (GWs) are produced.They immediately decouple from the matters and are wandering in the universe today, which may meet the GW detectors in the sky such as the LISA [7,8], Big Bang Observer (BBO) [9][10][11][12][13], DECIGO [13][14][15] and TianQin [16], where the sensitivities of frequencies cover a wide region from 10 −3 mHz to kHz; and the sensitivities can be improved in orders of magnitude over the next three decades.
To quantify the GW prediction, working on phenomenological models of Ginzburg-Landau type modeling the confinement-deconfinement PT is helpful.Since the PT is subject to the nonperturbative aspect of the PYM theory, the output from lattice simulations is often utilized to extract important clues in constructing the phenomenological models.Lattice simulations present the equilibrium thermodynamic properties of the SU (N ) plasma at finite T , in particular for N > 3 [18][19][20][21][22][23].The pressure p(T ) is directly related to the potential of the model in terms of some fundamental degrees of freedom such as the Polyakov loop (PL) [24], thus it is of special importance: • The pressure is found to scale as the ideal Stefan-Boltzmann (SB) limit when T > 4T c , namely p(T ) ∼ (N 2 −1)T 4 .
On the contrary, the pressure tends to vanish below T c .
introduction of a nonminimally-Z N charged PL variable (part A).With these preliminaries at hand, we employ a widely used conventional PLM, called the Haar-type PLM, and show that the model turns out to have incompatibility in a sense of the large N scaling for the latent heat and the other thermodynamic quantities at around the criticality for the confinement-deconfinement PT (part B).Then, we propose a set of PLMs which properly realize the thermodynamic and the large N properties at around the criticality, that is what we call the 4-6 and 4-8 PLMs, with explanation on the derivation for those highly nonminimal coupling terms.The models are shown to yield the good-fitness with the SU (N ) lattice data on the thermodynamic quantities, as well as the latent heat, regarding the large N scaling around the criticality (part C).In Sec.III, we move on to discussion on the generation mechanism of the GW spectra by the FOPT based on the 4-6 and 4-8 PLMs.The dark SU (N ) PYM is assumed to be secluded-dark gluonic plasma (part A).Passing discussion on a couple of the well-known GW sources (part B), we numerically evaluate the GW spectra, in comparison with the prospected GW detector sensitivities (part C).Our conclusion and several issues still left with us are presented in Sec.IV.

II. ROAD TO A PROPER PLM A. PLs and ZN center symmetry
The PL plays a center role in understanding the confinement-deconfinement PT in the PYM theory [24].When a immovable test quark is placed at the position ⃗ x, its color-averaged free energy F q (⃗ x) can be expressed as e −Fq(⃗ x)/T = ⟨tr c L(⃗ x)⟩ with L(⃗ x) the thermal Wilson line at a fixed spatial position ⃗ x extending in the temporal direction: where g denotes the PYM gauge coupling, and the path integration interval is closed due to periodicity of the gauge field A µ : A 4 (⃗ x, x 4 ) = A 4 (⃗ x, x 4 + 1/T ).Here A 4 ≡ A a 4 t a with t a (a = 1, • • • , N 2 − 1) being the generators of SU (N ) in the fundamental representation.Under the local SU (N ) transformation, L transforms as the adjoint representation.Therefore, taking trace in the color space, one obtains a gauge invariant operator, that is the PL: Quantum mechanically, l(⃗ x) is a local operator, while classically it is a complex scalar field with the field value bounded by one: |l(⃗ x)| ≤ 1.Note that l(⃗ x) is dimensionless by construction, and one should make up it by some dimensionful quantity in addressing the cosmic/thermal evolution of the GWs; we will come back to this point in Sec.III (around Eq.( 39)).Due to the compactified temporal direction, gauge fields should respect the periodicity shown before, however, gauge transformations V (⃗ x, x 4 ) can be non-periodic.Instead, to keep the PYM action gauge invariant, they are just required to satisfy the twisted boundary condition (the center twisted gauge transformations): V (⃗ x, β) = z k V (⃗ x, 0), where z k = e i2kπ/N with k = 1, • • • , N belonging to the discrete group Z N ; the transformation satisfying the periodicity corresponds to k = N .In particular, the gauge transformation with respect only to x 4 obviously obeys the twisted boundary condition; it is dubbed the center symmetry [38] since it is simply the power of the elements of z k = z k 1 N ×N ∈ Z N , the center of SU (N ).Under the action of V k (x 4 ), l is not invariant and changes as l → z k l although the Euclidean action is invariant.In this sense, l is charged under the global center Z N .Now, it is seen that the gauge invariant operator l can serve as a good order parameter of the confinement-deconfinement PT: • If ⟨tr c L(⃗ x)⟩ = 0 and Z N is not broken, then F q → ∞, thus quarks are confined.
• If ⟨tr c L(⃗ x)⟩ ̸ = 0 and Z N is broken, then F q = finite, thus quarks are deconfined.
The exact Z N charge of l is not determined from the first principle, and one may specify its charge as 1 and dub it as the fundamental PL, l 1 ≡ 1 N tr(L).Other traced PLs with different Z N charges can be introduced [37,39].Or, we introduce traced PLs in a different representation R under SU (N ), normalized by the dimension d R .But different representations may have the same Z N charge.For instance, in the large N -limit, both PLs in the representations of d R,± = (N 2 ± N )/2 have charge 2. In the large N limit, according to group theory [40], they can be decomposed as where trL 2 means the Wilson line wrapped around the imaginary time forward twice, and it is negligible in the large N limit.In general, PLs with Z N charge p + − p − (module N ) can be expanded as: with l1 ≡ tr(L † )/N , the PL in the anti-fundamental representation.To determine c n,m (N ) is nontrivial [40,41], but irrelevant to our discussions.In constructing PLM incorporating additional PLs, we just treat them as independent degrees of freedom, and require that the couplings among them are constrained by Z N .The higher-Z N charged PLs are assumed to be heavy and can be integrated out, which leads to the polynomial potential that we will use in Sec.II (See Eqs.( 29) and (30), and discussions around there.).
B. The Haar-type PLMs for general SU (N ) theory

The models
The Haar-type PLM is constructed based on the two-parameter model proposed by Fukushima [29,42] for QCD with N = 3.It contains the nearest-neighbor coupling among PLs, so-called the kinetic term for the traced PL l, and as well as the Haar measure term H N [L] which is crucial to trigger the FOPT.The Haar measure incorporates (N − 2) degrees of freedom other than l, which causes the complexity for a larger N .In the Polyakov gauge where A 4 is diagonal, with i q i = 0, then H N [L] is given by The Haar-type PLM, generalized from the original model for N = 3 [29,42] to any color [32], thus goes like where the form of the kinetic term has been fixed inspired by the strong coupling expansion; C 1 (N ) and C 2 (N ) are two parameters needed to be determined by lattice data.Although it successfully explains the confinement-deconfinment PT of the first order, it can not reproduce proper thermodynamics observed in lattice simulations, as will be clarified later.
Soon later we will show that it is true even in another Haar-type PLM (for N = 3) studied in Ref. [43], which has a kinetic term including more parameters, generalized to any color number: The potential parameters a 0 , a 1 , a 2 and b 3 are constrained and fixed by imposing some physical conditions and fitting this potential to the lattice data on the thermodynamic quantities.

Constraining the potential parameters from first-order deconfinement PT
As to the required thermodynamic constraints, there are two.One is set by requiring that as noted in the Introduction, at high temperature the PL potential should realize the SB limit in the deconfinement vacuum l = l c , which must go to 1 as T → ∞ (because in this limit L → 1 by definition, see Eq. ( 2)) and therefore It fixes a 0 in Eq.( 10): The other condition is given by assuming that the FOPT takes place at the critical temperature T c , namely that two minima should be degenerate at T = T c .This condition gives a set of equations where the labels c and 0 attached on the Ls respectively stand for the values measured in the deconfinement phase (T > T c ) and the confinement phase (T < T c ).The last condition in Eq.( 12) gives the stationary point in the confinement phase, which has to be realized at L 0 = l 0 = 0 so as to reflect the Z N center symmetry in the confinement phase.Note from the potential form in Eq.( 10) that this condition is controlled solely by the Haar measure term, no matter what values independent of the potential parameters a 0 , a 1 , a 2 , b 3 take.The remaining first two conditions in Eq. ( 12) read i.e., Thus, given the Haar measure term H N , the field value l c can be numerically determined by solving the second relation in Eq.( 14), as is independent of the potential parameters.With l c , we combine the first and second relations in Eq.( 14), and obtain the relation between the potential parameters

Fitting the potential parameters from thermodynamics
The potential parameters should also be fitted to the thermodynamic quantities measured in lattice simulations, such as pressure P , energy density ϵ and the entropy density s.From the potential in Eq. (10), those thermodynamic quantities can be evaluated as In light of GW signals, of particular importance among the thermodynamic properties is the latent heat for the confinement-deconfinement PT of the first order.Given the PL potential as in Eq.( 10), the latent heat L N for SU (N ) PYM theory is defined as where ∆V (T ) = V (L 0 , T ) − V (L c , T ).Taking into account the first relation in Eq.( 12), we have To this L N , lattice simulations tell us [18][19][20][21][22][23] L which implies log This is the potential parameter relation for the Haar-type PLM, which must be satisfied to be consistent with the vicinity of the FOPT in the SU (N ) PYM theory.Actually, using b 3 given in Eq. ( 15), the above equation turns out to be a relation among a i : where a 0 has been fixed as in Eq. (11).It turns out, however, that the potential parameters in the Haar-type PLMs are inconsistent with the lattice data on thermodynamic quantities in SU (N ) gauge theory, in terms of the large N scaling.We come to this conclusion with great help from the work of Kubo [32], and let us explain the details in the following.From the numerical result [32] that no matter how large N is taken, l c is found to be around 0.5.Taking into account this result and a 0 = 2(N 2 −1)π 2 45 in Eq.( 11), we have Thus the parameters a 0 , a 1 and a 2 can scale with (N 2 − 1) for a sufficiently large N .On the contrary, actually, the parameter b 3 cannot do because of nontrivial N dependence of the Haar measure part present as a prefactor in Eq. (15).As references, we know [32] that for large enough N 1 and N 2 .Therefore, the parameter b 3 cannot simply scale with (N 2 − 1).This is in contradiction with the lattice result [22] showing the universal large N scaling for all thermal quantities ∝ (N 2 − 1), which indicates that all Polyakov loop potential parameters should be proportional to N 2 − 1 as well.This incompatibility arises from the Haar measure part included in the potential, driving the nontrivial N dependence into the thermodynamics in SU (N ) gauge theory.As a concrete example, let us take a look at the case with N = 4. Then we have log Hence Eq. ( 20) reads On the other hand, from Eq.( 15), we have Combining them gives Substituting these conditions into the fitting program for the thermodynamic quantities in SU (4) gauge theory, with data available from [22], we have found the fit does not match the lattice result at all5 .Thus, it is concluded that the Haar type PL potential in Eq.( 10) for generic SU (N ) gauge theory possesses incompatibility between the latent heat, i.e., confinement-deconfinement phase transition feature, and the overall large N scaling for the thermodynamic quantities.Now that we have demonstrated the Haar type PLM in Eq. (10) to have incompatibility between the latent heat nature of the FOPT and the large N scaling, we need to consider other PLMs.For N = 3, there is a very successful potential structure of polynomial form with powers of l not exceeding 4 (renormalizablility) [44], where the FOPT is triggered by the well-known cubic term l 3 + l * 3 by virtue of the underlying Z 3 [38].Nevertheless, this minimal polynomial form obviously fails for a generic Z N .On the other hand, the large N scaling property revealed by lattice data strongly indicates that there must be a simple way to construct the PLM, irrespective of N not less than 3.
Thus we propose a simple PLM containing three terms, All the terms in Eq. ( 28) are allowed by Z N symmetry, and even keep a global U (1) symmetry.Actually, it turns out that this enhanced U (1) symmetry allows us to derive a unified large N scaling property of the PLM.Hence we will ignore other terms that may be allowed by the Z N symmetry, for instance l 4 + l * 4 for N = 4, and so on.
Comments on the potential parameters in Eq.( 28) are in orders.It is assumed that a(T ) = a 0 + a 1 ( Tc T ) + a 2 ( Tc T ) 2 + a 3 ( Tc T ) 3 as usual, but b and c are T -independent.In order to fit the sQGP region 1.2T c < T < 4T c , the pressure may take the simple form shown in Eq. ( 1), and consequently the a 1 -term namely the T 3 term (in terms of the potential V (l, T )) is not important while the a 2 -term is important.However, the a 3 -term namely the T term is crucial to fit the region very close to T c , which will be confirmed by fitting to the lattice data later (See Figs. 1, 3 and Table I).A different viewpoint has been proposed by Fukushima in Ref. [42]: The PL potential is merely dominant to the contributions to the pressure of the system below about 2T c , and it should give way to the transverse gluons at higher temperatures; therefore, the PLM giving a pressure below the lattice data for T ≳ 2T c is acceptable.We will investigate both scenarios (with or without a 2 ) in the lattice fitting procedure, and find that, as expected, they are almost "degenerate" PLMs in a view of the GW signals (See Fig. 7).
Our PLM is characterized by two features: 1) The quartic term with a negative sign b < 0; 2) a nonrenormalizable term with a positive sign c > 0 to stabilize the potential.These two terms create the deconfined vacuum at l ̸ = 0.In the mean field theory, the nonrenormalizable terms usually are not taken into account, but its legitimacy, along with the sign of the quartic term, may be justified in the presence of additional and non-minimally charged PLs defined in Eq. ( 5).This potential form frees from extra large N scaling, other than the overall factor of (N 2 − 1) for all the potential parameters at around the SU (N ) gauge criticality, as will be clearly seen below.
Let us briefly explain how the 4-6 structure of the PLM in Eq.( 28) can arise following the spirit of Pisarski outlined in Ref. [37].The key is the presence of charge-two PL, l 2 , which always allows a cubic coupling with the fundamental PL, l, like l 2 l * 2 .The most general two-field model is given by with taking a 12 to be real.Assuming m 2 to be heavy while taking m to be light around T c , then l 2 can be integrated out.Among others, this procedure produces four-point interactions |l 1 | 4 with a negative coupling −a 2 12 /m 2 2 ≡ − λ 2 < 0. The negativeness of the induced quartic coupling is robust, because, in terms of the ll → ll scattering process, it arises from the l 2 -scalar exchange involving the trilinear a 12 interaction, which gives rise to the repulsive potential among the two l probes (in the nonrelativistic limit).The total effective PLM below the m 2 scale reads Hence, if the bare coupling λ 1 is properly chosen to make λ 1 < λ 2 , the total quartic coupling is naturally negative, triggering the FOPT.Moreover, when we turn off either λ 2 or λ 12 for simplicity or minimality, then the above model respectively gives rise to the 4-6 (as in Eq.( 28)) or 4-8 structure as desired -which we shall call the 4-6 and 4-8 PLMs, respectively.The case N = 2 is an exception since the confinement-deconfinement PT becomes second order, which may be accidentally due to the smallness of λ 2 .Actually, in this case l 2 is the "quarkless-baryon" loop [37], which is a singlet thus supposed to develop the vacuum expectation value any temperature.This vacuum expectation value should be very small, otherwise it would contribute significantly to the pressure of the system, contradicting with the lattice data.This argument is applicable to any other quarkless-baryon loops.

Lattice fitting of the potential parameters in the 4-6 PLM
The potential parameters are constrained as follows.First, the trace of PL field (l) needs to asymptotically get close to 1 when one approaches the classical limit where T → ∞.This condition can be realized by taking or assuming the SB limit value, as done in Eq. (11).The second condition comes from the FOPT nature, as in Eq.( 12): We have particular interest in the vicinity of the confinement-deconfinement PT at T = T c and consistency with the large N scaling for the thermodynamic quantities around there.In that case, it turns out that the parameters a 1 and a 2 are fairly insensitive to this FOPT criticality: Near the critical temperature T c , the potential can be simplified as V = 1 2 3 i=0 a i |l| 2 + b|l| 4 + c|l| 6 , where among the parameters a i , the a 3 term is most dominant because of its lowest power contribution by T c T around T = T c in the potential.Hence the potential around T = T c should be almost controlled by the a 3 term, so should the pressure.Taking into account the classical limit satisfied by Eq.( 31), we see that the a 0 term is also relevant.
The relevant parameter space is then constrained by the conditions in Eqs.( 31) and ( 32), together with the desired large N scaling of the latent heat in Eq.( 19) along with Eq.( 17).In the case with a 1 = a 2 = 0, those constraints read Thus, all the potential parameters can scale with (N 2 −1) in the large N limit, consistently with the lattice observation.Figure 1 shows that the 4-6 PLM fits to the lattice SU (N ) thermodynamics with data from [22], with focusing on the confinement-deconfinement FOPT at around T = T c .The left panel has been drawn by assuming a 1 = a 2 = 0, while the right panel made with assuming nonzero a 2 .Comparison of both panels indeed implies the insensitivity of a 1 and a 2 to the thermodynamics around the criticality.See also Table I, which more precisely shows the best fit to the lattice data points for the N = 4 4-6 PLM without those parameters.Of importance is to note that the 4-6 PLM with N = 4, which has been fixed by fitting to the lattice data as above, can be used to construct another 4-6 PLM with N ≥ 5: Given certain a 3 for some SU (N ) PLM, using the large N scaling a 3 → M 2 −1 N 2 −1 a 3 from N to M and Eq. ( 33), one can immediately get a 4-6 PLM for SU (M ).This is possible because the last equation in Eq. ( 33) tells us that, in the large N limit, the parameter a 3 is approximately proportional to N 2 − 1 and other equations do not contain direct N dependence.Note also that the lattice data tell us thermodynamic quantities ∝ N 2 − 1 and the large N scaling of the latent heat ≈ N 2 − 1.Therefore, if N is large enough then all relevant parameters for the SU (M ) 4-6 PLM should scale like and similar scaling laws for a 1 and a 2 .Note here that in performing the fitting procedure, quoting the available lattice results has been restricted to those which exhibit both the thermodynamics quantities and the latent heat on the same lattice setting.Otherwise, we will create a "prospected" data from the available ones, by the large N scaling.Currently, only the N = 3, 4, 6 lattice simulations have so far been worked out for both those two [22] 6 .Since the reported errors on the data for the thermodynamics quantities are quite small [22], we work on the least square (LS) fit, and evaluate the size of the systematic error by the root-mean square-deviation (RMSD).
Figure 2 and Table I show how the large N scaling law works for SU (5) and SU (6) in Fig. 1 (and also the same Table I), where the desired data are available only for N = 6, from [22].
Thus the large N scaling procedure allows us to save the fitting steps.The procedure will work better if the starting color number is set to bigger: Let a from N = m 2,3 | N =m be the naively scaled a 2,3 at N = n from a 2,3 at N = m.These scaled parameters are compared to a 2,3 | N =n at N = n determined fully by the fitting procedure, as follows: One can thus see that as N gets larger, the naively rescaled values get closer to the ones determined by the fitting procedure.FIG.2: Check on the large N scaling for the 4-6 PLM.Left panel: Fitting the 4-6 PLM for SU (5) to the thermodynamic quantities "measured" in the lattice simulation, where the lattice data for N = 5 have been read off by using the large N scaling for the N = 6 data in [22], as noted in the text.Right panel: The same for SU (6).For both N = 5 and N = 6 cases, we have selected a3 as the scaling parameter (scaled from the N = 4 case), and the other parameters have been fixed by the criticality and the latent heat constraints given in Eq. (33).See Table I (4-6 PLM with N = 5 and N = 61) for the values of the scaled a3 and the other fixed parameters for SU (5) and SU (6) cases, and the goodness-of-fit.In both panels the LS fits have been performed for 1 ≤ T /Tc ≤ 1.5 focusing only around the criticality, as zoomed-in the close-up windows.

Lattice fitting of the potential parameters in the 4-8 PLM
Following the same method as done in discussing the 4-6 PLM above, we can analyze the 4-8 PLM.The potential parameter relation then takes the form similar to Eq. ( 33): Thus the 4-8 PLM can also realize the desired large N scaling for the thermodynamic quantities in SU (N ) PYM theory consistently with the large N property of the latent heat.The result for fitting to the lattice data is shown in Fig. 3, and the best-fit parameters are listed in Table I.
Similarly to the 4-6 PLM, this 4-8 PLM can also create a set of simple large N copies because of the manifest large N scaling: staring with the 4-8 PLM model having color N with the potential parameters fixed, one can easily get the corresponding set of the potential parameters for another model with color M , just by accessing the simple large N scaling for a 3 and fitting to the (rescaled) lattice data, to determine the other potential parameters.

A. Dark confinement-deconfinement FOPT in a secluded gluonic plasma
The GW associated with the confinement-deconfinement PT would be generated from the bubble nucleation of the confinement vacuum.The bubble nucleation rate per spacetime, Γ, contains two contributions.The first one comes from its thermal fluctuation, Γ(T ), and the second one generated by the quantum effects.In the present study, we ignore the second contribution because the quantum effects would dominate only in supercooling.
The thermal contribution to the bubble nucleation rate per volume per time is generically given by  I.The LS fit has been worked out only for 1 ≤ T /Tc ≤ 1.5 (focusing on the criticality regime) in the left panel, as zoomed-in the sub-small window, while in the right panel for 1 ≤ T /Tc ≤ 4. The fitting parameters for the 4-6 and 4-8 PLMs.The two groups on the list for the cases with N = 4, 6, 8 separate two categories: one is the case where only a3 is used as the fitting parameter (corresponding to N1), and the other case treats all a1, a2, and a3 as the fitting parameters (N2).The n denotes the number of data and the LS fits for the N = 4 PLMs have been performed based on the data for T ≤ 1.2Tc.As to the cases with N = 5 and N ≥ 7, which are directly unavailable from [22], the prospected data have been read off by using the large N scaling from N = 4 or 6 data, as noted in the text, and the fitting has been done by using the method described in the text.Though their RMSDs (=LS/ √ n) are relatively poorer, the relative errors (RE), defined by RMSD 2 /(N 2 − 1), remain small.This implies that the fitting still works well.
where A is supposed to be of O(1), and S 3 (T ) is the O(3) symmetric action in three-dimensional Euclidean space: with R = |⃗ x|.V (ϕ, T ) denotes an effective potential for two vacua (false vacuum and true vacuum), and ϕ(R) is the bounce solution satisfying the Euclidean equation of motion with V ′ ≡ ∂V (ϕ, T )/∂ϕ, and the boundary conditions lim R→∞ ϕ(R) = 0 (at the false vacuum position) and dϕ(R) dR | R→0 = 0.
One thing to note is that the generic bounce solution has mass dimension 1, so that the exponent S 3 (T )/T in Eq.( 37) is dimensionless.On the contrary, in the PYM case, the role of the bounce solution is played by the trace of the PL field, l, which is dimensionless.To match the case with the generic thermal evolution in Eq.( 37), we need to rescale l so as to give mass dimension 1 to l, like l → lΛ, and fix the scale parameter as Λ = T .This scaling would be reasonable because l, or L in Eq.( 2), involves only the temperature T as the dimensionful parameter.Similar prescription has been applied in the literature [39].
Another thing to note is the distinction between the temperature in the dark gluonic bath and that in the bath made of the standard model (SM) particles.We label them respectively as T (for the dark PYM sector) and T (for the SM sector).Since PYM theory incorporates no matters, the two baths should decouple and thermalize individually.In Γ(T ), Eq.( 37), the prefactor T 4 originates from the spacetime volume evolved by the scale factor a of the universe as V = a 3 t, so should be identified simply with neither T nor T .However, the difference in T and T actually does not matter so much since it merely gives a logarithmic correction to the nucleation temperature T n .Whereas both temperatures in S 3 (T )/T (the argument of S 3 and the normalization factor for the exponential, 1/T ) should be associated with the dark PYM-sector: In S 3 (T ), T is supplied from a given PLM, and the origin of 1/T is traced back to the imaginary-time formalism in the thermal field theory, where T should be the temperature of the thermal bath created in the thermalized dark PYM.
In the temperature region of interest (say, before the big-bang nucleosynthesis (BBN)) where both baths are in the form of radiation7 , then with T 1 = T and T 2 = T , and similar for others.Here ρ i and g * i are the energy density and relativistic degrees of freedom in the bath, respectively.Further assuming that the entropy is conserved, then both temperatures scale as 1/a, and are expanding with the Hubble parameter As a consequence, in the thermal epoch under consideration, the ratio ζ ≡ T /T keeps a constant, determined by the production mechanism of the two sectors.To the end of having a correct relic density of dark glueball dark matter, ζ ≫ 1 is favored [6].But it will significantly suppress the amplitude of GW.Thus in this article we just take ζ to be merely a parameter, and suppose a situation that the dark glueball is not the main component of dark matter, or not the dark matter candidate at all.Such a setup justifies our treatment later, when the predicted GWs are discussed in comparison between the case with ζ = 0 and ζ = 1 (Figs. 5 and 6).Anyway, in what follows we will use T only, and the SM-sector temperature is simply ζT .We use the Mathematica program AnyBubble [45] to numerically solve Eq. ( 39), and evaluate S 3 (T )/T and the bubble nucleation rate Γ(T ) in Eq. (37).The criterion for bubble nucleation is that at the nucleation temperature T n , a single bubble is nucleated within one Hubble horizon volume; it roughly amounts to the condition Γ(T n ) ∼ H(T n ).In the radiation dominated universe, this condition reads S 3 (T n )/T n ∼ 140, and see Ref. [46] for more detailed discussions.Through this equation we determine the bubble nucleation temperature T n .
With the bounce action and T n at hand, we can calculate two key parameters characterizing FOPT, which almost completely determine GW spectra.One is the α parameter, which measures the strength of the FOPT via the ratio between the latent heat and the total energy density of the universe at The other one is β, the inverse time scale of PT duration in terms of the Hubble time scale, with t n being the cosmic time corresponding to T n .It is convenient to replace t with T , via the relation dt dT = −1 HT , and then β goes like Now, both α and β are calculated as a function of the dark -PYM sector temperature T (given certain ζ).The resultant numbers for the parameters relevant to the GW spectra are give in Table II.We see that the typical β is as huge as 10 5 , namely a very fast PT.That short duration of PT suppresses the production of GW.This result turns out to be fairly irrespective to whether ζ is zero or nonzero (See Figs. 5 and 7, and related discussions around there).
Let us end up with this subsection by stressing an interesting observation of β.That is, it shows a minimum at some color number N m .This feature has been pointed out in Refs.[36,47], but we would like to stress that the value of N m depends on the model.For the model in Ref. [47], N m = 6, while N m = 4 in Ref. [36] which adopts a very different model, the matrix model.In our model, we find that N m = 7(9) for the 4 − 6(8) PLM; see Table .II.
The above feature may be partially understood.We can prove that β decreases with N by the large N scaling.For a sufficiently large N , Eq. ( 39) can be solved once and for all, due to the scaling behavior of the effective potential V N = N 2 U with U independent of N .We assume that temperature is fixed, hence shall drop its dependence.Writing the bounce solution for the potential U as ϕ(R), we see that the bounce solution for V N is simply given by ϕ N (R) = ϕ(N R).With this solution, the three-dimensional Euclidean action S 3 for color N can be rewritten in the following form: Now it is evident that S 3 simply scales as 1/N .Note that our T n is always very close to T c , which is determined by U universal for the large N .Thus the nucleation condition S 3 /T n ≃ 140 leads to a slightly larger T n when increasing N .Moreover, the slope of S 3 /T sharply increases when the temperature approaches T c .Therefore, β, which is proportional to the slope, increases with N .Taking into account the fact (not proved, just a numerical observation) that β decreases with N for a small N , we deduce that β reaches the minimum at some moderately large N .

B. GW sources from the fast PT
In the cosmic FOPT, there are three kinds of sources for GWs, i.e., bubble collision, sound wave and magnetohydrodynamic (MHD) turbulence.The latent heat released during the FOPT distributes into these sources, but the fractions are not precisely determined yet.Usually, it is believed that the fraction of the latent heat transferred to the bubble collision is negligible, κ col ≪ 1, provided that α does not become extremely large [48].The GW sources are dominated by the two bulk motions of the plasma.
One bulk motion is the sound wave propagating in the plasma after the percolation, and the latent heat that goes into it is estimated to be [49] suppressed by the small α.The GW peak frequency at T n is f sw, * = 2(8π) 1/3 /[ √ 3(v w − c s )R * ] with R * being the average bubble separation at the collision and the speed of sound in the plasma c s = 1/ √ 3.This separation is related to the typical time scale of PT via R * = (8π) 1/3 v w /β n , given the exponential approximation of Γ(T ) around T n .The observed GW spectra peak at f sw = f sw, * a 0 /a(T n ), as a result of redshift from the GW production time t n to today.Then the peak frequency is parameterized as [48] f sw = 1.65 × 10 −5 T n 100GeV g * 100 T n 100GeV g * 100 1 6 Hz.(46) The amplitude of the GW spectrum for the sound wave is then given by [48,50] v w S sw (f ), (47) where the shape of the spectrum takes the form The amplitude is proportional to τ sw , the duration of the sound wave source, defined as with U f being the root-mean square -fluid velocity [51] A smooth expression for H * τ sw has been derived in Ref. [52], and it is simply reduced to H * R * /U f in the limit H * R * /U f ≪ 1.This is true for almost all confinement-deconfinement PTs, which will generate huge β.Consequently, the sound wave source is considerably suppressed by such a factor.So, in this case, the GW from the sound wave becomes Regarding the bubble wall velocity, we first note that the background field l does not interact with the SM particles, and therefore the bubble expansion seems not to be hampered by the SM plasma.Even without the SM plasma, the dark gluonic plasma itself would provide friction during the bubble expansion, which should be too nonperturbative to simply be quantified.In the present analysis, we will take the bubble wall velocity v w = 1 to shed light only on the maximized amplitude.We will come back to this point later in the discussion section.
When the sound wave period ends and the fluid flow becomes non-linear, then the other fluid bulk motion, the MHD turbulence is generated.If τ sw is relatively long, at least one Hubble time scale, then the MHD turbulence is suppressed, having an efficiency factor κ turb ∼ 0.05κ sw [51].However, as shown above, τ sw ≪ 1/H * and consequently κ turb may be significantly enhanced.But the quantitative enhancement is unknown yet owing to the fact that a part of the fluid motion not consumed by the sound speed wave can convert into heat.Optimistically, we assume that all of the energy transfers to turbulence, giving rise to the GW spectrum [53], with the shape function which, compared to S sw (f ), shows a moderately large suppression ∼ O (10) in the high frequency region.The peak frequency is similar to that of the sound wave source [48], T n 100GeV g * 100  Predicted GW spectra Using the formulae above combined with the 4-6 PLM in the previous section, we plot the GW spectra predicted from the dark SU (N ) PYM theory in Fig. 4 for N = 4 marked as "4-6 PLM", in comparison with the prospected GW interferometer sensitivities in the future (LISA8 [7,8], BBO [9][10][11][12][13], DECIGO [13][14][15] and TianQin [16]).In the figure, as a reference, we set ζ = 0 (i.e.no SM contributions), which can ideally be realized in the large N limit.The critical temperature T c for the confinement-deconfinement PT has been set to 270 MeV, just as a reference point inspired by the PYM limit for QCD 9 .In Fig. 4, we have also made comparison with the Haar-type PLM in Eq. ( 10) (labeled by "Haar-type"), though it is inconsistent with lattice data on thermodynamic quantities in terms of the large N scaling as well as the confinement-deconfinement criticality scaling in T .From Fig. 4 we see that the predicted GW signal from the N = 4 4-6 PLM can have high sensitivity reach to be probed by the prospected BBO detector.Note also that the 4-6 PLM's signal is significantly different and gets larger than the Haar-type PLM's.Varying the input T c value, we also plot the predicted GW signals for N = 4, in Fig. 5.In the left panel, we have taken ζ = 0 as in Fig. 4, while in the right panel ζ = 1.As clearly seen from the figure, different critical temperatures will only shift the peak frequencies of GW and hardly affect the peak amplitudes.We also find that the SM contribution give no significant contribution to the GW spectra from the dark PYM.In Fig. 6, similar plots with ζ = 1 have been displayed for N = 4, 5, 6.   4) SU( 5) SU( 6 II. In Fig. 7 we show the N = 4 4-8 PLM prediction to the GW spectra without the SM term (ζ = 0).The left panel displays the two scenarios with a 2 = 0, i.e. the category N = 4 1 , or with a 2 ̸ = 0, i.e. the category N = 4 2 , while in the right panel, comparison between the 4-6 PLM and 4-8 PLM has been shown for the case with a 2 = 0, i.e. the category N = 4 1 .In the figure we have taken ζ = 0 and T c =270 MeV.We see that both 4-6 and 4-8 PLMs would be possible to be probed by BBO.
Taking larger N , the 4-8 PLM sensitivity gets larger and reaches the maximum at N = 9, as expected from the minimum point of β in Table II.The left panel of Fig. 8 shows the predicted signals for N = 7, 8, 9 with the SM contribution included (ζ = 1) and T c = 10 GeV fixed, and the right panel displays the T c dependence for the N = 9 signals.The figure implies that the GW spectra from the 4-8 PLM with larger N can be probed at earliest by BBO and DECIGO.

IV. CONCLUSION AND DISCUSSIONS
A pure SU (N ) Yang-Mills (PYM) sector, for instance, predicted by the string theory, may be an active part of particle physics beyond the standard model (SM).Maybe, such a PYM sector is completely secluded to the SM sector, and then we should peek into this dark sector by means of its gravitational wave (GW) signals, which are generated during the first-order (FO) confinement-deconfinement phase transition (PT).
The gauge invariant object, the Polyakov loop L, associated with the global Z N symmetry of SU (N ) (called the center symmetry), is a well defined order parameter of this PT.And a proper effective model based on L, namely the Polyakov loop model (PLM), successfully describes the confinement-deconfinement PT.In this article we first looked into the widely used Haar-type PLMs, and demonstrated that they fail to produce the large N scaling for the thermodynamic quantities and the latent heat at around the criticality of the PT reported from the lattice simulations.We then proposed a couple of PLMs with polynomial terms, dubbed the 4-6 PLM and 4-8 PLM, with a negative quartic term to trigger the FOPT.Such a structure may naturally arise in the presence of a heavy PL with the center charge 2, given that the fundamental or normal PL having charge 1.We showed that those models give the desired  thermodynamic and large N properties at around the criticality.The predicted GW spectra were shown to have high enough sensitivity to be probed in the future prospected interferometers such as BBO and DECIGO.
The dark PYM opens novel interesting possibilities in the early universe, and here are some open questions need to be addressed elsewhere: • In discussing the bubble velocity, in Sec.III, we have argued that the bubble is expanding relativistically.However, as was mentioned in the text, the dark gluonic plasma should presumably be nonperturbative, so may provide nontrivial friction during the bubble expansion.In this sense the present analysis might merely give an ideal prediction to the created GW signals having its maximal size.More precise estimates might involve a quasi particle description for the dark gluonic plasma.
• The universe may experience a stage of PYM dominance.In this scenario, the dark glueball has a relatively short lifetime, decaying away before the onset of BBN due to the proper higher dimensional operators linking the two sectors.
• As can be read off from Table II, the parameter β decreases with the number of color N , implying that the PT in our models tends to happen more slowly, as N gets larger.This tendency is still operative even starting from N = 3.This decreasing β in N is a characteristic feature of our PLM, which is compared to the PT nature predicted from the matrix model [36] exhibiting a growing scaling of the β in N for N > 4, but a decreasing scaling for N < 4. At N = 4 the parameters α and β have the same order of magnitude for both the matrix model and ours, so both models coincide at this N = 4 and GW signals do as well.As one can see from Fig. 6, the peak strength of GW gets smaller as N decreases, hence it will be impossible to detect for N = 3 by the prospected detectors.Thus the GW detection for the larger N signals would be a good discriminator between the matrix model and ours, where the latter's signals will be much stronger because of the much slower FOPTs.More detailed analysis along this line would be worth pursuing.

4 .
No Go in a View of Large N Scaling and the SU(N) PYM Criticality

3 FIG. 1 :
FIG.1: Fitting the 4-6 PLM to the thermodynamic quantities of SU (4) PYM theory observed in the lattice simulations[22], at around the critical temperature (Tc) for the confinement-deconfinement phase transition of the first order.In the left panel, several thermodynamic quantities are plotted as a function of T /Tc with the best fit parameters a0 = 4.95, a1 = 0, a2 = 0, a3 = −6, b = −2.19 and c = 2.28, while those with a0 = 6.36, a1 = 0.04, a2 = −2.68,a3 = −4.72,b = −2.28 and c = 2.58 are displayed in the right panel.The LS fit has been worked out only for 1 ≤ T /Tc ≤ 1.5 (focusing on the criticality regime) in the left panel, as zoomed-in the close-up small window, while in the right panel for 1 ≤ T /Tc ≤ 4.

3 FIG. 3 :
FIG. 3: Fitting the N = 4 4-8 PLM to the thermodynamic quantities of SU (4) PYM theory.Left panel: Thermodynamic quantities as a function of T /Tc with the best fit parameters listed in Table I corresponding to the case N = 41 for the 4-8 PLM.Right panel: The same for the case N = 42 with the best-fit parameters in the list of the 4-8 PLM, in TableI.The LS fit has been worked out only for 1 ≤ T /Tc ≤ 1.5 (focusing on the criticality regime) in the left panel, as zoomed-in the sub-small window, while in the right panel for 1 ≤ T /Tc ≤ 4.

FIG. 4 :
FIG.4: GW spectra predicted from different models for N = 4 under an extreme scenario that the secluded dark gluon plasma to highly dominate over the SM plasma (i.e.ζ = 0).The critical temperature has been set as Tc = 270 MeV.For more details, see the text.

FIG. 5 :
FIG. 5: GW spectra with the same condition as in Fig. 4, but with varying Tc in several orders.The right panel includes the SM contribution by ζ = 1, while the left panel doe not (ζ = 0).

FIG. 6 :
FIG. 6: GW signals from the 4-6 PLM for N = 4, 5, 6 with Tc = 100 GeV (left panel) and Tc = 100 MeV (right panel).The SM contribution has been set as ζ = 1.The signals with N = 7, 8 for 4 − 6 PLMs will look like nearly overlapping with the N = 6 signal in this plot, because of almost degenerate βs as read off from TableII.

FIG. 8 :
FIG. 8: Left panel: GW signals of 4 − 8 PLM PTs with N = 7, 8, 9 at Tc = 10 GeV.Right panel: The dependence of the critical temperature Tc on the GW signals from the N = 9 4-8 PLM.Both panels have taken into account the SM contribution (ζ = 1).

TABLE II :
The reference values for the parameters relevant to the GW spectra, predicted from the 4-6 and 4-8 PLMs with N = 4 − 8, and N = 4 − 11, respectively.The subscript attached on 4 of N = 4 denotes the same category as introduced in TableI.The SM contribution parameter ζ has been set to 1.