# Bubble wall velocity during electroweak phase transition in the inert
doublet model

Siyu Jiang  Fa Peng Huang  Corresponding Author<EMAIL_ADDRESS>Xiao Wang  MOE Key Laboratory of TianQin Mission, TianQin Research Center for
Gravitational Physics & School of Physics and Astronomy, Frontiers Science
Center for TianQin, Gravitational Wave Research Center of CNSA, Sun Yat-sen
University (Zhuhai Campus), Zhuhai 519082, China

###### Abstract

The bubble wall velocity is essential for the spectra of phase transition
gravitational wave, electroweak baryogenesis and the dynamical dark matter
mechanism. We perform detailed calculations of the bubble wall velocity in the
well-motivated inert doublet model using the microphysical approach with some
recent new methods.

##  I Introduction

Many fundamental problems, like the electroweak baryogenesis Cline and
Kainulainen (2020), the dark matter (DM) formation mechanism Baker et al.
(2020); Chway et al. (2020); Huang and Li (2017) and the strength or energy
budget Espinosa et al. (2010); Giese et al. (2020, 2021); Wang et al. (2021a);
Leitao and Megevand (2011, 2015) of the phase transition gravitational wave
(GW) strongly depend on the precise value of the bubble wall velocity. For
traditional electroweak baryogenesis, bubble wall velocities lower than the
sound speed are favored to guarantee sufficient diffusion time for the
baryons. The net baryon-to-photon ratio observed by the cosmic microwave
background or big bang nucleosynthesis may be sensitive to the bubble wall
velocity. New mechanisms for baryogenesis with large bubble wall velocity were
recently proposed in Refs. Laurent and Cline (2020); De Curtis et al. (2022);
Laurent and Cline (2022). Thus, to understand the origin of the matter-
antimatter asymmetry of our Universe, precise calculation of the bubble wall
velocity is crucial for a given electroweak baryogenesis model. Recently,
motivated by the current situation of DM detection, a new mechanism is
proposed to produce DM by the bubble dynamics through a strong first-order
phase transition (SFOPT) process in the early Universe Baker et al. (2020);
Chway et al. (2020); Huang and Li (2017). For this new DM mechanism, the
bubble wall velocity is also important as shown in Ref. Huang and Li (2017).
The DM relic density might depends on the bubble wall velocity. Both the
electroweak baryogenesis and the new DM mechanisms could be probed by phase
transition GW signals, whose spectra strongly depend on the bubble wall
velocity. The GW spectra for subsonic bubble wall velocity and supersonic
bubble wall velocity might have several orders of hierarchy. For example, the
energy budget of the phase transition GW is sensitive to the bubble wall
velocity. Therefore, from the perspective of GW experiments, bubble wall
velocity would be the first dynamical parameter of a SFOPT to be confirmed if
the phase transition GW is observed at future GW experiments like LISA Amaro-
Seoane et al. (2017), TianQin Luo et al. (2016); Mei et al. (2021), and Taiji
Hu and Wu (2017).

The bubble wall velocity is very important for both theoretical study and
experimental search. However, the precise calculation of bubble wall velocity
is complicated and difficult. The bubble wall velocity of a given new physics
model is generally determined by the hydrodynamic and microphysical processes
which basically involve the particles scattering processes between the
relevant particles of the new model at the vicinity of the bubble wall.
Besides the complicated hydrodynamic effects, the bubble wall velocity might
be model dependent based on the fact that it is complicated to quantify these
particle scattering processes precisely in the thermal plasma due to the
thermal effects and the infrared behavior. These issues make it difficult to
precisely predict the bubble wall velocity. Most of the previous studies on
electroweak baryogenesis, phase transition GW, and phase transition related DM
mechanism just take the bubble wall velocity as an input parameter. Therefore,
there exists large theoretical uncertainties for these predictions due to the
model-dependent bubble wall velocity. Some pioneering works Dine et al.
(1992); Liu et al. (1992); Ignatius et al. (1994); Moore and Prokopec (1995a,
b); Moore (2000a) show us how to calculate the wall velocity under certain
assumptions. The bubble wall velocity depends on the friction force acting on
the expanding bubbles, which is determined by the deviation of the massive
particle populations from thermal equilibrium. Basically, under semiclassical
approximation, the bubble wall velocity can be obtained by simultaneously
solving the equation of motion for the Higgs field (or the order-parameter
scalar field for the phase transition process) and the Boltzmann equations of
the massive particle species. To solve the Boltzmann equations, it is crucial
to calculate the collision terms, which quantify the particle scattering
processes in the vicinity of bubble wall.

The bubble wall velocity was first calculated microphysically in Refs. Moore
and Prokopec (1995b, a) by Moore and Prokopec for the Standard Model (SM) case
and then in Minimal Supersymmetric Standard Model John and Schmidt (2001). To
avoid their expansive calculations, some works modeled the friction by using
phenomenological method Espinosa et al. (2010); Huber and Sopena (2013);
Mégevand (2013); Megevand and Sanchez (2009) and they also considered the
hydrodynamic effects in detail to the velocity calculation. As the method of
Moore and Prokopec is problematic when the wall velocity is approaching the
speed of sound, Refs. Laurent and Cline (2020); De Curtis et al. (2022);
Lewicki et al. (2022); Laurent and Cline (2022) proposed some new ansatz, and
Refs. Dorsch et al. (2021, 2022) did higher-order calculations. Recently, it
is also interesting to consider bubble wall velocity in local equilibrium
Barroso Mancha et al. (2021); Balaji et al. (2021); Ai et al. (2022); Wang and
Yuwen (2022).

Using the microphysical approach with the recent new methods, for the first
time we calculate the bubble wall velocity in the well-motivated inert doublet
model (IDM), which could help to greatly reduce the large uncertainties in
calculating the phase transition GW and electroweak baryogenesis from the
unknown bubble wall velocity. We show the basic procedure to calculate the
bubble wall velocity in a given new physics model with a SFOPT in Fig. 1.

This paper is organized as follows. We briefly discuss the IDM and choose the
benchmark parameters in Sec. II, and describe the equation of motion for the
Higgs background field in Sec. III. Then we derive the hydrodynamic effects in
the calculation of bubble wall velocity in Sec. IV. In Sec. V, taking the IDM
model as concrete model, we introduce the basic method to calculate the bubble
wall velocity. After that, we calculate the collision terms of the IDM in Sec.
VI. Therefore, we obtain the final bubble wall velocity in Sec. VII. The
conclusion is given in Sec. VIII.

Figure 1:  Schematic procedure to calculate the bubble wall velocity during a
cosmological phase transition.

##  II A case study: inert doublet model

Firstly, we briefly introduce the well-studied IDM, which could improve the
naturalness problem and provide a natural DM candidate Barbieri et al. (2006);
Lopez Honorez and Yaguna (2010). Meanwhile, the IDM could trigger a SFOPT
Chowdhury et al. (2012) in the early Universe around 100 GeV during the
electroweak spontaneously symmetry breaking process. Then, one could use the
phase transition GW to explore the IDM. However, the GW spectra strongly
depend on the bubble wall velocity, which has not been calculated before in
the IDM. In this work, we perform detailed calculations on the bubble wall
velocity. We begin our discussions from the tree-level scalar potential at
zero temperature,

$\displaystyle V_{0}=$
$\displaystyle\mu^{2}_{1}|\Phi|^{2}+\mu^{2}_{2}|\eta|^{2}+\frac{1}{2}\lambda_{1}|\Phi|^{4}+\frac{1}{2}\lambda_{2}|\eta|^{4}$
$\displaystyle+\lambda_{3}|\Phi|^{2}|\eta|^{2}+\lambda_{4}|\Phi^{\dagger}\eta|^{2}+\frac{1}{2}\\{\lambda_{5}(\Phi^{\dagger}\eta)^{2}+\mathrm{H.c.}\\}\,\,,$
(1)

where $\Phi$ is the SM Higgs doublet and $\eta$ represents the inert doublet.
The vacuum stability requires $\lambda_{1}>0$, $\lambda_{2}>0$,
$\sqrt{\lambda_{1}+\lambda_{2}}+\lambda_{3}>0$,
$\lambda_{3}+\lambda_{4}\pm|\lambda_{5}|>0$ Barbieri et al. (2006); Lopez
Honorez and Yaguna (2010). At zero temperature, the two doublet scalar fields
could be written as

$\Phi=\begin{pmatrix}G^{+}\\\ \frac{1}{\sqrt{2}}(h+v+iG^{0})\end{pmatrix},\
\eta=\begin{pmatrix}H^{+}\\\ \frac{1}{\sqrt{2}}(H+iA)\end{pmatrix},$ (2)

where the mass of SM Higgs boson $h$ is 125 GeV and the vacuum expectation
value (VEV) $v=246$ GeV. $G^{\pm}$ and $G^{0}$ are the charged and neutral
Nambu-Goldstone bosons respectively. $H$ and $A$ are the $CP$-even and
$CP$-odd inert scalars respectively. $H^{\pm}$ are the charged inert scalars.
At zero temperature, we show the inert scalar masses beyond the SM

$\displaystyle\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H}^{2}$
$\displaystyle=\mu_{2}^{2}+\frac{1}{2}(\lambda_{3}+\lambda_{4}+\lambda_{5})v^{2},$
(3) $\displaystyle\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{A}^{2}$
$\displaystyle=\mu_{2}^{2}+\frac{1}{2}(\lambda_{3}+\lambda_{4}-\lambda_{5})v^{2},$
(4) $\displaystyle\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}^{2}_{H^{\pm}}$
$\displaystyle=\mu_{2}^{2}+\frac{1}{2}\lambda_{3}v^{2}\,\,.$ (5)

These new inert scalars could modify $T$ parameter with the deviation $\Delta
T$ approximated as

$\displaystyle\Delta T\simeq\frac{1}{6\pi e^{2}v^{2}}(\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H^{\pm}}-\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H})(\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H^{\pm}}-\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{A})\,\,.$
(6)

From the above expression, there exists a simple and obvious parameter region
that $\Delta T\simeq 0$ if we assume $\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{A}^{2}=\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}^{2}_{H^{\pm}}$.
This assumption corresponds to $\lambda_{4}=\lambda_{5}<0,~{}~{}\lambda_{3}>0$
, and could satisfy all the constraints from electroweak precise measurements,
DM direct searches and the collider data Chowdhury et al. (2012). In other
words, we assume degenerate $CP$-odd and charged scalar masses

$\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{A}^{2}=\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}^{2}_{H^{\pm}}=\mu_{2}^{2}+\frac{1}{2}\lambda_{3}v^{2}\,\,.$
(7)

This makes $\hbox to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H}^{2}=\mu_{2}^{2}+\lambda_{L}v^{2}$
the lightest particle, which becomes the natural DM candidate Barbieri et al.
(2006); Lopez Honorez and Yaguna (2010) with the tiny DM-Higgs boson coupling
$\lambda_{L}=(\lambda_{3}+\lambda_{4}+\lambda_{5})/2$. It should be very small
to be consistent with the DM direct search. Using the package micrOMEGAs
Barducci et al. (2018), we include the resonant effects and require the DM
relic abundance Aghanim et al. (2020) $\Omega_{\text{DM}}h^{2}=0.11933\pm
0.00091$. Then we could obtain the constraint Wang et al. (2021b)
$\lambda_{L}\lesssim 0.003$. These constraints almost reach the blind spots of
the IDM, which are difficult for DM direct searches. Future lepton colliders
in synergy with phase transition GW might help to explore the DM blind spots
Huang and Yu (2018); Wang et al. (2021b). When $\lambda_{3}$, $\lambda_{4}$,
$\lambda_{5}$ are $\mathcal{O}(1)$, a SFOPT with associated phase transition
GW could be produced Chowdhury et al. (2012); Borah and Cline (2012); Gil et
al. (2012); Cline and Kainulainen (2013); AbdusSalam and Chowdhury (2014);
Blinov et al. (2015); Cao et al. (2018); Huang and Yu (2018); Laine et al.
(2017); Senaha (2019); Huang and Senaha (2019); Kainulainen et al. (2019);
Wang et al. (2021b).

When we discuss the phase transition dynamics of the SFOPT, we denote the
background of the Higgs field as $\phi$, namely,

$\Phi=\begin{pmatrix}0\\\ \frac{1}{\sqrt{2}}\phi\end{pmatrix}\,\,.$ (8)

Then we could have the following new field-dependent masses beyond the SM,

$\displaystyle m_{H}^{2}(\phi)$
$\displaystyle=\mu_{2}^{2}+\frac{1}{2}(\lambda_{3}+\lambda_{4}+\lambda_{5})\phi^{2},$
(9) $\displaystyle m_{A}^{2}(\phi)$
$\displaystyle=\mu_{2}^{2}+\frac{1}{2}(\lambda_{3}+\lambda_{4}-\lambda_{5})\phi^{2},$
(10) $\displaystyle m^{2}_{H^{\pm}}(\phi)$
$\displaystyle=\mu_{2}^{2}+\frac{1}{2}\lambda_{3}\phi^{2}\,\,.$ (11)

We should obtain the finite-temperature effective potential Niemi et al.
(2021); Schicho et al. (2022) after taking the one-loop quantum and thermal
corrections with daisy resummation into account. Namely,

$V_{\text{eff}}(\phi,T)=V_{0}(\phi)+V_{\text{CW}}({\phi})+V_{\text{T}}(\phi,T)\,\,.$
(12)

The one-loop quantum correction of the potential in the on-shell
renormalizaiton scheme is

$V_{\text{CW}}({\phi},T=0)=\sum_{i}\frac{n_{i}}{64\pi^{2}}\left[{m}_{i}^{4}(\phi)\left(\ln\frac{{m}_{i}^{2}(\phi)}{\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{i}^{2}}-\frac{3}{2}\right)+2\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{i}^{2}\bar{m}_{i}^{2}(\phi)\right],$
(13)

where $n_{i}$ is the degree of freedom for each massive particle
$m_{i}(\phi)$. For the daisy resummation, we use the scheme proposed by Dolan
and Jackiw Dolan and Jackiw (1974) where one only needs to consider the
thermal correction in the thermal function $I_{b}$. Namely, the one-loop
thermal correction $V_{\text{T}}$ should be improved as

$V_{\text{T}}(\phi,T>0)=\sum_{i}n_{i}\frac{T^{4}}{2\pi^{2}}I_{b}\left(\frac{{M}_{i}^{2}}{T^{2}}\right),$
(14)

where ${M}_{i}^{2}={m}_{i}^{2}(\phi)+\Pi_{i}(T)$ and
$I_{b,f}\left(a^{2}\right)=\int_{0}^{\infty}dx~{}x^{2}\ln\left(1\mp
e^{-\sqrt{x^{2}+a^{2}}}\right)$. $\Pi_{i}(T)$ represents the thermal
correction for scalar bosons and the longitudinal components of gauge bosons.
Scalars with $\mu_{1}^{2}$ and $\mu_{2}^{2}$ appearing in their field-
dependent masses can be simply replaced by $\mu_{1}^{2}+\Pi_{\Phi}(T)$ and
$\mu_{2}^{2}+\Pi_{\eta}(T)$, respectively, where

$\displaystyle\Pi_{\Phi}$
$\displaystyle=\frac{T^{2}}{12}\left[3\lambda_{1}+2\lambda_{3}+\lambda_{4}+\frac{3}{4}(3g_{w}^{2}+g_{Y}^{2})+3y_{t}^{2}\right],$
(15) $\displaystyle\Pi_{\eta}$
$\displaystyle=\frac{T^{2}}{12}\left[3\lambda_{2}+2\lambda_{3}+\lambda_{4}+\frac{3}{4}(3g_{w}^{2}+g_{Y}^{2})\right].$
(16)

For the gauge bosons, only their longitudinal components need to include the
thermal corrections in the thermal function $I_{b}$ as the following:

$\displaystyle\Pi_{W}(T)=\Pi_{W}^{\Phi}(T)+\Pi_{W}^{\eta}(T)=\left[\frac{11}{6}+\frac{1}{6}\right]g_{w}^{2}T^{2}=2g_{w}^{2}T^{2}\,\,,$
(17)
$\displaystyle\Pi_{B}(T)=\Pi_{B}^{\Phi}(T)+\Pi_{B}^{\eta}(T)=\left[\frac{11}{6}+\frac{1}{6}\right]g_{Y}^{2}T^{2}=2g_{Y}^{2}T^{2}\,\,,$

where $g_{w}$ and $g_{Y}$ are the gauge coupling of $SU(2)_{L}$ and
$U(1)_{Y}$, respectively. Hence, for the longitudinal components of the gauge
bosons, their physical masses are eigenvalues of the following matrix

${M}_{\text{L}}^{2}=\left(\begin{array}[]{cccc}m_{1}^{2}+\Pi_{W}&0&0&0\\\
0&m_{1}^{2}+\Pi_{W}&0&0\\\ 0&0&m_{1}^{2}+\Pi_{W}&m_{12}^{2}\\\
0&0&m_{12}^{2}&m_{2}^{2}+\Pi_{B}\\\ \end{array}\right)\,\,,$ (18)

where $m_{1}^{2}=g_{w}^{2}\phi^{2}/4$, $m_{2}^{2}=g_{Y}^{2}\phi^{2}/4$ and
$m_{12}^{2}=-g_{w}g_{Y}\phi^{2}/4$.

To focus on the calculation of bubble wall velocity during the SFOPT in the
IDM, we simply choose three sets of benchmark model parameters based on our
numerical calculations where we use the package micrOMEGAs Barducci et al.
(2018) for DM relic abundance, DM direct search, collider constraints Belyaev
et al. (2018), and use CosmoTransitions Wainwright (2012) to calculate the
phase transition dynamics. Taking all the above discussions into
consideration, we choose the following benchmark parameters where $\hbox
to0.0pt{\hskip
1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{h}=125$
GeV, $\lambda_{2}=0.2$. Other parameters are shown in Table 1. For these
benchmark point sets, the DM relic density, the DM direct search, the collider
constraints and condition of the SFOPT could be satisfied simultaneously. We
perform detailed calculations of the bubble wall velocity based on Benchmark A
and give the final results of bubble wall velocity for Benchmarks B and C.

|  $\hbox to0.0pt{\hskip 1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H}$ [GeV] |  $\hbox to0.0pt{\hskip 1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{A}=\hbox to0.0pt{\hskip 1.89006pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$m$}}_{H^{\pm}}$ [GeV] | $\lambda_{L}$ |  $T_{c}$ [GeV] |  $T_{N}$ [GeV]
---|---|---|---|---|---
Benchmark A | 62.66 | 300 | 0.0015 | 118.3 | 117.1
Benchmark B | 65.00 | 300 | 0.0015 | 118.6 | 117.5
Benchmark C | 63.00 | 295 | 0.0015 | 119.4 | 118.4
Table 1: Three sets of benchmark model parameters.

##  III Standard method for solving bubble wall velocity

After bubbles of the broken phase are nucleated, they will expand at a steady
bubble wall velocity under some certain circumstances. To obtain the reliable
bubble wall velocity in a given new physics model, we firstly need to know the
bubble dynamics which are described by the equation of motion (EOM) of the
background field for the SFOPT process. We can obtain the EOM by using energy-
momentum conservation of the scalar-plasma system. The energy-momentum tensor
of the field is

$T^{\mu\nu}_{\phi}=\partial^{\mu}\phi\partial^{\nu}\phi-g^{\mu\nu}\left(\frac{1}{2}\partial_{\alpha}\phi\partial^{\alpha}\phi-
V_{T=0}(\phi)\right)\,\,,$ (19)

where $\phi$ is the background Higgs field in the IDM,
$V_{T=0}(\phi)=V_{0}(\phi)+V_{\rm CW}(\phi)$ is the effective potential at
zero temperature. The energy-momentum tensor of the plasma is

$T_{\mathrm{pl}}^{\mu\nu}=\sum_{i}\int\frac{\mathrm{d}^{3}{p}}{(2\pi)^{3}}\frac{p^{\mu}p^{\nu}}{E_{i}}f_{i}(x,p)\,\,,$
(20)

where $f_{i}(x,p)$ is the distribution function of the particle and
$E_{i}=\sqrt{p^{2}+m_{i}^{2}}$. The sum is over all particle species. Using
energy-momentum conservation condition

$\nabla_{\mu}(T_{\phi}^{\mu\nu}+T_{\rm{pl}}^{\mu\nu})=0\,\,,$ (21)

we can derive the EOM of the background field Moore and Prokopec (1995b, a);
Espinosa et al. (2010); Konstandin et al. (2014),

$\Box\phi+\frac{\partial
V_{T=0}(\phi)}{\partial\phi}+\sum_{i}\frac{dm_{i}^{2}}{d\phi}\int\frac{d^{3}p}{(2\pi)^{3}2E_{i}}f_{i}(x,p)=0\,\,,$
(22)

During a SFOPT, some massive particles will be in an out-of-equilibrium state
near the bubble wall. Hence the distribution function could be approximated as
a thermal equilibrium part plus an out-of-equilibrium part, namely,

$f_{i}\equiv f_{0,i}+\delta f_{i}\,\,.$

The equilibrium distribution function for fermions and bosons in the plasma
frame are given by $f_{0,i}=\frac{1}{\exp(E_{i}/T)\pm 1}$, respectively. The
integral of the equilibrium part of the distribution functions gives the
thermal-correction part of the effective potential $V_{\rm T}(\phi,\text{T})$.
Therefore, the EOM for the background field could be further simplified to

$\Box\phi+\frac{\partial V_{\rm
eff}(\phi,T)}{\partial\phi}+\underbrace{\sum_{i}\frac{dm_{i}^{2}}{d\phi}\int\frac{d^{3}p}{(2\pi)^{3}2E_{i}}\delta
f_{i}(x,p)}_{\mathclap{\text{friction term}}}=0\,\,,$ (23)

where
$V_{\text{eff}}(\phi,T)=V_{0}(\phi)+V_{\text{CW}}({\phi})+V_{\text{T}}(\phi,T)$
is the thermal effective potential in a specific model. The second term is the
driving term that accelerates the bubble wall. The third term in Eq. (23)
behaves as the friction force acting on the bubble wall. We can see that the
contributions mainly come from “heavy” particles with non-negligible field-
dependent mass. “Heavy” means that the particle could get field-dependent mass
comparable with the temperature. Thus, in the IDM, we could only consider the
contributions from $t,W^{\pm},Z,A,H^{\pm}$. For the stationary wall, choosing
the $-z$ as the propagating direction as shown in Fig. 2, all quantities $Q$
are functions of $z+v_{w}t$, so we can set $z+v_{w}t\rightarrow z$ and hence
$\partial_{t}Q\rightarrow v_{w}Q^{\prime}$, $\partial_{z}Q\rightarrow
Q^{\prime}$. $v_{w}$ is the bubble wall velocity and prime means derivative
with respect to $z$. Then in the _plasma frame_ , the EOM with bubble wall
velocity could be further expressed as in Refs. Moore and Prokopec (1995a, b)

$(1-v_{w}^{2})\phi^{\prime\prime}+\frac{\partial V_{\rm
eff}(\phi,T)}{\partial\phi}+\sum_{i}\frac{dm_{i}^{2}}{d\phi}\int\frac{d^{3}p}{(2\pi)^{3}2E_{i}}\delta
f_{i}(x,p)=0\,\,.$ (24)

##  IV Hydrodynamic treatment

Having the above EOM, we begin to consider the hydrodynamic effects where the
temperature might play important roles. For example, Eq. (24) should be
evaluated at a specific temperature. Then we should consider the hydrodynamic
effects which cause the temperature and velocity profiles varying across the
bubble wall Wang et al. (2021a); Espinosa et al. (2010); Leitao and Megevand
(2011, 2015); Giese et al. (2020, 2021). This could be seen from the second
term of Eq. (24) in the wall frame where we multiply it by $\phi^{\prime}$ and
integrate over $z$,

$\displaystyle\int dz\phi^{\prime}\frac{\partial V_{\rm
eff}(\phi,T)}{\partial\phi}$ $\displaystyle=$ $\displaystyle\int
dz\frac{dV_{\rm eff}(\phi,T)}{dz}-\int dz\frac{\partial V_{\rm
eff}(\phi,T)}{\partial T}\frac{dT}{dz}$ (25) $\displaystyle\simeq$
$\displaystyle V_{\rm eff}(\phi_{-},T_{-})-V_{\rm eff}(\phi_{+},T_{+})+\langle
s\rangle(T_{-}-T_{+})\,\,,$ (26)

where $V_{\rm eff}(\phi_{-},T_{-})-V_{\rm eff}(\phi_{+},T_{+})$ is the
pressure difference that push the bubble wall out. Notice that
$\phi_{-}=\phi(z>0)$ and $T_{-}=T(z>0)$ because we let wall move along the
negative $z$ direction as shown in Fig. 2 . Entropy density
$s(\phi,T)\equiv-\partial V_{\rm eff}(\phi,T)/\partial T$ and it gives non-
negligible contribution to the pressure. $\langle
s\rangle\simeq\frac{s(\phi_{+},T_{+})+s(\phi_{-},T_{-})}{2}$ is the average
entropy density. Notice that the presence of $T_{+}$ here in stead of $T_{N}$
is due to the hydrodynamic heating effects in front of the bubble wall.
Besides, as we consider the hydrodynamic effects, temperature in front of the
bubble wall $T_{+}$ is usually different from the one behind the wall $T_{-}$.
Therefore, it is important to calculate the temperature profile around the
bubble wall. All the perturbations are evaluated in front of the wall where
$T=T_{+}$.

The energy-momentum tensor of the equilibrium part of plasma and scalar field
can be combined in the form of ideal fluid,

$T_{f}^{\mu\nu}=(e_{f}+p_{f})u^{\mu}u^{\nu}+p_{f}g^{\mu\nu}=\omega_{f}u^{\mu}u^{\nu}+p_{f}g^{\mu\nu}$
(27)

where $e_{f}$ is the energy density, $p_{f}$ is the pressure and $\omega_{f}$
is the enthalpy of the fluid. $u^{\mu}$ is the fluid four-velocity
$u^{\mu}=\gamma(v)(1,\vec{v})$ in the background plasma frame with the Lorentz
factor $\gamma(v)=1/\sqrt{1-v^{2}}$.

We can write the hydrodynamic equations in the rest frame of Universe Espinosa
et al. (2010),

$\displaystyle\begin{gathered}(\xi-v)\frac{\partial_{\xi}e_{f}}{\omega_{f}}=2\frac{v}{\xi}+\left[1-\gamma^{2}v(\xi-v)\right]\partial_{\xi}v,\\\
(1-v\xi)\frac{\partial_{\xi}p_{f}}{\omega_{f}}=\gamma^{2}(\xi-v)\partial_{\xi}v\,\,,\end{gathered}$
(30)

where $\xi=r/t$ is the self-similar variable of the above equations and has
unit of velocity. $r$ is the distance to the bubble center and $t$ is the time
since the bubble nucleation. “Self-similar” means that there is no
characteristic length or time scale for one steadily expanding bubble.
$v\left(\xi_{w}\right)$ is the fluid velocity at the location of the bubble
wall $\xi_{w}$. Notice that $\xi_{w}=v_{w}$.

We can obtain the differential equation for the temperature from Eq. (30) by
using the following identity

$\frac{\partial p_{f}}{\partial
T}=\partial_{\xi}p_{f}\frac{\partial\xi}{\partial T}\,\,,$ (31)

then the enthalpy can be written as

$\omega_{f}\equiv T\frac{\partial p_{f}}{\partial
T}=T\partial_{\xi}p_{f}\left(\partial_{\xi}T\right)^{-1}\,\,.$ (32)

Substituting it into Eq. (30) one obtains

$\frac{\partial_{\xi}T}{T}=\gamma^{2}\mu\partial_{\xi}v\,\,,$ (33)

where $\mu(\xi,v)=(\xi-v)/(1-\xi v)$ is Lorentz transformation of the
velocity. Given speed of sound in the plasma,
$c_{s}^{2}\equiv(dp_{f}/dT)/(de_{f}/dT)$, from Eq. (30) we can get the
equation of the velocity profile,

$2\frac{v}{\xi}=\gamma^{2}(1-v\xi)\left[\frac{\mu^{2}}{c_{s}^{2}}-1\right]\partial_{\xi}v\,\,.$
(34)

Figure 2: Dynamical process for a bubble wall in the deflagration mode. The
deflagration wall (w) is moving along the negative $z$ direction (as indicated
by arrows) with a shock front (s) propagating in front of it. $v_{+},T_{+}$
are the velocity and temperature in front of the wall, and $v_{-},T_{-}$ are
the velocity and temperature behind the wall. All $v_{-},T_{-},v_{+},T_{+}$
are depicted in the bubble wall frame. $v_{s+},T_{s+}$ are the velocity and
temperature in front of the shock front and $v_{s-},T_{s-}$ the velocity and
temperature behind the shock front. All $v_{s-},T_{s-},v_{s+},T_{s+}$ are
depicted in the shock frame.

In order to solve the hydrodynamic fluid equations, one needs to consider the
boundary conditions in the system by using the conservation of the energy
momentum tensor $\nabla_{\mu}(T_{\phi}^{\mu\nu}+T_{\rm
pl}^{\mu\nu})=\nabla_{\mu}T_{f}^{\mu\nu}=0$ across the bubble wall,

$\omega_{+}\gamma_{+}^{2}v_{+}^{2}+p_{+}=\omega_{-}\gamma_{-}^{2}v_{-}^{2}+p_{-},\quad\omega_{+}\gamma_{+}^{2}v_{+}=\omega_{-}\gamma_{-}^{2}v_{-}\,\,,$
(35)

which can be transformed as

$v_{+}v_{-}=\frac{p_{+}-p_{-}}{e_{+}-e_{-}},\quad\frac{v_{+}}{v_{-}}=\frac{e_{-}+p_{+}}{e_{+}+p_{-}}\,\,,$
(36)

where $+(-)$ means in front of (behind) the bubble wall. These quantities can
be estimated by using the so-called bag model. In the bag model, the pressure
and the energy density can be written as

$\displaystyle\begin{aligned}
&p_{+}=\frac{1}{3}a_{+}T_{+}^{4}-\epsilon_{+},\quad
e_{+}=a_{+}T_{+}^{4}+\epsilon_{+},\quad\text{where}\quad\epsilon_{+}\equiv
V_{T=0}\left(\phi_{+}\right);\\\
&p_{-}=\frac{1}{3}a_{-}T_{-}^{4}-\epsilon_{-},\quad
e_{-}=a_{-}T_{-}^{4}+\epsilon_{-},\quad\text{where}\quad\epsilon_{-}\equiv
V_{T=0}\left(\phi_{-}\right)\text{. }\end{aligned}$ (37)

Notice that these conditions are derived in the rest frame of the bubble wall.
We have pointed out that in order to solve the bubble wall velocity, one
important quantity is the temperature in front of the wall $T_{+}$, where
massive particles meet the bubble wall. Actually, $v_{+},T_{+}$ are important
not only for solving wall velocity but also for calculating electroweak
baryogenesis. In bag model we can see the first two terms in Eq. (26),

$\displaystyle V_{\rm eff}(\phi_{-},T_{-})-V_{\rm eff}(\phi_{+},T_{+})$
$\displaystyle=$ $\displaystyle
p_{+}-p_{-}=\epsilon_{-}-\epsilon_{+}+\frac{1}{3}a_{+}T_{+}^{4}-\frac{1}{3}a_{-}T_{-}^{4}$
(38) $\displaystyle=$
$\displaystyle\epsilon_{-}-\epsilon_{+}+\frac{1}{3}a_{-}(T_{+}^{4}-T_{-}^{4})+\frac{1}{3}(a_{+}-a_{-})T_{+}^{4}\,\,.$

For weak phase transitions $T_{+}\simeq T_{-}$ the third term is approximately
zero and the fourth term accounts for entropy variation which may behave as
friction in _local equilibrium_ as discussed in Refs. Konstandin and No
(2011); Barroso Mancha et al. (2021); Balaji et al. (2021); Ai et al. (2022);
Wang and Yuwen (2022).

Having assumed that the system could be approximated by the bag equation of
state, one could obtain a relation between the plasma velocities in front of
and behind the bubble wall,

$v_{+}=\frac{1}{1+\alpha_{+}}\left[\left(\frac{v_{-}}{2}+\frac{1}{6v_{-}}\right)\pm\sqrt{\left(\frac{v_{-}}{2}+\frac{1}{6v_{-}}\right)^{2}+\alpha_{+}^{2}+\frac{2}{3}\alpha_{+}-\frac{1}{3}}\right]\,\,,$
(39)

where $\alpha_{+}=(\epsilon_{+}-\epsilon_{-})/(a_{+}T_{+}^{4})$ is the phase
transition strength. There would be three possible stable solutions to the
hydrodynamic equations Steinhardt (1982): 1) Deflagration solution
corresponding to the minus sign of above formula has a subsonic wall velocity,
and $v_{+}<v_{-}$. In order to fulfill the requirement that the velocity far
away from the wall should be zero, there should be a shock front in front of
the bubble wall, as we can see in Fig. 2. Note that in this case, the fluid
behind the wall is at rest so we have $v_{-}=v_{w}$; 2) For the detonation
solution, the wall velocity is supersonic and the fluid in front of the wall
is at rest. In this case the solution corresponds to the positive sign of Eq.
(39) and we have $v_{+}=v_{w}$; 3)The hybrid expansion occurs when
$v_{+}<v_{-}=c_{s}$. This is a hybrid state that there should be a shock front
in front of the wall and a rarefaction wave behind the wall.

In this work, we concentrate on deflagration mode where a shock front is
formed in front of the bubble wall as shown in Fig. 2. We can clearly see that
there exist two different interface boundaries, which lead to discontinuous
thermodynamic quantities and two different reference frames. Therefore, we
should carefully mention the designated reference frame when we deal with the
velocity. For convenience, we denote the fluid velocities as measured with
respect to the wall frame as $v_{+}$ if the fluid is in front of the wall and
as $v_{-}$ if the fluid is behind the wall. And we denote the fluid velocities
as measured with respect to the shock frame as $v_{s+}$ if the fluid is in
front of the shock front and as $v_{s-}$ if the fluid is behind the shock
front.

For the shock front, the fluid on both sides of it is in the symmetric phase
and it has the same matching conditions as bubble wall in Eq. (35), then the
velocities on both sides can be written as

$v_{s+}v_{s-}=\frac{1}{3},\quad\frac{v_{s+}}{v_{s-}}=\frac{T_{s+}^{4}+3T_{s-}^{4}}{T_{s-}^{4}+3T_{s+}^{4}}\,\,.$
(40)

We use tildes $\tilde{v}_{\pm}$ to refer to fluid velocities in the fluid
reference frame. This is also the frame of the universe or the reference frame
of the center of the bubble. In this frame the fluid far behind the bubble
wall and far ahead of the shock front is at rest. The fluid velocities in wall
frame or shock frame and fluid frame are related by Lorentz transformations,
i.e.,

$\tilde{v}_{\pm}=\frac{v_{w}-v_{\pm}}{1-v_{w}v_{\pm}}\,\,,\quad\tilde{v}_{s\pm}=\frac{v_{sh}-v_{s\pm}}{1-v_{\omega}v_{s\pm}}\,\,,$
(41)

where $v_{w}$ is just the bubble wall velocity, i.e., the relative velocities
between the wall and the center of the bubble, $v_{sh}$ is the relative
velocity between the shock front and the center of the bubble.

Figure 3: Velocity and temperature profiles in the deflagration mode for
different bubble wall velocity $v_{w}=0.1,0.2,0.3,0.4$. Here we choose
$\alpha_{N}=0.005$. Left panel is the velocity profile and right panel is the
temperature profile.

Since we focus on the subsonic bubble wall velocity, we should give a more
precise definition of a deflagration. In this case the fluid velocity behind
the wall is zero in the fluid frame, i.e,

$\tilde{v}_{\boldsymbol{-}}=0\quad\longrightarrow\quad
v_{w}=v_{\boldsymbol{-}}\,\,.$ (42)

Similarly, we require that the fluid far away should be at rest as in the
presence of the shock front. Then we will get

$\tilde{v}_{s+}=0\quad\longrightarrow\quad v_{sh}=v_{s+}\,\,.$ (43)

In order to obtain the fluid velocity and temperature profile in front of the
bubble wall, we need to solve their corresponding hydrodynamic equations (33)
and (34). Given $\alpha_{+}$ and some value for the wall velocity
$v_{w}<c_{s}$, noting we have $\xi_{w}=v_{w}=v_{-}$, the boundary condition
for the velocity equation in front of bubble wall is

$\tilde{v}_{+}\equiv\frac{\xi_{w}-v_{+}}{1-\xi_{w}v_{+}}=v\left(\xi_{w}\right),$

with $v_{+}$ which is determined by Eq. (39). From Eq. (34), we can see that
the solution would have a singularity when $\mu(\xi,v)=c_{s}$. However, the
solution will reach the shock front position $\xi_{sh}$ before the
singularity. The shock front position satisfies the second boundary condition,

$\tilde{v}_{s-}\equiv\frac{\xi_{sh}-v_{s-}}{1-\xi_{sh}v_{s-}}=v\left(\xi_{sh}\right)\,\,.$
(44)

Notice that $\xi_{sh}=v_{sh}=v_{s+}$. However, we still do not know the
position of the shock front $\xi_{sh}$. One can use the inverse transformation
of Eq. (44) and then substitute it into Eq. (40) which describes the
discontinuities of velocities across the shock front, then we will get

$\mu\left(\xi_{sh},v\left(\xi_{sh}\right)\right)\xi_{sh}=\frac{1}{3}=c_{s}^{2}\,\,,$
(45)

which determines the position of the shock front. Having the velocity profile,
we can substitute it into Eq. (33) to get the temperature profile. The
boundary values are the values of the temperature outside the wall $T_{+}$ and
the temperature inside the shock front $T_{s-}$. In Fig. 3 we show an
illustration of velocity profile and temperature profile where we choose
$\alpha_{N}=0.005$. We could integrate Eq. (33) and get

$\frac{T_{+}}{T_{s-}}=\exp\left[\int_{\xi_{sh}}^{\xi_{w}}d\xi\gamma^{2}\mu\partial_{\xi}v\right]=\exp\left[\int_{\xi_{sh}}^{\xi_{w}}d\xi\frac{2c_{s}^{2}v(\xi-v)}{\xi\left((\xi-v)^{2}-c_{s}^{2}(1-v\xi)^{2}\right)}\right]\,\,.$
(46)

Notice that the temperature in front of the shock front corresponds to the
nucleation temperature of the transition, i.e., $T_{s+}=T_{N}$. Using Eq. (40)
and eliminate $v_{s+}$ we can find

$\frac{T_{s-}^{4}}{T_{N}^{4}}=\frac{T_{s-}^{4}}{T_{s+}^{4}}=\frac{3\left(1-v_{s-}^{2}\right)}{9v_{s-}^{2}-1}\,\,.$
(47)

Using Eqs. (46) and (47) the temperature in front of the bubble wall $T_{+}$
is related to the nucleation temperature $T_{N}$ as

$\frac{T_{+}}{T_{N}}=\left(\frac{3\left(1-v_{s-}^{2}\right)}{9v_{s-}^{2}-1}\right)^{1/4}\exp\left[\int_{\xi_{sh}}^{\xi_{w}}d\xi\frac{2c_{s}^{2}v(\xi-v)}{\xi\left((\xi-v)^{2}-c_{s}^{2}(1-v\xi)^{2}\right)}\right]\,\,.$
(48)

Then we can use

$\frac{a_{+}T_{+}^{4}}{a_{-}T_{-}^{4}}=\frac{\omega_{+}}{\omega_{-}}\,\,,$
(49)

to determine the temperature behind the wall $T_{-}$.

In Fig.4 we show an example of the heating effect with $T_{c}=118.3$ GeV and
$T_{N}=117.1$ GeV. It can be seen that there is one value of $v_{w}$ at which
$T_{+}$ will be larger than $T_{c}$ (However, in reality we think that the
$T_{+}$ could not be larger than $T_{c}$. It could only approach $T_{c}$
because when $T_{+}=T_{c}$ the potential difference is almost zero and the
wall could not be accelerated anymore.). We can give an approximation of the
upper limit of the wall velocity Konstandin and No (2011),

$v_{w}\simeq\left(\frac{{\rm
log}\frac{T_{c}}{T_{N}}}{6\alpha_{c}}\right)^{1/2}.$ (50)

Actually, the heating effect can effectively be viewed as hydrodynamic
backreaction because it effectively reduces the potential difference in front
of and behind the wall Konstandin and No (2011); Wang and Yuwen (2022).

Figure 4: An illustration of temperature variation in front of the wall. We
have $T_{c}=118.3$ GeV and $T_{N}=117.1$ GeV for the benchmark model
parameters in the IDM.

##  V Basic method of calculating bubble wall velocity

Besides heavy particles in the standard model (top quark, $W$ and $Z$ bosons),
in IDM we have extra heavy particle species, the $CP$-odd scalar $A$ and two
charged scalars $H^{\pm}$ of the inert doublet $\eta$. As we can see in Eqs.
(10) and (11), these three scalars have the same coupling such that same
field-dependent mass and same interactions. Therefore, in this work one can
treat them as the same species $A$. For the same reason, since $W^{\pm}$ and
$Z$ have nearly degenerate mass, we can also treat them as the same species
$W$. Other massive species that obtain light field-dependent mass should be
almost in thermal equilibrium and we can treat them as background.

In order to evaluate the friction term, we adopt the flow ansatz. We must know
the out-of-equilibrium distribution part $\delta f$ for each massive particle
population. In principle, the distribution function $f$ for the microscopic
particles are described by the quantum Liouville equation. However, when the
de-Broglie thermal wavelength of particles in the system is smaller than the
bubble wall thickness, the Wentzel-Kramers-Brillouin (WKB) condition $p\gg
1/L$ ($p$ is the momentum of concerned particle and $L$ is the bubble wall
thickness) is satisfied. Then the background field varies slowly and hence the
distribution function for each particle can be approximated by the following
Boltzmann equation

$\frac{d}{dt}f=\left(\frac{\partial}{\partial
t}+\dot{z}\frac{\partial}{\partial z}+\dot{p_{z}}\frac{\partial}{\partial
p_{z}}\right)f=-C[f]\,\,,$ (51)

where $\dot{z}=p_{z}/E$ and
$\dot{p_{z}}=-\partial_{z}E=-(m^{2})^{\prime}/(2E)$. The distribution function
for one species $f$, which deviates from its equilibrium form, can be
expressed as

$f=\frac{1}{e^{(E+\delta)/T}\pm 1}\,\,,$ (52)

where $+$ is for fermions, $-$ is for bosons, and $\delta$ is the
perturbations ($\delta\ll 1$) that can be parametrized as Moore and Prokopec
(1995a, b)

$\delta=-\mu-\mu_{bg}-\frac{E}{T}(\delta T+\delta T_{bg})-p_{z}(\delta
v+\delta v_{bg})\,\,.$ (53)

And we take the background chemical potential perturbation as zero Moore and
Prokopec (1995a, b). We treat particles and antiparticles as one species
neglecting $\mathrm{CP}$ violation.111Noting that this description is not
appropriate for Infrared (IR) excitation with momenta $p\ll T$. These
contributions can be important for boson species Moore (2000b); Kozaczuk
(2015). $C[f]$ is the collision term, we will discuss the detailed
calculations in next section. The collision term and the EOM of the background
field are model dependent. Then the Boltzmann equation can be written as

$\displaystyle\begin{aligned}
&\left(-f_{0}^{\prime}\right)\left(\frac{p_{z}}{E}\left[\partial_{z}\mu+\frac{E}{T}\partial_{z}\left(\delta
T+\delta T_{bg}\right)+p_{z}\partial_{z}\left(\delta v+\delta
v_{bg}\right)\right]+\partial_{t}\mu\right.\\\
&\left.+\frac{E}{T}\partial_{t}\left(\delta T+\delta
T_{bg}\right)+p_{z}\partial_{t}\left(\delta v+\delta
v_{bg}\right)\right)+TC[\mu,\delta T,\delta
v]=\left(-f_{0}^{\prime}\right)\frac{\partial_{t}\left(m^{2}\right)}{2E}\,\,,\end{aligned}$
(54)

where $f_{0}$ is the equilibrium distribution function and
$f_{0}^{\prime}=-\exp(E/T)/[\exp(E/T)\pm 1]^{2}$. It is difficult to solve the
full Boltzmann equation. However, one can truncate the full Boltzmann equation
with three moments for an approximate solution. The three moments are chosen
as $\int d^{3}p/(2\pi)^{3},\int Ed^{3}p/(2\pi)^{2}$, and $\int
p_{z}d^{3}p/(2\pi)^{3}$ Moore and Prokopec (1995a, b). Then, after the
integration and keeping in mind that $\partial_{t}Q\rightarrow
v_{w}Q^{\prime}$, $\partial_{z}Q\rightarrow Q^{\prime}$ for every heavy
species $i=t,W,A$, we obtain

$\displaystyle v_{w}c_{2}^{i}\left(\mu_{i}^{\prime}+\mu_{bg}^{\prime}\right)$
$\displaystyle+v_{w}c_{3}^{i}\left(\delta T_{i}^{\prime}+\delta
T_{bg}^{\prime}\right)+\frac{c_{3}^{i}T}{3}\left(\delta v_{i}^{\prime}+\delta
v_{bg}^{\prime}\right)+\mu_{i}\Gamma_{\mu 1,i}+\delta
T_{i}\Gamma_{T1,i}=\frac{v_{w}c_{1}^{i}}{2T}\left(m_{i}^{2}\right)^{\prime}\,\,,$
$\displaystyle v_{w}c_{3}^{i}\left(\mu_{i}^{\prime}+\mu_{bg}^{\prime}\right)$
$\displaystyle+v_{w}c_{4}^{i}\left(\delta T_{i}^{\prime}+\delta
T_{bg}^{\prime}\right)+\frac{c_{4}^{i}T}{3}\left(\delta v_{i}^{\prime}+\delta
v_{bg}^{\prime}\right)+\mu_{i}\Gamma_{\mu 2,i}+\delta
T_{i}\Gamma_{T2,i}=\frac{v_{w}c_{2}^{i}}{2T}\left(m_{i}^{2}\right)^{\prime}\,\,,$
$\displaystyle\frac{c_{3}^{i}}{3}$
$\displaystyle\left(\mu_{i}^{\prime}+\mu_{bg}^{\prime}\right)+\frac{c_{4}^{i}}{3}\left(\delta
T_{i}^{\prime}+\delta
T_{bg}^{\prime}\right)+\frac{v_{w}c_{4}^{i}T}{3}\left(\delta
v_{i}^{\prime}+\delta v_{bg}^{\prime}\right)+\delta v_{i}T\Gamma_{v,i}=0\,\,,$
(55)

with different integration coefficients of the collision term

$\displaystyle\begin{split}&\int\frac{d^{3}p}{(2\pi)^{3}T^{2}}C[f_{i}]=\mu_{i}\Gamma_{\mu
1,i}+\delta T_{i}~{}\Gamma_{T1,i}\,\,,\\\
&\int\frac{d^{3}p}{(2\pi)^{3}T^{3}}EC[f_{i}]=\mu_{i}\Gamma_{\mu 2,i}+\delta
T_{i}~{}\Gamma_{T2,i}\,\,,\\\
&\int\frac{d^{3}p}{(2\pi)^{3}T^{3}}p_{z}C[f_{i}]=T\delta
v_{i}~{}\Gamma_{v,i}\,\,.\end{split}$ (56)

In Eqs. (V), the constants $c_{j}^{b/f}$ are defined as

$c_{j}^{b/f}T^{j+1}=\int
E^{j-2}(-f_{0}^{\prime})\frac{d^{3}p}{(2\pi)^{3}}\,\,.$ (57)

At lowest order in $m/T$, we have

$c_{1}^{b}=\frac{\log(2T/m_{b})}{2\pi^{2}}\,\,,\quad
c_{2}^{b}=\frac{1}{6}\,\,,\quad c_{3}^{b}=\frac{3\zeta(3)}{\pi^{2}}\,\,,\quad
c_{4}^{b}=\frac{2\pi^{2}}{15}\,\,,$ (58)

for bosons, whereas for fermions we have

$c_{1}^{f}=\frac{\log(2)}{2\pi^{2}}\,\,,\quad c_{2}^{f}=\frac{1}{12}\,\,,\quad
c_{3}^{f}=\frac{9\zeta(3)}{4\pi^{2}}\,\,,\quad
c_{4}^{f}=\frac{7\pi^{2}}{60}\,\,.$ (59)

$\zeta(x)$ is the Riemann zeta function.

Eqs. (V) are not enough to solve for six quantities $\mu$, $\delta T$, $\delta
v$, $\mu_{bg}$, $\delta T_{bg}$, $\delta v_{bg}$. Therefore, we need three
more equations which describe the evolution of the background fluid Moore and
Prokopec (1995b). The background fluid is the sum of all light particles which
have negligible mass variation across the bubble wall. All light particle
species are treated as being at the same temperature $T+\delta T_{bg}$ and
same velocity $\delta v_{bg}$ and the annihilation rates are fast enough that
the background is in chemical equilibrium $\mu_{bg}=0$. Then we obtain three
background equations which have a similar form as the heavy particles in Eqs.
(V):

$\tilde{c}_{4}\left(v_{w}\delta T_{bg}^{\prime}+\frac{\delta
v_{bg}^{\prime}}{3}T\right)=N_{t}(\mu_{t}\Gamma_{\mu 2,t}+\delta
T_{t}\Gamma_{T2,t})+\sum_{\rm bosons}N_{b}(\mu_{b}\Gamma_{\mu 2,b}+\delta
T_{b}\Gamma_{T2,b})\,\,,$ $\frac{\tilde{c}_{4}}{3}(\delta
T_{bg}^{\prime}+v_{w}T\delta v_{bg}^{\prime})=N_{t}T\delta
v_{t}\Gamma_{v,t}+\sum_{\rm bosons}N_{b}T\delta
v_{b}\Gamma_{v,b}\,\,,\quad\mu_{bg}=0\,\,,$ (60)

where we have used the fact that the collision term between massive particle
species and light particle species in the fluid equations with opposite sign.
Here the heat capacity of the light degrees of freedom
$\tilde{c}_{4}=78c_{4}^{f}+19c_{4}^{b}$.

From Eqs. (60) the derivative of the background temperature and velocity can
be written as:

$\displaystyle\begin{split}&\delta T_{bg}^{\prime}=\frac{-v_{w}(\hbox
to0.0pt{\hskip
1.25pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$A$}}+\hbox
to0.0pt{\hskip
1.5434pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$B$}})+(\hbox
to0.0pt{\hskip
1.43123pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$C$}}+\hbox
to0.0pt{\hskip
1.77846pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$D$}})}{\tilde{c}_{4}\left(1/3-v_{w}^{2}\right)}\,\,,\\\
&\delta v_{bg}^{\prime}=\frac{-3v_{w}(\hbox to0.0pt{\hskip
1.43123pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$C$}}+\hbox
to0.0pt{\hskip
1.77846pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$D$}})+(\hbox
to0.0pt{\hskip
1.25pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$A$}}+\hbox
to0.0pt{\hskip
1.5434pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$B$}})}{T\tilde{c}_{4}\left(1/3-v_{w}^{2}\right)}\,\,,\end{split}$
(61)

where

$\displaystyle\begin{split}&\hbox to0.0pt{\hskip
1.25pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$A$}}=N_{t}(\mu_{t}\Gamma_{\mu
2,t}+\delta T_{t}\Gamma_{T2,t})\,\,,\\\ &\hbox to0.0pt{\hskip
1.5434pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$B$}}=\sum_{\rm
bosons}N_{b}(\mu_{b}\Gamma_{\mu 2,b}+\delta T_{b}\Gamma_{T2,b})\,\,,\\\ &\hbox
to0.0pt{\hskip
1.43123pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$C$}}=N_{t}T\delta
v_{t}\Gamma_{v,t}\,\,,\\\ &\hbox to0.0pt{\hskip
1.77846pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$D$}}=\sum_{\rm
bosons}N_{b}T\delta v_{b}\Gamma_{v,b}\,\,.\end{split}$ (62)

After substituting Eqs. (LABEL:background) and $\mu_{bg}=0$ into Eqs. (V), one
can simplify the Boltzmann equation into the following compact matrix form

$\hat{A}\delta^{\prime}+\Gamma\delta=\Sigma\,\,,$ (63)

in our case the vector of perturbations comprise top quark $t$, vector bosons
$W$ and scalar $A$,

$\mathbf{\delta}=(\mu_{t},\delta T_{t},T\delta v_{t},\mu_{W},\delta
T_{W},T\delta v_{W},\mu_{A},\delta T_{A},T\delta v_{A})\,\,,$ (64)

the source term

$\Sigma=\frac{v_{w}}{2T}\left(c_{1}^{t}(m_{t}^{2})^{\prime},c_{2}^{t}(m_{t}^{2})^{\prime},0,c_{1}^{W}(m_{W}^{2})^{\prime},c_{2}^{W}(m_{W}^{2})^{\prime},0,c_{1}^{A}(m_{A}^{2})^{\prime},c_{2}^{A}(m_{A}^{2})^{\prime},0\right)\,\,,$
(65)

and

$\Gamma=\Gamma_{0}+\frac{1}{\tilde{c}_{4}}\mathbb{M}$ (66)
$\hat{A}=\left(\begin{array}[]{ccc}\hat{A}_{t}&0&0\\\ 0&\hat{A}_{W}&0\\\
0&0&\hat{A}_{A}\end{array}\right),\quad\text{where}\quad\hat{A}_{i}=\left(\begin{array}[]{ccc}v_{w}c_{2}^{i}&v_{w}c_{3}^{i}&\frac{1}{3}c_{3}^{i}\\\
v_{w}c_{3}^{i}&v_{w}c_{4}^{i}&\frac{1}{3}c_{4}^{i}\\\
\frac{1}{3}c_{3}^{i}&\frac{1}{3}c_{4}^{i}&\frac{1}{3}v_{w}c_{4}^{i}\\\
\end{array}\right),$ (67)
$\Gamma_{0}=\left(\begin{array}[]{ccc}\Gamma_{t}&0&0\\\ 0&\Gamma_{W}&0\\\
0&0&\Gamma_{A}\end{array}\right),\quad\text{where}\quad\Gamma_{i}=\left(\begin{array}[]{ccc}\Gamma_{\mu
1,i}&\Gamma_{T1,i}&0\\\ \Gamma_{\mu 2,i}&\Gamma_{T2,i}&0\\\
0&0&\Gamma_{v,i}\end{array}\right)\,\,.$ (68)

Here $\mathbb{M}$ is

$\mathbb{M}=\left(\begin{array}[]{ccc}M_{tt}&M_{tW}&M_{tA}\\\
M_{Wt}&M_{WW}&M_{WA}\\\
M_{At}&M_{AW}&M_{AA}\end{array}\right),\quad\text{where}\quad
M_{ij}=N_{j}\left(\begin{array}[]{ccc}c_{3}^{i}\Gamma_{\mu
2,j}&c_{3}^{i}\Gamma_{T2,j}&0\\\ c_{4}^{i}\Gamma_{\mu
2,j}&c_{4}^{i}\Gamma_{T2,j}&0\\\
0&0&c_{4}^{i}\Gamma_{v,j}\end{array}\right)\,\,,$ (69)

where $N_{j}$ is the total degree of freedom for the massive species
($N_{t}=12$, $N_{W}=9$ and $N_{A}=3$).

We should note that when $m_{b}\gg T$ the approximation of $c_{1}^{b}$ in Eqs.
(58) break down and we have

$c_{1}^{b}=(m_{b}/T)^{1/2}{\rm exp}(-m_{b}/T)/(2\pi)^{3/2}\,\,,$ (70)

so that the contribution of the heavy boson particles to the right-hand side
term in Eq. (65) is suppressed by a Boltzmann factor. Their contribution to
friction would then be suppressed. For an extremely strong phase transition,
this Boltzmann factor should be taken into account and tends to increase the
predicted velocity.

##  VI Collision terms

The collision terms are very important for calculating bubble wall velocity.
If they are large enough, we expect that massive particles are only slightly
apart from equilibrium and that they can not produce enough friction to
prevent the bubble expansion. Then the wall velocity may be very large. On the
other hand, if the collision terms are too small, our method that introduce
perturbations may be problematic; thus, it is essential to calculate collision
terms accurately.

In most of the previous studies, only two massive species are considered;
namely, the top quark and the $W$ bosons. Here, in our work we include extra
massive particle species, the new scalars $H^{\pm}$ and $A$. We treat them as
same species when calculating wall velocity so we can only calculate their
collision terms once.

The collision term of the Boltzmann equation for species $i$ can be expressed
as

$\displaystyle\begin{aligned} C[f_{i}]=\frac{1}{2\hbox to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{i}}\sum\frac{1}{2E_{p}}\int\frac{d^{3}kd^{3}p^{\prime}d^{3}k^{\prime}}{(2\pi)^{9}2E_{k}2E_{p^{\prime}}2E_{k^{\prime}}}&\left|\mathcal{M}_{i}\left(p,k;p^{\prime},k^{\prime}\right)\right|^{2}(2\pi)^{4}\\\
&\times\delta^{4}\left(p+k-p^{\prime}-k^{\prime}\right)\mathcal{P}\left[f_{i}(p)\right]\,\,,\end{aligned}$
(71)

where the sum is over all 4-body processes. The matrix elements are summed
over helicities and colors of all four external quasiparticles, then divided
by the number of degrees of freedom $\hbox to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{i}$
corresponding to species $i$,

$\displaystyle\hbox to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{t}=\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{\bar{t}}=6,\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{W^{+}}=\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{W^{-}}=\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{Z}=3,\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{A}=\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{H^{+}}=\hbox
to0.0pt{\hskip
2.06247pt\leavevmode\hbox{$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{$N$}}_{H^{-}}=1\,\,.$
(72)

$p$ is the momentum of the concerned heavy particle species, $k$ represents
the momentum of the other incoming particle. $p^{\prime}$ and $k^{\prime}$
denote the momenta of the outgoing particles. Then the Mandelstam variables
$s$, $t$, $u$ are defined as $s=(p+k)^{2}=(p^{\prime}+k^{\prime})^{2}$,
$t=(p-p^{\prime})^{2}=(k-k^{\prime})^{2}$, and
$u=(p-k^{\prime})^{2}=(k-p^{\prime})^{2}$. The population factor for process
$ij\rightarrow mn$ is

$\mathcal{P}\left[f_{i}(p)\right]\equiv f_{i}(p)f_{j}(k)\left(1\pm
f_{m}(p^{\prime})\right)\left(1\pm
f_{n}(k^{\prime})\right)-f_{m}(p^{\prime})f_{n}(k^{\prime})\left(1\pm
f_{i}(p)\right)\left(1\pm f_{j}(k)\right)$ (73)

with the upper (lower) signs corresponding to bosons (fermions) and $f_{i}$
the appropriate perturbed distribution function for particle $i$, which we
assume to take the form

$f_{i}=\left(e^{\left(E+\delta_{i}\right)/T}\pm 1\right)^{-1}.$ (74)

We now analyze the distribution function at first order. We can write
$f_{i}=1/(\exp a_{i}\pm 1)$ with $a_{i}=(E_{i}-\delta_{i})/T$ in the plasma
frame such that $1\pm f_{i}=f_{i}\exp a_{i}$. Then Eq. (73) can be expressed
as

$\mathcal{P}[f_{i}]=(e^{a_{m}+a_{n}}-e^{a_{i}+a_{j}})f_{i}f_{j}f_{m}f_{n}\,\,,$
(75)

and we have

$\displaystyle\begin{split}\exp(a_{i}+a_{j})&=\exp\left[(E_{i}+E_{j})/T\right]\times\exp(-\delta_{i}/T-\delta_{j}/T)\\\
&\simeq\exp\left[(E_{i}+E_{j})/T\right]\times(1-\delta_{i}/T-\delta_{j}/T)\,\,.\end{split}$
(76)

Then to the first order in $\delta$, we have

$\mathcal{P}[f_{i}]\simeq\frac{\delta_{i}+\delta_{j}-\delta_{m}-\delta_{n}}{T}f_{0,i}f_{0,j}(1\pm
f_{0,m})(1\pm f_{0,n})\,\,,$ (77)

where $f_{0,i}=(e^{E_{i}/T}\pm 1)^{-1}$ is the equilibrium distribution for
species $i$ and we use the fact that
$f_{i}e^{E_{i}/T}=f_{0,i}e^{E_{i}/T}+\mathcal{O}\left(\delta_{i}\right)=1\pm
f_{0,i}+\mathcal{O}\left(\delta_{i}\right)$. After integration we get the
collision terms with the forms in Eqs. (56).

Then we need to consider the dominant scattering processes of the massive
particle species. We calculate the collision terms by using the _leading-log_
approximation Moore and Prokopec (1995a): a) Neglecting masses of all the
external particles; b) Neglecting $s$-channel contributions because they are
not logarithmic; c) The logarithmic IR divergences are regularized by the the
thermal mass of the mediator. We use propagators of the forms
$1/(t-m_{i,T}^{2})$ or $1/(u-m_{i,T}^{2})$, where $m_{i,T}$ is thermal mass of
the mediator. The thermal masses for quarks and gluons are
$m_{q,T}^{2}=g_{s}^{2}T^{2}/6$ and $m_{g,T}^{2}=2g_{s}^{2}T^{2}$,
respectively. $g_{s}$ is the strong coupling constant. The thermal mass of $W$
boson is $m_{W,T}^{2}=11g_{w}^{2}T^{2}/6$. And for the new scalars,
$m_{A,T}^{2}=\lambda_{3}T^{2}/24$. To see where the “leading-log” comes from,
we show an example of the $t$-channel annihilation process
$t\bar{t}\rightarrow gg$ with matrix element
$-(64/9)g_{s}^{4}st/(t-m_{t,T}^{2})^{2}$. The integral about $p^{\prime}$ and
$k^{\prime}$ gives (in the center-of-mass frame),

$\displaystyle\int\frac{p^{\prime
2}dp^{\prime}d\Omega_{p^{\prime}}}{(2\pi)^{3}4E_{p^{\prime}}E_{k^{\prime}}}2\pi\delta\left(2E_{p}-2E_{p^{\prime}}\right)\frac{(2p\cdot
k)2pp^{\prime}\left(1-\cos\theta^{\prime}\right)}{\left[2pp^{\prime}\left(1-\cos\theta^{\prime}\right)+m_{t,T}^{2}\right]^{2}}\simeq\frac{1}{8\pi}\log\frac{2p\cdot
k}{m_{t,T}^{2}}\,\,,$ (78)

where we have used the leading-log approximation that
$E_{k^{\prime}}=E_{p^{\prime}}=\left|\vec{p}^{\prime}\right|=|\vec{k}^{\prime}|=p^{\prime}\text{
and }E_{k}=E_{p}=|\vec{p}|=|\vec{k}|=p$. Here, $2p\cdot k=2pk(1-$
$\cos\theta)+O\left(m_{t,T}^{2}\right)$, with $\theta$ being the plasma-frame
angle between $\vec{p}$ and $\vec{k}$. Then the remaining integrals about
$\theta$ in the plasma frame contain

$\displaystyle\int
d\cos\theta\frac{1}{2}\log\left(\frac{2pk(1-\cos\theta)}{m_{t,T}^{2}}\right)=-1+\log\frac{4pk}{m_{t,T}^{2}}\,\,.$
(79)

The remaining integrals can be done numerically.

However, in order to get more accurate results, we do the numerical
integration directly like $\int\frac{d^{3}p}{(2\pi)^{3}T^{2}}C[f_{i}]$ by
using the phase space parametrization as discussed in Ref. John and Schmidt
(2001). The Monte Carlo package VEGAS Hahn (2005) is used in our numerical
integration calculations. Then we can use Eqs. (56) to extract $\Gamma_{\mu
1}$, $\Gamma_{T1}$, $\Gamma_{\mu 2}$, $\Gamma_{T2}$ and $\Gamma_{v}$ for each
particle species. We use another approximation in that we neglect the
collisions between different massive species since we expect their
contributions are subdominant compared with other processes. This is implied
in Eq. (68). The results of the matrix elements are shown in Table 2. We only
include processes at order of $g_{s}^{4}$, $y_{t}^{2}g_{s}^{2}$,
$g_{w}^{2}g_{s}^{2}$, $g_{w}^{4}$, $\lambda_{3}^{4}$ as the dominant
contributions. The matrix elements are summed over the helicities and colors
of all four external states.

We perform numerical integration and get, for the top quarks,

$\displaystyle\begin{aligned} &\Gamma_{\mu 1,t}\simeq\left(5.0\times
10^{-4}g_{s}^{4}+5.8\times 10^{-4}g_{s}^{2}y_{t}^{2}\right)T\,\,,\\\
&\Gamma_{T1,t}\simeq\Gamma_{\mu 2,t}\simeq\left(1.1\times
10^{-3}g_{s}^{4}+1.3\times 10^{-3}g_{s}^{2}y_{t}^{2}\right)T\,\,,\\\
&\Gamma_{T2,t}\simeq\left(1.1\times 10^{-2}g_{s}^{4}+4.0\times
10^{-3}g_{s}^{2}y_{t}^{2}\right)T\,\,,\\\ &\Gamma_{v,t}\simeq\left(2.0\times
10^{-2}g_{s}^{4}+1.8\times
10^{-3}g_{s}^{2}y_{t}^{2}\right)T\,\,,\end{aligned}$ (80)

for the W bosons,

$\displaystyle\begin{aligned} &\Gamma_{\mu 1,W}\simeq\left(2.3\times
10^{-3}g_{s}^{2}g_{w}^{2}+2.0\times 10^{-3}g_{w}^{4}\right)T\,\,,\\\
&\Gamma_{T1,W}\simeq\Gamma_{\mu 2,W}\simeq\left(4.7\times
10^{-3}g_{s}^{2}g_{w}^{2}+4.1\times 10^{-3}g_{w}^{4}\right)T\,\,,\\\
&\Gamma_{T2,W}\simeq\left(1.5\times 10^{-2}g_{s}^{2}g_{w}^{2}+1.5\times
10^{-2}g_{w}^{4}\right)T\,\,,\\\ &\Gamma_{v,W}\simeq\left(5.7\times
10^{-2}g_{s}^{2}g_{w}^{2}+1.5\times
10^{-2}g_{w}^{4}\right)T\,\,,\end{aligned}$ (81)

and for the new bosons $A$,

$\displaystyle\begin{aligned} &\Gamma_{\mu 1,A}\simeq 1.0\times
10^{-2}\lambda_{3}^{4}T\,\,,\\\ &\Gamma_{T1,A}\simeq\Gamma_{\mu 2,A}\simeq
4.9\times 10^{-3}\lambda_{3}^{4}T\,\,,\\\ &\Gamma_{T2,A}\simeq 5.1\times
10^{-3}\lambda_{3}^{4}T\,\,,\\\ &\Gamma_{v,A}\simeq 1.8\times
10^{-3}\lambda_{3}^{4}T\,\,.\end{aligned}$ (82)

Process | $|\mathcal{M}_{i}|^{2}$
---|---
$\mathcal{O}\left(g_{s}^{4}\right):$ |
$t\bar{t}\leftrightarrow gg:$ | $\frac{128}{3}g_{s}^{4}\left(\frac{u}{t-m_{t,T}^{2}}+\frac{t}{u-m_{t,T}^{2}}\right)$
$tg\leftrightarrow tg:$ | $-\frac{128}{3}g_{s}^{4}\frac{su}{\left(u-m_{g,T}^{2}\right)^{2}}+96g_{s}^{4}\frac{s^{2}+u^{2}}{\left(t-m_{t,T}^{2}\right)^{2}}$
$tq(\bar{q})\leftrightarrow tq(\bar{q}):$ | $160g_{s}^{4}\frac{u^{2}+s^{2}}{\left(t-m_{t,T}^{2}\right)^{2}}$
$\mathcal{O}\left(y_{t}^{2}g_{s}^{2}\right):$ |
$t\bar{t}\leftrightarrow hg,G^{0}g:$ | $8y_{t}^{2}g_{s}^{2}\left(\frac{u}{t-m_{t,T}^{2}}+\frac{t}{u-m_{t,T}^{2}}\right)$
$t\bar{b}\leftrightarrow hG^{+}:$ | $8y_{t}^{2}g_{s}^{2}\left(\frac{u}{t-m_{t,T}^{2}}+\frac{t}{u-m_{b,T}^{2}}\right)$
$tg\leftrightarrow th,tG^{0}:$ | $-8y_{t}^{2}g_{s}^{2}\frac{s}{t-m_{t,T}^{2}}$
$tg\leftrightarrow bG^{+}:$ | $-8y_{t}^{2}g_{s}^{2}\frac{s}{t-m_{b,T}^{2}}$
$tG^{-}\leftrightarrow bg:$ | $-8y_{t}^{2}g_{s}^{2}\frac{s}{t-m_{t,T}^{2}}$
$\mathcal{O}\left(g_{w}^{2}g_{s}^{2}\right):$ |
$Wq\leftrightarrow qg:$ | $-72g_{s}^{2}g_{w}^{2}\frac{s}{t-m_{q,T}^{2}}$
$Wg\leftrightarrow q\bar{q}:$ | $-72g_{s}^{2}g_{w}^{2}\frac{s}{t-m_{q,T}^{2}}$
$\mathcal{O}\left(g_{w}^{4}\right):$ |
$WW\leftrightarrow f\bar{f}:$ | $-\frac{27}{2}g_{w}^{4}\left(\frac{3s}{t-m_{q,T}^{2}}+\frac{s}{t-m_{l,T}^{2}}\right)$
$Wf\leftrightarrow Wf:$ | $360g_{w}^{4}\frac{u^{2}}{\left(t-m_{W,T}^{2}\right)^{2}}-\frac{27}{2}g_{w}^{4}\left(\frac{3s}{u-m_{q,T}^{2}}+\frac{s}{u-m_{l,T}^{2}}\right)$
$\mathcal{O}\left(\lambda_{3}^{4}\right):$ |
$AA\leftrightarrow hh:$ | $\frac{\lambda_{3}^{4}v^{4}}{2}\left[\frac{1}{\left(t-m_{A,T}^{2}\right)^{2}}+\frac{1}{\left(u-m_{A,T}^{2}\right)^{2}}\right]$
$Ah\leftrightarrow hA:$ | $\frac{\lambda_{3}^{4}v^{4}}{2}\frac{1}{\left(t-m_{A,T}^{2}\right)^{2}}$
Table 2:  Relevant 4-body processes and their corresponding matrix elements
for massive particles in the leading-log approximation. The matrix elements
are summed over the helicities and colors of all four external states. We only
include processes at the order of $g_{s}^{4}$, $y_{t}^{2}g_{s}^{2}$,
$g_{w}^{2}g_{s}^{2}$, $g_{w}^{4}$, $\lambda_{3}^{4}$ as the dominant
contributions. Here $q$, $b$ and $l$ represent quarks, bottom quark and
leptons respectively.

Figure 5: Chemical, temperature and velocity perturbations of $W$, $t$ and $A$
for $v_{w}=0.1$ and $L=0.1~{}\rm GeV^{-1}$. Left: perturbations for top quark
(solid lines) and W boson (dashed lines). Right: perturbations for new
scalars.

##  VII Bubble wall velocity

Given the above results, the perturbations could be derived from Eq. (63) by
using Green’s function method,

$(\hat{A}^{-1}\Gamma)_{ij}\chi_{jk}=\chi_{jk}\rho_{k}\,\,,$ (83)

where $\rho_{k}$ are the eigenvalues of $\hat{A}^{-1}\Gamma$ and $\chi$ is the
matrix constituting of the eigenvectors. Note that there is no sum on $k$. It
is then straightforward to write down the Green’s function

$G_{i}(z,y)=\mathrm{sgn}(\rho_{i})e^{-\rho_{i}(z-y)}\Theta[\mathrm{sgn}(\rho_{i})(z-y)].$
(84)

The Heaviside function means the boundary condition
$\delta_{i}(z\rightarrow\pm\infty)=0$. Then the perturbation $\delta_{i}$ will
be given by

$\displaystyle\delta_{i}(z)=\chi_{ij}\int_{-\infty}^{\infty}dyG_{j}(z,y)[\chi^{-1}\hat{A}^{-1}\Sigma(y)]_{j}\,\,,$
(85)

where $\Sigma$ is given by Eq. (65). However, in order to evaluate this
integration, we must know the explicit form of $m_{i}^{2}(\phi(z))$. To
proceed further, we choose the ansatz of the bubble wall profile as

$\phi(z)=\frac{\phi_{-}}{2}\left(1+\tanh\frac{z}{L}\right)\,\,,$ (86)

where $\phi_{-}$ is the VEV of the Higgs boson in the broken phase at $T_{-}$,
$\phi_{-}=\phi(T_{-})$. $L$ is the bubble wall thickness. This bubble wall
profile might work well for relatively weak first-order phase
transition.222For a SFOPT or the ultra supercooling phase transition, this
profile ansatz might not be appropriate Wang et al. (2020a). In this work, for
the chosen benchmark parameters in the IDM, we just have a relatively weak
first-order phase transition. Thus, it is reasonable to assume this shape in
the following analysis.

Numerically solving the truncated Boltzmann equation, we could obtain the
perturbations of all massive species as shown in Fig. 5 and background
perturbations in left panel of Fig. 6, where we have chosen wall velocity
$v_{w}=0.1$ and wall thickness $L=0.1~{}\rm GeV^{-1}$. We can give brief
discussions on the underlying physics of these numerical results. One can see
that the values of perturbations for $t$ and $W$ are larger than $A$’s. This
is because the magnitudes of the collision rates of $A$ are larger than the
rates of $t$ and $W$, such that $A$ is closer to equilibrium than $t$ and $W$.
Also, for $t$ and $W$, their collision rates with respect to the chemical
potential are smaller than those with respect to temperature, so that $t$ and
$W$ are easier to deviate from chemical equilibrium than thermal equilibrium.
Then the extent of their decoupling from chemical equilibrium should be larger
than thermal equilibrium. Such that in Fig. 5, we can see that the
perturbations of the chemical potential are much larger than temperature and
velocity. It is the same for $A$, whose particle number-changing collision
rate is larger than temperature-exchanging rate so the magnitude of chemical
perturbations of $A$ is smaller than temperature and velocity. As the wall is
thick enough respect to the de-Broglie wavelength of massive particles, we
expect that the shape of the dominant perturbations $\mu/T$ for $t$ and $W$
and $\delta T/T$ for $A$ will be closely proportional to $mm^{\prime}$ and
$\delta T_{bg}/T$ be proportional to $m^{2}$, which is consistent with our
results.

Substituting the results obtained from Boltzmann equation under semiclassical
and fluid approximations into Eq. (24) we get that at leading-order
perturbations, the EOM of Higgs can be approximated as

$\displaystyle S_{\mathrm{EOM}}\equiv(1$
$\displaystyle\left.-v_{w}^{2}\right)\phi^{\prime\prime}+\frac{\partial
V_{\mathrm{eff}}\left(\phi,T_{+}\right)}{\partial\phi}+\frac{N_{t}T_{+}}{2}\frac{dm_{t}^{2}}{d\phi}\times\left(c_{1}^{t}\mu_{t}+c_{2}^{t}\left(\delta
T_{t}+\delta T_{bg}\right)\right)$ (87)
$\displaystyle+\sum_{b}\frac{N_{b}T_{+}}{2}\frac{dm_{b}^{2}}{d\phi}\left(c_{1}^{b}\mu_{b}+c_{2}^{b}\left(\delta
T_{b}+\delta T_{bg}\right)\right)=0\,\,,$

where $N_{i}$ is the degree of freedom of particle $i$. Here $T_{+}$ is
temperature just in front of the bubble wall and can be solved with the
hydrodynamic treatment of expanding bubble which we have discussed previously.

Given specific bubble wall profile in Eq. (86), it is still difficult to fully
solve the EOM in Eq. (87). In practice, we could obtain an approximate
solution when the following two constraints are satisfied Konstandin et al.
(2014)

$M_{1}=\int S_{\mathrm{EOM}}\phi^{\prime}dz=0,\quad M_{2}=\int
S_{\mathrm{EOM}}(2\phi-\phi_{-})\phi^{\prime}dz=0\,\,.$ (88)

The first equation $M_{1}=0$ in Eq. (88) means the total pressure on the
bubble wall vanishes in the steady velocity regime. And the second equation
$M_{2}=0$ indicates the bubble wall thickness should not change anymore, which
means the total pressure gradient should be zero for the steady bubble wall
velocity.

The equilibrium, constant temperature part of Eq. (88) can be easily
performed,

$\displaystyle\int\left[(1-v_{w}^{2})\phi^{\prime\prime}+\frac{\partial V_{\rm
eff}(\phi,T_{+})}{\partial\phi}\right]\phi^{\prime}dz$ $\displaystyle=V_{\rm
eff}(\phi_{-},T_{+})-V_{\rm eff}(\phi_{+},T_{+})\,\,,$ (89a)
$\displaystyle\int\left[(1-v_{w}^{2})\phi^{\prime\prime}+\frac{\partial V_{\rm
eff}(\phi,T_{+})}{\partial\phi}\right](2\phi-\phi_{-})\phi^{\prime}dz$
$\displaystyle=\frac{2(1-v_{w}^{2})\phi_{-}^{3}}{15L^{2}}+\int_{\phi_{+}}^{\phi_{-}}\frac{\partial
V_{\rm eff}(\phi,T_{+})}{\partial\phi}(2\phi-\phi_{-})d\phi\,\,,$ (89b)

The kinetic term $\phi^{\prime\prime}$ tends to stretch the wall (increasing
$L$) while $V_{\rm eff}$ tends to accelerate and compress the wall.

Before further calculations, one should emphasize that the vacuum value
$\phi_{-}$ using here is its value at $z\rightarrow\infty$. More precisely, it
is actually $\phi_{-}=\phi(T_{-})$ that evaluated at $T_{-}$ instead of
$\phi(T_{+})$.333Generally $\phi(T_{+})\neq\phi_{+}$ where $\phi(T_{+})$ is
the VEV at $T_{+}$ but $\phi_{+}$ is the field value in front of the bubble
wall. This temperature jump is actually from the hydrodynamic effects
discussed in Sec.IV. However, as we can see from Eq. (87),
$\phi^{\prime\prime}$, the perturbations $\mu_{i}$, $\delta T_{i}$ and $v_{i}$
should vanish at $z\rightarrow\infty$. On the contrary, $\delta T_{bg}$ is not
zero there. We will find that $\partial V/\partial\phi$ and the friction term
coming from $\delta T_{bg}$ do not exactly cancel with each other. This
behavior appears for the reason that the temperature jump is nonlinear so
$T_{-}\neq T_{+}+\delta T_{bg}$ and the vacuum value $\phi_{0}$ is very
sensitive to temperature.

We can see this more clearly in the following. Since we already had

$\frac{\partial V_{\rm
T}(\phi,T)}{\partial\phi}=\sum_{i}\frac{dm_{i}^{2}}{d\phi}\int\frac{d^{3}p}{(2\pi)^{3}2E_{i}}f_{0,i}(p,T)\,\,,$
(90)

we simply expect that

$0\equiv\frac{\partial V_{\rm
eff}(\phi_{-},T_{-})}{\partial\phi_{-}}\simeq\frac{\partial V_{\rm
eff}(\phi_{-},T_{+})}{\partial\phi_{-}}+\frac{\partial^{2}V_{\rm
T}(\phi_{-},T_{+})}{\partial\phi_{-}\partial T_{+}}\delta T_{bg+}\,\,,$ (91)

where $\delta T_{bg+}=\delta T_{bg}(z\rightarrow\infty)$ is the background
temperature perturbation at large positive $z$. The second term of Eq. (91) is
just the friction term coming from $\delta T_{bg}$. However, in reality, the
temperature varies nonlinearly in front and behind the wall, so this equality
is too crude. In other words, Eq. (87) is not the real EOM at large positive
$z$.

Actually, we can use $\partial V_{\rm T}/\partial T=-s$ such that

$\frac{\partial^{2}V_{\rm T}(\phi_{,}T_{+})}{\partial\phi\partial T_{+}}\delta
T_{bg}(z)=-\frac{\partial s(\phi,T_{+})}{\partial\phi}\delta T_{bg}(z)\,\,,$
(92)

which can be multiplied by $\phi^{\prime}$ and integrated,

$\displaystyle\int dz\phi^{\prime}\left(\frac{\partial V_{\rm
eff}(\phi,T_{+})}{\partial\phi}-\frac{\partial
s(\phi,T_{+})}{\partial\phi}\delta T_{bg}(z)\right)$ $\displaystyle\simeq$
$\displaystyle V_{\rm eff}(\phi_{-},T_{+})-V_{\rm
eff}(\phi_{+},T_{+})-s(\phi_{-},T_{+})\delta T_{bg+}$ (93) $\displaystyle+$
$\displaystyle\langle s\rangle\delta T_{bg+}$ $\displaystyle\simeq$
$\displaystyle V_{\rm eff}(\phi_{-},T_{-})-V_{\rm eff}(\phi_{+},T_{+})+\langle
s\rangle\delta T_{bg+}\,\,,$

where the last term $\langle s\rangle\delta T_{bg+}$ can be approximated as
$\langle s\rangle(T_{-}-T_{+})$ at first order which is consistent with Eq.
(26).

So in order to minimize this inconsistency, we should find a new vacuum value
$\phi_{-}$ to make EOM hold at infinity Friedlander et al. (2021),

$\left.\left(\frac{\partial
V_{\mathrm{eff}}\left(\phi,T_{+}\right)}{\partial\phi}+\left(\frac{N_{t}T_{+}}{2}\frac{dm_{t}^{2}}{d\phi}c_{2}^{t}+\sum_{b}\frac{N_{b}T_{+}}{2}\frac{dm_{b}^{2}}{d\phi}c_{2}^{b}\right)\delta
T_{bg}\right)\right|_{\phi=\phi_{-},z\rightarrow\infty}=0\,\,.$ (94)

The new $\phi_{-}$ does not minimize the potential but it cancels the friction
term. Then we should recompute the perturbations and iterate.

Ref. Friedlander et al. (2021) also introduced another method that one can
find an $\mathcal{O}(1)$ parameter $y$ so that the modified EOM

$\displaystyle\begin{aligned} (1-v_{w}^{2})\phi^{\prime\prime}+\frac{\partial
V_{\rm
eff}(\phi,T_{+})}{\partial\phi}&+\frac{N_{t}T_{+}}{2}\frac{dm_{t}^{2}}{d\phi}\left(c_{1}^{t}\mu_{t}+c_{2}^{t}(\delta
T_{t}+y\delta T_{bg})\right)\\\
&+\sum_{b}\frac{N_{b}T_{+}}{2}\frac{dm_{b}^{2}}{d\phi}\left(c_{1}^{b}\mu_{b}+c_{2}^{b}(\delta
T_{b}+y\delta T_{bg})\right)=0\end{aligned}$ (95)

is satisfied for larger positive values of $z$. In this work we will use the
former method that redefines the VEV.

The last term of Eq. (87), multiplied by $\phi^{\prime}$, gives the friction
term

$F(z)=\frac{N_{t}T_{+}}{2}\frac{dm_{t}^{2}}{dz}\left(c_{1}^{t}\mu_{t}+c_{2}^{t}(\delta
T_{t}+\delta
T_{bg})\right)+\sum_{b}\frac{N_{b}T_{+}}{2}\frac{dm_{b}^{2}}{dz}\left(c_{1}^{b}\mu_{b}+c_{2}^{b}(\delta
T_{b}+\delta T_{bg})\right)\,\,.$ (96)

Then the $M_{1}$ of Eq. (88) gives the condition of steady wall expansion,

$V_{\rm eff}(\phi_{-},T_{+})-V_{\rm eff}(\phi_{+},T_{+})=\int dzF(z)\,\,.$
(97)

We can see the behavior of $\int dzF(z)$ in Fig. 6. We find that when $v_{w}$
approaches $c_{s}$ the friction force will approach negative infinity which is
due to the singularity of $\delta T_{bg}$, as seen in Eq. (LABEL:background).
This singularity may be caused by the selection of ansatz Laurent and Cline
(2020); De Curtis et al. (2022) or some physical reasons Dorsch et al. (2021,
2022); Laurent and Cline (2022) but these works have shown that in the low-
velocity regime our treatment still works well. Fortunately, we found that the
velocity in the IDM is indeed small enough.

Figure 6: Left panel: An illustration for the perturbation of background
velocity and background temperature for $v_{w}=0.1,L=0.1~{}\rm GeV^{-1}$.
Right panel: $\int dzF(z)/T_{N}^{4}$ for different species as a function of
wall velocity $v_{w}$.

To do other integrals, it is much simpler to work in momentum space. However,
it should be noted that $m_{b}=m_{b}(z)$ so the $c_{1}^{b}=c_{1}^{b}(z)$ is
function of $z$. But we have noticed that the source term of perturbations is
proportional to

$\phi\phi^{\prime}=\frac{\phi_{-}^{2}}{4}{\rm
sech}^{2}(\frac{z}{L})\left(1+{\rm tanh}(\frac{z}{L})\right)\,\,,$ (98)

which has a peak at $z=\frac{\rm{ln}2}{2}L$. Therefore, the main contributions
should come from $z$ around the peak. So it is reasonable to make an
approximation that $c_{1}^{b}\simeq c_{1}^{b}(z=\frac{\rm{ln}2}{2}L)$.

Then the Boltzmann equations can be written in Fourier space as

$\displaystyle
ik\tilde{\delta_{i}}+(\hat{A}^{-1}\Gamma)_{ij}\tilde{\delta_{j}}$
$\displaystyle=$
$\displaystyle(\hat{A}^{-1}K)_{i}\widetilde{\phi\phi^{\prime}}\,\,,$ (99)
$\displaystyle ik\widetilde{\delta T}_{bg}$ $\displaystyle=$ $\displaystyle
R_{i}\tilde{\delta_{i}}\,\,,$ (100)

where $K$ and $R$ can be read from Eqs. (65) and (LABEL:background),
respectively.

$\displaystyle
K=\frac{v_{w}}{2T}\left(c_{1}^{t}y_{t}^{2},c_{2}^{t}y_{t}^{2},0,\frac{c_{1}^{W}g_{w}^{2}}{2},\frac{c_{2}^{W}g_{w}^{2}}{2},0,c_{1}^{A}\lambda_{3},c_{2}^{A}\lambda_{3},0\right)\,\,,$
(101) $\displaystyle
R=\frac{1}{\tilde{c}_{4}\left(1/3-v_{w}^{2}\right)}\left(-v_{w}N_{t}\Gamma_{\mu
2,t},-v_{w}N_{t}\Gamma_{T2,t},N_{t}\Gamma_{v2,t},-v_{w}N_{W}\Gamma_{\mu
2,W},-v_{w}N_{W}\Gamma_{T2,W},\right.$
$\displaystyle\left.N_{W}\Gamma_{v2,W},-v_{w}N_{A}\Gamma_{\mu
2,A},-v_{w}N_{A}\Gamma_{T2,A},N_{A}\Gamma_{v2,A}\right)\,\,.$ (102)

Substituting these into Eq. (88), we get the contributions of perturbations to
$M_{1}$ and $M_{2}$,

$\displaystyle\Omega_{1}=\int\frac{dk}{2\pi}\left[f_{i}\chi_{ij}\frac{S_{j}}{\rho_{j}+ik}+\Upsilon
R_{i}\chi_{ij}\left(\frac{S_{j}}{ik(\rho_{j}+ik)}-\frac{S_{j}\pi}{\rho_{j}}\delta(k)\right)\right]\widetilde{\phi\phi^{\prime}}(k)\widetilde{\phi\phi^{\prime}}(-k)\,\,,$
(103)
$\displaystyle\Omega_{2}=\int\frac{dk}{2\pi}\left[f_{i}\chi_{ij}\frac{S_{j}}{\rho_{j}+ik}+\Upsilon
R_{i}\chi_{ij}\left(\frac{S_{j}}{ik(\rho_{j}+ik)}-\frac{S_{j}\pi}{\rho_{j}}\delta(k)\right)\right]\widetilde{\phi\phi^{\prime}}(k)\widetilde{2\phi^{2}\phi^{\prime}}(-k)-\phi_{-}\Omega_{1}\,\,,$
(104)

where $S=\chi^{-1}\hat{A}^{-1}K$, $f$ and $\Upsilon$ can be read from Eq.
(87),

$f=T\left(\frac{c_{1}^{t}N_{t}y_{t}^{2}}{2},\frac{c_{2}^{t}N_{t}y_{t}^{2}}{2},0,\frac{c_{1}^{W}N_{W}g_{w}^{2}}{4},\frac{c_{2}^{W}N_{W}g_{w}^{2}}{4},0,\frac{c_{1}^{A}N_{A}\lambda_{3}}{2},\frac{c_{2}^{A}N_{A}\lambda_{3}}{2},0\right)\,\,,$
(105)
$\Upsilon=\frac{c_{2}^{t}N_{t}y_{t}^{2}}{2}+\frac{c_{2}^{w}N_{W}g_{w}^{2}}{4}+\frac{c_{2}^{A}N_{A}\lambda_{3}}{2}\,\,.$
(106)

It should be noted that we introduce Dirac delta function $\delta(k)$ in Eqs.
(103) and (104) because Eq. (100) can only determine $\delta T_{bg}$ up to a
constant of integration. In order to satisfy the boundary condition $\delta
T_{bg}(z\rightarrow-\infty)=0$, we must add this into the equation.

Using

$\displaystyle\widetilde{\phi\phi^{\prime}}(k)=\frac{\phi_{-}^{2}}{2}(1-ikL/2)\frac{kL\pi}{2}\mathrm{csch}\frac{kL\pi}{2}\,\,,$
(107)
$\displaystyle\widetilde{2\phi^{2}\phi^{\prime}}(k)=\frac{\phi_{-}^{3}}{12}(8-6ikL-k^{2}L^{2})\frac{kL\pi}{2}\mathrm{csch}\frac{kL\pi}{2}\,\,,$
(108)

Eqs. (103) and (104) can be solved by using the following integrals:

$\displaystyle\int\frac{dk}{2\pi}\frac{1}{\rho+ik}\widetilde{\phi\phi^{\prime}}(k)\widetilde{\phi\phi^{\prime}}(-k)$
$\displaystyle=\frac{\phi_{-}^{4}}{16}\left[\left(\rho L-\frac{(\rho
L)^{3}}{4}\right)I_{1}\left(\frac{\rho L\pi}{2}\right)+\frac{\rho
L}{3}\right]\,\,,$ (109a)
$\displaystyle\int\frac{dk}{2\pi}\frac{1}{ik(\rho+ik)}\widetilde{\phi\phi^{\prime}}(k)\widetilde{\phi\phi^{\prime}}(-k)$
$\displaystyle=-\frac{1}{\rho}\frac{\phi_{-}^{4}}{16}\left[\left(\rho
L-\frac{(\rho L)^{3}}{4}\right)I_{1}\left(\frac{\rho
L\pi}{2}\right)+\frac{\rho L}{3}\right]\,\,,$ (109b)
$\displaystyle\int\frac{dk}{2\pi}\frac{1}{\rho+ik}\widetilde{\phi\phi^{\prime}}(k)\widetilde{2\phi^{2}\phi^{\prime}}(-k)$
$\displaystyle=\frac{\phi_{-}^{5}}{12}\left[\left(\frac{\rho^{4}L^{4}}{16}-\frac{\rho^{3}L^{3}}{4}-\frac{\rho^{2}L^{2}}{4}+\rho
L\right)I_{1}(\frac{\rho L\pi}{2})\right.$
$\displaystyle\left.-\left(\frac{\rho^{2}L^{2}}{12}-\frac{\rho
L}{3}-\frac{2}{5}\right)\right]\,\,,$ (109c)
$\displaystyle\int\frac{dk}{2\pi}\frac{1}{ik(\rho+ik)}\widetilde{\phi\phi^{\prime}}(k)\widetilde{2\phi^{2}\phi^{\prime}}(-k)$
$\displaystyle=-\frac{1}{\rho}\frac{\phi_{-}^{5}}{12}\left[\left(\frac{\rho^{4}L^{4}}{16}-\frac{\rho^{3}L^{3}}{4}-\frac{\rho^{2}L^{2}}{4}+\rho
L\right)I_{1}(\frac{\rho L\pi}{2})\right.$
$\displaystyle\left.-\left(\frac{\rho^{2}L^{2}}{12}-\frac{\rho
L}{3}\right)\right]\,\,,$ (109d)

where

$I_{1}(a)=\int_{-\infty}^{\infty}dx\frac{x^{2}\mathrm{csch}^{2}x}{x^{2}+a^{2}}\,\,,$
(110)

which can be evaluated by contour integrals in the complex plane,

$I_{1}(a)=\frac{\pi|a|}{\mathrm{sin}^{2}a}-2-\sum_{n=1}^{\infty}\frac{n(2\pi
a)^{2}}{[(n\pi)^{2}-a^{2}]^{2}}\,\,.$ (111)

Figure 7: Illustration of the two constraints $M_{1}/T_{N}^{4}$ and
$M_{2}/T_{N}^{4}$. The red dots represent the final bubble wall velocity and
bubble wall thickness for the benchmark parameters in IDM.

In summary, the main procedure of calculating the bubble wall velocity is
shown in Fig. 1. For a given set of model parameters (in this case, the
Benchmark A of IDM), we firstly calculate the thermodynamic quantities like
the nucleation temperature $T_{N}$ and the phase transition strength
$\alpha_{N}$ of the nucleation process with the ComoTransitions package. Then
we perform a grid scan in the $(v_{w},L)$ space. At each point in the grid we
perform the next steps:

1\. Solve the hydrodynamic equations of the SFOPT. This allows us to obtain
the thermodynamic parameters evaluated in front of and behind the wall, that
is $T_{+},\alpha_{+},v_{+},\phi_{+}$, $T_{-}$, and $v_{-}$. Then we modify the
field value by using the correct minimization conditions with the effective
potential which should be evaluated at the temperature inside the bubble, with

$\left.\frac{\partial
V_{\mathrm{eff}}\left(\phi,T_{-}\right)}{\partial\phi}\right|_{\phi=\phi_{-}}=0\,\,.$
(112)

3\. Solve the truncated Boltzmann equations for concerned perturbations, which
determines the friction force from the massive particles.

4\. Choose new vacuum value $\phi_{-}$ such that it can satisfy the condition

$\left.\left(\frac{\partial
V_{\mathrm{eff}}\left(\phi,T_{+}\right)}{\partial\phi}+\sum_{i}N_{i}\frac{dm_{i}^{2}}{d\phi}\int\frac{d^{3}p}{(2\pi)^{3}2E_{i}}\delta
f_{i}(p,z)\right)\right|_{\phi=\phi_{-},z\rightarrow\infty}=0\,\,.$ (113)

5\. Recompute the perturbations and iterate. Then calculate $M_{1}$ and
$M_{2}$.

In Fig. 7, we show the two constraints $M_{1}$ and $M_{2}$ as functions of
$v_{w}$ and $L$. For the two red dots, $M_{1}=0$ and $M_{2}=0$ are satisfied
simultaneously. The two dots correspond to the same value of bubble wall
velocity and thickness. Therefore, the solution of two constraints is
$v_{w}\simeq 0.165,L\simeq 0.084~{}\rm GeV^{-1}$ which implies that the wall
velocity is subsonic in IDM . This solution is consistent with our assumption
that $LT_{N}\gg 1$ such the WKB approximation holds in this case.

In Table 3, we also give the results of Benchmark B and C. We can see that in
the allowed parameter spaces of IDM, the bubble wall velocity varies slightly
around $v_{w}\simeq 0.165$.

There are some uncertainties in our calculation. First, we parametrize the
nonequilibrium with leading-order perturbations which may be problematic.
Recently there were some modified schemes Laurent and Cline (2020); De Curtis
et al. (2022); Dorsch et al. (2021, 2022); Laurent and Cline (2022) but they
have shown that in low-velocity regime the result of new scheme is almost the
same as in the old scheme. Second, we calculated collision terms by using
leading-log approximation. Actually, if the collision term is larger, we
expect that the velocity would be increased because the massive particles
would be more closely to equilibrium. Then the friction would be decreased.

|  $T_{c}$ [GeV] |  $T_{N}$ [GeV] | $v_{w}$ |  $L$ [$\rm GeV^{-1}$]
---|---|---|---|---
Benchmark A | 118.3 | 117.1 | 0.165 | 0.084
Benchmark B | 118.6 | 117.5 | 0.164 | 0.085
Benchmark C | 119.4 | 118.4 | 0.164 | 0.088
Table 3: Bubble wall velocity, bubble wall thickness and phase transition
parameters for different benchmark points.

##  VIII Conclusion

It is crucial and complicated to precisely calculate the bubble wall velocity,
which is essential to the phase transition dynamics, the electroweak
baryogenesis, the new dark matter mechanism from the bubbles, and the spectra
of the phase transition gravitational wave. And more gravitational wave
experiments (LISA, TianQin, Taiji,…) are proposed and they need accurate
theoretic predictions on the gravitational wave signals. Thus, it becomes
necessary and useful to calculate the model-dependent bubble wall velocity
both for the fundamental problems and the gravitational wave experiments
Azatov and Vanvlasselaer (2021); Azatov et al. (2021).

In this work, we have systematically calculated the bubble wall velocity in
the well-motivated inert doublet model for the first time. The contribution
from the heavy inert scalar bosons are taken into account. All the thermal
masses are included in our numerical calculations. One difficulty is to
consistently consider the hydrodynamic effects including the heating effects.
To obtain the correct bubble wall velocity, it is important to figure out the
correct temperature and vacuum value. We consider the hydrodynamic effects to
obtain more reliable results of the temperature and vacuum value. Another
difficulty is the precise calculations of the collision terms, which represent
various particle scattering processes at the vicinity of the bubble wall so we
have used the Monte Carlo algorithm to numerically obtain more precise
collision terms. After having more rigorous hydrodynamic effects and collision
terms, we get relatively more reliable nonequilibrium perturbations of various
massive particles, which are essential for the friction force. Finally, by
scanning the grid of wall velocity and wall thickness, we obtain the
relatively precise bubble wall velocity in the inert doublet model for the
first time. It is obviously smaller than the speed of sound and is favored by
the traditional electroweak baryogenesis. This result could also help us to
greatly reduce the uncertainty of the gravitational wave spectra from the
strong first-order phase transition.

Our results might be useful for other similar models. This procedure is also
appropriate for other standard model-like models that have weak first-order
phase transitions.444For strong enough phase transitions, it is necessary to
introduce new methods to solve Boltzmann equations Laurent and Cline (2020);
De Curtis et al. (2022); Laurent and Cline (2022).  Regarding the precise
prediction of bubble wall velocity for a given new physics model, there are
still some uncertainties from the collision terms and fluid approximations.
More precise calculations on the collision terms Arnold et al. (2000, 2003a,
2003b); Bodeker and Moore (2017); Wang et al. (2020b) including the
resummation over the large logarithmic terms Höche et al. (2021); Gouttenoire
et al. (2022) and more accurate fluid ansatzes are left for the future study.

###### Acknowledgements.

We acknowledge valuable discussions with Eibun Senaha on the renormalization
and daisy resummation schemes. The authors thank David Tucker-Smith, Marek
Lewicki, Thomas Konstandin, Glauber Carvalho Dorsch, Benoit Laurent, Denis
Werth and Carlos Tamarit for helpful correspondence. This work is supported by
the National Natural Science Foundation of China (NNSFC) under Grant No.
12205387. X.W. is supported by the China Postdoctoral Science Foundation under
Grant No. 2022M713642. This work is supported in part by by Guangdong Major
Project of Basic and Applied Basic Research (Grant No. 2019B030302001).

## References

  * Cline and Kainulainen (2020) J. M. Cline and K. Kainulainen, Phys. Rev. D 101, 063525 (2020), eprint 2001.00568.
  * Baker et al. (2020) M. J. Baker, J. Kopp, and A. J. Long, Phys. Rev. Lett. 125, 151102 (2020), eprint 1912.02830.
  * Chway et al. (2020) D. Chway, T. H. Jung, and C. S. Shin, Phys. Rev. D 101, 095019 (2020), eprint 1912.04238.
  * Huang and Li (2017) F. P. Huang and C. S. Li, Phys. Rev. D 96, 095028 (2017), eprint 1709.09691.
  * Espinosa et al. (2010) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 06, 028 (2010), eprint 1004.4187.
  * Giese et al. (2020) F. Giese, T. Konstandin, and J. van de Vis, JCAP 07, 057 (2020), eprint 2004.06995.
  * Giese et al. (2021) F. Giese, T. Konstandin, K. Schmitz, and J. van de Vis, JCAP 01, 072 (2021), eprint 2010.09744.
  * Wang et al. (2021a) X. Wang, F. P. Huang, and X. Zhang, Phys. Rev. D 103, 103520 (2021a), eprint 2010.13770.
  * Leitao and Megevand (2011) L. Leitao and A. Megevand, Nucl. Phys. B 844, 450 (2011), eprint 1010.2134.
  * Leitao and Megevand (2015) L. Leitao and A. Megevand, Nucl. Phys. B 891, 159 (2015), eprint 1410.3875.
  * Laurent and Cline (2020) B. Laurent and J. M. Cline, Phys. Rev. D 102, 063516 (2020), eprint 2007.10935.
  * De Curtis et al. (2022) S. De Curtis, L. D. Rose, A. Guiggiani, A. G. Muyor, and G. Panico, JHEP 03, 163 (2022), eprint 2201.08220.
  * Laurent and Cline (2022) B. Laurent and J. M. Cline, Phys. Rev. D 106, 023501 (2022), eprint 2204.13120.
  * Amaro-Seoane et al. (2017) P. Amaro-Seoane et al. (LISA) (2017), eprint 1702.00786.
  * Luo et al. (2016) J. Luo et al. (TianQin), Class. Quant. Grav. 33, 035010 (2016), eprint 1512.02076.
  * Mei et al. (2021) J. Mei et al. (TianQin), PTEP 2021, 05A107 (2021), eprint 2008.10332.
  * Hu and Wu (2017) W.-R. Hu and Y.-L. Wu, Natl. Sci. Rev. 4, 685 (2017).
  * Dine et al. (1992) M. Dine, R. G. Leigh, P. Y. Huet, A. D. Linde, and D. A. Linde, Phys. Rev. D 46, 550 (1992), eprint hep-ph/9203203.
  * Liu et al. (1992) B.-H. Liu, L. D. McLerran, and N. Turok, Phys. Rev. D 46, 2668 (1992).
  * Ignatius et al. (1994) J. Ignatius, K. Kajantie, H. Kurki-Suonio, and M. Laine, Phys. Rev. D 49, 3854 (1994), eprint astro-ph/9309059.
  * Moore and Prokopec (1995a) G. D. Moore and T. Prokopec, Phys. Rev. D 52, 7182 (1995a), eprint hep-ph/9506475.
  * Moore and Prokopec (1995b) G. D. Moore and T. Prokopec, Phys. Rev. Lett. 75, 777 (1995b), eprint hep-ph/9503296.
  * Moore (2000a) G. D. Moore, JHEP 03, 006 (2000a), eprint hep-ph/0001274.
  * John and Schmidt (2001) P. John and M. G. Schmidt, Nucl. Phys. B 598, 291 (2001), [Erratum: Nucl.Phys.B 648, 449–452 (2003)], eprint hep-ph/0002050.
  * Huber and Sopena (2013) S. J. Huber and M. Sopena (2013), eprint 1302.1044.
  * Mégevand (2013) A. Mégevand, JCAP 07, 045 (2013), eprint 1303.4233.
  * Megevand and Sanchez (2009) A. Megevand and A. D. Sanchez, Nucl. Phys. B 820, 47 (2009), eprint 0904.1753.
  * Lewicki et al. (2022) M. Lewicki, M. Merchand, and M. Zych, JHEP 02, 017 (2022), eprint 2111.02393.
  * Dorsch et al. (2021) G. C. Dorsch, S. J. Huber, and T. Konstandin, JCAP 08, 020 (2021), eprint 2106.06547.
  * Dorsch et al. (2022) G. C. Dorsch, S. J. Huber, and T. Konstandin, JCAP 04, 010 (2022), eprint 2112.12548.
  * Barroso Mancha et al. (2021) M. Barroso Mancha, T. Prokopec, and B. Swiezewska, JHEP 01, 070 (2021), eprint 2005.10875.
  * Balaji et al. (2021) S. Balaji, M. Spannowsky, and C. Tamarit, JCAP 03, 051 (2021), eprint 2010.08013.
  * Ai et al. (2022) W.-Y. Ai, B. Garbrecht, and C. Tamarit, JCAP 03, 015 (2022), eprint 2109.13710.
  * Wang and Yuwen (2022) S.-J. Wang and Z.-Y. Yuwen (2022), eprint 2205.02492.
  * Barbieri et al. (2006) R. Barbieri, L. J. Hall, and V. S. Rychkov, Phys. Rev. D 74, 015007 (2006), eprint hep-ph/0603188.
  * Lopez Honorez and Yaguna (2010) L. Lopez Honorez and C. E. Yaguna, JHEP 09, 046 (2010), eprint 1003.3125.
  * Chowdhury et al. (2012) T. A. Chowdhury, M. Nemevsek, G. Senjanovic, and Y. Zhang, JCAP 02, 029 (2012), eprint 1110.5334.
  * Barducci et al. (2018) D. Barducci, G. Belanger, J. Bernon, F. Boudjema, J. Da Silva, S. Kraml, U. Laa, and A. Pukhov, Comput. Phys. Commun. 222, 327 (2018), eprint 1606.03834.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Wang et al. (2021b) Y. Wang, C. S. Li, and F. P. Huang, Phys. Rev. D 104, 053004 (2021b), eprint 2012.03920.
  * Huang and Yu (2018) F. P. Huang and J.-H. Yu, Phys. Rev. D 98, 095022 (2018), eprint 1704.04201.
  * Borah and Cline (2012) D. Borah and J. M. Cline, Phys. Rev. D 86, 055001 (2012), eprint 1204.4722.
  * Gil et al. (2012) G. Gil, P. Chankowski, and M. Krawczyk, Phys. Lett. B 717, 396 (2012), eprint 1207.0084.
  * Cline and Kainulainen (2013) J. M. Cline and K. Kainulainen, Phys. Rev. D 87, 071701 (2013), eprint 1302.2614.
  * AbdusSalam and Chowdhury (2014) S. S. AbdusSalam and T. A. Chowdhury, JCAP 05, 026 (2014), eprint 1310.8152.
  * Blinov et al. (2015) N. Blinov, S. Profumo, and T. Stefaniak, JCAP 07, 028 (2015), eprint 1504.05949.
  * Cao et al. (2018) Q.-H. Cao, F. P. Huang, K.-P. Xie, and X. Zhang, Chin. Phys. C 42, 023103 (2018), eprint 1708.04737.
  * Laine et al. (2017) M. Laine, M. Meyer, and G. Nardini, Nucl. Phys. B 920, 565 (2017), eprint 1702.07479.
  * Senaha (2019) E. Senaha, Phys. Rev. D 100, 055034 (2019), eprint 1811.00336.
  * Huang and Senaha (2019) F. P. Huang and E. Senaha, Phys. Rev. D 100, 035014 (2019), eprint 1905.10283.
  * Kainulainen et al. (2019) K. Kainulainen, V. Keus, L. Niemi, K. Rummukainen, T. V. I. Tenkanen, and V. Vaskonen, JHEP 06, 075 (2019), eprint 1904.01329.
  * Niemi et al. (2021) L. Niemi, P. Schicho, and T. V. I. Tenkanen, Phys. Rev. D 103, 115035 (2021), eprint 2103.07467.
  * Schicho et al. (2022) P. Schicho, T. V. I. Tenkanen, and G. White, JHEP 11, 047 (2022), eprint 2203.04284.
  * Dolan and Jackiw (1974) L. Dolan and R. Jackiw, Phys. Rev. D 9, 3320 (1974).
  * Belyaev et al. (2018) A. Belyaev, G. Cacciapaglia, I. P. Ivanov, F. Rojas-Abatte, and M. Thomas, Phys. Rev. D 97, 035011 (2018), eprint 1612.00511.
  * Wainwright (2012) C. L. Wainwright, Comput. Phys. Commun. 183, 2006 (2012), eprint 1109.4189.
  * Konstandin et al. (2014) T. Konstandin, G. Nardini, and I. Rues, JCAP 09, 028 (2014), eprint 1407.3132.
  * Konstandin and No (2011) T. Konstandin and J. M. No, JCAP 02, 008 (2011), eprint 1011.3735.
  * Steinhardt (1982) P. J. Steinhardt, Phys. Rev. D 25, 2074 (1982).
  * Moore (2000b) G. D. Moore, Phys. Rev. D 62, 085011 (2000b), eprint hep-ph/0001216.
  * Kozaczuk (2015) J. Kozaczuk, JHEP 10, 135 (2015), eprint 1506.04741.
  * Hahn (2005) T. Hahn, Comput. Phys. Commun. 168, 78 (2005), eprint hep-ph/0404043.
  * Wang et al. (2020a) X. Wang, F. P. Huang, and X. Zhang, JCAP 05, 045 (2020a), eprint 2003.08892.
  * Friedlander et al. (2021) A. Friedlander, I. Banta, J. M. Cline, and D. Tucker-Smith, Phys. Rev. D 103, 055020 (2021), eprint 2009.14295.
  * Azatov and Vanvlasselaer (2021) A. Azatov and M. Vanvlasselaer, JCAP 01, 058 (2021), eprint 2010.02590.
  * Azatov et al. (2021) A. Azatov, M. Vanvlasselaer, and W. Yin, JHEP 03, 288 (2021), eprint 2101.05721.
  * Arnold et al. (2000) P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 11, 001 (2000), eprint hep-ph/0010177.
  * Arnold et al. (2003a) P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 01, 030 (2003a), eprint hep-ph/0209353.
  * Arnold et al. (2003b) P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 05, 051 (2003b), eprint hep-ph/0302165.
  * Bodeker and Moore (2017) D. Bodeker and G. D. Moore, JCAP 05, 025 (2017), eprint 1703.08215.
  * Wang et al. (2020b) X. Wang, F. P. Huang, and X. Zhang (2020b), eprint 2011.12903.
  * Höche et al. (2021) S. Höche, J. Kozaczuk, A. J. Long, J. Turner, and Y. Wang, JCAP 03, 009 (2021), eprint 2007.10343.
  * Gouttenoire et al. (2022) Y. Gouttenoire, R. Jinno, and F. Sala, JHEP 05, 004 (2022), eprint 2112.07686.