SISSA 12/2022/FISI

TU-1157

Ultra$-$relativistic bubbles from the simplest Higgs portal and their
cosmological consequences

Aleksandr Azatova,b,c,1, Giulio Barnia,b, Sabyasachi Chakrabortya,b,c, Miguel
Vanvlasselaera,b,c,2, Wen Yind,3,

a SISSA International School for Advanced Studies, Via Bonomea 265, 34136,
Trieste, Italy
b INFN - Sezione di Trieste, Via Bonomea 265, 34136, Trieste, Italy
c IFPU, Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy
d Department of Physics, Tohoku University, Sendai, Miyagi 980-8578, Japan

Abstract

> We analyze the phase transitions in the minimal extension of the SM with a
> real singlet scalar field. The novelty of our study is that we identify and
> analyze in details the region of parameter space where the first order phase
> transition can occur and in particular when the bubbles with true vacuum can
> reach relativistic velocities. This region is interesting since it can lead
> to the new recently discussed baryogenesis and Dark Matter production
> mechanisms. We fully analyze different models for the production of Dark
> Matter and baryogenesis as well as the possibilities of discovery at the
> current and future experiments.

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

The origin of Dark Matter (DM) and matter antimatter asymmetry are one of the
most important unresolved puzzles of the early Universe cosmology. In this
paper we will investigate in detail the novel mechanisms proposed in [1, 2, 3]
for DM production and baryon asymmetry generation during the first order phase
transitions (FOPT).

The FOPT are known to be very useful for constructing baryogenesis models
since the process is out-of-equilibrium as required by one of the Sakharov’s
conditions [4]. Applications of this mechanism to the electroweak phase
transition lead to the seminal scenarios of electroweak baryogenesis [5, 6].
In the Standard Model (SM), electroweak phase transition is expected to be a
smooth crossover [7, 8]. However, this is not the case for various Beyond
Standard Model (BSM) frameworks where successful scenarios of electroweak
baryogenesis can be built (see for a recent review [9]). The connection
between DM and phase transitions is less direct. However, there are class of
models where these two phenomenon are strongly related (see for example
Ref[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28]).

The focus of this paper will be electroweak FOPT with very relativistic
bubbles. In this context, collision between bubbles and particles in the
plasma, as was shown recently [29], can lead to the production of new states
with mass scale significantly larger than the scale of the phase
transition111Note that this effect is different from the mechanism presented
in[25, 30] where the heavy fields are produced by the collisions between
different bubble walls.. These heavy states can serve as DM candidates [2] or
their production and subsequent decay [1, 3], if accompanied with $C$, $CP$
and baryon number violating interactions can lead to baryon asymmetry
generation. The successful realization of this scenario require bubble
expansion with ultra-relativistic velocities. Such bubble wall motion is known
to be possible and is expected within certain classes of the potentials ([31,
32, 33, 34, 35]). However in the case of electroweak phase transition, the
requirement of ultra relativistic velocities leads to very non-trivial
constraints on the effective potential. In this paper we analyze in detail the
simplest extension of the SM that can lead to a FOPT: real singlet scalar
field. There have been numerous studies of phase transition in this type of
scenario [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. The Lagrangian can be
restricted further by considering $Z_{2}$ [47] or $Z_{3}$ symmetries[48, 49].
Interestingly such a real singlet scalar field can appear in composite Higgs
models [50] (see also [51, 52] for studies pertaining to phase transition).
However, most of the above mentioned works did not analyze in detail the
region of parameter space with relativistic bubbles since slow bubble wall
velocities are required [5, 6] for the usual electroweak baryogenesis. We fill
this gap by analyzing in detail phase transitions in the real singlet
extension of SM with an approximate $Z_{2}$ discrete symmetry, focusing on the
parameter space which results in very fast bubble expansions. Anticipating the
results of this paper we found that the most promising scenario for
relativistic bubbles is the case where the phase transition occurs in two
steps: the first transition is related to the $Z_{2}$ symmetry breaking while
keeping EW symmetry unbroken and in the second one the Higgs field acquires
its vacuum expectation value (vev) which comes together with $Z_{2}$ symmetry
restoration.

This will permit us to evaluate the region of parameter space where
Baryogenesis via relativistic bubble walls is viable and to put constraints on
the mass scale of the dark sector responsible for $CP$ and $B$ violation.
Implications for non-thermal DM production from bubble wall plasma particle
collisions will also be considered.

The article is organised as follows: in section 2, we review the singlet
extension of the SM and write the two field potential with thermal
corrections. In section 3, we remind the reader of the basics of FOPTs in
early Universe and present the computation of the terminal velocity of the
bubble wall in the ultra-relativistic regime, for the generic and for SM. In
section 4, we qualitatively discuss the different patterns of phase
transitions and argue that only a two-steps PT can yield ultra relativistic
bubble wall motion. In section 5 we present a numerical study of the two-step
PT analyzing in detail the region with fast bubble motion. In section 6, we
draw the consequences of our previous results for the production of heavy DM
and Baryogenesis and in section 7, we comment on the GW signal induced by such
strong transitions. Finally, in section 8, we summarize and conclude.

##  2 Review of the singlet extension of the SM

Let us start by reviewing the effective potential of the SM with a real scalar
field ($s$). The scalar potential will be given by

$\displaystyle
V(\mathcal{H},s)=-\frac{m_{h}^{2}}{2}(\mathcal{H}^{\dagger}\mathcal{H})+\lambda(\mathcal{H}^{\dagger}\mathcal{H})^{2}-\frac{m_{s}^{2}}{4}s^{2}+\frac{\lambda_{hs}}{2}s^{2}(\mathcal{H}^{\dagger}\mathcal{H})+\frac{\lambda_{s}}{4}s^{4}\;,$
(1)

where $\mathcal{H}$ is the SM Higgs doublet and $m_{h}\approx 125$ GeV is the
physical mass of the Higgs. For simplicity, we will impose $Z_{2}$ symmetry on
the potential to avoid the terms with odd powers of ${s}$ field. As a result,
when $\langle{s}\rangle=0$, we avoid any scalar mixing terms which are
constrained by the recent Higgs signal strength measurements [53]. The Higgs
doublet can be decomposed as usual

$\displaystyle\mathcal{H}^{T}=\bigg{(}G^{+},\frac{h+iG^{0}}{\sqrt{2}}\bigg{)}\;,$
(2)

where $h$ is the usual Higgs boson getting a vev which is given by
$v_{EW}=\sqrt{m_{h}^{2}/2\lambda}\approx 246$ GeV if $\langle{s}\rangle\equiv
v_{s}=0$.

###  2.1 Coleman-Weinberg potential

Next, we take into account 1-loop corrections which are encapsulated by the
Coleman-Weinberg (CW) potential [54]

$\displaystyle
V_{CW}=\sum_{i=Z,h,W,t}\frac{n_{i}(-1)^{F}}{64\pi^{2}}\left[M_{i}^{4}\bigg{(}\log\frac{M_{i}^{2}}{M^{2}_{i0}}-\frac{3}{2}\bigg{)}+2M^{2}_{i}M^{2}_{i0}\right]\;.$
(3)

In this expression, $M_{i}$ stand for the masses depending on the Higgs and
singlet fields values, $M_{i}\equiv M_{i}(h,s)$, and $M_{i0}$ are the field
values from the tree-level vev, $M_{i0}\equiv M_{i}(v_{EW},0)$. Then one can
easily check that this potential corresponds to the following renormalization
conditions

$\displaystyle V_{\rm eff}=V_{0}+V_{CW},$ $\displaystyle\frac{dV_{\rm
eff}}{dh}\bigg{|}_{(h,s)=(v_{EW},0)}=0,\qquad\frac{d^{2}V_{\rm
eff}}{dh^{2}}\bigg{|}_{(h,s)=(v_{EW},0)}=m^{2}_{h}\;.$ (4)

Eq.(3) cannot be used directly for contributions of the Goldstone bosons,
since their masses vanish in the true vacuum and the CW potential is IR
divergent. However the solution to this issue was suggested in Ref. [55] which
emphasized that the physical Higgs mass is defined at $p^{2}=m^{2}_{h}$ and
the effective potential at $p^{2}=0$ (see for details Ref. [55]). Thus it
would be better to use modified renormalization condition

$\displaystyle\frac{d^{2}V_{\rm
eff}}{dh^{2}}\bigg{|}_{(h,s)=(v_{EW},0)}=m_{h}^{2}-\Sigma\left(p^{2}=m_{h}^{2}\right)+\Sigma\left(0\right)\;,$
(5)

where the differences of self-energies are taking into account the running of
self-energy from $p^{2}=0$ to $p^{2}=m_{h}^{2}$. In this case the IR
divergences in $\Sigma(0)$ due to the virtual Goldstone bosons are cancelled
against the IR divergences of $V_{\rm eff}$. In practice, at the end of this
renormalization procedure, the contribution of the Goldstone bosons is given
by :

$\displaystyle
V^{GB}_{CW}(h)=\sum_{i=G}\frac{n_{G}}{64\pi^{2}}M_{G^{\pm,0}}^{4}(h)\bigg{[}\log\frac{M_{G}^{2}}{m_{h}^{2}}-\frac{3}{2}\bigg{]}\;.$
(6)

The number of d.o.f, masses of various particles and the scalar mass matrix as
a functions of $h,s$ are given by:

$\displaystyle\qquad\qquad\qquad\qquad n_{W^{\pm}}=6,\qquad n_{Z}=3,\qquad
n_{G^{\pm,0}}=3,\qquad n_{t}=-12\;,$ $\displaystyle
M^{2}_{W^{\pm}}(h)=\frac{g^{2}h^{2}}{4},\quad
M^{2}_{Z}(h)=\frac{g^{2}+g^{\prime 2}}{4}h^{2},\quad
M^{2}_{G^{\pm,0}}(h)=-\frac{m_{h}^{2}}{2}+\lambda h^{2},\quad
M^{2}_{t}(h)=\frac{y^{2}_{t}h^{2}}{2}$ $\displaystyle
M^{2}(h,s)=\begin{pmatrix}3\lambda
h^{2}+\frac{\lambda_{hs}}{2}s^{2}-\frac{m_{h}^{2}}{2}&\lambda_{hs}sh\\\
\lambda_{hs}sh&3\lambda_{s}s^{2}+\frac{\lambda_{hs}}{2}h^{2}-\frac{m_{s}^{2}}{2}\;\end{pmatrix}\;.$
(7)

As a side remark, in the region where the Higgs $h\to 0$ and $s\sim
O(v_{EW})$, there will be two scales involved in the problem: the value of the
singlet field $s$ and the masses of the SM particles $M_{W,Z,t}(h\to 0,s)\to
0\ll s$. This type of two scale potential has been studied in the past [56],
by using two different renormalisation scales. It was concluded that
resummation is needed when the $\log\left(M_{i}(0,s)/s\right)$ is large enough
to cancel the loop suppression. Although we have two largely separated scales,
we have checked that for our region of the parameter space, such a resummation
is not necessary.

###  2.2 Finite temperature potential

The temperature and the density effects can be taken into account by
complementing the zero temperature potential with thermal corrections (see for
example [57, 58]),

$\displaystyle V(T,M_{i})=V^{T=0}_{\rm eff}(M_{i})+V_{T}(M_{i})\;.$ (8)

In Eq.(8), $V^{T=0}_{\rm eff}(M_{i})$ is the potential we computed in the
previous subsection and the thermal potential $V_{T}(M_{i})$ is given by

$\displaystyle V_{T}(M_{i}(h,s))$ $\displaystyle=\sum_{i\in
B}\frac{n_{i}}{2\pi^{2}}T^{4}J_{B}\bigg{(}\frac{M_{i}^{2}(h,s)}{T^{2}}\bigg{)}-\sum_{i\in
F}\frac{n_{i}}{2\pi^{2}}T^{4}J_{F}\bigg{(}\frac{M_{i}^{2}(h,s)}{T^{2}}\bigg{)}\;,$
$\displaystyle J_{B/F}(y^{2})$ $\displaystyle=\int\limits_{0}^{\infty}dx\
x^{2}\log\Big{[}1\mp\exp{(-\sqrt{x^{2}+y^{2}})}\Big{]}\;.$ (9)

However, to save computation time, we use the approximate expansion of the
function $J_{B/F}(y^{2})$ as given in Ref. [57]:

$\displaystyle J_{B}(y^{2})$
$\displaystyle=\begin{cases}-\dfrac{\pi^{4}}{45}+\dfrac{\pi^{2}}{12}y^{2}-\dfrac{\pi}{6}y^{3}-\dfrac{y^{4}}{32}\log\bigg{[}\dfrac{y^{2}}{16\pi^{2}\exp[3/2-2\gamma]}\bigg{]}\;,&y^{2}\ll
1\;,\\\ -{\sum}_{n=1}^{m>3}\dfrac{1}{n^{2}}y^{2}K_{2}(yn)\;,&y^{2}\gg
1\;,\end{cases}$ $\displaystyle J_{F}(y^{2})$
$\displaystyle=\begin{cases}\dfrac{7\pi^{4}}{360}-\dfrac{\pi^{2}}{24}y^{2}-\dfrac{y^{4}}{32}\log\bigg{[}\dfrac{y^{2}}{\pi^{2}\exp[3/2-2\gamma]}\bigg{]}\;,&y^{2}\ll
1\;,\\\ -{\sum}_{n=1}^{m>3}\dfrac{(-1)^{n}}{n^{2}}y^{2}K_{2}(yn)\;,&y^{2}\gg
1\;.\end{cases}$ (10)

where $\gamma\approx 0.577$ is the _Euler constant_ , and $K_{2}(z)$ are the
second-kind Bessel function. To account for dangerously divergent higher loops
due to the Daisy diagrams at finite temperature, we follow the so-called
“Truncated-Full-Dressing” procedure [57]. Doing so, the full one-loop
potential becomes

$\displaystyle V(h,s,T)=V_{\rm tree}(h,s)+$
$\displaystyle\sum_{i}\bigg{[}V_{CW}\bigg{(}M_{i}^{2}(h,s)+\Pi_{i}(h,T)\bigg{)}+V_{T}\bigg{(}M_{i}^{2}(h,s)+\Pi_{i}(h,T)\bigg{)}\bigg{]},$
(11)

where $\Pi_{i}(T)$ are the thermal masses of various degrees of freedom. In
the real singlet extension of the SM[57], the expressions of the thermal
masses read

$\displaystyle\text{Scalar:}\qquad\Pi_{h}(T)=T^{2}\bigg{(}\frac{3g^{2}}{16}+\frac{g^{\prime
2}}{16}+\frac{\lambda}{2}+\frac{y_{t}^{2}}{4}+\frac{\lambda_{hs}}{24}\bigg{)}\;,\qquad\Pi_{s}(T)=T^{2}\bigg{(}\frac{\lambda_{hs}}{6}+\frac{\lambda_{s}}{4}\bigg{)}\;,$
(12)
$\displaystyle\text{Gauge:}\qquad\Pi^{L}_{g}(T)=T^{2}\text{diag}\bigg{(}\frac{11}{6}g^{2},\frac{11}{6}(g^{2}+g^{\prime
2})\bigg{)}\;,\qquad\Pi^{T}_{g}(T)=0\;,$ (13)

where $\Pi_{g}^{L}(T)$ denote the thermal mass of the _longitudinal_ mode of
the gauge bosons, while transverse modes $\Pi^{T}_{g}(T)$ are protected by
gauge invariance and thus do not receive a mass at leading order in
perturbation theory.

##  3 First order phase transitions with relativistic bubbles

Starting with the effective potential derived in the previous section, we can
proceed to the analysis of the phase transition. FOPT happens when the minima
of the potential corresponding to different phases are separated by a
potential barrier and the transition happens via bubble nucleation. The
probability of nucleation of a bubble is given by [59, 60, 61]

$\displaystyle\Gamma(T)\simeq\Gamma_{3}+\Gamma_{4}=T^{4}\bigg{(}\frac{S_{3}}{2\pi
T}\bigg{)}^{3/2}e^{-S_{3}(T)/T}+\frac{1}{R_{0}^{4}}\left(\frac{S_{4}}{2\pi}\right)e^{-S_{4}}\;,$
(14)

where $S_{3},S_{4}$ are $O(3,4)$ bounce actions and $R_{0}$ is the initial
bubble radius. Bubble nucleation is characterised by a critical temperature
$T_{\rm crit}$, defined as the point when the two phases of the system have
vacua with the same energy. Below $T_{\rm crit}$ phase transition becomes
energetically possible. The probability to find a specific point of the
Universe to be in the false vacuum is given by [34, 62, 63]:

$\displaystyle P_{f}[T]=\exp[-I(T)]\;,\qquad
I(T)\equiv\frac{4\pi}{3}\int^{T_{c}}_{T}\frac{dT_{1}\Gamma(T_{1})v_{w}^{3}}{T_{1}^{4}H(T_{1})}\bigg{[}\int_{T}^{T_{1}}\frac{dT_{2}}{H(T_{2})}\bigg{]}^{3}.$
(15)

In Eq.(14) the strongest dependence on the temperature comes from the
$\Gamma(T)\propto\exp{(-S_{3}/T)}$, so that the quantity $I(T)$ is mostly
controlled by the ratio $\Gamma\left(T\right)/H^{4}\left(T\right)$ and an
order one fraction of the volume of the Universe will be in the true vacuum
when $\Gamma[T]\sim H^{4}[T]$. The temperature that satisfies this condition
is coined as the nucleation temperature $T_{\rm nuc}$ or, to phrase
differently, the nucleation temperature $T_{\rm nuc}$ corresponds to the
appearance of roughly one bubble of true vacuum per Hubble volume. Assuming
the scale of the phase transition to be $O(100)$ GeV we get

$\displaystyle\Gamma\sim H^{4},~{}~{}~{}H(T)^{2}=\frac{1}{3M^{2}_{\rm
pl}}\big{(}\rho_{\rm rad}+\rho_{\rm vac}+\rho_{\rm
wall}\big{)},\quad\rho_{w}\approx 0,\quad\rho_{\rm
rad}=\frac{30g_{\star}}{\pi^{2}}T^{4},$ (16) $\displaystyle\quad\rho_{\rm
vac}=\Delta V,\quad\Rightarrow\quad\frac{S_{3}}{T_{\rm nuc}}={3\over
2}\log\left(\frac{S_{3}}{2\pi T_{\rm nuc}}\right)+4\log\left(\frac{T_{\rm
nuc}}{H}\right)\;.$ (17)

In case we want to be more precise about the temperature when the phase
transition completes, we can define the so-called _percolation temperature_
$T_{\rm per}$ as the temperature when around $\sim$30 % of the space has been
converted to the true phase

$\displaystyle I(T\equiv T_{\rm per})=0.34\;,\qquad\text{(percolation
temperature)}.$ (18)

If the condition in Eq.(18) is not fulfilled, $I(T)<0.34$, then the bubbles do
not percolate.

At last the energy released during the PT is traditionally quantified by the
ratio between the energy stored in the vacuum and in the plasma at the moment
of the transition

$\displaystyle\alpha\equiv\frac{\Delta V}{\rho_{\rm rad}}\qquad\text{(strength
parameter)}.$ (19)

###  3.1 Velocity of the EW bubble

Let us proceed to the computation of the bubble wall velocity $v_{w}$. The
dynamics of the bubble wall motion is controlled by the driving force due to
the potential differences between the false and true vacuum

$\displaystyle\Delta V\equiv V_{\rm false}-V_{\rm true}\;,$ (20)

and the friction due to finite temperature effects. The calculation of the
friction is generically a very complicated problem, however for the ultra-
relativistic bubble motion at leading order (LO - tree level) very simple
expressions have been obtained[64, 65] for the pressure force from friction,

$\displaystyle\Delta\mathcal{P}_{\text{LO}}\simeq\sum_{i}g_{i}c_{i}\frac{\Delta
M_{i}^{2}}{24}T_{\text{nuc}}^{2}\;,$ (21)

where $\Delta M_{i}^{2}$ is the change in the mass of the particle $i$ during
the PT, $c_{i}=1\left(1/2\right)$ for bosons (fermions) and $g_{i}$ is the
number of d.o.f of the incoming particle. Eq.(21) assumes that the masses of
the particles outside of the bubble are less than the temperature, otherwise
the friction will have an additional Boltzmann suppression
$\propto\exp[-M_{\rm false}/T]$. Interestingly the production of the heavy
particles can also contribute to the friction at the same order [29]

$\displaystyle\Delta\mathcal{P}^{\rm mixing}_{\text{LO}}\propto
v^{2}T^{2}\Theta(\gamma_{w}T_{\rm nuc}-M^{2}_{\rm heavy}/v)\;,$ (22)

where $v$ is the vev of the Higgs field and we review the heavy particle
production later in the section 6.1. One can see that the friction (pressure
from plasma on the bubble wall) becomes velocity independent so that permanent
accelerating (runaway) behavior of the bubble expansion becomes possible.
However for theories where the gauge bosons receive a mass during the phase
transition, this is not the case and the effect with multiple gauge boson
emissions, leads to the additional contribution (NLO- Next to Leading order)
to the pressure which scales as [66, 29, 67, 20]222There is a claim that this
friction scales as $\gamma^{2}T^{4}$[68] see [29, 67] for criticism.

$\displaystyle\Delta\mathcal{P}_{\text{NLO}}$ $\displaystyle\propto$
$\displaystyle\sum_{a}g_{a}g_{gauge}^{3}\gamma_{w}T_{\text{nuc}}^{3}v\;.$ (23)

At this point we can see that bubbles will keep accelerating till the moment
when NLO friction becomes large enough to balance the driving force

$\displaystyle\Delta V=\Delta{\cal P}_{\rm LO}+\Delta{\cal P}^{\rm
mixing}_{\rm LO}+\Delta{\cal P}_{\rm NLO},$ (24)

which will set the $\gamma_{w}^{\rm terminal}$ the bubble can reach,

$\displaystyle\gamma_{w}^{\text{terminal}}\sim\left(\frac{\Delta V-\Delta{\cal
P}_{\rm LO}-\Delta{\cal P}^{\rm mixing}_{\rm LO}}{T_{\rm nuc}^{3}v}\right).$
(25)

Note that we can be as well in the situation where the percolation (bubble
collision) starts before the terminal $\gamma_{w}^{\rm terminal}$ is reached.

###  3.2 Friction forces during the electroweak phase transitions

Let us apply the discussion of the previous section to the EW phase transition
(EWPT). Our main interest will be the possibility of heavy particle production
which can later source baryogenesis and DM production following the ideas in
[29, 2, 1, 3]. The heavy particle production (see also the discussion in
section 6.1) happens due to the collision between the plasma particles and the
bubble wall. The typical center of mass energy for such process will be
roughly $\sim\sqrt{\gamma_{w}Tv_{EW}}$, thus $\gamma_{w}$ will be controlling
the maximal mass of the new fields that can be produced. To calculate
$\gamma_{w}$ in the context of the EWPT, we need to know the forces acting on
the bubble wall.

The LO friction from the top, $Z$, and $W$ takes the form:

$\displaystyle\Delta\mathcal{P}^{SM}_{\rm LO}\approx
T_{\text{nuc}}^{2}v_{EW}^{2}\bigg{(}\frac{y_{t}^{2}}{8}+\frac{g^{2}+g^{\prime
2}}{32}+\frac{g^{2}}{16}\bigg{)}\approx 0.17\;T_{\text{nuc}}^{2}v_{EW}^{2}.$
(26)

There is also a contribution to the LO friction from the singlet and Higgs
scalars, however it depends on the masses of these fields in the false vacuum
and we find it to be numerically subleading compared to the estimate in
Eq.(26)333These contributions are smaller due to the number of d.o.f. and
possible Boltzmann suppression factors $\propto\exp[-M_{\rm false}/T].$. In
our numerical calculation we took this additional contribution into account,
however, for the current discussion it is sufficient to use Eq.(26). This
gives a rough condition on the nucleation temperature for the transition to
become ultra-relativistic

$\displaystyle\Delta
V>0.17\;T_{\text{nuc}}^{2}v_{EW}^{2}\qquad\text{(relativistic wall
condition)}.$ (27)

The computation of the NLO friction in the SM has been carried out in [67] and
the following approximate expression has been derived:

$\displaystyle\Delta\mathcal{P}^{SM}_{\rm
NLO}\approx\bigg{[}\sum_{abc}\nu_{a}g_{a}\beta_{c}C_{abc}\bigg{]}\frac{\kappa\zeta(3)}{\pi^{3}}\alpha
M_{Z}(v_{EW})\gamma_{w}\log\frac{M_{Z}(v_{EW})}{\mu}T_{\text{nuc}}^{3}\;,$
(28)

where $\nu_{a}=1(3/4)$ for $a$ a boson (fermion), $g_{a}$ is the number of
degrees of freedom of $a$ and $M_{g,i}(v_{EW})$ is the mass of the gauge boson
inside the bubble. $C_{abc}$ represents the coupling which appears in the
vertex (see appendix C), $\beta_{c=Z^{0}}=1$, and
$\beta_{c=W^{\pm}}=\cos\theta_{W}=M_{W}/M_{Z}$, $\alpha=e^{2}(v_{EW})/4\pi\sim
1/128$ is the electromagnetic fine structure constant. The renormalization
scale $\mu$ has to be understood as the lower cut-off for the integration over
soft momentum, typically the thermal mass
$\mu\propto\alpha_{i}^{1/2}T_{\text{nuc}}$[29]. The $\kappa$ factor is
introduced [67] to account for the contributions of the both reflected and
transmitted bosons and is approximately equal to $\kappa\sim 4$. The sum in
the square brackets is approximately equal to

$\displaystyle\bigg{[}\sum_{abc}\nu_{a}g_{a}\beta_{c}C_{abc}\bigg{]}\approx
157\;,$ (29)

see appendix C for details. At this point we can compute the terminal wall
velocity by balancing the pressure against the driving force

$\displaystyle\Delta V-\Delta\mathcal{P}^{SM}_{\rm
LO}=\Delta\mathcal{P}^{SM}_{\rm NLO}(\gamma_{w}=\gamma^{\rm
terminal}_{w})\quad\text{(Terminal velocity criterion)}$ (30)
$\displaystyle\Rightarrow\gamma^{\rm terminal}\sim 6\times\left(\frac{\Delta
V-\Delta\mathcal{P}_{\rm
LO}^{SM}}{(100{\rm~{}GeV})^{4}}\right)\left(\frac{100{\rm GeV}}{T_{\rm
nuc}}\right)^{3}\frac{1}{\log\frac{M_{z}}{gT}}.$ (31)

Taking into account that $\left(\Delta V-\Delta{\cal P}_{\rm
LO}\right)/(100{\rm~{}GeV})^{4}\lesssim\mathcal{O}(1)$, we can see that the
bubbles will become ultra-relativistic, i.e. $\gamma^{\rm terminal}\gg 1$ only
if $T_{\rm nuc}$ is significantly lower than the scale of the phase transition
$\sim 100$ GeV.

##  4 Phase transition in the singlet extension

After the preparatory discussion in the previous sections 2 and 3, we can
proceed to the analysis of phase transition in the model with the singlet
field, Eq.(1). Our analysis will be focused on the region of parameter space
with relativistic bubble expansion, for other studies of phase transition in
SM plus $Z_{2}$ real singlet scalar, see Refs. [69, 47, 70, 71].

In the previous section 3, we have seen that the velocity of the bubble
expansion is fixed by the balance between the friction from the plasma and the
driving force. At low temperatures the friction is suppressed (see
Eq.(26)-(28)) so that we expect the bubbles to become relativistic (large
Lorentz $\gamma^{\rm terminal}_{w}$ factor).

Let us check whether low nucleation temperatures are feasible in the singlet
extension. We will assume that $Z_{2}$ remains unbroken in the true vacuum in
order to avoid to constraints from the Higgs-scalar mixing (see for example
[72]). Then in the model with $Z_{2}$ odd singlet, the phase transition can
occur in two ways: one-step
$(\langle{h}\rangle=0,\langle{s}\rangle=0)\to(v_{EW},0)$ and two-steps444We
will see later that at temperatures $T\ll T_{c}$ Coleman-Weinberg potential
can shift a little bit the false vacuum position to $(\delta v_{h},v_{s})$
where $\delta v_{h}\ll v_{s},v_{h}$.
$(\langle{h}\rangle=0,\langle{s}\rangle=0)\to(0,\langle{s}\rangle|_{\neq
0})\to(v_{EW},0)$ [44, 70, 47], and each of these phase transitions can be
first or the second order. We review both of these scenarios of phase
transitions in order to understand in what case it is possible to obtain
relativistic bubbles.

###  4.1  One-step phase transition:

This case has been largely studied in the literature [70, 47] and we will not
provide a complete description of it. In this scenario the singlet never gets
a vev, and all of its effect reduce to the additional contributions to the
Higgs potential from Coleman-Weinberg terms and thermal corrections. However
it turns out that relativistic bubbles are very unlikely for such phase
transitions (see also results in [47]). In the limit when $|m_{s}|\ll T_{c}$
we can show analytically that this is indeed the case. Near the origin $h\to
0$, the potential in $h$ direction is dominated by the $\propto h^{2}$ terms

$\displaystyle V_{\text{dominant}}(h\to 0,s=0,T)$ $\displaystyle=$
$\displaystyle\frac{m_{\rm eff}^{2}(T)}{2}h^{2}+\dots$ (32)

The effective mass $m_{\rm eff}^{2}$ include the tree-level terms and the
leading thermal contributions and is approximately equal to

$\displaystyle m_{\rm eff}^{2}(T)\simeq-
m_{h}^{2}+T^{2}\bigg{(}\frac{m_{h}^{2}}{4v_{EW}^{2}}+\frac{y_{t}^{2}}{4}+\frac{g^{2}+g^{\prime
2}}{16}+\frac{g^{2}}{8}+\frac{\lambda_{hs}}{24}\bigg{)}>0\;.$ (33)

The FOPT can happen only if $m_{\rm eff}^{2}(T)>0$. Then assuming perturbative
values of the coupling $\lambda_{hs}$ we can estimate the lowest temperature
where the FOPT can occur to be $T^{nuc}_{\text{min}}\gtrsim 100$ GeV.
Comparing this value with the discussion in the section 3.2 we can see that
the bubble wall velocities will always satisfy $\gamma_{w}\lesssim 10$. As
mentioned before, we are interested in the expansions with much larger
$\gamma_{w}$ factors, so that we do not discuss one step phase transition
further.

###  4.2 Two step FOPT with relativistic bubbles

Two-steps realisations of the EWPT have already been studied in many works,
see for example[73, 44, 74, 70, 75, 76, 47, 77]. The novelty of our study is
that we will be focusing on the parameter space with relativistic bubbles
which was previously ignored. We organize the discussion as follows: In
section 4.2, we show qualitative results based on approximate treatment of the
potential and then in the section 5 we present the exact numerical results
obtained with our code. The two step phase transition

$\displaystyle(0,0)\to(0,\langle{s}\rangle)\to(v_{EW},0)\;,$ (34)

can happen if the $m_{s}^{2}$ parameter of Eq.(1) is positive. In this case it
is convenient to parameterize the Lagrangian in the following way

$\displaystyle V_{\rm
tree}(h,s)=-\frac{m_{h}^{2}}{4}h^{2}+\frac{m_{h}^{2}}{8v_{EW}^{2}}h^{4}-\frac{m_{s}^{2}}{4}s^{2}+\frac{\lambda_{hs}}{4}s^{2}h^{2}+\frac{m_{s}^{2}}{8v_{s}^{2}}s^{4}\;.$
(35)

where $v_{EW}=\sqrt{m_{h}^{2}/2\lambda}$ GeV and
$v_{s}=\sqrt{m_{s}^{2}/2\lambda_{s}}$ correspond to the local minima at
$(\langle{h}\rangle=v_{EW},{\langle{s}\rangle=0)}$ and
$(\langle{h}\rangle=0,{\langle{s}\rangle=v_{s})}$ respectively. The origin of
two-step PT can be intuitively understood from the following considerations.
For simplicity let us ignore the Coleman-Weinberg potential and restrict the
discussion by considering only the thermal masses. Then the potential will be
given by

$\displaystyle V({{\cal H},s})$ $\displaystyle\approx$ $\displaystyle V_{\rm
tree}({\cal H},s)+\frac{T^{2}}{24}\left[\sum_{bosons}n_{i}M_{i}^{2}({\cal
H},s))+\frac{1}{2}\sum_{fermions}n_{F}M_{F}^{2}({\cal H},s)\right]\;,$
$\displaystyle=$ $\displaystyle V_{\rm tree}({\cal H},s)$ $\displaystyle+$
$\displaystyle T^{2}\left[h^{2}\left(\frac{g^{\prime
2}}{32}+\frac{3g^{2}}{32}+\frac{m_{h}^{2}}{8v_{EW}^{2}}+\frac{y_{t}^{2}}{8}+\frac{\lambda_{hs}}{48}\right)+s^{2}\left(\frac{m_{s}^{2}}{16v_{s}^{2}}+\frac{\lambda_{hs}}{12}\right)\right].$

From this expression we can clearly see that the temperatures when the minima
with non zero vevs appear for the Higgs and singlet fields can be different.
Then it can happen that the $Z_{2}$ breaking phase transition occurs before
the EW one. This means that there will be first a phase transition from
$(0,0)\to(0,v_{s})$. After this phase transition the Universe keeps cooling
down and the minimum with $\langle h\rangle\neq 0$ will be generated. Choosing
the appropriate values of masses and couplings we can make sure that the
minimum with $(v_{EW},0)$ is the true minimum of the system. The transition
$(0,v_{s})\rightarrow(v_{EW},0)$ will be of the first order if there is a
potential barrier in the between two minima. One of the necessary condition in
this case will be $\partial^{2}V/\partial h^{2}|_{s=v_{s},h\to 0}>0$, which
using Eq.(4.2) we get

$\displaystyle\frac{-m_{h}^{2}}{4}+\frac{\lambda_{hs}v_{s}^{2}}{4}+T^{2}\left(\frac{g^{\prime
2}}{32}+\frac{3g^{2}}{32}+\frac{m_{h}^{2}}{8v^{2}}+\frac{y_{t}^{2}}{8}+\frac{\lambda_{hs}}{48}\right)>0\;.$
(37)

From this discussion we can see that there are qualitatively two different
cases depending on whether the potential barrier between two minima remains or
disappears at zero temperature, i.e. when
$m_{h}^{2}\gtrless\lambda_{hs}v_{s}^{2}$. In the first case the phase
transition will necessarily happen before the “No Barrier”(NB) temperature

$\displaystyle T^{\rm
NB}=\sqrt{\frac{m_{h}^{2}-\lambda_{hs}v_{s}^{2}}{\frac{g^{\prime
2}}{8}+\frac{3g^{2}}{8}+\frac{m_{h}^{2}}{2v^{2}}+\frac{y_{t}^{2}}{2}+\frac{\lambda_{hs}}{12}}}\sim\frac{v_{s}\sqrt{-\lambda_{hs}+\big{(}\frac{m_{h}}{v_{s}}\big{)}^{2}}}{0.8}\;,$
(38)

since the bounce action drops once the potential barrier between two minima
reduces.

In the other case, when the barrier remains even for zero temperature
potential, the phase transition is obviously of the first order. However it
might happen that the tunneling rate is too slow and the system remains stuck
in the false vacuum.

From Eq.(37) we can see that size of the potential barrier between two minima
will be controlled by the coupling $\lambda_{hs}$. Increasing this parameter
will enhance the potential barrier and will lead to the reduction of $T^{\rm
NB}$ so that $(v_{s},0)$ remains a local minimum even at zero temperature. At
some point we expect the potential barrier to become so large that the system
will remain trapped in the false vacuum forever. From this discussion we can
expect that the lowest temperatures will be achieved at the boundary of the
region where no PT can occur. The lowest temperature will correspond to the
(global) minimum of $S_{3}/T$ controlling the tunneling rate, as explained
later in section 5. If this minimum happens at temperatures much lower than
$T_{\rm crit}$, then there will be super-cooling ($T_{\rm nuc}\ll T_{\rm
crit}$). Even though this discussion was made by considering only the thermal
masses in the potential, numerically we find that for the effective potential
with truncated full dressing [57] the qualitative behaviour does not change,
and only explicit conditions for $m_{h}^{2}$ and $T^{\rm nuc}_{\rm min}$ are
modified.

At last we always check whether the condition for the EWSB minimum to be the
global minimum at zero temperature is satisfied:

$\displaystyle m_{s}^{2}v_{s}^{2}<m_{h}^{2}v_{EW}^{2}\qquad\text{(EWSB is
global minimum)}\;,$ (39)

where above equation is valid up to small loop-level corrections.

##  5 Numerical results

After the qualitative discussion, let us proceed to the numerical
calculations. The bounce action for $(0,v_{s})\to(v_{EW},0)$ was computed
using our own dedicated code (cross checked against FindBounce[78]) and we
relegate the details and methods of this calculation to Appendix A. We
parameterize our model in terms of ($m_{s},\lambda_{hs},v_{s}$) parameters
(see Eq.(35)). Instead of analyzing all the possible values of $m_{s}$, in
this section we report the results (see Fig.1 left panel) by fixing
$m_{s}=125$ GeV (there is no particular reason for this value of $m_{s}$ and
the results for the other values of $m_{s}$ are similar and are shown in the
appendix B). In the plane (see Fig.1) $v_{s}-\lambda_{hs}$, we identify four
regions with different behaviours under the phase transition. The blue region
shows second order transition when there is no barrier between the two
separated vacuum. Next to it there is a light red region where the transition
is of the first order. We indicate separately (dark red) the region, where the
transition is of the first order and the bubbles are relativistic. In
particular the boundary between the regions with relativistic and non-
relativistic bubbles is defined by the criteria of Eq.(27), i.e. when the LO
pressure for relativistically expanding bubbles is less than the driving
force. At last, there is NO PT region, where the system remains stuck in the
false vacuum since the tunnelling rate is too small.

The structure of the diagram (on the Fig.1 top$-$left panel) can be easily
understood from the qualitative discussion in the previous section. Indeed,
keeping $v_{s}$ fixed, the size the potential barrier is controlled by the
coupling $\lambda_{hs}$. Moving from left to right the size of the potential
barrier increases and we are gradually moving from the region of second order
phase transition $\to$ FOPT $\to$ FOPT with relativistic bubbles $\to$ no PT
region. Similarly moving up (increasing $v_{s}$ for fixed values of
$\lambda_{hs}$) also corresponds to the increase of the potential barrier as
we pass though the regions with different PT in the same order. On the upper
axis we report physical mass of the singlet in the true vacuum
$M_{s}(v_{EW},0)$ and exclude the constrained region where $h\to ss$ is
kinematically allowed (the gray meshed region).

Next we plot the values of $T_{\rm nuc}$ and $\gamma^{\rm terminal}_{w}$ as a
function of $(\lambda_{hs},v_{s})$ (Fig.1). As discussed in the previous
section, the region with the smallest values of the nucleation temperature
(thus the fastest bubbles) is located near the “NO PT” region, i.e. where the
system remains trapped in the false vacuum. The blue dot on the $T_{\rm nuc}$
plot indicates the last point we have found before the system enters the
regime of no phase transition (NO PT). In the bottom right we plot $M_{\rm
max}=\sqrt{\gamma_{w}T_{\rm nuc}v_{EW}}$ quantity which indicates the maximal
energy in the c.o.m frame for the plasma particle$-$bubble wall collision,
which corresponds to the largest mass of the heavy particles we can produce
(see discussion in sec 6.1).

Figure 1:  Top-Left: Scan of the parameter space in the plane
$\lambda_{hs}-v_{s}$ for $m_{s}=125$ GeV. The four regions are as follows: i)
white - NO PT, the region where the transition never completes because the
barrier remains at zero temperature and the function $S_{3}/T$ never passes
below the nucleation condition, ii) light and dark red are the regions where
the FOPT happens. Dark (light) red corresponds to the region with relativistic
(non-relativistic) bubble expansion. The boundary between two regions is given
by Eq.(27) iv) blue - the phase transition is of the second order. The gray
meshed region is the one in which $M_{s}(v_{EW},0)<m_{h}/2$, that is
constrained from collider experiments. Top-Right and Bottom: the dependencies
of $T_{\rm nuc},\gamma_{w},\ M_{max}$ on the $\lambda_{hs}$ coupling for
$m_{s}=125$ GeV. The blue dot on the Top-Right plot designates the end of the
curve, when the tuning becomes $10^{-6}$ and the red dot signals the
appearance of a barrier at zero temperature (all the points above the red dot
have a barrier at $T=0$). For the last three plots we varied the
$\lambda_{hs}$ parameter with the steps of $10^{-6}$. The value of $v_{s}$ is
encoded in color according to the bottom-right plot.

In order to better understand the dependencies of $T_{\rm nuc}$ and
$\gamma_{w}$ on the parameters of the model, it is useful to separate further
the parameter space depending on whether the potential barrier disappears at
zero temperature or not. This is important since the bounce action in such
cases has very different dependencies on the temperature.

Figure 2: Left: $S_{3}/T$ function with $v_{s}=170$ GeV and $m_{s}=125$ GeV.
We observe that the nucleation temperature saturates around $T_{\rm
nuc}\approx 52$ GeV. The horizontal gray line satisfy the nucleation condition
$S_{3}/T_{\rm nuc}={3\over 2}\log\left(S_{3}/2\pi T_{\rm
nuc}\right)+4\log\left(T_{\rm nuc}/H\right)$. Right: Zooming in on the region
with lower nucleation temperatures.

See for example Fig.2, where we have plotted the $S_{3}/T$ for
$v_{s}=170,m_{s}=125$ GeV for the various couplings $\lambda_{hs}$. In the
case when the potential barrier disappears at some temperature $T^{\rm NB}$
the function $S_{3}/T$ drops to zero, but if the barrier remains even at zero
temperature $S_{3}/T$ has a global minimum for $T\neq 0$ which will be
controlling the lowest nucleation temperature possible.

###  5.1 No potential barrier at zero temperature

Let us start by defining the region where there is no potential barrier at
zero temperature. On Fig.1, we demonstrate the curves where the barrier
disappears for various values of the temperatures ($T^{\rm NB}$). For $T^{\rm
NB}=0$ case, approximate curve can be obtained analytically by looking at the
leading terms in the zero temperature CW potential

$\displaystyle\lambda_{hs}\lesssim\frac{m_{h}^{2}}{v_{s}^{2}}-\frac{n_{t}y_{t}^{4}}{32\pi^{2}}\frac{v_{EW}^{2}}{v_{s}^{2}}\;.$
(40)

The agreement between this equation and exact $T^{\rm NB}=0$ curve is at the
level of a few permille discrepancies. To the right of $T^{\rm NB}=0$ curve,
the potential barrier between the two minima remains even at zero
temperatures. For the values $v_{s}\gtrsim 200$ GeV, we find that the line
$T^{\rm NB}=0$ approximately coincides with the boundary of no phase
transition region (where system remains stuck in the false vacuum) but
obviously the boundary of “NO PT” is always to the right of $T^{\rm NB}=0$
curve. The size of this narrow strip is of the order $10^{-4}$ in
$\lambda_{hs}$ values. One can see it from the $T_{\rm nuc}$ panel of Fig.1
where we have indicated the value of $\lambda_{hs}$ when $T_{\rm NB}=0$ by
vertical thin line and red dot (for intersection) and the position of the blue
dot which is the last point where the transition is of the first order before
we enter NO PT region. The boundaries of this region were obtained by
numerical calculations where we have scanned $\lambda_{hs}$ parameter with a
step $10^{-6}$. We postpone the discussion of the FOPT in this narrow region
to the next section 5.2.

In this section we restrict our discussion only on the region to the left of
$T^{\rm NB}=0$ curve. Then the phase transition will be always completed
before the universe cools down to $T^{\rm NB}$, i.e. $T^{\rm NB}<T_{\rm nuc}$,
which provides a lower bound for the nucleation temperature. At the same time
the velocity of the bubbles become largest for the smallest possible values
for the nucleation temperature. So that the fastest bubbles will be near
$T^{\rm NB}=0$ curve. Looking at Fig.1 we can see that the largest
$\gamma_{w}$ (Lorentz boost factor) and lowest nucleation temperatures happen
for $v_{s}\gtrsim 200$ GeV, where the $T^{\rm NB}=0$ curve passes very close
to the NO PT boundary. The shape of the lines in Fig.1 clearly indicate the
necessity of tuning in order to obtain low nucleation temperatures (large
$\gamma_{w}$). In particular for the values of $v_{s}\gtrsim 200$ GeV we can
see that the nucleation temperature drops by choosing $\lambda_{hs}$ close to
the NB value (similarly $\gamma_{w}$ becomes maximal see Fig.1).

Figure 3: Tuning of the coupling $\lambda_{hs}$ as a function of the
nucleation temperature. The dashed green line represent the naive tuning
$\sim(T_{\rm nuc}/m_{h})^{2}$. We observe that this naive estimation for the
tuning is rather precise at large nucleation temperature but can underestimate
the tuning by one order of magnitude for very low nucleation temperature.

We can estimate this tuning by looking at
$\partial\log\lambda_{hs}/\partial\log T_{\rm nuc}$ quantity (analogue of
Giudice-Barbieri [79] measure of the tuning) as a function of $T_{\rm nuc}$.
This result agrees with our expectation from the steepness of the curves in
the Fig.3 and with the naive tuning expectation which scales as
$\sim\big{(}T_{\rm nuc}/m_{h}\big{)}^{2}$555This expression follows Eq.(37) if
we require the cancellation between the terms independent of temperature..

At last we would like to remind that the discussion in this section always
assumed that the phase transition completes before the potential barrier
disappears. We have checked numerically that this is always the case. Indeed
the time of the phase transition is given approximately by the bubble radius
at the moment of percolation [80, 81]

$\displaystyle
R_{\star}\equiv\frac{(8\pi)^{1/3}}{\tilde{\beta}}=\left(\int_{T_{\rm
per}}^{T_{\rm crit}}\frac{dT}{T}\frac{\Gamma(T)}{H(T)}\left(\frac{T_{\rm
per}}{T}\right)^{3}\right)^{-1/3}.$ (41)

This radius is related to the
$\beta\equiv-\frac{d}{dt}\frac{S_{3}}{T}=HT\frac{d}{dT}\frac{S_{3}}{T}$
parameter by an approximate relation [81]

$\displaystyle{\beta\over H}\biggl{|}_{T_{\rm
per}}\simeq\frac{\tilde{\beta}}{H}=\frac{(8\pi)^{1/3}}{R_{\star}H}\;,$ (42)

where we find $R_{\star}^{-1}\sim\beta_{\rm typical}\sim(10-10^{4})H.$ At this
point the temperature drops during the bubble expansion will scale as

$\displaystyle\Delta T\sim T_{\rm nuc}(H\Delta t)\sim\frac{T_{\rm
nuc}H}{\beta}\;.$ (43)

Due to the large value of $\beta/H$ we find numerically that this drop of the
temperature is not enough for the barrier to disappear or in other words

$\displaystyle T_{\rm nuc}-\Delta T>T^{\rm NB}\;.$ (44)

Such behaviour can be understood from the following consideration near $T_{\rm
NB}$ the bounce action drops very quickly and so that the tunneling becomes
very efficient almost instantaneous and typical bubble radiuses are much
smaller than the Hubble scale. This leads to another prediction that GW signal
will be suppressed as well since it is controlled by the ($\beta/H$) quantity
Eq.(42). As we will see, even with this suppression the GW signal is efficient
enough to be detected in the future.

###  5.2 Tunneling with potential barrier at zero temperature

Let us proceed to the analysis of the case when the potential barrier does not
disappear at zero temperature. The parameter space with the lowest nucleation
temperatures (fastest bubbles) will be located again near the “NO PT”
boundary. However in this case the nucleation temperature will be controlled
by the local minima of the $S_{3}/T$ function. (see Fig.2). At least a minimum
is expected since the potential at low temperature becomes fully temperature
independent and $S_{3}(T\to 0)\to\text{const.}$, so that $S_{3}/T$ necessarily
starts to grow for $T\to 0$.

Numerically (see Fig.1) for the value of $m_{s}=125$ GeV we find that for
$v_{s}\lesssim 180$ GeV entire region with the fast bubbles has a potential
barrier at zero temperature. On Fig.2, we present the euclidean action for
$v_{s}=170$ GeV and $m_{s}=125$ GeV. Going back to Fig.1, we see that for
those values the “NO PT” curve and the “$T_{\rm NB}=0$” curve are largely
separated. This is not a surprise since in this region of parameter space the
bounce action $S_{3}/T\sim O(10^{2})$ is small enough to guarantee the
successful tunneling even when the barrier remains at zero temperature. In the
range of $\lambda_{hs}$ from $0.544$ to $0.58$, numerically we find that
nucleation temperatures are between $65-44$ GeV, with a clear saturation at
$T^{\rm sat}_{\rm nuc}\approx 44$ GeV, and the corresponding Lorentz factor
for the velocities of the bubble expansion in the ranges of $\lesssim 10$.
Interestingly we find that in this case the bubble radius
$R_{\star}\sim\frac{(8\pi)^{1/3}v_{w}}{\beta}$ are a little bit larger than
the ones discussed in the section 5.1, corresponding to a bit smaller values
of $\beta/H$ parameter.

##### Super fine-tuned region

We finally comment on the parameter space with $v_{s}\gtrsim 200$ GeV (again
we are fixing $m_{s}=125$ GeV), where the curves “NO PT” and “$T_{\rm NB}=0$”
almost superimpose (the region between red and blue dots on the $T_{\rm nuc}$
panel of the Fig.1). There will be a very narrow strip between the curves “NO
PT” and “$T_{\rm NB}=0$” regions, where the tunnelling will happen even though
the barrier remains at zero temperature. We find (see Fig.5) that the region
corresponds to the variations of the $\lambda_{hs}$ parameter of the order
$\delta\lambda_{hs}\sim O(10^{-4})$, i.e. two order of magnitude smaller than
the full region with relativistic bubbles. In this very small region various
additional effects can start playing a role. For example let us look at the
Fig.5 we can see that the bounce action $S_{3}/T$ has a local maximum and a
deeper (global) minimum with respect to the standard scenario. Such a
behaviour of the action is coming from the cancellations of various terms in
the effective potential. For simplicity let us look at $T=0$ case. Then there
is a region of parameter space where purely polynomial potential has no local
minimum at $(0,v_{s})$, but the effects of the
$\frac{-3M_{t}^{4}(h)}{8\pi^{2}}\big{(}\log\frac{M_{t}(h)}{M_{t}(v_{EW})}\big{)}$
terms in Coleman-Weinberg contribution lead to the appearance of the local
minimum at $(\delta v_{h},v_{s})$. On Fig.4 we plot the contributions of the
various terms in the effective potential leading to the appearance of this
local minimum and the trajectory of the typical bounce solution in this case.
As a result the distance in the fields space between the two minima decreases
and the tunneling becomes faster, which leads to the appearance of the second
(global) minimum in $S_{3}/T$.

Figure 4: Here is presented an explicit example of displacement of the false
vacuum for $\\{\lambda_{hs},v_{s},T\\}=\\{0.36784,220\ {\rm GeV},3\ {\rm
GeV}\\}$. Left: we show the 2D potential where the blue line corresponds to
the part of the potential plotted in the right panel and the (purple) dot is
the position of the displaced false minimum, in both the plots. Red line
indicates the bounce trajectory. Right: plot of the different contributions to
the potential. We see that a displaced minimum can be generated balancing the
tree level and the CW potential of the top quark, for low enough temperature,
in such a way all the other particles, that are massless in the false vacuum,
have a negligible contribution. It can be shown that they cause, as the
temperature increases, the shift of the local minimum towards $h=0$.

####  5.2.1 Benchmark points

Figure 5: Left: Plot of $S_{3}(T)/T$ as a function of the temperature, for
different values of $\lambda_{hs}$ in the case where $m_{s}=125$ GeV and
$v_{s}=$ 220 GeV. As we increase the value of the coupling $\lambda_{hs}$, the
disappearance of the potential barrier happens later, allowing for longer
supercooling, until it is large enough to remain even at zero temperature. For
the first four curves from the top, we observe a second drop in the function
at very low temperature. This second drop corresponds to the displacement on
the false minimum that we describe in this section. Right: Same plot as in the
left panel, but with a lower value for $v_{s}$. The pattern we found is the
same, but lowering $v_{s}$ causes a lowering of the curves and the
displacement of the false minimum is less pronounced.

In Table 1 and 2, we give typical values of the nucleation temperature, the
Lorentz factor $\gamma_{w}$, the $\beta/H$ factor, and we indicate if the
barrier remains at zero temperature, applying the criterion in Eq.(44). We can
see that the largest bubble radius at the collision (smallest $\beta$)
correspond to the case when $S_{3}/T$ is monotonic and very flat near the
tunneling temperature (c.f. the right panel of Fig. 5).

$m_{s}=125$ GeV, $v_{s}=170$ GeV |
---|---
$\lambda_{hs}$ | $\frac{T_{\rm reh}}{100\text{GeV}}$ | $\frac{T_{\rm nuc}}{100\text{GeV}}$ | $\frac{T_{\rm per}}{100\text{GeV}}$ | $\gamma_{w}$ | ${\tilde{\beta}\over H}={(8\pi)^{1/3}\over R_{\star}H}$ | $m_{H}^{\rm False}/\text{GeV}$ |  FM0
$0.56$ | $0.880$ | $0.877$ | $0.850$ | $-$ | $434$ | $35.4$ | No
$0.58$ | $0.855$ | $0.851$ | $0.822$ | $-$ | $355$ | $37.3$ | No
$0.6$ | $0.829$ | $0.824$ | $0.790$ | $-$ | $296$ | $39.2$ | No
$0.62$ | $0.800$ | $0.795$ | $0.762$ | $1.2$ | $209$ | $40.1$ | No
$0.64$ | $0.769$ | $0.762$ | $0.714$ | $2.4$ | $158$ | $42.5$ | No
$0.66$ | $0.729$ | $0.720$ | $0.661$ | $4$ | $108$ | $43.4$ | No
$0.68$ | $0.678$ | $0.666$ | $0.582$ | $6.6$ | $51$ | $44.8$ | No
$0.69$ | $0.642$ | $0.627$ | $0.506$ | $8.8$ | $18$ | $45.0$ | No
$0.695$ | $0.612$ | $0.594$ | $0.412$ | $11$ | $5$ | $44.7$ | No
$0.6997$ | $0.566$ | $0.542$ | $0.237$ | $15$ | $1.4$ | $43.8$ | No
Table 1: We report for Fig.2, $m_{s}=125$ GeV and $v_{s}=170$ GeV, reheating, nucleation and percolation temperatures, respectively, for different values of $\lambda_{hs}$ and $\gamma_{w}$ reached by the expanding walls as well as the parameter $\tilde{\beta}/H$ computed using Eq.(42). We also show the effective Higgs mass in the false vacuum at the nucleation temperature defined as $(m_{H}^{\rm False})^{2}={\partial^{2}V\over\partial h^{2}}\big{|}_{T=T_{\rm nuc}}$, relevant for DM production explained in section 6.2. In the last column, FM0 concerns the displacement of the false minimum. No if it is at the (Higgs) origin, yes if it moved. In this case, the minimum is always at the origin.  $m_{s}=125$ GeV, $v_{s}=205$ GeV |
---|---
$\lambda_{hs}$ | $\frac{T_{\rm reh}}{100\text{GeV}}$ | $\frac{T_{\rm nuc}}{100\text{GeV}}$ | $\frac{T_{\rm per}}{100\text{GeV}}$ | $\gamma_{w}$ | ${\tilde{\beta}\over H}={(8\pi)^{1/3}\over R_{\star}H}$ | $m_{H}^{\rm False}/\text{GeV}$ |  FM0
$0.397$ | $0.577$ | $0.564$ | $0.544$ | $4$ | $371$ | $19.1$ | No
$0.405$ | $0.530$ | $0.512$ | $0.488$ | $8$ | $268$ | $19.1$ | No
$0.4155$ | $0.448$ | $0.412$ | $0.379$ | $18$ | $130$ | $17.7$ | No
$0.42$ | $0.393$ | $0.330$ | $0.290$ | $37$ | $72$ | $15.2$ | No
$0.423$ | $0.339$ | $0.161$ | $0.124$ | $270$ | $66$ | $7.1$ | No
$0.4234$ | $0.335$ | $0.107$ | $0.095$ | $805$ | $109$ | $3.9$ | No
$0.424$ | $0.335$ | $0.051$ | $0.051$ | $5.7\cdot 10^{3}$ | $3.3\cdot 10^{3}$ | $0.7$ | No
$0.4242$ | $0.335$ | $0.0337$ | $0.0337$ | $1.8\cdot 10^{4}$ | $3.2\cdot 10^{4}$ | $0.25$ | No
$0.42424$ | $0.335$ | $0.028$ | $0.0279$ | $3.0\cdot 10^{4}$ | $1.8\cdot 10^{3}$ | $4.4$ | No
$0.424266$ | $0.335$ | $0.018$ | $0.017$ | $1.0\cdot 10^{5}$ | $99$ | $6.2$ | Yes
$0.424267$ | $0.335$ | $0.016$ | $0.014$ | $1.3\cdot 10^{5}$ | $44$ | $6.3$ | Yes
Table 2: Same as Table 1, but for Fig.5 and with $v_{s}=205$ GeV. We observe
that the last two points display a displacement of the false minimum.

##  6 Consequences for production of dark matter and Baryogenesis

One of the motivation for the study of a model of EWPT with relativistic
bubbles is the relation between relativistic expansion and the out-of-
equilibrium production of heavy states presented for the first time in [29],
when the field undergoing the PT (here the Higgs) is coupled to some heavy
dark sector at typical mass $M_{N}$. In this section, we remind the principle
of the production mechanism and we study the scenario of the production of
Dark Matter[2] and Baryogenesis[1], that were previously agnostic about the
EWPT realisation.

First of all, the strong FOPT involves a supercooling represented by a
dilution factor,

$\displaystyle D\equiv\frac{g^{\rm sym}_{\star s}(T_{\rm nuc})}{g_{\star
s}(T_{\rm reh})}\left(\frac{T_{\rm nuc}}{T_{\rm reh}}\right)^{3},$ (45)

with $g_{\star s}(T)$ ($g^{\rm sym}_{\star s}(T)$) being the number of
relativistic degrees of freedom of the entropy in the broken (symmetric)
phase. This means that with $D\ll 1$, any type of dark matter production or
Baryogenesis mechanism that happens much earlier than the PT should provide
values denser than the conventional estimation by a factor of $1/D$ (see for
example [19, 82]). For instance, the WIMP cross section should be
$\sigma\sim\frac{D}{10^{-3}}10^{-29}\rm cm^{3}/s$ to produce a correct dark
matter abundance. This is the case when freeze-out happens at temperatures
much higher than the reheating.

###  6.1 Production of heavy states during ultra-relativistic expansion

There are few mechanisms which can lead to heavy particle production during
FOPT. This can happen if the incoming massless particle in the unbroken phase
gets a very large mass from the Higgs vev [3, 64] (mass gain), or due to the
bubble- bubble collision [25, 30] or due to the plasma particle$-$bubble
collision [29]. Our study will be focused on the later one.

Let us assume that FOPT happens and the bubble expansion is indeed
relativistic with $\gamma_{w}^{\rm terminal}\gg 1$. The simplest model where
the production of heavy particles during plasma$-$bubble wall collision can be
realized, is described by the following Lagrangian [29]:

$\displaystyle\mathcal{L}={1\over
2}(\partial_{\mu}h)^{2}+i\bar{q}\mathord{\not\mathrel{{\mathrel{\partial}}}}q+i\bar{N}\mathord{\not\mathrel{{\mathrel{\partial}}}}N-M_{N}\bar{N}N-Yh\bar{N}q-V(h)\;,$
(46)

where $q$ is a massless particle in the symmetric phase and $N$ is a heavy
field with large vev$-$independent mass $M_{N}\gg v_{EW}$. $h$ is the Higgs
field undergoing a FOPT. With no loss of generality, we go to the basis where
fermion masses are real. Before the strong phase transition starts, the
abundance of heavy states $N$ in the plasma is strongly Boltzmann suppressed
and they would naively seem irrelevant for the dynamics of the transition. In
an homogeneous vacuum, the transition from light to heavy state $q\to N$ is
obviously forbidden by the conservation of momentum. However, in the presence
of the bubble wall, the conservation of momentum along the $z$ direction is
broken (assuming a planar wall expanding in the $x-y$ plane) and a computation
using WKB phases for the $q$ and $N$ fields demonstrates that the probability
$\mathcal{P}(q\to N)$ is non-vanishing [29] and is given by

$\displaystyle\mathcal{P}(q\to
N)\approx\frac{Y^{2}v_{EW}^{2}}{M_{N}^{2}}\Theta(\gamma_{w}T_{\text{nuc}}-M_{N}^{2}L_{w})\;,$
(47)

with $L_{w}\sim 1/v_{EW}$ the width of the wall. Behind the bubble wall a
large abundance of $N$ and $\bar{N}$, $n_{N}^{BE}$ is produced. Let us
emphasize that this abundance is much larger than its equilibrium value.

Another possibility of the heavy particle production can be realized for the
following interaction

$\displaystyle\Delta\mathcal{L}\supset\frac{\lambda_{h\phi}}{2}\phi^{2}h^{2}+\frac{1}{2}M^{2}_{\phi}\phi^{2}.$
(48)

In this case $\phi$ is a heavy scalar field with mass $M_{\phi}\gg v_{EW}$,
then in the vicinity of the wall, the process $h\to\phi\phi$ has the
probability [29]

$\displaystyle\qquad\mathcal{P}(h\to\phi^{2})\approx\bigg{(}\frac{\lambda_{h\phi}v_{EW}}{M_{\phi}}\bigg{)}^{2}\frac{1}{24\pi^{2}}\Theta(\gamma_{w}T_{\text{nuc}}-M_{\phi}^{2}L^{h}_{w}).$
(49)

The results in Sec. 3.1 on the terminal velocity in the singlet extension of
SM allow us to compute the maximal mass of the particles which can be produced
during the electroweak FOPT in the singlet extension. Indeed saturating the
step function in the above equation and assuming the $L_{w}\sim 1/v_{EW}$ we
get approximately:

$\displaystyle
M^{MAX}\approx\frac{{400}{\rm~{}GeV}}{\log^{1/2}\frac{M_{z}}{gT_{\rm
nuc}}}\left(\frac{\Delta V-\Delta P_{\rm
LO}}{(100{\rm~{}GeV})^{4}}\right)^{1/2}\left(\frac{100{\rm GeV}}{T_{\rm
nuc}}\right).$ (50)

Numerical results for the maximal mass $M^{\rm MAX}$ are reported in Fig.1. We
can see that the maximal mass we can produce is roughly $\sim 10$ TeV scale.

We would like to note that our results can be easily applied for the mass gain
mechanism of the heavy state production [3]. Indeed in this case the maximal
mass will be $M_{\rm mass~{}gain}\simeq\gamma_{w}T$, and can be read off from
the bottom right plot of the Fig.1 by noting that it will scale as $M_{\rm
mass~{}gain}\sim M_{\rm MAX}^{2}/v_{EW}$. Since the mass of the heavy field
comes from the vev of the Higgs, it will additionally be bounded by the
unitarity considerations to be below $\lesssim 2$ TeV.

###  6.2 Dark Matter production

In this section we will apply the results for the velocity of the bubble
expansion for DM model building.

####  6.2.1 Scalar DM coupled to the Higgs portal

We assume a heavy scalar $\phi$ coupled to the SM via the traditional Higgs
portal

$\displaystyle{\cal
L}_{DM}=\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{M_{\phi}^{2}\phi^{2}}{2}-{\frac{\lambda_{\phi
h}}{2}}h^{2}\phi^{2}\;.$ (51)

The DM ($\phi$) field is stabilized by some additional $Z_{2}^{\phi}$ (we use
this subscript to differentiate it from $Z_{2}$ of the singlet potential).
After the Higgs transition, the abundance of massive $\phi$, $n^{\rm
BE}_{\phi}$, behind the wall is given by

$\displaystyle n_{\phi}^{\text{BE}}$ $\displaystyle\approx$
$\displaystyle\frac{2}{\gamma_{w}v_{w}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p_{z}}{p_{0}}\mathcal{P}(h\to\phi^{2})\times
f_{h}(p,T_{\text{nuc}})\;.$ (52)

We can see that DM production during the bubble expansion is strongly
dependent on the density of the Higgs field available at the nucleation
temperature $f_{h}(p,T_{\text{nuc}})$. The relevant parameter for the
discussion is the ratio

$\displaystyle\sqrt{\frac{d^{2}V}{dh^{2}}\bigg{|}_{\rm
fv}}\frac{1}{T}\equiv\frac{m_{H}^{\rm False}}{T}\;,$ (53)

where fv denotes the position of the false vacuum and
$\sqrt{d^{2}V/dh^{2}\big{|}_{\rm fv}}$ is the mass of the Higgs $m_{H}^{\rm
False}$ in the false vacuum. As soon as this quantity becomes larger than 1,
we expect exponential suppression of the Higgs abundance

$\displaystyle C_{\rm
eff}\frac{\zeta(3)T_{\text{nuc}}^{3}}{\pi^{2}}\equiv\int\frac{d^{3}p}{(2\pi)^{3}}f_{h}(p,T_{\text{nuc}})\approx\begin{cases}\frac{\zeta(3)T_{\text{nuc}}^{3}}{\pi^{2}}\quad\text{if}\quad
m_{H}^{\rm False}<T\;,\\\ \bigg{(}\frac{m_{H}^{\rm
False}T_{\text{nuc}}}{2\pi}\bigg{)}^{3/2}e^{-m_{H}^{\rm
False}/T_{\text{nuc}}}\quad\text{if}\quad m_{H}^{\rm False}>T\;.\end{cases}$
(54)

Here we define $C_{\rm eff}$ to take into account the Boltzmann suppression.
After redshifting to today, the stable produced abundance takes the form

$\displaystyle\Omega^{\text{today}}_{\phi,\text{BE}}h^{2}$
$\displaystyle=\frac{M_{\phi}n_{\phi}^{\text{BE}}}{\rho_{c}/h^{2}}\frac{g_{\star
S0}T_{0}^{3}}{g_{\star S}(T_{\text{reh}})T_{\text{reh}}^{3}}\approx 6.3\times
10^{8}\;\frac{M_{\phi}n_{\phi}^{\text{BE}}}{\text{GeV}}\frac{1}{g_{\star
S}(T_{\text{reh}})T_{\text{reh}}^{3}}\;,$ $\displaystyle\approx{5.4}\times
10^{5}\times\bigg{(}\frac{C_{\rm
eff}\lambda_{h\phi}^{2}v_{EW}}{M_{\phi}g_{\star
S}(T_{\text{reh}})}\bigg{)}\bigg{(}\frac{v_{EW}}{\text{GeV}}\bigg{)}\bigg{(}\frac{T_{\text{nuc}}}{T_{\text{reh}}}\bigg{)}^{3}\times
e^{-\frac{M_{\phi}^{2}}{2\gamma_{w}v_{EW}T_{\rm nuc}}}\;.$ (55)

This expression has to be supplemented with the _freeze-out_(FO) contribution
which is produced _before_ the phase transition

$\displaystyle\Omega^{\text{today}}_{\phi,\text{FO}}h^{2}$
$\displaystyle\approx$ $\displaystyle
0.1\times\bigg{(}\frac{T_{\text{nuc}}}{T_{\text{reh}}}\bigg{)}^{3}\times\bigg{(}\frac{0.03}{\lambda_{\phi
h}}\bigg{)}^{2}\bigg{(}\frac{M_{\phi}}{100\text{ GeV}}\bigg{)}^{2}\;,$
$\displaystyle\Omega^{\text{today}}_{\phi,\text{tot}}h^{2}$ $\displaystyle=$
$\displaystyle\Omega^{\text{today}}_{\phi,\text{BE}}h^{2}+\Omega^{\text{today}}_{\phi,\text{FO}}h^{2}\;.$
(56)

Note that the FO contribution is suppressed by the factor
$\big{(}T_{\text{nuc}}/T_{\text{reh}}\big{)}^{3}$ due the brief stage of
inflation during the phase transition. Obviously the prediction for relic
density must match the experimental observations:
$\Omega^{\text{today}}_{\phi,\text{tot}}h^{2}\approx 0.1$. We can see from
Eqs.(55),(6.2.1) that for small values of the portal coupling
$\lambda_{h\phi}$, DM production will be dominated by the freeze out mechanism
while bubble expansion takes over for larger values of $\lambda_{h\phi}$.

Next we can check whether this mechanism for DM production can lead to viable
phenomenology, given the results on bubble dynamics in section 5. Instead of
making a scan of the parameter space, we will just focus on a few
representative benchmark points.

Figure 6: Left: DM abundance in the parameter space $\lambda_{\phi
h}-M_{\phi}$ for different values of $\lambda_{hs}$ by fixing $v_{s}=205$ GeV
and $m_{s}=125$ GeV (that fixes the values of $T_{\rm nuc},T_{\rm reh}$ and
$\gamma_{w}^{\rm terminal}$). The solid lines represent correct DM abundance,
while underproduced inside and overproduced outside. The lower part of each
contour is dominated by freeze-out and the upper part via bubble expansion.
The connecting vertical line (independent of the portal) comes from thermal
production after the reheating of the transition. The magenta shaded region is
excluded by XENON1T while the dotted green and blue lines are projected limits
from XENONnT and DARWIN respectively. Right: same plot for $v_{s}=175$ GeV and
$m_{s}=150$ GeV. As expected, increasing the tuning from red to magenta
increases the amplitude of the curve. The values used are extracted from Table
3.

For $m_{s}=125,v_{s}=205$ GeV, we show in Fig.6 the isocontours reproducing
the correct relic density for three reference values of $\lambda_{hs}$
(corresponding nucleation temperatures can be found in the Table 1). Firstly,
$C_{\rm eff}\simeq 1$ for all three reference points. For $\lambda_{hs}=0.424$
the upper red curve corresponds to the case when DM production is dominated by
the BE (bubble expansion) and the lower curve by FO. The steepness of the
upper red curve (BE) comes from the fact that we are always in the region of
parameter space where $\exp\big{[}{-M_{\phi}^{2}/(2\gamma_{w}v_{EW}T_{\rm
nuc})}\big{]}\ll 1$, leading to a very strong sensitivity on $M_{\phi}$ mass.
Physically this means that the model generically predicts large overproduction
of DM in BE process unless the Boltzmann suppression
$\exp\big{[}{-M_{\phi}^{2}/(2\gamma_{w}v_{EW}T_{\rm nuc})}\big{]}$ is playing
a role. For the other two reference points $\lambda_{hs}=0.4242,0.42424$ we
can see that there is an additional part of parameter space for the DM masses
$M_{\phi}\sim 1-4$ TeV, which corresponds to the region without the Boltzmann
suppression $\exp\big{[}{-M_{\phi}^{2}/(2\gamma_{w}v_{EW}T_{\rm
nuc})}\big{]}\sim 1$. This is related to larger values of $M^{\rm
MAX}\sim\sqrt{\gamma_{w}v_{EW}T_{\rm nuc}}$ and smaller values of the
nucleation temperature, reducing the excess of the DM abundance. On the right
panel of Fig.6, we report similar plots for $v_{s}=175,m_{s}=150$ GeV.

Finally, before closing this section, we comment about the possibility of
considering the singlet $s$ itself, in the limit of very precise $Z_{2}$, as
DM. After the phase transition $T\sim 40$ GeV, the singlet is in thermal
equilibrium and we can apply straightforwardly the freeze-out expression:

$\displaystyle\Omega^{\text{today}}_{s,\text{ FO}}h^{2}\approx
0.1\bigg{(}\frac{0.06}{\lambda_{hs}}\bigg{)}^{2}\bigg{(}\frac{M_{s}(v_{EW},0)}{100\text{
GeV}}\bigg{)}^{2}\;.$ (57)

From this estimate of the FO abundance for $s$ and recalling that we
considered $\lambda_{hs}\sim 0.3-0.6$ and $M_{s}(v_{EW},0)\sim 100$ GeV, we
conclude that the abundance of $s$ produced in this fashion, today, is
underproduced by one or two orders of magnitude to fit the observed amount of
DM $\Omega^{\text{today}}_{s,\text{ FO}}h^{2}\approx 0.1$. Even in this
underproduced case, there are severe bounds from the direct detection
experiments except for the resonant region, where (57) is over-estimated.
However, as we will discuss in the Appendix. D, we will have a $Z_{2}$
explicit breaking which makes $s$ decay much before today.

####  6.2.2 Singlet portal DM

Figure 7: DM production for the singlet portal model. In this context, the
Boltzmann suppression $\sim e^{-m_{s}/T_{\rm nuc}}$ plays a strong role and
allows for bubble expansion produced DM with much higher nucleation
temperature $T_{\rm nuc}\sim 15$ GeV.

In this section we mention an alternate possibility of coupling the DM
($\phi$) to the singlet field $s$ via “singlet portal”

$\displaystyle\Delta\mathcal{L}\supset\frac{\lambda_{s\phi}}{2}\phi^{2}s^{2}+\frac{1}{2}M^{2}_{\phi}\phi^{2}\;,\quad\mathcal{P}(s\to\phi^{2})\approx\
\bigg{(}\frac{\lambda_{s\phi}v_{s}}{M_{\phi}}\bigg{)}^{2}\frac{1}{24\pi^{2}}\Theta(\gamma_{w}T_{\text{nuc}}-M_{\phi}^{2}L^{s}_{w})\;,$
(58)

where the width of the singlet wall is similar to the length of the Higgs wall
$L_{w}$. Interestingly even though FOPT is from $(0,v_{s})\to(v_{EW},0)$, the
singlet scattering of the wall can lead to the production of the $\phi$ field.
Phenomenology of DM production is very similar to the Higgs portal case
discussed in the previous section, with one main difference: in the false
vacuum, the mass of the singlet is not small and the factor $C_{\rm eff}$
introduced in the Eq.(54) plays an important role. The results are shown in
Fig.7. For example, if we compare the curves for
$\lambda_{hs}=0.424,m_{s}=125,v_{s}=205$ in Fig.7 and in Fig.6, we can see
that, for the Higgs portal DM, the isocontour has the same shape, but with the
larger values of DM $M_{\phi}$ masses. As shown in Eq.(55), this is due to the
proportionality between the DM relic abundance produced during the bubble
expansion with $\propto C_{\rm eff}/M_{\phi}$.

##### Singlet portal with additional field

A slight modification of this scenario is to further introduce a light scalar
$\tilde{s}$. Then we can have,

$\displaystyle{\cal L}={\cal
L}_{SM}-\tilde{\lambda}\tilde{s}s\phi^{2}-\frac{M_{\phi}^{2}}{2}\phi^{2},$
(59)

where again $\phi$ is the DM and we did not write down the mass terms for
simplicity of notation. We assume that $\tilde{s}$ is in the thermal bath
before the PT. Then due to the field change of $s$ in the bubble wall, the
momentum conservation violating process $\tilde{s}\to\phi\phi$ can occur
($h\to\phi\phi$ may also occur if there is the $h^{2}\phi^{2}$ term.) In this
model $s\to\tilde{s}+SM$ particles happen via the DM loop.

####  6.2.3 Fermion-mediated Dark Matter

In the previous section we noticed that the DM production during bubble
expansion strongly depends on the mass of the incoming particle in the
symmetric phase, due to the Boltzmann suppression factor. In case of an
incoming scalar, generically this effect is relevant and crucially modifies
the phenomenology, as we have seen in the section 6.2.2. In this section we
construct a model where the incoming particle is a massless fermion in the
symmetric phase, so that $C_{\rm eff}\equiv 1$ by definition. The model
consists of a vector-like neutral fermion $N$ which is a singlet under SM and
a couple of $Z_{2}$ odd fields $\phi$ and $\chi$:

$\displaystyle{\cal L}={\cal
L}_{SM}+Y_{*}\bar{L}HN+M_{N}\bar{N}N+Y_{DM}\bar{N}\chi\phi\;.$ (60)

Here, $L,H$ are SM lepton and Higgs doublets, respectively. The production
mechanism works as follows: the heavy field $N$ is produced during the phase
transition $L\to N$ and it will subsequently decay into $N\to\chi\phi,N\to
LH$. The field $N$ can be Majorana or Dirac (in the former case there will be
a relation to neutrino masses and in the later it will be completely
independent from neutrinos).

In this model, heavy $N$ are produced via $L\to N$ with a probability

$\displaystyle\mathcal{P}^{\rm tree}(L\to
N)\approx\frac{Y_{\star}^{2}v_{EW}^{2}}{M_{N}^{2}}\Theta(\gamma_{w}T_{\text{nuc}}-M_{N}^{2}L_{w})\;.$
(61)

As a consequence, unstable heavy $N$ accumulate behind the wall with initial
density given by

$\displaystyle n_{N}^{\text{BE}}$ $\displaystyle\approx$
$\displaystyle\frac{Y_{\star}^{2}v_{EW}^{2}}{M_{N}^{2}\gamma_{w}v_{w}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p_{z}}{p_{0}}\times
f_{L}(p,T_{\text{nuc}})\Theta(p_{z}-M_{N}^{2}/v_{EW})\;,$ (62)
$\displaystyle\approx$
$\displaystyle\frac{Y_{\star}^{2}v_{EW}^{2}T_{\text{nuc}}^{3}}{2\pi^{2}M_{N}^{2}}e^{-\frac{M_{N}^{2}}{2v_{EW}T_{\text{nuc}}\gamma_{w}}}+\mathcal{O}(1/\gamma_{w})\;,$

where $v_{w}=\sqrt{1-1/\gamma_{w}^{2}}$, we expanded for large $\gamma_{w}$
and approximated the Fermi-Dirac distribution as a Boltzmann distribution.
Compared to the original proposal in Ref.[2], the density of the heavy fields
inside the bubble will be additionally enhanced by $\sim 16\pi^{2}$ factor
since $1\to 1$ transitions are more effective than $1\to 2$. Let us assume
that $M_{\phi}<M_{\chi}$ so that $\phi$ is the DM candidate, then DM
production will happen via the following chain of processes:

$\displaystyle\chi\underset{\text{via PT}}{\rightarrow}N\underset{\text{via
decay}}{\rightarrow}\phi\chi\to\phi\phi+\text{SM}\;.$ (63)

However, the heavy $N$ has two channels of decay: toward the heavy dark sector
$\phi,\chi$ and back to the light $L$. The abundance of $\phi,\chi$ after the
transition is thus suppressed and given by

$\displaystyle n_{\chi}\approx
n_{\phi}\approx\frac{Y_{DM}^{2}Y_{\star}^{2}}{Y^{2}_{DM}+Y^{2}_{\star}}\frac{v_{EW}^{2}T_{\text{nuc}}^{3}}{2\pi^{2}M_{N}^{2}}e^{-\frac{M_{N}^{2}}{2v_{EW}T_{\text{nuc}}\gamma_{w}}}+\mathcal{O}(1/\gamma_{w})\;,$
(64)

and the final relic abundance redshifted to today thus reads

$\displaystyle\Omega^{\text{today}}_{\phi,\text{BE}}h^{2}$
$\displaystyle\approx$ $\displaystyle{1.5}\times
10^{8}\times\frac{Y^{2}_{\star}Y_{DM}^{2}}{Y_{\star}^{2}+Y^{2}_{DM}}\frac{2M_{\phi}}{M_{N}}\bigg{(}\frac{v_{EW}}{M_{N}}\bigg{)}\bigg{(}\frac{v_{EW}}{246{\rm
GeV}}\bigg{)}\bigg{(}\frac{T_{\text{nuc}}}{T_{\text{reh}}}\bigg{)}^{3}e^{-\frac{M_{N}^{2}}{2v_{EW}T_{\text{nuc}}\gamma_{w}}}\;.$

For the freeze-out process in the symmetric phase, we have: $\phi\phi\to LHLH$
by neglecting co-annihilation. The cross-section is highly phase space
suppressed (closing a loop for a 2 to 2 annihilation gives a similar scaling):
$\sigma_{\phi\phi\to(LH)^{*}LH}\sim\frac{M_{\chi}^{2}(Y_{\rm
DM}Y_{*})^{2}}{(16\pi^{2})^{2}4\pi M_{N}^{4}}.$ The abundance by taking
account the supercooling is

$\displaystyle\Omega^{\text{today}}_{\phi,\text{FO}}h^{2}=10^{3}\left(\frac{T_{\text{nuc}}}{T_{\text{reh}}}\right)^{3}\frac{M_{N}^{4}/M_{\chi}^{2}}{(6{\rm
TeV})^{2}}\frac{10}{(Y_{\rm DM}Y_{*})^{4}}\;.$ (66)

The total DM density today will be given by the sum of
Eq.(LABEL:eq:relic_ab_alt)-(66). Therefore, this scenario leads to the over-
production of DM unless $M_{\phi},M_{\chi}\lesssim 10$ GeV. However, note that
these equations are valid only for the heavy DM candidates which do not go
back to equilibrium after the phase transition. Otherwise, the final density
will be given by Eq.(66) only without
$\big{(}T_{\text{nuc}}/T_{\text{reh}}\big{)}^{3}$ and we are going back to the
normal freeze-out scenario.

Let us now investigate the regime $M_{\phi}\simeq M_{\chi}$, precisely
$|M_{\phi}-M_{\chi}|\lesssim M_{\phi}/20$, where the co-annihilation takes
place. In this case we have the channel $\phi\chi\to H{\bar{L}}$ to decrease
the abundance of $\phi$ as well as $\chi$. The cross-section is
$\sigma_{\phi\chi\to H{\bar{L}}}\sim\frac{(Y_{\rm DM}Y_{*})^{2}}{4\pi
M_{N}^{2}}.$ Therefore, we have

$\displaystyle\Omega^{\text{today}}_{\phi,\text{FO,Co}}h^{2}\sim
0.1\left(\frac{T_{\text{nuc}}}{T_{\text{reh}}}\right)^{3}\times\frac{M_{N}^{2}}{(10{\rm
TeV})^{2}}\frac{1}{(Y_{\rm DM}Y_{*})^{2}}\;.$ (67)

Summing this estimate with the $\Omega^{\text{today}}_{\phi,\text{BE}}h^{2}$
in Eq.(LABEL:eq:relic_ab_alt) we find that it becomes possible to reproduce
the observed DM abundance. However we see that bubble expansion tends to
overproduce the DM and the relic abundance in BE can be reproduced if only the
factor $\exp[-M_{N}^{2}/(2v_{EW}T_{\rm nuc}\gamma_{w})]$ starts playing a role
in suppressing DM relic density (left boundary of Fig.8 is almost vertical).

Figure 8: DM production in the fermionic portal model. We fix $M_{\rm
DM}=M_{N}/5$, and $Y_{\rm DM}=3$. Only the co-annihilation regime where
$|M_{\phi}-M_{\chi}|\lesssim M_{\phi}/20$ appears to be viable and shown by
dashed lines.

###  6.3 Baryogenesis mechanism

Now we remind the scenario of Baryogenesis with relativistic bubble walls that
was proposed in [1]. As a prototype, we worked with the following model,
reminiscent of the toy model of Eq.(46), (omitting the kinetic terms)

$\displaystyle{\cal L}=$ $\displaystyle{\cal
L}_{SM}+m_{\eta}^{2}|\eta|^{2}+\sum_{I=1,2}M_{I}\bar{B}_{I}B_{I}$
$\displaystyle+$
$\displaystyle\left(\sum_{I=1,2}Y_{I}(\bar{B}_{I}H)P_{L}Q+{y_{I}\eta^{*}\bar{B}_{I}P_{R}\chi}+\kappa\eta^{c}du+\frac{1}{2}m_{\chi}\bar{\chi^{c}}\chi+h.c.\right).$
(68)

Thus, additionally to the SM and the singlet sector, the model contains a
Majorana field $\chi$ and two vector-like $B$ quarks with the masses
$M_{1,2}\sim m_{\chi}$. $\eta$ is a scalar field in the fundamental
representation of QCD and with electric charge $Q(\eta)=1/3$. We defined
$Q,u,d$ as the SM quark doublet and singlets respectively, and ignored the
flavour indices. $H$ is the SM Higgs.

We first proved that the production mechanism, when computed up to one loop-
level, indeed transforms a CP-violating phase into a chiral asymmetry in the
abundances produced, in a fashion very similar to usual leptogenesis decay.
Let us now sketch the mechanism in itself. First, $b-$quarks collide with the
relativistic wall and produce $B_{I},B_{I}^{c}$ via the mechanism explained
above. Thus inside the bubble we have

$\displaystyle
n_{B_{I}}-n_{B^{c}_{I}}=-\bigg{(}\frac{Y_{I}v_{EW}}{M_{I}}\bigg{)}^{2}\epsilon_{I}n_{b}^{\rm
inc}\;,\qquad
n_{b}-n_{b^{c}}=\sum_{I}\bigg{(}\frac{Y_{I}v_{EW}}{M_{I}}\bigg{)}^{2}\epsilon_{I}n_{b}^{\rm
inc}\;,$ (69)

where $n_{b}^{\rm inc}$ is the number density of the bottom-type quark
colliding with the wall from outside, $n_{b}$ and $n_{B}$ are the abundances
inside of the bubble and $\epsilon$ is a loop suppressed coefficient which
parametrizes the $CP$ violating phase and the resonance between one-loop and
tree-level diagrams. After the passage through the wall, the following
asymmetric abundances are then generated

$\displaystyle\sum_{I}\left(n_{B_{I}}-n_{B_{I}^{c}}\right)=-(n_{b}-n_{b^{c}})\;.$
(70)

From this expression we can read an apparent asymmetry in the bottom quark
abundances. However, if the heavy $B$ freshly produced were to decay back into
$b$, the asymmetry would be washed out. This is however not the case if they
decay in a dark sector containing $\chi,\eta$, where the asymmetry is enhanced
by the presence of a Majorana mass for $\chi$. The final unsuppressed produced
asymmetry is given by

$\displaystyle\frac{\Delta n_{\rm Baryon}}{s}\approx$
$\displaystyle\frac{135\zeta(3)}{8\pi^{4}}\sum_{I,{J}}\bigg{(}\frac{Y_{I}\langle{H}\rangle}{M_{I}}\bigg{)}^{2}\frac{|y_{I}|^{2}}{|y_{I}|^{2}+|Y_{I}|^{2}}\times\frac{g_{b}}{g_{\star}}\bigg{(}\frac{T_{\rm
nuc}}{T_{reh}}\bigg{)}^{3}$ (71) $\displaystyle\times\;{\rm
Im}(Y_{I}Y^{*}_{J}y_{I}^{*}y_{J})\left(-\frac{2{\rm
Im}[f_{B}^{IJ}]}{|Y_{I}|^{2}}+\frac{4{\rm Im}[f_{B}^{IJ}]|_{m_{\chi,\eta}\to
0}}{|y_{I}|^{2}}\right)\;,$

where $s$ is the entropy at the moment of the production, $g_{b}$ is number of
degrees of freedom of the bottom and $g_{\star}$ the number of relativistic
degrees of freedom. The loop functions $f_{B}^{IJ}$ have been computed in [1]
and are controlled by the CP violating sector. Absence of strong wash-out
conditions

$\displaystyle M_{B,\chi,\eta}\gtrsim 30T_{\rm reh}\sim 10^{3}\text{
GeV}\;,\qquad\text{(suppressed wash-out)}$ (72)

as well as experimental signatures (direct production in colliders, flavor
violation, neutron oscillations) pushed the heavy particles to be

$\displaystyle M_{B,\chi,\eta}\gtrsim 2\times 10^{3}\text{ GeV}\;.$ (73)

In the context of singlet extension with $Z_{2}$ that we studied, this opens
up the range

$\displaystyle M_{B,\chi,\eta}\in[2,10]\text{ TeV}\;,$ (74)

where the Baryogenesis mechanism proposed above is operative.

###  6.4 Impact of the heavy sector on the phase transition

The models we are considering by construction have new heavy fields coupled to
the Higgs boson. These will lead to the finite corrections to the scalar
parameters of the form (assuming a Yukawa type connection $yBHb$)

$\displaystyle\delta
m_{h}^{2}\sim\frac{-8g_{N}M_{N}^{2}y^{2}}{64\pi^{2}}\bigg{(}\log\frac{M_{N}^{2}}{v_{EW}^{2}}-\frac{3}{2}\bigg{)},\qquad\delta\lambda\sim
4g_{N}\frac{y^{4}}{64\pi^{2}}\bigg{(}\log\frac{M_{N}^{2}}{v_{EW}^{2}}-\frac{3}{2}\bigg{)},$
(75)

where $g_{N}$ is the number of heavy degrees of freedom and $M_{N}$ is the
typical mass of the heavy sector. One can wonder how these corrections can
effect the tuning of the Higgs potential. However note that in our model the
Higgs mass hierarchy problem is not addressed and generically we expect the
size of $m_{h}^{2}$ to be of the order of the cut off scale ($M_{\rm pl}$ in
SM). So the corrections in Eq.(75) do not make the tuning worse.

In case the Higgs hierarchy problem is solved at the scale of the heavy fields
in Eq.(51),(6.3) the tuning in the Higgs potential will be roughly,

$\displaystyle{\rm tuning}\sim\frac{m_{h}^{2}}{\Delta_{\rm
Heavy~{}physics}m_{h}^{2}}\sim\frac{8\pi^{2}m_{h}^{2}}{y^{2}M_{N}^{2}}.$ (76)

We can combine this estimate with a tuning for low nucleation temperatures
(see discussion in section 5.1) which are necessary for the heavy field
production and the tuning estimate becomes:

$\displaystyle{\rm tuning}\sim\frac{T^{2}_{\rm nuc}}{\Delta_{\rm
Heavy~{}physics}m_{h}^{2}}\sim\frac{8\pi^{2}T_{\rm nuc}^{2}}{y^{2}M_{N}^{2}}.$
(77)

Using the estimates of the maximal values of $\gamma_{w}$ and the maximal mass
of heavy particles which can be produced during the bubble-plasma collisions
(see Eq.(31) and Eq.(61)) we get the following estimate for the maximal tuning
in the model

$\displaystyle{\rm tuning}^{MAX}\sim\left(\frac{T_{\rm nuc}}{20{\rm
GeV}}\right)^{4},$ (78)

where we remind the reader that this estimate is valid only if the Higgs
hierarchy problem is solved at the heavy fields scale.

##  7 Gravitational waves emitted

It is well known that strong gravitational waves background will be emitted,
with peak frequency around the mHz if the EWPT happens to be strongly first
order. This is the optimal range of sensitivity of the forthcoming LISA
detector [83, 84] and also GW detectors such as eLISA[85], LIGO[86, 87],
BBO[88, 89], DECIGO[90, 91, 92], ET[93, 94, 95], AION[96], AEDGE[97].

The signal produced at the moment of the transition can be separated into
different contributions: the _bubble collision_[98] contribution, the _plasma
sound waves_[84] and finally the _turbulence_. Only the two first sources of
GW are well understood. Another nice feature of those two sources is that they
are expected to dominate in different physical situations; the bubble
collision would dominate in case of runaway wall and the sound waves if the
wall reaches a terminal velocity. We have already mentioned that the EWPT, if
first order, will always happen in the regime of terminal velocity, because of
the large number of strongly coupled vector bosons666A possible exception
would be the case of extreme cooling, as hinted in [67] where even the
pressure from gauge bosons cannot stop the acceleration of the wall. However,
in our study, we do not find such a situation.. For GW produced by plasma
sound wave, the peak frequency and amplitude are given by

$\displaystyle\Omega_{\text{plasma}}^{\text{peak}}h^{2}$ $\displaystyle\approx
0.7\times
10^{-5}\bigg{(}\frac{100}{g_{\star}}\bigg{)}^{1/3}\bigg{(}\frac{\kappa_{sw}\alpha}{1+\alpha}\bigg{)}^{2}(H_{\text{reh}}R_{\star}),$
$\displaystyle f_{\text{peak}}$ $\displaystyle\approx 2.6\times
10^{-5}\bigg{(}\frac{1}{H_{\text{reh}}R_{\star}}\bigg{)}\bigg{(}\frac{z_{p}}{10}\bigg{)}\bigg{(}\frac{T_{\text{reh}}}{100\text{
GeV}}\bigg{)}\bigg{(}\frac{g_{\star}}{100}\bigg{)}^{1/6}\text{ Hz}$ (79)

with $z_{p}\sim 10$, $\kappa_{sw}$ is the efficiency factor for the production
of sound waves in the plasma[99],

$\displaystyle\kappa_{sw}\approx\frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha},$
(80)

$\alpha$ and $R_{\star}$ have been defined in Eqs.(19) and (41) respectively,
$R_{\star}\sim\mathcal{O}(10^{-1}-10^{-4})H^{-1}$ is the approximate size of
the bubble at collision and all quantities $(T,H,g_{\star})$ have to be
evaluated at _reheating_.

As we have seen in sections 6.3-6.2, for the baryogenesis and DM production we
need relativistic walls with relatively low nucleation temperature $\lesssim
10$ GeV. In this context, $\alpha\gg 1$ and $\kappa_{sw}\to 1$. The peak
frequency and the signal amplitude are only function of the size of the
bubbles at collision, which are reported in Table 2 and 3. We can observe that
in this range $\beta/H$ spans the value between $[50,10^{4}]$, with a
preference for lower values. Going back to Eq.(79), emitted amplitude and
frequencies will be of the order

$\displaystyle\Omega_{\text{plasma}}^{\text{peak}}h^{2}\in[5\times
10^{-7},2\times 10^{-9}],\qquad
f_{\text{peak}}\in[10^{-4},0.03]\quad\text{Hz}$ (81)

where we set $z_{p}=10,g_{\star}=100$. This range of frequencies and amplitude
are largely in the expected sensitivity of the coming observer LISA[85, 84],
as expected for this class of models[34]. We thus conclude that strong GW
signal in the LISA with spectrum controlled by the plasma sound waves is a
generic prediction of Baryogenesis with relativistic bubble walls. This is in
sharp opposition with the general expectation that usual EWBG demands slow
walls, and thus suppressed signals.

As a final comment, it should however be noticed that the current simulations
do not directly provide a solutions for the regime of large $\alpha$, and we
only have an extrapolation of the numerical result. Thus, the conclusion above
should be taken with a grain of salt.

##  8 Conclusion

In this study we have presented the first explicit realization of the
baryogenesis and DM production during electroweak phase transition for ultra-
relativistic bubble expansion. The work is based on the proposals in [29, 2,
1, 3] where new heavy particles are produced in plasma$-$bubble wall
collisions. We have shown that the model with SM extended by a real singlet
with a $Z_{2}$ symmetry can indeed lead to ultra-relativistic bubbles, where
the Lorentz factor $\gamma_{w}$ can reach the values $\sim 10^{5-6}$. Such
fast bubbles can appear if the symmetry breaking occurs in two steps: first
discrete $Z_{2}$ is spontaneously broken and in the second step electroweak
symmetry breaking is accompanied by $Z_{2}$ restoration
$(0,0)\xrightarrow{SOPT}(0,v_{s})\xrightarrow{FOPT}(v_{EW},0)$. We find that
there exists a region of parameter space where the nucleation temperature can
become as low as $1-2$ GeV and the collision of the bubble wall with the
plasma particles can lead to the non-thermal heavy particle production with
the masses up to $\sim 10$ TeV. Interestingly we find that the mechanism is
most efficient for relatively low masses of the singlet field
${M_{s}(v_{EW},0)}\sim 70-100\ {\rm GeV}$, close to the region excluded by the
Higgs invisible decays. Subsequently, this region of parameter space will be
probed by HL-LHC ([70, 100]) in the singlet production mediated by off-shell
Higgs boson. By noting the slight $Z_{2}$ breaking, $s$, if produced, can
decay into $b\bar{b}$ in collider experiments. Depending on the size of the
breaking displaced vertices of $b\bar{b}$ may be probed. We find the typical
bubble radius parameter of the order of $R_{\star}\sim(10^{-4}-1)H^{-1}$ so
that stochastic gravitational background signal becomes observable at GW
experiments like LISA[85, 84].

The model necessarily requires tuning $\propto(T_{\rm nuc}/m_{h})^{2}$ which
numerically turns out to be of the order of $10^{-4}-10^{-2}$ (using Giudice-
Barbieri measure) for successful baryogenesis and DM production mechanism. In
spite of this we believe it can provide a useful guidance for more appealing
models where these hierarchies can appear naturally.

## Acknowledgements

AA and SC were supported by the MIUR contract 2017L5W2PT. WY was supported by
JSPS KAKENHI Grant Nos.20H05851, 21K20364, 22K14029, and 22H01215. MV thanks
the DESY institute for hospitality during the writing of this manuscript as
well as Alberto Mariotti and Simone Blasi for fruitful discussions in the late
stages of the completion of the paper.

##  Appendix A The bounce in two dimensions

In this paper we studied numerically the phase transition from the minimum
$(0,v_{s})$ (or in the vicinity of it) to $(v_{EW},0)$. The bounce computation
can be done using existing codes for example FindBounce or CosmoTransition.
However we have found that in the regime of long supercooling where the
potential around the false vacuum is very flat and, the existing codes are
often not stable and lead to numerical errors. Thus we have developed our own
code (more stable for the flat potentials), following the procedure described
in [101], while cross-checking the available values with FindBounce.

###  A.1 Computation of the bounce profile

In this Appendix, we briefly review the standard computation of the bounce
action with only one field before going to describe the algorithm we used for
the same computation but for the case of two fields PT. In order to compute
the vacuum tunneling probability from the false vacuum to the true one in $d$
dimensions, we need to minimize the Euclidean action given by

$\displaystyle S_{E}=\int
d^{d}x\bigg{[}\frac{1}{2}(\partial_{\mu}\phi)^{2}+V[\phi]\bigg{]}.$ (82)

It is known that the field configurations leading to the minimal action are
the ones that exhibit an $O(d)$ spherical symmetry, then the so-called bounce
solution is the solution of the following Cauchy problem

$\displaystyle\frac{d^{2}\phi}{dr^{2}}+\frac{d-1}{r}\frac{d\phi}{dr}=\frac{dV[\phi]}{d\phi},\qquad\lim_{r\to\infty}\phi(r)=0,\qquad\frac{d\phi}{dr}\bigg{|}_{r=0}=0,$
(83)

where we have chosen the false minimum to be at $\phi=0$. If we interpret the
parameter $r$ as a time and $\phi$ as a position, this problem becomes
formally equivalent to the evolution of a mechanical ball in a potential
$-V[\phi]$ undergoing a friction $\frac{d-1}{r}\frac{d\phi}{dr}$, released
with vanishing velocity and stopping its evolution for $r\to\infty$ at
$\phi=0$. It is well known that this problem can be solved by applying
numerically an overshoot/undershoot method on the position of the released
point. Releasing the ball too close to the true vacuum would induce an
overshoot configuration (the ball would continue after crossing $\phi=0$), we
would thus shift the release point toward the false vacuum, while releasing it
too close would end up in an undershoot configuration (the ball would never
reach $\phi=0$ and starts oscillate around the minimum of $-V[\phi]$) and we
correct it by shifting the release point farther from the false vacuum.
Iterating between those two situations, we are able to find the correct
release point and obtain the _bounce solution_.

It is well known that the case of a PT triggered by temperature fluctuation at
temperature $T$ is formally equivalent to imposing a periodicity $T^{-1}$ in
the imaginary time $t_{E}$, which imposes the following constraint on the
field

$\displaystyle\phi[t_{E},\vec{x}]=\phi[t_{E}+T^{-1},\vec{x}]$ (84)

and the computation of thermally induced phase transition thus amounts to take
$d=3$ in the above equations.

###  A.2 Bounce action in two dimensions and path deformation

The problem complexifies when the transition involves many fields. Here there
is no straightforward intuition for the path followed by the fields in field
space during the tunneling. One can think that a straight line, connecting the
two minima, could be a reasonable guess, but it turns out that it cannot be
considered as a good approximation of the euclidean action777Let us emphasize
that in the region of the parameter space we studied, the straight line
between the false and the true vacuum gives an Euclidean action which is often
wrong by _orders of magnitudes_ , as the path is often very far from the
straight line, as a consequence, we cannot dispense from the effort of
studying the exact 2D path in field space.. Here we thus describe the
algorithm [101] to find the right path in field space. In a multi-field case,
the Eq.(83) becomes

Figure 9: On the left, the potential, along the path, experienced by the field
$x$ and the corresponding escape point $x_{\star}$. On the right the
potential, along the path, projected on the $h$ direction.

$\displaystyle\frac{d^{2}\vec{\phi}}{dr^{2}}+\frac{d-1}{r}\frac{d\vec{\phi}}{dr}=\vec{\nabla}V[\vec{\phi}],\qquad\lim_{r\to\infty}\vec{\phi}(r)=0,\qquad\frac{d\vec{\phi}}{dr}\bigg{|}_{r=0}=0.$
(85)

Since in this case an overshoot/undershoot procedure cannot be easily applied,
the idea is to reduce the problem to one dimensional tunneling. In order to do
so we start guessing the path, $\vec{\phi}_{g}(x)$, where $x$ is now to be
understood as the parameter that measure the distance along the path, i.e. the
so-called curvilinear abscissa. For the present case, if we parametrize the
path in the field space as $(h(t),s(t))=(t,f(t))\equiv(h,s(h))$ it is defined
as

$\displaystyle
x(h)=\int^{h}_{h_{fm}}\sqrt{1+\left(\frac{ds(h^{\prime})}{dh^{\prime}}\right)^{2}}dh^{\prime},$
(86)

where $h_{fm}$ is the value of the Higgs field in the false minimum. With this
choice of field coordinate to parametrize the path, the condition
$\big{|}\frac{d\vec{\phi}(x)}{dx}\big{|}^{2}=1$ is satisfied, and the
Euclidean equation of motion in Eq. (85) can be rewritten along the parallel
and the perpendicular direction

$\displaystyle\frac{d^{2}x}{dr^{2}}+\frac{d-1}{r}\frac{dx}{dr}=\partial_{x}V[\vec{\phi}_{g}(x)],$
$\displaystyle\frac{d^{2}\vec{\phi}_{g}(x)}{dx^{2}}\bigg{(}\frac{dx}{dr}\bigg{)}^{2}=\vec{\nabla}_{\perp}V[\vec{\phi}_{g}(x)].$
(87)

Here, we have been able to separate the dynamics along the parallel and
perpendicular direction in such a way the first equation defines a new
undershoot/overshoot problem, that we solve to obtain the value of the _escape
point_ , $\vec{\phi}_{0}(x_{\star})$, and the Euclidean action corresponding
to the potential along the path considered $\vec{\phi}_{g}$, as in Fig.9. On
the other hand, the second equation can be seen as a condition that the bounce
solution has to satisfy and can be thought as a force field acting on the
path, defined as following

$\displaystyle\vec{N}\equiv\frac{d^{2}\vec{\phi}_{g}(x)}{dx^{2}}\bigg{(}\frac{dx}{dr}\bigg{)}^{2}-\vec{\nabla}_{\perp}V[\vec{\phi}_{g}(x)].$
(88)

The right path will be the one where $\vec{N}$ is vanishing. The algorithm
proceeds iteratively: first we guess a path, the straight line connecting the
two minima, then we find the bounce solution along this path, we compute the
normal force and deform the guessed path according to it. In practice, to
define the path at the step $n$, $\vec{\phi}_{n}$, we need to solve for the
bounce profile for the path at $\vec{\phi}_{n-1}$, extract the escape point
$x_{\star,n-1}$, that is $(h(x_{\star,n-1}),s(x_{\star,n-1}))$ in field space,
we then discretize the path in the interval $x\in[0,x_{\star,n-1}]$, creating
a series $(\vec{\phi}_{n-1})_{j}$, for $j=1,...,N$ and a series of values for
the normal force $(\vec{N}_{n-1})_{j}$. We then shift each point of the
discretized path by

$\displaystyle(\vec{\phi}_{n})_{j}=(\vec{\phi}_{n-1})_{j}+\rho(\vec{N}_{n-1})_{j}\qquad
j=1,\dots,N.$ (89)

In the end, we fit a path $\vec{\phi}_{n}$ along the shifted points from
$(\vec{\phi}_{n-1})_{j}$. The procedure of deformation of the path will
produce a series of paths $\vec{\phi}_{i}[x]$, over which we compute the
Euclidean action according to Eq.(82) at each step of the deformation, like in
Fig.10. The algorithm stops when the difference in the bounce action, $S_{3}$,
between two successive iterations is below some imposed precision. At a
definite temperature $T$, we start by identifying the two minima, the false
and the true ones

Figure 10: Left: iterative procedure for the correct path, starting from the
straight line connecting the two minima and then modified according to the
field $\vec{N}$. Right: bounce profile of the fields (black for $h(r)$ and red
for $s(r)$) on the correct path.

$\displaystyle(\langle h\rangle,\langle s\rangle)_{\rm
fm}=(v(T),v_{s}(T))\to(v_{EW},0),$ (90)

and will keep the false minimum fixed during the whole procedure of
deformation. Generally, especially when we have a sizable amount of
supercooling, the escape point is just behind the barrier, so the escape point
$(v_{\star}(T_{n}),v_{s,\star}(T_{n}))$ will be different from the, zero-
temperature, EWSB vacuum, but when the tunneling happens the system will
classically roll down towards the global minimum, as we can see from Fig.10.
We do not track the evolution of the fields profile after the tunneling.

##  Appendix B Supplemental numerical results

Figure 11: Here is presented the same results found in Fig. 6 of Ref. [47]. It
has to be noted that these results are obtained with only the thermal
potential and without Daisy resummation, i.e. without thermal masses. The
relation with our parameters is $(\lambda_{s},\ \lambda_{hs})=(\eta,\
2\kappa)$ and $M_{s}(v_{EW},0)=300$ GeV.

Figure 12: Similar plots than in Fig.1 and Fig.3 for the value of $m_{s}=150$ GeV. Figure 13: On this plot we show the regime of transition. We observe that several points display a disappearance of the barrier which is typical of the regime of no barrier at $T=0$. However, the nucleation temperature is controlled by the first minimum of $S_{3}/T$, which is typical of the regime with a barrier at $T=0$. $m_{s}=150$ GeV, $v_{s}=175$ GeV |
---|---
$\lambda_{hs}$ | $\frac{T_{\rm reh}}{100\text{GeV}}$ | $\frac{T_{\rm nuc}}{100\text{GeV}}$ | $\frac{T_{\rm per}}{100\text{GeV}}$ | $\gamma_{w}$ | ${\tilde{\beta}\over H}={(8\pi)^{1/3}\over R_{\star}H}$ | $m_{H}^{\rm False}/\text{GeV}$ |  FM0
$0.53$ | $0.616$ | $0.608$ | $0.592$ | $-$ | $580$ | $17.4$ | No
$0.54$ | $0.581$ | $0.570$ | $0.552$ | $3$ | $420$ | $17.8$ | No
$0.56$ | $0.492$ | $0.470$ | $0.444$ | $9$ | $228$ | $17.7$ | No
$0.58$ | $0.329$ | $0.141$ | $0.130$ | $348$ | $160$ | $4.7$ | No
$0.582$ | $0.327$ | $0.051$ | $0.0508$ | $5.3\cdot 10^{3}$ | $9.5\cdot 10^{3}$ | $0.4$ | No
$0.582262$ | $0.327$ | $0.025$ | $0.0236$ | $3.7\cdot 10^{4}$ | $194$ | $8.2$ | Yes
$0.582264$ | $0.327$ | $0.024$ | $0.0219$ | $4.3\cdot 10^{4}$ | $130$ | $8.3$ | Yes
$0.582266$ | $0.327$ | $0.021$ | $0.017$ | $5.7\cdot 10^{4}$ | $24$ | $8.5$ | Yes
Table 3: Same as Tables 1 and 2, but with $m_{s}=150\ {\rm GeV}$ and
$v_{s}=175$ GeV.

In this Appendix, we present all our supplemental numerical results. First,
tough we focused mostly on a more weakly coupled part of the parameter space,
we would like to compare our findings with the ones in the Ref.[47] and argue
that we observed only small changes, due to the inclusion of loop-corrections
and Daisy resummation. On Fig.11 we make a reproduction of the scan of the
Fig.6 of [47] using our potential and emphasize the close similarities. The
relations between the parameters $\kappa,\eta$ and the couplings in the Eq.(1)
is as follows

$\displaystyle\kappa_{\hbox{\cite[cite]{[\@@bibref{}{Kurup:2017dzf}{}{}]}}}=\frac{\lambda_{hs}}{2},~{}~{}~{}\eta_{\hbox{\cite[cite]{[\@@bibref{}{Kurup:2017dzf}{}{}]}
}}=\lambda_{s}.$ (91)

In the main text, we studied specifically the case where the parameter
$m_{s}=125$ GeV, we observed that for this value, the region of deep
supercooling displayed small masses of the singlet in the real vacuum, being
on the verge of detection due to $h\to ss$ at $M_{s}\lesssim 62$ GeV. We also
concluded in section 5 that this region was closing around $M_{s}\approx 75$
GeV. We could wonder if this conclusion would change if we modify the value of
the parameter $m_{s}$, and if so in which direction. On Fig.12 we show similar
plots than in Fig.1 and 3 for the case of $m_{s}=150$ GeV. Thus, increase the
value of $m_{s}$ pushes the deep supercooling region to $M_{s}\approx 90$ GeV,
at the price of increasing the portal coupling $\lambda_{hs}$. However, we can
observe on the last plot of Fig.12 that the typical tuning remains roughly the
same and that we can still trust our naive $(T_{\rm nuc}/m_{h})^{2}$ for an
order-of-magnitude estimate of the tuning.

On the other hand, we also observed that decreasing the parameter $m_{s}$ to
$\approx 100$ GeV was pushing all the deep supercooling region inside
$M_{s}\lesssim 62$ GeV, which is thus strongly disfavored by colliders. We
hope that this trend can be extrapolated to larger values of $m_{s}$, until we
hit perturbativity bounds for $\lambda_{hs}$.

We could also wonder what happens at the _upper_ boundary of the deep
supercooling region, as we have observed on Fig.1 a sharp decrease in the
supercooling allowed around $v_{s}\lesssim 200$ GeV (for $m_{s}=125$ GeV).
This transition regime can be understood if we plot the explicit $S_{3}/T$
functions on Fig.13. Comparing the plot in Fig.13 with the one in Fig.5, we
see that as we decrease $v_{s}$, the full pattern of $S_{3}/T$ is shifted
toward smaller values. At some critical point around $v_{s}\approx 200$ GeV,
the nucleation becomes controlled by the first minimum in the function
$S_{3}/T$ and not by the disappearance of the barrier. This largely suppresses
the possibility for large supercooling.

Finally, in Table 3, we provide the value of the velocity, reheating and
nucleation temperature for $m_{s}=150$ GeV and $v_{s}=175$ GeV that was used
in Fig.6.

##  Appendix C The coefficient of NLO pressure

In this appendix we will review the calculation of the friction coefficient
for the NLO pressure for EW phase transition. We will follow closely the
discussion in [67] and report the quantity

$\displaystyle\bigg{[}\sum_{abc}\nu_{a}g_{a}\beta_{c}C_{abc}\bigg{]}$ (92)

where $\nu_{a}=1(3/4)$ for $a$ a boson (fermion),
$\beta_{c}\equiv\frac{M_{c}}{M_{Z}}$and $C_{abc}$ stands for the couplings
appearing in the vertex. Normalization of the $C_{abc}$ coefficient is the
following: for a chiral fermion coupled to the vector field the amplitude for
the process $\psi\to\psi A_{soft}$ is equal to

$\displaystyle
g_{\psi}\bar{\psi}_{L}A_{\mu}\psi_{L}\Rightarrow~{}~{}~{}C_{\psi\psi
A}=\frac{g_{\psi}^{2}}{4\pi\alpha_{em}}.$ (93)

Where in the relation Eq.(93) is written only for one polarization of the
vector field. Similarly for the scalar field

$\displaystyle
ig_{\phi}(\phi^{*}\partial_{\mu}\phi-\partial_{\mu}\phi^{*}\phi)A^{\mu}\Rightarrow~{}~{}~{}C_{\phi\phi
A}=\frac{g_{\psi}^{2}}{4\pi\alpha_{em}},$ (94)

and the vector fields

$\displaystyle
g_{V^{1}V^{2}A^{3}}\left(V^{1}_{\mu\nu}V^{2}_{\mu}A_{\nu}+V^{2}_{\mu\nu}A_{\mu}V^{1}_{\nu}+A_{\mu\nu}V^{1}_{\mu}V^{2}_{\nu}\right)\Rightarrow
C_{V_{1}V_{2}A}=\frac{g_{V^{1}V^{2}A}^{2}}{4\pi\alpha_{em}},$ (95)

where in all of these formulas $C_{abc}$ coefficients are reported only for
one polarization of the vector fields both in the initial and the final
states. Summing all of these contributions and taking care of the
multiplicities of the initial and final states we find

$\displaystyle\bigg{[}\sum_{abc}\nu_{a}g_{a}\beta_{c}C_{abc}\bigg{]}={2}\left(\frac{7+14c_{w}}{s_{w}^{2}}-\frac{7-15s_{w}^{2}}{c_{w}^{2}}\right)\simeq
157.$ (96)

For the interested reader we refer various individual contributions in the
Table 4. If in the false vacuum the Higgs doublet ${\cal H}$ is too heavy its
contribution must be subtracted and the sum in Eq.(96) reduces to

$\displaystyle\bigg{[}\sum_{abc}\nu_{a}g_{a}\beta_{c}C_{abc}\bigg{]}_{\rm
No~{}Higgs}\simeq 145.$ (97)

At last we would like to emphasize that these results include only the
transverse polarizations of the vector fields. NLO effects of the longitudinal
polarizations are not fully established and we omit them here, however these
cannot qualitatively modify the results.

Process | $\sum g_{a}C_{abc}$ | $\beta$ | $\nu$ | Result
---|---|---|---|---
$\psi\to W^{\pm}\psi$ | $\frac{24}{s_{w}^{2}}$ | $c_{w}$ | $\frac{3}{4}$ | $\frac{18c_{w}}{s_{w}^{2}}$
$\psi\to Z\psi$ | $\frac{4(3-6s_{w}^{2}+8s_{w}^{4})}{s_{w}^{2}c_{w}^{2}}$ | 1 | $\frac{3}{4}$ | $\frac{3(3-6s_{w}^{2}+8s_{w}^{4})}{s_{w}^{2}c_{w}^{2}}$
${\cal H}\to W{\cal H}$ | $\frac{2}{s_{w}^{2}}$ | $c_{w}$ | $1$ | $\frac{2c_{w}}{s_{w}^{2}}$
${\cal H}\to Z{\cal H}$ | $\frac{1-2s_{w}^{2}c_{w}^{2}+c_{w}^{4}+s_{w}^{4}}{2s_{w}^{2}c_{w}^{2}}$ | $1$ | $1$ | $\frac{1-2s_{w}^{2}c_{w}^{2}+c_{w}^{4}+s_{w}^{4}}{2s_{w}^{2}c_{w}^{2}}$
$A\to W_{soft}W$ $\&$ $W\to W_{soft}A$ | $8$ | $c_{w}$ | $1$ | $8c_{w}$
$Z\to W_{soft}W$ $\&$ $W\to W_{soft}Z$ | $\frac{8c_{w}^{2}}{s_{w}^{2}}$ | $c_{w}$ | $1$ | $\frac{8c_{w}^{3}}{s_{w}^{2}}$
$W\to Z_{soft}W$ | $\frac{4c_{w}^{2}}{s_{w}^{2}}$ | $1$ | $1$ | $\frac{4c_{w}^{2}}{s_{w}^{2}}$
Total: |  |  |  | $2\left(\frac{7+14c_{w}}{s_{w}^{2}}-\frac{7-15s_{w}^{2}}{c_{w}^{2}}\right)\simeq 157$
Table 4: Different contributions to the sum in Eq.(28).

##  Appendix D Domain wall collapse

Our main discussion was focused on the two step phase transition
$(0,0)\to(0,v_{s})\to(v_{h},0)$ where the first phase transition is $Z_{2}$
breaking. Obviously during such a phase transition domain walls will be formed
which can drastically modify the cosmology of the system. We can avoid the
stable domain walls if we assume some small $Z_{2}$ breaking, however in this
case the question rises about the timescale for the stability of the domain
walls. This is particularly important since recently it was shown [102] that
for singlet extension of the SM the domain walls (if still present) will
become seeds of the secondary phase transition $(0,v_{s})\to(v_{h},0)$ and
will dominate the phase transition. We will follow closely the discussion in
the section 4.2 using only the tree level potential and the thermal
corrections to the masses. Then the $Z_{2}$ breaking phase transition will
occur at the temperatures

$\displaystyle T_{Z_{2}}\simeq
v_{s}\left(\frac{12}{3+4\lambda_{hs}\frac{v_{s}^{2}}{m_{s}^{2}}}\right)^{1/2}\simeq
200-300~{}~{}\rm GeV,$ (98)

which is a temperature of the domain wall formation. The domain wall mediated
transition will happen at the temperature $T_{w}$ which is found to be order
one different from $T_{Z_{2}}$. The exact mechanism of the transition depends
on the values of the couplings and can proceed either with the classical
rolling or 2D bounces localized on the domain wall. The temperature when the
classical rolling can start is reported in Ref.[102] and is equal to

$\displaystyle T_{w}^{\rm rolling}\simeq
T_{Z_{2}}\left[\frac{4m_{h}^{2}+m_{s}^{2}\bigg{(}1-\sqrt{1+\frac{8v_{s}^{2}\lambda_{hs}}{m_{s}^{2}}}\bigg{)}}{8\Pi_{h}(T_{Z_{2}})+m_{s}^{2}\bigg{(}1-\sqrt{1+\frac{8v_{s}^{2}\lambda_{hs}}{m_{s}^{2}}}\bigg{)}}\right]^{1/2}.$
(99)

The nucleation temperature $(T_{w}^{2D})$ of $2D$ bounces should be found
numerically (Ref.[102]) however it will be obviously smaller than $T_{\rm
crit}$ (of $(v_{s},0)\to(v_{h},0)$ phase transition). At this point we can
safely ignore the seeded phase transition effects if all of the domain walls
annihilate in the interval of temperatures

$\displaystyle\left[T_{*},T_{Z_{2}}\right],~{}~{}~{}T_{*}<T_{\rm crit},$ (100)

where $T_{*}$ is the temperature when the seeded phase transition will be
completed and it is obviously less than $T_{\rm crit}$ of EW phase transition.
Let us estimate how strong should be the bias $\Delta V$ between the potential
energies of the two minima of $Z_{2}$ potential so that all of the walls can
disappear. For these estimates it is sufficient to assume that there is order
one difference between $T_{*}$ and $T_{Z_{2}}$, which is generically the case.
The critical radius (above which) areas with true vacuum will start to expand
is roughly

$\displaystyle R_{c}\sim\frac{\sigma}{\Delta V},$ (101)

where $\sigma$ is surface energy density of the wall. So the domain walls will
exist on the time scale of

$\displaystyle\Delta t_{w}\sim\frac{R_{c}}{u}\sim\frac{\sigma}{u\Delta V},$
(102)

where $u$ is velocity of the wall motion. The change of the temperature during
the wall annihilation will be roughly

$\displaystyle\Delta T\sim TH\Delta t_{w}.$ (103)

So that if $\Delta t_{w}H\ll 1\Rightarrow\Delta T\ll T$ the wall annihilation
happens almost instantaneously. Assuming $\sigma\sim T_{\rm crit}^{3}$ and
$H\sim\frac{T_{\rm crit}^{2}}{M_{p}}$ we get

$\displaystyle\frac{\Delta V}{T_{\rm crit}^{4}}\gg\frac{T_{\rm
crit}}{uM_{pl}}.$ (104)

Balancing the pressure against the friction forces $\Delta V\sim uT_{\rm
crit}^{4}$ we can estimate the velocity and then the condition for the quick
wall annihilation becomes

$\displaystyle\frac{\Delta V}{T_{\rm crit}^{4}}\gg\sqrt{\frac{T_{\rm
crit}}{M_{pl}}}\sim 10^{-8},$ (105)

which is not restrictive at all.

## References

  * [1] A. Azatov, M. Vanvlasselaer, and W. Yin JHEP 10 (2021) 043, [arXiv:2106.14913].
  * [2] A. Azatov, M. Vanvlasselaer, and W. Yin JHEP 03 (2021) 288, [arXiv:2101.05721].
  * [3] I. Baldes, S. Blasi, A. Mariotti, A. Sevrin, and K. Turbang arXiv:2106.15602.
  * [4] A. D. Sakharov Pisma Zh. Eksp. Teor. Fiz. 5 (1967) 32–35. [Usp. Fiz. Nauk161,no.5,61(1991)].
  * [5] V. Kuzmin, V. Rubakov, and M. Shaposhnikov Phys. Lett. B 155 (1985) 36\.
  * [6] M. Shaposhnikov JETP Lett. 44 (1986) 465–468.
  * [7] K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov Nucl. Phys. B 493 (1997) 413–438, [hep-lat/9612006].
  * [8] M. D’Onofrio and K. Rummukainen Phys. Rev. D 93 (2016), no. 2 025003, [arXiv:1508.07161].
  * [9] D. E. Morrissey and M. J. Ramsey-Musolf New J. Phys. 14 (2012) 125003, [arXiv:1206.2942].
  * [10] D. N. Schramm Nucl. Phys. B 252 (1985) 53–71.
  * [11] M. J. Baker and J. Kopp Phys. Rev. Lett. 119 (2017), no. 6 061801, [arXiv:1608.07578].
  * [12] M. J. Baker and L. Mittnacht JHEP 05 (2019) 070, [arXiv:1811.03101].
  * [13] L. Heurtier and H. Partouche Phys. Rev. D 101 (2020), no. 4 043527, [arXiv:1912.02828].
  * [14] T. Cohen, D. E. Morrissey, and A. Pierce Phys. Rev. D 78 (2008) 111701, [arXiv:0808.3994].
  * [15] L. Bian and Y.-L. Tang JHEP 12 (2018) 006, [arXiv:1810.03172].
  * [16] L. Hui and E. D. Stewart Phys. Rev. D 60 (1999) 023518, [hep-ph/9812345].
  * [17] D. Chung, A. Long, and L.-T. Wang Phys. Rev. D 84 (2011) 043523, [arXiv:1104.5034].
  * [18] D. J. H. Chung and A. J. Long Phys. Rev. D 84 (2011) 103513, [arXiv:1108.5193].
  * [19] T. Hambye, A. Strumia, and D. Teresi JHEP 08 (2018) 188, [arXiv:1805.01473].
  * [20] I. Baldes, Y. Gouttenoire, and F. Sala JHEP 04 (2021) 278, [arXiv:2007.08440].
  * [21] J. Shelton and K. M. Zurek Phys. Rev. D 82 (2010) 123512, [arXiv:1008.1997].
  * [22] K. Petraki, M. Trodden, and R. R. Volkas JCAP 02 (2012) 044, [arXiv:1111.4786].
  * [23] I. Baldes JCAP 05 (2017) 028, [arXiv:1702.02117].
  * [24] E. Hall, T. Konstandin, R. McGehee, and H. Murayama arXiv:1911.12342.
  * [25] A. Falkowski and J. M. No JHEP 02 (2013) 034, [arXiv:1211.5615].
  * [26] M. J. Baker, J. Kopp, and A. J. Long Phys. Rev. Lett. 125 (Oct, 2020) 151102.
  * [27] D. Chway, T. H. Jung, and C. S. Shin Phys. Rev. D 101 (2020), no. 9 095019, [arXiv:1912.04238].
  * [28] D. Marfatia and P.-Y. Tseng arXiv:2006.07313.
  * [29] A. Azatov and M. Vanvlasselaer JCAP 01 (2021) 058, [arXiv:2010.02590].
  * [30] A. Katz and A. Riotto JCAP 11 (2016) 011, [arXiv:1608.00583].
  * [31] P. Baratella, A. Pomarol, and F. Rompineve JHEP 03 (2019) 100, [arXiv:1812.06996].
  * [32] L. Delle Rose, G. Panico, M. Redi, and A. Tesi JHEP 04 (2020) 025, [arXiv:1912.06139].
  * [33] N. Craig, N. Levi, A. Mariotti, and D. Redigolo JHEP 21 (2020) 184, [arXiv:2011.13949].
  * [34] J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen JCAP 1906 (2019), no. 06 024, [arXiv:1903.09642].
  * [35] A. Azatov and M. Vanvlasselaer JHEP 09 (2020) 085, [arXiv:2003.10265].
  * [36] J. R. Espinosa and M. Quiros Phys. Rev. D 76 (2007) 076004, [hep-ph/0701145].
  * [37] A. Noble and M. Perelstein Phys. Rev. D 78 (2008) 063518, [arXiv:0711.3018].
  * [38] A. Ashoorioon and T. Konstandin JHEP 07 (2009) 086, [arXiv:0904.0353].
  * [39] P. H. Damgaard, D. O’Connell, T. C. Petersen, and A. Tranberg Phys. Rev. Lett. 111 (2013), no. 22 221804, [arXiv:1305.4362].
  * [40] T. Alanne, K. Tuominen, and V. Vaskonen Nucl. Phys. B 889 (2014) 692–711, [arXiv:1407.0688].
  * [41] P. Huang, A. J. Long, and L.-T. Wang Phys. Rev. D 94 (2016), no. 7 075008, [arXiv:1608.06619].
  * [42] A. Beniwal, M. Lewicki, J. D. Wells, M. White, and A. G. Williams JHEP 08 (2017) 108, [arXiv:1702.06124].
  * [43] V. Barger, P. Langacker, M. McCaskey, M. J. Ramsey-Musolf, and G. Shaughnessy Phys. Rev. D 77 (2008) 035005, [arXiv:0706.4311].
  * [44] J. R. Espinosa, T. Konstandin, and F. Riva Nucl. Phys. B 854 (2012) 592–630, [arXiv:1107.5441].
  * [45] J. Kozaczuk, M. J. Ramsey-Musolf, and J. Shelton Phys. Rev. D 101 (2020), no. 11 115035, [arXiv:1911.10210].
  * [46] Z. Wang, X. Zhu, E. E. Khoda, S.-C. Hsu, N. Konstantinidis, K. Li, S. Li, M. J. Ramsey-Musolf, Y. Wu, and Y. E. Zhang, Study of Electroweak Phase Transition in Exotic Higgs Decays at the CEPC, 3, 2022.  arXiv:2203.10184.
  * [47] G. Kurup and M. Perelstein Phys. Rev. D 96 (2017), no. 1 015036, [arXiv:1704.03381].
  * [48] Z. Kang, P. Ko, and T. Matsui JHEP 02 (2018) 115, [arXiv:1706.09721].
  * [49] L. Niemi, P. Schicho, and T. V. I. Tenkanen Phys. Rev. D 103 (2021), no. 11 115035, [arXiv:2103.07467].
  * [50] B. Gripaios, A. Pomarol, F. Riva, and J. Serra JHEP 04 (2009) 070, [arXiv:0902.1483].
  * [51] J. R. Espinosa, B. Gripaios, T. Konstandin, and F. Riva JCAP 01 (2012) 012, [arXiv:1110.2876].
  * [52] S. De Curtis, L. Delle Rose, and G. Panico JHEP 12 (2019) 149, [arXiv:1909.07894].
  * [53] CMS Collaboration, A. M. Sirunyan et al. Phys. Lett. B 805 (2020) 135425, [arXiv:2002.06398].
  * [54] E. J. Weinberg, Radiative corrections as the origin of spontaneous symmetry breaking.  PhD thesis, Harvard U., 1973.  hep-th/0507214.
  * [55] C. Delaunay, C. Grojean, and J. D. Wells JHEP 04 (2008) 029, [arXiv:0711.2511].
  * [56] M. B. Einhorn and D. R. T. Jones Nucl. Phys. B 230 (1984) 261–272.
  * [57] D. Curtin, P. Meade, and H. Ramani Eur. Phys. J. C78 (2018), no. 9 787, [arXiv:1612.00466].
  * [58] 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.
  * [59] S. R. Coleman Phys. Rev. D15 (1977) 2929–2936. [Erratum: Phys. Rev.D16,1248(1977)].
  * [60] A. D. Linde Phys. Lett. 100B (1981) 37–40.
  * [61] A. D. Linde Nucl. Phys. B216 (1983) 421. [Erratum: Nucl. Phys.B223,544(1983)].
  * [62] A. H. Guth and S. H. H. Tye Phys. Rev. Lett. 44 (Apr, 1980) 963–963.
  * [63] A. H. Guth and E. J. Weinberg.
  * [64] D. Bodeker and G. D. Moore JCAP 0905 (2009) 009, [arXiv:0903.4099].
  * [65] M. Dine, R. G. Leigh, P. Y. Huet, A. D. Linde, and D. A. Linde Phys. Rev. D46 (1992) 550–571, [hep-ph/9203203].
  * [66] D. Bodeker and G. D. Moore JCAP 1705 (2017), no. 05 025, [arXiv:1703.08215].
  * [67] Y. Gouttenoire, R. Jinno, and F. Sala arXiv:2112.07686.
  * [68] S. Höche, J. Kozaczuk, A. J. Long, J. Turner, and Y. Wang JCAP 03 (2021) 009, [arXiv:2007.10343].
  * [69] J. M. Cline and K. Kainulainen JCAP 01 (2013) 012, [arXiv:1210.4196].
  * [70] D. Curtin, P. Meade, and C.-T. Yu JHEP 11 (2014) 127, [arXiv:1409.0005].
  * [71] V. Vaskonen Phys. Rev. D 95 (2017), no. 12 123515, [arXiv:1611.02073].
  * [72] D. Buttazzo, F. Sala, and A. Tesi JHEP 11 (2015) 158, [arXiv:1505.05488].
  * [73] S. Profumo, M. J. Ramsey-Musolf, and G. Shaughnessy JHEP 08 (2007) 010, [arXiv:0705.2425].
  * [74] H. H. Patel and M. J. Ramsey-Musolf Phys. Rev. D 88 (2013) 035013, [arXiv:1212.5652].
  * [75] W. Huang, Z. Kang, J. Shu, P. Wu, and J. M. Yang Phys. Rev. D 91 (2015), no. 2 025006, [arXiv:1405.1152].
  * [76] M. Jiang, L. Bian, W. Huang, and J. Shu Phys. Rev. D 93 (2016), no. 6 065032, [arXiv:1502.07574].
  * [77] C.-W. Chiang, M. J. Ramsey-Musolf, and E. Senaha Phys. Rev. D 97 (2018), no. 1 015005, [arXiv:1707.09960].
  * [78] V. Guada, M. Nemevšek, and M. Pintar Comput. Phys. Commun. 256 (2020) 107480, [arXiv:2002.00881].
  * [79] R. Barbieri and G. F. Giudice Nucl. Phys. B 306 (1988) 63–76.
  * [80] J. Ellis, M. Lewicki, and J. M. No arXiv:1809.08242. [JCAP1904,003(2019)].
  * [81] K. Enqvist, J. Ignatius, K. Kajantie, and K. Rummukainen Phys. Rev. D45 (1992) 3415–3428.
  * [82] I. Baldes, Y. Gouttenoire, F. Sala, and G. Servant arXiv:2110.13926.
  * [83] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, J. Camp, C. Caprini, V. Cardoso, M. Colpi, J. Conklin, N. Cornish, C. Cutler, K. Danzmann, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, J. Gair, L. G. Bote, D. Giardini, F. Gibert, C. Grimani, H. Halloin, G. Heinzel, T. Hertog, M. Hewitson, K. Holley-Bockelmann, D. Hollington, M. Hueller, H. Inchauspe, P. Jetzer, N. Karnesis, C. Killow, A. Klein, B. Klipstein, N. Korsakova, S. L. Larson, J. Livas, I. Lloro, N. Man, D. Mance, J. Martino, I. Mateos, K. McKenzie, S. T. McWilliams, C. Miller, G. Mueller, G. Nardini, G. Nelemans, M. Nofrarias, A. Petiteau, P. Pivato, E. Plagnol, E. Porter, J. Reiche, D. Robertson, N. Robertson, E. Rossi, G. Russano, B. Schutz, A. Sesana, D. Shoemaker, J. Slutsky, C. F. Sopuerta, T. Sumner, N. Tamanini, I. Thorpe, M. Troebs, M. Vallisneri, A. Vecchio, D. Vetrugno, S. Vitale, M. Volonteri, G. Wanner, H. Ward, P. Wass, W. Weber, J. Ziemer, and P. Zweifel, Laser interferometer space antenna, 2017.
  * [84] C. Caprini et al. arXiv:1910.13125.
  * [85] C. Caprini et al. JCAP 1604 (2016), no. 04 001, [arXiv:1512.06239].
  * [86] B. Von Harling, A. Pomarol, O. Pujolàs, and F. Rompineve JHEP 04 (2020) 195, [arXiv:1912.07587].
  * [87] V. Brdar, A. J. Helmboldt, and J. Kubo JCAP 02 (2019) 021, [arXiv:1810.12306].
  * [88] V. Corbin and N. J. Cornish Class. Quant. Grav. 23 (2006) 2435–2446, [gr-qc/0512039].
  * [89] J. Crowder and N. J. Cornish Phys. Rev. D 72 (2005) 083005, [gr-qc/0506015].
  * [90] N. Seto, S. Kawamura, and T. Nakamura Phys. Rev. Lett. 87 (2001) 221103, [astro-ph/0108011].
  * [91] K. Yagi and N. Seto Phys. Rev. D 83 (2011) 044011, [arXiv:1101.3940]. [Erratum: Phys.Rev.D 95, 109901 (2017)].
  * [92] S. Isoyama, H. Nakano, and T. Nakamura PTEP 2018 (2018), no. 7 073E01, [arXiv:1802.06977].
  * [93] S. Hild et al. Class. Quant. Grav. 28 (2011) 094013, [arXiv:1012.0908].
  * [94] B. Sathyaprakash et al. Class. Quant. Grav. 29 (2012) 124013, [arXiv:1206.0331]. [Erratum: Class. Quant. Grav.30,079501(2013)].
  * [95] M. Maggiore et al. JCAP 03 (2020) 050, [arXiv:1912.02622].
  * [96] L. Badurina et al. JCAP 05 (2020) 011, [arXiv:1911.11755].
  * [97] AEDGE Collaboration, Y. A. El-Neaj et al. EPJ Quant. Technol. 7 (2020) 6, [arXiv:1908.00802].
  * [98] D. Cutting, M. Hindmarsh, and D. J. Weir Phys. Rev. D97 (2018), no. 12 123513, [arXiv:1802.05712].
  * [99] J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant JCAP 1006 (2010) 028, [arXiv:1004.4187].
  * [100] S. Argyropoulos, O. Brandt, and U. Haisch Symmetry 2021 (9, 2021) 13, [arXiv:2109.13597].
  * [101] S. Bruggisser, B. Von Harling, O. Matsedonskyi, and G. Servant JHEP 12 (2018) 099, [arXiv:1804.07314].
  * [102] S. Blasi and A. Mariotti arXiv:2203.16450.