Model Discrimination in Gravitational Wave spectra from Dark Phase Transitions

In anticipation of upcoming gravitational wave experiments, we provide a comprehensive overview of the spectra predicted by phase transitions triggered by states from a large variety of dark sector models. Such spectra are functions of the quantum numbers and (self-) couplings of the scalar that triggers the dark phase transition. We classify dark sectors that give rise to a first order phase transition and perform a numerical scan over the thermal parameter space. We then characterize scenarios in which a measurement of a new source of gravitational waves could allow us to discriminate between models with differing particle content.


Introduction
The detection of gravitational waves (GW) [1] established a new and independent probe of New Physics. It has already been suggested that the data from resolvable events such as from binary mergers could help constrain interacting dark matter [2,3] or exotic compact objects [4][5][6]. The observation also implies that we may anticipate the detection of a stochastic GW background at both current and future detectors. This may be a rare probe of the Cosmic Dark Ages and the first observational window onto cosmic phase transitions (PTs). Such cosmic phase transitions leave behind a characteristic broken power-law gravitational wave spectrum. The relic spectral shape depends on the strength of the transition, the speed of the transition, the bubble wall velocity and the temperature of the transition.
In terms of the LISA-inverse problem, much attention has been focused towards either arguing for a new scale of physics [10,25], or for relic backgrounds from certain well motivated extensions of the standard model, assuming the reheating temperature is sufficiently high [8, 9, 11-24, 31, 32]. Little work to date has focused on the question of model discrimination [18]. In this work we endeavour to see how much model discrimination is in principle possible from the frequency spectrum of a future stochastic gravitational wave signal. In particular, we consider renormalizable and non-renomalizable effective field theories of interacting hidden sectors, in which a gauge symmetry is spontaneously broken. We also consider the effect of fermions that couple to the scalar.
Simulations of gravitational wave backgrounds from cosmic phase transitions indicate that there are three spectral contributions: the collision spectrum is the direct effect of bubbles of true vacuum colliding, the sound wave spectrum is the result of the fluid dynamics after such collisions, and the turbulence spectrum, which is usually subdominant. It has been realized recently that the sound wave contribution dominates in most relevant scenarios. In particular, this is true in all cases that do not display "runaway" behaviour, and such runaway is blocked by any gauge bosons acquiring a mass in the transition [33].
All spectra are controlled by four thermal parameters: the velocity of the bubble wall, v w , the ratio of the free energy density difference between the true and false vacuum and the total energy density, ξ, the speed of the phase transition β/H and the nucleation temperature T N . In the special case in which two peaks are visible, the four thermal parameters can in principle be reconstructed.
In this paper we focus on the thermal parameters T N , β/H and ξ. The first determines the scale of the phase transition, and the latter two are most powerful at model discrimination.
To study the thermal parameters in a general context, we observe that first order phase transitions are realized in (effective) double well potentials, from the interaction of terms with alternating signs. As such, we will study multiple models within two limiting scenarios, where all coefficients are positive at the time of transition. Most phase transitions can be mapped onto these effective scenarios. In particular, the EWPT in a Higgs+singlet model is an example of (1.2), upon integrating out the heavy singlet (up to dimension-6 operators). The thermal parameters are strongly dependent on the nature of the thermal corrections. These thermal corrections are functions of the bosonic and fermionic degrees of freedom coupling to the scalar. The bosonic degrees of freedom are given by the gauge structure of the theory. Here we will consider models within the following scenarios, corresponding to the limiting cases (1.1) and (1.2): 1. A dark Higgs -SU (N ) breaking into SU (N − 1). In this case the barrier between the true and false vacuum during the transition is caused by dark gauge bosons that provide an effective cubic term.

2.
A dark Higgs -SU (N ) breaking into SU (N − 1) -with significant non-renormalizable operators. In this case the barrier between the true and false vacuum is caused by the quartic dark Higgs coupling being negative and the vacuum being stabilized by the positive Wilson coefficient of the sextet interaction.
Such scenarios may arise for example in the context of Composite Higgs models of cosmology [34,35], where the Dark Higgs would be represented by a pseudo-Goldstone boson state (a generalization of a QCD pion) whose interesting potential is due to explicit breaking of the global symmetry via SU (N ) gauge and Yukawa couplings.
In each case we consider gauge groups of different ranks as well as models with and without a thermal mass produced by dark fermions. For all cases ξ is independent of the scale of the potential and β/H has a weak logarithmic dependence whereas both thermal parameters are controlled by the ratio of the vev with the scale of the potential x ≡ v/Λ. Therefore the renormalizable potentials are a 2(3) parameter problems for each model and the non-renormalizable potentials are a 3(4) parameter problem without (with) the addition of a dark fermion.
We find that non-renormalizable operators dramatically improve the visibility of gravitational wave spectra, whereas adding a dark fermions N f and increasing the rank of the group N provide a more modest boost, which becomes reasonably large in the limit of large N f or N . The boosts to visibility in each case are non-degenerate. In the renormalizable case (1.1), we find that both the effect of a larger gauge group (SU (N ) → SU (N + 1)), and the effect of increasing the number of fermions (with significant thermal mass) are essentially to shift the thermal parameter space, and increase the detection prospects. Of course, there is a degeneracy of predictions for specific models. It has been suggested that anisotropy measurements could break this degeneracy, for example by a cross-correlation with the CMB data [36].
The structure of this paper is as follows. In section 2 we summarize the models we are attempting to discriminate. In section 4 we review the spectra of gravitational waves from a cosmic phase transition and in section 5 we present our results. In section 6 we relate our results to studies of dark matter, before concluding with a discussion and an outlook to future work in the final section.

Scenarios for a dark first order Phase Transition
A first order phase transition may occur for a potential with three competing terms, with alternating signs, such that it has a double well separated by a barrier. Moreover, the vacuum energy corresponding to these minima will be temperature-dependent, such that the ground state changes as the Universe cools. The first order phase transition may then happen if the potential barrier is present at the critical temperature T c , when the minima are degenerate.
We will consider two limiting cases of such potentials. In the renormalizable case, the potential barrier is generated effectively at finite temperatures, but does not exist at zero temperature. As we will see, the zero-temperature masses and self-couplings, the quantum numbers of the scalar, and the couplings to fermions crucially determine the thermal parameters of the phase transition. For all the models we are considering the part of the Lagrangian relevant to phase transitions can be written We will consider potentials of the form (1.1) and (1.2).

SU (N )/SU (N − 1) models with renormalizable operators
The first case has a double well generated from the quadratic, cubic, and quartic interactions at finite temperature. We parametrize the potential such that the overal scale (Λ) and the zero temperature vacuum expectation value (v) are inputs. This implies the following redefinitions for zero temperature parameters in the potential (1.1), As we will see below, we find that some thermal parameters are only functions of the ratio of the zero temperature vev and the scale of the potential (v/Λ). Using this parametrization, the finite temperature potential is given by, 5) where N G is the number of gauge bosons coupling to the scalar sector with coupling constant g, N GB is the number of Goldstone degrees of freedom, and N f is the number of fermions with Yukawa coupling y. For simplicity, we consider degenerate Yukawa couplings, as the gravitational waves produced by (y, N F ) and ({y i }, N F ) are related by y 2 × N F = N F y 2 i . 1 In the second line we have applied a high temperature expansion, All field dependent masses which enter into the effective potential are provided in the appendix. 2

SU (N )/SU (N − 1) models with non-renormalizable operators
The second limiting case has the double well resulting from the interplay between the quadratic, quartic, and sextic terms. We again choose a parametrization of the potential such that the scale of the potential Λ and the zero temperature vacuum expectation v value are inputs. This will leave us with one free parameter α, which parameterizes the difference in vacuum energy of the two minima at zero temperature. In the high temperature expansion (2.6), the potential becomes, It is seen that at zero temperature, the potential has minima at h D = 0 and h D = v respectively, overall scale Λ, and the (dimensionless) non-renormalizable coupling is α. That is, we have made the following redefinitions in Eq. (1.2): Up to a small change in the number of relativistic degrees of freedom g * . Since the gravitational wave spectra has a very weak dependence on g * , making this simplification is at little cost to generality. 2 Note that the use of perturbation theory introduces some theoretical uncertainty as perturbativity at finite temperature breaks down above the critical temperature [37,38], a fact that can be delayed somewhat by the inclusion of "daisy terms" [39] although in reality one requires a lattice simulation for a robust treatment. In spite of this theoretical uncertainty we expect our results to be indicative of the overall thermal parameter space including its overall scope and dependence on the model. Finally note that the most important points in our scan are where a lot of supercooling occurs and TC is significantly higher than TN meaning that these are the points where perturbation theory is most valid. At zero temperature, a value of α = 1/2 corresponds to degenerate minima, and the upper limit α = 2/3 corresponds to the value for which there is no zero temperature barrier between the vacua (as the zero temperature mass term changes sign), see Fig. 1. Of course, finite temperature corrections may allow for a higher value of α, as positive corrections to the mass term may reintroduce the barrier. Note we have once again assumed degenerate Yukawas with little loss of generality as explained in the previous section.
For operators up to dimension-6, models for the electroweak phase transition (EWPT) can be captured effectively by a special case of the above, with where we have defined, Here Λ 6 is the scale associated with the dimension-6 operators which arise from integrating out BSM physics, such as a singlet scalar. 3 . We will also consider the EWPT with nonrenormalizable operators for the sake of comparison later. Finally, note that if one rewrites Eq. 2.7 in terms of implicitly defined temperature dependent parameters one can follow the process in [31,40] and fit the action to the function for the range α(T ) ∈ [0.51, 0.65] In this section we give examples of hidden sector models which can be mapped onto our general framework given above. Of course we are not completely general as we do not consider for example the case where multiple scalars acquire a vev at the same time (such as a multi dark higgs doublet model) or more complicated gauge group structures SU(N)×SU(N ) where both gauge couplings are large. However, scenarios which can be mapped onto our framework are ubiquitous including, Pati-Salem symmetry breaking 4 [41], colour breaking intermediate phase transitions [42,43], atomic dark matter [44], asymmetric dark matter [45] and compositeness [46] to give a non-exhaustive list. We give more details of three of these examples and how they map to the various models we consider below.

Generalized baryon number
As was suggested in [45], the dark sector relic abundance and the baryon asymmetry in the SM can have a common origin in models with a generative symmetry breaking. In such models, there is a generative gauge group G, for example SU (2) G which is broken spontaneously through a first-order phase transition in the early universe. The asymmetry generated in this phase transition is communicated to the dark and visible sectors through a mixed Yukawa term. The degenerative scalar has tree-level zero temperature potential, and quartic mixing terms with the SM Higgs, B-L breaking scalar σ, and dark scalar χ. For small mixing, such as is the case in various supersymmetric models, the mass contributions are small. For non-supersymmetric models, the mixing can be significant, and contribute to the thermalization and decay properties of the various sectors. The mass hierarchies are small, such that the scalar ϕ can have a mass at the electroweak scale. In this case there are significant cosmological and astrophysical constraints as discussed in [45]. The first order PT can be induced when one includes an effective dimension-6 operator, which can arise at the one loop level from the mixed quartic interactions [21] from which it is seen that this an example of a model within the scenario given by (2.7). 4 This phase transition is more likely to occur at a scale visible to aLIGO than LISA

Atomic dark matter
A further possibility is that the dark sector contains a confining group, as well as fermions charged under an unbroken U (1) . Then, dark atoms can be formed [44]. The strongest constraint on atomic dark matter comes from the self scattering bound [47,48], where m χ is the heavier particle, which forms the nucleus of dark atoms. The mass of m χ can be heavier than a TeV [49] in which case the constraint on the gauge coupling is very modest (α D ∼ 0.1, implying g ∼ O(1)). A simple example is an SU(4) gauge group, which breaks into SU (3) ×U (1), allowing for the formation of nuclei during dark BBN [50].

Composite Dark Matter models
A final example is a dark matter candidate as the lightest bound state of a confining gauge group SU (N ), such as has been discussed in [32]. The spontaneous symmetry breaking of an approximate global symmetry, which is only partially gauged, gives rise to pseudo-Goldstone bosons. These light states are sensitive to an effective scalar potential at the 1-loop level, which in turn initiates a further breaking. A particularly interesting possilibility has the SM Higgs and the dark matter candidate both as pseudo-Goldstone bosons of the same symmetry breaking [46]. Various symmetry breaking cosets have been studied in the literature, with scalar potentials of the form (2.5) or (2.7). The couplings in such scenarios correspond to 1-loop integrals in the UV theory. The GW spectra for benchmarks of thermal parameters for the breaking SU (3) and SU (4) dark gauge symmetries were previously considered in [32], where it was argued that scalar DM bound states and dark quarks (carrying EW quantum numbers) are most relevant for detection at LISA.

Thermal parameters
The dynamics of the phase transition are controlled by a bounce solution φ c (r, T ), which is a spherically symmetric classical solution to the Euclidean equations of motion [40,51,52] We compute the bounce solutions with potentials in the previous section. The thermal parameters of the phase transition can then be computed from the bounce solution. First, the nucleation temperature of bubbles of the new vacuum T N is conventionally defined as the temperature for which a volume fraction e −1 is in the true vacuum state. This corresponds approximately to where p(t) is the nucleation probability per unit time per unit volume, and where t N is the nucleation time. The nucleation probability can be calculated from the bounce solution as, where S E is the Euclidean action evaluated on the bounce. We assume a radiation dominated universe to relate the nucleation temperature and time. The speed of the phase transition is controlled by the parameter β, which can also be related to the bounce action, Last, the latent heat parameter is given by, Where ∆ indicates that the quantity should be evaluated on both sides of the bubble wall, and where ρ N = π 2 g * T 4 N /30 is the equilibrium energy density at T N .

Gravitational wave spectrum and the LISA inverse problem
The gravitational wave profiles can be related to the thermal parameters. We will adopt a parametrization introduced by [53], but our analysis can be adapted when future models become available. In principle, there are three contributions to the power spectrum, Ω GW = Ω col + Ω sw + Ω turb (4.6) Where the first term corresponds to the spectrum from bubble collisions, the second is a spectrum due to sound waves in the fluid after collisions, and the third a turbulence term. As realized last year [33], in any model in which gauge bosons gain a mass in the transition, the bubble wall velocity approaches a finite limit. Therefore, the sound wave contribution [28] is typically dominant in all of the cases we consider in this work. Its power spectrum can be expressed as [53], where Γ ∼ 4/3 is the adiabatic index, andŪ 2 f ∼ (3/4)κ f ξ is the rms fluid velocity. For v w → 1, the efficiency parameter is well approximated by [54] κ f ∼ ξ 0.73 + 0.083 √ ξ + ξ (4.8) For v w ≈ 0.5, we use [54] κ f ∼ ξ 2/5 0.017 + (0.997 + ξ) 2/5 (4.9) Figure 2: (Schematically) the LISA inverse problem. In the above, the subscript x refers to the dominant peak of the GW spectra (collision, sound wave, or turbulence). As described in the text, for most models the sound wave contribution is dominant. The thermal parameters of the PT can be calculated by solving the bounce EOM (4.1), and then related to the GW spectra using (4.7) and (4.11). This paper finds general relations between the GW spectra and the Lagrangian. and the spectral shape is given by From this we notice that the amplitude of the signal is a function of the parameters β/H, the wall velocity v w , and the latent heat ξ; whereas the position of the peak depends on β/H and T N . We will use this insight in the next section, to compare the predictions of the different models (2.5) and (2.7). This effort can be summarized by the LISA inverse problem, in fig.  2. We should mention some previous work towards solving the LISA inverse problem. The link between gravitational waves detection of collision and turbulence peaks and the thermal parameters has previously been summarized in ref. [25], which highlighted visible regions in the thermal parameter space. On the link between the Lagrangian and thermal parameters some thorough work has been done in the case of the EWPT with extended scalar sectors, [23,[55][56][57]. The aim of this paper is to compliment these previous works by studying the general case of a (single) scalar, with couplings to different numbers of fermions and gauge bosons, as well as other scalars separated in mass.

Spectra from models
We compute the thermal parameters for scenarios (2.5) and (2.7), for different dark sectors. We are specifically interested in light scalar sectors, with masses around the EW scale. For comparison, we also study the SMEFT case, in which the electroweak phase transition is catalyzed by a non-renormalizable H 6 effective operator. The SMEFT case is then well approximated by a dark SU(2) with three dark fermions.
We find bounce solutions using two techniques to ensure accuracy: a numerical finitedifference algorithm, where we discretized the radial direction r and the analytic technique described in section 2. The thermal parameters are then found by substituting the bounce solution into the Euclidean action S E as described in the previous paragraph.
In both the renormalizable and non-renormalizable models the thermal parameter set (ξ, H/β) governs the peak amplitude. We find that these results are essentially independent of the scale of the potential Λ. Specifically, ξ is independent, whereas β/H has a weak Logrthmic dependence. The nucleation temperature by contrast scales linearly with Λ. In the case where we have only renormalizable operators (2.5), we scan over (g, v/Λ), with scan ranges g ∈ (0.1, 1), and v/Λ ∈ (0.5, 4) In the non-renormalizable case (2.7), we fix g = (0.5, 1) and scan over (α, v/Λ), where we fix Λ = 200 GeV. The scan ranges are v/Λ ∈ (0.5, 4) and α ∈ (0.55, 1.5). We assume that the fermions are massless before the PT. The parameter that enters the scan is then N f × y χ . For convenience, we have set y χ = 1 in the figures.
We summarize the results for the peak amplitude and peak frequency in Figs. 3 and 4 respectively where in the spirit of reference [25] we include visibility curves for LISA and plot the (ξ, β/H) and (β/H, T n ) planes. We check explicitly that the high temperature expansion is valid for the results of our scan, by ensuring that 2m 2 i < T 2 with i = h, GB for the gauge boson and Higgs mass at the critical vev and temperature. 5 The effect of excluding points for which this check fails is to mildly trim the very tip of the peaks of the thermal parameter space in Figs. 3. The fact that the trimming occurs for low dark Higgs mass can be understood in direct analogy with early studies of the EWPT (before the Higgs mass was known). In this model one finds that for fixed vev, the strength of the phase transition grows inversely with the Higgs mass. In the limit of small Higgs mass, the gauge boson masses (which scale with v(T n )) become large, invalidating the high temperature expansion (m G /T < 1) to be valid.
The different shape of the results for the potential (2.5) with fermions can be understood as the fermions contribute only to the mass term. Therefore the potential barrier is no longer just a function of the gauge coupling, which we scan over, and the zero temperature mass. The reader will also notice that the results for the different potentials (2.5) and (2.7) have different zero temperature mass ranges. This can be understood by considering the contribution of the dimension-6 term to the latter.  .7), where we have chosen vw = 0.5 in the left plot, and vw = 1 in the right plot (with the corresponding efficiencies from [54]), as motivated using the conditions in [58]. The upper thicker contour corresponds to the LISA 1-year peak sensitivity [59]. The lower thicker dashed contour corresponds to LISA for a power-law spectrum (integrated over frequency), taken from [60]. The width of the contours is found from varying the zero-temperature potential parameters. Left: unless otherwise indicated, the number of Yukawa couplings is taken to be zero. If present, the Yukawa couplings are set to yχ = 1. Right: unless otherwise indicated, g = 1. The light blue dashed line corresponds to the predictions from the EWPT.
From the results for the non-renormalizable operators, it would naively seem that gauge bosons and fermions change the zero temperature mass of the scalar. However, the more accurate statement is that the presence of fermions and the rank of the gauge group determines which zero temperature masses lead to a strong first order PT, and are not disallowed by supercooling. Furthermore, for the case where g = 0.5 rather than g = 1, the high temperature expansion is valid for lower dark Higgs masses, before it is rendered invalid by large gauge boson masses.
In the right panels, we compare our result to the predictions from the EWPT up to dimension-6 operators (2.11), with the dashed blue line. We find that the results in Fig. 3 overlap, demonstrating that these results are insensitive to the scale Λ (but sensitive to the ratio v/Λ). As expected, the predictions for the peak frequency (Fig. 4) do not overlap, as T N scales with Λ.
Some qualitative features can lead to model discrimination, which we list below: 1. The thermal parameter space available for SU (N ) is essentially the same as that of SU (2) apart from a shift in log ξ by an amount Figure 4: Thermal parameters from the PT described by Eqs (2.5) and (2.7) respectively. The dashed contours in the plots correspond to the sound wave peak fsw (4.11), where we have chosen the wall velocities as in Fig. 3. The thicker dashed contour corresponds to the LISA frequency peak [60]. Note that the EWPT results do not overlap with our scans, since the nucleation temperature T N is sensitive to the scale Λ.
where the coefficient A(y χ × N f ) depends on y χ × N f and is around 2.4 for y χ × N F ∼ 0 and decreases to about 1.8 for y × N f = 10. Note that in general increasing the rank of the gauge group improves visibility although one has diminishing returns for large N which we show in Fig. 3.

2.
Adding fermions qualitatively changes the available thermal parameter space slightly. Comparing N f × y χ > 0 and N F × y χ = 0, we notice a shift and a slight change in shape. For 1 < y × N f < 10 we find that the thermal parameter space merely shifts according to where we find that C(2) ∼ −0.35 for SU (2)  That is ξ is shifted in the direction of greater visibility whereas β/H is shifted in a direction of weaker visibility. Since the amplitude is more sensitive to ξ this overall means that adding fermions increases the visibility of the transition which we show in both Fig. 5 and Fig. 3. The increase in β/H is due to T c − T n reducing in magnitude as one adds strongly coupled fermions. For ξ there is a competition between two effects: the reduction in T c − T n which tends to reduce ξ and an increase in dV /dT which increases ξ. Its the latter that wins. 3. The presence of nonrenormalizable operators boosts H/β by orders of magnitude compared to what is possible in the renormalizable case. This is a striking signal suggesting that a large H/β indicates the presence of more than one new scale of physics. In this case the effect of adding extra fermions is to shift and slightly rotate the thermal paramater space (ξ, H/β), this time in the ∆ log ξ direction although the relationship is less clean than the case of renormalizable operators. In contrast the effect of increasing the rank of the group is to both shift and somewhat contract the parameter space. The shift in both cases is in a direction of increased visibility.

Relic abundance example
The scenarios discussed in the previous sections constitute hidden sectors, which may explain the present relic abundance of Dark Matter (DM). As an example, we discuss the contribution to DM relic abundance from the coupling to a single Dirac fermion to the scalar responsible for the PT. We will also assume the region m h D < m χ , which corresponds to the majority of the scenarios we covered in the last section. The fermionic DM may not have tree-level couplings to the SM, just Yukawa interactions with the Dark Higgs, and thermalize at a dark temperature, which in principle could be different from the SM evolution, T D = T SM . But provided that there was thermal equilibrium between the SM and hidden sector at some scale (above the weak scale), one can assume that at freeze-out of the χ particles T D ∼ T SM . This scenario can explain the observed DM relic abundance [61], which is mostly determined by the internal dynamics of the hidden sector.
In particular, the annihilationχχ → h D h D sets the relic abundance of χ particles.
To avoid over-closure, the h D scalar is expected to have a decay channel to the SM, such as via Yukawa couplings to the SM fermions, via a mixing θ with the Higgs, of magnitude g f = (m f /v) sin θ, where y 2 χ sin 2 θ 2 × 10 −13 [62]. This coupling is small enough such that the SM fermions are not expected to play a significant role in the h D phase transition.
Under these assumptions, the dominant annihilation cross section is p-wave, and an approximate expression for the relic abundance is then given by [63] Ω DM h 2 2.1 × 10 8 GeV −1 where the fermion masses are m χ = y χ v/ √ 2, x F = m χ /T F 20 and b = (3/128π)y 4 χ /m 2 χ . To illustrate the possible interplay between DM observations and the discovery of a new source of gravitational waves, we explore the region of correct relic abundance in the model (1.1) with Λ = 200 GeV. The results are shown in Fig. 6.
These results are based on a toy model for DM, and many other scenarios could be considered. In particular, one could explore non-thermal production of DM and its relation with the scalar potential responsible for the PT. An alternative scenario has the heavy gauge bosons of the broken symmetry as the most important component of the dark matter relic abundance. Such a scenario was considered in [64] for the symmetry breaking pattern SU (3)/SU (2), and is sensitive to additional cosmological constraints from structure formation.

Discussion
In this work we have considered the relic gravitational wave spectra from phase transitions in a hidden sector. These spectra can be related to the thermal parameters of the transition, which can be computed from first principles: β/H, the speed of the transition, Υ, the latent heat, and T N , nucleation temperature. We have distinguished between two limiting cases, with potentials (1.1) and (1.2), which effectively capture the main classes of models. Furthermore, we have studied the effect of varying the quantum numbers of the scalar, the gauge coupling, and the number of coupled fermions. The results of these studies are summarized in Figs. 3 and 4, and some general conclusions are derived in section 5. We find that although there is some degeneracy in the predictions, a level of model discrimination is possible. This is due to the fact that increasing the number of strongly coupled fermions, the rank of the group, or the number of scales involved all increase the visibility of a gravitational wave signal. Moreover the changes in thermal parameters due to each of these model changes are qualitatively different. In section 6, we comment on the relic abundance of hidden sectors that could be constrained through their GW spectra.
A few caveats to our work. First, the renormalizable potential (2.5) does not have a zero-temperature potential barrier, such as could be the case for a singlet scalar with a cubic self-interaction. Phase transitions resulting from such a potential are qualitatively different, and the thermal corrections may restore the vacuum to a unique field value in such a way that no first order phase transition is expected. If a first order phase transition does occur, it may exhibit runaway behaviour, such that the GW spectrum from bubble shell collisions becomes relevant. This would lead to a different spectral shape, which in principle may be distinguishable in future experiments, for T N around the weak scale. A detailed analysis of such a scenario is beyond the scope of the present work.
Second, in the present work we have employed a high temperature expansion, Eq. (2.6), which has a limited range of validity. Phase transitions not captured by this approximation may also give observable spectra; this is most noticeable in the results from the renormalizable potential (2.5). In future work, it will be interesting to explore the models using the full thermal functions. Another possible extension is the inclusion of higher dimensional potential corrections to the two limiting potential forms considered here. An analysis with the inclusion of such operators will be presented in a future paper.
Finally, we have not calculated the wall velocities v w in the phase transitions, instead making conservative assumptions to calculate the spectra. Calculating the bubble wall velocity for a general model with general parameters is a highly non-trivial task which we leave to future work. However, we can briefly comment on how measuring the bubble wall velocity can lead to further model discrimination. The wall velocity can be estimated in the limit that the departure from each particle in the plasma's equilibrium distributions are slowly varying near the bubble wall [65][66][67]. In this case the bubble wall velocity solves for boundary conditions h D (−∞) = 0 and h D (∞) = v(T n ), that is, the value of the non trivial minimum at the nucleation temperature. Only a particular value of the combined friction term on the right hand side will satisfy the boundary conditions, and since η is determined by particle physics, the problem reduces to choosing an appropriate value for v w γ (where γ is the Lorentz factor). In the above η can be written as a matrix product G T Γ −1 F where G and F are vectors, whose components scale as g 2 or y 2 χ . The matrix of coefficients scale either as g 2 y 2 χ , g 4 or y 4 χ . Therefore the bubble wall velocity can give more information on the size of both the gauge coupling and the fermion couplings, if present.
Future work may also include further analysis of the internal hidden sector dynamics, including a thorough calculation of the thermal histories and relic abundances of hidden sector degrees of freedom. In this work we have chosen to focus on a decoupled hidden sector, but it is in principle straightforward to extend the results presented here to sectors with significant portal couplings.
The next order in the expansion is given by a logarithm, which is cancelled by the zero temperature one-loop Coleman Weinberg potential. Note we have also ignored the constant term. We find numerically that the high temperature expansion is valid almost exactly for m 2 < 2 × T 2 . Our values of n i are n H = 1 n G = 2N − 1 n GB = 3 × (2N − 1) n f = 2 × N × N f (7.7) where N f is the number of fermions and N is the rank of the group. Note that we follow the standard practice of ignoring the second term in the high temperature expansion for Goldstones and Higgs ∼ m 3 T , such that the only cubic self-interaction comes from the gauge bosons.

Non-renormalizable potential
For the nonrenomalizable potential (2.7) we proceed as before, but here we assume the cubic corrections due to gauge bosons are subdominant compared to the zero temperature terms with alternating signs. This corresponds to only taking the first term in the high temperature expansion for every species.