Turn up the volume: Listening to phase transitions in hot dark sectors

Stochastic gravitational wave (GW) backgrounds from first-order phase transitions are an exciting target for future GW observatories and may enable us to study dark sectors with very weak couplings to the Standard Model. In this work we show that such signals may be significantly enhanced for hot dark sectors with a temperature larger than the one of the SM thermal bath. The need to transfer the entropy from the dark sector to the SM after the phase transition can however lead to a substantial dilution of the GW signal. We study this dilution in detail, including the effect of number-changing processes in the dark sector (so-called cannibalism), and show that in large regions of parameter space a net enhancement remains. We apply our findings to a specific example of a dark sector containing a dark Higgs boson and a dark photon and find excellent detection prospects for LISA and the Einstein telescope.


Introduction
The observation of gravitational wave (GW) signals from binary mergers has opened up a new window to the universe [1][2][3]. Not only can we expect major breakthroughs in our understanding of compact astrophysical objects, but future GW observatories will enable us to explore the early universe beyond the Cosmic Microwave Background. One of the most exciting prospects is the detection of a stochastic GW background of cosmological origin, which may arise for example from inflation, cosmic strings or a first-order phase transition [4,5].
A particularly attractive possibility is that future GW observatories may search for GWs produced by phase transitions within a dark sector [6][7][8][9]. The presence of such dark sectors is well-motivated by the need to explain the role of dark matter in the early universe and during structure formation. Since any attempts to discover dark matter in the laboratory have so far been unsuccessful, it is a plausible possibility that the dark sector only interacts very feebly with Standard Model particles. In such a case, gravitational wave signals may offer unique opportunities to study the structure and dynamics of dark sectors.
If the dark sector is not in thermal contact with the Standard Model, its temperature T DS may be different from the temperature of the thermal bath of SM particles T SM . Indeed, the larger the ratio ξ = T DS /T SM the more energy is stored in the dark sector and can be released during a phase transition in the form of gravitational wave signals [10]. In the present work we therefore focus on GW signals from hot dark sectors with ξ > 1. Such a temperature difference could be a direct result of the details of reheating [11], but it could also be generated much later, for example if heavy particles in the dark sector annihilate and transfer their entropy to the lighter degrees of freedom.
However, an often overlooked problem is what happens to the energy density of a decoupled dark sector after the end of the phase transition. If the dark sector contains any light or massless states, measurements of the number of relativistic degrees of freedom during Big Bang Nucleosynthesis (BBN) and recombination place strong bounds on the temperature ratio ξ (see e.g. [12]). In the presence of massive stable states, on the other hand, the universe would typically enter matter domination much earlier than observed. Based on this line of reasoning, it was argued in Ref. [10] that the dark sector should be much colder than the SM, which in turn places a strong bound on the magnitude of any gravitational wave signal (see also Ref. [13]).
Here we consider an alternative possibility, namely that the energy of the dark sector is transferred to the SM via out-of-equilibrium decays of the lightest dark sector particle. Such decays inject entropy into the SM thermal bath and thereby alter the expansion history of the universe. Indeed, such decays have been studied as a possibility of decreasing the dark matter relic abundance after these particles have decoupled from the thermal bath [14][15][16][17][18][19][20][21]. In a similar fashion, entropy injection leads to additional red-shifting, i. e. an effective dilution, of stochastic GW backgrounds [22].
To investigate these effects in detail, we perform a model-independent calculation of the effect of out-of-equilibrium decays on GW signals in terms of the properties and the abundance of the decaying particles. We improve upon previous studies by including the effect of number-changing processes (so-called cannibalism [23][24][25][26][27][28]) when the lightest dark sector particle becomes non-relativistic. We then apply our results to a specific dark sector model of a dark photon coupled to a dark Higgs field, which develops a non-zero vacuum expectation value (vev) through a first-order phase transition. In the parameter region where the phase transition is strongest, the lightest dark sector particle is the dark Higgs boson, which can then decay for example via a tiny mixing with the SM Higgs boson.
We estimate the resulting GW signals and the corresponding signal-to-noise ratios and find that it is possible in this set-up to produce observable signals in planned GW observatories such as LISA [29,30] and the Einstein Telescope (ET) [31]. These signals can be enhanced for temperature ratios ξ > 1 and an excessive dilution of the signal can be avoided if the dark Higgs boson decays sufficiently quickly after the phase transition. This is illustrated in figure 1, which shows an example of a GW signal produced by a dark sector phase transition, as well as its dependence on the temperature ratio ξ and the lifetime τ of the dark Higgs boson. In this example, if ξ is sufficiently large and τ is sufficiently small, the GW signal can be substantially enhanced and can lie within the projected sensitivity of LISA. These conclusions are very general and apply to other types of dark sectors as well as more refined calculations of GW signals from first-order phase transitions.
The remainder of this work is structured as follows. In section 2 we introduce the general formalism for describing a hot dark sector that is not in thermal equilibrium with the SM thermal bath. We also review the calculation of the effective potential and of the stochastic GW background arising from strong first-order phase transitions. In section 3 we then consider the subsequent evolution of the dark sector and how the transfer of entropy to the SM bath leads to a dilution of GW signals. Finally, we apply this general discussion to a specific dark sector  Figure 1. Example for the effects of the dark sector temperature ratio ξ and the dark Higgs lifetime τ on the stochastic GW background spectrum. An increase in ξ increases the transition strength and thereby amplifies the signal. A long-lived dark Higgs however injects a considerable amount of entropy into the SM bath, which dilutes the signal. The specific scenario considered here can be tested by LISA for sufficiently large temperature ratios and small lifetimes as indicated by the gray shaded power-law integrated sensitivity curve. model in section 4 and obtain the predicted signal-to-noise ratios in future GW observatories. The code used to obtain our results, described in detail in appendix B, is publicly available as TransitionListener at https://github.com/tasicarl/TransitionListener.

General formalism
We begin this section by reviewing some relevant concepts from cosmology and introducing our notation for describing dark sectors. We then briefly present the calculation of the effective potential for a dark Higgs field at finite temperatures and how this can give rise to a first-order phase transition. Finally, we discuss the resulting stochastic GW background and summarize the approximations made in the present work.

Dark sector cosmology
In a flat Friedmann-Lemaître-Robertson-Walker universe, the Friedmann equations read ,ρ tot (t) + 3 H(t) [ρ tot (t) + P tot (t)] = 0 . (2.1) Here and in the following, m Pl = (8 π G) −1/2 2·10 18 GeV denotes the reduced Planck mass. The Hubble rate H(t) works as a measure for the expansion rate of the universe and can be calculated using the total energy density ρ tot (t) of the primordial plasma. The time evolution of ρ tot (t) in an expanding universe is described by the second Friedmann equation, where P tot (t) denotes the pressure of the primordial plasma. The total energy density and pressure can be obtained by summing over the contributions from all individual particle species x, that is ρ tot (t) = x ρ x (t) and P tot (t) = x P x (t). We further introduce the volume heating rate [32]q 2) for a given particle species x. The second Friedmann equation therefore states that the total heat is conserved: To make use of the Friedmann equations, the time dependences of ρ tot (t) and P tot (t), and therefore of all individual ρ x (t) and P x (t), have to be known. The full evolution of these thermodynamical quantities for a given particle species x is encoded in its distribution function f x (t, p). If particles of the species x scatter frequently enough, they will thermalize and f x (t, p) will follow a Bose-Einstein or Fermi-Dirac distribution. The thermal bath is then completely determined by its temperature T x (t) and the chemical potential µ x (t). In the case of this so-called "local thermal equilibrium", the time-temperature relation T x (t) can be inverted and used to replace the time dependence in the previous functions by a temperature dependence.
The distribution function f x (t, p) can further be used to obtain the comoving entropy density S x (t) of a generic particle species x. If x follows a local thermal equilibrium, one finds that the second law of thermodynamics holds individually for x, that is [33] where N x (t) = n x (t) a 3 (t) is the comoving particle number density of x. Hence, the comoving entropy density S x is conserved, ifQ x = 0 and µ x (t)Ṅ x (t) = 0. This is the case, when no heat is transferred between x and other particle species and when either µ x (t) = 0 oṙ N x (t) =ṅ x + 3 H(t) n x (t) = 0. Moreover, eq. (2.3) can be used to define the entropy density s x (t) = S x (t)/a 3 (t), which implies [33] T in local thermal equilibrium.
If µ x T x , the thermal distribution functions f x (t, p) can therefore be integrated to obtain Here the substitutions u x = m 2 x + p 2 /T x and z x = m x /T x have been employed and a + (−) sign refers to a fermionic (bosonic) species x.
A handy feature of these equations is that one can introduce effective relativistic degrees of freedom, which can be used to elegantly express energy and entropy densities of thermal baths consisting of multiple particle species. Dividing ρ x (T x ) by ρ rel bos (T x ) g=1 = π 2 30 T 4 x and P x (T x ) by P rel bos (T x ) g=1 = π 2 90 T 4 x , one can define [34] In the present work, we will consider a bath of SM particles and a separate thermal bath that we call the dark sector. As the interactions between the SM particles and the dark sector are assumed to be too feeble for the two sectors to thermalize, the two baths will in general have distinct temperatures T SM and T DS [10]. Introducing the ratio ξ = T DS /T SM of these temperatures, we find that the total energy and entropy densities of the primordial plasma is given by where g SM eff,ρ (g SM eff,s ) denotes the effective energy (entropy) degrees of freedom for the photon bath. The corresponding functions for the dark sector are denoted by the index "DS". Due to the high powers of ξ that enter into eqs. (2.7a) and (2.7b), the contributions from dark sectors that are only slightly hotter than the SM bath (i. e. ξ > 1) can have a large influence on ρ tot and s tot . Hot dark sectors can therefore significantly modify the thermal history of the early universe. This effect is shown in figure 2 for two hot dark sector species (a dark Higgs boson and a dark photon, see section 4) in addition to the particles of the SM bath for ξ = 3. In this plot and the following work, we used the data for the SM effective degrees of freedom given in the ancillary material of Ref. [35].
To describe the overall evolution of the combined system comprising the dark sector and the SM bath, we have to calculate the time (i. e. temperature) dependence of ξ(T SM ). For this purpose, we can use the fact that entropy is conserved individually in the two decoupled baths. Thus, S SM = 2π 2 45 g SM eff,s T 3 SM a 3 and S DS = 2π 2 45 g DS eff,s T 3 DS a 3 are both constant and hence Here, the quantityT SM specifies the temperature of the SM bath at a point in time where the two sectors have already been decoupled. Note that g DS eff,s (T SM ) depends implicitly on ξ(T SM ) and therefore eq. (2.8) must in general be solved numerically. The general result is that the temperature ratio ξ increases when the dark sector degrees of freedom g DS eff,s decrease and that it decreases when the SM degrees of freedom g SM eff,s decrease.  The temperature evolution of the effective degrees of freedom of the SM (left) and the total system, including also a dark sector (right), assuming that the dark sector particle species follow their equilibrium distributions for all times. This dark sector consists of a dark photon with mass m A = 10 6 GeV (and three internal degrees of freedom) and a dark Higgs with mass m φ = 10 4 GeV (and one internal degree of freedom). The temperature ratio between the two thermal baths was fixed to ξ = 3 to show that a dark sector slightly hotter than the SM bath can already yield interesting new dynamics. A possible temperature dependence of ξ(T SM ) as it would arise from the reheating of either sector was ignored here.
To conclude this discussion we emphasize that while the description we have presented above is valid for large parts of the thermal history, we will also encounter out-of-equilibrium processes as soon as only the lightest species remains in the dark sector. In section 3 we will discuss in detail the evolution of this stage as well as its consequences for the observable signals of a dark phase transition.

Dark scalar effective potential
In this work we are interested in thermal phase transitions within the dark sector. These transitions can occur when a scalar field has a temperature-dependent vev [36,37]. The vev of a given field is dictated by the principle of stationary action, where the action is calculated using only the static components of the field. Setting all kinetic terms in the corresponding Lagrangian to zero, the principle of stationary action reduces to the minimization of the field's potential energy density. Since the scalar field of interest is a quantum field in a finitetemperature environment, several corrections have to be added to its tree-level potential. Including finite-temperature effects up to 1-loop order and daisy diagram contributions to the vacuum energy, the effective potential reads [38] (2.9) The first term on the right-hand side is the tree-level potential of the scalar field φ, V CW (φ) + V ct (φ) is the Coleman-Weinberg contribution and the corresponding counterterm, V T is the finite-temperature contribution to the vacuum energy density, and V daisy encodes the contributions from the resummation of the Matsubara-zero modes of bosonic ring diagrams.
Treating the ultraviolet divergences in V CW with dimensional regularization in the MS renormalization scheme and the infrared divergences of the boson modes with vanishing Matsubara frequency in the thermal corrections to V tree with the Arnold-Espinoza method [39], the contributions read Here, n x are the degrees of freedom of the fields coupled to φ, n L b are their longitudinal boson components, η x is +1 (−1) for bosons (fermions), Λ is the renormalization scale, which will be set to the tree-level vev v of φ, and C x = 3/2 are the renormalization constants for scalars and fermions, while C x = 5/6 holds for gauge bosons. Goldstone modes have to be counted in addition to the longitudinal gauge boson degrees of freedom (see Ref. [38]) and the expression m 2 (φ) + Π(T ) has to be understood as the b-th eigenvalue of the temperature-dependent mass matrix. The functions Π(T ) denote the hard Debye masses of the longitudinal gauge boson components. For a quartic tree-level potential V tree (φ) = − µ 2 2 φ 2 + λ 4 φ 4 , as we will consider in section 4, the counterterm potential is given by where the counter-mass δµ 2 and the counter-coupling δλ can be calculated using [40] (2.12b)

First-order phase transitions in the dark sector
The transition of the real part of the dark Higgs field to different vevs can occur in two different fashions: continuously or discontinuously. In the first case, the global minimum of the effective potential shifts continuously with decreasing temperature, while in the opposite case competing minima in field space occur, to which φ has to tunnel to minimize its action. As was shown in Ref. [41], the euclidean tunneling action is given by if the field is embedded in a sufficiently hot thermal bath. Imposing stationarity of the action and considering O(3)-symmetric solutions, this yields the so-called bounce equation with the boundary conditions φ(r → ∞) → 0 and φ (r = 0) = 0. The euclidean distance measure r = |x| = √ R 2 + c 2 t 2 can be understood as the radius of a single expanding bubble that reaches luminal bubble wall velocities after having nucleated with an initial radius R.
By solving the bounce equation for a given temperature, one obtains a bubble profile φ T (r). Plugging this solution into eq. For simplicity, numerical factors of O(1) for the conversion of the dark sector temperature to a SM temperature are ignored in the second logarithm, here. This equation can be solved iteratively for the dark sector nucleation temperature T n DS . We emphasize that this procedure is an numerically expensive task, as it requires the computation of the bounce action S(T DS ) at each iteration step.
For the effective potential introduced above, one can show that the hard thermal loops of the scalar field cancel the thermally induced potential barrier [10]. The same holds for the longitudinally polarized modes of gauge bosons coupled to φ. The transversal modes of coupled gauge bosons however do not obtain thermal masses and can therefore still contribute to the thermally induced barrier. As we are interested in thermally induced first-oder phase transitions, we will in the following consider the case that the scalar is complex and charged under a gauge group. In this way, once the scalar field undergoes spontaneous symmetry breaking, it will give rise to a mass of the gauge boson it is coupled to, with the longitudinal polarization corresponding to the degree of freedom encoded in the angular mode of φ. As one can always project the vev of the complex scalar field to lie on its real axis, the above discussion of the effective potential including only a single scalar component φ still applies. Since the coupled gauge boson will become massive in the phase transition, we will refer to the complex scalar as a dark Higgs field in the following, in analogy to the SM Higgs field. Accordingly, the massive φ bosons will be referred to as dark Higgs bosons.

Gravitational waves from a dark first-order phase transition
The emission of GWs in a first-order phase transition is a result of the collision of bubbles in which φ already obtained its new vev. Additional contributions to the GW spectrum come from the excitation of the primordial plasma during the collision of bubbles in the form of sound waves and magnetohydrodynamic turbulence. The GW spectrum at its emission Ω em GW (f ) can be calculated as described in eq. (19) in Ref. [10]. The quantity α is a measure of the strength of a first-order phase transition. It is proportional to the amount of energy and pressure liberated in the phase transition, which can be characterized by the difference ∆θ of the trace of the energy momentum tensor between the broken and unbroken phase [8], that is Here, ∆V eff (T DS ) < 0 denotes the difference in potential energy density of the two minima of V eff , between which φ tunnels in the phase transition. To receive the dimensionless parameter that quantifies the strength of the phase transition, ∆θ is normalized to the total energy density of the surrounding plasma of relativistic species ρ n tot = π 2 30 g tot,n eff,ρ (T n SM ) 4 , such that α ≡ ∆θ ρ n tot . (2.17) As ρ n tot scales approximately with ξ −4 n for a fixed dark sector nucleation temperature T n DS and assuming g SM,n eff,ρ g DS,n eff,ρ , the transition strength scales as α ∝ ξ 4 n as was first shown in Ref. [10]. This is the reason why first-order phase transitions from dark sectors already slightly hotter than the SM bath can potentially emit strong GW signals.
Another important quantity in our analysis arises when one instead normalizes ∆θ to the energy density of only the dark sector species, that is This strength parameter will be used to calculate the energy budget of the contributions from the individual GW sources, which should only depend on the hydrodynamics of the dark sector after the phase transition but not on the decoupled SM bath. The parameter α DS is therefore independent of the temperature ratio ξ n for a fixed dark sector nucleation temperature T n DS as opposed to the transition strength α. The inverse time-scale β/H of the phase transition can be computed using the derivative of the bounce action S(T DS ) at the nucleation:

(2.19)
A fast transition happens on a short time scale and thus leads to a large β/H, which damps the resulting spectrum. This damping is due to the almost simultaneous production of many bubbles in a fast transition, which will collide while still being relatively small. In the opposite case of a slow transition, the bubble nucleation rate is low, leading to the eventual collision of larger, more energetic bubbles. Since the small bubbles after a fast transition collide more frequently than in the opposite case, the corresponding spectrum will have its peak at a higher frequency. The ratio β/H is almost independent of ξ n , as was shown in Ref. [10]. The bubble wall velocity v w is the most intricate parameter, since its calculation requires knowledge of the diverse, highly model-dependent particle processes that can happen at the accelerating bubble wall, see e. g. Refs. [42][43][44][45][46]. In general, the collision of particles in the plasma with an expanding bubble exerts a non-negligible pressure on its moving wall. Additionally, next to the mere change of momentum of particles being reflected, there is also an additional friction term due to transition radiation by gauge bosons, which is likely to dominate over the friction from particles colliding with the wall [47]. A detailed analysis of the processes happening at the bubble wall requires the solution of Boltzmann-like equations and is a subject of current research. However, for sufficiently strong first-order phase transitions, the bubble walls will quickly reach luminal velocities. Since strong phase transitions are favorable for detectable stochastic GW backgrounds, we will focus on very strong transitions for which v w ∼ 1, and neglect the details of the bubble wall dynamics.
If α DS exceeds a certain threshold strength α ∞ , the bubble walls can accelerate continuously (the "runaway regime"), while in the opposite case there will be a terminal velocity (the "non-runaway regime"). One can show that the friction exerted by particles getting (more) massive in a first-order phase transition is approximately given by [48] The condition > P fric for runaway bubbles is thus equivalent to α DS > α ∞ ≡ P fric /ρ DS (T n DS ), which will be used as a definition for the threshold transition strength. Note that the expression in eq. (2.20) is meant to sum over all particles gaining masses in the phase transition and thus excludes Goldstone bosons [7].
The efficiency factors κ of the three contributions entering Ω em GW (f ) depend on both α DS as well as the coupling between the plasma and the bubble wall. If the bubble walls have a terminal velocity, the latent heat of the transition is rather converted into kinetic energy of the plasma than used to accelerate the bubble walls. Consequently, the contribution from bubble collisions to the GW signal are negligible in the non-runaway bubble scenario, that is κ φ = 0. The efficiency factor κ sw of sound wave contributions then follows a function κ(α DS ) that can be approximated for luminal bubble wall velocities v w ∼ 1 as [7] κ Conversely, when the bubble walls can accelerate continuously, bubble collisions are nonnegligible sources of GWs and contribute with κ φ = 1 − α ∞ /α DS . The efficiency of sound wave contributions to the GW signal in the runaway regime then reads In either case, a fraction turb of the bulk motion energy is converted into turbulence, such that κ turb = turb κ sw . Following Ref. [10] we employ the optimistic estimate turb 10 %. An accurate prediction of the spectrum of GWs generated in a first-order phase transition in the early universe for a given model can be very challenging, requiring advanced methods for the calculation of the effective potential [49,50], for the onset and duration of the phase transition, as well as for a description of bubble walls [51]. The focus of the present work is however not on the detailed dynamics of the phase transition itself, but on the cosmological evolution of the dark sector that gives rise to such a phase transition subsequent to the generation of a GW signal. We will therefore limit ourselves to the simplified calculation of the stochastic GW background shown above with the understanding that all the effects discussed in the remainder of this work would equally apply to a more refined approach. We note however that the most important uncertainties of our calculations are either negligible or lead to a conservative estimate of the resulting signals of the analyzed first-order phase transitions as we explain in the following.
As was described in Ref. [8], it is not clear whether the equations used to compute the contributions to Ω em GW (f ) from sound waves and magnetohydrodynamic turbulence in the primordial plasma can be applied to the case of two decoupled sectors. In our analysis, we will focus on transitions deep in the runaway bubble regime for which α DS α ∞ , such that κ φ κ sw and κ φ κ turb . Thus, virtually no latent heat gets transferred to bulk plasma motion such that bubble collisions provide the dominant source of the emitted stochastic GW background. Therefore, the uncertainties connected with the sound wave and turbulence production of a stochastic GW background are negligible for the presented analysis.
The description of contributions from bubble collisions to Ω em GW (f ) however relies heavily on previous semi-analytical work that utilized the envelope approximation. In this approach, it is assumed that the scalar field's stress-energy is located in an infinitesimally thin shell around the bubble wall, which vanishes when two bubbles collide. Lately, this approximation was shown to only yield slightly larger signal strengths than predicted by fully numerical simulations of bubble collisions for vacuum transitions in the runaway bubble regime and to outperform the alternative bulk flow model [52]. However, Refs. [53][54][55] find substantial deviations from the envelope approximation for strongly supercooled phase transitions. Nevertheless, these uncertainties do not affect the main findings of our work and are therefore neglected in the following.
A general review of more refined computational approaches than the one used in this work can be found in Ref. [51]. Most importantly, we used the nucleation temperature as a reference scale to define thermodynamic quantities. New findings suggest to rather employ the "percolation temperature" instead, corresponding to the point in time when a significant fraction of the Universe is already filled with bubbles of the broken phase. This especially makes a difference for very strong transitions with α 1, as in these cases the percolation temperature is significantly lower than the nucleation temperature [56]. Then, also the energy density ρ tot that is used as a normalization in the definition of α in eq. (2.17) is lower, such that α increases effectively. Our approach hence yields a good approximation for transition strengths up to about 1 and underestimates the expected signal strength for even stronger transitions.
Redshift of the GW background after its emission. After its generation, the stochastic GW background propagates freely and undisturbed until today, effectively being a form of dark radiation. 1 The expansion of the Universe, however, redshifts both its amplitude and its frequency such that today's power spectrum, Ω GW (f ), can be expressed as where f denotes the spectrum's frequency in today's units, which is shifted from its value at nucleation by multiplication with the scale factor ratio a 0 /a n , where a n (a 0 ) is the scale factor at nucleation (today) [10]. The amplitude of the spectrum redshifts like a −4 H −2 , since the energy density of radiation scales with a −4 , while the critical energy density ρ c = 3 m 2 Pl H 2 with which ρ GW has been normalized to obtain Ω GW scales with H 2 . The prefactor R is thus defined as R ≡ a n a 0 In the following chapter 3, we will see how the dilution effect by the out-of-equilibrium decay of a dark sector into SM particles can contribute to the ratio a 0 /a n and the quantity R.

Evolution and decay of a hot dark sector
In this section we consider the evolution of the dark sector after the phase transition has ended and a stochastic GW background has been produced. Since we are interested in temperature ratios ξ > 1, large amounts of energy are stored in the dark sector and need to be transferred to the visible sector in order to satisfy observational constraints. This transfer of energy implies that the entropy in each sector is no longer conserved, i. e. it is no longer possible to directly calculate the scale factor a n of bubble nucleation in terms of the SM temperature T n SM . As a result, the present-day GW signal depends not only on the details of the phase transition itself, but also on the subsequent energy transfer from the dark to the visible sector.
In principle, there are a number of different ways for depleting the energy density of the dark sector. Here, we will focus on the case that the lightest particle in the dark sector (henceforth referred to as the mediator) is unstable and eventually decays into SM particles. If the mediator has a very long lifetime, it will come to completely dominate the energy density of the universe leading to an early period of matter domination (or, if number-changing processes are efficient, of cannibal domination). When these particles eventually decay, the resulting entropy injection into the SM sector may then lead to a significant dilution of any previously produced GW signal. We define the dilution factor as the ratio of the entropy in the SM sector before and after the mediator decays: Most of the discussion in this section is very general in the sense that the dilution factors that we calculate apply to a wide range of dark sectors and any type of GW signal (or, in fact, any form of fully decoupled matter or radiation). Nevertheless, we will assume for concreteness that the mediator is a scalar boson, for example a dark Higgs boson originating from spontaneous symmetry breaking in the dark sector. In this case, 3 → 2 processes play an important role in the evolution of the dark sector energy density once the temperature of the dark sector drops below the mass of the lightest state.
For simplicity, we assume that decays of the mediator happen after this particle has become non-relativistic. This assumption implies in particular that all heavier dark sector particles have annihilated away and transferred their energy to the lightest states. We denote the temperature at which the heavier states decouple by T cd DS . Moreover, if the mediator particles are non-relativistic when they decay, we can neglect effects related to inverse decays from SM states. On the other hand, the requirement that BBN proceeds as in standard cosmology places an upper bound on the lifetime of the lightest dark sector state, which is approximately given by τ < 1 s [32].
We will begin our discussion by considering the evolution of a dark sector away from thermal equilibrium in section 3.1. The equations describing the mediator decays will be derived and solved numerically in section 3.2. The overall effect of these decays on GW signals will be investigated in section 3.3.

Evolution of the mediator energy density
Once the mediator is the only particle species remaining in the dark sector, the only mechanisms that can change its energy density are the expansion of the universe, possible numberchanging processes and the eventual decays of the mediator into SM particles. We will first consider the case that number-changing processes are negligible and then extend our discussion to include mediator cannibalism.
In principle, the Boltzmann equation for the phase space density of a decaying mediator with no other interactions can be directly integrated given the time dependence of the scale factor a(t). The energy density ρ med (t) is then obtained by integrating the distribution function over momentum space [32]. Since this procedure is numerically rather expensive, we present below a simple approximation that can be used to describe the thermal history of the mediator species from its chemical decoupling until its decay.
Following the approach presented in Ref. [19] and extending it to allow for relativistic mediators, we can approximately write where τ is the lifetime of the mediator species and ζ(t) is an appropriate function (to be discussed in more detail below) such that ζ(t) = 4/3 for relativistic mediators and ζ(t) = 1 for non-relativistic mediators. We point out that eq. (3.2) implicitly assumes that the decays of the mediator only become relevant after the mediator has become non-relativistic, because we have not included a temperature-dependent Lorentz factor to account for time dilation in the final term [32].
As initial condition we assume that ρ med is given by an equilibrium distribution with vanishing chemical potential at sufficiently early times. In practice, we start our calculation at the time when all heavier states in the dark sector have decoupled, which we denote by t cd (for chemical decoupling), assuming t cd τ . To simplify notation, we introduce the dimensionless quantities and denote derivatives with respect to θ by a prime. Eq. (3.2) then becomes In figure 3, we compare the solution of this equation to the result obtained by integrating the full Boltzmann equation, assuming radiation domination for simplicity. We find that a good overall agreement between the two curves is achieved when performing the transition from ζ = 4/3 to ζ = 1 when T DS ≈ 0.38 m med , which can be written as 2 θ ≈ θ nr ≡ 7.0 θ cd T cd DS /m med 2 . In this approximation, the mediator behaves as a relativistic particle species with ρ med ∝ a −4 for θ < θ nr , while for θ > θ nr it behaves as nonrelativistic matter and its energy density scales as ρ med ∝ a −3 . For θ 1, mediator decays become relevant such that the energy density decreases exponentially as ρ med a 3 ∝ exp(−θ).
The assumption that the mediator is non-relativistic when it decays therefore translates to the requirement θ nr 1. In this case, a precise description of the mediator energy density around θ nr is irrelevant for its subsequent evolution and hence the approximation introduced above is sufficient for our purposes.
So far, we have assumed that the only processes changing the comoving number density are decays of the mediator into SM particles. As was first argued in Ref. [23], however, this description is incomplete, because a secluded particle species can perform number-changing processes like 3 → 2 or 4 → 2, thereby reducing its comoving number density while conserving 2 We point out that θnr does not depend on T cd DS . This can be seen explicitly in the cancellation of T cd   its entropy. This leads to an unusual relationship between the energy density and the scale factor until the number-changing processes become inefficient. Since in this process the species consumes itself to keep warm (i. e. to prevent becoming non-relativistic), it is casually referred to as "cannibalism" [25]. 3 If number-changing processes are efficient, i. e. their rate exceeds the Hubble rate, the chemical potential of the particle species vanishes. This is typically the case at sufficiently early times and high temperatures when the number densities are large. In this case the equilibrium energy density ρ med and the equilibrium entropy density s med only depend on z med ≡ m med /T DS . By eliminating z med we can therefore obtain a functions med (ρ med ), wherē s med ≡ 2 π 2 s med /(g med T 3 DS ) andρ med ≡ 2 π 2 ρ med /(g med T 4 DS ). Since we know that, as long as the decay rate of the mediator is negligible, the dark sector entropy is conserved (s med a 3 = const), we can use this function to calculate the evolution of the energy density of the mediator species: We note that the function d lnρ d lns (ρ med ) is close to 4/3 for large energy densities, corresponding to high temperatures and relativistic species, and approaches 1 for low energy densities ρ med , corresponding to non-relativistic species. Hence, for the case that number-changing processes are efficient, d lnρ d lns replaces the function ζ(θ) introduced above to describe the transition from the relativistic to the non-relativistic scaling.
As the number density of the mediator species decreases, number-changing processes eventually become inefficient and the chemical potential µ med can no longer be neglected. This transition is typically quite sudden, meaning that number-changing processes are either sufficient to keep the mediators in chemical equilibrium, or completely negligible [28]. Hence, as soon as the rate of number-changing processes Γ nc drops below the Hubble rate, we can revert to the description for non-interacting mediators from above. To combine both of these phases, we can therefore define which can be used in eq. (3.5) to include the effect of cannibalism. The only remaining task is then to calculate Γ nc for a given particle physics model. In the absence of a stabilising symmetry, the dominant contribution typically arises from the 3 → 2 rate Γ 32 ≈ σ 32 v 2 n 2 med ≈ σ 32 v 2 ρ 2 med /m 2 med [25], where the thermally averaged cross section of the 3 → 2 process can be written as for a scalar mediator. If the scalar potential is written as , the effective 3 → 2 coupling is given by [25,27] (4 π α 32 ) 3 ≡ If the mediator acquires its mass through the spontaneous breaking of a symmetric potential with quartic interactions, κ 3 = √ 3 κ 4 m med holds after the phase transition. The resulting time evolution of ρ medā 3 and s medā 3 is shown in figure 4. In the violet shaded area (phase I) and the red shaded area (phase II), Γ nc ≥ H holds, such that numberchanging processes are efficient, the chemical potential is negligible and ζ(θ) follows the gradual decrease from 4/3 to 1 as described above. The transition between phase I and II corresponds to θ = θ nr . For larger θ, the mediator energy density scales as ρ medā 3 ∝ 1/ lnā [23]. The end of the cannibalism period, i. e. the transition from phase II to III is given by the condition Γ nc = H. From this point onward, the expected behavior for a non-relativistically decaying species, ρ med ∝ a −3 exp(−θ), is recovered. As one can see on the right-hand side of figure 4, entropy is conserved throughout the first two phases and only decreases when the mediators start decaying.
Comparing the left panel of figure 4 to figure 3, for which the same benchmark point (but without including number-changing processes) was considered, shows that a cannibalism phase can substantially reduce the energy density stored in the dark sector. Eventually, this reduces the amount of energy injected into the SM bath and hence results in less reheating of the SM bath, as will be shown in the next section.

Entropy injection into the Standard Model
For the discussion above we have assumed that the scale factor is proportional to √ t, corresponding to a radiation dominated era. However, this is not necessarily the case in our set-up, because the non-relativistic mediators can come to dominate the energy density of the universe, leading to an early era of matter domination. The period only ends when the mediators decay, leading to the injection of a considerable amount of entropy into the SM bath [14,33]. It is well known that such an entropy injection can have a profound impact on the abundance of a frozen-out dark matter component [19]. As we will see below, the same is true for a stochastic GW background produced before the era of matter domination. Hence, it is essential for our purposes to obtain an accurate description of the relevant effects.

Differential equations governing entropy injection
As we are interested in the process of entropy transfer from one sector to another through an exchange of energy, we need to simultaneously consider the evolution of the SM and dark sector energy densities, the SM degrees of freedom and the scale factor. If the energy density of the mediators ρ med (θ) is non-negligible compared to the energy density of SM radiation ρ rad , the first Friedmann equation becomes where we have introduced the dimensionless variable r(θ) ≡ ρ med (θ)/ρ cd med and the constants θ H ≡ 3 m 2 Pl /(τ 2 ρ cd med ) and r cd rad ≡ ρ cd rad /ρ cd med . Furthermore, we have introduced the functions , which encode the change in the SM entropy degrees of freedom and the total entropy of the SM sector, respectively. Note that we have implicitly assumed that any other form of nonrelativistic matter gives a negligible contribution to the energy density and that g SM eff,ρ (θ) ≈ g SM eff,s (θ), which are both excellent approximation for the temperature range that we will be interested in.
As discussed above, the evolution of the mediator energy density is given by eq. (3.5), which can be written as The change of SM entropy due to mediator decays is directly proportional to the mediator energy density and is given by Using the relation we can write the change of the SM degrees of freedom as where we have introducedĜ see appendix A for details. Note that, equivalently, the latter equation could be replaced by a differential equation for T SM . The differential eqs. (3.10)-(3.14) can be easily solved with initial conditions given bȳ a cd = S cd = r cd = G cd = 1. The actual particle physics properties of the dark sector are hidden in the various quantities defined above, specifically r cd rad , T cd SM and θ H , as well as θ nr and α 32 , which enter in the definition of ζ(θ) in eq. (3.7). The evolution of the mediator species as well as the scale factor and the SM entropy ratio generated during the decay are therefore fully described by these five parameters.
For a more intuitive interpretation of our results, we would like to express the parameters r cd rad , θ H and θ nr in terms of the mediator lifetime τ , the mediator mass m med , the temperature ratio ξ cd as well as the temperature T cd SM . To do so, we calculate the initial mediator energy density ρ cd med by setting T cd DS = ξ cd T cd SM and assuming that the mediator has an equilibrium distribution with vanishing chemical potential at chemical decoupling. As mentioned above, we limit ourselves to the case that the degrees of freedom of the mediator are given by g med = 1. The initial ratio of energy densities r cd rad can then be calculated directly from T cd SM , using the appropriate number of relativistic degrees of freedom. The age of the Universe at chemical decoupling t cd in units of τ is given by θ cd = (2 H cd τ ) −1 , where H 2 cd = π 2 90 g tot,cd

Numerical solution
Let us consider an example to highlight the various evolutionary stages encoded in eqs. (3.10)-(3.14). Figure 5 shows an overview of the evolution of the energy densities ρ med and ρ rad , the normalized scale factorā, the temperature of the SM bath T SM , and the amount of injected entropy S SM /S cd SM into the SM bath as functions of the dimensionless time parameter θ = t/τ . The five physical input parameters are T cd SM = 1 TeV, m med = 100 GeV, τ = 0.1 s, ξ cd = 1, and α 32 = 0.01.
We can identify six distinct stages in the evolution.
Phase I: Relativistic mediator. At the initial temperature T cd SM = 1 TeV the mediators are still relativistic and the universe is dominated by SM radiation such that the scale factor and temperature obey the well-known relations of a ∝ √ t and T SM ∝ 1/a. Since the assumed mediator mass is only slightly smaller than the initial dark sector temperature, this phase only lasts for a short period of time.
Phase II: Cannibalism. Once the mediators become non-relativistic, cannibalism processes become relevant and lead to a ρ medā 3 ∝ 1/ lnā behaviour. As discussed in section 3.1, this scaling implies a more efficient depletion of the energy density compared to the case without any cannibalistic effects.
Phase III: Non-relativistic mediator. At some point 3 → 2 processes cease to be efficient while decays are not yet relevant, such that the comoving mediator energy density stays constant. Since the energy density of the SM radiation still scales with ρ rad ∝ a −4 , the contribution of the mediators becomes increasingly important for the expansion history of the universe.
Phase IV: Early matter domination. Once the mediator density dominates the universe's energy content, we find ourselves in a phase of early matter domination. This means that the scale factor no longer scales as a ∝ t 1/2 , but rather with a ∝ t 2/3 . Therefore, we can see thatā starts to deviate from its initial time evolution (marked in light violet in the upper-right panel). As a result, the temperature of the SM bath falls off slightly more quickly (T SM ∝ t −2/3 ) than predicted in ΛCDM (T SM ∝ t −1/2 ).
Phase V: Entropy injection. As soon as the decay term becomes non-negligible, the injection of entropy into the SM bath becomes relevant, influencing the entropy ratio S SM /S cd SM from θ 2 · 10 −2 onwards. The energy density of the mediator slowly starts to decrease as a result of its decay, reheating the SM bath due to T SM ∝ S 1/4 , see eq. (3.13). Note, however, that there is no literal "reheating" but rather a decrease of the SM temperature that is less steep than T SM ∝ a −1 . As shown in reference [19], the scaling in this period is T SM ∝ a −3/8 . With the decrease of mediator energy density, the radiation energy density increases, such that the universe converges towards a radiation dominated era again.
Phase VI: Decay. The decays of the mediator continue after the end of the early matter domination. At θ ∼ 1, the two energy densities considered here are again equal and the temperature has almost reached its ΛCDM evolution, as described by the curve in light red in the bottom-left panel. This curve was calculated using eq. (3.13), but therein setting a ∝ √ t (as in radiation domination) and S = 1 (for no entropy injection). Once the decaying mediators have injected most of their entropy into the SM bath, the curve for S SM /S cd SM saturates and the subsequent evolution of the universe follows the usual picture for radiation domination.
To conclude this discussion, we note that not all phases described above are present for all parameter points. For example, if α 32 is very small, there may never be a cannibalism phase. Conversely, if α 32 is very large, the universe enters a period of cannibal domination rather than matter domination, which only ends when the mediator decays.

Dilution of gravitational waves
The out-of-equilibrium decay of the mediator can result in a considerable injection of entropy and energy into the SM bath, as shown in figure 5, where the comoving entropy of the SM bath after the dark sector decay is more than two orders of magnitude larger than before. While the standard cosmological evolution will be recovered after the mediators have decayed, there is nevertheless one main consequence of the entropy injection: the dilution of frozen-out abundances.
Here we focus on the dilution effect on the generated stochastic GW background, which can be interpreted as a form of dark radiation [22]. The relevant quantity to quantify the dilution is the scale factor ratio between bubble nucleation and today, see eq. (2.22). Using the usual entropy-scale factor relation, we obtain a n a 0 = 1 where we have used that the comoving SM entropy is separately conserved before and after the mediator decay. We therefore observe once again that an increase of the SM entropy (D SM > 1) corresponds to a larger scale factor today. The effects on the frequency spectrum of the GW are then determined by the relation corresponding to a frequency shift towards smaller values when dilution effects are present. 4 Importantly, another redshift contribution affects the amplitude of the GW signal directly, which can be included in the R factor from eq. (2.23): which takes into account the contribution of the dark sector to the total energy density (see figure 2). Using this definition, eq. (3.20) The advantage of this expression is that for temperature ratios ξ 1 we expect g tot,n eff,ρ ∝ ξ 4 and g tot,n eff,s ∝ ξ 3 , see eqs. (2.7a) and (2.7b). Hence, the ξ dependence in the degrees of freedom cancels and all effects are encoded entirely in D. Moreover, as shown in Ref. [19], the dilution factor D also saturates for large ξ cd , such that in this limit R h 2 becomes independent of the temperature ratio between the two sectors.
We are now in the position to describe the effect of the mediator decay on the stochastic GW background by investigating the dependence of the dilution factor D SM on the input parameters defined at the end of section 3.2.1. For a first impression of the impact of the different quantities that specify the dark sector and its decay, we show in figure 6 the results from scans over different planes in the resulting parameter space. More specifically, we show the dependence of D SM on the temperature of the SM bath at decoupling T cd SM , the mediator mass m med , the temperature ratio at chemical decoupling ξ cd , the effective 3 → 2 coupling α 32 and the mediator lifetime τ . 5 The parameters not varied explicitly in each panel are fixed to m med = 110 GeV, T cd SM = 175 GeV, α 32 = 3.6 · 10 −3 and ξ cd = 2. Note that we exclude parameter regions where the mediator decays are already efficient during chemical decoupling (θ cd ≥ 1) and where the approximation of non-relativistic mediator decays breaks down (θ nr ≥ 1).
In the top-left panel we consider the dependence of the dilution factor D SM on the chemical decoupling temperature T cd SM and the mediator lifetime τ . We find that the entropy injection (and hence the dilution) only becomes sizeable for a sufficiently long-lived mediator, such that θ nr 1 and there is a substantial period of early matter domination. For T cd DS > m med = 110 GeV, corresponding to T cd SM > m med /ξ cd = 55 GeV, the dependence of D SM on the chemical decoupling temperature is very mild and results only from changes in the number of relativistic degrees of freedom. For smaller values of T cd SM , on the other hand, D SM decreases rapidly. This decrease is a direct consequence of our assumption that the mediator abundance is given by an equilibrium distribution at chemical decoupling and therefore becomes Boltzmann suppressed at small decoupling temperatures. We note, however, that such a Boltzmann suppression is difficult to achieve in realistic models, and that the examples that we will consider in section 4 always correspond to T cd DS > m med . In the top-right panel we focus on the dependence of D SM on the mediator mass. We find that smaller mediator masses lead to a decrease of the dilution factor. The reason is that lighter mediators experience a longer period of relativistic and cannibalistic evolution, such that the duration of the early matter domination and thus the amount of entropy injection decreases. As before, when m med > T cd DS , the mediator is Boltzmann suppressed at decoupling, thus reducing the dilution.
The bottom-left panel shows the effect of varying the temperature ratio ξ at chemical decoupling. Since an increase in ξ cd corresponds to an increase in the energy and entropy stored in the dark sector, it is clear that D SM grows with increasing ξ cd . We note that the dilution factor changes more rapidly with τ for large lifetimes, corresponding to mediator decays during matter domination, than for smaller lifetimes, corresponding to mediator decays during cannibal domination.
The general effect of the cannibalistic era can be observed in the bottom-right panel, which describes the dependence of D SM on α 32 . A large effective 3 → 2 coupling means that the number-changing processes stay efficient for a longer period of time. During such a cannibalistic era the comoving mediator energy density ρ med a 3 decreases, such that the universe enters into the phase of early matter domination at a later point for larger α 32 . We emphasize that this is potentially a large effect: Compared to the case where cannibalism is negligible (α 32 = 10 −4 ) the dilution factor can be suppressed by a factor of a few if α 32 is large. For large α 32 , the dilution factor becomes essentially independent of α 32 , as the mediator decays before cannibalism ends.
To conclude this section, we remind the reader that large dilution factors correspond to small GW signals in the present universe. We have therefore identified two competing effects: Increasing the temperature ratio ξ increases the stochastic GW background produced during a first-order phase transition (see section 2.4) but also leads to larger dilution factors. However, it should be clear from figure 6 that ξ is not the only relevant parameter. In particular, we expect the mediator lifetime to play a decisive role in determining whether increasing ξ leads to an overall enhancement or suppression of GW signals. We will study this question in a more concrete setting in the following section.

Example: Hot dark Higgs bosons
We are now equipped with all the necessary tools to investigate a particular dark sector model, determine its phase transitions and the resulting stochastic GW background, calculate the decay of a mediator species and the consequent entropy injection, and finally obtain the present-day GW signal. In section 4.1, we will provide a description of the model that we study in the rest of this section and identify the five parameters relevant to the entire phenomenological discussion. A study of the effects occurring in the different regions of the available parameter space will be presented in section 4.2. Finally, we analyze the expected signal-to-noise ratio in our model for LISA and ET for two benchmark points in section 4.3.

Model definition
We consider a simple dark sector model given by a dark photon that arises from a new U (1) D gauge group under which a complex Higgs field Φ = (φ + i ϕ) / √ 2 is charged, while the entire SM field content is neutral. The radial mode φ is the dark Higgs boson and ϕ is the Goldstone boson contributing to the longitudinal mode of the dark photon after symmetry breaking. The underlying Lagrangian can be written as [10] with the covariant derivative and the field-strength tensor X µν = ∂ µ X ν − ∂ ν X µ , where X stands for A (dark photon) or B (SM hypercharge gauge boson). The scalar potential involving the dark Higgs field Φ and the SM Higgs field H is given by [59] where µ 2 i and λ i are the different quadratic and quartic couplings. We assume that both and λ p are sufficiently small that they do not lead to the thermalisation of the dark sector with the SM thermal bath, which also implies that they play a negligible role for the phase transition. We can therefore treat the dark and visible sectors independently at early times and study the dark Higgs potential in terms of the couplings λ, g and µ. The various mass parameters read [10] At high temperatures, the effective potential given in eq. (2.9) has a minimum at φ = 0, such that the dark photon is massless. At zero temperature, on the other hand, the minimum of the potential for φ lies at v = µ/ √ λ, leading to m A = g v and m φ = √ 2 λ v. Furthermore, the coupling relevant for number-changing processes is given by α 32 = 9/(2 1/3 π) λ ≈ 2.3 λ.
We observe that, for g > √ 2 λ, the dark Higgs boson is lighter than the dark photon and therefore the lightest particle in the dark sector after the phase transition. As we are going to see, this will be the case for all of the interesting regions of parameter space. The dark Higgs boson therefore takes on the role of the decaying mediator particle discussed in chapter 3, which is responsible for the energy transfer from the dark sector to the SM. In the following, we will treat the dark Higgs lifetime τ as an independent parameter that may be determined by λ p or some other unspecified mechanism and assume that the kinetic mixing parameter is sufficiently small to neglect dark photon decays into the SM. In this set-up the chemical decoupling temperature is approximately given by the temperature when the dark photons become Boltzmann-suppressed: T cd DS = m A > m φ . 6 Another important quantity is the temperature ratio ξ n , which is fixed at the time of nucleation and evolve to later times by tracking the relevant degrees of freedom, see eq. (2.8).
While we remain agnostic about the precise process that leads to different temperatures of the dark and visible sector, simple possibilities would be a difference of the number of degrees of freedom much earlier in the universe and/or the decay of a heavy particle species that heats the dark sector. In total, our set-up is therefore characterised by five parameters: the dark Higgs quartic coupling λ, the U(1) D gauge coupling g, the vev v, the dark Higgs lifetime τ and the temperature ratio ξ n at the nucleation time of the first-order phase transition.
In practice, we first compute the details of the first-order phase transition, i. e. we determine the nucleation temperature T n DS , the inverse time scale β/H and the transition strength α. In a second step, the dilution factor D SM is computed as described in the previous chapter with the initial temperature for the evolution given by the time of decoupling of dark photons. As a final step, we then calculate D according to eq. (3.19). Having determined all these quantities, we can then compute the relevant redshift factors given in eqs. (3.17) and (3.20) which enter the final GW spectrum in eq. (2.22), linking the entire GW evolution from the time of bubble collisions up until the present day. These calculations are performed with a modified version of CosmoTransitions [60] described in more detail in appendix B. This code is publicly available as TransitionListener at https://github.com/tasicarl/ TransitionListener.
We emphasize that a subtletly arises in this computation. For some parameter points the nucleation temperature is found to be smaller than the dark photon mass in the broken phase. In such a case the number of relativistic degrees of freedom in the dark sector changes discontinuously from four before the phase transition to one after the phase transition. 7 Whenever this happens, we set the temperature of chemical decoupling T cd DS equal to the nucleation temperature T n DS (rather than to m A ) and assume four degrees of freedom for the evaluation of the nucleation criterion in eq. (2.15). This approximation accounts for our ignorance of out-of-equilibrium effects at the bubble wall. Indeed, the naively expected particle abundances may be modified by two additional processes known as "bubble filtering" [42,43] and "bubble expansion production" of heavy particles [45,46].

Exploration of model parameter space
In order to understand the parameter dependences of our set-up, it is important to realize that three out of the five aforementioned parameters determine the requirement for bubble nucleation: g, λ and v, with the first two having the strongest influence. In the effective potential the barrier height separating the true and false vacua is set by the gauge coupling g, the quartic coupling λ determines the depth of the tree-level minimum and, finally, the vev v sets the overall temperature scale for the phase transition. Moreover, the value of v In the white area above, the potential barrier is too high to be overcome, such that bubbles cannot nucleate. In the white area below, a smooth crossover occurs in which no bubbles are emitted either. The tree-level vev was fixed to v = 2 TeV and the temperature ratio was set to ξ n = 1 to generate this figure. enters the calculations through the number of effective degrees of freedom at nucleation (see below for a more detailed discussion). The dependence of the transition parameters α and β/H is visualized in figure 7 (see also Ref. [10]). As expected, we find that increasing g corresponds to stronger phase transitions, until eventually g is so large that the tunneling rate to the new vacuum receives a great suppression and the universe remains trapped inside of the false vacuum [61]. For very small values of g, conversely, the barrier becomes so small that the phase transition happens smoothly. In figure 7 the tree-level vev was fixed to v = 2 TeV, but very similar results would be obtained for somewhat different values. The temperature ratio was set to ξ n = 1.
For the purpose of the following discussion, we are mainly interested in strong phase transitions with large α, which also corresponds to an overall slow process with small β/H. Based on the insight from figure 7, the benchmark point λ = 1.5·10 −3 and g = 0.5 fulfills this criterion and will be the focus of the subsequent discussion. We point out, however, that any other point along the upper border of the coloured band would give rise to a qualitatively similar phenomenological discussion presented in the following.
In order to consider the influence of the temperature ratio on the phase transition parameters, we show the dependence of α and β/H on ξ n in figure 8 (left) for v = 2 TeV. While β/H is insensitive to ξ n [10], we observe that even a mild increase from ξ n = 1 to ξ n = 2 boosts α and therefore the GW spectrum by more than an order of magnitudesee the discussion below eq. (2.17). A further increase in the temperature ratio (for fixed T n DS ) does not have a large influence on the GW spectrum, because the dark sector begins to dominate the total energy density. Consequently ρ n tot and therefore α become independent of ξ n . This feature is further illustrated in the right panel of figure 8, which shows that the energetic effective degrees of freedom g tot,n eff,ρ grow with ξ 4 n for sufficiently large temperature ratios such that the total energy density becomes independent of the SM temperature, see eq. (2.7a). In a similar manner, g tot,n eff,s 4/3 grows with ξ 4 n for large enough ξ n as can be deduced from eq. (2.7b), in line with our reasoning below eq. (3.20).
In figure 9 we explore the effects of varying the vev v, the lifetime τ and the temperature ratio ξ n on the GW spectra. For comparison we show the expected power-law integrated sensitivity curves for LISA, ET and BBO [10]. As expected, the main effect of varying the vev (indicated by the different colours) is to change the nucleation temperature and hence the peak frequency. A change in the temperature ratio ξ n for fixed v amounts to altering the GW spectra as visualized by the dotted, dashed, and (dark) solid lines. As already observed in figure 8, increasing the temperature ratio by a factor of two from ξ n = 1 to ξ n = 2 is already sufficient to boost the GW spectrum significantly by one to two orders of magnitude, much to the benefit of experimental prospects.
The influence of the lifetime τ and hence the dilution effects can be observed in figure 9 by comparing the solid lines of different shading, with dilution effects increasing from darker to lighter shading for fixed v and ξ n = 10. As discussed above, increasing the lifetime of the dark Higgs boson corresponds to a larger entropy injection, diluting the signal towards both smaller signal strengths and frequencies, which eventually leads to a reduced experimental sensitivity in these scenarios. Note that we consider different lifetimes of the dark Higgs boson for the different vevs. The reason is that the quantity that sets the relevant scale for τ in the calculation of D is the Hubble parameter, which scales with T 2 DS . Hence, to obtain a comparable value for the signal strength and the dilution factor, smaller values of τ are needed for larger values of v, corresponding to larger nucleation temperatures. For the darkest (uppermost) curves, we have chosen τ such that the dilution factor is of order unity for all three cases.
Let us finally have a closer look at the dependence of the GW spectra on the vev. We observe that for v = 1 GeV (ξ n = 1) a peak signal strength of h 2 Ω GW ≈ 10 −13 is obtained, whereas the corresponding signals for v = 1 TeV and v = 1 PeV lie at significantly smaller values. The underlying reason for this is that larger values of v correspond to more relativistic 1 P e V τ = 1 ms, ξ n = 1 τ = 1 ms, ξ n = 2 τ = 1 ms, ξ n = 10 τ = 1 s, ξ n = 10 τ = 1 ks, ξ n = 10 τ = 1 ns, ξ n = 1 τ = 1 ns, ξ n = 2 τ = 1 ns, ξ n = 10 τ = 1 µs, ξ n = 10 τ = 1 ms, ξ n = 10 τ = 1 fs, ξ n = 1 τ = 1 fs, ξ n = 2 τ = 1 fs, ξ n = 10 τ = 1 ps, ξ n = 10 τ = 1 ns, ξ n = 10 Figure 9. An overview plot for the different possible GW spectra that can be provided by our model for a strong first-order phase transition (λ = 1.5 · 10 −3 , g = 0.5), compared to the expected power-law integrated sensitivity curves for LISA, ET and BBO [10] (see also Ref. [62,63]). The plot shows the resulting spectra of the phase transition of a dark Higgs acquiring its tree-level vev v = 1 GeV (blue), v = 1 TeV (purple), or v = 1 TeV (red). Dotted lines refer to the case when ξ n = 1, whereas dashed and solid lines indicate ξ n = 2 and ξ n = 10, respectively. The dependence of the spectrum on the dark Higgs lifetime τ is indicated by lighter colors. The main result is that increasing the temperature ratio ξ n leads to a strong enhancement of the signal strength when the dark Higgs decays sufficiently fast. The tree-level vev shifts the signal to different frequencies and has a mild impact on the signal strength for v ≤ 100 GeV.
degrees of freedom and hence a larger energy density in the SM at nucleation, leading to smaller α. However, eventually the degrees of freedom reach their maximal value in the SM at about v 100 GeV, beyond which α does not change considerably with the vev anymore. For our example, we have α = 0.51 for v = 1 GeV and α = 0.11 for v = 1 TeV, 1 PeV when considering an equally hot dark and visible sector (ξ n = 1). Given the scaling h 2 Ω GW ∝ α 2 /(1 + α) 2 [10,64] we are also able to deduce that the enhancement of the GW signal with increasing temperature ratio is stronger the weaker the initial spectrum, i. e. if the initial value of α is rather small. Therefore, the effect of increasing ξ n is larger for v = 1 TeV, 1 PeV than for v = 1 GeV. This consideration also explains the fact that all spectra converge towards a comparable peak signal strength with large α for ξ n = 10 despite different underlying values for the vev.
An important conclusion that can be drawn from figure 9 is that TeV-and PeV-scale vevs are most favorable for LISA and ET, respectively. We emphasize, however, that farfuture GW interferometers such as BBO [65], which would operate in an intermediate frequency range with very high sensitivity, would also be sensitive to these benchmark scenarios. For the following section, we will take a closer look at the two benchmark values v 1 = 2 TeV and v 2 = 10 PeV for λ = 1.5·10 −3 and g = 0.5. In these cases, the nucleation temperatures lie at T n DS,1 = 175 GeV and T n DS,2 = 851 TeV, respectively, while the dark photon mass is given by m A ,1 = 1 TeV and m A ,2 = 5 PeV. Consequently, we encounter the scenario mentioned before: the dark photon decouples on a very short timescale once the universe enters the new phase. As discussed above, we therefore set T cd DS = T n DS and identify the temperature ratio at nucleation ξ n with the one at chemical decoupling ξ cd .

Sensitivity for gravitational waves
In this section, we discuss the prospects for LISA and ET for the two benchmark points described above. While figure 9 allows to assess the observability of a GW signal in a qualitative manner, an analysis based on signal-to-noise ratios is in fact better suited to quantify the experimental sensitivities [10]. In this discussion, we aim to study the two main competing effects in particular: the enhancement of GWs at production through a large temperature ratio ξ n and the subsequent dilution of the spectrum determined primarily by the dark Higgs lifetime τ . In the following, we therefore analyse the signal-to-noise ratios in terms of these two quantities for fixed values of the other three parameters.
For the calculation of signal-to-noise ratios, we employ the auto-correlated optimal-filter measure as derived in Ref. [10]. The observed signal strength Ω GW (see eq. (2.22)) is normalized to the effective noise energy density spectrum Ω eff , integrated within the frequency band [f min , f max ] in which the detector is sensitive, and weighted with the duration of observation t obs . The noise spectrum not only encompasses instrumental noises but also noise from unresolved galactic binaries. We will refer to the signals as being detectable by LISA and the Einstein Telescope for signal-to-noise ratios exceeding the respective threshold SNR values of 10 and 5. A detailed overview of the calculation of signal-to-noise ratios as well as the used noise curves can be obtained from Ref. [10].
In figure 10 we visualize our main results for the LISA and ET benchmark points. In the left column, the dependence of the dilution factor D on ξ n and τ is shown, while the right-hand side focuses on the signal-to-noise ratios in the same parameter plane. For both LISA and ET we obtain qualitatively similar situations. For equally hot dark and visible sectors (ξ n = 1), the GWs from the first-order phase transition are not observable irrespective of what the dark Higgs lifetime is. Only when increasing the temperature ratio to ξ n 2, amounting to more energy stored in the GWs, we reach parameter regions for which both experiments become sensitive, thus validating our initial motivation that hot dark sectors greatly increase the observability of GWs. However, we also observe that for larger lifetimes, eventually, the signals get weaker again in spite of large ξ n . This is the case because dilution effects become increasingly important the larger the dark Higgs lifetime is, redshifting the GW spectrum and decreasing the experimental reach for such scenarios. The sensitivity loss for larger τ is also independent of ξ n for large enough temperature ratios since both α and D saturate eventually, leading to a nearly horizontal turnover of the LISA and ET sensitivity. This interplay between enhancing and diluting the GWs ultimately results in a rectangular shape of the relevant signal-to-noise ratio region in the right column.
We emphasize that in the most interesting parameter regions the lifetime of the dark Higgs boson is sufficiently small to not interfere with BBN. However, given our assumption that the dark Higgs boson decays non-relativistically, the lifetime cannot be arbitrarily small. This is indicated by the blue shaded region in the right column, which corresponds to θ nr ≥ 1  Figure 10. Dependence of the dilution factor (left) and the expected signal-to-noise ratios in future GW observatories (right) as a function of the temperature ratio at bubble nucleation and the dark Higgs lifetime. In the top row we consider a Higgs vev v = 2 TeV, corresponding to potentially observable signals in LISA, while in the bottom row we choose v = 10 PeV and consider ET. In both cases we find that large temperature ratios significantly enhance the predicted strength of the GW signal, provided the dark Higgs lifetime is short enough to avoid significant dilution.
(see section 3). Based on the left column of figure 10, however, we do not expect dilution effects to become relevant in this parameter region, which would correspond to the case that the universe never enters a period of early matter domination. Taking all relevant effects and conditions into account, we are therefore able to find sizeable parameter regions in ξ n and τ that predict sufficiently high signal-to-noise ratios to offer attractive prospects for future GW observatories. For LISA the signal-to-noise ratio can be as large as O(100), while the corresponding values for ET are somewhat smaller due to its steeper sensitivity loss for small frequencies. This finding demonstrates clearly that GW signatures, which may seem out of reach, can in fact be enhanced as soon as the dark sector is hot, enabling one to distinguish GWs emerging in the early universe from noise affecting the experimental measurement. Further progress in testing the GW scenarios discussed here can be expected from BBO, which will also be able to cover vevs within the TeV-PeV range and also be sensitive to scenarios with stronger dilution effects (see figure 9).
To conclude this discussion, we emphasize once more that the enhancement and dilution effects considered in this work apply to a wide range of GW signals from dark sector phase transitions. In particular, our findings can be directly applied to different values of λ and g, which would lead to smaller α and larger β/H, corresponding to overall weaker GW signals. In this case a larger vev would be required to achieve a similar peak frequency and to counterbalance the change in the nucleation temperature. At the same time, the effect from increasing ξ n would be even larger, given that the saturation of h 2 Ω GW ∝ α 2 /(1 + α) 2 would be delayed. Hence, for sufficiently large temperature ratios ξ n even such comparably weak GW signals may be rendered observable.

Conclusions
In this work we have considered the exciting prospect that future gravitational wave (GW) observatories will be able to measure the stochastic GW background from first-order phase transitions. Such first-order phase transitions arise frequently in extensions of the Standard Model (SM) that feature a dark sector, i. e. a collection of new states that interact with each other but only very weakly with the SM. A crucial property of such a dark sector is that it may have a temperature different from the temperature of the SM thermal bath. The hotter the dark sector, the larger its contribution to the total energy density and therefore the stronger the GW signals that can be produced. In this work we have focused on the case that the dark sector has a larger temperature than the visible sector, which can significantly boost GW signals that would otherwise be unobservable.
Such a set-up however faces a great challenge: In order to recover the standard cosmological evolution at low temperatures, the entropy stored in the dark sector needs to be transferred to the visible sector. To achieve this goal we have considered the case that the lightest dark sector particle (called the mediator) is unstable against decays into SM particles. These decays must be sufficiently slow that they do not bring the two sectors into thermal equilibrium at early times. On the other hand, if these decays happen too late, the universe enters a phase of early matter domination and the eventual transfer of entropy to the visible sector leads to a strong dilution of the GW signal.
We have explored the dependence of this dilution effect on the properties of the mediator (specifically its mass, lifetime, decoupling temperature and self-interactions) and the dark sector temperature. We have extended previous works on the topic by considering a period of cannibalism, during which the energy in the dark sector decreases through number-changing processes. Such cannibalism can significantly reduce the resulting dilution factors and thereby extend the range of mediator lifetimes, for which GW signals may be observable.
To apply our findings to a realistic scenario, we have considered a dark sector describing the spontaneous breaking of a new U (1) D gauge symmetry with a dark Higgs field and a dark photon as field content. In the parameter regions where a strong first-order phase transition is predicted, the lightest dark sector particle is the dark Higgs boson, while the dark photon obtains a mass larger than the dark sector temperature at bubble nucleation. If the dark sector temperature is equal to the SM temperature, the predicted GW signals are below the sensitivity of next-generation GW observatories such as LISA or the Einstein Telescope. However, temperature ratios of order 2 are sufficient to boost the GW signal above the expected level of noise. The subsequent dilution of the signal due to entropy injection remains small provided the dark Higgs bosons decay sufficiently quickly after the phase transition. We find that the interesting regions of parameter space span several orders of magnitude in the dark Higgs lifetime. Further regions of parameter space are expected to open up when extending our analysis to the case that the dark Higgs bosons decay before they become non-relativistic. The code used to obtain our results is publicly available as TransitionListener at https://github.com/tasicarl/TransitionListener. We emphasize that our calculation of the actual GW signals and the resulting signalto-noise ratios is rather simplified and could be improved in a number of ways, for example by considering the percolation temperature instead of the nucleation temperature or by including effects such as bubble filtering or the production of heavy particles during bubble expansion. Nevertheless, our central findings are independent of these approximations: The enhancements that we find for large dark sector temperatures and the dilution factors that we calculate are independent of the details of the phase transition and can equally be applied to more refined calculations. Our largely model-independent treatment facilitates the transfer of our results to different settings.
For the specific dark sector model that we consider it will be interesting to further explore the connection to dark matter. To do so, it will be essential to specify in detail the couplings of the dark sector to the SM. For example, if kinetic mixing between the dark photon and hypercharge gauge bosons is absent, the dark photon itself could be stable and a viable dark matter candidate. Alternatively, the dark sector could contain additional fermions that freeze out before the phase transition. Such a set-up might furthermore provide an explanation for the difference in temperature between the dark and the visible sector that we have assumed from the beginning. Observatories such as LISA and the Einstein Telescope may therefore have a unique possibility to combine both dark sector and GW physics, with the experimental sensitivity enhanced greatly when hot dark sectors are involved.
of any decoupled non-relativistic species, such as a frozen-out dark matter component, which therefore scales as ρ mat = ρ cd matā −3 . The Friedmann equation thus reads a (θ) = τā(θ) ρ med (θ) + ρ mat (θ) + ρ rad ( where all ofā, r, g SM eff,ρ , g SM eff,s and S SM depend implicitly on θ. The last term in eq. (A.6) describes the amount of entropy that has been injected into the SM bath since T SM = T cd SM , thus increasing the radiation energy density therein.
In order to derive a second differential equation quantifying the entropy injection, we point out the relation The time evolution of the scale factor, the SM entropy, the mediator energy density and the effective degrees of freedom in the SM bath can therefore be described by the following set of coupled differential equations: The equations used in the main text are then obtained by neglecting the contribution from non-relativistic matter (r cd mat = 0) and focusing on temperatures above the MeV-scale, for which γ = 1.

B Modifications of CosmoTransitions
To perform our analysis of possible phase transitions in dark sectors, we used a customized version of CosmoTransitions [60]. CosmoTransitions comes with the necessary tools to trace the global and local minima of a given effective potential of one or multiple scalar fields. Moreover, it allows to identify the possible phase transitions between these minima. CosmoTransitions is often used as a benchmark code in the literature [8,66], as it is sufficiently stable and fast. Apart from the identification of first-order phase transitions, also the calculation of bounce actions and bubble profiles is possible with CosmoTransitions.
We first updated the individual modules of the program to work with Python 3 and extended it by an accurate nucleation criterion for first-order phase transitions in dark sectors with a temperature different from the SM bath. Next, we added the code necessary to compute the important phase transition parameters α and β/H. This requires a model file, in which the effective potential V 1−loop eff (φ) and the mass spectrum of the dark sector are given, and an additional module to calculate the effective degrees of freedom. Furthermore, we added a module for the calculation of dilution factors D SM , in which the decay of the dark sector is modeled. Another module for the calculation of stochastic gravitational wave spectra and the signal-to-noise ratios has been added to interpret the observability of the generated signals. This set of modules is controlled by an interface, which itself is executed by a small scan file, which defines the region of parameter space that one wishes to analyze. In addition, there are a few parameters for adjusting the settings, such as the accuracy of scans and the grid for the scan.