#  Particle dark matter density and entropy production in the early universe

Arnab Chaudhuri a111e-mail<EMAIL_ADDRESS>*  |, Maxim Yu. Khlopovb222e-mail<EMAIL_ADDRESS>*  | and Shiladitya Poreya *  |
a Novosibirsk State University
Pirogova ul., 2, 630090 Novosibirsk, Russia
b CNRS, Astroparticule et Cosmologie, Université de Paris, F-75013 Paris,
France;
Institute of Physics, Southern Federal University, 344090 Rostov on Don,
Russia
Center for Cosmopartilce Physics “Cosmion”
National Research Nuclear University “MEPhI” (Moscow Engineering Physics
Institute)
115409 Moscow, Russia  corresponding author

###### Abstract

Dark Matter (DM) density is reduced if entropy production takes place after DM
particles abundance is frozen out in the early universe. We study a
possibility of such reduction due to entropy production in the electroweak
phase transition (EWPT). We compare scenarios of entropy production in the
standard model (SM) and its simplest extension, the two-Higgs doublet model
(2HDM). Assuming the EWPT is of second order in the SM scenario and the first
order in the 2HDM, we calculate the entropy release in these scenarios and the
corresponding dilution of preexisting DM density in the early universe. We
find the effect of dilution in EWPT significant for confrontation with
observations of any form of possible DM (including primordial black holes
(PBHs)), which is frozen out, decoupled, frozen in, or formed before EWPT.

Keywords— electroweak phase transition; dark matter; entropy.

##  1 Introduction

The early universe is associated with the high temperature and density of
particles. In accordance with the Big Bang theory, the temperature during the
very early stages of expansion ranged from few $eV$ to Planck scale $\sim
M_{Pl}$ $\sim 10^{19}$ $GeV$ (more precisely to reheating temperature in
inflationary models). The physical basis of now standard cosmological paradigm
of the inflationary universe with baryosynthesis and Dark Matter (DM)/energy
involves extensions of the particle Standard Model (SM) and depending on the
symmetry breaking pattern of the Beyond the Standard model (BSM) physics may
include several phase transitions in the cosmological scenario. However, any
scenario based on BSM physics should inevitably include phase transitions
predicted by SM, and first of all, related to breaking of the electroweak
symmetry.

The period of the universe before breaking of electroweak symmetry, when
electromagnetic and weak force remained merged into a single electroweak force
is generally known as the electroweak epoch [1, 2]. The electroweak epoch
ended with electroweak phase transition (EWPT). Depending on the parameters of
the BSM model, the nature of the phase transition can be of the first order,
second order, or it can be a smooth crossover. The first order EWPT can be a
source of gravitational waves [3, 4]. EWPT, in general, cannot only be a
potential candidate for baryogenesis (if the three Sakharov’s conditions [5]
are satisfied) but also it can be a source to dilute the preexisting baryon
asymmetry as well as of any form of DM, which is frozen out, decoupled, frozen
in or formed before the EWPT.

With the proviso that thermal equilibrium is sustained, the entropy of the
universe should remain conserved according to the second law of
thermodynamics. In this article we are engrossed in the earliest epochs of the
universe such that Bosonic and fermionic species were massless and with
negligible chemical potential, as mentioned in the textbooks [6, 7], relative
to the temperature of the universe $T(t)$. The above assumption is accurate
just before the EWPT and this with the condition of conservation of energy
density of plasma lead to entropy conservation law:

$\displaystyle s=\frac{\mathcal{P}+\rho}{T}a^{3}=const.$ (1)

where $s$ is the entropy density per comoving volume, $a(t)$ is the scale
factor, $T(t)$ is the temperature of the fluid (or plasma), $\rho$ and $P$ are
the plasma energy density and pressure, respectively.

To attain the state of local thermal equilibrium and the maximum entropy
state, the interaction rate $\Gamma$ among the ultra-relativistic particles in
the hot dense plasma must be greater than Hubble expansion rate, ${\cal H}$.
$\Gamma=\sigma nv$, the product of cross section among particles, their number
density and relative velocity ($v\sim 1$ for photons) and thus it can be
stated that for ${\cal O}(1)$ estimation, $\Gamma\simeq\alpha^{n}T(t)$ ($n=1$
for decay and $2$ for two-body scattering). Under radiation domination, the
${\cal H}$ varies as ${\cal H}\sim T(t)^{2}/M_{Pl}$, $M_{Pl}$ being the
reduced Planck mass. The condition $\Gamma>{\cal H}$ implies
$T<\alpha^{n}M_{Pl}$. This relation was verifiable in the very early universe
because $\alpha\sim 10^{-2}$ and varies logarithmically with $T(t)$ [8].

In thermal equilibrium, the distribution function of any species is defined by
their chemical potential ($\mu_{j}$) and the temperature ($T$) of the plasma.
The temperature of the plasma is similar to all the species. The equilibrium
distribution function of the species are represented by the Bose-Einstein
distribution and Fermi-Dirac distribution for bosons and fermions,
respectively, with $\mu_{j}=0$ and has the following form:

$\displaystyle f(E)=\frac{1}{exp(E/T\pm 1)}.$ (2)

Entropy density is of our special interest because it helps us to estimate the
number density of any particle per comoving volume and more entropy influx
implies dilution of the number density of any particle. The number density of
the frozen out species are conserved in the early universe but there can be
some periods when their dilution is possible due to the influx of entropy into
the plasma. If the energy density of the universe was dominated by that of the
primordial black holes of sufficient masses such that they evaporate in the
corresponding early matter-dominated epoch, the influx of entropy can
sufficiently reduce the pre-existing baryon asymmetry and DM density, if they
were created before this epoch (see [9, 10]). Quantum chromodynamics (QCD)
phase transition at $T\sim$ 150–200 $GeV$ can also lead to such dilution as
well [11].

According to the inflationary cosmology with baryogenesis and DM/energy,
physics beyond the Standard model (BSM) underlying the basis elements of
modern cosmology can provide many scenarios of entropy production (see, e.g.,
[12, 13, 14, 15] for review and references). In this article we have
considered the entropy production within SM and the minimal extension of SM,
namely, two Higgs Doublet Model (2HDM) and discuss the effects of this
mechanism.

According to the $\Lambda-$CDM (Lambda cold dark matter) model of the universe
at the present time, the $\Omega-$parameter of the universe for the cold DM is
$\Omega_{c}\simeq 0.265(7)$ i.e., around $27\%$ of the total matter-energy of
the today’s universe is added from DM sector. DM only interacts by
gravitational force with the SM members and thus it makes the circumstances
really complex for the cosmologists to claim with certainty about the
characteristics and constituents of DM sector. Due to this reason, numerous
models of DM have been proposed for the last a few decades. All those models
can be categorized into three groups - cold dark matter, warm dark matter (see
[16]), and hot dark matter (see p.178 of [7]). There are several proposed
members of these three groups. Among them the promising candidates are
priomordial black holes [17, 18], Axion [19] or Fuzzy DM or Axion like
particles [20], Weakly interacting massive particles [21], 4th generation DM
[22], neutrinos (both cold and hot DM possible) [23, 24], fermionic DM [25],
gravitino [26] etc. Howbeit, the calculation of this script is independent of
the particular DM model as long as their production or freezing out or
freezing in is not connected with the EWPT mechanisms we considered or affect
their parameters and frozen out event has been completed before the beginning
of EWPT. If the DM are produced during or later stage (e.g., [27]), the
dilution will be change due to increase in number density of DM particles. The
scenario of production of DM before EWPT is already well known in literature
(example [20, 28]). But the temperature scale of their decoupling from
radiation depends their mass and cross section of their reaction with the SM
particles.

In this article, entropy density implies the entropy of the relativistic
particles (photons) which dominates the entropy of the universe because there
are significantly more photons than baryons and non-baryons in our universe.
For non-relativistic matter, chemical potential $\mu_{\rm chm.pot.}$ is not
negligible and conservation of entropy can be expressed as continuity equation
of entropy current. If the dark matter is in the form of cold dark matter, its
entropy density behaves as non-relativistic. If it is in the form of hot dark
matter, then entropy density of DM behaves as relativistic and follow Equation
(1). For warm dark matter, which is an intermediate particle between hot and
cold dark matter then it can behave as both. Disquisition about the entropy of
DMs and non-relativistic particles is beyond the scope of this work.

The largest dilution of DM as a result of entropy influx during the course of
universe expansion took place during EWPT. According to electroweak (EW)
theory when $T>T_{c}$, where $T_{c}$ is the critical temperature of EWPT,
$\expectationvalue{\phi}=0$ of the Higgs field, see [29] but for $T<T_{c}$,
$\expectationvalue{\phi}$ increases to the second vacuum expectation value
$\eta$. Such a nonequilibrium state violates the entropy conservation law
releasing entropy into the primeval plasma.

In the SM, EWPT can be of second order or a smooth crossover. In this work, we
have considered EWPT as a second order phase transition. The influx of entropy
can significantly increase if EWPT is of the first order, which happens in
multi Higgs extension of the standard model, like two-Higgs doublet model
(2HDM).

For $T\gg T_{c}$, the universe was in thermal equilibrium, and relativistic
particles dominated the universe. The contribution to the overall energy
density of the universe from those who were already massive (e.g., decoupled
DM) was also insignificant. To this extend, the entropy density per unit
comoving volume follows the conservation law, as mentioned in Equation (1). In
this scenario, the sum of the energy and the pressure density can lead up to

$\rho_{r}+P_{r}\sim g_{*}T^{4}.$ (3)

where $g_{*}$ is the effective number of particle species at or near the EWPT.

It is to be noted that $g_{*}$ is not constant but varies over time. It
depends on the components of the primordial plasma. Equations (1) and (3)
imply:

$T\sim a^{-1}.$ (4)

If at a later stage the thermal equilibrium is broken, $s$ and thus
$g_{*}\left(T\right)a^{3}T^{3}$ increases as entropy should always increase or
remain unchanged. The biggest contribution comes from the heaviest particle
whose mass ($m(T)$) is closest to the temperature $T$, i.e., $m(T)<T$. Here
the product of the effective number of particles to the temperature is
unchanged and the rise in entropy is dominated by the change in the scale
factor, i.e., $a^{3}T^{3}$, suggesting the influx of entropy happening over
$g_{*}a^{3}T^{3}$. And hence the net entropy release is given by

$\displaystyle\frac{\delta
s}{s}=\frac{(a_{c}\,T_{c})^{3}-(a\,T)^{3}}{(a_{c}\,T_{c})^{3}}.$ (5)

In this work, we show that the dilution happens in EWPT in both SM and 2HDM
for all the DM candidates, created, frozen in or frozen out before EWPT. The
paper is organized as follows. In the next two Sections 2 and 3 we calculate
the entropy influx due to SM and 2HDM into the primeval plasma and show the
dilution of pre-EWPT DM density. These are followed by a generic conclusion
that though this effect cannot influence predicted DM density of such popular
DM candidates as WIMPs or QCD axions, which are frozen out or formed after
EWPT, the dilution of the pre-EWPT DM density of sterile neutrinos, superheavy
gravitinos or Primordial Black hole (PBHs) (see, e.g., [12] for review and
references) is significant and should be taken into account in the detailed
analysis of the corresponding cosmological scenarios.

##  2 Entropy Production and Dilution of Frozen Out DM Density within the
Framework of the SM

The Lagrangian of the EWPT in SM framework containing only a single Higgs
field is given as:

$\displaystyle\mathcal{L}=g^{\xi\tau}\partial_{\xi}\phi^{\dagger}\partial_{\tau}\phi-
U_{\phi}(\phi)+\sum_{j}i\left[g^{\xi\tau}\bar{\chi_{j}}\gamma_{\xi}\partial_{\tau}\chi_{j}-U_{j}(\chi_{j})\right]+\mathcal{L}_{int}$
(6)

where $\phi$ is the Higgs field with potential $U_{\phi}(\phi)$. $\chi_{j}$
represents the fermionic fields which gain masses during or immediately after
EWPT, $U_{j}(\chi_{j})$ is the potential of the field $\chi_{j}$, and the
metric is flat FRW background with the signature of the metric
$g_{\xi\tau}=(+,-,-,-)$ with the Greek indices in sub- or superscripts running
from $0$ to $3$. The interaction part of the above Lagrangian density is

$\displaystyle\mathcal{L}_{int}=\phi\sum_{j}g_{j}\bar{\chi_{j}}{\chi}_{j},$
(7)

where $g_{j}$ are the Yukawa coupling constants, and the summation is over all
fields $\chi_{j}$.

In Equation (6), we have used simplistic Lagrangian density instead of the
original form of the Lagrangian density described in [30]. The reason behind
this is that the detailed form of the Gauge bosons and fermionic sector, as we
will see below, is not needed for our calculations, and the reader, thus, can
get the main working principle at ease. However, for completeness, we have
used the exact form of the Lagrangian density in Section 3 for the case of
2HDM.

The homogeneous and isotropic classical field $\phi(t)$ follows the equation
of motion as

$\displaystyle\ddot{\phi}+3\mathcal{H}\dot{\phi}+U^{\prime}_{\phi}(\phi)-\sum_{\chi}g_{j}\bar{\chi_{j}}\chi_{j}=0,$
(8)

where the dot represents derivative w.r.t. time $t$, $\mathcal{H}=\dot{a}/a$
is the Hubble parameter for the cosmological scale factor $a(t)$ and
$U^{\prime}_{\phi}=dU_{\phi}/d\phi$. .

We have considered the following SM potential

$\displaystyle
U_{\phi}(\phi)=\frac{\lambda}{4}(\phi^{2}-\eta^{2})^{2}+\frac{T^{2}\phi^{2}}{2}\,\sum_{j}\,h_{j}\left(\frac{m_{j}(T)}{T}\right),$
(9)

Here, $T$ and $m_{j}(T)$ are the temperature of the relativistic fluid and
mass of the $\chi_{j}$-particle at temperature $T$, respectively.

The last term depends on the temperature and arises from (7) after thermal
averaging and involves the contributions of $\phi$ itself and of all particles
$\chi_{j}$. The summation is made over all particles and the function
$h_{j}(m_{j}/T)\propto g_{j}^{2}$ and is positive always. For $T>m_{j}(T)$ the
function takes a noticeable value and it is negligible for $T<m_{j}(T)$.

According to the experiment, the vacuum expectation value (vev) of the Higgs
field $\phi$ is equal to $v_{\rm sm}\equiv\eta=246$ GeV and the quartic self-
coupling of $\phi$ is $\lambda=0.13$, see [31].

Within the framework of SM the fermionic contribution is important and their
masses at zero temperature is given by $m_{f}=g_{f}\eta$, $g_{f}$ is the
Yukawa coupling between the Higgs and fermionic members. The thermal masses of
the fermions are proportional to the temperature which depends on the value of
$\phi$ at the minima of the potential (9).

$\phi^{2}_{min}(T)=\eta^{2}-(T^{2}/\lambda)\sum_{j}h_{j}\left(\frac{m_{j}(T)}{T}\right)$
(10)

and analogously

$m^{2}_{f}(T)=g^{2}_{f}\phi^{2}_{min}(T)=g^{2}_{f}\left[\eta^{2}-(T^{2}/\lambda)\sum_{j}h_{j}\left(\frac{m_{j}(T)}{T}\right)\right].$
(11)

The dilution of preexisting baryon asymmetry and the dilution of DM density
can be computed from Equation (1) and to do this we need $\mathcal{P}$ and
$\rho$ which can be obtained from the energy–momentum tensor:

$\displaystyle T_{\mu\nu}$ $\displaystyle=$
$\displaystyle\partial_{\mu}\phi\,\partial_{\nu}\phi-
g_{\mu\nu}\left(g^{\alpha\beta}\partial_{\alpha}\phi\,\partial_{\beta}\phi-
U_{\phi}(\phi)\right)$ $\displaystyle+$
$\displaystyle\sum_{j}\left[\partial_{\mu}\bar{\chi_{j}}\,\partial_{\nu}\chi_{j}+\partial_{\nu}\bar{\chi_{j}}\,\partial_{\mu}\chi_{j}-g_{\mu\nu}\left(g^{\alpha\beta}\partial_{\alpha}\bar{\chi_{j}}\,\partial_{\beta}\chi_{j}-U_{j}(\chi_{j})+2\mathcal{L}_{int}\right)\right],$

where $\mathcal{L}_{int}$ is already specified in Equation (7).

The energy and pressure density can be calculated from Equation (2) and is
given by:

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\dot{\phi}^{2}/2+U_{\phi}(\phi)+\sum_{j}\left[\dot{\bar{\chi_{j}}}\dot{\chi}_{j}+\partial_{l}\bar{\chi_{j}}\,\partial_{l}\chi_{j}/a^{2}+U_{j}(\chi_{j})\right]-\mathcal{L}_{int};$
(13) $\displaystyle\mathcal{P}$ $\displaystyle=$
$\displaystyle\dot{\phi}^{2}/2-U_{\phi}(\phi)+\sum_{j}\left[\dot{\bar{\chi_{j}}}\dot{\chi}_{j}-(1/3)\partial_{l}\bar{\chi_{j}}\,\partial_{l}\chi_{j}/a^{2}-U_{j}(\chi_{j})\right]+\mathcal{L}_{int},$
(14)

where $\partial_{i}\chi$ is the derivative with respect to spatial coordinate.

To simplify the calculation, it is assumed that the energy density is the
summation of the energy density of $\phi(t)$ which sits at the minimum of the
potential and relativistic particles, hence $\rho$ becomes:

$\rho\approx
U_{\phi}(\phi_{min})+\frac{\dot{\phi}^{2}}{2}+\frac{\pi^{2}g_{*}}{30}T^{4},$
(15)

along with

${\mathcal{P}}+\rho\approx\dot{\phi}^{2}+\frac{4}{3}\,\frac{\pi^{2}g_{*}}{30}T^{4},$
(16)

At or around the electroweak phase transition, $g_{*}\sim 10^{2}$ is the
effective number of relativistic particles. It is a temperature-dependent
function that decreases as the cosmological cooling process progresses. In the
limit of instant thermalization, Equation (15) holds true.

Under this assumption, the oscillations of $\phi$ around $\phi_{min}$ follows
fast damping, thus we take $\dot{\phi}=\dot{\phi}_{min}$ and neglect higher
order terms as the evolution of $\phi_{min}$ is induced by the universe
expansion which is very slow. A single differential equation of the
temperature (or scale factor) is thus obtained. Finally the equation of energy
conservation can be rewritten as

$\frac{\dot{T}}{T}\left[h_{tot}(m)\eta^{2}T^{2}\left(1-\frac{T^{2}}{T^{2}_{c}}\frac{h_{tot}(m)}{h_{tot}(0)}\right)+\frac{4\pi^{2}g_{*}}{30}T^{4}\right]=-4\mathcal{H}\frac{\pi^{2}g_{*}}{30}T^{4}.$
(17)

where $\sum h_{j}(m_{j}/T)\equiv h_{tot}(m)$ and $\sum h_{j}(0)\equiv
h_{tot}(0)$.

Following the above assumptions, the relative influx of entropy which leads to
the corresponding dilution of DM is given by:

$\frac{\delta
s}{s}=\sum_{f}x_{j,min}^{6\nu_{f}}\exp\left[3\nu_{f}\left(\frac{1}{x_{f,min}^{2}}-1\right)\right]-1,$
(18)

where $x=T/T_{c}$, $\nu_{f}=\frac{\kappa\eta^{2}}{2T_{c}^{2}}$ and
$\kappa=\frac{30h_{tot}(m)}{4\pi^{2}g_{*}}$. As a result of entropy inflow
into the primordial plasma, the total dilution of frozen-out DM particles is
shown in the Figure 1 .

Figure 1: The black line shows the entropy released into the plasma, and the
blue line shows the corresponding dilution of DM density. The net dilution is
$\sim$13%.

Within the framework of the SM, the EWPT is assumed to be of second order.
Even though the dilution is quite small, it is noticeable. The $t$-quark (top-
quark), the heaviest SM particle contributes mostly to the influx of the
entropy. In the extended versions of SM with multi Higgs, the influx and in
return the dilution is quite larger as seen in Section 3.

##  3 Entropy Production and Dilution of Frozen Out DM Density within the
Framework of 2HDM

In this section, we explore entropy generation and the dilution of frozen out
DM density within the context of the two-Higgs doublet model (2HDM). It is the
most popular extension of the scalar sector of the SM. An additional scalar
doublet is necessary for the supersymmetric models. Non-supersymmetric of this
model includes a successful electroweak baryogenesis [32]. There are various
sectors of 2HDM like Type-I 2HDM, Type-II 2HDM, etc. A detailed review of 2HDM
describing all the sectors can be found in [33]. In this work, we confine
ourselves to real Type-I 2HDM where the charged fermions only couple to second
doublet.

The Lagrangian for the 2HDM scenario is given in the following form:

$\mathcal{L}=\mathcal{L}_{f}+\mathcal{L}_{\rm Yuk}+\mathcal{L}_{\rm
gauge,kin}+\mathcal{L}_{\rm Higgs}.$ (19)

and each term on the R.H.S. is defined as follows: The first term is the
kinetic term originating due to the fermionic fields:

$\displaystyle\mathcal{L}_{f}$ $\displaystyle=$
$\displaystyle\sum_{j}i(\bar{\Psi}^{(j)}_{L}\not{D}\Psi^{(j)}_{L}+\bar{\Psi}^{(j)}_{R}\not{D}\Psi^{(j)}_{R})$
(21) $\displaystyle=$ $\displaystyle
i\bar{\Psi}_{L}\gamma^{\xi}(\partial_{\mu}+igW_{\mu}+ig^{\prime}Y_{L}B_{\mu})\Psi_{L}$
$\displaystyle+i\bar{\Psi}_{R}\gamma^{\mu}(\partial_{\mu}+igW_{\mu}+ig^{\prime}Y_{R}B_{\mu})\Psi_{R}.$

Here, $L$ and $R$ are the left and right chiral field of that fermion,
$\not{D}$ is the the gauge covariant derivative in Feynman notation (Feynman
notation is defined as $\not{A}=\gamma^{\mu}A_{\mu}$, where $\gamma$ are the
gamma matrices) and $j$ runs over all fermionic species (the field
$\Psi_{j}$). In addition, $g$ is the coupling constant. The partial derivative
follows from the previous section, e.g., $\partial_{0}=d/dt$ and $\mu$ runs
from $0$ to $3$, [34, 35].

The second term of Equation (19), the Yukawa interaction term is given by

$\displaystyle\mathcal{L}_{\rm
Yuk}=-\left[y_{e}\bar{e_{R}}\Phi_{a}^{\dagger}L_{L}+y_{e}^{*}\bar{L_{L}}\Phi_{a}^{\dagger}e_{R}+\cdots\right],$
(22)

where $y_{e}$ is a dimensionless constant which is complex in nature,
$\Phi_{a}$ ($a=1,2$) is a $SU(2)_{L}$ doublet, and since the Lagrangian needs
to be gauge invariant, it is coupled with another $SU(2)_{L}$ fermion $L_{L}$.
$e_{R}$ follows from above, is the right chiral electron field, and so are the
other fermions, e.g., quarks, neutrinos, etc.

The third term $\mathcal{L}_{\rm gauge,\,kin}$ represents $U(1)$ invariant
kinetic term of four gauge bosons ($W^{i},\,i=1,2,3$, and $B$). It can be
written in the following manner:

$\displaystyle\mathcal{L}_{\rm
gauge,\,kin}=-\frac{1}{4}G^{i}_{\mu\nu}{G^{i}}^{\mu\nu}-\frac{1}{4}F^{B}_{\mu\nu}{F^{B}}^{\mu\nu},$
(23)

where
$G^{i}_{\mu\nu}=\partial_{\mu}W^{i}_{\nu}-\partial_{\nu}W^{i}_{\mu}-g\epsilon^{ijk}W_{\mu}^{j}W_{\nu}^{k}$
and $F^{B}_{\mu\nu}=\partial_{\mu}B_{\nu}-\partial_{\nu}B_{\mu}$ with
$\epsilon_{ijk}$ is the Levi-Civita symbol.

The fourth term is the Lagrangian for the Higgs boson doublets:

$\displaystyle\mathcal{L}_{\rm Higgs}$ $\displaystyle=$
$\displaystyle(D^{\mu}\Phi_{1})^{\dagger}(D_{\mu}\Phi_{2})+(D^{\mu}\Phi_{1})^{\dagger}(D_{\mu}\Phi_{2})-V_{\rm
tot}(\Phi_{1},\Phi_{2})$ $\displaystyle=$
$\displaystyle\\{(\partial_{\mu}+igT^{i}W^{i}_{\mu}+ig^{\prime}YB_{\mu})\Phi_{1}\\}^{\dagger}\\{(\partial_{\mu}+igT^{i}W^{i}_{\mu}+ig^{\prime}YB_{\mu})\Phi_{1}\\}$
$\displaystyle+\\{(\partial_{\mu}+igT^{i}W^{i}_{\mu}+ig^{\prime}YB_{\mu})\Phi_{2}\\}^{\dagger}\\{(\partial_{\mu}+igT^{i}W^{i}_{\mu}+ig^{\prime}YB_{\mu})\Phi_{2}\\}-V_{\rm
tot}(\Phi_{1},\Phi_{2},T).$

For simplification, we define
$\mathcal{W}_{\mu}=gT^{i}W^{i}_{\mu}+g^{\prime}YB_{\mu}$, we’ve from Equation
(3):

$\displaystyle\mathcal{L}_{\rm Higgs,\,kin}$ $\displaystyle=$
$\displaystyle(\partial^{\mu}{\Phi_{a}}^{\dagger})(\partial_{\mu}\Phi_{a})-i({\mathcal{W}^{\mu}}\Phi_{a})^{\dagger}(\partial_{\mu}\Phi_{a})+$
$\displaystyle
i(\partial^{\mu}{\Phi_{a}}^{\dagger})\mathcal{W}_{\mu}\Phi_{a}+(\mathcal{W}^{\mu}\Phi_{a})^{\dagger}\mathcal{W}_{\mu}\Phi_{a}.$

The type-I 2HDM potential and correction terms are as follows:

$V(\Phi_{1},\Phi_{2},T)=V_{tree}(\Phi_{1},\Phi_{2})+V_{CW}(\Phi_{1},\Phi_{2})+V_{T}(T){+V_{\rm
daisy}(T)},$ (26)

The tree level potential is given by:

$\displaystyle V_{\rm tree}(\Phi_{1},\Phi_{2})=$ $\displaystyle
m_{11}^{2}\Phi_{1}^{\dagger}\Phi_{1}+m_{22}^{2}\Phi_{2}^{\dagger}\Phi_{2}-\left[m_{12}^{2}\Phi_{1}^{\dagger}\Phi_{2}+m_{12}^{*}\Phi_{2}^{\dagger}\Phi_{1}\right]+\frac{1}{2}\lambda_{1}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)^{2}$
$\displaystyle+\frac{1}{2}\lambda_{2}\left(\Phi_{2}^{\dagger}\Phi_{2}\right)^{2}+\lambda_{3}\left(\Phi_{1}^{\dagger}\Phi_{1}\right)\left(\Phi_{2}^{\dagger}\Phi_{2}\right)+\lambda_{4}\left(\Phi_{1}^{\dagger}\Phi_{2}\right)\left(\Phi_{2}^{\dagger}\Phi_{1}\right)$
(27)
$\displaystyle+\left[\frac{1}{2}\lambda_{5}\left(\Phi_{1}^{\dagger}\Phi_{2}\right)^{2}+\frac{1}{2}\lambda_{5}^{*}\left(\Phi_{2}^{\dagger}\Phi_{1}\right)^{2}\right].$

Here, $m_{12}$ is a function of the mixing angle $\beta$ which follows
$tan\beta=v2/v1$, where $v_{1}$ and $v_{2}$ are the expectation values of
$\Phi_{1}$ and $\Phi_{2}$, respectively. Despite the fact that $\tan\beta$ is
a physical parameter, it is still feasible to redefine the two doublets and
arrive at a basis where the full vev is totally contained within one of the
two doublets, which is named as Higgs basis, see [36]. The quartic coupling
constants are real for our scenario and hence $\lambda_{5}^{*}=\lambda_{5}$.
The explicit expressions for the couplings are given by:

$\displaystyle\lambda_{1}$ $\displaystyle=\frac{1}{v_{\rm
sm}^{2}\cos^{2}\beta}\left(-\mu^{2}\tan\beta+m_{h}^{2}\sin^{2}\alpha+m_{H}^{2}\cos^{2}\alpha\right)$
(28) $\displaystyle\lambda_{2}$ $\displaystyle=\frac{1}{v_{\rm
sm}^{2}\sin^{2}\beta}\left(-\frac{\mu^{2}}{\tan\beta}+m_{h}^{2}\cos^{2}\alpha+m_{H}^{2}\sin^{2}\alpha\right)$
(29) $\displaystyle\lambda_{3}$ $\displaystyle=\frac{1}{v_{\rm
sm}^{2}}\left(-\frac{2\mu^{2}}{\sin
2\beta}+2m_{H_{\pm}}^{2}+\left(m_{H}^{2}-m_{h}^{2}\right)\frac{\sin
2\alpha}{\sin 2\beta}\right)$ (30) $\displaystyle\lambda_{4}$
$\displaystyle=\frac{1}{v_{\rm sm}^{2}}\left(\frac{2\mu^{2}}{\sin
2\beta}+m_{A}^{2}-2m_{H_{\pm}}^{2}\right)$ (31) $\displaystyle\lambda_{5}$
$\displaystyle=\frac{1}{v_{\rm sm}^{2}}\left(\frac{2\mu^{2}}{\sin
2\beta}-m_{A}^{2}\right)$ (32)

All the masses (i.e., $m_{h},m_{H_{\pm}},m_{H},m_{A}$ and other particles) are
dependent on the second minimum of the potential. For detailed calculation of
the masses, see [8]. Here, $\alpha$.

The other terms in Equation (26) are given below:

$\displaystyle V_{\rm
CW}\left(v_{1}+v_{2}\right)=\sum_{j}\frac{n_{j}}{64\pi^{2}}(-1)^{2s_{j}}m_{j}^{4}\left(v_{1},v_{2}\right)\left[\log\left(\frac{m_{j}^{2}\left(v_{1},v_{2}\right)}{\mu^{2}}\right)-c_{j}\right]$
(33) $\displaystyle V_{T}=\frac{T^{4}}{2\pi^{2}}\left(\sum_{j={\rm
bosons}}n_{j}J_{B}\left[\frac{m_{j}^{2}(v_{1},v_{2})}{T^{2}}\right]+\sum_{j={\rm
fermions}}n_{j}J_{F}\left[\frac{m_{j}^{2}(v_{1},v_{2})}{T^{2}}\right]\right)$
(34)

where $s_{j}$ is the spin of the $j-$th particle, and $\mu$ is the
renormalisation scale, whose value is takes as $246$ GeV. Detailed description
of the parameters are done in [8].

$J_{B}$ and $J_{F}$ are approximated in series expansion at high temperature
in the Landau gauge up to the leading orders:

$\displaystyle T^{4}J_{B}\left[\frac{m^{2}}{T}\right]$
$\displaystyle=-\frac{\pi^{4}T^{4}}{45}+\frac{\pi^{2}}{12}T^{2}m^{2}-\frac{\pi}{6}T(m^{2})^{3/2}-\frac{1}{32}m^{4}\ln\frac{m^{2}}{a_{b}T^{2}}+\cdots,$
(35) $\displaystyle T^{4}J_{F}\left[\frac{m^{2}}{T}\right]$
$\displaystyle=\frac{7\pi^{4}T^{4}}{360}-\frac{\pi^{2}}{24}T^{2}m^{2}-\frac{1}{32}m^{4}\ln\frac{m^{2}}{a_{f}T^{2}}+\cdots,$
(36)

where $a_{b}=16a_{f}=16\pi^{2}\exp(3/2-2\gamma_{E})$ with $\gamma_{E}$ being
the Euler–Mascheroni constant.

The daisy correction term is given by:

$\displaystyle
V_{\text{daisy}}(T)=-\frac{T}{12\pi}\left[\sum_{j=1}^{n_{\text{Higgs}}}\left(\bar{(m}^{2}_{j})^{3/2}-(m_{j}^{2})^{3/2}\right)+\sum_{j=1}^{n_{\text{gauge}}}\left(\bar{(m}^{2}_{j})^{3/2}-(m_{j}^{2})^{3/2}\right)\right]$
(37)

As the universe cools down, a second local minimum appears. As the temperature
of the plasma drops to $T_{c}$ the second local minimum takes the form
$(\left<\Phi_{1}\right>=v_{1},\left<\Phi_{2}\right>=v_{2})$ and become
degenerate with the global one at
$(\left<\Phi_{1}\right>=\left<\Phi_{2}\right>=0)$. Thus the critical
temperature is obtained by the following expression:

$V_{\rm tot}\left(\Phi_{1}=0,\Phi_{2}=0,T_{c}\right)=V_{\rm
tot}\left(\Phi_{1}=v_{1},\Phi_{2}=v_{2},T_{c}\right).$ (38)

The energy density of the homogeneous classical field $\Phi$:

$\displaystyle\rho=$
$\displaystyle\partial^{0}\Phi_{a}^{\dagger}\partial^{0}\Phi_{a}-(\mathcal{W}^{0}\Phi_{a})^{\dagger}\mathcal{W}^{0}\Phi_{a}-(\mathcal{W}^{j}\Phi_{a})^{\dagger}\mathcal{W}_{j}\Phi_{a}$
$\displaystyle+\left[V_{\rm tot}(\Phi_{1},\Phi_{2},T)-\mathcal{L}_{\rm
gauge,kin}-\mathcal{L}_{f}-\mathcal{L}_{\rm Yuk}\right]$ (39)

where, the spatial derivatives of Higgs fields vanish due to the condition of
homogeneity and isotropy. Similarly,

$\displaystyle\rho+P=2\partial^{0}\Phi_{a}\partial^{0}\Phi_{a}^{\dagger}-i({\mathcal{W}^{0}}\Phi_{a})^{\dagger}(\partial_{0}\Phi_{a})+i(\partial^{0}{\Phi_{a}}^{\dagger})\mathcal{W}_{0}\Phi_{a}.$
(40)

Following the assumption from Section 2 (before Equation (17)), the
oscillations follow fast damping and hence higher order terms are neglected.
And hence we simplify Equation (40) to the following expression:

$\displaystyle\rho=\dot{\Phi}_{a,{\rm min}}^{2}+V_{\rm
tot}(\Phi_{1},\Phi_{2},T)+\frac{g_{*}\pi^{2}}{30}T^{4}.$ (41)

The endmost term on the right hand side of Equation (41) is the energy density
of the relativistic particles that have not obtained mass until EWPT. The
Yukawa interaction between fermions and Higgs bosons, as well as the energy
density of fermions, gauge bosons, and the Higgs-gauge boson interaction
contribute to this term.

Entropy is conserved, and the density of the frozen out DM species remains
constant as long as the primeval plasma, which is assumed to be an ideal
fluid, remains constant. But below the critical temperature $T_{c}$, the
universe goes into the state of non-equilibrium. This situation is rather
useful since 2HDM can successfully explain electroweak baryogenesis and
deviation from thermal equilibrium is one of the primary conditions for
baryogenesis. Due to this non-equilibrium state, the entropy conservation law
breaks and the frozen out DM density starts diluting due to the influx of
entropy into the plasma.

In order to determine the dilution of the frozen out species it is important
to solve the evolution equation for the conservation of the energy density,

$\dot{\rho}=-3\mathcal{H}(\rho+P).$ (42)

Unlike the previous section, it is not easy to proceed forward with analytical
calculations. Hence numerical calculations are done from henceforth. For this
purpose we use the public code: Beyond the Standard Model Phase Transitions
(BSMPT) [38, 37], which is which is written in C++ programming language and it
can be employed for calculating the strength of the electroweak phase
transition in extended Higgs sectors. This tool takes care the corrected terms
of the tree-level potential, including loop-correction and daisy re-summation
of the bosonic masses.

The analytical calculations for ${\delta s}/{s}$ were possible in the previous
section because primarily because of the single Higgs field and the simplified
Lagrangian that was used. However, in the case of 2HDM, there are two Higgs
fields, and the situation gets more complex as we are taking the exact
Lagrangian, and hence numerical calculations are needed.

As we have seen in the previous section, near $T_{c}$ the Higgs potential gets
a degenerate minimum and as the temperature falls off to the mass of any
relativistic component, that component acquires mass, and becomes non-
relativistic. Assuming this process is instantaneous and there is a change in
$g_{*}$, entropy density $s$ increases. With the evolution of scale factor
(and with time actually), $T$ changes, and with the temperature the value of
$\phi_{\rm min}$ where the potential gets second minima. The variation of
$\phi_{\rm min}$ with $T$ can be obtained from BSMPT, as we will see below. To
find the variation of $T$ with $a$, we need Equation (42).

The program is used to calculate the temperature dependent vev ($v$) and
$T_{c}$. Here, we considered only the real sector of Type-I CP-conserving 2HDM
potential and numerically computed the vevs, $T_{c}$ and the effective value
of the potential $V_{eff}(T)$ for all the benchmark points in Table 1 with the
constraining condition $vev/T_{c}>0.01$. The differential Equation (42) was
computed numerically for all the benchmark points in Table 1, and the dilution
factor, as well entropy production for some benchmark points, are shown in
Figure 2.

Table 1: 2HDM Benchmark points for entropy production.

|  $m_{h}$ [GeV] |  $m_{H}$ [GeV] |  $m_{H^{\pm}}$ [GeV] |  $m_{A}$ [GeV] | $\tan\beta$ | $\cos\left(\beta-\alpha\right)$ | $m_{12}^{2}~{}\text{GeV}^{2}$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{3}$ | $\lambda_{4}$ | $\lambda_{5}$ | $T_{c}$ | $vev/T_{c}$ | $\delta s/s[\%]$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
BM1 | $125$ | $500$ | $500$ | $500$ | $2$ | $0$ | $10^{5}$ | $0.258$ | $0.258$ | $0.258$ | $0$ | $0$ | $161.36$ | 1.4 | $57$
BM2 |  "‡ | " | " | " | " | $0.06$ | " | $1.14$ | $0.037$ | $0.63$ | $0$ | $0$ | $167.95$ | 1.6 | $59$
BM3 | " | " | " | " | $10$ | $0$ | 24,752.5 | $0.258$ | $0.258$ | $0.258$ | $0$ | $0$ | $161.02$ | 1.4 | $56$
BM4 | " | " | " | " | $10$ | $0.1$ | 24,752.5 | $4.13$ | $0.22$ | $4.15$ | $0$ | $0$ | $255.71$ | 1.9 | $73$
BM5 | " | " | " | $485$ | $2$ | $0.0$ | $10^{5}$ | $0.26$ | $0.26$ | $0.26$ | $-0.244$ | $0.244$ | $161.53$ | 1.4 | $57$
BM6 | " | " | " | $485$ | $2$ | $0.07$ | $10^{5}$ | $1.28$ | $0.002$ | $0.7$ | $-0.244$ | $0.244$ | $169.81$ | 1.7 | $60$
BM7 | " | " | " | $477$ | $2$ | $0.07$ | $10^{5}$ | $1.28$ | $0.002$ | $0.7$ | $-0.37$ | $0.37$ | $169.53$ | 1.7 | $60$
BM8 | " | " | " | $485$ | $10$ | $0.0$ | 24,752.5 | $0.258$ | $0.258$ | $0.258$ | $-0.244$ | $0.244$ | $160.76$ | 1.3 | $56$
BM9 | " | " | " | $350$ | $10$ | $0.1$ | 24,752.5 | $4.13$ | $0.22$ | $4.15$ | $-2.1$ | $2.1$ | $209.87$ | 1.8 | $68$
BM10 | " | " | $485$ | $500$ | $2$ | $0.00$ | $10^{5}$ | $0.258$ | $0.258$ | $-0.23$ | $0.49$ | $0$ | $153.27$ | 1.25 | $53$
BM11 | " | " | $485$ | $500$ | $2$ | $0.07$ | $10^{5}$ | $1.28$ | $0.002$ | $0.21$ | $0.49$ | $0$ | $169.28$ | 1.7 | $60$
BM12 | " | " | $485$ | $500$ | $10$ | $0.0$ | 24,752.5 | $0.258$ | $0.258$ | $-0.23$ | $0.49$ | $0$ | $160.51$ | 1.3 | $56$
BM13 | " | " | $485$ | $500$ | $10$ | $0.1$ | 24,752.5 | $4.13$ | $0.22$ | $3.66$ | $0.49$ | $0$ | $241.75$ | 1.88 | $70$
BM14 | " | " | $485$ | $485$ | $2$ | $0$ | $10^{5}$ | $0.258$ | $0.258$ | $-0.23$ | $0.24$ | $0.24$ | ${159.76}$ | 1.3 | $59$
BM15 | " | " | $485$ | $485$ | $2$ | $0.07$ | $10^{5}$ | $1.28$ | $0.002$ | $0.21$ | $0.244$ | $0.244$ | ${168.61}$ | 1.7 | $59$
BM16 | " | " | $485$ | $485$ | $10$ | $0$ | 24,752.5 | $0.258$ | $0.258$ | $-0.23$ | $0.244$ | $0.244$ | ${160.19}$ | 1.3 | $56$
BM17 | " | $485$ | $485$ | $485$ | $2$ | $0.0$ | 94,090 | $0.258$ | $0.258$ | $0.258$ | $0$ | $0$ | ${161.31}$ | 1.4 | $57$
BM18 | " | $485$ | $485$ | $485$ | $2$ | $0.07$ | 94,090 | $1.22$ | $0.02$ | $0.67$ | $0$ | $0$ | ${169.7}$ | 1.7 | $60$
BM19 | " | $485$ | $485$ | $485$ | $10$ | $0$ | 23,289.6 | $0.258$ | $0.258$ | $0.258$ | $0$ | $0$ | ${160.96}$ | 1.3 | $57$
BM20 | " | $485$ | $485$ | $485$ | $10$ | $0.1$ | 23,289.6 | $3.9$ | $0.22$ | $3.9$ | $0$ | $0$ | ${230.18}$ | 1.86 | $70$
BM21 | " | $485$ | $485$ | $500$ | $2$ | $0$ | 94,090 | $0.258$ | $0.258$ | $0.258$ | $0$ | $0$ | ${161.31}$ | 1.4 | $57$
BM22 | " | $90$ | $200$ | $300$ | $2$ | $0$ | $3240$ | $0.258$ | $0.258$ | $1.31$ | $0.3$ | $-1.35$ | ${150.76}$ | 1.2 | $51$
BM23 | " | $90$ | $200$ | $300$ | $10$ | $0$ | $801.98$ | $0.258$ | $0.258$ | $1.31$ | $0.3$ | $-1.35$ | ${135.38}$ | 1.06 | $37$
BM24 | " | $90$ | $200$ | $300$ | $10$ | $0.2$ | $801.98$ | $0.263$ | $0.258$ | $1.06$ | $0.3$ | $-1.35$ | ${141.06}$ | 1.1 | $42$

$\ddagger$ It means the present value follows the previous value in the same
column.

Figure 2: For $6$ benchmark points the dilution factor,i.e., the amount of
dilution of the frozen out DM particle density (blue lines) and the entropy
influx factor (other coloured lines) are shown. Each plot corresponds to a
particular benchamrk point which is mentioned within the plot itself.

As mentioned above, and now evident from the Table 1 and Figure 2 the dilution
is more intense in the scenario when multiple higgs fields are involved.

##  4 Conclusions

As shown in this paper, the preexisting baryon asymmetry, and the frozen out
DM density of the universe can be diluted significantly due to electroweak
phase transition. In particular, such dilution can take place for DM particles
with the mass $m\gg 5$ TeV, particles interacting with plasma with cross-
section $\sigma\leq(T_{c}M_{Pl})^{-1}$ and primordial black holes (PBH) with
mass (see, e.g., [18] for review and references)

$M<M_{Pl}\left(\frac{M_{Pl}}{T_{c}}\right)^{2}\approx 2\times
10^{28}\gamma\approx 10^{-5}M_{\odot}.$

($5.6\times 10^{28}$ $gram$ or $3\times 10^{-5}M_{\odot}$ for SM and $3\times
10^{28}$ $gram$ or $1.5\times 10^{-5}$ $M_{\odot}$ for 2HDM
($T_{c}=245~{}$$GeV$)). These types of DM are frozen out, decoupled, frozen
in, or formed before EWPT. Their dilution is achieved with the influx of
entropy into the primordial plasma, which was in the state of thermal
equilibrium before the onset of EWPT. We have considered two scenarios. First,
in the SM scenario, we have assumed EWPT is of second order and we see the
dilution factor is $\sim$13%. Second, we have considered the simplest
extension of SM, namely 2HDM. From Table 1, we see a significant difference
between the two scenarios. Figures 1 and 2 also show major differences. This
is because more entropy influx happens for first order phase transition which
is the case for 2HDM.

One should note that such popular DM candidates, as WIMPs, QCD axion or Axion
Like Particles (ALP) may avoid the considered effect of entropy inflow at the
EWPT. Indeed, WIMPs with mass, less than $fT_{c}$ (where the numerical factor
$f\sim 20$) freeze out after EWPT, while QCD axion condensate is formed at QCD
phase transition and the condensate of ALP with masses of order $10^{-20}$ eV
even much later. The existence of these popular DM candidates is, however,
still not experimentally verified and turning to such DM candidates as sterile
neutrinos, which decouple before EWPT or superheavy gravitino, created just
after inflation (see, e.g., [12] for review and references), we find that
effect of entropy inflow at EWPT would be significant for detailed
descriptions of the corresponding cosmological scenarios. We may conclude that
the entropy inflow effect at EWPT may play important role in all the scenarios
with superheavy and super-weakly interacting DM candidates.

DM candidates, which can be diluted in EWPT, are predicted in much wider
extensions of SM, than we have considered here. It would imply with necessity
self-consistent treatment of effect of dilution with the account for all the
specific model-dependent consequences of physics of early Universe, which are
predicted in each corresponding scenario. However, the results of the present
paper have demonstrated that even in the case of minimal extension of SM
dilution in the EWPT of primordial DM or PBH density is significant and
deserves special study for the confrontation of DM models with observations.

## Acknowledgement

The work of S.P. and A.C. is funded by RSF Grant 19-42-02004. The research by
M.K. was supported by the Ministry of Science and Higher Education of the
Russian Federation under Project “Fundamental problems of cosmic rays and dark
matter”, No. 0723-2020-0040.

## References

  * [1] Ryden, B. Introduction to Cosmology; Addison-Wesley: San Francisco, CA, USA, 2003; ISBN 0-8053-8912-1. https://doi.org/10.1017/9781316651087
  * [2] Jonathan, A. Quarks, Leptons, and the Big Bang; CRC Press: Boca Raton, FL, USA, 2016; ISBN 978-0-7503-0806-9. https://doi.org/10.1201/9781315381367
  * [3] Caprini, C.; Chala, M.; Dorsch, G.C.; Hindmarsh, M.; Huber, S.J.; Konstandin, T.; Kozaczuk, J.; Nardini, G.; No, J.M.; Rummukainen, K.; et al. Detecting gravitational waves from cosmological phase transitions with LISA: An update. J. Cosmol. Astropart. Phys. 2020, 3, 024, https://doi.org/10.1088/1475-7516/2020/03/024.
  * [4] Ghiglieri, J.; Jackson, G.; Laine, M.; Zhu, Y. Gravitational wave background from Standard Model physics: Complete leading order. J. High Energy Phys. 2020, 7, 092. https://doi.org/10.1007/JHEP07(2020)092.
  * [5] Sakharov, A.D. Violation of CP Invariance, C Asymmetry, and Baryon Asymmetry of the Universe. J. Exp. Theor. Phys. Lett. 1967, 5, 24–27. https://doi.org/10.1070/PU1991v034n05ABEH002497
  * [6] Gorbunov, D.S.; Rubakov, V.A. Introduction to the Theory of the Early Universe: Cosmological Perturbations and Inflationary Theory; World Scientific: Singapore, 2011. https://doi.org/10.1142/7874.
  * [7] Bambi, C.; Dolgov, A.D. Introduction to Particle Cosmology; Springer: Berlin/Heidelberg, Germany, 2015. https://doi.org/10.1007/978-3-662-48078-6.
  * [8] Chaudhuri, A.; Khlopov, M.Y. Entropy production due to electroweak phase transition in the framework of two Higgs doublet model. Physics 2021, 3, 275. https://doi.org/10.3390/physics3020020.
  * [9] Dolgov, A.D.; Naselsky, P.D.; Novikov, I.D. Gravitational waves, baryogenesis, and dark matter from primordial black holes. arXiv 2000, arXiv:astro-ph/0009407.
  * [10] Chaudhuri, A.; Dolgov, A. PBH evaporation, baryon asymmetry, and dark matter. arXiv 2020, arXiv:2001.11219.
  * [11] Boeckel, T.; Schettler, S.; Schaffner-Bielich, J. The Cosmological QCD Phase Transition Revisited. Prog. Part. Nucl. Phys. 2011, 66, 266–270. https://doi.org/10.1016/j.ppnp.2011.01.017
  * [12] Khlopov, M.Y. What comes after the Standard model? Prog. Part. Nucl. Phys. 2021, 116 103824. https://doi.org/10.1016/j.ppnp.2020.103824
  * [13] Chaudhuri, A.; Khlopov, M.Y. Balancing Asymmetric Dark Matter with Baryon Asymmetry and Dilution of Frozen Dark Matter by Sphaleron Transition. Universe 2021, 7, 275. https://doi.org/10.3390/universe7080275.
  * [14] Chaudhuri, A.; Khlopov, M.Y.; Porey, S. Effects of 2HDM in Electroweak Phase Transition. Galaxies 2021, 9, 45. https://doi.org/10.3390/galaxies9020045
  * [15] Chaudhuri, A.; Khlopov, M.Y.; Porey, S. Entopy release in Electroweak Phase Transition in 2HDM. arXiv 2021, arXiv:2111.00139.
  * [16] Viel, M.; Becker, G.D.; Bolton, J.S.; Haehnelt, M.G. Warm dark matter as a solution to the small scale crisis: New constraints from high redshift Lyman-$\alpha$ forest data. Phys. Rev. D 2013, 88, 043502. https://doi.org/10.1103/PhysRevD.88.043502.
  * [17] Carr, B.; Kühnel, F. Primordial black holes as dark matter: Recent developments. Ann. Rev. Nucl. Part. Sci. 2020, 70, 355–394. https://doi.org/10.1146/annurev-nucl-050520-125911
  * [18] Khlopov, M.Y. Primordial black holes. Res. Astron. Astrophys. 2010, 10, 495–528. https://doi.org/10.1088/1674-4527/10/6/001.
  * [19] Chadha-Day, F.; Ellis, J.; Marsh, D.J. Axion Dark Matter: What is it and Why Now? arXiv 2021, arXiv:2105.01406.
  * [20] Visinelli, L. Light axion-like dark matter must be present during inflation. Phys. Rev. D 2017, 96, 023013. https://doi.org/10.1103/PhysRevD.96.023013.
  * [21] Schumann, M. Direct detection of WIMP dark matter: Concepts and status. J. Phys. G Nucl. Part. Phys. 2019, 46, 103003. https://doi.org/10.1088/1361-6471/ab2ea5.
  * [22] Khlopov, M.Y. Composite dark matter from 4th generation. Pisma Zh. Eksp. Teor. Fiz. 2006, 83, 3–6. https://doi.org/10.1134/ S0021364006010012.
  * [23] Ibarra, A. Neutrinos and dark matter. AIP Conf. Proc. 2015, 1666, 140004. https://doi.org/10.1063/1.4915588.
  * [24] Pandey, S.; Karmakar S.; Rakshit S. Interactions of astrophysical neutrinos with dark matter: A model building perspective.  J. High Energy Phys. 2019, 1, 95. https://doi.org/10.1007/JHEP11(2021)215.
  * [25] Lee, K.Y.; Kim, Y.G.; Shin, S. Singlet fermionic dark matter. J. High Energy Phys. 2008, 5, 100. https://doi.org/10.1088/1126-6708/2008/05/100.
  * [26] Steffen, F.D. Gravitino dark matter and cosmological constraints. J. Cosmol. Astropart. Phys. 2006, 9, 1. https://doi.org/10.1088/1475-7516/2006/09/001.
  * [27] Falkowski, A.; No, J.M. Non-thermal dark matter production from the electroweak phase transition: Multi-TeV WIMPs and “baby-zillas”. J. High Energy Phys. 2013, 2, 34. https://doi.org/10.1007/JHEP02(2013)034.
  * [28] Peebles, P.J.E. Fluid dark matter. Astrophys. J. Lett. 2000, 534, L127. https://doi.org/10.1086/312677.
  * [29] Bochkarev, A.I.; Kuzmin, S.V.; Shaposhnikov, M.E. Electroweak baryogenesis and the Higgs boson mass problem. Phys. Lett. B 1990, 244, 275–278. https://doi.org/10.1016/0370-2693(90)90069-I.
  * [30] Logan, H.E. TASI 2013 lectures on Higgs physics within and beyond the Standard Model. arXiv 2013, arXiv:1406.1786.
  * [31] Melo, I. Higgs potential and fundamental physics. Eur. J. Phys. 2017, 38, 065404. https://doi.org/10.1088/1361-6404/aa8c3d.
  * [32] Trodden, M. Electroweak baryogenesis. Rev. Mod. Phys. 1999, 71, 1463. https://link.aps.org/doi/10.1103/RevModPhys.71.1463
  * [33] Branco, G.C.; Ferreira, P.M.; Lavoura, L.; Rebelo, M.N.; Sher, M.; Silva, Joao, P. Theory and phenomenology of two-Higgs-doublet models. Phys. Rep. 2021, 516, 1–102. https://doi.org/10.1016/j.physrep.2012.02.002
  * [34] Karmakar, S.; Rakshit, S. Effective Route Beyond the Extended Scalar Sectors of the Standard Model. Doctoral Dissertation, Discipline of Physics, IIT Indore, Indian Institute of Technology Indore, Indore, India, 2020.
  * [35] Karmakar S. Relaxed Constraints on Masses of New Scalars in 2HDM. Springer Proc. Phys. 2020, 248, 193–198. https://doi.org/10.1007/978-981-15-6292-1_23.
  * [36] Bernon, J.; Bian, L.; Jiang, Y. A new insight into the phase transition in the early Universe with two Higgs doublets. J. High Energy Phys. 2018, 5, 151. https://doi.org/10.1007/JHEP05(2018)151
  * [37] Basler, P.; Mühlleitner, M. BSMPT (Beyond the Standard Model Phase Transitions): A tool for the electroweak phase transition in extended Higgs sectors. Comput. Phys. Commun. 2019, 237, 62–85. https://doi.org/10.1016/j.cpc.2018.11.006
  * [38] Basler, P.; Mühlleitner, M.; Müller, J. BSMPT v2 A Tool for the Electroweak Phase Transition and the Baryon Asymmetry of the Universe in Extended Higgs Sectors. arXiv 2020, arXiv:2007.01725.