Cosmic Evolution of Viscous QCD Epoch in Causal Eckart Frame

Abstract

It is conjectured that in cosmological applications the particle current is not modified but finite heat or energy flow. Therefore, comoving Eckart frame is a suitable choice, as it merely ceases the charge and particle diffusion and conserves charges and particles. The cosmic evolution of viscous hadron and parton epochs in casual and non-casual Eckart frame is analyzed. By proposing equations of state deduced from recent lattice QCD simulations including pressure p, energy density $\rho$, and temperature T, the Friedmann equations are solved. We introduce expressions for the temporal evolution of the Hubble parameter $\dot{H}$, the cosmic energy density $\dot{\rho }$, and the share $\dot{\eta }$ and the bulk viscous coefficient $\dot{\zeta }$. We also suggest how the bulk viscous pressure $\mathsf{\Pi }$ could be related to H. We conclude that the relativistic theory of fluids, the Eckart frame, and the finite viscous coefficients play essential roles in the cosmic evolution, especially in the hadron and parton epochs.

1. Introduction

Giving a reliable description for matter and/or radiation filling in the cosmological background geometry, the space, is an essential ingredient to determine the temporal evolution of the Universe. A realistic picture of the cosmic evolution can only be drawn, if reliable equations of state could be proposed. Since no observational evidence of the cosmic evolution, especially during the very early epochs of the Universe, was made so far, we first have to do the best of all in order to best describe the cosmic matter and/or radiation [1].

Recent perturbative and non-perturbative lattice simulations with almost all quark flavors and thermal contributions from charged leptons, electroweak particles, scalar Higgs bosons, and photons have calculated various thermodynamic quantities like pressure, energy density, bulk viscosity, relaxation time, and temperature up to the TeV-scale [2]. Therefore, an access to the relativistic dissipative cosmic fluid covering hadron, QGP, and electroweak (EW) epochs is likely gained. On the other hand, the shear viscosity could be originated to curvature effects through Ricci tensor or Ricci scalar of a given spacetime, for example, some radiation processes from an accretion disk may attribute to shear viscosity and heat flow. The exponential expansion of the Universe significantly contributes to the curvature effects in the bulk viscosity.

An approach for relativistic hydrodynamics has found applications in various disciplines, especially in physics of the early Universe [3,4,5]. For ideal fluids, the theory of hydrodynamics was first formulated by Euler [6]. This was much later extended to non-ideal fluids by Navier [7] and Stokes [8] and utilized in determining dissipative quantities like viscosity and heat flow.

About eight decades ago, the relativistic theory of fluids was introduced by Eckart [9] and Landau–Lifshitz [10], in which the relativistic version of Euler ideal and Navier-Stokes non-ideal fluids was developed. The latter is characterized by first-order gradients of the hydrodynamical variable, i.e., the dissipative quantities are constructed from the first-order gradients of the fluid’s four-velocity, temperature, and chemical potential. But it was ascertained that the dissipative perturbation has infinite speed [11,12,13,14]. A linear relation between the bulk viscous pressure and the expansion in the early Universe was also suggested [5].

Hiscock and Lindblom pointed out that under linear perturbations the relativistic theory for non-ideal fluids suffers from instabilities and allows a non-casual signal propagation [15]. It was also found that the dissipative quantities seem to instantaneously response to the first-order changes in Eckart and Landau–Lifshitz theory. When repairing this in light of the theory of general relativity, Mueller [16] and Israel–Stewart [17,18] came up with a causal theory which is stable under linear perturbations and the dissipative perturbation has finite speed [15]. Here, the second-order gradients are dominant [11,12,13].

To address phenomenology of relativistic heavy-ion collisions assuming vanishing conserved charges in Landau frame, a systematic expansion of the hydrodynamical gradients has been done in context of conformal and non-conformal theories and relativistic viscous fluids could be investigated [19,20]. Under the assumption of conserved charges and finite heat flux, Lahiri started with the second-order gradients expansion in the Eckart frame and developed a causal theory of relativistic non-ideal fluids for various astrophysical applications [14]. In other words, general relativistic causal hydrodynamical equations in the Eckart frame like of bulk viscosity, shear viscosity tensor, and heat flow are obtained up to second-order gradients.

In astrophysical and cosmological applications, where the particle current is not modified but the heat or energy flow is non-negligible, comoving Eckart frame is a suitable choice. As it merely ceases the charge and particle diffusion (conserved charges and particles of the system n are given from $N^{\mu }=n{u}^{\mu }$, where $N^{\mu }$ is the conserved charge/particle flux and $u^{\mu }$ is the four-velocity), while in the Landau frame there is no energy or heat flow ($T^{\mu \nu }{u}_{\nu }=\rho u^{\mu }$, where $\rho$ is the energy density and no dissipation of energy appears), where $u_{\mu }$ stands for four-velocity of particle motion and/or energy transport [21]. To summarize, in comoving frame of reference, the Eckart frame conserves $N^i=o$, while the Landau frame conserves $T^{i0}=0$.

The present work aims at parameterizing various thermodynamical and transport quantities for the general relativistic causal hydrodynamical equations in the Eckart frame. For the first time, a cosmological implication of causal Eckart frame is presented. We focus on recent lattice QCD simulations [22]. The temporal evolution of cosmological parameters like the scale factor, the Hubble parameter, etc. is presented.

The paper is organized as follows. The Friedmann equations are reviewed in Section 2. The equations of state of viscous QCD epoch and the parameterizations for different thermodynamic properties are presented in Section 3. The resulting cosmic evolution is elaborated in Section 5. Section 6 is devoted to the conclusions.

2. Friedmann Equations in Eckart Frame

We assume that the bulk viscous fluid shapes the background geometry, the space, of the early Universe and limit to discussion to the QCD epoch. For spatially flat background geometry, the Friedmann–Lemaiture–Robertson–Walker (FLRW) line element reads [5]

$$\begin{array}{c}\hfill d{s}^2=d{t}^2-a^2\left(t\right)d{r}^2+a^2\left(t\right)\left[r^{2}d{\theta }^2+r^2{sin}^2\theta d{\varphi }^2\right],\end{array}$$

(1)

where $a\left(t\right)$ is the scale factor. The temporal evolution becomes accessible when integrating this in the theory of general relativity. Then, the Hubble parameter $H=\dot{a}\left(t\right)/a\left(t\right)$, where $\dot{a}\left(t\right)$ is the time evolution of the scale factor, can be estimated from the Einstein field Equations [5]

$$\begin{array}{ccc}\hfill H^2& =& \frac{1}{3m_{Pl}^2}\rho ,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \dot{H}+H^2& =& -\frac{1}{6m_{Pl}^2}\left(3p_{eff}+\rho \right),\hfill \end{array}$$

(3)

where $m_{Pl}$ is known as the reduced Planck mass $1/m_{Pl}^2=1.687\times {10}^{-43}$ MeV${}^{-2}$ [5], $p_{eff}$ is the effective thermodynamic pressure $p_{eff}=p+\mathsf{\Pi }$ and $\rho$ is the energy density. Here, we use natural units with $c=ħ=k_B=1$, so that $8\pi G=1/m_{Pl}^2=1.687\times {10}^{-43}{\mathrm{MeV}}^{-2}$. The reduced Planck time is given as $t_{Pl}=1/m_{Pl}=4.0\times 1152\times {10}^{-22}$ MeV${}^{-1}$. Therefore, the Hubble parameter reads

$$\begin{array}{ccc}\hfill H& =& {\left(\frac{1}{3m_{Pl}^2}\rho \right)}^{1/2},\hfill \end{array}$$

(4)

and its time evolution can be deduced from Equations (2) and (3)

$$\begin{array}{ccc}\hfill \dot{H}& =& -\frac{1}{2m_{Pl}^2}\left(p_{eff}+\rho \right).\hfill \end{array}$$

(5)

The temporal evolution of the energy density can be related to H [5]

$$\begin{array}{c}\hfill \dot{\rho }=-3(p_{eff}+\rho )H,\end{array}$$

(6)

To solve the field equations, various baratropic EoS including the bulk viscous pressure $\mathsf{\Pi }$ have to be substituted. Relating $\mathsf{\Pi }$ to $\rho$ or indirectly to H depends the features of causality and stability on the corresponding relativistic fluid theory. As discussed in introduction, the Eckart theory, which is a non-casual and unstable theory, is the relativistic version of non-viscous and viscous fluid theories introduced by Euler and Navier–Stokes, respectively.

In astrophysics and cosmology, the particle current-in contrary to heat or energy flow-likely remains unmodified. This makes the Eckart frame the suitable choice. With a second-order gradients expansion in the Eckart frame, non-causal theory of relativistic non-ideal fluids was developed for astrophysical applications [14]. This is similar to the causal stable Mueller–Israel–Stewart theory but explicitly in the Eckart frame. Various general relativistic causal hydrodynamical equations like bulk viscosity, shear viscosity tensor, and heat flow up to second-order gradients have been derived. The four-velocity normalization conditions reads $u^{\nu }{u}_{\nu }=-1$ and $u_{\nu }{\Delta }^{\nu \mu }=0$, where the projection tensor in curved spacetime is given as ${\Delta }^{\nu \mu }=g^{\nu \mu }+u^{\nu }{u}^{\mu }$ [14]. To summarize, the proposed Eckart and Mueller–Israel–Stewart theories are causal. While the latter considers up to the second-order derivations of the thermodynamic quantities, the earlier is limited to first order, only.

Comments of features of non-casual and causal Eckart Frame are now in order.

• Non-casual Eckart frame: The entropy current can be given as $S^{\nu }=sn{u}^{\nu }$ but it isn’t conserved all the time, where s is the entropy density and n is the number density. It was found that Eckart theory violates the second law of thermodynamics as a result of the divergence of the entropy current $T{S}^{\nu }=-3H\mathsf{\Pi }$ [5]. There is a linear relation between the Hubble parameter H and the bulk viscous pressure $\mathsf{\Pi }$ which results form the first-order derivations of the equilibrium states [5]

  $$\begin{array}{ccc}\hfill \mathsf{\Pi }& =& -3\zeta H,\hfill \end{array}$$

  (7)

  where $\zeta$ is the bulk viscosity coefficient.
• Causal Eckart frame: To avoid violating the causality principle, higher-order deviations from the equilibrium states are taken into consideration, mostly the second-order derivations of the thermodynamic quantities [14]. Using gradient expansion scheme of fluid hydrodynamical variables, the causality preserving property of the theory could be defined. The Mueller–Israel–Stewart theory, which is second-order theory in gradients of fluid hydrodynamical variables gives rise to causal and stable theory of viscous hydrodynamics. By using gradient expansion scheme up to second order and keeping the spirit of Mueller–Israel–Stewart theory, non-causal stable theory of relativistic non-ideal fluids in the Eckart frame was developed [14]. By solving the continuity equation, the modified Friedmann equation and the evolution equation of the bulk viscosity, the effects of chemical potential can be studied. Accordingly, the structure of bulk viscosity up to second-order gradients is given as

  $$\begin{array}{ccc}\hfill \tau \dot{\mathsf{\Pi }}+\mathsf{\Pi }& =& -{\xi }_1(\nabla .u)+{\xi }_2({\mathcal{D}}_{\perp }^{\alpha }lnT){\mathcal{D}}_{\alpha }^{\perp }lnT\hfill \\ & +& {\xi }_3{(\nabla .u)}^2+{\xi }_{4}R+{\xi }_5{u}^{\alpha }{u}^{\beta }{R}_{\alpha \beta }+{\xi }_6{\sigma }^{\alpha \beta }{\sigma }_{\alpha \beta }+{\xi }_7{\Omega }^{\alpha \beta }{\Omega }_{\alpha \beta }\hfill \\ & +& {\xi }_8\left({\mathcal{D}}_{\perp }^{\alpha }\mu \right){\mathcal{D}}_{\alpha }^{\perp }\mu +{\xi }_9\left({\mathcal{D}}_{\alpha }^{\perp }\mu \right){\mathcal{D}}_{\perp }^{\alpha }\mu +{\xi }_{10}\left({\mathcal{D}}_{\perp }^{\alpha }\mu \right){\mathcal{D}}_{\alpha }^{\perp }lnT,\hfill \end{array}$$

  (8)
  The first line in rhs gives the first-order gradient terms. The other lines list out the possible second-order gradient terms for scalars, vectors and tensors constructing the dissipative flux quantities. T is temperature, $\mu$ is chemical potential, $-\nabla .u={\partial }_{t}lns$, and ${\mathcal{D}}_{\perp }^{\alpha }={\partial }_i$.

3. Equations of State of Viscous QCD Matter

The lattice QCD simulations are based on non-perturbative solution of the quantum chromodynamics (QCD) theory of the strong interactions. In high energy physics, the colored parton phase, the quark-gluon plasma, which is likely created in the relativistic collisions [23,24], is conjectured to expand with a fixed number of both baryon number and entropy, so the ratio entropy/baryon number remains constant, $S/N_B=51$. For each ratio, the phase diagram, T-${\mu }_B$, is obtained and also an EoS can be deduced. For the lattice simulations 4-stout improved staggered fermions with temporal lattice sizes 10, 12, and 16 with physical quark masses. are used. An isospin symmetry is assumed. Describing the state of the system of interest relies on EoS proposed; thermodynamic properties like pressure p, density $\rho$, and temperature T.

From recent lattice QCD calculations, at non-vanishing baryon chemical potential [22], different EoS could be derived. As presented in Figure 1, we parameterize the results as follows.

• For pressure p vs. temperature T

  $$\begin{array}{c}\hfill p\left(T\right)=a_1+a_{2}T+a_3{T}^2+a_4{T}^3,\end{array}$$

  (9)

  where $a_1=-0.851\pm 0.1672$, $a_2=21.094\pm 2.675$, $a_3=-176.893\pm 14.01$, and $a_4=506.068\pm 24.04$, Figure 1.
• For T vs. energy density $\rho$

  $$\begin{array}{c}\hfill T=b_1{\rho }^{b_2}+b_3\rho ,\end{array}$$

  (10)

  where $b_1=0.174\pm 0.001$, $b_2=0.131\pm 0.002$, and $b_3=0.006\pm 0.003$, Figure 2.
• For $\rho$ vs. T

  $$\begin{array}{c}\hfill \rho =c_1{T}^{c_2}+c_{3}T,\end{array}$$

  (11)

  where $c_1=2470.27\pm 7.597$, $c_2=4.316\pm 0.046$, and $c_3=-2.59\pm 0.021$, as represented in Figure 2.
• Then, the baratropic EoS can be given as

  $$p\left(\rho \right)=d_1{\rho }^{d_2},$$

  (12)

  where $d_1=0.162\pm 0.001$ and $d_2=1.202\pm 0.004$, Figure 3.
• For speed of sound squared $c_s^2$ vs. T

  $$\begin{array}{c}\hfill c_s^2=f_1{T}^{f_2}+f_3,\end{array}$$

  (13)

  where $f_1=1.97\pm 0.234$, $f_2=0.171\pm 0.029$, and $f_3=-1.25\pm 0.26$. The speed of sound squared can be derived as $c_s^2=\frac{\partial p}{\partial \rho }=0.196{\rho }^{0.208}$.
• For $c_s^2$ vs. ratio of bulk to shear viscosity

  $$\begin{array}{c}\hfill \frac{\zeta }{\eta }=15{\left(\frac{1}{3}-c_s^2\right)}^2,\end{array}$$

  (14)

  which is based on kinetic model [25], Figure 4. Then, the temperature dependence of the ratio of shear to bulk viscosity $\zeta /\eta$ can be parameterized as

  $$\begin{array}{c}\hfill \frac{\zeta }{\eta }=g_1{T}^{g_2}+g_3,\end{array}$$

  (15)

  where $g_1=-4.653\pm 0.103$, $g_2=0.9885\pm 0.034$ and $g_3=2.154\pm 0.031$.
• For entropy density s vs. T

  $$\begin{array}{c}\hfill s\left(T\right)=\frac{\partial p}{\partial T}=a_2+2a_{3}T+3a_4{T}^2.\end{array}$$

  (16)
  At vanishing chemical potential, $s/T^3=p/T^4+\rho /T^4$ [5]. Also, $s\left(T\right)$ can be obtained from differentiation of Equation (9) wrt T. Figure 5 plots $s\left(T\right)$ as a function of T. From the lattice QCD calculations and Equation (16), a simple parameterization can be deduced

  $$\begin{array}{c}\hfill s\left(T\right)=h_1{T}^{h_2}+h_3,\end{array}$$

  (17)

  where $h_1=358.89\pm 663.4$, $h_2=3.451\pm 0.143$, and $h_3=-2.85\pm 0.542$, Figure 5.
• For bulk viscosity normalized to entropy $\zeta /s$ vs. T

  $$\begin{array}{c}\hfill \frac{\zeta }{s}=15{\sigma }_0{\left(\frac{1}{3}-c_s^2\right)}^2,\end{array}$$

  (18)

  where ${\sigma }_0=0.01$ [25]. The results are depicted in Figure 6. The results can be best fitted as

  $$\begin{array}{c}\hfill \frac{\zeta }{s}=j_1{T}^{j_2}+j_3,\end{array}$$

  (19)

  where $j_1=-0.046\pm 0.001$, $j_2=0.99\pm 0.034$, and $j_3=0.022\pm 0.003$.
• Then, for $\zeta$ vs. T

  $$\begin{array}{c}\hfill \zeta =\left[a_2+2a_{3}T+3a_4{T}^2\right]\times \left[15{\sigma }_0{\left(\frac{1}{3}-c_s^2\right)}^2\right],\end{array}$$

  (20)

  where the entropy density is substituted from Equation (16) into Equation (18). Figure 7 shows our calculations, Equation (20), as a solid curve. The results can be fitted,

  $$\begin{array}{c}\hfill \zeta =I_1+I_{2}T+I_3{T}^2+I_4{T}^3,\end{array}$$

  (21)

  where $I_1=0.445\pm 0.008$, $I_2=-8.451\pm 0.125$, $I_3=48.415\pm 0.655$, and $I_4=-69.804\pm 1.124$.
• For $\eta$ vs. T

  $$\begin{array}{c}\hfill \eta =\frac{I_1+I_{2}T+I_3{T}^2+I_4{T}^3}{15{\left(\frac{1}{3}-c_s^2\right)}^2},\end{array}$$

  (22)

  where Equations (14) and (21) are substantiated. The results from Equation (22) are depicted as a solid curve in Figure 7 right. A generic expression for the dependence of the shear viscosity $\eta$ on T could be obtained

  $$\begin{array}{c}\hfill \eta =k_1+k_{2}T+k_3{T}^2+k_4{T}^3,\end{array}$$

  (23)

  where $k_1=0.211\pm 2.832\times {10}^{-5}$, $k_2=-3.538\pm 0.001$, $k_3=15.186\pm 0.002$, and $k_4=-0.008\pm 0.004$.

It it worth highlighting that the analytical fits as given by Equations (9)–(23) agree well with the relevant results reported in ref. [22].

4. Evolution of the Cosmic Parameters

At Bag pressure $B^{1/4}=200$ MeV, we plot in Figure 8 data taken from ref. [5] for the time evolution of the temperature. This can be parameterized as

$$\begin{array}{c}\hfill T=l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3,\end{array}$$

(24)

where $l_1=-0.062\pm 0.002$, $l_2=0.423\pm 0.009$, and $l_3=0.274\pm 0.002$. An expression for t vs. T is suggested as

$$\begin{array}{c}\hfill t=\frac{1}{t_{Pl}}{\left(\frac{T-l_3}{l_1}\right)}^{1/l_2},\end{array}$$

(25)

where $t_{Pl}$ is the reduced Planck time $t_{Pl}=4.608\times {10}^{-16}$ GeV${}^{-1}$.

• In order to obtain an expression for the dependence of H on T, we substitute $\rho$ taken from the lattice QCD calculations [22] into Equation (4), Figure 9

  $$\begin{array}{c}\hfill H\left(T\right)=(m_1{T}^{m_2}+m_3)t_{Pl},\end{array}$$

  (26)

  where $m_1=0.023\pm 0.001$, $m_2=1.821\pm 0.0323$, and $m_3=-0.0004\pm 2.064\times {10}^{-5}$.
• Then, the scale factor a can be calculated using Equation (4),

  $$\begin{array}{ccc}\hfill \frac{1}{a}da& =& {\left(\frac{1}{3m_{Pl}^2}\rho \right)}^{1/2}dt.\hfill \end{array}$$

  (27)
  When integrating over the comoving time t, we get

  $$\begin{array}{ccc}\hfill a\left(t\right)& =& exp\left(t\sqrt{\frac{1}{3m_{Pl}^2}\rho }\right).\hfill \end{array}$$

  (28)
  The results are depicted as symbols in Figure 9. These results can then be fitted as (solid curve in Figure 10)

  $$\begin{array}{c}\hfill a\left(t\right)=n_1+n_2\left(\frac{t}{t_{Pl}}\right)+n_3{\left(\frac{t}{t_{Pl}}\right)}^2+n_4{\left(\frac{t}{t_{Pl}}\right)}^3,\end{array}$$

  (29)

  where $n_1=1.0\pm 1.129\times {10}^{-5}$, $n_2=0.0011\pm 1.861\times {10}^{-5}$, $n_3=-0.0002\pm 7.077\times {10}^{-6}$, and $n_4=1.309\times {10}^{-5}\pm 7.192\times {10}^{-7}$.
• The temperature dependence of the scale factor can be derived from Equations (25) and (28)

  $$\begin{array}{ccc}\hfill a\left(T\right)& =& exp\left(\sqrt{\frac{1}{3m_{Pl}^2}\rho }\left[\frac{1}{t_{Pl}}{\left(\frac{T-l_3}{l_1}\right)}^{1/l_2}\right]\right).\hfill \end{array}$$

  (30)
  The calculations from this expression are depicted as symbols in Figure 10. The solid curve refers to the proposed expression

  $$\begin{array}{c}\hfill a\left(T\right)=q_1+q_{2}T+q_3{T}^2+q_4{T}^3,\end{array}$$

  (31)

  where $q_1=0.978\pm 0.001$, $q_2=0.368\pm 0.012$, $q_3=-1.819\pm 0.064$, and $q_4=2.803\pm 0.11$.
• The time evolution of the Hubble parameter can be deduced as follows. First we substitute T from Equation (24) into Equation (26)

  $$\begin{array}{c}\hfill H\left(t\right)=m_1{\left[l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3\right]}^{m_2}{t}_{Pl}+m_3{t}_{Pl}.\end{array}$$

  (32)
  Second, we differentiate the previous expression wrt t

  $$\begin{array}{c}\dot{H}=\left(m_1{m}_2{\left[l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3\right]}^{m_2-1}\left[l_1{l}_2{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2-1}\right]\right)t_{Pl}^2\hfill \\ =0.001T^{0.821}{\left(t{t}_{Pl}\right)}^{-0.577}{t}_{Pl}^2.\hfill \end{array}$$

  (33)
  The results are shown as solid curve in Figure 11.
• For the time evolution of the energy density, we first substitute T from Equation (24) into Equation (11)

  $$\begin{array}{c}\hfill \rho =c_1{\left[l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3\right]}^{c_2}+c_3\left[l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3\right].\end{array}$$

  (34)
  Then, we differentiate this equation wrt t

  $$\begin{array}{c}\dot{\rho }=c_1{c}_2\left({\left[l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3\right]}^{c_2-1}+c_3\right)\left\{\left(\frac{l_1{l}_2}{t_{Pl}^{-l_2}}\right)t^{l_2-1}\right\},\hfill \\ =-279.35\left\{-2.58+{\left[l_1{\left(\frac{t}{t_{Pl}^{-1}}\right)}^{l_2}+l_3\right]}^{3.316}\right\}{t}^{-0.577}{t}_{Pl}^{0.423}.\hfill \end{array}$$

  (35)
  The resulting $\dot{\rho }$ are depicted as a solid curve in Figure 11.
• For $\dot{\zeta }$, we first substitute T from Equation (24) into Equation (21) and then differentiate the obtained expression wrt to t,

  $$\begin{array}{c}\dot{\zeta }=3l_2{l}_1^3{I}_4{t}_{Pl}^{3l_2}{t}^{3l_2-1}+2l_1^2{l}_2\left(I_3+3l_3{I}_4\right)t_{Pl}^{2l_2}{t}^{2l_2-1}\hfill \\ +l_1{l}_2\left(I_2+2l_3{I}_3+3I_4{l}_3^2\right)t_{Pl}^{l_2}{t}^{l_2-1},\hfill \end{array}$$

  (36)

  which can be simplified as

  $$\begin{array}{ccc}\hfill \dot{\zeta }& =& [p_1{t}^{(3l_2-1)}+p_2{t}^{(2l_2-1)}+p_3{t}^{(l_2-1)}]t_{Pl},\hfill \end{array}$$

  (37)

  where $p_1=7.286\times {10}^{-22}$, $p_2=-3.126\times {10}^{-15}$, and $p_3=-2.01\times {10}^{-8}$, Figure 11.
• For $\dot{\eta }$, we first substitute T from Equation (24) into Equation (23) and then differentiate the expression obtained wrt to t,

  $$\begin{array}{c}\dot{\eta }=3l_2{l}_1^3{j}_k{t}_{Pl}^{3l_2}{t}^{3l_2-1}+2l_1^2{l}_2\left(j_k+3l_3{k}_4\right)t_{Pl}^{2l_2}{t}^{2l_2-1}\hfill \\ +l_1{l}_2\left(k_2+2l_3{k}_3+3k_4{l}_3^2\right)t_{Pl}^{l_2}{t}^{l_2-1},\hfill \end{array}$$

  (38)

  which can be simplified as

  $$\begin{array}{ccc}\hfill \dot{\eta }& =& [y_1{t}^{(3l_2-1)}+y_2{t}^{(2l_2-1)}+y_3{t}^{(l_2-1)}]t_{Pl},\hfill \end{array}$$

  (39)

  where $y_1=8.351\times {10}^{-26}$, $y_2=5.233\times {10}^{-15}$, and $y_3=-4.096\times {10}^{-8}$, Figure 11.
• Finally, we propose an expression for $\mathsf{\Pi }$ in dependence on H by using Equation (6) and $p_{eff}=p+\mathsf{\Pi }$,

  $$\begin{array}{c}\dot{\rho }=-3(p+\mathsf{\Pi }+\rho )H,\\ \mathsf{\Pi }=\left(\frac{\dot{\rho }}{-3H}\right)-\rho -p.\end{array}$$

  (40)
  We calculate $\mathsf{\Pi }$ by substituting $\dot{\rho }$, H, $\rho$, and p into Equation (40). The results are plotted as squares in Figure 12. $\mathsf{\Pi }$ in the non-casual Eckart frame, Equation (7), is presented as circles. From statistical fit, we propose

  $$\begin{array}{ccc}\hfill \mathsf{\Pi }& =& z_1{\left(\frac{H}{t_{Pl}}\right)}^2+z_2\left(\frac{H}{t_{Pl}}\right)+z_3,\hfill \end{array}$$

  (41)

  which is presented as a solid curve in Figure 10, where $z_1=(-5.977\pm 0.125)\times {10}^{11}$, $z_2=(8.726\pm 0.149)\times {10}^8$, and $z_3=(-51.186\pm 1.75)\times {10}^4$.

5. Results

Figure 1 shows the dependence of the pressure p on the temperature T. The symbols refer to the lattice QCD calculations [22] and the solid curve gives to the proposed expression, Equation (9). There is an excellent agreement. We find that increasing T leads to a rapid increase in p. This characterizes a phase transition from confined hadron to deconfined parton phase [23,24]. The condition for the transition from hadron to parton epoch is the temperature, where the density is conjectured to remain small.

The left-hand panel of Figure 2 depicts the temperature T as a function of the energy density $\rho$. The right-hand panel presents the dependence of the energy density $\rho$ on the temperature T. The symbols refer to lattice QCD calculations [22] and the solid curve stand for the proposed expressions, Equations (10) and (11), respectively. We find that T largely increases with $\rho$, and $\rho$ also increases with T. As concluded for p, $\rho$ can be taken as an order parameter. It should be noted that the hadron-parton transition is not discontinuous. It is known as a rapid crossover.

Figure 3 presents the dependence of the pressure p on the energy density $\rho$. The symbols with errorbars are the lattice QCD calculations along trajectories of constant ratio of entropy to particle number, $S/N=51$ [22]. The solid curve refer to the proposed EoS, Equation (12). We find that p increases as $\rho$ increases. Such a dependence slightly varies when moving from hadron to parton phase. The results on the speed of sound squared, which also expresses the EoS obtained, are depicted in Figure 13.

Figure 4 presents the ratio of bulk to shear viscosity $\zeta /\eta$ as a function of T. The symbols refer to calculations from Equation (14) [25], while the solid curve shows the fitting according to Equation (15). We find that $\zeta /\eta$ is inversely linearly proportional to T.

To determine the ratio $\zeta /s$ in dependence on T, we first determine the dependence of the entropy density s on T. Figure 5 depicts s as a function of T. The symbols depict the lattice QCD calculations evaluated from Equation (16). The solid curve refers to the proposed expression, Equation (17). We notice that s increases with the increase in T. Similar to p and $\rho$, s acts as an order parameter differentiating between confined (hadron) and deconfined (parton) phases. Accordingly, we conclude that the local law of thermodynamics is not violated, as s is always positive. This apparently refers to casual Eckart theory [5].

Figure 6 shows the ratio $\zeta /s$ as a function of T. The symbols stand for the calculations from Equation (18) which are based on Ref. [25], while the solid curve represents Equation (19). We notice that $\zeta /s$ is inversely linearly proportional to T.

Figure 7 depicts the dependence of the both types of viscosity on the temperature. The left-hand panel shows the bulk viscosity $\zeta$ as a function of T. Our calculations is depicted as a solid curve using Equation (20), while the dashed curve gives the statistical fit, Equation (21). The right-hand panel presents the shear viscosity $\eta$ as a function of T. Our calculations is depicted as the solid curve using Equation (22), while the dashed curve shows the statistical fit, Equation (23). We find that both viscosity coefficients rapidly increase with the increase in T, Equations (21) and (23). From the dependence of time t on temperature T, Figure 8, we can conclude that $\zeta$ and $\eta$ are inversely proportional to t. This agrees well with the conclusion drawn in Ref. [5].

Figure 8 shows the time dependence of the temperature T, at the bag constant $B^{1/4}=200$ MeV. The symbols refer to results taken from Ref. [5], while the solid curve presents their statistical fit, Equation (24). This expression is crucial for further calculations. Figure 9 depicts the temperature dependence of the Hubble parameter H, Equation (4). The solid curve represents the statistical fit to Equation (26).

The left–hand panel of Figure 10 shows the time evolution of the scale factor $a\left(t\right)$ during the quark-hadron phase transition, at bag constant $B^{1/4}=200$ MeV. The numerical values of the scale factors raise with increasing B, Equation (28). The solid curve presents Equation (29). We notice that $a\left(t\right)$ rapidly increases with the increase in t. Then, $a\left(t\right)$ begins to decrease with increasing t. At small t, $a\left(t\right)$ rapidly increases. This means that the viscous Universe, during QCD epoch, was likely rapidly expanding. But, with the increase in t, the bulk viscosity $\zeta$ increases. This shrinks $a\left(t\right)$ [5]. The right-hand panel of Figure 10 depicts the scale factor as a function of the temperature T, Equation (30). The solid curve presents our fit to Equation (31). We find that $a\left[t\right(T\left)\right]$ increases with the increase in the temperature. Then, at $T\approx 230$ MeV, $a\left[t\right(T\left)\right]$ begins to decrease with the increasing in T. This results from the high increasing in $\zeta$ with the increase in T, which slows down the Universe expansion. Therefore, the scale factor $a\left(t\right)$ likely decreases, at high temperature.

Figure 11 depicts the time evolution of the cosmological parameters $\dot{H}$, $\dot{\rho }$, $\dot{\zeta }$, and $\dot{\eta }$. The solid curves refer to our calculations using Equations (6), (33), (36) and (38), respectively. In top panels, we find that the time evolution of $\dot{H}$ and $\dot{\rho }$ exponentially decrease with time t. This is slowing down of the Universe expansion. In bottom panels, $\dot{\zeta }$ and $\dot{\eta }$ increase with the increasing t. This ranges from negative to positive values or is an increase in viscosity of the Universe which slows down the expansion of the Universe.

Figure 12 shows the dependence of $\mathsf{\Pi }$ on H. The circle symbols refer to calculations using non-casual Eckart theory [5], while squared symbols represent the calculations using casual Eckart theory. The solid curve refers to Equation (41). We notice that in non-casual Eckart theory [5] $\mathsf{\Pi }$ decreases with the increase in H, while in the casual Eckart theory the increase in $\mathsf{\Pi }$ with H is then followed by a smooth decrease forming a parabolic shape. The maximum $\mathsf{\Pi }$ is positioned at $t\sim 6\times {10}^{-4}$.

6. Conclusions

The present work aims at solving the field equations in the hadron and the parton epochs and proposing barotropic equation for the bulk viscous pressure $\mathsf{\Pi }$. We solve the Friedmann equations by proposing various equations of state using recent lattice QCD simulations in the Eckart frame. While the non-casual Eckart theory violates the second law of thermodynamics so that $\mathsf{\Pi }=-3\zeta H$ which is defined as the covariant entropy current stemming from Gibbs equation, the causal Eckart theory-similar the the causal stable Mueller–Israel–Stewart theory-leads to general relativistic causal hydrodynamical equations in the Eckart frame such as bulk viscosity, shear viscosity tensor, and heat flow obtained up to second-order gradients.

Using recent lattice QCD simulations, we propose barotropic equations of state for hadron and parton matter. By means of statistical fits, we propose expressions for various thermodynamic quantities and cosmological parameters including $\dot{H}$, $\dot{\rho }$, $\dot{\zeta }$, and $\dot{\eta }$. A relation between $\mathsf{\Pi }$ and H was introduced. The validation of the causal theory could be guaranteed. In non-casual Eckart theory, there is a linear dependence of the Hubble parameter H and the bulk viscous pressure $\mathsf{\Pi }$, which results form the first-order derivations of the equilibrium state, while in causal Eckart theory, in which bulk viscosity, shear viscosity tensor, and heat flow are obtained up to second-order gradients, the causality is guaranteed.

The scale factor starts with a rapid increase with the time. Then, it begins to decrease. We conclude that the viscous Universe, during QCD epoch, was likely rapidly expanding. But, with the time, the scale factor shrinks as the bulk viscosity increases.

We also conclude that $\dot{H}$ and $\dot{\rho }$ exponentially decrease with time, i.e., slowing down of the Universe expansion. Both $\dot{\zeta }$ and $\dot{\eta }$ increase with the time. As this ranges from negative to positive values, an increase in viscosity seems to slow down the expansion of the Universe. The dependence of shear and bulk viscosity on the expansion rate shall be investigated in a future work.

In the causal relativistic theory which conserves the charges and the particles and simultaneously allows finite heat and energy flow, the viscous pressure plays an essential role in drawing a reliable picture on the cosmic evolution of the early Universe Cosmological implications of the proposed causal Eckart frame is presented, for the first time.