1 Introduction

The amount of theoretical and observational results strongly suggest that our universe had undergone one accelerated phase of expansion during its very early phase of evolution, known as inflation [1, 2], and presently it is again undergoing an accelerating phase of expansion, known as late time cosmic acceleration, dominated either by some hypothetical dark energy fluid adjusted in the Einstein’s General Relativity (GR) [3] or by some extra geometrical terms appearing due to modifications of GR or due to new gravitational theories [4, 5]. So, one accelerating phase (inflation) occurred in the high energy regime and the second acceleration is occurring in the low energy regime. The difference of the energy scales between the early and late accelerating phases of the universe is significantly high. Therefore, understanding the evolution of the universe in both these regimes has been one of the longstanding issues in modern cosmology. Despite of many investigations performed in these sectors, the actual nature is still unrevealed. The unification of the gravitational theory both at classical and quantum levels has therefore remained as an open problem in modern cosmology.

Concerning the early evolution of the universe, considerable attention on the nature of the gravitational theory at the quantum level has been paid by many investigators. The existing literature demands that the theory of quantum gravity [6,7,8,9,10,11,12,13,14] played a very crucial role to understand the very early evolution of the universe. However, the approach to reveal the quantum nature of the gravity is not unique. Various approaches were proposed over the years that include string theory (ST) [15], doubly special relativity (DSR) [16,17,18,19,20,21,22] , black hole physics (BHP) [23,24,25,26] etc. All of them reported the existence of a minimum length scale of the order of the Planck length \(l_{\mathrm {pl}}\) (equivalently, there exists a maximum energy scale in nature). This actually motivated to generalize the Heisenberg’s uncertainty principle to some generalized uncertainty principle (GUP) in quantum gravity [27].

The generalized uncertainty principle has been found to be very effective to explain several issues related to the dynamics of the universe. For instance, the origin of the magnetic fields with microgauss strength [28, 29] which is observed in the intergalactic regimes of the universe acquires some explanation if GUP is taken into account [30]. The potential application of GUP can be found in the context of black hole mechanics. Following the Bekenstein’s entropy relation [23,24,25,26] and the Hawking temperature [31, 32], the small black holes in the universe should emit black body radiation. Due to this radiation, these small black holes must loss their mass and consequently becomes hotter. This emission of radiation continues until they are completely evaporated. However, the above formulation was based under the assumption of a classical black hole metric together with the consideration that the emitting radiation of the black hole was much smaller (can be ignored) compared to the rest mass energy of the black hole. Now, during this emission process, the black holes naturally become smaller and lighter, and certainly within this length scale, the assumption of classical metric and the ignorance of the radiation seem to be invalid. Thus, whether the black holes will completely evaporate to either photons, ordinary matter particles or vacuum, or whatever it be (some remnant for instance) is highly questionable due to the lack of a definite quantum gravitational theory (see the discussions [33,34,35,36,37]). Here, GUP finds an answer towards this discrepancy [38]. As argued in Ref. [38] a small black hole having temperature greater than the ambient temperature will keep radiating photons and ordinary particles until it reaches the Planck scale. When the black hole will reach the Planck size, it will stop radiating and its entropy reaches zero while its effective temperature will reach the maximum limit, and finally, it reduces to an inert remnant having no radiating power, but only gravitational interactions may exist. It is not necessary that the created remnant will enjoy the horizon structure as in classical black hole [35]. Such remnants might be consider as the dark matter candidate [39]. This is a very surprising result because GUP indicates the origin of one of the heavy fluids of the universe. Naturally, the effects of GUP in various cosmological theories can be studied in order to see whether such effects give rise to some interesting directions related the dynamical history of the universe. Since dark energy is another heavy resource of the universe playing the leading role behind its accelerated expansion, so how GUP modified dark energy models behave, will be an interesting direction of research. As the minimum length effects are quite universal at any stage, so, there is no reason to exclude the GUP modified models, rather, one can generalize such cosmic theories. Thus, in this work shall focus on the cosmological models that are modified by the GUP. The work has been organized in the following way.

In Sect. 2 provide with a brief description on the GUP and the corresponding algebra. In Sect. 3 we discuss the quintessence scalar field model modified by the GUP. Then in Sect. 4 we discuss the gravitational field equations of the modified quintessence scalar field in the Friedmann–Lemaître–Robertson–Walker universe. In Sect. 5 we define the dimensionless dynamical variables and perform the stability analysis. After that in Sect. 6 we introduce the kinematical quantities and compare with the present model. Finally, in Sect. 7 we close the work with a brief summary of all the findings.

2 Generalized uncertainty principle

The GUP is related with the existence of a minimum measurable length, such that the Heisenberg uncertainty has to be modified as

$$\begin{aligned} \Delta X_{i}\Delta P_{j}\geqslant \frac{\hbar }{2}[\delta _{ij}(1+\beta P^{2})+2\beta P_{i}P_{j}]. \end{aligned}$$
(1)

where \(\beta \) is defined as \(\beta ={\beta }_{0}\ell _{Pl}^{2}/2\hbar ^{2}~\) and \(\beta _{0}~\) is the deformation parameter, which is also known as the quadratic generalized GUP. From (1) the generalized deformation of the Heisenberg algebra in a four-dimensional Minkowski spacetime of signature [40,41,42,43] follows

$$\begin{aligned}{}[X_{\mu },P_{\nu }]=-i\hbar [(1-\beta (\eta ^{\alpha \gamma }P_{\alpha }P_{\gamma }))\eta _{\mu \nu }-2\beta P_{\mu }P_{\nu }]. \end{aligned}$$
(2)

From the algebra (2) it leads to the deformation of the coordinate representation of the momentum operator [44]. The deformation parameter is usually defined to be positive but in the literature it has been defined also as a negative value parameter [45, 46]. Some theoretical constraints on the deformation parameter have been proposed before in [47, 48]. Specifically, in [47] the authors proposed the computation of the deformation parameter by using the leading quantum correction of the Newtonian potential. On the other hand, in [48] it has been proposed a theoretical calculation of the deformation parameter from the maximal acceleration theory. Furthermore, in [49] it has been found that the deformation parameter may be constrained by some phenomenological aspects of the astrophysical objects such as the binary systems and gravitational wave emissions.

By using the latter commutation relation we can write the coordinate representation of the operators \(X_{\mu },P_{\nu }\). We select \(X_{\mu }\) underformed, that is \(X_{\mu }=x_{\mu }\), and the momentum operator to modified as [50]

$$\begin{aligned} P_{\mu }=p_{\mu }(1-\beta (\eta ^{\alpha \gamma }p_{\alpha }p_{\gamma } ))~, \end{aligned}$$
(3)

in which \(p^{\mu }\) are defined as usual, that is, \(p^{\mu }=i\hbar \frac{\partial }{\partial x_{\mu }},\) and \([x_{\mu },p_{\nu }]=-i\hbar \eta _{\mu \nu }\). 

In the context of the GUP, the Klein–Gordon equation which describes a spin-0 particle, described by the wave function \(\Psi ,\) as \(\eta ^{\mu \nu }P_{\mu }P_{\nu }\Psi -\left( mc\right) ^{2}\Psi =0\), it is now a fourth-order partial differential equation as given

$$\begin{aligned} \Delta \Psi +2\beta \hbar \Delta \left( \Delta \Psi \right) +\left( \frac{mc}{\hbar }\right) ^{2}\Psi +O\left( \beta ^{2}\right) =0, \end{aligned}$$
(4)

where \(\Delta \) is the Laplace operator defined as \(\Delta =\square \) when the underlying space is \(\eta _{\mu \nu }\), or in general \(\Delta =\frac{1}{\sqrt{-g} }\partial _{\mu }\left( g^{\mu \nu }\sqrt{-g}\partial _{\mu }\right) \) for any Riemannian space of Lorentzian signature with metric \(g_{\mu \nu }\). Equation (4) is a singular perturbative equation.

An important observation that helps us to generalize the quintessence scalar field model in GUP, is that the fourth-order differential equation (4) follows from the variation of the action integral

$$\begin{aligned} S_{SF}=\int dx^{4}\sqrt{-g}L\left( \Psi ,\mathcal {D}_{\sigma }\Psi \right) , \end{aligned}$$
(5)

where \(L\left( \Psi ,\mathcal {D}_{\sigma }\Psi \right) \) is the usual Lagrangian for the Klein–Gordon equation

$$\begin{aligned} L\left( \Psi ,\mathcal {D}_{\sigma }\Psi \right) =\frac{1}{2}g^{\mu \nu }\mathcal {D}_{\mu }\Psi \mathcal {D}_{\nu }\Psi -\frac{1}{2}\left( \frac{mc}{\hbar }\right) ^{2}\Psi ^{2}, \end{aligned}$$
(6)

in which the new operator \(\mathcal {D}_{\mu }\) is defined as

$$\begin{aligned} \mathcal {D}_{\mu }=\nabla _{\mu }+\beta \hbar ^{2}\nabla _{\mu }\left( \Delta \right) , \end{aligned}$$
(7)

where \(\nabla _{\mu }\) is the covariant derivative for the metric tensor \(g_{\mu \nu }\).

We introduce the new field \(\Phi =\Delta \Psi \), and the Lagrange multiplier \(\lambda .\) Variation with respect to the Lagrange multiplier gives, \(\frac{\delta S}{\delta \lambda }=0\) from where it follows \(\lambda =-2\beta \hbar ^{2}\Phi \), such that the action integral (5) to be written as

$$\begin{aligned}&S=\int dx^{4}\sqrt{-g}\left( \frac{1}{2}g^{\mu \nu }\nabla _{\mu }\Psi \nabla _{\nu }\Psi +2\beta \hbar ^{2}g^{\mu \nu }\nabla _{\mu }\Psi \nabla _{\nu }\Phi \right. \nonumber \\&\qquad \left. +\beta \hbar ^{2}\Phi ^{2}-\frac{1}{2}V_{0}\Psi ^{2}\right) . \end{aligned}$$
(8)

Therefore, the new Lagrangian is

$$\begin{aligned}&L\left( \Psi ,\Psi _{;\mu },\Phi ,\Phi _{;\mu }\right) \nonumber \\&\quad =\sqrt{-g}\left( \frac{1}{2}g^{\mu \nu }\Psi _{;\mu }\Psi _{;\nu }+2\beta \hbar ^{2}g^{\mu \nu }\Psi _{;\mu }\Phi _{;\nu }\right) \nonumber \\&\qquad -\sqrt{-g}\left( \frac{1}{2}V_{0}\Psi ^{2}-\beta \hbar ^{2}\Phi ^{2}\right) \end{aligned}$$
(9)

Hence, the Euler–Lagrange equations for Lagrangian (9) are

$$\begin{aligned} g^{\mu \nu }\Psi _{,\mu \nu }-\Gamma ^{\mu }\Psi _{,\mu }-\Phi =0 \end{aligned}$$
(10)
$$\begin{aligned} 2\beta \hbar ^{2}\left( g^{\mu \nu }\Phi _{,\mu \nu }-\Gamma ^{\mu }\Phi _{,\mu }\right) +\left( V_{0}\Psi +\Phi \right) =0 \end{aligned}$$
(11)

where Eq. (10) is the constraint \(\Phi =\Delta _{g}\Psi \), and system (10), (11) is equivalent with the fourth-order partial differential equation (4).

3 Quintessence modified by the GUP

Quintessence is one of the simple dark energy models. Specifically, a canonical scalar field \(\phi \left( x^{\mu }\right) \) is responsible for the acceleration era of the universe. The late time acceleration or the early acceleration phase, where the scalar field plays the role of the inflaton [51].

In the case of a four-dimensional Riemannian manifold described by the metric tensor \(g_{\mu \nu }\), with Ricci scalar R, the Action integral which describes the gravitational theory for the quintessence model is

$$\begin{aligned} S_{Q}=\int dx^{4}\sqrt{-g}\left( -\frac{R}{2}+\frac{1}{2}g^{\mu \nu } \nabla _{\mu }\phi \nabla _{\nu }\phi -V(\phi )\right) . \end{aligned}$$
(12)

However, by using the operator \(\mathcal {D}_{\mu }\) defined in (7), we can generalize the latter Action integral as follows

$$\begin{aligned} S_{Q}^{GUP}=\int dx^{4}\sqrt{-g}\left( -\frac{R}{2}+\frac{1}{2}g^{\mu \nu }\mathcal {D}_{\mu }\phi \mathcal {D}_{\nu }\phi -V(\phi )\right) , \end{aligned}$$
(13)

which reduces to (12) when \(\beta =0\).

However, we can always define a new field \(\psi =\Delta \phi \), as before and write the Action Integral (13) like that of a second-order theory, that is,

$$\begin{aligned}&S_{Q}^{GUP}=\int dx^{4}\sqrt{-g}\left( -\frac{R}{2}+\frac{1}{2}g^{\mu \nu }\nabla _{\mu }\phi \nabla _{\nu }\phi -V(\phi )\right. \nonumber \\&\quad \left. +\beta \hbar ^{2}\left( 2g^{\mu \nu }\nabla _{\mu }\phi \nabla _{\nu }\psi +\psi ^{2}\right) \right) . \end{aligned}$$
(14)

Variation with respect to the metric tensor of (14) leads to the gravitational field equations, which are,

$$\begin{aligned} G_{\mu \nu }=T_{\mu \nu }, \end{aligned}$$
(15)

where \(G_{\mu \nu }\) is the Einstein tensor and \(T_{\mu \nu }\) is the energy momentum tensor for the two fields defined as

$$\begin{aligned} T_{\mu \nu }&=\left( \frac{1}{2}\nabla _{\mu }\phi \nabla _{\nu }\phi +2\beta \hbar ^{2}\nabla _{\mu }\phi \nabla _{\nu }\psi \right) \nonumber \\&\quad -g_{\mu \nu }\left( \frac{1}{2}g^{\alpha \beta }\nabla _{\alpha }\phi \nabla _{\beta }\phi +2\beta \hbar ^{2}g^{\alpha \beta }\nabla _{\alpha }\phi \nabla _{\beta }\psi \right. \nonumber \\&\quad \left. -V(\phi )+\beta \hbar ^{2}\psi ^{2}\right) , \end{aligned}$$
(16)

Furthermore, variation of (14) with respect to the scalar fields \(\phi ,~\psi \) leads to the equation of motions for the two fields

$$\begin{aligned}&\Delta \phi -\psi = 0, \end{aligned}$$
(17)
$$\begin{aligned}&\beta \hbar ^{2}\Delta \psi +\frac{1}{2}\left( \psi +V_{,\phi }\right) =0. \end{aligned}$$
(18)

4 FLRW background space

In the large scales our universe is almost homogeneous and isotropic and this description is well described by the Friedmann–Lemaître–Robertson–Walker (FLRW) line element with zero spatial curvature given by,

$$\begin{aligned} ds^{2}=dt^{2}-a^{2}\left( t\right) \left( dx^{2}+dy^{2}+dz^{2}\right) , \end{aligned}$$
(19)

where \(a\left( t\right) \) is the expansion scale factor of the universe, and \(H=\dot{a}/a\) is the Hubble function. For the line element (19) the gravitational field equations (15) are

$$\begin{aligned}&3H^{2}=\left( \frac{1}{2}\dot{\phi }^{2}+V\left( \phi \right) \right) +\beta \hbar ^{2}\left( 2\dot{\phi }\dot{\psi }-\psi ^{2}\right) +\rho _{m} \nonumber \\ \end{aligned}$$
(20)
$$\begin{aligned}&2\dot{H}+3H^{2}=-\left( \frac{1}{2}\dot{\phi }^{2}-V\left( \phi \right) +\beta \hbar ^{2}\left( 2\dot{\phi }\dot{\psi }+\psi ^{2}\right) \right) ,\nonumber \\ \end{aligned}$$
(21)

where \(\rho _{m}\) is the energy density for the dust fluid source, i.e. \(T_{\mu \nu }^{(m)}=\rho _{m}u_{\mu }u_{\nu }\), in which \(u^{\mu }=\delta _{t}^{\mu }\). For the fields \(\phi ,~\psi \) the equation of motions are

$$\begin{aligned} \ddot{\phi }+\frac{3}{a}\dot{a}\dot{\phi }-\psi&=0, \end{aligned}$$
(22)
$$\begin{aligned} \varepsilon \left( \ddot{\psi }+\frac{3}{a}\dot{a}\dot{\psi }\right) +\frac{1}{2}\left( \psi +V_{,\phi }\right)&=0, \end{aligned}$$
(23)

while for the dust fluid source, the Bianchi identity \(T_{~\ ~~~~;\nu }^{(m)~\mu \nu }=0\) gives \(\dot{\rho }_{m}+3H\rho _{m}=0,\) that is \(\rho _{m} =\rho _{m0}a^{-3}\), where \(\rho _{m0}\) is the present value of \(\rho _{m}\).

An equivalent way to write the field equations is to define the new variables

$$\begin{aligned}&3H^{2}=\rho _{GUP}+\rho _{m}, \end{aligned}$$
(24)
$$\begin{aligned}&2\dot{H}+3H^{2}=-p_{GUP}, \end{aligned}$$
(25)

where

$$\begin{aligned} \rho _{GUP}=\rho _{\phi }+\rho _{\psi },~p_{GUP}=p_{\phi }+p_{\psi }. \end{aligned}$$

The notations \(\rho _{\phi }\), \(~p_{\phi }\) are respectively the energy density and pressure of the quintessence field taking the explicit expressions

$$\begin{aligned} \rho _{\phi }=\frac{1}{2}\dot{\phi }^{2}+V\left( \phi \right) ,~~p_{\phi } =\frac{1}{2}\dot{\phi }^{2}-V\left( \phi \right) , \end{aligned}$$
(26)

and \(\rho _{\psi }\), \(p_{\psi }\) are respectively the energy and pressure for the second interacting field \(\psi \), that is,

$$\begin{aligned} \rho _{\psi }=\beta \hbar ^{2}( 2\dot{\phi }\dot{\psi }-\psi ^{2}) ~,~p_{\psi }=\beta \hbar ^{2}( 2\dot{\phi }\dot{\psi }+\psi ^{2}) . \end{aligned}$$
(27)

At this point we remark that the second field \(\psi \) is not a physical field, but it is introduced by the quantum corrections of GUP and it describes the new geometrodynamical degrees of freedom of GUP.

As far as concerns the parameter for the equation of state of the scalar field, defined as

$$\begin{aligned} w_{GUP}=\frac{p_{GUP}}{\rho _{GUP}}=\frac{\left( \frac{1}{2}\dot{\phi } ^{2}-V\left( \phi \right) \right) +\beta \hbar ^{2}\left( 2\dot{\phi } \dot{\psi }+\psi ^{2}\right) }{\left( \frac{1}{2}\dot{\phi }^{2}+V\left( \phi \right) \right) +\beta \hbar ^{2}\left( 2\dot{\phi }\dot{\psi }-\psi ^{2}\right) }. \end{aligned}$$
(28)

where expanding around \(\beta \hbar ^{2}\rightarrow 0\), it follows

$$\begin{aligned} w_{GUP}\simeq w_{\phi }+\frac{1}{\rho _{\phi }}\left( p_{\psi }-w_{\phi } \rho _{\psi }\right) +O( ( \beta \hbar ^{2}) ^{2}) \end{aligned}$$
(29)

while when the scalar field potential \(V\left( \phi \right) \), dominates, that is \(\dot{\phi }^{2}\ll V\), it follows

$$\begin{aligned} w_{GUP}\simeq -1+\frac{4\beta \hbar ^{2}}{V}\dot{\phi }\dot{\psi } \end{aligned}$$
(30)

which can cross the phantom divide line when \(~\dot{\phi }\dot{\psi }\le 0\). An important observation here is that the field equations form a singular perturbation system, also known as the slow-fast dynamical system. That means it is possible for \(~\left| \dot{\psi }\right| \) to be large enough such that the term \(\left| \beta \hbar ^{2}\dot{\psi }\right| \ \) will not be negligible. This is a main characteristic of the slow-fast dynamical systems.

In the following sections we study the global evolution of the dynamics, in order to understand the effects of GUP in the quintessence field. This work extends and generalizes the previous work on the specific theory [52].

5 Dynamical systems formulation

We continue our analysis by define the new dimensionless variables

$$\begin{aligned} x_{1}= & {} \frac{\dot{\phi }}{\sqrt{6}H}~,~y_{1}=\sqrt{\frac{V}{3H^{2}}}, x_{2}=\beta \hbar ^{2}\frac{2\sqrt{2}\dot{\psi }}{\sqrt{3}H} \nonumber \\ y_{2}= & {} \frac{\beta \hbar ^{2}\psi ^{2}}{3H^{2}},~\Omega _{m}=\frac{\rho _{m}}{3H^{2} }, \end{aligned}$$
(31)

such that the gravitational field equations to be written as

$$\begin{aligned} \frac{dx_{1}}{d\tau }= & {} \frac{1}{4}\left( 6x_{1}\left( x_{1}\left( x_{1} +x_{2}\right) -y_{1}^{2}-1\right) +y_{2}\left( 6x_{1}-\sqrt{6}\mu \right) \right) , \end{aligned}$$
(32)
$$\begin{aligned} \frac{dy_{1}}{d\tau }= & {} \frac{1}{2}y_{1}\left( 3\left( 1-y_{1}^{2} +y_{2}\right) +x_{1}\left( 3\left( x_{1}+x_{2}\right) -\sqrt{6} \lambda \right) \right) , \end{aligned}$$
(33)
$$\begin{aligned} \frac{dx_{2}}{d\tau }= & {} \frac{1}{2}\left( 3x_{1}x_{2}^{2}\!-\!3x_{2}\left( 1\!-\!x_{1}^{2}+y_{1}^{2}-y_{2}\right) \!+\!\sqrt{6}\left( 2\lambda y_{1}^{2}\!+\!\mu y_{2}\right) \right) \end{aligned}$$
(34)
$$\begin{aligned} \frac{dy_{2}}{d\tau }= & {} \frac{1}{4}y_{2}\left( 12x_{1}\left( x_{1} +x_{2}\right) +12\left( 1-y_{1}^{2}-y_{2}\right) -\sqrt{6}x_{2}\mu \right) \end{aligned}$$
(35)
$$\begin{aligned} \frac{d\mu }{d\tau }= & {} \frac{1}{4}\sqrt{\frac{3}{2}}x_{2}\mu ^{2}~,~\frac{d\lambda }{d\tau }=\sqrt{6}x_{1}\lambda ^{2}\left( \Gamma \left( \lambda \right) -1\right) \end{aligned}$$
(36)

where the new independent variable \(\tau =\ln a\), and variables \(\lambda ,\mu \) are defined as

$$\begin{aligned} \lambda =-\frac{V_{,\phi }}{V}~,~\beta \hbar ^{2}\mu =-\frac{2}{\psi } \end{aligned}$$
(37)

while \(\Gamma \left( \lambda \right) =\frac{V_{,\phi \phi }V}{V_{,\phi }^{2}}\). Moreover, the constraint equation is written as

$$\begin{aligned} \Omega _{m}\left( \mathbf {x,y}\right) =1-x_{1}\left( x_{1}+x_{2}\right) -y_{1}^{2}+y_{2}, \end{aligned}$$
(38)

from where it follows that the parameters are constraint as \(\,0\le \Omega _{m}\le 1.\)

The parameter for the equation of state for the effective fluid is derived to be

$$\begin{aligned} w_{tot}\left( \mathbf {x,y}\right) =x_{1}\left( x_{1}+x_{2}\right) -y_{1}^{2}+y_{2}, \end{aligned}$$
(39)

while we define the physical variables \(\Omega _{\phi },~\Omega _{\psi }\) such as, \(\Omega _{\phi }\left( \mathbf {x,y}\right) =x_{1}^{2}+y_{1}^{2}~,~\Omega _{\psi }\left( \mathbf {x,y}\right) =x_{1}x_{2}-y_{2},~\) where the constraint equation is written as \(\Omega _{m}\left( \mathbf {x,y}\right) =1-\Omega _{\phi }\left( \mathbf {x,y}\right) -\Omega _{\psi }\left( \mathbf {x,y}\right) \).

Finally, the parameter for the equation of state for the scalar field is expressed in terms of the new variables as

$$\begin{aligned} w_{\phi }\left( \mathbf {x,y}\right) =\frac{\left( x_{1}^{2}-y_{1} ^{2}\right) +x_{1}x_{2}+y_{2}}{\left( x_{1}^{2}+y_{1}^{2}\right) +x_{1}x_{2}-y_{2}}. \end{aligned}$$
(40)

5.1 Equilibrium points for exponential potential

Consider now that the scalar field potential is exponential, that is \(V\left( \phi \right) =V_{0}e^{-\lambda \phi }\), such that \(\Gamma \left( \lambda \right) =1\) and \(\frac{d\lambda }{d\tau }=0\). For that potential the dynamical system (32)–(36) is reduced by one-dimension. Therefore, the equilibrium points of the dynamical system are the points with coordinates \(P=\left( x_{1}\left( P\right) , x_{2}\left( P\right) , y_{1}\left( P\right) , y_{2}\left( P\right) , \mu \left( P\right) \right) \), where the right hand side of (32)–(36) are zero.

The equilibrium points are calculated to be

$$\begin{aligned} P_{0}= & {} \left( 0,0,0,0,0,\mu \right) ,~P_{1}^{\pm }=\left( \pm 1,0,0,0,\mu \right) , \end{aligned}$$
(41)
$$\begin{aligned} ~P_{2}= & {} \left( x_{1},\frac{1}{x_{1}}-x_{1},0,0,0\right) ,\nonumber \\ P_{3}= & {} \left( 0,\sqrt{\frac{2}{3}}\lambda y_{1}^{2},y_{1},y_{1}^{2}-1,0\right) . \end{aligned}$$
(42)

In the following lines we discuss the physical properties of the points as also their stability.

Point \(P_{0}\,\ \)describes a universe dominated by the dust fluid source, where \(\Omega _{m}\left( P_{0}\right) =1\), while the parameter for the equation of state \(w_{tot}\) is derived to be \(w_{tot}\left( P_{0}\right) =0\). The eigenvalues of the linearized system are calculated to be \(e_{1}\left( P_{0}\right) =3~,~e_{2}\left( P_{0}\right) =-\frac{3}{2}~,~e_{3}\left( P_{0}\right) =-\frac{3}{2},~e_{4}\left( P_{0}\right) =\frac{3}{2}\) and \(e_{5}\left( P_{0}\right) =0\), from where we infer that the point is a saddle point and the exact solution at the point is always unstable.

Points \(P_{1}^{\pm }\) describe universes dominated by the kinetic term of the scalar field, i.e. \(\Omega _{m}\left( P_{1}^{\pm }\right) =0\), and \(w_{tot}\left( P_{1}^{\pm }\right) =1.\) The eigenvalues are calculated to be \(e_{1}\left( P_{1}^{\pm }\right) =6,~e_{2}\left( P_{1}^{\pm }\right) =3,~e_{3}\left( P_{1}^{\pm }\right) =\frac{1}{2}\left( 6\mp \sqrt{6} \lambda \right) ,~e_{4}\left( P_{1}^{\pm }\right) =0\) and \(e_{5}\left( P_{1}^{\pm }\right) =0\), from where we infer that the exact solutions at the equilibrium points are always unstable. Point \(P_{1}^{+}\) is a saddle point when \(\lambda >\sqrt{6}\) while point \(P_{1}^{-}\) is a saddle point when \(\lambda <-\sqrt{6}\).

Point \(P_{2}\) is actually family of points that lie on a two-dimensional surface. The physical properties of this points are similar with that of points \(P_{1}^{\pm }\). The eigenvalues of the linearized system are derived \(e_{1}\left( P_{2}\right) =6,~e_{2}\left( P_{2}\right) =3,~e_{3}\left( P_{2}\right) =\frac{1}{2}\left( 6-\sqrt{6}x_{1}\lambda \right) ,~e_{4}\left( P_{2}\right) =0\) and \(e_{5}\left( P_{2}\right) =0\), from where we infer that the exact solutions at the family of points are always unstable.

The family of points on the surface with coordinates \(P_{3}\) describe de Sitter universes where \(\Omega _{m}\left( P_{3}\right) =0\) and \(w_{tot} \left( P_{3}\right) =-1\). The eigenvalues of the linearized system are determined to be \(e_{1}\left( P_{3}\right) =-3~,~e_{2}\left( P_{3}\right) =-3~,~e_{3}\left( P_{3}\right) =-3,~e_{4}\left( P_{3}\right) =0\) and \(e_{5}\left( P_{2}\right) =0\), hence the Center Manifold Theorem (CMT) should be applied to find the manifold where the de Sitter universes are attractors.

We observe that no scaling solutions or a tracking solutions exist in this specific model like in the quintessence theory. However, the critical points which describe the de Sitter solution do not exist in the case of quintessence for the exponential potential; these are the new equilibrium points provided by the new terms given by GUP.

We solve numerically the field equations (32)–(36) for three different sets of initial conditions, and we present the evolution of the physical parameters \(\Omega _{m},~\Omega _{\phi }=1-\Omega _{m}\) and \(w_{tot},~w_{\phi }\) if Fig. 1. For the three different sets of the initial conditions the final attractor is the de Sitter universe, while for the parameters of the equation of state we observe that they can cross the phantom divide line. In order to understand that behaviour we present the phase space diagram for the dynamical system in the plane \(\left\{ x_{1}-y_{1}\right\} \) from where it is clear that \(P_{3}\) can be a local attractor.

To analyze \(P_{3}\) in the following we proceed with the application of the CMT.

Fig. 1
figure 1

Qualitative evolution of the physical variables \(\Omega _{m},~\Omega _{\phi }\) (Left figs.) and \(w_{tot},~w_{\phi }\) (Right figs) for three sets of initial conditions. The plots of the first row is for initial conditions near to the point \(P_{0}\), the second row is for initial condition near to the point \(P_{1}^{+}\), while the plots of the third row are for initial conditions near to the point \(P_{1}^{-}\). For the energy density we observe that at late times the scalar field dominates \(\Omega _{\phi }\rightarrow 1\) and \(\Omega _{m}\rightarrow 0\), while the parameter for the equation of state \(w_{tot}\) have the limits \(w_{tot}\rightarrow -1\). Thus, we observe that \(w_{\phi }\) and \(w_{tot}\) can cross the phantom divide line and take values smaller than minus one

Full size image
Fig. 2
figure 2

Phase space diagram in the plane \(\left\{ x_{1}-y_{1}\right\} \) where for the rest of the parameters we selected \(\lambda =-1\), \(y_{2}=0.01\), \(x_{2}=0.01\) and \(\mu =0.1\). From the diagram we observe the attractor \(P_{3}\)

Full size image

5.1.1 Center manifold theorem for line of points \(P_{3}\)

Assuming \(y_{1c}\notin \left\{ 0, \frac{\sqrt{2}}{2}, 1\right\} , \lambda \ne 0\), and introducing the new variables

$$\begin{aligned} u_{1}&= \frac{1}{3} {y_{1c}} \left( \left( {y_{1c}}^{2}-1\right) \left( {y_{1c}} \left( \sqrt{6} \lambda {x_{1}}\right. \right. \right. \nonumber \\&\quad \left. \left. \left. +\lambda \mu \left( {y_{1c}}^{2}-1\right) +3\right) -6 {y_{1}}\right) +3 {y_{1c}} {y_{2}}\right) , \end{aligned}$$
(43)
$$\begin{aligned} u_{2}&= -\frac{1}{2} \lambda \mu {y_{1c}}^{2} \left( 2 {y_{1c}}^{4}-3 {y_{1c}}^{2}+1\right) , \end{aligned}$$
(44)
$$\begin{aligned} v_{1}&= \frac{1}{18} \Big \{6 \lambda {y_{1c}}^{2} \left( -2 \lambda {x_{1}} \left( {y_{1c}}^{2}-1\right) -\sqrt{6} \left( {y_{1c}}^{2}+{y_{2}}\right) \right) \nonumber \\&\quad +18 {x_{2}}+12 \sqrt{6} \lambda {y_{1}} \left( {y_{1c}}^{2}-1\right) {y_{1c}}\nonumber \\&\quad -\sqrt{6} \mu \left( {y_{1c}}^{2}-1\right) \left( \lambda ^{2} {y_{1c}}^{2} \left( 3 {y_{1c}}^{2}-2\right) +3\right) \Big \}, \end{aligned}$$
(45)
$$\begin{aligned} v_{2}&= \frac{1}{6} \lambda {y_{1c}}^{2} \left( -\lambda \left( 6 {x_{1}} \left( {y_{1c}}^{2}-2\right) +\sqrt{6} \mu \left( {y_{1c}}^{2}-1\right) ^{2}\right) \right. \nonumber \\&\quad \left. -3 \sqrt{6} \left( -2 {y_{1}} {y_{1c}}+{y_{1c}}^{2}+{y_{2} }+1\right) \right) , \end{aligned}$$
(46)
$$\begin{aligned} v_{3}&= -\frac{1}{3} \left( {y_{1c}}^{2}-1\right) \left( {y_{1c}}^{2} \left( \sqrt{6} \lambda {x_{1}}+\lambda \mu \left( {y_{1c}}^{2}-1\right) +3\right) \right. \nonumber \\&\quad \left. -6 {y_{1}} {y_{1c}}+3 ({y_{2}}+1)\right) , \end{aligned}$$
(47)

a particular point \(P_{3c}=\left( 0,\sqrt{\frac{2}{3}}\lambda y_{1c} ^{2},y_{1c},y_{1c}^{2}-1,0\right) \) on line of points \(P_{3}\) is translated to the origin and the linearization of the system (31), (32), (33), (34), (35) is transformed to its real Jordan canonical form:

$$\begin{aligned} \left( \begin{array} [c]{c} v_{1}^{\prime }\\ v_{2}^{\prime }\\ v_{3}^{\prime }\\ u_{1}^{\prime }\\ u_{2}^{\prime }\\ \end{array} \right) = \left( \begin{array} [c]{ccccc} -3 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} -3 &{} 1 &{} 0 &{} 0\\ 0 &{} 0 &{} -3 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 1\\ 0 &{} 0 &{} 0 &{} 0 &{} 0\\ &{} &{} &{} &{} \end{array} \right) \left( \begin{array} [c]{c} v_{1}\\ v_{2}\\ v_{3}\\ u_{1}\\ u_{2}\\ \end{array} \right) . \end{aligned}$$
(48)

Assuming \(y_{1c}\notin \left\{ 0, \frac{1}{\sqrt{2}}, 1\right\} , \lambda \ne 0\), and applying the CMT we obtain that the center manifold is given locally (up to second order) by the graph

$$\begin{aligned}&\Bigg \{(u_{1},u_{2}, v_{1}, v_{2}, v_{3})\in \mathbb {R}^{5}:\nonumber \\&\quad v_{1} = a_{1} u_{1}^{2} + a_{2} u_{1} u_{2} + a_{3} u_{2}^{2} + \mathcal {O}(3),\nonumber \\&\qquad v_{2} = b_{1} u_{1}^{2} + b_{2} u_{1} u_{2} + b_{3} u_{2}^{2} +\mathcal {O}(3),\nonumber \\&\qquad v_{3} = c_{1} u_{1}^{2} + c_{2} u_{1} u_{2} + c_{3} u_{2}^{2} +\mathcal {O}(3)\Bigg \}, \end{aligned}$$
(49)

where

$$\begin{aligned} a_{1}&= -\frac{\lambda \left( {y_{1c}}^{2}-1\right) }{2 \sqrt{6} {y_{1c}} ^{2}}, a_{2}\!=\! \frac{\lambda ^{2} \left( -11 {y_{1c}}^{4}+10 {y_{1c}} ^{2}-1\right) -6}{3 \sqrt{6} \lambda {y_{1c}}^{2} \left( 2 {y_{1c}}^{4}-3 {y_{1c}}^{2}+1\right) }, \nonumber \\ a_{3}&=\frac{6 \left( 4 {y_{1c}}^{2}-3\right) +\lambda ^{2} \left( 40 {y_{1c}}^{6}-57 {y_{1c}}^{4}+20 {y_{1c}}^{2}-1\right) }{18 \sqrt{6} \lambda \left( {y_{1c}}^{2}-1\right) \left( {y_{1c}}-2 {y_{1c} }^{3}\right) ^{2}}, \end{aligned}$$
(50)
$$\begin{aligned} b_{1}&=-\frac{1}{4} \sqrt{\frac{3}{2}} \lambda , b_{2}=\frac{5 \lambda }{\sqrt{6} \left( 2-4 {y_{1c}}^{2}\right) }, \nonumber \\ b_{3}&= \frac{\lambda ^{2} {y_{1c} }^{2} \left( 12 {y_{1c}}^{2}-5\right) +6}{12 \sqrt{6} \lambda \left( {y_{1c}}-2 {y_{1c}}^{3}\right) ^{2}}, \end{aligned}$$
(51)
$$\begin{aligned} c_{1}&= \frac{1}{4} \left( \frac{1}{{y_{1c}}^{2}}-1\right) , c_{2}=\frac{1}{2-4 {y_{1c}}^{2}}-\frac{1}{6 {y_{1c}}^{2}}, \nonumber \\ c_{3}&= \frac{{y_{1c}}^{2} \left( \lambda ^{2} \left( 12 {y_{1c}}^{4}-9 {y_{1c}}^{2}+1\right) +6\right) -6}{36 {y_{1c}}^{4} \left( \lambda -2 \lambda {y_{1c}}^{2}\right) ^{2}}. \end{aligned}$$
(52)

Hence, the dynamics on the center manifold is given locally (up to second order) by

$$\begin{aligned} u_{1}^{\prime }&=\frac{{u_{1}} {u_{2}} \left( 2 {y_{1c}}^{6}+7 {y_{1c}}^{4}-8 {y_{1c}}^{2}+1\right) }{4 {y_{1c}}^{6}-6 {y_{1c}}^{4}+2 {y_{1c}}^{2}}\nonumber \\&\quad +\frac{{u_{2}}^{2} ( \lambda ^{2} ( 2 {y_{1c}}^{4}-3 {y_{1c}} ^{2}+1) -6) }{6 \lambda ^{2} {y_{1c}}^{2} \left( 1-2 {y_{1c}} ^{2}\right) ^{2}}+{u_{2}} + \mathcal {O}(3), \end{aligned}$$
(53)
$$\begin{aligned} u_{2}^{\prime }&=-\frac{{u_{1}} {u_{2}}^{2}}{4 {y_{1c}}^{6}-6 {y_{1c}}^{4}+2 {y_{1c}}^{2}}\nonumber \\&\quad +\frac{{u_{2}}^{3} \left( \lambda ^{2} {y_{1c}}^{4}+3\right) }{6 {y_{1c}}^{4} \left( {y_{1c}}^{2}-1\right) \left( \lambda -2 \lambda {y_{1c} }^{2}\right) ^{2}}\nonumber \\&\quad -\frac{{u_{2}}^{2}}{4 {y_{1c}}^{4}-6 {y_{1c}}^{2}+2} + \mathcal {O}(3). \end{aligned}$$
(54)
Fig. 3
figure 3

Dynamics on the center manifold is given locally (up to second order) by (53), (54) for several choices of the parameters

Full size image

Defining the function

$$\begin{aligned} U(u_{1}, u_{2})=\frac{{u_{1}}^{2}}{2}+{u_{2}}^{2}, \end{aligned}$$
(55)

we have

$$\begin{aligned} U^{\prime }(\tau )&=\nabla U \cdot (u_{1}^{\prime }, u_{2}^{\prime })=u_{1} u_{2}\nonumber \\&\quad + \frac{{u_{1}}^{2} {u_{2}} \left( 2 {y_{1c}}^{6}+7 {y_{1c}}^{4}-8 {y_{1c} }^{2}+1\right) }{4 {y_{1c}}^{6}-6 {y_{1c}}^{4}+2 {y_{1c}}^{2}}\nonumber \\&\quad +{u_{1}} \left( \frac{{u_{2}}^{2} \left( \lambda ^{2} \left( 2 {y_{1c}}^{4}-3 {y_{1c}}^{2}+1\right) -6\right) }{6 \lambda ^{2} {y_{1c}}^{2} \left( 1-2 {y_{1c}}^{2}\right) ^{2}}\right. \nonumber \\&\quad \left. -\frac{{u_{2}}^{3}}{2 {y_{1c}}^{6}-3 {y_{1c}} ^{4}+{y_{1c}}^{2}}\right) \nonumber \\&\quad +\frac{{u_{2}}^{4} \left( 2 \lambda ^{2} {y_{1c}}^{4}+6\right) }{6 {y_{1c} }^{4} \left( {y_{1c}}^{2}-1\right) \left( \lambda -2 \lambda {y_{1c}} ^{2}\right) ^{2}}\nonumber \\&\quad +\frac{{u_{2}}^{3}}{-2 {y_{1c}}^{4}+3 {y_{1c}}^{2}-1}\sim u_{1} u_{2} + \mathcal {O}(3). \end{aligned}$$
(56)

Hence, for small \(u_{1}, u_{2}\) such as \(u_{1} u_{2} < 0\) the origin is locally stable, whereas, for small \(u_{1}, u_{2}\) such as \(u_{1} u_{2} >0\) it is unstable.

In the Fig. 3 is it presented the dynamics on the center manifold is given locally (up to second order) by (53), (54) for several choices of the parameters. In Fig. 4 it is presented \(U^{\prime }(\tau )\) for the same choices of the parameters, showing that for small \(u_{1}, u_{2}\) such as \(u_{1} u_{2} < 0\) the origin is locally stable, whereas, for small \(u_{1}, u_{2}\) such as \(u_{1} u_{2} >0\) it is unstable.

Special case \(y_{1c}=\frac{\sqrt{2}}{2}\). Using the transformation

$$\begin{aligned}&{u_{1}}= \mu , \end{aligned}$$
(57)
$$\begin{aligned}&{u_{2}}= \frac{1}{24} \left( \lambda \mu -2 \sqrt{6} \lambda {x_{1}}+12 \sqrt{2} {y_{1}}+12 {y_{2}}-6\right) , \end{aligned}$$
(58)
$$\begin{aligned}&{v_{1}}= \frac{1}{144} \left( -\sqrt{6} ( \lambda ^{2}-12) \mu -12 \lambda ( -2 \lambda {x_{1}}+4 \sqrt{3} {y_{1}}\right. \nonumber \\&\qquad \left. +2 \sqrt{6} {y_{2} }+\sqrt{6}) \right) +{x_{2}}, \end{aligned}$$
(59)
$$\begin{aligned}&{v_{2}}= -\frac{1}{48} \lambda \left( \sqrt{6} (\lambda \mu +18)\!-\!36 \lambda {x_{1}}\!-\!24 \sqrt{3} {y_{1}}\!+\!12 \sqrt{6} {y_{2}}\right) , \end{aligned}$$
(60)
$$\begin{aligned}&{v_{3}}= \frac{1}{24} \left( -\lambda \mu +2 \sqrt{6} \lambda {x_{1}}-12 \sqrt{2} {y_{1}}+12 {y_{2}}+18\right) , \end{aligned}$$
(61)

the linearization transforms to

$$\begin{aligned} \left( \begin{array} [c]{c} v_{1}^{\prime }\\ v_{2}^{\prime }\\ v_{3}^{\prime }\\ u_{1}^{\prime }\\ u_{2}^{\prime }\\ \end{array} \right) = \left( \begin{array} [c]{ccccc} -3 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} -3 &{} 1 &{} 0 &{} 0\\ 0 &{} 0 &{} -3 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0\\ &{} &{} &{} &{} \end{array} \right) \left( \begin{array} [c]{c} v_{1}\\ v_{2}\\ v_{3}\\ u_{1}\\ u_{2}\\ \end{array} \right) . \end{aligned}$$
(62)

Applying the CMT we obtain that the center manifold is given locally (up to second order) by the graph

$$\begin{aligned}&\Bigg \{(u_{1},u_{2}, v_{1}, v_{2}, v_{3})\in \mathbb {R}^{5}:\nonumber \\&\quad v_{1} = a_{1} u_{1}^{2} + a_{2} u_{1} u_{2} + a_{3} u_{2}^{2} + \mathcal {O} (3),\nonumber \\&\quad v_{2} = b_{1} u_{1}^{2} + b_{2} u_{1} u_{2} + b_{3} u_{2}^{2} + \mathcal {O}(3),\nonumber \\&\quad v_{3} = c_{1} u_{1}^{2} + c_{2} u_{1} u_{2} + c_{3} u_{2}^{2} + \mathcal {O}(3)\Bigg \}, \end{aligned}$$
(63)

where

$$\begin{aligned} a_{1}&= 0, a_{2}= \frac{1}{\sqrt{6}}, a_{3}= 0, b_{1}=0, b_{2}= 0, b_{3}= 0, \nonumber \\ c_{1}&=-\frac{1}{192}, c_{2}= 0, c_{3}= \frac{1}{4}. \end{aligned}$$
(64)

Hence, the dynamics on the center manifold is given locally (up to second order) by

$$\begin{aligned}&u_{1}^{\prime }\!=\! \frac{1}{192} {u_{1}}^{2} ( {u_{1}} ( \!-\!\lambda ^{2}+24 {u_{2}}-12) \!+\!24 (\lambda +2 \lambda {u_{2}})) \!+\! \mathcal {O}(3), \end{aligned}$$
(65)
$$\begin{aligned}&u_{2}^{\prime }=-\frac{1}{32} {u_{1}} ({u_{1}}-8 \lambda {u_{2}}) + \mathcal {O}(3) . \end{aligned}$$
(66)

According to Fig. 5, the origin is a saddle.

Fig. 4
figure 4

\(U^{\prime }(\tau )\) for several choices of the parameters, showing that for small \(u_{1}, u_{2}\) such as \(u_{1} u_{2} < 0\) the origin is locally stable, whereas, for small \(u_{1}, u_{2}\) such as \(u_{1} u_{2} >0\) it is unstable

Full size image
Fig. 5
figure 5

Dynamics on the center manifold is given locally (up to second order) by (65), (66) for \(\lambda =-1\), (resp. \(\lambda =1\)). The dynamics is qualitative the same for other values of \(\lambda <0\) (resp. \(\lambda >0\))

Full size image

Special case \(y_{1c}=1\). Using the transformation

$$\begin{aligned}&u_{1} = \mu ,u_{2}= y_{2}, \nonumber \\&v_{1}= \frac{1}{2} \left( -2 \sqrt{6} \lambda +2 \lambda ^{2} x_{1}+2 \sqrt{6} \lambda y_{1}-\sqrt{6} \lambda y_{2}\right) , \nonumber \\&v_{2} = \frac{1}{2} (2 y_{1}-y_{2}-2),\nonumber \\&v_{3}= \frac{1}{3} \left( -\sqrt{6} \lambda +3 x_{2}-\sqrt{6} \lambda y_{2}\right) , \end{aligned}$$
(67)

the linearization transforms to

$$\begin{aligned} \left( \begin{array} [c]{c} v_{1}^{\prime }\\ v_{2}^{\prime }\\ v_{3}^{\prime }\\ u_{1}^{\prime }\\ u_{2}^{\prime }\\ \end{array} \right) = \left( \begin{array} [c]{ccccc} -3 &{} 1 &{} 0 &{} 0 &{} 0\\ 0 &{} -3 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} -3 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0\\ &{} &{} &{} &{} \end{array} \right) \left( \begin{array} [c]{c} v_{1}\\ v_{2}\\ v_{3}\\ u_{1}\\ u_{2}\\ \end{array} \right) . \end{aligned}$$
(68)

Applying the CMT we obtain that the center manifold is given locally (up to second order) by the graph

$$\begin{aligned}&\Bigg \{(u_{1},u_{2}, v_{1}, v_{2}, v_{3})\in \mathbb {R}^{5}: \nonumber \\&\quad v_{1} = a_{1} u_{1}^{2} + a_{2} u_{1} u_{2} + a_{3} u_{2}^{2} + \mathcal {O} (3),\nonumber \\&\quad v_{2} = b_{1} u_{1}^{2} + b_{2} u_{1} u_{2} + b_{3} u_{2}^{2} + \mathcal {O}(3),\nonumber \\&\quad v_{3} = c_{1} u_{1}^{2} + c_{2} u_{1} u_{2} + c_{3} u_{2}^{2} + \mathcal {O}(3)\Bigg \}, \end{aligned}$$
(69)

where

$$\begin{aligned}&a_{1}=0, a_{2}=0, a_{3}=-\frac{1}{4} \sqrt{\frac{3}{2}} \lambda , b_{1}=0, b_{2}= \frac{\lambda }{12},\nonumber \\&\quad b_{3}= -\frac{1}{8}, c_{1}=0, c_{2}=\frac{\lambda ^{2}+3}{3 \sqrt{6}}, c_{3}=0. \end{aligned}$$
(70)

Hence, the dynamics on the center manifold is given locally (up to second order) by

$$\begin{aligned}&u_{1}^{\prime }=\frac{1}{24} u_{1}^{2} (( \lambda ^{2}+3) u_{1} u_{2}+6 \lambda (u_{2}+1)) + \mathcal {O}(3), \end{aligned}$$
(71)
$$\begin{aligned}&u_{2}^{\prime }=-\frac{1}{2} \lambda u_{1} u_{2} (3 u_{2}+1) + \mathcal {O}(3) . \end{aligned}$$
(72)

According to Fig. 6, the origin is a saddle.

Fig. 6
figure 6

Dynamics on the center manifold is given locally (up to second order) by (71), (72) for \(\lambda =-1\), (resp. \(\lambda =1\)). The dynamics is qualitative the same for other values of \(\lambda <0\) (resp. \(\lambda >0\))

Full size image

5.2 Equilibrium points for arbitrary potential

Consider now an arbitrary potential \(V\left( \phi \right) \), such that \(\Gamma \left( \lambda \right) \) to be arbitrary. In that consideration, the dynamical system of our consideration has dimension six, where the equilibrium points have coordinates \(\bar{P}=\left( x_{1}\left( \bar{P}\right) ,x_{2}\left( \bar{P}\right) , y_{1}\left( \bar{P}\right) ,\right. \left. y_{2}\left( \bar{P}\right) ,\mu \left( \bar{P}\right) ,\lambda \left( \bar{P}\right) \right) \).

In particular the equilibrium points are \(\bar{P}_{A}=\left( x_{1}\left( \bar{P}_{A}\right) ,\right. \left. x_{2}\left( \bar{P}_{A}\right) ,y_{1}\left( \bar{P}_{A}\right) ,y_{2}\left( \bar{P}_{A}\right) ,\mu \left( \bar{P}_{A}\right) ,\lambda \left( \bar{P}_{A}\right) \right) \), where \(A=1,2,3\) and \(\lambda \left( P_{A}\right) \) solves the algebraic equation \(\lambda ^{2}\left( \Gamma \left( \lambda \right) -1\right) =0\). While there are more equilibrium points because of the different values of \(\lambda \left( P_{A}\right) \), the physical physical properties of the equilibrium points are exactly the same with the exponential potential. That is an interesting result, since any kind of potential has a similar physical behaviour.

Since the physical properties on the additional equilibrium points do not change, we omit the presentation of the stability analysis.

Fig. 7
figure 7

Qualitative evolution of the cosmographic parameters q,  j and s, for the same initial conditions as used for Fig. 2. Right Figs are subplots of the Left Figs in the range near the de Sitter attractor

Full size image

6 Evolution of the cosmographic parameters

There are usually two known ways through which one tries to understand the dynamical history of our universe. One is the model independent approach in which one deals with the kinematical quantities that are extracted directly from the space-time metric and secondly, one needs to fix either the matter sector of the universe or the underlying gravitational sector and then explore the dynamics of the universe. Indeed, the kinematical quantities are more flexible and novel compared to the model dependent approach since within the kinematical approach the dynamics of the universe is solely dependent on the geometrical quantities and hence this is more robust. The study of the kinematical quantities derived from a homogeneous and isotropic universe is called the cosmography [53]. Using the variation of the kinematical quantities, one can measure the viability of a proposed cosmological scenario.

Under the background of the homogeneous and isotropic universe characterized by the FLRW universe, one can define the kinematical quantities as [54, 55]

$$\begin{aligned}&H(t) = \frac{1}{a}\frac{da}{dt},\\&q(t) = -\frac{1}{a}\frac{d^{2}a}{dt^{2}}\left[ \frac{1}{a}\frac{da}{dt}\right] ^{-2},\\&j(t) = \frac{1}{a}\frac{d^{3}a}{dt^{3}}\left[ \frac{1}{a}\frac{da}{dt}\right] ^{-3},\\&s(t) = \frac{1}{a}\frac{d^{4}a}{dt^{4}}\left[ \frac{1}{a}\frac{da}{dt}\right] ^{-4}, \end{aligned}$$

where H, q, j, s are respectively called the Hubble, deceleration, jerk, snap parameters and except for the Hubble parameter, the last three parameters are dimensionless. With these terminology, the expansion scale factor a(t) can be written as

$$\begin{aligned} a(t)&= a_{0} \left[ 1 + H_{0} (t-t_{0}) - \frac{1}{2} q_{0}^{2} (t-t_{0})^{2} + \frac{1}{3!} j_{0} (t-t_{0})^{3}\right. \nonumber \\&\quad \left. + \frac{1}{4!} s_{0} (t-t_{0})^{4} + \mathcal {O} [(t-t_{0})^{5}] \right] , \end{aligned}$$
(73)

in which any sub-index attached to any quantity refers to its present value. We note that one can extend this series by allowing higher order derivatives of the scale factor beyond order four. Here we are interested to investigate the deceleration parameter, and its next two hierarchies, namely, jerk and the snap parameter which in terms of the Hubble parameter can alternatively be rewritten as

$$\begin{aligned} q&=-1-\frac{\dot{H}}{H^{2}} \end{aligned}$$
(74)
$$\begin{aligned} j&=\frac{\ddot{H}}{H^{3}}-3q-2 \end{aligned}$$
(75)
$$\begin{aligned} s&=\frac{\dddot{H}}{H^{4}}+4j+3q\left( q+4\right) +6, \end{aligned}$$
(76)

where an overhead dot represents the cosmic time derivative; therefore, two or three overhead dots represents the two- or three-times derivative with respect to the cosmic time. Thus, from the evolution of q, j, s, one can understand the expansion of the universe (accelerating or decelerating), the rate of acceleration (j) and its next derivative. In Fig. 7 we present the qualitative evolution of the cosmographic parameters.

7 Conclusions

The dynamics of the universe both at its very early and late evolution are two mysterious issues in cosmology. Various approaches have been proposed to understand both the phases of the universe but none of them are found to be extremely satisfactory. The classical theory of gravity does not work out during the very early evolution of the universe but the quantum gravity does a bit, according to the available quantum gravitational theories. The quantum gravity theory predicts the existence of a minimum length scale of the order of the Planck length \(l_{\mathrm {pl}}\) and consequently, this motivated to generalize the Heisenberg’s uncertainty principle to some generalized uncertainty principle (GUP) in quantum gravity. The generalized uncertainty principle has been found to be very effective to explain several important issues that are directly related to our universe and its evolution. In particular, in the context of black hole evaporation, GUP has a significant role. Moreover, GUP has a explanation towards the existence of the dark matter fluid in the universe’s sector. Thus, it will be interesting whether the GUP modified cosmological models could give rise to some interesting features of the universe. Hence, in the present work we study the modified quintessence model and investigate the stability of the model by performing its dynamical analysis.

We first investigate the modified quintessence model for the exponential potential which gives rise to various possibilities. For various critical points, we get different solutions such as the matter dominated universe and the universe dominated by the kinetic term only, however, all of them are found to be unstable in nature. Only the family of points on the surface with coordinates \(P_{3}\) describe the de Sitter universes and they all are attractors. It is interesting to note that although no scaling solutions or a tracking solution exist in this specific model contrary to the quintessence theory, however, the critical points describing the de Sitter solution do not exist in the case of a quintessence model with an the exponential potential. This is a new result in this work. In order to understand the behavior of the family of critical points, \(P_{3}\), we applied the Center Manifold Theorem.

Furthermore, we considered an arbitrary potential and performed the dynamical system analysis similar to the exponential potential. We found that the physical properties of the equilibrium points for the general potential remain exactly same as in the case with the exponential potential. That is a very striking and exciting result in this work, because this clearly indicates that the investigations with the exponential potential is enough for the modified quintessence models with GUP.

Finally, we observe that the parameter for the equation of state for the effective cosmological fluid can cross the phantom divine line and it was observed by the numerical simulations without the appearance of ghosts, while the late time attractor is the de Sitter universe. Of course, for the initial conditions of our numerical simulations we have selected those where point \(P_{3}\) is an attractor as it is presented by the Center Manifold Theorem. However, it should be noted here and we also stress that, keeping aside the initial conditions for the numerical simulations, the ghost issue should be investigated further for a definite answer towards this direction.