A Generalized Solution of Bianchi Type-V Models with Time-Dependent G and Λ

Abstract

We study the homogeneous but anisotropic Bianchi type-V cosmological model with time-dependent gravitational and cosmological “constants”. Exact solutions of the Einstein field equations (EFEs) are presented in terms of adjustable parameters of quantum field theory in a spatially curved and expanding background. It has been found that the general solution of the average scale factor a as a function of time involved the hypergeometric function. Two cosmological models are obtained from the general solution of the hypergeometric function and the Emden–Fowler equation. The analysis of the models shows that, for a particular choice of parameters in our first model, the cosmological “constant” decreases whereas the Newtonian gravitational “constant” increases with time, and for another choice of parameters, the opposite behaviour is observed. The models become isotropic at late times for all parameter choices of the first model. In the second model of the general solution, both the cosmological and gravitational “constants” decrease while the model becomes more anisotropic over time. The exact dynamical and kinematical quantities have been calculated analytically for each model.

1. Introduction

The Bianchi models are described as classes of non-standard cosmological models that are in principle spatially homogeneous but anisotropic. Alternatively, they can be considered as a generalization of the well-known standard Friedman–Lemaître–Robertson–Walker (FLRW) models of cosmology. Originally, these models date back to the work of Bianchi [1,2], who classified them according to their construction of homogeneous surfaces in space-time. These surfaces are constructed by the action of a three-dimensional group of isometrics $G_3$ upon the space-like 3−surfaces [3]. From the cosmological point of view, these models are of great interest because they provide a way of studying the anisotropy at an early period of our universe’s expansion history [4].

One of the most puzzling and unsolved problems in physics today is the Cosmological Constant Problem [5] in the Einstein field equations (EFEs). In cosmology, it is regarded as a matter field with negative pressure (or as a vacuum energy density) that drives the accelerated expansion of the universe. The calculations of quantum field theory predict the value of vacuum to be of the order of ${10}^{130}$ s${}^{-2}$, whereas the cosmological constant value is of the order of ${10}^{-20}$ s${}^{-2}$; hence, a new thought is required to explain or resolve this puzzle. The idea of considering the cosmological models as varying vacuum energy density $\mathsf{\Lambda }$ has extensively been adapted by several authors [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. Chen [23] considered $\mathsf{\Lambda }$ proportional to $1/a^2$, and a generalized form $\mathsf{\Lambda }=\alpha /a^2+\beta H^2$, which depends on adjustable parameters $\alpha$ and $\beta$ of the quantum field on a curved and expanding background, the Hubble parameter H and the average scale factor of the universe a was also introduced [24].

Since the Newtonian constant G is a coupling “constant” between the geometry of spacetime and energy in the general theory of relativity (GR), and the universe is evolving with time, it is natural to assume that G varies with time. This consideration was first introduced by Dirac [25]. Based on this idea, there have been many attempts to modify the theory of GR, but, unfortunately, one of these efforts has yet to be universally accepted or studied extensively. Recently, there has been an increasing interest in studying modifications of GR with variable cosmological and Newtonian “constants” [26,27,28,29,30,31,32,33]. Apart from these investigations that include the variation in both G and $\mathsf{\Lambda }$ within the limit of GR, there are also studies on cosmological models with viscous fluids in the presence of G and $\mathsf{\Lambda }$ [17,34,35].

Many authors have studied the solutions of EFEs for homogeneous and anisotropic Bianchi type models [27,36,37,38,39,40,41,42,43,44,45,46]. For instance, Dwivedi [47] has solved the EFEs for a Bianchi Type-V model with variable cosmological and Newtonian “constants“ for a stiff perfect fluid. In his solution, the physical and kinematical parameters have been fully described, and he has found that the model has a singularity point, and all the physical parameters decrease as time increases, with the model isotropizing at late times. In addition to that, the universe described by such a model expands at a constant rate (i.e., the deceleration parameter equals zero). On the other hand, Yadav [48] has solved the same problem and generated two types of solution for the average scale factor, namely power-law type and exponential-type solutions. He found that the “cosmological constant” decreases with time and it reaches a small positive value at late and early times, respectively, results supported by supernova Type Ia observations. It is now well established from a variety of studies [47] that all the previous solutions are for special cases where the values of $\alpha$ and $\beta$ are chosen from the beginning of the solution process, and no previous study has investigated the general case.

The main aim of this paper is to find the general solution of the EFEs for Bianchi type-V models for a stiff perfect fluid with variable $\mathsf{\Lambda }$ and G without making any constraints on the value of $\alpha$ and $\beta$ in the $\mathsf{\Lambda }$ term, and to describe the behavior of the physical and kinematical parameters of the models. The remaining part of the paper is structured as follows: Section 2 is dedicated for a short description of the Bianchi type-V model and their general mathematical solutions for the EFEs. In Section 3, we present our models for restricted values of either $\alpha$ or $\beta$, and considering the other parameter as free. Finally, we conclude our results in Section 4.

We follow the Misner–Wheeler–Thorne notation where $x^a=(-x^0,x^1,x^2,x^3)$, Latin letters for the metric and geometric units with $c=1$.

2. Bianchi Type-V Cosmology

We consider the spatially homogeneous and anisotropic Bianchi type-V space-time that is represented by the following line-element:

$$d{s}^2=-d{t}^2+A^2(t)d{x}^2+e^{2x}\left[B^2(t)d{y}^2+C^2(t)d{z}^2\right],$$

(1)

where $A(t),B(t)$ and $C(t)$ are the components of the fundamental metric tensor. We assume that the cosmic matter is a perfect fluid that is represented by the following energy-momentum tensor:

$$T_{ij}=(p+\rho )u_i{u}_j+p{g}_{ij},$$

(2)

where $\rho$ is matter density, $u^i={\delta }_t^i=(-1,0,0,0)$ is the normalized fluid four-velocity, which is a time-like quantity such that $u^i{u}_i=-1$, and p is the fluid’s isotropic pressure. $\rho$ and p are related through the barotropic equation of state

$$p=w\rho ,0\le w\le 1,$$

(3)

where w is the equation-of-state (EoS) parameter. The EFEs with time-dependent $\mathsf{\Lambda }$ and G are given by

$$R_{ij}-\frac{1}{2}{g}_{ij}R=-8\pi G(t)T_{ij}+g_{ij}\mathsf{\Lambda }(t),$$

(4)

where $R_{ij}$ and $g_{ij}$ are the Ricci and metric tensors, respectively, and R is the Ricci scalar. Substituting Equations (1) and (2) into Equation (4), we get the field equations

$$\begin{array}{cc}\hfill \frac{\dot{A}\dot{B}}{AB}+\frac{\dot{B}\dot{C}}{BC}+\frac{\dot{A}\dot{C}}{AC}-\frac{3}{A^2}& =8\pi G(t)\rho +\mathsf{\Lambda }(t),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{\ddot{B}}{B}+\frac{\ddot{C}}{C}+\frac{\dot{B}\dot{C}}{BC}-\frac{1}{A^2}& =-8\pi G(t)p+\mathsf{\Lambda }(t),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\ddot{C}}{C}+\frac{\ddot{A}}{A}+\frac{\dot{C}\dot{A}}{CA}-\frac{1}{A^2}& =-8\pi G(t)p+\mathsf{\Lambda }(t),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{\ddot{A}}{A}+\frac{\ddot{B}}{B}+\frac{\dot{A}\dot{B}}{AB}-\frac{1}{A^2}& =-8\pi G(t)p+\mathsf{\Lambda }(t),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \frac{2\dot{A}}{A}-\frac{\dot{B}}{B}-\frac{\dot{C}}{C}& =0,\hfill \end{array}$$

(9)

with the dot over a letter representing differentiation with respect to time. The covariant divegence of the left hand side of Equation (4) produces

$$8\pi G\left[\dot{\rho }+(\rho +p)\left(\frac{\dot{A}}{A}+\frac{\dot{B}}{B}+\frac{\dot{C}}{C}\right)\right]+8\pi \dot{G}+\dot{\mathsf{\Lambda }}=0,$$

(10)

while the conservation of the usual energy-momentum tensor $T^{ij}$ (i.e., ${\nabla }_j{T}^{ij}=0$) yields

$$\dot{\rho }+(\rho +p)\left(\frac{\dot{A}}{A}+\frac{\dot{B}}{B}+\frac{\dot{C}}{C}\right)=0.$$

(11)

Substituting Equation (11) into (10), we get

$$8\pi \rho \dot{G}+\dot{\mathsf{\Lambda }}=0.$$

(12)

This equation shows how $\mathsf{\Lambda }$ and G evolve with time. It also shows that the two “constants” form a coupled system and, therefore, do not evolve independently of each other. The average scale factor $a=a(t)$ for Bianchi-V models is defined to be

$$a={(ABC)}^{1/3},$$

(13)

and the generalized Hubble parameter H is defined as in [47]

$$H=\frac{\dot{a}}{a}=\frac{1}{3}\left(\frac{\dot{A}}{A}+\frac{\dot{B}}{B}+\frac{\dot{C}}{C}\right)=\frac{1}{3}(H_1+H_2+H_3),$$

(14)

where $H_1,H_2$ and $H_3$ are directional Hubble’s parameters along x, y and z directions, respectively. The deceleration parameter q follows the usual definition:

$$q=-\frac{\ddot{a}a}{{\dot{a}}^2}=-1-\frac{\dot{H}}{H^2}.$$

(15)

The volume expansion parameter $\theta$, the average anisotropy parameter $A_p$, and shear modules $\sigma$ are defined as in [47]

$$\begin{array}{c}\theta ={\nabla }_i{u}^i=3H,\hfill \end{array}$$

(16)

$$A_p=\frac{1}{3}\sum _{i=1}^3{\left(\frac{H_i-H}{H}\right)}^2,$$

(17)

$${\sigma }^2=\frac{1}{2}{\sigma }_{ij}{\sigma }^{ij}=\frac{1}{2}\left(\frac{{\dot{A}}^2}{A^2}+\frac{{\dot{B}}^2}{B^2}+\frac{{\dot{C}}^2}{C^2}\right)-\frac{{\theta }^2}{6}$$

(18)

$$\begin{array}{c}=\frac{1}{3}\left(\frac{{\dot{A}}^2}{A^2}+\frac{{\dot{B}}^2}{B^2}+\frac{{\dot{C}}^2}{C^2}\right)-\frac{1}{3}\left(\frac{\dot{A}\dot{B}}{AB}+\frac{\dot{B}\dot{C}}{BC}+\frac{\dot{A}\dot{C}}{AC}\right),\hfill \end{array}$$

(19)

where the term ${\sigma }^{ij}$ represents the shear tensor. For this model, its scalar quantity comes out to be

$$\sigma =\frac{K}{a^3},$$

(20)

where K is a positive constant that is related to the anisotropy of the model. Having introduced these quantities, we can re-express the field Equations (5)–(8) and (11) in terms of a, H, q and $\sigma$ as

$$\begin{array}{cc}\hfill & 8\pi Gp-\mathsf{\Lambda }=(2q-1)H^2-{\sigma }^2+\frac{1}{A^2},\hfill \end{array}$$

(21)

$$\begin{array}{c}8\pi G\rho +\mathsf{\Lambda }=3H^2-{\sigma }^2-\frac{3}{A^2},\hfill \end{array}$$

(22)

so that, when subtracting Equation (21) from Equation (22), we eliminate the $\sigma$ term, obtaining

$$\frac{\ddot{a}}{a}+2\frac{{\dot{a}}^2}{a^2}-\frac{2}{A^2}=4\pi G(t)\rho (1-w)+\mathsf{\Lambda }(t).$$

(23)

We see that this equation cannot be integrated as it is because of the unknown functions A, G, $\rho$ and the parameter w. Now, integrating Equation (9) and absorbing the integration constant into A or B, gives

$$A^2=BC,$$

(24)

and substituting it into Equation (13), we see that the average scale factor reads

$$a=A.$$

(25)

Subtracting Equation (6) from (7) gives

$$\begin{array}{c}\hfill \frac{\ddot{B}}{B}+\frac{\dot{B}\dot{C}}{BC}=\frac{\ddot{A}}{A}+\frac{\dot{C}\dot{A}}{CA}⟶\frac{C\ddot{B}+\dot{B}\dot{C}}{BC}=\frac{C\ddot{A}+\dot{C}\dot{A}}{CA},\end{array}$$

(26)

which is equivalent to

$$\begin{array}{cc}\hfill \frac{1}{BC}\frac{d}{dt}\left(C\dot{B}\right)& =\frac{1}{AC}\frac{d}{dt}\left(C\dot{A}\right),\hfill \end{array}$$

(27)

and, multiplying this equation by $ABC$, we obtain

$$\begin{array}{cc}\hfill A\frac{d}{dt}\left(C\dot{B}\right)& =B\frac{d}{dt}\left(C\dot{A}\right).\hfill \end{array}$$

(28)

Both sides of this equation are a factorization of a product rule of the following

$$\begin{array}{cc}\hfill \frac{d\left(A(C\dot{B})\right)}{dt}-\dot{A}\dot{B}C& =\frac{d\left(B(C\dot{A})\right)}{dt}-\dot{A}\dot{B}C.\hfill \end{array}$$

(29)

Eliminating $\dot{A}\dot{B}C$ from both sides and integrating with respect to t yields

$$\begin{array}{cc}\hfill A(C\dot{B})& =B(C\dot{A})+k_1,\hfill \end{array}$$

(30)

and dividing throughout by $ABC$ yields

$$\begin{array}{cc}\hfill \frac{\dot{B}}{B}& =\frac{\dot{A}}{A}+\frac{k_1}{ABC},\hfill \end{array}$$

(31)

where $k_1$ is a constant of integration. Recalling that $A=a$, we finally get a first-order coupled differential equation of A and B as follows:

$$\begin{array}{cc}\hfill \frac{\dot{B}}{B}& =\frac{\dot{A}}{A}+\frac{k_1}{a^3}.\hfill \end{array}$$

(32)

Similarly, the subtraction of (6) from (8), and (7) from (8), respectively, gives

$$\begin{array}{cc}\hfill \frac{\dot{C}}{C}& =\frac{\dot{B}}{B}+\frac{k_3}{a^3},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \frac{\dot{C}}{C}& =\frac{\dot{A}}{A}+\frac{k_2}{a^3},\hfill \end{array}$$

(34)

where $k_2$ and $k_3$ are constants of integration. Again integrating the above set of Equations (32)–(33), we get

$$\begin{array}{cc}\hfill \frac{B}{A}& =d_1\exp \left(k_1\int \frac{dt}{a^3}\right),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \frac{C}{B}& =d_2\exp \left(k_2\int \frac{dt}{a^3}\right),\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \frac{C}{A}& =d_3\exp \left(k_3\int \frac{dt}{a^3}\right).\hfill \end{array}$$

(37)

These equations can be combined to give

$$\begin{array}{cc}\hfill B& =m_{1}a\exp \left(k\int \frac{dt}{a^3}\right),\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill C& =m_{2}a\exp \left(-k\int \frac{dt}{a^3}\right),\hfill \end{array}$$

(39)

where $m_1,m_2$ and k are constant values that depend on the undetermined constants $d_1,d_2$ and $d_3$, see [49] for more details. It is clear that Equations (25), (38) and (39) are totally dependent on the value of a, which is yet to be determined. To complete the solution process, we substitute Equation (25) into Equation (23) so that

$$\frac{\ddot{a}}{a}+2\frac{{\dot{a}}^2}{a^2}-\frac{2}{a^2}=4\pi G(t)\rho (1-w)+\mathsf{\Lambda }(t).$$

(40)

Now, following [24] for the parametrization of $\mathsf{\Lambda }$, we assume

$$\begin{array}{c}\hfill \mathsf{\Lambda }(t)=\frac{\alpha }{a^2}+\beta H^2.\end{array}$$

(41)

Setting $w=1$ for a stiff cosmological fluid and substituting Equation (41) into Equation (40) produces

$$\frac{\ddot{a}}{a}+(2-\beta )\frac{{\dot{a}}^2}{a^2}-\frac{1}{a^2}(2+\alpha )=0,$$

(42)

which can also be re-arranged as

$$a\ddot{a}+(2-\beta ){\dot{a}}^2-(\alpha +2)=0.$$

(43)

This is the generalized Friedman equation of Bianchi type-V models. The first term $a\ddot{a}$ represents force per unit mass times distance or work of a system that is equivalent to a potential energy, the second term $(2-\beta ){\dot{a}}^2$ is a kinetic energy and $(2-\beta )$ is a kind of mass, and the last term is the total energy of the system.

As we have mentioned that $\alpha$ and $\beta$ are adjustable constants of quantum field theory in an expanding curved background, the challenging problem here is to solve Equation (43) for general values of $\alpha$ and $\beta$ without making any constraints from the beginning as some authors did [47,48]. It is obvious that, for $\alpha \ne -2,\beta \ne 2$, Equation (43) is a nonlinear second-order differential equation, which is not easy to solve unless we transform it into a set of first-order equations. Before we move to the general solution, however, we remark that, for the special cases of $\alpha =-2,\beta =2$, the model has a scale factor solution that grows linearly with cosmic time, i.e.,

$$a(t)=C_{1}t+a_0,$$

(44)

where $C_1$ and $a_0$ are integration constants. This solution has a constant expansion and an initially increasing $\mathsf{\Lambda }$ that asymptotes to a constant value at late times, as well as a decreasing G, as can be seen in Figure 1. Moreover, we can see from Figure 2 that the model describes an expanding universe with an overall increasing volume and asymptotically approaches isotropy at late times.

On the other hand, for $\beta \ne 2,3$ and $\alpha =-2$, the model provides a solution of the form

$$a(t)={\left[(3-\beta )\left(C_{2}t+C_3\right)\right]}^{\frac{1}{3-\beta }}$$

(45)

for integration constants $C_2$ and $C_3$. Unsurprisingly, this solution reduces to the linear expansion solution above when we fix $\beta =2$. We clearly notice that the behaviours of $\mathsf{\Lambda }$ and G depend on the choice of $\beta$. For example, while fixing $\beta =1$ results in Figure 3, a different behaviour is observed when $\beta =-0.5$, as shown in Figure 4 where $\mathsf{\Lambda }$ asymptotically attains a constant value. It is worthwhile noticing, however, that the general behaviour of V and $\sigma$ remains the same as shown in Figure 5 (infinitely increasing and asymptotically vanishing at time infinity, respectively).

Now, let us introduce an intermediate variable P such that

$$P=\dot{a}\Rightarrow \ddot{a}=\frac{dP}{da}\frac{da}{dt}=P\frac{dP}{da},$$

(46)

and use this in Equation (43) to obtain

$$aP\frac{dP}{da}+(2-\beta )P^2-(\alpha +2)=0.$$

(47)

Separating the variables, we can rewrite this as

$$\frac{PdP}{(\beta -2)P^2+(\alpha +2)}=\frac{da}{a},$$

(48)

and integrating both sides gives

$$\begin{array}{cc}\hfill \ln \left[\frac{(\beta -2)}{\gamma }{P}^2+\frac{(\alpha +2)}{\gamma }\right]& =\ln \left[a^{2(\beta -2)}\right],\hfill \end{array}$$

(49)

where $\gamma$ is an integration constant. Substituting the value of P from (46) and taking the anti-logarithm for both sides gives

$$\begin{array}{c}\hfill {\dot{a}}^2-\frac{\gamma }{(\beta -2)}{a}^{2(\beta -2)}=\left(\frac{\alpha +2}{2-\beta }\right).\end{array}$$

(50)

Integrating Equation (50), we can see that

$$\int \frac{da}{\sqrt{\frac{\gamma }{(\beta -2)}{a}^{2(\beta -2)}+\left(\frac{\alpha +2}{2-\beta }\right)}}=\frac{a}{\sqrt{\frac{\alpha +2}{2-\beta }}}{}_2{F}_1\left(\frac{1}{2},\frac{1}{2(\beta -2)},1+\frac{1}{2(\beta -2)};\frac{\gamma }{(\alpha +2)}{a}^{2(\beta -2)}\right)=\pm \tau ,$$

(51)

which can be rewritten as

$$a{}_2{F}_1\left(\frac{1}{2},\frac{1}{2(\beta -2)},1+\frac{1}{2(\beta -2)};\frac{\gamma }{(\alpha +2)}{a}^{2(\beta -2)}\right)=\pm \sqrt{\frac{\alpha +2}{2-\beta }}\tau ,$$

(52)

where $\beta \ne 2$ and $\alpha \ne -2$, $\tau =t-t_0$ and ${}_2{F}_1$ is the hypergeometric function. It is obvious that Equation (52) is a general form of the solution, and cannot be simplified into an analytical expression for a as an explicit function of time unless we have to make some gauge choices on the values of $\alpha$ and $\beta$.

Our approach is different from the other previous works of some authors, where they firstly fixed the value of $\alpha$ and $\beta$ in the $\mathsf{\Lambda }$ term, then solved the generalized Friedmann equation, but we do not assume a priori the value or functional dependence of these constants or any other parameters in our models.

3. Some Specific Models of the General Solution

In this section, we will present our cosmological models that emerge from the choice of suitable values of $\alpha$ and $\beta$.

3.1. Model I: $\beta =1$

In this model, we choose $\beta =1$ such that Equation (52) gives a power solution of a as

$$\begin{array}{c}\hfill a=\sqrt{\alpha +2}{\left[{\tau }^2+\frac{\gamma }{{(\alpha +2)}^2}\right]}^{\frac{1}{2}}.\end{array}$$

(53)

Having known the value of a in this model, then the kinematical and dynamical parameters $H,\theta ,\sigma ,q$, and $\rho$, the Newtonian gravitational “constant” G, and the metric variables $A,B,C$ are calculated as

$$\begin{array}{cc}\hfill H& =\frac{\dot{a}}{a}=\frac{\tau }{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}},\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill \theta & =3H=\frac{3\tau }{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}},\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill \sigma & =\frac{K}{a^3}=\frac{K}{{\left[\sqrt{\alpha +2}\sqrt{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}}\right]}^3},\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill \rho & =\frac{{\rho }_0}{a^6}=\frac{{\rho }_0}{{\left[(\alpha +2){\tau }^2+{\displaystyle \frac{\gamma }{\alpha +2}}\right]}^3},\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill G& =-\int \frac{\dot{\mathsf{\Lambda }}}{8\pi \rho }dt=\frac{{\alpha }^3+5{\alpha }^2+8\alpha +4}{8\pi {\rho }_0}{\tau }^4-\frac{\gamma }{4\pi {\rho }_0}{\tau }^2+G_0,\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill q& =-\frac{\ddot{a}a}{{\dot{a}}^2}=-\frac{\gamma }{{(\alpha +2)}^2{\tau }^2},\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill A& =\sqrt{\alpha +2}\sqrt{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}},\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill B& =B_0\sqrt{\alpha +2}\sqrt{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}}\exp \left[\frac{k}{\sqrt{\alpha +2}}\frac{\tau }{\sqrt{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}}}\right],\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill C& =C_0\sqrt{\alpha +2}\sqrt{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}}\exp \left[-\frac{k}{\sqrt{\alpha +2}}\frac{\tau }{\sqrt{{\tau }^2+{\displaystyle \frac{\gamma }{{(\alpha +2)}^2}}}}\right],\hfill \end{array}$$

(62)

where $B_0=m_1{e}^{\mathrm{const}},C_0=m_2{e}^{\mathrm{const}}$ and $G_0$ are constants of integration. This model has an initial point of singularity at $t=t_0$ and $\gamma =0$. Furthermore, although early-time behaviour, shown in Figure 6, Figure 7, Figure 8 and Figure 9, differs slightly for different choices of $\alpha$ and $\gamma$, it can be shown that, as $t\to \infty$, $H,\theta ,\sigma ,\mathsf{\Lambda }$, the magnitude of q and $\rho$ all decrease with time, whereas G, the volume V and and the metric variables $A,B$ and C are increasing functions of time. It is worth mentioning here that although the numerical values of $\mathsf{\Lambda }$ and q decrease asymptotically towards constant values, whether the model describes an accelerated or decelerated expansion solely depends on the sign of $\gamma$. For example, a positive choice of $\gamma$ describes an early accelerated expansion that eventually slows down to an asymptotically constant expansion expansion, whereas a negative $\gamma$ describes an early decelerated expansion that eventually asymptotes to a constant expansion at late times. In addition, $|\sigma /\theta |\to 0$ as $t\to \infty$, thus indicating that the model approaches isotropy for large values of t, as observed in [47] as well.

In redshift space, we can show that the deceleration parameter for the model can be given by

$$q(z)=\frac{\gamma (\alpha +2){(1+z)}^2}{\gamma {(1+z)}^2-(\alpha +2)},$$

(63)

where $1+z\equiv \frac{a_0}{a}$. In Figure 7, we show in redshift space the transition from an early deceleration epoch (at large redshifts) to late-time acceleration (at small redshifts).

Here, a mathematical singularity occurs at the value of z for which $\gamma {(1+z)}^2-(\alpha +2)=0$, and hence the singular point can be shifted either way by choosing appropriate values of $\alpha$ and $\gamma$.

3.2. Model II: The Emden–Fowler Approach

This model naturally arises when comparing the nonlinear Equation (50) to the first integral of the Emden–Fowler equation [50]

$$\ddot{a}=r{t}^n{a}^m,$$

(64)

which, for $n=0$, is given by

$$\begin{array}{c}\hfill {\dot{a}}^2-\frac{2r}{m+1}{a}^{m+1}=s,\end{array}$$

(65)

where $m,n,r$ and s are constants. This equation has a particular solution of the form

$$\begin{array}{c}\hfill a={\left[\frac{2(m+1)}{r{(m-1)}^2}\right]}^{\frac{1}{m-1}}{t}^{\frac{2}{1-m}}.\end{array}$$

(66)

Thus, comparing Equations (65) and (50), we find that $m=2\beta -5$ and $r=\gamma$, and a reads

$$\begin{array}{c}\hfill a={\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{1}{3-\beta }}.\end{array}$$

(67)

Although this is a particular solution, this still is a direct and exact generalization result of all work that has been done on Bianchi-V cosmological models with varying $\mathsf{\Lambda }$ and G. Having known the values of a, the dynamical and kinematical parameters of the model can be computed as follows. Using Equation (67), the Hubble parameter is

$$\begin{array}{cc}\hfill H& =\frac{\dot{a}}{a}=\frac{1}{(3-\beta )}\frac{1}{t}.\hfill \end{array}$$

(68)

Integrating Equation (11) gives an expression for the density as

$$\begin{array}{cc}\hfill \rho & =\frac{{\rho }_0}{a^6}=\frac{{\rho }_0}{{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{6}{3-\beta }}}.\hfill \end{array}$$

(69)

Using Equation (68) in Equation (19), the expansion scalar is

$$\begin{array}{cc}\hfill \theta & =3H=\frac{3}{(3-\beta )}\frac{1}{t},\hfill \end{array}$$

(70)

and, from Equations (67) and (16), the shear scalar reads

$$\begin{array}{cc}\hfill \sigma & =\frac{K}{a^3}=\frac{K}{{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{3}{3-\beta }}}.\hfill \end{array}$$

(71)

From Equation (67), we can compute $\dot{a}$ and $\ddot{a}$ and substitute them into Equation (15) to compute the constant deceleration parameter as

$$\begin{array}{cc}\hfill q& =-\frac{\ddot{a}a}{\dot{a^2}}=2-\beta .\hfill \end{array}$$

(72)

Thus, this model can have $q>0$ (deceleration) if $\beta <2$, $q<0$ (acceleration) if $\beta >2$ and $q=0$ (constant expansion) if $\beta =2$.

The direct substitution of a in Equation (41) gives

$$\begin{array}{cc}\hfill \mathsf{\Lambda }& =\alpha {\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{2}{\beta -3}}+\frac{\beta }{{(3-\beta )}^2}\frac{1}{t^2},\hfill \end{array}$$

(73)

and differentiating Equation (73) with respect to t gives

$$\begin{array}{cc}\hfill \dot{\mathsf{\Lambda }}& =-2\alpha \sqrt{\frac{\gamma }{\beta -2}}{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{5-\beta }{\beta -3}}-\frac{2\beta }{{(3-\beta )}^2}\frac{1}{t^3}.\hfill \end{array}$$

(74)

We use this result together with Equation (69) to compute the Newtonian “constant" as

$$G=-\int \frac{\dot{\mathsf{\Lambda }}}{8\pi \rho }dt=G_0-\frac{{\gamma }^2{(\beta -3)}^4}{8\pi {\rho }_0{(\beta -2)}^3}\left[\gamma (\beta -3)t^{\frac{2\beta }{3-\beta }}-\frac{\alpha (\beta -2)}{2}{t}^{\frac{4}{3-\beta }}\right],$$

(75)

and, from Equations (25), (38) and (39), the metric variables can be computed as

$$\begin{array}{cc}\hfill A& ={\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{1}{3-\beta }},\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill B& =B_0{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{1}{3-\beta }}\exp \left\{\frac{k(\beta -3)}{3}{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}\right]}^{\frac{3}{\beta -3}}{t}^{\frac{\beta }{3-\beta }}\right\},\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill C& =C_0{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}t\right]}^{\frac{1}{3-\beta }}\exp \left\{-\frac{k(\beta -3)}{3}{\left[(\beta -3)\sqrt{\frac{\gamma }{\beta -2}}\right]}^{\frac{3}{\beta -3}}{t}^{\frac{\beta }{3-\beta }}\right\},\hfill \end{array}$$

(78)

where $B_0=m_1{e}^{\mathrm{const}},C_0=m_2{e}^{\mathrm{const}},{\rho }_0$ and $G_0$ are constants of integration whose actual values can be chosen using appropriate initial conditions.

The above results show that the model has a point of initial singularity or a Big Bang at $t=0$. This result improves on the result obtained in [47], where the singularity occurs at a negative point in time. As Figure 10, Figure 11 and Figure 12 depict, this model represents a contracting universe solution where V and both $\mathsf{\Lambda }$ and G simultaneously decrease, while H, $\theta$, $\sigma$ and $\rho$ all increase as time increases, with H asymptotically approaching zero from below. It is worth mentioning here that the model produces a constant deceleration of the expansion for acceptable values of $\beta$. In this model, as the expansion rate $\theta$ of the universe slows down as time increases, with the expansion eventually stopping, the ratio of $|\sigma /\theta |\to \infty$ as $t\to \infty$, thus predicting that the universe in this model becomes anisotropic at late times. These anisotropic predictions contradict the results reported by Dwivedi and Tiwari [47]. Moreover, as $t\to \infty$, the metric components A, B and C approach zero rapidly, and the average scale factor a tends to zero, which, with an even more rapidly increasing energy density, potentially results in a Big Crunch.

4. Conclusions

The main objective of this paper was to find general exact solutions for Bianchi Type-V cosmological models for a stiff perfect fluid with time-varying cosmological and gravitational “constants”. In practice, this is a generalization process of previous works without prior choice of the quantum field theoretically adjustable parameters $\alpha$ and $\beta$ that define the cosmological constant as $\mathsf{\Lambda }=\alpha /a^2+\beta H^2$.

In our solution process, we have followed two steps: (i) the EFEs were derived and reduced to a coupled system of differential equations, plus arbitrary constants of integration. This provides a way to integrate for the values of A, B, C, $\rho$, $\mathsf{\Lambda }$ and G, but it is not a complete solution because of the unknown value of a; (ii) The EFEs have been transformed to a nonlinear second-order DE in a, whose solution has been obtained using the hypergeometric function. Throughout the solution steps, we did not assume a priori the values of $\alpha$ and $\beta$, unlike previous works in the literature. We have therefore tried to solve the system in a more general setting. We have shown that the hypergeometric function is controlled by the value of $\beta$ and makes it very difficult to have a closed form for a as an explicit function of t.

Finally, two cosmological models are obtained with a choice of suitably fixed values of $\beta$ and through the transformation of the generalized Friedman equation into a special case of the Emden–Fowler equation. The dynamical and kinematical parameters of each model are exactly computed and clearly discussed. It has been shown that, while one of these models results in a universe that asymptotically isotropizes at late times, the other becomes increasingly anisotropic.

As more precise data become available, it will, in principle, be possible to constrain the integration constants that we chose arbitrarily in this study to get a better picture of this class of cosmological models. We leave this particular task to a future consideration.