Interacting quintessence in light of generalized uncertainty principle: cosmological perturbations and dynamics

Abstract

We consider a cosmological scenario endowed with an interaction between the universe’s dark components – dark matter and dark energy. Specifically, we assume the dark matter component to be a pressure-less fluid, while the dark energy component is a quintessence scalar field with Lagrangian function modified by the quadratic Generalized Uncertainty Principle. The latter modification introduces new higher-order terms of fourth-derivative due to quantum corrections in the scalar field’s equation of motion. Then, we investigate asymptotic dynamics and general behaviour of solutions of the field equations for some interacting models of special interests in the literature. At the background level, the present interacting model exhibits the matter-dominated and de Sitter solutions which are absent in the corresponding quintessence model. Furthermore, to boost the background analysis, we study cosmological linear perturbations in the Newtonian gauge where we show how perturbations are modified by quantum corrected terms from the quadratic Generalized Uncertainty Principle. Depending on the coupling parameters, scalar perturbations show a wide range of behavior.

1 Introduction

The physics of the dark sector of our Universe is obscure, at least according to the up-to-date observational shreds of evidence. This dark sector comprises of a dark matter (DM) which is almost pressure-less or cold, and a dark energy (DE) fluid, and they jointly contribute nearly 96% ($\sim$ 68% for DE and $\sim$ 28% for DM) of the entire energy content of the Universe [1]. The most mathematically simplest cosmological model that explains this dark sector is the $\Lambda$-Cold Dark Matter ($\Lambda$CDM) in which $\Lambda >0$, the cosmological constant, plays the role of DE. Even though the $\Lambda$CDM model has been found to fit excellently to most of the available astrophysical and cosmological probes, the most significant assumption within such a cosmological picture is the independent evolution of both DM and DE, that means, in other words, we have a non-interacting cosmological scenario. However, despite the tremendous success of the $\Lambda$CDM model, there are already known severe issues associated with it. We can reconcile the existing problems related to the $\Lambda$CDM model by an alternative description of the Universe along with the excellent fits to the available observational data. The list of the models is quite extensive, including a variety of cosmological theories and models; see for instance [2]. Amongst many attractive cosmological models, in this article, we focus on a very generalized cosmological theory where DM and DE are freely allowed to interact with one other, known as Interacting DE (IDE) or coupled DE.

The IDE theory which was formally proposed by Amendola [3] received tremendous attention in the scientific community for its ability to offer a possible explanation to the cosmic coincidence problem [4,5,6,7,8,9,10] where a direct interaction between DM and DE (in a non-gravitational way) can play the game. However, we recall an earlier work by Wetterich [11] where the author first argued the possibility of an interaction in the Universe where a coupling or interaction between a scalar field and gravity can lead to a time-dependent cosmological ‘constant’ which asymptotically becomes constant and consequently a possible explanation to the cosmological constant problem was placed. With such motivations, IDE theory was investigated widely in the literature. Subsequently it was proved that interaction between DM and DE could offer some more interesting possibilities, for example, the crossing of the phantom divide line without any need to scalar field [4, 12,13,14,15], and most importantly, it can explain the recent $H_0$ tension [16,17,18] and the $\sigma _8$ tension [18,19,20]. We refer to a recent review on the $H_o$ tension where the observational data report that IDE rocks (see Tables B1 and B2 of [21]). Nevertheless, the questions regarding the derivation of the interaction function from some fundamental action have not been truly understood yet. Attempts to derive the interaction function have been made by various investigators with some reasonable answers, see for instance [22,23,24,25,26,27,28,29], however, the development of this section is still in progress.

In this article, we consider an interacting scheme where DM is pressure-less, and DE is a scalar field. Even though we are considering an interacting scalar field model which reminds us the earlier works in this direction [3, 4, 30], however, the framework and the formulation of the present interacting scalar field scenario is different from them [3, 4, 30] because here the scalar field sector which is interacting with the pressure-less DM has been modified following the generalized uncertainty principle (GUP). The inclusion of GUP adds a novel feature to this work. According to the past historical records, GUP has many astrophysical and cosmological implications, for instance, the origin of the magnetic fields in the Universe sector [31], black hole thermodynamics [32] and some others (see [33, 34] for a detailed description). Since GUP plays a vital role in the context of the early physics of the Universe, therefore, its effects on the late Universe physics is a topical issue that should be investigated, aiming to understand the complete dynamical picture of the Universe. It is always very appealing to encounter the new physical theories modified by the GUP’s inclusion to understand the minimum length scale effects in the late Universe physics. Following this motivation, the authors of [34] first investigated the dynamics of the GUP modified quintessence scalar field with two different potentials; one is exponential, and the other is an arbitrary. The resulted dynamics offered some interesting possibilities that are absent in the usual quintessence scalar field model. For instance, the appearance of the critical points describing the de Sitter Universe in the GUP modified quintessence is a new piece of result, and this does not depend on the choice of the potential. In this article, we extended the previous work [34] by allowing an interaction between the GUP modified quintessence scalar field and the pressure-less DM to examine the dynamical features of this generalized cosmic scenario.

To improve the background analysis, we also investigate cosmological linear perturbations. More specifically, we study the dynamics of scalar perturbations during the matter-dominated era, where the effect of GUP might imprint exciting features during the growth of structure. Such analysis is crucial to identify specific signatures of the interacting scalar field theories in the light of GUP to look for by cosmological observations. Therefore, the present study will serve as a preliminary investigation of the modified interacting IDE model’s cosmic viability.

The article is structured as follows: In Sect. 2 we present the basic equations of the interacting scenario. For the dark energy, we assume that it is described by a minimally coupled scalar field where the GUP modifies the Action Integral. The new set of equations are presented in Sect. 3. We continue our analysis to study the asymptotic solutions for the field equations for a linear interacting model. The detailed analysis of the asymptotic solutions and their physical properties are given in Sect. 4. Furthermore, in Sect. 5, we derive the scalar equations of the cosmological linear perturbations. We find that the perturbed system is also a singular perturbation system due to the new terms introduced by GUP. Therefore, the qualitative evolution of the perturbed equations in the slow-fast manifolds is studied in the matter-dominated era from where we find new growing modes that follow from GUP corrections. Finally, in Sect. 6 we summarize the main results of this article.

2 Interacting dark energy

According to cosmological principle, our Universe in the large scales is almost homogeneous and isotropic. Therefore, we consider that the background physical space is described by Friedmann–Lemaître–Robertson–Walker (FLRW) line-element with zero curvature given by

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

(1)

where $a\left( t\right)$ is the expansion scale factor and (t, x, y, z) are the co-moving coordinates. In the context of Einstein’s General Relativity (GR) the gravitational field equations are (in the units where $8 \pi G =1$):

$$\begin{aligned} G_{~\nu }^{\mu }=T_{~\nu }^{\mu },{\mu }, \nu = 0,1,2,3 \end{aligned}$$

(2)

where $G_{~\nu }^{\mu }=R_{~\nu }^{\mu }-\frac{1}{2}Rg_{~\nu }^{\mu }$ is the Einstein tensor and $T_{~\nu }^{\mu }$ is the effective momentum tensor

$$\begin{aligned} T_{\mu \nu }=T_{\mu \nu }^{\left( m\right) }+T_{\mu \nu }^{\left( d\right) }+T_{\mu \nu }^{\left( r\right) }+T_{\mu \nu }^{\left( b\right) }, \end{aligned}$$

(3)

in which $T_{\mu \nu }^{\left( b\right) }=\rho _{b}u_{\mu }u_{\nu }$ describes baryonic matter, $T_{\mu \nu }^{\left( r\right) }=\left( \rho _{r} +p_{r}\right) u_{\mu }u_{\nu }$ corresponds to radiation, $T_{\mu \nu }^{\left( m\right) }=\rho _{m}u_{\mu }u_{\nu }$ describes dark matter component and $T_{\mu \nu }^{\left( d\right) }=\left( \rho _{d}+p_{d}\right) u_{\mu }u_{\nu }+p_{d}g_{\mu \nu }$ is the energy-momentum tensor for DE fluid source which has a negative equation of state parameter $w_{d}= p_{d}/\rho _{d}<0$ such that the late-time accelerating phase of the Universe is realized [35]. Here $u^{\mu }$ describes a comoving observer $u^{\mu }=\delta _{t}^{\mu },~u^{\mu }u_{\mu }=-1$ and $\rho _i$ ($i= b, r, m, d$) stands for the energy density of i-th fluid; $p_r$ and $p_d$ denote the pressure of radiation and DE fluid. We note that baryons and DM are assumed to be pressureless.

For spatially flat FLRW background space the gravitational field equations (2) read

$$\begin{aligned}&3H^{2}=\rho _{m}+\rho _{d}+\rho _{r}+\rho _{b}, \end{aligned}$$

(4)

$$\begin{aligned}&-2{\dot{H}}-3H^{2}=p_{d}+p_{r}, \end{aligned}$$

(5)

where an overhead dot represents the derivative with respect to the cosmic time and $H={\dot{a}}/a$ is the Hubble function. From Bianchi identity $\nabla _{\nu }G^{\mu \nu }=0$, we get $\nabla _{\nu }T^{\mu \nu }=0,$ which means,

$$\begin{aligned} \nabla _{\nu }T^{\left( m\right) \mu \nu }+\nabla _{\nu }T^{\left( d\right) \mu \nu }+\nabla _{\nu }T^{\left( r\right) \mu \nu }+\nabla _{\nu }T^{\left( b\right) \mu \nu }=0. \end{aligned}$$

(6)

Often, it is assumed that the fluid components do not interact with each other, that means, $\nabla _{\nu }T^{\left( m\right) \mu \nu }=0$ , $\nabla _{\nu }T^{\left( d\right) \mu \nu }=0$ , $\nabla _{\nu }T^{\left( r\right) \mu \nu }=0$ and $\nabla _{\nu }T^{\left( b\right) \mu \nu }=0$. However, over the last several years, cosmological models allowing an interaction between DM and DE have drawn a significant attention. In this case the conservation equation (6) gives the following conditions $~\nabla _{\nu }T^{\left( r\right) \mu \nu }=0$ ,$~\nabla _{\nu }T^{\left( b\right) \mu \nu }=0$ and $\nabla _{\nu }T^{\left( m\right) \mu \nu }+\nabla _{\nu }T^{\left( d\right) \mu \nu }=0$ where the last relation between DM and DE can be decoupled into two separated equations by introducing an interaction term Q as follows

$$\begin{aligned} \nabla _{\nu }T^{\left( m\right) \mu \nu }=Q,\quad \nabla _{\nu }T^{\left( d\right) \mu \nu }=-Q. \end{aligned}$$

(7)

We note that the nature of interaction function Q is unknown and there is not a unique theoretical model which describes its origin [36,37,38,39]. In terms of phenomenology, there are various approaches in the literature which have shown that a nonzero function Q is supported by the cosmological observations [40,41,42,43,44,45,46,47].

For the background space (1), conservation equations (7) take the forms

$$\begin{aligned}&{\dot{\rho }}_{m}+3H\rho _{m} = Q, \end{aligned}$$

(8)

$$\begin{aligned}&{\dot{\rho }}_{d}+3H\left( \rho _{d}+p_{d}\right) =-Q. \end{aligned}$$

(9)

Thus, once an interaction function Q is given, one can solve the conservation equations either analytically or numerically (depending on the nature of the interaction function) and finally using the Friedmann equation (4), the dynamics of interacting Universe can be explored. Let us make an important comment that is essential to understand the dynamics in the following sections. In the total energy density of the Universe (see Eq. (4)) the contributions of baryonic matter and radiation fluid are very small compared to that of DM and DE, therefore, the total dynamics of the Universe is not influenced by their joint contribution. Hence, we shall omit their contribution from now on and we shall focus on the dynamics of the Universe driven mainly by DM and DE. That means, the effective cosmological fluid responsible for the Universe’s dynamics is identified with the energy-momentum tensor $T_{\mu \nu }^\mathrm{eff}=T_{\mu \nu }^{\left( m\right) }+T_{\mu \nu }^{\left( d\right) }$.

As mentioned earlier, the DE fluid is quintessence scalar field, one of the simplest time-varying DE models [48]. In quintessence DE, energy density and pressure terms of$~T_{\mu \nu }^{\left( d\right) }$ are written as

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

(10)

where $V (\phi )$ is the potential of quintessence scalar field. With the above expressions for energy density and pressure of the quintessence DE, the set of Eqs. (8), (9) can now be explicitly written as

$$\begin{aligned}&{\dot{\rho }}_{m}+3H\rho _{m} =Q, \end{aligned}$$

(11)

$$\begin{aligned}&{\dot{\phi }}\left( \ddot{\phi }+3H{\dot{\phi }}+V_{,\phi }\right) =-Q. \end{aligned}$$

(12)

The latter system has been the case of study of various analysis in the literature [49,50,51,52,53,54,55,56,57,58,59,60] where $Q=Q\left( \rho _{m},\rho _{d}\right)$ has been assumed to be either a linear function $Q=\alpha H\rho _{m}+\kappa ~H\rho _{d}$ or a nonlinear function such as $Q=\alpha \frac{H\rho _{m}}{{\dot{\phi }}^{2}\rho _{\phi }},~Q=\alpha \frac{H{\dot{\phi }}^{2}\rho _{\phi }}{\rho _{m}}$, $Q=-\frac{1}{2}(4-3 \gamma ) \rho _m {\dot{\phi }}\frac{\chi '(\phi )}{\chi (\phi )}$ [46,47,48,49, 61,62,63,64,65,66] and many others. Additionally, with the above choices for the coupling function, various functional forms of the scalar field potential function $V\left( \phi \right)$ can also be chosen [33, 35, 50,51,52, 60, 67,68,69,70,71,72].

In this work we are interested to extend the analysis of IDE models for a generalization of the quintessence scalar field model inspired by GUP.

3 Quintessence modified by the GUP

In this section we offer a brief description about GUP and the modifications of quintessence scalar field due to this principle. The quadratic GUP is based on the introduction of a minimum measurable length in which the Heisenberg uncertainty principle is 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}$$

(13)

where $\beta =\beta _{0}\ell _{Pl}^{2}/2\hbar ^{2}$ and $\beta _{0}~$is the deformation parameter [73,74,75,76]. The deformation parameter can be positive or negative [77, 78]. Moreover, there are various proposed theoretical constraints in the literature from quantum systems [79, 80] and also from gravitational systems [81, 82].

The main characteristic of GUP is that it modifies the second-order Klein–Gordon equation for a scalar field into a fourth-order. Inspired by this observation in Ref. [33] a generalized quintessence scalar field model modified by the GUP has been proposed; in which Klein–Gordon equation is modified.

The Action Integral for the scalar field has been proposed to be [33]

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

(14)

in which$~{\mathcal {D}}_{\mu }=\nabla _{\mu }+\beta \hbar ^{2}\nabla _{\mu }\left( \Delta \right) ~$and $\Delta$ is the Laplace operator. For zero deformation parameter Action Integral (14) reduces to that of quintessence. For $\beta _{0}\ne 0$, with the use of a Lagrange multiplier the Action Integral (14) is written in the equivalent form

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

(15)

where the new scalar field $\psi$ attributes the degrees of freedom provided by higher-order derivatives.

Therefore, in this context DE components read

$$\begin{aligned} \rho _{d}=\rho _{\phi }+\rho _{\psi },\quad p_{d}=p_{\phi }+p_{\psi }, \end{aligned}$$

(16)

where $\rho _{\phi },~p_{\phi }$ are usual terms of energy density and pressure for quintessence scalar field given by (10). The quantities $\rho _{\psi }$, $p_{\psi }$ are respectively the energy density and pressure terms which correspond to higher-order derivatives and are described by kinematic and dynamic quantities for the second field $\psi$ as follows

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

(17)

Hence, using (16) and (17) the interacting conservation equations (11), (12) are modified as

$$\begin{aligned}&{\dot{\rho }}_{m}+3H\rho _{m} =Q, \end{aligned}$$

(18)

$$\begin{aligned}&2\beta \hbar ^{2}{\dot{\phi }}\left( \ddot{\psi }+3H{\dot{\psi }}\right) +\dot{\phi }\left( \psi +V_{,\phi }\right) =-Q, \end{aligned}$$

(19)

with constraint equation

$$\begin{aligned} \ddot{\phi }+3H{\dot{\phi }}-\psi =0, \end{aligned}$$

(20)

which follows from definition of Lagrange multiplier, where the dot means derivative with respect to cosmic time t.

We continue our analyses by studying the dynamics and asymptotic behaviour of the solutions of the gravitational field equations for various interacting functions of the form $Q=Q\left( \rho _{m},\rho _{d}\right)$ in the light of GUP. In particular, we classify stationary points according to their stability and determine generic features irrespective of the initial conditions. This is achieved using dynamical systems tools. Concerning the choice of interaction function we consider its linear form:

$$\begin{aligned} Q_{A} =\frac{{\dot{\phi }}}{\sqrt{6}}\left( \alpha \rho _{m}+\kappa ~\rho _{d}\right) , \end{aligned}$$

(21)

where $\alpha$ and $\kappa$ are coupling parameters of the interaction function. Notice that for different values of $\alpha$ and $\kappa$, one can realize a number of interaction functions. For instance, $\alpha =0$ reduces to $Q_{A} =\frac{{\dot{\phi }}}{\sqrt{6}}\; \kappa \rho _{d}$. Similarly, one can also get another two versions, namely, $Q_{A} =\frac{{\dot{\phi }}}{\sqrt{6}}\; \alpha \rho _{m}$ for $\kappa =0$ and $Q_{A} =\frac{{\dot{\phi }}}{\sqrt{6}} \delta \left( \rho _{m}+\rho _{d}\right)$ for $\alpha = \kappa = \delta$. These kind of interactions were motivated from scalar-tensor theories [11, 83]. In particular, model $Q_{A} =\frac{{\dot{\phi }}}{\sqrt{6}}\; \alpha \rho _{m}$ is conformally equivalent to power-law potential model of Brans–Dicke theory. As far as the scalar field potential is concerned we consider the exponential potential $V\left( \phi \right) =V_{0}e^{-\lambda \phi }$, which as it was found in [34] can describe different kinds of asymptotic solutions with physical interest.

4 Asymptotic solutions and stability

To study asymptotic evolution of interacting gravitational model we define dimensionless variables following [34]:

$$\begin{aligned}&x_{1}=\frac{{\dot{\phi }}}{\sqrt{6}H},\quad y_{1}=\sqrt{\frac{V}{3H^{2}}} ,\quad 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}},\quad \Omega _{m}=\frac{\rho _{m}}{3H^{2}}. \end{aligned}$$

(22)

Using the new variables, the field equations (4), (5), (11), (12) reduce to an algebraic-differential system of first-order. We determine the stationary points of the new system and investigate their stability. Every stationary point describes an asymptotic solution for the cosmological field equations.

4.1 Interaction $Q_{A}$

For interaction model $Q_{A}$, the field equations reduce to the equivalent system

$$\begin{aligned} \frac{dx_{1}}{d\tau }&= \frac{1}{4}\left( 6x_{1}\left( x_1(x_1+x_2)-y_1^{2}-1\right) +y_2\left( 6x_1-\sqrt{6}\mu \right) \right) , \end{aligned}$$

(23)

$$\begin{aligned} \frac{dy_{1}}{d\tau }&= \frac{1}{2}\,y_1\left( x_1\,\left( -\sqrt{6}\lambda +3\,x_1+3x_2\right) -3y_1^{2}+3y_2+3\right) , \end{aligned}$$

(24)

$$\begin{aligned} \frac{d x_2}{d\tau }&= \frac{1}{2}\left( -2\alpha \left( -x_{1}\left( x_1+x_2\,\right) -y_{1}^{2}+y_{2}+1\right) \right. \nonumber \\&\quad +x_1(3x_2-2\kappa )\left( x_1+x_2\right) +y_{1}^{2}\left( 2\sqrt{6}\lambda -2\kappa -3 x_2 \right) \nonumber \\&\quad \left. +\,y_{2}\left( \sqrt{6}\mu +2\kappa + 3x_2\right) -3 x_2 \right) , \end{aligned}$$

(25)

$$\begin{aligned} \frac{dy_{2}}{d\tau }&=\frac{1}{4}y_2 \left( 12 x_1^{2}+12\, x_{1}x_{2}-\sqrt{6}\mu \,x_2 +12\left( -y_1^{2}+y_2+1\right) \right) , \end{aligned}$$

(26)

$$\begin{aligned} \frac{d\mu }{d\tau }&=\frac{1}{4}\sqrt{\frac{3}{2}}\mu ^{2}\,x_2, \end{aligned}$$

(27)

with the constraint equation

$$\begin{aligned} \Omega _{m}=1-( x_{1}^{2}+y_{1}^{2}) -( x_{1}x_{2} -y_{2}) , \end{aligned}$$

(28)

where the new phase space variable $\mu$ is defined as $\mu =-2\left( \beta \hbar ^{2} \psi \right) ^{-1}$. From the requirement that the dimensionless energy densities of matter fluids must be non-negative, i.e., $\Omega _m\ge 0$ and $\Omega _d:= \left( x_{1}^{2}+y_{1}^{2}\right) +\left( x_{1}x_{2}-y_{2}\right) \ge 0$ we have the physical part of the phase space is $0 \le \left( x_{1}^{2}+y_{1}^{2}\right) +\left( x_{1}x_{2} -y_{2}\right) \le 1$. However, it is important to mention here, that from expression (28) the phase space variables are not bounded, but the set of initial conditions in which parameters $x_2, y_2$ reach infinity are not physically acceptable (because $x_2, y_2$ are proportional to $\beta$ and $\beta \ll 1$ from the model’s construction). For this reason one can omit the analysis when $x_2$ or $y_2$ tends to infinity.

In the new variables the effective equation of state parameter $w_\mathrm{eff} =-1-\frac{2}{3}\frac{{\dot{H}}}{H^{2}}$ is written as

$$\begin{aligned} w_\mathrm{eff}\left( x_{1},x_{2},y_{1},y_{2}\right) =x_{1}^{2}-y_{1}^{2}+x_{1} x_{2}+y_{2}. \end{aligned}$$

(29)

Note that at a stationary point the scale factor is expressed as $a\left( t\right) =a_{0}t^{\frac{2}{3\left( 1+w_\mathrm{eff}\right) }},$ for $w_\mathrm{eff} \ne -1$ and $a\left( t\right) =a_{0}e^{H_{0}t}$ for $w_\mathrm{eff}=-1$. On the other hand, according to (17), we have

$$\begin{aligned} w_{\psi }= \frac{p_{\psi }}{\rho _{\psi }}= \frac{\left( 2{\dot{\phi }}{\dot{\psi }}-\psi ^{2}\right) }{\left( 2{\dot{\phi }}{\dot{\psi }}+\psi ^{2}\right) }= \frac{x_1 x_2-y_2}{x_1 x_2+y_2}. \end{aligned}$$

(30)

The deceleration parameter $q =-1-\frac{{\dot{H}}}{H^{2}}$ is written as

$$\begin{aligned} q= \frac{1}{2} (3 w_\mathrm{eff}+1)= \frac{1}{2} \left( 3 x_1 (x_1+x_2)-3 y_1^2+3 y_2+1\right) . \end{aligned}$$

(31)

It is useful to define $c_{d}^2$, the effective sound speed of dark energy perturbations (the corresponding quantity for matter is zero in the dust case) as

$$\begin{aligned}&c_{d}^2 = w_{d} - \frac{w_{d}'}{3(1+w_{d})}, \end{aligned}$$

(32)

where $w_d=\frac{p_d}{\rho _d}$ and prime denotes derivative with respect to $\tau$. In terms of phase-space variables we have

$$\begin{aligned} c_{d}^2&= \frac{x_1 x_2^2 \left( 12 x_1-\sqrt{6} \mu y_2\right) +x_2 y_2 \left( 12 x_1+\sqrt{6} \mu y_2\right) -2 \sqrt{6} x_1 y_2 \left( 2 \lambda y_1^2+\mu y_2\right) }{12 x_1 x_2 (x_1 x_2+y_2)}\nonumber \\&\quad -\frac{\alpha y_2 \left( x_1 (x_1+x_2)+y_1^2-y_2-1\right) }{3 x_2 (x_1 x_2+y_2)}+\frac{\kappa y_2 \left( x_1 (x_1+x_2)+y_1^2-y_2\right) }{3 x_2 (x_1 x_2+y_2)}. \end{aligned}$$

(33)

Before studying the general case for arbitrary values of the coupling parameters $\left\{ \alpha ,\kappa \right\}$ in the following sections we consider the special cases where $\kappa =0$ or $\alpha =0$.

Fig. 1

Qualitative evolution of $w_\mathrm{eff}\left( a\right)$ (left figure) and $\Omega _{m}\left( a\right)$ (right figure) for the interacting model $Q_A$ with $\kappa =0$. Solid line $\left( \alpha =5\right)$ and dashed line $\left( \alpha =-5\right)$ are for the initial conditions $\left( x_{1}^{0},x_{2} ^{0},y_{1}^{0},y_{2}^{0},\mu ^{0}\right) =\left( 0,-\frac{2}{3} \alpha ,0,-0.01.0.2\right)$, while the dotted line is for $\alpha =-5$ with initial conditions $\left( x_{1}^{0},x_{2}^{0},y_{1}^{0},y_{2}^{0},\mu ^{0}\right) =\left( 0.5,0,0,0.01,0.2\right)$. The unique attractor of the dynamical system is point $A_{2}$ which describes the de Sitter universe

4.1.1 Interaction $Q_{A}$ with $\kappa =0$

For $\kappa =0$, the dynamical system consisted by the differential equations (23)–(27) admits the stationary points ${\mathbf {A}} ={\mathbf {A}}\left( \mathbf {x,y,}\mu ;\lambda ,\alpha \right)$, $\mathbf {x=} \left( x_{1},x_{2}\right) ,~{\mathbf {y}}=\left( y_{1},y_{2}\right)$ with coordinates

$$\begin{aligned} A_{0}^{\pm }&=\left( \pm 1,0,0,0,\mu \right) ,\quad A_{1} =\left( x_{1},\frac{1}{x_{1}}-x_{1},0,0,0\right) , \\ A_{2}&=\left( 0,\sqrt{\frac{2}{3}}\lambda y_{1}^{2},y_{1},y_{1}^{2}-1,0\right) ,\\ A_{3}&=\left( 0,-\frac{2}{3}\alpha ,0,0,0\right) ,\quad A_{4}=\left( \frac{3}{\alpha },0,0,-1-\frac{9}{\alpha ^{2}},\frac{6\sqrt{6}\alpha }{9+\alpha ^{2}}\right) . \end{aligned}$$

While the family of points $A_0^{\pm }$ constitutes a line, $A_{1},~A_{2}$ describe surfaces in the phase space. We continue with the presentation of physical properties of the stationary points and investigate asymptotic solutions’ stability properties. To determine the stability behavior, we examine the nature of eigenvalues corresponding to a perturbed matrix of the system (23)–(27) for each point.

The family of points $A_{0}^{\pm }$ describe asymptotic solutions where only the kinetic parts of scalar field contribute to the effective cosmological fluid, that is, $\Omega _{m}\left( A_{0}^\pm \right) =0$ and $w_\mathrm{eff}\left( A_{0}^\pm \right) =1$. Hence, the asymptotic solution for the scale factor is $a\left( t\right) =a_{0}t^{\frac{1}{3}}$. The eigenvalues of the linearized system around each stationary point are $e_{1}\left( A_{0}^\pm \right) =6$, $e_{2}\left( A_{0}^\pm \right) =3\pm \alpha$, $e_{3}\left( A_{0}^\pm \right) =\frac{1}{2}\left( 6\mp \sqrt{6}\lambda \right)$, $e_{4}\left( A_{0}^\pm \right) =0$ and $e_{5}\left( A_{0}^\pm \right) =0$. Therefore, the family of points $A_{0}^{+}$ is unstable when $\alpha \,>-3$ and $\lambda \,<\sqrt{6}$, otherwise it is saddle. Also, the family of points $A_{0}^{-}$ is unstable when $\alpha \,<3$ and $\lambda \,>-\sqrt{6}$, otherwise it is saddle.

The family of points $A_{1}$ describe the same physical properties as that of $A_0^\pm$. The family of points $A_1$ and $A_{0}^\pm$ intersect at points $(\pm 1,0,0,0,0)$. The eigenvalues of the linearized system around the stationary point are $e_{1}\left( A_{1}\right) =6$, $e_{2}\left( A_{1}\right) =3+\alpha x_{1}$, $e_{3}\left( A_{1}\right) =\frac{1}{2}\left( 6-\sqrt{6}\lambda x_{1}\right)$, $e_{4}\left( A_{1}\right) =0$ and $e_{5}\left( A_{1}\right) =0$. Looking from the nature of eigenvalues, the family of points $A_{1}$ is unstable when $\alpha \,x_1>-3$ and $\lambda \,x_1<\sqrt{6}$, otherwise it is saddle as one of the eigenvalues is always positive. Therefore, the asymptotic solutions at $A_{1}$ are not stable.

The family of points $A_{2}$ describes de Sitter asymptotic solution where the dynamical part of scalar field dominates in the field equations, that is, $\Omega _{m}\left( A_{2}\right) =0,~w_\mathrm{eff}\left( A_{2}\right) =-1$ and $a\left( t\right) =a_{0}e^{H_{0} t}$. The eigenvalues of the linearized system are $e_{1}\left( A_{2}\right) =-3$, $e_{2}\left( A_{2}\right) =-3$, $e_{3}\left( A_{2}\right) =-3$, $e_{4}\left( A_{2}\right) =0$, $e_{5}\left( A_{2}\right) =0$. Therefore, we apply the center manifold theorem (CMT) to analyze the stability of $A_2$.

Point $A_{3}$ describes a universe dominated by DM fluid $\Omega _{m}\left( A_{3}\right) =1~w_\mathrm{eff}\left( A_{3}\right) =0$ with scale factor $a\left( t\right) =a_{0}t^{\frac{2}{3}}$. The eigenvalues of the linearized system are $e_{1}\left( A_{3}\right) =3$, $e_{2}\left( A_{3}\right) =-\frac{3}{2}$, $e_{3}\left( A_{3}\right) =-\frac{3}{2}$, $e_{4}\left( A_{3}\right) =\frac{3}{2}$, $e_{5}\left( A_{3}\right) =0$. Therefore, the asymptotic solution is not stable and the point is a saddle.

Finally, for point $A_{4}$ we find $\Omega _{m}\left( A_{4}\right) =-\frac{18}{\alpha ^{2}}<0$ and so the point is not physically acceptable. Therefore it is not necessary to investigate its stability.

Compared with the corresponding model of usual quintessence theory [84], the present model does not contain any accelerated scaling solution. However, the present model exhibits new stationary points describing late time de Sitter solution and matter-dominated solution.

In Fig. 1 we present the qualitative evolution of the physical parameters $w_\mathrm{eff}\left( a\right)$ and $\Omega _{m}\left( a\right)$ for various sets of initial conditions near to the matter dominated era. Note that we have chosen initial conditions where the final attractor is the de Sitter solution described by the point $A_{2}$.

4.1.2 Interaction $Q_{A}$ with $\alpha =0$

For arbitrary parameter $\kappa$ and $\alpha =0$, the stationary points of the dynamical system (23)–(27) are

$$\begin{aligned} {\bar{A}}_{1}&=\left( 0,0,0,0,\mu \right) ,\\ {\bar{A}}_{2}&=\left( 0,\frac{\sqrt{6}\lambda y_{1}^{2}-\kappa }{3},y_{1},y_{1}^{2}-1,0\right) ,\\ {\bar{A}}_{3}&=\left( \frac{\sqrt{6}}{\lambda },\frac{\left( \lambda ^{2}-12\right) \kappa -\sqrt{6}\lambda \left( \lambda ^{2}-6\right) }{2\lambda \left( \sqrt{6}\kappa -3\lambda \right) }, \right. \\&\quad \left. \sqrt{\frac{\kappa \left( 9\lambda -2\sqrt{6}\kappa \right) }{2\sqrt{6}\kappa ^{2}-30\kappa \lambda +9\sqrt{6}\lambda ^{2}}},0,0\right) ,\\ {\bar{A}}_{4}^{\left( \pm \right) }&=\left( -\frac{3\pm \sqrt{9-2\kappa ^{2}}}{2\kappa },0,0,-\frac{9+\kappa ^{2}\pm 3\sqrt{9-2\kappa ^{2}}}{2\kappa ^{2} }, \right. \\&\quad \left. -\frac{2\sqrt{6}\kappa \left( 9\pm \sqrt{9-2\kappa ^{2}}\right) }{36+\kappa ^{2}}\right) . \end{aligned}$$

The family of stationary points ${\bar{A}}_{1}$ describes a universe dominated by dark matter fluid, that is, $\Omega _{m}\left( {\bar{A}}_{1}\right) =1,~w_\mathrm{eff} \left( {\bar{A}}_{1}\right) =0$ and $a\left( t\right) =a_{0}t^{\frac{2}{3}}$. The eigenvalues of the linearized system are $e_{1}\left( {\bar{A}} _{1}\right) =3$, $e_{2}\left( {\bar{A}}_{1}\right) =-\frac{3}{2}$, $e_{3}\left( {\bar{A}}_{1}\right) =-\frac{3}{2}$, $e_{4}\left( {\bar{A}} _{1}\right) =\frac{3}{2}$, $e_{5}\left( {\bar{A}}_{1}\right) =0$ which means that the point is a saddle point and the asymptotic solution is not stable.

${\bar{A}}_{2}$ describes a family of stationary points of de Sitter asymptotic solutions, with physical properties $\Omega _{m}\left( {\bar{A}}_{2}\right) =0,~w_\mathrm{eff}\left( {\bar{A}}_{2}\right) =-1$ and $a\left( t\right) =a_{0}e^{H_{0}t}$. The eigenvalues are exactly that of $A_{2}$, therefore we shall apply CMT.

For point ${\bar{A}}_{3}$, we find $\Omega _{m}\left( {\bar{A}}_{3}\right) =\frac{2\kappa }{2\kappa -\sqrt{6}\lambda },~w_\mathrm{eff}\left( {\bar{A}}_{3}\right) =1$. Therefore, for the point to be real and physically acceptable $\left\{ y_{1}\left( {\bar{A}}_{3}\right) \ge 0,~0\le \Omega _{m}\left( {\bar{A}}_{3}\right) \le 1\right\}$, we must have $\kappa =0$ i.e. there is no interaction. Hence, we omit the study of its stability.

Finally, for the point ${\bar{A}}_{4}^{\left( \pm \right) }$, we obtain $\Omega _{m}\left( {\bar{A}}_{4}^{\left( \pm \right) }\right) =1-\frac{9}{\kappa ^{2} }\mp \frac{3}{\kappa ^{2}}\sqrt{9-2\kappa ^{2}}$, $w_\mathrm{eff}\left( {\bar{A}} _{4}^{\left( \pm \right) }\right) =-1$. However, the point is not physically acceptable as there is no value of $\kappa$ satisfying the physical conditions $\left\{ 0\le \Omega _{m}\left( {\bar{A}}_{4}^{\left( \pm \right) }\right) \le 1,9-2\kappa ^{2}\ge 0\,\right\}$.

4.1.3 Interaction $Q_{A}$ with arbitrary $\alpha$ and $\kappa$

In the general case where the coupling parameters $\alpha$ and $\kappa$ are arbitrary we find the stationary points

$$\begin{aligned} A_{1}^{\prime }&=\left( 0,-\frac{2\alpha }{3},0,0,0\right) ,\\ A_{2} ^{\prime }&=\left( 0,\frac{\sqrt{6}\lambda y_{1}^{2}-\kappa }{3},y_{1},y_{1} ^{2}-1,0\right) ,\\ A_{3}^{\prime }&=\left( \frac{\sqrt{6}}{\lambda },\frac{\left( \lambda ^{2}-6\right) \left( 2\alpha +\sqrt{6}\lambda \right) -\left( \lambda ^{2}-12\right) \kappa }{2\lambda \left( \sqrt{6}\left( \alpha -\kappa \right) +3\lambda \right) }, \right. \\&\quad \left. \sqrt{\frac{\kappa \left( 2\sqrt{6}\left( \alpha -\kappa \right) +9\lambda \right) }{4\sqrt{6}\left( \alpha -\kappa \right) ^{2}+3\lambda \left( 10\left( \alpha -\kappa \right) +3\sqrt{6}\lambda \right) }},0,0\right) ,\\ A_{4}^{\prime \left( \pm \right) }&=\left( \frac{3\pm K_{4}}{2\left( \alpha -\kappa \right) },0,0,-\frac{9+2\alpha ^{2}-3\alpha \kappa +\kappa ^{2} \pm 3K_{4}}{2\left( \alpha -\kappa \right) ^{2}}, \right. \\&\quad \left. -\frac{2\sqrt{6} \left( \left( \kappa -2\,\alpha \right) \left( 3\pm K_4 \right) +6\,\kappa \right) }{ \left( 2\,\alpha -\kappa \right) ^{2}+36}\right) . \end{aligned}$$

where $K_{4}=\sqrt{9+2\kappa \left( \alpha -\kappa \right) }$.

The stationary point $A_{1}^{\prime }$ describes a matter dominated solution with $\Omega _{m}\left( A_{1}^{\prime }\right) =1,~w_\mathrm{eff}\left( A_{1}^{\prime }\right) =0$ and asymptotic solution $a\left( t\right) =a_{0}t^{\frac{2}{3}}$. The linearized system around the stationary point admits eigenvalues $e_{1}\left( A_{1}^{\prime }\right) =3$, $e_{2}\left( A_{1}^{\prime }\right) =-\frac{3}{2}$, $e_{3}\left( A_{1}^{\prime }\right) =-\frac{3}{2}$, $e_{4}\left( A_{1}^{\prime }\right) =\frac{3}{2}$, $e_{5}\left( A_{1}^{\prime }\right) =0$ from where we infer that $A_{1}^{\prime }$ is a saddle point.

The family of stationary points described by $A_{2}^{\prime }$ have the same physical properties and eigenvalues as that of family of points $A_{2}$ and ${\bar{A}}_{2}$. Hence, we shall apply the CMT to investigate the stability of the points.

Point $A_{3}^{^{\prime }}$ describes an asymptotic solution with $\Omega _{m}\left( A_{3}^{\prime }\right) =-\frac{2\kappa }{\sqrt{6} \lambda +2(\alpha -\kappa )}$ and $w_\mathrm{eff}\left( A_{3}^{\prime }\right) =1$. We can easily find that the asymptotic solution at the stationary point is physically acceptable only when $\kappa =0$, which belongs to a family of points $A_{1}$ for $x_1=\frac{\sqrt{6}}{\lambda }$.

For a stationary point $A_{4}^{\prime \left( \pm \right) }$, we find that $\Omega _{m}\left( A_{4}^{\prime \left( \pm \right) }\right) =\frac{\kappa ^{2}-a\kappa -3\left( 3\pm K_{4}\right) }{\left( \alpha -\kappa \right) ^{2}}$, $w_\mathrm{eff}\left( A_{4}^{\prime \left( \pm \right) }\right) =-1$. However, this point is not physically acceptable as there are no real values for $\kappa$ and $\alpha$ such that $0\le \Omega _m\le 1$. In the next subsection, we shall proceed with the analysis of a non-hyperbolic equilibrium point $A_{2}^{^{\prime }}$.

4.1.4 CMT for $A_{2}^{\prime }$

We use the Center Manifold Theorem to analyze the stability of the stationary point $A_{2}^{^{\prime }}$ of system (23)–(27) and we discuss the specific cases $\kappa =0,~$or $\alpha =0$ which correspond to the points $A_{2}$ and ${\bar{A}}_{3}$ respectively.

Firstly, we assume $(\sqrt{6} \kappa + 6 (1 - 2 y_c^2) \lambda )\ne 0, y_c\ne 0, y_c^2\ne 1$. Then, evaluating the Jacobian matrix of system (23)–(27) at $A_{2}^{^{\prime }}$ we obtain

$$\begin{aligned} J(A_{2}^{^{\prime }})= \scriptscriptstyle \left( \begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} -3 &{} 0 &{} 0 &{} 0 &{} -\frac{1}{2} \sqrt{\frac{3}{2}} \left( y_1^2-1\right) \\ \frac{1}{2} y_1 \left( \sqrt{6} \lambda \left( y_1^2-1\right) -\kappa \right) &{} -3 y_1^2 &{} 0 &{} \frac{3 y_1}{2} &{} 0 \\ -\frac{1}{6} \left( \kappa -\sqrt{6} \lambda y_1^2\right) \left( 2 \alpha -3 \kappa +\sqrt{6} \lambda y_1^2\right) &{} -y_1 \left( -2 \alpha +\kappa +\sqrt{6} \lambda \left( y_1^2-2\right) \right) &{} -3 &{} \frac{1}{2} \left( -2 \alpha +\kappa +\sqrt{6} \lambda y_1^2\right) &{} \sqrt{\frac{3}{2}} \left( y_1^2-1\right) \\ \left( y_1^2-1\right) \left( \sqrt{6} \lambda y_1^2-\kappa \right) &{} -6 y_1 \left( y_1^2-1\right) &{} 0 &{} 3 \left( y_1^2-1\right) &{} \frac{1}{12} \left( y_1^2-1\right) \left( \sqrt{6} \kappa -6 \lambda y_1^2\right) \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ \end{array} \right) \,,\nonumber \\ \end{aligned}$$

(34)

which can be reduced to its canonical Jordan Form

$$\begin{aligned} j(A_{2}^{^{\prime }})= \left( \begin{array}{ccccc} -3 &{} 1 &{} 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) , \end{aligned}$$

(35)

defining a similarity matrix s which satisfies $J(A_{2}^{^{\prime }})= s j(A_{2}^{^{\prime }}) s^{-1}$. The eigenvalues of $J(A_{2}^{^{\prime }})$ and $j(A_{2}^{^{\prime }})$ are $-3,-3,-3,0,0$. Defining new variables

$$\begin{aligned}&( v_1, v_2, v_3, u_1, u_2)\nonumber \\&\quad = s^{-1}\cdot \left( x_1,y_1-y_c,x_2+\frac{1}{3} \left( \kappa -\sqrt{6} \lambda y_c^2\right) ,y_2-y_c^2+1,\mu \right) ,\nonumber \\ \end{aligned}$$

(36)

which translates a fixed point on the curve $A_{2}^{^{\prime }}$ labelled by the value of $y_c$ to the origin, that is,

$$\begin{aligned} v_1&=\frac{1}{216} \Bigg (12 \Big (3 \sqrt{6} \mu +y_c^2 \Big (-\left( 3 \sqrt{6} \mu +\lambda \left( \lambda \left( y_c^2-1\right) \right. \right. \nonumber \\&\quad \left. \left. \times \left( 12 x_1+\sqrt{6} \mu \left( 3 y_c^2-2\right) \right) +6 \sqrt{6} \left( y_2+y_c^2\right) \right) \right) \Big ) \nonumber \\&\quad +12 \sqrt{6} \lambda y_1 \left( y_c^2-1\right) y_c\Big ) +\kappa \left( 72-2 \mu \left( y_c^2-1\right) \right. \nonumber \\&\quad \left. \times \left( \sqrt{6} \alpha +3 \lambda y_c^2\right) \right) +3 \sqrt{6} \kappa ^2 \mu \left( y_c^2-1\right) \Bigg )+x_2, \end{aligned}$$

(37a)

$$\begin{aligned} v_2&= \frac{x_1 \left( \kappa -\sqrt{6} \lambda \left( y_c^2-1\right) \right) \left( 2 \alpha -\kappa +\sqrt{6} \lambda y_c^2\right) ^2}{6 \left( \kappa (\kappa -2 \alpha )+\sqrt{6} \lambda \left( \kappa +2 \alpha \left( y_c^2-1\right) -2 \kappa y_c^2\right) +6 \lambda ^2 y_c^2 \left( y_c^2-1\right) \right) ^2} \nonumber \\&\quad \times \Bigg (\kappa ^2 (3 \kappa -2 \alpha )+6 \lambda ^2 y_c^2 \left( \kappa \left( 3 y_c^2-2\right) -2 \alpha \left( y_c^2-1\right) \right) \nonumber \\&\quad +\sqrt{6} \kappa \lambda \left( \alpha \left( 4 y_c^2-2\right) +\kappa \left( 3-7 y_c^2\right) \right) \nonumber \\&\quad +6 \sqrt{6} \lambda ^3 y_c^2 \left( y_c^4-3 y_c^2+2\right) \Bigg ) \nonumber \\&\quad +y_c (y_1-y_c) \left( 2 \alpha -\kappa +\sqrt{6} \lambda y_c^2\right) \nonumber \\&\quad +\frac{1}{2} \left( y_2-y_c^2+1\right) \left( -2 \alpha +\kappa -\sqrt{6} \lambda y_c^2\right) \nonumber \\&\quad -\frac{\mu \left( y_c^2-1\right) \left( 2 \alpha -\kappa +\sqrt{6} \lambda y_c^2\right) }{72 \left( \kappa (\kappa -2 \alpha )+\sqrt{6} \lambda \left( \kappa +2 \alpha \left( y_c^2-1\right) -2 \kappa y_c^2\right) +6 \lambda ^2 y_c^2 \left( y_c^2-1\right) \right) ^2} \nonumber \\&\quad \times \Bigg (3 \sqrt{6} \kappa ^5-24 \kappa ^2 \lambda \left( y_c^2-1\right) \nonumber \\&\quad \times \left( 2 \alpha ^2+\sqrt{6} \alpha \lambda \left( 5 y_c^2-2\right) +3 \lambda ^2 y_c^2 \left( 3 y_c^2-1\right) \right) \nonumber \\&\quad +12 \sqrt{6} \kappa \lambda ^2 \left( y_c^2-1\right) ^2 \left( 2 \alpha ^2+3 \lambda ^2 y_c^2 \left( 2-3 y_c^2\right) \right) \nonumber \\&\quad -4 \kappa ^4 \left( 2 \sqrt{6} \alpha +9 \lambda \left( 2 y_c^2-1\right) \right) \nonumber \\&\quad +144 \lambda ^4 y_c^2 \left( y_c^2-1\right) ^3 \left( \sqrt{6} \alpha +3 \lambda y_c^2\right) +2 \kappa ^3 \nonumber \\&\quad \times \left( 2 \sqrt{6} \alpha ^2+24 \alpha \lambda \left( 3 y_c^2-2\right) +3 \sqrt{6} \lambda ^2 \left( 16 y_c^4-16 y_c^2+3\right) \right) \Bigg ), \end{aligned}$$

(37b)

$$\begin{aligned} v_3&= \frac{1}{24} \kappa \left( 12 x_1+\sqrt{6} \mu \left( y_c^2-1\right) \right) \left( -2 \alpha +\kappa -\sqrt{6} \lambda y_c^2\right) , \end{aligned}$$

(37c)

$$\begin{aligned} u_1&= \frac{1}{3} y_c \left( y_c \left( \left( y_c^2-1\right) \left( \sqrt{6} \lambda x_1+\lambda \mu \left( y_c^2-1\right) +3\right) +3 y_2\right) \nonumber \right. \\&\quad \left. -6 y_1 \left( y_c^2-1\right) \right) , \end{aligned}$$

(37d)

$$\begin{aligned} u_2&=\frac{1}{12} \mu y_c^2 \left( y_c^2-1\right) \left( \sqrt{6} \kappa +6 \lambda \left( 1-2 y_c^2\right) \right) , \end{aligned}$$

(37e)

we obtain evolution equations $u_1', u_2', v_1', v_2'$ and $v_3'$ from the equations that results from taking derivative with respect to $\tau$ of (37) and substituting $x_1', y_1', x_2', y_2'$ and $\mu '$ from equations (23)–(27) and taking inverse transformation of (37). The linear part of the resulting system is written as

$$\begin{aligned} \left( \begin{array}{ccccc} v_1' \\ v_2' \\ v_3' \\ u_1' \\ u_2' \\ \end{array} \right) = \left( \begin{array}{ccccc} -3 &{} 1 &{} 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}{ccccc} v_1 \\ v_2 \\ v_3 \\ u_1 \\ u_2 \\ \end{array} \right) . \end{aligned}$$

(38)

Therefore, the center manifold is given by the graph

$$\begin{aligned} W^{c}_{loc}({\mathbf {0}})=&\Bigg \{(v_1, v_2, v_3,u_1, u_2)\in {\mathbb {R}}^5: v_i=h_i(u_1, u_2), \nonumber \\&h_i(0,0)=0, \nonumber \\&\frac{\partial h_i}{\partial u_j} \Bigg {|}_{(u_1,u_2)=(0,0)}= 0, i=1,2,3, j=1,2, ||(u_1, u_2)||<\delta \Bigg \}, \end{aligned}$$

(39)

for some $\delta >0$.

The $h_i$’s satisfy a set of partial quasi-linear differential equations than can be written symbolically by

$$\begin{aligned}&-u_2' \frac{\partial h_1}{\partial u_2}(u_1,u_2)-u_1' \frac{\partial h_1}{\partial u_1}(u_1,u_2)+v_1'=0, \end{aligned}$$

(40a)

$$\begin{aligned}&-u_2' \frac{\partial h_2}{\partial u_2}(u_1,u_2)-u_1' \frac{\partial h_2}{\partial u_1}(u_1,u_2)+v_2',\end{aligned}$$

(40b)

$$\begin{aligned}&-u_2' \frac{\partial h_3}{\partial u_2}(u_1,u_2)-u_1' \frac{\partial h_3}{\partial u_1}(u_1,u_2)+v_3'=0. \end{aligned}$$

(40c)

where the prime means derivative with respect to $\tau$. These equations are obtained by substituting the expressions $u_1', u_2', v_1', v_2'$ and $v_3'$ from the equations that results from taking derivative with respect to $\tau$ of (37) and substituting $x_1', y_1', x_2', y_2'$ and $\mu '$ from equations (23)–(27) and taking inverse transformation of (37). Next, replacing $v_1\mapsto h_1(u_1, u_2), v_2\mapsto h_2(u_1, u_2), v_3\mapsto h_3(u_1, u_2)$ we obtain the final equations for $h_1(u_1, u_2), \; h_2(u_1, u_2), \; h_3(u_1, u_2)$. In general, we solve these equations by using Taylor expansions in the independent variables $u_1, u_2$ that are of order greater or equal than two. Therefore, we propose the expansion in series

$$\begin{aligned} h_i(u_1, u_2)= & {} \sum _{n=2}^{N}\sum _{k=0}^{n} a^{[i]}_{n k} u_1^{n-k} u_2^{k} \nonumber \\&+ {\mathcal {O}}(\Vert (u_1,u_2)^{N+1}\Vert ),\quad i=1,2, 3. \end{aligned}$$

(41)

Fig. 2

Phase plot of system (42) for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point

Setting $N=2$, we obtain the coefficients $a^{[1]}_{20}= -\frac{\left( y_c^2-1\right) \lambda }{2 \sqrt{6} y_c^2}, \; a^{[2]}_{20}= \frac{-\sqrt{6} \lambda y_c^2-2 \alpha +\kappa }{8 y_c^2}, \; a^{[3]}_{20}= 0, \; a^{[3]}_{21}= -\frac{\kappa \left( -6 \lambda y_c^2-2 \sqrt{6} \alpha +\sqrt{6} \kappa \right) }{2 y_c^2 \left( y_c^2-1\right) \left( \sqrt{6} \kappa +6 \left( 1-2 y_c^2\right) \lambda \right) },$ and $a^{[1]}_{21}, a^{[2]}_{21}, a^{[1]}_{22}, a^{[2]}_{22}$ and $a^{[3]}_{22}$ are complicated expressions of the parameters.

Finally, the dynamics on the center manifold up to third order is given by a system of the form

$$\begin{aligned} u_1'&= \alpha _1 {u_2^2}+ \alpha _2 {u_1 u_2}+u_2 + {\mathcal {O}}(\Vert (u_1,u_2)^{3}\Vert ), \end{aligned}$$

(42a)

$$\begin{aligned} u_2'&= -\frac{u_2^2 \left( \sqrt{6} \kappa -6 y_c^2 \lambda \right) }{2 y_c^2 \left( y_c^2-1\right) \left( \sqrt{6} \kappa +6 \left( 1-2 y_c^2\right) \lambda \right) } \nonumber \\&\quad + {\mathcal {O}}(\Vert (u_1,u_2)^{3}\Vert ), \end{aligned}$$

(42b)

where $\alpha _1$ and $\alpha _2$ are complicated functions of the parameters.

From Eq. (42b) the center is unstable to perturbations along the $u_2$-axis (saddle behavior). This is due to (42b) which is a gradient like equation

$$\begin{aligned}&u_2'=-U'(u_2), \quad \nonumber \\&U(u_2)= \frac{u_2^3 \left( \sqrt{6} \kappa -6 y_c^2 \lambda \right) }{6 y_c^2 \left( y_c^2-1\right) \left( \sqrt{6} \kappa +6 \left( 1-2 y_c^2\right) \lambda \right) }, \end{aligned}$$

(43)

such that $u_2=0$ is an inflection point of $U(u_2)$. That is, the solution either depart or approach the origin along the $u_2$-axis, depending on the sign of the initial value of $u_2$, $u_{20}$ and the sign of $\frac{\left( \sqrt{6} \kappa -6 y_c^2 \lambda \right) }{ y_c^2 \left( y_c^2-1\right) \left( \sqrt{6} \kappa +6 \left( 1-2 y_c^2\right) \lambda \right) }$. The system also admits an invariant line of stationary points which at second order is approximated by $\alpha _1 u_2 + \alpha _2 u_1 +1=0$ which behaves as saddle.

We represent the phase plot of the system (42) for several choices of parameters in Fig. 2 where the origin is a saddle point.

Now, we study the particular cases $(\sqrt{6} \kappa + 6 (1 - 2 y_c^2) \lambda )= 0$, or $y_c=0$, or $y_c^2=1$.

(a) Subcase $\lambda =\frac{\kappa }{\sqrt{6} \left( 2 y_c^2-1\right) }$, $y_c^2\ne \frac{1}{2}$. In this subcase the similarity matrix is

$$\begin{aligned} s= \left( \begin{array}{ccccc} -\frac{y_c^2-1}{2 \sqrt{6}} &{} 0 &{} 0 &{} 0 &{} \frac{1}{72 \kappa \left( 2 y_c^2-1\right) ^3 \left( \alpha \left( 2-4 y_c^2\right) +\kappa \left( y_c^2-1\right) \right) } \\ \frac{\kappa y_c \left( y_c^2-1\right) }{12 \sqrt{6} \left( 2 y_c^2-1\right) } &{} \frac{1}{2 y_c} &{} 0 &{} \frac{y_c}{12 \left( 2 y_c^2-1\right) \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right) } &{} \frac{\kappa y_c \left( y_c^2-1\right) ^2}{216 \left( 1-2 y_c^2\right) ^4 \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right) ^2} \\ \frac{\left( y_c^2-1\right) \left( 2 \alpha \kappa +(y_c-1) (y_c+1) \left( 3 \kappa ^2+8 y_c^2 (\kappa (\alpha -\kappa )+18)\right) +36\right) }{36 \sqrt{6} \left( 1-2 y_c^2\right) ^2} &{} \frac{\kappa }{6 y_c^2-3} &{} 1 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} \frac{y_c^2-1}{6 \left( 2 y_c^2-1\right) \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right) } &{} \frac{\left( y_c^2-1\right) \left( \alpha \left( 4 y_c^2-2\right) +\kappa \left( 2 y_c^4-5 y_c^2+3\right) \right) }{216 \left( 1-2 y_c^2\right) ^4 \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right) ^2} \\ 1 &{} 0 &{} 0 &{} 0 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(44)

and the real Jordan form of the matrix is

$$\begin{aligned} j(A_{2}^{^{\prime }})= \left( \begin{array}{ccccc} 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} -3 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} -3 &{} 1 \\ 0 &{} 0 &{} 0 &{} 0 &{} -3 \\ \end{array} \right) . \end{aligned}$$

(45)

Defining

$$\begin{aligned}&(u_1, u_2, v_1, v_2, v_3 )\\&\quad = s^{-1} \cdot \left( x_1,y_1-y_c,x_2+\frac{1}{3} \left( \kappa -\frac{\kappa y_c^2}{2 y_c^2-1}\right) ,\right. \\&\left. \qquad y_2-y_c^2+1,\mu \right) , \end{aligned}$$

and using the CMT we obtain that the center manifold is given by (39), where the $h_i$ satisfy (40).

Using the ansatz (41) with $N=2$, we obtain the coefficients

$a^{[1]}_{20}= -\frac{\kappa \left( y_c^2-1\right) ^2 \left( -4 \alpha ^2+12 \alpha \kappa -9 \kappa ^2+8 y_c^6 \left( 4 \alpha ^2-11 \alpha \kappa +7 \kappa ^2-72\right) +y_c^4 \left( -48 \alpha ^2+136 \alpha \kappa -101 \kappa ^2+720\right) +y_c^2 \left( 24 \alpha ^2-70 \alpha \kappa +54 \kappa ^2-288\right) +36\right) }{2592 \left( 2 y_c^2-1\right) ^3}$,

$a^{[2]}_{20}= \frac{\left( y_c^2-1\right) ^2 \left( 144 y_c^4+\left( \kappa ^2-144\right) y_c^2+36\right) \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right) }{144 \left( 2 y_c^2-1\right) }, a^{[3]}_{20}= 0, a^{[1]}_{21}= \frac{2 \alpha \kappa -\kappa ^2+y_c^4 \left( 8 \alpha \kappa +3 \left( \kappa ^2+48\right) \right) -2 y_c^2 \left( 4 \alpha \kappa +\kappa ^2+72\right) +36}{36 \sqrt{6} \left( 1-2 y_c^2\right) ^2}$,

$a^{[2]}_{21} = \frac{\kappa \left( \alpha \left( 20 y_c^4-22 y_c^2+6\right) +\kappa \left( -9 y_c^4+16 y_c^2-7\right) \right) }{2 \sqrt{6}}, a^{[3]}_{21} =6 \sqrt{6} \kappa \left( 2 y_c^2\right. \left. -1\right) ^3 \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right)$,

$a^{[1]}_{22}= -\frac{\kappa -\kappa y_c^2}{12 y_c^2-24 y_c^4}, a^{[2]}_{22}= -\frac{3 \left( 2 y_c^2-1\right) \left( \kappa +\alpha \left( 4 y_c^2-2\right) -\kappa y_c^2\right) }{2 y_c^2}, a^{[3]}_{22}= 0$.

Fig. 3

Phase plot of system (46) for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point

Therefore, the dynamics on the center manifold is governed up to third order by

$$\begin{aligned} u_1'&= \frac{\sqrt{6} \kappa u_1^2 \left( y_c^2-1\right) }{24(1-2 y_c^2)}, \end{aligned}$$

(46a)

$$\begin{aligned} u_2'&=-\frac{u_1^2 \left( y_c^2-1\right) ^2 y_c^2 \left( -\alpha \kappa +2 \left( \kappa ^2-9\right) +2 y_c^2 (\kappa (\alpha -\kappa )+18)\right) }{72 \left( 2 y_c^2-1\right) } \nonumber \\&\quad -\frac{\kappa u_1 u_2 \left( y_c^2-1\right) y_c^2}{\sqrt{6} \left( 2 y_c^2-1\right) }. \end{aligned}$$

(46b)

From Eq. (46a) the center is unstable to perturbations along the $u_1$-axis (saddle behavior). This is due to (46a) which is a gradient like equation

$$\begin{aligned} u_1'=-U'(u_1), \quad U(u_1)=- \frac{\sqrt{6} \kappa u_1^3 \left( y_c^2-1\right) }{72(1-2 y_c^2)}, \end{aligned}$$

(47)

such that $u_1=0$ is an inflection point of $U(u_1)$. That is, the solution either depart or approach the origin along the $u_1$-axis, depending on the sign of the initial value of $u_1$, $u_{10}$ and the sign of $\frac{\kappa \left( y_c^2-1\right) }{(1-2 y_c^2)}$.

In Fig. 3 a phase plot of system (46) is presented for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point.

(b) Subcase $\kappa =0, y_c^2= \frac{1}{2}$.

For $y_c= \frac{\sqrt{2}}{2}$ the similarity matrix is

$$\begin{aligned} s=\left( \begin{array}{ccccc} \frac{1}{\frac{\lambda ^2}{4 \alpha +\sqrt{6} \lambda }+\frac{\lambda }{\sqrt{6}}} &{} 0 &{} \frac{1}{\sqrt{\frac{2}{3}} \alpha \lambda +\lambda ^2} &{} \frac{1}{4 \sqrt{6}} &{} 0 \\ -\frac{3 \lambda }{4 \sqrt{3} \alpha +6 \sqrt{2} \lambda } &{} 0 &{} \frac{1}{2 \sqrt{2} \alpha +2 \sqrt{3} \lambda } &{} -\frac{\lambda }{24 \sqrt{2}} &{} \frac{1}{\sqrt{2}} \\ 0 &{} 1 &{} 0 &{} -\frac{\lambda ^2+12}{24 \sqrt{6}} &{} \sqrt{\frac{2}{3}} \lambda \\ 1 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(48)

and the real Jordan matrix is

$$\begin{aligned} j(A_{2}^{^{\prime }})= \left( \begin{array}{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) . \end{aligned}$$

(49)

In this case we define

$$\begin{aligned} v_1&=\frac{1}{24} \left( -\lambda \mu +2 \sqrt{6} \lambda x_1-12 \sqrt{2} y_1+12 y_2+18\right) , \\ v_2&= \frac{1}{144} \left( -\sqrt{6} \left( \lambda ^2-12\right) \mu \right. \\&\quad \left. -12 \lambda \left( -2 \lambda x_1+4 \sqrt{3} y_1+2 \sqrt{6} y_2+\sqrt{6}\right) \right) +x_2,\\ v_3&=\alpha \left( \frac{\lambda x_1}{\sqrt{6}}+\sqrt{2} y_1-y_2-\frac{3}{2}\right) \\&\quad -\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) ,\\ u_1&= \mu ,\\ 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}$$

On the other hand, for $y_c=- \frac{\sqrt{2}}{2}$ the similarity matrix changes to

$$\begin{aligned} s=\left( \begin{array}{ccccc} \frac{1}{\frac{\lambda ^2}{4 \alpha +\sqrt{6} \lambda }+\frac{\lambda }{\sqrt{6}}} &{} 0 &{} \frac{1}{\sqrt{\frac{2}{3}} \alpha \lambda +\lambda ^2} &{} \frac{1}{4 \sqrt{6}} &{} 0 \\ \frac{3 \lambda }{4 \sqrt{3} \alpha +6 \sqrt{2} \lambda } &{} 0 &{} \frac{1}{-2 \sqrt{2} \alpha -2 \sqrt{3} \lambda } &{} \frac{\lambda }{24 \sqrt{2}} &{} -\frac{1}{\sqrt{2}} \\ 0 &{} 1 &{} 0 &{} -\frac{\lambda ^2+12}{24 \sqrt{6}} &{} \sqrt{\frac{2}{3}} \lambda \\ 1 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(51)

and the real Jordan matrix is (49). In this case we define

$$\begin{aligned} v_1&=\frac{1}{24} \left( -\lambda \mu +2 \sqrt{6} \lambda x_1+12 \sqrt{2} y_1+12 y_2+18\right) ,\\ v_2&= \frac{1}{144} \left( -\sqrt{6} \left( \lambda ^2-12\right) \mu \right. \\&\quad \left. -12 \lambda \left( -2 \lambda x_1-4 \sqrt{3} y_1+2 \sqrt{6} y_2+\sqrt{6}\right) \right) +x_2,\\ v_3&= -\frac{1}{6} \alpha \left( -\sqrt{6} \lambda x_1+6 \sqrt{2} y_1+6 y_2+9\right) \\&\quad -\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) ,\\ u_1&= \mu , \\ 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}$$

The real Jordan matrix is the same as before. In both cases the CMT gives the same result. That is, the center manifold is given by (39), where the $h_i$ satisfy (40).

Fig. 4

Phase plot of system (53) for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point

Using the ansatz (41) with $N=2$, we obtain the coefficients $a^{[1]}_{20}= \frac{-\lambda ^2-12}{2304}, a^{[2]}_{20}=\frac{\lambda \left( \lambda ^2+24\right) }{1152 \sqrt{6}}, a^{[3]}_{20}= \frac{\left( \lambda ^2+12\right) \left( 4 \alpha +\sqrt{6} \lambda \right) }{4608}, a^{[1]}_{21}= -\frac{\lambda }{16}, a^{[2]}_{21}= \frac{24-5 \lambda ^2}{24 \sqrt{6}}, a^{[3]}_{21}= -\frac{1}{96} \lambda \left( 4 \alpha +5 \sqrt{6} \lambda \right) , a^{[1]}_{22}=\frac{1}{4} , a^{[2]}_{22}=\frac{\lambda }{2 \sqrt{6}}, a^{[3]}_{22}=\frac{1}{8} \left( -4 \alpha -\sqrt{6} \lambda \right)$.

Therefore, the dynamics on the center manifold is governed up to third order by

$$\begin{aligned}&u_1'= \frac{\lambda u_1^2}{8}, \end{aligned}$$

(53a)

$$\begin{aligned}&u_2'=\frac{\lambda u_1 u_2}{4}-\frac{u_1^2}{32}. \end{aligned}$$

(53b)

From Eq. (53a)which the center is unstable to perturbations along the $u_1$-axis (saddle behavior). This is due to (53a) which is a gradient like equation

$$\begin{aligned} u_1'=-U'(u_1), \quad U(u_1)= -\frac{\lambda u_1^3}{24}, \end{aligned}$$

(54)

such that $u_1=0$ is an inflection point of $U(u_1)$. That is, the solution either depart or approach the origin along the $u_1$-axis, depending on the sign of the initial value of $u_1$, $u_{10}$ and the sign of $\lambda$.

In Fig. 4 a phase plot of system (53) is presented for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point.

(c) Subcase $y_c=0$, $\kappa \ne 2\alpha$. In this subcase the similarity matrix is

$$\begin{aligned} s= \left( \begin{array}{ccccc} 0 &{} 0 &{} -\frac{2}{2 \alpha \kappa -\kappa ^2} &{} \frac{1}{2 \sqrt{6}} &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 1 &{} 0 &{} 0 &{} -\frac{2 \alpha \kappa -3 \kappa ^2+36}{36 \sqrt{6}} &{} 0 \\ 0 &{} \frac{2}{\kappa -2 \alpha } &{} \frac{4 \alpha -6 \kappa }{3 (\kappa -2 \alpha )^2} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(55)

and the real Jordan matrix is

$$\begin{aligned} j(A_{2}^{^{\prime }})= \left( \begin{array}{ccccc} -3 &{} 1 &{} 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) . \end{aligned}$$

(56)

Defining

$$\begin{aligned} v_1&=\frac{\mu \left( 2 \alpha \kappa -3 \kappa ^2+36\right) }{36 \sqrt{6}}+\frac{\kappa }{3}+x_2, \\ v_2&= \frac{1}{72} \left( 3 \kappa \left( -\sqrt{6} \kappa \mu +12 \kappa x_1+12 y_2+12\right) \right. \\&\quad \left. -2 \alpha \left( -\sqrt{6} \kappa \mu +12 \kappa x_1+36 y_2+36\right) \right) ,\\ v_3&= -\frac{1}{24} \kappa (2 \alpha -\kappa ) \left( 12 x_1-\sqrt{6} \mu \right) , \quad u_1= \mu , \quad u_2=y_1, \end{aligned}$$

and using the CMT we obtain that the center manifold is given by (39), where the $h_i$ satisfy (40).

Using the ansatz (41) with $N=2$, we obtain the coefficients $a^{[1]}_{20}= \frac{\kappa \left( -4 \alpha ^2+12 \alpha \kappa -9 \kappa ^2+36\right) }{2592}, a^{[2]}_{20}= \frac{1}{48} (2 \alpha -\kappa ), a^{[3]}_{20}=0, a^{[1]}_{21}=0, a^{[2]}_{21}=0, a^{[3]}_{21}= 0, a^{[1]}_{22}= \sqrt{\frac{2}{3}} \lambda , a^{[2]}_{22}=\frac{1}{2} (\kappa -2 \alpha ), a^{[3]}_{22}=0$. Therefore, the dynamics on the center manifold is governed up to third order by

$$\begin{aligned}&u_1'= -\frac{\kappa u_1^2}{4 \sqrt{6}}, \end{aligned}$$

(58a)

$$\begin{aligned}&u_2'= \frac{1}{24} u_1 u_2 \left( -\sqrt{6} \kappa -6 \lambda \right) . \end{aligned}$$

(58b)

From Eq. (58a) the center is unstable to perturbations along the $u_1$-axis (saddle behavior). This is due to (58a) which is a gradient like equation

$$\begin{aligned} u_1'=-U'(u_1), \quad U(u_1)=\frac{\kappa u_1^3}{12 \sqrt{6}}, \end{aligned}$$

(59)

such that $u_1=0$ is an inflection point of $U(u_1)$. That is, the solution either depart or approach the origin along the $u_1$-axis, depending on the sign of the initial value of $u_1$, $u_{10}$ and the sign of $\kappa$.

Fig. 5

Phase plot of system (58) for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point

In Fig. 5 a phase plot of system (58) is presented for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point.

(d) Subcase $y_c=0$, $\kappa = 2\alpha$.

In this case, the similarity matrix is

$$\begin{aligned} s= \left( \begin{array}{ccccc} 0 &{} 0 &{} 1 &{} \frac{1}{2 \sqrt{6}} &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ -\frac{\kappa }{3} &{} \frac{\kappa ^2}{3} &{} 0 &{} \frac{\kappa ^2-18}{18 \sqrt{6}} &{} 0 \\ 0 &{} \kappa &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(60)

and the real Jordan matrix is (49).

Defining

$$\begin{aligned}&v_1= \frac{\sqrt{6} \left( \kappa ^2-18\right) \mu -108 x_2+36 \kappa y_2}{36 \kappa }, \\&v_2= \frac{y_2+1}{\kappa }, \quad v_3= x_1-\frac{\mu }{2 \sqrt{6}}, \quad u_1= \mu , \quad u_2= y_1 \end{aligned}$$

and using the CMT we obtain that the center manifold is given by (39), where the $h_i$ satisfy (40).

Using the ansatz (41) with $N=2$, we obtain the coefficients $a^{[1]}_{20}= \frac{1}{216} \left( \kappa ^2-18\right) , a^{[2]}_{20}= -\frac{1}{24 \kappa }, a^{[3]}_{20}= 0,a^{[1]}_{21}= 0, a^{[2]}_{21}= 0, a^{[3]}_{21}= 0, a^{[1]}_{22}= 1-\frac{\sqrt{6} \lambda }{\kappa }, a^{[2]}_{22}= \frac{1}{\kappa }, a^{[3]}_{22}= 0$. Hence, we obtain the system governing the dynamics on the center manifold is (58). Then, it follows the same result as before. That is, $A_{2}^{^{\prime }}$ is a saddle point.

(e) Subcase $y_c^2=1$, $2 \alpha -\kappa +\sqrt{6} \lambda \ne 0$. In this subcase the similarity matrix is

$$\begin{aligned} s= \left( \begin{array}{ccccc} 0 &{} 0 &{} \frac{2}{\kappa \left( -2 \alpha +\kappa -\sqrt{6} \lambda \right) } &{} 0 &{} 0 \\ 0 &{} \frac{\epsilon }{2 \alpha -\kappa +\sqrt{6} \lambda } &{} \frac{\epsilon \left( \kappa -\sqrt{6} \lambda \right) \left( -2 \alpha +3 \kappa -\sqrt{6} \lambda \right) }{3 \kappa \left( 2 \alpha -\kappa +\sqrt{6} \lambda \right) ^2} &{} 0 &{} \frac{\epsilon }{2} \\ 1 &{} 0 &{} 0 &{} 0 &{} \sqrt{\frac{2}{3}} \lambda \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(62)

where $\epsilon =\pm 1$ is the value of $y_c$ and the real Jordan matrix is the same as (56).

Defining

$$\begin{aligned}&v_1=\frac{1}{3} \left( \kappa +3 x_2-\sqrt{6} \lambda (y_2+1)\right) , \\&v_2= \frac{1}{6} \left( 6 \lambda ^2 x_1-2 \alpha \left( \kappa x_1-\sqrt{6} \lambda x_1-6 y_1 \epsilon +3 y_2+6\right) \right. \\&\quad \left. -\sqrt{6} \lambda (4 \kappa x_1-6 y_1 \epsilon +3 y_2+6)+3 \kappa (\kappa x_1-2 y_1 \epsilon +y_2+2)\right) ,\\&v_3= \frac{1}{2} \kappa x_1 \left( -2 \alpha +\kappa -\sqrt{6} \lambda \right) ,\quad u_1= \mu , \quad u_2=y_2, \end{aligned}$$

and using the CMT we obtain that the center manifold is given by (39), where the $h_i$ satisfy (40).

Using the ansatz (41) with $N=2$, we obtain the coefficients $a^{[1]}_{20}= 0, a^{[2]}_{20}= 0, a^{[3]}_{20}= 0, a^{[1]}_{21}= \frac{1}{216} \left( 2 \sqrt{6} \alpha \kappa -3 \sqrt{6} \kappa ^2+6 \kappa \lambda +12 \sqrt{6} \left( \lambda ^2+3\right) \right) , a^{[2]}_{21}= \frac{1}{72} \kappa \left( 2 \sqrt{6} \alpha -3 \sqrt{6} \kappa +18 \lambda \right) , a^{[3]}_{21}= \frac{1}{24} \kappa \left( 2 \sqrt{6} \alpha -\sqrt{6} \kappa \right. \left. +6 \lambda \right) , a^{[1]}_{22}=0,a^{[2]}_{22}= \frac{1}{8} \left( -2 \alpha +\kappa -\sqrt{6} \lambda \right) , a^{[3]}_{22}= 0$. Therefore, the dynamics on the center manifold is governed up to third order by

$$\begin{aligned}&u_1'= -\frac{1}{24} u_1^2 \left( \sqrt{6} \kappa -6 \lambda \right) , \end{aligned}$$

(64a)

$$\begin{aligned}&u_2'= \frac{1}{12} u_1 u_2 \left( \sqrt{6} \kappa -6 \lambda \right) . \end{aligned}$$

(64b)

Fig. 6

Phase plot of system (64) for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point

From Eq. (64a) the center is unstable to perturbations along the $u_1$-axis. This is due to (64a) which is a gradient like equation (saddle behavior)

$$\begin{aligned} u_1'=-U'(u_1), \quad U(u_1)=\frac{1}{72} u_1^2 \left( \sqrt{6} \kappa -6 \lambda \right) , \end{aligned}$$

(65)

such that $u_1=0$ is an inflection point of $U(u_1)$. That is, the solution either depart or approach the origin along the $u_1$-axis, depending on the sign of the initial value of $u_1$, $u_{10}$ and the sign of $\sqrt{6} \kappa -6 \lambda$.

In Fig. 6 a phase plot of system (58) is presented for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point.

(f) Subcase $y_c^2=1$, $2 \alpha -\kappa +\sqrt{6} \lambda =0, \alpha \ne 0$. In this subcase the similarity matrix is

$$\begin{aligned} s= \left( \begin{array}{ccccc} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ \frac{3 \epsilon }{2 \left( \kappa -\sqrt{6} \lambda \right) } &{} -\frac{\kappa \epsilon }{2} &{} 0 &{} 0 &{} \frac{\epsilon }{2} \\ 0 &{} \frac{1}{3} \kappa \left( \kappa -\sqrt{6} \lambda \right) &{} 0 &{} 0 &{} \sqrt{\frac{2}{3}} \lambda \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(66)

where $\epsilon =\pm 1$ is the value of $y_c$ and the real Jordan matrix is the same as (49).

Defining

$$\begin{aligned}&v_1= \frac{1}{3} \left( 3 x_2-\kappa (-2 y_1 \epsilon +y_2+1)+\sqrt{6} \lambda (1-2 y_1 \epsilon )\right) ,\\&v_2= \frac{\kappa +3 x_2-\sqrt{6} \lambda (y_2+1)}{\kappa \left( \kappa -\sqrt{6} \lambda \right) }, \quad v_3=x_1, \quad u_1= \mu , \quad u_2= y_2, \end{aligned}$$

and using the CMT we obtain that the center manifold is given by (39), where the $h_i$ satisfy (40).

Using the ansatz (41) with $N=2$, we obtain the nonzero coefficients $a^{[1]}_{20}= 0, a^{[2]}_{20}= 0, a^{[3]}_{20}= 0, a^{[1]}_{21}= \frac{1}{\sqrt{6}}-\frac{1}{108} \kappa \left( \sqrt{6} \kappa -6 \lambda \right) , a^{[2]}_{21}= \frac{-\kappa ^2+6 \lambda ^2+18}{6 \sqrt{6} \kappa ^2-36 \kappa \lambda }, a^{[3]}_{21}= -\frac{1}{2 \sqrt{6}}, a^{[1]}_{22} = \frac{1}{12} \left( \sqrt{6} \lambda -\kappa \right) , a^{[2]}_{22}= 0, a^{[3]}_{22}= 0$. Therefore, the dynamics on the center manifold is governed up to third order by (64) and we have the same result as before that $A_{2}^{^{\prime }}$ is a saddle point.

(g) Subcase $y_c^2=1$, $\kappa =\sqrt{6} \lambda , \alpha = 0$.

Fig. 7

Phase plot of system (70) for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point

In this subcase the similarity matrix is

$$\begin{aligned} s= \left( \begin{array}{ccccc} 0 &{} 0 &{} -\frac{\sqrt{\frac{2}{3}} \epsilon }{\lambda } &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 &{} \frac{\epsilon }{2} \\ 1 &{} 0 &{} 0 &{} 0 &{} \sqrt{\frac{2}{3}} \lambda \\ 0 &{} 0 &{} 0 &{} 0 &{} 1 \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 \\ \end{array} \right) , \end{aligned}$$

(68)

where $\epsilon =\pm 1$ is the value of $y_c$ and the real Jordan matrix is the same as (49).

Defining

$$\begin{aligned}&v_1=x_2-\sqrt{\frac{2}{3}} \lambda y_2, \quad v_2= y_1-\frac{1}{2} (y_2+2) \epsilon , \\&v_3=-\sqrt{\frac{3}{2}} \lambda x_1 \epsilon ,\quad u_1=\mu , \quad u_2= y_2, \end{aligned}$$

and using the CMT we obtain that the center manifold is given by (39), where the $h_i$ satisfy (40).

Using the ansatz (41) with $N=3$, we obtain the coefficients $a^{[1]}_{20}=0, a^{[2]}_{20}=0, a^{[3]}_{20}= 0, a^{[1]}_{21}=\frac{1}{\sqrt{6}}, a^{[2]}_{21}= \frac{\lambda \epsilon }{12}, a^{[3]}_{21}=\frac{\lambda \epsilon }{4}, a^{[1]}_{22}= 0, a^{[2]}_{22}= -\frac{\epsilon }{8}, a^{[1]}_{22}= 0, a^{[1]}_{30}=0, a^{[2]}_{30}= 0, a^{[3]}_{30}= 0, a^{[1]}_{31}= 0, a^{[2]}_{31}=0, a^{[3]}_{31}= 0, a^{[1]}_{32}= \frac{5 \lambda ^2}{6 \sqrt{6}}, a^{[2]}_{32}= \frac{\lambda \epsilon }{24}, a^{[3]}_{32}=0, a^{[1]}_{33}=0, a^{[2]}_{33}=\frac{\epsilon ^3}{16}, a^{[3]}_{33}= 0$.

Therefore, the dynamics on the center manifold is governed up to fourth order by

$$\begin{aligned} u_1'= \frac{1}{4} \lambda u_1^2 u_2, \quad u_2'=- \lambda u_1 u_2^2. \end{aligned}$$

(70)

In Fig. 7 a phase plot of system (70) is presented for some values of the parameters showing that $A_{2}^{^{\prime }}$ is a saddle point.

5 Cosmological linear perturbations

We continue our analysis by investigating the modified scalar field Lagrangian’s effects in the cosmological linear perturbation theory. For this purpose, we consider the simple linear perturbation theory in the Newtonian gauge, where the spacetime line element is described as [85]

$$\begin{aligned} ds^{2}= & {} a^{2}\left( \eta \right) \left( -\left( 1+2\Phi \left( \eta ,x,y,z\right) \right) d\eta ^{2} \right. \nonumber \\&\left. +\left( 1-2\Phi \left( \eta ,x,y,z\right) \right) \left( dx^{2}+dy^{2}+dz^{2}\right) \right) ,\nonumber \\ \end{aligned}$$

(71)

where $\eta$  ($= \int a^{-1} dt$) is the conformal time .

Before presenting the perturbed equations, let us define the perturbed physical quantities of the usual scalar field $\phi$, the scalar field $\psi$ due to the effect of GUP, the matter-energy density, the DE density, the dark energy pressure as follows:

$$\begin{aligned} \phi +\delta \phi ,\quad \psi +\delta \psi ,\quad \rho _m+\delta \rho _m,\quad \rho _d+\delta \rho _d,\quad p_d+\delta p_d.\nonumber \\ \end{aligned}$$

(72)

We note here that $\phi , \psi , \rho _m, \rho _d, p_d$ denote the corresponding background quantities.

For the line element (71), the perturbed terms of the Einstein tensor are [85]

$$\begin{aligned} \delta G_{0}^{0}&=-\frac{2}{a^{2}}\nabla ^{2}\Phi +\frac{6}{a} {\mathcal {H}}\left( {\dot{\Phi }}+a{\mathcal {H}}\Phi \right) , \end{aligned}$$

(73)

$$\begin{aligned} \delta G_{i}^{0}&=-\frac{2}{a^{2}}\nabla _{i}\left( \dot{\Phi }+a{\mathcal {H}}\Phi \right) , \end{aligned}$$

(74)

$$\begin{aligned} \delta G_{j}^{i}&=2\Phi \left( \frac{2}{a}\mathcal {{\dot{H}}}+3{\mathcal {H}} ^{2}\right) +\frac{2}{a^{2}}\left( \ddot{\Phi }+3a{\mathcal {H}}\dot{\Phi }\right) , \end{aligned}$$

(75)

where now dot means the derivative with respect to the variable $\eta$, and ${\mathcal {H}}=\frac{{\dot{a}}}{a^{2}}$ is the Hubble parameter in the new frame.

We assume the comoving observer $u^{\mu }=\left( \frac{1-\Phi }{a},\delta u^{i}\right)$ such that the energy momentum tensor for the dark matter to have nonzero perturbed components [85]

$$\begin{aligned} \delta T_{0}^{\left( m\right) 0}=-\rho _{m}\delta _{m},\quad \delta T_{0} ^{i}=a\rho _{m}\delta u^{i}, \end{aligned}$$

(76)

where $\delta _{m}=\delta \rho _{m}/\rho _{m}$. Similarly for the scalar field it follows

$$\begin{aligned} \delta T_{0}^{\left( d\right) 0}= & {} -\frac{1}{a^{2}}\left( {\dot{\phi }} \dot{\left( \delta \phi \right) } -{\dot{\phi }}^{2}\Phi \right) -V_{,\phi }\delta \phi \nonumber \\&-\frac{2\beta \hbar ^{2}}{a^{2}}\left( \dot{\left( \delta \phi \right) } {\dot{\psi }}-2{\dot{\phi }}{\dot{\psi }}\Phi \right) \nonumber \\&-\frac{2\beta \hbar ^{2}}{a^{2} }{\dot{\phi }} \dot{\left( \delta \psi \right) } +2\beta \hbar ^{2}\psi \delta \psi , \end{aligned}$$

(77)

$$\begin{aligned} \delta T_{i}^{\left( d\right) i}= & {} \frac{1}{a^{2}}\left( {\dot{\phi }}\dot{\left( \delta \phi \right) } -{\dot{\phi }}^{2}\Phi \right) -V_{,\phi }\delta \phi \nonumber \\&+\frac{2\beta \hbar ^{2}}{a^{2}}\left( \dot{\left( \delta \phi \right) } {\dot{\psi }}-2{\dot{\phi }}{\dot{\psi }}\Phi \right) \nonumber \\&+\frac{2\beta \hbar ^{2}}{a^{2} }{\dot{\phi }} \dot{\left( \delta \psi \right) } +2\beta \hbar ^{2}\psi \delta \psi , \end{aligned}$$

(78)

$$\begin{aligned} \delta T_{i}^{\left( d\right) 0}= & {} -\frac{{\dot{\phi }}}{a^{2}}\left( \delta \phi +4\beta \hbar ^{2}\delta \psi \right) _{,i}, \end{aligned}$$

(79)

where we have replaced $\phi \rightarrow \phi +\delta \phi$ and $\psi \rightarrow \psi +\delta \psi$. However, $\delta \phi ~$and $\delta \psi$ are not independent, they are related with the constraint condition (20). At this point, someone could omit the terms $\beta \hbar ^{2}\delta \phi$ ,$~\beta \hbar ^{2}\delta \psi$, etc., because they are small. However, it is important to remark that in terms of $\beta \hbar ^{2}$, the related dynamical system is a singular perturbation system, and as it has been shown before in [33] that the singular perturbation term can be relatively large such that it can drive the dynamics in the background space. That is why we found the existence of de Sitter asymptotic universe for the background space in the previous section.

We now proceed with our analysis by considering the interacting model of previous section with $\kappa =0$, that is $Q_{A}=\frac{\alpha }{\sqrt{6}a}\left( {\dot{\phi }}\delta _{m}+\rho _{m}\delta {\dot{\phi }}\right)$. Therefore, the perturbative equations for the matter components are

$$\begin{aligned} {\dot{\delta }}_{m}-3{\dot{\Phi }}+a( \delta u^{i}) _{,i}-\frac{\alpha }{\sqrt{6}}\left( \Phi {\dot{\phi }}+\delta {\dot{\phi }}\right) =0, \end{aligned}$$

(80)

where we have replaced $Q=Q_{0}+\delta Q$, in which $Q_{0}$ is the unperturbed term and $\delta Q$ is the linear perturbation term.

For the scalar field it follows

$$\begin{aligned}&\frac{1}{a^{2}}\left( \ddot{(\delta \phi )}-\nabla ^{2}\left( \delta \phi \right) -4{\dot{\Phi }}{\dot{\phi }}\right) +\frac{2}{a}{\mathcal {H}}\delta {\dot{\phi }} -2\Phi \psi -\delta \psi =0,\nonumber \\ \end{aligned}$$

(81)

$$\begin{aligned}&\frac{2\beta \hbar ^{2}}{a^{2}}\left( \left( \delta \ddot{\psi }\right) -\nabla ^{2}\left( \delta \psi \right) -{\dot{\psi }}{\dot{\Phi }}+2a{\mathcal {H}} \dot{\left( \delta \psi \right) }+2a{\mathcal {H}}{\dot{\psi }}\Phi \right) \nonumber \\&\quad +\left( \delta \psi +V_{,\phi \phi }\delta \phi \right) +2\left( \psi +V_{,\phi }\right) \Phi +\frac{2\alpha }{\sqrt{6}}\Phi \rho _{m} \nonumber \\&\quad +\frac{\alpha }{\sqrt{6}}\rho _{m}\delta _{m}=0\,, \end{aligned}$$

(82)

where we have omitted the presentation of the equations for the velocity divergence.

From Eqs. (81) and (82) it is clear that for $\beta \hbar ^{2}=0$, the equations for quintessence are recovered [49]. However, from Eq. (82), we observe that it is singular in terms of $\beta \hbar ^{2},$ which means that the singular perturbative term $\left( \left( \delta \ddot{\psi }\right) -\nabla ^{2}\left( \delta \psi \right) +\frac{2}{a}{\mathcal {H}}{\dot{\psi }}\Phi -4{\dot{\Phi }}\dot{\psi }\right)$ can be significantly large such that it can dominate and drive the dynamics.

We continue by define the new independent variable to be $\tau =\ln a$, that is $\frac{1}{a}\frac{d}{d\eta }={\mathcal {H}}\frac{d}{d\ln a}$. Hence, by writing in Fourier form and setting $\delta u^{i}=0$, for the exponential potential, the perturbative equations (80)–(82) are written as follows

$$\begin{aligned}&0={\delta }_{m}^{\prime }-3\Phi ^{\prime }-\frac{\alpha }{\sqrt{6}}\left( \delta \phi ^{\prime }+2\sqrt{6}x_{1}\Phi \right) , \end{aligned}$$

(83)

$$\begin{aligned}&0=\delta \phi ^{\prime \prime }+\left( \frac{{\mathcal {H}}^{\prime }}{{\mathcal {H}}}+4\right) \delta \phi ^{\prime } \nonumber \\&\quad -4\sqrt{6}x_{1}\Phi ^{\prime } +\frac{k^{2}}{a^2 {\mathcal {H}}^{2}}\delta \phi - \frac{\delta \psi +3\mu y_{2}\Phi }{{\mathcal {H}}^{2}}, \end{aligned}$$

(84)

$$\begin{aligned}&0=2\beta \hbar ^{2}\left( \delta \psi ^{\prime \prime }+\left( \frac{{\mathcal {H}}^{\prime }}{{\mathcal {H}}}+4\right) \delta \psi ^{\prime }+\frac{k^{2} }{a^{2}{\mathcal {H}}^{2}}\delta \psi \right) \nonumber \\&\quad +\frac{\delta \psi +3\mu y_{2}\Phi }{{\mathcal {H}}^{2}}+3y_{1}^{2}\lambda \left( \lambda \delta \phi -2\Phi \right) \nonumber \\&\quad +2\sqrt{6}x_{2}\left( \Phi -\Phi ^{\prime }\right) +\frac{3}{\sqrt{6}} \alpha \left( \delta _{m}+2\Phi \right) , \end{aligned}$$

(85)

where a prime means derivative with respect to $\tau =\ln a.$ Recalling

$$\begin{aligned} \frac{{\mathcal {H}}^{\prime }}{{\mathcal {H}}}= -\frac{3}{2}(1+ w_\mathrm{eff}), \end{aligned}$$

(86)

where $w_\mathrm{eff}$ is given by (29), one could write,

$$\begin{aligned} \frac{{\mathcal {H}}^{\prime }}{{\mathcal {H}}}+4= \frac{1}{2} \left( 5-3 x_1 (x_1+x_2)+3 y_1^2-3 y_2\right) . \end{aligned}$$

(87)

Note that the above equations can be simplified by considering the Poisson equation, which in sub-horizon scales ($k/a\gg {\mathcal {H}}$) becomes [86]:

$$\begin{aligned} -\frac{k^2}{a^2 {\mathcal {H}}^2 }\Phi = \frac{3}{2} \left( \Omega _m \delta _m + (1+ 3 c_{d}^2)\Omega _{d} \delta _{d}\right) , \end{aligned}$$

(88)

where $c_{d}^2$ is the effective sound speed of DE perturbations (the corresponding quantity for matter is zero in the dust case),

$$\begin{aligned}&c_{d}^2 = w_{d} - \frac{w_{d}'}{3(1+w_{d})}, \end{aligned}$$

(89)

$$\begin{aligned}&\Omega _{d}= ( x_{1}^{2}+y_{1}^{2}) + (x_{1}x_{2} -y_{2}), \end{aligned}$$

(90)

$$\begin{aligned}&\Omega _{m}=1-( x_{1}^{2}+y_{1}^{2}) -(x_{1}x_{2} -y_{2}), \end{aligned}$$

(91)

$$\begin{aligned}&\Omega _{d} \delta _{d} \equiv - \Omega _{d}\delta T_{0}^{\left( d\right) 0}/ \rho _{d}= -\frac{\delta T_{0}^{\left( d\right) 0}}{3 {\mathcal {H}}^2}, \end{aligned}$$

(92)

where we have used $\delta T_{0}^{\left( d\right) 0}=-\rho _{d}\delta _{d}$ with $\delta T_{0}^{\left( d\right) 0}$ given by (77). That is,

$$\begin{aligned} \Omega _{d} \delta _{d}&= \frac{1}{3 {\mathcal {H}}^2 a^{2}}\left( {\dot{\phi }} \dot{\left( \delta \phi \right) } -{\dot{\phi }}^{2}\Phi \right) \nonumber \\&\quad +\frac{V_{,\phi }\delta \phi }{3 {\mathcal {H}}^2}+\frac{2\beta \hbar ^{2}}{3 {\mathcal {H}}^2 a^{2}}\left( \dot{\left( \delta \phi \right) } {\dot{\psi }}-2{\dot{\phi }}{\dot{\psi }}\Phi \right) \nonumber \\&\quad +\frac{2\beta \hbar ^{2}}{3 {\mathcal {H}}^2 a^{2} }{\dot{\phi }} \dot{\left( \delta \psi \right) } -2\beta \hbar ^{2} \frac{\psi \delta \psi }{3 {\mathcal {H}}^2}. \end{aligned}$$

(93)

Using $\dot{\left( \delta \phi \right) }= a {\mathcal {H}} {\delta \phi }', {\dot{\phi }}= a {\mathcal {H}} \phi ', \dot{\left( \delta \psi \right) }= a {\mathcal {H}} {\delta \psi }', {\dot{\psi }}= a {\mathcal {H}} \psi '$, where the dot denotes derivative with respect to the conformal time $\eta$, and $\phi '(\tau )=\sqrt{6} x_1, V(\phi )=3 y_1^2 {\mathcal {H}}^2, \psi '(\tau )=\frac{\sqrt{\frac{3}{2}} x_2}{2 \beta \hbar ^2},\psi = -\frac{3}{2} {\mathcal {H}}^2 \mu y_2$, we obtain

$$\begin{aligned} \Omega _{d} \delta _{d}&= \left( {\phi }' {\delta \phi }' -{\phi '}^{2}\Phi \right) +\frac{2\beta \hbar ^{2}}{3 }\left( {\delta \phi }' {\psi }'-2 {\phi }' {\psi }' \Phi \right) \nonumber \\&\quad +\frac{2\beta \hbar ^{2}}{3 } {\phi }' { \delta \psi }' +\frac{V_{,\phi }\delta \phi }{3 {\mathcal {H}}^2} -2\beta \hbar ^{2} \frac{\psi \delta \psi }{3 {\mathcal {H}}^2} \nonumber \\&= \left( {\phi }' {\delta \phi }' -{\phi '}^{2}\Phi \right) +\frac{2\beta \hbar ^{2}}{3 }\left( {\delta \phi }' {\psi }'-2 {\phi }' {\psi }' \Phi \right) \nonumber \\&\quad +\frac{2\beta \hbar ^{2}}{3 } {\phi }' { \delta \psi }' -\lambda \frac{V\delta \phi }{3 {\mathcal {H}}^2} -2\beta \hbar ^{2} \frac{\psi \delta \psi }{3 {\mathcal {H}}^2} \nonumber \\&= - 2 \Phi x_1^2 (4 \beta \hbar ^2 +3) +\ 2 \sqrt{\frac{2}{3}} \beta x_1 \hbar ^2 {\delta \psi }' \nonumber \\&\quad +\frac{(6 x_1+x_2) \delta \phi '}{\sqrt{6}}-\lambda y_1^2 {\delta \phi }+\beta \mu y_2 \hbar ^2 {\delta \psi }. \end{aligned}$$

(94)

The effective sound speed of the dark energy perturbations $c_{d}^2$ for $\kappa =0$ can be expressed as

$$\begin{aligned} c_{d}^2&= \frac{x_1 x_2^2 \left( 12 x_1-\sqrt{6} \mu y_2\right) +x_2 y_2 \left( 12 x_1+\sqrt{6} \mu y_2\right) -2 \sqrt{6} x_1 y_2 \left( 2 \lambda y_1^2+\mu y_2\right) }{12 x_1 x_2 (x_1 x_2+y_2)} \nonumber \\&\quad -\frac{\alpha y_2 \left( x_1 (x_1+x_2)+y_1^2-y_2-1\right) }{3 x_2 (x_1 x_2+y_2)}. \end{aligned}$$

(95)

In what follows, we are going to explore how scalar perturbations behave in the matter dominated solutions obtained from the background analysis.

5.1 Matter era

Consider now that we are near the stationary point $A_{3}$, then from the linear perturbation terms of the Einstein field equations and (83)–(85) it follows

$$\begin{aligned}&\delta _{m}^{\prime \prime }+\frac{7}{2}\delta _{m}^{\prime }+\frac{\alpha }{\sqrt{6}}\left( 2\delta \phi ^{\prime }-\delta _{\psi }\right) =0, \end{aligned}$$

(96)

$$\begin{aligned}&\delta \phi ^{\prime \prime }+\left( \frac{{\mathcal {H}}^{\prime }}{{\mathcal {H}} }+4\right) \delta \phi ^{\prime }-\delta _{\psi } =0\,, \end{aligned}$$

(97)

$$\begin{aligned}&2\beta \hbar ^{2}{\mathcal {H}}^{2}\left( \delta _{\psi }^{\prime \prime }+\left( \frac{{\mathcal {H}}^{\prime }}{{\mathcal {H}}}+4\right) \delta \psi ^{\prime } +\frac{33}{2}\delta _{\psi }\right) +\delta _{\psi } \nonumber \\&\quad +\frac{2\alpha }{3\sqrt{6} }\left( 6\delta _{m}+5\delta _{m}^{\prime }-6\zeta \delta \phi ^{\prime }\right) =0, \end{aligned}$$

(98)

in which $\delta \psi ={\mathcal {H}}^{2}\delta _{\psi }$, ${\mathcal {H}}^{2}\left( \tau \right) ={\mathcal {H}}_{0}^{2}e^{-3\tau }$ and for simplicity, we have assumed $k=0$. The latter system is a singular perturbation system that possesses a slow invariant manifold. Therefore, we continue our analysis by studying the evolution in the slow-fast manifolds [87, 88].

Now in the fast manifold we do the change of independent variable $\tau =2\beta \hbar ^{2}s,$ in Eq. (98) the dominated terms are $\left( \frac{\partial ^{2}\delta _{\psi }}{\partial s^{2}}-\frac{3}{2}\frac{\partial \delta _{\psi }}{\partial s}\right) \simeq 0,~$which provides $\delta _{\psi }\left( a\right) =\delta _{\psi }^{0}+\delta _{\psi }^{1} e^{\frac{3}{2}s}$, recall that $\ln a=2\beta \hbar ^{2}s$ and we assume that $\delta _{\psi }^{1}$ is large enough. Therefore, in the fast manifold from (96), (97) we find

$$\begin{aligned} \frac{\partial ^{2}\delta \phi }{\partial s^{2}}-\frac{3}{2}\frac{\partial \delta \phi }{\partial s}-\delta _{\psi }^{1}e^{\frac{3}{2}s}\simeq 0, \end{aligned}$$

(99)

and

$$\begin{aligned} \frac{\partial ^{2}\delta _{m}}{\partial s^{2}}-\frac{\alpha }{\sqrt{6}}\left( {\bar{\delta }}_{\psi }^{1}e^{\frac{3}{2}s}\right) =0. \end{aligned}$$

(100)

Fig. 8

The dependency of the exponents of the linear system (103), (104) with the coupling parameter $\alpha$ is displayed in the plots

Therefore,

$$\begin{aligned} \delta \phi \left( s\right)&=\frac{2}{9}e^{\frac{3}{2}s}\left( \left( 3s-2\right) \delta _{\psi }^{1}+\delta \phi _{1}\right) +\delta \phi _{2}, \end{aligned}$$

(101)

$$\begin{aligned} \delta _{m}\left( s\right)&=\frac{4}{9}\frac{\alpha }{\sqrt{6}}\left( {\bar{\delta }}_{\psi }^{1}e^{\frac{3}{2}s}\right) +s\delta _{m}^{0}+\delta _{m}^{0}, \end{aligned}$$

(102)

from where we observe that a growing mode exists in the fast manifold because of the higher-order derivative terms of the GUP. Hence, we see that the GUP component’s presence enhances scalar perturbations’ growth compared to the usual quintessence case.

On the other hand, in the slow manifold where $2\beta \hbar ^{2}\rightarrow 0$, we find the system

$$\begin{aligned}&\delta _{m}^{\prime \prime }+\frac{\alpha }{9}\left( 6\alpha \delta _{m} +5\alpha \delta _{m}^{\prime }-\sqrt{6}\left( \alpha ^{2}-3\right) \delta \phi ^{\prime }\right) +\frac{7}{2}\delta _{m} =0, \end{aligned}$$

(103)

$$\begin{aligned}&\delta \phi ^{\prime \prime }+\frac{\alpha }{9}\left( \sqrt{6}\left( 6\delta _{m}+5\delta _{m}^{\prime }\right) -6\alpha \delta \phi ^{\prime }\right) +\frac{5}{2}\delta \phi ^{\prime } =0.\nonumber \\ \end{aligned}$$

(104)

We replace $\delta _{m}\left( \tau \right) =\delta _{m}^{0}e^{A\tau }$ and $\delta \phi \left( \tau \right) =\delta \phi _{0}e^{A\tau }$ and the values of the exponent A depend of the coupling parameter $\alpha$, as they are presented in Fig. 8. The behavior of A concerning coupling parameter $\alpha$ exhibits two interesting results. First, the evolution of perturbations is independent of the sign of $\alpha$. Second, unlike the usual quintessence model, there is a decay of scalar perturbations (i.e., $A<0$) in the present model for the uncoupled scenario. However, scalar perturbations either decay or grow or describe an oscillatory solution for the coupling case. It is worth mentioning that in the case where scalar perturbations decay, the decay rate is slower compared to the uncoupled case. From Fig. 9, it is interesting to see that while for $\alpha =0,~1$ there is a decay of scalar perturbations, however, for $\alpha =3$ there is a growth of scalar perturbations. On assuming a negligible oscillating behavior of $\delta \phi$, the results reduce to that of [49]. However, in general, such an assumption cannot be applied as $\delta \phi$ may have non-oscillating terms since it is a singular perturbation system.

We have presented a qualitative analysis for the evolution of the perturbations in the slow-fast manifolds. From the above results, it is clear that the modification of the scalar field Action Integral by the GUP indeed modifies the evolution at the perturbation level.

Fig. 9

The evolution of matter density perturbation for $\alpha =0$ (solid curve), $\alpha =1$ (dotted-dashed curve) and $\alpha =3$ (dashed curve)

6 Conclusions

The Heisenberg Uncertainty Principle needs to be modified whenever we consider gravitational interaction. The GUP, which emerges from one such modification, has given rise to a host of exciting results. For instance, the quantum effects due to GUP modifies Hawking temperature and the Bekenstein entropy. Consequently, we get a modified Bekenstein system that critically affects Black hole evaporation. Similarly, in cosmology, the GUP is also very effective in explaining various issues related to the universe’s evolution. The cosmological background dynamics of the GUP modified quintessence scalar field analyzed in [34] provides engaging scenarios which are absent in the usual quintessence scalar field model. These exciting results motivate us to extend the work of [34] to the GUP modified interacting quintessence scenario at the background and perturbation levels.

More elaborately, here, we investigated a very generalized cosmological scenario of recent interest. We consider a quintessence scalar field cosmological model in a spatially flat FLRW background space in which the Lagrangian of the scalar field has been modified according to the quadratic GUP. Besides, we assumed that the scalar field is coupled to the dark matter component of the cosmological fluid. We considered three models already proposed in the literature for the interacting functions between the two fluid sources. For each cosmological model, we investigated the evolution of the global dynamics and studied the asymptotic behaviour.

First, we analyze the interacting DE model’s dynamics at the background level using the dynamical analysis tools. Indeed, we found that the models’ cosmological evolution is different in the presence of the minimum length or its absence. One of the present model’s interesting features compared to the corresponding interacting model of usual quintessence is that it exhibits matter-dominated and de Sitter solutions. While the former is crucial to explain the structure formation, the latter describes the late time acceleration at the background level. It is worth mentioning that these two solutions are due to the GUP modification and are independent of the choice of potential or coupling functions. Unlike the usual quintessence theory, the GUP modified interacting model cannot provide a solution that alleviates the coincidence problem.

To better understand background dynamics, we explored the model’s behavior at the linear perturbation level. More precisely, we focused on the scalar perturbations. We derived the dynamical equations for linear cosmological perturbation, where we found that the resulting linear equations form a singular perturbation system. Therefore, in contrast to the usual quintessence, we can explicitly write the perturbed equations’ solution in the fast and slow manifolds. In the fast manifold, the GUP components enhance the scalar perturbations’ growth due to the higher-order derivative terms of the GUP. However, in the slow manifold, the scalar perturbations either decay or grow or describe an oscillatory solution. Consequently, in the presence of the minimum length, the perturbation equations are also affected. It means that the present model differs from the usual quintessence scalar field theory.

In Ref. [89], the possibility to solve the cosmological constant problem with the introduction of the minimum length and specifically with the quadratic GUP has been investigated. However, the predicted value of the cosmological constant has been found to be $\Lambda \sim \frac{1}{\beta ^{2}}$ which cannot solve the smallness problem for the cosmological constant. A similar result was found recently in [90] by introducing a more general expression for the GUP where the authors found that $\Lambda \sim \frac{1}{\beta ^{4}}$. On the other hand, quintessence has been introduced as a theoretical model to introduce a time-varying effects of the DE fluid. However, for the simple exponential potential for the quintessence model the de Sitter universe is not recovered [2]. Thus, this is not the case where the minimum length is introduced to modify the scalar field Lagrangian. As we investigated in this study, the presence of the minimum length modifies the field equations in such a way so that the quintessence may reach the limit of the cosmological constant. This is an interesting mechanism in order the field equation to reach the de Sitter limit. The de Sitter expansion era can solve the “flatness”, “horizon” and monopole problems [91, 92] as it is explained by the support the cosmic “no-hair” conjecture [93, 94]. Furthermore, because the deformation parameter $\beta$ is found to be multiplied always with the derivatives of the scalar field, the new degrees of freedom can solve the smallness problem of the cosmological constant. For instance, considering the de Sitter point $A_{2}$ of the interaction model $Q_{A}$, it follows

$$\begin{aligned}&\sqrt{\frac{2}{3}}\lambda y_{1}^{2}=\beta \hbar ^{2}\frac{2\sqrt{2}{\dot{\psi }} }{\sqrt{3}H_{0}},\quad \sqrt{3}y_{1}=\frac{\sqrt{V}}{3H_{0}}, \\&y_{1}^{2} -1=\frac{\beta \hbar ^{2}\psi ^{2}}{3H_{0}^{2}}\quad \text {with }H_{0}=\sqrt{\frac{2}{3}\Lambda }. \end{aligned}$$

Hence, it is clear that $\Lambda$ does not depend explicitly on $\beta .$ However, if we assume $\lambda =0$ and $y_{1} =0$, which means that the initial scalar field model is the stiff fluid, from the latter expressions, it follows $\Lambda =-\frac{\beta \hbar ^{2}\psi ^{2}}{2}$. Hence, either from a simple quintessence model which does not provide any inflationary era, the modification of the Lagrangian function by the GUP provides an inflationary epoch driven by the higher-order derivatives of the scalar field.

This work contributes to the study of the existence of the minimum length in cosmological models. It would be interesting to investigate the constraint of the model concerning the upcoming precise cosmological observations in future works.