# Hunting WIMPs with LISA: Correlating dark matter and gravitational wave
signals

Torsten Bringmann  Tomás E. Gonzalo  Felix Kahlhoefer  Jonas Matuszak  and
Carlo Tasillo

###### Abstract

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)^{\prime}$ 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.

##  1 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
$\Omega_{\text{DM}}h^{2}\simeq 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]. 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
$\Omega_{\text{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
$\chi$ charged under a new $U(1)^{\prime}$ gauge group that is spontaneously
broken by the vacuum expectation value (vev) $v_{\phi}$ of a new dark Higgs
field. It is well known that strong PTs can occur in this model for a
sufficiently large gauge coupling [23, 24]. All newly introduced particles are
massless before symmetry breaking and acquire a mass proportional to
$v_{\phi}$ afterwards. The dark gauge boson $A^{\prime}$ (a.k.a. dark photon)
and the dark Higgs boson $\phi$ are generally unstable against decays into SM
particles, but $\chi$ 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]

$\displaystyle f_{\rm peak}\simeq
10\,\text{mHz}\left(\frac{\beta/H}{100}\right)\left(\frac{T_{\text{p}}}{1\,\text{TeV}}\right).$
(1.1)

Here $\beta/H$ denotes the speed of the PT and $T_{\text{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 $\beta/H\sim 100$ and $T_{\text{p}}\sim v_{\phi}$.

The relic density from thermal freeze-out, on the other hand, can in leading-
order approximation be written as [25]

$\displaystyle\Omega_{\text{DM}}\simeq
0.1\,\frac{10^{-8}\,\text{GeV}^{-2}}{\langle\sigma_{\text{ann}}v\rangle}\,,$
(1.2)

with ${\langle\sigma_{\text{ann}}v\rangle}$ the thermally averaged DM
annihilation cross section. If the DM particles dominantly annihilate into the
dark Higgs bosons $\phi$, arising from the same dark Higgs field that
generates the DM mass, it is parametrically of the form

$\langle\sigma_{\text{ann}}v\rangle\sim\frac{y^{4}}{m_{\text{DM}}^{2}}\sim\frac{y^{2}}{v_{\phi}^{2}}\,,$
(1.3)

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 $\lambda$ must be sizeable, which implies that
the dark Higgs boson mass $m_{\phi}$ cannot be much smaller than $v_{\phi}$.
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
[26]). This, in turn, implies that $y$ cannot be much smaller than unity, and
hence $v_{\phi}\sim\mathrm{TeV}$ once we require $\Omega_{\text{DM}}h^{2}\sim
0.1$. Combining this with the conclusion from eq. (1.1), we thus expect a peak
frequency of $f_{\rm peak}\sim 10\,\text{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 [27]. 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$, $\lambda$ and $y$ and the dark Higgs vev $v_{\phi}$, 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, 28] – 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

###  2.1 Dark sector model

The model we study in this work is an extension of the models considered in
refs. [29, 11] and consists of a complex scalar $\Phi$ charged under the
$U(1)^{\prime}$ gauge symmetry, the associated gauge boson $A^{\prime}_{\mu}$,
and two chiral fermion DM candidates, $\chi_{L}$ and $\chi_{R}$. The
Lagrangian describing the model is

$\displaystyle\mathcal{L}=$
$\displaystyle|D_{\mu}\Phi|^{2}-\frac{1}{4}A^{\prime}_{\mu\nu}A^{{}^{\prime}\mu\nu}+\mu^{2}\Phi^{*}\Phi-\lambda(\Phi^{*}\Phi)^{2}$
$\displaystyle+\chi_{L}^{\dagger}i\not{D}\chi_{L}+\chi_{R}^{\dagger}i\not{D}\chi_{R}-y\Phi\chi_{L}^{\dagger}\chi_{R}-y\Phi^{*}\chi_{R}^{\dagger}\chi_{L}\,,$
(2.1)

where $A^{\prime}_{\mu\nu}$ is the field strength tensor of
$A^{\prime}_{\mu}$, $\lambda$ and $y$ are dimensionless couplings and $\mu$ is
a bare mass parameter for $\Phi$. The complex scalar and fermions are charged
under the $U(1)^{\prime}$ group as $Q_{\Phi}=+1$, $Q_{\chi_{L}}=+1/2$ and
$Q_{\chi_{R}}=-1/2$. The kinetic mixing of the dark photon with SM hypercharge
and the portal coupling of the scalar field $\Phi$ with the SM Higgs field are
assumed to be small enough that they satisfy experimental bounds (see, e.g.,
refs. [30, 31] for bounds on dark photons and ref. [32] 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_{\phi}=\pm\sqrt{\mu^{2}/\lambda}$. One can hence expand the complex field
as $\Phi=(v_{\phi}+\phi+i\varphi)/\sqrt{2}$, where $\phi$ and $\varphi$ are
real scalar fields. In addition, the chiral fermions $\chi_{L}$ and $\chi_{R}$
can be written as a Dirac fermion $\chi$. The Lagragian in eq. (2.1) can thus
be re-written as

$\displaystyle\mathcal{L}=$
$\displaystyle\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi+\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi-\frac{1}{4}A^{\prime}_{\mu\nu}A^{{}^{\prime}\mu\nu}-\frac{1}{2}m_{\phi}^{2}\phi^{2}+\frac{1}{2}m_{A^{\prime}}^{2}A^{{}^{\prime}2}_{\mu}$
$\displaystyle-
gA^{\prime}_{\mu}[\varphi\partial^{\mu}\phi-\phi\partial^{\mu}\varphi-
v_{\phi}\partial^{\mu}\varphi]+\frac{g^{2}}{2}\phi^{2}A^{{}^{\prime}2}_{\mu}+\frac{g^{2}}{2}\varphi^{2}A^{{}^{\prime}2}_{\mu}+g^{2}v_{\phi}\phi
A^{{}^{\prime}2}_{\mu}$ $\displaystyle-\lambda v_{\phi}\phi^{3}-\lambda
v_{\phi}\varphi^{2}\phi-\frac{\lambda}{4}\phi^{2}\varphi^{2}-\frac{\lambda}{4}\phi^{4}-\frac{\lambda}{4}\varphi^{4}$
$\displaystyle+i\bar{\chi}\not{\partial}\chi-
m_{\chi}\bar{\chi}\chi+\frac{g}{2}\bar{\chi}\not{A}^{\prime}\gamma^{5}\chi-\frac{y}{\sqrt{2}}\phi\bar{\chi}\chi+i\frac{y}{\sqrt{2}}\varphi\bar{\chi}\gamma^{5}\chi\,,$
(2.2)

where $g$ is the gauge coupling associated with the $U(1)^{\prime}$ symmetry,
and the bare masses of the various fields depend on the vacuum expectation
value (vev) $v_{\phi}$ as

$m_{\phi}^{2}=-\mu^{2}+3\lambda v_{\phi}^{2}=2\lambda v_{\phi}^{2},\quad
m_{\varphi}^{2}=0,\quad m_{A^{\prime}}^{2}=g^{2}v_{\phi}^{2},\quad
m_{\chi}^{2}=\frac{y^{2}}{2}v_{\phi}^{2}.$ (2.3)

###  2.2 The effective potential

The properties of the PT can be computed from the temperature-dependent
effective potential $V_{\rm 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_{\rm CW}$ and the finite temperature potential
$V_{T}$ [33]. The CW potential in our model, computed using the background
field ($\phi_{b}$) method [34], in the Landau gauge and the
$\overline{\text{MS}}$ renormalisation scheme, has the form [35, 36]

$V_{\rm CW}(\phi_{b})=\sum_{a=\phi,\varphi,A^{\prime},\chi}\pm
g_{a}\frac{m_{a}^{4}(\phi_{b})}{64\pi^{2}}\left(\log\frac{m_{a}^{2}(\phi_{b})}{\Lambda^{2}}-k_{a}\right),$
(2.4)

where the $+$ ($-$) sign applies to bosons (fermions), $g_{a}$ are the degrees
of freedom of particle $a$, $m_{a}(\phi_{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 $\Lambda$ to be
the vev at zero temperature, $v_{\phi}(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_{\phi}\to\phi_{b}$ and with
the addition of relevant counterterms [33]. 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

$V_{T}(\phi_{b})=\frac{T^{4}}{2\pi^{2}}\sum_{a}\pm
g_{a}J_{\text{b}/\text{f}}\left(\frac{m_{a}^{2}(\phi_{b})}{T^{2}}\right),$
(2.5)

with the thermal functions $J_{\text{b}/\text{f}}$ for bosons/fermions defined
in ref. [33]. The finite temperature potential in eq. (2.5) suffers from
infrared divergences for bosonic modes when $T\gg m$, leading to a breakdown
of finite-temperature perturbation theory at high temperatures [36]. The
resummation of these modes is commonly done by “daisy resummation” [37]. Here
we follow the Espinosa method [38], which adds an additional contribution to
the effective potential of the form

$V_{\rm
daisy}(\phi_{b})=-\frac{T}{12\pi}\sum_{a=\phi,\varphi,A_{\text{L}}^{\prime}}g_{a}\left[(m_{a}^{2}(\phi_{b})+\Pi_{a}(T))^{3/2}-(m_{a}^{2}(\phi_{b}))^{3/2}\right],$
(2.6)

where $a$ now only runs over the scalar fields and the longitudinal component
of the dark photon. The one-loop thermal masses $\Pi_{a}$ in our model are
given by

$\displaystyle\Pi_{\phi}$
$\displaystyle=\left(\frac{\lambda}{3}+\frac{y^{2}}{12}+\frac{g^{2}}{4}\right)T^{2}\,,$
$\displaystyle\Pi_{\varphi}$
$\displaystyle=\left(\frac{\lambda}{3}+\frac{y^{2}}{12}+\frac{g^{2}}{4}\right)T^{2}\,,$
$\displaystyle\Pi_{A^{\prime}}$
$\displaystyle=\left(\frac{1}{12}+\frac{2}{3}\right)g^{2}T^{2}\,.$ (2.7)

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

$V_{\rm eff}(\phi_{b},T)=V_{0}(\phi_{b})+V_{\rm
CW}(\phi_{b})+V_{T}(\phi_{b},T)+V_{\rm daisy}(\phi_{b},T)\,.$ (2.8)

For illustration, we plot this potential in figure 1 for a choice of
parameters $g=0.67,\lambda=0.0035$, and varying values of $y$.

Figure 1: Total effective potential (solid) and finite temperature part (dot-
dashed) for $g=0.67,\lambda=0.0035,v_{\phi}=1000\,$GeV and $T=240\,$GeV, for
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 $\lambda$ and $g$ the potential can
become unbounded from below [27]. 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_{\text{eff}}(0)>V_{\text{eff}}(v_{\phi})$.
Hence, it is sufficient to check whether there exist vacua with lower
potential energy for large field values, i.e. whether $V_{\rm
eff}(\phi_{b})<V_{\rm eff}(v_{\phi})$ for $\phi_{b}\gg v_{\phi}$. 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 [39]. 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 [40]. In order to validate our simpler approach, we therefore also
implemented our model in DRalgo [41], 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\gg m_{\phi}$) 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).

###  2.3 Properties of the phase transition

The PT in our model occurs when the dark Higgs $\phi$ acquires a non-vanishing
vev at the minimum of the effective potential $V_{\rm eff}$, thereby breaking
the $U(1)^{\prime}$ symmetry. At the nucleation temperature $T_{\text{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 $\Gamma(T)$, can be computed from the effective potential in a
semiclassical formalism by solving the bounce equation [42, 43, 44], which
takes into account thermal tunneling through the potential barrier. We use
TransitionListener [11] for this purpose, an extension of CosmoTransitions
[45], 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

$\xi\equiv\frac{T_{\text{DS}}}{T}=1\,.$ (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 $\xi\neq 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 [46,
47, 7], 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\left[-I(T)\right]$ with [36]

$\displaystyle
I(T)=\frac{4\pi}{3}v_{w}^{3}\int^{T_{\text{c}}}_{T}\mathop{}\\!\mathrm{d}T^{\prime}\,\frac{\Gamma(T^{\prime})}{T^{\prime\,4}\,H(T^{\prime})}\left(\int_{T}^{T^{\prime}}\frac{\mathop{}\\!\mathrm{d}T^{\prime\prime}}{H(T^{\prime\prime})}\right)^{3},$
(2.10)

where $T$ is the common temperature of the dark sector (DS) and the SM bath,
$\Gamma(T)$ is the bubble nucleation rate, $H(T)$ is the Hubble parameter,
$v_{w}$ is the wall velocity and $T_{\text{c}}$ is the ‘critical’ temperature
where the minimum of the effective potential with non-vanishing vev becomes a
global minimum [48]. Hence, the temperature at which percolation occurs can be
computed by solving $I(T_{\text{p}})=0.34$ [49] with the approximation
$v_{w}\rightarrow 1$. A discussion about the bubble wall velocity in our model
can be found in appendix A.

Figure 2: The percolation temperature $T_{\text{p}}$ (top) and the ratio of
percolation temperature and critical temperature $T_{\text{p}}/T_{\text{c}}$
(bottom) in the $\lambda-g$ plane for Yukawa couplings of $y=0.0$ (left) and
$y=0.5$ (right) and a vev of $v_{\phi}=1$ TeV. The coloured band shows the
parameter region where a first-order PT is possible.

The dependence of the percolation temperature $T_{\text{p}}$ on the model
parameters can be seen in figure 2. In the top panels $T_{\text{p}}$ is shown
as a function of the quartic coupling $\lambda$ 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 $\lambda$ and $g$ that produce a
first-order PT, with lower values of $T_{\text{p}}$ in the top-left end of the
allowed band, and higher values of $T_{\text{p}}$ in the bottom-right. The
disallowed areas correspond to parameter regions where the transition is not
first-order or it does not occur at all. These effects are better illustrated
in the bottom panels, where the color scale indicates the ratio
$T_{\text{p}}/T_{\text{c}}$ in the same parameter plane. The amount of
supercooling of the transition is largest when $T_{\text{p}}$ is much lower
than $T_{\text{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_{\text{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 $\theta\equiv g_{\mu\nu}T^{\mu\nu}$
of the energy momentum tensors between the two phases [50],

$\displaystyle\Delta\theta=\left[-4\Delta V_{\text{eff}}+T\frac{\partial\Delta
V_{\text{eff}}}{\partial T}\right]_{T_{\text{p}}}>0\,,$ (2.11)

where $\Delta V_{\text{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 $\Delta\theta$ by the total energy
density of the SM and dark sector plasma,

$\displaystyle\alpha\equiv\frac{\Delta\theta/4}{\rho_{\text{SM},\text{p}}+\rho_{\text{DS},\text{p}}}\,,$
(2.12)

where $\rho_{\text{SM},\text{p}}$ ($\rho_{\text{DS},\text{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

$\displaystyle\beta/H=T_{\text{p}}\left.\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}T}\frac{S_{3}(T)}{T}\right|_{T=T_{\text{p}}}\,,$
(2.13)

where $S_{3}$ is the $O(3)$-symmetric bounce action for thermal tunneling.

Figure 3: The transition strength $\alpha$ (top) and speed $\beta/H$ (bottom)
of the PT, in the $\lambda-g$ plane for $y=0$ (left) and $y=0.5$ (right), with
$v=1$ TeV and $\xi=1$.

In figure 3 we show how the transition strength $\alpha$ and transition speed
$\beta/H$ depend on the model parameters, $\lambda$, $g$ and $y$. The PT is
relatively strong for most of the allowed region $\alpha\in(10^{-2},10^{2})$
and it is particularly strong close to the supercooled limit, where
percolation is delayed ($T_{\text{p}}\ll T_{\text{c}}$). On the other hand,
the speed of the PT $\beta/H$ becomes smaller in the supercooling limit,
reaching values of $\beta/H\approx 10^{2}-10^{3}$.

###  2.4 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 [51, 52, 53]. 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 [54], 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
$\left\\{\alpha,\beta/H,T_{\text{p}}\right\\}$. 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]

$\displaystyle h^{2}\Omega_{\text{GW}}(f)$
$\displaystyle=\mathcal{R}h^{2}\,\tilde{\Omega}\left(\frac{\kappa_{\text{sw}}\,\alpha}{\alpha+1}\right)^{2}\left(\frac{\beta}{H}\right)^{-1}\mathcal{Y}\,S(f)\,.$
(2.14)

In this expression, the prefactor $\mathcal{R}h^{2}$ quantifies the redshift
of the amplitude of the emitted GW signal and is given by [11, 12]

$\displaystyle\mathcal{R}h^{2}$
$\displaystyle=\Omega_{\gamma}h^{2}\left(\frac{h_{\text{SM},0}}{h_{\text{tot},\text{p}}}\right)^{4/3}\left(\frac{g_{\text{tot},\text{p}}}{g_{\gamma,0}}\right)=1.653\cdot
10^{-5}\left(\frac{100}{h_{\text{tot},\text{p}}}\right)^{4/3}\left(\frac{g_{\text{tot},\text{p}}}{100}\right)\,,$
(2.15)

where $\Omega_{\gamma}h^{2}=2.473\cdot 10^{-5}$ is the present radiation
energy density [2]; $h_{\text{SM},0}=3.93$ and $g_{\gamma,0}=2$ are the
entropy and energy degrees of freedom of the SM bath today [55]. The total
degrees of freedom at percolation $g_{\text{tot},\text{p}}$ and
$h_{\text{tot},\text{p}}$ are fixed through
$g_{\text{tot},\text{p}}=g_{\text{SM},\text{p}}+g_{\text{DS},\text{p}}$ and
$h_{\text{tot},\text{p}}=h_{\text{SM},\text{p}}+h_{\text{DS},\text{p}}$,
respectively. The normalization of the signal in eq. (2.14) is given by
$\tilde{\Omega}=3\times 0.012\times 0.687\times(8\pi)^{1/3}=0.07$, with the
first two factors coming as normalization constants from ref. [56], 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
$\beta/H$. The efficiency $\kappa_{\text{sw}}$ of converting sound waves into
GWs is calculated from $\Delta\theta/(4\,\rho_{\text{DS},\text{p}})$, cf. eq.
(2.11), using the high-$v_{w}$ approximation from ref. [57]. The factor
$\mathcal{Y}$ takes into account the lifetime of the sources of GWs and is
given by [58, 7]

$\displaystyle\mathcal{Y}=\min\left[1,\tau_{\text{sh}}H\right]\simeq\min\left[1,\frac{3.38}{\beta/H}\sqrt{\frac{1+\alpha}{\kappa_{\text{sw}}\,\alpha}}\right].$
(2.16)

The spectral shape of the GWs is given by [5]

$\displaystyle
S(f)=\left(\frac{f}{f_{\text{peak}}}\right)^{3}\left(\frac{7}{4+3\,(f/f_{\text{peak}})^{2}}\right)^{7/2},$
(2.17)

with peak frequency [12]

$f_{\text{peak}}=8.9\,\text{mHz}\left(\frac{T_{\text{p}}}{100\text{
GeV}}\right)\left(\frac{\beta/H}{1000}\right)\left(\frac{g_{\text{tot},\text{p}}}{100}\right)^{1/2}\left(\frac{100}{h_{\text{tot},\text{p}}}\right)^{1/3}\,,$
(2.18)

where we used that the peak frequency of the sound wave signal lies at
$0.53\,\beta$ at the time of its emission [56].

We will refer to a GW signal $h^{2}\Omega_{\text{GW}}(f)$ as being observable
by LISA, the Einstein Telescope [59, 60] or pulsar timing arrays [61], if it
achieves a sufficiently large signal-to-noise ratio. This means that the
signal will be observable if its peak amplitude
$h^{2}\Omega_{\text{GW}}^{\text{peak}}\equiv
h^{2}\Omega_{\text{GW}}(f_{\text{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].

###  2.5 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
$\lambda$, the Yukawa coupling of the dark fermion $y$ and the vacuum
expectation value of the dark Higgs $v_{\phi}$. 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 $\alpha$ and values of
$\beta/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
$\lambda$. Too small values of these two couplings would imply very large
values of $\beta/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 $\lambda$. As will be
seen below in section 3, successfully producing the right DM relic density
requires $m_{\chi}\gtrsim m_{\phi}$, which implies a lower limit
$y>2\sqrt{\lambda}$. Lastly, the vev $v_{\phi}$ 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\leq g\leq 1$,
$10^{-4}\leq\lambda\leq 10^{-2}$, $0.01\leq y\leq 0.7$ and $10^{-3}\text{
GeV}\leq v_{\phi}\leq 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_{\phi}$) 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. [62]), 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), $\lambda>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.

|  $g$ |  $\lambda$ |  $y$ |  $v_{\phi}$ |  $m_{\chi}$ |  $m_{\phi}$ |  $m_{A^{\prime}}$
---|---|---|---|---|---|---|---
Benchmark point |  $0.67$ |  $0.0035$ |  $0.62$ |  $430\,\text{GeV}$ |  $189\,\text{GeV}$ |  $36\,\text{GeV}$ |  $288\,\text{GeV}$
Table 1: Benchmark point used for discussing the thermalisation of visible and
dark sector.

##  3 Dark sector relic density

Figure 4:  The upper (lower) panels show the ratio of the dark Higgs boson
mass $m_{\phi}$ (the dark photon mass $m_{A^{\prime}}$) to the mass of the DM
fermion $m_{\chi}$ for $y=0.1$ (left) and $y=0.5$ (right), as a functions of
the gauge coupling $g$ and the self-interaction $\lambda$. Note that these
ratios are independent of the dark Higgs vev.

During the PT, the dark sector particles $\chi$, $\phi$ and $A^{\prime}$ all
obtain masses proportional to the dark Higgs vev $v_{\phi}$. In the parameter
regions of interest for a strong first-order PT, we generally find
$g>\sqrt{2\lambda}$ and $g>y/\sqrt{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 $\chi$ is thus determined through a freeze-out
process [63] 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: $\chi\chi\to\phi\phi$, $\chi\chi\to\phi
A^{\prime}$ amd $\chi\chi\to A^{\prime}A^{\prime}$. If the DM fermion is the
lightest particle in the dark sector, its annihilation cross section into
other dark sector states becomes exponentially suppressed at low temperatures;
a relic abundance in accordance with observations in this ‘forbidden’ regime
[26] thus implies 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_{\chi}\gtrsim m_{\phi}+m_{A^{\prime}}$, the
annihilation channel $\chi\chi\to\phi A^{\prime}$ 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,
$\chi\chi\to\phi\phi$, 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 [64],
which calculates the thermal averages and solves the full Boltzmann equation
[65]. While DarkSUSY allows precision calculations of the relic density in a
fully secluded dark sector with a varying temperature ratio $\xi$ between the
dark and the SM sector, cf. ref. [66], we will set $\xi=1$ for the purpose of
this section. We will revisit this assumption of thermal equilibrium between
the two sectors in section 4.

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

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 $\Omega_{\text{DM}}h^{2}$, the peak frequency $f_{\text{peak}}$
well as the peak amplitude $\Omega_{\text{GW}}^{\text{peak}}h^{2}$. One can
immediately see that $\Omega_{\text{DM}}$ and
$\Omega_{\text{GW}}^{\text{peak}}$ 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_{\phi}$ (indicated by the
colour of each point). A smaller value of $v_{\phi}$ 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_{\phi}$ 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_{\phi}$. 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
$\Omega_{\text{GW}}^{\text{peak}}h^{2}\approx 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
$\Omega_{\text{DM}}h^{2}\approx 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_{\phi}$ and the lower bound on $y$.

Figure 6: Scatter plots and 1D distributions as in figure 5, but with the
additional constraint $0.06\,\leq\Omega_{\text{DM}}h^{2}\leq 0.12$. Here, the
color scale does not encode the vev, but the ratio of percolation temperature
$T_{\text{p}}$ to the vev $v_{\phi}$, thus indicating the amount of
supercooling.

In figure 6 we show the result of sharpening the relic density requirement by
requiring that $0.06\,\leq\Omega_{\text{DM}}h^{2}\leq 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\,\mathrm{mHz}$ and $100\,\mathrm{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_{\rm peak}$, in particular, correspond to parameter
points in the ‘forbidden’ regime, $m_{\chi}<m_{\phi}$, where DM annihilations
are exponentially suppressed at small temperatures. Smaller values of
$v_{\phi}$ (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_{\rm 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_{\phi}$. 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.

##  4 Thermalisation of the two sectors

In this section we revisit the assumption that the temperature ratio of the
dark and visible sectors is $\xi=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. We then 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 allows 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 $\xi=1$, the parameters of this point lead to
$\alpha=0.258$, $\beta/H=874$, $T_{\text{n}}=39.7\,$GeV,
$T_{\text{p}}=39.1\,$GeV, $f_{\text{peak}}=3\,\text{mHz}$,
$\Omega_{\text{DM}}h^{2}=0.117$, and
$\Omega_{\text{GW}}^{\text{peak}}h^{2}=3\cdot 10^{-13}$. The rationale behind
choosing this benchmark point is that i) the observed DM relic abundance is
reproduced (for $\xi=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.

###  4.1 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
$\Omega_{\text{GW}}^{\text{peak}}h^{2}<10^{-10}$ and can therefore safely be
ignored. We also neglect the effect of bubble filtering [16, 67], 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_{\text{DS}}^{\text{br}}$ immediately after the PT. This temperature is
obtained from the temperature of the symmetric dark sector phase
$T_{\text{DS}}^{\text{sym}}$ using energy conservation:

$\rho_{\text{vac}}(T_{\text{DS}}^{\text{br}})+\rho_{\text{DS}}(T_{\text{DS}}^{\text{br}})=\rho_{\text{vac}}(T_{\text{DS}}^{\text{sym}})+\rho_{\text{DS}}(T_{\text{DS}}^{\text{sym}})\,,$
(4.1)

where $T_{\text{DS}}^{\text{br}}$ denotes the temperature in the broken phase
and

$\rho_{\text{DS}}(T_{\text{DS}}^{\text{br}})=\frac{\pi^{2}}{30}\,g_{\text{DS}}(T_{\text{DS}}^{\text{br}})\,\left(T_{\text{DS}}^{\text{br}}\right)^{4},$
(4.2)

where $g_{\text{DS}}(T)$ takes into account the $T$-dependence stemming both
from thermal (field-dependent) masses and the minimum of the effective
potential. Eq. (4.1) can easily be solved numerically for
$T_{\text{DS}}^{\text{br}}$ for a given $T_{\text{DS}}^{\text{sym}}$ .

In practice, we find that slightly different temperatures
$T_{\text{DS}}^{\text{br}}$ of the broken phase are obtained when solving the
equation taking $T_{\text{DS}}^{\text{sym}}=T_{\text{p}}$ or
$T_{\text{DS}}^{\text{sym}}=T_{\text{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 [68], 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_{\text{DS}}^{\text{sym}}=T_{\text{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.

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
$\xi=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’).

###  4.2 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 $\tau$ 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:

$\Gamma_{X}=\sum_{Y}\langle\sigma_{XY}v\rangle n_{Y}^{\rm eq}\,,$ (4.3)

where the sum is over all dark sector states $Y$, $\sigma_{XY}$ denotes the
total interaction cross section of $X$ and $Y$, brackets denote thermal
averaging (for simplicity calculated by assuming Boltzmann distributions) and
$n_{Y}^{\rm eq}$ 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 $\Gamma_{X}$ is only a few orders of magnitude smaller than $m_{X}$.
For example, dark Higgs bosons can thermalise via self-scattering, i.e.
$\phi\phi\leftrightarrow\phi\phi$, for which the scattering cross section is
$9\lambda^{2}/(8\pi s)$. For temperatures comparable to the dark Higgs boson
mass, we have $s\approx 4m_{\phi}^{2}$ and
$n_{\phi}\approx\zeta(3)m_{\phi}^{3}/\pi^{2}$, such that $\Gamma_{\phi}\sim
10^{-7}m_{\phi}$ for the benchmark point. Interactions of the dark Higgs
bosons with dark fermions or dark photons benefit from the larger couplings
$y,g\gg\lambda$, but suffer from a Boltzmann suppression if
$T_{\text{p}}<m_{\chi},m_{A^{\prime}}$.

A rough estimate of the thermalisation timescale is then obtained via

$\tau=\max_{X}\,\Gamma_{X}^{-1}\,.$ (4.4)

Given the timescale $\tau$ 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_{\text{t}}(t)=1-P(t)$. Its time derivative $\dot{P}_{\text{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

$F(t)\equiv P(t-\tau)-P(t)>0\,,$ (4.5)

which can hence be interpreted as the volume fraction
$\dot{P}_{\text{t}}\,\Delta t$ that just transitioned to the broken phase
within the small thermalization time scale $\Delta t=\tau$, cf. figure 7.
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_{\text{p}}$ for our benchmark point. As expected, we find that $F$
takes its maximal value close to the percolation temperature $T_{\text{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_{\text{p}}$ at the percolation
temperature for more points from our parameter scan, varying both $g$ and
$\lambda$. Similar conclusions are reached in all cases.

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
$\lambda$ and $g$ for $y=0.62$ and $v_{\phi}=430$ GeV.

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\phi\to 2\phi$, 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 [69, 70]. 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\approx
2\,\text{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.

###  4.3 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 [71, 72, 73, 74].111A second
possibility would be to consider kinetic mixing between the dark photon and SM
hypercharge. However, given that the dark photon is typically the heaviest
state in the dark sector, it will be strongly Boltzmann-suppressed at low
temperatures, and therefore cannot efficiently keep the two sectors in
equilibrium.  Such a mixing arises from an additional term in the scalar
potential:

$V_{\text{mix}}(H,\Phi)=\lambda_{h\phi}|H|^{2}|\Phi|^{2}\,.$ (4.6)

As long as the vevs of both Higgs bosons vanish, the dominant process
connecting the two sectors is $HH\to\Phi\Phi$. As soon as one of the two Higgs
bosons acquires a vev, it can decay into the other one, i.e. $h\to\Phi\Phi$
(if the electroweak symmetry breaks first) or $\phi\to HH$ (if the dark
symmetry breaks first). If kinematically allowed, these decay processes
typically dominate over the $2\to 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

$\theta=\frac{\lambda_{h\phi}v_{\phi}v_{h}}{m_{h}^{2}-m_{\phi}^{2}}\,,$ (4.7)

where we have assumed $\theta\ll 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 $\theta$. 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

$\Gamma_{\phi\to b\bar{b}}=\frac{3m_{\phi}m_{b}^{2}\sin^{2}\theta}{8\pi
v_{h}^{2}}\sqrt{1-\frac{4m_{b}^{2}}{m_{\phi}^{2}}}\,.$ (4.8)

To calculate the entropy transfer between the dark and visible sector, we
define the heat transfer rate

$\dot{q}\equiv\dot{\rho}+3H(\rho+P)\,,$ (4.9)

which is related to the change of entropy density via

$\dot{s}=-3Hs+\frac{\dot{q}}{T}\;.$ (4.10)

At the same time, the first moment of the Boltzmann equation gives

$\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}E\,L[f]=\dot{q}=\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}E\,C[f]\,.$
(4.11)

The general expression for the first moment of the collision operator for
decays (including relativistic corrections and quantum effects) was derived in
refs. [75, 76]. It was shown there that the leading relativistic effects
(namely a time dilation of the decay proportional to $1/\gamma$ and an
increase of the injected energy by a factor $\gamma$) 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

$\dot{q}\simeq m_{\phi}\left(\dot{n}+3Hn\right)\,,$ (4.12)

and the evolution of the number density is given by the usual Boltzmann
equation

$\dot{n}_{\phi}+3Hn_{\phi}=-\Gamma_{\phi\to
b\bar{b}}\,n^{\text{eq}}_{\phi}\left(\frac{n_{\phi}}{n^{\text{eq}}_{\phi}}-1\right).$
(4.13)

Here $n^{\text{eq}}$ denotes the equilibrium number density for
$T_{\text{DS}}=T_{\text{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

$\dot{s}_{\text{DS}}=-3Hs_{\text{DS}}-\frac{m_{\phi}}{T_{\text{DS}}}\,\Gamma_{\phi\to
b\bar{b}}\,n^{\text{eq}}_{\phi}\left(\frac{n_{\phi}}{n^{\text{eq}}_{\phi}}-1\right).$
(4.14)

A completely analogous equation holds for the evolution of the SM entropy
density $s_{\text{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 $\dot{a}=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\to\phi\phi$ and $\phi\to hh$ if they are
kinematically allowed. We do not, however, include $2\to 2$ processes of the
form $q\bar{q}\to g\phi$ or $q\phi\to qg$, which may give a non-negligible
contribution for light dark Higgs bosons [77, 78]. 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.

Figure 9: Evolution of the dark sector temperature (solid) and the SM
temperature (dashed) for different values of the initial temperature ratio
$\xi_{\text{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
non-relativistic ($T_{\text{DS}}=m_{\phi}/3$) and approximately when the DM
particles freezes out ($T_{\text{DS}}=m_{\chi}/25$). The different panels
correspond to different values of the portal coupling $\lambda_{h\phi}$.

In the following, we will take the temperature ratio of the two sectors at the
electroweak phase transition (EWPT) as initial value, $\xi_{\text{EWPT}}$, and
calculate the evolution of $\xi$ during the subsequent cosmological stages
(see appendix B for details). In figure 9 we show the visible and dark sector
temperatures as a function of scale factor for different initial values of
$\xi_{\text{EWPT}}$, defined at $T_{\text{SM}}=150\,\text{GeV}$. In the four
panels the portal coupling was set to the representative values
$\lambda_{h\phi}=10^{-6}$ (top left), $3\cdot 10^{-7}$ (top right), $10^{-7}$
(bottom left) and $3\cdot 10^{-8}$ (bottom right). In each panel, we indicate
the moment of percolation by a vertical dot-dashed line; the approximate
temperature of the dark Higgs becoming non-relativistic
($m_{\phi}/T_{\text{DS}}=3$) and the DM fermion freezing out
($m_{\chi}/T_{\text{DS}}=25$) are indicated by horizontal dotted and dashed
lines, respectively. We make the following observations:

  * •

For $\lambda_{h\phi}=10^{-6}$, the two sectors thermalise efficiently already
before percolation. The initial value of $\xi_{\text{EWPT}}$ is therefore
inconsequential for the subsequent evolution, and we obtain the same results
for all cases.

  * •

For $\lambda_{h\phi}=3\times 10^{-7}$, the two sectors do exchange energy and
entropy already in the unbroken phase, but do not fully thermalise before
percolation, such that $\xi_{\text{p}}$ depends on $\xi_{\text{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 $\xi_{\text{EWPT}}$.

  * •

For $\lambda_{h\phi}=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 $\lambda_{h\phi}=3\times 10^{-8}$, the two sectors do not quickly
thermalize after the PT, and the universe enters an early period of cannibal
domination.222For the purpose of this plot, we assume that number changing
processes in the dark sector remain efficient throughout, such that the
chemical potential of the dark Higgs boson vanishes and (in the absence of
decays) the dark sector temperature grows relative to the SM temperature.

In general the value of $\lambda_{h\phi}$ needed to ensure thermalisation
depends somewhat on the dark Higgs vev, since for $m_{\phi}\ll m_{h}$ the
mixing angle scales as $\theta\propto\lambda_{h\phi}v_{\phi}$. 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 $\xi=1$ made
in section 3 is well justified for $\lambda_{h\phi}$ 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
$\lambda_{h\phi}\gtrsim 10^{-3}$ [32].

##  5 Gravitational waves from hot dark sectors

While the analysis of the dark sector PT is simplest for
$\lambda_{h\phi}>10^{-6}$, it is phenomenologically interesting to also
consider smaller values of $\lambda_{h\phi}$, such that the temperature ratio
of the two sectors before the PT may differ from unity. The reason is that
$\xi_{\text{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 $\lambda_{h\phi}=3\times 10^{-8}$ the temperature
ratio $\xi$ still differs from unity when the lightest dark sector particle
becomes non-relativistic. 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. [79]). 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 $\xi=1$) to be
valid.

###  5.1 Dilution of gravitational waves

Figure 10: Dilution factor $D$, cf. eq. (5.5), as a function of the portal
coupling $\lambda_{h\phi}$, and for different values of the initial
temperature ratio $\xi_{\text{EWPT}}$ as indicated.

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 [80]

$\displaystyle g_{\text{tot},\text{p}}$
$\displaystyle=g_{\text{SM},\text{p}}+g_{\text{DS},\text{p}}\,\xi_{\text{p}}^{4}\,,$
(5.1) $\displaystyle h_{\text{tot},\text{p}}$
$\displaystyle=h_{\text{SM},\text{p}}+h_{\text{DS},\text{p}}\,\xi_{\text{p}}^{3}\,,$
(5.2)

where the temperature ratio at percolation $\xi_{\text{p}}$ is derived from
the initial temperature ratio $\xi$ 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 $\Omega_{\text{GW}}$ observed today can be calculated from
the emitted spectrum $\Omega_{\text{GW}}^{\text{em}}$ as

$\displaystyle\Omega_{\text{GW}}(f)$
$\displaystyle=\left(\frac{a_{\text{p}}}{a_{0}}\right)^{4}\left(\frac{H_{\text{p}}}{H_{0}}\right)^{2}\Omega_{\text{GW}}^{\text{em}}\left(\frac{a_{0}}{a_{\text{p}}}\,f\right)\equiv\mathcal{R}\,\Omega_{\text{GW}}^{\text{em}}\left(\frac{a_{0}}{a_{\text{p}}}\,f\right).$
(5.3)

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 [81]. 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

$\displaystyle\frac{a_{0}}{a_{\text{p}}}=\left(\frac{h_{\text{SM},\text{p}}}{h_{\text{SM},0}}\right)^{1/3}\left(\frac{T_{\text{SM},\text{p}}}{T_{\text{SM},0}}\right).$
(5.4)

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_{\text{p}}$ in eq. (5.4) needs to be corrected by a
factor of $D_{\text{SM}}^{1/3}$, where
$D_{\text{SM}}=S_{\text{SM},0}/S_{\text{SM},\text{p}}$ [82, 11]. We hence
find, more generally, that

$\displaystyle\mathcal{R}h^{2}$
$\displaystyle=\frac{\Omega_{\text{rad}}h^{2}}{D_{\text{SM}}^{4/3}}\,\left(\frac{h_{\text{SM},0}}{h_{\text{SM},\text{p}}}\right)^{4/3}\frac{g_{\text{tot},\text{p}}}{2}=\frac{\Omega_{\text{rad}}h^{2}}{D^{4/3}}\,\left(\frac{h_{\text{SM},0}}{h_{\text{tot},\text{p}}}\right)^{4/3}\frac{g_{\text{tot},\text{p}}}{2}\,,$
(5.5)

where we introduced the more convenient dilution factor $D\equiv
D_{\text{SM}}\,h_{\text{SM},\text{p}}/h_{\text{tot},\text{p}}$ [11, 82].
Analogously, the peak frequency obtains an additional redshift for the case of
a non-conserved SM entropy and hence reads

$\displaystyle
f_{\text{peak}}=\frac{8.9\,\text{mHz}}{D^{1/3}}\left(\frac{T_{\text{SM},\text{p}}}{100\text{
GeV}}\right)\left(\frac{\beta/H}{1000}\right)\left(\frac{g_{\text{tot},\text{p}}}{100}\right)^{1/2}\left(\frac{100}{h_{\text{tot},\text{p}}}\right)^{1/3}\,.$
(5.6)

We show in figure 10 the dilution factor $D$, as a function of the portal
coupling $\lambda_{h\phi}$. 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 $\xi_{\text{EWPT}}$, defined at
$T_{\text{SM}}=150\,\text{GeV}$. As expected, for sufficiently large
$\lambda_{h\phi}$ there is no significant dilution333In the limit of large
portal couplings $\lambda_{h\phi}$, the dilution factor $D$ approaches a value
slightly larger than 1. This is an expected feature, indicating a negligible
dilution effect that is entirely due to the additional degrees of freedom in
the combined thermal bath of SM and dark sector particles and not a
consequence of additional entropy injection (with respect to $\Lambda$CDM)
into the SM bath after the PT. , as entropy is conserved and only radiation
degrees of freedom contribute to the energy content of the universe. However,
with decreasing $\lambda_{h\phi}$ this is no longer the case and the dilution
factor grows, becoming sizeable for $\lambda_{h\phi}\ll 10^{-7}$.

###  5.2 Results

Figure 11: Gravitational wave signal for the same scenarios as considered in
figure 9, i.e. for various values of the portal coupling $\lambda_{h\phi}$
(different panels) and the initial temperature ratio $\xi_{\text{EWPT}}$
(different colours).

We are finally in the position to calculate the GW signal for $\xi\neq 1$ and
$\lambda_{h\phi}<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
$\xi_{\text{EWPT}}>1$, provided that $\lambda_{h\phi}$ is sufficiently small
for the sectors to not equilibrate before the PT (cf. figure 9). The
enhancement saturates for $\xi_{\text{EWPT}}\gtrsim 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
$\lambda_{h\phi}\lesssim 3\times 10^{-8}$, on the other hand, the effect of
dilution becomes relevant and the GW signal starts to become suppressed.
Crucially, changing $\xi$ and $\lambda_{h\phi}$ does not significantly affect
the peak frequency, such that the GW signal remains within the LISA
sensitivity range.

Figure 12: Same as figure 6, but without the assumption of thermal equilibrium
between the two sectors. Specifically, we consider $\xi_{\text{EWPT}}=2$ and
$\lambda_{h\phi}=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.

We can now test the robustness of our results from section 3, where we assumed
$\xi=1$ and $\lambda_{h\phi}>10^{-6}$, by allowing larger values of $\xi$ and
smaller values of $\lambda_{h\phi}$. In figure 12 we show the same result as
in figure 6, but now for $\xi_{\text{EWPT}}=2$ and $\lambda_{h\phi}=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 $\xi_{\text{nr}}>1.1$ at
$T_{\text{DS},\text{nr}}=m_{\phi}/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_{\text{peak}}$ and
$\Omega_{\text{GW}}^{\text{peak}}h^{2}$ look very similar with and without
this additional requirement.

The parameter combination $\xi_{\text{EWPT}}=2$ and $\lambda_{h\phi}=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
one-dimensional 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 observable events
increases from 35% to 69%. We summarize our findings in table 2.

| Fraction of parameter points observable by LISA
---|---
|  $\xi_{\text{EWPT}}\\!=\\!1$, $\lambda_{h\phi}\\!=\\!10^{-6}$ |  $\xi_{\text{EWPT}}\\!=\\!2$, $\lambda_{h\phi}\\!=\\!10^{-7}$
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.

We emphasize that this large increase is a result of fixing $\xi$ and
$\lambda_{h\phi}$ to particular values. If we instead vary $\xi$ and
$\lambda_{h\phi}$ as part of the scan, most parameter combinations will either
give very similar results to the case $\xi=1$ considered in section 3 or lead
to an extended period where the dark sector energy density dominates.

##  6 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)^{\prime}$ 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
$\lambda$, the dark Yukawa coupling $y$ and the dark Higgs vev $v_{\phi}$. 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 $\alpha$) and not too fast ($\beta/H\sim 10^{2}$–$10^{4}$) PT,
corresponding to potentially observable GW signals. We showed in particular
that large GW signals require sizeable couplings $g$ and $\lambda$ 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_{\phi}\sim v_{\phi}$, which in turn implies a sizable Yukawa coupling and a
DM mass that is comparable to the dark Higgs vev $v_{\phi}$. 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_{\text{p}}/T_{\text{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 $\lambda_{h\phi}$. 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 $\lambda_{h\phi}>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
$\xi_{\text{EWPT}}$ as a free parameter. While small values of
$\lambda_{h\phi}$ combined with large initial values of $\xi_{\text{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
$\lambda_{h\phi}\approx 10^{-7}$ a net enhancement remains, demonstrating that
it is possible to have $\xi>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 $\xi=2$ and
$\lambda_{h\phi}=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 [83, 84,
85]. 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 non-relativistic 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 [86] such as CTA [87]. 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.

## Acknowledgments

We thank Thomas Konstandin and Kai Schmidt-Hoberg for helpful comments on the
manuscript and Peter Athron, Csaba Balazs, Andreas Ekstedt, Fatih Ertas, Henda
Mansour, Lachlan Morris and Margarethe Mühlleitner for discussions. CT thanks
FNAL for its hospitality. This project has received funding from the European
Union’s Horizon Europe research and innovation programme under the Marie
Skłodowska-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.

##  Appendix A Dark bubble walls

In addition to the strength $\alpha$ of the PT, its speed $\beta/H$, the
percolation temperature $T_{\text{p}}$ and the temperature ratio
$\xi_{\text{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_{\text{eff}}(\phi,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 [88, 89, 51, 90, 91]. 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 [53].

In this paper we take the approach suggested in ref. [92], 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 [89], which relates the velocity-
independent leading-order (LO) bubble wall friction
$\mathcal{P}_{\mathrm{LO}}$ to the amount of liberated vacuum energy density
$\Delta V_{\text{eff}}$:

$\displaystyle\text{Bödeker-Moore criterion:}\quad\begin{cases}\Delta
V_{\text{eff}}>\mathcal{P}_{\mathrm{LO}}\quad\text{Relativistic bubble
walls}\\\ \Delta V_{\text{eff}}<\mathcal{P}_{\mathrm{LO}}\quad\text{Non-
relativistic bubble walls}\end{cases}\,.$ (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 $\mathcal{P}_{\mathrm{NLO}}$, which scales with powers of
$\gamma_{\text{w}}=1/\sqrt{1-v_{w}^{2}}$ [51, 93, 52]. Bubbles can run away
only if $\Delta
V_{\text{eff}}>\mathcal{P}_{\mathrm{LO}}+\mathcal{P}_{\mathrm{NLO}}$ for all
$v_{w}$. Otherwise they will reach a relativistic terminal velocity given by
the equilibrium of forces, $\Delta
V_{\text{eff}}=\mathcal{P}_{\mathrm{LO}}+\mathcal{P}_{\mathrm{NLO}}$.

The LO friction due to particles acquiring a mass when crossing the bubble
wall is given by [89, 57]

$\displaystyle\mathcal{P}_{\mathrm{LO}}\simeq\sum_{i}g_{i}\,c_{i}\frac{\Delta
m_{i}^{2}}{24}T_{\text{n}}^{2}\,,$ (A.2)

where $\Delta m_{i}^{2}$ 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 $\gamma_{\text{w}}T_{\text{n}}$ can also add to the
LO friction [90]. 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)^{\prime}$ with a fermionic species, the LO
friction reads [57]

$\displaystyle\mathcal{P}_{\text{LO}}\simeq\left(3m_{A^{\prime}}^{2}+m_{\phi}^{2}+2m_{\chi}^{2}\right)\,\frac{{T_{\text{p}}}^{2}}{24}\approx\frac{m_{A^{\prime}}^{2}\,{T_{\text{p}}}^{2}}{8}\,,$
(A.3)

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 $\Delta
V_{\text{eff}}$ can be estimated from the zero-temperature potential, which
gives

$\Delta
V_{\text{eff}}\approx\frac{\lambda}{4}v_{\phi}^{4}=\frac{m_{\phi}^{2}}{8}v_{\phi}^{2}\,.$
(A.4)

Hence we find that the Bödeker-Moore criterion for relativistic bubble walls
is satisfied if $m_{\phi}/m_{A^{\prime}}>{T_{\text{p}}}/v_{\phi}$, which is
the case for the parameter regions that give strongly supercooled PTs
corresponding to observable GW signals ($T_{\text{p}}\lesssim 0.1\,v_{\phi}$).
In these parameter regions, the bubble walls are therefore expected to be
relativistic, $v_{w}\rightarrow 1$. This finding also implies that we can
neglect the effect of bubble filtering, which is only relevant in the (deeply)
non-relativistic regime $v_{w}\ll 1$ [16, 67]. 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 $\gamma_{\text{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 $\gamma_{\text{w}}^{\text{terminal}}$
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. [94, 95, 96]. In ref. [94] 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.

Figure 13: The leading-order friction $\mathcal{P}_{\text{LO}}$ over the
difference $\Delta V_{\text{eff}}$ in potential energy between the true and
false vacuum phases as a function of the quartic coupling $\lambda$ 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.

##  Appendix 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.7) 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. [75] 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 $\Gamma_{h}$ of the
Higgs resonance. As shown in ref. [78], including the width leads to an
effective mixing angle given by

$\displaystyle\theta_{\mathrm{eff}}^{2}(T)=\frac{\left(\lambda_{p}v_{h}(T)v_{\phi}(T)\right)^{2}}{(m_{\phi}^{2}(T)-m_{h}^{2}(T))^{2}+(m_{\phi}(T)\Gamma_{h})^{2}}\,.$
(B.1)

###  B.2 The Boltzmann equation for entropy transfer

In our analysis we specify the initial conditions, i.e. the temperature ratio
$\xi_{\text{EWPT}}$ of the dark sector and SM bath at the EWPT, for which we
take $T_{\text{EWPT}}=150\,$GeV. Which processes contribute to the
thermalisation of the two sectors depends on whether or not the
$U(1)^{\prime}$ 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\leftrightarrow\Phi\Phi$
and the decay of the SM Higgs $h\rightarrow\Phi\Phi$. 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\leftrightarrow\phi\phi$, the decays of the SM Higgs into dark Higgs
$h\rightarrow\phi\phi$ if it is kinematically allowed, and the decays of the
dark Higgs into SM particles $\phi\rightarrow\mathrm{SM}\ \mathrm{SM}$ through
Higgs mixing.

In the following we give the relevant expressions contributing to the entropy
transfer.

####  B.2.1 $2\to 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\to 2$ process
induced by the portal coupling $\Phi\Phi\leftrightarrow 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. [97, 98], which we briefly sketch here. The heat transfer for $2\to 2$
processes can be expressed as

$\displaystyle\dot{q}\big{|}_{2\to 2}$
$\displaystyle=\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}EC[f]$
$\displaystyle=\int\frac{\mathrm{d}^{3}p_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}p_{2}}{(2\pi)^{3}2E_{2}}\frac{\mathrm{d}^{3}k_{1}}{(2\pi)^{3}2E_{k_{1}}}\frac{\mathrm{d}^{3}k_{2}}{(2\pi)^{3}2E_{k_{2}}}|\mathcal{M}|^{2}(2\pi)^{4}\delta(\Sigma
p)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
f(p_{1})f(p_{2})(1+f(k_{1})(1+f(k_{2}))$
$\displaystyle=\int\frac{\mathrm{d}^{3}p_{1}}{(2\pi)^{3}2E_{1}}\frac{\mathrm{d}^{3}p_{2}}{(2\pi)^{3}2E_{2}}8E_{1}F(p_{1},p_{2})\sigma(p_{1},p_{2})f(p_{1})f(p_{2})\,.$
(B.2)

Here the final state statistical factors are absorbed into the cross section
and $F(p_{1},p_{2})=\sqrt{(p_{1}\cdot p_{2})^{2}-m^{2}_{1}m_{2}^{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 $\Lambda$ 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 $\eta$ and two angles $\theta$ and $\varphi$;
for details see ref. [97]. The phase space distribution becomes

$\displaystyle f(k)=\frac{1}{e^{u\cdot k/T}\mp
1}\overset{\Lambda}{\rightarrow}f_{\Lambda}(k)=\frac{1}{e^{\left(k_{0}\cosh\eta+k_{1}\sinh\eta\right)/T}\mp
1}\,.$ (B.3)

With this, we can rewrite the center of mass cross section as

$\displaystyle\sigma_{\mathrm{CM}}(p_{1},p_{2})=\frac{1}{(2\pi)^{2}16F(p_{1},p_{2})}\int\mathrm{d}\varphi\mathrm{d}\cos\theta\left|\mathcal{M}\right|^{2}\frac{\sqrt{E^{2}-m_{f}^{2}}}{E}\left((1+f_{\Lambda}(k_{1}))(1+f_{\Lambda}(k_{2}))\right)\,.$
(B.4)

The matrix element for the processes $\Phi\Phi\to hh$ is at tree level simply
given by $\mathcal{M}=i\lambda_{h\phi}$. For this angle-independent transition
amplitude, we can integrate over the solid angle and obtain

$\displaystyle\sigma_{\mathrm{CM}}=\frac{|\mathcal{M}|^{2}}{64\pi
E^{2}}\frac{T}{\sqrt{E^{2}-m_{\Phi}^{2}}\sinh{\eta}}\frac{1}{1-e^{-\frac{2E}{T}\cosh\eta}}\ln\underbrace{\left[\frac{\sinh\left(\frac{E\cosh\eta+\sqrt{E^{2}-m_{h}^{2}}\sinh\eta}{2T}\right)}{\sinh\left(\frac{E\cosh\eta-\sqrt{E^{2}-m_{h}^{2}}\sinh\eta}{2T}\right)}\right]}_{\equiv\lambda(E,\eta,T,m_{h})}\,.$
(B.5)

In the case of initial and final states with respectively equal masses, eq.
(B.2.1) reduces to

$\displaystyle I_{\Phi\Phi\rightarrow
hh}\equiv\frac{2T_{\mathrm{DS}}}{\pi^{4}}\int_{m_{\Phi}(T)}^{\infty}\mathrm{d}E\sqrt{E^{2}-m_{\Phi}^{2}}E^{4}\int_{0}^{\infty}\mathrm{d}\eta\frac{\sinh\eta\cosh\eta}{e^{2E\cosh\eta/T_{\mathrm{DS}}}-1}\ln\lambda(E,\eta,T_{\mathrm{DS}},m_{\Phi})~{}\sigma_{\mathrm{CM}}\,.$
(B.6)

This expression can now be efficiently evalutated numerically. An analogous
expression is obtained for the process $hh\rightarrow\phi\phi$.

We note that in principle there are additional $2\to 2$ processes that may
contribute to thermalisation. The process $\phi\phi\to t\bar{t}$ via an off-
shell SM Higgs boson is strongly Boltzmann-suppressed below the EWPT [75] and
does not give a relevant contribution in the temperature range that we
consider. However, processes such as $\phi+q\to g+q$ with a quark in the
$t$-channel can give a relevant contribution if the decay $\phi\to b\bar{b}$
is kinematically forbidden [77, 78]. 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 pair-annihilation of dark Higgs bosons into
predominantly bottom quarks: $\phi\phi\to h\to b\bar{b}$. In thermal
equilibrium, the rate of this process can be related to the inverse process,
which is the decay $h\to\phi\phi$ with partial width given by

$\Gamma_{h\to\phi\phi}=\frac{(m_{h}^{2}+2\,m_{\phi}^{2})^{2}\,\sin^{2}2\theta_{\mathrm{eff}}}{128\pi\,m_{h}}\left(1-\frac{4\,m_{\phi}^{2}}{m_{h}^{2}}\right)^{1/2}\left(\frac{1}{v_{\phi}}\cos\theta_{\mathrm{eff}}+\frac{1}{v_{h}}\sin\theta_{\mathrm{eff}}\right)^{2}\,.$
(B.7)

This gives the additional term in the heat transfer rate

$\displaystyle\dot{q}\big{|}_{h\to\phi\phi}=-m_{h}\Gamma_{h\to\phi\phi}n^{\text{eq}}_{h}\left[1-\left(\frac{n_{\phi}}{n^{\text{eq}}_{\phi}}\right)^{2}\right]\;.$
(B.8)

Before the dark sector symmetry breaking, we also have the decay
$h\rightarrow\Phi\Phi$ with the decay rate

$\displaystyle\Gamma_{h\to\Phi\Phi}=\frac{\lambda_{h\phi}^{2}v_{h}^{2}}{128\pi\,m_{h}}\,,$
(B.9)

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

$\displaystyle\dot{s}_{\mathrm{DS}}+3Hs_{\mathrm{DS}}=-\frac{m_{h}}{T_{\mathrm{DS}}}\Gamma_{h\rightarrow\Phi\Phi}n_{h}^{\mathrm{eq}}\left[\left(\frac{n_{\phi}}{n_{\phi}^{\mathrm{eq}}}\right)^{2}-1\right]-\frac{1}{T_{\text{DS}}}I_{\Phi\Phi\to
hh}+\frac{1}{T_{\text{DS}}}I_{hh\to\Phi\Phi}\,.$ (B.10)

The equation for the case of the PT occuring before the EWPT follows
analogously. After both symmetries are broken the full Boltzmann equation
reads

$\displaystyle\dot{s}_{\mathrm{DS}}+3Hs_{\mathrm{DS}}$
$\displaystyle=-\frac{m_{\phi}}{T_{\mathrm{DS}}}\Gamma_{\phi\to\mathrm{SM}\mathrm{SM}}n^{\text{eq}}_{\phi}\left(\frac{n_{\phi}}{n^{\text{eq}}_{\phi}}-1\right)+\frac{m_{h}}{T_{\mathrm{DS}}}\Gamma_{h\to\phi\phi}n^{\text{eq}}_{h}\left[1-\left(\frac{n_{\phi}}{n^{\text{eq}}_{\phi}}\right)^{2}\right]$
$\displaystyle-\frac{1}{T_{\text{DS}}}I_{\phi\phi\to
hh}+\frac{1}{T_{\text{DS}}}I_{hh\to\phi\phi}\,.$ (B.11)

where we use the decay widths of the SM Higgs into other SM particle pairs
from ref. [32]. 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):

$\displaystyle(\sigma v)_{\chi\chi\to A^{\prime}A^{\prime}}=$
$\displaystyle\frac{m_{A^{\prime}}^{4}(m_{\chi}^{2}-m_{A^{\prime}}^{2})^{3/2}}{64\pi\,v_{\phi}^{4}\,m_{\chi}(m_{A^{\prime}}^{2}-2m_{\chi}^{2})^{2}}$
$\displaystyle+\frac{v^{2}\sqrt{m_{\chi}^{2}-m_{A^{\prime}}^{2}}}{384\pi\,v_{\phi}^{4}\,m_{\chi}(m_{\phi}^{2}-4m_{\chi}^{2})^{2}(m_{A^{\prime}}^{2}-2m_{\chi}^{2})^{4}}$
$\displaystyle\quad\times\Big{[}144m_{A^{\prime}}^{12}m_{\chi}^{2}+2m_{A^{\prime}}^{10}(7m_{\phi}^{4}-88m_{\phi}^{2}m_{\chi}^{2}-432m_{\chi}^{4})$
$\displaystyle\qquad+128m_{\chi}^{10}(m_{\phi}^{4}+8m_{\chi}^{4})-64m_{A^{\prime}}^{2}m_{\chi}^{8}(3m_{\phi}^{4}+16m_{\phi}^{2}m_{\chi}^{2}+32m_{\chi}^{4})$
$\displaystyle\qquad+4m_{A^{\prime}}^{4}m_{\chi}^{6}(17m_{\phi}^{4}+600m_{\phi}^{2}m_{\chi}^{2}+592m_{\chi}^{4})$
$\displaystyle\qquad+m_{A^{\prime}}^{8}(-73m_{\phi}^{4}m_{\chi}^{2}+1128m_{\phi}^{2}m_{\chi}^{4}+1840m_{\chi}^{6})$
$\displaystyle\qquad+4m_{A^{\prime}}^{6}(25m_{\phi}^{4}m_{\chi}^{4}-648m_{\phi}^{2}m_{\chi}^{6}-496m_{\chi}^{8})\Big{]}$
$\displaystyle+\mathcal{O}\left(v^{4}\right)$ (B.12) $\displaystyle(\sigma
v)_{\chi\chi\to\phi\phi}=$
$\displaystyle\frac{v^{2}m_{\chi}\sqrt{m_{\chi}^{2}-m_{\phi}^{2}}}{192\pi\,v_{\phi}^{4}(m_{\phi}^{2}-4m_{\chi}^{2})^{2}(m_{\phi}^{2}-2m_{\chi}^{2})^{4}}\left(3m_{\phi}^{4}-8m_{\phi}^{2}m_{\chi}^{2}+8m_{\chi}^{4}\right)$
$\displaystyle\quad\times\left(9m_{\phi}^{8}-64m_{\phi}^{6}m_{\chi}^{2}+200m_{\phi}^{4}m_{\chi}^{4}-352m_{\phi}^{2}m_{\chi}^{6}+288m_{\chi}^{8}\right)$
$\displaystyle+\mathcal{O}\left(v^{4}\right)$ (B.13) $\displaystyle(\sigma
v)_{\chi\chi\to A^{\prime}\phi}=$
$\displaystyle\frac{(m_{A^{\prime}}^{4}-2m_{A^{\prime}}^{2}m_{\phi}^{2}+m_{\phi}^{4}-8m_{A^{\prime}}^{2}m_{\chi}^{2}-8m_{\phi}^{2}m_{\chi}^{2}+16m_{\chi}^{4})^{3/2}}{2048\pi
v^{4}m_{\chi}^{4}}$
$\displaystyle+\frac{v^{2}\sqrt{m_{A^{\prime}}^{4}+(m_{\phi}^{2}-4m_{\chi}^{2})^{2}-2m_{A^{\prime}}^{2}(m_{\phi}^{2}+4m_{\chi}^{2})}}{12288\pi\,v_{\phi}^{4}m_{\chi}^{4}(m_{A^{\prime}}^{2}-4m_{\chi}^{2})^{2}(m_{A^{\prime}}^{2}+m_{\phi}^{2}-4m_{\chi}^{2})^{4}}$
$\displaystyle\quad\times\Big{[}-2m_{A^{\prime}}^{14}(11m_{\phi}^{2}-228m_{\chi}^{2})-16m_{\chi}^{4}(m_{\phi}^{2}-4m_{\chi}^{2})^{4}(15m_{\phi}^{4}-80m_{\phi}^{2}m_{\chi}^{2}+16m_{\chi}^{4})$
$\displaystyle\qquad-11m_{A^{\prime}}^{16}-m_{A^{\prime}}^{12}(-11m_{\phi}^{4}-1200m_{\phi}^{2}m_{\chi}^{2}+7264m_{\chi}^{4})$
$\displaystyle\qquad+4m_{A^{\prime}}^{10}(11m_{\phi}^{6}+266m_{\phi}^{4}m_{\chi}^{2}-4328m_{\phi}^{2}m_{\chi}^{4}+15136m_{\chi}^{6})$
$\displaystyle\qquad+8m_{A^{\prime}}^{2}m_{\chi}^{2}(m_{\phi}^{2}-4m_{\chi}^{2})^{2}(15m_{\phi}^{8}-244m_{\phi}^{6}m_{\chi}^{2}+1280m_{\phi}^{4}m_{\chi}^{4}-2880m_{\phi}^{2}m_{\chi}^{6}+6912m_{\chi}^{8})$
$\displaystyle\qquad-
m_{A^{\prime}}^{4}(m_{\phi}^{2}-4m_{\chi}^{2})^{2}(11m_{\phi}^{8}-344m_{\phi}^{6}m_{\chi}^{2}+1808m_{\phi}^{4}m_{\chi}^{4}-8704m_{\phi}^{2}m_{\chi}^{6}+80384m_{\chi}^{8})$
$\displaystyle\qquad-
m_{A^{\prime}}^{8}(-11m_{\phi}^{8}-32m_{\phi}^{6}m_{\chi}^{2}+12464m_{\phi}^{4}m_{\chi}^{4}-114688m_{\phi}^{2}m_{\chi}^{6}+291840m_{\chi}^{8})$
$\displaystyle\qquad-2m_{A^{\prime}}^{6}(11m_{\phi}^{10}-140m_{\phi}^{8}m_{\chi}^{2}+1216m_{\phi}^{6}m_{\chi}^{4}-29056m_{\phi}^{4}m_{\chi}^{6}+201984m_{\phi}^{2}m_{\chi}^{8}-412672m_{\chi}^{10})\Big{]}$
$\displaystyle+\mathcal{O}\left(v^{4}\right)$ (B.14)

## References

  * [1] Particle Data Group Collaboration, R. L. Workman et al., Review of Particle Physics,  PTEP 2022 (2022) 083C01.
  * [2] Planck Collaboration, N. Aghanim et al., Planck 2018 results. vi. cosmological parameters,  1807.06209.
  * [3] M. Pospelov, A. Ritz, and M. B. Voloshin, Secluded WIMP Dark Matter,  Phys. Lett. B 662 (2008) 53–61, [0711.4866].
  * [4] LISA Collaboration, P. Amaro-Seoane et al., Laser Interferometer Space Antenna, 1702.00786.
  * [5] C. Caprini et al., Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions,  JCAP 04 (2016) 001, [1512.06239].
  * [6] C. Caprini and D. G. Figueroa, Cosmological Backgrounds of Gravitational Waves,  Class. Quant. Grav. 35 (2018), no. 16 163001, [1801.04268].
  * [7] C. Caprini et al., Detecting gravitational waves from cosmological phase transitions with LISA: an update,  JCAP 03 (2020) 024, [1910.13125].
  * [8] J. R. Espinosa, T. Konstandin, J. M. No, and M. Quiros, Some Cosmological Implications of Hidden Sectors,  Phys. Rev. D 78 (2008) 123528, [0809.3215].
  * [9] P. Schwaller, Gravitational Waves from a Dark Phase Transition,  Phys. Rev. Lett. 115 (2015), no. 18 181101, [1504.07263].
  * [10] M. Breitbach, J. Kopp, E. Madge, T. Opferkuch, and P. Schwaller, Dark, Cold, and Noisy: Constraining Secluded Hidden Sectors with Gravitational Waves,  JCAP 07 (2019) 007, [1811.11175].
  * [11] F. Ertas, F. Kahlhoefer, and C. Tasillo, Turn up the volume: listening to phase transitions in hot dark sectors,  JCAP 02 (2022), no. 02 014, [2109.06208].
  * [12] T. Bringmann, P. F. Depta, T. Konstandin, K. Schmidt-Hoberg, and C. Tasillo, Does NANOGrav observe a dark sector phase transition?, 2306.09411.
  * [13] T. Alanne, K. Tuominen, and V. Vaskonen, Strong phase transition, dark matter and vacuum stability from simple hidden sectors,  Nucl. Phys. B 889 (2014) 692–711, [1407.0688].
  * [14] I. Baldes, Gravitational waves from the asymmetric-dark-matter generating phase transition,  JCAP 05 (2017) 028, [1702.02117].
  * [15] E. Hall, T. Konstandin, R. McGehee, and H. Murayama, Asymmetric matter from a dark first-order phase transition,  Phys. Rev. D 107 (2023), no. 5 055011, [1911.12342].
  * [16] M. J. Baker, J. Kopp, and A. J. Long, Filtered Dark Matter at a First Order Phase Transition,  Phys. Rev. Lett. 125 (2020), no. 15 151102, [1912.02830].
  * [17] A. Azatov, M. Vanvlasselaer, and W. Yin, Dark Matter production from relativistic bubble walls,  JHEP 03 (2021) 288, [2101.05721].
  * [18] M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht, Primordial Black Holes from First-Order Cosmological Phase Transitions, 2105.07481.
  * [19] M. Kierkla, A. Karam, and B. Swiezewska, Conformal model for gravitational waves and dark matter: a status update,  JHEP 03 (2023) 007, [2210.07075].
  * [20] M. Lewicki, K. Müürsepp, J. Pata, M. Vasar, V. Vaskonen, et al., Dynamics of false vacuum bubbles with trapped particles,  Phys. Rev. D 108 (2023), no. 3 036023, [2305.07702].
  * [21] M. Lewicki, P. Toczek, and V. Vaskonen, Primordial black holes from strong first-order phase transitions,  JHEP 09 (2023) 092, [2305.04924].
  * [22] S. Kanemura and S.-P. Li, Dark phase transition from WIMP: complementary tests from gravitational waves and colliders, 2308.16390.
  * [23] D. Croon, T. E. Gonzalo, and G. White, Gravitational Waves from a Pati-Salam Phase Transition,  JHEP 02 (2019) 083, [1812.02747].
  * [24] P. Athron, C. Balázs, T. E. Gonzalo, and M. Pearce, Falsifying Pati-Salam models with LIGO, 2307.02544.
  * [25] E. W. Kolb and M. S. Turner, The Early Universe, vol. 69.  1990\.
  * [26] R. T. D’Agnolo and J. T. Ruderman, Light Dark Matter from Forbidden Channels,  Phys. Rev. Lett. 115 (2015), no. 6 061301, [1505.07107].
  * [27] S. Weinberg, Mass of the Higgs Boson,  Phys. Rev. Lett. 36 (1976) 294–296.
  * [28] M. Fairbairn, E. Hardy, and A. Wickens, Hearing without seeing: gravitational waves from hot and cold hidden sectors,  JHEP 07 (2019) 044, [1901.11038].
  * [29] B. Holdom, Two U(1)’s and $\epsilon$ charge shifts,  Physics Letters B 166 (Jan., 1986) 196–198.
  * [30] P. Ilten, Y. Soreq, M. Williams, and W. Xue, Serendipity in dark photon searches,  JHEP 06 (2018) 004, [1801.04847].
  * [31] M. Bauer, P. Foldenauer, and J. Jaeckel, Hunting All the Hidden Photons,  JHEP 07 (2018) 094, [1803.05466].
  * [32] T. Ferber, A. Grohsjean, and F. Kahlhoefer, Dark Higgs Bosons at Colliders, 2305.16169.
  * [33] M. Quiros, Finite temperature field theory and phase transitions, in ICTP Summer School in High-Energy Physics and Cosmology, pp. 187–259, 1, 1999.  hep-ph/9901312.
  * [34] S. R. Coleman and E. J. Weinberg, Radiative Corrections as the Origin of Spontaneous Symmetry Breaking,  Phys. Rev. D 7 (1973) 1888–1910.
  * [35] P. Basler and M. Mühlleitner, BSMPT (Beyond the Standard Model Phase Transitions): A tool for the electroweak phase transition in extended Higgs sectors,  Comput. Phys. Commun. 237 (2019) 62–85, [1803.02846].
  * [36] P. Athron, C. Balázs, A. Fowlie, L. Morris, and L. Wu, Cosmological phase transitions: from perturbative particle physics to gravitational waves, 2305.02357.
  * [37] M. E. Carrington, The Effective potential at finite temperature in the Standard Model,  Phys. Rev. D 45 (1992) 2933–2944.
  * [38] P. Arnold and O. Espinosa, The Effective Potential and First-Order Phase Transitions: Beyond Leading Order,  Physical Review D 50 (Nov., 1994) 6662–6662, [hep-ph/9212235].
  * [39] D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, and G. White, Theoretical uncertainties for cosmological first-order phase transitions,  JHEP 04 (2021) 055, [2009.10080].
  * [40] P. Ginsparg, First and second order phase transitions in gauge theories at finite temperature,  Nuclear Physics B 170 (1980), no. 3 388–408.
  * [41] A. Ekstedt, P. Schicho, and T. V. I. Tenkanen, DRalgo: A package for effective field theory approach for thermal phase transitions,  Comput. Phys. Commun. 288 (2023) 108725, [2205.08815].
  * [42] S. R. Coleman, The Fate of the False Vacuum. 1. Semiclassical Theory,  Phys. Rev. D 15 (1977) 2929–2936. [Erratum: Phys.Rev.D 16, 1248 (1977)].
  * [43] A. D. Linde, Decay of the False Vacuum at Finite Temperature,  Nucl. Phys. B 216 (1983) 421. [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * [44] A. D. Linde, Fate of the False Vacuum at Finite Temperature: Theory and Applications,  Phys. Lett. B 100 (1981) 37–40.
  * [45] C. L. Wainwright, CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields,  Comput. Phys. Commun. 183 (2012) 2006–2013, [1109.4189].
  * [46] H.-K. Guo, K. Sinha, D. Vagie, and G. White, The Benefits of Diligence: How Precise are Predicted Gravitational Wave Spectra in Models with Phase Transitions?,  Journal of High Energy Physics 2021 (June, 2021) 164, [2103.06933].
  * [47] P. Athron, L. Morris, and Z. Xu, How robust are gravitational wave predictions from cosmological phase transitions?, 2309.05474.
  * [48] A. H. Guth and E. J. Weinberg, Cosmological consequences of a first-order phase transition in the s${\mathrm{u}}_{5}$ grand unified model,  Phys. Rev. D 23 (Feb, 1981) 876–885.
  * [49] J. Ellis, M. Lewicki, and J. M. No, On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal,  JCAP 04 (2019) 003, [1809.08242].
  * [50] F. Giese, T. Konstandin, and J. van de Vis, Model-independent energy budget of cosmological first-order phase transitions—A sound argument to go beyond the bag model,  JCAP 07 (2020), no. 07 057, [2004.06995].
  * [51] D. Bodeker and G. D. Moore, Electroweak Bubble Wall Speed Limit,  JCAP 05 (2017) 025, [1703.08215].
  * [52] Y. Gouttenoire, R. Jinno, and F. Sala, Friction pressure on relativistic bubble walls,  JHEP 05 (2022) 004, [2112.07686].
  * [53] A. Azatov, G. Barni, R. Petrossian-Byrne, and M. Vanvlasselaer, Quantisation Across Bubble Walls and Friction, 2310.06972.
  * [54] A. Roper Pol, S. Mandal, A. Brandenburg, T. Kahniashvili, and A. Kosowsky, Numerical simulations of gravitational waves from early-universe turbulence,  Phys. Rev. D 102 (2020), no. 8 083512, [1903.08585].
  * [55] K. Saikawa and S. Shirai, Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model,  JCAP 05 (2018) 035, [1803.01038].
  * [56] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Shape of the acoustic gravitational wave power spectrum from a first order phase transition,  Phys. Rev. D 96 (2017), no. 10 103520, [1704.05871]. [Erratum: Phys.Rev.D 101, 089902 (2020)].
  * [57] J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, Energy Budget of Cosmological First-order Phase Transitions,  JCAP 06 (2010) 028, [1004.4187].
  * [58] J. Ellis, M. Lewicki, and J. M. No, Gravitational waves from first-order cosmological phase transitions: lifetime of the sound wave source,  JCAP 07 (2020) 050, [2003.07360].
  * [59] M. Punturo et al., The Einstein Telescope: A third-generation gravitational wave observatory,  Class. Quant. Grav. 27 (2010) 194002.
  * [60] M. Maggiore et al., Science Case for the Einstein Telescope,  JCAP 03 (2020) 050, [1912.02622].
  * [61] J. Antoniadis et al., The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background,  Mon. Not. Roy. Astron. Soc. 510 (2022), no. 4 4873–4887, [2201.03980].
  * [62] S. S. AbdusSalam et al., Simple and statistically sound recommendations for analysing physical theories,  Rept. Prog. Phys. 85 (2022), no. 5 052201, [2012.09874].
  * [63] B. W. Lee and S. Weinberg, Cosmological Lower Bound on Heavy Neutrino Masses,  Phys. Rev. Lett. 39 (1977) 165–168.
  * [64] T. Bringmann, J. Edsjö, P. Gondolo, P. Ullio, and L. Bergström, DarkSUSY 6 : An Advanced Tool to Compute Dark Matter Properties Numerically,  JCAP 07 (2018) 033, [1802.03399].
  * [65] P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: Improved analysis,  Nucl. Phys. B 360 (1991) 145–179.
  * [66] T. Bringmann, P. F. Depta, M. Hufnagel, and K. Schmidt-Hoberg, Precise dark matter relic abundance in decoupled sectors,  Phys. Lett. B 817 (2021) 136341, [2007.03696].
  * [67] D. Chway, T. H. Jung, and C. S. Shin, Dark matter filtering-out effect during a first-order phase transition,  Phys. Rev. D 101 (2020), no. 9 095019, [1912.04238].
  * [68] Y. Gouttenoire and T. Volansky, Primordial Black Holes from Supercooled Phase Transitions, 2305.04942.
  * [69] D. Pappadopulo, J. T. Ruderman, and G. Trevisan, Dark matter freeze-out in a nonrelativistic sector,  Phys. Rev. D 94 (2016), no. 3 035005, [1602.04219].
  * [70] M. Farina, D. Pappadopulo, J. T. Ruderman, and G. Trevisan, Phases of Cannibal Dark Matter,  JHEP 12 (2016) 039, [1607.03108].
  * [71] R. M. Schabinger and J. D. Wells, A Minimal spontaneously broken hidden sector and its impact on Higgs boson physics at the large hadron collider,  Phys. Rev. D 72 (2005) 093007, [hep-ph/0509209].
  * [72] B. Patt and F. Wilczek, Higgs-field portal into hidden sectors, hep-ph/0605188.
  * [73] E. Weihs and J. Zurita, Dark Higgs Models at the 7 TeV LHC,  JHEP 02 (2012) 041, [1110.5909].
  * [74] M. Duerr, F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz, and S. Vogl, How to save the WIMP: global analysis of a dark matter model with two s-channel mediators,  JHEP 09 (2016) 042, [1606.07609].
  * [75] T. Bringmann, S. Heeba, F. Kahlhoefer, and K. Vangsnes, Freezing-in a hot bath: resonances, medium effects and phase transitions,  JHEP 02 (2022) 110, [2111.14871].
  * [76] E. Bernreuther, F. Kahlhoefer, M. Lucente, and A. Morandini, Searching for dark radiation at the LHC,  JHEP 08 (2022) 169, [2204.01759]. [Erratum: JHEP 06, 071 (2023)].
  * [77] J. A. Evans, S. Gori, and J. Shelton, Looking for the WIMP Next Door,  JHEP 02 (2018) 100, [1712.03974].
  * [78] A. Fradette, M. Pospelov, J. Pradler, and A. Ritz, Cosmological beam dump: constraints on dark scalars mixed with the Higgs boson,  Phys. Rev. D 99 (2019), no. 7 075004, [1812.07585].
  * [79] A. Berlin, D. Hooper, and G. Krnjaic, Thermal Dark Matter From A Highly Decoupled Sector,  Phys. Rev. D 94 (2016), no. 9 095019, [1609.02555].
  * [80] L. Husdal, On Effective Degrees of Freedom in the Early Universe,  Galaxies 4 (2016), no. 4 78, [1609.04979].
  * [81] G. Barenboim and W.-I. Park, Gravitational waves from first order phase transitions as a probe of an early matter domination era and its inverse problem,  Phys. Lett. B 759 (2016) 430–438, [1605.03781].
  * [82] M. Cirelli, Y. Gouttenoire, K. Petraki, and F. Sala, Homeopathic Dark Matter, or how diluted heavy substances produce high energy cosmic rays,  JCAP 02 (2019) 014, [1811.03608].
  * [83] NANOGrav Collaboration, G. Agazie et al., The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background,  Astrophys. J. Lett. 951 (2023), no. 1 L8, [2306.16213].
  * [84] EPTA Collaboration, J. Antoniadis et al., The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals,  Astron. Astrophys. 678 (2023) A50, [2306.16214].
  * [85] D. J. Reardon et al., Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,  Astrophys. J. Lett. 951 (2023), no. 1 L6, [2306.16215].
  * [86] N. F. Bell, Y. Cai, and R. K. Leane, Dark Forces in the Sky: Signals from Z’ and the Dark Higgs,  JCAP 08 (2016) 001, [1605.09382].
  * [87] CTA Collaboration, A. Acharyya et al., Sensitivity of the Cherenkov Telescope Array to a dark matter signal from the Galactic centre,  JCAP 01 (2021) 057, [2007.16129].
  * [88] B.-H. Liu, L. D. McLerran, and N. Turok, Bubble nucleation and growth at a baryon number producing electroweak phase transition,  Phys. Rev. D 46 (1992) 2668–2688.
  * [89] D. Bodeker and G. D. Moore, Can electroweak bubble walls run away?,  JCAP 05 (2009) 009, [0903.4099].
  * [90] A. Azatov and M. Vanvlasselaer, Bubble wall velocity: heavy physics effects,  JCAP 01 (2021) 058, [2010.02590].
  * [91] I. Garcia Garcia, G. Koszegi, and R. Petrossian-Byrne, Reflections on bubble walls,  JHEP 09 (2023) 013, [2212.10572].
  * [92] M. E. A. Vanvlasselaer, Dynamics of phase transitions in the early universe and cosmological consequences.  PhD thesis, SISSA, Trieste, 2022.
  * [93] S. Höche, J. Kozaczuk, A. J. Long, J. Turner, and Y. Wang, Towards an all-orders calculation of the electroweak bubble wall velocity,  JCAP 03 (2021) 009, [2007.10343].
  * [94] M. Lewicki and V. Vaskonen, Gravitational waves from bubble collisions and fluid motion in strongly supercooled phase transitions,  Eur. Phys. J. C 83 (2023), no. 2 109, [2208.11697].
  * [95] J. Ellis, M. Lewicki, and V. Vaskonen, Updated predictions for gravitational waves produced in a strongly supercooled phase transition,  JCAP 11 (2020) 020, [2007.15586].
  * [96] J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, Gravitational wave energy budget in strongly supercooled phase transitions,  JCAP 06 (2019) 024, [1903.09642].
  * [97] G. Arcadi, O. Lebedev, S. Pokorski, and T. Toma, Real Scalar Dark Matter: Relativistic Treatment,  Journal of High Energy Physics 2019 (Aug., 2019) 50, [1906.07659].
  * [98] V. De Romeri, D. Karamitros, O. Lebedev, and T. Toma, Neutrino dark matter and the Higgs portal: Improved freeze-in analysis,  Journal of High Energy Physics 2020 (Oct., 2020) 137, [2003.12606].