# Phase transition gravitational waves from pseudo-Nambu-Goldstone dark matter
and two Higgs doublets

Zhao Zhang  Chengfeng Cai  Xue-Min Jiang  Yi-Lei Tang
<EMAIL_ADDRESS>Zhao-Huan Yu<EMAIL_ADDRESS>Hong-Hao
Zhang<EMAIL_ADDRESS>School of Physics, Sun Yat-Sen University,
Guangzhou 510275, China

###### Abstract

We investigate the potential stochastic gravitational waves from first-order
electroweak phase transitions in a model with pseudo-Nambu-Goldstone dark
matter and two Higgs doublets. The dark matter candidate can naturally evade
direct detection bounds, and can achieve the observed relic abundance via the
thermal mechanism. Three scalar fields in the model obtain vacuum expectation
values, related to phase transitions at the early Universe. We search for the
parameter points that can cause first-order phase transitions, taking into
account the existed experimental constraints. The resulting gravitational wave
spectra are further evaluated. Some parameter points are found to induce
strong gravitational wave signals, which have the opportunity to be detected
in future space-based interferometer experiments LISA, Taiji, and TianQin.

###### Contents

  1. I Introduction
  2. II The model
  3. III Experimental bounds
  4. IV Effective potential
  5. V Phase transitions
  6. VI Gravitational wave spectra
  7. VII Numerical analyses
  8. VIII Summary

##  I Introduction

It is conventionally believed that dark matter (DM) originates from thermal
production at the early Universe Bertone:2004pz ; Feng:2010gw ; Young:2016ala
. Thus, the DM relic abundance would be determined by the annihilation cross
section at the freeze-out epoch. The relic abundance observation suggests that
the natural strength of the DM couplings to standard model (SM) particles
should be close to the weak interaction strength. This motivates the worldwide
establishment of various direct detection experiments searching for nuclear
recoil signals induced by DM scattering. Nonetheless, no DM signal is robustly
found in these experiments so far, leading to stringent constraints on the DM-
nucleon scattering cross section Akerib:2016vxi ; Cui:2017nnn ; Aprile:2018dbl
. Therefore, the thermal production paradigm faces a serious challenge.

Such a situation can be circumvented if one can effectively suppress DM-
nucleon scattering at zero momentum transfer without reducing DM annihilation
at the freeze-out epoch. An appealing approach to achieve this is provided by
Higgs-portal pseudo-Nambu-Goldstone boson (pNGB) DM models Gross:2017dan ;
Azevedo:2018exj ; Ishiwata:2018sdi ; Huitu:2018gbc ; Alanne:2018zjm ;
Kannike:2019wsn ; Karamitros:2019ewv ; Cline:2019okt ; Jiang:2019soj ;
Arina:2019tib ; Abe:2020iph ; Okada:2020zxo ; Glaus:2020ihj ; Abe:2020ldj ,
where the DM candidate is a pNGB protected by a global symmetry which is
softly broken by quadratic mass terms. The pNGB nature makes the tree-level
DM-nucleon scattering amplitude vanish in the zero momentum transfer limit
Gross:2017dan . Although loop corrections give rise to a nonzero scattering
cross section, one-loop calculations show that near future direct detection
experiments would not be able to probe the pNGB DM Azevedo:2018exj ;
Ishiwata:2018sdi ; Glaus:2020ihj . Therefore, other experimental approaches
are crucial for exploring these models.

There are various ways to experimentally test DM models, including direct and
indirect DM detection, collider searches, etc. The discovery of gravitational
waves (GWs) by LIGO Abbott:2016blz  provides a new path. DM fields could be
relevant to strong first-order electroweak phase transitions (EWPTs) that
produce detectable stochastic GW signals in the proposed future GW experiments
Jiang:2015cwa ; Chala:2016ykx ; Chao:2017vrq ; Beniwal:2017eik ; Huang:2017rzf
; Huang:2017kzu ; Hektor:2018esx ; Baldes:2018emh ; Madge:2018gfl ;
Beniwal:2018hyi ; Bian:2018mkl ; Bian:2018bxr ; Shajiee:2018jdq ;
Mohamadnejad:2019vzg ; Kannike:2019mzk ; Paul:2019pgt ; Chen:2019ebq ;
Barman:2019oda ; Chiang:2019oms ; Borah:2020wut ; Kang:2020jeg ;
Pandey:2020hoq ; Han:2020ekm ; Alanne:2020jwx ; Wang:2020wrk ; Ghosh:2020ipy ;
Huang:2020mso ; Chao:2020adk . Such stochastic GWs typically peak around the
mHz frequency band Mazumdar:2018dfl , to which ground-based laser
interferometers are not sensitive. Nonetheless, future space-based GW
interferometer plans, e.g., LISA Audley:2017drz , TianQin Luo:2015ght ;
Hu:2018yqb ; Mei:2020lrl , Taiji Hu:2017mde ; Guo:2018npi , DECIGO Seto:2001qf
; Kudoh:2005as , and BBO Ungarelli:2005qb ; Cutler:2005qq , are able to probe
sub-Hz bands and look for the stochastic GW signals.

The minimal setup of the pNGB DM involves a complex scalar singlet with a
global $\mathrm{U}(1)$ symmetry and a quadratic term that softly breaks
$\mathrm{U}(1)$ into $Z_{2}$ Gross:2017dan . The singlet and the SM Higgs
doublet together could induce two-step phase transitions. Nevertheless, a
study Kannike:2019wsn  showed that such phase transitions can only be of
second order and impossible to produce stochastic GWs. Further studies tried
to introduce extra terms to break the $\mathrm{U}(1)$ symmetry, e.g., the soft
cubic terms Kannike:2019mzk , or the most general breaking terms
Alanne:2020jwx . These efforts successfully achieved first-order phase
transitions (FOPTs) and stochastic GWs, but the essential merit of the
vanishing tree-level DM-nucleon scattering in the zero momentum transfer limit
is sacrificed.

It would be rather interesting if we can find out a pNGB DM setup that allows
both the vanishing DM-nucleon scattering and detectable stochastic GW signals.
Inspired by the notable GW signals from the strong FOPTs obtained in the two-
Higgs-doublet models Dorsch:2016nrg ; Wang:2019pet ; Zhou:2020xqi , we study
the possibility of extending the minimal pNGB DM setup with an additional
Higgs doublet Jiang:2019soj , which is expected to involve more phase
transition patterns. The corresponding phase transitions and stochastic GWs
will be studied in this paper in detail.

In the following Sec. II, we briefly introduce the model and the particle
masses. Existed experimental bounds are described in Sec. III. The effective
potential at finite temperature is constructed in Sec. IV. We analyze key
properties of the EWPT in Sec. V, which are relevant to the GW spectra
discussed in Sec. VI. Sec. VII give numerical analyses of GW signals based on
random parameter scans. We summarize the paper in Sec. VIII.

##  II The model

In this section, we briefly describe the model we are interested in. More
details can be found in Ref. Jiang:2019soj . This model involves two
$\mathrm{SU}(2)_{\mathrm{L}}$ Higgs doublets $\Phi_{1}$ and $\Phi_{2}$, both
carrying hypercharge $1/2$, and a complex scalar $S$ which is a SM gauge
singlet. The Lagrangian respects a global $\mathrm{U}(1)$ symmetry $S\to
e^{i\alpha}S$ explicitly violated into a $Z_{2}$ symmetry $S\to-S$ by soft
breaking quadratic terms. The $\mathrm{U}(1)$ symmetry is further
spontaneously broken after $S$ develops a vacuum expectation value (VEV). Then
the imaginary part of $S$ becomes a pNGB, acting as a DM candidate.

As in the simplified versions of the two-Higgs-doublet models Branco:2011iw ,
we assume that the scalar potential respects a $Z_{2}$ symmetry
$\Phi_{1}\to-\Phi_{1}$ or $\Phi_{2}\to-\Phi_{2}$ which is only softly broken
by quadratic terms. Moreover, $CP$ conservation is assumed in the scalar
sector, leading to only real coefficients. The potential satisfying these two
assumptions and the global $\mathrm{U}(1)$ symmetry reads

$\displaystyle V_{\Phi_{i},S}$
$\displaystyle=m_{11}^{2}|{\Phi_{1}}|^{2}+m_{22}^{2}|{\Phi_{2}}|^{2}-m_{12}^{2}(\Phi_{1}^{\dagger}{\Phi_{2}}+\Phi_{2}^{\dagger}{\Phi_{1}})+\frac{\lambda_{1}}{2}|{\Phi_{1}}{|^{4}}+\frac{\lambda_{2}}{2}|{\Phi_{2}}{|^{4}}$
$\displaystyle\quad+{\lambda_{3}}|{\Phi_{1}}|^{2}|{\Phi_{2}}|^{2}+{\lambda_{4}}|\Phi_{1}^{\dagger}{\Phi_{2}}|^{2}+\frac{{{\lambda_{5}}}}{2}[{(\Phi_{1}^{\dagger}{\Phi_{2}})^{2}}+{(\Phi_{2}^{\dagger}{\Phi_{1}})^{2}}]$
$\displaystyle\quad-
m_{S}^{2}|S|^{2}+\frac{{{\lambda_{S}}}}{2}|S{|^{4}}+{\kappa_{1}}|{\Phi_{1}}|^{2}|S|^{2}+{\kappa_{2}}|{\Phi_{2}}|^{2}|S|^{2}.$
(1)

The $\mathrm{U}(1)$ soft breaking terms

$V_{\mathrm{soft}}=-\frac{m^{\prime 2}_{S}}{4}S^{2}+\mathrm{H.c.}$ (2)

are further introduced in the potential. Thus, the total scalar potential is
$V=V_{\Phi_{i},S}+V_{\mathrm{soft}}$.

We can always make the soft breaking parameter $m^{\prime 2}_{S}$ real and
positive through a phase redefinition of $S$. Consequently, the potential
respects a dark $CP$ symmetry $S\to S^{*}$. For $m^{\prime 2}_{S}>0$, the VEV
of $S$ developed must be real, and the dark $CP$ symmetry remain unbroken,
ensuring that the imaginary part of $S$ acts as a stable DM candidate
Gross:2017dan .

If the charged component of $\Phi_{1}$ or $\Phi_{2}$ gains a nonzero VEV, the
photon would become massive, and the theory is unacceptable. If the neutral
component of $\Phi_{1}$ or $\Phi_{2}$ develops an imaginary VEV, $CP$ would be
spontaneously broken. Detailed discussions on vacuum configurations and
parameter relations in general two-Higgs-doublet models can be found in Ref.
Ginzburg:2004vp . Here we are particularly interested in the case that only
the neutral real parts of $\Phi_{1}$, $\Phi_{2}$, and $S$ develop nonzero VEVs
$v_{1}$, $v_{2}$, and $v_{s}$, respectively. Thus, at the zero temperature,
these scalar fields can be expanded as

$\displaystyle\Phi_{1}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}{\phi_{1}^{+}}\\\
{(v_{1}+{\rho_{1}}+i{\eta_{1}})/\sqrt{2}}\\\ \end{pmatrix},$ (3)
$\displaystyle\Phi_{2}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}{\phi_{2}^{+}}\\\
{(v_{2}+{\rho_{2}}+i{\eta_{2}})/\sqrt{2}}\\\ \end{pmatrix},$ (4)
$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{v_{s}+s+i\chi}{\sqrt{2}}.$ (5)

The potential is minimized at $(v_{1},v_{2},v_{s})$, leading to three
stationary point conditions,

$\displaystyle m_{11}^{2}$ $\displaystyle=$
$\displaystyle{m}_{12}^{2}\tan\beta-\frac{1}{2}({\lambda_{1}}v_{1}^{2}+\lambda_{345}v_{2}^{2}+{\kappa_{1}}v_{s}^{2}),$
(6) $\displaystyle m_{22}^{2}$ $\displaystyle=$
$\displaystyle{m}_{12}^{2}\cot\beta-\frac{1}{2}({\lambda_{2}}v_{2}^{2}+\lambda_{345}v_{1}^{2}+{\kappa_{2}}v_{s}^{2}),$
(7) $\displaystyle m_{S}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}({\kappa_{1}}v_{1}^{2}+{\kappa_{2}}v_{2}^{2}+{\lambda_{S}}v_{s}^{2}-m^{\prime
2}_{S}),$ (8)

where

$\displaystyle\beta$ $\displaystyle\equiv$
$\displaystyle\arctan\frac{v_{2}}{v_{1}},$ (9) $\displaystyle\lambda_{345}$
$\displaystyle\equiv$
$\displaystyle{\lambda_{3}}+{\lambda_{4}}+{\lambda_{5}}.$ (10)

The VEVs contribute to a $3\times 3$ mass-squared matrix for the $CP$-even
neutral scalars $\rho_{1}$, $\rho_{2}$, and $s$. The eigenvalues
$m_{h_{1}}^{2}$, $m_{h_{2}}^{2}$, and $m_{h_{3}}^{2}$ of this matrix are the
masses squared for the mass eigenstates $h_{1}$, $h_{2}$, and $h_{3}$,
respectively. One of $h_{i}$ must behave as a SM-like Higgs boson with a mass
of $\sim 125~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, satisfying the experimental
observations. After rotations with the angle $\beta$, the $CP$-odd neutral
scalars $\eta_{1}$ and $\eta_{2}$ are transformed into the mass eigenstates
$G^{0}$ and $a$, while the charged scalars $\phi_{1}^{+}$ and $\phi_{2}^{+}$
are transformed into the mass eigenstates $G^{+}$ and $H^{+}$. $G^{0}$ and
$G^{\pm}$ are the Nambu-Goldstone bosons eaten by the $Z$ and $W^{\pm}$ gauge
bosons. $a$ and $H^{\pm}$ are extra Higgs bosons, whose masses squared are
given by

$\displaystyle m_{a}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sin\beta\cos\beta}\left(m_{12}^{2}-\lambda_{5}v_{1}v_{2}\right),$
(11) $\displaystyle m_{H^{\pm}}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sin\beta\cos\beta}\left[m_{12}^{2}-\frac{1}{2}(\lambda_{4}+\lambda_{5})v_{1}v_{2}\right].$
(12)

For $m^{\prime 2}_{S}=0$, the neutral boson $\chi$ is a massless Nambu-
Goldstone boson due to the global $\mathrm{U}(1)$ symmetry. The soft breaking
terms endow the pNGB $\chi$ with a mass of

$m_{\chi}=m^{\prime}_{S}.$ (13)

Besides, $\chi$ only appears in pairs in the interaction terms, guaranteeing
its stability to become a DM candidate. The pNGB feature also eliminates the
tree-level $\chi$-nucleon scattering amplitude in the zero momentum transfer
limit without any parameter tuning Gross:2017dan ; Jiang:2019soj . Thus, this
model is hardly constrained by DM direct detection experiments.

The masses of the $W$ and $Z$ gauge bosons are given by

$m_{W}=\frac{gv}{2},\quad m_{Z}=\frac{v}{2}\sqrt{g^{2}+g^{\prime 2}},$ (14)

where $v\equiv\sqrt{v_{1}^{2}+v_{2}^{2}}$, and $g$ and $g^{\prime}$ denote the
$\mathrm{SU}(2)_{\mathrm{L}}$ and $\mathrm{U}(1)_{\mathrm{Y}}$ gauge
couplings, respectively. Thus, we observe that $v$ is equivalent to the Higgs
VEV in the SM and can be expressed as $v=(\sqrt{2}G_{\mathrm{F}})^{-1/2}$,
where $G_{\mathrm{F}}$ is the Fermi constant.

For the two Higgs doublets, four types of Yukawa couplings without tree-level
flavor-changing neutral currents (FCNCs) can be constructed Glashow:1976nt ;
Paschos:1976ay ; Branco:2011iw . In this paper, we only focus on the type-I
and type-II Yukawa couplings, whose Lagrangians are respectively given by

$\displaystyle{\mathcal{L}_{{\mathrm{Y,I}}}}$ $\displaystyle=$
$\displaystyle-\tilde{y}_{d}^{ij}{{\bar{Q}}_{i{\mathrm{L}}}}{d^{\prime}_{j{\mathrm{R}}}}{\Phi_{2}}-\tilde{y}_{u}^{ij}{{\bar{Q}}_{i{\mathrm{L}}}}{u^{\prime}_{j{\mathrm{R}}}}{{\tilde{\Phi}}_{2}}-{\tilde{y}_{{\ell_{i}}}}{{\bar{L}}_{i{\mathrm{L}}}}{\ell_{i{\mathrm{R}}}}{\Phi_{2}}+\mathrm{H.c.},$
(15) $\displaystyle{\mathcal{L}_{{\mathrm{Y,II}}}}$ $\displaystyle=$
$\displaystyle-\tilde{y}_{d}^{ij}{{\bar{Q}}_{i{\mathrm{L}}}}{d^{\prime}_{j{\mathrm{R}}}}{\Phi_{1}}-\tilde{y}_{u}^{ij}{{\bar{Q}}_{i{\mathrm{L}}}}{u^{\prime}_{j{\mathrm{R}}}}{{\tilde{\Phi}}_{2}}-{\tilde{y}_{{\ell_{i}}}}{{\bar{L}}_{i{\mathrm{L}}}}{\ell_{i{\mathrm{R}}}}{\Phi_{1}}+\mathrm{H.c.},$
(16)

where
$L_{i\mathrm{L}}\equiv(\nu_{i\mathrm{L}},\ell_{i\mathrm{L}})^{\mathrm{T}}$,
$Q_{i\mathrm{L}}\equiv(u^{\prime}_{i\mathrm{L}},d^{\prime}_{i\mathrm{L}})^{\mathrm{T}}$,
${\tilde{\Phi}}_{2}\equiv i\sigma^{2}\Phi_{2}^{*}$, and $i,j=1,2,3$. The
Yukawa coupling matrices $\tilde{y}_{d}^{ij}$ and $\tilde{y}_{u}^{ij}$ can be
diagonalized by unitary matrices, which transform the gauge eigenstates
$u^{\prime}_{i}$ and $d^{\prime}_{i}$ into the mass eigenstates $u_{i}$ and
$d_{i}$. We remark that due to the similarity of the Yukawa couplings in the
quark sectors and the smallness of the ones in the leptonic sectors, many of
the following analyses for the type-I (type-II) case can be cast to the lepton
specific (flipped) case.

##  III Experimental bounds

In our analyses, we carry out random scans in the parameter space. The
following 12 parameters are adopted as the free parameters:

$\lambda_{1},~{}~{}\lambda_{2},~{}~{}\lambda_{3},~{}~{}\lambda_{4},~{}~{}\lambda_{5},~{}~{}\lambda_{S},~{}~{}\kappa_{1},~{}~{}\kappa_{2},~{}~{}\tan\beta,~{}~{}m_{\chi},~{}~{}v_{s},~{}~{}m_{12}^{2}.$
(17)

Each parameter point in the scans should be tested by existed experimental
bounds.

Firstly, we require that $m_{h_{i}}^{2}$ ($i=1,2,3$), $m_{a}^{2}$, and
$m_{H^{+}}^{2}$ should be positive to guarantee physical scalar masses.
Moreover, in order to ensure that the scalar potential is bounded from below,
the following conditions from copositivity criteria Klimenko:1984qx ;
Kannike:2012pe  should be satisfied:

$\displaystyle\lambda_{1}\geq 0,\quad\lambda_{2}\geq 0,\quad\lambda_{S}\geq
0,$ (18) $\displaystyle
a_{12}\equiv\lambda_{3}+\sqrt{\lambda_{1}\lambda_{2}}\geq 0,\quad
a^{\prime}_{12}\equiv\lambda_{3}+\lambda_{4}-|\lambda_{5}|+\sqrt{\lambda_{1}\lambda_{2}}\geq
0,$ (19) $\displaystyle\quad
a_{13}\equiv\kappa_{1}+\sqrt{\lambda_{1}\lambda_{S}}\geq 0,\quad
a_{23}\equiv\kappa_{2}+\sqrt{\lambda_{2}\lambda_{S}}\geq 0,$ (20)
$\displaystyle\sqrt{\lambda_{1}\lambda_{2}\lambda_{S}}+\lambda_{3}\sqrt{\lambda_{S}}+\kappa_{1}\sqrt{\lambda_{2}}+\kappa_{2}\sqrt{\lambda_{1}}+\sqrt{2a_{12}a_{13}a_{23}}\geq
0,$ (21)
$\displaystyle\sqrt{\lambda_{1}\lambda_{2}\lambda_{S}}+(\lambda_{3}+\lambda_{4}-|\lambda_{5}|)\sqrt{\lambda_{S}}+\kappa_{1}\sqrt{\lambda_{2}}+\kappa_{2}\sqrt{\lambda_{1}}+\sqrt{2a^{\prime}_{12}a_{13}a_{23}}\geq
0.$ (22)

Furthermore, we require one of $h_{i}$ acting as the SM-like Higgs boson with
a mass within the $3\sigma$ range of the measured value $m_{h}=125.18\pm
0.16~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$ Tanabashi:2018oca . The numerical
tool Lilith 2 Bernon:2015hsa ; Kraml:2019sis  is used to test whether the SM-
like Higgs boson is consistent with LHC run 1 and run 2 Higgs measurements
from ATLAS and CMS. Parameter points excluded by the data at $95\%$ confidence
level (C.L.) are abandoned.

Although FCNCs have been forbidden at tree level, they can arise from loop
corrections. In particular, the loops involving the charged Higgs boson
$H^{\pm}$ significantly contribute to the FCNC $B$-meson decays, depending on
$m_{H^{\pm}}$ and $\tan\beta$. The analysis by the Gfitter Group
Haller:2018nnx  shows that the strongest constraint on the type-I (type-II)
Yukawa couplings comes from the measurement of the FCNC decay
$B_{d}\to\mu^{+}\mu^{-}$ ($B_{s}\to\mu^{+}\mu^{-}$ and $B\to X_{s}\gamma$). We
further reject the parameter points that are excluded at 95% C.L. by these
flavor physics constraints.

Then we impose the constraints from DM phenomenology. We utilize FeynRules 2
Alloul:2013bka  and the MadGraph5_aMC@NLO Alwall:2014hca  plugin MadDM 3
Ambrogi:2018jqj  to calculate the prediction of the DM relic abundance. The
observed value of the relic abundance from the Planck experiment is given by
$\Omega_{\mathrm{DM}}h^{2}=0.1200\pm 0.0012$ Aghanim:2018eyx , where
$\Omega_{\mathrm{DM}}$ is the ratio of the DM energy density to the critical
density of the Universe and $h$ is the Hubble constant in unit of
$100~{}$\mathrm{k}\mathrm{m}~{}\mathrm{s}^{-1}~{}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$.
Only the parameter points predicting the observed relic abundance are
preserved. MadDM 3 is also used to compute the DM annihilation cross section
$\left<\sigma_{\mathrm{ann}}v\right>_{\mathrm{d}}$ with an average velocity of
$2\text{\times}{10}^{-5}$, which is corresponding to DM annihilation processes
at dwarf spheroidal galaxies. The 95% C.L. upper limits on
$\left<\sigma_{\mathrm{ann}}v\right>$ in the $b\bar{b}$ channel from the
$\gamma$-ray observations of dwarf galaxies by the Fermi-LAT satellite
experiment and the MAGIC Cherenkov telescopes Ahnen:2016qkx  are employed to
test the parameter points.

Below, we study the effective potential, cosmological phase transitions, and
gravitational waves for the parameter points surviving from all the
experimental bounds above.

##  IV Effective potential

In order to investigate the cosmological phase transitions in the model, we
need to construct the effective potential. We assume that only the $CP$-even
neutral scalar fields $\rho_{1}$, $\rho_{2}$, and $s$ can develop VEVs in the
cosmological history. The effective potential is then expressed as a function
of the classical background fields $\tilde{\rho}_{1}$, $\tilde{\rho}_{2}$, and
$\tilde{s}$.

The tree-level effective potential in terms of the classical fields derived
from Eqs. (II) and (2) is

$\displaystyle V_{\mathrm{0}}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})$
$\displaystyle=\frac{m_{11}^{2}}{2}\tilde{\rho}_{1}^{2}+\frac{m_{22}^{2}}{2}\tilde{\rho}_{2}^{2}-\frac{2m_{S}^{2}+m^{\prime
2}_{S}}{4}\tilde{s}^{2}-m_{12}^{2}\tilde{\rho}_{1}\tilde{\rho}_{2}$
$\displaystyle\quad+\frac{\lambda_{1}}{8}\tilde{\rho}_{1}^{4}+\frac{\lambda_{2}}{8}\tilde{\rho}_{2}^{4}+\frac{\lambda_{S}}{8}\tilde{s}^{4}+\frac{\lambda_{345}}{4}\tilde{\rho}_{1}^{2}\tilde{\rho}_{2}^{2}+\frac{\kappa_{1}}{4}\tilde{\rho}_{1}^{2}\tilde{s}^{2}+\frac{\kappa_{2}}{4}\tilde{\rho}_{2}^{2}\tilde{s}^{2}.$
(23)

Here, $m_{11}^{2}$, $m_{22}^{2}$, and $m_{S}^{2}$ should be expressed as in
Eqs. (6), (7), and (8), respectively.

At zero temperature, the one-loop effective potential $V_{1}$ receives the
Coleman-Weinberg terms Coleman:1973jx  in the $\overline{\mathrm{MS}}$
renormalization scheme Quiros:1999jp ,

$V_{1}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=\frac{1}{64\pi^{2}}\sum\limits_{i}{{n_{i}}\tilde{m}_{i}^{4}\left({\ln\frac{\tilde{m}_{i}^{2}}{\mu^{2}}-C_{i}}\right)},$
(24)

where the sum runs over all the particles $i$ coupling to the classical
fields, and $\tilde{m}_{i}^{2}$ are the corresponding particle masses squared
in terms of the classical fields. For the SM fermions, we only take into
account the top and bottom quark contributions, and neglect all the other much
smaller Yukawa couplings. Hence, all the particles we include in the
calculations are

$h_{1},~{}h_{2},~{}h_{3},~{}a,~{}H^{\pm},~{}G^{0},~{}G^{\pm},~{}\chi,~{}W^{\pm},~{}Z,~{}\gamma,~{}t,~{}b.$
(25)

Although the photon $\gamma$ would not contribute to Eq. (24), its
longitudinal mode can contribute to the daisy potential $V_{\mathrm{D}}$,
which will be discussed below. $\mu$ is the renormalization scale. For
transverse gauge bosons, $C_{i}=1/2$, while for longitudinal gauge bosons,
scalar bosons and fermions, $C_{i}=3/2$. $n_{i}$ count the degrees of freedom
of the particles, given by

$\displaystyle n_{h_{i}}$ $\displaystyle=$ $\displaystyle
n_{a}=n_{G^{0}}=n_{\chi}=n_{Z_{\mathrm{L}}}=n_{\gamma_{\mathrm{L}}}=1,$ (26)
$\displaystyle\quad n_{H^{\pm}}$ $\displaystyle=$ $\displaystyle
n_{G^{\pm}}=n_{Z_{\mathrm{T}}}=n_{\gamma_{\mathrm{T}}}=n_{W^{\pm}_{\mathrm{L}}}=2,$
(27) $\displaystyle n_{W^{\pm}_{\mathrm{T}}}$ $\displaystyle=$ $\displaystyle
4,\quad n_{t}=n_{b}=-12,$ (28)

where the minus signs for $n_{t}$ and $n_{b}$ characterize the feature of
fermion loops. The subscripts L and T denote the longitudinal and transverse
polarizations of the gauge bosons.

In terms of the classical background fields $\tilde{\rho}_{1}$,
$\tilde{\rho}_{2}$, and $\tilde{s}$, the elements of the symmetric mass-
squared matrix $\tilde{M}_{h}^{2}$ for the $CP$-even neutral scalar bosons are
derived as

$\displaystyle\tilde{M}_{h,11}^{2}$ $\displaystyle=$ $\displaystyle
m_{11}^{2}+\frac{1}{2}(3\lambda_{1}\tilde{\rho}_{1}^{2}+\lambda_{345}\tilde{\rho}_{2}^{2}+\kappa_{1}\tilde{s}^{2}),$
(29) $\displaystyle\tilde{M}_{h,22}^{2}$ $\displaystyle=$ $\displaystyle
m_{22}^{2}+\frac{1}{2}(\lambda_{345}\tilde{\rho}_{1}^{2}+3\lambda_{2}\tilde{\rho}_{2}^{2}+\kappa_{2}\tilde{s}^{2}),$
(30) $\displaystyle\tilde{M}_{h,33}^{2}$ $\displaystyle=$ $\displaystyle-
m_{S}^{2}-\frac{1}{2}(m^{\prime
2}_{S}-\kappa_{1}\tilde{\rho}_{1}^{2}-\kappa_{2}\tilde{\rho}_{2}^{2}-3\lambda_{S}\tilde{s}^{2}),$
(31) $\displaystyle\tilde{M}_{h,12}^{2}$ $\displaystyle=$ $\displaystyle-
m_{12}^{2}+\lambda_{345}\tilde{\rho}_{1}\tilde{\rho}_{2},\quad\tilde{M}_{h,13}^{2}=\kappa_{1}\tilde{\rho}_{1}\tilde{s},\quad\tilde{M}_{h,23}^{2}=\kappa_{2}\tilde{\rho}_{2}\tilde{s}.$
(32)

The mass-squared matrix for the $CP$-odd neutral scalar bosons is

$\tilde{M}_{0}^{2}=\begin{pmatrix}m_{11}^{2}+(\lambda_{1}\tilde{\rho}_{1}^{2}+\hat{\lambda}_{345}\tilde{\rho}_{2}^{2}+\kappa_{1}\tilde{s}^{2})/2&-m_{12}^{2}+\lambda_{5}\tilde{\rho}_{1}\tilde{\rho}_{2}\\\
-m_{12}^{2}+\lambda_{5}\tilde{\rho}_{1}\tilde{\rho}_{2}&m_{22}^{2}+(\lambda_{2}\tilde{\rho}_{2}^{2}+\hat{\lambda}_{345}\tilde{\rho}_{1}^{2}+\kappa_{2}\tilde{s}^{2})/2\end{pmatrix},$
(33)

with $\hat{\lambda}_{345}\equiv\lambda_{3}+\lambda_{4}-\lambda_{5}$, while the
mass-squared matrix for the charged scalar bosons is

$\tilde{M}_{+}^{2}=\begin{pmatrix}m_{11}^{2}+(\lambda_{1}\tilde{\rho}_{1}^{2}+\lambda_{3}\tilde{\rho}_{2}^{2}+\kappa_{1}\tilde{s}^{2})/2&-m_{12}^{2}+(\lambda_{4}+\lambda_{5})\tilde{\rho}_{1}\tilde{\rho}_{2}/2\\\
-m_{12}^{2}+(\lambda_{4}+\lambda_{5})\tilde{\rho}_{1}\tilde{\rho}_{2}/2&m_{22}^{2}+(\lambda_{2}\tilde{\rho}_{2}^{2}+\lambda_{3}\tilde{\rho}_{1}^{2}+\kappa_{2}\tilde{s}^{2})/2\end{pmatrix}.$
(34)

The eigenvalues of these matrices give the masses squared of the scalar
bosons, i.e.,

$\displaystyle\mathrm{eigenvalues}(\tilde{M}_{h}^{2})$ $\displaystyle=$
$\displaystyle\left\\{\tilde{m}_{h_{1}}^{2},~{}\tilde{m}_{h_{2}}^{2},~{}\tilde{m}_{h_{3}}^{2}\right\\},$
(35) $\displaystyle\mathrm{eigenvalues}(\tilde{M}_{0}^{2})$ $\displaystyle=$
$\displaystyle\left\\{\tilde{m}_{G^{0}}^{2},~{}\tilde{m}_{a}^{2}\right\\},$
(36) $\displaystyle\mathrm{eigenvalues}(\tilde{M}_{+}^{2})$ $\displaystyle=$
$\displaystyle\left\\{\tilde{m}_{G^{\pm}}^{2},~{}\tilde{m}_{H^{\pm}}^{2}\right\\}.$
(37)

The masses squared of the DM candidate $\chi$ is obtained as

$\tilde{m}_{\chi}^{2}=-m_{S}^{2}+\frac{1}{2}(m^{\prime
2}_{S}+\kappa_{1}\tilde{\rho}_{1}^{2}+\kappa_{2}\tilde{\rho}_{2}^{2}+\lambda_{S}\tilde{s}^{2}).$
(38)

The mass squared of the $W^{\pm}$ boson is given by

$\tilde{m}_{W^{\pm}_{\mathrm{L}}}^{2}=\tilde{m}_{W^{\pm}_{\mathrm{T}}}^{2}=\frac{g^{2}}{4}(\tilde{\rho}_{1}^{2}+\tilde{\rho}_{2}^{2}).$
(39)

The mass-squared matrix of the $B$ and $W^{3}$ gauge fields is

$\tilde{M}^{2}_{W^{3},B}=\frac{1}{4}(\tilde{\rho}_{1}^{2}+\tilde{\rho}_{2}^{2})\begin{pmatrix}g^{2}&-gg^{\prime}\\\
-gg^{\prime}&g^{\prime 2}\end{pmatrix}.$ (40)

After diagonalization, the masses squared of the $Z$ boson and the photon are

$\displaystyle\tilde{m}_{Z_{\mathrm{L}}}^{2}$ $\displaystyle=$
$\displaystyle\tilde{m}_{Z_{\mathrm{T}}}^{2}=\frac{1}{4}(g^{2}+g^{\prime
2})(\tilde{\rho}_{1}^{2}+\tilde{\rho}_{2}^{2}),$ (41)
$\displaystyle\tilde{m}_{\gamma_{\mathrm{T}}}^{2}$ $\displaystyle=$
$\displaystyle\tilde{m}_{\gamma_{\mathrm{L}}}^{2}=0.$ (42)

For the type-I Yukawa couplings, the masses squared of the top and bottom
quarks are

$\tilde{m}_{t}^{2}=\frac{y_{t}^{2}}{2\sin^{2}\beta}\,\tilde{\rho}_{2}^{2},\quad\tilde{m}_{b}^{2}=\frac{y_{b}^{2}}{2\sin^{2}\beta}\,\tilde{\rho}_{2}^{2},$
(43)

where the couplings $y_{t}=\sqrt{2}m_{t}/v$ and $y_{b}=\sqrt{2}m_{b}/v$ are
defined the same as in the SM. For the type-II Yukawa couplings, the masses
squared become

$\tilde{m}_{t}^{2}=\frac{y_{t}^{2}}{2\sin^{2}\beta}\,\tilde{\rho}_{2}^{2},\quad\tilde{m}_{b}^{2}=\frac{y_{b}^{2}}{2\cos^{2}\beta}\,\tilde{\rho}_{1}^{2}.$
(44)

Notice that loop corrections generally shift the values of the VEVs as well as
the renormalized mass-squared matrix of the $CP$-even neutral scalar bosons.
To keep them intact, we introduce the following counterterms Cline:2011mm ;
Basler:2016obg ,

$\displaystyle V_{\mathrm{CT}}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})$
$\displaystyle=\delta m_{1}^{2}\tilde{\rho}_{1}^{2}+\delta
m_{2}^{2}\tilde{\rho}_{2}^{2}+\delta
m_{s}^{2}{\tilde{s}^{2}}+\delta{\lambda_{1}}\tilde{\rho}_{1}^{4}+\delta{\lambda_{2}}\tilde{\rho}_{2}^{4}+\delta{\lambda_{s}}{\tilde{s}^{4}}$
$\displaystyle\quad+\delta\lambda_{12}\tilde{\rho}_{1}^{2}\tilde{\rho}_{2}^{2}+\delta\lambda_{1s}\tilde{\rho}_{1}^{2}\tilde{s}^{2}+\delta\lambda_{2s}\tilde{\rho}_{2}^{2}\tilde{s}^{2}.$
(45)

The nine counterterm coefficients are determined by the following nine
equations at
$(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=(v_{1},v_{2},v_{s})$,

$\displaystyle\frac{\partial V_{\mathrm{CT}}}{\partial\tilde{\rho}_{1}}$
$\displaystyle=-\frac{\partial V_{1}}{\partial\tilde{\rho}_{1}},\quad$
$\displaystyle\frac{\partial V_{\mathrm{CT}}}{\partial\tilde{\rho}_{2}}$
$\displaystyle=-\frac{\partial V_{1}}{\partial\tilde{\rho}_{2}},\quad$
$\displaystyle\frac{\partial V_{1}}{\partial\tilde{s}}$
$\displaystyle=-\frac{\partial V_{1}}{\partial\tilde{s}},$ (46)
$\displaystyle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\tilde{\rho}_{1}^{2}}$
$\displaystyle=-\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{1}^{2}},\quad$
$\displaystyle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\tilde{\rho}_{2}^{2}}$
$\displaystyle=-\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}^{2}},\quad$
$\displaystyle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\tilde{s}^{2}}$
$\displaystyle=-\frac{\partial^{2}V_{1}}{\partial\tilde{s}^{2}},\quad$ (47)
$\displaystyle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\tilde{\rho}_{2}\partial\tilde{\rho}_{1}}$
$\displaystyle=-\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}\partial\tilde{\rho}_{1}},\quad$
$\displaystyle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\tilde{s}\partial\tilde{\rho}_{1}}$
$\displaystyle=-\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{1}},\quad$
$\displaystyle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\tilde{s}\partial\tilde{\rho}_{2}}$
$\displaystyle=-\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{2}}.$
(48)

The masses of the Nambu-Goldstone bosons $G^{0}$ and $G^{\pm}$ vanish at
$(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=(v_{1},v_{2},v_{s})$ in the
Landau gauge, inducing logarithmic IR divergence terms in Eqs. (47) and (48)
proportional to

$\frac{\partial\tilde{m}_{G}^{2}}{\partial\phi_{i}}\frac{\partial\tilde{m}_{G}^{2}}{\partial\phi_{j}}\ln\frac{\tilde{m}_{G}^{2}}{\mu^{2}},\quad\phi_{i}=\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s}.$
(49)

This problem is due to the ill-defined renormalized Higgs boson masses at
$p^{2}=0$ with massless Nambu-Goldstone modes, and one can fix it by setting
the momenta of the Higgs bosons on shell Cline:1996mga ; Casas:1994us .
Similar problems exist in the effective potential with higher loops, and more
details can be found in Refs. Elias-Miro:2014pca ; Martin:2014bca . An
approximate treatment is to give an IR cutoff $\Lambda_{\mathrm{IR}}$ to the
Nambu-Goldstone boson masses Cline:2011mm , i.e., to set
$\tilde{m}_{G^{0}}^{2}=\tilde{m}_{G^{\pm}}^{2}=\Lambda_{\mathrm{IR}}^{2}$ in
the logarithms at
$(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=(v_{1},v_{2},v_{s})$. Here, we
take $\Lambda_{\mathrm{IR}}$ to be the mass of the SM-like Higgs boson.
Solving Eqs. (46)–(48), we obtain

$\displaystyle\delta m_{1}^{2}$ $\displaystyle=$
$\displaystyle-\frac{3}{4v_{1}}\frac{\partial
V_{1}}{\partial\tilde{\rho}_{1}}+\frac{1}{4}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{1}^{2}}+\frac{v_{2}}{4v_{1}}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}\partial\tilde{\rho}_{1}}+\frac{v_{s}}{4v_{1}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{1}},$
(50) $\displaystyle\delta m_{2}^{2}$ $\displaystyle=$
$\displaystyle-\frac{3}{4v_{2}}\frac{\partial
V_{1}}{\partial\tilde{\rho}_{2}}+\frac{1}{4}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}^{2}}+\frac{v_{1}}{4v_{2}}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}\partial\tilde{\rho}_{1}}+\frac{v_{s}}{4v_{2}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{2}},$
(51) $\displaystyle\delta m_{s}^{2}$ $\displaystyle=$
$\displaystyle-\frac{3}{4v_{s}}\frac{\partial
V_{1}}{\partial\tilde{s}}+\frac{1}{4}\frac{\partial^{2}V_{1}}{\partial\tilde{s}^{2}}+\frac{v_{1}}{4v_{s}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{1}}+\frac{v_{2}}{4v_{s}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{2}},$
(52) $\displaystyle\delta\lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{8v_{1}^{3}}\frac{\partial
V_{1}}{\partial\tilde{\rho}_{1}}-\frac{1}{8v_{1}^{2}}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{1}^{2}},\quad\delta\lambda_{2}=\frac{1}{8v_{2}^{3}}\frac{\partial
V_{1}}{\partial\tilde{\rho}_{2}}-\frac{1}{8v_{2}^{2}}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}^{2}},$
(53) $\displaystyle\delta\lambda_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{8v_{s}^{3}}\frac{\partial
V_{1}}{\partial\tilde{s}}-\frac{1}{8v_{s}^{2}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}^{2}},\quad\delta\lambda_{12}=-\frac{1}{4v_{1}v_{2}}\frac{\partial^{2}V_{1}}{\partial\tilde{\rho}_{2}\partial\tilde{\rho}_{1}},$
(54) $\displaystyle\delta\lambda_{1s}$ $\displaystyle=$
$\displaystyle-\frac{1}{4v_{1}v_{s}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{1}},\quad\delta\lambda_{2s}=-\frac{1}{4v_{2}v_{s}}\frac{\partial^{2}V_{1}}{\partial\tilde{s}\partial\tilde{\rho}_{2}},$
(55)

at $(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=(v_{1},v_{2},v_{s})$.

Thermal corrections to the effective potential are crucial for studying the
EWPT. The one-loop finite-temperature effective potential Dolan:1973qd  can be
expressed as

$V_{\mathrm{1T}}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s},T)=\frac{T^{4}}{2\pi^{2}}\left[{\sum\limits_{i=\mathrm{bosons}}{n_{i}}{J_{\mathrm{B}}}\left(\frac{\tilde{m}_{i}^{2}}{T^{2}}\right)}+\sum\limits_{i={t,b}}{n_{i}}{J_{\mathrm{F}}}\left(\frac{\tilde{m}_{i}^{2}}{T^{2}}\right)\right],$
(56)

where $T$ is the temperature and the functions $J_{\mathrm{B}}$ and
$J_{\mathrm{F}}$ are defined as

$\displaystyle J_{\mathrm{B}}(x)$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}{{y^{2}}\ln\left({1-e^{-\sqrt{y^{2}+x}}}\right)dy},$
(57) $\displaystyle J_{\mathrm{F}}(x)$ $\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}{{y^{2}}\ln\left({1+e^{-\sqrt{y^{2}+x}}}\right)dy}.$
(58)

We also consider the daisy diagrams, which can be significant. The
corressponding contribution to the effective potential can be estimated by
Carrington:1991hz ; Arnold:1992rz

$V_{\mathrm{D}}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s},T)=-\frac{T}{12\pi}\sum\limits_{i=\mathrm{bosons}}{n_{i}}\left[(\bar{m}_{i}^{2})^{3/2}-(\tilde{m}_{i}^{2})^{3/2}\right].$
(59)

$\bar{m}_{i}^{2}$ are the field-dependent boson masses squared with thermal
corrections in the high-temperature limit and can be derived by

$\bar{m}_{i}^{2}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s},T)=\mathrm{eigenvalues}[\tilde{M}_{X}^{2}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})+\Pi_{X}(T)],$
(60)

where $\tilde{M}_{X}^{2}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})$
represents the mass-squared matrices or masses squared in terms of the
classical fields, and $\Pi_{X}(T)$ denotes the thermal corrections to
$\tilde{M}_{X}^{2}$. The subleading off-diagonal elements of $\Pi_{X}(T)$ can
be neglected Carrington:1991hz ; Blinov:2015vma . The diagonal elements of
$\Pi_{X}(T)$ for the scalar bosons are derived as

$\displaystyle\Pi_{h,11}$ $\displaystyle=$
$\displaystyle\Pi_{0,11}=\Pi_{+,11}=\frac{T^{2}}{48}(9{g^{2}}+3{g^{\prime}}^{2}+12\lambda_{1}+8\lambda_{3}+4\lambda_{4}+4\kappa_{1}+y_{1}),$
(61) $\displaystyle\Pi_{h,22}$ $\displaystyle=$
$\displaystyle\Pi_{0,22}=\Pi_{+,22}=\frac{T^{2}}{48}(9{g^{2}}+3{g^{\prime}}^{2}+12\lambda_{2}+8\lambda_{3}+4\lambda_{4}+4\kappa_{2}+y_{2}),$
(62) $\displaystyle\Pi_{h,33}$ $\displaystyle=$
$\displaystyle\Pi_{\chi}=\frac{T^{2}}{6}(\lambda_{S}+\kappa_{1}+\kappa_{2}).$
(63)

Here, $y_{1}$ and $y_{2}$ are the contributions from the Yukawa couplings. For
the type-I and -II cases, they are given by

Type I: $\displaystyle\quad y_{1}=0,\quad
y_{2}=\frac{12(y_{t}^{2}+y_{b}^{2})}{\sin^{2}\beta},$ (64) Type II:
$\displaystyle\quad y_{1}=\frac{12y_{b}^{2}}{\cos^{2}\beta},\quad
y_{2}=\frac{12y_{t}^{2}}{\sin^{2}\beta}.$ (65)

The thermal corrections to the electroweak gauge bosons are

$\displaystyle\Pi_{W^{\pm}_{\mathrm{L}}}$ $\displaystyle=$
$\displaystyle\Pi_{W^{3}_{\mathrm{L}}}=2g^{2}T^{2},$ (66)
$\displaystyle\Pi_{B_{\mathrm{L}}}$ $\displaystyle=$ $\displaystyle
2{g^{\prime}}^{2}T^{2},$ (67) $\displaystyle\Pi_{W^{\pm}_{\mathrm{T}}}$
$\displaystyle=$
$\displaystyle\Pi_{Z_{\mathrm{T}}}=\Pi_{\gamma_{\mathrm{T}}}=0.$ (68)

Note that $\Pi_{W^{3}_{\mathrm{L}}}$ and $\Pi_{B_{\mathrm{L}}}$ are the
corrections to the diagonal elements of $\tilde{M}^{2}_{W^{3},B}$ in Eq. (40).

Finally, we obtain the total effective potential111Discussions on theoretical
uncertainties in perturbative calculations of the effective potential can be
found in Ref. Croon:2020cgk .

$V_{\mathrm{eff}}(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s},T)=V_{0}+V_{1}+V_{\mathrm{CT}}+V_{\mathrm{1T}}+V_{\mathrm{D}}.$
(69)

##  V Phase transitions

Based on the effective potential constructed in the previous section, we can
study its evolution with temperature. At sufficiently high temperatures, the
effective potential is minimized at the origin
$(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=(0,0,0)$, implying the
restoration of the electroweak gauge symmetry. As the Universe cools down,
extra minima appear. In particular, if there are two coexisted minima
separated by a high barrier, strong FOPT could take place and result in a
stochastic GW background. We utilize the numerical package CosmoTransitions
Wainwright:2011kj  to analyze the phase transitions. For each parameter point
in the random scans, we verify whether or not the minimum
$(\tilde{\rho}_{1},\tilde{\rho}_{2},\tilde{s})=(v_{1},v_{2},v_{s})$ is the
global one of the zero-temperature effective potential. The parameter points
that fail this test are rejected. Then we use CosmoTransitions to trace the
temperature evolution of the local minima.

In this model, the three classical $CP$-even neutral scalar fields would
develop VEVs, typically leading to multi-step cosmological phase transitions.
In Fig. 1, we demonstrate the temperature evolution of multiple phases for a
benchmark point (BP), whose parameters can be found in the BP3 column of Table
1 in Sec. VII. In the plots, $v_{1}(T)$, $v_{2}(T)$, and $v_{s}(T)$ are the
$T$-dependent values of the classical fields $\tilde{\rho}_{1}$,
$\tilde{\rho}_{2}$, and $\tilde{s}$ at the local minima of the effective
potential. The red, green, and blue lines indicate the positions of three
local minima.

(a) Minima in the $T$-$v_{1}(T)$ plane. (b) Minima in the $T$-$v_{2}(T)$
plane. (c) Minima in the $T$-$v_{s}(T)$ plane.

Figure 1: Temperature evolution of the positions of the minima in the axes
$v_{1}(T)$, $v_{2}(T)$, and $v_{s}(T)$ for BP3. The red, green, and blue lines
denote three local minima. The vertical dashed lines indicate the critical,
nucleation, and percolation temperatures $T_{\mathrm{c}}$, $T_{\mathrm{n}}$,
and $T_{\mathrm{p}}$.

At $T\gtrsim 460~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, the system stays at the
red minimum with $\big{(}v_{1}(T),v_{2}(T),v_{s}(T)\big{)}=(0,0,0)$,
respecting the electroweak gauge symmetry. At $T\simeq
460~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, a second-order phase transition
occurs and the system turns into the green minimum, where $\tilde{s}$ develops
a nonzero VEV. At $T\simeq 148~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, the blue
minimum appears, accompanied with a barrier that separates it from the green
minimum. These two minima coexist till the zero temperature.

The effective potential at the blue minimum is higher than at the green
minimum until the critical temperature $T_{\mathrm{c}}\simeq
119~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$. Below $T_{\mathrm{c}}$, the green
minimum becomes a metastable state, i.e., a “false vacuum”. The system finally
undergoes a FOPT through quantum tunneling and turns into the blue minimum, or
the “true vacuum”. Such a FOPT nucleates bubbles, inside which the system is
trapped at the true vacuum. In this FOPT, $v_{1}(T)$ and $v_{2}(T)$ increase
from zero to $\mathcal{O}(100)~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, while
$v_{s}(T)$ slightly decreases. At zero temperature, the true vacuum satisfies
$\big{(}v_{1}(0),v_{2}(0),v_{s}(0)\big{)}=(v_{1},v_{2},v_{s})$.

In our parameter scans, we usually find that $v_{s}(T)$ does not evolve
synchronously with $v_{1}(T)$ and $v_{2}(T)$, probably due to the less
couplings of the singlet field to other fields compared with the two Higgs
doublets. Typically, $v_{s}(T)$ becomes nonzero much earlier than the
conventional EWPT epoch via a second-order or first-order phase transition.
$\tilde{\rho}_{1}$ and $\tilde{\rho}_{2}$ then gain VEVs in an subsequent
phase transition, which could be a strong FOPT similar to those in the
conventional two-Higgs-doublet models Dorsch:2017nza ; Bernon:2017jgv .

Below we discuss the dynamics of the FOPTs. The bubble nucleation rate per
unit time and unit volume is given by Linde:1980tt ; Linde:1981zj

$\Gamma\sim AT^{4}e^{-S},$ (70)

where $A$ is an $\mathcal{O}(1)$ constant and
$S=\mathrm{min}\\{S_{4},S_{3}/T\\}$. $S_{4}$ and $S_{3}$ are the Euclidean
actions of the scalar fields for $O(4)$\- and $O(3)$-symmetric bubbles,
respectively. The three-dimensional action $S_{3}$ can be simplified to

$S_{3}=4\pi\int_{0}^{\infty}dr\,r^{2}\left[{\frac{1}{2}{\frac{d\phi_{i}}{dr}}{\frac{d\phi_{i}}{dr}}+V_{\mathrm{eff}}(\phi_{i},T)}\right],$
(71)

where $r$ is the radius of the bubble.
$\phi_{i}(r)=\big{(}\tilde{\rho}_{1}(r),\tilde{\rho}_{2}(r),\tilde{s}(r)\big{)}$
is given by the bounce solution of the equations of motion

$\frac{d^{2}\phi_{i}}{dr^{2}}+\frac{2}{r}\frac{d\phi_{i}}{dr}=\frac{\partial
V_{\mathrm{eff}}}{\partial\phi_{i}}$ (72)

with boundary conditions

$\left.\frac{d\phi_{i}}{dr}\right|_{r=0}=0,\quad\phi_{i}(\infty)=\phi_{i}^{\mathrm{false}},$
(73)

where $\phi_{i}^{\mathrm{false}}$ is the field configuration of the false
vacuum.

(a) $S_{3}/T$ and $S_{4}$. (b) Nucleation rate $\Gamma$.

Figure 2: The actions $S_{3}/T$ and $S_{4}$ (a) and the nucleation rate
$\Gamma$ (b) as functions of the temperature $T$ for BP3. The dashed lines
denote the critical, nucleation, and percolation temperatures
$T_{\mathrm{c}}$, $T_{\mathrm{n}}$, and $T_{\mathrm{p}}$.

We present $S_{4}$ and $S_{3}/T$ as functions of the temperature for BP3 in
Fig. 2(a), as well as the corresponding nucleation rate $\Gamma$ in Fig. 2(b).
We find that $S_{3}/T$ is the smaller one until temperatures below $\sim
5~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$. The minimal of $S_{3}/T$ is reached at
$T\sim 40~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$. As the Universe cools down,
below the critical temperature $T_{\mathrm{c}}$, the nucleation rate increases
before the peak around $T\sim 40~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, and then
decreases.

The bubbles are actually nucleated at the nucleation temperature
$T_{\mathrm{n}}$, where the nucleation probability for a single bubble within
a Hubble volume reaches $\mathcal{O}(1)$. Thus, $T_{\mathrm{n}}$ can be
estimated by Moreno:1998bq

$\int_{t_{\mathrm{c}}}^{t_{\mathrm{n}}}{dt\,\frac{\Gamma}{H^{3}}}=\int_{T_{\mathrm{n}}}^{T_{\mathrm{c}}}{dT\,\frac{\Gamma}{H^{4}T}}=1,$
(74)

where $H$ is the Hubble rate, and $t_{\mathrm{c}}$ and $t_{\mathrm{n}}$ denote
the critical and nucleation times, respectively. Note that the differential
relation between the time $t$ and the temperature $T$ is $dt=-(HT)^{-1}dT$ in
the radiation-dominated epoch.

Below the nucleation temperature $T_{\mathrm{n}}$, an increasing number of
bubbles thrive and collide with each other. The maximum of bubble collisions
that remarkably produces stochastic GWs is expected to be reached when
percolation occurs Leitao:2012tx . In order to evaluate the percolation time
$t_{\mathrm{p}}$, we need to estimate the fraction of space that still remains
in the false vacuum at time $t$, which can be computed by Guth:1979bh ;
Guth:1981uk

$P(t)=\exp\left[-\frac{4\pi}{3}\int_{t_{\mathrm{c}}}^{t}{dt^{\prime}\,\Gamma(t^{\prime})\,a^{3}(t^{\prime})\,r^{3}(t,t^{\prime})}\right],$
(75)

where $a(t^{\prime})$ is the scale factor. $r(t,t^{\prime})$ is the comoving
radius of a bubble growing from $t^{\prime}$ to $t$, given by

$r(t,t^{\prime})=\int_{t^{\prime}}^{t}{d\tau}\,\frac{v_{\mathrm{w}}}{a(\tau)},$
(76)

where $v_{\mathrm{w}}$ is the velocity of the bubble wall. For randomly
distributed spherical bubbles with equal size in the three-dimensional space,
the percolation threshold is reached when the fraction of space converted to
the true vacuum, $1-P(t)$, increases to $\sim 0.29$ Shante_1971 ; Rintoul_1997
. Thus, the percolation time $t_{\mathrm{p}}$ can be derived by requiring
$P(t_{\mathrm{p}})\simeq 0.71$, with the corresponding temperature
$T_{\mathrm{p}}$ characterizing GW production from FOPTs Leitao:2012tx ;
Kobakhidze:2017mru ; Ellis:2018mja ; Wang:2020jrd .

FOPTs are able to release latent heat from the vacuum energy, which drives the
expansion of the bubbles and also converts into the thermal and bulk kinetic
energies of the plasma Steinhardt:1981ct ; Kamionkowski:1993fg ;
Espinosa:2010hh . The density of the released vacuum energy is given by
Enqvist:1991xw

$\rho_{\mathrm{vac}}=V_{\mathrm{eff}}(\phi_{i}^{\mathrm{false}},T)-V_{\mathrm{eff}}(\phi_{i}^{\mathrm{true}},T)-T\frac{\partial}{\partial
T}[V_{\mathrm{eff}}(\phi_{i}^{\mathrm{false}},T)-V_{\mathrm{eff}}(\phi_{i}^{\mathrm{true}},T)],$
(77)

where $\phi_{i}^{\mathrm{true}}$ is the field configuration of the true
vacuum. It is useful to define a dimensionless strength parameter

$\alpha\equiv\frac{\rho_{\mathrm{vac}}}{\rho_{\mathrm{rad}}},$ (78)

with $\rho_{\mathrm{rad}}={{\pi^{2}}{g_{*}}T^{4}}/{30}$ the radiation energy
density in the plasma. $g_{*}$ is the effective relativistic degrees of
freedom in the plasma.

The expansion of the bubbles depends on the interactions between the bubble
walls and the plasma, analogous to chemical combustion in a relativistic fluid
Steinhardt:1981ct . Hydrodynamic analyses show that bubble propagation have
diverse modes, including Jouguet detonations, weak detonations, subsonic
deflagrations, supersonic deflagrations (hybrid), and runway bubble walls
Espinosa:2010hh . Thus, it is difficult to completely work out the bubble wall
velocity $v_{\mathrm{w}}$. For Jouguet detonations, the Chapman-Jouguet
condition leads to a wall velocity of Steinhardt:1981ct

$v_{\mathrm{CJ}}=\frac{1+\sqrt{3\alpha^{2}+2\alpha}}{\sqrt{3}(1+\alpha)}.$
(79)

This is a typical assumption when evaluating GW signals.

Expanding the action $S$ around the time $t^{\prime}=t_{\mathrm{n}}$ or
$t^{\prime}=t_{\mathrm{p}}$, we have

$S(t)\simeq
S(t^{\prime})-\beta(t-t^{\prime})+\mathcal{O}[(t-t^{\prime})^{2}],$ (80)

where

$\beta\equiv\left.-\frac{dS}{dt}\right|_{t=t^{\prime}}=\left.\left(HT\frac{dS}{dT}\right)\right|_{T=T^{\prime}}$
(81)

can be roughly understood as the inverse time duration of the phase transition
Kosowsky:1991ua . For the electroweak FOPTs in which we are interested, the
derivative $dS/dT$ is positive at $T=T^{\prime}$, leading to positive $\beta$.
In addition, $S_{3}/T$ is typically smaller than $S_{4}$ at $T=T^{\prime}$. In
order to conveniently compare the phase transition time scale $\beta^{-1}$ and
the cosmological expansion time scale $H^{-1}$, we define a dimensionless
quantity

$\tilde{\beta}(T^{\prime})\equiv\frac{\beta(T^{\prime})}{H(T^{\prime})}.$ (82)

Based on Eqs. (74) and (75), further calculations show that the nucleation and
percolation temperatures $T_{\mathrm{n}}$ and $T_{\mathrm{p}}$ can be
approximately determined by Huber:2007vva

$\displaystyle\frac{S_{3}(T_{\mathrm{n}})}{T_{\mathrm{n}}}$
$\displaystyle\simeq$ $\displaystyle
141.5-2\ln\frac{g_{*}}{100}-4\ln\frac{T_{\mathrm{n}}}{100~{}$\mathrm{G}\mathrm{e}\mathrm{V}$}-\ln\frac{\tilde{\beta}(T_{\mathrm{n}})}{100},$
(83) $\displaystyle\frac{S_{3}(T_{\mathrm{p}})}{T_{\mathrm{p}}}$
$\displaystyle\simeq$ $\displaystyle
132.0-2\ln\frac{g_{*}}{100}-4\ln\frac{T_{\mathrm{p}}}{100~{}$\mathrm{G}\mathrm{e}\mathrm{V}$}-4\ln\frac{\tilde{\beta}(T_{\mathrm{p}})}{100}+3\ln
v_{\mathrm{w}}.$ (84)

For BP3, the nucleation temperature is $T_{\mathrm{n}}\simeq
64~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, while the percolation temperature
assuming $v_{\mathrm{w}}=v_{\mathrm{CJ}}$ is slightly lower,
$T_{\mathrm{p}}\simeq 60~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, as denoted in
Figs. 1 and 2.

##  VI Gravitational wave spectra

Electroweak FOPTs could induce significant perturbations of the Friedmann-
Robertson-Walker metric and produce stochastic GWs around the mHz band. Two
key parameters relevant to the relic GW spectrum are $\alpha$ and
$\tilde{\beta}$ evaluated at the time $t_{*}$ when GWs are produced. There are
three coexisting GW sources at a FOPT, namely bubble collisions, sound waves,
and magnetohydrodynamic (MHD) turbulence Binetruy:2012ze ; Caprini:2015zlo ;
Caprini:2019egz ; Hindmarsh:2020hop . Denoting $\Omega_{\mathrm{GW}}$ to be
the present GW energy density per logarithmic frequency interval divided by
the critical density, we separate the contributions from the three sources as

$\Omega_{\mathrm{GW}}h^{2}=\Omega_{\mathrm{col}}h^{2}+\Omega_{\mathrm{sw}}h^{2}+\Omega_{\mathrm{turb}}h^{2}.$
(85)

(a) Bubble collisions

The nucleated bubbles expand and finally collide with each other. Their
collisions break the spherical symmetry and generate gravitational waves
Kosowsky:1991ua . This process can be well described by the envelope
approximation Kosowsky:1991ua ; Kosowsky:1992rz ; Kosowsky:1992vn . Numerical
simulations for bubble collisions in the thermal plasma Kamionkowski:1993fg ;
Huber:2008hg  show that the resulting GW spectrum at present can be
approximated by

$\Omega_{\mathrm{col}}h^{2}=1.67\times 10^{-5}\
\frac{0.11v_{\mathrm{w}}^{3}}{(0.42+v_{\mathrm{w}}^{2})\tilde{\beta}^{2}}\left(\frac{\kappa_{\phi}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}\frac{3.8(f/f_{\mathrm{col}})^{2.8}}{1+2.8(f/f_{\mathrm{col}})^{3.8}},$
(86)

where $g_{*}$ is evaluated at $T=T_{*}$, the temperature corresponding to
$t=t_{*}$. The peak frequency of the spectrum can be modeled as Huber:2008hg

$f_{\mathrm{col}}=\frac{0.62\,\tilde{\beta}h_{*}}{1.8-0.1v_{\mathrm{w}}+v_{\mathrm{w}}^{2}}.$
(87)

The redshift of the frequency has been taken into account by the factor

$h_{*}=\frac{a(t_{*})H(t_{*})}{a(t_{0})}=1.65\times 10^{-5}\ \mathrm{Hz}\
\frac{T_{*}}{100\ \mathrm{GeV}}\left(\frac{g_{*}}{100}\right)^{1/6},$ (88)

which is the inverse Hubble time at $t=t_{*}$ redshifted to today ($t=t_{0}$).
The efficiency factor $\kappa_{\phi}$ characterizes the fraction of the
available vacuum energy converted into the gradient energy of the scalar
fields.

(b) Sound waves

The explosive bubble expansion in the plasma induces a sound shell around the
bubble wall. After the bubble collisions, the sound shells propagate into the
fluid as sound waves, which become a significant GW source Hindmarsh:2013xza ;
Hindmarsh:2015qta ; Hindmarsh:2017gnf . This source lasts until the sound
waves are disrupted by the development of nonlinear shocks and turbulence
Hindmarsh:2017gnf ; Ellis:2018mja ; Ellis:2019oqb ; Ellis:2020awk . Therefore,
the duration of the sound wave source can be determined by the nonlinearity
timescale estimated as Ellis:2020awk

$\tau_{\mathrm{nl}}\sim\frac{(8\pi)^{1/3}v_{\mathrm{w}}}{\tilde{\beta}H_{*}}\sqrt{\frac{4(1+\alpha)}{3\kappa_{v}\alpha}},$
(89)

where $H_{*}\equiv H(t_{*})$ is the Hubble rate at $t=t_{*}$ and $\kappa_{v}$
is the fraction of the available vacuum energy converted into the kinetic
energy of the fluid bulk motion. For Jouguet detonations,
$v_{\mathrm{w}}=v_{\mathrm{CJ}}$, and $\kappa_{v}$ can be approximated by
Espinosa:2010hh

$\kappa_{v}^{\mathrm{CJ}}=\frac{\sqrt{\alpha}}{0.135+\sqrt{0.98+\alpha}}.$
(90)

In an expanding radiation-dominated Universe, the finite duration of the sound
wave source leads to a suppression factor Guo:2020grp

$\Upsilon=1-\frac{1}{\sqrt{1+2\tau_{\mathrm{nl}}H_{*}}}.$ (91)

Thus, the GW spectrum contributed by the sound waves is given by
Hindmarsh:2017gnf ; Guo:2020grp

$\Omega_{\mathrm{sw}}h^{2}=1.17\times 10^{-6}\ \frac{\Upsilon
v_{\mathrm{w}}}{\tilde{\beta}}\left(\frac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{1/3}\left(\frac{f}{f_{\mathrm{sw}}}\right)^{3}\left(\frac{7}{4+3f^{2}/f_{\mathrm{sw}}^{2}}\right)^{7/2},$
(92)

where the peak frequency is estimated to be Hindmarsh:2017gnf

$f_{\mathrm{sw}}=\frac{0.54\tilde{\beta}h_{*}}{v_{\mathrm{w}}}.$ (93)

(c) MHD turbulence

Bubble collisions can stir up turbulence in the fluid, as the energy injection
to the plasma results in an extremely high Reynolds number Kamionkowski:1993fg
. Since the plasma is fully ionized, the magnetic field, along with the
velocity field, should be considered, leading to MHD turbulence
Kosowsky:2001xp . It takes several Hubble times for the MHD turbulence to
decay, and the stochastic GWs arise continuously during this period
Caprini:2009yp . The corresponding GW spectrum can be fitted as Caprini:2009yp
; Caprini:2015zlo

$\Omega_{\mathrm{turb}}h^{2}=3.35\times 10^{-4}\
\frac{v_{\mathrm{w}}}{\tilde{\beta}}\left(\frac{\kappa_{\mathrm{turb}}\alpha}{1+\alpha}\right)^{3/2}\left(\frac{100}{g_{*}}\right)^{1/3}\frac{(f/f_{\mathrm{turb}})^{3}}{(1+f/f_{\mathrm{turb}})^{11/3}(1+8\pi
f/{h_{*}})},$ (94)

with

$f_{\mathrm{turb}}=\frac{3.5\tilde{\beta}h_{*}}{2v_{\mathrm{w}}}.$ (95)

Based on the suggestion from simulations, we optimistically set
$\kappa_{\mathrm{turb}}\simeq 0.1\kappa_{v}$ Hindmarsh:2015qta ;
Caprini:2015zlo .

In general, the contribution from the sound waves dominates in the GW spectrum
Hindmarsh:2020hop . Moreover, $\kappa_{\phi}$ is typically negligible, except
for runaway bubble walls Espinosa:2010hh ; Caprini:2015zlo ; Ellis:2019oqb .
Thus, we omit the contribution from the bubble collisions in the following
calculations.

##  VII Numerical analyses

We perform random scans with the model parameters in the ranges of
$v_{s}\in[10,~{}1000]~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$,
$m_{\chi}\in[58,~{}800]~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$,
$|m^{2}_{12}|\in[1,~{}500^{2}]~{}$\mathrm{G}\mathrm{e}\mathrm{V}^{2}$$,
$\tan\beta\in[0.5,~{}20]$, $\lambda_{1}$, $\lambda_{2}$, $|\lambda_{3}|$,
$|\lambda_{4}|$, $|\lambda_{5}|$, $\lambda_{S}\in[0.8,~{}8]$, and
$|\kappa_{1}|$, $|\kappa_{2}|\in[0.01,~{}8]$ for the two types of Yukawa
couplings. We assume that the prior probabilities for the random parameters
follow uniform distributions in the logarithmic scale. The parameter points
are required to pass all the experimental constraints described in Sec. III,
as well as to cause a FOPT.

Note that positive $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{S}$ are
required to satisfy the bounded-from-below conditions (18). A negative $v_{s}$
would be totally equivalent to a positive one due to the $Z_{2}$ symmetry
$S\to-S$. Besides, a parameter point with $\tan\beta$ and $m_{12}^{2}$ is
equivalent to one with $-\tan\beta$ and $-m_{12}^{2}$, since the potential
respects the $Z_{2}$ symmetry $\Phi_{1}\to-\Phi_{1}$ or $\Phi_{2}\to-\Phi_{2}$
expect for the soft breaking quadratic terms with $m_{12}^{2}$. Thus, we can
just take positive $v_{s}$ and $\tan\beta$ in the scans, while $m_{12}^{2}$,
$\lambda_{3}$, $\lambda_{4}$, $\lambda_{5}$, $\kappa_{1}$, and $\kappa_{2}$
can be either positive or negative. In addition, the $v_{s}$ range of
$10~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$ to
$1~{}$\mathrm{T}\mathrm{e}\mathrm{V}$$ ensures that $S$ has a VEV near the
electroweak scale, and the interplay between $S$ and the two Higgs doublets
could be important.

Figure 3: Contours of the peak amplitudes of the GW spectra in the
$\alpha$-$\tilde{\beta}^{-1}$ plane assuming Jouguet detonations. Purple and
green points denote the parameter points for type-I and type-II Yukawa
couplings, respectively. Four BPs are also indicated.

The strength of the stochastic GW signals from the FOPT depend on $\alpha$ and
$\beta$. Larger $\alpha$ implies a stronger FOPT, while smaller $\beta$
corresponds to a longer FOPT time duration. Consequently, larger $\alpha$ and
$\tilde{\beta}^{-1}$ lead to stronger GW signals, as implied in Eqs. (86),
(92) and (94). For the surviving parameter points, we calculate the resulting
values of $\alpha$ and $\tilde{\beta}^{-1}$, and then project the points in
the $\alpha$-$\tilde{\beta}^{-1}$ plane, as presented in Fig. 3. The purple
and green points are corresponding to type-I and type-II Yukawa couplings,
respectively. The parameter points lie in the ranges of
$10^{-4}\lesssim\alpha\lesssim 0.3$ and
$10^{-7}\lesssim\tilde{\beta}^{-1}\lesssim 0.02$. The relic GW spectra for the
parameter points are further evaluated, assuming Jouguet detonations.

We introduce $\hat{\Omega}_{\mathrm{GW}}h^{2}$ to denote the peak amplitudes
of the GW spectra. The contours of $\hat{\Omega}_{\mathrm{GW}}h^{2}$ are
demonstrated in Fig. 3 and we can easily read off the GW signal strengths of
the parameter points from the plot. The strongest GW signal we find reaches up
to $\hat{\Omega}_{\mathrm{GW}}h^{2}\sim 10^{-11}$.

(a) Type-I Yukawa couplings. (b) Type-II Yukawa couplings.

Figure 4: Peak amplitudes of the total GW spectra versus frequency for the
parameter points with type-I (a) and type-II (b) Yukawa couplings assuming
Jouguet detonations. Sensitivity curves for the future space-based GW
interferometers LISA Audley:2017drz , Tianqin Mei:2020lrl , Taiji Guo:2018npi
, BBO Cutler:2005qq , and ultimate DECIGO Kudoh:2005as  are also plotted. The
color axes denote the LISA signal-to-noise ratio
$\mathrm{SNR}_{\mathrm{LISA}}$ for the parameter points with
$\mathrm{SNR}_{\mathrm{LISA}}>10$. The gray points yield
$\mathrm{SNR}_{\mathrm{LISA}}<10$.

Fig. 4 illustrates $\hat{\Omega}_{\mathrm{GW}}h^{2}$ versus the peak frequency
$f$ for the parameter points. For comparison, we also plot the sensitivity
curves for the future space-based interferometers LISA Audley:2017drz ,
TianQin Mei:2020lrl , Taiji Guo:2018npi , BBO Cutler:2005qq , and DECIGO
Kudoh:2005as . Some of the curves are converted from the sensitivity on
amplitude spectral density or characteristic strain. The conversions of the
related quantities can be found in, e.g., Ref. Moore:2014lga . The DECIGO
curve we adopt here is the ultimate sensitivity that is only limited by
quantum noises, and it can be regarded as an observational limitation
Kudoh:2005as .

The GWs produced by FOPTs become an isotropic and stochastic background in the
present Universe. The detectability of the GW signals in the space-based
interferometers increases with the practical observation time $\mathcal{T}$.
The signal-to-noise ratio can be defined as Thrane:2013oya ; Caprini:2015zlo

$\mathrm{SNR}\equiv\sqrt{\mathcal{T}\int_{f_{\min}}^{f_{\max}}{\frac{\Omega_{\mathrm{GW}}^{2}(f)}{\Omega_{\mathrm{sens}}^{2}(f)}\,df}},$
(96)

where $\Omega_{\mathrm{sens}}(f)$ is the sensitivity of the experiment. Below
we take the practical observation time $\mathcal{T}=9.46\times
10^{7}~{}\text{s}$ (3 years) for LISA Caprini:2019egz , Taiji, and TianQin.
The signal $\Omega_{\mathrm{GW}}(f)$ can be detected if the corresponding
$\mathrm{SNR}$ is larger than a signal-to-noise ratio threshold
$\mathrm{SNR}_{\mathrm{thr}}$. For the six (four) link configuration of LISA,
the threshold is $\mathrm{SNR}_{\mathrm{thr}}=10~{}(50)$ Caprini:2015zlo . We
find that some parameter points yield the LISA signal-to-noise ratio
$\mathrm{SNR}_{\mathrm{LISA}}>10$ and could be probed by LISA. We denote
$\mathrm{SNR}_{\mathrm{LISA}}$ for them as the color axes in Fig. 4, with the
remaining gray points corresponding to $\mathrm{SNR}_{\mathrm{LISA}}<10$. The
next-generation plans aiming at
$f\sim\mathcal{O}(0.1)~{}$\mathrm{H}\mathrm{z}$$, like BBO and DECIGO, may
probe much more parameter points.

Table 1: Detailed information for four benchmark points. | BP1 | BP2 | BP3 | BP4
---|---|---|---|---
Type | I | I | II | II
$v_{s}~{}(\mathrm{GeV})$ | $542.40$ | $384.26$ | $64.987$ | $138.82$
$m_{\chi}~{}(\mathrm{GeV})$ | $117.88$ | $78.191$ | $134.03$ | $76.678$
$m^{2}_{12}~{}(\mathrm{GeV}^{2})$ | $2.0210\times 10^{4}$ | $1.5876\times 10^{2}$ | $1.7696\times 10^{5}$ | $1.5042\times 10^{5}$
$\tan\beta$ | $2.8616$ | $3.2654$ | $0.91655$ | $1.1732$
$\lambda_{1}$ | $2.1496$ | $2.1882$ | $1.5297$ | $0.87839$
$\lambda_{2}$ | $0.80887$ | $0.85479$ | $1.2074$ | $0.80222$
$\lambda_{3}$ | $2.3925$ | $2.2628$ | $1.5741$ | $2.8002$
$\lambda_{4}$ | $3.0027$ | $1.4715$ | $5.3967$ | $4.4643$
$\lambda_{5}$ | $-6.2187$ | $-4.0567$ | $-7.8556$ | $-7.5755$
$\lambda_{S}$ | $3.4048$ | $2.5502$ | $6.0689$ | $4.8644$
$\kappa_{1}$ | $-1.4852$ | $1.0295$ | $0.80378$ | $-0.38075$
$\kappa_{2}$ | $1.1727$ | $-1.2142$ | $-0.83745$ | $-0.14591$
$m_{h_{1}}~{}(\mathrm{GeV})$ | $125.11$ | $91.459$ | $125.38$ | $124.87$
$m_{h_{2}}~{}(\mathrm{GeV})$ | $282.02$ | $124.77$ | $158.83$ | $307.56$
$m_{h_{3}}~{}(\mathrm{GeV})$ | $1014.5$ | $641.83$ | $650.98$ | $582.08$
$m_{a}~{}(\mathrm{GeV})$ | $664.75$ | $496.49$ | $911.87$ | $874.04$
$m_{H^{\pm}}~{}(\mathrm{GeV})$ | $402.96$ | $280.94$ | $655.60$ | $631.66$
$\left<\sigma_{\mathrm{ann}}v\right>_{\mathrm{d}}~{}(\mathrm{cm}^{3}/\mathrm{s})$ | $1.30\times 10^{-26}$ | $3.68\times 10^{-27}$ | $1.72\times 10^{-26}$ | $6.82\times 10^{-27}$
$\alpha$ | $0.240$ | $0.160$ | $0.181$ | $0.346$
$\tilde{\beta}^{-1}$ | $1.33\times 10^{-2}$ | $4.02\times 10^{-3}$ | $7.71\times 10^{-3}$ | $2.15\times 10^{-2}$
$T_{\mathrm{p}}~{}(\mathrm{GeV})$ | $55.3$ | $74.9$ | $60.2$ | $47.2$
$\mathrm{SNR_{LISA}}$ | $96.6$ | $37.7$ | $60.1$ | $120$
$\mathrm{SNR_{Taiji}}$ | $83.3$ | $23.9$ | $42.3$ | $155$
$\mathrm{SNR_{TianQin}}$ | $5.50$ | $2.39$ | $3.07$ | $9.20$

For a closer look at the results, we choose four benchmark points, whose
detailed information is listed in Table 1. BP1 and BP2 (BP3 and BP4)
correspond to the type-I (type-II) Yukawa couplings. In these BPs, the masses
of the Higgs bosons $h_{1,2,3}$, $a$, and $H^{\pm}$ are all below 1 TeV, while
the mass of the DM candidate is less than 140 GeV. The SM-like Higgs boson is
$h_{2}$ in BP2, while it is $h_{1}$ in the rest BPs. The DM annihilation cross
sections $\left<\sigma_{\mathrm{ann}}v\right>_{\mathrm{d}}$ at dwarf galaxies
predicted by the BPs are below $2\text{\times}{10}^{-26}$
$\mathrm{c}\mathrm{m}^{3}\mathrm{/}\mathrm{s}$, beyond the reach of Fermi-LAT
and MAGIC Ahnen:2016qkx .

Figure 5: GW spectra for four benchmark points assuming Jouguet detonations.

For the four BPs, percolation of the FOPT occurs in
$47~{}$\mathrm{G}\mathrm{e}\mathrm{V}$\lesssim T_{\mathrm{p}}\lesssim
75~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$, with $\alpha$ ranging from 0.16 to
0.35 and $\tilde{\beta}^{-1}$ ranging from $4\text{\times}{10}^{-3}$ to
$2.2\text{\times}{10}^{-2}$. Assuming Jouguet detonations, we derive the GW
spectra for these BPs, as presented in Fig. 5. The BPs are also indicated in
Figs. 3 and 4. We find that the GW signal strengths decrease according to the
order of BP4, BP1, BP3, and BP2, reflecting the descending order of $\alpha$.
In Table 1, we also list the signal-to-noise ratios $\mathrm{SNR_{LISA}}$,
$\mathrm{SNR_{Taiji}}$, and $\mathrm{SNR_{TianQin}}$ for the LISA, Taiji, and
TianQin experiments, respectively. LISA and Taiji look promising to detect all
BPs, while TianQin may probe BP4 with a sightly longer observation time.

(a) Type-I Yukawa couplings. (b) Type-II Yukawa couplings.

Figure 6: Parameter points projected in the $\tan\beta$-$m_{H}^{\pm}$ plane
for type-I (a) and type-II (b) Yukawa couplings. Yellow regions are excluded
at 95% C.L. by the FCNC bounds from the Gfitter global fit Haller:2018nnx .
The color axes denote the peak amplitude of the GW spectrum
$\hat{\Omega}_{\mathrm{GW}}h^{2}$ for Jouguet detonations.

In order to show the most important flavor constraints mentioned in Sec. III,
we plot our parameter points confronting the FCNC bounds. We have adopted the
data from the Gfitter global fit Haller:2018nnx  to reject parameter points.
In Fig. 6, the parameter points for the two types of Yukawa couplings are
projected in the $\tan\beta$-$m_{H}^{\pm}$ plane, with the color axis
indicating the peak amplitude of the GW spectrum,
$\hat{\Omega}_{\mathrm{GW}}h^{2}$.

For the type-I case in Fig. 6(a), the most stringent FCNC bound comes from the
LHCb and CMS measurements of $B_{d}\to\mu^{+}\mu^{-}$ CMS:2014xfa ;
Aaij:2017vad , excluding a region with $\tan\beta\lesssim 3$. For the type-II
case in Fig. 6(b), the bounds from the observations of $B\to X_{s}\gamma$
Amhis:2016xyh ; Misiak:2006zs ; Misiak:2015xwa  and $B_{s}\to\mu^{+}\mu^{-}$
CMS:2014xfa ; Aaij:2017vad  exclude a region with $m_{H^{\pm}}\lesssim
750~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$. Thus, the FCNC constraints remove
small $\tan\beta$ and light $H^{\pm}$ for type-I and type-II Yukawa couplings,
respectively. We remark that strong GW signals typically favor small
$m_{H^{\pm}}$. There are many parameter points with $m_{H^{\pm}}\lesssim
600~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$ in the type-I case leading to large
$\hat{\Omega}_{\mathrm{GW}}h^{2}$. In the type-II case, the
$m_{H^{\pm}}\lesssim 600~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$ region is
basically excluded by the FCNC constraints, and hence it is more difficult to
achieve strong GW signals.

(a) Type-I Yukawa couplings. (b) Type-II Yukawa couplings.

Figure 7: Parameter points projected in the
$m_{\chi}$-$\left<\sigma_{\mathrm{ann}}v\right>_{\mathrm{d}}$ plane for type-I
(a) and type-II (b) Yukawa couplings. Dashed lines denote the 95% C.L. upper
limits on the DM annihilation cross section from the $\gamma$-ray observations
of dwarf galaxies by Fermi-LAT and MAGIC Ahnen:2016qkx .

In Fig. 7, we project the parameter points in the
$m_{\chi}$-$\left<\sigma_{\mathrm{ann}}v\right>_{\mathrm{d}}$ plane. Although
all these parameter points are required to predict the observed relic
abundance, $\left<\sigma_{\mathrm{ann}}v\right>_{\mathrm{d}}$ can deviate from
the canonical annihilation cross section $3\times
10^{-26}~{}$\mathrm{c}\mathrm{m}^{3}\mathrm{/}\mathrm{s}$$, due to the
velocity dependence of $\left<\sigma_{\mathrm{ann}}v\right>$. One reason is
that some annihilation channels kinematically forbidden at low velocities
could be opened at the freeze-out epoch, and another reason is the resonance
effects Griest:1990kh . The pile-up of points around $m_{\chi}\sim
78~{}$\mathrm{G}\mathrm{e}\mathrm{V}$$ is mainly related to the annihilation
channel $\chi\chi\to W^{+}W^{-}$.

In the above numerical analyses, we have assumed the bubble propagation mode
to be Jouguet detonations with bubble wall velocity
$v_{\mathrm{w}}=v_{\mathrm{CJ}}$. Below, we study the effect of various bubble
propagation modes. In general, the dependence of $\kappa_{v}$ on
$v_{\mathrm{w}}$ and $\alpha$ can be found in Ref. Espinosa:2010hh . The sound
speed $c_{\mathrm{s}}$ in the relativistic plasma is very close to
$1/\sqrt{3}$. For $v_{\mathrm{w}}\ll c_{\mathrm{s}}$,
$v_{\mathrm{w}}=c_{\mathrm{s}}$, and $v_{\mathrm{w}}=1$, the vacuum energy
fraction converted into the bulk kinetic energy of the fluid $\kappa_{v}$ has
the following analytic approximations, based on fit.

$\displaystyle v_{\mathrm{w}}\ll c_{\mathrm{s}}:$
$\displaystyle\quad\kappa_{v}^{\mathrm{A}}=\frac{{6.9\alpha
v_{\mathrm{w}}^{6/5}}}{{1.36-0.037\sqrt{\alpha}+\alpha}}.$ (97) $\displaystyle
v_{\mathrm{w}}=c_{\mathrm{s}}:$
$\displaystyle\quad\kappa_{v}^{\mathrm{B}}=\frac{{{\alpha^{2/5}}}}{{0.017+{(0.997+\alpha)^{2/5}}}}.$
(98) $\displaystyle v_{\mathrm{w}}=1:$
$\displaystyle\quad\kappa_{v}^{\mathrm{D}}=\frac{\alpha}{{0.73+0.083\sqrt{\alpha}+\alpha}}.$
(99)

Furthermore, for subsonic deflagrations (${v_{\mathrm{w}}}<{c_{\mathrm{s}}}$),
supersonic deflagrations
(${c_{\mathrm{s}}}<{v_{\mathrm{w}}}<{v_{{\mathrm{CJ}}}}$), and detonations
(${v_{\mathrm{w}}}\gtrsim{v_{{\mathrm{CJ}}}}$), $\kappa_{v}$ is roughly given
by

$\displaystyle{\kappa_{v}}({v_{\mathrm{w}}}<{c_{\mathrm{s}}})$
$\displaystyle=$
$\displaystyle\frac{{c_{\mathrm{s}}^{11/5}\kappa_{v}^{\mathrm{A}}\kappa_{v}^{\mathrm{B}}}}{{(c_{\mathrm{s}}^{11/5}-v_{\mathrm{w}}^{11/5})\kappa_{v}^{\mathrm{B}}+{v_{\mathrm{w}}}c_{\mathrm{s}}^{6/5}\kappa_{v}^{\mathrm{A}}}},$
(100)
$\displaystyle{\kappa_{v}}({c_{\mathrm{s}}}<{v_{\mathrm{w}}}<{v_{{\mathrm{CJ}}}})$
$\displaystyle=$
$\displaystyle\kappa_{v}^{\mathrm{B}}+({v_{\mathrm{w}}}-{c_{\mathrm{s}}})\delta\kappa+{\left({\frac{{{v_{\mathrm{w}}}-{c_{\mathrm{s}}}}}{{{v_{{\mathrm{CJ}}}}-{c_{\mathrm{s}}}}}}\right)^{3}}\left[\kappa_{v}^{\mathrm{CJ}}-\kappa_{v}^{\mathrm{B}}-({v_{{\mathrm{CJ}}}}-{c_{\mathrm{s}}})\delta\kappa\right],$
(101) $\displaystyle{\kappa_{v}}({v_{\mathrm{w}}}\gtrsim{v_{{\mathrm{CJ}}}})$
$\displaystyle=$
$\displaystyle\frac{{{{({v_{{\mathrm{CJ}}}}-1)}^{3}}{{({v_{{\mathrm{CJ}}}}/{v_{\mathrm{w}}})}^{5/2}}\kappa_{v}^{\mathrm{CJ}}\kappa_{v}^{\mathrm{D}}}}{{[{{({v_{{\mathrm{CJ}}}}-1)}^{3}}-{{({v_{\mathrm{w}}}-1)}^{3}}]v_{{\mathrm{CJ}}}^{5/2}\kappa_{v}^{\mathrm{CJ}}+{{({v_{\mathrm{w}}}-1)}^{3}}\kappa_{v}^{\mathrm{D}}}},$
(102)

where $\delta\kappa=-0.9\ln[\sqrt{\alpha}/(1+\sqrt{\alpha})]$.

Figure 8: GW spectra for BP4 with several assumptions of $v_{\mathrm{w}}$.

According to these expressions, we derive the GW spectra for BP4 assuming the
bubble wall velocity $v_{\mathrm{w}}=0.05$, $0.2$, $0.72$, and $1$. These
spectra are demonstrated in Fig. 8, along with the previously obtained BP4
spectrum for $v_{\mathrm{w}}=v_{\mathrm{CJ}}=0.87$. Compared to Jouguet
detonations with $v_{\mathrm{w}}=v_{\mathrm{CJ}}$, supersonic deflagrations
with $v_{\mathrm{w}}=0.72$ lead to a stronger GW signal, which could be
properly tested by TianQin with $\mathrm{SNR}_{\mathrm{TianQin}}=15.8$. On the
other hand, subsonic deflagrations with $v_{\mathrm{w}}=0.05$ and $0.2$ give
much weaker GW signals.

##  VIII Summary

In this paper, we have studied the stochastic GW signals from electroweak
FOPTs in the model comprising the pNGB dark matter framework and two Higgs
doublets with type-I or type-II Yukawa couplings. The DM candidate is a pNGB
whose tree-level scattering off nucleons vanishes at zero momentum transfer,
evading the constraints from direct detection experiments. The three scalar
fields in the model have nonzero VEVs at zero temperature, which should be
developed from EWPTs at the early Universe. If such EWPTs are of strongly
first order, stochastic GWs could be effectively produced. The effective
potential has been carefully constructed with one-loop corrections at zero
temperature as well as thermal corrections, allowing us to carry out accurate
analyses on EWPTs.

We have performed random scans in the 12-dimensional parameter space, taking
into account the constraints from bounded from below conditions, LHC run 1 and
run 2 measurements of the 125 GeV Higgs boson, FCNC $B$-meson decays, the
Planck observation of the DM relic abundance, and the $\gamma$-ray
observations of dwarf galaxies by Fermi-LAT and MAGIC. The surviving parameter
points are also required to induce an electroweak FOPT. We have further
analyzed the characteristic temperatures, the phase transition strength, and
the characteristic time duration of the FOPT. Based on such information, the
resulting relic GW spectra from sound waves and MHD turbulence have been
evaluated.

Assuming that the bubble propagation mode is Jouguet detonations with
$v_{\mathrm{w}}=v_{\mathrm{CJ}}$, we have found that the FOPTs of some
parameter points could induce peak amplitudes of the GW spectra around
$10^{-13}\text{--}10^{-11}$, which could be well detected by the future space-
based GW interferometers LISA and Taiji. The next-generation GW
interferometers BBO and DECIGO are capable of probing much more parameter
points. We have noticed that a lighter charged Higgs boson $H^{\pm}$ in this
model is more probable to induce a strong GW signal. Since the FCNC
constraints on $m_{H^{\pm}}$ for type-II Yukawa couplings at large $\tan\beta$
are more stringent than those for type-I Yukawa couplings, the type-I case
typically leads to stronger GW signals.

We have also investigated the effects of different bubble propagation modes
with several values of the bubble wall velocity $v_{\mathrm{w}}$. For the
benchmark point BP4, supersonic deflagrations with $v_{\mathrm{w}}=0.72$ can
induce a stronger GW signal than Jouguet detonations with
$v_{\mathrm{w}}=v_{\mathrm{CJ}}$. In this optimistic case, BP4 could be well
tested by LISA, Taiji, and TianQin. Detonations with $v_{\mathrm{w}}=1$ lead
to a slightly weaker GW signal, while subsonic deflagrations with
$v_{\mathrm{w}}=0.2$ and $0.05$ result in much weaker signals.

###### Acknowledgements.

We thank Fa Peng Huang and Ligong Bian for helpful discussions. This work is
supported in part by the National Natural Science Foundation of China under
Grants No. 11805288, No. 11875327, No. 11905300, and No. 12005312, the China
Postdoctoral Science Foundation under Grant No. 2018M643282, the Fundamental
Research Funds for the Central Universities, and the Sun Yat-Sen University
Science Foundation.

## References

  * (1) G. Bertone, D. Hooper, and J. Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept. 405 (2005) 279–390, arXiv:hep-ph/0404175 [hep-ph].
  * (2) J. L. Feng, “Dark Matter Candidates from Particle Physics and Methods of Detection,” Ann. Rev. Astron. Astrophys. 48 (2010) 495–545, arXiv:1003.0904 [astro-ph.CO].
  * (3) B.-L. Young, “A survey of dark matter and related topics in cosmology,” Front. Phys.(Beijing) 12 (2017) 121201. [Erratum: Front. Phys.(Beijing)12,no.2,121202(2017)].
  * (4) LUX Collaboration, D. S. Akerib et al., “Results from a search for dark matter in the complete LUX exposure,” Phys. Rev. Lett. 118 (2017) 021303, arXiv:1608.07648 [astro-ph.CO].
  * (5) PandaX-II Collaboration, X. Cui et al., “Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment,” Phys. Rev. Lett. 119 (2017) 181302, arXiv:1708.06917 [astro-ph.CO].
  * (6) XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * (7) C. Gross, O. Lebedev, and T. Toma, “Cancellation Mechanism for Dark-Matter–Nucleon Interaction,” Phys. Rev. Lett. 119 (2017) 191801, arXiv:1708.02253 [hep-ph].
  * (8) D. Azevedo, M. Duch, B. Grzadkowski, D. Huang, M. Iglicki, and R. Santos, “One-loop contribution to dark-matter-nucleon scattering in the pseudo-scalar dark matter model,” JHEP 01 (2019) 138, arXiv:1810.06105 [hep-ph].
  * (9) K. Ishiwata and T. Toma, “Probing pseudo Nambu-Goldstone boson dark matter at loop level,” JHEP 12 (2018) 089, arXiv:1810.08139 [hep-ph].
  * (10) K. Huitu, N. Koivunen, O. Lebedev, S. Mondal, and T. Toma, “Probing pseudo-Goldstone dark matter at the LHC,” Phys. Rev. D 100 (2019) 015009, arXiv:1812.05952 [hep-ph].
  * (11) T. Alanne, M. Heikinheimo, V. Keus, N. Koivunen, and K. Tuominen, “Direct and indirect probes of Goldstone dark matter,” Phys. Rev. D 99 (2019) 075028, arXiv:1812.05996 [hep-ph].
  * (12) K. Kannike and M. Raidal, “Phase Transitions and Gravitational Wave Tests of Pseudo-Goldstone Dark Matter in the Softly Broken U(1) Scalar Singlet Model,” Phys. Rev. D 99 (2019) 115010, arXiv:1901.03333 [hep-ph].
  * (13) D. Karamitros, “Pseudo Nambu-Goldstone Dark Matter: Examples of Vanishing Direct Detection Cross Section,” Phys. Rev. D 99 (2019) 095036, arXiv:1901.09751 [hep-ph].
  * (14) J. M. Cline and T. Toma, “Pseudo-Goldstone dark matter confronts cosmic ray and collider anomalies,” Phys. Rev. D 100 (2019) 035023, arXiv:1906.02175 [hep-ph].
  * (15) X.-M. Jiang, C. Cai, Z.-H. Yu, Y.-P. Zeng, and H.-H. Zhang, “Pseudo-Nambu-Goldstone dark matter and two-Higgs-doublet models,” Phys. Rev. D 100 (2019) 075011, arXiv:1907.09684 [hep-ph].
  * (16) C. Arina, A. Beniwal, C. Degrande, J. Heisig, and A. Scaffidi, “Global fit of pseudo-Nambu-Goldstone Dark Matter,” JHEP 04 (2020) 015, arXiv:1912.04008 [hep-ph].
  * (17) Y. Abe, T. Toma, and K. Tsumura, “Pseudo-Nambu-Goldstone dark matter from gauged $U(1)_{B-L}$ symmetry,” JHEP 05 (2020) 057, arXiv:2001.03954 [hep-ph].
  * (18) N. Okada, D. Raut, and Q. Shafi, “Pseudo-Goldstone Dark Matter in gauged $B-L$ extended Standard Model,” arXiv:2001.05910 [hep-ph].
  * (19) S. Glaus, M. Mühlleitner, J. Müller, S. Patel, T. Römer, and R. Santos, “Electroweak Corrections in a Pseudo-Nambu Goldstone Dark Matter Model Revisited,” JHEP 12 (2020) 034, arXiv:2008.12985 [hep-ph].
  * (20) Y. Abe, T. Toma, and K. Yoshioka, “Non-thermal Production of PNGB Dark Matter and Inflation,” arXiv:2012.10286 [hep-ph].
  * (21) LIGO Scientific, Virgo Collaboration, B. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116 (2016) 061102, arXiv:1602.03837 [gr-qc].
  * (22) M. Jiang, L. Bian, W. Huang, and J. Shu, “Impact of a complex singlet: Electroweak baryogenesis and dark matter,” Phys. Rev. D 93 (2016) 065032, arXiv:1502.07574 [hep-ph].
  * (23) M. Chala, G. Nardini, and I. Sobolev, “Unified explanation for dark matter and electroweak baryogenesis with direct detection and gravitational wave signatures,” Phys. Rev. D 94 (2016) 055006, arXiv:1605.08663 [hep-ph].
  * (24) W. Chao, H.-K. Guo, and J. Shu, “Gravitational Wave Signals of Electroweak Phase Transition Triggered by Dark Matter,” JCAP 09 (2017) 009, arXiv:1702.02698 [hep-ph].
  * (25) A. Beniwal, M. Lewicki, J. D. Wells, M. White, and A. G. Williams, “Gravitational wave, collider and dark matter signals from a scalar singlet electroweak baryogenesis,” JHEP 08 (2017) 108, arXiv:1702.06124 [hep-ph].
  * (26) F. P. Huang and J.-H. Yu, “Exploring inert dark matter blind spots with gravitational wave signatures,” Phys. Rev. D 98 (2018) 095022, arXiv:1704.04201 [hep-ph].
  * (27) F. P. Huang and C. S. Li, “Probing the baryogenesis and dark matter relaxed in phase transition by gravitational waves and colliders,” Phys. Rev. D 96 (2017) 095028, arXiv:1709.09691 [hep-ph].
  * (28) A. Hektor, K. Kannike, and V. Vaskonen, “Modifying dark matter indirect detection signals by thermal effects at freeze-out,” Phys. Rev. D 98 (2018) 015032, arXiv:1801.06184 [hep-ph].
  * (29) I. Baldes and C. Garcia-Cely, “Strong gravitational radiation from a simple dark matter model,” JHEP 05 (2019) 190, arXiv:1809.01198 [hep-ph].
  * (30) E. Madge and P. Schwaller, “Leptophilic dark matter from gauged lepton number: Phenomenology and gravitational wave signatures,” JHEP 02 (2019) 048, arXiv:1809.09110 [hep-ph].
  * (31) A. Beniwal, M. Lewicki, M. White, and A. G. Williams, “Gravitational waves and electroweak baryogenesis in a global study of the extended scalar singlet model,” JHEP 02 (2019) 183, arXiv:1810.02380 [hep-ph].
  * (32) L. Bian and Y.-L. Tang, “Thermally modified sterile neutrino portal dark matter and gravitational waves from phase transition: The Freeze-in case,” JHEP 12 (2018) 006, arXiv:1810.03172 [hep-ph].
  * (33) L. Bian and X. Liu, “Two-step strongly first-order electroweak phase transition modified FIMP dark matter, gravitational wave signals, and the neutrino mass,” Phys. Rev. D 99 (2019) 055003, arXiv:1811.03279 [hep-ph].
  * (34) V. R. Shajiee and A. Tofighi, “Electroweak Phase Transition, Gravitational Waves and Dark Matter in Two Scalar Singlet Extension of The Standard Model,” Eur. Phys. J. C 79 (2019) 360, arXiv:1811.09807 [hep-ph].
  * (35) A. Mohamadnejad, “Gravitational waves from scale-invariant vector dark matter model: Probing below the neutrino-floor,” Eur. Phys. J. C 80 (2020) 197, arXiv:1907.08899 [hep-ph].
  * (36) K. Kannike, K. Loos, and M. Raidal, “Gravitational wave signals of pseudo-Goldstone dark matter in the $\mathbb{Z}_{3}$ complex singlet model,” Phys. Rev. D 101 (2020) 035001, arXiv:1907.13136 [hep-ph].
  * (37) A. Paul, B. Banerjee, and D. Majumdar, “Gravitational wave signatures from an extended inert doublet dark matter model,” JCAP 10 (2019) 062, arXiv:1908.00829 [hep-ph].
  * (38) N. Chen, T. Li, Y. Wu, and L. Bian, “Complementarity of the future $e^{+}e^{-}$ colliders and gravitational waves in the probe of complex singlet extension to the standard model,” Phys. Rev. D 101 (2020) 075047, arXiv:1911.05579 [hep-ph].
  * (39) B. Barman, A. Dutta Banik, and A. Paul, “Singlet-doublet fermionic dark matter and gravitational waves in a two-Higgs-doublet extension of the Standard Model,” Phys. Rev. D 101 (2020) 055028, arXiv:1912.12899 [hep-ph].
  * (40) C.-W. Chiang and B.-Q. Lu, “First-order electroweak phase transition in a complex singlet model with $\mathbb{Z}_{3}$ symmetry,” JHEP 07 (2020) 082, arXiv:1912.12634 [hep-ph].
  * (41) D. Borah, A. Dasgupta, K. Fujikura, S. K. Kang, and D. Mahanta, “Observable Gravitational Waves in Minimal Scotogenic Model,” JCAP 08 (2020) 046, arXiv:2003.02276 [hep-ph].
  * (42) Z. Kang and J. Zhu, “Scale-genesis by Dark Matter and Its Gravitational Wave Signal,” Phys. Rev. D 102 (2020) 053011, arXiv:2003.02465 [hep-ph].
  * (43) M. Pandey and A. Paul, “Gravitational Wave Emissions from First Order Phase Transitions with Two Component FIMP Dark Matter,” arXiv:2003.08828 [hep-ph].
  * (44) X.-F. Han, L. Wang, and Y. Zhang, “Dark matter, electroweak phase transition and gravitational wave in the type-II two-Higgs-doublet model with a singlet scalar field,” arXiv:2010.03730 [hep-ph].
  * (45) T. Alanne, N. Benincasa, M. Heikinheimo, K. Kannike, V. Keus, N. Koivunen, and K. Tuominen, “Pseudo-Goldstone dark matter: gravitational waves and direct-detection blind spots,” JHEP 10 (2020) 080, arXiv:2008.09605 [hep-ph].
  * (46) Y. Wang, C. S. Li, and F. P. Huang, “Complementary probe of dark matter blind spots by lepton colliders and gravitational waves,” arXiv:2012.03920 [hep-ph].
  * (47) T. Ghosh, H.-K. Guo, T. Han, and H. Liu, “Electroweak Phase Transition with an SU(2) Dark Sector,” arXiv:2012.09758 [hep-ph].
  * (48) W.-C. Huang, M. Reichert, F. Sannino, and Z.-W. Wang, “Testing the Dark Confined Landscape: From Lattice to Gravitational Waves,” arXiv:2012.11614 [hep-ph].
  * (49) W. Chao, X.-F. Li, and L. Wang, “Filtered pseudo-scalar dark matter and gravitational waves from first order phase transition,” arXiv:2012.15113 [hep-ph].
  * (50) A. Mazumdar and G. White, “Review of cosmic phase transitions: their significance and experimental signatures,” Rept. Prog. Phys. 82 (2019) 076901, arXiv:1811.01948 [hep-ph].
  * (51) LISA Collaboration, P. Amaro-Seoane et al., “Laser Interferometer Space Antenna,” arXiv:1702.00786 [astro-ph.IM].
  * (52) TianQin Collaboration, J. Luo et al., “TianQin: a space-borne gravitational wave detector,” Class. Quant. Grav. 33 (2016) 035010, arXiv:1512.02076 [astro-ph.IM].
  * (53) X.-C. Hu, X.-H. Li, Y. Wang, W.-F. Feng, M.-Y. Zhou, Y.-M. Hu, S.-C. Hu, J.-W. Mei, and C.-G. Shao, “Fundamentals of the orbit and response for TianQin,” Class. Quant. Grav. 35 (2018) 095008, arXiv:1803.03368 [gr-qc].
  * (54) J. Mei et al., “The TianQin project: current progress on science and technology,” PTEP 2020 (2020) , arXiv:2008.10332 [gr-qc].
  * (55) W.-R. Hu and Y.-L. Wu, “The Taiji Program in Space for gravitational wave physics and the nature of gravity,” Natl. Sci. Rev. 4 (2017) 685–686.
  * (56) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, “Taiji program: Gravitational-wave sources,” Int. J. Mod. Phys. A 35 (2020) 2050075, arXiv:1807.09495 [gr-qc].
  * (57) N. Seto, S. Kawamura, and T. Nakamura, “Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space,” Phys. Rev. Lett. 87 (2001) 221103, arXiv:astro-ph/0108011.
  * (58) H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, “Detecting a gravitational-wave background with next-generation space interferometers,” Phys. Rev. D 73 (2006) 064006, arXiv:gr-qc/0511145.
  * (59) C. Ungarelli, P. Corasaniti, R. Mercer, and A. Vecchio, “Gravitational waves, inflation and the cosmic microwave background: Towards testing the slow-roll paradigm,” Class. Quant. Grav. 22 (2005) S955–S964, arXiv:astro-ph/0504294.
  * (60) C. Cutler and J. Harms, “BBO and the neutron-star-binary subtraction problem,” Phys. Rev. D 73 (2006) 042001, arXiv:gr-qc/0511092.
  * (61) G. Dorsch, S. Huber, T. Konstandin, and J. No, “A Second Higgs Doublet in the Early Universe: Baryogenesis and Gravitational Waves,” JCAP 05 (2017) 052, arXiv:1611.05874 [hep-ph].
  * (62) X. Wang, F. P. Huang, and X. Zhang, “Gravitational wave and collider signals in complex two-Higgs doublet model with dynamical CP-violation at finite temperature,” Phys. Rev. D 101 (2020) 015015, arXiv:1909.02978 [hep-ph].
  * (63) R. Zhou and L. Bian, “Baryon asymmetry and detectable Gravitational Waves from Electroweak phase transition,” arXiv:2001.01237 [hep-ph].
  * (64) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept. 516 (2012) 1–102, arXiv:1106.0034 [hep-ph].
  * (65) I. F. Ginzburg and M. Krawczyk, “Symmetries of two Higgs doublet model and CP violation,” Phys. Rev. D 72 (2005) 115013, arXiv:hep-ph/0408011.
  * (66) S. L. Glashow and S. Weinberg, “Natural Conservation Laws for Neutral Currents,” Phys. Rev. D 15 (1977) 1958.
  * (67) E. Paschos, “Diagonal Neutral Currents,” Phys. Rev. D 15 (1977) 1966.
  * (68) K. G. Klimenko, “On Necessary and Sufficient Conditions for Some Higgs Potentials to Be Bounded From Below,” Theor. Math. Phys. 62 (1985) 58–65.
  * (69) K. Kannike, “Vacuum Stability Conditions From Copositivity Criteria,” Eur. Phys. J. C 72 (2012) 2093, arXiv:1205.3781 [hep-ph].
  * (70) Particle Data Group Collaboration, M. Tanabashi et al., “Review of Particle Physics,” Phys. Rev. D 98 (2018) 030001.
  * (71) J. Bernon and B. Dumont, “Lilith: a tool for constraining new physics from Higgs measurements,” Eur. Phys. J. C 75 (2015) 440, arXiv:1502.04138 [hep-ph].
  * (72) S. Kraml, T. Q. Loc, D. T. Nhung, and L. D. Ninh, “Constraining new physics from Higgs measurements with Lilith: update to LHC Run 2 results,” SciPost Phys. 7 (2019) 052, arXiv:1908.03952 [hep-ph].
  * (73) J. Haller, A. Hoecker, R. Kogler, K. Mönig, T. Peiffer, and J. Stelzer, “Update of the global electroweak fit and constraints on two-Higgs-doublet models,” Eur. Phys. J. C 78 (2018) 675, arXiv:1803.01853 [hep-ph].
  * (74) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, “FeynRules 2.0 - A complete toolbox for tree-level phenomenology,” Comput. Phys. Commun. 185 (2014) 2250–2300, arXiv:1310.1921 [hep-ph].
  * (75) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * (76) F. Ambrogi, C. Arina, M. Backovic, J. Heisig, F. Maltoni, L. Mantani, O. Mattelaer, and G. Mohlabeng, “MadDM v.3.0: a Comprehensive Tool for Dark Matter Studies,” Phys. Dark Univ. 24 (2019) 100249, arXiv:1804.00044 [hep-ph].
  * (77) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO].
  * (78) MAGIC, Fermi-LAT Collaboration, M. L. Ahnen et al., “Limits to Dark Matter Annihilation Cross-Section from a Combined Analysis of MAGIC and Fermi-LAT Observations of Dwarf Satellite Galaxies,” JCAP 1602 (2016) 039, arXiv:1601.06590 [astro-ph.HE].
  * (79) S. R. Coleman and E. J. Weinberg, “Radiative Corrections as the Origin of Spontaneous Symmetry Breaking,” Phys. Rev. D 7 (1973) 1888–1910.
  * (80) M. Quiros, “Finite temperature field theory and phase transitions,” in ICTP Summer School in High-Energy Physics and Cosmology, pp. 187–259.  1, 1999, arXiv:hep-ph/9901312.
  * (81) J. M. Cline, K. Kainulainen, and M. Trott, “Electroweak Baryogenesis in Two Higgs Doublet Models and B meson anomalies,” JHEP 11 (2011) 089, arXiv:1107.3559 [hep-ph].
  * (82) P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt, and A. Wlotzka, “Strong First Order Electroweak Phase Transition in the CP-Conserving 2HDM Revisited,” JHEP 02 (2017) 121, arXiv:1612.04086 [hep-ph].
  * (83) J. M. Cline and P.-A. Lemieux, “Electroweak phase transition in two Higgs doublet models,” Phys. Rev. D 55 (1997) 3873–3881, arXiv:hep-ph/9609240.
  * (84) J. Casas, J. Espinosa, M. Quiros, and A. Riotto, “The Lightest Higgs boson mass in the minimal supersymmetric standard model,” Nucl. Phys. B 436 (1995) 3–29, arXiv:hep-ph/9407389. [Erratum: Nucl.Phys.B 439, 466–468 (1995)].
  * (85) J. Elias-Miro, J. Espinosa, and T. Konstandin, “Taming Infrared Divergences in the Effective Potential,” JHEP 08 (2014) 034, arXiv:1406.2652 [hep-ph].
  * (86) S. P. Martin, “Taming the Goldstone contributions to the effective potential,” Phys. Rev. D 90 (2014) 016013, arXiv:1406.2355 [hep-ph].
  * (87) L. Dolan and R. Jackiw, “Symmetry Behavior at Finite Temperature,” Phys. Rev. D 9 (1974) 3320–3341.
  * (88) M. Carrington, “The Effective potential at finite temperature in the Standard Model,” Phys. Rev. D 45 (1992) 2933–2944.
  * (89) P. B. Arnold and O. Espinosa, “The Effective potential and first order phase transitions: Beyond leading-order,” Phys. Rev. D 47 (1993) 3546, arXiv:hep-ph/9212235. [Erratum: Phys.Rev.D 50, 6662 (1994)].
  * (90) N. Blinov, S. Profumo, and T. Stefaniak, “The Electroweak Phase Transition in the Inert Doublet Model,” JCAP 07 (2015) 028, arXiv:1504.05949 [hep-ph].
  * (91) D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, and G. White, “Theoretical uncertainties for cosmological first-order phase transitions,” arXiv:2009.10080 [hep-ph].
  * (92) C. L. Wainwright, “CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields,” Comput. Phys. Commun. 183 (2012) 2006–2013, arXiv:1109.4189 [hep-ph].
  * (93) G. Dorsch, S. Huber, K. Mimasu, and J. No, “The Higgs Vacuum Uplifted: Revisiting the Electroweak Phase Transition with a Second Higgs Doublet,” JHEP 12 (2017) 086, arXiv:1705.09186 [hep-ph].
  * (94) J. Bernon, L. Bian, and Y. Jiang, “A new insight into the phase transition in the early Universe with two Higgs doublets,” JHEP 05 (2018) 151, arXiv:1712.08430 [hep-ph].
  * (95) A. D. Linde, “Fate of the False Vacuum at Finite Temperature: Theory and Applications,” Phys. Lett. B 100 (1981) 37–40.
  * (96) A. D. Linde, “Decay of the False Vacuum at Finite Temperature,” Nucl. Phys. B 216 (1983) 421. [Erratum: Nucl.Phys.B 223, 544 (1983)].
  * (97) J. Moreno, M. Quiros, and M. Seco, “Bubbles in the supersymmetric standard model,” Nucl. Phys. B 526 (1998) 489–500, arXiv:hep-ph/9801272.
  * (98) L. Leitao, A. Megevand, and A. D. Sanchez, “Gravitational waves from the electroweak phase transition,” JCAP 10 (2012) 024, arXiv:1205.3070 [astro-ph.CO].
  * (99) A. H. Guth and S. Tye, “Phase Transitions and Magnetic Monopole Production in the Very Early Universe,” Phys. Rev. Lett. 44 (1980) 631. [Erratum: Phys.Rev.Lett. 44, 963 (1980)].
  * (100) A. H. Guth and E. J. Weinberg, “Cosmological Consequences of a First Order Phase Transition in the SU(5) Grand Unified Model,” Phys. Rev. D 23 (1981) 876.
  * (101) V. K. S. Shante and S. Kirkpatrick, “An introduction to percolation theory,” Advances in Physics 20 (1971) 325–357.
  * (102) M. D. Rintoul and S. Torquato, “Precise determination of the critical threshold and exponents in a three-dimensional continuum percolation model,” J. Phys. A: Math. Gen. 30 (1997) L585–L592.
  * (103) A. Kobakhidze, C. Lagger, A. Manning, and J. Yue, “Gravitational waves from a supercooled electroweak phase transition and their detection with pulsar timing arrays,” Eur. Phys. J. C 77 (2017) 570, arXiv:1703.06552 [hep-ph].
  * (104) J. Ellis, M. Lewicki, and J. M. No, “On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal,” JCAP 04 (2019) 003, arXiv:1809.08242 [hep-ph].
  * (105) X. Wang, F. P. Huang, and X. Zhang, “Phase transition dynamics and gravitational wave spectra of strong first-order phase transition in supercooled universe,” JCAP 05 (2020) 045, arXiv:2003.08892 [hep-ph].
  * (106) P. J. Steinhardt, “Relativistic Detonation Waves and Bubble Growth in False Vacuum Decay,” Phys. Rev. D 25 (1982) 2074.
  * (107) M. Kamionkowski, A. Kosowsky, and M. S. Turner, “Gravitational radiation from first order phase transitions,” Phys. Rev. D 49 (1994) 2837–2851, arXiv:astro-ph/9310044.
  * (108) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, “Energy Budget of Cosmological First-order Phase Transitions,” JCAP 06 (2010) 028, arXiv:1004.4187 [hep-ph].
  * (109) K. Enqvist, J. Ignatius, K. Kajantie, and K. Rummukainen, “Nucleation and bubble growth in a first order cosmological electroweak phase transition,” Phys. Rev. D 45 (1992) 3415–3428.
  * (110) A. Kosowsky, M. S. Turner, and R. Watkins, “Gravitational radiation from colliding vacuum bubbles,” Phys. Rev. D 45 (1992) 4514–4535.
  * (111) S. J. Huber and T. Konstandin, “Production of gravitational waves in the nMSSM,” JCAP 05 (2008) 017, arXiv:0709.2091 [hep-ph].
  * (112) P. Binetruy, A. Bohe, C. Caprini, and J.-F. Dufaux, “Cosmological Backgrounds of Gravitational Waves and eLISA/NGO: Phase Transitions, Cosmic Strings and Other Sources,” JCAP 06 (2012) 027, arXiv:1201.0983 [gr-qc].
  * (113) C. Caprini et al., “Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions,” JCAP 04 (2016) 001, arXiv:1512.06239 [astro-ph.CO].
  * (114) C. Caprini et al., “Detecting gravitational waves from cosmological phase transitions with LISA: an update,” JCAP 03 (2020) 024, arXiv:1910.13125 [astro-ph.CO].
  * (115) M. B. Hindmarsh, M. Lüben, J. Lumma, and M. Pauly, “Phase transitions in the early universe,” arXiv:2008.09136 [astro-ph.CO].
  * (116) A. Kosowsky, M. S. Turner, and R. Watkins, “Gravitational waves from first order cosmological phase transitions,” Phys. Rev. Lett. 69 (1992) 2026–2029.
  * (117) A. Kosowsky and M. S. Turner, “Gravitational radiation from colliding vacuum bubbles: envelope approximation to many bubble collisions,” Phys. Rev. D 47 (1993) 4372–4391, arXiv:astro-ph/9211004.
  * (118) S. J. Huber and T. Konstandin, “Gravitational Wave Production by Collisions: More Bubbles,” JCAP 09 (2008) 022, arXiv:0806.1828 [hep-ph].
  * (119) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Gravitational waves from the sound of a first order phase transition,” Phys. Rev. Lett. 112 (2014) 041301, arXiv:1304.2433 [hep-ph].
  * (120) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Numerical simulations of acoustically generated gravitational waves at a first order phase transition,” Phys. Rev. D 92 (2015) 123009, arXiv:1504.03291 [astro-ph.CO].
  * (121) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, “Shape of the acoustic gravitational wave power spectrum from a first order phase transition,” Phys. Rev. D 96 (2017) 103520, arXiv:1704.05871 [astro-ph.CO]. [Erratum: Phys.Rev.D 101, 089902 (2020)].
  * (122) J. Ellis, M. Lewicki, J. M. No, and V. Vaskonen, “Gravitational wave energy budget in strongly supercooled phase transitions,” JCAP 06 (2019) 024, arXiv:1903.09642 [hep-ph].
  * (123) J. Ellis, M. Lewicki, and J. M. No, “Gravitational waves from first-order cosmological phase transitions: lifetime of the sound wave source,” JCAP 07 (2020) 050, arXiv:2003.07360 [hep-ph].
  * (124) H.-K. Guo, K. Sinha, D. Vagie, and G. White, “Phase Transitions in an Expanding Universe: Stochastic Gravitational Waves in Standard and Non-Standard Histories,” JCAP 01 (2021) 001, arXiv:2007.08537 [hep-ph].
  * (125) A. Kosowsky, A. Mack, and T. Kahniashvili, “Gravitational radiation from cosmological turbulence,” Phys. Rev. D 66 (2002) 024030, arXiv:astro-ph/0111483.
  * (126) C. Caprini, R. Durrer, and G. Servant, “The stochastic gravitational wave background from turbulence and magnetic fields generated by a first-order phase transition,” JCAP 12 (2009) 024, arXiv:0909.0622 [astro-ph.CO].
  * (127) C. J. Moore, R. H. Cole, and C. P. L. Berry, “Gravitational-wave sensitivity curves,” Class. Quant. Grav. 32 (2015) 015014, arXiv:1408.0740 [gr-qc].
  * (128) E. Thrane and J. D. Romano, “Sensitivity curves for searches for gravitational-wave backgrounds,” Phys. Rev. D 88 (2013) 124032, arXiv:1310.5300 [astro-ph.IM].
  * (129) CMS, LHCb Collaboration, V. Khachatryan et al., “Observation of the rare $B^{0}_{s}\to\mu^{+}\mu^{-}$ decay from the combined analysis of CMS and LHCb data,” Nature 522 (2015) 68–72, arXiv:1411.4413 [hep-ex].
  * (130) LHCb Collaboration, R. Aaij et al., “Measurement of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ branching fraction and effective lifetime and search for $B^{0}\to\mu^{+}\mu^{-}$ decays,” Phys. Rev. Lett. 118 (2017) 191801, arXiv:1703.05747 [hep-ex].
  * (131) HFLAV Collaboration, Y. Amhis et al., “Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of summer 2016,” Eur. Phys. J. C 77 (2017) 895, arXiv:1612.07233 [hep-ex].
  * (132) M. Misiak et al., “Estimate of $\mathcal{B}(\bar{B}\to X_{s}\gamma)$ at $O(\alpha_{s}^{2})$,” Phys. Rev. Lett. 98 (2007) 022002, arXiv:hep-ph/0609232.
  * (133) M. Misiak et al., “Updated NNLO QCD predictions for the weak radiative B-meson decays,” Phys. Rev. Lett. 114 (2015) 221801, arXiv:1503.01789 [hep-ph].
  * (134) K. Griest and D. Seckel, “Three exceptions in the calculation of relic abundances,” Phys. Rev. D 43 (1991) 3191–3203.