###### Abstract

Based on recent perturbative and non-perturbative lattice calculations with
almost quark flavors and the thermal contributions from photons, neutrinos,
leptons, electroweak particles, and scalar Higgs bosons, various thermodynamic
quantities, at vanishing net-baryon densities, such as pressure, energy
density, bulk viscosity, relaxation time, and temperature have been calculated
up to the TeV-scale, i.e. covering hadron, QGP and electroweak (EW) phases in
the early Universe. This remarkable progress motivated the present study to
determine the possible influence of the bulk viscosity in the early Universe
and to understand how this would vary from epoch to epoch. We have taken into
consideration first- (Eckart) and second-order (Israel-Stewart) theories for
the relativistic cosmic fluid and integrated viscous equations of state in
Friedmann equations. Nonlinear nonhomogeneous differential equations are
obtained as analytical solutions. For Israel-Stewart, the differential
equations are very sophisticated to be solved. They are outlined here as road-
maps for future studies. For Eckart theory, the only possible solution is the
functionality, $H(a(t))$, where $H(t)$ is the Hubble parameter and $a(t)$ is
the scale factor, but none of them so far could to be directly expressed in
terms of either proper or cosmic time $t$. For Eckart-type viscous background,
especially at finite cosmological constant, non-singular $H(t)$ and $a(t)$ are
obtained, where $H(t)$ diverges for QCD/EW and asymptotic EoS. For non-viscous
background, the dependence of $H(a(t))$ is monotonic. The same conclusion can
be drawn for an ideal EoS. We also conclude that the rate of decreasing
$H(a(t))$ with increasing $a(t)$ varies from epoch to epoch, at vanishing and
finite cosmological constant. These results obviously help in improving our
understanding of the nucleosynthesis and the cosmological large-scale
structure.

###### keywords:

Viscous Cosmology, Particle-theory and field-theory models of the early
Universe, Mathematical and relativistic aspects of cosmology, Thermodynamic
functions and equations of state

1 1 0 https://doi.org/ Early Universe Thermodynamics and Evolution in
Nonviscous and Viscous Strong and Electroweak epochs: Possible Analytical
Solutions Title Abdel Nasser Tawfik1,2†,‡ and Carsten Greiner 2,‡
Correspondence<EMAIL_ADDRESS>In the early Universe both QCD and
EW eras play an essential role in laying seeds for nucleosynthesis and even
dictating the cosmological large-scale structure. Taking advantage of recent
developments in ultrarelativistic nuclear experiments and nonperturbativ and
perturbative lattice simulations, various thermodynamic quantities including
pressure, energy density, bulk viscosity, relaxation time, and temperature
have been calculated up to the TeV-scale, in which the possible influence of
finite bulk viscosity is characterized for the first time and the analytical
dependence of Hubble parameter on the scale factor is also introduced.
98.80.-k, 98.80.Cq, 98.80.Jk, 05.70.Ce

##  1 Introduction

The current thorough knowledge on the cosmic evolution is primarily based on
the standard model of cosmology (SMC), which introduces a generic hypothesis
that the cosmic background is isotropically and homogeneously filled up with
an exclusively ideal fluid. After all, we simply realize that this is an
abstraction, i.e. a general description that isn’t based on real physical
situation. Apart from approaches, models, and theories, the real situation
could be arose out from recent high-energy experiments Heinz (2001); Tawfik
and Ganssauge (2000); Gyulassy and McLerran (2005); Heinz et al. (2012);
Adamczyk et al. (2014); Ryu et al. (2018); Bzdak et al. (2020) and cosmic
observations Komatsu et al. (2011); Ade et al. (2016); Aghanim et al. (2020);
Akrami et al. (2020); Aghanim et al. (2020); Akrami et al. (2020). Over years,
it was assumed that the impacts of the viscosity coefficients on cosmology
should be weak or at least subdominant so that the inclusion of viscous
concepts in the macroscopic theory of the cosmic fluid appeared as most
natural improvement. It was first assumed that the influence of viscosity in
the early Universe would be the largest at the the end of the lepton era, i.e.
during the neutrino decoupling era, at temperature $\simeq 10^{10}~{}$K.
Viscous coefficients connected with particle physics have been also proposed
by Misner Misner (1968); Zeldovich and Novikov (1983). Recent studies reveal
that the impact of viscosity likely sets on during the very early stages of
the Universe Tawfik and Mishustin (2019). The present study suggests extending
SMC to beyond SMC. In bCMS, the cosmic background geometry is filled with
viscous matter whatever its constituents are, so that isotropicity and
homogeneity are generalized.

For the inclusion of the viscous properties, one would like to start with
small perturbations from the thermal equilibrium. The suitable theoretical
framework for this is the first-, Sec. 4.2, and the second-order cosmic
relativitic fluid, Sec. 4.3. Both viscosity coefficients, the bulk viscosity
$\zeta$ and the shear viscosity $\eta$ can be determined. From SMC
considerations that the Universe is spatially homogeneously expands, and
cosmological observations Komatsu et al. (2011); Ade et al. (2016); Aghanim et
al. (2020); Akrami et al. (2020); Aghanim et al. (2020); Akrami et al. (2020),
$\zeta$ would be taken as a dominant component, while $\eta$ would be
neglected. In bCMS, this assumption could be also generalized. The motivations
for viscous theories in cosmology have are diversified. For instance, over the
last three decades, various attempts have been reported in literature Gron
(1990); Maartens (1995); Tawfik et al. (2010); Tawfik and Wahba (2010); Tawfik
et al. (2010); Adamczyk et al. (2019); Tawfik and Mishustin (2019). A direct
implementation of the equations of state (EoS) deduced from recent lattice
quantum chromodynamic calculations and/or heavy-ion collisions on physics of
the early Universe was initiated in various studies conducted by one of the
authors Tawfik and Harko (2012); Tawfik (2011); Tawfik and Magdy (2012);
Tawfik et al. (2011, 2010); Tawfik (2010); Tawfik and Greiner (2019); Tawfik
and Mishustin (2019). The present paper resumes these studies, especially in
light of the recent progress enabled us to explore the very early epochs of
the evolution of the Universe Laine and Schroder (2006); Laine and Meyer
(2015); D’Onofrio and Rummukainen (2016); Borsanyi et al. (2016); Tawfik and
Mishustin (2019). The procedure goes as follows. The viscous EoS introduced in
Section 4.1 and taken from Ref. Tawfik and Mishustin (2019) shall be
substituted in the Friedmann equations. This leads to sophisticated
differential equations. Their analytical solutions turn into a very
challenging mathematical task. By finding unambitious analytical solutions,
bSMC becomes a feasible approach. In this paper, we introduce and discuss the
possible analytical solutions; the ones expressing the Hubble parameter in
dependence on the scale factor, i.e. functionality $H(a(t))$, where both
quantities are also functions of the cosmic time $t$. We also introduce a
road-map for future studies based on bSMC.

The present script is organized as follows. The cosmic geometry and the field
equations will be reviewed in section 2. The cosmic evolution in non-viscous
background geometry classified into different epochs will be discussed in
section 3. The cosmic evolution in viscous background geometry, section 4, is
based on viscous EoS introduced in section 4.1, where the background fluid is
described by first-order Eckart theory, section 4.2 and second-order Israel-
Stewart theory, Section 4.3. The results on the possible analytical solutions,
i.e. functionality $H(a(t))$, where both $H$ and $a$ are also functions of the
cosmic time $t$. will be elaborated in Section 5. Section 6 is devoted to draw
the final conclusions.

##  2 Geometry and field equations

In curved cosmic geometry under the assumptions of SMC (homogeneity and
isotropy) for cosmic space and matter, the Friedmann-Lemaitre-Robertson-Walker
(FLRW) metric reads

$ds^{2}=dt^{2}-a(t)^{2}\left[\frac{dr^{2}}{1-kr^{2}}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\right],$ (1)

where $a(t)$ is the dimensionless scale factor, which describes the expansion
of the Universe. $k$ characterizes elliptical, flat (Euclidean), and
hyperbolic cosmic space, where $k=\\{-1,0,+1\\}$ stands for negative, flat,
and positive curvature, respectively. It should be noticed that if $r$ is
taken dimensionless, $a(t)$ shall be given in a unit of length. In Eq. (1) and
to simplify the cosmology notation, we use the natural units $c=G=1$. So-far,
the theory of general relativity doesn’t inter the play. It certainly does,
when the evolution of $s(t)$, the temporal evolution of the line element,
should be tackled. Towards this end, the theory of general relativity should
be combined with the matter/energy content of the space-time within the cosmic
geometry.

The Einstein gravitational fields with finite cosmological constant are given
as

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\,R+\Lambda_{\mu\nu}=\frac{8\pi}{3}\,T_{\mu\nu},$
(2)

where the indices $\mu$, $\nu$ take discrete values $0$, $1$, $2$, $3$. The
energy-momentum tensor of the bulk viscous cosmological fluid filling the very
early Universe can be expressed as Maartens (1995)

$T_{\mu\nu}=\left(\rho+p+\Pi\right)\,u_{\mu}u^{\nu}-\left(p+\Pi\right)g_{\mu\nu},$
(3)

where $\rho$ is the energy density, $p$ is the thermodynamic pressure, $\Pi$
is the bulk viscous pressure, and $u_{\mu}$ is the four velocity satisfying
the normalization condition $u_{\mu}u^{\mu}=1$. The bulk pressure $\Pi$ can
formally be included in the thermodynamic pressure $p_{\mathtt{eff}}=p+\Pi$.
We shall discuss on how to evaluate $\Pi$, concretely the bulk viscous
pressure, in framework of Eckart (first-order), section 4.2, and Israel-
Stewart (second-order) theories, section 4.3, for relativistic viscous cosmic
fluid.

For number density $n$, specific entropy $s$, finite temperature $T$, bulk
viscosity coefficient $\zeta$, and relaxation time $\tau$, the particle and
entropy fluxes are to be related to each other as $N^{i}=\nu^{i}$ and
$S^{i}=sN^{i}-\left(\tau\Pi^{2}/2\zeta T\right)u^{i}$, respectively. It should
be emphasized that the evolution of the cosmological fluid is subject to the
dynamical laws of particle number conservation $N_{\;;i}^{i}=0$ and Gibbs’
equation $Td\rho=d\left(\rho/n\right)+pd\left(1/n\right)$ Maartens (1995). In
what follows, we assume that the energy-momentum tensor of the cosmological
fluid is locally conserved, i.e. $T_{i;k}^{k}=0$, where ; denotes the
covariant derivative with respect to the line metric.

In the proper frame, i.e. the inertial frame of reference comoving with the
the fluid, the components $T_{0}^{0}=\rho$,
$T_{1}^{1}=T_{2}^{2}=T_{3}^{3}=-p_{\mathtt{eff}}$. For the isotropic and
homogeneous metric given in Eq. (1), the Einstein field equations in natural
units read

$\displaystyle H(t)^{2}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{3}\;\rho(t)-\frac{k}{a(t)^{2}}+\frac{\Lambda}{3},$
(4) $\displaystyle\dot{H}(t)+H(t)^{2}$ $\displaystyle=$
$\displaystyle-\frac{4\pi}{3}\;\left[\rho(t)+3p_{\mathtt{eff}(t)}\right]+\frac{\Lambda}{3},$
(5)

where the dot refers to differentiation with respect to the cosmic time $t$,
and $H(t)=\dot{a}(t)/a(t)$ is the Hubble parameter. From the expressions (4)
and (5), the time evolution of Hubble parameter can be deduced as

$\displaystyle\dot{H}(t)$ $\displaystyle=$
$\displaystyle-4\pi\,\left[\rho(t)+p_{\mathtt{eff}}(t)\right]+\frac{k}{a(t)^{2}}.$
(6)

From the local conservation of the energy-momentum tensor in the Universe,
following equation has been proposed by McCrea and Milde and by Peebles with
vanishing McCrea and Milne (1934, 2000) and finite pressure Peebles (1993),
respectively, as being equivalent to Newtonian mechanics,

$\dot{\rho}(t)+3\left[\rho(t)+p_{\mathtt{eff}}(t)\right]H(t)=0.$ (7)

This means that the decrease in the energy content of a cube with side $a(t)$
equals the energy budget due to the expansion of the Universe and the work
done by the pressure on the surface.

To obtain a closed system of equations, we need to propose EoS relating $p$ to
$\rho$. Depending on the approach we are applying, we might also need to
propose a reliable estimation for $\Pi$, as well. In the section that follows,
we introduce solution for the Friedmann equation based on various types of
EoS. We first assume vanishing bulk viscosity, Section 3. Then, we discuss on
the extension to finite bulk viscosity, Section 4, which then required
barotropic equations for the pressure, Eq. (8), (9), (10), the temperature,
Eq. (37), the bulk viscosity coefficient, Eqs (34), (35), (36) and the
relaxation time, Eq. (39) (40), (41). These are examples on novel
contributions presented by the present script.

##  3 Cosmic evolution in non-viscous approach

By combining recent non-perturbative and perturbative calculations with other
degrees of freedom (dof), such as photons, neutrinos, leptons, electroweak
particles, and Higgs bosons, various thermodynamic quantities for almost net-
baryon-free cosmic matter have been calculated up to the TeV-scale, i.e.
covering quantum chromodynamic (QCD) and electroweak (EW) eras of the early
Universe Tawfik and Mishustin (2019). It was found that while the EoS relating
the pressure $p$ to the energy density $\rho$ for the hadronic matter is
simple, the one for QCD and that EW matter are rather complicated. It is worth
highlighting that in the cosmological context, the various thermodynamic
quantities should be translated into time-depending quantities. When the
cosmic time elapses, the spacial dimensions of the Universe expand, and
accordingly thermodynamic quantities characterizing the background geometry
vary. The EoS proposed in Fig. 1 Tawfik and Mishustin (2019) are preliminary
depending on the energy density, whose decrease might be - for simplicity -
taken as a scale for increasing cosmic time and vice versa.

$\displaystyle\mathtt{Hadron:}\qquad\qquad p(t)$ $\displaystyle=$
$\displaystyle\alpha_{1}+\beta_{1}\rho(t),$ (8)
$\displaystyle\mathtt{QCD/EW:}\qquad\qquad p(t)$ $\displaystyle=$
$\displaystyle\alpha_{2}+\beta_{2}\rho(t)+\gamma_{2}\rho(t)^{\delta_{2}},$ (9)

where $\alpha_{1}=0.0034\pm 0.0023$, $\beta_{1}=0.1991\pm 0.0022$,
$\alpha_{2}=0.0484\pm 0.0164$, $\beta_{2}=0.3162\pm 0.031$,
$\gamma_{2}=-0.21\pm 0.014$, and $\delta_{2}=-0.576\pm 0.034$. At very large
$\rho(t)$, the asymptotic behavior becomes very close to that of an ideal gas
limit,

$\displaystyle\mathtt{Asymp.:}\qquad\qquad p(t)$ $\displaystyle=$
$\displaystyle\gamma_{3}\rho(t),$ (10)

where $\gamma_{3}=0.3304\pm 0.0236$. It should be noticed that the EoS (8), in
the hadronic phase, with its positive parameters $\alpha_{1}$ and $\beta_{1}$
could easily be - due to its large uncertainty - reexpressed with vanishing
$\alpha_{1}$. Nevertheless, in the present calculations, we keep $\alpha_{1}$
finite. In section 3, we present solutions for the Friedmann equations, Eqs.
(4), (5), (6), with the various Eqs. (9)-(10), which as mentioned characterize
various types of cosmic backgrounds corresponding various epochs of the early
Universe. The results obtained for the dependence of the Hubble parameter on
the scale factor are presented in Fig. 4.

At vanishing bulk viscosity, $\Pi(t)=0$, the effective pressure
$p_{\mathtt{eff}}(t)=p(t)+\Pi(t)$ can be simplified as the thermodynamic
pressure $p(t)$. Then, the Friedmann equation (6) can be rewritten as

$\displaystyle\ddot{a}(t)\,a(t)-\dot{a}(t)^{2}+4\pi\left[\rho(t)+p(t)\right]{a}(t)^{2}-k$
$\displaystyle=$ $\displaystyle 0,$ (11)

which can be solved if combined with set of closed equations, such as Eq. (4)
and suitable EoS. Accordingly, we have various solutions characterizing the
various eras in the early Universe.

Fig. 1: The pressure is depicted as a function of the energy density. Both
quantities are given in GeV/fm3 units. The dashed lines present the various
parameterizations (see text).

###  3.1 Hadronic Era

By substituting Eq. (4) and Eq. (8) into Eq. (11), we get

$\displaystyle\ddot{a}(t)\,a(t)+C_{1}\;\dot{a}(t)^{2}+C_{2}\;{a}(t)^{2}+C_{1}\;k$
$\displaystyle=$ $\displaystyle 0,$ (12)

where the variables $C_{1}$ and $C_{2}$ are given in Tab. 1. $C_{1}$ and
$C_{2}$ are functions of the coefficients obtained in the parameterized EoS,
which in tern vary from epoch to epoch. For the sake of simplicity, the
coefficients are conjectured remaining constant within each epoch. We notice
that the value of $C_{1}$ \- in the hadron era - is finite but not necessarily
unity. This assures that $k$, the curvature contact, remains finite.

  * •

At vanishing $k$, which is the case at $\beta_{1}=-1/3$, we have analytical
solutions. Then, Eq. (12) can be solved as,

$\displaystyle a(t)$ $\displaystyle=$ $\displaystyle
c_{2}\cosh\left[\sqrt{C_{2}(1-C_{1})}(t+c_{1})\right]^{1/(1-C_{1})},$ (13)

where $c_{1}$ and $c_{2}$ are integration constants which can be fixed at
boundary and initial conditions. For instance, at $t=0$, $H(t)=0$, Eq. (14),
then, $c_{1}=-t$. In general $c_{1}$ has the dimension of the cosmic time $t$
and therefore varies with the evolution of the Universe. Hence, $c_{1}$ is
finite and the cosmological parameters of the hadron epoch ($\sim
10^{-6}-10^{-5}~{}$s or111assuming that $6.58\times
10^{-16}\mathtt{s}=\mathtt{eV}^{-1}$ $\sim 15.197\times 10^{14}~{}$MeV-1), for
instance Eq. (14) can be estimated, numerically. On the other hand, for the
scale factor, $a(t)$, we still need to estimate the other integration constant
$c_{2}$. Having an analytical expression for the scale factor, Eq. (13), then
the Hubble parameter can then be obtained,

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{\dot{a}(t)}{a(t)}=\sqrt{\frac{C_{2}}{1+C_{1}}}\;\tanh\left[\sqrt{C_{2}(1+C_{1})}(t+c_{1})\right].$
(14)

  * •

At non-vanishing $k$, there is no direct analytical solution for $a(t)$. But
when assuming that $u=\dot{a}^{2}(t)$ and substituting this into Eq. (12),

$\displaystyle\frac{du(a(t))}{da(t)}+2C_{1}\,\frac{u(a(t))}{a(t)}+2C_{2}a(t)+2C_{1}\frac{k}{a(t)}$
$\displaystyle=$ $\displaystyle 0,$ (15)

a solution for $\dot{a}(t)$ can be proposed

$\displaystyle u(a(t))$ $\displaystyle=$ $\displaystyle\dot{a}(t)^{2}$ (16)
$\displaystyle=$ $\displaystyle
c_{1}a(t)^{-2C_{1}}-\frac{C_{2}}{1+C_{1}}a(t)^{2}-\frac{k}{C_{1}+1}a(t)^{2}.$

The physical solution is the one assuring that,

$\displaystyle a(t)$ $\displaystyle<$
$\displaystyle\left[\frac{c_{3}}{C_{2}}(1+C_{1})\right]^{\frac{1}{2(1-C_{1})}}.$
(17)

Hence, the Hubble parameter can be deduced as

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\left\\{c_{1}a(t)^{-2C_{1}}-\frac{C_{2}+k}{C_{1}+1}a(t)^{2}\right\\}^{1/2}\frac{1}{a(t)}.$
(18)

By solving the second-order differential equation (16), an analytical
expression for the scale factor $a(t)$ can also be deduced,

$\displaystyle a(t)$ $\displaystyle=$
$\displaystyle\left[\sqrt{\frac{C_{2}+k}{c_{1}(C_{1}+1)}}\right]^{\frac{1}{-(1+C_{1})}},$
(19)

whose time dependence is given by the time dependence of the corresponding
EoS, namely $C_{1}$ and $C_{2}$, which are listed in Tab. 1. As discussed they
have an indirect time dependence through the coefficients of the corresponding
EoS. But within one era, they are conjectured to remain constant. The latter
might be the only way possible to gain an analytical solution. Then, $H(t)$,
Eq. (18), can be rewritten as

$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle
c_{1}\left[\left[\sqrt{\frac{C_{2}+k}{c_{1}(C_{1}+1)}}\right]^{-\frac{1}{1+C_{1}}}\right]^{-(1+2C_{1})}+\frac{C_{2}+k}{C_{1}+1}\left[\sqrt{\frac{C_{2}+k}{c_{1}(C_{1}+1)}}\right]^{\frac{1}{-(1+C_{1})}}.\hskip
22.76219pt$ (20)

Apparently, all coefficients involved in can be determined.

###  3.2 QCD and EW Era

Based on the EoS outlined in Eq. (9), and the dependence of energy density on
the Hubble parameter, Eq. (4), then Eq. (11) can be reexpressed as

$\displaystyle\ddot{a}(t)\,a(t)+C_{1}\;\dot{a}(t)^{2}+C_{2}\;{a}(t)^{2}+C_{1}\;k+4\pi\gamma_{2}a(t)^{2}\rho^{\delta_{2}}$
$\displaystyle=$ $\displaystyle 0,$ (21)

where $C_{1}$ and $C_{2}$ are variables depending on the corresponding EoS,
Tab. 1. The last term in lhs of this expression gives another difference with
Eq. (12). The other terms remaining can be estimated as outlined in the
previous section. Now we focus on the the contributions added by this term,

$\displaystyle
4\pi\gamma_{2}a(t)^{2}\rho^{\delta_{2}}=4\pi\gamma_{2}a(t)^{2(1-\delta_{2})}\left(\frac{3}{8\pi
k}\right)^{\delta_{2}}\left[1-\left(\frac{\Lambda
a(t)^{2}}{3k}-\frac{\dot{a}(t)^{2}}{k}\right)\right]^{\delta_{2}}.$ (22)

For $\delta_{2}=-0.576\simeq-0.5$, the square bracket can be expressed as a
binomial expansion,

$\displaystyle\left[1-\left(\frac{\Lambda
a(t)^{2}}{3k}-\frac{\dot{a}(t)^{2}}{k}\right)\right]^{-1/2}=1-\frac{\Lambda
a(t)^{2}}{6k}-\frac{\dot{a}(t)^{2}}{2k}+\cdots.$ (23)

Thus, we might approximate the entire bracket to first terms outlined. This
result can also be obtained when assuming that the exponent $\delta_{2}$
approaches unity. Then, Eq. (21) becomes

$\displaystyle\ddot{a}(t)\,a(t)+C_{1}\dot{a}(t)^{2}+C_{2}a(t)^{2}+C_{1}k+C_{3}a(t)^{2}\left(1-C_{4}\frac{a(t)^{2}}{k}-\frac{\dot{a}(t)^{2}}{3k}\right)$
$\displaystyle=$ $\displaystyle 0.$ (24)

As done while solving Eq. (15), we assume that $u(a(t))=\dot{a}(t)^{2}$. Then,
the approximated Eq. (24) becomes

$\displaystyle\frac{du(a(t))}{da(t)}+2C_{1}\frac{u(a(t))}{a(t)}+2C_{2}a(t)+2C_{1}k\,a(t)^{-1}+$
$\displaystyle
2C_{3}a(t)\left(1-C_{4}\frac{a(t)^{2}}{k}-\frac{u(a(t))}{3k}\right)=0,\hskip
14.22636pt$ (25)

which can be solved as

$\displaystyle u(t)$ $\displaystyle=$ $\displaystyle\dot{a}(t)^{2}$ (26)
$\displaystyle=$
$\displaystyle\frac{1}{C_{3}}\left\\{3k\left[C_{2}+C_{3}-3\left(1+C_{1}\right)C_{4}\right]-3C_{3}C_{4}a(t)^{2}+c_{1}a(t)^{-2C_{1}}e^{\frac{C_{3}}{3k}a(t)^{2}}\right.$
$\displaystyle\left.-C_{1}k\left[-3C_{2}-4C_{3}+9\left(1+C_{1}\right)C_{4}\right]e^{\frac{C_{3}}{3k}a(t)^{2}}Ein_{1-C_{1}}\left(\frac{C_{3}}{3k}a(t)^{2}\right)\right\\},\hskip
22.76219pt$

where $c_{1}$ is another integration constant to be fixed for boundary
conditions and the exponential integral represents a special case of the
incomplete gamma function

$\displaystyle
Ein_{n}\left(x\right)=\left(x\right)^{n-1}\Gamma\left[1-n,y\right].$ (27)

For equation (26), there is no analytical solution. Nevertheless, the Hubble
parameter, $H(t)=\dot{a}(t)/a(t)$, can be constructed as

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{1}{C_{3}^{1/2}a(t)}\left\\{3k\left[C_{2}+C_{3}-3\left(1+C_{1}\right)C_{4}\right]-3C_{3}C_{4}a(t)^{2}+c_{4}a(t)^{-2C_{1}}e^{\frac{C_{3}}{3k}a(t)^{2}}\right.$
(28)
$\displaystyle\left.-C_{1}k\left[-3C_{2}-4C_{3}+9\left(1+C_{1}\right)C_{4}\right]e^{\frac{C_{3}}{3k}a(t)^{2}}Ein_{1-C_{1}}\left(\frac{C_{3}}{3k}a(t)^{2}\right)\right\\}^{1/2}.$

$C_{3}$ and $C_{4}$ are given in Tab. 1. The results obtained for the Hubble
parameter as a function of the scale factor shall be presented in Fig. 4, in
which the exponential function, Eq. (27), is - for the sake of simplicity -
assigned to the unity.

###  3.3 Asymptotic Limit

Again when substantiating Eq. (4) and Eq. (10) into Eq. (11), we get

$\displaystyle\ddot{a}(t)\,a(t)+C_{1}\,\dot{a}(t)^{2}+C_{2}\,a(t)^{2}+C_{1}\,k$
$\displaystyle=$ $\displaystyle 0,$ (29)

which apparently looks almost identical to Eq. (12) in section 3.1. The
possible analytical solution reads

$\displaystyle u(t)$ $\displaystyle=$ $\displaystyle\dot{a}(t)^{2}$ (30)
$\displaystyle=$ $\displaystyle
c_{4}a(t)^{-2C_{1}}-\left[\frac{C_{2}}{C_{1}+1}+k\right]a(t)^{2},$

for which the Hubble parameter can be given as

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\left\\{c_{4}a(t)^{-2C_{1}}-\left[\frac{C_{2}}{C_{1}+1}+k\right]a(t)^{2}\right\\}^{1/2}\frac{1}{a(t)}.$
(31)

The results of this expression are given in Fig. 4. By integrating (30), an
expression for the scale factor can be obtained

$\displaystyle a(t)$ $\displaystyle=$
$\displaystyle\left[\sqrt{\frac{C_{2}+k(1+C_{1})}{c_{1}(1+C_{1})}}\right]^{\frac{1}{-(1+C_{1})}},$
(32)

which helps in constructing the corresponding Hubble parameter

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\left[\frac{C_{2}+k(1+C_{1})}{c_{1}(1+C_{1})}\right]^{\frac{1}{1+C_{1}}}\left\\{-\frac{C_{2}+k(1+C_{1})}{1+C_{1}}\left[\sqrt{\frac{C_{2}+k(1+C_{1})}{c_{1}(1+C_{1})}}\right]^{-\frac{2}{1+C_{1}}}\right.$
(33) $\displaystyle+$
$\displaystyle\left.\left[\left[\sqrt{\frac{C_{2}+k(1+C_{1})}{c_{1}(1+C_{1})}}\right]^{-\frac{1}{1+C_{1}}}\right]^{-2C1}\right\\}^{1/2}.$

The various coefficients characterizing the various EoS and also combining
cosmological constant, $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$, are listed in Tab.
1. Accordingly, analytical solutions similar to Eq. (13) are obtained. The
cosmological constant $\Lambda$ is conjectured to count for the dark energy
component Shalyt-Margolin (2010). For Eq. (29), expressions for Hubble
parameter similar to (14) can then be derived. Also, with the variable change
$u=\dot{a}(t)^{2}$, (18) can be obtained, as well.

| $\mathtt{Section}$ | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}$
---|---|---|---|---|---
$\mathtt{Hadron}$ | 3.1 | $\frac{3}{2}(1+\beta_{1})-1$ | $4\pi\alpha_{1}-\frac{1}{2}(1-\beta_{1})\Lambda$ |  |
$\mathtt{QCD/EW}$ | 3.2 | $\frac{3}{2}(1+\beta_{2})-1$ | $4\pi\alpha_{2}-\frac{1}{2}(1-\beta_{2})\Lambda$ | $4\pi\left(\frac{3}{8\pi k}\right)^{\delta_{2}}\gamma_{2}$ | $\frac{\Lambda}{6}$
$\mathtt{Asymp.}$ | 3.3 | $\frac{3}{2}(1+\gamma_{3})-1$ | $-\frac{1}{2}(1+\gamma_{3})\Lambda$ |  |
Tab. 1: The parameters defining different solutions for Hadron, QCD/EW and
Asympt. phases corresponding to various EoS, (12), (21), and (29),
respectively.

##  4 Cosmic evolution in viscous approaches

###  4.1 Viscous equations of state

The recent results for the bulk viscosity are based on non-perturbative and
perturbative calculations with as much quark flavors as possible. By combining
these calculations with additional dof, such as photons, neutrinos, leptons,
electroweak particles, and Higgs bosons, various thermodynamic quantities
including bulk viscosity, for almost net-baryon-free cosmic matter, have been
calculated up to the TeV-scale Tawfik and Mishustin (2019). The dependence of
the bulk viscosity on the energy density Tawfik and Greiner (2019) is depicted
in the top panel of Fig. 2, in which both quantities are given in the physical
units. As discussed, such a barotropic dependence straightforwardly allows for
direct cosmological implications Tawfik and Harko (2012); Tawfik et al. (2010,
2011, 2010), where $\rho(t)$ can be directly substituted by $H(t)$, Eq. (4).
These wide values of $\rho(t)$ which are accompanied by a wide range of
temperatures cover quantum chromodynamic (QCD) (Hadron and QGP) and
electroweak (EW) phases in the early Universe. Accordingly, the dependence of
the bulk viscosity on the energy density can be parameterized.

$\displaystyle\mathtt{Hadron-QGP:}$
$\displaystyle\zeta(t)=d_{1}+d_{2}\rho(t)+d_{3}\rho(t)^{d_{4}},$ (34)
$\displaystyle\mathtt{QCD:}$
$\displaystyle\zeta(t)=e_{1}+e_{2}\rho(t)^{e_{3}},$ (35)
$\displaystyle\mathtt{EW:}$
$\displaystyle\zeta(t)=f_{1}+f_{2}\rho(t)^{f_{3}}.$ (36)

The various fit parameters are given as follows. For Hadron-QCD:
$d_{1}=-9.336\pm 4.152$, $d_{2}=0.232\pm 0.003$, $d_{3}=11.962\pm 4.172$, and
$d_{4}=0.087\pm 0.029$. For QCD: $e_{1}=8.042\pm 0.056$, $e_{2}=0.301\pm
0.002$, and $e_{3}=0.945\pm 0.0001$. For EW: $f_{1}=0.350\pm 0.065$,
$f_{2}=10.019\pm 0.934$, and $f_{3}=0.929\pm 8.898\times 10^{-5}$.

While $\zeta(t)$ vs. $\rho(t)$ is much structured in the hadron era, there are
three domains to be emphasized (from low to large energy density).

  * •

The first one is the hadron-QGP domain (Hadron-QGP), which spans over
$\rho(t)\lessapprox 100~{}$GeV/fm3. At the beginning, there is a rapid
increase in $\zeta(t)$, i.e. $\zeta\eqsim 1~{}$GeV3, at $\rho(t)\simeq
1~{}$GeV/fm3, which is then followed by a slight increase in $\zeta(t)$. For
example, at $\rho(t)\simeq 100~{}$GeV/fm3, $zeta(t)$ reaches $\sim
130~{}$GeV3. It is apparent that the hadron-parton phase transition seems to
take place at $\rho(t)\lessapprox 0.5~{}$GeV/fm3 Tawfik (2005a, b). At this
value, $\zeta(t)\lessapprox 0.5~{}$GeV3.

  * •

The second domain, the QGP epoch, seems to be formed, at
$0.5\lessapprox\rho(t)\lessapprox 100~{}$GeV/fm3, i.e. a much wider $\rho(t)$
than that of the hadron domain. Thus, we could conclude that over this wide
range of $\rho(t)$, the bulk viscosity is obviously not only finite but rather
largely supporting the RHIC discovery of strongly correlated viscous QGP Ryu
et al. (2018); Heinz et al. (2012); Gyulassy and McLerran (2005). At higher
$\rho(t)$, we observe a tendency of a linear increase in $\zeta(t)$ with
further increasing $\rho(t)$. Thus, the second domain is the one where
$100\gtrapprox\rho(t)\lessapprox 5\times 10^{7}~{}$GeV/fm3 and
$80\gtrapprox\zeta(t)\lessapprox 10^{6}~{}$GeV3. In light on this observation,
we conclude that the phase transition from QCD to EW domain is very smooth.

  * •

The third domain is also characterized by an almost linear increase in
$\zeta(t)$ with increasing $\rho(t)$. For $10^{8}\lessapprox\rho(t)\lessapprox
10^{15}~{}$GeV/fm3, there is a nearly steady increase in $\zeta(t)$ from
$10^{8}$ to $10^{14}~{}$GeV3.

For temperatures ranging from a few MeV to TeV and energy densities up to
$10^{16}~{}$GeV/fm3, we have taken into consideration almost all possible
contributions to the bulk viscosity. With these we mean the thermodynamic
quantities calculated in non-perturbation and perturbation QCD with up, down,
strange, charm, and bottom quark flavors. The second type of contributions is
the guage bosons, the entire gluonic sector. We have also included photons,
$W^{\pm}$, and $Z^{0}$, charged leptons (neutrino, electron, muon, and tau),
and scalar Higgs particle. The third type of contributions is the vacuum and
thermal condensations. We have included condensations for up, down, strange
and charm quarks. We merely still miss the vacuum and the thermal bottom quark
condensates, besides the entire gravitational, the neutral leptons, and the
top quark sector to compile the entire standard model for elementary
particles.

Fig. 2: Top panel depicts the energy-density dependence of the bulk viscosity.
Bottom panel illustrates the temperature as a function of energy density. The
parameterizations are depicted as curves.

The dependence of temperature $T(t)$ on the energy density $\rho(t)$, the
barotropic equation of state, is depicted in bottom panel of Fig. 2. We notice
that $T(t)$ almost linearly depends on $\rho(t)$. A best parametrization reads

$\displaystyle T(t)=\alpha_{4}+\beta_{4}\rho(t)^{\gamma_{4}},$ (37)

where $\alpha_{4}=0.048\pm 0.001$, $\beta_{4}=0.13\pm 2\times 10^{-4}$, and
$\gamma_{4}=0.25\pm 8\times 10^{-5}$. We notice that at low $\rho(t)$ the
temperature looks a little bit structured. While with increasing $\rho(t)$,
the temperature goes almost linearly with increasing $\rho(t)$, especially at
very high temperatures, where $\rho(t)$ becomes related in $T^{4}$, i.e. ideal
gas.

The third quantity, for which we need to propose a barotropic EoS, is the
relaxation time, $\tau(t)$. We assume to apply the phenomenological model
presented in refs. Maartens (1995); Pun et al. (2008); Tawfik et al. (2010,
2011), which is based on dissipative relativistic fluid. This model was
assumed to characterize the evolution of the Universe with a flat homogeneous
isotropic Friedmann-Robertson-Walker geometry filled with viscous cosmic
fluid, but still valid for other types of curvature and cosmic backgrounds.
Accordingly, we have

$\displaystyle\tau(t)$ $\displaystyle=$
$\displaystyle\frac{\zeta(t)}{\rho(t)}.$ (38)

Fig. 3: The energy-density dependence of the relaxation time. The
parameterizations, Eqs. (39), (40), (41), are depicted as curves.

Figure 3 shows the energy-density dependence of the relaxation time $\tau(t)$.
Bearing in mind the linear dependence of the energy density $\rho(t)$ on the
temperature $T^{4}$, bottom panel of Fig. 2 and Eq. (37), the temperature
dependence $\tau(t)$ can be almost straightforwardly estimated. As done in the
present script, we would like to distinguish between hadron-QGP (squares), QCD
(circles) and electroweak eras (diamonds). In the hadron-QGP era, there is a
very rapid decrease in $\tau(t)$ with increasing $rho(t)$. The QCD epoch is
characterized by a slower decline in $\tau(t)$ with increasing $rho(t)$. The
relaxation time within the electroweak epoch starts and ends with a slow
decrease, while in the middle EW era, $\tau(t)$ rapidly decreases with
increasing $\rho(t)$. This region likely characterizes the electroweak phase
transition222In this case, $\tau(t)$ is conjectured to play the role of an
order parameter. For the electroweak phase transition, other thermodynamic
order parameters should be proposed and then analyzed.. The dependence of
$\tau(t)$ on $rho(t)$ is proposed as follows.

$\displaystyle\mathtt{Hadron-QGP:}$
$\displaystyle\tau(t)=g_{1}+g_{2}\exp(-g_{3}\rho(t)^{g_{4}}),$ (39)
$\displaystyle\mathtt{QCD:}$
$\displaystyle\tau(t)=h_{1}+\frac{h_{2}}{h_{3}+\log(h_{4}\rho(t))},$ (40)
$\displaystyle\mathtt{EW:}$
$\displaystyle\tau(t)=k_{1}\rho(t)^{k_{2}}\log(k_{3}\rho(t)),$ (41)

where $g_{1}=0.0002\pm 1.945\times 10^{-5}$, $g_{2}=0.008\pm 0.001$,
$g_{3}=1.671\pm 0.097$, $g_{4}=0.312\pm 0.0226$, $h_{1}=-1.605\times
10^{-7}\pm 2.556\times 10^{-6}$, $h_{2}=0.0015\pm 0.0013$, $h_{3}=0.935\pm
1.084\times 10^{5}$, $h_{3}=10.524\pm 1.141\times 10^{6}$, $k_{1}=9.582\times
10^{-4}\pm 9.504\times 10^{-5}$, $k_{2}=0.216\pm 0.0035$, and $k_{3}=5.9\times
10^{-7}\pm 9.165\times 10^{-8}$.

In the sections that follow, we apply well-know theories for relativistic
dissipative fluid to the cosmic background. We start with the Eckart
relativistic theory of a simple dissipative fluid, which is used to simplify
the nuclear motion arising in the second Born–Oppenheimer approximation. The
cosmic relevance of this theory is remarkable because it introduces the so-
called Eckart frame, which is a frame of orthonormal vectors following a
vibro-rotating object. The orientation of this frame is governed by the so-
called Eckart conditions assuring minimal Coriolis interaction. Second we
apply the Israel-Stewart theory as this theory is conjectured to solve Eckart
theory’s lack of causality and its obvious instabilities by introducing a
second-order term to the entropy.

###  4.2 Eckart relativistic viscous fluid

For cosmological context, the first theory of relativistic dissipative fluid
has been presented by Eckart Eckart (1940) and Landau and Lifshitz Landau and
Lifshitz (1987). It was pointed out that regardless the choice of EoS, the
equilibrium states of this theory are found unstable Osada (2012). In this
theory, only the first-order deviation from the equilibrium is taken into
consideration. But this leads to the superluminal velocities of the
dissipative signals, i.e. signals propagate through the relativistic
dissipative fluid with velocities exceeding the speed of light $c$ and hence
the theory violates the causality principle Piattella et al. (2011). Moreover,
it was shown that the resulting equilibrium states are unstable Israel (1976).
All these severe problems are originated from the fact that the Eckart theory
merely considers first-order deviations from the equilibrium leading to
parabolic differential equations, Eq. (45). The applicability of this theory
can only be thought for quasi-stationary phenomena, i.e. temporally and
spatially slowly varying, which are characterized by mean free-path and mean
collision-time.

The Eckart theory introduces a linear relationship between the bulk viscous
pressure and the rate of expansion of the Universe Coley and van den Hoogen
(1995). Obviously, this feature - despite the severe contains - makes it
possible to work out an analytical method for the cosmic parameters in the
expanding Universe. For bulk viscous cosmic fluid, whose energy-momentum
tensors are given as

$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle(\rho+p+\Pi)\,u_{\mu}\,u^{\nu}-(p+\Pi)\delta_{\mu}^{\nu},$ (42)

the line element in flat homogeneous isotropic Friedmann-Robertson-Walter
metric reads

$\displaystyle ds^{2}$ $\displaystyle=$ $\displaystyle
dt^{2}+a(t)^{2}\left[dx^{2}+dy^{2}+dz^{2}\right].$ (43)

When applying Eckart theory on modeling such a cosmic fluid, we assume
averaged 4-velocity fields $u^{\alpha}$ with $u^{\alpha}u_{\alpha}=1$ and
vector number density $n^{\alpha}=n\,u^{\alpha}$. For unbalanced
creation/annihilation processes in gravitational fields,
$n^{\alpha}_{;\alpha}=0$,

$\dot{n}+3\,H\,n=0,$ (44)

with Hubble parameter $H=u^{\alpha}_{;\alpha}$. In this theory, the entropy
current is given as

$\displaystyle S^{\alpha}$ $\displaystyle=$ $\displaystyle s\,n\,u^{\alpha}.$
(45)

As discussed, this is a non-conserved quantity. The covariant form of second
law of thermodynamics reads $S^{\alpha}_{;\alpha}\geq 0$. The divergence of
this quantity is given as $TS^{\alpha}_{;\alpha}=-3H\Pi$.

With respect to the proposed cosmic fluid, the temporal evolution can be
related to the dynamical laws of the particle number conservation
$N_{;i}^{i}=0$. Gibbs equation implies that
$Td\rho=d\left(\rho/n\right)+pd\left(1/n\right)$. Then, the covariant entropy
current assures a linear first-order relationship between the thermodynamical
flux $\Pi(t)$ and the corresponding $H(t)$. It is worth highlighting that
$H(t)$ in this context plays the role of a force Tawfik et al. (2013)

$\Pi(t)=-3\,\zeta(t)\,H(t).$ (46)

Having an estimation of the bulk viscous pressure $\Pi$, we can now substitute
Eq. (46) in Eq. (6),

$\displaystyle\dot{H}(t)$ $\displaystyle=$
$\displaystyle-4\,\pi\,\left[\rho(t)+\,p(t)-3\,\zeta(t)\,H(t)\right]+\frac{k}{a(t)^{2}},$
(47)

which can be rewritten as

$\displaystyle\ddot{a}(t)\,a(t)-\left[1+12\pi\zeta(t)\right]\dot{a}(t)^{2}+4\pi\,\left[\rho(t)+p(t)\right]a(t)^{2}-k$
$\displaystyle=$ $\displaystyle 0.$ (48)

For the present calculations, we start with Eq. (48), in which we substitute
with the barotropic EoS of the pressure and the bulk viscosity. They can then
be related to the Hubble parameter, Eq. (4). Due to the various barotropic EoS
introduced, Eq. (8), Eq. (9), Eq. (10) and section 4.1, the evolution of the
various cosmological parameters obviously differ from epoch to epoch. Equation
(48) combines extended assumptions and ingredients of SMC. This is viscous
cosmic background.

The sections that follow elaborate details on the dependence of $H(t)$ on
$a(t)$, where both quantities are functions of the cosmic time $t$. Such a
limitation is merely based on the currently available analytical solutions.
When lifting such mathematical limitations, bSMC emerges as a proper
cosmological approach.

####  4.2.1 Hadron-QGP era

Substituting with the pressure, Eq. (8), and the bulk viscosity, Eq. (34) into
Eq. (48) leads to a second-order differential equation

$\displaystyle\ddot{a}(t)a(t)-D_{1}\dot{a}(t)^{2}+D_{2}a(t)^{2}+D_{3}k$
$\displaystyle-$ $\displaystyle
D_{4}\left[1-D_{5}a(t)^{2}-\frac{\dot{a}(t)^{2}}{2k}\right]-\frac{9}{2}d_{2}a(t)^{-2}\dot{a}(t)^{4}-\frac{9}{2}d_{2}a(t)^{-2}\dot{a}(t)^{2}$
$\displaystyle=$ $\displaystyle 0,$ (49)

where the various coefficients of the different EoS are included in the
following parameters

$\displaystyle D_{1}$ $\displaystyle=$ $\displaystyle 1+12\pi
d_{1}-\frac{3}{2}(1+\beta_{1})-\frac{3}{2}d_{2}\Lambda,$ $\displaystyle D_{2}$
$\displaystyle=$ $\displaystyle
4\pi\alpha_{1}-\frac{\Lambda}{2}(1+\beta_{1}),$ $\displaystyle D_{3}$
$\displaystyle=$ $\displaystyle\frac{3}{2}(1+\beta_{1})-1,$ $\displaystyle
D_{4}$ $\displaystyle=$ $\displaystyle
2^{2-3d_{4}}3^{1+d_{4}}\pi^{1-d_{4}}d_{3}k^{-d_{4}},$ $\displaystyle D_{5}$
$\displaystyle=$ $\displaystyle\frac{\Lambda}{6k}.$

For Eq. (49), there is no analytical solution. But if assuming that
$u(a(t))=\dot{a}(t)^{2}$, this can be reduced to

$\displaystyle
u^{{}^{\prime}}(a(t))-2D_{1}\frac{u(a(t))}{a(t)}+2D_{2}a(t)+2D_{3}ka(t)^{-1}$
$\displaystyle-$ $\displaystyle
2D_{4}\left[1+\left(\frac{u(a(t))}{k}-\frac{\Lambda}{3k}\right)\right]^{d_{4}}a(t)^{-1}-9d_{2}u(a(t))^{2}a(t)^{-3}-9d_{2}ku(a(t))a(t)^{-3}$
$\displaystyle=$ $\displaystyle 0,\hskip 22.76219pt$ (50)

where $u^{{}^{\prime}}(a(t))=du(a(t))/da(t)$. With binomial expansion, the
squared bracket can be approximated to unity. Such an assumption333This
reduces the certainty of the proposed solution but this seems to be the only
approximation possible! The inclusion of higher terms simply prevents any
analytical solution. leads to an analytical solution,

$\displaystyle u(a(t))$ $\displaystyle=$ $\displaystyle\dot{a}(t)^{2}$ (51)
$\displaystyle=$ $\displaystyle-\frac{3^{\frac{{\cal
A}}{k}}(d_{2}k)^{\frac{{\cal A}}{2k}}}{{\cal
K}_{1}(t)}\left\\{\frac{18d_{2}k^{2}{\cal
M}\left[\frac{5D_{4}-2k(D_{1}+2D_{3}-3)+{\cal A}}{4k},2+\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]}{[5D_{4}-2k(D_{1}+2D_{3}-1)+{\cal
A}]^{-1}(2k+{\cal A})}\right\\}a(t)^{-\frac{{\cal A}}{k}}$ $\displaystyle+$
$\displaystyle\frac{2\,3^{\frac{{\cal A}}{k}}(d_{2}k)^{\frac{{\cal
A}}{2k}}}{{\cal K}_{1}(t)}\left\\{(D_{4}+2k-2kD_{1}+{\cal A}){\cal
M}\left[\frac{5D_{4}-2k(D_{1}+2D_{3}-1)+{\cal A}}{4k},1+\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]\right\\}a(t)^{2-\frac{{\cal A}}{k}}$
$\displaystyle-$ $\displaystyle\frac{c_{4}2^{\frac{{\cal A}}{2k}}}{{\cal
K}_{1}(t)}\left\\{\frac{18d_{2}k^{2}{\cal
M}\left[-\frac{-5D_{4}+2k(D_{1}+2D_{3}-3)+{\cal A}}{4k},2-\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]}{[-5D_{4}+2k(D_{1}+2D_{3}-1)+{\cal
A}]^{-1}({\cal A}-2k)}\right\\}$ $\displaystyle-$
$\displaystyle\frac{2c_{4}2^{\frac{{\cal A}}{2k}}}{{\cal
K}_{1}(t)}[2k(D_{1}-1)-D_{4}+{\cal A}]{\cal
M}\left[-\frac{-5D_{4}+2k(D_{1}+2D_{3}-1)+{\cal A}}{4k},1-\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]a(t)^{2},$

where

$\displaystyle{\cal K}_{1}(t)$ $\displaystyle=$ $\displaystyle
36d_{2}k\left\\{3^{\frac{{\cal A}}{k}}(d_{2}k)^{\frac{{\cal A}}{2k}}{\cal
M}\left[\frac{1}{4k}\left(5D_{4}-2k(D_{1}+2D_{3}-1)+{\cal
A}\right),1+\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]a(t)^{-\frac{{\cal A}}{k}}\right.$
$\displaystyle\left.\hskip 36.98857pt+c_{4}2^{\frac{{\cal A}}{2k}}{\cal
M}\left[\frac{-1}{4k}\left(-5D_{4}+2k(D_{1}+2D_{3}-1)+{\cal
A}\right),1-\frac{{\cal A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]\right\\},$
$\displaystyle{\cal A}$ $\displaystyle=$
$\displaystyle\left[72d_{2}(D_{2}+D_{4}D_{5})k^{2}+[D_{4}-2(D_{1}-1)k]^{2}\right]^{1/2}$
(52)

Even for the resultant differential equation (51) there is no analytical
solution in terms of the cosmic time $t$. The only possible solution is
$\dot{a}(t)$ as a function of $a(t)$. Such a solution [impeded in Eq. (53)]
leads to an analytical expression for $H(t)$ as a function of $a(t)$, i.e.
functionality, which in turn depends on regularized confluent hypergeometric
function whose asymptotic limit reads ${\cal
M}(\\{a,b,z\\})\sim\Gamma(b)(e^{z}z^{a-b}+(-)^{-a}/\Gamma(b-a))$ Abramowitz
and Stegun (1964). Two of the three regular singularities of ${\cal
M}(\\{a,b,z\\})$ are conjectured to merge into an irregular singularity and
therefrom the conjugate "confluent" emerges. The Hubble parameter reads

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\left\\{-\frac{3^{\frac{{\cal A}}{k}}(d_{2}k)^{\frac{{\cal
A}}{2k}}}{{\cal K}_{1}(t)}\left[\frac{18d_{2}k^{2}{\cal
M}\left[\frac{5D_{4}-2k(D_{1}+2D_{3}-3)+{\cal A}}{4k},2+\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]}{[5D_{4}-2k(D_{1}+2D_{3}-1)+{\cal
A}]^{-1}(2k+{\cal A})}\right]a(t)^{-2-\frac{{\cal A}}{k}}\right.$ (53)
$\displaystyle+$ $\displaystyle\left.\frac{2\,3^{\frac{{\cal
A}}{k}}(d_{2}k)^{\frac{{\cal A}}{2k}}}{{\cal
K}_{1}(t)}\left[(D_{4}+2k-2kD_{1}+{\cal A}){\cal
M}\left[\frac{5D_{4}-2k(D_{1}+2D_{3}-1)+{\cal A}}{4k},1+\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]\right]a(t)^{-\frac{{\cal
A}}{k}}\right.$ $\displaystyle-$ $\displaystyle\left.\frac{c_{4}2^{\frac{{\cal
A}}{2k}}}{{\cal K}_{1}(t)}\left[\frac{18d_{2}k^{2}{\cal
M}\left[-\frac{-5D_{4}+2k(D_{1}+2D_{3}-3)+{\cal A}}{4k},2-\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]}{[-5D_{4}+2k(D_{1}+2D_{3}-1)+{\cal
A}]^{-1}({\cal A}-2k)}\right]a(t)^{-2}\right.$ $\displaystyle-$
$\displaystyle\left.\frac{2c_{4}2^{\frac{{\cal A}}{2k}}}{{\cal
K}_{1}(t)}[2k(D_{1}-1)-D_{4}+{\cal A}]{\cal
M}\left[-\frac{-5D_{4}+2k(D_{1}+2D_{3}-1)+{\cal A}}{4k},1-\frac{{\cal
A}}{2k},-\frac{9kd_{2}}{2a(t)^{2}}\right]\right\\}^{1/2}.$

It is worth highlighting that Kummer confluent hypergeometric functions, for
instance, which are common standard forms of the confluent hypergeometric
functions ${\cal M}$, have a regular singular point, at
$z\equiv-9kd_{2}/2a(t)^{2}=0$ and an irregular singular point at
$z\equiv-9kd_{2}/2a(t)^{2}=\infty$. Thus, the curvature parameter $k$ and the
scale factor $a(t)$ define whether regular or irregular singular point
appears. At vanishing and finite cosmological constant, the results of $H(t)$
vs. $a(t)$ are shown in Fig. 5.

####  4.2.2 QCD-EW era

When substituting with the barotropic equation for the pressure, Eq. (9), and
the bulk viscosity, Eq. (35), in Eq. (48), we get

$\displaystyle\ddot{a}(t)a(t)-E_{1}\dot{a}(t)^{2}+E_{2}a(t)^{2}+E_{3}k-E_{4}\left[1+\left(\frac{\dot{a}(t)^{2}}{k}-\frac{\Lambda}{3k}\right)\right]^{e_{3}}\frac{\dot{a}(t)^{2}}{a(t)^{2e_{3}}}+$
$\displaystyle
E_{5}\left[1+\left(\frac{\dot{a}(t)^{2}}{k}-\frac{\Lambda}{3k}\right)\right]^{\delta_{2}}a(t)^{2-2e_{3}}=0,$
(54)

with the coefficients

$\displaystyle E_{1}$ $\displaystyle=$ $\displaystyle 1+12\pi
e_{1}-\frac{3}{2}(1+\beta_{2}),$ $\displaystyle E_{2}$ $\displaystyle=$
$\displaystyle 4\pi\alpha_{2}-\frac{\Lambda}{2}(1+\beta_{2}),$ $\displaystyle
E_{3}$ $\displaystyle=$ $\displaystyle\frac{3}{2}(1+\beta_{2})-1,$
$\displaystyle E_{4}$ $\displaystyle=$ $\displaystyle
2^{2-3e_{3}}3^{1+e_{3}}\pi^{1-e_{3}}e_{2}k^{-e_{3}},$ $\displaystyle E_{5}$
$\displaystyle=$ $\displaystyle
2^{2-3\delta_{2}}3^{\delta_{2}}\pi^{1-\delta_{2}}\gamma_{2}k^{-\delta_{2}}.$

Assuming that $u(a(t))=\dot{a}^{2}(t)$ and applying the same approximation
given in Eq. (23), the previous differential equation can be reduced to

$\displaystyle
u^{\prime}(a(t))-2E_{1}\frac{u(a(t))}{a(t)}+2E_{2}a(t)-2k\frac{E_{3}}{a(t)}$
$\displaystyle+$ $\displaystyle
2E_{4}\left[1-C_{4}\frac{a(t)^{2}}{k}-\frac{u(a(t))}{3k}\right]\frac{u(a(t))}{a(t)^{1+2e_{3}}}+2E_{5}\left[1-C_{4}\frac{a(t)^{2}}{k}-\frac{u(a(t))}{3k}\right]a(t)^{1-3\delta_{3}}$
$\displaystyle=$ $\displaystyle 0.\hskip 22.76219pt$ (55)

An analytical solution is only possible when both squared brackets are
replaced by unity

$\displaystyle u(a(t))$ $\displaystyle=$ $\displaystyle\dot{a}(t)^{2}$ (56)
$\displaystyle=$
$\displaystyle\left\\{\left[c_{4}\,e_{3}\,a(t)^{2E_{1}}-E_{3}\,k\;{\cal
L}_{1-\frac{E_{1}}{e_{3}}}\left(-\frac{E_{4}}{e_{3}}a(t)^{-2e_{3}}\right)\right.\right.$
$\displaystyle-$ $\displaystyle\left.\left.E_{2}\;{\cal
L}_{\frac{1-E_{1}+e_{3}}{e_{3}}}\left(-\frac{E_{4}}{e_{3}}a(t)^{2}\right)a(t)^{2}\right]a(t)^{2\delta_{2}}\right.$
$\displaystyle-$ $\displaystyle\left.E_{5}\;{\cal
L}_{1-\frac{E_{1}-1+\frac{3}{2}\delta_{2}}{e_{3}}}\left(-\frac{E_{4}}{e_{3}}a(t)^{-2e_{3}}\right)a(t)^{2}\right\\}\frac{e^{-\frac{E_{4}}{e_{3}}a(t)^{-2e_{3}}}}{e_{3}}a(t)^{-2\delta_{2}},$

where

$\displaystyle{\cal L}_{\nu}(z)$ $\displaystyle=$
$\displaystyle\int_{1}^{\infty}\frac{e^{-zt}}{t^{\nu}}dt.$ (57)

Thus, the corresponding Hubble parameter reads

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{1}{a(t)}\left\\{\left[\left[c_{4}e_{3}a(t)^{2E_{1}}-E_{3}k\;{\cal
L}_{1-\frac{E_{1}}{e_{3}}}\left(-\frac{E_{4}}{e_{3}}a(t)^{-2e_{3}}\right)\right.\right.\right.$
(58) $\displaystyle-$ $\displaystyle\left.\left.\left.E_{2}\;{\cal
L}_{\frac{1-E_{1}+e_{3}}{e_{3}}}\left(-\frac{E_{4}}{e_{3}}a(t)^{2}\right)a(t)^{2}\right]a(t)^{2\delta_{2}}\right.\right.$
$\displaystyle-$ $\displaystyle\left.\left.E_{5}\;{\cal
L}_{1-\frac{E_{1}-1+\delta_{2}}{e_{3}}}\left(-\frac{E_{4}}{e_{3}}a(t)^{-2e_{3}}\right)a(t)^{2}\right]\frac{e^{-\frac{E_{4}}{e_{3}}a(t)^{-2e_{3}}}}{e_{3}}\;a(t)^{-2\delta_{2}}\right\\}^{1/2}.$

The results of $H(t)$ vs. $a(t)$ are depicted in Fig. 5.

####  4.2.3 EW (asymptotic) era

For pressure, Eq. (10), and bulk viscosity, Eq. (36), Eq. (48) leads to

$\displaystyle\ddot{a}(t)a(t)-F_{1}\dot{a}^{2}(t)-F_{2}a^{2}(t)+F_{3}k-F_{4}\left[1\left(\frac{\dot{a}^{2}(t)}{k}-\frac{\Lambda}{3k}a^{2}(t)\right)\right]^{f_{3}}a^{-2f_{3}}(t)\,\dot{a}^{2}(t)$
$\displaystyle=$ $\displaystyle 0,\hskip 22.76219pt$ (59)

where

$\displaystyle F_{1}$ $\displaystyle=$ $\displaystyle 1+12\pi
f_{1}-\frac{3}{2}(1+\gamma_{3}),$ $\displaystyle F_{2}$ $\displaystyle=$
$\displaystyle\frac{\Lambda}{2}(1+\gamma_{3}),$ $\displaystyle F_{3}$
$\displaystyle=$ $\displaystyle\frac{3}{2}(1+\gamma_{3})-1,$ $\displaystyle
F_{4}$ $\displaystyle=$ $\displaystyle
2^{2-3f_{3}}3^{1+f_{3}}\pi^{1-f_{3}}f_{2}k^{-f_{3}}.$

Applying the same substitution, $u(a(t))=\dot{a}(t)^{2}$, and taking into
account the first term of the binomial expansion as unity, we get an
analytical solution, functionality,

$\displaystyle u(a(t))$ $\displaystyle=$ $\displaystyle\dot{a}(t)^{2}$ (60)
$\displaystyle=$
$\displaystyle\frac{e^{-\frac{F_{4}}{f_{3}}a(t)^{-2f_{3}}}}{f_{3}}\left\\{c_{4}f_{3}a(t)^{2F_{1}}-F_{3}k{\cal
L}_{1-\frac{F_{1}}{f_{3}}}\left(-\frac{F_{4}}{f_{3}}a(t)^{-2f_{3}}\right)\right.$
$\displaystyle+$ $\displaystyle\left.F_{2}{\cal
L}_{\frac{1-F_{1}+f_{3}}{f_{3}}}\left(-\frac{F_{4}}{f_{3}}a(t)^{-2f_{3}}\right)a(t)^{2}\right\\}.\hskip
19.91692pt$

Accordingly, the Hubble parameters is given as

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{1}{a(t)}\left[\frac{e^{-\frac{F_{4}}{f_{3}}a(t)^{-2f_{3}}}}{f_{3}}\left\\{c_{4}f_{3}a(t)^{2F_{1}}-F_{3}k{\cal
L}_{1-\frac{F_{1}}{f_{3}}}\left(-\frac{F_{4}}{f_{3}}a(t)^{-2f_{3}}\right)\right.\right.$
(61) $\displaystyle+$ $\displaystyle\left.\left.F_{2}{\cal
L}_{\frac{1-F_{1}+f_{3}}{f_{3}}}\left(-\frac{F_{4}}{f_{3}}a(t)^{-2f_{3}}\right)a(t)^{2}\right\\}\right]^{1/2}.\hskip
19.91692pt$

The dependence of $H(t)$ on $a(t)$ is presented in Fig. 5.

###  4.3 Israel-Stewart relativistic viscous fluid

In order to fix the acausality and instability problem of Eckart theory,
Israel and Stewart have introduced a relativistic second-order theory for
relativistic fluid Israel (1976); Israel and Stewart (1976). With extended
irreversible thermodynamics, this theory was then developed by Hiscock and
Lindblom Hiscock and Lindblom (1983). This theory is also characterized by a
deviation from equilibrium as defined by Eckart theory. Quantities such as
bulk stress, heat flow, and shear stress are treated as independent dynamical
variables. Accordingly, $14$ dynamical fluid variables have to be estimated.
The role that this type of causal thermodynamics would play in the general
theory of relativity was reported in ref. Maartens (1995). A general algebraic
form for $S^{\alpha}$ including a second-order term in the dissipative
thermodynamical flux $\Pi$ Israel (1976); Israel and Stewart (1976);
Chattopadhyay (2016) reads

$S^{\alpha}=s\,n\,u^{\alpha}+\frac{\tau}{\zeta}\,\Pi^{2}\,\frac{u^{\alpha}}{2T},$
(62)

where $\tau$ is the relaxation time. Similar to Eckart theory, the
corresponding number flux could be given as

$\displaystyle N^{\alpha}$ $\displaystyle=$ $\displaystyle N\,u^{\alpha}.$
(63)

For the evolution of the bulk viscous pressure, we adopt the causal evolution
equation in the simplest way, i.e. linear in $\Pi$ satisfying the $H$-theorem
Maartens (1995). Accordingly, the entropy production remains nonnegative,
$S_{;i}^{i}=\Pi^{2}/\zeta T\geq 0$ Israel (1976); Israel and Stewart (1976)).
The causal transport equation of the bulk viscous pressure reads Maartens
(1995)

$\tau\dot{\Pi}+\Pi=-3\,\zeta\,H-\frac{\epsilon}{2}\
\tau\,\Pi\,\left(3\,H+\frac{\dot{\tau}}{\tau}-\frac{\dot{\zeta}}{\zeta}-\frac{\dot{T}}{T}\right),$
(64)

where $\epsilon$ is a parameter controlling the type of considered theory.
$\epsilon=1$ assures full theory, while $\epsilon=0$ a truncated one. It is
obvious that the non-causal Eckart theory can be retrieved, Eq. (46), at
$\tau=0$. In order to have a closed system from Eqs. (4), (7) and (64), we
have to introduce EoS for the pressure $p(t)$, the temperature $T(t)$, bulk
viscosity coefficient $\zeta$(t), and the relaxation time $\tau(t)$,
respectively, section 4.1.

In the sections that follow, we elaborate the consequences of the various
barotropic EoS for $p(t)$, $T(t)$, $\zeta(t)$, and $\tau(t)$ in strong of
electroweak epochs of the early Universe. We get sophisticated differential
equations. We hope that this concrete mathematical problem finds resonances
among mathematicians. Despite their apparent challenging analytical solutions,
we separately derive them in the appendices. A future work shall be devoted in
order to propose numerical solutions for all these differential equations.

##  5 Results

The present section summarizes the results of the possible analytical
solutions outlined in sections 3 and 4.2. They are only limited to the Hubble
parameter in dependence on the scale factor for non-viscous, section 3, and
Eckart-type viscous cosmic backgrounds, section 4.2. As introduced, from the
corresponding EoS, we could differentiate between the various epochs of the
early Universe. Nevertheless, we did not emphasize when each epoch starts or
when ends, i.e. in terms of the cosmic time. Such a concrete limitation
becomes only possible when the initial and the final conditions are precisely
determined. This is not precisely available. As alternatives, we would be able
to propose for each epoch an interval of cosmic energy densities, which in
turn could be related to an interval of the Hubble parameter. The latter can
hen be given as functions of the scale factor; the proposed analytical
solutions. On the other hand, such a concrete limitation would be only
urgently needed, when a complete or an inter-epochal picture is to be drawn.
The results discussed in the sections that follow are not limiting the
evolution of the Hubble parameter within the successive epochs. They cover a
wider range than than of the corresponding epoch. Accordingly, we conclude
that the evolution during the successive epochs characterized by electroweak
and strong interactions would not be monotonic.

###  5.1 Non-viscous fluid

Fig. 4: The dependence of the Hubble parameter on the scale factor in non-
viscous cosmic background is depicted for finite (top) and vanishing
cosmological constant (bottom panel). The equations of state for hadron, QCD-
EW and asymptotic limit are presented as dashed, dotted, long dashed curves,
respectively.

Figure 4 depicts the dependence of the Hubble parameter on the scale factor at
finite (top panel) and vanishing cosmological constant (top panel). The
results for the equations of state characterizing hadron, QCD-EW, and
asymptotic limit are presented as dashed, dotted, long dashed curves,
respectively. We also draw the ideal gas results as tiny dashed curves. There
is a rapid decrease in $H(t)$ with increasing $a(t)$. The various epochs (the
different EoS) show miscellaneous rates. Relative to the ideal gas EoS, hadron
and asymptotic EoS look very similar, especially, at finite cosmological
constant (top panel). At large $a(t)$, both hadron and asymptotic EoS become
almost identical. The QCD/EW EoS shows a slightly different behavior,
especially at large $a(t)$, where $H(t)$ diminishes.

At vanishing cosmological constant (bottom panel), the rate of decreasing
$H(t)$ with the increase in $a(t)$ is larger than the one observed in the top
panel. Here, only QCD/EW EoS looks similar to the ideal gas EoS, while both
hadron and asymptotic epochs look almost identical. At small $a(t)$, both have
a similar decrease as the one of ideal and QCD/EW EoS, while, at large $a(t)$,
their corresponding $H(t)$ vanishes.

As discussed, each EoS should be restrictively utilized within a specific
interval of the cosmic time characterizing the corresponding cosmic epoch. Due
to the mathematical difficulties associated with the resulting differential
equations so that the proposed analytical solutions are restricted to the
functionality $H(a(t))$ but not in terms of the cosmic time $t$, directly, we
are left with a unique alternative. This is relating the various epochs to an
interval of energy densities, as introduced in ref. Tawfik and Greiner (2019);
Tawfik and Mishustin (2019) and section 4.1. It is obvious that even this
option is an approximation. Thus, we leave the results drawn in Fig. 4
unchanged. The conclusion which could be drawn here is that the cosmic
evolution [$H(t)$ vs. $a(t)$] seems not monotonic, especially over the
successive asymptotic, EW-QCD and hadron epochs.

###  5.2 Eckart-type viscous fluid

Fig. 5: The same as in Fig. 4 but here for viscousv cosmic geometry (Eckart
theory).

Figure 5 presents the same as in in Fig. 4 but here for viscous cosmic
geometry (Eckart theory). Comparing with the results depicted in Fig. 4, the
dependence of $H(t)$ vs. $a(t)$ in viscous background geometry looks very
different. Hadronic EoS is associated with non-singularity. Finite
cosmological constant likely assures non-singularity in both quantities, while
vanishing cosmological constant is associated with non-singular Hubble
parameter. The QCD-EW EoS results in diverging Hubble parameter within a tiny
range of the scale parameter. At lower $a(t)$, we find that $H(t)$ remains
almost vanishing. For the asymptotic EoS, at finite cosmological constant and
low $a(t)$, $H(t)$ vanishes. Then, $H(t)$ gets positive small values. At
higher $a(t)$, $H(t)$ becomes non-physical. Again, within the short range of
$a(t)$, $H(t)$ diverges. For the asymptotic EoS, at vanishing cosmological
constant and low $a(t)$, $H(t)$ vanishes. Then increasing $a(t)$, the
resulting $H(t)$ very slightly linearly decreases. Then, $H(t)$ diverges
within the short range of $a(t)$.

##  6 Conclusions

Based on recent progress achieved, especially the numerical and experimental
studies on hadron, parton, and EW matter, the main conclusion of the present
study is that the analytical solutions for EoS, in which as much as possible
contributions from both standard model for elementary particles and standard
model for cosmology are taken into consideration, are sophisticated. The only
possible analytical solutions are the ones relating the Hubble parameter to
the scale factor, functionality, in non-viscous and Eckart-type-viscous cosmic
backgrounds. For Israel-Stewart-viscosity, the resulting differential
equations are challenging tasks for mathematicians. We have outlined these
differential equations as road-maps for future studies.

Recent non-perturbative and perturbative calculations with as much as possible
quark flavors at almost physical masses have been combined with the thermal
contributions from photons, charged neutrinos, leptons, electroweak particles
($W^{\pm}$ and $Z^{0}$ bosons), and the scalar Higgs bosons. Various
thermodynamic quantities, including pressure, energy density, bulk viscosity,
relaxation time, and temperature for almost net-baryon-free cosmic matter
could be calculated up to the TeV-scale, i.e. covering hadron, QGP and
electroweak (EW) phases.

In equivalence with Newtonian mechanics and based on Friedman solutions and
the conservation of the energy-momentum tensor, McCrea and Milde and by
Peebles derived the temporal evolution of the energy density, i.e. an equation
of motion, with vanishing and finite pressure, that dictates that the decrease
in the energy content of the Universe is given by the energy budget due to the
expansion and the work done by the pressure. We have followed the same
procedure and in order to have a closed set of equations, we have integrated
with various equations of state, such as pressure vs. energy density. For the
present study, we have introduced a reliable estimation for the bulk pressure,
for which we have taken into consideration Eckart (first order) and Israel-
Stewart (second order) theories for relativistic fluid. For the latter, we
found that the resulting differential equations are higher-ordered nonlinear
nonhomogeneous so that no analytical solution could be proposed, so far. For
the earlier, the only possible solutions relates the Hubble parameter with the
scale factor, but none of them could be directly given in terms of the cosmic
time.

The present study has a potential to be extended to cover new standard
inflationary cosmology with baryosynthesis and dark matter, for which reliable
barotropic EoS are unfortunately missing. Taking into consideration the
possible influence processes of the beyond standard model is also conditioned
to reliable barotropic EoS. Despite the observational constraints on the
cosmological evolution at earlier stages are still challenging, another
extension to cover light element abundance with BBN predictions could subject
to a future study. We would like to suggest concrete predictions and/or
observable features of the effects of bulk viscosity in the early cosmological
evolution. A framework of new standard cosmology would be rather the ultimate
goal. The present script is designed to pave a path towards these goals.

###### Acknowledgements.

The work of AT was supported by the ExtreMe Matter Institute (EMMI) at the GSI
Helmholtz Centre for Heavy Ion Research, Visiting Professor 2019. The authors
declare no conflict of interest. yea

##  Appendix A Relativistic viscous fluid in the expanding early Universe

###  A.1 Israel-Stewart second-order theory

####  A.1.1 Hadron epoch

When starting with the continuity equation, Eq. (7), which can be rewritten as

$\displaystyle\Pi(t)$ $\displaystyle=$
$\displaystyle-\frac{\dot{\rho}(t)}{3H(t)}-\left[\rho(t)+p(t)\right],$ (65)

and substituting with the corresponding EoS, Eq. (8), where $\rho(t)$ can be
replaced as in Eq. (4), we obtain an expression for the viscous stress tensor,

$\displaystyle\Pi(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\dot{H}(t)-\left(1+\beta_{1}\right)\frac{3}{8\pi}H^{2}(t)$
(66)
$\displaystyle+\left[1-\frac{3}{2}(1+\beta_{1})\right]\frac{k}{4\pi}a^{-2}(t)+\left(1+\beta_{1}\right)\frac{\Lambda}{8\pi}-\alpha_{1},$

which is valid for Eckart as well as for Israel-Stewart theories. For the
latter, we take into consideration the second-order entropy fulfilling
causality and stability conditions. Then, the time derivative of the viscous
stress tensor reads

$\displaystyle\dot{\Pi}(t)$ $\displaystyle=$
$\displaystyle-\left[1-\frac{3}{2}(1+\beta_{1})\right]\frac{k}{2\pi}a^{-2}(t)H(t)$
(67)
$\displaystyle-(1+\beta_{1})\frac{3}{4\pi}H(t)\dot{H}(t)-\frac{1}{4\pi}\ddot{H}(t).$

Having both expressions, Eqs. (66) and (67), we still need additional EoS for
$\tau(t)$, $\zeta(t)$ and $T(t)$ and their temporal evolutions to solve the
differential equation resulting from Eq. (64). As $\tau(t)$, $\zeta(t)$ and
$T(t)$ are expressed in dependence on $\rho(t)$, given Eq. (7), their temporal
evolutions shall be depending on $\dot{\rho}(t)$, which in turn can be
expressed in dependence on the scale factor $a(t)$. Then Eq. (64) leads to a
sophisticated third-order inhomogeneous differential equation

$\displaystyle
3\left[d_{1}-\frac{d_{2}\Lambda}{8\pi}+\frac{3d_{2}}{8\pi}\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}+8^{-d_{3}}d_{3}\left(-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right)^{d_{3}}\right]\frac{\dot{a}(t)}{a(t)}+$
$\displaystyle\frac{3\epsilon}{16\pi}\left[\left(-8\pi\alpha_{1}+\Lambda(1+\beta_{1})\right)a(t)^{2}-\frac{k+\dot{a}(t)^{2}}{(1+3\beta_{1})^{-1}}-2a(t)\ddot{a}(t)\right]\left[g_{1}+g_{2}e^{-8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}\right]\frac{\dot{a}(t)}{a(t)^{3}}$
$\displaystyle\left\\{1-\frac{k+\dot{a}(t)-a(t)\ddot{a}(t)}{\Lambda
a(t)^{2}-3(k+\dot{a}(t)^{2})}\left[\frac{2^{1-3g_{4}}g_{2}g_{3}g_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}{g_{1}+g_{2}e^{-8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}}-\frac{2\beta_{4}\gamma_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{\gamma_{4}}}{8^{\gamma_{4}}\alpha_{4}+\beta_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{\gamma_{4}}}-\right.\right.$
$\displaystyle\left.\left.\frac{2\left[3\times
8^{d_{4}}d_{2}(k+\dot{a}(t)^{2})+a(t)^{2}\left[-8^{d_{4}}d_{2}\Lambda+8\pi
d_{3}d_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{d_{4}}\right]\right]}{3\times
8^{d_{4}}d_{2}(k+\dot{a}(t)^{2})+a(t)^{2}\left[8^{1+d_{4}}\pi
d_{1}-8^{d_{4}}d_{2}\Lambda+8\pi
d_{3}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{d_{4}}\right]}\right]\right\\}+$
$\displaystyle\frac{e^{-8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}}{8\pi
a(t)^{3}}\left\\{\left[-8\pi\alpha_{1}+\Lambda(1+\beta_{1})\right]a(t)^{3}e^{8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}+\right.$
$\displaystyle\left.2\left[k+\dot{a}(t)^{2}\right](1+3\beta_{1})\left[g_{2}+g_{1}e^{8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}\right]\dot{a}(t)-\right.$
$\displaystyle\left.\left[e^{8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}\frac{k+\dot{a}(t)^{2}}{(1+3\beta_{1})^{-1}}+6\beta_{1}\left[g_{2}+g_{1}e^{8^{-g_{4}}g_{3}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}\right]\dot{a}(t)\ddot{a}(t)\right]a(t)-\right.$
$\displaystyle\left.2\left[e^{-8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}\ddot{a}(t)+\left[g_{2}+g_{1}e^{8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{g_{4}}}\right]\dddot{a}(t)\right]a(t)^{2}\right\\}=0.\hskip
19.91692pt$ (68)

Despite the obvious assessment that there is no analytical solution to be
proposed, we state this expression here and hope that interested
mathematicians become interested in such highly complicated physical problems.

####  A.1.2 QGP epoch

As discussed, for cosmic relativistic fluid, no matter whether Eckart and
Israel-Stewart theories are applied, the viscous stress tensor $\Pi(t)$ can be
deduced for a given EoS, Eq. (9), which in turn can be expressed in terms of
the Hubble parameter, Eq. (4)

$\displaystyle\Pi(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\dot{H}(t)-\left(1+\beta_{2}\right)\frac{3}{8\pi}H^{2}(t)+\left(1-\frac{3}{2}(1+\beta_{2})\right)\frac{k}{4\pi}a^{-2}(t)$
(69) $\displaystyle+$
$\displaystyle\left(1+\beta_{2}\right)\frac{\Lambda}{8\pi}-\alpha_{2}-\gamma_{2}\left[\frac{3}{8\pi}\left(H^{2}(t)+\frac{k}{2a^{2}(t)}-\frac{\Lambda}{3}\right)\right]^{\delta_{2}},$

Then, the time derivative of $\Pi(t)$ is given as

$\displaystyle\dot{\Pi}(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\ddot{H}(t)-\left(1+\beta_{2}\right)\frac{3}{4\pi}H(t)\dot{H}(t)-\left(1-\frac{3}{2}(1+\beta_{2})\right)\frac{k}{2\pi}a^{-2}(t)H(t)$
(70) $\displaystyle-$ $\displaystyle
2\gamma_{2}\delta_{2}\left[\frac{3}{8\pi}\left(H^{2}(t)+\frac{k}{2a^{2}(t)}-\frac{\Lambda}{3}\right)\right]^{\delta_{2}-1}\left(H(t)\dot{H}(t)-ka^{-2}(t)H(t)\right).$

This results in a highly sophisticated third-order inhomogeneous differential
equation

$\displaystyle-\alpha_{2}-8^{-\delta_{2}}\gamma_{2}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{\delta_{2}}+\frac{\Lambda
a(t)^{2}-3[k+\dot{a}(t)^{2}]}{8\pi
a(t)^{2}}(1+\beta_{2})+\frac{k+\dot{a}(t)\ddot{a}(t)}{4\pi a(t)^{2}}+$
$\displaystyle
3\left[d_{1}-\frac{d_{2}\Lambda}{8\pi}+\frac{3d_{2}}{8\pi}\frac{k+\dot{a}(t)^{2}}{a(T)^{2}}+8^{-d_{4}}d_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{d_{4}}\right]\frac{\dot{a}(t)}{a(t)}+$
$\displaystyle\epsilon\frac{3\times 2^{-4-3\delta_{2}}}{\pi
a(t)^{3}}\left[g_{1}+g_{2}e^{-8^{g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{g_{4}}}\right]\left[-8^{\delta_{2}}\frac{k+\dot{a}(t)^{2}}{(1+3\beta_{2})^{-1}}+a(t)^{2}\left[\frac{1+\beta_{2}}{8^{-\delta_{2}}}\Lambda-\right.\right.$
$\displaystyle\left.\left.8\pi\left(8^{\delta_{2}}\alpha_{2}+\gamma_{2}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{\delta_{2}}\right)\right]-2^{1+3\delta_{2}}a(t)\ddot{a}(t)\right]\left\\{1-\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{\Lambda
a(t)^{2}-3\left[k+\dot{a}(t)^{2}\right]}\right.$
$\displaystyle\left.\left[\frac{2^{1-3g_{4}}g_{2}g_{3}g_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{g_{4}}}{g_{2}+g_{1}e^{8^{-g_{4}}g_{3}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{g_{4}}}}-\frac{2\beta_{4}\gamma_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{\gamma_{4}}}{8^{\gamma_{4}}\alpha_{4}+\beta_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{\gamma_{4}}}-\right.\right.$
$\displaystyle\left.\left.\frac{2\left[\frac{3\times
8^{d_{4}}d_{2}}{\left[k+\dot{a}(t)^{2}\right]^{-1}}+a(t)^{2}\left(-8^{d_{4}}d_{2}\Lambda+8d_{3}d_{4}\pi\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}\right)\right]^{d_{4}}\right)\right]}{\left[\frac{3\times
8^{d_{4}}d_{2}}{\left[k+\dot{a}(t)^{2}\right]^{-1}}+a(t)^{2}\left(8^{1+d_{4}}d_{1}\pi-8^{d_{4}}d_{2}\Lambda+8d_{3}\pi\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{d_{4}}\right)\right]}\right]\right\\}\dot{a}(t)+$
$\displaystyle\frac{g_{1}+g_{2}e^{-8^{g_{4}}g_{3}\left[-\frac{\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(T)^{2}}}{\pi}\right]^{g_{4}}}}{4\pi
a(t)^{3}\dot{a}(t)^{2}}\left\\{\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{\left[4\dot{a}(t)^{3}\right]^{-1}}+\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{\left[3\beta_{2}\dot{a}(t)^{3}\right]^{-1}}+\right.$
$\displaystyle\left.\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{\left[2^{\delta_{2}}8^{1-\delta_{2}}\gamma_{2}\delta_{2}\dot{a}(t)^{3}\right]^{-1}}\left[\frac{-\frac{\Lambda}{3}+\frac{k+\dot{a}(t)^{2}}{a(T)^{2}}}{\pi}\right]^{\delta_{2}-1}-\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{\left[a(t)\dot{a}(t)\ddot{a}(t)\right]^{-1}}+\right.$
$\displaystyle\left.\left[3\dot{a}(t)^{4}+\dot{a}(t)^{2}\left[3k-5a(t)\ddot{a}(t)\right]+\frac{a(t)\ddot{a}(t)-k}{\left[a(t)\ddot{a}(t)\right]^{-1}}+a(t)^{2}\dot{a}(t)\dddot{a}(t)\right]\dot{a}(t)\right\\}=0,\hskip
14.22636pt$ (71)

for which the analytical solution a very challenging task.

####  A.1.3 QCD-EW epoch

In this era, we assume that the EoS, Eq. (9), so that an equation very similar
to (69) and (70) shall be obtained. For the time derivative of the bulk
stress, we apply with related barotropic EoS, (37), (36), and (41) for $T(t)$,
$\zeta(t)$ and $\tau(t)$, respectively,

$\displaystyle-\alpha_{2}+\frac{\Lambda a(t)^{2}-3(k+\dot{a}(t)^{2})}{8\pi
a(t)^{2}}(1+\beta_{2})-8^{\delta_{2}}\gamma_{2}\left[-\frac{1}{\pi}\left(\Lambda-\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{\delta_{2}}+\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{4\pi
a(t)^{2}}+$ $\displaystyle
3\left[e_{1}+8^{-e_{3}}e_{2}\left(-\frac{1}{\pi}\left[\Lambda-\frac{3(k+\dot{a}(t)^{2})}{a(t)}\right]\right)^{e_{3}}\right]\frac{\dot{a}(t)}{a(t)}+\epsilon\frac{3\times
2^{-4-3\delta_{2}}}{\pi
a(t)^{3}}\left[h_{1}+\frac{h_{2}}{h_{3}+\log\left(\frac{h_{4}\left[-\Lambda
a(t)^{2}+3(k+\dot{a}(t)^{2}\right]}{8\pi a(t)^{2}}\right)}\right]$
$\displaystyle\left[-8^{\delta_{2}}\frac{k+\dot{a}(t)^{2}}{(1+3\beta_{2})^{-1}}+a(t)^{2}\left[8^{\delta_{2}}(1+\beta_{2})\Lambda-8\pi\left(8^{\delta_{2}}\alpha_{2}+\gamma_{2}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2})}{a(t)^{2}}\right)\right]^{\delta_{2}}\right)\right]-2^{1+3\delta_{2}}a(t)\ddot{a}(t)\right]$
$\displaystyle\left\\{1-\frac{2h_{2}\left[k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)\right]}{\left[h_{3}+\log\left(\frac{h_{4}\left[-\Lambda
a(t)^{2}+3(k+\dot{a}(t)^{2})\right]}{8\pi
a(t)^{2}}\right)\right]\left[h_{2}+h_{1}h_{3}+h_{1}\log\left(\frac{h_{4}\left[-\Lambda
a(t)^{2}+3(k+\dot{a}(t)^{2})\right]}{8\pi a(t)^{2}}\right)\right]\left[\Lambda
a(t)^{2}-3\left[k+\dot{a}(t)^{2}\right]\right]}-\right.$
$\displaystyle\left.\frac{2e_{2}e_{3}\left[k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)\right]\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)}{a(t)^{2}}\right)\right]^{e_{3}}}{\left[\Lambda
a(t)^{2}-3(k+\dot{a}(t)^{2}\right]\left(8^{e_{3}}e_{1}+e_{2}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{e_{3}}\right)}-\right.$
$\displaystyle\left.\frac{2\beta_{4}\gamma_{4}\left[k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)\right]\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)}{a(t)^{2}}\right)\right]^{\gamma_{4}}}{\left[\Lambda
a(t)^{2}-a(k+\dot{a}(t)^{2}\right]\left(8^{\gamma_{4}}\alpha_{4}+\beta_{4}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{\gamma_{4}}\right)}\right\\}\dot{a}(t)-$
$\displaystyle\frac{\left[h_{3}+\log\left(\frac{h_{4}\left[-\Lambda
a(t)^{2}+3(k+\dot{a}(t)^{2}\right]}{8\pi a(t)^{2}}\right)\right]^{-1}}{4\pi
a(t)^{3}\left[\Lambda
a(t)^{2}-3(k+\dot{a}(t)^{2}\right]}\left\\{\left[h_{2}+h_{1}h_{3}+h_{1}\log\left(\frac{h_{4}\left[-\Lambda
a(t)^{2}+3(k+\dot{a}(t)^{2}\right]}{8\pi
a(t)^{2}}\right)\right]\left[\frac{3\dot{a}(t)^{5}}{(1+3\beta_{2})^{-1}}+\right.\right.$
$\displaystyle\left.\left.8^{\delta_{2}}\dot{a}(t)^{3}\left(\frac{3\times
2^{1+3\delta_{2}}}{(1+3\beta_{2})^{-1}}-a(t)^{2}\left[\frac{8^{\delta_{2}\Lambda}}{(1+3\beta_{2})^{-1}}-24\pi\gamma_{2}\delta_{2}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{\delta_{2}}\right]-3^{2}8^{\delta_{2}}\beta_{2}a(t)\ddot{a}(t)\right)+\right.\right.$
$\displaystyle\left.\left.8^{-\delta_{2}}\dot{a}(t)\left(\frac{3\times
8^{\delta_{2}}k^{2}}{(1+3\beta_{2})^{-1}}-ka(t)^{2}\left[\frac{8^{\delta_{1}}\Lambda}{(1+3\beta_{2})^{-1}}-24\pi\gamma_{2}\delta_{1}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{\delta_{2}}\right]-3^{2}8^{\delta_{2}}k\beta_{2}a(t)\ddot{a}(t)+\right.\right.\right.$
$\displaystyle\left.\left.\left.3a(t)^{3}\left[8^{\delta_{2}}\beta_{2}\Lambda-8\pi\gamma_{2}\delta_{2}\left[\frac{1}{\pi}\left(-\Lambda+\frac{3(k+\dot{a}(t)^{2}}{a(t)^{2}}\right)\right]^{\delta_{2}}\right]\ddot{a}(t)\right)+\frac{a(t)^{2}\dddot{a}(t)}{(-3k+\Lambda
a(t)^{2})^{-1}}-3a(t)^{2}\dot{a}(t)^{2}\dddot{a}(t)\right]\right\\}=0.\hskip
14.22636pt$ (72)

####  A.1.4 EW (asymptotic) epoch

Substituting with $\rho(t)$, Eq. (4), and the corresponding EoSe, Eq. (10), in
Eq. (7), the viscous stress tensor and its temporal evolution, respectively,
become

$\displaystyle\Pi(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\dot{H}(t)-\left(1+\gamma_{3}\right)\frac{3}{8\pi}H^{2}(t)$
(73)
$\displaystyle+\left[1-\frac{3}{2}(+\gamma_{3})\right]\frac{3k}{8\pi}a^{-2}(t)+\left(1+\gamma_{3}\right)\frac{\Lambda}{8\pi},$
$\displaystyle\dot{\Pi}(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\ddot{H}(t)-\left(1+\gamma_{3}\right)\frac{3}{4\pi}H(t)\dot{H}(t)$
(74)
$\displaystyle-\left[1-\frac{3}{2}(1+\gamma_{3})\right]\frac{3}{4\pi}a^{-2}(t)H(t).$

Then, the cosmic evolution, Eq. (64), can be expressed as

$\displaystyle
3\left[f_{1}+8^{-f_{3}}f_{2}\left(\frac{1}{\pi}\left[\Lambda-\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{f_{3}}\right]\frac{\dot{a}(t)}{a(t)}+\epsilon\frac{3\times
a^{-4-3k_{2}}k_{1}}{\pi
a(t)^{2}}\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{k_{2}}$
$\displaystyle\log\left[\frac{k_{3}}{8\pi a(t)^{2}}\left[-\Lambda
a(t)^{2}+3\left[k+\dot{a}(t)^{2}\right]\right]\right]\left[\frac{\Lambda
a(t)^{2}}{(1+\gamma_{3})^{-1}}-\frac{k+\dot{a}(t)^{2}}{(1+3\gamma_{3})^{-1}}-2a(t)\ddot{a}(t)\right]\,\dot{a}(t)$
$\displaystyle\left\\{\frac{1}{a(t)}+\frac{2k_{2}\left[k+\dot{a}(t)-a(t)\ddot{a}(t)\right]}{\Lambda
a(t)^{3}-3a(t)\left[k+\dot{a}(t)^{2}\right]}+\frac{2\left[k+\dot{a}(t)-a(t)\ddot{a}(t)\right]}{\Lambda
a(t)^{2}-3\left[a(t)+\dot{a}(t)*2\right]}\frac{1}{a(t)\log\left[\frac{k_{3}\left[-\Lambda
a(t)^{2}+3\left[k+\dot{a}(t)^{2}\right]\right]}{8\pi
a(t)^{2}}\right]}-\right.$
$\displaystyle\left.\frac{k+\dot{a}(t)^{2}-a(t)\ddot{a}(t)}{\Lambda
a(t)^{2}-3\left[k+\dot{a}(t)^{2}\right]}\frac{1}{a(t)}\left[\frac{2f_{2}f_{3}\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{f_{3}}}{8^{f_{3}}f_{1}+f_{2}\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{f_{3}}}-\frac{2\beta_{4}\gamma_{4}\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{\gamma_{4}}}{8^{\gamma_{4}}\alpha_{4}+\beta_{4}\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{\gamma_{4}}}\right]\right\\}+$
$\displaystyle\frac{8^{-1-k_{2}}}{\pi a(t)^{3}}\left\\{\frac{8^{k_{2}}\Lambda
a(t)^{3}}{(1+\gamma_{3})^{-1}}+\frac{2k_{1}\dot{a}(t)}{(1+3\gamma_{3})^{-1}}\log\left[\frac{k_{3}}{8\pi
a(t)^{2}}\left[-\Lambda
a(t)^{2}+3\left[k+\dot{a}(t)^{2}\right]\right]\right]\frac{\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{k_{2}}}{\left[k+\dot{a}(t)^{2}\right]^{-1}}-\right.$
$\displaystyle\left.\left[\frac{8^{k_{2}}k}{(1+3\gamma_{3})^{-1}}+\frac{8^{k_{2}}\dot{a}(t)^{2}}{(1+3\gamma_{3})^{-1}}+6k_{1}\gamma_{3}\log\left[\frac{k_{3}\left[-\Lambda
a(t)^{2}+3\left[k+\dot{a}(t)^{2}\right]\right]}{8\pi
a(t)^{2}}\right]\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{k_{2}}\dot{a}(t)\ddot{a}(t)\right]a(t)-\right.$
$\displaystyle\left.2\left[8^{k_{2}}\ddot{a}(t)+k_{1}\log\left[\frac{k_{3}\left[-\Lambda
a(t)^{2}+3\left[k+\dot{a}(t)^{2}\right]\right]}{8\pi
a(t)^{2}}\right]\left(\frac{1}{\pi}\left[-\Lambda+\frac{k+\dot{a}(t)^{2}}{a(t)^{2}}\right]\right)^{k_{2}}\dddot{a}(t)\right]a(t)^{2}\right\\}=0.\hskip
14.22636pt$ (75)

References

yes

## References

  * Heinz (2001) Heinz, U.W.  The Little bang: Searching for quark gluon matter in relativistic heavy ion collisions.  Nucl. Phys. 2001, A685, 414–431, [arXiv:hep-ph/hep-ph/0009170].  doi:black10.1016/S0375-9474(01)00558-9.
  * Tawfik and Ganssauge (2000) Tawfik, A.M.; Ganssauge, E.  Levy stable law description of the intermittent behavior in Pb + Pb collisions at 158/A-GeV.  Acta Phys. Hung. 2000, A12, 53, [arXiv:hep-ph/hep-ph/0012008].
  * Gyulassy and McLerran (2005) Gyulassy, M.; McLerran, L.  New forms of QCD matter discovered at RHIC.  Nucl. Phys. 2005, A750, 30–63, [arXiv:nucl-th/nucl-th/0405013].  doi:black10.1016/j.nuclphysa.2004.10.034.
  * Heinz et al. (2012) Heinz, U.; Shen, C.; Song, H.  The viscosity of quark-gluon plasma at RHIC and the LHC.  AIP Conf. Proc. 2012, 1441, 766–770, [arXiv:nucl-th/1108.5323].  doi:black10.1063/1.3700674.
  * Adamczyk et al. (2014) Adamczyk, L.; others.  Energy Dependence of Moments of Net-proton Multiplicity Distributions at RHIC.  Phys. Rev. Lett. 2014, 112, 032302, [arXiv:nucl-ex/1309.5681].  doi:black10.1103/PhysRevLett.112.032302.
  * Ryu et al. (2018) Ryu, S.; Paquet, J.F.; Shen, C.; Denicol, G.; Schenke, B.; Jeon, S.; Gale, C.  Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider.  Phys. Rev. 2018, C97, 034910, [arXiv:nucl-th/1704.04216].  doi:black10.1103/PhysRevC.97.034910.
  * Bzdak et al. (2020) Bzdak, A.; Esumi, S.; Koch, V.; Liao, J.; Stephanov, M.; Xu, N.  Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan.  Phys. Rept. 2020, 853, 1–87, [arXiv:nucl-th/1906.00936].  doi:black10.1016/j.physrep.2020.01.005.
  * Komatsu et al. (2011) Komatsu, E.; others.  Seven-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation.  Astrophys. J. Suppl. 2011, 192, 18, [arXiv:astro-ph.CO/1001.4538].  doi:black10.1088/0067-0049/192/2/18.
  * Ade et al. (2016) Ade, P.A.R.; others.  Planck 2015 results. XIII. Cosmological parameters.  Astron. Astrophys. 2016, 594, A13, [arXiv:astro-ph.CO/1502.01589].  doi:black10.1051/0004-6361/201525830.
  * Aghanim et al. (2020) Aghanim, N.; others.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys. 2020, 641, A6, [arXiv:astro-ph.CO/1807.06209].  doi:black10.1051/0004-6361/201833910.
  * Akrami et al. (2020) Akrami, Y.; others.  Planck 2018 results. X. Constraints on inflation.  Astron. Astrophys. 2020, 641, A10, [arXiv:astro-ph.CO/1807.06211].  doi:black10.1051/0004-6361/201833887.
  * Aghanim et al. (2020) Aghanim, N.; others.  Planck 2018 results. V. CMB power spectra and likelihoods.  Astron. Astrophys. 2020, 641, A5, [arXiv:astro-ph.CO/1907.12875].  doi:black10.1051/0004-6361/201936386.
  * Akrami et al. (2020) Akrami, Y.; others.  Planck 2018 results. VII. Isotropy and Statistics of the CMB.  Astron. Astrophys. 2020, 641, A7, [arXiv:astro-ph.CO/1906.02552].  doi:black10.1051/0004-6361/201935201.
  * Misner (1968) Misner, C.W.  The isotropy of the universe.  The Astrophysical Journal 1968, 151, 431.
  * Zeldovich and Novikov (1983) Zeldovich, Ya.B.; Novikov, I.D.  RELATIVISTIC ASTROPHYSICS. VOL. 2. THE STRUCTURE AND EVOLUTION OF THE UNIVERSE; 1983.
  * Tawfik and Mishustin (2019) Tawfik, A.N.; Mishustin, I.  Equation of State for Cosmological Matter at and beyond QCD and Electroweak Eras.  J. Phys. 2019, G46, 125201, [arXiv:hep-ph/1903.00063].  doi:black10.1088/1361-6471/ab46d4.
  * Gron (1990) Gron, O.  Viscous inflationary universe models.  Astrophys. Space Sci. 1990, 173, 191–225.  doi:black10.1007/BF00643930.
  * Maartens (1995) Maartens, R.  Dissipative cosmology.  Class. Quant. Grav. 1995, 12, 1455–1465.  doi:black10.1088/0264-9381/12/6/011.
  * Tawfik et al. (2010) Tawfik, A.; Harko, T.; Mansour, H.; Wahba, M.  Dissipative Processes in the Early Universe: Bulk Viscosity.  Uzbek J. Phys. 2010, 12, 316–321, [arXiv:gr-qc/0911.4105].
  * Tawfik and Wahba (2010) Tawfik, A.; Wahba, M.  Bulk and Shear Viscosity in Hagedorn Fluid.  Annalen Phys. 2010, 522, 849–856, [arXiv:hep-ph/1005.3946].  doi:black10.1002/andp.201000056.
  * Tawfik et al. (2010) Tawfik, A.; Wahba, M.; Mansour, H.; Harko, T.  Hubble Parameter in QCD Universe for finite Bulk Viscosity.  Annalen Phys. 2010, 522, 912–923, [arXiv:gr-qc/1008.0971].  doi:black10.1002/andp.201000103.
  * Adamczyk et al. (2019) Adamczyk, L.; others.  Constraining the initial conditions and temperature dependent viscosity with three-particle correlations in Au+Au collisions.  Phys. Lett. B 2019, 790, 81–88, [arXiv:nucl-ex/1701.06497].  doi:black10.1016/j.physletb.2018.10.075.
  * Tawfik and Harko (2012) Tawfik, A.; Harko, T.  Quark-Hadron Phase Transitions in Viscous Early Universe.  Phys. Rev. 2012, D85, 084032, [arXiv:astro-ph.CO/1108.5697].  doi:black10.1103/PhysRevD.85.084032.
  * Tawfik (2011) Tawfik, A.  The Hubble parameter in the early universe with viscous QCD matter and finite cosmological constant.  Annalen Phys. 2011, 523, 423–434, [arXiv:gr-qc/1102.2626].  doi:black10.1002/andp.201100038.
  * Tawfik and Magdy (2012) Tawfik, A.; Magdy, H.  Thermodynamics of viscous Matter and Radiation in the Early Universe.  Can. J. Phys. 2012, 90, 433–440, [arXiv:gr-qc/1109.6469].  doi:black10.1139/p2012-037.
  * Tawfik et al. (2011) Tawfik, A.; Wahba, M.; Mansour, H.; Harko, T.  Viscous Quark-Gluon Plasma in the Early Universe.  Annalen Phys. 2011, 523, 194–207, [arXiv:gr-qc/1001.2814].  doi:black10.1002/andp.201000052.
  * Tawfik (2010) Tawfik, A.  Thermodynamics in the Viscous Early Universe.  Can. J. Phys. 2010, 88, 825–831, [arXiv:gr-qc/1002.0296].  doi:black10.1139/P10-058.
  * Tawfik and Greiner (2019) Tawfik, A.N.; Greiner, C.  Bulk viscosity at high temperatures and energy densities 2019.  [arXiv:hep-ph/1911.02797].
  * Laine and Schroder (2006) Laine, M.; Schroder, Y.  Quark mass thresholds in QCD thermodynamics.  Phys. Rev. 2006, D73, 085009, [arXiv:hep-ph/hep-ph/0603048].  doi:black10.1103/PhysRevD.73.085009.
  * Laine and Meyer (2015) Laine, M.; Meyer, M.  Standard Model thermodynamics across the electroweak crossover.  JCAP 2015, 1507, 035, [arXiv:hep-ph/1503.04935].  doi:black10.1088/1475-7516/2015/07/035.
  * D’Onofrio and Rummukainen (2016) D’Onofrio, M.; Rummukainen, K.  Standard model cross-over on the lattice.  Phys. Rev. 2016, D93, 025003, [arXiv:hep-ph/1508.07161].  doi:black10.1103/PhysRevD.93.025003.
  * Borsanyi et al. (2016) Borsanyi, S.; others.  Calculation of the axion mass based on high-temperature lattice quantum chromodynamics.  Nature 2016, 539, 69–71, [arXiv:hep-lat/1606.07494].  doi:black10.1038/nature20115.
  * McCrea and Milne (1934) McCrea, W.H.; Milne, E.A.  Newtonian universes and the curvature of space.  The quarterly journal of mathematics 1934, pp. 73–80.
  * McCrea and Milne (2000) McCrea, W.H.; Milne, E.A.  Newtonian Universes and the Curvature of Space.  General Relativity and Gravitation 2000, 32, 1949–1958.  doi:black10.1023/A:1001949117817.
  * Peebles (1993) Peebles, P.J.E.  Principles of Physical Cosmology; 1993.
  * Shalyt-Margolin (2010) Shalyt-Margolin, A.E.  Entropy In The Present And Early Universe, New Small Parameters And Dark Energy Problem.  Entropy 2010, 12, 932–952, [arXiv:gr-qc/0911.5597].  doi:black10.3390/e12040932.
  * Tawfik (2005a) Tawfik, A.  QCD phase diagram: A Comparison of lattice and hadron resonance gas model calculations.  Phys. Rev. 2005, D71, 054502, [arXiv:hep-ph/hep-ph/0412336].  doi:black10.1103/PhysRevD.71.054502.
  * Tawfik (2005b) Tawfik, A.  The Influence of strange quarks on QCD phase diagram and chemical freeze-out: Results from the hadron resonance gas model.  J. Phys. 2005, G31, S1105–S1110, [arXiv:hep-ph/hep-ph/0410329].  doi:black10.1088/0954-3899/31/6/068.
  * Pun et al. (2008) Pun, C.S.J.; Gergely, L.A.; Mak, M.K.; Kovacs, Z.; Szabo, G.M.; Harko, T.  Viscous dissipative Chaplygin gas dominated homogenous and isotropic cosmological models.  Phys. Rev. 2008, D77, 063528, [arXiv:gr-qc/0801.2008].  doi:black10.1103/PhysRevD.77.063528.
  * Eckart (1940) Eckart, C.  The Thermodynamics of Irreversible Processes. 1. The Simple Fluid.  Phys. Rev. 1940, 58, 267–269.  doi:black10.1103/PhysRev.58.267.
  * Landau and Lifshitz (1987) Landau, L.; Lifshitz, E.  Fluid Mechanics; Number v. 6, Elsevier Science, 1987.
  * Osada (2012) Osada, T.  Modification of Eckart theory of relativistic dissipative fluid by introducing extended matching conditions.  Phys. Rev. 2012, C85, 014906, [arXiv:nucl-th/1111.1276].  doi:black10.1103/PhysRevC.85.014906.
  * Piattella et al. (2011) Piattella, O.F.; Fabris, J.C.; Zimdahl, W.  Bulk viscous cosmology with causal transport theory.  JCAP 2011, 1105, 029, [arXiv:astro-ph.CO/1103.1328].  doi:black10.1088/1475-7516/2011/05/029.
  * Israel (1976) Israel, W.  Thermo field dynamics of black holes.  Phys. Lett. 1976, A57, 107–110.  doi:black10.1016/0375-9601(76)90178-X.
  * Coley and van den Hoogen (1995) Coley, A.A.; van den Hoogen, R.J.  Qualitative analysis of viscous fluid cosmological models satisfying the Israel-Stewart theory of irreversible thermodynamics.  Class. Quant. Grav. 1995, 12, 1977–1994, [arXiv:gr-qc/gr-qc/9605061].  doi:black10.1088/0264-9381/12/8/015.
  * Tawfik et al. (2013) Tawfik, A.; Magdy, H.; Ali, A.F.  Effects of quantum gravity on the inflationary parameters and thermodynamics of the early universe.  Gen. Rel. Grav. 2013, 45, 1227–1246, [arXiv:gr-qc/1208.5655].  doi:black10.1007/s10714-013-1522-0.
  * Abramowitz and Stegun (1964) Abramowitz, M.; Stegun, I.A.  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, ninth dover printing, tenth gpo printing ed.; Dover: New York, 1964.
  * Israel (1976) Israel, W.  Nonstationary irreversible thermodynamics: A Causal relativistic theory.  Annals Phys. 1976, 100, 310–331.  doi:black10.1016/0003-4916(76)90064-6.
  * Israel and Stewart (1976) Israel, W.; Stewart, J.M.  Thermodynamics of nonstationary and transient effects in a relativistic gas.  Physics Letters A 1976, 58, 213–215.  doi:black10.1016/0375-9601(76)90075-X.
  * Hiscock and Lindblom (1983) Hiscock, W.A.; Lindblom, L.  Stability and causality in dissipative relativistic fluids.  Annals of Physics 1983, 151, 466–496.  doi:black10.1016/0003-4916(83)90288-9.
  * Chattopadhyay (2016) Chattopadhyay, S.  Israel-Stewart approach to viscous dissipative extended holographic Ricci dark energy dominated universe.  Adv. High Energy Phys. 2016, 2016, 8515967, [arXiv:gr-qc/1604.05297].  doi:black10.1155/2016/8515967.