Speed of sound in cosmological phase transitions and effect on gravitational waves

The energy budget for gravitational waves of a cosmological first order phase transitions depends on the speed of sound in the thermal plasma in both phases around the bubble wall. Working in the real-singlet augmented Standard Model, which admits a strong two-step electroweak phase transition, we compute higher order corrections to the pressure and sound speed. We compare our result to lower-order approximations to the sound speed and the energy budget and investigate the impact on the gravitational wave signal. We find that deviations in the speed of sound from $c_s^2 = 1/3$ are enhanced up to $\mathcal O(5\%)$ in our higher-order computation. This results in a suppression in the energy budget of up to $\mathcal O (50\%)$ compared to approximations assuming $c_s^2 = 1/3$. The effect is most significant for hybrid and detonation solutions. We generalise our discussion to the case of multiple inert scalars and the case of a reduced number of fermion families in order to mimic hypothetical dark sector phase transitions. In this sector with modified field content, the sound speed can receive significant suppression, with potential order-of-magnitude impact on the gravitational wave amplitude.


Introduction
Gravitational waves from first order phase transitions A stochastic gravitational wave (GW) background generated by a cosmological phase transition can offer a window to the early universe predating the recombination era. In a first order phase transition, a barrier separates the two degenerate minima of the free energy of the thermal plasma at the critical temperature T c . When the temperature (T ) of the plasma drops below the critical temperature, the system transitions from the high temperature to the lower temperature phase via thermal or quantum fluctuations. Bubbles of the low-T phase nucleate, expand due to the pressure difference between the phases inside and outside the bubble, and collide eventually, filling the whole universe with the low-T phase. A fraction of the latent heat of the transition gets converted into kinetic energy in the form of sound waves within the plasma, which act as a source of shear stress producing gravitational waves, resulting in a stochastic background [1][2][3][4][5][6]. A detection of such a hypothetical background is one of the science goals of the space-based GW observatory LISA [7] and other similar next-generation GW experiments.
The Standard Model (SM) of particle physics does not predict any cosmological first order phase transitions, as both the electroweak (EW) and QCD phase transitions occur via a smooth crossover [8][9][10][11]. Hence, the existence of a cosmic GW background would point to physics beyond the Standard Model (BSM). The GW signal would offer a probe complementary to high energy colliders [12][13][14]. Furthermore, a first order EW phase transition is a crucial ingredient in EW baryogenesis, a dynamical mechanism aiming to explain the matter/antimatter imbalance of the present day universe [15][16][17][18][19]. Multiple BSM theories have been proposed that are believed to incorporate cosmic phase transitions, see e.g. the references listed in [14,20].
Modelling the signal -importance of the sound speed The gravitational wave signal of a cosmological phase transition receives contributions from the kinetic energy stored in the bubble walls [21][22][23][24], sound waves in the plasma, and turbulence [25][26][27][28][29]. In the case of weak to moderate phase transitions, the sound wave contribution is expected to be dominant [2,3,6], and it will be the main focus of this work.
The gravitational wave signal from sound waves has been determined in a range of numerical simulations [2,4,5,30]. The results of [4,5] were used to obtain a fitting function [6] which approximates the gravitational wave signal as a function of the nucleation temperature (T n ), the inverse transition duration (β), the bubble wall speed (v w ), the sound speed (c s ) -which is obtained from the pressure, and usually assumed to be c s ∼ 1/ √ 3 -and the kinetic energy fraction, or energy budget, K. The value of K can be determined from the fluid velocity profile of an isolated bubble, but is usually obtained from a fit provided by [31], which gives the value of K as a function of the phase transition strength α and v w .
The results of [31] were obtained in the 'bag model', which assumes that both the broken and the symmetric phase are composed of radiation degrees of freedom only. The same assumption underlies the choice of c s = 1/ √ 3. In a more realistic model, these assumptions do not hold, as some degrees of freedom are usually massive around the phase transition temperature. The procedure sketched above thus neglects some of the model-dependence of the signal, by describing it fully in terms of T n , α, β and v w . This leads to an inaccuracy in the estimated gravitational wave signal, but also possibly to an unnecessary degeneracy between different models.
Recently, it was demonstrated in [32,33] that further model-dependence of K can, to a good approximation, be captured by a modification of the phase transition strength (see Eq.(2.12) below) and the values of the speed of sound of both phases of the phase transition. A code snippet provided by [33] serves as an alternative to the fit based on the bag model by [31]. It was demonstrated for a set of toy models, that the code snippet approximates the numerically obtained energy budget with an error of 5%, compared to an error of 50% for a mapping onto the bag model.
Although some studies of the effect of the sound speed on the shape and strength of the gravitational wave signal now exist [32][33][34][35], an accurate computation of the value of the speed of sound in a first order phase transition is still lacking. Careful studies of the sound speed thus far have been limited to QCD [36,37] and the SM electroweak crossover [38]. The goal of this work is to improve the accuracy of the computation of the pressure -for a representative BSM model that accommodates a strong electroweak or dark sector phase transition [39] -and in particularly to determine what are realistic values of the sound speed. To understand the importance of our analysis, we will compare the resulting amplitude of the gravitational wave spectrum with the result obtained in the approximation c s = 1/ √ 3. In particular, we will study the SM extended by N singlets with O(N ) symmetry as a proxy for first order electroweak phase transitions with an extended Higgs sector. We will also investigate the effect of the removal of some of the light SM fermions, mimicking a dark sector with a smaller number of degrees of freedom.
The reason that we include multiple inert singlets is that one can expect deviations in the speed of sound when additional massive particles are present. In particular we expect that the deviation from c 2 s = 1/3 increases with increasing N . In addition, we also manually remove some of the light degrees of freedom by reducing number of fermion families (N f ) to one instead of three. In this case, one can expect the sound speed to drop, because the many light degrees of freedom all push the speed of sound (squared) towards 1/3. We indeed observe such behaviour, and can conclude that deviations in the speed of sound can be relevant for dark sector phase transitions.
BSM model While the features that we discuss are generic for thermal phase transitions in BSM scalar extensions, practically we work in the model where the SM is augmented with an O(N ) symmetric scalar singlet [40]. The Langrangian reads where the singlet decomposes as S = (S 1 , ..., S N ) T such that each component S i is real. We will focus on a two-step phase transition. We assume that in the first step of the phase transition only one of the components of the singlet gets a non-zero vev and participates dynamically in the phase transition, while the other components remain inert. Therefore, for N = 1 this setup matches the Z 2 -symmetric xSM studied in [20,[41][42][43] in the context of dimensionally reduced effective field theories. We generalise the results of these works by including the inert singlet contributions for general N . As another novelty, we also compute the coefficient of the unit operator that describes the hard and soft scale contributions to the symmetric phase pressure. These contributions were not discussed in earlier singlet dimensional reduction literature as therein the intrest was purely in the pressure difference between different phases, and corrections to the speed of sound were not determined.
Outline of this work This article is organised as follows. In Sec. 2 we review the setup for gravitational wave production, paying attention to the role of the sound speed. In Sec. 3 we discuss the computation of the pressure and inclusion of several higher order corrections. In Sec. 4 we present our results for quantities important for the GW wave power spectrum. In Sec. 5 we discuss how our findings could generalise to phase transitions in a dark sector. In Sec. 6 we discuss our conclusions and their impact. Appendix A collects multiple technical details of our analysis.
2. Model dependence of the gravitational wave signal: pressure and speed of sound For a BSM model with a first order phase transition, the gravitational wave signal from sound waves can be estimated as [6] dΩ gw,0 where F gw,0 accounts for the redshift from the generation of the signal until now, H * is the Hubble parameter at the phase transition temperature, R * is the mean bubble separation at the moment of collision, given by R * ∼ (8π) 1/3 v w /β. c s is the speed of sound in the plasma, Ω gw is a numerical factor, C is a function that determines the shape of the spectrum and f p,0 is the peak frequency. τ s , the lifetime of the source, is estimated as [44] where τ sh is the time of shock formation. 1 When shock formation is relatively slow, the gravitational wave spectrum scales with K 2 , and for fast-developing shocks the spectrum scales as K 3/2 . The kinetic energy budget K is the fraction of energy that is available for the production of gravitational waves. It can be determined by solving the hydrodynamic equations for the plasma around a single expanding bubble. We will now give a summary of this computation, for a more complete description, see e.g. [31,[46][47][48]. The hydrodynamics follow from the energy-momentum tensor of the plasma: where u µ is the four-velocity of the plasma and p and w are the pressure and enthalpy respectively. The pressure is equal to the negative of the free energy of the plasma, and the enthalpy and energy density e can be obtained from the pressure via The fluid equations are simply the continuity equations ∂ µ T µν . By projecting these along the directions parallel and perpendicular to the fluid flow, and by introducing the radial coordinate ξ = r/t, with r the distance from the center of the bubble and t the time since nucleation, the following hydrodynamic equations are found where v(ξ) is the fluid velocity, γ the Lorentz factor and µ(ξ, v) the boosted velocity The speed of sound can be derived from the pressure via The shock formation time can be estimated as [5,45] whereŪ f is the enthalpy-weighted root-mean-square fluid velocity and Γ ∼ 4/3 is the adiabatic index.
We compute the sound speed separately for each metastable phase. 2 The boundary conditions at the bubble wall are given by where the +(−) sign denotes a quantity right in front of (behind) the bubble wall. The subscript s(b) refers to the equation of state of the symmetric (broken) phase. The boundary condition for w is that w = w n in the region in front of the bubble where the fluid is at rest. We use the subscript n to denote quantities evaluated at the nucleation temperature T n . The kinetic energy fraction K is obtained from the solution of the hydrodynamic equations via where v w is the bubble wall velocity. Although both the fluid equations (2.5) and (2.6) and the boundary conditions (2.9) and (2.10) depend on the equation of state, Refs. [32,33] show that this dependence can, to a good approximation, be captured by the following parameters only: the speed of sound in the broken phase c s,b , the speed of sound in the symmetric phase c s,s (both evaluated at the nucleation temperature) and the phase transition strength and Dθ(T n ) = θ s (T n ) − θ b (T n ). (2.12) Consequently, instead of solving the full hydrodynamic equations for each BSM model, K can be obtained as a function of c s,b , c s,s , v w and αθ in a simplified model using the code snippet provided by [33]. This procedure is a generalisation of the fit in terms of the phase transition strength and v w provided by Ref. [31], using the bag equation of state. The new approximation works particularly well when the temperature-dependence of the sound speed is weak, but approximates the full numerical solution to an accuracy of < 5% for all toy models considered in [32,33]. To quantify the importance of an accurate computation of the sound speed, in Sec. 4 we will determine the kinetic energy fraction for our benchmark point using three different methods (these methods are identical to methods 1-3 in [33]): • Method 1: full solution of Eqs. (2.5) and (2.6), using the complete, temperaturedependent speed of sound.
• Method 3: assuming the speed of sound is c s = 1/ √ 3. This corresponds to a mapping onto the bag equation of state. Note that this is the most common method in the literature, and it was also used in Ref. [50], where the uncertainty in the GW spectrum associated to the other thermodynamic quantities was estimated.
The difference between Method 1 and Method 3 gives a measure of the importance of a careful computation of the speed of sound. The difference between Method 1 and Method 2 gives a measure of the temperature-dependence of the speed of sound.
In this work we will assume that the gravitational wave spectrum is fully described by the broken power law of Eq. (2.1), and only investigate the effect of our improved sound speed computation on the overall amplitude. According to Eq. (2.1), the position of the peak is determined by the size of the bubbles at the moment of collision. However, results of the sound shell model [51,52] and hybrid simulations [53] suggest that the width of the sound shells also gets imprinted onto the spectrum, resulting in a doubly broken power law shape. The results of [51][52][53] were all obtained under the assumption that the sound speed is c s = 1/ √ 3. Since deviations on the sound speed affect the shape of the fluid profile [33,54] it is expected that the shape of the gravitational wave spectrum is also affected by variations in the sound speed. This effect was demonstrated for the sound shell model in [35], but will not be considered further in our analysis.

Computation of the pressure and sound speed
In this section we review how to compute the pressure perturbatively in high temperature field theory. We present this discussion for the SM augmented with N singlets with O(N ) symmetry, but this discussion readily generalises to other BSM theories. In fact, many of the computations presented in this work can be obtained with the automated DRalgo package for generic, user defined models [55]. Many details of our computation in the model in question are relegated to Appendix A.
Using dimensional reduction to determine the phase transition parameters In order to improve the accuracy of the computation of the sound speed, we will use the method of high temperature dimensional reduction [56,57]. Therein, the perturbative computation is organised in the effective field theory language [58,59], based on a chain of scale hierarchies at high temperature. In the rest of this article, we will refer to the different scales as hard (πT ), soft (gT ) and ultrasoft (g 2 T ), where g is a weak coupling, often identified with the weak gauge coupling in EW theories. 3 This method was recently applied in [50,60] to obtain higher order corrections of T n , α and β. Ref. [50] reported an alarming leftover renormalisation scale dependence plaguing the conventional analysis based on a mere one-loop determination of the effective potential -that describes the free energy and pressure of the plasma -which in general signals slow convergence of perturbation theory. In particular, there could be a multiple order-of-magnitude uncertainty in the predicted amplitude of the GW power spectrum, location of the peak frequency, and the consequent signal-to-noise (SNR) ratio for experiments such as LISA. In [60], it was reviewed and discussed in detail why conventional one-loop analyses fail to provide better accuracy: due to slower convergence of the perturbative expansion at high temperatures, several two-loop contributions have to be included to enable leading renormalisation group (RG) improvement. After such RG improvement, order-by-order in perturbation theory, the renormalisation scale-dependence cancels between the running of lower order contributions and explicit logarithms at higher orders. The inclusion of these two-loop level contributions is standard in a three-dimensional effective field theory (3d EFT) description of the phase transition thermodynamics, from which equilibrium properties of the transition, such as the pressure, and subsequent quantities of interest such as the critical temperature, latent heat, and the speed of sound can be derived.
In this work at hand, we focus on the computation of the pressure and the speed of sound, and investigate the importance of its higher order corrections to the energy budget K, and hence to the GW power spectrum. In practice, we compute the pressure by following [38,61,62], and extend the SM computation of these works to a BSM setup. For a similar, related discussion of QCD and the speed of sound in a quark-gluon plasma, see e.g [36,37] (c.f. also the review [63]).
Although a necessary improvement in the accuracy of equilibrium quantities can be made by including higher order corrections to the pressure, or the effective potential, it was concluded in [60] that the most limiting source of theoretical uncertainty in the phase transition parameters is the bubble nucleation rate. This rate is used to derive T n and all other parameters have to be evaluated at this temperature. In [60], the bubble nucleation rate was determined only at leading order within the 3d EFT. However, recent work has described how to extend such a computation to higher orders [64][65][66][67] (also c.f. [68][69][70]). In the work at hand, we will not focus on such an improvement for the bubble nucleation rate, but instead treat T n as an input parameter that we vary; once its computation at higher orders is available, our results for the pressure, sound speed and the kinetic energy fraction can be used to determine the actual effect on GW power spectra.
Pressure in perturbation theory In perturbation theory, the pressure admits the following form and subsequent formal expansion 4 in the high temperature limit: where p sym describes the pressure in the symmetric phase where all scalar background fields vanish and V eff is the effective potential as a function of the background fields, and a-f represent constants. The MS mass parameters are assumed to scale as m 2 ∼ (gT ) 2 .
We comment below on our conventions for the normalisations of both terms. A formal expansion in g follows when a power counting rule is associated to each coupling. In particular we assume that all quartic couplings λ scale as g 2 , and we assume that there is a barrier present in the leading order effective potential; for radiatively generated transitions, the corresponding expansion is different [82]. 5 Compared to the perturbative expansion of the effective potential at zero temperature, the structure of the perturbative expansion of the pressure is inflicted by several finite-T peculiarities: odd powers of g appear, logarithmic terms include ratios of various thermal mass scales, and crucially the O(g 6 ) term is not even attainable in perturbation theory, due to the Linde problem [83]. In addition, the coupling expansion misaligns with the loop expansion due to the enhancement of IR contributions.
We emphasise, however, that in practice the pressure is not computed directly order-byorder in the above expansion in g. The reason is, that the pressure is not directly computable in such an expansion, since the direct use of perturbation theory is inhibited by various IR singularities. These have to be tackled by thermal resummations, based on a scale hierarchy of the different mass scales (hard/soft/ultrasoft) at high temperature. The perturbative result will include logarithms of all three mass scales in the problem, and it is not possible to choose a UV cutoff that would simultaneously remove all large logarithms. In particular, these logarithms are of the form ln(µ/T ), ln(T /m E ) and ln(m E /m M ), where m E ∼ gT and m M ∼ g 2 T are the soft and ultrasoft -or electric (E) and magnetic (M) -mass scales. These logarithms lead formally to orders g 4 ln(g) and g 6 ln(g 6 ), where logarithms are large for g → 0.
Dimensional reduction to a chain of EFTs resolves this problem systematically: the thermal contributions of the hard scale are resummed as contributions to the parameters of the soft and ultrasoft scales. The running within the EFT introduces new RG scales (µ 3 andμ 3 ), which can be chosen independently, allowing to avoid large logarithms. In our computation, we associate p sym with the coefficient of the unit operator [59] in dimensional reduction, and this contribution describes the hard and soft contributions to the symmetric phase pressure.
Contributions of the ultrasoft scale are, by our convention, encoded in V eff , for both symmetric and broken phases.
Let us take a closer look at how different order terms in the pressure arise, and how we label them: Note that terms that are formally of higher order are resummed inside the 3d EFT parameters in different contributions and, in accord with the EFT construction, we will keep them untruncated in the final result [36,84]. We remind that the O(g 6 )-contribution is non-perturbative, although the O(g 6 ln g) logarithmic term is still computable at 4-loop level [79,80].
The LO result for the sound speed, c 2 s = 1/3, follows from the LO result for the pressure. Typical phase transition analyses, based on the one-loop daisy resummed thermal effective potential (c.f. Appendix A.6), produce a pressure that is correct to NNLO at O(g 3 ) for the broken phase contributions, but are lacking two-loop O(g 2 ) terms in the symmetric phase. Also, such analyses include some O(g 4 ) corrections, but are not complete at N 3 LO. In our computation in this article at hand, we provide complete accuracy at N 3 LO: this requires a three-loop computation for the hard contributions, and a two-loop computation of the soft contributions and the broken phase effective potential; technical details can be found in Appendix A. At N 4 LO the O(g 5 ) correction would require a three-loop computation in the 3d EFT, which we do not include here, but leave to future work. It is tempting to pursue this correction, as it is the last term that can still be obtained in perturbation theory, c.f. [79].
Numerical study For our numerical analysis, we select a single, representative benchmark point with a two-step phase transition [85][86][87][88][89] with a strong second transition to the Higgs phase.
BM: (m S , λ S , λ m ) = (160 GeV, 1.0, 1.6), (3.2) where m S is the singlet mass. The order parameters of the different phases are depicted in Fig. 1. The Higgs (singlet) background field is denoted by v (s). The pressure and sound speed -computed from Eq. (2.8)  Figure 1: Minima of the real part of the effective potential as a function of the temperature, determined in a one-loop approximation (left) and at two-loop within the 3d EFT with NLO dimensional reduction (right). Blue (red) denotes the Higgs (singlet) direction, and dots (crosses) present global (local) minima. Local minima describe the metastable phases. The different dark and light colours present two choices for the 4d RG scale µ: 0.5πT and 2πT . The vertical lines depict the critical temperatures for the second transition (T c,φ ), and the band between them the theoretical uncertainty in their determination. We note that the critical temperatures for the first transition (T c,S ) is very sensitive to two-loop corrections, while T c,φ is not.
-are shown in Figs. 2 and 3, respectively. In all these plots, we present a comparison of our full N 3 LO computation with the frequently used, sole one-loop approximation [90] (cf. Appendix A.6). In addition, we vary the RG scale µ from 0.5 πT to 2 πT , to monitor convergence, and we find in all the aforementioned figures that in the full N 3 LO computation the sensitivity to the 4d theory RG scale is milder, as expected. For the 3d RG scales we have Figs. 2 and 3. We note that our use of minima of the effective potential as order parameters leads to an undesired gauge-dependence that signals an inconsistent perturbative expansion. We carefully comment on this issue in Appendix A.5 and argue why this does not compromise our discussion regarding the sound speed.
We find that the deviations in the speed of sound squared from c 2 s = 1/3 are small in the singlet phase, where they remain subpercent over the entire temperature range for the one-loop result. Deviations are somewhat more apparent in the N 3 LO analysis, of the order of 1%. Note that the regime where the deviations are largest is the regime with most supercooling, where the high-temperature expansion might not be reliable. In the Higgs phase, the one-loop result deviates from the relativistic value by ∼ 2% over a large range of temperatures. Again, the deviations in the N 3 LO analysis are larger, ∼ 3 − 4%. In both cases, the speed of sound exceeds the relativistic value in the less trustworthy regime of strong supercooling.

Impact on gravitational wave predictions
We have seen in Sec. 2 that the gravitational wave spectrum depends on the sound speed through the kinetic energy fraction K and through an explicit appearance in Eq. (2.1). Fig. 3 demonstrates that, for N = 1, the deviations from c s = 1/ √ 3 are moderate, so we do not expect major deviations in the gravitational wave spectrum from the explicit dependence on c s . As demonstrated in [33] however, a small deviation in the sound speed can already lead to a significant effect on the kinetic energy budget. Bear in mind that the gravitational wave amplitude scales with K 2 or K 3/2 , which amplifies the dependence. We will quantify the effect of the accurate computation of the sound speed on the gravitational wave spectrum by comparing the kinetic energy budget computed in the bag model (Method 3) with the computations accounting for the deviations in the sound speed (Method 1 and 2). In all following plots, we demonstrate the results obtained in the 3d EFT with fixed RG scales . It should be noted that an accurate prediction of the gravitational wave spectrum depends on the nucleation temperature, which enters in H * , R * , f p,0 and K. A self-consistent determination of the nucleation temperature requires a computation of the bubble nucleation rate within the 3d nucleation EFT [65] (for applications, see [64,69]). This goes beyond the scope of this work. Instead we will present our results for K as a function of an undetermined nucleation temperature.
One distinguishes three types of hydrodynamic solutions, characterised by the wall velocity: (subsonic) deflagrations, in which the fluid is at rest inside the bubble, and forms a shock wave in front of the bubble wall; (supersonic) detonations, for which the fluid in front of the bubble is at rest and a rarefaction wave forms inside; hybrids, or supersonic deflagrations, which consist of a rarefaction wave and a shock. We treat the wall velocity as an external parameter, but in reality, the computation of its value is challenging. Recently, it has been argued that solutions likely fall into one of two categories: deflagrations and ultrarelativistic detonations, but that no solutions exist between the Jouguet velocity and wall velocities with γ w 10 [91,92].
We have determined the ratio of the kinetic energy fraction K computed in Method 2, K 2 and Method 3, K 3 , to the result in Method 1, K 1 for our benchmark point. The results are demonstrated in Fig. 4 for a hybrid and detonation solution. The temperature has been rescaled by the critical temperature of the transition to the singlet to the Higgs phase. We found that for the deflagration solution, the differences between the three methods are very small (< O(1%)) for the entire temperature range. This implies that the gravitational wave spectrum can be estimated accurately with c s = 1/ √ 3 and a mapping onto the bag equation of state. For a hybrid, the story is different. Especially when the nucleation temperature is relatively close to the critical temperature, 6 the difference between K 3 and the other methods is more pronounced, approaching O(40%). The estimate based on the bag equation of state typically overestimates the kinetic energy fraction. Over the entire temperature range, K 2 reproduces Method 1 very well, meaning that the temperature dependence of the sound speed plays an insignificant role, and a full numerical solution of the hydrodynamics is not required. For detonations, the effect of the sound speed is again smaller, with a maximum of O(8%) for T n ∼ T c,φ . For significant supercooling, we see that K 2 starts to deviate slightly from K 1 , implying that the temperature-dependence of the sound speed mildly affects the result. Fig. 5 showcases the absolute value of K for a hybrid solution as a function of the nucleation temperature, for the N 3 LO (black) and sole one-loop result (grey). The solid line displays K 1 , and the dashed line K 3 . We see that the difference between K 1 and K 3 is larger in the N 3 LO analysis, which could be expected from Fig. 3. We see that, for a given nucleation temperature, the N 3 LO computation yields a smaller value of K than the one-loop analysis (but of course the two approaches do not yield the same nucleation temperature). In addition, our accurate treatment of the speed of sound suppresses the kinetic energy fraction compared to the treatment with c 2 s = 1/3. For a given temperature, the determination of the pressure at one loop and a determination of K via the bag model together overestimate K by 3% in the strong supercooled regime up to 50% in the regime around 0.9T c,φ .

Speed of sound at a dark sector
We have seen in the previous section that the sound speed only deviates from c 2 s = 1/3 by a few percent. The main reason for this is that the light SM fermions give the dominant contribution to the pressure. In this section, we study how the the sound speed behaves when the particle content deviates further from the SM. We investigate two different effects: • Increasing the number of singlets N .
• Reducing the number of fermions by setting N f = 1. This modification mimics a dark sector with a smaller amount of fermions than the SM. of singlets suppresses the sound speed. Again, the deviations are largest in the broken phase, and approach 5%. The dotted lines demonstrate the sound speed in the dark sector with only 1/3 of the fermionic content of the SM. We see that this leads to an even stronger suppression in the sound speed, approaching 8% in the singlet phase and 15% in the Higgs phase. This is a significant result, as deviations of that size can suppress the gravitational wave signal by an order of magnitude [33,35] and also affect its shape. An even stronger suppression is expected if the particle content gets reduced further. From this simple exercise we thus conclude that especially in dark sectors with a particle content that strongly deviates from the SM, a computation of the gravitational wave spectrum requires a careful computation of the sound speed.

Discussion
In this work, we have computed the pressure and the corresponding sound speed at N 3 LO in formal expansion in g (c.f. Eq. (3.1) and discussion below it) in the Standard Model augmented by N singlets with O(N ) symmetry. The speed of sound enters in the computation of the energy budget of gravitational waves. Whenever the sound speed deviates from the often-assumed value c 2 s = 1/3, the amplitude and shape of the gravitational wave spectrum are modified. In order to obtain N 3 LO accuracy, we used a dimensionally reduced EFT. We worked with a benchmark point defined in Eq. (3.2) giving rise to a two-step phase transition. We expect our qualitative results to hold for other parameter choices, and other BSM models that accommodate similar strong two-step phase transitions, relevant for observable GW wave signatures.
As known from earlier comparisons between the mere one-loop analysis and higher order analysis in the 3d EFT, thermodynamic quantities such as the critical temperature and latent heat differ significantly between the two computations, potentially leading to multiple order of magnitude differences in the gravitational wave signal. The one-loop result suffers from a disturbing leftover RG-scale dependence, which gets significantly reduced in the result obtained from the 3d EFT. We observe the same reduction in the RG-scale dependence in the pressure and the sound speed. In the N 3 LO analysis, suppressions of the sound speed compared to the LO value c 2 s = 1/3 are more significant than in the mere one-loop computation. This suppression affects the energy budget, which sets the amplitude of the gravitational wave signal: in a phase transition with c 2 s = 1/3 the energy budget should be computed with the methods of [32,33] instead of the commonly-used bag model. The error associated with using the bag model or an inaccurately computed c 2 s is not as large as the error associated to some of the other parameters of the phase transition [50,60]. Nevertheless, especially in the case of hybrid and detonation solutions, the error is non-negligible, and can be as high as O(50%) for the energy budget.
In Sec. 5 we further modified the particle content, by increasing the number of singlets, and by removing 2/3 of the fermionic content. These modifications allowed us to envision the possible behaviour of the sound speed in phase transition in dark sectors. We found that both an increase in N , as well as the removal of the 2/3 of the fermions cause larger deviations in the pressure compared to its LO value. As a result, the sound speed decreases with N , and especially small values are obtained in the fewer fermion case. The effect was strongest in the Higgs phase, where the sound speed could decrease by almost 15%. Such strong deviations from c 2 s = 1/3 lead to a significant suppression in the gravitational wave signal, by possibly an order of magnitude, and also significantly modify its shape. We thus conclude that computations of the gravitational wave signal from phase transitions in a dark sector require an accurate computation of the sound speed in order to even obtain the right order of magnitude of the signal.

A. Dimensional reduction and pressure with O(N ) singlet
In this appendix we collect multiple details of our computation, in particular the dimensional reduction matching relations for the SM + O(N ) singlet, as well as the computation of the pressure. For the SM parameters, we use the same power counting in terms of g as in [20,41], and in analogy to these references we use The detailed structure and exact definitions of the 3d EFTs at the soft and ultrasoft scales are discussed in [20,41,43]. We do not repeat the discussion here, as the only new aspect is the singlet having N components instead of one. Parameters of the 3d EFTs will be denoted with an additional subscript 3, and parameters of the ultrasoft theory with an additional bar. We emphasise, that compared to [62] in which the SM pressure is computed, we use a different power counting for the 3d Higgs mass parameter. In this reference, in the context of pure SM, the proper counting is µ 2 h,3 ∼ g 3 3 , which corresponds to a radiatively generated barrier (see for a similar discussion the recent work [69,82]), whereas here we use µ 2 h,3 ∼ g 2 3 which corresponds to the case with a tree-level barrier in a two-step phase transition of the Higgs and singlet.

A.1. Renormalisation
We start with the zero T renormalisation. All the relevant counterterms are listed in [41] and [20], where the case of SM with a single additional singlet was considered. Here, we consider the case of N singlets, and list the counterterms that get modified by the presence of these additional singlets. We use the modified minimal subtraction scheme (MS).
with λ h the Higgs quartic coupling and µ h its mass parameter, g and g the SU(2) and U(1) gauge couplings, g Y the Yukawa coupling and ξ 2 and ξ 1 the SU(2) and U(1) gauge fixing parameters, in general covariant gauge, respectively.
The one-loop RG-equations, which encode the running of parameters as a function of the RG-scale µ to order O(g 4 ), read For N = 1, they agree with the RG-equations of [41], but be aware of the different convention for the sign of µ 2 h and µ 2 S . All MS parameters in the Lagrangian can be related to pole masses and other physical parameters as in [43,58,93], but in this work at hand, we merely use leading, tree-level relations without one-loop corrections.

A.2. Dimensional reduction
For the dimensional reduction procedure, the Feynman rules of the symmetric phase are listed in [41]. Here we only list the Feynman rules for the singlet. The singlet propagator reads S a (P )S b (Q) = δ abδ (P + Q) where P, Q denote Euclidean four-momentum andδ(K) where the indices a, b, c, d label singlet components, for which δ aa = N , while i, j are fundamental SU(2) indices for the Higgs. The 3d matching relations, or short distance coefficients, can be obtained following [20, 41-43, 58, 59]. 7 We list the matching relations that involve the singlet fields here at NLO, i.e. at O(g 4 ): where µ 3 denotes the RG-scale of the 3d EFT. Here we have used the shorthand notation For the mass parameter we have included the O( T 2 g 2 ) level pieces, as these contribute to the ultrasoft pressure p M2 when multiplied with 1/ poles therein. The coefficients β λm2 , γ λm2 and γ λs2 read Similar pure SM contributions to the Higgs mass parameter can be found in [62] in Appendix B, there denoted by β A2 , β B2 , β λ2 , β Y 2 . The 3d gauge parameters do not receive any singlet contributions at NLO and can be read off from [41,58]. The Debye masses for the temporal scalars read The SM parts at two-loop order can be found in Refs. [61,62]. The singlet interacts with the temporal scalar fields through [20] However, since these interactions only start at O(g 4 ), they do not contribute to any of our further expressions.

A.3. Coefficient of the unit operator, or pressure, in the symmetric phase
In addition to the matching of the couplings and masses, we also need the matching of the unit operator (c.f. [59] Sec. 3A), which amounts to computing vacuum contributions to the pressure from the hard and soft modes, as in [61]. Note that by convention, we include the similar ultrasoft pieces of the symmetric phase in the effective potential.
Leading order pressure At leading order, the pressure is given in terms of the sumintegrals [77] arising from hard mode contributions of one-loop vacuum bubble diagrams. For the definition of the sum-integration measure, see [77]. Following [61], different fields contribute to the total pressure as To label different contributions, we have used the following notation of [61] for the constants: N f = 3 for the number of fermion families, N c = 3 for the SU(3) colour group, d F = 2 for the dimensionality of the fundamental representation of SU(2) and d A = 3 the dimensionality of the adjoint representation of SU (2). In total, This is the leading order pressure in the symmetric phase.
Higher order corrections The full pressure in the symmetric phase has the form Here, different contributions are organised in a manner similar to [61,62]: p E collects the hard mode contributions, and p M1 the soft mode contributions from temporal scalars. The subscripts denote electric and magnetic contributions, as in hot EQCD [77]. We note that by convention, the leading order QCD contributions are included in p E , while higher order QCD corrections are collected in p QCD [75][76][77][78][79]. We do not include these QCD corrections in our computation, as consecutive orders are known to fluctuate around the ideal gas pressure -due to the large values of the strong coupling g s and g Y -unless the temperature is asymptotically large, and here our focus is solely on the EW corrections. In addition, for p sym there would also be an ultrasoft, or additional magnetic contribution p M2 (T ) [61,62] that arises from one-and two-loop diagrams with ultrasoft scalars. However, by convention we include these contributions in the effective potential part of the pressure, i.e. p M2 (T ) = −T V 3d eff (0, 0) (c.f. Sec. A.4), where both Higgs and singlet background fields vanish in the symmetric phase.
To reach the N 3 LO, or O(g 4 ) precision for the vacuum diagrams requires a computation up to -and including -three-loop topologies, for the hard scale contributions. These can be split into where the SM contribution can be read off from Eq. (20) and Appendix A in [61]. 8 In this work, we include the singlet contributions where we used notation similar to [61,62].  As in [62], we compute p M2 (T ) at two-loop order, providing O(g 4 ) accuracy. Finally, we comment that the divergent 1/ poles in p E , p M1 and p M2 cancel when all these terms are summed together, i.e. the cancellation happens between contributions from different scales. 9 Note that p singlet M2 is finite and p M1 does not have divergences related to the singlet contributions. Therefore all leftover divergent singlet pieces from the hard modes are cancelled by the ultrasoft Higgs-gauge field contribution that is proportional toμ 2 h,3 and that encodes the contribution from the singlet portal coupling.
Vacuum Feynman diagrams at the hard scale with the singlet The computation of the singlet contributions depicted in Fig. 7 is relatively straightforward -even at three-loop level -compared to the computation of the SM contributions: gauge, ghost and fermion topology diagrams do not appear at all at two-loop, and even at three-loop there are only three of them in total. Furthermore, all these reduce trivially to lower order master integrals, since the momentum in the singlet one-loop bubble separates from the other momenta. All other three-loop topologies with the singlet are pure scalar diagrams of the singlets and Higgs, and therefore straightforward to compute. As in Appendix D of Ref. [61] for the SM contributions, we list the explicit midstage results for all different new diagrams containing singlets. Note that here we show also pure Higgs diagrams that include counterterms with singlet couplings.
The results and definitions for the master integrals I 1,2 , I 1 and M 0,0 can be found in Appendix D of Ref. [61] and Appendix A of Ref. [77]. For diagrams with gauge fields, we have used the Fermi gauge, or general covariant gauge, with the identical choice for the SU(2) and U(1) sector gauge fixing parameters ξ 2 = ξ 1 = ξ. It can be seen how the gauge fixing parameter ξ cancels for the sum of diagrams (m)+(n). Note that the SM parts in [61] were computed in Feynman gauge, and the authors did not explicitly check the gauge invariance of their final result for the pressure. We note that it would be a valuable crosscheck of the correctness of the computation to confirm the gauge invariance of the p SM E (T ), and for this the automated computation developed and used in Refs. [20,50,94] could be utilised.
From our computation in Eqs. (A.42) -(A.56), we find the following results for the different contributions from the hard modes to Eq. (A.38) We have used similar notation to that of [61], and L b was defined in Eq. (A.22). In analogy to [61], the normalisation p(T = 0) = 0 in the symmetric phase is accounted for by α Eµ 4 S . In practice, we have subracted the zero temperature Coleman-Weinberg part from the hard mode contribution Here, the first term is the contribution from the hard modes and the second term is the T = 0 contribution, c.f. e.g. Eqs. (A.22) and (2.6) in [60]. Note in particular that the dependence on the RG-scale vanishes in α Eµ 4

S
. The divergence related to µ 4 S is removed by imposing a vacuum counterterm, that reads, together with the similar contribution of the Higgs This completes our computation of p singlet E .
A.4. The effective potential, or pressure, in the broken phase The pressure in the broken phases, with non-zero vacuum expectation values for the Higgs or singlet, can be decomposed as where the effective potential is evaluated at its minima (v min , s min ). Note that according to our EFT construction, the effective potential is computed in the final ultrasoft scale EFT. The parameters, background fields and the ultrasoft RG scale (μ 3 ) of this EFT are denoted by bars. For the xSM the effective potential has been computed in [43], in Landau gauge, and here we generalise this computation for N singlets used in this work at hand. We compute the effective potential assuming only one singlet component gets a vev at high temperature, i.e. S → (S 1 +x, ..., S N ) T . S 1 mixes with the Higgs field, and they constitute two mass eigenstates h 1 and h 2 as in the xSM case, where N = 1. The mixing angle is defined as The other N − 1 inert singlets do not couple to gauge fields, and they merely have interaction with h 1 , h 2 , Goldstones and among themselves. All inert singlets have mass (squared) eigenvaluem In total, we compose the effective potential as where the indices denote loop-order. The tree-level potential is the same as in the xSM. At one-loop, the inert singlets add a contribution to the effective potential At two-loop, the inert singlets add where the master integrals can be found in the supplementary material of [88] and the vertex coefficients read where we denote cθ ≡ cosθ and sθ ≡ sinθ. The 3d EFT counterterms read  Figure 8: Schematic illustration of order parameters in perturbation theory, as a function of temperature. Left: minimum of the effective potential (v). At high T , v is zero, but non-zero at low T . If there is a discontinuity, the transition between the two phases is of first order (A). If v is continuous, but there is a discontinuity in the first derivative, the transition is of second order (B). A smooth crossover (C), where all order derivatives are continuous, is not possible. Furthermore, v itself is not a gauge invariant quantity. Right: the scalar condensate φ † φ has a non-vanishing value even in the high T phase, and hence a smooth interpolation between the two phases is possible, leading to a crossover between the phases. Condensates can be computed in a gauge invariant manner, even in perturbation theory.
Note that the vacuum counterterm δV 3 does not include contributions from the singlet, but since it was not listed in [43], we add it here. The magnetic, or ultrasoft, contribution to the symmetric phase pressure is encoded in the effective potential as p M2 (T ) = T V 3d eff (0, 0).

A.5. Order parameters in perturbation theory and gauge invariance
Determining the character of a phase transition -whether it is of first or second order, or crossover -is a non-trivial endeavour. In a crossover type transition the system transitions smoothly from one phase to another, and there is no discontinuity in the order parameter at any order in its derivatives. In conventional analyses in which the minima of the potential are treated as order parameters, such a transition can not be recovered: in the high temperature symmetric phase such minima vanish, but they are non-zero and vary smoothly at lower temperatures, so there cannot be a smooth interpolation between the phases. Furthermore, these minima are notoriously gauge dependent, 10 and hence cannot serve as realistic orderparameters. More realistic order parameters can be defined in terms of condensates, such as φ † φ [95]. They have a non-zero value even in the symmetric phase and can be computed in a gauge-invariant manner in perturbation theory as , (A.94) 10 Only the value of the effective potential at the minima is gauge invariant, and corresponds to the pressure.
where the effective potential is evaluated in an expansion around its leading order minimum. Even better, the condensates can be determined in a lattice Monte Carlo simulation, which does not require gauge fixing. We point out, however, that the computation of these condensates is dependent on their UV renormalisation, and hence their direct numerical values do not have a physical meaning. 11 These order parameter-like quantities are schematically illustrated in Fig. 8.
In this article at hand, we are interested in the sound speed for a strong phase transition, and we have focused solely on the second transition of a two-stepper. In such a case, the barrier separating the phases exists already at tree-level, and the phase transition can hence safely be assumed to be of first order, even without a dedicated non-perturbative study. Regardless, the computation in perturbation theory should be arranged in a gauge-invariant manner, in order to study the physical thermodynamic properties of the plasma. Despite this, however, in this article at hand we have resorted to a naive analysis in terms of a direct minimisation of the real part of the potential, and treat the location of the minima as a first estimate of the order parameters, based on the arguments presented in [43]. With this decision, none of our results are properly gauge invariant. Furthermore, in [96] it has been demonstrated that the use of either minima in Landau gauge or condensates does not make a qualitative difference. Hence, we expect our results in Landau gauge to match the gauge invariant analysis, qualitatively. We argue, that since our main interest is to get a handle on higher order corrections to the speed of sound and estimate their effect on the GW power spectrum, the limitation of a gauge dependent analysis does not compromise the conclusion of this work. For future work, that would also include a determination of the bubble nucleation rate and the nucleation temperature, and would hence provide a more complete analysis of the predicted GW signature, we envision an upgrade of our current computation of the pressure, 12 in terms of a gauge invariant analysis in perturbation theory, c.f. recent [38,50,60,96] but also [97,98]. In addition, in such an analysis, a consistent perturbative expansion around the leading order minima, imaginary parts do not show up. The effective potential itself has a spurious imaginary part, that we have simply disregarded in our direct minimisation of the potential.

A.6. Relation to lower order computations
For readers not familiar with the 3d EFT approach, we describe here the commonly adopted one-loop computation, in order to provide a comparison. For the remainder of this section, we set N = 1, i.e. we discuss the xSM. In a one-loop computation, the EFT technology is not 11 Physical quantities such as the latent heat released in the transition, can be related to these condensates though, see e.g. [50]. 12 In our computation, the dimensional reduction step is gauge invariant [20], as well as the symmetric phase pressure [61,62], but the use of the effective potential is not. as apparent as in computations required for higher order corrections. Regardless, even the one-loop effective potential based on daisy resummation by Arnold and Espinosa [99] utilises the EFT picture for the thermal scale hierarchy, as only the soft modes are resummed, corresponding to screening of the hard scale. However, the underlying EFT picture is not apparent in [99] and the computation is performed directly in the 4d parent theory; a formal discussion in the EFT language is formulated in [58,59]. At leading order, the parameters of the EFT have trivial relations to the 4d parameters, i.e. λ 3 ∼ T λ and v 2 3 ∼ v 2 /T , and only the masses receive thermal corrections, at one-loop level. Let us start by taking a closer look at the all-order resummation of the leading daisy diagrams.
Leading order daisy resummation In Fig. 9 we illustrate schematically the premise of daisy resummation. Fig. 9 (a) showcases the division of the one-loop vacuum bubble into soft (dashed line) and hard contributions (double lines). While the hard modes are regulated in the IR by non-zero Matsubara frequencies, the soft part is IR-sensitive. Physically, this is related to macroscopic collective phenomena, i.e. screening of different scales in the plasma. Technically, a problem appears in loops including the zero-mode, since the propagator has the same form as at T = 0, while the integration measure has one less power of momentum: this leads to worse behaviour in the IR compared to the T = 0 case. This is problem is demonstrated in Fig. 9 (b): consider a k-loop daisy diagram where the inner loop has a soft loop-momentum with k propagators with unresummed mass m 2 ∼ (gT ) 2 , and k loops with hard loop momenta P . When each vertex contributes with g 2 this diagram has the form [100] Results for the used one-loop integrals can be found in e.g. [77]. Here (x) n denotes the falling factorial, but the overall numerical constant is irrelevant for the observation, that regardless of k, any such diagram contributes at O(g 3 ), since m ∼ gT . Even worse, for k ≥ 2 and m → 0 all these contributions are IR-divergent. However, for massive fields -i.e. scalars and temporal scalar fields, but not for massless, spatial gauge fields -there is a salvation by means of a resummation. In Fig. 9 (c) the soft scalar propagator is resummed as a geometric Dyson series. Formally, the scalar 2-point Green's function G can be witten in terms of the tree-level propagator G 0 = 1/(p 2 + m 2 ) and the one-loop thermal correction Π ∼ g 2 T 2 as which leads to a thermally corrected mass for the three-dimensional zero-mode m 2 3 ≡ m 2 + Π, which is resummed to all-orders, by the one-loop hard correction. When all different daisy diagrams in Fig. 9 (b) are summed together, the result for the one-loop soft vacuum bubble reads When this contribution is combined with the tree-level terms and thermally corrected mass at O(g 2 ), the effective potential is consistent at O(g 3 ). In Fig. 9 (a) in addition to the O(g 3 ) soft contribution, there are one-loop hard mode contributions -that include the T = 0 Coleman-Weinberg potential (see end of this section) -which appear at order O(g 4 ). Alas, for a complete computation at this order, also higher order resummations are required. Such higher order resummations are those depicted in they are resummed to all-orders in dimensional reduction to 3d EFT, in analogy to daisy resummation that resums one-loop thermal mass to all orders.
Comparison of 3d and 4d effective potentials Schematically, the relation between the lower-order 4d effective potential and the higher-order 3d effective potential is 13 This relation is only an approximate equality, since V 3d eff (v 3 , s 3 ) contains higher-order contributions. Technically, the one-loop V 4d eff with daisy resummation is correct at O(g 3 ), but incomplete at O(g 4 ). By comparing the above expressions, we have for the soft terms In practice one uses O(g 4 ) accurate relations in the 3d EFT, providing higher accuracy. We emphasise, that this is exactly the term that corresponds to leading daisy resummation. For the hard terms we have Again, equality holds exactly if the expressions are computed to the same order in the coupling expansion. In the 3d picture, hard contributions are encoded in the parameters of the EFT.
Here we emphasise, that in the 4d approach, the leading contributions that depend on the external momenta -corresponding to a screening of the soft field by the heat bath -are not included. In the 3d computation these terms are included, and necessary for a consistent computation at O(g 4 ) and RG-improvement, see Ref. [60,99,101]. As discussed in these references, at O(g 4 ) also two-loop contributions to the thermal masses are required, as well as a two-loop effective potential for the soft or 3d terms. For concreteness, below we construct again the one-loop effective potential using the background field method. We need the mass squared eigenvalues for fluctuating quantum fields where m ± are the two neutral scalar mass eigenvalues, and the Goldstone mass eigenvalue m 2 χ is triply-degenerate. For the gauge field and top quark pieces, we need the mass eigenvalues The mass squared eigenvalues related to the temporal scalar fields read [58,101] where M 2 3 is doubly-degenerate. Here all quantities are those of the 3d EFT, since temporal scalars live at the soft scale. 14 Coefficients h 3 , h 3 , h 3 are couplings between Higgs and temporal scalars, see [20,58].
14 Temporal scalars A a 0 and B0 are mixing through a h 3 term, but in practise one can linearise this mixing.
Above, the mass eigenvalues are input for the one-loop master integral V 4d 1-loop ∼ where we have used the high-T expansion in m 2 /T 2 for the last line. Concretely where we give the result for both bosonic and fermionic fields. Typically only the first three terms are kept. Higher order terms would result in higher dimensional, marginal, operators in the 3d EFT [58].
In total, the one-loop effective potential reads (in Landau gauge)  The counterterms cancel the 1/ poles in J hard . In the soft pieces, the conventional daisyresummed approach corresponds to utilising the dimensional reduction matching relations at leading order, i.e. one-loop for the masses and tree-level for the couplings This completes the construction of the one-loop effective potential. Note that the zero temperature Coleman-Weinberg potential is implicitly included above, as can be evaluated numerically without using the high-T expansion, while the resummation of the zero mode is encoded in J daisy (m 3 , m). However, the zero mode is resummed due to screening by hard non-zero modes, and still assumes the high temperature scale hierarchy. Therefore, in this formulation it is not straightforward to move to a low-temperature limit even if the thermal function is handled numerically, without high temperature expansion. where a b = (4π) 2 Exp( 3 2 −2γ) and a f = a b /16. In the high-T expansion it is straightforward to check that indeed J CW + J T = J hard , i.e. the zero-T Coleman-Weinberg potential is included in J hard .
Finally, we comment that above we did not yet include #T 4 contributions from gluons and light fermions in the effective potential, but this can be done as in Sec. A.3 for the leading order pressure.