Gravitational Waves from Composite Dark Sectors

We study under which conditions a first-order phase transition in a composite dark sector can yield an observable stochastic gravitational-wave signal. To this end, we employ the Linear-Sigma model featuring $N_f=3,4,5$ flavours and perform a Cornwall-Jackiw-Tomboulis computation also accounting for the effects of the Polyakov loop. The model allows us to investigate the chiral phase transition in regimes that can mimic QCD-like theories incorporating in addition composite dynamics associated with the effects of confinement-deconfinement phase transition. A further benefit of this approach is that it allows to study the limit in which the effective interactions are weak. We show that strong first-order phase transitions occur for weak effective couplings of the composite sector leading to gravitational-wave signals potentially detectable at future experimental facilities.


I. INTRODUCTION
Present and future gravitational-wave (GW) detectors add extra dimensions to the way we observe and explore the world around us [1,2].The added dimension for beyond the Standard Model physics is the possibility to investigate a number of physical scenarios in a complementary way to direct observations, for example, at colliders [3] (for a recent review, see e.g.[4]).The existing and planned GW observatories such as NANOGrav [5,6] and Laser Interferometer Space Antenna (LISA) [7,8] can detect the stochastic GW background sourced by violent cosmological events such as first-order phase transitions in the early universe.Such events are typically very sensitive to the presence of new particles and interactions determining the shape of the effective potential at finite temperatures.For an overview of the existing new physics realisations that feature the strongly first-order phase transitions, see for instance Ref. [9], and for possible new physics implications of the recent NANOGrav measurement of the stochastic GW background, see e.g.Ref. [10] and references therein.
Among such possibilities, an exciting one is that our universe may feature new composite dynamics.Nussinov first considered this possibility in [11] within models of dynamical electroweak symmetry breaking in which dark matter emerges as dark baryons.It was, however, later understood that new composite asymmetric pions [12,13] are interesting candidates for dark matter trig-gering the first dedicated lattice study of composite dark matter [14].Having constructed concrete underlying extensions of the Standard Model in terms of new fundamental composite dynamics it became clear that their thermal history could be investigated via gravitational wave detection [15].
More recent applications of composite dynamics to the dark sector of the universe can be found in [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] literature.Our most recent contributions to the field started with asking the simple question: Can the simplest model of dark composite dynamics be detectable by present and future GWs observatories [35,36]?Arguably, the simplest model for dark composite dynamics is the one stemming from a purely dark gluonic sector interacting only with itself and gravity.Although strong dynamics are not yet fully solved the inevitable dark confinement phase transition occurring during the universe's evolution leads to a primordial background of GWs.In our work, we combined symmetries, effective approaches and state-of-the-art lattice results to arrive at precise predictions [35] (see also [36][37][38][39][40][41][42][43][44][45][46][47][48][49]).While much is still left to be understood even for such a simple model, it is a fact that composite dynamics is much richer and besides dark gluonic degrees of freedom can feature new light matter.The latter greatly alters the phase diagram of the theory adding to the confinement phase transition, typically also new ones such as the chiral phase transition [36,50,51], and even a quantum one when modifying the number of matter fields to achieve a conformal phase.
It is a fact that for such a wide spectrum of composite theories, we lack a full understanding because strongly coupled dynamics is yet unsolved.Among the various methodologies we could use to tackle the dynamics the lattice field theory would be ideal but, unfortunately, lattice methods are expensive and often still too immature to investigate the full dynamics of these theories.Other arXiv:2309.16755v1[hep-ph] 28 Sep 2023 methods make use of toy models of holographic nature [52,53].
We resort, here, to the time-honoured approach of effective theories such as the Linear Sigma Model (LSM).The reason for such a choice is two-fold.The first is that there is a great body of literature that allows us to test our results and the second is that in certain coupling regimes when the theory becomes weakly coupled, the computations allow us to extend our results also to perturbative models.The latter feature is, for example, absent in holographic approaches which by construction can only access some features of the strongly interacting composite dynamics.Specifically, in this work, we use the Polyakov quark meson model (PQM) also known as Polyakov-loop improved Linear Sigma Model (PLSM) [54][55][56] as an effective theory to study the firstorder chiral phase transition for a generic number of Dirac matter fields N f .As a finite-temperature computational framework, we use the Cornwall-Jackiw-Tomboulis (known as CJT) method which allows us to go beyond the mean-field approximation [57][58][59].As advertised above, the (P)LSM with the CJT method compared to other model computations such as holography [52,53,[60][61][62] and the Polyakov-loop improved Nambu-Jona-Lasinio (PNJL) model [50,51,63,64], can bridge perturbative and non-perturbative regimes of the effective theory corresponding to distinct underlying dynamical realizations.
Our main findings can be summarised as follows.Stronger first-order phase transitions appear for less strongly coupled composite models and light sigma masses expressed in units of the pion decay constant.When considering the phase transition as a function of the number of flavours we find that stronger phase transitions appear when increasing the number of flavours.However, exceptions to this behaviour occur depending on the interplay, for example, between the anomalyinduced terms and the values assumed by the effective parameters of the theory.Specifically, as the sigma mass in units of the critical temperature decreases the order of the phase transition increases for different N f .As mentioned above, intriguingly for N f = 3 such a dependence on the lightness of the sigma state suddenly amplifies when the mass drops around twice the critical temperature in part due to the presence of a relevant anomalyinduced operator in the action.
We then translate the information about the type and strength of the studied phase transitions into potential GW-related signals.We show that these will be observable at the proposed future BBO [65][66][67][68][69] and DECIGO [69][70][71][72] GW experiments but fall short of the LISA expected sensitivity [73][74][75][76].A summary of the sensitivity curves for the above measurements is given in Ref. [77].
The paper is organised as follows.In Sec.II, we present the theoretical formulation of the PLSM incorporating finite temperature effects via the CJT approach.In Sec.III, we briefly overview the standard picture of the first-order phase transition dynamics and characteristics of the resulting GW power spectrum.In Sec.IV, we present the main results of our analysis including the flavour dependence, the effect of the Polyakov loop sector as well as the impact of the determinantal term.The physics implications of our results are discussed in Sec.V. Finally, in Sec.VI we summarise our main conclusions.

II. THEORETICAL SETUP
In this section, we briefly summarize the effective actions we use to investigate the order and strength of the phase transition as well as the computational method used to extract the finite temperature corrections relevant to obtain the desired information.

A. Polyakov-Loop Improved Linear Sigma Model
The model we use combines information about confinement contained in the Polyakov loop with the physics of chiral symmetry encoded in the LSM.The resulting effective model has different names in the literature from the PQM to the PLSM.An alternative approach is the PNJL model [54-56, 63, 64], which prefers to couple the Polyakov loop to the non-renormalizable NJL model.It is fair to say that the approaches are complementary and feature their own set of well-known shortcomings.
The Lagrangian of the PLSM where mesons couple to a spatially constant temporal background gauge field A 0 reads where PLM +V LSM +V medium denote, respectively, the Polyakov-loop potential, the LSM potential and the interaction terms between the Polyakov loop and meson fields.All three potentials are explained in detail below.
In the pure gluon sector of a SU (N ) Yang-Mills theory, a centre symmetry Z N is used to distinguish the confinement and deconfinement phases.The corresponding order parameter is the colour trace of the Polyakov loop where F represents fundamental representation of the gauge group, and L F is the Polyakov loop defined as Here, A 4 = A a 4 T a and T a are the generators in the defining (fundamental) representation of the gauge group.More generically, we can define L R similarly with generators T R in a general representation R. In what follows, we focus on the fundamental representation and denote, for simplicity, L ≡ L F and ℓ ≡ ℓ F .
The thermal effective potential of the gluon sector which preserves the Z N symmetry in the polynomial form is described by the Polyakov loop as [78,79]: The above "• • • " represent any required lower-dimension operator than ℓ N i.e. (ℓℓ * ) k = |ℓ| 2k with 2k < N .In [80] it was first realized how to transfer the information about confinement encoded in the somewhat mathematical Polyakov loop model to the hadronic quasiparticle states of the theory and vice versa.The approach was later generalized to take into account an interplay between confinement and chiral symmetry breaking [81,82].
For the LSM we assume the symmetry to be SU (N f )× SU (N f ) for a given number of flavours N f and the potential reads [83] where the meson field Φ is a N f × N f matrix defined as where I denotes the identity matrix and T a is the corresponding SU (N f ) generators with a = 1, 2 The determinant term in (4) breaks the U (1) A symmetry and thus the potential (4) is invariant under Furthermore, the σ field in (5) acquires a vacuum expectation value (VEV) triggering the spontaneous chiral symmetry breaking i.e.
where f π is the pion decay constant.The interaction between the Polyakov loop and meson fields is described by the medium potential V medium ⟨ ψψ⟩, ℓ, ℓ * evaluated as follows [84] V medium ⟨ ψψ⟩, ℓ, ℓ in which S R and S † R are defined for a representation R as Here, µ is the chemical potential that is taken to be zero, and E p = ⃗ p 2 + m 2 q denotes the quark energy where is the constituent quark mass.In the PNJL model, the constituent quark masses contain in general, a linear as well as higher powers of σ stemming, for example, from the anomaly fermion operator.In this model, this information is now carried by the scalar potential term.In the fundamental representation with N = 3 colours, we have where, for simplicity, we have implemented the reality condition i.e. ℓ = ℓ * .For N > 3, it requires an additional traced Polyakov loop in S F expressions.
In a simultaneous expansion in powers of ℓ and in the number of meson fields while keeping chiral symmetry manifest, we expect operators of the type T 2 (g 1 ℓ + g 2 ℓ 2 )Tr Φ † Φ to naturally emerge [81] from the first principle computations.However, given the level of approximations made here, one generates only a mixing between ℓ and the scalar field σ which might be sufficient for a phenomenological analysis of the thermal phase transitions.
To summarize, the total effective potential is In the next section, we will use the Cornwall-Jackiw-Tomboulis (CJT) method to obtain the finite temperature effective potential of the LSM sector.To avoid double counting the pure gluon sector of the theory is encoded directly in (V PLM ) and the interplay with the scalar mesonic sector by the medium potential.The CJT treatment will be reserved for the purely mesonic sector of the theory.

B. CJT Method
Cornwall, Jackiw and Tomboulis first proposed a generalized version of effective action Γ (ϕ, G) of composite operators [57], where the effective action not only depends on ϕ(x) but also on the propagator G(x, y).This approach is known as the CJT method.In this generalized version, the effective action becomes the generating functional of the two-particle irreducible (2PI) vacuum graphs rather than the conventional 1PI diagrams.The zero temperature version of CJT method has been generalized to finite temperatures by Amelino-Camelia and Pi in [58,59].According to [59], the CJT method is equivalent to summing up the infinite class of "daisy" and "super daisy" graphs and is thus useful in studying such strongly coupled models beyond the mean-field approximation.
Below, we will use the imaginary time Matsubara formalism.We evaluate the momentum space integrals with the time component k 4 replaced by a summation over discrete frequencies because the temporal coordinate is compactified on a circle.Here, β is the inverse of the temperature T .This leads to the following relationship: where f is a generic function and we have used k 4 = 2πinT in the last step.
According to the CJT formalism [59], the finite temperature effective potential with a generic scalar field ϕ is given by: where the sum runs over all meson species and V 0 , D −1 (ϕ; k), V 2 (ϕ, G) denote, respectively, the tree level potential, tree level propagator as well as the infinite sum of the two-particle irreducible vacuum graphs.In this work, we use the Hartree approximation according to which V 2 (ϕ, G) is simplified to a one "double bubble" diagram.In the simplest one-meson case, it is proportional We therefore obtain a gap equation by minimizing the above effective potential with respect to the dressed propagator G i (ϕ; k): where δV 2 (ϕ, G)/δG i (ϕ; k) denotes the self energy operator.It is convenient to introduce the full field-dependent propagator: where M ≡ M (ϕ, k) is the effective mass which can be interpreted as tree-level mass dressed by the tadpole contributions.By using the full propagator (15), the gap equation ( 14) yields a set of equations leading to effective temperature-dependent masses (the detailed expressions for general N f are given in App.A).We can easily check that under the Hartree-Fock approximation the effective mass M becomes momentum-independent.It is because in (14), the momentum k in the dressed propagator G and tree level propagator D cancel each other while δV2(ϕ,G) δGi(ϕ;k) in the Hartree-Fock approximation leads to a term proportional to β G (ϕ; k) where the momentum k is integrated out.
Our final results for the CJT improved finitetemperature effective potential of the LSM sector can be written as: with where Here, J B (R 2 ) comes from the logarithmic integral term 13) while I B (R 2 ) comes from the β G (ϕ; k) term.Also, we have implemented the definitions r i ≡ m i (σ)/T and R i ≡ M i (σ, T )/T where m i (σ) are the tree-level meson masses while M i (σ, T ) are the thermal masses (or dressed masses) defined in (15).We list the tree-level meson masses below: In the chiral limit, i.e., m π = 0, (19) reduces to

C. Model Parameters and Observables
In this section, we summarize all the key model parameters and observables in this work.
The Polyakov-loop part of the model consists of the Polyakov-loop potential and the interaction term of the Polyakov loop with the LSM.The free parameters in the Polyakov-loop potential (3) have been fitted in [35] to lattice data from [85] for the pure Yang-Mills theory.The introduction of quarks affects these coefficients.Here, we consider only the impact generated from the medium potential (7).We expect the medium potential to capture the main features of the fully non-perturbative results.The medium potential describes the interaction between the Polyakov loop and the LSM and introduces one more free parameter, i.e. the coupling g between the scalar field and the fermions.
For the purely LSM part of the model, for fixed N f , we have four input parameters m 2 , λ σ , λ a , c which we fix via the four observables (f π , m σ , m η , m a ).We are not considering the pion mass since we are working in the chiral limit and therefore the pion mass vanishes, m π = 0.It should be noted here that the above observables are defined at zero temperature, see ( 6) and (20).

III. PHASE TRANSITION AND BUBBLE NUCLEATION A. Bubble nucleation
In the case of a first-order phase transition, the latter occurs via the nucleation and expansion of vacuum bubbles.Its dynamics is captured by the bubble nucleation rate.One starts with a computation of the temperaturedependent tunnelling rate between the meta-stable and stable vacuum related to the three-dimensional Euclidean action S 3 (T ) [86][87][88][89] via 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 scalar field Φ (see (5)).Once we choose a vacuum configuration with VEV along the σ direction, we can focus on the effective potential only depending on ℓ and σ.In the cases considered here (i.e.N f ≥ 3), we always have a first-order chiral phase transition and the Polyakov loop acts as a spectator.Since we are investigating matter in the fundamental representation of the underling gauge theory there is no proper confinement transition justifying why we do not investigate the latter transition [81].In this respect, the Polyakov loop works indeed as a spectator field for the chiral transition, even though its properties, such as the spatial two-point function, will be affected by the chiral transition itself [81].
The chiral phase transition is described by the scalar field σ, representing the chiral condensate.When the Polyakov loop is included, we use a mean-field approximation in the Polyakov loop ℓ.This means that we evaluate the Polyakov loop at the minimum of the effective potential for given values of σ and T .Thus, the potential becomes a function of only the sigma field, (σ, T ), V eff (σ, T ) = V eff (σ, T, ℓ min (σ, T )).This is a good approximation as long as the minimum value of the Polyakov loop does not strongly depend on the sigma field, which is indeed the case in our computations.Based on the above argument, the three-dimensional Euclidean action simplifies to The bubble profile is obtained by solving the equation of motion of the action in (23) and is given by with the associated boundary conditions To solve ( 24) and ( 25), we have implemented the overshooting/undershooting method and employed the Python package CosmoTransitions [90].We substitute the solved bubble profile σ(r, T ) into the threedimensional Euclidean action (23) and, after integrating over r, S 3 depends only on T .
B. Gravitational-wave parameters

Inverse duration time
We start with the decay rate of the false vacuum to the true vacuum due to thermal effects (while the decay rate due to quantum corrections is strongly suppressed).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.One of the three key parameters in the GW spectrum, the inverse duration time β is defined as: It is often convenient to introduce a dimensionless version β which is defined relative to the Hubble parameter H * at the characteristic time t * as where we used that dT /dt = −H(T )T .Note that here, for simplicity, we have assumed that the temperature in the hidden and visible sectors are the same, T d = T v .But in a more general case, these two temperatures can be different.Introducing a difference between the hidden and visible temperature can in some cases avoid phenomenological constraints [34,91].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 of being in the false vacuum is about 0.7.It is apparent that the percolation temperature corresponds to a point when the phase transition has proceeded further compared with that at the nucleation temperature, and thus T p ≤ T n .We use the percolation temperature throughout this work.To find it, we first write the false-vacuum probability as [92,93] with the weight function [94] The percolation temperature is defined by I(T p ) = 0.34, corresponding to P (T p ) = 0.7 [95].Using T * = T p in (28) yields the dimensionless inverse duration time.

Energy budget
For the second important parameter the phase transition strength α, there are two different definitions.In the literature, α is often defined as the latent heat during the phase transition per d.o.f.In this work, we define the strength parameter α by using the trace of the energymomentum 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).This definition quantifies the jump of the energy density and the pressure across the phase boundary weighted by the enthalpy and thus it measures the strength of the first-order phase transition.The enthalpy quantifies the total d.o.f.'s of the system which participates in the phase transition and 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 [96][97][98].With (32) and (33), α is given by 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 Sec.III C.

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 [99][100][101][102][103][104].In our case, we have a strongly-coupled system and most likely a fully nonperturbative analysis would be necessary to determine the friction.Some initial works towards this direction by using holographic methods can be found in [62,[105][106][107].
We display results for two different wall velocities, once the Chapman-Jouguet (CJ) detonation velocity, which is given by and once for v w = 0.1.Using the CJ velocity is an optimistic assumption since it leads to the largest efficiency factor for the production of GWs, see next section, and therefore also to the largest GW signal.The results with the CJ velocity can be interpreted as an upper bound for the true GW signal.Note however, that as long as the wall velocity is larger than the speed of sound, v w ≥ c s = 1/ √ 3, the wall velocity does not strongly impact 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 [108].

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 to the phase transitions considered here.The efficiency factor for the sound waves κ sw consists of the factor κ v [109] as well as an additional suppression due to the length of the sound-wave period τ sw [110][111][112] In our notation, τ sw is dimensionless and measured in units of the Hubble time.It is given by [112] τ sw = 1 − 1/ 1 + 2 (8π) and for β ≫ 1, τ sw can be simplified to It is apparent that for larger β, τ sw will be suppressed.
Ūf is the root-mean-square fluid velocity [110,113] Ū We follow [109] for κ v where it was numerically fitted to simulation results.The factor κ v depends on α and v w , and, for example, at the CJ velocity it reads For a phase transition with α ∼ O(10 −2 ), we have κ v ∼ 0.1.
The GW spectrum sourced by sound waves in the plasma is given by with the peak frequency f peak ≃ 1.9 • 10 −5 Hz g * 100 and the peak amplitude Here, h = H/(100km/s/Mpc) is the dimensionless Hubble parameter and g * is the effective number of relativistic d.o.f.'s, including the SM d.o.f.'s g * ,SM = 106.75 and the dark sector ones, which is g * ,dark = 16 + 2N 2 f , where the 16 comes from the SU (3) gluon d.o.f.'s while 2N 2  f comes from the meson field Φ (see (5) where N 2 − 1 d.o.f.'s, respectively, for a field and π field and 2 extra d.o.f.'s for σ and η ′ ).
The factor Ω 2 dark in (43) accounts for the dilution of the GWs by the visible SM matter which does not participate in the phase transition.The factor reads It is also possible to absorb the above dilution factor into the redefinition of the strength parameter (see e.g.[35,91]).The decisive quantity that determines the detectability of a GW signal is the signal-to-noise ratio (SNR) [135,136] given by Here, h 2 Ω GW is the GW spectrum given by (41), h 2 Ω det is the sensitivity curve of the detector, and T is the observation time, for which we assume T = 3 years.We compute the SNR of the GW signals for the future GW observatories LISA [73][74][75], BBO [65][66][67][68][69], and DECIGO [69][70][71][72].The sensitivity curves of these detectors are nicely summarised and provided in [77].

IV. RESULTS
We are now ready to present the results starting with the LSM as a benchmark model for N f = 3, 4, and 5.For N f = 3, we compare the PLSM and the LSM while for N f = 5 we use LSM alone but consider separately the cases in which the determinant operator is turned on or off in (4).
For all these models, we randomly choose values for the meson masses and the pion decay constant.Afterwards, the ratios are adjusted so the phase transition temperature is at T c = 1 GeV.The values of the mass parameters are chosen from the range  Note that we are working in the chiral limit and thus m π = 0.After these values are chosen, we relate them to the Lagrangian parameters of the theory (m 2 , λ σ , λ a , c) using (20).We then perform some sanity checks on these theories: firstly, we demand that the vacuum is stable.The vacuum stability condition reads [137] Secondly, we demand that the theory has a broken chiral symmetry at zero temperature.Lastly, we limit the magnitude of the couplings to be, in absolute value, within the range λ a , λ σ < 16π.If these demands are not met we discard the given theory point.Around 60k values are randomly chosen for each model in the ranges (46), and all relevant phase transition parameters are computed.In this manner, we obtain plots such as the left panel in Fig. 1, where the phase transition parameters α and β are displayed for the LSM with N f = 3, 4, 5.Each point in the left panel of Fig. 1 represents a theory obtained for a given choice of input parameters.Due to the large amount of randomly chosen parameter values, we obtain a thorough picture of the entire (α, β) LSM landscape.Note that for readability reasons, we have reduced the number of points shown to three thousand per given number of flavours in the left panel of Fig. 1.

While the plots give us a thorough picture of what a
given underlying theory can achieve, they are not immediate to interpret.Therefore, we create more accessible plots by averaging over the points.For example, we sort the results in α and then group them in sets of roughly 2k points per set.We compute the mean and assume for the variance a Gaußian distribution above and below the mean.This leads to the plot in the right panel of    features.For example, the α − β-data displays a characteristic 'edge' around (α, β) = (0.033, 3 • 10 4 ) for N f = 3 and (α, β) = (0.048, 3 • 10 4 ) for N f = 4.This edge is also visible in the averaged data in Fig. 1.We made sure that the averaging process did not wash out important physics features in the results that we displayed.

A. Flavour dependence
In Fig. 1, we show the averaged plots of the β vs α plane for N f = 3, 4, 5.The strongest first-order phase transitions are in the bottom right corner of these plots, where α is large and β is small.We find that on average the first-order phase transitions are stronger for larger N f due to the larger α values (the curves are shifted towards the right with increasing N f ).Note however, that the strongest first-order phase transitions occur in the N f = 3 case: there are a few sparse blue points in the left panel of Fig. 1 with β ∼ 10 3 .These happen due to a zero-temperature barrier in the potential which is induced from the determinantal term going ∼ σ 3 for N f = 3.Note that this cannot happen for other values of N f .To support the statement of Fig. 1, we show in Fig. 2 the GW parameters α and β separately as a function of the sigma meson mass m σ .The sigma meson mass is the best parameter to plot against as we will show below.In Fig. 2, we can clearly see that α gets larger (left panel) and β gets smaller (right panel) with increasing N f , both leading to a stronger GW signal.

B. Dependence on physical and theory parameters
We now want to understand what are the dependencies on the physical parameters (masses and pion decay constant) as well as on the theoretical model parameters.For this, we compute the peak amplitude, (43), with the assumption that wall velocity is given by the Chapman-Jouguet detonation velocity (40).This is an optimistic assumption as it maximises the peak amplitude and most likely the strong gluon dynamics lead to significant friction and a slower wall velocity [31,32,[105][106][107].Nonetheless, we use this as a benchmark and discuss the impact of the wall velocity later in Sec.IV E. Overall, the peak amplitude gives us the most direct re-   lation to the strength of the first-order phase transition, and therefore we plot it as a function of the physical and theory parameters.In Fig. 3, we show the averaged plots of the peak amplitude Ω peak as a function of the pion decay constant f π and the masses m σ , m η , m a .The latter are all measured in units of T c .We observe that the strongest correlation are between Ω peak and f π as well as m σ : Ω peak increases with increasing f π and with decreasing m σ .In particular, the width of the curves is the smallest for m σ indicating that the distribution of theories is most clearly sorted with this parameter.Note, that for N f = 3 and small m σ , the largest GW amplitudes are generated, reaching up to Ω peak = O(10 −14 ), which would be almost in reach of the sensitivity of LISA.Due to the clear correlation between Ω peak and m σ , we will use m σ as the preferred plotting parameter in the upcoming sections.
On the other hand, we only see a mild dependence of the peak amplitude on the masses m η and m a .In particular, there is almost no dependence on m η , while we observe a mild increase of Ω peak with m a .The N f = 3 curves are oscillating strongly as a function of m η and especially as a function of m a .The reason is simply that the curves have not converged yet despite us having in-cluded 60k theory points.This is due to the very wide distributions as a function of these parameters.
In Fig. 4, we show the averaged plots of the peak amplitude Ω peak as a function of the couplings m 2 , λ σ , λ a , and c.The parameter m 2 is measured in units of N f T c and c is measured in units of T c if it is dimensionful.
We observe a very similar picture as with the physical parameters: we have two parameters that display a strong correlation (m 2 and λ σ ) and two parameters that display almost no correlation (λ a and c).The peak amplitude increases with a decreasing m 2 and a decreasing λ σ .The strongest is between Ω peak and λ σ , which is in straight analogy to the correlation with m σ .This can be understood directly at the hand of the relation between λ σ and m σ , as we will discuss in Sec.V.
Overall, our observations are consistent with the results discussed in Sec.IV A since the GW peak amplitude generically increases with N f , which is most easily visible in the top two panels of both, Figs. 3 and 4.

C. PLSM vs LSM
In this section, we investigate the effect of the Polyakov loop on the GW spectrum.This implies that we include the Polyakov loop potential (3) and the medium potential (7) to the LSM.This results in one new free parameter, the coupling g in the constituent quark mass (9), which we set to g = 1 for this investigation.We also restrict ourselves to N f = 3 but we expect that the inclusion of the Polyakov loop has similar overall effects for all N f .We note that the inclusion of the Polyakov loop increases the complexity of the computation and therefore we use fewer statistics for the PLSM (14k points versus 60k points), which results in a less converges and less monotonic plot for the PLSM.
In Fig. 5, we show the GW parameters α and β as well as the peak amplitude Ω peak as a function of the sigmameson mass m σ in the PLSM and LSM for N f = 3.We observe that the Polyakov loop has little effect on the parameter β (central panel) while it increases the strength parameter α compared to the LSM (left panel).In summary, this leads to an increased GW peak am-plitude in the PLSM (right panel).Note, however, that the strongest GW signals are still present for the LSM at very small values of the sigma meson mass.

D. Determinant operator
In this section, we investigate the influence of the determinant operator.This is in particular relevant for larger N f since the mass dimension of the determinant term is [c] = 4 − N f , which means that it is a relevant operator for N f = 3, marginal for N f = 4, and non-renormalisable for N f = 5.Neglecting the determinant term leads to a vanishing mass for the eta-meson, m η = 0. To ensure a fair comparison between the theories, we included the determinant term for N f = 5 in the previous sections.Since for N f = 5, the term is non-renormalisable this leads to an un-or meta-stable tree-level potential and we interpret the theory to be valid below a given cutoff below.The instabilities of the potential are above this given cutoff scale.We discard theories that have a cutoff scale which is too close to the critical temperature.To Figure 7. Signal-to-noise ratio at BBO as a function of the sigma-meson mass for the optimistic case that the wall velocity is vw = vJ (α) (left panel) and for the wall velocity vw = 0.1 (central panel).We choose the critical temperature of the phase transition such that the peak frequency falls on the maximal sensitivity of the BBO detector.The critical temperature that maximises the signal-to-noise ratio is shown in the right panel.
ensure that this is a valid treatment, we compare in this section the results to the case of a vanishing determinant term.
In Fig. 6, we compare the LSM at N f = 5 with and without the determinant term.We show the GW parameters α (left panel) and β (central panel) as well as the peak amplitude Ω peak (right panel) as a function of the sigma-meson mass m σ .We observe that on average curves agree very well with each other, and we conclude that the determinant term does not influence the strength of the GW signal.This justifies our treatment in the previous sections and also agrees with our observations in Figs. 3 and 4 where we have observed that the GW peak amplitude neither depends on m η nor on c.This also holds true for N f = 3 and 4. The latter is indeed surprising since the determinant term is a relevant operator for N f = 3 and it is also responsible for the strongest first-order phase transitions due to a zero temperature barrier in the potential, as shown in Fig. 3 for small m σ and as discussed in Sec.IV A.

E. Signal-to-noise ratio
The signal-to-noise ratio (45) is the key quantity to decide whether a GW signal is detectable or not.We focus on BBO since it has the greatest sensitivity of the detectors and assume that the signal is detectable for SNR > 1.We further choose the critical temperature of the phase transition such that the peak frequency falls on the maximal sensitivity of the detector.This implies for BBO a critical temperature of the order of T c = O(100) GeV for v w = v J and T c = O(10) GeV for v w = 0.1, which we show in the right panel of Fig. 7.
Here we already indicated that the result depends on the terminal wall velocity, which we use as an external parameter since it is difficult to compute in a nonperturbative setting.There have been some indications that the wall velocity is rather small in these kinds of phase transitions since the strong gluon dynamics could provide a lot of friction that slows down the bubble wall.Here, we exemplify our results for two different wall velocities: the Chapman-Jouguet detonation velocity (40), which leads to the strongest GW signal and is an optimistic best-case scenario, and secondly v w = 0.1.
The results are displayed in Fig. 7.In the left panel, we show the SNR for N f = 3, 4, 5 as a function of the sigmameson mass in the best-case scenario with v w = v J .Here, we find a detectable GW signal for m σ /T c < 1.96 (N f = 3), m σ /T c < 2.65 (N f = 4), and m σ /T c < 4.19 (N f = 5).For smaller wall velocities, v w = 0.1, see central panel, none of the GW signals are detectable, with the exception of the very strong first-order phase transitions in the N f = 3 case for small m σ .

V. DISCUSSIONS: CONDITIONS FOR A STRONGER FIRST-ORDER PHASE TRANSITION
In the previous section, we have elucidated the phase structure of the (P)LSM and have identified corners of the theories where strong first-order phase transitions take place, some measurable at future GW detectors.We found the strongest first-order phase transitions for small sigma masses m σ , large pion decay constants f π , small m 2 , and small λ σ .Now we want to discuss how these conditions are interrelated.
Naively, a smaller λ σ leads to a larger m σ , which is contradictory to our above results.This is most easily seen at the hand of the tree-level expression, see (19), and utilising the minimum condition of the tree-level effective potential.For N f = 3, we find which shows that a smaller λ σ leads to a larger m σ .Similar expressions hold for larger N f .Opposing this naive view, our full data shows a different correlation: we chose the physical parameters ran-  domly in the ranges given in (46) from which we can directly compute the couplings of the theory.In this data set, we observe the correlation that larger λ σ leads to a larger m σ .
The reason that the naive expectation from equation (48) fails, is that it is based on an analysis for fixed values of m 2 and c.For a decreasing λ σ at fixed m 2 and c, the other masses (m a and m η ) would increase as well and leave the range given in (46).
It is more useful to express (48) in terms of other physical observables.For example, the equation can be recast in the form (N f = 3) from which it becomes apparent that small m 2 corresponds to small m 2 σ for fixed m 2 η .The latter is a valid assumption since m η is dominantly determined by c, see (46).Similarly, the relation (N f = 3) can be misleading, since it does not imply a larger f π for smaller c.Instead, a smaller c corresponds to a smaller m 2 η , which together keeps f π constant.More straightforwardly, a large f π corresponds to a large vev, since Together with the minimum condition ( which entails that a large v corresponds to small λ σ , we arrive at the desired conclusion that small λ σ corresponds to large f π .We infer the last key relation from (20), which reads For fixed m 2 η , which is justified as explained before, m 2 σ is directly related to λ σ v 2 .From ( 52), we know that a large v corresponds to small λ σ but the explicit λ σ dependence in ( 53) is stronger and overall, we arrive at correlation that a larger λ σ corresponds to a larger m σ , contrary to the naive statement in (48).
In summary, we have the following correlations: • Smaller m leads to smaller m σ and thus stronger GW signals, see the top left panel of Fig. 4.
• Smaller λ σ leads to smaller m σ and thus stronger GW signals.This is clearly shown in the top right panel of Fig. 4.
• Smaller c leads to smaller m σ and thus stronger GW signals.See lower right panel of Fig. 4.
These are the observations on the physical and theoretical parameters.From the physical nature, we observe the following interesting relations • Weak coupling regime: Our model allows us to study the transition from an effective strong to an effective weak coupling regime.The strongest firstorder phase transitions take place for small λ σ and m 2 couplings, where the latter is measured in units of T c , see Fig. 3. Note, however, that this does not hold for the other couplings, λ a and c.This supports our previous findings [35,39,51] where we observed that strongly coupled systems have a strongly suppressed GW signal.
• Light sigma mass: In the limit of a small sigma meson mass, the state shares some similarities with the dilaton often associated with the dynamics of a near-conformal field theory.Near conformal theories are known to have a very strong first-order phase transition [111,[138][139][140][141][142][143][144], which we also observe here at small sigma masses.
• Supercooling: The dominant reason for a strong first-order phase transition is typically a large temperature difference between the critical and the phase-transition temperature, also known as supercooling.We analysed this relation in Fig. 8 where we show the supercooling parameter 1 − T p /T c as a function of the sigma mass m σ (left panel), the GW parameters β (central panel) and α (right panel).We clearly observe the stronger phase transition with an increasing supercooling parameter.Note, the sharp edges for N f = 3, 4 in the α plot (right panel), which are related to the edges already observed in Fig. 1.These edges limit the maximal strength of the phase transition since α increases very slowly despite a strongly increasing supercooling parameter.The exception is N f = 5 where α is continuously increasing with the supercooling parameter.

VI. CONCLUSIONS
In this work, we explored the possibility that a strongly coupled dark sector may produce an observable stochastic gravitational wave signal.We employed the Polyakov Linear Sigma Model as an effective theory to study the first-order chiral phase transitions with N f = 3, 4, 5 flavours.We further implemented the well-established Cornwall-Jackiw-Tomboulis (known as CJT) method to bridge perturbative and non-perturbative regimes of the theory and thus are able to study both strongly-and weakly-coupled systems at once.
We observed that stronger first-order phase transitions and thus larger GW signals generically occur when the system features a light sigma meson m σ and/or is weakly coupled, corresponding to larger f π .The strongest phase transitions with β = O(10 3 ) and α = O(10 −1 ) are present for N f = 3 and small sigma mass.This is due to a zero temperature barrier in the potential induced by the cubic determinant term.It provides GW signals which are just outside the sensitivity of LISA, but that can be observable via the BBO and DECIGO detectors assuming that the confinement temperature is such that the peak frequency falls into the maximal sensitivity range of the detectors.We have also observed that the strength of the first-order phase transition generically increases with N f due to an increase of the latent heat (strength parameter α).
The main features of our results are expected to be sufficiently general, robust and largely independent of the details of the model computations.This is so since all the physical underlying mechanisms at play have been embedded in the effective model description.
Furthermore, we introduce the tensors F, G, and H describing the tensor structures in the potential, see (28) of [145].We also use S for the full quantum propagator of the scalar particles (σ and a) and P for the full quantum propagator of the pseudoscalar particles (π and η).With these definitions, the relevant structures for the thermal masses are given by where R i ≡ M i (σ, T )/T and the I B is given by (18).The expressions are in straight analogy for the contractions with the pseudoscalars, P. We furthermore have two contractions that are proportional to the determinant term, which are only relevant for the N f = 4 case With these results and returning to our standard coupling conventions with (A1), the thermal masses are given by These thermal masses are used in (17) and thereby provide the finite temperature part of the potential.

)N f = 3 N f = 4 N f = 5 Figure 1 .
Figure 1.We show the range of α and β values of the LSM for N f = 3, 4, 5.In the left panel, we show the actual distribution of theory points, while in the right panel, we display the averaged values.On average, the LSM produces stronger GW signals with increasing N f due to the larger α values.Nonetheless, the strongest GW signals are achieved with the LSM for N f = 3, corresponding to the sparse blue dots at small β in the left panel.

Figure 2 .
Figure 2. We show the averaged values of α (left panel) and β (right panel) in the LSM for N f = 3, 4, 5 as a function of the sigma meson mass in units of Tc.Smaller values of mσ correspond to larger values of α and smaller values of β leading to a stronger first-order phase transition.

Fig. 1 .
While this process washes out some features of the full data it still displays the important qualitative

Figure 3 .
Figure 3.We show the averaged values of the peak amplitude Ω peak as a function of physical observables fπ (top left), mσ (top right), mη (bottom left), and ma (bottom right) all in units of Tc in the LSM for N f = 3, 4, 5.The pion decay constant and the sigma meson mass have the strongest correlation with the peak amplitude: larger values of fπ and smaller values of mσ lead to a larger Ω peak .

Figure 4 .
Figure 4. We show the averaged values of the peak amplitude Ω peak as a function of the model parameters m 2 (top left), λσ (top right), λa (bottom left), and c (bottom right) in the LSM for N f = 3, 4, 5.The strongest correlation is found between Ω peak and m 2 as well as λσ: smaller values of m 2 and λσ lead to a larger Ω peak .

Figure 5 .Figure 6 .
Figure 5.We show the impact of the Polyakov loop on the GW parameters α (left panel), β (central panel), and on the peak amplitude Ω peak as a function of the sigma meson mass mσ at the example of N f = 3 and the coupling g = 1 in the constituent quark mass,(9).While the Polyakov loop has no effect on β, it increases the strength of the phase transition by increasing the value of α.

5 Figure 8 .
Figure 8. Supercooling parameter 1 − Tp/Tc as a function of the sigma mass mσ in units of Tc (left panel), and the GW parameters β (central panel) and α (right panel).