Viscosity in Modified Gravity

Abstract

A bulk viscosity is introduced in the formalism of modified gravity. It is shown that, based on a natural scaling law for the viscosity, a simple solution can be found for quantities such as the Hubble parameter and the energy density. These solutions may incorporate a viscosity-induced Big Rip singularity. By introducing a phase transition in the cosmic fluid, the future singularity can nevertheless in principle be avoided.

1. Introduction

Modified gravity has become an active branch of modern cosmology, attempting to give a unified description of the early (inflationary) epoch of the universe and at the same time intending to account for the accelerated expansion at the later stages. Useful reviews on modified gravity theories can be found in [1,2,3].

Most treatises on modified gravity, as well as on standard gravity, assume the cosmic fluid to be ideal, i.e., non-viscous. From a hydrodynamicist’s point of view this is somewhat surprising, since there are several situations in fluid mechanics, even in homogeneous space without boundaries, where the two viscosity coefficients, the shear coefficient η and the bulk coefficient, ζ come into play. This means a deviation from thermal equilibrium to the first order. Such a theory in effect means acceptance of the Eckart 1940 theory [4]. An important property of the Eckart assumption is that the theory becomes non-causal. By taking into account second order deviations from thermal equilibrium, one can obtain a causal theory respecting special relativity. Pioneering articles on causal fluid mechanics are those of Müller [5], Israel [6], and Israel and Stewart [7]. A recent review can be found in [8]. Because of the assumption about spatial isotropy, the shear coefficient is usually omitted.

Our purpose in the following will be to include the bulk viscosity ζ in the modified gravity formalism. We consider the case when ζ is satisfying a scaling law, reducing in the Einstein case to a form proportional to the Hubble parameter. It turns out that this scaling law is quite useful. We survey first earlier developments along this line, extracting material largely from our earlier papers [9,10,11]. Thereafter, as a novel development, we investigate how the occurrence of a phase transition can change the development of the universe, especially in the later stages approaching the future singularity. (It may here appear natural to relate such a phase transition with the onset of a turbulent state of motion.) It is shown that such a transition may in principle be enough to prevent the singularity to occur at all. This part of the paper, covered in Section 4, is a generalization of the viscous/turbulent theory for standard cosmology recently given in [12,13].

2. Fundamental Formalism

The action in modified gravity is conventionally written in the general form

$$S=\int d^{4}x\sqrt{-g}\left[\frac{F(R)}{2{\kappa }^2}+{\mathcal{L}}_{\mathrm{matter}}\right]$$

(1)

where ${\kappa }^2=8\pi G$, and where ${\mathcal{L}}_{\mathrm{matter}}$ is the matter Lagrangian. The equations of motion are

$$-\frac{1}{2}{g}_{\mu \nu }F(R)+R_{\mu \nu }{F}^{\prime }(R)-{\nabla }_{\mu }{\nabla }_{\nu }{F}^{\prime }(R)+g_{\mu \nu }\square F^{\prime }(R)={\kappa }^2{T}_{\mu \nu }^{\mathrm{matter}}$$

(2)

where $T_{\mu \nu }^{\mathrm{matter}}$ is the energy-momentum tensor corresponding to ${\mathcal{L}}_{\mathrm{matter}}$.

We shall however in the following not consider the general case, but limit ourselves to the special form where

$$F(R)=f_0{R}^{\alpha }$$

(3)

with $f_0$ and α being constants. This model has been used before, by Abdalla et al. [14] and others. The case of Einstein gravity corresponds to $f_0=1$ and $\alpha =1$. This choice appears to be natural from a mathematical viewpoint, and in our case it will play an important role in connection with the scaling law for the bulk viscosity; cf. Equation (25) below.

We assume the spatially flat FRW metric

$$d{s}^2=-d{t}^2+a^2(t)d{\mathbf{x}}^2$$

(4)

and put the cosmological constant $\Lambda =0$. In comoving coordinates, the components of the four-velocity $U^{\mu }$ are $U^0=1,U^i=0$. Introducing the projection tensor $h_{\mu \nu }=g_{\mu \nu }+U_{\mu }{U}_{\nu }$, we have for the energy-momentum tensor

$$T_{\mu \nu }=\rho U_{\mu }{U}_{\nu }+\tilde{p}{h}_{\mu \nu }$$

(5)

where $\tilde{p}$ is the effective pressure

$$\tilde{p}=p-3H\zeta$$

(6)

The scalar expansion is $\theta =3\dot{a}/a=3H$, with H the Hubble parameter. The shear viscosity is here omitted.

The equations of motion following from the above action are

$$-\frac{1}{2}{f}_0{g}_{\mu \nu }{R}^{\alpha }+\alpha f_0{R}_{\mu \nu }{R}^{\alpha -1}-\alpha f_0{\nabla }_{\mu }{\nabla }_{\nu }{R}^{\alpha -1}+\alpha f_0{g}_{\mu \nu }\square R^{\alpha -1}={\kappa }^2{T}_{\mu \nu }^{\mathrm{matter}}$$

(7)

The equation of state for the fluid is written as

$$p=w\rho \equiv (\gamma -1)\rho$$

(8)

If $w=-1$ or $p=-\rho$, the fluid is a vacuum fluid with strange thermodynamical properties such as negative entropies (cf., for instance, [15]). Recent observations indicate that $w=-1.{04}_{-0.10}^{+0.09}$ [16,17]. It has been conjectured that w is a function varying with time, perhaps even oscillatory, and that w might have been around 0 at redshift z of order unity [18]. The quintessence region $-1<w<-1/3$ and the phantom region region $w<-1$ are both of physical interest. Both quintessence and phantom fluids imply the inequality $\rho +3p\le 0$, thus breaking the strong energy condition.

We now consider the (00) component of Equation (7), observing that $R_{00}=-3\ddot{a}/a$ and $R=6(\dot{H}+2H^2)$. With $T_{00}^{\mathrm{matter}}=\rho$ we obtain

$$\frac{1}{2}{f}_0{R}^{\alpha }-3\alpha f_0(\dot{H}+H^2)R^{\alpha -1}+3\alpha (\alpha -1)f_{0}H{R}^{\alpha -2}\dot{R}={\kappa }^2\rho$$

(9)

An important property of (9) is that the four-divergence of the LHS is equal to zero, ${\nabla }^{\nu }{T}_{\mu \nu }^{\mathrm{matter}}=0$ [19]. This is as in Einstein’s gravity, meaning that conservation of energy-momentum follows from the field equations. The energy conservation equation becomes

$$\dot{\rho }+(\rho +p)3H=9\zeta H^2$$

(10)

Differentiating Equation (9) with respect to t and inserting $\dot{\rho }$ from Equation (10), we get

$$\frac{3}{2}\gamma f_0{R}^{\alpha }+3\alpha f_0[2\dot{H}-3\gamma (\dot{H}+H^2)]R^{\alpha -1}+3\alpha (\alpha -1)f_0[(3\gamma -1)H\dot{R}+\ddot{R}]R^{\alpha -2}+3\alpha (\alpha -1)(\alpha -2)f_0{\dot{R}}^2{R}^{\alpha -3}=9{\kappa }^2\zeta H$$

(11)

Inserting $R=6(\dot{H}+2H^2)$, we see that this equation for $H(t)$ is quite complicated. We shall be interested in solutions related to the future singularity, and make therefore the ansatz

$$H=\frac{H_0}{X},\mathrm{where}X=1-B{H}_{0}t$$

(12)

where $H_0$ is the Hubble parameter at present time, and B a non-dimensional constant. If a future singularity is to happen, B must be positive.

Before closing this section, it is desirable to comment on stability issues for our ansatz Equation (3) for the modified Lagrangian. A theory of modified gravity should admit an asymptotically flat, static spherically symmetric solution. Now, we expect that the expression (2) for the complete action, with Equation (3) inserted, will not be the full solution. It is reasonable to expect that the modified part will contain also other terms so that Equation (2) makes up only a part of the complete action. Nevertheless, it is of interest to ask to what extent Equation (3), when take separately, will behave with respect to the stability requirements. In the Solar system, far from the sources, it is known that $R\approx {10}^{-61}$ eV${}^2$; it corresponds to one hydrogen atom per cubic centimeter. [Note that 1 eV = $5.068\times {10}^4$ cm${}^{-1}$.] On a planet, $R=R_b\approx {10}^{-38}$ eV${}^2$, whereas the average curvature in the universe is $R\approx {10}^{-66}$ eV${}^2$. According to the stability analysis of Elizalde et al. [20], the stability condition for matter is

$$F^{″}(R_b)>0,\mathrm{where}{R}_b\approx {10}^{-38}{\mathrm{eV}}^2$$

(13)

In our case, this means merely that the exponent α in the expression (3) has to be greater than one. The stability condition on α is quite modest.

3. Special Cases

It is now mathematically simplifying, and physically instructive, to focus on special cases.

3.1. Einstein Action

As mentioned above, Einstein’s gravity corresponds to $f_0=1$ and $\alpha =1$. It means that the Lagrangian is linear in R. It is natural to consider this case as a reference case before embarking on the nonlinear general situation.

We first have to adopt a definite form for the bulk viscosity. The simplest choice would be to put $\zeta =$ constant. There are however reasons to assume a slightly more complicated form, namely to put ζ proportional to the Hubble parameter H. This is physically natural, in view of the large fluid velocities expected near the future singularity. Such violent conditions should correspond to an increased value of ζ. We shall take ζ to be proportional to the scalar expansion, $\theta =3H$,

$$\zeta =3{\tau }_{E}H$$

(14)

with ${\tau }_E$ being the proportionality constant in the Einstein theory. An important property of this particular form, shown in [21], is that if ${\tau }_E$ is sufficiently large to satisfy the condition

$$\chi \equiv -\gamma +3{\kappa }^2{\tau }_E>0$$

(15)

then a Big Rip singularity is encountered after a finite time t. Even if the universe starts out from the quintessence region $(-1<w<-1/3)$ or $(0<\gamma <2/3)$, the presence of a sufficiently large bulk viscosity will drive it into the phantom region $(w<-1)$ and thereafter inevitably into the Big Rip singularity.

From the governing equations we get

$$B=\frac{3}{2}\chi ,H_0=\sqrt{\frac{1}{3}{\kappa }^2{\rho }_0}$$

(16)

where ${\rho }_0$ is the present ($t=0$) value of the energy density. The time dependent value ${\rho }_E$ for the energy density according to the Einstein theory becomes

$${\rho }_E=\frac{{\rho }_0}{X^2}$$

(17)

Here we ought to mention that other forms for the bulk viscosity, more complicated than the form Equation (14) above, have been suggested. One possibility is that ζ, in addition to the term proportional to H, contains also a term proportional to $\ddot{a}/a$. See further discussions on this topic in [22,23].

3.2. Modified Gravity Action

Consider now the modified gravity fluid for which $f_0$ and α are arbitrary constants. As before, we look for solutions satisfying the ansatz Equation (12). It turns out that such solutions exist, if we model the bulk viscosity ${\zeta }_{\alpha }$ according to the following scaling law [10,24],

$${\zeta }_{\alpha }={\tau }_{\alpha }{\theta }^{2\alpha -1}={\tau }_{\alpha }{(3H)}^{2\alpha -1}$$

(18)

We see that this scaling fits nicely with our results from the preceding subsection: if $\alpha =1$, our previous form Equation (14) follows. The time-dependent factors in Equation (11) drop out, and we get the following equation determining B,

$${(B+2)}^{\alpha -1}\left\{9(2-\alpha )\gamma +3[\alpha +3\gamma +\alpha (2\alpha -3)(3\gamma -1)]B+6\alpha (\alpha -1)(2\alpha -1)B^2\right\}=\frac{18{\kappa }^2}{f_0}{\left(\frac{3}{2}\right)}^{\alpha }{\tau }_{\alpha }$$

(19)

This equation is in general complicated. Let us consider $\alpha =2,\gamma =0$ as a typical example (recall that $\gamma =0$ corresponds to a vacuum fluid). Then we obtain from Equation (19) $({\tau }_{\alpha }\to {\tau }_2)$,

$$B^3+2B^2-\frac{9{\kappa }^2{\tau }_2}{8f_0}=0$$

(20)

If the LHS is drawn as a function of B, it is seen that there is a local maximum at $B=-4/3$ and a local negative minimum at $B=0$, irrespective of the value of ${\tau }_2$. For all positive ${\tau }_2$ there is thus one single positive root. This root is viscosity-induced and leads to the Big Rip singularity. When ${\tau }_2$ increases from zero, there is a parameter region in which there are three real roots. Assume this region, and introduce an angle $\varphi \in [0,{180}^0]$ such that

$$cos\varphi =-\left(1-\frac{243}{128}\frac{{\kappa }^2{\tau }_2}{f_0}\right)$$

(21)

then the actual value of the root can be expressed as

$$B=-\frac{2}{3}+\frac{4}{3}cos\left(\frac{\varphi }{3}+{240}^0\right)$$

(22)

For instance, if we choose $\varphi ={120}^0$, the positive solution becomes $B=0.3547$. According to Equation (12), this gives the following Big Rip time

$$t_{BR}=\frac{1}{B}\frac{1}{H_0}=\frac{2.819}{H_0}$$

(23)

We may also note the general relation for B following from the energy conservation Equation (10), when $\rho \to {\rho }_{\alpha },p\to p_{\alpha },\zeta \to {\zeta }_{\alpha }$,

$$B=-\frac{3\gamma }{2\alpha }+\frac{3{\tau }_{\alpha }}{2\alpha }\frac{{(3H_0)}^{2\alpha }}{{\rho }_0}$$

(24)

Here we used

$${\zeta }_{\alpha }={\tau }_{\alpha }{\left(\frac{3H_0}{X}\right)}^{2\alpha -1},{\rho }_{\alpha }=\frac{{\rho }_0}{X^{2\alpha }}$$

(25)

and for simplicity we used the same initial conditions at $t=0$ for the modified fluid as for the Einstein fluid, ${\rho }_{0\alpha }={\rho }_{0E}\equiv {\rho }_0$, and $H_{0\alpha }=H_{0E}\equiv H_0$.

4. On the Possibility of a Phase Transition in the Late Universe

In the preceding, we have surveyed bulk viscosity-induced generalizations of modified gravity, following essentially the earlier treatments in [10,11,24]. Our intention in the following, as a new contribution, will be to discuss the flexibility that the above model possesses with respect to sudden changes in the time development (we will refer to it as phase transitions) in the late universe. The main point is the different solutions for B in the governing Equation (11) that are possible when the scaling ansatz Equation (18) is inserted. We obtain the following algebraic equation for B, for definiteness still assuming $\alpha =2$,

$$B^3+(2+\frac{3}{4}\gamma )B^2+\frac{3}{2}\gamma B-\frac{9}{8}\frac{{\kappa }^2{\tau }_2}{f_0}=0$$

(26)

This equation generalizes Equation (20) to the case of non-vanishing γ.

Consider the following scenario: the universe starts out from present time $t=0$ and follows the equations of modified gravity, with a ${\tau }_2$-induced bulk viscosity corresponding to a positive value of B. That means the universe develops according to

$$H=\frac{H_0}{X},{\zeta }_2={\tau }_2{\left(\frac{3H_0}{X}\right)}^3,{\rho }_2=\frac{{\rho }_0}{X^4}$$

(27)

with $X=1-B{H}_{0}t$. The universe thus enfaces a future singularity at large times. Let now, at a fixed time that we shall call $t_{*}$, there be a phase transition in the cosmic fluid implying that the effect from ${\tau }_2$ goes away. It means that the further development of the fluid will be determined by the γ-dependent roots of Equation (26) when ${\tau }_2=0$. There are three roots:

(1) The first is $B=0$. This is the de Sitter case, corresponding to

$$H=H_{*},{\rho }_2={\rho }_{*}$$

(28)

where $H_{*}$ and ${\rho }_{*}$ follow from Equation (27) when $t=t_{*}$. By assuming that $\left|\gamma \right|\ll 1$, which is of main physical interest, we see that it is easy to determine the remaining two roots. One of them is

(2) $B=-2$. This means

$$H=\frac{H_{*}}{1+2H_{*}(t-t_{*})},\rho =\frac{{\rho }_{*}}{{[1+2H_{*}(t-t_{*})]}^2}$$

(29)

The accelerated expansion is accordingly reversed at $t=t_{*}$, and the density goes smoothly to zero at large times.

(3) The third root is $B=-3\gamma /4$, which yields

$$H=\frac{H_{*}}{1+\frac{3}{4}\gamma H_{*}(t-t_{*})},\rho =\frac{{\rho }_{*}}{{[1+\frac{3}{4}\gamma H_{*}(t-t_{*})]}^2}$$

(30)

The sign of γ is important here. If the equation-of-state parameter w lies in the quintessence region, $w>-1$ ($\gamma >0$), then the density of the universe will go to zero for large times, like for the case (2) above. By contrast, in the phantom region $w<-1$ ($\gamma <0$), the universe will actually move towards a Big Rip again, although very weakly so.

Finally, it is of interest to compare the above results with those obtained in ordinary viscous cosmology when the universe, similarly as above, is thought to undergo a phase transition at a definite time $t_{*}$. Such an investigation was recently carried out in [13] (the one-component case treated in Section VI). Consider the following model: the universe starts from $t=0$ as an ordinary viscous fluid with a constant bulk viscosity,

$$\zeta =\mathrm{constant}\equiv {\zeta }_0$$

(31)

and develops according to the Friedmann equations. Assume that the universe is in the phantom region, $\gamma <0$. It follows that in the initial period $0<t<t_{*}$,

$$H=\frac{H_0{e}^{t/t_c}}{1-\frac{3}{2}\left|\gamma \right|H_0{t}_c(e^{t/t_c}-1)}$$

(32)

$$\rho =\frac{{\rho }_0{e}^{2t/t_c}}{[1-\frac{3}{2}\left|\gamma \right|H_0{t}_c(e^{t/t_c}-1){]}^2}$$

(33)

where $t_c$ means the “viscosity time”,

$$t_c={\left(\frac{3}{2}{\kappa }^2{\zeta }_0\right)}^{-1}$$

(34)

According to these equations, the universe develops towards a Big Rip. Now, after $t=t_{*}$ we imagine an era for which ${\gamma }_{\mathrm{turb}}=1+w_{\mathrm{turb}}>0$ and an equation of state of the form

$$p_{\mathrm{turb}}=w_{\mathrm{turb}}{\rho }_{\mathrm{turb}}$$

(35)

Here the subscript “turb” refers to our association in [13] of the transition at $t=t_{*}$ into an era dominated by turbulence.

Then, for $t>t_{*}$,

$$H=\frac{H_{*}}{1+\frac{3}{2}{\gamma }_{\mathrm{turb}}{H}_{*}(t-t_{*})}$$

(36)

$$\rho =\frac{{\rho }_{*}}{{[1+\frac{3}{2}{\gamma }_{\mathrm{turb}}{H}_{*}(t-t_{*})]}^2}$$

(37)

This means a dilution of the density again, at large times. The Big Rip may thus be avoided, as a result of a phase transition in the cosmic fluid. We see that in this sense the behavior is similar in the two cases, modified or ordinary, gravity.

It ought to be made clear that we do not at present have a specific model of the phase transition suggested at $t=t_{*}$. Our association with a turbulent state of motion is however quite natural, on the basis of the following consideration: In states of violent local motions of the cosmic fluid near a future singularity, the transition into a turbulent kind of motion seems physically inevitable, as the local Reynolds number becomes then very high. That brings the shear viscosity concept back into consideration, now not in a macroscopic but in a microscopic (local) sense. We expect that there is an established distribution of eddies over the wave number spectrum. Most likely this distribution can be taken to be approximately isotropic, implying the existence of an inertial subrange in which the energy density is $E(k)={\alpha }_K{ϵ}^{2/3}{k}^{-5/3}$, where ${\alpha }_K$ is the Kolmogorov constant, ϵ is the mean energy dissipation, and k is the wave number. Ultimately, when the magnitude of k reaches the inverse Kolmogorov length $1/{\eta }_K={(ϵ/{\nu }^3)}^{1/4}$ with ν being the kinematic viscosity, the local Reynolds number becomes of order unity and heat dissipation occurs. What we have done above is to denote the post-transition period $t>t_{*}$ conceptually as a turbulent region, where the influence from the bulk viscosity has essentially gone away and where ${\gamma }_{\mathrm{turb}}$ has become positive, without going into further detail regarding the underlying physical transition process.

5. Conclusions

Starting from the modified gravity action integral Equation (5), we solved the (00) component of the governing Equation (7) inserting the scaling relation in Equation (12), $H=H_0/X$. The bulk viscosity ζ was assumed in the form Equation (18), generalizing the Einstein-case value Equation (14) frequently used in the literature. It turned out that the mentioned ansatz Equation (18) permitted solutions in the form of Equation (25), corresponding to future Big Rip singularities. This is thus a Big Rip scenario induced by the bulk viscosity.

In Section 4 it was discussed how the future singularity can nevertheless in principle be avoided, if one allows for a future phase transition in the cosmic fluid where the influence from viscosity goes to zero. Modified or conventional cosmology behaves in this sense essentially in the same way.

Finally, one may ask to what extent the above theory can be generalized to more complicated forms of the modified Lagrangian term than the simple power-law form given in Equation (3). Such a general expression for $F(R)$ would then have to be inserted into the field Equation (2), and thereafter to be combined with the energy conservation Equation (10). The general case seems difficult to handle, but there may be special cases that are mathematically tractable and of physical interest. Such an investigation is however outside the scope of the present paper. At present, the scaling laws like Equation (25) appears to be closely linked to our basic ansatz Equation (3).