Viscous Cosmologies

TL;DR

This paper investigates the effects of non-zero viscosities, especially shear viscosity, in cosmological models filled with Quark Gluon Plasma, demonstrating conditions for non-singular universe evolution.

Contribution

It introduces a causal viscous cosmology model incorporating shear viscosity from QGP, extending previous bulk viscosity models, and analyzes conditions for non-singular universe evolution.

Findings

Bulk viscosity leads to a non-singular universe with scale factor zero only asymptotically.

Shear viscosity can be compatible with FLRW metric under specific conditions.

The model with shear viscosity also results in a non-singular cosmological evolution.

Abstract

We probe into universes filled with Quark Gluon Plasma with non-zero viscosities. In particular, we study the evolution of a universe with non-zero shear viscosity motivated by the theoretical result of a non-vanishing shear viscosity in the Quark Gluon Plasma due to quantum-mechanical effects. We first review the consequences of a non-zero bulk viscosity and show explicitly the non-singular nature of the bulk-viscosity-universe by calculating the cosmological scale factor $R(t)$ which goes to zero only asymptotically. We further extend the model of bulk viscosity to include a Cosmological Constant. We contrast the previous results with the cosmology of universes with non-zero shear viscosity. We first clarify under which conditions shear viscosity terms are compatible with the Friedmann-Lama\^itre-Robertson-Walker metric. To this end we use a version of the energy-momentum tensor from the M\"uller-Isreal-Stewart theory which leads to causal Navier-Stoke equations. We then derive the corresponding Friedmann equations and show under which conditions the universe emerges non-singular.

Equations

$$K={\mathcal{R}}_{\mu \nu \alpha \beta }{\mathcal{R}}^{\mu \nu \alpha \beta }$$

$$K=12{\left(\frac{\ddot{a}}{a}\right)}^2+{\left(\frac{\dot{a}}{a}\right)}^4,$$

$$\frac{\eta }{s}\ge \frac{1}{4\pi }$$

$$\Delta T^{\alpha \beta }$$

$$h^{\alpha \beta }={\eta }^{\alpha \beta }+u^{\alpha }{u}^{\beta }.$$

$$\Delta T^{ij}$$

$$\Delta T^{i0}$$

$$\Delta T^{00}$$

$$\Delta T^{ij}=-\eta \left(\text{accentset}\circ {\nabla }^j{u}^i+\text{accentset}\circ {\nabla }^i{u}^j-\frac{2}{3}{\nabla }_{\mu }{u}^{\mu }{g}^{ij}\right)-\zeta \text{accentset}\circ {\nabla }_{\mu }{u}^{\mu }{g}^{ij}$$

$$d{s}^2=-d{t}^2+R^2(t)(d{r}^2+r^{2}d{\theta }^2+r^2{\sin }^2\theta d{\varphi }^2),$$

$$g_{ij}=R^2{\tilde{g}}_{ij}$$

$${\tilde{g}}_{ij}=\text{mbox}diag(1,r^2,r^2{\sin }^2\theta ).$$

$$\text{accentset}\circ {\nabla }_{\mu }{u}^{\mu }=\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left(\sqrt{-g}{u}^{\mu }\right)=\frac{1}{\sqrt{-g}}({\partial }_{\mu }\sqrt{-g})u^{\mu }+{\partial }_{\mu }{u}^{\mu }$$

$$\sqrt{-g}=r^2{\sin }^2\theta R^3(t)$$

$$({\partial }_{\mu }\sqrt{-g})u^{\mu }=({\partial }_0\sqrt{-g})u^0=3r^2{\sin }^2\theta R^2(t)\dot{R}(t)$$

$$\frac{1}{\sqrt{-g}}{\partial }_0(\sqrt{-g})=3\frac{\dot{R}(t)}{R(t)}=3H.$$

$$\text{accentset}\circ {\nabla }_{\mu }{u}^{\mu }=3H$$

$$\text{accentset}\circ {\nabla }_j{u}_i=\left({\partial }_j{u}_i-\text{accentset}\circ {\Gamma }_{ji}^{\lambda }{u}_{\lambda }\right)=-\text{accentset}\circ {\Gamma }_{ji}^0{u}_0=\text{accentset}\circ {\Gamma }_{ji}^0$$

$$\text{accentset}\circ {\Gamma }_{ij}^0=R(t)\dot{R}(t){\tilde{g}}_{ij}.$$

$$T_{ij}=\left[p{R}^2-3\zeta H{R}^2\right]{\tilde{g}}_{ij}.$$

$$-2\frac{\ddot{R}}{R}-\frac{{\dot{R}}^2}{R^2}=\kappa \left[p-3\zeta H\right].$$

$$p\to p^{\prime }=p-3\zeta H.$$

$$d{s}^2=-d{t}^2+\sum _{i=1}^3{a}_i^2(t)d{x}_i^2,$$

$$\eta$$

$$\zeta$$

$$p^{\prime }=\left(\gamma -1-3\alpha \frac{\dot{a}}{a}\right)\rho .$$

$$\frac{\ddot{a}}{a}=-\frac{4\pi G_N}{3}(\rho +3p^{{}^{\prime }}),H^2=\frac{8\pi G_N}{3}\rho .$$

$$\dot{H}=\frac{9}{2}\alpha H^3-\frac{3}{2}\gamma H^2.$$

$$\ln {\left[\frac{1}{H}\left|H-\frac{\gamma }{3\alpha }\right|\right]}^{\frac{3\alpha }{\gamma }}+\frac{1}{H}=\frac{3}{2}\gamma (t-t_0)+\frac{1}{H_0}+\ln {\left[\frac{1}{H_0}\left|H_0-\frac{\gamma }{3\alpha }\right|\right]}^{\frac{3\alpha }{\gamma }}.$$

$$H_{>}(t)=\frac{\gamma }{3\alpha \left[1+W\left(-e^{\frac{{\gamma }^2}{2\alpha }(t-t_0)}+f(H_0)\right)\right]},f(H_0)=\frac{\gamma }{3\alpha H_0}+\ln \left[\frac{1}{H_0}\left(H_0-\frac{\gamma }{3\alpha }\right)\right].$$

$$H_{<}(t)=\frac{\gamma }{3\alpha \left[1+W\left(e^{\frac{{\gamma }^2}{2\alpha }(t-t_0)}+g(H_0)\right)\right]},g(H_0)=\frac{\gamma }{3\alpha H_0}+\ln \left[\frac{1}{H_0}\left(\frac{\gamma }{3\alpha }-H_0\right)\right].$$

$$\dot{\rho }+3H\rho (\gamma -3\alpha H)=0.$$

Full text

Cosmology, Viscosity, Shear, Quantum Effects in Gravity

I Introduction

The phenomenological choice of a standard perfect fluid energy–momentum tensor for the standard cosmology has yielded appropriate results matching observations. Nevertheless, the existence of unresolved issues in cosmology has made cosmologists wonder if this could be due to a modification of the geometric part of the Einstein Field Equations (EFE) or of the energy-momentum tensor. Consider for instance the present accelerated stage of the universe U1 ; U2 . While the simplest way to reproduce it is to include a positive cosmological constant in the EFE (with all the consequences Boehmer ; U3 ; PositiveLambda ; ScalesLambda ), its present value interpreted as vacuum energy seems to be incompatible with contributions arising from the standard Quantum Field Theory. Therefore, some modifications of the Einstein equations or the energy-momentum tensor were suggested Chaplygin ; Quintessence ; Harko . Apart from the problem of the actual accelerated stage of the universe, there are two other problems associated with the early universe. One of them is the initial singularity when the cosmological scale factor goes to zero. The invariants calculated at this value of $a$ indicate that $a=0$ is a true singularity. For instance the Kretschmann scalar given as

[TABLE]

in Friedmann-Lamaître-Robertson-Walker (FLRW) cosmology yields

[TABLE]

which indeed tends to $\infty$ as $a\rightarrow 0$ and thus this sets a true singularity at $a=0$ KretschmannST . The other problem in the early universe is the choice of the inflationary scenario at the beginning of the universe. Many solutions for both problems have been suggested including scalar fields LindeScalar and higher order gravity Starobinsky ; Odinstov ; Odinstov2 for the inflationary mechanism and versions of quantum gravity LoopQG ; Bojowald ; NonC or quantum corrected cosmology NC for the first problem.

One scenario is particularly interesting as it requires only the modification of the energy-momentum tensor whose origin would be the Quark Gluon Plasma at the early stage of the universe. The appealing aspect of this scenario is the fact that one would indeed expect a deconfining phase of Quantum Chromo Dynamics at an early stage of cosmological evolution QGPNature ; QGPRafelski . The ingredients necessary to make the universe non-singular in this model could be the two possible viscosity terms in the energy-momentum tensor: the bulk and the shear viscosity. It has been shown by Murphy Murphy that bulk viscosity (plus additional assumptions about the spatial change of the fluid velocity) leads to a non-singular universe. A number of authors have also examined the consequences of the bulk universe Barrow1 ; Barrow2 ; Barrow3 ; Lima ; Barbosa ; viscosity1 ; viscosity2 ; Wilson2007 ; Pourhassan2013-1 ; Pourhassan2013-2 ; Brevik2017 , while recently there has been a rising interest in the study of self interacting viscous dark matter VDM1 ; VDM2 ; Blas ; VDM3 . The shear viscosity has been treated so far in a quite different context arguing that its presence is incompatible with the Friedman-Lamaître-Walker-Robertson (FLRW) metric requiring considerations of different Bianchi types of cosmologies Banerjee1985 ; Belinskii1975 ; Barrow4 ; Zimdahl ; Huang1990 ; Gron ; Banerjee1990 ; Singh2007 . The main obstacle to treat shear viscosity seems to have been terms with spatial derivatives of the velocities. Indeed, putting these derivatives to zero in the special relativistic energy-momentum tensor would make all viscosity terms (bulk and shear) vanish. However, as we will show below, going over to the general relativistic context there will be a non-zero residual effect for the bulk viscosity even if the derivatives of velocities are zero in the co-moving frame of the cosmological fluid. Consider, e.g., a total divergence $V\indices{{}^{\mu}_{;\mu}}=g^{-1/2}\partial_{\mu}(g^{1/2}V^{\mu})$ which will be non-zero even if one insists on $\partial_{\mu}V^{\mu}=0$. The bulk viscosity is proportional to a total derivative of velocities and its non-zero effect in the FLRW metric goes back to exactly the above result. It will become clear that shear viscosity cannot be treated in an analogous way. Paving its entry into homogeneous and isotropic universes is, however, possible by extending the standard viscous energy-momentum tensor. This extension is one of the possible versions of the Müller-Israel-Stewart theory Mueller ; Israel ; Stewart which contains a new parameter $\tau_{\pi}$ of dimension of time proportional to the shear viscosity. This new time scale makes the resulting Navier-Stokes equations causal NS . Several possible new energy–momentum tensors are possible and we refer the reader to NS ; Romatschkenew for more details. We shall come back to this subject in section IV and we mention here that our choice of a new causal shear energy–momentum tensor is motivated by simplicity. We use this new source (shear viscous energy–momentum tensor) in cosmology and show that it is compatible with the FLRW metric, i.e. with isotropy and homogeneity. We this consider as a main result of the paper since shear viscosity was treated in the past mostly in Bianchi types of cosmological models (see the discussion at the end of section II).

There are two reasons why such a result is of interest. For one, bulk viscosity is proportional to $[1/3-\partial p/\partial\rho]$ which for relativistic matter would be zero unless some unknown (quantum) effects would correct this null result. Secondly, and in contrast to the bulk viscosity case, it appears that the shear viscosity cannot reach an absolute zero. An equivalent version of the uncertainty principle in the Quark Gluon Plasma was studied to suggest that the shear viscosity in such a plasma could never reach zero. This was obtained by using string theory methods in strongly interacting Quantum Field Theories Son-2007 . Later using the AdS/CFT correspondence, this value was suggested as a lower bound for viscosity on such systems. The shear viscosity $\eta$ has a quantum mechanical lower bound Kovtun2005 of the form

[TABLE]

where $s$ is the entropy density. Although laboratory experiments on Quark Gluon Plasma have not shown the existence of any viscous terms, their sensitivity has not yet reached yet the quantum mechanical limit Heinz2011 ; Dress ; Teaney2009 . However, it is also worth noting that some models suggest the violation of such a limit under particular conditions, see for instance Pourhassan1 ; Pourhassan2 . Motivated by this result in QGP, and further driven by the possibility that the early universe was filled with QGP, we shall explore viscous QGP Policastro2001 ; Denielewicz1985 ; Chen2013 as a quantum effect in the cosmological evolution.

We first discuss the special and general relativistic versions of the energy-momentum tensor with shear and bulk viscosity in a frame co-moving with the (cosmological) fluid. Although it is usually assumed that bulk viscosity is compatible with the FLRW metric and the shear one not, it will become apparent that in the general relativistic context both can be treated on an equal footing. We then review briefly the bulk viscosity early universe calculating explicitly the cosmological scale factor. The results will depend on the constant $\rho_{0}/\rho_{Bulk}$ where $\rho_{0}$ is the initial value for the density and $\rho_{Bulk}$ is inversely proportional to the bulk viscosity. Taking upon the shear viscosity we show that it could also lead to non-singular universes provided one of the parameters we introduce as an initial value in the solutions of a differential equation is negative. However, the way they avoid the initial singularity is quite different from the bulk case.

Our main motivation is to study the fate of the initial singularity of the very early universe under the assumption of the extended energy–momentum tensor for shear viscosity. We therefore will always assume the almost-relativistic equation of state $p=\left(\frac{1}{3}+\epsilon\right)\rho$ allowing small deviations by including a small parameter $\epsilon$.

II Shear and Bulk Viscosity in General Relativity

Let us first establish the energy-momentum tensor with viscosities in a special relativistic context. We split the energy-momentum tensor into a sum $T_{\mu\nu}={\cal T}_{\mu\nu}+\Delta T_{\mu\nu}$ where ${\cal T}_{\mu\nu}$ is the energy-momentum tensor of a perfect fluid, i.e., ${\cal T}_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}$, and $\Delta T_{\mu\nu}$ is the part containing the viscosities. The full form of the latter has been studied by many authors (see for instance HawkingEllis ). Here we follow reference WeinbergApJ which gives it as

[TABLE]

with

[TABLE]

. Here $\eta^{\alpha\beta}$ is the Minkowski metric, $\eta$ the shear viscosity coefficient, $\zeta$ the bulk viscosity coefficient, $\xi$ the heat conduction coefficient and $\chi$ is due to a purely relativistic effect. It is worth pointing out that in the special relativistic context putting $\frac{\partial v_{i}}{\partial x_{j}}=0$ results in vanishing bulk and shear viscosities. This might have led some authors to the statement that the FLRW metric being isotropic and homogeneous is incompatible with shear viscosity. However, in the general relativistic context where $\eta_{\mu\nu}$ is replaced by $g_{\mu\nu}(x)$ and $\partial_{\mu}$ by the covariant $\nabla_{\mu}$ there remains a term proportional to the shear viscosity $\eta$ even if we put the partial derivatives of the velocity to zero. We think that this is what Murphy called “the motion of pure expansion” Murphy . To see this point let us start with Weinberg’s WeinbergApJ expression for $\Delta T_{\mu\nu}$ in a locally co-moving fluid with $u^{i}=0$ and $u^{0}=1$. It reads

[TABLE]

Going to general relativity (GR) one replaces $\delta_{ij}\rightarrow g_{ij}$, $\frac{\partial}{\partial x^{i}}\rightarrow\nabla_{i}$ and $\mathbf{\nabla}\cdot\mathbf{u}\rightarrow\accentset{\circ}{\nabla}_{\mu}u^{\mu}$ (for the sake of a compact notation in section III we reserve $\accentset{\circ}{\nabla}$ for the standard covariant derivate with the Christoffel connection denoted by $\accentset{\circ}{\Gamma}$) and obtains

[TABLE]

We shall take the cosmological FLRW metric which reads

[TABLE]

as well as the local system where $u^{i}=0$ (i.e. in the co-moving fluid frame). The spatial components of the metric can be written in the form

[TABLE]

with

[TABLE]

Let us examine the expressions above in more detail. We evaluate first the total divergence by writing

[TABLE]

where

[TABLE]

and

[TABLE]

so that

[TABLE]

The term $\partial_{\mu}u^{\mu}$ is just zero if the fluid is incompressible, namely $\vec{\nabla}\cdot\vec{v}=0$. Finally, we arrive at

[TABLE]

With the spatial covariant derivatives being

[TABLE]

we can also evaluate the explicit derivatives of the velocities in (5). Given that $u^{i}=0$ and $u^{0}=1$, the Chrisfoffel symbols $\accentset{\circ}{\Gamma}_{\mu\nu}^{\alpha}$ in FLRW of interest in our case are

[TABLE]

It is evident from the above that to get a non-zero effect of the bulk viscosity compatible with the FLRW metric we had to put the divergence of the velocities to zero. Similarly, requiring $\frac{\partial u^{i}}{\partial x_{j}}=0$ simplifies the above expression for $\Delta T_{ij}$. We finally obtain, adding $\Delta T$ to the perfect fluid energy-momentum tensor, the full expression as

[TABLE]

Einstein’s equations $G_{ij}\equiv{\cal R}_{ij}-(g_{ij}/2){\cal R}=\kappa T_{ij}$ (the $0-0$ components of Einstein’s equations give simply $\dot{R}^{2}=(\kappa/3)\rho R^{2}$ with $\kappa=8\pi G$ and they do not receive any contribution from the viscosities) in the FLRW metric give the following result

[TABLE]

This means that the effect of the shear and bulk viscosity coefficients in FLRW is to modify the pressure. This modification is explicitly given by

[TABLE]

This is not exactly the same expression mentioned by Murphy Murphy following the energy-momentum tensor of Landau and Lifshitz Landau1959 , but given there Murphy without derivation. Murphy’s expression reads $p\rightarrow p^{\prime}=p+(4\eta-3\zeta)H$. As far as the bulk viscosity is concerned we confirm the result in Murphy , but obtain a zero effect of the shear viscosity. We think that it is worthwhile to derive the effects of the viscosities to demonstrate the following conclusions explicitly. In the special relativistic framework putting the divergence of the velocity to zero implies no effect of the bulk viscosity. However, in the general relativistic framework we can have at the same time $\mathbf{\nabla}\cdot\mathbf{u}=0$ and a residual effect of the bulk viscosity (compatible with the FLWR metric). The special relativistic energy-momentum tensor with viscosities would vanish if all spatial partial derivatives of the velocity were zero. But, in the case of bulk viscosity in general relativity there will remain a non-zero contribution even if all derivatives of the velocity are zero in the co-moving cosmological fluid.

Regarding the zero effect of the shear viscosity we think that we are in agreement with literature here. It is claimed that working with the energy-momentum tensor outlined in the beginning of section II a non-zero shear viscosity in cosmology would require to work in Bianchi types of models. Consider for instance the Bianchi type I of cosmological models where the metric is taken in the form

[TABLE]

which has been considered by some authors Banerjee1985 ; Belinskii1975 ; Huang1990 . In such a case, in order to obtain explicit solutions we would need full information on the derivation of the velocities which are given by the Navier-Stokes equations. These Navier-Stokes equations would be the general relativistic version of the special relativistic version NS ; NSStrickland by the replacement $\eta_{\mu\nu}\to\ g_{\mu\nu}$ and $\partial_{\alpha}\to\accentset{\circ}{\nabla}_{\alpha}$. This is a much more challenging undertaking as the non-relativistic Navier-Stokes is already a complicated system. Going to the relativistic one complicates the matter and adding gravity (by doing the replacements mentioned) makes it highly complex. Fortunately, there might exist an alternative to make the shear viscosity compatible with the FLRW metric. The crucial point is to realize that the shear viscosity energy-momentum tensor which we used so far is not complete. Indeed, the Navier-Stokes equations based on this tensor are acausal and need to be remedied by adding new terms NS . We will come to this point in section IV after having discussed the cosmological implications of the bulk viscosity.

III Cosmology with Bulk viscosity

In order to have a comparison of the effects of the bulk versus shear viscosity let us briefly review the state of art of the bulk viscosity cosmology. In this context we mention two important results which relate the bulk and the shear viscosity with the energy density in a material medium with very short mean free times $\tau$. These relations are given as Weinberg ; Misner ; WeinbergApJ .

[TABLE]

where $\bar{a}$ is the Stefan-Boltzmann constant and $T$ the temperature. We can take the above expressions as proportional to the energy density $\rho$. In the case of non-zero bulk viscosity Murphy names the proportionality constant $\alpha$, namely $\zeta=\alpha\rho$. We will follow this convention. Some discussion of equation (21) is in order. If we take strictly the relativistic equation of state, $p=\frac{1}{3}\rho$, the bulk viscosity $\zeta$ will come out zero. What apparently is meant by $\zeta=\alpha\rho$ Murphy ; Brevik2017 ; Belinskii1975 ; Zimdahl is to allow a small deviation in the form $p=\left(\frac{1}{3}+\epsilon\right)\rho$ or to consider a small correction to (21). Either way we can speculate that such corrections might exist due to quantum mechanics which also, as an example, correct the classical equation for the ideal gas. In such a case $\zeta$ will be, albeit small, proportional to the energy density. Even if $\zeta$ is small the mere fact that it is non-zero can affect the nature of the initial singularity as shown below. This means we can rewrite the pressure as $p^{\prime}=p-3\alpha\rho(\dot{a}/a)$ with $p=(\gamma-1)\rho$ and end up with the expression

[TABLE]

Then, the Friedmann equations with $\Lambda=0$, $k=0$ and $H=\dot{a}/a$ become

[TABLE]

where we used the notation $a\equiv R/R_{0}$ with $R_{0}=R(t_{0})$ being an initial value. It is not difficult to verify that $H$ must satisfy the Abel equation of the first kind

[TABLE]

Note that for $\alpha=0$ we recover the case without viscosity. We first make use of the Abel equation since it is given completely in terms of the Hubble parameter $H$. Solving this equation under the initial condition $H(t_{0})=H_{0}$ yields

[TABLE]

We can solve analytically for $H$ by using the Lambert $W$ function Lambert . In the case $H_{0}>\gamma/(3/\alpha)$ we find

[TABLE]

If $H_{0}<\gamma/(3/\alpha)$ we obtain instead

[TABLE]

Taking the alternative path to arrive at a solution (i.e., determining first $\rho=\rho(a)$ via the continuity equation and using the latter in the second Friedmann equation in (23) $H^{2}\propto\rho$) we note that

[TABLE]

Using the second Friedmann equation in (23) and assuming radiation domination (the equation of state is of the standard form $p=(\gamma-1)\rho$ with $\gamma=4/3$ for relativistic matter) we can write the above equation in the following way

[TABLE]

where the $\pm$ signs come from the fact that $H=\pm\sqrt{\kappa\rho/3}$. The solution takes the form

[TABLE]

or, alternatively, solving the above equation for $\rho(a)$

[TABLE]

where we have set $a_{0}=1$ as it should be. Together with this result the second Friedmann equation in (23) can be used to infer the behavior of $a(t)$. From

[TABLE]

we obtain

[TABLE]

Here, we have defined the density $\rho_{Bulk}$ as

[TABLE]

which is a critical value for $\rho$, since $\rho=\rho_{Bulk}$ implies $\dot{\rho}=0$ as can be seen in the continuity equation (28). Note that in the standard case of zero viscosity, i.e, $\alpha=0$ one has $\rho_{Bulk}^{-\frac{1}{2}}\rightarrow 0$ and thus one recovers the expected solution of the form $a\propto\sqrt{t}$. The solution above can be further simplified if we define

[TABLE]

and we take the dimensionless variable

[TABLE]

making equation (33) take the form

[TABLE]

The particular case in which $\xi_{0}=1$ has a very special behavior. We find that for the first sign, i.e. positive $H$, we have

[TABLE]

while for negative $H$ we find

[TABLE]

The exponential solution makes it clear that the early universe has an inflationary expansion. It goes hand in hand with $H=const=4/(9\alpha)$ consistent with (24). Indeed, note that the choice of the parameter $\xi_{0}$ turns out to be rather important in general and it is related to parameter $C$ in Murphy . To see that we plot $a(\tilde{\tau})$ with $\tilde{\tau}=\tau-\tau_{0}$ in Figures 1-3 for the upper sign of equations (37) choosing $\xi_{0}$ below, equal and above 1. The condition to avoid the initial singularity has to do with $\xi_{0}$. For $\xi_{0}\leq 1$ the scale factor $a$ of the universe approaches zero asymptotically (as $\tilde{\tau}\to-\infty$) and hence avoids the singularity. The smaller the $\xi_{0}$ the more the universe displays a “coasting” non-singular character.

Moreover, from the first Friedmann equation in (23) we see that whenever

[TABLE]

there is an accelerated expansion. This is true all the time for the model in Figure 2 since $\rho^{1/2}=\rho_{Bulk}^{1/2}=\frac{4\sqrt{3}}{9\alpha\sqrt{\kappa}}$ implies that $H$ is also constant, while for the model plotted in Figure 2 this happens only for a finite amount of time. The models with initial values $\xi_{0}>1$ appear to have a similar behavior towards the past (i.e. $\tilde{\tau}\rightarrow-\infty$) also avoiding the singularity, but towards higher values of $\tilde{\tau}$, at some point this Universe ‘stops evolving’, or as Murphy writes this Universe “runs out of time” Murphy .

This becomes more transparent if we find the critical value $\tilde{\tau}_{crit}$ at which $\frac{d\tilde{\tau}}{da}=0$. We start by first finding its corresponding value $a_{crit}$ from equation (37) which yields

[TABLE]

where one finds that this value is a maximum for $a$. Note that this value is not always real for the expanding Universe. We use then equation (37) to obtain the $\tilde{\tau}_{crit}$ that corresponds to it:

[TABLE]

We have included this value explicitly in the plots.

This type of model will also have an accelerated expansion which will never reach its end and this is the main reason why it exhibits a non-invertible behavior at late times. This can be easily checked using equation (31) and noting that the density will always be above the value $\rho^{1/2}=\frac{2\sqrt{3}}{9\sqrt{\kappa}\alpha}$.

The lower sign in our solution of the Friedmann equations describes collapsing universes which abruptly end at a certain scale factor or expanding universes starting with a flat non-singular evolution at the beginning and also ending abruptly (see Figures 5, 5 and 6).

III.1 The case with $\Lambda\neq 0$ and $k=0$

In this case the Friedmann equations are

[TABLE]

Since $\rho\geq 0$ we have $H^{2}\geq\Lambda/3$. It can be easily verified that also in this case $H$ satisfies an Abel equation of the first kind, namely

[TABLE]

We mention here that the same steady state solution for $H$, namely $H=4/(9\alpha)$ which we encountered above with the bulk viscosity non-zero and $\Lambda=0$ is also possible for $\Lambda\neq 0$ and $\gamma=4/3$ (radiation) as a direct check reveals. The resulting exponential inflation does not contain the cosmological constant which enters only the expression for the constant density.

To get more insight into the consequences of the introduction of a positive cosmological constant we will employ different semi-analytical methods. We start by studying some global properties of the solution that can be evinced from the theory of autonomous differential equations. If we introduce the new dependent variable $h=\sqrt[3]{\alpha/\Lambda}H$ and the independent variable $\tau=(9\sqrt[3]{\alpha\Lambda^{2}}/2)t$, the above differential equation can be brought into the form

[TABLE]

Imposing the initial condition $h(\tau_{0})=h_{0}>0$ the general integral of the differential equation above reads

[TABLE]

The roots of the cubic polynomial appearing in the denominator of the integral are

[TABLE]

We always have one negative and two positive roots. Observe that $y_{2}<y_{3}$ whenever $\alpha\sqrt{3\Lambda}<\gamma$ and $y_{2}>y_{3}$ for $\alpha\sqrt{3\Lambda}>\gamma$. From the theory of autonomous differential equations $y_{1}$, $y_{2}$ and $y_{3}$ will represent the equilibrium solutions. Furthermore, $p(h)=h^{3}+a_{2}h^{2}+a_{1}h+a_{0}$ has local maxima and minima at

[TABLE]

with $p(h_{min})<0$, $p(h_{max})>a_{0}$, and $y_{2}<h_{min}<y_{3}$. Taking into account that we are only interested in the case of positive $h$, it is not difficult to check that $h$ increases whenever $0<h_{0}<y_{2}$, or $h_{0}>y_{3}$. On the other hand, $h$ decreases if $y_{2}<h_{0}<y_{3}$. We have the following allowed scenarios

$0<h(\tau)<y_{2}$ for $h_{0}\in(0,y_{2})$: this case corresponds to $H^{2}<\Lambda/3$ and is, by the second Friedmann equation in (43), unphysical. 2. 2.

$y_{2}<h(\tau)<y_{3}$ for $h_{0}\in(y_{2},y_{3})$: the Universe decreases in size and approaches a lower bound given by the horizontal asymptote $y_{2}$. 3. 3.

$h(\tau)>y_{3}$ for $h_{0}>y_{2}$: in this case the Universe grows without bound.

For an expanding Universe we are left with case 3. Either we have $h>y_{2}>y_{3}$ or $h>y_{3}>y_{2}$. This relates the bulk viscosity to the (positive) Cosmological Constant. In our Universe (with $\sqrt{\Lambda}\sim 10^{-42}$GeV2 Planck ) this puts a very weak bound on $\alpha$, i.e. $\alpha<\frac{4}{3\sqrt{3}}\frac{1}{\sqrt{\Lambda}}$.

Integrating (46) in the case $h_{0}>y_{3}$ yields the solution

[TABLE]

with

[TABLE]

If $y_{2}<h_{0}<y_{3}$, the solution reads

[TABLE]

Finally, if $0<h_{0}<y_{2}$, we get

[TABLE]

In examining the solutions obtained above, it is clear that above we assumed that explicit solutions are possible. If we drop this assumption and are satisfied with an implicit solution we can integrate (44) by standard methods. First of all, observe that (44) can be rewritten as

[TABLE]

Since from the second equation in (43) $H^{2}>\frac{\Lambda}{3}$, it follows that the solution of (53) will increase or decrease depending whether $H>\frac{\gamma}{3\alpha}$ or $H<\frac{\gamma}{3\alpha}$, respectively. The above equation can be trivially solved by the method of separation of variables. Taking into account the partial fraction expansion

[TABLE]

and an initial condition $H_{0}=H(t_{0})$ with $H_{0}^{2}>\Lambda/3$ the corresponding solution of the initial value problem associated to (53) can be obtained by solving

[TABLE]

Since for $H>\sqrt{\frac{\Lambda}{3}}$ we have the indefinite integral

[TABLE]

we find that the solution of our initial value problem becomes

[TABLE]

from which we conclude that

[TABLE]

The cosmological Abel equation (44) was also obtained in Tawfik as an approximation for the Müller-Israel-Stewart cosmological model Mueller ; Israel ; Stewart . However, in Tawfik this equation is solved using a $\tanh^{-1}$ function in place of $\coth^{-1}$ and an inappropriate argument in the logarithm (both violating the inequality $H^{2}\geq\Lambda/3$). Note that for $\alpha=0$ we obtain the correct solution for a universe with a positive cosmological constant. This suggests that the solution given by Tawfik should actually be replaced by (54).

It is interesting to study yet another road to the solution starting from the continuity equation. In order to avoid repetition we will focus below on the functional form of the energy density. We have

[TABLE]

Relying on the relation between $H$ and $\rho$ from the Friedmann equations means that we can cast the above equation in the form

[TABLE]

Solving this for the upper sign gives

[TABLE]

Note that as expected the previous results are recovered as $\Lambda\rightarrow 0$. Moreover, from this form of $a$ in terms of $\rho$ we may also observe that the critical density obtained above would be affected by the cosmological constant. In the case of an expanding Universe we would obtain the following modified critical density

[TABLE]

where $\rho_{\rm vac}=\frac{\Lambda}{8\pi G}$. From the continuity equation we get $\dot{\rho}<0$ (expected in an expanding universe) only if $0<H<4/(9\alpha)$. The Friedmann equation which relates $H$ to $\rho$ implies then $\rho<\rho_{Bulk,\Lambda}$. The effect of the Cosmological Constant is to diminish the critical density and change the condition for an accelerated Universe $\ddot{a}>0$ which will now be

[TABLE]

Equations (58) and (59) might be considered as the effects of $\Lambda$ in an early Universe with bulk viscosity. The other interesting characteristic of such a model is that $\Lambda$ would still account for the late time acceleration of the Universe, since we would expect viscosity to be negligible later. It is worth noticing that $\alpha$ can be estimated in the lattice gauge of Quantum Chromodynamics as suggested in Tawfik1 ; QCDL . Here we just give the end result which is in the form $\alpha=\frac{1}{9\omega_{0}}\frac{9\gamma^{2}-24\gamma+16}{\gamma-1}$, where $\gamma\simeq 1.183$ and $\omega_{0}\simeq 0.5-1.5$GeV. This gives a range for $\alpha$, namely $0.0823$ GeV*-1* $<$ $\alpha$ $<$ 0.247 GeV*-1*.

Implicit information on the density can also be obtained if we replace in equation (55) all terms with the Hubble parameter by $H=\pm\sqrt{\frac{\kappa}{3}\rho+\frac{\Lambda}{3}}$ . We arrive then at

[TABLE]

whose implicit solutions can be parametrized by $\Lambda=\kappa\rho_{\rm vac}$ and

[TABLE]

[TABLE]

For the the upper sign (the expanding Universe) the solution takes then the form

[TABLE]

It is of some interest to have a closer look at the global behavior of the density in an expanding/contracting universe with bulk viscosity. Take, e.g., a contracting universe with $H<0$. Then equation (28) guarantees that the density is decreasing and hence with $\dot{\rho}>0$ we get $\xi_{0}>1$. This almost self-evident fact is not true any more for an expanding universe with $H>0$. We have $\dot{\rho}<0$ if $1-(9\alpha H/4)>0$ (for the radiation case in (28)). It follows that $H<4/(9\alpha)$ which leads to $H=\sqrt{\kappa\rho/3}<4/(9\alpha)$. The immediate conclusion is that $\xi_{0}<1$. Expanding Universes with bulk viscosity and $\xi_{0}>1$ might have the paradoxical behavior (at least from the mathematical point of view) that their densities increase. We might exclude them on interpretative grounds. We will find a similar situation in the shear cosmology, but with reversed roles of expanding and collapsing.

IV Shear Viscosity

We further explore models with viscosity, motivated by the notion that the shear viscosity in the QGP may have a minimum bound as mentioned above. We put forward the question here as to what would happen if we had a fluid with shear viscosity instead of bulk. This undertaking is motivated further by the fact that bulk viscosity as given in equation (21) would vanish in a radiation dominated era (unless unknown quantum effects prevent this from happening).

The ideal energy-momentum fluid tensor ${\cal T}$ does not contain any dissipative (bulk or shear) terms. The dissipation is introduced by the term $\Delta T$ at the beginning of section III. We found that the only non-zero effect of the dissipation compatible with the FLRW metric is the bulk viscosity. As far as the bulk term is concerned we can leave the matter as it is and focus exclusively on the shear viscosity from now on. The special relativistic Navier-Stokes equations based on the energy-momentum tensor $T={\cal T}+\Delta T$ are acausal (we refer the reader to literature S1 ; S2 ; S3 ; S4 for details). Indeed, the speed of diffusion exceeds the velocity of light. There seem to be different ways to regulate the theory by introducing the relaxation time $\tau_{\pi}$ (and other transport coefficients) which for the sake of mathematical simplicity we keep constant proportional to the shear viscosity parameter $\eta$. This should suffice to answer the question how shear viscosity affects the early universe, especially its initial singularity. We note that different versions of density dependent $\tau_{\pi}$ exist in the literature Romatschkenew and we will come back to these versions in future investigations. We will explore here the most simple remedy in the form of Müller-Israel-Stewart theory based on the second law of thermodynamics. There exists also a version based on kinetic equations and a higher order theory based on conformal symmetry Baier . In other words, several versions of a causal shear viscous energy–momentum tensors are possible. Put simply, the classification scheme to distinguish them relies on conformal symmetry (or the lack of it) or the inclusion of geometric tensor like the Ricci and Riemann tensors, i.e. we can have conformal and non-conformal versions as well as types of shear energy–momentum tensors which themselves will contain the geometric tensors usually associated with the left hand side of the Einstein equations. Another tool to classify the causal shear energy–momentum tensor would be the case of equilibrium versus out–of–equilibrium. For details, we refer the reader to NS ; Romatschkenew . Our guiding principle is based on simplicity since the structure of the new energy–momentum tensor is usually complicated. We opt for the simplest viscosity case given in NS (see equation 67 below) called there the Boltzmann gas which we treat in two versions: traceless and non-traceless. We follow the reference NS which as most of the article uses the signature $(+,-,-,-)$. We therefore will also from now on use this convention. To be specific we have

[TABLE]

with its inverse given by

[TABLE]

such that $g_{\mu\nu}g^{\mu\nu}=\delta_{\alpha}^{\alpha}$. With this metric all the non-zero Christoffel symbols are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

As before we will be working in the comoving frame $u^{\mu}=(1,0,0,0)$. Following NS we introduce first some convenient abbreviated notation

[TABLE]

[TABLE]

[TABLE]

where $\accentset{\circ}{\nabla}_{\mu}$ corresponds to the covariant derivative with respect to the Christoffel symbols. The full energy momentum tensor will be given by

[TABLE]

where the behavior of the additional term $\pi_{\mu\nu}$ is given by the following equation

[TABLE]

Neglecting the term proportional to $\tau_{\pi}$ gives us back the viscosity contribution discussed in section III. Indeed, with the FLRW metric the term is zero as can be easily shown. We show that the first term on the right hand side is zero for $\mu=\nu=0$ and for $\mu=\nu=i$. Expanding our first order expression gives

[TABLE]

such that for the $\mu=\nu=0$ case we easily establish

[TABLE]

For the spatial components we have by virtue of $u^{i}=0$

[TABLE]

and a similar equation for $\accentset{\circ}{\nabla}_{\gamma}^{\alpha}$. Hence we conclude again that

[TABLE]

Together with $\nabla^{<0}u^{i>}=0$, which can be readily shown to hold, this confirms the results obtained in section II. A non-zero contribution is possible only through the term proportional to $\tau_{\pi}$ in (67). The compatibility with the FLWR metric requires the $\pi_{\mu\nu}$ part to be diagonal. In view of that we would like to have

[TABLE]

[TABLE]

which we will satisfy by choosing appropriate initial conditions. Starting with the $0-0$ component of equation (67)

[TABLE]

we notice that

[TABLE]

since the Christoffel symbols are zero. On the other hand we also have

[TABLE]

Putting the results together gives a differential equation for $\pi_{00}$

[TABLE]

which can be readily integrated in the form

[TABLE]

where $\pi^{00}(t_{0})\equiv\pi_{0}$.

For the case $\mu=i$, $\nu=j$ we can perform a similar analysis. Indeed, we can write

[TABLE]

It is easy to see that

[TABLE]

The full differential equation then reads

[TABLE]

Its solution is given by

[TABLE]

The condition that $\pi_{\mu\nu}$ be diagonal can be implemented by choosing appropriate initial conditions. Indeed, proceeding along the same lines as above one can also show that

[TABLE]

which yields the solution

[TABLE]

guaranting that by choosing appropriate initial conditions one can make $\pi_{\mu\nu}$ diagonal.

A few comments on the transport coefficient $\tau_{\pi}$ are in order. There exist different theories/estimates to calculate $\tau_{\pi}$ and $\eta$. Among them are the gauge/gravity duality Baier , BKG BGKc (Boltzmann equation in the form used by Bhatnagar, Gross and Krook BGK ) perturbative QCD pQCD and lattice QCD lQCD . In some of the approaches $\eta$ and $\tau_{\pi}$ are density dependent. In the BGK theory, used by many authors, one takes $\tau_{\pi}=\tau_{R}$ where $\tau_{R}$ is the relaxation time in the Boltzmann equation Romatschkenew . The relaxation time $\tau_{R}$ can be estimated by the mean free path MFP or by taking it constant BGKc . For instance in BGKc $\tau_{R}$ is taken as $\tau_{R}=0.5$ fm, which is much bigger than the Planck length.

IV.1 Conservation Laws

Let us recall that the Einstein equations will be given in the form $G_{\mu\nu}=-\kappa T_{\mu\nu}$. Since the left hand side satisfies $\accentset{\circ}{\nabla}_{\mu}G^{\mu\nu}=0$ it is stringent that our energy-momentum tensor fullfills

[TABLE]

We know from standard cosmology that $\accentset{\circ}{\nabla}_{\mu}{\cal T}^{\mu\nu}$ is zero for $\nu=i$, and that for $\nu=0$ it gives

[TABLE]

Thus we concentrate on the term related to $\pi^{\mu\nu}$

[TABLE]

Specializing first on $\nu=0$. i.e.,

[TABLE]

[TABLE]

we infer that

[TABLE]

On the other hand for $\nu=i$ we have

[TABLE]

from which it follows that

[TABLE]

We get different relations for each value of $i$. The $i=1$ case gives

[TABLE]

which in short form reads

[TABLE]

For $i=2$ we have

[TABLE]

from which we conclude that

[TABLE]

or alternatively

[TABLE]

[TABLE]

This way we can relate the two other components of the diagonal $\pi^{ij}$ to one component only

[TABLE]

[TABLE]

Now, since

[TABLE]

the modified continuity equation is

[TABLE]

This suggests a modification of the density and pressure in the following way

[TABLE]

IV.2 Friedmann Equations

Proceeding to calculate the Friedmann equations from the Einstein Field equations, namely $G_{\mu\nu}=-\kappa T_{\mu\nu}$, we need the $0-0$ and $i-i$ components of the Einstein and the full energy-momentum tensor. Starting with the $0-0$ components we see that

[TABLE]

and the Ricci scalar comes out to be

[TABLE]

This gives us the $0-0$ component of the Einstein tensor

[TABLE]

Thus the first Friedmann equation reads

[TABLE]

For the second Friedmann equation we take the $i-j$ with $i=j$ components of the Einstein field equations $G_{ij}=-\kappa({\cal T}_{ij}+\pi_{ij})$. It can be easily seen that

[TABLE]

On the other hand we have

[TABLE]

and

[TABLE]

leading to

[TABLE]

where we have used the form of the $\pi$-tensor from eq. (100). The second Friedmann equation can now be cast into the

[TABLE]

The two Friedmann equations can be shown to be consistent with the continuity equation we derived before. We see that all of the above suggests that the quantity $\pi^{\mu\nu}$ may be written in a similar way as ${\cal T}^{\mu\nu}$, namely

[TABLE]

so that

[TABLE]

The apparent formal steady state solution $\dot{R}=0$ of the Friedmann equation is possible if $\rho=-\pi^{00}$ and $\pi^{11}=-\frac{p}{R^{2}}$. If we insist on a non-empty Universe, this, however, leads to a contradiction. Since the velocities $\partial_{\beta}u^{\alpha}$ are zero, the covariant derivative takes the form $\accentset{\circ}{\nabla}_{\lambda}u^{\mu}=\accentset{\circ}{\Gamma}_{\lambda 0}^{\mu}$. Because of $\accentset{\circ}{\Gamma}_{0j}^{i}=\frac{\dot{R}}{R}\delta_{j}^{i}$ and $\accentset{\circ}{\Gamma}_{ij}^{0}=R\dot{R}\tilde{g}_{ij}$, the solution $R=$constant implies that all covariant derivatives are zero and hence also $\pi^{\mu\nu}$ (see eq. (67)). We end up with $\rho=p=0$.

IV.3 Traceless Case

We know that for radiation, the standard perfect fluid energy-momentum tensor is traceless, i.e., ${\cal T}\indices{{}_{\mu}^{\mu}}=0$. We can impose the same condition on the trace of $\pi^{\mu\nu}$ and study the consequences. Putting the trace to zero

[TABLE]

[TABLE]

gives

[TABLE]

which looks similar to the equation of state for radiation in which $p=\frac{1}{3}\rho$. We note that using the time dependence of $\pi^{00}$ and $\pi^{11}$ we also find a relation between the initial values, namely,

[TABLE]

The $\pi$-contribution to the energy-momentum tensor now reads

[TABLE]

while the Friedmann equations simplify to

[TABLE]

[TABLE]

From the continuity equation

[TABLE]

we infer the solution for the density

[TABLE]

Putting this into the first Friedmann we obtain

[TABLE]

the cosmological scale factor comes out as

[TABLE]

with $R_{0}=R(t_{0})\neq 0$. This corresponds to the standard case with the replacement $\rho_{0}\rightarrow\rho_{0}+\pi_{0}$. It is evident that we do not avoid the initial singularity here, as $R(t)$ becomes zero at some finite time. The only effect of the traceless viscosity energy-momentum tensor, is apart form the modification of $\rho_{0}$, the behavior of the density at early times

[TABLE]

which at later times goes over to the standard expression. It appears that, at least formally $H^{2}$ could be zero if $\pi_{0}<0$, leading to a possible bounce. However equation (123) says that in such a case

[TABLE]

and thus will lead to $R$ being constant, discarding the bounce possibility.

IV.4 The non-traceless case

In a more general case the conservation law

[TABLE]

together with the Friedman equations

[TABLE]

[TABLE]

are the determining cosmological equations. Noting that $\pi^{11}$ may be written as

[TABLE]

and taking $\rho$ from the first Friedmann equation, we arrive at

[TABLE]

Finally, we can write this in terms of $R$ and $t$ defining

[TABLE]

The differential equation for the scale factor $R(t)$ is given by

[TABLE]

where

[TABLE]

In terms of $a=(R/R_{0})$ this reads as

[TABLE]

This equation is non-linear and difficult to solve analytically, at least in an explicit form. However, with the following ansatz

[TABLE]

we can get a special solution. We see that (137) is indeed a special solution, provided we satisfy

[TABLE]

If we consider this to be a physical solution we have to respect also $a_{0}=1$ which makes (138) a relation between $\kappa\xi\pi_{0}$ and $\tau_{\pi}$. Furthermore, it follows that $\xi\pi_{0}$ has to be negative ($\xi=0$ brings us back to the traceless case). Mathematically, we can distinguish the two possibilities

[TABLE]

The first one sets a maximally possible value for $R_{0}$ in the form

[TABLE]

provided $\pi_{1}\equiv\pi_{11}(t_{0})$ is positive. The second possibility gives us a minimal value of $R_{0}$ in the form

[TABLE]

if $\pi_{1}$ is also negative. Since we do not know anything about $\pi_{1}$ and $\pi_{0}$ we can interpret the above results as a bound on $\pi_{1}/\pi_{0}$ or equivalently as a relation between these two. With the special solution of a collapsing universe at hand we can also obtain an expression for $\rho$ by noting that $H=-b$ and then putting this into equation (126). In this case we arrive at

[TABLE]

which corresponds to a constant density (provided we impose a condition to keep it bigger than or equal zero) as well as a constant $\pi_{00}$ in time

[TABLE]

It is convenient to rewrite (136) in dimensionless form by defining $\eta\equiv(t-t_{0})/\tau_{\pi}$ and introducing a dimensionless parameter $\alpha$ motivated by (137), such that

[TABLE]

Equation (136) is equivalent to

[TABLE]

The special solution corresponds to $sgn(\xi)sgn(\pi_{0})<0$, $\alpha=1$, $a_{0}=a(\eta=0)=1$ and $a^{\prime}(\eta=0)=-1/4$. We can then generalize our initial conditions as

[TABLE]

with $\beta$ a real number and $\alpha>0$. It is possible to determine under which conditions on the physical parameters and for which choice of the initial conditions the universe undergoes an initial acceleration/deceleration. For this purpose we look for a solution of (145) subject to the initial conditions (146) in the form of a MacLaurin series

[TABLE]

Using (145) gives immediately

[TABLE]

Since $\alpha>0$ there will be an initial acceleration when $\xi$ and $\pi_{0}$ have opposite signs and $\alpha>8\beta^{2}$. If $\xi$ and $\pi_{0}$ have the same sign, the universe will initially decelerate. We also note that the density $\rho$ can be written in this dimensionless form by first using equation (129) to obtain

[TABLE]

and in order to leave it fully dimensionless we write

[TABLE]

where we have used equation (144) to obtain $\alpha$. We have plotted $a(\eta)$ and $\sigma(\eta)$ as numerical solutions for different values of $\alpha$ and $\beta$ as well as different signs of $\pi_{0}$ and $\xi$.

With different choices of the initial value and signs, different universes emerge. In Figure 7 we have plotted the scale factor $a$ for the special solution describing a collapsing universe and its constant density (in Figure 8). Most importantly there are bouncing universes whose scale factor can be seen in Figures 9, 15 and 17 with the corresponding densities depicted in 10, 16 and 18. It is clear from the first Friedmann equation (107) that a bouncing universe has to have a negative $\pi_{00}$ and hence a negative $\pi_{0}$ which comes out correctly when we plot the corresponding densities. Singular universes, expanding, contracting or recollapsing, emerge whenever we choose sgn$(\xi\pi_{0})$ to be positive as shown in Figures 11, 12 (the recollapsing case), 13 and 14.

The choice of the FRLW metric and the Einstein equations fix the behavior of the universe. The advantage of working with dimensionless parameters enables us to examine the global behavior of the universe. One can pose the valid question addressing the fate of the viscosity during and after the phase transition to hadrons (or what will happen in the matter dominated universe). As for now we have to leave the answer open. However, there exist theories in which the dark matter also exhibits a viscous component VDM1 ; VDM2 ; Blas ; VDM3 . We notice that the energy–momentum tensor used in Blas in the context of dark matter contains only a subset of possible terms. This could be a starting point to use a more general viscous energy–momentum tensor (e.g. equation (67) in the study of dark matter in the matter dominated epoch of the universe).

V Alternative Version of $\pi^{\mu\nu}$ for a fluid with Shear Viscosity

In this section we will briefly touch upon one of the different versions of the energy-momentum tensor and present some preliminary results. This alternative version of the behavior of a fluid with shear viscosity NS reads

[TABLE]

where

[TABLE]

and in passing we note that yet another, third version can be consistently derived. We explore the possibility given in (151) in order to make a first comparison with our previous model. We will leave, however, the details for a future publication. We use the same notation as in section IV. We concentrate first on the case with $\mu=\nu=0$ where we already found that $\nabla^{<0}u^{0>}=0$. It is straightforward to calculate the other terms appearing in (151). We obtain for the first two terms in the square brackets in (151)

[TABLE]

[TABLE]

and for the vorticity term,

[TABLE]

Using the previous definition for $\nabla^{<\mu}u^{\nu>}$ we have

[TABLE]

The full equation then reads

[TABLE]

We proceed to check the $\mu=i$ and $\nu=j$ case. As before we have $\nabla^{<i}u^{j>}=0$. Continuing with the next terms we can write

[TABLE]

[TABLE]

For the $\Omega$ term we get

[TABLE]

For the above to be non-zero we require that $\alpha=j$ in the first term within the brackets while $\alpha=i$ in the second term within the brackets. It is worth to recall that

[TABLE]

[TABLE]

With this at hand we deduce

[TABLE]

With

[TABLE]

[TABLE]

[TABLE]

we arrive at the result

[TABLE]

The last relevant term takes the form

[TABLE]

and the full equation now reads

[TABLE]

For $\mu=i$ and $\nu=0$ we calculate the following expressions

[TABLE]

[TABLE]

[TABLE]

This enables us to write

[TABLE]

V.1 Conservation Laws

As done in the previous case, whatever the modification of the energy-momentum tensor we expect it to fulfill its conservation given by the Bianchi identities for the Einstein tensor. This conservation is given by

[TABLE]

From the standard energy momentum tensor we know that when $\nu=i$, we get $\accentset{\circ}{\nabla}_{\mu}\mathcal{T}^{\mu i}=0$ and thus we should have

[TABLE]

[TABLE]

Under the assumption that the $\pi^{\mu\nu}$ are not dependent upon spatial coordinates we have

[TABLE]

The above equation sets some conditions for the different $i$ cases, namely,

[TABLE]

[TABLE]

[TABLE]

One can show that the assumption that the tensor $\pi^{\mu\nu}$ be diagonal does not lead to any contradiction. Let us take the equation for $\pi^{i0}$ , namely equation (175) and let us assume that $\pi$ is diagonal. We obtain

[TABLE]

Both sides are equal to zero under the assumption of diagonality and hence consistent. For the $\pi^{ij}$ equation (169) taking $i\neq j$ we get something equivalent after noting that $g^{ij}$ is zero for $i\neq j$. We have

[TABLE]

which again yields zero on both sides, thus the diagonal assumption of $\pi^{\mu\nu}$ is consistent. We note that the continuity equation will also set the conditions already found in the previous version of $\pi^{\mu\nu}$, namely, those given in equations (99) and (100). With this consistency in mind from now on we work with a diagonal $\pi$ tensor so that the equation above for $\pi^{00}$ actually reads

[TABLE]

Solving for $\pi^{00}$ we may write

[TABLE]

In contrast to the previous case, this is an algebraic equation, while the differential equations will give

[TABLE]

where we have used

[TABLE]

We note that the differential equations are all equivalent when this condition is applied. The relevant differential equation reads

[TABLE]

With $H=\dot{R}/R$ we can re-write it in the following form

[TABLE]

We will address the possible solutions in an upcoming work, but we notice that in comparison to our previous case this equation is inhomogeneous. Furthermore, it is straightforward to note that the Friedmann equations obtained in the previous section will be the same for this case, only that the behavior of $\pi^{\mu\nu}$ is now determined by equation (151). This implies that the first Friedmann equation may be written as

[TABLE]

Recalling that the bounce is likely to happen if at some point $H=0$, we see that in this case this possibility can be realized as long as the density takes the value

[TABLE]

The additional condition for the bounce to be possible is that $\ddot{R}>0$. With the second Friedmann equation in the form

[TABLE]

this implies that for the possible bounce at $t_{B}$ we have

[TABLE]

The bounce happens at an early time during the evolution of the Universe. We can therefore again make use of the equation of state $p=\frac{1}{3}\rho$ leading to

[TABLE]

[TABLE]

or

[TABLE]

This last condition is particularly interesting, given that one usually expects a minimum value for $R$ at the moment of the bounce. We also find another result of interest from the time derivative of $\pi^{00}$,

[TABLE]

which in the condition for the bounce implies

[TABLE]

which suggests that $\pi^{00}$ is decreasing at the moment of the bounce.

We also note that the Friedmann equation may be written as

[TABLE]

from which one can obtain $H(\rho)$ or $\rho(H)$ which may be useful for an eventual solution of the equations obtained. The full exploration of this version is beyond the scope of this paper and we postpone the details to a future publication. The role of viscosity has been investigated in a similar fashion in Belinski1 ; Belinski2 , however this was done, as far as we can see, with yet another version of the extended energy–momentum tensor.

VI Conclusions

Whereas there seems to be a general agreement that bulk viscosity is still compatible with the FLRW metric, one often finds the statement that this is not true for the shear part. As we have shown, both can enter a homogeneous, isotropic cosmology and the reason for this can be traced back to the difference of the energy-momentum tensor in general relativity as compared to the analogous expression in the special relativistic context. In the latter case the vanishing of the divergence of velocities results in the vanishing of the bulk viscosity and the vanishing of partial derivatives of velocities leads to a null effect of the shear viscosity. With this mind it appears that the bulk and shear viscosities are not compatible with a homogeneous, isotropic spacetime. This conclusion is erroneous in the context of general relativity when derivatives are replaced by covariant counterparts. Even with constant velocities in the co-moving frame, there will be a non-zero effect of the bulk viscosity given by (18). We agree here with Murphy Murphy who mentioned this result without an explicit derivation. Using this we reviewed the main results of cosmology with bulk viscosity. The bulk universe avoids the initial singularity by following an asymptotic behavior at initial times (i.e., the scale factor approaches zero asymptotically).

However, we cannot confirm Murphy’s result regarding the shear viscosity. Indeed, with the simplest choice of the energy-momentum tensor which contains the shear viscosity we obtain a zero result when trying to make the shear case compatible with FLRW universes. The above mentioned simplest choice of the energy-momentum tensor leads to acausal Navier-Stokes equations at the special relativistic level and needs additional terms to become consistent. One such improved version has been used in the present paper introducing a new time scale $\tau_{\pi}$. It is exactly this part of the new energy-momentum tensor which is compatible with the FRLW metric leading to a panoply of universes depending on the signs and parameters chosen. It would also be a worthwhile undertaking to examine the new causal shear viscous tensor in the context of non-isotropic metrics, i.e. Bianchi type of cosmology metrics.

We have briefly touched upon a second possible version of the energy-momentum tensor which leads to different differential equations. As compared to the first case these equations are inhomogeneous. It would be interesting to see what consequences this has on the cosmology. Finally, we mention that there exists also a third possible version of a viscous energy-momentum tensor which includes the Riemann tensor. Needless to say it would be also of interest to explore this model further. We intend to come back to these problems in a future publication. We have not probed into the question of the hadronic phase transition from QGP to hadrons. This is an important issue in the cosmological context and it is partly examined in Tawfik3 .

Acknowledgements.

We would like to thank Paul Romatschke and Michael Strickland for very useful discussions and comments.

Appendix A Reduction of (136) to an Abel equation and its solution

It is possible to obtain an implicit analytical solution of (136) which we briefly sketch below. Even though we do not make use of it, it might prove useful in future research.

By means of the particular solution to (136)

[TABLE]

it is possible to construct an ansatz

[TABLE]

which leads after rescaling $t=\tau/b$ to the following differential equation

[TABLE]

By means of the transformation

[TABLE]

it is not difficult to verify that (205) can be cast into the form of an Abel equation of the $2^{nd}$ kind, namely,

[TABLE]

Observe that in the case we succeed to find the unknown function $w$, the transformation (206) will give $\tau$ as a function of $\varphi$ and at this step it is not clear at all if it will be possible to invert this relation. Equation (207) can be brought into a simpler form by the substitution

[TABLE]

which yields

[TABLE]

Finally, by means of the transformation

[TABLE]

equation (209) can be cast into its canonical form, more precisely

[TABLE]

We follow Greek to construct exact analytic solutions to (211) for $z>0$. To this purpose we introduce the transformation

[TABLE]

where $h$ and $V$ are unknown differentiable functions of $z$ on $(0,\infty)$. Then, (211) becomes

[TABLE]

The introduction of a further unknown differentiable function $U=U(z)$ followed by addition and subtraction of the term $UdV/dz$ in (213) gives

[TABLE]

By means of a further unknown function $G=G(z)$, it is possible to split (214) into the following Abel’s equations

[TABLE]

According to Theorem $1$ in Bou if there exists a constant $\lambda$ such that

[TABLE]

then equation (215) admits a solution

[TABLE]

Similarly, if it is possible to find a constant $\widetilde{\lambda}$ such that

[TABLE]

then equation (216) has a solution

[TABLE]

Combining (217) and (219) yields an integral equation for $h$ whose solution is

[TABLE]

with $c$ an arbitrary integration constant. Substituting (221) into (219) and using (217) to eliminate the dependence on $U$ gives the relation $\lambda=-\widetilde{\lambda}$. According to (217) the simplest choice for $\lambda$ is $\lambda=2$. Then, we have

[TABLE]

It remains to determine the subsidiary function $G$ and the function $V$. To find $V$ we need to require that the quadratic equations (218) and (220) which can be cast in the form

[TABLE]

have one common root. Let

[TABLE]

Then, the quadratic equations in (223) will have a common root whenever

[TABLE]

Squaring (225) followed by differentiation with respect to $z$ gives a cubic equation for $\sqrt{1+\Psi(z)}$, namely

[TABLE]

By means of the transformation

[TABLE]

we can bring (226) into its normal form

[TABLE]

with

[TABLE]

Depending on the sign of the discriminant

[TABLE]

there are the following scenarios Bron :

For $D>0$ there are only one real and two complex solutions. More precisely, if $D>0$ and $p<0$ the only real solution is

[TABLE]

where the sign of $r$ must be chosen so that it coincides with the sign of $q$. If both $D$ and $p$ are positive, we have

[TABLE]

and $r$ is defined as in (231). 2. 2.

In the case $D<0$ and $p<0$ there are three distinct real roots given by

[TABLE]

with $\cos{\widetilde{\varphi}}=q/(2r^{3})$. 3. 3.

If $D=0$ there is one real solution with algebraic multiplicity three if $p=q=0$ or two real solutions (one with algebraic multiplicity one the other having algebraic multiplicity two) whenever $p$ and $q$ do not vanish at the same time. Note that the case $p=0$ and $q=0$ is never satisfied because these two conditions give rise to two distinct subsidiary functions. Finally, if $p$ and $q$ do not vanish but $D=0$, we find

[TABLE]

The solution of (211) reads

[TABLE]

with $Z$ given as in (231), (233), or (234). A general formula to determine the subsidiary function has been provided by Greek . Here, we limit us to provide the main results, namely

[TABLE]

with $\omega=\ln{|z+2c|}$ and $\mbox{ci}(\cdot)$ denoting the cosine integral given by Grad

[TABLE]

where $\gamma$ is the Euler-Mascheroni constant. By means of (208) and (210) we find that

[TABLE]

Finally, (206) yields

[TABLE]

with $K$ an arbitrary integration constant.