Dark confinement and chiral phase transitions: gravitational waves vs matter representations

We study the gravitational-wave signal stemming from strongly coupled models featuring both, dark chiral and confinement phase transitions. We therefore identify strongly coupled theories that can feature a first-order phase transition. Employing the Polyakov-Nambu-Jona-Lasinio model, we focus our attention on SU(3) Yang-Mills theories featuring fermions in fundamental, adjoint, and two-index symmetric representations. We discover that for the gravitational-wave signals analysis, there are significant differences between the various representations. Interestingly we also observe that the two-index symmetric representation leads to the strongest first-order phase transition and therefore to a higher chance of being detected by the Big Bang Observer experiment. Our study of the confinement and chiral phase transitions is further applicable to extensions of the Standard Model featuring composite dynamics.


Introduction
Future gravitational-wave (GW) observatories provide new opportunities to investigate the existence of dark sectors that are currently inaccessible. Because the experiments will be sensitive to strong first-order phase transitions in the early universe, it is paramount to understand the landscape of theories that can lead to GW signals. Within the landscape of theories, asymptotically free gauge-fermion systems are privileged, being already chosen by nature to constitute the backbone of the Standard Model (SM). It is therefore reasonable to expect them also to appear in the dark sector [1][2][3][4][5][6][7][8][9][10][11][12][13][14].
In [15], we embarked on a systematic investigation of the composite landscape by providing the first comprehensive study of the dark confinement phase transition stemming from pure gluonic theories. There, we first extended the state-of-the-art knowledge of the confinement phase transition to arbitrary number of colours by combining effective approaches with lattice results and then determined their GW imprints.
Using our work [15] on pure gauge dynamics as a stepping stone, we now investigate the dark dynamics stemming from the confinement and chiral phase transitions arising when adding fermions in different matter representations to gauge theories. For previous works and complementary approaches on chiral and confinement phase transitions see e.g. [16][17][18][19][20][21][22][23][24][25][26][27][28]. We highlight below the main findings of our work: 1. We have systematically investigated the chiral and confinement phase transitions including their interplay for SU(3) Yang-Mills theory with fermions in fundamental, adjoint, and two-index symmetric representations via the Polyakov-Nambu-Jona-Lasinio (PNJL) model.

2.
Representations matter: fermions in the two-index symmetric representation increase the strength of the first-order confinement phase transition while the fermions in the adjoint representation decrease it. Notably, for the adjoint and two-index symmetric representation cases with one Dirac flavour, the chiral phase transition is a secondorder phase transition.
3. For all representations (chiral and confinement phase transitions), the inverse duration of the phase transition is large,β ∼ O(10 4 ). We discuss the large value ofβ in the thin-wall approximation where it is given by a competition of the surface tension with the latent heat. This sheds light on generic features of non-abelian gauge-fermion systems and helps finding models with stronger gravitational-wave signals. 4. The confinement phase transition with fermions in the two-index symmetric representation has the best chance of being detected by the Big Bang Observer (BBO) with a potential signal-to-noise ratio SNR ∼ O(10).

5.
We further provide an outlook of which strongly coupled models can potentially lead to a first-order chiral phase transition. The methods used in this work can be readily applied to some of these models as well as other dark and bright extensions of the Standard Model featuring new composite dynamics.

JHEP01(2022)003
This work is structured as follows. In section 2, we introduce the PNJL model as an effective low energy model which we use here to describe the dynamics in the dark gauge-fermion sector. In section 3, we discuss the details and interplay of the chiral and confinement phase transitions, including the bubble nucleation and the GW parameters and spectrum. In section 4, we present our results which includes a detailed analysis of the GW spectrum for the various models and the corresponding signal-to-noise ratios. In section 5, we discuss our results in the light of the thin-wall approximation and analyse which further models are likely to have a first-order phase transition. In section 6, we offer our conclusions.

Basic considerations
A first-order phase transition in the early universe generates a GW signal that might be detectable by future observations [29][30][31][32][33][34][35][36][37][38][39][40], for recent reviews see, e.g., [41,42]. This offers unprecedented detection possibilities for dark sectors that are otherwise (almost) decoupled from the visible SM sector. In this section, we concentrate on the description of such dark phase transitions, which fall into two categories: 1. Confinement phase transition characterized by the restoration of the centre symmetry at low temperatures [43]. The order parameter is the traced Polyakov loop [44].
2. Chiral phase transition characterized by the spontaneous breaking of chiral symmetry at low temperatures. The order parameter is the chiral condensate.
These phenomena occur in the strongly coupled regime of the gauge dynamics and cannot be described by perturbation theory. To make further progress, three approaches can be envisioned: 1. Universality analysis [45][46][47]. This is useful for investigating the order of the phase transition but does not provide a quantitative way to compute thermodynamic observables near a first-order phase transition.
2. First-principle non-perturbative approaches, e.g. lattice gauge theory [48] and functional renormalization group [49]. Unlike pure-gauge theories, first-principle results for thermodynamic observables are very limited for gauge-fermion theories in the chiral limit (i.e. zero fermion mass limit).
3. Effective theories [50,51]. These are constructed by including the relevant degrees of freedom and enforcing symmetry principles. They allow for a quantitative framework to compute thermodynamic observables near a first-order phase transition. However, there is the possibility that they do not provide a faithful modelling of the phase transition dynamics. The results obtained from effective theories provide valuable hints about possible dynamics of the underlying gauge-fermion theory, but should always be interpreted with care.

JHEP01(2022)003
Model name Gauge group Fermion irrep Reality N f KMT term Centre symmetry Chiral symmetry breaking pattern Table 1.
The three dark gauge-fermion models studied here, see section 2.1 for a detailed explanation.
In this work, we employ effective theories to obtain a quantitative description of the phasetransition dynamics in dark gauge-fermion sectors. For the description of the chiral phase transition at finite temperature, chiral effective theories such as the Nambu-Jona-Lasinio (NJL) model [52,53], see [54][55][56][57] for reviews, and the Quark-Meson (QM) model [58,59] are frequently adopted in the literature. In order to account also for the confinement phase transition, these models have been generalized to include the Polyakov-loop dynamics, with quarks propagating in a constant temporal background gauge field. The resulting models are called Polyakov-Nambu-Jona-Lasinio (PNJL) model [60,61] and Polyakov-Quark-Meson (PQM) model [62,63], respectively. See also [64] for a study of the Polyakovextended linear sigma model. Here, we focus on the PNJL approach and leave the PQM approach for future work.
We study three dark gauge-fermion models as shown in table 1. In the "Fermion irrep" column, F, Adj, and S 2 denotes the fundamental, adjoint, and two-index symmetric representation, respectively. The "Reality" column refers to the reality property of the fermion representation, which is complex (C) for F and S 2 , and real (R) for Adj. The "KMT term" column indicates the number of fermions needed to form a Kobayashi-Maskawa-'t Hooft term [65,66] (i.e. the 't Hooft determinantal term), e.g. 6F means the KMT term is a six-fermion interaction. For the 3G1 and 3S1 models the KMT terms are 12-fermion and 10-fermion terms respectively, and in our PNJL treatment, their effects are ignored since they correspond to operators of very high dimensions. The chiral symmetry breaking patterns for 3G1 and 3S1 models are also marked with an asterisk to indicate that the effects of KMT term are not considered (so the U(1) A part is included in the chiral symmetry).
We restrict our attention to the SU(3) gauge group for simplicity. The treatment of the PNJL grand potential in gauge groups of larger rank requires introducing two or more independent Polyakov-loop variables and is left for future work. Three smallest fermion representations F, Adj, and S 2 are then the natural targets for investigation. For the fundamental representation case, we consider the three Dirac flavour case, which is most likely to exhibit a first-order chiral phase transition (see discussion on universality below). The phase transition and GW signature of the 3F3 model have been also studied in [19], using several effective theories including the PNJL model. Compared to [19], we employ a different regularisation as discussed below. For the adjoint and two-index symmetric representations, we consider one Dirac flavour since the case of two Dirac flavours is believed JHEP01(2022)003 to be close to the lower boundary of the conformal window [67][68][69]. 1 We work in the chiral limit, i.e. setting all current quark masses to zero for simplicity. In a more generic setup, nonzero current quark masses could be introduced in such cases, leading to pseudo-Nambu-Goldstone bosons (pNGB). We leave this more generic case for future study.

The PNJL model
The PNJL Lagrangian can be generically written as [60,61] where L pure-gauge is the pure-gauge part of the Lagrangian whose effect is to contribute as the Polyakov-loop potential in the full grand potential, to be described below. L 4F and L 6F are the multi-fermion interaction terms which exist in the NJL model. They read Here the same notation ψ for the dark quark fields, and G S for the four-fermion coupling are adopted in all three models to avoid proliferation of new symbols. Note that ψ is a colour triplet in 3F3, a colour octet in 3G1, and a colour sextet in 3S1. It is understood that in the above equations, the colour indices of ψ are contracted to form singlets inside the fermion bilinear it resides in. For example,ψ and ψ are contracted with a Kronecker delta in colour space, while ψ C and ψ are contracted with a unitary symmetric matrix P * in the colour space in the 3G1 case. 2 In the 3F3 model case, ψ is also a three-component vector in the flavour space and we write ψ = (u, d, s) T when we want to make its individual components explicit. Note that these are dark quarks, not to be confused with the SM quarks. The λ a are 3 × 3 matrices in flavour space with λ 0 ≡ 2 3 1 and λ a , with a = 1, . . . , 8, are the usual Gell-Mann matrices written in the flavour space, normalized as Tr F (λ a λ b ) = 2δ ab (Tr F denotes trace in the flavour space only). Finally, L k is the covariant kinetic term for the quark field [60,61] with A µ = δ µ 0 A 0 being the temporal background gauge field living in the corresponding fermion representation. In the Polyakov gauge [60,72], A 0 can be taken to be diagonal and 1 If the two Dirac flavour case is inside the conformal window, then there is no confinement and chiral symmetry breaking at low temperature. It is also possible that the two Dirac flavour case is just below the lower boundary of the conformal window, associated with a weakly first-order phase transition in the sense discussed in [70]. Some nice discussions also appear in [71].
2 ψ C is defined as ψ C ≡ Cψ T as usual, with C being the charge conjugation matrix. The existence of the unitary symmetric matrix P is guaranteed in the 3G1 case because the adjoint representation is real.

JHEP01(2022)003
static. In the mean-field approximation to be introduced later, A 0 is also taken to be spatially homogeneous so it acts as a constant imaginary chemical potential [60] (A 0 = −iA 4 when continued to Euclidean spacetime). This way of coupling the gauge field to the quark field allows us to investigate the interplay between confinement and chiral dynamics in a convenient manner. However, only temporal gauge fields play a role in the modelling here. It is expected that for high temperatures (a few times the confinement phase-transition temperature), the transverse gluons are also important and the PNJL modelling should be revised accordingly [61,73]. A few remarks are in order regarding the construction of L 4F and L 6F in the NJL model. In principle one should write down to a given order (e.g. four-fermion level) all operators that are compatible with the full symmetry (spacetime, colour, and flavour) of the theory, with each independent operator carrying an independent coefficient [74,75]. A further complication arises due to the fact that in computing fermion loops, both the direct term and the exchange term may contribute. One can achieve simplifications by considering a Fierz-invariant Lagrangian from the beginning, which can be obtained by adding to the original Lagrangian its Fierz transformation [55,75]. With the Fierz-invariant Lagrangian, one only needs to consider the direct terms in the computation [55]. Fortunately, in many cases (including the present work) we do not need to carry out the exercise of writing down the full Fierz-invariant Lagrangian compatible with the symmetry of the theory. This is because one may work in a mean-field approximation and only care about the condensate in some but not all channels. For example in the mean-field approximation here, we are only concerned with the condensate ψ ψ since we work at finite temperature but zero chemical potential. We therefore only need to retain four-fermion terms related to this particular channel and express them in a form that preserves the full symmetry of the theory. Adding the Fierz transformation amounts to a redefinition of the couplings in the terms relevant for calculation. Since the coupling G S and G D are left arbitrary at this stage, we may also stick to the Fierz-non-invariant Lagrangian and take into account only the direct terms in the computation, without loss of generality [57].
The chiral symmetry of the L 4F term in the 3F3 case can be made manifest by introducing composite fields Φ ij ≡ψ jR ψ iL , (Φ † ) ij ≡ψ jL ψ iR with i, j = 1, 2, 3 being flavour indices and making use of the following identity for the four-fermion term [56,76] which can be proven by brute force using the explicit form of the λ a matrices. Φ transforms as (3,3) under the chiral symmetry group SU(3) L × SU(3) R , and thus Tr The six-fermion term L 6F in the 3F3 case is a parametrization of the U(1) A -breaking instanton effect. It is also chirally-invariant because it can be written as G D (detΦ + h.c.) and the SU(3) L and SU(3) R transformation matrices have unit determinants. In the 3G1 model, the L 4F term can be rewritten using the Nambu spinor

JHEP01(2022)003
in which Ψ C ≡ P ψ C L (P ψ C R ) C and the 2 × 2 Σ matrices are defined by Σ 1 = 1, Σ 2 = iσ 1 , Σ 3 = iσ 3 . The unitary matrix P is a symmetric matrix in colour space such that P ψ C R has the same gauge transformation properties under the gauge group as ψ L . This rewriting of L 4F using Nambu spinors makes the SU(2) × U(1) A symmetry manifest in the 3G1 case, since Ψ transforms as a doublet in the Nambu space, while Ψ C ΣΨ can be shown to transform as a complex three-vector in the Nambu space. 3 Finally, in the 3S1 case, the L 4F term in eq. (2.4) is already manifestly invariant under the chiral U(1) L × U(1) R symmetry.
The finite-temperature grand potential of the PNJL models can be generically written as (2.8) Here and * denote the traced Polyakov loop and its conjugate (to be defined more precisely below). V PLM [ , * ] is the Polyakov-loop potential describing the glue sector where all the potential can be fully determined by the existing lattice results (see e.g. [15]). V cond ψ ψ represents the condensate energy while V zero ψ ψ denotes the fermion zeropoint energy. The medium potential V medium ψ ψ , , * encodes the interactions between the chiral and gauge sector which arises from an integration over the quark fields coupled to a background gauge field. Each part is described in detail in the following sections. Before moving on it is important to note that PNJL models are non-renormalizable due to the multi-fermion interactions. We truncate the models to the six-fermion operators at most and use a 3D cutoff to obtain meaningful predictions from divergent momentum integrals. The cutoff Λ should be understood as a model parameter [55]. This regularization scheme is most convenient for finite-temperature computation and has been widely used in the NJL literature.

Polyakov-loop potential
Pisarski proposed the Polyakov-Loop Model (PLM) in [78,79] as an effective field theory to describe the confinement-deconfinement phase transition of SU(N ) gauge theories. The fundamental (traced) Polyakov loop plays the role of the order parameter (Tr c denotes the trace in colour space) where

JHEP01(2022)003
is the thermal Wilson line at temperature T , P denotes the path ordering along the time direction, and A 4 is the Euclidean temporal component of the gauge field in the fundamental representation (with the gauge coupling absorbed). 4 The symbols x and τ denote a spatial point and the Euclidean time, respectively. Note that an ordinary gauge transformation of the gauge field reads with the transformation matrix V (x) ∈ SU(N c ). A centre symmetry transformation is defined to be a transformation in the form of eq. (2.11) with V (x) satisfying a twisted boundary condition (β ≡ 1/T ) [51,80] V ( The connection between the centre symmetry and confinement is due to the fact that the free energy of a single static quark in the fundamental representation F q is related to the thermal average of the traced fundamental Polyakov loop [51,80] exp(−βF q ) ∝ . (2.14) So an unbroken centre symmetry implies the vanishing of , which in turn implies F q is infinite, and vice versa. Note the above discussion of centre symmetry applies to the pure-glue sector. The coupling to dynamical fermions may or may not explicitly break the centre symmetry, depending on the quark representation [80]. For the three models in table 1, fermions in the fundamental and two-index symmetric representation of the SU(3) gauge group leads to explicit breaking of the centre symmetry, while fermions in the adjoint representation preserve the centre symmetry. This is indicated in the "centre symmetry" column of table 1.
We adopt the Polyakov gauge [60,72] in which A 4 is diagonal and static. Also, in the spirit of mean-field approximation, we consider A 4 to be spatially homogeneous [81]. Then is independent of x. Introducing also the conjugate traced Polyakov loop * and the notation | | 2 ≡ * , the simplest effective potential preserving the Z Nc symmetry in the polynomial form is given by (2.16)

JHEP01(2022)003
Parametrization We have chosen the coefficients b 3 and b 4 to be temperature independent following the treatment in [15,51,61], and also neglected higher orders in | |. For the SU(3) case, there is also an alternative logarithmic parameterization which includes the information on the Haar measure, 5 see e.g. [51,84], given by PLM have been determined with a dedicated fit in [15] to available pure-glue lattice data [85]. The results are shown in table 2.
The temperature T 0 in eqs. (2.16), (2.18) is identical to the critical temperature of the confinement phase transition T c , in the pure-glue case. In the presence of fermions, the relation of T 0 to T c is slightly modified, which we will discuss later.

Condensate energy
The condensate energy V cond ψ ψ can be viewed as a tree-level contribution from the chiral condensate ψ ψ to the grand potential V PNJL . It can be derived in a self-consistent mean-field approximation [19,56,86] in which one introduces auxiliary fields for the condensate and splits the original Lagrangian into a mean-field part L MFA and a residual interacting part L res so that in the Bogoliubov-Valatin (BV) vacuum defined by a set of self-consistent conditions (SCC), L res has vanishing expectation value, and the SCC coincide with the equations of motion for the auxiliary fields derived from L MFA . For computing the condensate energy, the procedure can be simplified to a linearization of the fermion bilinears around the condensate [57,87]. Moreover, as explained in section 2.2, we only need to consider direct terms.
(2.19) 5 The Haar measure in the context of confinement physics is the Jacobian of the variable transformation from the gauge potential to the Polyakov loop [82,83].

JHEP01(2022)003
Then the trick is to rewrite (ūu) 2 as (ūu) 2 = [(ūu − ūu ) + ūu ] 2 = (ūu − ūu ) 2 + 2 ūu (ūu − ūu ) + ūu 2 − ūu 2 + 2 ūu ūu , (2.20) where in the last step the (ūu − ūu ) 2 term is dropped in the spirit of the mean-field approximation. In the remaining terms, the 2 ūu ūu term contributes to the constituent quark mass of u which plays an important role in determining the zero-point energy and medium part of the potential. The − ūu 2 term leads to a contribution to the condensate energy. Similar procedures can be applied to (dd) 2 and (ss) 2 , and to L 6F . In the chiral limit, the condensates should exhibit flavour universality, therefore we introduce, in the 3F3 case The condensate energy is then which is consistent with [19,88] after conversion of notations and conventions. The σ 3 term that originates from the 't Hooft determinental interaction turns out to be an important driving force for a first-order chiral phase transition. In the 3G1 and 3S1 cases, we define 3G1 and 3S1 : σ ≡ ψ ψ , (2.23) and their condensate energies are found to be In these cases, the 't Hooft determinantal interaction is associated with some very highdimensional operator which we neglect in the current approximation, and thus there is no σ 3 term.

Zero-point energy and medium potential
In the mean-field approximation outlined in section 2.4, the effects of the multi-fermion interaction terms L 4F and L 6F boil down to a contribution to the condensate energy shown in eqs. (2.22) and (2.24), and a contribution to the constituent quark mass as discussed below eq. (2.20). With the inclusion of the constituent quark mass, the covariant kinetic term for the quark field shown in eq. (2.5) becomes

JHEP01(2022)003
with A µ = δ µ 0 A 0 and the constituent quark mass M is given by (2.28) and the result for the 3F3 case is again found to be consistent with [19,88] after conversion of notations and conventions. In the Polyakov gauge and the mean-field approximation, A 0 = −iA 4 is diagonal and constant, acting as an imaginary chemical potential. The contribution of L k to the grand potential can be readily evaluated by a functional integration over the fermion at finite temperature and imaginary chemical potential as L k is quadratic in the fermion field. The calculation can be found in standard textbooks in thermal field theory [89]. The resulting contribution to the grand potential can be decomposed into a temperature-independent zero-point energy contribution V zero ψ ψ and a temperature-dependent thermal quark energy (called medium potential) contribution V medium ψ ψ , , * . The expression for the zero-point energy is given by [51,57] is the energy of a free quark with constituent mass M and three-momentum p, dim(R) is the dimension of the quark representation R, and N f is the number of Dirac quark flavours. The momentum integral is understood to be regularized by a sharp three-momentum cutoff Λ, which enters the expression for observables and is thus also a parameter of the theory. The integration can be carried out analytically and the result is [50] The spontaneous chiral symmetry breaking at zero temperature is the result of the interplay between the negative contribution from V zero ψ ψ which favours large values of M , and the positive contribution from V cond ψ ψ which favours small values of M [57].
The medium potential V medium ψ ψ , , * is evaluated to be [51] V medium ψ ψ , , in which S R and S † R are defined as

JHEP01(2022)003
Here µ is the chemical potential that we take to be zero. L R is the Polyakov loop matrix in the quark representation R, defined in a similar manner to eq. (2.10) with the Euclidean temporal gauge field A 4 now taken to be in the representation R. Accordingly, we also define the normalized traced Polyakov loop in the representation R .
We now proceed to evaluate S R and S † R at zero chemical potential. In the spirit of the meanfield approximation, we utilize the properties of the traced Polyakov loops that are satisfied at the saddle point of the grand potential. Especially, we note that the traced Polyakov loop R at the saddle point is always real, 6 and becomes equal to its conjugate * R at zero chemical potential [51]. With this in mind, in the grand potential we set R = * R [19,90], and the L R in the fundamental representation of SU(3) can be parameterized as assuming Polyakov gauge and spatial homogeneity. The phase θ is then related to R in the fundamental representation as S R and S † R in the fundamental representation can now be easily evaluated using the parametrization in eq. (2.36), and the result is To evaluate S R and S † R in higher representations, we note that S R is invariant under a similarity transformation in colour R-representation space, thus only the eigenvalues of L R matter for the calculation. L R can always be brought into a diagonal form via a similarity transformation, and its diagonal entries become pure phase , are weights of the irreducible representation R, and they can be obtained from the weights (i.e. eigenvalue phases) of the fundamental representation which we parameterize. For example, in the adjoint representation case, L Adj can be parameterized as

JHEP01(2022)003
with θ defined in the parameterization of L F in eq. (2.36). S R in the adjoint representation is then computed to be For the two-index symmetric representation, L S 2 can be parameterized as with the same θ introduced in eq. (2.36). S R in the two-index symmetric representation is then computed to be The expressions we obtained for S R in eqs. (2.38), (2.40), and (2.42) agree with [90].
If we consider a gauge group of larger rank, then the fundamental traced Polyakov loop F alone is not sufficient to describe the medium potential. For example, in the SU(4) case, we need two angles to parameterize L F L (4) where the superscript "(4)" is used to indicate the quantities associated with the SU(4) gauge group. The θ 1 , θ 2 angles can be traded for two independent traced Polyakov loops, in the fundamental and two-index antisymmetric representations, respectively A 2 simultaneously. This is true even for fermions in the fundamental representation S (4) and similarly for S R in higher representations (with more complicated expressions). This suggests treating the confinement dynamics and the interaction between the quark and gluon sectors not in terms of a single traced Polyakov loop in the fundamental representation but in terms of eigenvalues of the Polyakov-loop matrix, which goes in the line of

JHEP01(2022)003
the matrix-model approach [22,81,91,92]. Moreover, for the convenience of studying the bubble nucleation, some method needs to be introduced to reduce the multi-variable problem to the tunnelling in a single dimension [22]. The extension of the current work to these cases will be left for future study. The medium potential V medium ψ ψ , , * does not contain ultraviolet (UV) divergence and there are different procedures on the market regarding the regularization of this contribution, even within a 3D momentum cutoff framework. For example, in [60], no momentum cutoff is imposed on the medium potential and the 3D momentum is integrated to infinity. On the other hand, in [93], a sharp 3D momentum cutoff has been employed everywhere, including the medium potential. The choice is motivated by the authors' wish to describe certain mesonic properties. When it comes to quarks in higher representations, [90] regulates the medium potential by introducing a momentum-dependent four-fermion coupling which implies that for three-momentum larger than Λ, the medium potential is not set to zero, but rather computed as if the quarks have zero constituent mass. It is found in [90] that such a regularization treatment is needed to obtain clearly separated confinement and chiral phase transition in the case of adjoint quarks. 7

Model parameters and observables
Apart from the coefficients in the Polyakov-loop potential, the PNJL models we are considering have only two parameters (G S and Λ) in the 3G1 and 3S1 cases, and three parameters (G S , G D , and Λ) in the 3F3 case. In principle, these parameters should be determined from observables (meson masses, decay constants) measured from experiments or predicted in lattice calculations. However, because we work in the chiral limit, even in the 3F3 case it is difficult to determine the parameters precisely. In [19], four benchmark points were chosen to study the 3F3 model in the chiral limit. We deem it reasonable in the 3F3 case to use the parameters corresponding to physical real-world values as a reference and then investigate variations away from the physical point by some amount. The values for such a choice can be found in [57], see also [94]. In the 3G1 and 3S1 cases, however, there are no clear guidelines to determine the parameter and we thus allow the parameters to vary in a larger range. Nonetheless, we provide formulae for a set of observables to gain more physical insight from the model parameters that we use, and also to facilitate future comparison of computations done in different approaches (because a meaningful comparison should be carried out at the same value of observables). The set of observables that we consider include the chiral condensate σ, the constituent quark mass M , the pion-decay constant f π , and the σ-meson mass m σ .

JHEP01(2022)003
The chiral condensate σ is determined from the saddle point equation at zero temperature and chemical potential The pion-decay constant is determined from the vacuum to one-pion axial-vector matrix element, and we use the normalization in [55]. In the 3D momentum cutoff scheme, it is given by Finally, the σ-meson mass m σ is the root of the 1PI σ-σ 2-point function Γ σσ whose expressions are lengthy and thus given in appendix A. It turns out that in the 3G1 and 3S1 cases the σ-meson mass is simply twice the constituent quark mass.

Comments on universality
We have collected various pieces of the PNJL grand potential for the three models of our interest. Before moving to the exploration of the phase transitions based on the expressions for the grand potential, let us comment on the relation between the universality argument and our analysis. There is a nice summary of the logic of the universality argument in [95]: 1. One first assumes the phase transition is continuous. The asymptotic critical behaviour must be associated with a 3D universality class with the same symmetry breaking pattern as the original theory.
2. The existence of such a 3D universality class can be investigated by considering the most general Landau-Ginzburg-Wilson Φ 4 theory compatible with the given symmetry breaking pattern.
3. The critical behaviour at a continuous transition is determined by the fixed points (FPs) of the renormalization group (RG) flow: the absence of a stable FP generally implies first-order transitions. However, if a stable FP exists, the phase transition can be of second-order or first-order (if the theory is outside the attractive domain of the FP).
RG predictions for the type of chiral phase transitions in SU(3) QCD theories with quarks in the complex or real representations have also been summarized in [95]. The universality analysis can also be carried out for the confinement phase transition, see [50] for a summary.
Among the three models we are considering, 3F3 is predicted to exhibit a first-order chiral phase transition. On the other hand, for 3G1 and 3S1, no definite prediction can be made. It is therefore well-motivated to explore what type of phase transitions these models exhibit using an effective theory approach like PNJL.

JHEP01(2022)003
For quarks in the fundamental representation, the universality argument predicts a first-order phase transition also for a number of Dirac flavours larger than 3 (but below the lower boundary of the conformal window). On the other hand, our study suggests that for N f = 4 the PNJL approach does not seem to exhibit a first-order chiral phase transition even if the model parameters are allowed to vary in a large range. This does not mean the universality argument is invalid. Rather, it is likely that this points to the possibility that the PNJL approach fails to model the phase-transition dynamics faithfully in such a case. This motivates further modelling of the phase transition using alternative effective theories in order to deliver a first-order chiral phase transition compatible with the prediction of universality.

Confinement and chiral phase transitions
In this section, we discuss the nature of the confinement and chiral phase transition for the models studied in this work. We start with some cosmological considerations and the discussion of order parameters, followed by a generic review of the bubble nucleation and the resulting GW spectrum.

Cosmological considerations
We briefly discuss here the cosmological constraints on the dark pions in our model. This can be divided into the two cases of massless and massive dark pions.
Massless dark pions. The dark pions play the role of dark radiation which is strongly constrained by the cosmic microwave background (CMB). Following the nice discussion in [96], there are two viable cases: 1. The dark sector and visible sector are thermalized in the very early universe but decouple prior to the electroweak scale T ew . In addition, the chiral phase transition should happen even before the decoupling. In this case, for the SU(3) gauge group with fermions in the fundamental and two-index symmetric representation, only 1 ≤ N f ≤ 3 Dirac flavours are viable while in the adjoint representation, only one Dirac flavour is viable. Thus our choice of N f = 3 and N f = 1 in respectively fundamental and adjoint (two-index symmetric) representations is valid.
2. The dark sector and visible sector never thermalize. In this case, the key parameter from the CMB constraint i.e. the ratio between the hidden and visible sector tem- can be made arbitrarily small and thus avoid the constraints from CMB. Note that the GW signal is suppressed if the dark temperature is much colder than the visible temperature [97,98] and thus we do not consider this case.
Massive dark pions. This scenario is less constrained and can be achieved by adding small explicit quark masses. We assume the dark gauge sector is in thermal equilibrium with the SM at early times before the Big-Bang Nucleosynthesis (BBN) and that the pseudo-Nambu-Goldstone bosons (pNGBs) have decay channels into lighter SM particles.

JHEP01(2022)003
As long as they have a mass larger than a few MeVs and can decay before BBN, these pNGBs do not cause conflicts with cosmological and laboratory constraints [99], which we explain in more detail in appendix B.
In this work, we consider the simplest case where at the phase transition the dark and visible sectors are in thermal equilibrium, T d = T v . This implies that the dark pions can be massless for T c > T ew , while they have to be massive for T c < T ew .

Order parameters
In this section, we focus on the interplay between chiral and confinement phase transitions for the models studied here, see table 1. The order parameter of the confinement phase transition is the Polyakov-loop expectation value, while the chiral condensate is the order parameter for the chiral phase transition. Fermions in the fundamental and two-index symmetric representation explicitly break the Z 3 centre symmetry of SU(3) and thus the Polyakov loop is no longer a rigorous order parameter for the confinement phase transition. Nevertheless, the Polyakov loop can still serve as an indicator of a crossover between confinement and deconfinement [61,100].
We In the pure-glue case, both potentials have the property that → 1 for T → ∞, however, with the addition of the medium potential eq. (2.32) this property is lost in the polynomial case. In principle, one would need to refit the coefficients to the lattice data with the refined constraint of → 1 for T → ∞ that includes the properties in the medium potential. We refrain from doing so since this effect only becomes relevant at large T and is not relevant for the dynamics of the phase transition. For the purpose of the figures, we implement the constraint ≤ 1 by hand. We also emphasise that the logarithmic potential naturally implements the constraint → 1 for T → ∞ and therefore might be better suited for the considerations in this work. N f = 3 (figure 1). The chiral condensate has a discontinuity at the critical temperature and thus we have a first-order chiral phase transition. The Polyakov loop expectation value undergoes a cross over (there is a small discontinuity at T c due to the discontinuity in σ). The confinement crossover happens roughly at the same temperature as the first-order chiral phase transition. N f = 1 (figure 2). The fermions in the adjoint representation do not break the Z 3 centre symmetry and thus the Polyakov loop expectation value remains a good order parameter for the confinement phase transition. We find a firstorder confinement phase transition, while the chiral phase transition is of second order and happens at much larger temperatures.   N f = 1 (figure 3). The twoindex symmetric representation is a very interesting case. In principle, the centre symmetry is explicitly broken by the fermions and thus there should be no confinement phase transition. However, it turns out that the centre symmetry is only weakly broken [90]. mass. The latter is rather large in the two-index symmetric representation, see table 3, and consequently, the centre symmetry is only softly broken. As in the adjoint case, the chiral phase transition is of second-order and happens at larger temperatures. Therefore the centre symmetry is almost restored at T c , and we observe a first-order confinement phase transition. The small negative dip of the Polyakov loop expectation value is precisely due to the breaking of the centre symmetry induced via the medium potential eq. (2.32).

Bubble nucleation
In case of a first-order phase transition, the transition occurs via bubble nucleation and it is essential for the understanding of the dynamics to compute the nucleation rate. 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 ) [101][102][103][104] The three-dimensional Euclidean action reads where ρ denotes a generic scalar field with mass dimension one, [ρ] = 1, and V eff denotes its effective potential. In our case, the effective potential depends on two scalar fields, the Polyakov loop and the chiral condensate σ. Which field takes the leading role depends on whether we have a first-order confinement or chiral phase transition and therefore we discuss them separately.

JHEP01(2022)003
Confinement phase transition. The phase transition is described by the Polyakov loop and it is a first-order phase transition in the adjoint and two-index symmetric case, see figure 2 and figure 3. In both cases, the second-order chiral phase transition is at significantly higher temperatures and has already been completed. Therefore, we can work in the approximation that σ is constant. Note also that is dimensionless while ρ in eq. (3.2) has mass dimension one. We therefore rewrite the scalar field as ρ = T and convert the radius into a dimensionless quantity r = r T . Thus, the action becomes with the associated boundary conditions To attain the solutions, we used the method of overshooting/undershooting and employ the Python package CosmoTransitions [105].

Chiral phase transition.
The chiral phase transition is described by the chiral condensate σ, see eqs. (2.21), (2.23). In the three models studied here, we only find a first-order chiral phase transition in the fundamental case, see figure 1. In order to have a field with mass dimension one, we defineσ We work in the mean-field approximation where we evaluate the Polyakov loop for given values ofσ and T at the minimum of the effective potential. Thus the potential becomes a function of only (σ, T ), V eff (σ, T ) = V eff (σ, T, min (σ, T )).
Sinceσ is not a fundamental field, we have to include its wave-function renormalization Z σ , see appendix A for more details. In figure 4, we display the wave-function renormalization as a function of the chiral condensate and the temperature. The three-dimensional Euclidean action is slightly modified [19] The bubble profile is obtained by solving the equation of motion of the action in eq. (3.7) and is given by

Inverse duration time
An important parameter for determining the GW spectrum is the rate at which the phase transition completes. For sufficiently fast phase transitions, the decay rate can be approximated by where t * is a characteristic time scale for the production of GWs to be specified below. The inverse duration time then follows as The dimensionless versionβ is defined relative to the Hubble parameter H * at the characteristic time t * β

JHEP01(2022)003
where we used that dT /dt = −H(T )T . Note that here we assumed that the temperature in the hidden and visible sectors are the same, T d = T v . The phase-transition temperature T * is often identified with 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). More accurately one can 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 in the vacuum decay rate Γ, see eq. (3.10), and consequently, we use the percolation temperature throughout this work. We write the false-vacuum probability as [106,107] P (T ) = e −I(T ) , (3.13) with the weight function [108] (3.14) The percolation temperature is defined by I(T p ) = 0.34, corresponding to P (T p ) = 0.7 [109].
Using T * = T p in eq. (3.12) yields the dimensionless inverse duration time. We will see that all phase transitions considered here have very fast rates,β ∼ O(10 4 ).

Energy budget
We define the strength parameter α from the trace of the energy-momentum tensor θ weighted by the enthalpy 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 relations between enthalpy w, pressure p, and energy e are given by These are hydrodynamic quantities and we work in the approximation where do not solve the hydrodynamic equations but instead extract them from the effective potential with This treatment should work well for the phase transitions considered here, see [110][111][112]. With eqs. (3.16), (3.17), α is given by In case of the confinement phase transition, we find that the contribution from ∆V eff is negligible since e + p + and therefore α ≈ 1/3. In the case of the chiral phase transition, we find smaller values, α ∼ O(10 −2 ), which relates to the fact that there are more relativistic d.o.f.s participating in the phase transition. Note that relativistic SM d.o.f.s do not contribute to our definition of α since they are fully decoupled from the phase transition. The dilution due to the SM d.o.f.s is included at a later stage, see section 3.5.

Bubble-wall velocity
We treat the bubble-wall velocity v w as a free parameter. A reliable estimate of the wall velocity would require a detailed analysis of the pressure and friction on the bubble wall. The latter is typically evaluated in an expansion of 1 → n processes [113][114][115][116][117][118]. In our case, we have a strongly coupled system and most likely a fully non-perturbative analysis would be necessary to determine the friction.
We make the assumption that the wall velocity is larger than the speed of sound, v w ≥ c s = 1/ √ 3. In this regime, the wall velocity does not have a strong impact on the GW peak amplitude. For wall velocities smaller than the speed of sound, the efficiency factor decreases rapidly and the generation of GW from sound waves is suppressed [119].

Efficiency factors
The efficiency factors determine which fraction of the energy budget is converted into GWs. In this work, we focus on the GWs from sound waves, which is the dominating contribution for the phase transitions considered here. The efficiency factor for the sound waves κ sw consist of the factor κ v [120] as well as an additional suppression due to the length of the sound-wave period τ sw [121][122][123] In our notation, τ sw is dimensionless and measured in units of the Hubble time. It is given by [123] τ sw = 1 − 1/ 1 + 2 (8π) 1 3 v w βŪ f . (3.20) whereŪ f is the root-mean-square fluid velocity [121,124] We follow [120] for κ v where it was numerically fitted to simulation results. The factor κ v depends α and v w , and, for example, at the Chapman-Jouguet detonation velocity it reads For the confinement phase transition with α ≈ 1/3, this leads to κ v ≈ 0.45, and for the chiral phase transition with α ∼ O(10 −2 ), we have κ v ∼ 0.1.

Results
With the setup of the PNJL model in section 2 and the tools for phase transition in section 3, we are now ready to display our results.  Table 3. Table of parameters of the PNJL model with resulting GW parameters and observables for the benchmark and best-case scenario for T c = 100 GeV.

Choice of PNJL parameters
The fundamental QCD-type model of our dark gauge-fermion sector has only one parameter, which can be for example the critical temperature T c or the value of the strong coupling at any arbitrary scale. However, we are working in the effective PNJL model which has the advantage of being well suited to describe the strong dynamics at the phase transition, but also has the drawback of many model parameters, including the couplings G S and G D , the cutoff Λ, and the a i , b i coefficients as well as the temperature T 0 in the Polyakov-loop potential. The a i , b i coefficients have been fitted against pure-glue lattice data, see table 2, and we make the approximation that the coefficients remain the same in the presence of quarks.
In the 3F3 model, we have further guidance from the lattice and we use values from [51] as a benchmark, and rescale them such that, for example, T c = 100 GeV, see table 3. We also scan the PNJL parameters in the vicinity of this benchmark point and search for a best-case scenario that gives us the strongest GW signal and thereby explore the lattice uncertainty on the PNJL model parameters. Note that we keep the ratio Λ/T 0 = 3.8 fixed such that it agrees with [88].
In the 3G1 and 3S1 models, there are no lattice data available 8 to fit to as in the 3F3 case and therefore the PNJL parameters are basically unconstrained. Nonetheless, we use the benchmark values 9 from [90], again rescaled to, e.g., T c = 100 GeV, see table 3. We again scan the parameters for the best-case scenario, however, in these models we also vary the ratio Λ/T 0 due to the lack of constraints from the lattice data. Remarkably, although the parameters have few constraints, we find clear predictions for these models. 8 There are quite a few lattice works in the direction of fermions with higher dimension representations.
However, most of these focus on studying the phase structure between chiral symmetry breaking and (near)conformal phase (see e.g. [161,162]) rather than producing the mass spectrums which are important in determining the PNJL parameters.

Parameter scan
Fundamental representation. We vary the couplings G S and G D , while we keep the ratio Λ/T 0 fixed, since the latter is well constrained by lattice data. We parameterize the couplings with G S = k G S · 4.6 GeV −2 and G D = k G D · (−743 GeV −5 ), where k G S = k G D = 1 is the benchmark point from [51], which we choose as a starting point of the scan. Below a certain value of k G S , chiral symmetry breaking does not occur. This value depends on k G D and, for example, for k G D = 1, we have k G S ,crit = 0.882.
The resulting GW parameters and PNJL observables are displayed in figure 5. For the purpose of the plot, we keep either k G S or k G D equal to one, while varying the other. We observe that the strength of the phase transition is increasing with a larger G D and decreasing with larger G S , i.e.,β is decreasing and α is increasing with a larger G D . The phase transition cannot become arbitrarily strong, instead the GW parameters approach asymptotic values which we estimate via extrapolation to beβ asymp ≈ 3.5·10 3 and α asymp ≈ 0.06. Consequently, we use a large value of G D as best-case scenario in table 3.
In terms of PNJL observables, we observe that the constituent mass M and the sigmameson mass m σ grow roughly linear with G S and G D , while the pion decay constant remains roughly constant. Also we observe the approximate relation 2M ≈ m σ .

JHEP01(2022)003
Pure glue Adjoint and two-index symmetric representation. We vary the coupling G S and the cutoff Λ while the temperature T 0 is determined by the choice of T c . We start from the benchmark values in [90], see table 3. The resulting GW parameters and PNJL observables are displayed in figure 6. For the purpose of the plot, we keep one coupling fixed to the benchmark value while varying the other.
Both, G S and Λ, have a lower bound below which chiral symmetry breaking does not occur. We denote these lower bounds by G S,crit and Λ crit . If the respectively other coupling is at the benchmark point, they take the values G S,crit,adj = 2.68 GeV −2 and Λ crit,adj /T 0 = 4.41 in the adjoint representation case as well as G S,crit,sym = 14.3 GeV −2 and Λ crit,sym /T 0 = 3.24 in the two-index symmetric representation case.
For Λ/Λ crit 1, the fermions are decoupling from the theory and the pure-glue result is approached, see figure 6. This is most easily understood in terms of the constituent mass, which is growing with Λ 3 and triggers this decoupling. The large constituent mass makes the medium potential trivial, see eq. (2.32), and we are left with the pure-glue dynamics.
For Λ/Λ crit → 1, we observe that the fermion representation makes a big difference for the GW parameterβ. Fermions in the adjoint representation lead to larger values ofβ (weaker phase transition), while fermions in the two-index symmetric representation lead to smaller values ofβ (stronger phase transition).

JHEP01(2022)003
Also when varying the parameter G S , we observe that the two-index symmetric representation case always has smaller values ofβ than the pure-glue case while the adjoint representation case always has larger values ofβ. In terms of the PNJL observables, we see that the constituent mass grows linear with G S while the pion decay constant is decreasing with G S . For G S /G S,crit 1, the GW parameterβ becomes constant but does not approach the pure-glue result, although the constituent mass is growing linear with G S . The difference to the large-Λ limit is that we implement a momentum-dependent four-fermion coupling, see eq. (2.46), and therefore the mediums potential still gives contributions for momenta larger than the cutoff.
Consequently, we choose for the best-case scenario a small value of the cutoff in the two-index symmetric representation case and a large value of the cutoff in the adjoint representation case, see table 3.

Gravitational-wave spectrum
With the GW parameters displayed in table 3, we compute the GW spectrum as described in section 3.5. We use both, the polynomial and the logarithmic parameterization of the Polyakov-loop potential, see eqs. (2.15), (2.17). The central value of both parameterizations is our main result, β = (β poly + β log )/2, and we estimate the error of the parameters via the difference between the parameterizations, δβ = |β poly − β log |/2. Furthermore we vary the wall velocity between the speed of sound and light speed, c s ≤ v w ≤ 1. With the latter treatment, we define the error bands for GW signal. This is compared to the power-law integrated sensitivity curves of LISA [149][150][151], BBO [152][153][154][155][156], and DECIGO [156][157][158][159]. The power-law integrated sensitivity curves provide a qualitative visualisation for the detectability of a GW signal. Since the GW signals considered here are not simple power-law signals within the frequency range of the detector, we refer to the SNR for the quantitative discussion of detectability of the GWs, which is discussed in section 4.4.

Fundamental representation.
We present the GW spectrum from the chiral phase transition with fermions in the fundamental representation in figure 7. The best-case scenario has a peak frequency of h 2 Ω peak GW ∼ O(10 −17 ) and might be detectable with the BBO. The benchmark scenario has a peak frequency of h 2 Ω peak GW ∼ O(10 −18 ) and is out of reach of any planned future detector. Adjoint representation. We present the GW spectrum from the confinement phase transition with fermions in the adjoint representation in figure 8. The GW signal is strongly suppressed and out of reach of any planned future detector. It has a peak frequency of h 2 Ω peak GW ∼ O(10 −18 ). Remarkably, the spectra from the benchmark and the best-case scenario are almost identical, which shows that in this case, the PNJL model is highly predictive.

Two-index symmetric representation.
We present the GW spectrum from the confinement phase transition with fermions in the two-index symmetric representation in figure 9. Similar to the GW spectrum of the fundamental case, the best-case scenario has a peak frequency of h 2  benchmark scenario has a peak frequency of h 2 Ω peak GW ∼ O(10 −18 ) and is out of reach of any planned future detector.
Dependence on T c . In figure 10, we show the dependence of the GW spectrum on the critical temperature at the example of the best-case scenario of the fundamental representation. As expected, the critical temperature simply shifts the peak frequency of the GW spectrum. For the fundamental representation, we have the biggest overlap with BBO for T c ∼ 100 GeV. For the two-index symmetric representation, this happens at T c ∼ 50 GeV, and for the adjoint representation even at T c ∼ 10 GeV, see also section 4.4.
Sound waves vs turbulence. In figure 11, we display the comparison between GWs from sound waves and from turbulence at the example of the best-case scenario of the two- index symmetric representation. The sound-wave contribution is clearly dominating over the turbulence with the exception of frequencies far above and below the peak frequency. This is related to the slower fall-off behaviour of the turbulence contribution compared to the sound-wave contribution. We have checked that the sound waves are dominating the GW spectra for all phase transitions considered here.

Signal-to-noise ratio
In figure 12 we present the signal-to-noise ratio for the detectors BBO and DECIGO as a function of the critical temperature for the best-case scenario of the three models, see If we assume that for a successful detection we need SNR > 1, 10 then BBO will test the chiral phase transition in the 3F3 model for 50 GeV T c 200 GeV as well as the confinement phase transition in the 3S1 model for 10 GeV T c 200 GeV. The SNR at DECIGO is slightly below one for these models for all T c . A detection of the chiral phase transition in the 3G1 model is out of reach for both detectors. 10 In general, it is difficult to determine the exact threshold from which a signal is detectable. It also depends on how well the astrophysical foreground such as gravitational radiation from inspiralling compact binaries can be subtracted from the signal, see, e.g., [163][164][165]. Here we make the probably optimistic assumption that a signal with SNR > 1 is detectable.

Thin-wall approximation
All phase transitions in the gauge-fermion systems studied here, as well as the phase transitions in SU(N ) gauge theories without fermions studied in [15], have one thing in common: the inverse duration is remarkably largeβ ∼ O(10 4 ). Here we want to provide an instructive argument via the thin-wall approximation to explain this feature in an intuitive picture. The advantage of the thin-wall approximation is that we can analytically calculate the decay rate of the false vacuum in terms of the latent heat and the surface tension. Eventually, we will relate the large inverse duration time to a competition between the latent heat and the surface tension.
The thin-wall approximation for the Euclidean action was derived in [104,166] and we briefly review it here. The three-dimensional Euclidean action is written as where p f and p t denote respectively the pressure in the false vacuum and true vacuum, σ 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 in the pressure between the false vacuum and true vacuum is also linked to the latent heat L via The thin-wall approximation works with the assumption that η is small, η 1. With eqs. (5.2), (5.3), the three-dimensional Euclidean action eq. (5.1) can be written as a function of the latent heat L and the surface tension σ where we have made explicit that the surface tension and latent heat are evaluated at T c . The ratio S 3 (T p )/T p is typically a number O(150) for phase transitions around the electroweak scale. From this we infer that , (5.5) and the inverse durationβ follows as

JHEP01(2022)003
Model name Gauge group Fermion irrep Reality Number of Weyl flavours Centre Symmetry Table 4. List of gauge-fermion models that lead to cubic terms in the condensate energy of the corresponding NJL model. see also [39]. We see thatβ stems from the competition between latent heat L and surface tension σ. However, already the prefactor is O(10 3 ) and therefore we would need a surface tension that is much larger than the latent heat to achieve a strong first-order phase transition. We caution that this analysis relies on η 1, which is true for the phase transitions investigated here but must not hold for all strongly coupled gauge-fermion systems.
In the case of pure-gluon dynamics, lattice results have provided fitting functions for the surface tension and the latent heat as functions of the number of colours and T c [85,167]. With those lattice results, we find the small surface tension can not compensate the large latent heat and thus leads toβ ∼ O(10 4 ) [15]. We expect that it is quite common that in strongly coupled gauge-fermion systems, the latent heat is big while the surface tension is not sufficiently large and thus many small bubbles are formed (nucleate very quickly everywhere in the space-time) rather than large and long-lasting big bubbles which are crucial to generate a strong GW signal.

Models with cubic terms in the condensate energy
In the three gauge-fermion models we explored, a first-order chiral phase transition is only found for the 3F3 case which features a cubic (σ 3 ) term in its condensate energy, see eq. 2.22. For the 3G1 and 3S1 case such a cubic term is absent and we do not discover a first-order chiral phase transition in the parameter space. The relevance of the cubic term for the triggering of a first-order chiral phase transition is well-known [50]. In the 3F3 case, such a cubic term originates from the 't Hooft determinantal interaction. It is therefore tempting to ask the following question: suppose we consider a gauge-fermion theory featuring a simple gauge group and one type of fermion under a single irreducible representation of the gauge group (for complex representations, its conjugate does not count), can we list all such gauge-fermion models that deliver a six-fermion 't Hooft determinantal term so that their NJL condensate energy contain cubic terms? The answer is positive. It turns out there is JHEP01(2022)003 a finite number of such models, characterized by the condition coming from an analysis of the discrete chiral symmetry that must be preserved by the 't Hooft determinantal term with N W f being the number of Weyl flavours, T (r) is the trace normalization factor for the fermion representation r. The resulting solutions are listed in table 4, in which we also list the remnant centre symmetry for each case. The 3F3 model belongs to the M1 category. The models listed in the table are the natural targets for the next-step investigation. They require extensions of the effective theory framework used in this work so that larger gauge groups can be dealt with, as commented in section 2.5.
The possibility of having both a first-order chiral phase transition and a first-order confinement phase transition is intriguing. If they are separated in scale, a double-peak feature in the GW signature is in principle possible. Table 4 gives some hints on where to find candidate models with this feature. Here actually the requirement of centre symmetry is not that strict. If a nontrivial centre symmetry remains, one may have an idea about the order of the confinement phase transition from universality argument [50]. However, as commented in section 2.7, only the absence of a FP requires a first-order phase transition, while the existence of a FP does not lead to definite predictions. Therefore, although an SU(2) confinement phase transition (with Z 2 centre) is known to be of second-order, it is not straightforward to exclude the possibility of first-order phase transitions for other cases with Z 2 centre listed in table 4. Moreover, even if the remnant centre symmetry is trivial, there still can be a first-order confinement phase transition. There can be two possibilities in such a case. The first is that the centre symmetry is only broken weakly by dynamical quarks. This is the case when the chiral phase transition occurs at a scale a few times larger than the confinement transition scale (which is typical in higher representations, see [168]). When the temperature drops below the chiral phase-transition temperature, the quarks obtain a large dynamical mass and are decoupled from the low-energy gauge dynamics. A first-order confinement phase transition can then occur as if the quarks are absent. In fact, this phenomenon already takes place in the 3S1 case we studied and has been found earlier in [90]. The second is that the confinement phase transition might only be directly related to the change of some order parameter, but not related to the centre symmetry. For example, the G 2 gauge group is known to have a trivial centre even in the absence of dynamical quarks, but lattice studies suggest it exhibits a first-order confinement phase transition at finite temperature [169,170]. Modelling the confinement phase transitions in such cases is an interesting issue which we leave for future work. Finally, we would like to comment that the model M9 corresponds to the supersymmetric SU(3) gauge theory. A recent lattice study [171] suggests that it exhibits a single first-order phase transition where chiral symmetry is restored and centre symmetry gets broken at high temperature.

Conclusions
We studied the GW signal from strongly coupled gauge-fermion systems featuring both, dark chiral and confinement phase transitions. We employed the PNJL model to systemat-

JHEP01(2022)003
ically study SU(3) Yang-Mills theory with fermions in fundamental, adjoint, and two-index symmetric representations. We studied in detail the interplay between chiral and confinement phase transitions.
We discovered that the representation of the fermions matters: the two-index symmetric representation case leads to the strongest first-order phase transition and has the highest chance of being detected by the Big Bang Observer experiment with a potential signal-to-noise ratio of SNR ∼ 10. Conversely, fermions in the adjoint representation lead to the weakest GW signal. However, for all models considered here, the inverse duration time is large,β ∼ O(10 4 ), and the GW signal is generically suppressed. We analyse this observation through the thin-wall approximation and show that the large rate stems from the competition between the small surface tension and the large latent heat.
Beyond gravitation waves, our study of the confinement and chiral phase transitions can be readily employed for different models of composite dynamics for beyond standard model physics. In the future, it will be intriguing to study the impact of larger numbers of colours as well as adding dark scalars and dilatons, and therefore study the near-conformal dynamics [172].

Acknowledgments
MR thanks M. Hindmarsh for helpful discussions and S. van der Woude for correspondence on [19]. MR acknowledges support by the Science and Technology Research Council

A Wave-function renormalization
The NJL model Lagrangians are constructed purely from fermion fields, as shown in eqs. (2.2), (2.3), (2.4). This raises the question of how to evaluate the kinetic term in the bounce action for the chiral phase transition, and more generally, also the question of how to interpret bosonic d.o.f. in such models. As pointed out in [19], the mesons should be viewed as non-propagating at tree level, but their kinetic terms can be induced by quantum corrections due to fermion loops. Therefore, unlike models in which the order parameter field is elementary, in NJL models the kinetic terms for the order parameter field need to be computed as a quantum effect. For the three gauge-fermion models of our interest, only the 3F3 model exhibits a first-order chiral phase transition, and accordingly, we need to compute the kinetic term (i.e. wave-function renormalization) for the σ meson for the 3F3 case. This exercise has been done in [19] using a 4D momentum cutoff scheme, however, some key steps in the derivation are hidden, and some subtle issues are not emphasized.
Here we present the derivation for the PNJL model in the 3D momentum cutoff scheme

JHEP01(2022)003
which is the regularization scheme adopted in this work and much of the (P)NJL literature. This derivation reveals a number of subtle issues which deserve further investigation. We believe the presentation of the derivation here will also be helpful for future studies of GW signatures from a first-order chiral phase transition in other gauge-fermion models using a (P)NJL approach.
We now definē in the 3F3 case. The following equations are expressed in terms of the barred couplings G,Ḡ D and the barred chiral condensateσ. These barred quantities correspond to the corresponding unbarred quantities in App. B of [19]. Note however we stick to the 3D momentum cutoff scheme with Λ denoting the 3D momentum cutoff, while in [19] Λ denotes the 4D momentum cutoff. The 3D Euclidean bounce action can be written as with V PNJL being the PNJL grand potential, see eq. (2.8). If we explicitly display the functional (and function) dependence of S 3 , Z −1 σ , and V PNJL , we have In the PNJL model, the quark propagator is modified in the temporal background gauge field, which enters into the computation of Z −1 σ . Eventually Z −1 σ also depends on the traced Polyakov loop as we will show. To simplify the treatment of the bounce action, we adopt the mean-field approximation for the Polyakov loop, in the sense that we set [19]  Then Z −1 σ and V PNJL become functions ofσ and T only, and the bounce action becomes a functional ofσ and a function of T : The wave-function renormalization factor Z −1 σ can be computed from the derivative of the polarization propagator at finite temperature with respect to the spatial momentum squared [173] Because at finite temperature Lorentz invariance is lost, in the above equation we treat the temporal and spatial components of the external four-momentum q separately. In the JHEP01(2022)003 evaluation of Z −1 σ the derivative should be taken with respect to the spatial momentum squared in light of correspondence to the 3D bounce action. The extra minus sign in eq. (A.6) comes about because we adopt the metric signature (+, −, −, −). Γ σσ (q 0 , q,σ) is the polarization propagator for σ at finite temperature. It is related to Feynman graphs via In the (P)NJL model,σ is non-propagating at classical level. To make sense of the computation of Γ σσ (q 0 , q,σ) we work in the framework of a self-consistent mean-field approximation [19,56,86]. 11 The mean-field approximation Lagrangian at finite chemical potential reads One can read off the Feynman rules from the MFA PNJL Lagrangian in eqs. (A.8), (A.9). The vertex rules for each quark colour and flavour read and the quark propagator reads quark propagator : .
Note that iA 4 acts as a constant imaginary chemical potential. There is no scalar propagator since scalars are non-propagating at tree level. Because σ has a nonzero expectation value (at zero and finite temperature), the twopoint 1PI σ-σ graphs should contain those in which extra σ's are inserted, which eventually take their expectation values. We adopt the random phase approximation [55] which implies that extra σ insertions in the middle of quark propagators need not be taken into account. In this approximation, their effects are assumed to have been absorbed into the constituent quark mass. With this in mind, we arrive at a finite set of two-point 1PI σ-σ graphs from which we derive an expression for Γ σσ (q 0 , q,σ) In the self-consistent mean-field approximation σ is introduced as an auxiliary field. The vacuum expectation value of σ is also denoted σ and it is required by the self-consistent condition to satisfy eq. (A.1).

JHEP01(2022)003
where N f = N c = 3 and we have performed the analytic continuation q 0 = iω E , and ω E takes value at bosonic Matsubara frequencies 2nπT, n ∈ Z. The two loop integrals I V and I S are defined by introducing a modified loop momentum so that (A.14) and Here tr denotes the trace in Dirac space only, and T is defined as 16) with ω n = (2n + 1)πT, n ∈ Z taking values at fermionic Matsubara frequencies.
For the computation of Z −1 σ , only the term containing I S (iω E , q,σ) in eq. (A.12) is relevant because only this term depends on the external momentum. By simple algebraic manipulations, I S (iω E , q,σ) can be decomposed as with q 2 = −ω 2 E − q 2 , and I(iω E , q,σ) is defined as .
With this decomposition, Z −1 σ can be computed as Here I(0) and I (0) are defined as Now the crucial task is the computation of I(iω E , q,σ). This can be achieved by the use of partial fraction and the summation formulae [93] 1

JHEP01(2022)003
So in the 3D cutoff scheme at zero temperature and chemical potential, Lorentz invariance is not restored even when we send the cutoff to infinity. Because we indeed expect the restoration of Lorentz invariance at zero temperature and chemical potential, the above calculation reveals a limitation of the 3D cutoff scheme when applied to the evaluation of the Z −1 σ factor. A negative Z −1 σ factor signals an unphysical instability which is not acceptable in the computational framework. Motivated by the consideration of Lorentz invariance at zero temperature and chemical potential, we propose to remedy the problem by using Z −1 σ [temporal] for the zero temperature contribution in Z −1 σ , while the finite temperature contribution (terms that depend on f or its derivative) is unchanged. This amounts to modifying eq. (A.31) into with the subscript m in Z −1 σ,m indicating this is the modified wave-function renormalization factor. Fortunately, this modified Z −1 σ,m factor exhibits the expected positivity in the parameter space of interest, and thus we use it for the computation of bubble nucleation and GW signature.
Finally, as a byproduct of computing Γ σσ , here we also provide the formulae for evaluating the σ meson mass in the 3F3, 3G1 and 3S1 models. In all three models the σ meson mass m σ is determined from the equation while the expression for Γ σσ (ω, 0,σ) in the three models are given by , (A.39) 3S1 : Γ σσ (ω, 0,σ) = − 2 G S + 6 π 2 Λ 0 dp .
(A. 40) in which we introduced the principal value integral (denoted by P.V. in order to get rid of the singularity in the integrand when ω 2 = 4E 2 p , as done in [87,174]. In fact, further JHEP01(2022)003 simplification occurs due to the gap equations in the 3G1 and 3S1 cases, This simple relation is also a consequence of the fact that we are considering the chiral limit. On the other hand in the 3F3 case due to a nonzero G D , the cancellation is not perfect and m σ deviates from 2M .

B Cosmological and laboratory constraints on massive dark pions
In this appendix, we show how massive dark pions can decay prior to BBN (i.e. with a lifetime τ π 1 s) while being compatible with laboratory constraints. For example, consider the dimension-five Hermitian gauge-invariant effective operators O BE = (ψ iL σ µν ψ jR +ψ jR σ µν ψ iL )B µν , In the above equations, ψ denotes the dark quark fields, with i, j being flavour indices. B µν is the hypercharge gauge field strength tensor, while H is the SM Higgs doublet. The subscript "B" or "H" indicates the hypercharge portal or the Higgs portal respectively, while the subscript "E" or "O" indicates the CP property (even or odd). Replacing B µν with the dual field strength does not lead to new independent operators due to the selfduality property of the σ µν matrix. For simplicity, let us consider the case in which the dark QCD sector conserves CP. CP conservation means that we can only add the CPeven operators O BE and O HE to the theory. Note that O HE can lead to the decay of CP-even (but not CP-odd) dark pions, via mixing with the SM Higgs (see ref. [98] for a lifetime estimate). In order to allow the decay of both CP-even and CP-odd dark pions, let us consider the effect of the operator O BE . With only one insertion of O BE it is not possible to make dark pions decay as the gauge-invariant effective operator ∂ µ ∂ ν πB µν vanishes identically. Nevertheless, with two insertions of O BE one may form gauge-invariant effective operators π E B µν B µν and π O B µνB µν , with π E , π O denoting CP-even and CP-odd dark pions respectively. The existence of these operators relies on a three-point function of composite operators in the dark QCD sector which is not protected to vanish by any symmetry. Then, both π E and π O are allowed to decay to diphoton: π E → γγ, π O → γγ.

JHEP01(2022)003
With the above in mind we add the CP-conserving hypercharge portal interaction with c BE a coefficient of O(1) and Λ BE a suppression scale. Suppose this leads to gaugeinvariant effective operators The coefficients a E , a O can be estimated from dimensional analysis, which gives roughly Since our analysis of phase transition and gravitational wave is performed in the massless quark limit, in the massive dark pion scenario we hope to make dark pions as light as possible (but still heavy enough to decay prior to BBN). Thus as a benchmark we take M π = 10 MeV. 12 The typical dark pion decay constant relevant for our analysis can be taken as f π = 100 GeV. The dark pion lifetime can then be estimated as (for c BE ≈ 1) (B.5) Requiring τ π 1 s with the above values of dark pion mass and decay constant leads to an upper bound on the suppression scale Since we are considering a CP-conserving dark sector, we expect no constraints from electron or neutron electric dipole moment measurements or other CP-violating observables. The precision electroweak measurements and collider direct production can currently probe new physics scale of a few TeV at most, thus there exists a window of values for Λ BE which allows the decay of dark pions prior to BBN while being compatible with laboratory constraints. We also note that for f π ≈ 100 GeV, dark baryons are expected to have masses around 1 TeV and have a relic abundance that is smaller than the required amount needed to act as dark matter (assuming zero dark baryon asymmetry) [98]. In such a case the main component of dark matter would be composed of other particles whose impact on the phase transition and gravitational wave signals can be made negligible. 13 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. 12 Using Gell-Mann-Oakes-Renner relation, a dark quark gets bare a bare mass 1 MeV. We expect its impact on the phase transition of our interest which occurs around 100 GeV to be small but a dedicated future study is needed to make a quantitative estimate. 13 The main component of dark matter can be some heavy particle outside the dark QCD sector. Its energy density scales as ρ DM ∝ a −3 (a denotes the scale factor) while the radiation energy density scales as ρ R ∝ a −4 so the dark matter energy density is highly suppressed in the early Universe.