1 Introduction

Cosmological and astrophysical data obtained thanks to the Supernovae Ia (SNeIa), the Cosmic Microwave Background (CMB) radiation anisotropies, the Large Scale Structure (LSS) and X-ray experiments provide strong evidence supporting a phase of accelerated expansion of the present day universe [1–13].

In order to find a suitable model for our universe, some possible reasons for this observed accelerated expansion have been investigated. Three main classes of models are usually proposed to describe this phenomenon: (a) a cosmological constant \(\Lambda \); (b) Dark Energy (DE) models; (c) modified theories of gravity.

The cosmological constant \(\Lambda \), which solves the equation of state (EoS) parameter \(\omega = -1\), represents the simplest candidate proposed to explain the accelerated expansion of the universe. However, \(\Lambda \) is well known to be related to the fine-tuning and the cosmic coincidence problems [14, 15]. According to the first, the vacuum energy density is about \(10^{123}\) times smaller than what is observed. Instead, according to the cosmic coincidence problem, the vacuum energy and Dark Matter (DM) models are nearly equal nowadays, although they have evolved independently and from different mass scales (which is a particular fortuity if no interaction exists between them). Many attempts have been made with the aim to propose a solution to the coincidence problem [16–19].

DE may explain the origin of the accelerated expansion of the universe [20, 21]. There are several models to describe the dark energy such as the models based on quintessence [22–24]. K-essence [25] and tachyonic models [26] are also two among many other ways to describe dark energy.

A successful model to describe DE is based on the Chaplygin gas (CG) equation of state [26, 27] and yields the Generalized Chaplygin Gas (GCG) model [28, 29]. It initially emerged in cosmology from string theory point of view [30–32]. One can indeed unify dark matter and dark energy using this model. It is also possible to study the effect of viscosity in GCG [33–35]. However, observational data ruled out such a proposal and the modified Chaplygin gas (MCG) model introduced [36]. Recently, viscous MCG has also been suggested and studied [37, 38]. A further extension of CG model is called the modified cosmic Chaplygin gas (MCCG), which was recently proposed [39–42]. Moreover, various Chaplygin gas models were studied from the holography point of view [43–45].

The MCG equation of state has two parts, the first term gives an ordinary fluid obeying a linear barotropic EoS, while the second term relates the pressure to some power of the inverse of energy density. Therefore, we are dealing with a two-fluid model. However, it is possible to consider a barotropic fluid with quadratic EoS or even with higher orders EoS [46, 47]. Therefore, it is interesting to extend the MCG EoS which recovers at least the barotropic fluid with quadratic EoS. We called the new version the extended Chaplygin gas (ECG) model [48, 49]. In order to get better agreement with observational data one can develop interesting models by varying the constants in EoS parameter.

The cosmic acceleration has also been accurately studied by imposing the concept of modification of gravity [50]. This new model of gravity (predicted by string/M theory) gives a very natural gravitational alternative for exotic matter. The explanation of the phantom, non-phantom, and quintom phases of the universe can be well described using modified gravity without the necessity of the introduction of a negative kinetic term in DE models. The cosmic acceleration is evinced by the straightforward fact that terms like \(1/R\) might become fundamental at small curvature. Furthermore, modified gravity models provide a natural way to join the early-time inflation and late-time acceleration. Such theories are also prime candidates for the explanation of the DE and DM, including for instance the anomalous galaxies rotation curves. The effective DE dominance may be assisted by the modification of gravity. Hence, the coincidence problem is solved in this case simply by the fact that the universe expands.

Modified gravity is also expected to be useful in high energy physics, explaining the hierarchy problem or unification of other forces with gravity [50]. Some of the most famous models of modified gravity are represented by braneworld models, \(f(T)\) gravity (where \(T\) indicates the torsion scalar), \(f(R)\) gravity (where \(R\) indicates the Ricci scalar curvature), \(f(G)\) gravity where

$$\begin{aligned} G=R^2-4R_{\mu \nu }R^{\mu \nu } + R_{\mu \nu \lambda \sigma }R^{\mu \nu \lambda \sigma } \end{aligned}$$

is the Gauss–Bonnet invariant, with \(R_{\mu \nu }\) representing the Ricci curvature tensor and \(R_{\mu \nu \lambda \sigma }\) representing the Riemann curvature tensor, \(f(R,T)\) gravity, DGP models, DBI models and Ho\({\check{\mathrm{r}}}\)ava–Lifshitz gravity [51–68].

In order to obtain a comprehensive model, we also add two modifications to the ordinary model. First, we consider a fluid which governs the background dynamics of the universe in a higher derivative theory of gravity. Second, we consider time-varying \(G\) and \(\Lambda \). As we know, the Einstein equations of general relativity do not permit any variations in the gravitational constant \(G\) and cosmological constant \(\Lambda \). Since the Einstein tensor has zero divergence and the energy conservation law is also zero, some modifications of Einstein equations are required. A similar study has been recently performed for another fluid model instead of the Chaplygin gas [69]. There are also several works on cosmological models with varying \(G\) and \(\Lambda \) [70, 71].

Therefore, in this paper, we study two different models of Chaplygin gas models in higher order gravity with varying \(G\) and \(\Lambda \).

This paper is organized as follows. In Sect. 2, we introduce our models. In Sect. 3, we study special cases corresponding to constant \(G\) and \(\Lambda \). In Sect. 4, we investigate the Hubble, the deceleration and the EoS parameters for the two models that were introduced. In Sect. 5, we consider the statefinder diagnostics for both models. In Sect. 6, we make a perturbation analysis for the Chaplygin gas. Finally, in Sect. 7, we draw the conclusions of this paper.

2 The models

We consider two different models for a fluid which govern the background dynamics of the universe in a higher derivative theory of gravity in the presence of time-varying \(G\) and \(\Lambda \). Within modified theories of gravity, we hope to solve the problems of the dark energy which originate with general relativity.

A gravitational action with higher order term in the Ricci scalar curvature \(R\) containing a variable gravitational constant \(G(t)\) is given by

$$\begin{aligned} I = - \int {\mathrm{d}^{4}x \sqrt{-g}\left[ \frac{1}{ 16 \pi G} f(R) +L_\mathrm{m} \right] }, \end{aligned}$$
(1)

where \(f(R)\) is a function of \(R\) and its higher power (including a variable \(\Lambda \)), \(g\) is the determinant of the four dimensional tensor metric \(g^{\mu \nu }\) and \(L_\mathrm{m}\) represents the matter Lagrangian. Considering second order gravity, we can take into account the following expression of \(f(R)\):

$$\begin{aligned} f(R)=R + \alpha R^{2}-2\Lambda , \end{aligned}$$
(2)

where \(\alpha \) is a constant parameter.

By using the following flat FRW metric:

$$\begin{aligned} \mathrm{d}s^2=-\mathrm{d}t^2+a(t)^2(\mathrm{d}r^{2}+r^{2}\mathrm{d}\Omega ^{2}), \end{aligned}$$
(3)

where \(\mathrm{d}\Omega ^{2}=\mathrm{d}\theta ^{2}+\sin ^{2}\theta \mathrm{d}\phi ^{2}\) represents the angular part of the metric and \(a(t)\) represents the scale factor (which gives information as regards the expansion of the universe), we get the following Friedmann equations [72]:

$$\begin{aligned}&H^{2}-6\alpha (2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2})=\frac{8\pi G}{3}\rho +\frac{\Lambda }{3},\end{aligned}$$
(4)
$$\begin{aligned}&\dot{H}+H^{2}-6\alpha (2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2})\nonumber \\&\quad =-\frac{4\pi G}{3}(3P+\rho )+\frac{\Lambda }{3} \end{aligned}$$
(5)

where an overdot and two overdots indicate, respectively, the first and the second derivative with respect to the cosmic time \(t\).

The energy conservation equation is given by the following relation:

$$\begin{aligned} \dot{\rho }+3H(\rho +P)=-\left( \frac{\dot{G}}{G}\rho +\frac{\dot{\Lambda }}{8\pi G}\right) , \end{aligned}$$
(6)

where \(\rho \) and \(P\) are the energy density and the pressure of the perfect fluid, respectively. Modified theories of gravity (like \(f(R)\) theories) give an opportunity to find a natural representation and introduction of the dark energy into the theory. Therefore, the type of the dark energy and dynamics of the universe depends on the form of \(f(R)\) which will be considered. The type of the work which we would like to consider in this paper assumes the existence of an effective fluid controlling the dynamics of the universe composed non-interacting dark energy (emerging from \(f(R)\)) and a fluid emerging from our assumptions.

Assuming there is not interaction, for matter we have the following continuity equation:

$$\begin{aligned} \dot{\rho }+3H(\rho +P)=0. \end{aligned}$$
(7)

Therefore, comparing Eqs. (6) and (7), we can easily derive the following relation for the dynamics of \(G\):

$$\begin{aligned} \dot{G} +\frac{\dot{\Lambda }}{8\pi \rho }=0. \end{aligned}$$
(8)

The energy density \(\rho \) can be assumed to originate by some kind of Chaplygin gas. In particular, we find that the MCG is described by the following EoS:

$$\begin{aligned} p=A\rho -\frac{B}{\rho ^{n}}, \end{aligned}$$
(9)

where \(A\) and \(B\) are two arbitrary constant parameters which may fitted using observational data. The special case corresponding to \(A = 0\) yields the GCG EoS, while the special case corresponding to \(A = 0\) and \(n = 1\) recovers the pure Chaplygin gas EoS. Moreover, the limiting case corresponding to \(n=0.5\) has been studied in [73]. The authors of Ref. [74] concluded that the best fitted parameters are \(A = -0.085\) and \(\alpha = 1.724\), while Constitution + CMB + BAO data suggests \(A = 0.061 \pm 0.079\), \(n = 0.053 \pm 0.089\), and Union + CMB + BAO results suggest \(A = 0.110 \pm 0.097\), \(n = 0.089 \pm 0.099\) [75]. Other observational constraints on the MCG model using Markov Chain Monte Carlo suggest \(A = 0.00189^{+0.00583}_{-0.00756}\), \(\alpha =0.1079^{+0.3397}_{-0.2539}\) at \(1\sigma \) level and \(A = 0.00189^{+0.00660}_{-0.00915}\), \(n=0.1079^{+0.4678}_{-0.2911}\) at \(2\sigma \) level [76]. There are also other constraints; for example those reported by Refs. [77–80].

It is possible to consider \(B\) as a variable (depending on the time \(t\)) instead of a constant. So, a time-varying MCG can be described by the following EoS:

$$\begin{aligned} P=A\rho -\frac{B(t)}{\rho ^{n}}, \end{aligned}$$
(10)

where

$$\begin{aligned} B(t)=\omega (t)a(t)^{-3(1-\omega (t))(1-n)}, \end{aligned}$$
(11)

with \(\omega (t)\) given by [81]:

$$\begin{aligned} \omega (t)=\omega _{0}+\omega _{1}t\frac{\dot{H}}{H}. \end{aligned}$$
(12)

where \(\omega _0\) and \(\omega _1\) are two constant parameters.

In the second model we consider ECG with the following EoS:

$$\begin{aligned} P=\sum _{k=1}^{m}{A_{k}\rho ^{k}}-\frac{B}{\rho ^{n}}, \end{aligned}$$
(13)

\(B\) and \(n\) are arbitrary constants, and \(A_{k}=1/k\) assumed in this paper. The ECG EoS reduces to the MCG EoS in the limiting case of \(k=1\), and it can recover the barotropic fluid with quadratic EoS by setting \(k=2\). Moreover, higher values of \(k\) may recover a higher order barotropic fluid, which is indeed our motivation to consider the ECG.

3 Numerical results with constant \(G\) and \(\Lambda \)

We start the analysis of the models considered in this paper with the case of constant \(\Lambda \) and \(G\) (we will also consider units of \(8\pi G=c=1\)). Therefore, we have the following two equations, which will describe the dynamics of the universe; in this case the dynamics of the Chaplygin gases described by the equations:

$$\begin{aligned} H^{2}-6\alpha (2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2})=\frac{\rho }{3}+\frac{\Lambda }{3}, \end{aligned}$$
(14)

and

$$\begin{aligned} \dot{\rho }+3H(\rho +P)=0. \end{aligned}$$
(15)

In the first model, given by the EoS (10), we obtain the behavior of the Hubble parameter and the EoS parameter \(\omega =P/\rho \) by using a numerical analysis which is represented in Fig. 1. The left plot of Fig. 1 shows a typical time evolution of the Hubble parameter. It is found that the value of \(H\) decreases with the increasing of the value of \(\Lambda \). For \(\Lambda <2.5\), the Hubble parameter is an increasing function of time. \(\Lambda =2.5\) yields the constant \(H\), while for \(\Lambda >2.5\) the Hubble parameter is a decreasing function of the time.

Fig. 1
figure 1

Behavior of \(H\) and \(\omega \) against \(t\)

Full size image

Since one would expect the Hubble parameter to decrease with time and become constant at the present epoch, this model is not in good agreement with observational data. However, the evolution of the EoS parameter (right plot of Fig. 1) agrees with the \(\Lambda \)CDM model where it is expected to behave as \(\omega \rightarrow -1\).

Using Fig. 1, we can obtain the following fit of the function of the Hubble parameter:

$$\begin{aligned} H=H_{0}+C(2.5-\Lambda )t, \end{aligned}$$
(16)

where \(H_{0}\) is the current value of the Hubble parameter and \(C\) is a constant. Using the expression of \(H\) obtained in Eq. (16), we can now investigate the deceleration parameter \(q\) via the following relation:

$$\begin{aligned} q=-1-\frac{\dot{H}}{H^{2}} = - 1 -\frac{C(2.5-\Lambda )}{[ H_{0}+C(2.5-\Lambda )t]^2}. \end{aligned}$$
(17)

Numerical analysis yields Fig. 2, which shows that, for appropriate values of the parameters, we obtain \(q=-1\), which is in agreement with the \(\Lambda \)CDM model. In Fig. 2, we also see deceleration to acceleration and acceleration to deceleration phase transitions. The curves of this figure were drawn for \(\Lambda \ge 2.5\); however, in the cases of \(\Lambda \le 2.5\) we find \(q\le -1\).

Fig. 2
figure 2

Behavior of \(q\) against \(t\). \(\Lambda =2.5\) (orange), \(\Lambda =3.5\) (red), \(\Lambda =4\) (green), \(\Lambda =5\) (black)

Full size image

The analysis of the second model given by Eq. (13) is based on the following EoS parameter:

$$\begin{aligned} \omega = \frac{P}{\rho } =\sum _{k=1}^{m}{A_{k}\rho ^{k-1}}-\frac{B}{\rho ^{n+1}}. \end{aligned}$$
(18)

In the plots of Fig. 3 we represent behavior of some cosmological parameters by numerical analysis up to third order (\(m = 3\)). We obtained similar results with the first model. Therefore, we should apply a modification to obtain results which are in agreement with the current observations. For this reason, we will consider time-varying \(G\) and \(\Lambda \) in the next section.

Fig. 3
figure 3

Behavior of \(H\), \(\omega \) and \(q\) against \(t\)

Full size image

4 Numerical results with varying \(G\) and \(\Lambda \)

To perform an analysis of the dynamics of the universe, we assume a specific model where effectively the form \(\Lambda \) is given by the following relation:

$$\begin{aligned} \Lambda =\gamma \rho , \end{aligned}$$
(19)

where \(\gamma \) is a positive constant. Therefore, using the expression given in Eq. (8), we have the following expression for the dynamics of \(G\):

$$\begin{aligned} \dot{G} +\frac{ \gamma \dot{\rho } }{8\pi \rho }=0. \end{aligned}$$
(20)

Equation (20) can easily be integrated leading to the following expression for \(G\):

$$\begin{aligned} G=G_{0}-\frac{\gamma }{8 \pi } \ln {\rho }, \end{aligned}$$
(21)

where \(G_{0}\) is an integration constant. In the following subsections, we give a numerical analysis of the two models.

4.1 Variable modified Chaplygin gas model

Within this subsection, we analyze our first model, which is the time-varying modified Chaplygin gas with EoS given by

$$\begin{aligned} \omega =A-\frac{B(t)}{\rho ^{n+1}}, \end{aligned}$$
(22)

where we used Eq. (10).

The plots of Fig. 4 show that the Hubble parameter in this model is a decreasing function of time The first plot shows that increasing \(\alpha \) increases the value of \(H\); therefore we find that higher order terms of gravity increase the value of the Hubble parameter. The second plot deals with the variation of the parameter \(A\). It has been shown that increasing \(A\) increases the value of \(H\). In the third plot we look at how the Hubble parameter changes with \(n\). It is clear that the variation of \(H\) depends on the time period. During the period when \(t < 5\), increasing \(n\) increases the value of \(H\); but for \(t > 5\), increasing \(n\) decreases the value of \(H\). Finally, in the last plot we can see the variation of \(H\) with \(\omega _0\) and \(\omega _1\). In the plots of Fig. 5 we see the evolution of the deceleration parameter with various values of \(\alpha \), \(A\), \(n\), \(\omega _0\), and \(\omega _1\). In all cases we see that \(q\) takes value between \(-1\) and 0, which is in agreement with current observational data. So in this case there is no acceleration to deceleration phase transition.

Fig. 4
figure 4

Behavior of \(H\) against \(t\). Model 1

Full size image
Fig. 5
figure 5

Behavior of \(q\) against \(t\). Model 1

Full size image

In Fig. 6 we draw the EoS parameter versus time. We see from the first plot that increasing \(\alpha \) decreases the value of \(\omega =P/\rho \). It is illustrated that higher values of \(\alpha \) make \(\omega \rightarrow -1\) faster than lower values of \(\alpha \). A similar situation happens by varying \(A\) (at least for the late time behavior). The third plot of Fig. 6 shows the variation of \(\omega \) with \(n\). The first and the second plots were drawn for \(n = 0.1\) and \(n = 0.3\), respectively, causing \(\omega \) to be a decreasing function of time and eventually asymptotically approach constant negative values. But, for different values of \(n\) and \(\omega \), the situation is different, as illustrated in the third and fourth plots of Fig. 6.

Fig. 6
figure 6

Behavior of \(\omega \) against \(t\). Model 1

Full size image

Finally, in Fig. 7 we draw the variation of \(\dot{G}/G\) versus time. According to the most observational data, we have [82]

$$\begin{aligned} \left| \frac{\dot{G}}{G}\right| \le 1.3\times 10^{-12}\,\mathrm{year}^{-1}. \end{aligned}$$
(23)

The plots of Fig. 7 show good agreement with the observational constraint given in Eq. (23).

Fig. 7
figure 7

Behavior of \(\dot{G}/G\) against \(t\) in unit of \(10^{-11}\,\mathrm{year}^{-1}\). Model 1

Full size image

Therefore, we can conclude that the first model agrees with observational data with the exception of the results of the deceleration parameter.

4.2 Extended Chaplygin gas model

In the second model we use the EoS parameter given in Eq. (18). The plots of Fig. 8 show the behavior of Hubble parameter as a function of time.

Fig. 8
figure 8

Behavior of \(H\) against \(t\). Model 2

Full size image

We can see that the Hubble parameter is a decreasing function of time. These plots suggest the following general behavior for \(H\):

$$\begin{aligned} H=H_{0}-Ct(1+t), \end{aligned}$$
(24)

where \(C\) is a positive constant. Using this fit for the Hubble parameter we can obtain the following expression for the deceleration parameter \(q\):

$$\begin{aligned} q= -1 + \frac{C(1+2t)}{[ H_{0}-Ct(1+t)]^2}. \end{aligned}$$
(25)

We plotted the behavior of the deceleration parameter obtained in Eq. (25) in Fig. 9. We can see that \(q\) is negative during the early universe as well as the late universe therefore the universe is accelerating during these two eras. Between these two epochs \(q\) is positive, hence the universe is decelerating. This behavior is in agreement with the \(\Lambda \)CDM model.

Fig. 9
figure 9

Behavior of \(q\) against \(t\). Model 2

Full size image

In Fig. 10, we plot the EoS parameter as a function of the time. The first plot shows that varying of \(\alpha \) does not have an effect on \(\omega \) in the presence of higher order terms. The second plot deals with the variation of the constant \(B\). For the cases of \(B < 3\) the EoS parameter is a decreasing function of time which tends to \(-1\), while for the case of \(B = 3\) the EoS parameter is an increasing function of time which also goes to the value of \(-1\). The third plot shows that increasing \(n\) decreases the value of \(\omega \). These plots are obtained for \(m = 3\). The last plot contains five curves corresponding to \(m = 1,\ldots , 5\) with increasing \(n\). The case of \(m = 1\), which corresponds to the modified Chaplygin gas, is illustrated by a violet line. In the plots of Fig. 11 we see the evolution of \(\dot{G}/G\) is in agreement with observational data.

Fig. 10
figure 10

Behavior of \(\omega \) against \(t\). Model 2

Full size image
Fig. 11
figure 11

Behavior of \(\dot{G}/G\) against \(t\). Model 1

Full size image

5 Statefinder diagnostics

In the framework of general relativity, dark energy can explain the present cosmic acceleration. Except the cosmological constant there are many other candidates of dark energy (quintom, quintessence, brane, modified gravity etc.). The property of dark energy is that it is model dependent and in order to differentiate different models of dark energy, a sensitive diagnostic tool is needed. The Hubble parameter H, and deceleration parameter \(q\) are very important quantities which can describe the geometric properties of the universe. Since \(\dot{a}>0\), \(H>0\) means that the universe expands. Also, \(ddot{a}>0\), which is \(q<0\), indicates the accelerated expansion of the universe. Since the various dark energy models give \(H>0\) and \(q<0\), one needs further evidence to differentiate general models of dark energy by investigating cosmological observational data more accurately. For this aim, we need higher order of time derivatives of the scale factor, a geometrical tool. In Ref. [83] the geometrical statefinder diagnostic tool was proposed, based on dimensionless parameters \((r, s)\) which are functions of the scale factor and its time derivative.

These parameters are defined as follows:

$$\begin{aligned} r=\frac{1}{H^{3}}\frac{\dddot{a}}{a} \quad s=\frac{r-1}{3\left( q-\frac{1}{2}4\right) }. \end{aligned}$$
(26)

Plots of Fig. 12 show results of our numerical analysis in both models for the cases of constant (\(G\), \(\Lambda \)), and varying (\(G\), \(\Lambda \)). We can see that the value of \(r\) in the cases of varying \(G\) and \(\Lambda \) are lower than the cases of constant \(G\) and \(\Lambda \). However, \((r, s)=(1, 0)\) is verified in all cases.

Fig. 12
figure 12

\(r\)–\(s\). Top panels correspond to the cases of constant \(G\) and \(\Lambda \) for model 1 and model 2, respectively. Bottom panels correspond to the cases of varying \(G\) and \(\Lambda \) for model 1 and model 2, respectively

Full size image

6 Perturbation analysis

One of the best ways to investigate the stability of a model is studying cosmological perturbations [84]. The first equation is obtained using the continuity equation given in Eq. (15):

$$\begin{aligned} \delta \dot{\rho }=-3\delta H(\rho +P)-3H(1+C_\mathrm{s}^{2})\delta \rho , \end{aligned}$$
(27)

where

$$\begin{aligned} C_\mathrm{s}^{2}=\frac{\delta P}{\delta \rho } \end{aligned}$$
(28)

is the squared speed of sound. On the other hand, using Eq. (4), we can obtain the following relation:

$$\begin{aligned} \delta \dot{H}&= -2H\delta H+6\alpha \delta (2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2}) \nonumber \\&+\frac{4\pi }{3}\left[ \frac{\gamma }{8\pi }\delta (1+3\omega )\rho -G(1+C_\mathrm{s}^{2})\delta \rho \right] +\frac{\gamma }{3}\delta \rho ,\nonumber \\ \end{aligned}$$
(29)

where we used Eqs. (19) and (21) together with \(P=\omega \rho \) and \(\delta =\frac{\delta \rho }{\rho }\). We can then rewrite Eq. (14) in first order in \(\delta H\), and differentiate it with respect to \(t\), obtaining

$$\begin{aligned} \delta \dot{H}&= \left[ 3\alpha \frac{2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2}}{H}-H+\frac{4\pi G}{3H}\rho \right] \ddot{\delta }\nonumber \\&+ \left[ 3\alpha \frac{\mathrm{d}}{\mathrm{d}t}(2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2})\right. \nonumber \\&\left. -\dot{H}+(1+C_\mathrm{s}^{2}+(1+\omega )\rho )(\gamma +4\pi G)\right] \dot{\delta }, \end{aligned}$$
(30)

where we used \(\delta =\frac{\delta H}{H}\) along with Eq. (27). We are interested in the case of the extended Chaplygin gas; so combining Eqs. (29) and (30) gives

$$\begin{aligned} X\ddot{\delta }+Y\dot{\delta }+Z\delta =0, \end{aligned}$$
(31)

where

$$\begin{aligned} X&\equiv \left( 5-2\frac{\dot{H}}{H^{2}}\right) \nonumber \\&\times \left[ 3\alpha \frac{2H\ddot{H}-\dot{H}^{2}+6\dot{H}H^{2}}{H}-H+\frac{4\pi G}{3H}\rho \right] ,\end{aligned}$$
(32)
$$\begin{aligned} Y&\equiv \left( 5-2\frac{\dot{H}}{H^{2}}\right) \left\{ 3\alpha \frac{\mathrm{d}}{\mathrm{d}t}(2H\ddot{H}-\dot{H}^{2} +6\dot{H}H^{2}) \right. \nonumber \\&\left. -\dot{H}+[1+C_\mathrm{s}^{2}+(1+\omega )\rho ](\gamma +4\pi G)\right\} ,\end{aligned}$$
(33)
$$\begin{aligned} Z&\equiv 12\dot{H}+2\frac{\ddot{H}}{H}\!+\!\frac{\gamma }{6\pi }(1+3\omega )-\frac{4\pi G}{3}(1\!+\!C_\mathrm{s}^{2}) \!+\!\frac{\gamma }{3}-2.\nonumber \\ \end{aligned}$$
(34)

Then, using Eqs. (4) and (24), we can obtain a time-dependent equation to investigate the perturbation evolution. We can see from Fig. 13 that perturbations will grow for \(m>1\). However, after long time, they lead to a constant. Moreover, the analysis of the squared speed of sound \(C_\mathrm{s}^{2}\) shows that the model is stable at all time. In Fig. 14 we have plotted the behavior of the squared speed of the sound for three different values of \(m\), in particular \(m=1\), \(m=2\), and \(m=3\).

Fig. 13
figure 13

Typical behavior of \(\delta \) in terms of \(t\) for the extended Chaplygin gas. \(m=1\) (solid), \(m=2\) (dotted), \(m=3\) (dashed)

Full size image
Fig. 14
figure 14

Typical behavior of squared sound speed in terms of \(t\) for the extended Chaplygin gas. a \(m=1\), b \(m=2\), c \(m=3\)

Full size image

In the case of \(m = 1\), which corresponds to MCG, we have constant squared sound speed at initial time increasing at the late time. For the case of \(m = 2\), i.e. a quadratic barotropic fluid, we see that the value of the squared sound speed increases dramatically at the early universe and is decreasing at later stages to a constant value. A similar behavior is observed for the case of \(m = 3\). Therefore, we can conclude that the squared sound speed is positive for these models, and hence stable.

7 Conclusion

In this work, we considered higher order \(f(R)\) gravity with time-dependent \(G\) and \(\Lambda \). We assumed the Chaplygin gas as a candidate for dark energy and considered two models to describe the evolution of our universe. The first one was the variable modified Chaplygin gas model (VMCG). The constant \(B\) in ordinary MCG is a time-dependent quantity in the context of VMCG. The second model was the extended Chaplygin gas. This is indeed an extended version of MCG recovering higher order barotropic fluid EoS.

We assumed the specific case where \(\Lambda \) is proportional to the energy density and analyzed the Hubble, the deceleration, and the EoS parameters. We found that varying \(G\) and \(\Lambda \) fit the observational data compared to the cases of constant \(G\) and \(\Lambda \).

We also found that the extended Chaplygin gas is a more appropriate model compared to the variable modified Chaplygin gas model. We can see also an acceleration to deceleration phase transition in the extended Chaplygin gas model.

Finally, we investigated the evolution of the density perturbations by looking at perturbed Friedmann equations. We also confirmed the stability of the extended Chaplygin gas by investigating the square of the speed of sound. In summary, we think that the higher order corrected extended Chaplygin gas with time-dependent \(G\) and \(\Lambda \) as a possible model enables one to describe the evolution of our universe.