Hunting WIMPs with LISA: correlating dark matter and gravitational wave signals

The thermal freeze-out mechanism in its classical form is tightly connected to physics beyond the Standard Model around the electroweak scale, which has been the target of enormous experimental efforts. In this work we study a dark matter model in which freeze-out is triggered by a strong first-order phase transition in a dark sector, and show that this phase transition must also happen close to the electroweak scale, i.e. in the temperature range relevant for gravitational wave searches with the LISA mission. Specifically, we consider the spontaneous breaking of a U(1)′ gauge symmetry through the vacuum expectation value of a scalar field, which generates the mass of a fermionic dark matter candidate that subsequently annihilates into dark Higgs and gauge bosons. In this set-up the peak frequency of the gravitational wave background is tightly correlated with the dark matter relic abundance, and imposing the observed value for the latter implies that the former must lie in the milli-Hertz range. A peculiar feature of our set-up is that the dark sector is not necessarily in thermal equilibrium with the Standard Model during the phase transition, and hence the temperatures of the two sectors evolve independently. Nevertheless, the requirement that the universe does not enter an extended period of matter domination after the phase transition, which would strongly dilute any gravitational wave signal, places a lower bound on the portal coupling that governs the entropy transfer between the two sectors. As a result, the predictions for the peak frequency of gravitational waves in the LISA band are robust, while the amplitude can change depending on the initial dark sector temperature.


Introduction
Dark matter (DM) is known to be the dominant form of matter in the universe, but it has so far evaded any attempt of detection in the laboratory or by other non-gravitational means [1].These null results have cast doubt on the so-called WIMP miracle, where DM is produced from the thermal bath of Standard Model (SM) particles in the early universe, and which for a long time has been used to motivate sizeable couplings between DM particles and the SM.Indeed, it has been shown that thermal freeze-out may happen entirely within an extended dark sector, such that the observed DM relic abundance Ω DM h 2 ≃ 0.12 [2] can be reproduced without the need for any sizeable couplings between the dark and the visible sector [3].These so-called secluded DM models pose a great challenge for laboratory searches due to their apparent lack of testable predictions.
At the same time, gravitational wave (GW) observatories have opened a completely new window into the universe, making it possible to observe objects and phenomena that affect visible matter only through gravity.The proposed LISA mission [4] will extend this window to the mHz frequency range, allowing in particular for the observation of a stochastic GW background that would be connected to a strong first-order phase transition (PT) close to the electroweak scale [5][6][7].LISA therefore raises new hopes to detect dark sectors that are otherwise unobservable.Over the past few years, first-order PTs in dark sectors have been studied in great detail [8][9][10][11][12], and various correlations between GW signals and the phenomenology of DM have been explored [13][14][15][16][17][18][19][20][21][22][23][24][25].The conclusion of these studies is that it is difficult to robustly predict the expected amplitude of the GW signal for a given DM model, because strong PTs often only happen in special regions of parameter space.In other words, it appears generally challenging to identify a strong correlation between the GW amplitude and the DM abundance.In this work, we instead focus on the peak frequency of the GW signal and show that it can be tightly correlated with the predicted DM relic abundance.Intriguingly, when imposing the observed value of Ω DM h 2 = 0.12 and focussing on GW signals strong enough to be potentially observable, we predict a GW peak frequency that falls right into the most sensitive range of LISA.
Before describing our analysis in detail, let us provide a rough sketch of the argument.We consider a dark sector comprised of a fermionic DM candidate χ charged under a new U (1) ′ gauge group that is spontaneously broken by the vacuum expectation value (vev) v ϕ of a new dark Higgs field.It is well known that strong PTs can occur in this model for a sufficiently large gauge coupling [26,27].All newly introduced particles are massless before symmetry breaking and acquire a mass proportional to v ϕ afterwards.The dark gauge boson A ′ (a.k.a.dark photon) and the dark Higgs boson ϕ are generally unstable against decays into SM particles, but χ is stable and may obtain a sizeable relic abundance through thermal freeze-out.If the spontaneous symmetry breaking occurs in a first-order PT, bubbles of the new phase will nucleate spontaneously, expand and collide.This process perturbs the dark plasma and leads to the emission of GWs, with a present-day peak frequency very roughly given by [7] f peak ≃ 10 mHz β/H 100 Here β/H denotes the speed of the PT and T p is the temperature of the SM heat bath at the time of percolation.For a not-too-strongly supercooled dark sector PT, which is what we consider here, one expects β/H ∼ 100 and T p ∼ v ϕ .The relic density from thermal freeze-out, on the other hand, can in leading-order approximation be written as [28] Ω DM ≃ 0.1 10 −8 GeV −2 ⟨σ ann v⟩ , with ⟨σ ann v⟩ the thermally averaged DM annihilation cross section.If the DM particles dominantly annihilate into the dark Higgs bosons ϕ, arising from the same dark Higgs field that generates the DM mass, it is parametrically of the form where y denotes the DM Yukawa coupling.At first sight, this coupling is arbitrary, and hence the freeze-out mechanism does not predict a specific dark sector mass scale.However, if we are interested in dark sectors that produce strong first-order PTs and large GW signals, the dark gauge coupling g and the dark Higgs quartic coupling λ must be sizeable, which implies that the dark Higgs boson mass m ϕ cannot be much smaller than v ϕ .At the same time, the observed DM relic abundance can only be obtained through dark sector freeze-out if the DM particle is not the lightest particle in the dark sector (or at least not much lighter than its annihilation products [29]).This, in turn, implies that y cannot be much smaller than unity, and hence v ϕ ∼ TeV once we require Ω DM h 2 ∼ 0.1.Combining this with the conclusion from eq. (1.1), we thus expect a peak frequency of f peak ∼ 10 mHz -which, as advocated, lies right within the LISA band.
A possible concern with the simplified reasoning above is that a large Yukawa coupling will affect the effective potential and may possibly prevent a first-order phase transition, or even destabilize the scalar potential [30].We have also neglected the impact of additional DM annihilation channels involving dark photons.In our full analysis, we explore the entire parameter space of the model, calculating in detail the effective potential, the thermodynamic quantities characterising the PT and the relic density from thermal freeze-out.We then identify viable combinations of the different dark sector couplings and show that the qualitative argument from above is confirmed by quantitative calculations.In order to further refine the analysis, we also perform parameter scans over all relevant model parameters -namely the three couplings g, λ and y and the dark Higgs vev v ϕ , and we identify parameter points for which the correct DM abundance is obtained.Interpreting the sampling distributions for the model parameters as prior probabilities thus enables us to define "typical" model predictions and quantify the probability (in the Bayesian sense) of a detectable signal.
A significant focus of our analysis is to extend the simple argument sketched above to situations where the couplings are so weak that the dark and visible sectors do not necessarily share a common temperature, which would be maintained through (inverse) decays of SM and dark Higgs bosons.Indeed, even if the two sectors have the same temperature initially, the first-order PT in the dark sector will change the temperature ratio, as the vacuum energy in the dark Higgs field is converted to rest mass and kinetic energy.This additional energy needs to be rapidly transferred to the SM in order to avoid a dilution of GW signals from late-time entropy injection [11,12].We calculate the dilution of the GW background and derive a lower bound on the portal coupling from the requirement that no significant dilution occurs.We show that the portal coupling required for this purpose is well below the sensitivity of laboratory experiments.Finally, we explore what happens if the initial temperature ratio of the two sectors differs from unity.In this case the amplitude of the GW signal will change [10,31,32] -but the peak frequency remains almost unaffected, such that the estimate from above remains robust even for portal couplings that are too small to quickly (re-)thermalize the sectors after the transition.This conclusion is only modified if the portal coupling is so weak that the energy density of the dark sector cannot be depleted and starts to dominate the energy density of the universe.
The remainder of this work is structured as follows.In section 2 we introduce the model under consideration and discuss the finite-temperature effective potential.We also briefly review the calculation of the temperature and strength of the PT, and we identify the interesting regions of parameter space.In section 3 we calculate the DM relic density under the assumption that the dark and SM sector remain in thermal equilibrium throughout their evolution, and explore the correlation between the relic density and the GW signal.We revisit this assumption in section 4, and discuss in detail the processes that thermalize the dark sector with itself and with the SM.In section 5 we finally calculate the effect of inefficient thermalisation on the GW signal.We consider the dilution due to entropy injection and show that for hot dark sectors a net enhancement of the GW amplitude can remain, while the peak frequency is essentially unaffected.We conclude in section 6 with a summary of our results and some remarks about their consequences.In two technical appendices, we provide details on the bubble wall velocity (appendix A) and on the Boltzmann equations for entropy transfer (appendix B).
2 Dark sector phase transition

Dark sector model
The model we study in this work is an extension of the models considered in refs.[11,33] and consists of a complex scalar Φ charged under the U (1) ′ gauge symmetry, the associated gauge boson A ′ µ , and two chiral fermion DM candidates, χ L and χ R .The Lagrangian describing the model is where A ′ µν is the field strength tensor of A ′ µ , λ and y are dimensionless couplings and µ is a bare mass parameter for Φ.The complex scalar and fermions are charged under the U (1) ′ group as The kinetic mixing of the dark photon with SM hypercharge and the portal coupling of the scalar field Φ with the SM Higgs field are assumed to be small enough that they satisfy experimental bounds (see, e.g., refs.[34,35] for bounds on dark photons and ref. [36] for a recent review of bounds on Higgs mixing) and do not play a role during the PT.We will return to these terms in section 4.2 when discussing the thermalisation of the dark sector with the SM bath.
The tree-level scalar potential of our model has a minimum at v ϕ = ± µ 2 /λ.One can hence expand the complex field as Φ = (v ϕ + ϕ + iφ)/ √ 2, where ϕ and φ are real scalar fields.In addition, the chiral fermions χ L and χ R can be written as a Dirac fermion χ.The Lagragian in eq.(2.1) can thus be re-written as where g is the gauge coupling associated with the U (1) ′ symmetry, and the bare masses of the various fields depend on the vacuum expectation value (vev) v ϕ as

The effective potential
The properties of the PT can be computed from the temperature-dependent effective potential V eff .This effective potential can be constructed out of the tree-level potential extracted from eq. (2.2), the 1-loop Coleman-Weinberg (CW) potential V CW and the finite temperature potential V T [37].The CW potential in our model, computed using the background field (ϕ b ) method [38], in the Landau gauge and the MS renormalisation scheme, has the form [39,40] V where the + (−) sign applies to bosons (fermions), g a are the degrees of freedom of particle a, m a (ϕ b ) are the background field-dependent masses, and k a = 3/2 for scalars and fermions and k a = 5/6 for gauge bosons.In this study we take the renormalisation scale Λ to be the vev at zero temperature, v ϕ (T = 0).The field-dependent masses in the CW potential are computed with renormalised parameters, i.e. using the expressions from eq. ( 2.3) with the replacement v ϕ → ϕ b and with the addition of relevant counterterms [37].We choose the renormalisation conditions such that the dark Higgs mass and vev are fixed to their tree-level values at zero temperature.
The finite-temperature part has the form with the thermal functions J b/f for bosons/fermions defined in ref. [37].The finite temperature potential in eq.(2.5) suffers from infrared divergences for bosonic modes when T ≫ m, leading to a breakdown of finite-temperature perturbation theory at high temperatures [40].The resummation of these modes is commonly done by "daisy resummation" [41].Here we follow the Espinosa method [42], which adds an additional contribution to the effective potential of the form where a now only runs over the scalar fields and the longitudinal component of the dark photon.The one-loop thermal masses Π a in our model are given by We only resum the Matsubara zero modes in the daisy potential, which are the ones causing the infrared divergences of eq.(2.5), such that the thermal mass of the DM fermion does not enter our calculations.The final form of the effective potential is then given by For illustration, we plot this potential in figure 1 for a choice of parameters g = 0.67, λ = 0.0035, and varying values of y.
In addition to encoding the properties of the PT, the effective potential also provides information about the stability of the true vacuum after the phase transition occurs.In fact, a new feature becoming important for non-zero Yukawa couplings is that for low values of λ and g the potential can become unbounded from below [30].To ensure vacuum stability we require that no deeper vacua are present at zero temperature.The requirement of a dark sector PT already implies that V eff (0) > V eff (v ϕ ).Hence, it is sufficient to check whether there exist vacua with lower potential energy for large field values, i.e. whether V eff (ϕ b ) < V eff (v ϕ ) for ϕ b ≫ v ϕ .In our analysis, we explicitly exclude such parameter points.
It is well known that the one-loop, daisy-resummed calculation of the effective potential can suffer from large theoretical uncertainties, foremost sourced by a large renomalization scale-dependence [43].A possibility to improve upon those uncertainties is to systematically resum higher orders of the thermal masses in the effective field theory framework of dimensional reduction [44].In order to validate our simpler approach, we therefore also implemented our model in DRalgo [45], which automates the task of dimensional reduction.We calculate the critical temperature in both our four-dimensional implementation and the reduced three-dimensional theory for the parameter space where we expect a first-order PT.In the regime where the effective field theory is valid (T ≫ m ϕ ) we find that the two results agree very well.We therefore conclude that we can take the computationally more economical approach of using the 1-loop, daisy-resummed effective potential stated in eq.(2.8).

Properties of the phase transition
The PT in our model occurs when the dark Higgs ϕ acquires a non-vanishing vev at the minimum of the effective potential V eff , thereby breaking the U (1) ′ symmetry.At the nucleation temperature T n the transition from the unbroken (false) vacuum to the broken (true) vacuum becomes energetically favourable and bubbles of the new phase start nucleating and expand rapidly to fill the entire Universe.The probability of transitioning from one phase to the other, parametrized by the bubble nucleation rate Γ(T ), can be computed from the effective potential in a semiclassical formalism by solving the bounce equation [46][47][48], which takes into account thermal tunneling through the potential barrier.We use TransitionListener [11] for this purpose, an extension of CosmoTransitions [49], which takes care of the computation of the bounce action, as well as phase tracing and the calculation of the thermodynamic properties of the phase transition.For the moment, we will assume equal temperatures for the SM bath and the dark sector, i.e. we adopt a temperature ratio of (2.9) The temperature that appears in the effective potential in eq.(2.8) and in the thermodynamic quantities discussed below can thus be identified with the temperature of SM photons.We discuss the general case with ξ ̸ = 1 in sections 4 and 5.
Gravitational waves from a PT are produced when the expanding bubbles or their sound shells collide (more in section 2.4).Therefore, the thermodynamic properties of the PT must be computed at a time when the rate of bubble collisions is maximised.This occurs at the time of percolation, when the Universe is permeated with a connected web of bubbles of the broken phase [7,50,51], which happens when approximately 70 % of the Universe is in the symmetric phase.Quantitatively, the fraction of the Universe remaining in the false vacuum after the PT is given by P (T ) = exp [−I(T )] with [40] where T is the common temperature of the dark sector (DS) and the SM bath, Γ(T ) is the bubble nucleation rate, H(T ) is the Hubble parameter, v w is the wall velocity and T c is the 'critical' temperature where the minimum of the effective potential with non-vanishing vev becomes a global minimum [52].Hence, the temperature at which percolation occurs can be computed by solving I(T p ) = 0.34 [53] with the approximation v w → 1.A discussion about the bubble wall velocity in our model can be found in appendix A. The dependence of the percolation temperature T p on the model parameters can be seen in figure 2. In the top panels T p is shown as a function of the quartic coupling λ and gauge coupling g for two values of the Yukawa coupling, y = 0 (left) and y = 0.5 (right).There is a strong correlation between the values of λ and g that produce a first-order PT, with lower values of T p in the top-left end of the allowed band, and higher values of T p in the bottomright.The disallowed areas correspond to parameter regions where the transition is not firstorder or it does not occur at all.These effects are better illustrated in the bottom panels, where the color scale indicates the ratio T p /T c in the same parameter plane.The amount of supercooling of the transition is largest when T p is much lower than T c and smallest when both temperatures almost coincide.For the points above the coloured contours, the potential barrier becomes so large that the bubble nucleation rate is too low for the transition to reach percolation; the region below instead indicates a smooth crossover transition in which no bubbles form since the potential does not develop a barrier between the phases.For non-zero values of the Yukawa coupling y, the enhanced thermal corrections in the effective potential cause a delay of the development of the true vacuum (cf.figure 1), thereby decreasing the value of T p .The vacuum also becomes deeper due to the Yukawa coupling, which increases the tunneling rate close to the supercooled region, and thus slightly larger values of g are within the allowed band.The grey shaded regions, finally, indicate parameter combinations where the potential is unstable (as discussed above in section 2.2).
The strength of the PT is quantified by the amount of released latent heat, described by the difference in the trace θ ≡ g µν T µν of the energy momentum tensors between the two phases [54], where ∆V eff is the difference in the effective potential between the vacua of the two phases.
A dimensionless quantity for the strength of the PT can be obtained by dividing ∆θ by the total energy density of the SM and dark sector plasma, ) where ρ SM,p (ρ DS,p ) is the energy density in the SM (DS) at percolation.The last missing piece for obtaining a GW spectrum is the speed of the PT, which can be determined through where S 3 is the O(3)-symmetric bounce action for thermal tunneling.
In figure 3 we show how the transition strength α and transition speed β/H depend on the model parameters, λ, g and y.The PT is relatively strong for most of the allowed region α ∈ (10 −2 , 10 2 ) and it is particularly strong close to the supercooled limit, where percolation is delayed (T p ≪ T c ).On the other hand, the speed of the PT β/H becomes smaller in the supercooling limit, reaching values of β/H ≈ 10 2 − 10 3 .

Gravitational wave spectrum
The spectrum of GWs in our scenario is produced dominantly through bulk fluid motion in the reheated plasma due to the large velocity-dependent friction from the emission of soft dark photons in the bubble wall, yielding a terminal bubble wall velocity [55][56][57].A discussion of this argument can be found in appendix A. As the case of runaway bubbles can hence be excluded, we neglect the contribution of bubble collisions to the GW signal.Since the onset of turbulence as a GW source is not yet understood well enough to make quantitative statements [7], and often requires complicated lattice simulations [58], we conservatively consider sound waves the only relevant source of GWs emitted during the PT.Therefore, the GW spectrum is exclusively determined by the set of parameters {α, β/H, T p }.We use semi-analytical approximations to compute the peak frequency and spectrum of GWs from sound shell collisions, based on simulations [7].Our computation of the GW spectrum includes some corrections taking into account that the transition happens in a dark sector and that the bubble wall dynamics are independent of the SM field content.The spectrum of gravitational waves is thus computed as [12] In this expression, the prefactor Rh 2 quantifies the redshift of the amplitude of the emitted GW signal and is given by [11,12] where Ω γ h 2 = 2.473 • 10 −5 is the present radiation energy density [2]; h SM,0 = 3.93 and g γ,0 = 2 are the entropy and energy degrees of freedom of the SM bath today [59].The total degrees of freedom at percolation g tot,p and h tot,p are fixed through g tot,p = g SM,p + g DS,p and h tot,p = h SM,p + h DS,p , respectively.The normalization of the signal in eq.(2.14) is given by Ω = 3 × 0.012 × 0.687 × (8π) 1/3 = 0.07, with the first two factors coming as normalization constants from ref. [60], the third one being the overall normalization of the spectrum S(f ) to unity [7] and the fourth one arising due to the conversion from mean bubble separations to β/H.The efficiency κ sw of converting sound waves into GWs is calculated from ∆θ/(4 ρ DS,p ), cf.eq.(2.11), using the high-v w approximation from ref. [61].The factor Y takes into account the lifetime of the sources of GWs and is given by [7,62] The spectral shape of the GWs is given by [5] with peak frequency [12] f peak = 8.9 mHz where we used that the peak frequency of the sound wave signal lies at 0.53 β at the time of its emission [60].
We will refer to a GW signal h 2 Ω GW (f ) as being observable by LISA, the Einstein Telescope [63,64] or pulsar timing arrays [65], if it achieves a sufficiently large signal-to-noise ratio.This means that the signal will be observable if its peak amplitude h 2 Ω peak GW ≡ h 2 Ω GW (f peak ) lies above the power-law-integrated (PLI) sensitivity curve at the peak frequency.The threshold signal-to-noise ratios and the PLI curves used in this work were derived in ref. [10].

Parameter ranges
The computation of properties and GW spectrum of the PT in section 2.3 and 2.4, respectively, depends only on the bare parameters of the Lagrangian in eq.(2.2), i.e. the gauge coupling g, the quartic coupling of the dark Higgs λ, the Yukawa coupling of the dark fermion y and the vacuum expectation value of the dark Higgs v ϕ .In our study we are specifically interested in exploring possible correlations between the DM relic density and a strong GW signal in our model.
A strong GW signal typically requires sizeable α and values of β/H that are not too large.From the discussion about figure 3, a strong first-order PT implies large but perturbative values of both g and λ.Too small values of these two couplings would imply very large values of β/H and correspondingly weak GW signals, and even cause issues of vacuum stability (for large y, cf.right panel of figure 3).This in turn induces an upper limit on the value of y as large values would cause an unstable vacuum for any perturbative value of g and λ.As will be seen below in section 3, successfully producing the right DM relic density requires m χ ≳ m ϕ , which implies a lower limit y > 2 √ λ.Lastly, the vev v ϕ is chosen in a range that produces GWs in the frequency range of near-future GW observatories, such as LISA.
Consequently, we randomly draw parameters from distributions that are logarithmically flat within the following ranges: 0.1 ≤ g ≤ 1, 10 −4 ≤ λ ≤ 10 −2 , 0.01 ≤ y ≤ 0.7 and 10 −3 GeV ≤ v ϕ ≤ 10 3 GeV.We then discard parameters that cause the vacuum to be unstable, that do not predict a first-order PT, or for which the PT is too supercooled and never percolates, thereby removing 82% of the points drawn.The remaining 18% of parameter points all feature a first-order PT with a corresponding GW signal.However, since the percolation temperature is very sensitive to small changes in the couplings, the PT is only strong enough to give an observable GW signal in certain small regions of parameter space.Indeed, only about 1% of paramter points from the original sample feature strong supercooling (T p /T c < 0.5).
We can quantify the fine-tuning required to obtain an observable GW signal by interpreting our parameter scan as a sample drawn from the prior distributions of the parameters.We then find that out of the parameter points that give a first order PT, only about 0.8% would be observable with LISA, whereas this number increases to 10% if we select parameter points that give a strongly supercooled PT.For the parameter ranges that we consider (in particular of v ϕ ) none would be observable with pulsar timing arrays or the Einstein Telescope.We note that these numbers do not correspond to rigorously calculated posterior probabilities, but rather rough estimates based on sampling densities.More precise estimates would require a different sampling strategy (see e.g.[66]), which is beyond the scope of this work.
We emphasize that this number is largely independent of the choice of priors as long as we select only parameter points that predict any kind of first-order PT.The probability to find parameter points that give a first-order PT does however depend sensitively on the choice of priors.If we were to extend the prior ranges for all parameters to lower couplings, the volume of parameter space without first-order PT would grow significantly.Choosing for example g > 0.01 (instead of 0.1), λ > 10 −5 (instead of 10 −4 ) and y > 10 −3 (instead of 0.01) would decrease the fraction of parameter points with a first-order PT from 18% to 6%.Out of these, 0.7% would be observable by LISA, which increases to 8.5% when considering only points with a strongly supercooled PT.As expected, our results are not very sensitive to different prior choices as we find that points that already have a first-order PT have a roughly equivalent probability of being visible at LISA regardless of the parameter ranges.In later sections, we will discuss how these numbers change when imposing additional constraints on the dark sector, such as the relic density requirement.Finally, when studying the effects of thermalisation in our model in section 4.2 it will be convenient to identify a benchmark scenario with the right properties for the PT and DM relic abundance.For reference the benchmark point is given in table 1.

Dark sector relic density
During the PT, the dark sector particles χ, ϕ and A ′ all obtain masses proportional to the dark Higgs vev v ϕ .In the parameter regions of interest for a strong first-order PT, we generally find g > √ 2λ and g > y/ √ 2 and hence the dark photon is usually the heaviest state in the dark sector, cf.eq. ( 2.3).Depending on the value of the Yukawa coupling y, the lightest dark sector particle will instead be either the DM fermion or the dark Higgs boson, as shown in figure 4. The dark sector equilibrises soon after the PT (see section 4 for a more detailed discussion).Typically, the heaviest particles will then first drop out of equilibrium as their number densities become strongly suppressed.The relic abundance of the dark fermions χ is thus determined through a freeze-out process [67] in the usual way.We assume that the dark photon is unstable, decaying for example through kinetic mixing, and therefore does not contribute to the DM relic density (unlike the case studied in ref. [22]).
In our model there are three possible DM annihilation processes that are relevant for setting the DM abundance: χχ → ϕϕ, χχ → ϕA ′ and χχ → A ′ A ′ .If the DM fermion is the lightest particle in the dark sector, annihilation into other dark sector states is kinematically forbidden for vanishing kinetic energy, such that the annihilation cross section becomes exponentially suppressed at low temperatures.In this so-called 'forbidden' regime [29], a relic abundance in accordance with observations requires that all mass scales must be correspondingly smaller, or the relevant couplings (much) larger.For the parameter values we are interested in here, it is therefore typically necessary for the DM particle to be heavier than the dark Higgs boson, which in turn requires a sizeable Yukawa coupling y.For even heavier DM, with 2m χ ≳ m ϕ + m A ′ , the annihilation channel χχ → ϕA ′ opens up.This process is a highly relevant contribution, once kinematically accessible, as it proceeds via an s-wave; the annihilation into a pair of dark Higgs bosons, χχ → ϕϕ, on the other hand, only proceeds via a p-wave.
To compute the DM relic density, we have calculated the amplitudes for all three processes, see appendix B.3, and implemented them in DarkSUSY [68], which calculates the thermal averages and solves the full Boltzmann equation [69].While DarkSUSY allows precision calculations of the relic density in a fully secluded dark sector with a varying temperature ratio ξ between the dark and the SM sector, cf.ref. [70], we will set ξ = 1 for the purpose of this section.We will revisit this assumption of thermal equilibrium between the two sectors in section 4.  We show the results from the parameter scan described in section 2.5 in figure 5.The three two-dimensional scatter plots show the correlation between the DM relic density Ω DM h 2 , the peak frequency f peak well as the peak amplitude Ω peak GW h 2 .One can immediately see that Ω DM and Ω peak GW are not tightly correlated (with a correlation coefficient of 0.20), while there exists a clear connection between the DM relic density and the peak frequency (with a correlation coefficient of 0.85).We can trace this correlation back to the fact that both quantities are determined by the dark Higgs vev v ϕ (indicated by the colour of each point).A smaller value of v ϕ implies a smaller DM mass and therefore a larger annihilation cross section, which in turn results in a smaller relic density.At the same time, a smaller v ϕ also implies a smaller percolation temperature, and hence a smaller peak frequency.The strength of the PT, on the other hand, depends on the details of the effective potential, and can vary over many orders of magnitude for any given value of v ϕ .We complement these scatter plots by showing distributions of the derived quantities, in the form of histograms based on our random scan described above.For example, one can infer that most samples drawn in  our setup correspond to a peak GW signal strength of Ω peak GW h 2 ≈ 10 −16 , i.e. a few orders of magnitude below the sensitivity of near-future GW observatories (indicated as grey shaded areas).Note also that the DM density caps at Ω DM h 2 ≈ 10, which would already correspond to an overclosed universe; even higher values are avoided by our prior choice, in particular the upper bound on v ϕ and the lower bound on y.
In figure 6 we show the result of sharpening the relic density requirement by requiring that 0.06 ≤ Ω DM h 2 ≤ 0.12.Demanding in this way that the fermionic DM candidate in our model constitutes the dominant form of DM, the predicted range of peak frequencies of the GW signal shrinks significantly -as expected from the discussion above.Interestingly, almost all viable parameter points now predict a peak frequency between 0.1 mHz and 100 mHz, largely overlapping with the frequency range to which LISA is sensitive.In fact, the peak frequencies for those parameter points that result in the strongest signal are the same as those where LISA is most sensitive.This striking correlation is a non-trivial feature of our model and constitutes one of our main results.Let us note that a few points remain that predict peak frequencies outside the LISA band.Much smaller values of f peak , in particular, correspond to parameter points in the 'forbidden' regime, m χ < m ϕ , where DM annihilations are exponentially suppressed at small temperatures.Smaller values of v ϕ (and hence smaller temperatures of the PT) can then still result in the correct DM relic abundance, but only at the cost of significant tuning between the various couplings (as reflected by the rareness of such parameter points, cf. the f peak histogram at the top of the plot).
In contrast to figure 5, where the colour coding of each point represents the dark Higgs vev, the points in figure 6 are coloured according to the percolation temperature of the PT, normalized to the vev v ϕ .Doing so allows us to confirm that the peak amplitude of the GW spectrum is determined primarily by the amount of supercooling.In other words, if the PT is delayed by a large potential barrier, the strength of the PT increases, yielding strong GW amplitudes (as expected from figures 2 and 3).As discussed in section 2.5, the predictions for the PT properties vary a lot with small changes of the model parameters, and thus only certain regions of the parameter space predict a strong first-order PT.For this reason, our model cannot in general guarantee a strong PT, and thus a GW signal that is visible with next-generation GW observatories.
We can make this statement more precise if we interpret the sampling distributions of the model parameters as prior probabilities (as we did in section 2.5), such that the density of points in figures 5 and 6 can be interpreted as probability distributions for the observables under consideration.As before, this makes it possible to quantify the amount of fine-tuning required to obtain a strong first-order PT, through the fraction of points with a first order PT that predict a signal observable with LISA.If we do not impose the relic density requirement (figure 5), only 0.8% of points with a first-order PT predict a GW signal visible at LISA, whereas this fraction increases to 3% once the relic density requirement is included (figure 6).If we restrict ourselves to parameter points with a strongly supercooled PT, the fraction of observable parameter points increases from 10% to 35%.Again, we have checked that this number is not very sensitive to our choice of parameter ranges.

Thermalisation of the two sectors
In this section we revisit the assumption that the temperature ratio of the dark and visible sectors is ξ = 1 throughout the PT.To do so, we first need to understand the evolution of the dark sector temperature during the PT, and convince ourselves that the dark sector quickly thermalizes with itself afterwards, such that the dark sector states remain in kinetic equilibrium with each other until after dark sector freeze-out (i.e. chemical decoupling).However, it is not necessarily the case that also the SM states are in kinetic equilibrium with the dark sector, such that their temperature may differ from the one of the dark sector both before and after the PT.We therefore discuss the various processes that allow for the exchange of energy and entropy between the dark and the SM sector, and the resulting Boltzmann equations.This enables us to identify the necessary portal couplings for efficient thermalisation.For the case of delayed thermalisation, after the end of the PT, we calculate the resulting dilution of the GW signal due to the injection of entropy into the SM thermal bath.
For the purpose of illustration, we will in this section consider a specific benchmark point that we selected from the random parameter scan discussed previously (see table 1).For ξ = 1, the parameters of this point lead to α = 0.258, β/H = 874, T n = 39.7 GeV, T p = 39.1 GeV, f peak = 3 mHz, Ω DM h 2 = 0.117, and Ω peak GW h 2 = 3•10 −13 .The rationale behind choosing this benchmark point is that i) the observed DM relic abundance is reproduced (for ξ = 1), and that ii) the PT is sufficiently strong in order to obtain an observable signal in LISA.We have explicitly checked that our choice is representative in the sense that other points fulfilling these two criteria lead to a very similar temperature evolution and resulting predictions.

The dark sector temperature
As the bubbles of the broken phase expand, more and more dark sector particles will pass through the bubble walls and enter the new phase.In the process, not only their rest masses but also their kinetic energies increase dramatically, by converting the vacuum energy of the dark Higgs field stored in the false vacuum.Here we neglect the small fraction of the energy density that is converted into GWs and assume that the bubble walls have already reached their terminal velocity, such that no energy is needed for their acceleration.As we have learned in the previous section, in particular, the energy density of GWs produced in the PT is bounded by Ω peak GW h 2 < 10 −10 and can therefore safely be ignored.We also neglect the effect of bubble filtering [16,71], i.e. we assume that all dark sector particles can enter the new phase.This is a good approximation for sufficiently fast bubble walls, see appendix A for details.
Since the different particle species in the dark sector were all relativistic before the PT, their number densities immediately after the PT will be comparable, even though their masses will now be very different.Indeed, for strongly supercooled PTs the dark photons (and possibly also the DM particles) will typically have a large mass compared to the temperature of the plasma, such that their equilibrium number density would be Boltzmann-suppressed.
In other words, right after the PT the dark sector finds itself far away from thermal equilibrium.Nevertheless, interactions between the different dark sector particles are rather strong, and hence the heavier particles are expected to annihilate rapidly into lighter ones, thereby restoring equilibrium.
As we will show below, the time required to reach equilibrium is negligible compared to the duration of the PT, such that we can to a very good approximation define a dark sector temperature of the broken phase T br DS immediately after the PT.This temperature is obtained from the temperature of the symmetric dark sector phase T sym DS using energy conservation: where T br DS denotes the temperature in the broken phase and where g DS (T ) takes into account the T -dependence stemming both from thermal (fielddependent) masses and the minimum of the effective potential.Eq. (4.1) can easily be solved numerically for T br DS for a given T sym DS .In practice, we find that slightly different temperatures T br DS of the broken phase are obtained when solving the equation taking T sym DS = T p or T sym DS = T n .This is because the energy density of the broken phase redshifts differently from the symmetric phase.We have therefore implemented a more detailed calculation, which tracks the temperature of the symmetric and broken phases from bubble nucleation to percolation and applies eq.(4.1) at each time step to the fraction of the universe entering the broken phase.Here the energy in the bubble walls, which for relativistic bubble wall velocities redshifts like radiation [72], is included in the energy in the symmetric phase.We find that this more careful treatment gives very similar results to simply applying eq. ( 4.1) at T sym DS = T p .We therefore use the latter prescription in the following when computing the temperature of the dark sector after the phase transition.
For our benchmark point, we find that the energy density of the dark sector before the PT is dominated by vacuum energy (see figure 7).In the broken phase, on the other hand, the vacuum energy is very small and quickly relaxes to a value very close to its zero temperature value as the temperature decreases further.This difference in vacuum energy leads to a substantial reheating of the dark sector, which, as a result, is hotter than the SM sector after the PT.For our benchmark point, we find that if the two sectors have equal temperature before the PT, in the broken phase the dark sector temperature will be larger by a factor of about 1.3.This reheating of the dark sector typically ensures that the dark Higgs bosons will be relativistic immediately after the PT.

Thermalisation within the dark sector
In the discussion above we have assumed that the dark sector can be characterised by a common temperature shortly after the PT.To justify this approach, we need an estimate of the time τ required to reach this equilibrium state and show that it is sufficiently small.For this purpose, we calculate the interaction rate for each dark sector state X in thermal equilibrium:  where the sum is over all dark sector states Y , σ XY denotes the total interaction cross section of X and Y , brackets denote thermal averaging (for simplicity calculated by assuming Boltzmann distributions) and n eq Y denotes the equilibrium number density of Y .A total of 20 different processes contribute to the thermalisation of the dark sector, the relative importance of which depends on the specific choice of parameters and the dark sector temperature.In the interest of brevity we refrain from stating the thermalisation rates explicitly.Broadly speaking, we find that Γ X is only a few orders of magnitude smaller than m X .For example, dark Higgs bosons can thermalise via self-scattering, i.e. ϕϕ ↔ ϕϕ, for which the scattering cross section is 9λ 2 /(8πs).For temperatures comparable to the dark Higgs boson mass, we have s ≈ 4m 2 ϕ and n ϕ ≈ ζ(3)m 3 ϕ /π 2 , such that Γ ϕ ∼ 10 −7 m ϕ for the benchmark point.Interactions of the dark Higgs bosons with dark fermions or dark photons benefit from the larger couplings y, g ≫ λ, but suffer from a Boltzmann suppression if T p < m χ , m A ′ .
A rough estimate of the thermalisation timescale is then obtained via Given the timescale τ we can estimate the out-of-equilibrium fraction of the universe F (t), which describes the relative volume of the universe that entered the broken phase so recently that it has not had enough time to reach thermal equilibrium.The related false vacuum fraction P (t), cf.eq.(2.10), describes the fraction of the universe which has not yet transitioned to the new phase, such that the true vacuum fraction is given by P t (t) = 1−P (t).Its time derivative Ṗt describes the rate with which the volume is transitioning to the true minimum of the potential for a given time t.We introduce the quantity  which can hence be interpreted as the volume fraction Ṗt ∆t that just transitioned to the broken phase within the small thermalization time scale ∆t = τ , cf. figure 7. The volume fraction F becomes small exactly when the thermalisation timescale τ is small compared to the transition timescale 1/β, as can be seen from the saddle point approximation of P (t) around the percolation time t p : ≈ βτ e β(t−tp) exp −0.34e β(t−tp) ≤ 0.37βτ .Here, the last term follows by inserting the time at which F (t) peaks, which is found to be t ≈ t p −1.08/β.Alternatively, one can interpret F as the volume fraction of a shell around the bubbles with the width of the mean free path of the particles that just entered the bubbles.In the left panel of figure 8 we show F as a function of T − T p for our benchmark point.As expected, we find that F takes its maximal value close to the percolation temperature T p .Nevertheless, even at the peak the value is tiny, implying that the fraction of the universe that is not in thermal equilibrium and therefore cannot be described by a temperature is completely negligible.In the right panel of figure 8 we show the maximum value F p at the percolation temperature for more points from our parameter scan, varying both g and λ.
Similarly we find that for all our cases the thermalisation is fast compared to the transition timescale.Even strongly supercooled points, where F p becomes close to unity, are expected to return to equilibrium before freeze out, since β/H ≫ 1.
So far, we have throughout assumed that chemical potentials in the dark sector can be neglected after thermalization has taken place.This is certainly a good approximation as long as at least the lightest state in the dark sector has a mass that is small compared to the dark sector temperature, which is typically the case shortly after the PT.However, as the universe continues to cool down, this assumption becomes increasingly critical.If we assume that the dark sector cannot transfer its entropy to the SM thermal bath, the subsequent evolution depends crucially on whether number-changing processes of the lightest state, such as 3ϕ → 2ϕ, are efficient enough to maintain chemical equilibrium in the dark sector.If this is the case, the dark sector temperature will decrease much more slowly than the SM temperature, and the universe will enter a period of 'cannibal' domination [73,74].If, on the other hand, number-changing processes are inefficient, the dark sector will develop large chemical potentials and the universe will eventually enter a period of matter domination.In both cases, the energy and entropy stored in the dark sector must later be transferred to the SM heat bath in order to recover radiation domination before neutrinos decouple at T ≈ 2 MeV, marking the onset of big bang nucleosynthesis [12].Neither of these scenarios is very desirable, as the GW signals from the PT will be strongly diluted in the process [11].
We are therefore more interested in the case where the dark and SM sector quickly equilibrate after the PT, such that their temperatures become equal, chemical potentials become negligible, and the universe evolves approximately as in radiation domination.In the following we discuss the processes that contribute to this process, and we derive the coupling strengths required for it to happen rapidly enough.

Thermalisation of the dark and visible sector
The conceptually simplest way for the dark and visible sector to exchange entropy and energy is via Higgs mixing [32,[75][76][77][78]. 1 Such a mixing arises from an additional term in the scalar potential: As long as the vevs of both Higgs bosons vanish, the dominant process connecting the two sectors is HH → ΦΦ.As soon as one of the two Higgs bosons acquires a vev, it can decay into the other one, i.e. h → ΦΦ (if the electroweak symmetry breaks first) or ϕ → HH (if the dark symmetry breaks first).If kinematically allowed, these decay processes typically dominate over the 2 → 2 process for non-relativistic particles.
After both electroweak symmetry and the dark gauge symmetry have been spontaneously broken, the Higgs mixing generates a non-diagonal mass term, which can be rotated away by introducing the mixing angle where we have assumed θ ≪ 1 both in order to satisfy experimental constraints on the properties of the SM-like Higgs boson and to ensure that thermal corrections from SM fields to the effective potential are negligible so that the dark sector PT can be treated separately from the EWPT.Both the masses and the vevs, and hence also the mixing angle, depend on the temperature.As a result of this mixing, the dark Higgs boson obtains couplings to SM fermions and gauge bosons proportional to θ.Of the greatest relevance for our discussion will be the decay of dark Higgs bosons into bottom quarks b, with a tree-level decay width given by To calculate the entropy transfer between the dark and visible sector, we define the heat transfer rate q ≡ ρ + 3H(ρ + P ) , (4.10) which is related to the change of entropy density via ṡ = −3Hs + q T .(4.11) At the same time, the first moment of the Boltzmann equation gives 12) The general expression for the first moment of the collision operator for decays (including relativistic corrections and quantum effects) was derived in refs.[79,80].It was shown there that the leading relativistic effects (namely a time dilation of the decay proportional to 1/γ and an increase of the injected energy by a factor γ) cancel, and it is therefore a good approximation to assume that the decaying particle is at rest.The integral of the collision operator thus gives q ≃ m ϕ ( ṅ + 3Hn) , A completely analogous equation holds for the evolution of the SM entropy density s SM .Since the Hubble rate H depends on the combined energy density of both sectors, both equations need to be solved simultaneously, together with the equation ȧ = Ha, which we will need below to calculate the evolution of the GW spectrum.
In practice, we also include decays into lighter quarks and leptons, which become relevant if decays into bottom quarks are kinematically forbidden.We further include the processes h → ϕϕ and ϕ → hh if they are kinematically allowed.We do not, however, include 2 → 2 processes of the form q q → gϕ or qϕ → qg, which may give a non-negligible contribution for light dark Higgs bosons [81,82].Additional details, and the relevant equations, can be found in appendix B. We note that analogous equations to the ones above can be derived for the case that only one Higgs boson has a vev and the case that both symmetries are unbroken.
As discussed above, it is not necessarily the case that the dark sector is in kinetic equilibrium with the SM at high temperatures.In the following, we will therefore take the temperature ratio of the two sectors at the electroweak phase transition (EWPT), ξ EWPT , as a free parameter and calculate the evolution of ξ during the subsequent cosmological stages • For λ hϕ = 10 −6 , the two sectors thermalise efficiently already before percolation.The initial value of ξ EWPT is therefore inconsequential for the subsequent evolution, and we obtain the same results for all cases.
• For λ hϕ = 3 × 10 −7 , the two sectors do exchange energy and entropy already in the unbroken phase, but do not fully thermalise before percolation, such that ξ p depends on ξ EWPT .After dark symmetry breaking, the two sectors thermalise rapidly, so that the subsequent evolution, and in particular the relic density calculation, do not depend on ξ EWPT .
• For λ hϕ = 10 −7 , the energy exchange before dark symmetry breaking is completely negligible.Even after dark symmetry breaking, it will take a while for the temperatures of the two sectors to approach each other.Nevertheless, the two sectors reach equilibrium before the dark Higgs bosons become non-relativistic.
• For λ hϕ = 3 × 10 −8 , the two sectors do not quickly thermalize after the PT, and the universe enters an early period of cannibal domination. 2n general the value of λ hϕ needed to ensure thermalisation depends somewhat on the dark Higgs vev, since for m ϕ ≪ m h the mixing angle scales as θ ∝ λ hϕ v ϕ .Moreover, for small dark Higgs boson masses, decays into bottom quarks are kinematically forbidden and thermalisation is less efficient.For the parameter points that reproduce the observed relic density we find that the assumption ξ = 1 made in section 3 is well justified for λ hϕ greater than 10 −6 -10 −5 .This value should be compared to the currently strongest bounds from direct detection experiments, which are only sensitive to λ hϕ ≳ 10 −3 [36].

Gravitational waves from hot dark sectors
While the analysis of the dark sector PT is simplest for λ hϕ > 10 −6 , it is phenomenologically interesting to also consider smaller values of λ hϕ , such that the temperature ratio of the two sectors before the PT may differ from unity.The reason is that ξ p > 1 leads to an enhancement of the GW signal, as a result of the larger total energy density in the dark sector compared to the SM thermal bath [10,11].In this case, however, we also need to consider what happens to the energy density of the dark sector after the PT.
If the transfer of energy from the dark to the visible sector after the PT is slow, the energy density of the universe will eventually be dominated by non-relativistic dark sector particles.This effect is already visible in figure 9, where for small values of λ hϕ = 3×10 −8 the temperature ratio ξ still differs from unity when the lightest dark sector particle becomes nonrelativistic.When the non-relativistic dark Higgs bosons eventually decay into SM particles, their entropy is transferred to the thermal bath; this modifies the expansion history of the universe and leads to a dilution of the GW signal.It is a priori unclear whether this dilution effect dominates over the enhancement with increasing dark sector temperature, or whether a net increase in the strength of GW signals remains.Moreover, the dilution effect also shifts the GW frequencies and might thereby spoil the correlation between peak frequency and relic density found in section 3.In the following we will investigate these effects in detail.
If the portal coupling is extremely small, in principle even the relic density calculation could be modified.If the dark Higgs bosons become non-relativistic after freeze-out, in particular, they may come to dominate the energy density of the universe at later times, and dilute not only the GW energy density but also the DM energy density through their decays (see e.g.ref. [83]).If, on the other hand, the dark Higgs bosons are non-relativistic already during freeze-out, inefficient thermalisation between the two sectors may additionally imply a non-trivial evolution of the dark sector temperature during freeze-out.While these are interesting scenarios in their own right, they are beyond the scope of this work.Instead, we will here focus exclusively on the case where the dark Higgs bosons decay sufficiently quickly for the standard freeze-out calculation (with temperature ratio ξ = 1) to be valid.

Dilution of gravitational waves
We described the GW spectrum in section 2.4 under the assumption of equal initial temperatures between the dark and SM sectors and the absence of any entropy injection into the SM bath.If the temperatures are unequal, the total degrees of freedom used in eqs.(2.15) and (2.18) are modified as [84] g tot,p = g SM,p + g DS,p ξ 4 p , where the temperature ratio at percolation ξ p is derived from the initial temperature ratio ξ by assuming entropy conservation in the dark sector, which is valid before the onset of thermalization of the two sectors (see section 4.2).Gravitational waves redshift like radiation, such that the spectrum Ω GW observed today can be calculated from the emitted spectrum Ω em GW as Here we assume that all modes redshift equally and thus ignore possible spectral features, which can for instance arise in the case of an early matter domination [85].This is reasonable as we are only interested in the GW spectra close to their peak.Assuming entropy conservation between the end of the phase transition and today, we have Together with eq. ( 5.3) this directly yields the previously used eq.(2.15).If entropy is however not conserved in the SM bath at early times, the ratio of scale factors a 0 /a p in eq. ( 5.4) needs to be corrected by a factor of D 1/3 SM , where D SM = S SM,0 /S SM,p [11,86].We hence find, more generally, that where we introduced the more convenient dilution factor D ≡ D SM h SM,p /h tot,p [11,86].Analogously, the peak frequency obtains an additional redshift for the case of a non-conserved SM entropy and hence reads ( We show in figure 10 the dilution factor D, as a function of the portal coupling λ hϕ .Here, we choose the same benchmark point as studied in section 4, and show the result for different values of the initial dark sector temperature ratio ξ EWPT , defined at T SM = 150 GeV.As expected, for sufficiently large λ hϕ there is no significant dilution3 , as entropy is conserved and only radiation degrees of freedom contribute to the energy content of the universe.However, with decreasing λ hϕ this is no longer the case and the dilution factor grows, becoming sizeable for λ hϕ ≪ 10 −7 .

Results
We are finally in the position to calculate the GW signal for ξ ̸ = 1 and λ hϕ < 10 −6 .We show the result in figure 11, for various initial values of the temperature ratio as well as portal couplings.As anticipated, a net enhancement of the GW signal is found for ξ EWPT > 1, provided that λ hϕ is sufficiently small for the sectors to not equilibrate before the PT (cf.figure 9).The enhancement saturates for ξ EWPT ≳ 2, see also the discussion in ref. [11], implying a dark sector energy density that initially dominates over that of the SM sector.For too small portal couplings λ hϕ ≲ 3 × 10 −8 , on the other hand, the effect of dilution becomes relevant and the GW signal starts to become suppressed.Crucially, changing ξ and λ hϕ does not significantly affect the peak frequency, such that the GW signal remains within the LISA sensitivity range.
We can now test the robustness of our results from section 3, where we assumed ξ = 1 and λ hϕ > 10 −6 , by allowing larger values of ξ and smaller values of λ hϕ .In figure 12 we show the same result as in figure 6, but now for ξ EWPT = 2 and λ hϕ = 10 −7 .For comparison we show the 1D distributions from figure 6 as red lines.In this plot we have removed points for which the dark sector temperature still differs significantly from the SM temperature when  the dark Higgs boson becomes non-relativistic, i.e. for which ξ nr > 1.1 at T DS,nr = m ϕ /3.The reason is that for such cases our final predictions depend on the details of chemical decoupling within the dark sector, which we do not study further in this work.While this requirement removes almost half of the parameter points, it does not introduce any significant bias, i.e. the distributions of f peak and Ω peak GW h 2 look very similar with and without this additional requirement.
The parameter combination ξ EWPT = 2 and λ hϕ = 10 −7 leads to a nearly maximal enhancement of the GW signal.As expected, we find that the peak position of the GW signal is not affected, such that the frequency range implied by the observed DM relic abundance remains within the LISA sensitivity window.The amplitude of the GW signal, on the other hand, is slightly enhanced, as can be seen from the comparison in the corresponding onedimensional histogram.
We can make this statement more precise by once again interpreting the density of points in the scatter plots as a probability distribution for the observables.Compared to figure 6, we find that the probability to obtain a signal observable with LISA increases from 3% to 8%.Limiting ourselves to parameter regions with strong supercooling, the fraction of Full sample 0.1% 0.5% First-order PT 0.8% 3% First-order PT + relic density 3% 8% Strong supercooling 10% 21% Strong supercooling + relic density 35% 69% Table 2: Fraction of parameter points that predict an observable GW signal for LISA after imposing various selection requirements on the sample of points drawn from the parameter ranges discussed in section 2.5.
observable events increases from 35% to 69%.We summarize our findings in table 2. We emphasize that this large increase is a result of fixing ξ and λ hϕ to particular values.If we instead vary ξ and λ hϕ as part of the scan, most parameter combinations will either give very similar results to the case ξ = 1 considered in section 3 or lead to an extended period where the dark sector energy density dominates.

Conclusions
In this work we explored correlations between the DM relic density and GW signals arising from a first-order PT that breaks a U (1) ′ gauge symmetry and gives rise to the mass of the fermionic DM particle.We demonstrated that, while the amplitude of the GW signal depends on the details of the effective potential and can vary over many orders of magnitude, the peak frequency is tightly constrained once we impose the observed value for the DM relic abundance.Intriguingly, the peak frequency is found to lie exactly in the milli-Hertz range, which will be explored by the LISA mission.
The dark sector considered in this work is characterised by four parameters: the dark gauge coupling g, the quartic coupling of the dark Higgs field λ, the dark Yukawa coupling y and the dark Higgs vev v ϕ .As a first step, we calculated the effective potential and the percolation temperature of the PT and identified the regions of parameter space that give a strong (large α) and not too fast (β/H ∼ 10 2 -10 4 ) PT, corresponding to potentially observable GW signals.We showed in particular that large GW signals require sizeable couplings g and λ and occur also for large values of y.The relic density of the dark sector is determined by annihilations of DM fermions into pairs of dark Higgs bosons.The requirement to match the observed DM relic density then requires that the DM fermion cannot be much lighter than the dark Higgs boson, with mass m ϕ ∼ v ϕ , which in turn implies a sizable Yukawa coupling and a DM mass that is comparable to the dark Higgs vev v ϕ .The dark Higgs vev, on the other hand, determines the percolation temperature and hence the peak frequency of the resulting GW signal.This connection leads to a tight correlation between relic density and GW peak frequency.Through comprehensive scans of the parameter space, we confirmed that this correlation is indeed highly generic in our model.
A rigorous statistical interpretation of our results is beyond the scope of this work, but some estimates of the significance of the samples can be performed.A rough measure of the fine-tuning required to have a visible signal at LISA can be obtained by assuming that the sampling distributions of the parameters act as prior probabilities, and that the sampling density of points hence indicates the posterior distributions of derived quantities.Indeed, the majority of points drawn from these distributions do not feature any first-order PT at all.Out of the points that feature a strongly-supercooled PT (T p /T c < 0.5), the probability of producing a visible signal at LISA in our model is around 10%.With the additional requirement that the observed DM relic abundance be reproduced, this probability increases to 35%, as a result of the strong correlation between the predicted relic density and the peak frequency of the GW signal (see figure 6).
We then studied two connected questions: How does the dark sector transfer its energy density to the Standard Model?And is it justified to assume the same temperature for both sectors?Indeed, the PT leads to an increase of the dark sector temperature, as vacuum energy is converted into rest mass and kinetic energy.Having confirmed that the dark sector itself thermalises immediately after the PT, we discussed in detail how the two sectors thermalise with each other for the specific case that the two Higgs bosons in the theory interact via the portal coupling λ hϕ .After both electroweak symmetry breaking and dark sector symmetry breaking, this interaction leads to mixing between the Higgs bosons, such that they can decay into each other as well as into fermions of both sectors.
We derived and solved the Boltzmann equations for the entropy transfer between the two sectors and showed that for λ hϕ > 10 −6 the assumption of equal temperature for both sectors is well justified.This portal coupling is small enough to be consistent with all laboratory constraints, in particular with direct detection experiments.For even smaller portal couplings, on the other hand, the temperatures of the two sectors may differ significantly, motivating us to consider the initial temperature ratio ξ EWPT as a free parameter.While small values of λ hϕ combined with large initial values of ξ EWPT may lead to an increase in the amplitude of the GW signal, the dark Higgs bosons will decay only slowly after the PT and may end up dominating the energy density of the universe.The resulting entropy injection would then lead to a substantial dilution of GW signals.We find that for λ hϕ ≈ 10 −7 a net enhancement remains, demonstrating that it is possible to have ξ > 1 during the dark sector PT while at the same time avoiding the dilution of the GW signal due to sufficiently rapid thermalisation afterwards.Rpeating our parameter scans for ξ = 2 and λ hϕ = 10 −7 , we find that the impact on the GW spectrum is only modest; while the typical amplitude is slightly enhanced, the peak frequency remains unchanged.In combination, these effects increase the fraction of points with a strongly supercooled PT that would be observable in LISA to 69%.Hence, our conclusion regarding the correlation between the DM relic abundance and the GW peak frequency applies also to dark sectors that thermalise only slowly with the Standard Model.
An interesting open question is how the relic density calculation would change for even smaller values of the portal coupling than what we consider.In this case, the energy density of the universe would be dominated by non-relativistic dark Higgs bosons, which may develop a chemical potential if number-changing processes are inefficient.The relic density calculation then requires solving a coupled set of Boltzmann equations with non-trivial evolution of the dark sector temperature.While the details of this calculation are beyond the scope of this work, the general expectation is that the DM relic abundance would be increased.This might open up the possibility to have a dark sector PT in the nano-Hertz frequency range and hence of interest in the context of recent results from pulsar timing arrays [87][88][89].Another promising avenue for further investigations opens up due to our findings regarding the bubble wall velocity presented in appendix A. There exists a potentially relevant part in the parameter space of our model in which the Bödeker-Moore criterion hints towards nonrelativistic bubble wall velocities.For such low wall speeds the effect of bubble filtering [16] may be non-negligible, and the calculation of the DM relic abundance is expected to be more involved than presented here.
Let us finally mention that even for tiny portal couplings there is a chance to actually detect the DM particles that we consider: unlike for annihilation into dark Higgs boson pairs, the mixed annihilation channel into one dark Higgs boson and one dark photon is not suppressed in the limit of small DM velocities.If kinematically allowed, it may thus lead to observable signals in indirect detection experiments [90] such as CTA [91].Such an observation would raise the possibility to explore in practice the correlations studied in this work and pin down the detailed structure of the dark sector.project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Sk lodowska-Curie Staff Exchange grant agreement No 101086085 -ASYMMETRY, and from the Deutsche Forschungsgemeinschaft (DFG) through the Emmy Noether Grant No. KA 4662/1-2 and grant 396021762 -TRR 257.

A Dark bubble walls
In addition to the strength α of the PT, its speed β/H, the percolation temperature T p and the temperature ratio ξ p , it is also necessary to know the speed of the bubble walls v w in order to determine the processes that dominate the GW signal from a dark sector PT.While the former parameters can be obtained from the effective potential V eff (ϕ, T ), the bubble wall velocity depends on plasma effects of the expanding bubble walls and therefore requires additional considerations.For bubbles expanding into the vacuum (i.e. if bubbles expand into a plasma that is not influenced by a change of the scalar vev), there is no source of friction, such that bubble walls can accelerate up to the point of their collision.For walls that interact with the surrounding plasma, on the other hand, several model-dependent sources of friction have been discussed [55,[92][93][94][95].If the friction increases with the bubble wall velocity, the acceleration of the bubble walls eventually stalls and a terminal velocity is reached.In this case, the bulk motion of the plasma dominates the GW spectrum [57].
In this paper we take the approach suggested in ref. [96], i.e. we show that v w is either expected to be non-relativistic in our model or that the bubbles are relativistic, but do not run away due to the emission of soft gauge bosons in the bubble walls.We conclude that the plasma motion is responsible for the dominant part of the GW signal and the contribution from bubble wall collisions is negligible.
To decide whether a bubble wall can accelerate up to relativistic velocities, we use the Bödeker-Moore criterion [93], which relates the velocity-independent leading-order (LO) bubble wall friction P LO to the amount of liberated vacuum energy density ∆V eff : Bödeker-Moore criterion: ∆V eff > P LO Relativistic bubble walls ∆V eff < P LO Non-relativistic bubble walls .(A.1) We emphasize that this criterion is insufficient to decide whether walls can run away (i.e.accelerate indefinitely), because of the next-to-leading-order (NLO) friction P NLO , which scales with powers of γ w = 1/ 1 − v 2 w [55, 56,97].Bubbles can run away only if ∆V eff > P LO + P NLO for all v w .Otherwise they will reach a relativistic terminal velocity given by the equilibrium of forces, ∆V eff = P LO + P NLO .
The LO friction due to particles acquiring a mass when crossing the bubble wall is given by [61,93] where ∆m 2 i is the positive change of the mass square of particle species i during the PT, g i is the corresponding number of degrees of freedom, and c i = 1 (1/2) for bosons (fermions).This expression assumes that the particle masses outside the bubble are below the dark sector temperature.For ultra-relativistic bubble walls, the production of heavy particles with mass up to γ w T n can also add to the LO friction [94].For the dark sector considered in this work, all particles are light before the PT and subsequently obtain a mass comparable to the scale of the PT, such that (A.2) is the only relevant contribution to the LO friction.
In the case of a broken U (1) ′ with a fermionic species, the LO friction reads [61]  where in the last step we have used that the dominant contribution comes from the heaviest state in the dark sector, which in the parameter regions of interest is the dark photon.The amount of released vacuum energy ∆V eff can be estimated from the zero-temperature potential, which gives Hence we find that the Bödeker-Moore criterion for relativistic bubble walls is satisfied if m ϕ /m A ′ > T p /v ϕ , which is the case for the parameter regions that give strongly supercooled PTs corresponding to observable GW signals (T p ≲ 0.1 v ϕ ).In these parameter regions, the bubble walls are therefore expected to be relativistic, v w → 1.This finding also implies that we can neglect the effect of bubble filtering, which is only relevant in the (deeply) nonrelativistic regime v w ≪ 1 [16,71].In the regions in which weaker GW signals are expected, the Bödeker-Moore criterion instead hints towards slower bubble walls, see figure 13.
The NLO friction created by the emission of soft dark photons into the broken phase quickly starts to grow with γ w .The bubble walls will therefore reach a terminal, asymptotic bubble wall velocity which is close to the speed of light.The precise value of γ terminal w is unnecessary for our purposes, as the existence of a terminal yet relativistic bubble wall velocity is sufficient to assume a dominant GW emission through bulk fluid motion.A more refined calculation of the respective energy budgets for the processes emitting gravitational radiation was performed in refs.[98][99][100].In ref. [98] it was shown that for sufficiently high terminal bubble wall velocities the fluid profiles are strongly peaked, such that the emitted GW spectral shapes are in fact indistinguishable from bubble collisions.We conclude that it is hence a reasonable approximation to work with a GW spectrum that is solely sourced through sound waves.

B Boltzmann equations for thermalisation and dark matter freeze-out
Here we discuss the processes we consider for the entropy transfer between the DS and the SM.

B.1 Thermal mixing angle
The Higgs mixing angle defined in eq.(4.8) depends on temperature through the masses and vevs of the two Higgs bosons.The temperature dependence of the dark Higgs boson can be directly obtained from the effective potential, whereas we follow ref.[79] to implement the temperature dependence of the SM Higgs boson.For large values of the dark Higgs vev, we sometimes encounter the situation that the SM Higgs and dark Higgs mass become similar or even cross, in which case the mixing angle apparently diverges.To regulate this unphysical divergence, we have to include the finite width Γ h of the Higgs resonance.As shown in ref. [82], including the width leads to an effective mixing angle given by

B.2 The Boltzmann equation for entropy transfer
In our analysis we specify the initial conditions, i.e. the temperature ratio ξ EWPT of the dark sector and SM bath at the EWPT, for which we take T EWPT = 150 GeV.Which processes contribute to the thermalisation of the two sectors depends on whether or not the U (1) ′ gauge symmetry is already broken at this point.We consider two different timelines: • The EW symmetry breaks before the dark sector PT.This is the case for the majority of our parameter space.Here, the thermalisation between the two sectors is initially determined by the processes hh ↔ ΦΦ and the decay of the SM Higgs h → ΦΦ.Additional processes can only contribute after the dark sector PT.
• The dark sector PT occurs before EW symmetry breaking.This can happen for parameter points with a large vev of the dark Higgs and not too strong supercooling.Once both symmetries are broken, thermalisation proceeds via hh ↔ ϕϕ, the decays of the SM Higgs into dark Higgs h → ϕϕ if it is kinematically allowed, and the decays of the dark Higgs into SM particles ϕ → SM SM through Higgs mixing.
In the following we give the relevant expressions contributing to the entropy transfer.

B.2.1 2 → 2 processes
In most regions of parameter space, the dark sector phase is unbroken immediately after the EWPT, such that there are no dark Higgs boson decays that can transfer entropy between the two sectors.Here the 2 → 2 process induced by the portal coupling ΦΦ ↔ hh become relevant.Since the particles' thermal masses are smaller than the temperature, we cannot take the usual non-relativistic approach.Instead we will follow the relativistic treatment of the calculation of the entropy transfer developed in refs.[101,102], which we briefly sketch here.The heat transfer for 2 → 2 processes can be expressed as q 2→2 = d 3 p (2π) .It is easiest to calculate this cross section in the center of mass frame.However, since the Bose-Einstein and Fermi-Dirac distributions are not Lorentz-invariant we have to apply a Lorentz boost Λ from the cosmic rest frame where u = (1, 0, 0, 0) T into the center of mass frame, which can be parameterised by the rapidity η and two angles θ and φ; for details see ref. [101].The phase space distribution becomes With this, we can rewrite the center of mass cross section as (B.4) The matrix element for the processes ΦΦ → hh is at tree level simply given by M = iλ hϕ .For this angle-independent transition amplitude, we can integrate over the solid angle and obtain This expression can now be efficiently evalutated numerically.An analogous expression is obtained for the process hh → ϕϕ.We note that in principle there are additional 2 → 2 processes that may contribute to thermalisation.The process ϕϕ → t t via an off-shell SM Higgs boson is strongly Boltzmannsuppressed below the EWPT [79] and does not give a relevant contribution in the temperature range that we consider.However, processes such as ϕ+q → g+q with a quark in the t-channel can give a relevant contribution if the decay ϕ → b b is kinematically forbidden [81,82].Since this is the case only in a very small fraction of the parameter space that we consider, we neglect these processes, thus giving a conservative estimate of the thermalisation rate.

B.2.2 Standard Model Higgs boson decays
After both symmetries are broken and for temperatures comparable to the SM Higgs boson mass, a second process of interest besides dark Higgs decays is the resonantly enhanced pairannihilation of dark Higgs bosons into predominantly bottom quarks: ϕϕ → h → b b.In thermal equilibrium, the rate of this process can be related to the inverse process, which is the decay h → ϕϕ with partial width given by Before the dark sector symmetry breaking, we also have the decay h → ΦΦ with the decay rate which can be treated in analogy to the case above.

B.2.3 Full Boltzmann equation
The full Boltzmann equation that we solve before the PT, in case it occurs after the EWPT then reads where we use the decay widths of the SM Higgs into other SM particle pairs from ref. [36].
The equations for the SM entropy follow analogously.

B.3 Annihilation cross sections
In the following we list the various DM annihilation cross sections in the non-relativistic limit, up to second order in the CMS velocity v (of each of the DM particles):

Figure 2 :
Figure 2: The percolation temperature T p (top) and the ratio of percolation temperature and critical temperature T p /T c (bottom) in the λ − g plane for Yukawa couplings of y = 0.0 (left) and y = 0.5 (right) and a vev of v ϕ = 1 TeV.The coloured band shows the parameter region where a first-order PT is possible.

y
.1, m /m s u p e r c o o l i n g c r o s s o v e r y = 0.5, m /m s u p e r c o o li n g c r o s s o v e r u n s t a b le p o t e n t ia l .1,m A /m s u p e r c o o l i n g c r o s s o v e = 0.5, m A /m s u p e r c o o li n g c r o s s o v e r u n s t a b le p o t e n t ia l

Figure 4 :
Figure 4: The upper (lower) panels show the ratio of the dark Higgs boson mass m ϕ (the dark photon mass m A ′ ) to the mass of the DM fermion m χ for y = 0.1 (left) and y = 0.5 (right), as a functions of the gauge coupling g and the self-interaction λ.Note that these ratios are independent of the dark Higgs vev.

Figure 5 :
Figure 5: Scatter plots and 1D distributions of the DM density Ω DM h 2 , the GW density Ω GW h 2 and the GW peak frequency f peak .For comparison, the dashed line shows the observed DM density, Ω DM h 2 = 0.12 [2]; grey shaded areas show the PLI sensitivities [10] of pulsar timing arrays, LISA and the Einstein Telescope, respectively.

Figure 6 :
Figure6: Scatter plots and 1D distributions as in figure5, but with the additional constraint 0.06 ≤ Ω DM h 2 ≤ 0.12.Here, the color scale does not encode the vev, but the ratio of percolation temperature T p to the vev v ϕ , thus indicating the amount of supercooling.

Figure 7 :
Figure 7: Contributions to the energy density around the PT for our benchmark point, as a function of the SM temperature and for a temperature ratio ξ = 1.The energy densities in the symmetric and broken phase of the dark sector have two contributions, namely the energy of the particles ('rad') and the potential energy of the scalar field ('vac').
e r c o o li n g c r o s s o v e r u n s t a b le p o t e n t ia l

Figure 8 :
Figure 8: Left: The black line shows the energy fraction of the dark sector that is out of equilibrium for our benchmark point, as a function of the SM photon temperature.For reference, the fraction of the dark sector in the broken phase is given by the orange line and the percolation temperature as the dashed vertical line.Right: Out-of-equilibrium fraction as a function of λ and g for y = 0.62 and v ϕ = 430 GeV.

14 )
(4.13)    and the evolution of the number density is given by the usual Boltzmann equation ṅϕ + 3Hn ϕ = −Γ ϕ→b b n eq ϕ Here n eq denotes the equilibrium number density for T DS = T SM , whereas the actual number density can be calculated from the dark sector temperature and the assumption of negligible chemical potential.Putting everything together, we obtain ṡDS = −3Hs DS − m ϕ T DS Γ ϕ→b b n eq ϕ

TDS = m / 3 TDS = m / 25 h = 3 10 7 TDS = m / 3 TDS 7 TDS = m / 3 TDS 8 TDS = m / 3 TDS = m / 25 Figure 9 :
Figure 9: Evolution of the dark sector temperature (solid) and the SM temperature (dashed) for different values of the initial temperature ratio ξ EWPT as a function of scale factor.The vertical lines indicate the scale factor at percolation, which depends on the temperature ratio.The horizontal lines indicate when the lightest state in the dark sector becomes nonrelativistic (T DS = m ϕ /3) and approximately when the DM particles freezes out (T DS = m χ /25).The different panels correspond to different values of the portal coupling λ hϕ .

Figure 10 :
Figure 10: Dilution factor D, cf.eq.(5.5), as a function of the portal coupling λ hϕ , and for different values of the initial temperature ratio ξ EWPT as indicated.

Figure 11 :
Figure 11: Gravitational wave signal for the same scenarios as considered in figure 9, i.e. for various values of the portal coupling λ hϕ (different panels) and the initial temperature ratio ξ EWPT (different colours).

Figure 12 :
Figure 12: Same as figure 6, but without the assumption of thermal equilibrium between the two sectors.Specifically, we consider ξ EWPT = 2 and λ hϕ = 10 −7 .Compared to the situation in figure 6 (indicated by the red lines in the 1D histograms), the GW amplitude is shifted to slightly larger values, while the peak position remains almost unaffected.

Figure 13 :
Figure13: The leading-order friction P LO over the difference ∆V eff in potential energy between the true and false vacuum phases as a function of the quartic coupling λ and the gauge coupling g for values of the Yukawa coupling y = 0 and y = 0.5.Values greater than one indicate that relativistic bubble wall velocities cannot be reached, cf.eq.(A.1).For values smaller than one, a relativistic terminal velocity is expected.

Table 1 :
Benchmark point used for discussing the thermalisation of visible and dark sector.