Observational prospects for gravitational waves from hidden or dark chiral phase transitions

We study the gravitational wave (GW) signature of first-order chiral phase transitions ($\chi$PT) in strongly interacting hidden or dark sectors. We do so using several effective models in order to reliably capture the relevant non-perturbative dynamics. This approach allows us to explicitly calculate key quantities characterizing the $\chi$PT without having to resort to rough estimates. Most importantly, we find that the transition's inverse duration $\beta$ normalized to the Hubble parameter $H$ is at least two orders of magnitude larger than typically assumed in comparable scenarios, namely $\beta/H\gtrsim\mathcal{O}(10^4)$. The obtained GW spectra then suggest that signals from hidden $\chi$PTs occurring at around 100 MeV can be in reach of LISA, while DECIGO and BBO may detect a stochastic GW background associated with transitions between roughly 1 GeV and 10 TeV. Signatures of transitions at higher temperatures are found to be outside the range of any currently proposed experiment. Even though predictions from different effective models are qualitatively similar, we find that they may vary considerably from a quantitative point of view, which highlights the need for true first-principle calculations such as lattice simulations.


I. INTRODUCTION
Since their first direct detection on Earth [1][2][3], gravitational waves (GW) provide a new and complementary approach towards observational astrophysics. Not only do they allow to study violent transient phenomena like black hole or neutron star mergers, but they also offer a unique way to probe the physics of the early Universe. In particular, it is well known that first-order cosmic phase transitions (PT) can create a stochastic background of gravitational waves [4][5][6][7][8][9]. Even though the Standard Model (SM) of particle physics does not feature such a PT [10][11][12][13][14], they are predicted by a plethora of beyond-the-SM (BSM) theories. For instance, successful electroweak baryogenesis typically requires a strong first-order electroweak phase transition, which can e.g. be obtained in models with an augmented scalar particle content [15]. Of course, gravitational waves may also be produced during phase transitions in dark or -more generally -hidden sectors. Correspondingly, studies investigating gravitational wave signatures of a wide range of new physics scenarios received a lot of attention in the recent years, see e.g. [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].
In the present article, we focus on BSM theories with a new confining gauge force SU(n c ) acting on n f flavors of hidden fermions in the fundamental representation. Similar to the situation in quantum chromodynamics (QCD), strong (i.e. non-perturbative) dynamics are anticipated to spontaneously break the hidden counterpart of chiral symmetry in such models. Unlike in realworld QCD, 1 the corresponding chiral phase transition can, however, be strong and of first order. It is the gravitational wave background associated with such a transition that we will be interested in, cf. also [16].
Hidden sectors of the above described type are particularly attractive for several reasons. They are, for instance, featured in various models of dark matter (DM) like strongly interacting massive particle (SIMP) DM [39,40], or composite DM as discussed in e.g. Refs. [41][42][43][44][45][46][47]. Similarly, they are introduced in Twin Higgs models [48,49], Hidden Valley theories [50] or vector-like confinement models [51]. Lastly, strongly coupled hidden sectors commonly occur in theories based on classical scale invariance, since they offer one of the two 2 known possibilities to dynamically generate a scale in a previously scaleless setup. Several realistic classically conformal SM extensions along those lines were suggested e.g. in Refs. [53][54][55][56][57][58]. Reflecting this rich variety of applications, well-motivated scenarios exist at vastly different scales ranging from 100 MeV to above 10 TeV, see Ref. [16] for an overview.
Given the presented appealing physics case for BSM theories with new confining gauge forces, it is desirable to quantify gravitational wave signatures from a potential hidden chiral phase transition, in order to assess the prospects for their observation. Due to the non-perturbative nature of the problem, the corresponding analyses generally prove difficult. Accordingly, no true first-principle calculations exist yet. There are, however, various possibilities how to still proceed. For a start, one may simply estimate the relevant quantities based on certain assumptions and/or general arguments. For instance, the inverse duration of the phase transition β, which is a key parameter in understanding the associated GW spectrum, is commonly estimated by [5,6] β where H is the Hubble parameter at the time of the transition, M Pl is the (reduced) Planck mass and the second relation holds for a wide range of transition temperatures T n . Analyses of GW signals based on the estimate of Eq. (1) can e.g. be found in Refs. [16,59,60].
Here, we will follow a different approach. Namely, we will model the dynamics in the strongly coupled hidden sector by employing various effective theories, which are known to work well in capturing non-perturbative aspects of the low-energy regime of real-world QCD. Thus we will be able to explicitly calculate relevant quantities characterizing the hidden chiral phase transition such as the bubble nucleation rate, the transition temperature or β/H. Interestingly, we find that all our calculations result in values for β/H which are much larger as compared to the rough estimate of Eq. (1), namely β/H O(10 4 ) . ( Since the energy density of the GW background decreases -depending on the nature of its sourcelinearly or even quadratically with β/H, and since the spectrum's peak frequency increases linearly with β/H (see e.g. [61]), conclusions about the observability of the GW signal change considerably, as we will demonstrate in the main text. The article is structured as follows. In Section II, we briefly introduce the low-energy effective models that we employ throughout the paper. Those are the Nambu-Jona-Lasinio (NJL) model [62,63] and its Polyakov-loop enhanced variant (PNJL) [64], as well as the linear sigma model (LSM) [65]. For definiteness, we will consistently focus on the case of QCD-like hidden sectors with n c = 3 colors and n f = 3 hidden fermion flavors, which we assume to be massless for simplicity. In Sections III and IV we review important aspects of first-order phase transitions and of the associated GW background spectra, respectively. Section V then contains this work's main results. Namely, we use the aforementioned effective models to explicitly calculate different parameters crucial in characterizing the hidden PT. Furthermore, we determine the corresponding GW signals and evaluate the prospects for their detectability with future observatories. We present our findings in a way that they can be applied to hidden or dark sectors over a wide range of inherent scales, i.e. transition temperatures. Finally, we conclude in Section VI.

II. LOW-ENERGY EFFECTIVE MODELS OF QCD-LIKE THEORIES
In the present work we are interested in strongly interacting hidden sectors which we assume to be describable by a theory similar to QCD. In particular, the hidden analogue of the chiral phase transition (χPT) will be triggered by non-perturbative effects, requiring a treatment beyond standard perturbation theory. To that end, we here employ different effective models for the purpose of reliably capturing the phase transition's dynamics. In order for such an approach to be successful, all these models necessarily have to share the relevant symmetry properties of the underlying hidden sector, which we briefly outline in the following.
The matter part of a QCD-like theory with n f fermion flavors q transforming in the fundamental representation of a hidden SU(n c ) gauge group reads where i is a flavor index. Color indices ranging from one to n c = 3 are suppressed for readability, but are to be understood as being summed over as well. In the chiral limit of vanishing fermion masses, m i = 0, the above Lagrangian classically exhibits a global chiral symmetry G := U(n f ) L × U(n f ) R . Separating the transformations which treat the chiral components q L and q R alike (vector subgroup) from those which transform them differently (axial-vector subgroup), as well as using the Lie algebra isomorphism U(N ) ∼ = SU(N ) × U(1), the group G can be conveniently written as ]   A  72  0  248  458  491   B  90  0  400  672  697   C  74  0  291  328  431   D  108  0  694  535  792   TABLE I. Benchmark points used throughout the paper. The model parameters which correspond to the shown mass spectra and pion decay constant are compiled in Appendix D.
It is, however, well known that the axial U(1) factor is explicitly broken by quantum effects due to the Adler-Bell-Jackiw anomaly [66,67]. The symmetry group of the quantum action is therefore Nevertheless, the theory's vacuum does not respect the full group G, which is spontaneously broken down to SU(n f ) V × U(1) V by a finite expectation value of the fermion condensate, q i q i = 0. However, finite-temperature effects are anticipated to restore the chiral symmetry group G at sufficiently high temperatures [68][69][70]. The properties of the associated chiral phase transition depend both on the number n f of fermion flavors and on their masses as summarized in the famous Columbia plot [71]. Since we are interested in the gravitational wave signatures from a first-order hidden χPT, we will focus on the scenario of n f = 3 massless flavors in the following. This choice corresponds to the minimal number of flavors needed to obtain a first-order χPT [71,72].
Apart from the correct symmetry structure, a successful low-energy effective model must, of course, also incorporate all degrees of freedom most relevant to the physics process of interest. In the case of the chiral phase transition in QCD-like theories, those will primarily be the hidden counterparts of the well-known mesons from real-world QCD, see e.g. [73]. In particular, the presence of a σ-like scalar flavor singlet is crucial since it carries the same quantum numbers as the aforementioned fermion condensate q i q i . Thus, its acquiring a finite vacuum expectation value (vev) v σ triggers the spontaneous breakdown of chiral symmetry as indicated above. In the course of the present work we will additionally take into account the pseudoscalar octet π, the pseudoscalar singlet η and the scalar octet a. Depending on the concrete model, further degrees of freedom may be taken into consideration.
The remainder of the present section briefly introduces the three effective theories considered throughout this article. In order to meaningfully compare predictions about gravitational waves as derived from the individual models, we choose their free parameters such that all models produce an identical meson spectrum, i.e. the same σ, η and a masses, as well as the same pion decay constant f π . Since we exclusively discuss the chiral limit, the pion octet is an exact Nambu-Goldstone boson and thus massless. The benchmark meson spectra used throughout the paper are listed in Table I. A. The Nambu-Jona-Lasinio (NJL) model The NJL model was originally proposed to describe the nucleus [62,63]. Nowadays it is primarily used as a low-energy effective model of strongly coupled theories like QCD. The NJL model only describes the quarks, its Lagrangian does not contain any of the mesons which are the physical degrees of freedom at low energies. As we will explicitly see the Lagrangian contains four-and six-fermion interactions which render the model perturbatively non-renormalizable. The theory therefore necessarily contains a UV cutoff Λ. Following Ref. [55], we employ a four-dimensional momentum cutoff scheme and fix Λ = 0.93 GeV throughout the paper. An extensive review of the NJL model is, for example, given in Ref. [74].
The Lagrangian of the three-flavor NJL model in the chiral limit is given by [55] where the q i are Dirac fields corresponding to the (hidden) quarks, and Ψ is a fermion bi-linear with flavor indices i, j running from 1 to n f = 3. Spinor and color indices are implicit. As we are only interested in the chiral limit, the above Lagrangian does not contain an explicit quark mass term. Importantly, the NJL model thus exhibits the global symmetries of massless QCD, see e.g. [75]. To be more precise, the first two terms are invariant under G = U(3) L × U(3) R , whereas the last term is the 't Hooft determinant which mimics the anomalous U(1) A breaking mentioned below Eq. (4) [67,76]. One can explicitly demonstrate that these symmetries are indeed present by realizing that the Dirac spinors q transform under U(1) V and U(1) A as q → e iϑ q and q → e iϑ 5 γ 5 q , respectively, for ϑ, ϑ 5 ∈ R. Under SU(3) L and SU(3) R the chiral components transform as respectively, where q = q L + q R . The operators U L and U R correspond to unitary matrices with determinant equal to one. The above transformations clearly show the G invariance of the kinetic term. Establishing the symmetry properties of the other two terms is slightly more involved. As a start, note that the bi-linear Ψ as defined in Eq. (6) can also be rewritten in terms of the chiral fermion fields, namely Ψ ij =q Rj q Li . Employing Eqs. (7) and (8), one then finds that Ψ transforms as follows under the (global) symmetries of massless QCD: These transformations show that also the second term of the Lagrangian is invariant under the full group G . The last term, the 't Hooft determinant, is trivially U(1) V invariant. It also exhibits an SU(3) L × SU(3) R symmetry, because the SU(3) operators U L and U R have unit determinant. In contrast, Eq. (9) reveals that Ψ transforms non-trivially under U(1) A , which is therefore broken by the 't Hooft determinant. Thus, the NJL Lagrangian of Eq. (6) has all the symmetry properties we require for the NJL model to be a suitable effective theory describing the low-energy regime of QCD.
To be able to perform explicit calculations in the NJL model we use the (self-consistent) mean field approximation (MFA) [75,[77][78][79]. Within the MFA one defines the expectation value of the fermion bi-linear Ψ as a sum over effective meson fields, namely −4G Ψ = (σ + iη )1 + 2(a a + iπ a )T a .
The SU(3) generators in the fundamental representation T a are normalized so as to satisfy Tr(T a T b ) = δ ab /2 and fulfill the appropriate SU(3) algebra relations with the usual antisymmetric and symmetric SU(3) structure constants f abc and d abc . Correspondingly, the indices a, b and c run from 1 to n 2 f − 1 = 8. The meson fields introduced into the theory via Eq. (10) are to be understood as bosonic auxiliary fields meaning that they do not have tree-level kinetic terms and are thus classically non-propagating. However, fermionic quantum fluctuations will eventually induce kinetic terms and hence render the fields dynamical. More details regarding the MFA can be found in Appendix B. In particular, the final form of the mean field Lagrangian in terms of the fermion and effective meson fields is given in Eq. (B4).
In order to determine the dynamics of the chiral phase transition the model's effective potential is required. Note that only the σ field acquires a finite vev during the χPT. 3 The effective potential will therefore only be a function of the corresponding classical fieldσ, which is defined as the expectation value of σ in the presence of some external source. From Appendix B the tree-level potential V 0 can be read off. Apart from the tree-level terms the one-loop vacuum and thermal contributions, V CW and V FT , are required and are readily obtained by integrating out the fermions. Since all mesons are non-propagating at tree-level they do not contribute to the one-loop effective potential -only the fermions contribute. The finite-temperature effective potential of the NJL model at one loop can thus be written as with Using the finite-temperature effective potential one can determine the nature of the χPT by computing the potential's global minimumσ min (T ) for each temperature, where v σ =σ min (T = 0). It can be shown that for suitable G and G D the transition is indeed of first order. As we will see in Section III, studying the dynamics of the χPT necessitates the calculation of the probability that the scalar field which drives the transition, σ, tunnels from the theory's false to its true vacuum. Recall, however, that all meson fields are classically non-propagating in the NJL model. Since tunneling processes are inherently dynamical, tunneling of the σ field can only be described in the current setup, once we compute the σ field's kinetic term as it is induced via fermionic quantum fluctuations. At this stage of the approximation, the auxiliary field σ is thus promoted to a propagating quantum field. In doing so, the crucial quantity is the (finite-temperature) wave-function renormalization, which is defined as where Γ σσ (p 2 , σ) is the σ field's one-loop propagator, cf. Eq. (B7a). Note that in order to determine Z −1 σ (σ) this propagator is to be evaluated at the given field value σ, not necessarily at 3 The effective scalar meson fields are defined in terms of their constituent fermions as σ ∼ q1q and aa ∼ qTaq . The isospin symmetry SU(3) V present in the Lagrangian is not spontaneously broken, so that all mesons whose definition contains generators with off-diagonal elements, i.e. terms which couple fermions of different flavor, cannot acquire a finite vev. The only mesons whose vev can be non-zero are therefore σ, a3 and a8. Furthermore, we assume all quark masses to be equal. The vacuum should also respect this symmetry, so that ūu = d d = ss must hold. Since σ is the only meson proportional to the identity matrix, it is the only one which can obtain a finite vev. Note that pseudoscalars cannot acquire a non-zero vev because the vacuum is parity even.
v σ . Additionally, the loop integrals in the expression for the propagator quoted in Appendix B are generally evaluated at finite temperature when calculating Z −1 σ , and at zero temperature to compute the meson pole masses (see below). The full expression for Z −1 σ can be found in Eq. (B8). Eventually, we want to compare results about gravitational waves as obtained from the different considered effective models, provided that they predict the same meson mass spectra as well as the same pion decay constant f π . In the NJL model, the masses of the mesons -in the broken phaseare derived from the roots of their propagators, which, in turn, are determined from the mean field Lagrangian L MFA , see Appendix B. The one-loop meson propagators can be found in Eqs. (B7a) to (B7d). The pion decay constant on the other hand is given by (see e.g. [74,80]) with the constituent quark mass M c = M (v σ ), which corresponds to the effective fermion mass defined in Eq. (13d) evaluated at the global minimum of the effective potential at zero temperature.
The integral I 0 (v σ ) can be found in Appendix A. Finally, note that the NJL model only has three free parameters in the chiral limit, namely Λ, G and G D . Since we are always fixing Λ = 0.93 GeV only two free parameters remain. Therefore not all of the physical parameters m σ , m a , m η and f π are independent of each other. Using f π and m σ as an input to determine G and G D , the other meson masses (m a and m η ) become predictions of the model. The pion mass m π is zero since we are working in the chiral limit.

B. The Polyakov loop enhanced NJL model
Real QCD does not only exhibit a chiral PT, but also a (de)confinement PT, whose dynamics is governed by the gluons. Both of the aforementioned phase transitions are known to occur at similar temperatures from lattice QCD, see e.g. Ref. [81]. Therefore, gluon dynamics might also be relevant to the physics of the chiral phase transition. The NJL model discussed above does not include any gauge dynamics and is consequently unable to describe confinement. This issue was remedied in 2004 by Fukushima [64]. He incorporated gluon dynamics into effective QCD models at finite temperature by adding the Polyakov loop. This resulted in the so-called Polyakov-enhanced NJL (PNJL) model. A recent review of the Polyakov loop can, for instance, be found in Ref. [81]. The Lagrangian of the PNJL model (again in the mean field approximation) is given by where L MFA NJL can be found in Eq. (B4) and V glue is the background gluon potential in terms of the deconfinement transition's (pseudo) order parameter L, which denotes the expectation value of the traced Polyakov loop [81]. 4 Unfortunately, the gluon potential cannot yet be directly calculated from first principles. Still, the impact of including gauge dynamics can be effectively captured by appropriately parameterizing V glue . Data from, for example, lattice QCD can then be used to fix the values of the thus introduced additional free parameters. A discussion on different parameterizations can be found in Ref. [82], where it is also shown that the possible choices for the gluon potential are equivalent for temperatures below two or three times the critical temperature of the deconfinement transition. The dynamics of the χPT is hence not expected to depend strongly on the choice of potential. In this paper we use the gluon potential first introduced in Ref. [83], namely with The potential is chosen such that the value of L is constrained to lie between 0 and 1, with the deconfined and confined phases being reached for L → 1 and L → 0, respectively. Furthermore, since we are always working at zero chemical potential, L is real and thus L =L [81]. The critical temperature T glue of the QCD deconfinement PT in the pure gauge limit, i.e. assuming infinitely heavy quarks, is known to be T glue = 270 MeV from lattice simulations [84]. However, when taking into account three mass-degenerate quark degrees of freedom as in our model, T glue is reduced to 178 MeV [84]. Lastly, the dimensionless parameters a i and b 3 in Eq. (18) are determined from lattice QCD calculations as well [83], namely Including the Polyakov loop modifies the thermal part of the NJL model's effective potential from Eq. (13c) in two ways. First, the background gluon potential V glue of Eq. (17) is to be added. Second, the quark coupling to gluons must be taken into account. In contrast, the vacuum part of the effective potential is unchanged with respect to the NJL model. In particular, the meson propagators and thus also the meson masses are therefore identical to before. Combining all of the above we obtain as the final expression for the PNJL model's thermal one-loop effective potential with where r ≡ r(σ, T ) = M (σ)/T with M denoting the effective fermion mass of the NJL model as introduced in Eq. (13d). The MFA was applied to arrive at the quoted expression for V PNJL FT [81]. Using the above potential, we found that the deconfinement transition described by the Polyakov loop L is a crossover, while the χPT with order parameterσ min is of first order for all considered benchmark points. To proceed we again need to determine how the σ field 5 tunnels from the theory's false to its true vacuum. In doing so, we choose L to minimize the effective potential for any given point (σ, T ), i.e. L := L min (σ, T ). This will reduce V PNJL eff to be only a function ofσ and T , thus allowing us to use the standard formalism for tunneling in one field dimension, which will be outlined in Section III. Note that our choice of L is similar to the one made in Ref. [85].

C. The linear sigma model
Just as the Nambu-Jona-Lasinio (NJL) model the linear sigma model (LSM), as first introduced by Gell-Mann and Levy [65], is a phenomenological theory to effectively describe certain nonperturbative aspects of the low-energy dynamics of QCD or related strongly coupled theories. Unlike the NJL model the LSM is a renormalizable, purely scalar theory with the scalar fields being identified with the well-known mesons of QCD. There exist several variants of the LSM which differ by the global internal symmetry that they are based on. Since we are interested in a QCD-like theory with n f = 3 massless fermion flavors, we will here discuss the chiral limit of the U(3) × U(3) LSM, whose Lagrangian reads as follows: where the tree-level potential V tree is given by (see e.g. [86]) The complex-valued bosonic field Φ is a 3 × 3 matrix defined as and thus collects both the scalar meson fields (σ, a) and the pseudoscalar ones (η , π). The SU (3) generators T a were introduced below Eq. (10), the index a runs again from 1 to n 2 so that Φ ij has the same quantum numbers as the fermion bilinearq Rj q Li , where q are the matter fields introduced in Eq. (3). It is straightforward to derive from Eqs. (24) and (25) [67,76]. Spontaneous chiral symmetry breaking G → SU(3) V ×U(1) V is finally triggered by the scalar singlet σ acquiring a finite vacuum expectation value, i.e. σ → v σ + σ, or in terms of the Φ field from Eq. (24) where the latter relation is e.g. demonstrated in Ref. [86]. It is well-known that the LSM at finite temperatures is plagued by infrared divergences which lead to the breakdown of standard perturbation theory [87]. Many-body resummation schemes are therefore necessary for the purpose of obtaining robust results. In the following, we will employ the composite operator or 2PI formalism due to Cornwall, Jackiw and Tomboulis (CJT) [88] to implement the so-called Hartree-Fock approximation [89,90]. Details of the corresponding calculations are relegated to Appendix C and can e.g. also be found in Refs. [91][92][93][94].
Here, we only quote the final result for the U(3) × U(3) LSM finite-temperature effective potential, which reads whereσ denotes the expectation value of the σ field in the presence of an external source, such thatσ → v σ as the source vanishes. The individual components in Eq. (26) are given by Note that we here neglect vacuum contributions, which was previously demonstrated to not change the outcome qualitatively [92,93]. A few further comments regarding Eq. (27) are in order. First, the sum in Eq. (27b) runs over all meson species introduced in Eq. (24), i.e. i ∈ {σ, η , π, a}. The corresponding statistical weights g i are g σ = g η = 1 and g π = g a = n 2 Here, the m i (σ) denote the field-dependent tree-level meson masses [86], while the M i (σ, T ) constitute effective meson masses, which are self-consistently calculated within the CJT formalism for given T andσ as described in Appendix C. Lastly, the thermal integrals J B and I B are defined in Eq. (A4) of Appendix A. Formally, the first term in Eq. (27b) reproduces the one-loop term of the conventional perturbative series for the effective potential at finite temperatures, see e.g. [95]. However, the tree-level masses m i usually appearing as an argument of the J B are replaced by the effective masses M i . The second term goes beyond standard perturbation theory, as well. In total, it can be shown that the CJT effective potential in the Hartree-Fock approximation accounts for the resummation of all daisy and super-daisy diagrams [96]. Let us finally mention that for the U(3) × U(3) LSM in the chiral limit, i.e. m π = 0, the number of free Lagrangian parameters (m 2 , λ σ , λ a , c) equals that of observable vacuum quantities (f π , m σ , m η , m a ). Hence, evaluating Eq. (28) at the physical vacuum, i.e. forσ = 3/2 f π , yields a system of equations that can be uniquely solved for the former set of parameters given the pion decay constant and a full set of meson masses.
We will use the finite-temperature effective potential of Eq. (26) for all further calculations related to the linear sigma model.

III. CHIRAL PHASE TRANSITION
As we have mentioned in the beginning of Section II, the chiral phase transition (PT) in QCDlike theories with three (or more) massless fermion flavors is predicted to be of first order on general grounds [72]. For the low-energy effective models introduced in the previous section, we can explicitly demonstrate the existence of such a first-order chiral PT by using the appropriate finite-temperature effective potentials in Eqs. (12), (20) and (26), respectively. To that end, we examine how the potential's global minimumσ min , which, in turn, determines the given theory's true groundstate, changes with temperature. Fig. 1 shows the functionσ min (T ) as calculated within the investigated low-energy models for one of the benchmark points collected in Table I. The appearance of two degenerate minima at some critical temperature T c clearly signals the occurrence of a first-order phase transition. As an aside, note that the figure also explicitly demonstrates the equivalence of the NJL and the PNJL model at zero temperature.
From a physics point of view, it is well known that a first-order PT proceeds via the nucleation and subsequent growth of bubbles of the true vacuum inside an expanding universe, which is still in the metastable symmetric phase. Hence, the transition's properties crucially depend both on the nucleation rate Γ of the aforementioned bubbles, and on the Hubble parameter H. In what follows, we will briefly introduce these two quantities in turn.  Table I. The discontinuity at some critical temperature T c is characteristic of a first-order phase transition.
First, the temperature-dependent bubble nucleation rate Γ, which is also referred to as the decay rate of the false vacuum, can be calculated as [4,5,97] In the above equation, the theory's three-dimensional Euclidean action S 3 is to be understood as being evaluated for the O(3)-symmetric bounce or tunneling solutionσ b (r), i.e.
The field configurationσ b (r), in turn, is determined as the solution of the scalar background field's equation of motion subject to the boundary conditionsσ → 0 as r → ∞ and dσ/dr = 0 at r = 0, where r is the radial coordinate of three-dimensional space. If theσ field represents a fundamental degree of freedom, as is e.g. the case in the LSM, the theory's Euclidean action can be written as so that the background field's equation of motion reads In contrast, theσ field is composite in the (P)NJL model (i.e. non-propagating at tree-level). The Euclidean action therefore has to be slightly modified with respect to Eq. (30) and is now given by where Z σ ≡ Z σ (σ) is the wave-function renormalization introduced in Eq. (14). The background field's equation of motion is then found to be Throughout the present work, we solve Eqs. (31) and (33) using the CosmoTransitions package [98] and an appropriately customized version thereof, respectively. The same code is employed to compute the quantity S 3 [σ b (r)] based on Eqs. (30) and (32). Next, let us discuss the expansion of the universe as it is governed by the temperature-dependent Hubble parameter H. In the absence of significant supercooling, the universe will be radiation dominated during the whole phase transition so that H(T ) is given by where M Pl is the reduced Planck mass, M Pl = 2.435 · 10 21 MeV, while the effective number of relativistic degrees of freedom in the thermal plasma is denoted as g .
Given the bubble nucleation rate and the Hubble parameter of Eqs. (29) and (34), the nucleation temperature T n is defined as the temperature, for which one bubble per Hubble volume is produced on average, i.e.
Typically, the above integral is dominated by temperatures very close to T n . Hence, the defining condition in Eq. (35) is approximately equivalent to where we used Eqs. (29) and (34) to arrive at the second relation, ignoring the slowly varying factor (S 3 /(2π 2 T )) 3/2 in the expression for the nucleation rate Γ. Notably, the right-hand side of Eq. (36) depends on the temperature only logarithmically, so that the nucleation condition roughly translates to S 3 /T ≈ O(100) for a wide range of temperatures [4,5]. The gravitational wave spectrum associated with a first-order PT, however, not only depends on the nucleation temperature, but also on various other quantities that characterize the transition's properties. First, the phase transition's inverse duration β compared to the expansion rate of the Universe at the time of the transition is a crucial parameter that can be calculated as 6 In the absence of a calculable model Eq. (37) cannot be evaluated without further assumptions on S 3 /T . One way is then to follow Refs. [4,5] and approximate dS 3 /dT S 3 /T . In that case β/H can be estimated as where the last relation holds for a wide range of transition temperatures T n . Since first-principle calculations of β/H are lacking for theories exhibiting a chiral phase transition in a strongly coupled hidden sector, the rough estimate in Eq. (38) was employed for a long time in the corresponding literature, see e.g. [16,59,60]. In contrast, the use of effective models in the present work permits us to determine the function S 3 (T )/T , so that we can evaluate Eq. (37) explicitly. A second important quantity in the calculation of gravitational wave spectra is the phase transition's latent heat normalized to the radiation energy density, i.e.
where ∆V eff (T ) := V eff (0, T )−V eff (σ min (T ), T ) withσ min (T ) being the temperature-dependent global minimum. The parameter α as defined above is a measure of the phase transition's strength.

IV. STOCHASTIC GRAVITATIONAL WAVE SIGNAL
The production of gravitational waves (GW) during violent astrophysical processes like black hole mergers is a well-known phenomenon, which is predicted by Einstein's theory of General Relativity [99]. The first GW signals of this kind have been measured recently by the LIGO and VIRGO collaborations [1][2][3]. Apart from these transient signals due to astrophysical sources, a stochastic background of GWs may exist as a result of, for instance, first-order cosmic phase transitions (PT) or an inflationary phase in the early universe [100]. Observing such a background would thus provide a unique opportunity to investigate primordial physics which would otherwise be impossible to study. As discussed in Section III, a first-order PT proceeds via the nucleation and subsequent growth of bubbles of the true groundstate. Gravitational waves are produced when the aforementioned bubbles collide via several different processes, namely collisions of bubble walls or, equivalently, scalar field shells [101], as well as sound waves [102] and magnetohydrodynamic turbulence [103] in the thermal plasma. The relative importance of these contributions is determined by the dynamical properties of the phase transition and hence depends on the underlying particle physics model. Still, it is possible to broadly distinguish between two different transition scenarios which are usually referred to as runaway and non-runaway, see e.g. [61].
In the former case the friction exerted by the thermal plasma on the nucleated bubbles is too small to slow down their accelerated expansion, which, in turn, is driven by the released latent heat. Consequently, the bubble wall velocity v b will continue to increase until it eventually reaches the speed of light c. In the non-runaway scenario on the other hand, the friction between bubbles and the surrounding plasma is large enough to counteract the bubbles' accelerated expansion. The bubble wall velocity hence approaches a terminal value, which is bounded by c but may still be relativistic.
In order to determine which of the above cases applies to the QCD-like hidden sector models discussed in the present paper, the following considerations are helpful: As outlined in Section II, the chiral phase transition in those models is driven by theσ field. The latter is expected to have sizable interactions with the theory's remaining degrees of freedom. The friction exerted by the plasma on the bubbles defined byσ is therefore anticipated to be non-negligible, which already strongly suggests a non-runaway scenario (i.e. Case I of Ref. [61]). Furthermore, we will see in the next section that our calculations predict rather small values for the phase transition strength α defined in Eq. (39), which also indicates a transition of this type. Accordingly, we will assume in the following that the χPT proceeds via non-runaway bubbles in all of the considered models and for all of our benchmark points. 7 In non-runaway scenarios, the fraction of the PT's latent heat that goes into kinetic energy of the scalar field is negligible, so that the contribution to GW production from scalar field shell collisions becomes irrelevant [61]. In contrast, and as a result of the large friction between bubbles and the surrounding thermal plasma, a substantial portion of the available energy is converted into bulk motion in the form of sound waves and magnetohydrodynamic turbulence. Typically, however, only a rather small part of bulk motion is turbulent and it can be shown that the corresponding contribution to GW production is only relevant for slow PTs [61], i.e. β/H 1. As we will explicitly see in Section V, the χPTs under consideration are throughout predicted to proceed very fast, β/H O(10 4 ). Consequently the associated GW background arises virtually exclusively from sound waves in the non-runaway case. The corresponding energy density is [61] Ω sw h 2 = 2.65 · 10 −6 H β where the spectral shape S sw (f ) and the spectrum's peak frequency f sw are given by respectively. As indicated above, α and β/H characterize strength and inverse duration of the PT (cf. Section III), while g is the effective number of relativistic degrees of freedom in the plasma and κ v denotes the fraction of latent heat converted into bulk motion. The actual terminal wall velocity v b of non-runaway bubbles is determined by the details of the underlying particle physics model. Calculating its specific value is, however, beyond the scope of the present work. In the following we will instead assume a certain range of different v b and investigate how the resulting GW spectra change. In doing so we will concentrate on highly relativistic bubble wall velocities, v b ≥ 0.75, since those lead to the strongest GW signals and are thus most interesting from an observational point of view. In the considered case of non-runaway bubbles with large wall velocity v b c the efficiency factor κ v is then given by [104] κ v ≈ α 0.73 + 0.083 √ α + α .

V. RESULTS
Having summarized the basics of first-order (chiral) phase transitions and the associated gravitational wave signal in the previous two sections, we will now apply the described formalism to derive predictions for the gravitational wave spectrum within the low-energy effective models introduced in Section II. Let us once more stress that using the aforementioned models puts us in the position to explicitly calculate the parameters characterizing the chiral phase transition (e.g. β/H) without having to resort to heuristic arguments or rough estimates. We will do so in the following, starting with the benchmark points introduced in Table I, for which the transition is expected to occur around the QCD scale. Building on these results we will, in a second step, also determine the gravitational wave signal originating from hidden chiral phase transitions at (much) higher temperatures.

A. Hidden chiral phase transitions at O(100 MeV)
Employing the finite-temperature effective potentials of Eqs. (12), (20) and (26), we can directly compute the bounce solutionσ b (r) for a given temperature T . Eqs. (30) and (32) Table II. The result of such a calculation is shown in Fig. 2 for one of our benchmark points. Interestingly, we observe that S 3 /T quite precisely follows a function of the simple form which was already anticipated by Hogan in Ref. [5]. Fitting our data points to the function in Eq. (43), the parameters γ and b can be determined, see the last two columns of Table II. Intriguingly, the exponent γ is very similar among all benchmark points and effective models. In contrast, the coefficient b is found to vary by several orders of magnitude. Using the best-fit values for γ and b, we plot the function of Eq. (43) alongside the data in Fig. 2 (solid lines).
For comparison, we also show natural cubic splines interpolating between the data points (dashed lines). All further calculations involving S 3 /T will, however, make use of the above described fit, since this is less sensitive to numerical errors. Given S 3 /T as a function of T we can now straightforwardly compute the nucleation temperature T n by solving Eq. (36). Most importantly, we always find T n T c implying that there exists virtually no supercooling, which a posteriori justifies our choice of the Hubble parameter in Eq. (34), as well as our assumption that the nucleation temperature approximately coincides with the transition temperature (cf. footnote 6 on page 12).
Next, the phase transition's inverse duration β normalized to the expansion rate of the Universe at the time of the transition can be determined from Eq. (37). The outcome is again compiled in Table II. Although the precise values for β/H vary between different models and benchmark points, we generally find them to be of order 10 4 or even larger. Note that these results are in stark contrast to the usual assumption 8 β/H ≈ O(100), which is being made in discussions of (hidden) chiral phase transitions throughout the literature, see e.g. [16,60]. 9 As we will see in more detail below, this has far-reaching consequences on the observational prospects for the associated gravitational wave signal. Lastly, the transition strength α introduced in Section III is obtained using Eq. (39), see Table II for our results. Let us note that, whereas T n and β/H depend only very mildly on the effective number of relativistic degrees of freedom, α is by definition rather sensitive to the precise value of g . By Eqs. (40) and (42) the same is then true for the energy density Ωh 2 of the GW background (but not for its peak frequency, cf. Eq. (41)). To be more precise, it is straightforward to see that increasing g effectively decreases Ωh 2 . For definiteness and unless explicitly stated otherwise, we have fixed g = 47.5 in all calculations, which corresponds to n f = 3 light hidden fermions and n 2 c − 1 = 8 hidden gluons. For models with different hidden sectors the presented results must be appropriately adapted.
With the quantities T n , β/H and α at hand, and assuming that the chiral phase transition proceeds via non-runaway bubbles with some terminal wall velocity v b , we are now able to compute the predicted gravitational wave spectra Ωh 2 as described in Section IV. Our findings for the benchmark points from Table I are displayed in Fig. 3. Recall that both position and height of the spectrum's peak are solely determined by the sound wave contribution to Ωh 2 . Importantly, Eq. (41) reveals that values of β/H ≈ O(10 4 ) imply a peak frequency at approximately 1 mHz for bubble nucleation at roughly 100 MeV. A first-order chiral phase transition at such temperatures is consequently predicted to produce stochastic GW background within the LISA frequency band. At the same time, large β/H suppresses the peak height according to Eq. (40). Unfortunately, Fig. 3 shows that our effective model analysis fails to conclusively answer whether the GW signal can be expected to be observed by LISA. More generally, comparing the predictions of different effective models for a given benchmark point (i.e. a fixed mass spectrum), we find that the parameters characterizing the chiral phase transition as well as the resulting gravitational wave spectrum may vary considerably. We therefore conclude that first-principle calculations, like lattice simulations, are required in order to get quantitatively robust predictions.  Table I together with the power-law integrated sensitivity curve [105] for the LISA experiment assuming five years of running and a threshold signal-to-noise ratio of 5. The strain noise power spectral density for LISA was adopted from Ref. [106]. The displayed signal bands were obtained for a fixed g = 47.5 by varying the bubble wall velocity between v b = 0.75 and v b = 1.

B. Hidden chiral phase transitions at higher temperatures
In the previous section, we exclusively studied benchmark points with an inherent scale of order 100 MeV. However, as we have argued in the introduction, there are equally well motivated BSM scenarios, where a hidden chiral phase transition is anticipated to occur at higher temperatures. In order to investigate those cases as well, we will now simply consider scaled-up versions of the benchmark points in Table I, obtained by multiplying all meson masses etc. by a common dimensionless factor ξ > 1, i.e.
Since it is fully determined by model parameters, the critical temperature T c will then scale with the same factor. For the same reason the action functional as a function of T /T c will remain unaltered, i.e.
In contrast, the nucleation temperature T n only approximately scales with ξ. An exact scaling is violated by the fact that Eq. (36) contains the Planck mass as an absolute energy scale, which is kept constant. Accordingly, the ratio T n /T c is observed to decrease mildly for growing ξ, while it stays always close to one. A small change in T n /T c can, however, have a sizable effect on β/H  Table I. The shown bands illustrate the dependence of the result on the effective number of relativistic degrees of freedom g , which we considered to range from g = 1 (upper edge) to g = 100 (lower edge).
for the rescaled point since S 3 /T has a pole at T = T c , close to which its derivative varies rapidly, see Fig. 2. The resulting scale dependence of β/H is displayed in Fig. 4. Note that β/H stays of order 10 4 for a wide range of critical temperatures in all considered models. Lastly, in line with its definition in Eq. (39) the parameter α was observed to not depend strongly on the precise value of T n /T c , so that it can be safely assumed as constant in ξ for any fixed choice of g . Following the above discussion and based on the results compiled in Table II, we can now compute the GW signal for any given ξ. Importantly, the spectrum's peak frequency is proportional to the nucleation temperature and will thus approximately scale linearly with ξ, cf. Eq. (41). For hidden chiral phase transitions occurring between roughly 1 GeV and 10 TeV (corresponding to ξ ranging from 10 to 10 5 ) the associated GW backgrounds are therefore anticipated to lie within the DECIGO/BBO frequency band. As Figs. 5 and 6 exemplarily demonstrate, our effective model study strongly suggests that the aforementioned experiments are, in fact, sufficiently sensitive for the GW signal to be observed.
The above described analysis can, of course, be repeated for larger scaling factors ξ in order to obtain results appropriate to hidden chiral phase transitions at even higher temperatures. Crucially, we find that for T n 10 TeV (or, equivalently, ξ 10 5 ) all considered models unanimously predict the GW signal to lie outside the sensitivity range of all planned gravitational wave observatories, cf. Fig. 7.

VI. SUMMARY AND CONCLUSIONS
In the present paper, we investigated the prospects for observing a stochastic gravitational wave background associated with a first-order chiral phase transition (PT) in a strongly coupled hidden or dark sector. For definiteness we studied a new confining SU(n c ) gauge force with n c = 3 colors acting on n f = 3 flavors of massless fermions. In order to reliably model the strong dynamics that drive the aforementioned transition, we employed various effective theories known to capture non-perturbative aspects of the low-energy regime of real-world QCD: the Nambu-Jona-Lasinio model, its Polyakov-loop enhanced variant and the linear sigma model. This approach allowed us to explicitly calculate the temperature-dependent action functional S 3 /T , which is, in turn, crucial to determine key characteristics of the hidden PT like its duration or the nucleation rate of critical bubbles of the true vacuum. We found that, in the vicinity of its pole at the critical temperature T c , the function S 3 (T )/T is of the simple form b(1 − T /T c ) −γ for some b, γ > 0. Intriguingly, the exponent γ is very similar among all considered benchmark points and effective models, namely γ 1.8 − 1.9. 10 As a next step, we then computed the phase transition's nucleation temperature T n , as well as its inverse duration normalized to the Hubble parameter β/H. Our results consistently predict the transition to proceed very fast with virtually no supercooling, T n T c . Crucially, our findings of β/H O(10 4 ) are in stark contrast to β/H ≈ 100 which is typically assumed in the literature whenever explicit calculations are unavailable. Lastly, we demonstrated that the observed enhancement of β/H by two orders of magnitude has significant impact on the predicted gravitational wave spectrum. The latter was determined for various transition temperatures appropriate for different well-motivated beyond-the-Standard Model scenarios featuring a new strongly coupled sector.
Our study suggests that a gravitational wave background associated with a hidden chiral PT at temperatures of order 100 MeV (1 GeV to 10 TeV) may be detected by LISA (DECIGO/BBO). In contrast, if the transition occurs at temperatures above roughly 10 TeV the corresponding GW signal is unlikely to be observed by any currently proposed experiment. Let us further remark that although all discussed low-energy effective models tend to give qualitatively similar results, the predicted gravitational wave spectra may vary considerably. This highlights the necessity of true first-principle calculations, like lattice simulations, for the purpose of reaching also quantitative precision.
While the present work focused on the the case of a new SU(3) gauge sector with three massless fermion flavors, our analysis can, in principle, likewise be applied to different combinations of (n c , n f ). Those are then, however, not automatically guaranteed to give a first-order PT (see Ref. [16] for an overview). We also do not expect our results to change qualitatively if the hidden fermions acquire a (sufficiently small) current mass. Thus, our and analogous studies can be employed to investigate the prospects for observing potential gravitational wave signals from a large variety of well-motivated new physics scenarios featuring strongly coupled hidden or dark sectors. where M c = M (v σ ), with M (σ) given by Eq. (13d), and v σ denotes the zero-temperature vacuum expectation value of the σ field. All of the above integrals correspond to diagrams which arise from integrating out the fermion fields. Similarly, the expression for the pion decay constant in the (P)NJL model contains Note that the integrals in Eqs. (A1) and (A2) are divergent and must therefore be regularized. We employ a hard four-dimensional Euclidean cutoff Λ to do so. Next, the integrals relevant for determining the wave-function renormalization Z σ in the (P)NJL model (cf. Eq. (B8)) are given by Lastly, the bosonic (B) and fermionic (F) thermal integrals needed to compute the one-loop thermal contributions to the effective potentials in Section II are defined as Appendix B: Details on the (P)NJL model The present appendix elaborates on certain aspects of the (P)NJL model and is meant to complement the corresponding discussions in Section II. Let us start with a slightly more detailed account on the self-consistent mean field approximation (MFA) [75,[77][78][79]. For convenience, we repeat here the Lagrangian for the three-flavor NJL model which is given by Using the MFA, the above Lagrangian can be written in terms of the fermion fields q i and auxiliary meson fields, cf. also Ref. [55]. This is done by going from the perturbative to the Bogoliubov-Valatin (BV) vacuum and using Wick's theorem for operator products normal ordered with respect to the BV vacuum. Consequently, one can split up the Lagrangian in an interacting part L int and a mean field term L MFA , such that L MFA contains all terms which are at most quadratic in the fermion fields. One can then rewrite L MFA in terms of Ψ and its average Ψ , which is defined as a sum over mesonic auxiliary fields or equivalently Working along the steps outlined above one can thus write the mean field Lagrangian in terms of the effective meson fields σ, a a , π a , η and the fermion field q. The NJL Lagrangian in the MFA is given by where repeated flavor indices a, running from 1 to n 2 f − 1 = 8, are summed over and we have defined π := 2π a T a and a := 2a a T a , which are matrices in flavor space. The effective fermion mass M was already introduced in Eq. (13d) and is given by The tree-level potential is V NJL tree = 1 8G 3σ 2 + 3η 2 + 2π a π a + 2a a a a − G D 16G 3 σ σ 2 + π a π a − 3η 2 − a a a a + 5a a π a η . (B6) Next, we briefly discuss the zero-temperature one-loop meson propagators, whose roots define the meson's effective pole masses. The propagators can be determined based on the mean field Lagrangian of Eq. (B4) by calculating all one-loop 1PI diagrams with two external φ lines, where φ stands for one of the meson fields. The correlator iΓ φφ is then given by the sum of all those diagrams. To find the relevant diagrams one has to keep in mind that only the fermions can run in the loop, because the mesons are represented by auxiliary fields in the NJL model and are thus non-propagating at tree-level. Furthermore, note that at zero temperature σ has a non-zero vev, v σ , which implies that also diagrams with v σ as an external source contribute to the meson propagators. Taking all these considerations into account, the one-loop meson propagators are computed to be Explicit formulas for I V , I S and I P can be found in Appendix A. Equivalent results are, for example, given in Ref. [113], where the meson propagators are determined for non-zero quark masses. Our results agree with theirs when taking the chiral limit of the latter. Finally, we provide the full expression for the σ field's wave-function renormalization Z σ as obtained within the MFA from its definition in Eq. (14): with r ≡ r(σ) = |M (σ)|/T and where M ≡ M (σ) is again the effective quark mass from Eq. (B5) and Λ denotes the fourdimensional hard momentum cutoff used to regulate the divergent vacuum parts of the occurring loop integrals. Expressions for the thermal parts i can be found in Eq. (A3) of Appendix A.

Appendix C: Basics of the CJT formalism
The results of the linear sigma model presented in this work's main part were derived within the composite operator formalism due to Cornwall, Jackiw and Tomboulis (CJT) [88]. Their formalism extends the approach to QFT based on the conventional effective action Γ[φ] by employing a generalized functional Γ[φ,Ḡ], the so-called two-particle irreducible (2PI) action. Here,φ and G are expectation values 11 of some quantum field Φ and of the corresponding two-point function, respectively. While the original CJT paper [88] only considered the vacuum case T = 0, the composite operator method was subsequently extended to finite-temperature field theory problems in Refs. [96,114], see also [115]. It was then employed to calculate the chiral phase transition in different versions of the linear sigma model e.g. in Refs. [91][92][93][94]. 12 In the present appendix, we will only give a brief overview of the CJT formalism's very basics. For more details we refer the interested reader to the given references.
In situations with translational invariance, however,φ becomes constant andḠ can be taken as a function of (x − y) only. All relevant information is then encoded in the 2PI effective potential V CJT [φ,Ḡ]. In the chiral limit of the U(3) × U(3) linear sigma model (LSM) discussed in Section II C, where mixing between the different meson species is absent, the aforementioned potential can be written as [93,94] V LSM CJT [σ,Ḡ(k)] = V 0 (σ) + 1 Here, the sum runs over all meson species and the corresponding multiplicities g i are given by g σ = g η = 1 and g π = g a = 8 for n f = 3. Finite-temperature momentum space integration is abbreviated as with ω n = 2πnT being the bosonic Matsubara frequencies. Furthermore, V 0 is the model's tree-level potential given in Eq. (27a), while the ∆ i represent the tree-level meson propagators, i.e.
with the field-dependent tree-level meson masses m i of Eq. (28). Finally, the quantity V 2 is computed as the sum of all two-particle irreducible vacuum graphs with each line corresponding to a full propagatorḠ, as well as generallyσ-dependent vertices. Importantly, the conventional effective potential V eff as a function of the classical fieldσ is then obtained from Eq. (C2) as being the scalar self-energy. Note that the well-known form of the one-loop one-particle irreducible effective potential V (1) eff will be exactly reproduced from Eq. (C5), if V 2 ≡ 0 is assumed. Consequently, any non-trivial improvement over V (1) eff requires a finite V 2 . On the other hand, an exact determination of V 2 is in general not feasible since it involves the calculation of infinitely many diagrams. In practice, one must therefore resort to truncations of V 2 .
Following Refs. [91][92][93][94], we will use the so-called Hartree-Fock approximation [89,90] throughout the present work. At two-loop level, this corresponds to only retaining contributions from doublebubble graphs (cf. Fig. 8a) in computing V 2 [88]. In contrast, sunset diagrams as the one displayed in Fig. 8b are consistently ignored. According to its definition in Eq. (C7) the scalar self-energy thus exclusively obtains contributions from graphs of the tadpole topology shown in Fig. 8c.  ) and (b)), as well as the associated one-loop contributions to the scalar self-energy (graphs (c) and (d)). All dashed lines correspond to full, field-dependent mesonic propagatorsḠ i , with i being one of σ, a, η or π. Vertices are obtained from the shifted LSM tree-level potential V LSM tree (Φ + Φ), where V LSM tree and the collective meson field Φ are given in Eqs. (23) and (24), respectively, whileΦ :=σ1/ √ 6. In the Hartree-Fock approximation used in the present paper, only diagrams of the double-bubble topology (a) are taken into account in calculating V 2 , so that only tadpole graphs (c) contribute to the self-energy.
For a given field value and at a given temperature, the system in Eq. (C9) can be solved for the four effective masses M i , which, in turn, define the full, field-dependent meson propagators according to the ansatz in Eq. (C8). The sought-after finite-temperature effective potential can then be computed via Eq. (C5), which eventually leads to the final expression quoted in Eq. (26).