1 Introduction

We have entered the age of precision cosmology. Late time cosmic acceleration is one of the important discoveries of this era. Various observational artifacts like measurement of luminosity distance and red-shift of type Ia supernova [1], baryon acoustic oscillation [2], Cosmic microwave background radiation [3] power spectrum show that the present universe is going through a phase of accelerated expansion. It is widely assumed that an esoteric entity called ‘Dark energy,’ is solely responsible for this late time acceleration. The origin and constituents of this entity are still a big mystery. To tackle the question related to the cause of late time acceleration the \(\Lambda \)-CDM [4] model has been introduced, where \(\Lambda \) corresponds to the cosmological constant and CDM refers to the cold dark matter. Despite the countless success, \(\Lambda \)-CDM suffers from fine-tuning and cosmic coincidence problems. Another important downside of the \(\Lambda \)-CDM model is, from the observational perspective, the discrepancy between the present observed value of Hubble’s constant and the predicted value of Hubble’s constant from theory. This statistically significant tension has arisen due to the estimation of the Hubble constant, \(H_0\), from the Cosmic Microwave Background (CMB) and the direct local distance ladder measurements. In \(\Lambda \)-CDM scenario, analyzed data of Planck-CMB gives \(H_0 = 67.4 \pm 0.5 \mathrm km s^{-1} Mpc^{-1} \), which is in \(4.4\sigma \) tension with the local measurement \(H_0 =74.03 \pm 1.42 \mathrm km s^{-1} Mpc^{-1}\) [5,6,7]. In addition to that, several other observations also suggest a higher value of the Hubble constant. These fundamental discrepancies motivate us to study various non-minimally coupled field-fluid models of dark sector.

Several theoretical ventures of field-theoretic models have already been introduced to investigate the behavior of dark energy. Among them, canonical [8,9,10,11,12,13,14,15,16,17,18,19], non-canonical [20,21,22,23,24,25,26,27,28,29] scalar field formulations, f(R) modified gravity theory [30, 31] and scalar-tensor theory [32, 33] are well known. Scenarios with a canonical scalar field and a relativistic fluid, with algebraic and derivative coupling terms, have already been explored in several pieces of literature [34,35,36,37,38,39,40,41,42]. Advancing this scenario, we are introducing a non-canonical scalar field with a relativistic fluid which are interacting with each other non-minimally. In the present paper the non-minimal field-fluid coupling involves the single order derivative of the \(k-\)essence scalar field. In the present paper the field sector is governed by a non-canonical k-essence scalar field. The fluid sector is comprised by an ideal relativistic hydrodynamic fluid. Using a variational approach, we will analyze this coupled system in a flat homogeneous Friedmann–Lemaitre–Robertson–Walker universe.

In this paper, we have taken the arbitrary form of the interaction between field-fluid as \(f(n,s,\phi ,X) U^{\mu } \partial {}_{\mu } \phi \). The interaction term f contains a non-canonical kinetic term X, which contains the spacetime derivative of the scalar field. Derivative of field component \(\phi \) has been introduced which couples with fluid 4-velocity \( U^{\mu } \) as \( U^{\mu } \partial {}_{\mu } \phi \). In the previous works [43, 44] the present authors introduced the topic of non-minimal coupling between the k-essence field and an ideal, relativistic hydrodynamic fluid. In those works the field-fluid interaction term was independent of the first derivative coupling \(U^{\mu } \partial {}_{\mu } \phi \). In this paper we see that inclusion of the first derivative coupling modifies the previous results and requires a new formulation. Due to the presence of \(U^{\mu } \partial {}_{\mu } \phi \) factor, the field gradient directly couples with the fluid 4-velocity. Due to the presence of this coupling, the background field equations get modified. The derivative dependent coupling term \(f(n,s,\phi ) J^{\mu }\partial _{\mu } \phi \) was initially introduced in the context of quintessence field coupled with dark matter in Refs. [36, 45]. In our approach, we have extended the interaction functional form to \(f(n,s, \phi , X)J^{\mu }\partial _{\mu } \phi \). Our ignorance about the proper phenomenology of the dark sector forces us to rely more on the formal understanding of the system. The derivative coupling is motivated by this formal urge. Except for the algebraic coupling, this is another form of non-minimal coupling which the dark sector can have in principle. Which one out of these two is physically realized can only be ascertained after more phenomenological results are obtained. The derivative coupling has its own particular properties which are distinct from the algebraic coupling. In the present model, the non-minimal coupling is suppressed when the scalar field attains an extrema or becomes constant in any phase of cosmological evolution. These properties are not present in the algebraic non-minimal coupling case. The resulting coupling of fluid velocity \(U^{\mu }\) with the space-time derivative of field results in more complicated dynamics (when compared to the algebraic coupling) at the background and perturbation level. Any variations in \( U^{\mu }\) would directly propagate to the field \(\phi \) and vice-versa. These kinds of models may play a relevant role in the theories related to linear growth of large-scale structures. It would be of particular interest to test models based on this type of interaction against recently published data [46, 47]. This interaction term also gives rise to nontrivial non-linear effects on the dynamics of the early phase of the universe.

Using the above mentioned field-fluid coupling in our cosmological model we will present the dynamical stability analysis by advancing two different kinds of interacting models (corresponding to two different choices of f) in the presence of inverse square law k-essence scalar field potential \( V(\phi ) \propto 1/\phi ^2 \). In recent times the stability analysis approach in dark energy models has become useful [36, 48,49,50,51,52]. In the present work we will show the behavior of the non-minimally coupled model. It will be shown that the presence of non-minimal coupling can produce different phases of the universe. The coupled model can produce both phantom and non-phantom types of dark energy.

The topics discussed in the paper are presented in the following manner. The basic model equations and analysis of the system have been presented in Sect. 2. In Sect. 3 we show the energy–momentum conservation of the field-fluid system. The interacting field-fluid theory has been formulated in the context of a flat FLRW universe in Sect. 4. This coupled system’s dynamical stability has been thoroughly investigated using two kinds of interacting models in Sect. 5. A generalized conclusion has been drawn in Sect. 6.

2 The basic formulation

The introduction of the k-essence theory does not only describe the early universe scenario; rather, this theory is essentially sound to explain the late-time cosmic behavior. Several articles [21, 22] have already encompassed the dual usefulness of this non-canonical scalar field theory. Main ingredients of this theory are a scalar field \(\phi \) and a non-canonical kinetic term, \(X=-\frac{1}{2} g^{\mu \nu } \nabla _\mu \phi \nabla _\nu \phi \, \). In this section we construct the interacting model from the Lagrangian approach where scalar field \(\phi \) couples with a relativistic fluid non-minimally. Hence the total action of this coupled sector can be written as follows:

$$\begin{aligned} S= & {} \int _\Omega d^4x \bigg [\sqrt{-g}\frac{R}{2\kappa ^2}-\sqrt{-g}\rho (n,s) \nonumber \\{} & {} + J^\mu (\varphi _{,\mu } + s\theta _{,\mu } + \beta _A\alpha ^A_{,\mu })-\sqrt{-g}{{\mathcal {L}}}(\phi ,X)\nonumber \\{} & {} + f(n,s,\phi ,X)J^{\mu }\partial {}_{\mu }\phi \bigg ]. \end{aligned}$$
(1)

Here the First term of this total action is the Einstein–Hilbert action, where R is the Ricci scalar, \( \kappa ^2 = 8 \pi G \) and g specifies the determinant of metric \( g_{{\mu \nu }} \). We have considered the metric

$$\begin{aligned} ds^2=-dt^2 + a(t)^2 d\textbf{x}^2, \end{aligned}$$

with an isotropic and homogeneous (FLRW) background throughout the paper. Here a(t) is the scale factor in the FLRW spacetime. The next two terms of the grand action represent the relativistic fluid introduced by Brown [34]. The action consists of energy density term, \(\rho (n,s)\) which depends on particle number density, n and entropy density per particle, s. Other parameters \(\varphi , \theta , \alpha ^A, \beta _A\) are all Lagrangian multipliers. The current density \( J^{\mu } \) is related with the particle density n as,

$$\begin{aligned} J^\mu= & {} \sqrt{-g}\, n\,U^\mu , \quad |J|=\sqrt{-g_{\mu \nu }J^\mu J^\nu }, \quad \nonumber \\ n= & {} \dfrac{|J|}{\sqrt{-g}}, \quad U^\mu U_\mu =-1. \end{aligned}$$
(2)

where \(U^{\mu } = (1, \varvec{0})\) specifies the 4-velocity of a fluid parcel. The action of k-essence scalar field is specified by the fourth term in the action and the final term signifies the coupling between the scalar field and the relativistic fluid. The coupling term f is an arbitrary function of \( (n,s,\phi , X) \) and it is accompanied by the first derivative of the field \(\phi \), which is contracted with fluid current density \( J^{\mu } \). Because of this type of construction, the coupling term f has mass dimension one.

To find the equation of motion of the grand action we vary the Lagrangian in Eq. (1) with respect to the dynamical variables \(g_{{\mu \nu }}, J^{\mu },s, \theta ,\varphi , \beta _{A}, \alpha ^{A}, \phi \). The variation of the action with respect to \( g^{{\mu \nu }} \) gives the energy momentum tensor

$$\begin{aligned} T_{\mu \nu } \equiv \dfrac{-2}{\sqrt{-g}}\dfrac{\delta S}{\delta g^{\mu \nu }}. \end{aligned}$$
(3)

Therefore the energy–momentum tensor of the relativistic fluid, k-essence field and interaction becomes

$$\begin{aligned} T^{\mu \nu }_{(M)}= & {} n\dfrac{\partial {}\rho }{\partial {}n} \left[ U^{\mu } U^{\nu } + g^{{\mu \nu }} \right] -g^{{\mu \nu }} \rho , \end{aligned}$$
(4)
$$\begin{aligned} T^{\mu \nu }_{(\phi )}= & {} - {{\mathcal {L}}}_{,X}\,(\partial ^\mu \phi )(\partial ^\nu \phi )-g^{\mu \nu }\,{{\mathcal {L}}}, \end{aligned}$$
(5)
$$\begin{aligned} T^{{\mu \nu }}_{(\mathrm int)}&{=}&\left[ - n^2 \dfrac{\partial {}f}{\partial {}n}\left[ U^{\mu } U^{\nu } + g^{{\mu \nu }} \right] {+}n \dfrac{\partial {}f}{\partial {}X} (\partial {}^\mu \phi \partial {}^\nu \phi ) \right] U^{\alpha }\partial {}_{\alpha }\phi . \nonumber \\ \end{aligned}$$
(6)

Comparing these energy momentum tensor with the perfect fluid \( T^{{\mu \nu }} = \rho U^\mu U^{\nu }+P (g^{{\mu \nu }} + U^{\mu } U^{\nu }) \), we get energy density and pressure of the fluid from Eq. (4)

$$\begin{aligned} \rho _{M} = \rho , \quad P_{M} = \left( n\dfrac{\partial {}\rho }{\partial {}n} - \rho \right) . \end{aligned}$$
(7)

The k-essence energy density and pressure from Eq. (5) as

$$\begin{aligned} \rho _{\phi }= -{{\mathcal {L}}}_{,X}({\dot{\phi }}^2) + {\mathcal {L}}, \quad P_{\phi } = - {\mathcal {L}} \end{aligned}$$
(8)

The interaction energy density and pressure can be written from Eq. (6) as,

$$\begin{aligned} \rho _\textrm{int} = n\dfrac{\partial {}f}{\partial {}X}{\dot{\phi }}^3 , \quad P_\textrm{int} = \left[ - {n^2} \dfrac{\partial {}f}{\partial {}n} \right] {\dot{\phi }} . \end{aligned}$$
(9)

These three individual sectors’ energy–momentum tensors jointly produce the total energy–momentum tensor of the coupled sector, which can be presented as,

$$\begin{aligned} T^{\mu \nu }= T^{\mu \nu }_{(M)} + T^{\mu \nu }_{(\phi )} + T^{\mu \nu }_{(\textrm{int})}. \end{aligned}$$
(10)

Varying the action in Eq. (1) with respect to fluid parameters gives

$$\begin{aligned} J^{\mu }:&\quad \mu U_{\mu } + \left( \varphi _{,\mu }+ s\, \theta _{,\mu } + \beta _A \alpha _{,\mu }^A\right) \nonumber \\&\quad + \left[ \left( - n\dfrac{\partial {}f}{\partial {}n} \right) U^{\nu }U_{\mu } + f \delta _{\mu }^{\nu } \right] \partial {}_{\nu }\phi = 0 \,, \end{aligned}$$
(11)
$$\begin{aligned} s:&\qquad -\dfrac{\partial {}\rho }{\partial {}s} + \dfrac{\partial {}f}{\partial {}s} n U^{\mu } \partial {}_{\mu } \phi + nU^{\mu } \theta _{,\mu } = 0\, \,, \end{aligned}$$
(12)
$$\begin{aligned} \varphi :&\qquad J^\mu {}_{,\mu }=0 \,,\end{aligned}$$
(13)
$$\begin{aligned} \theta :&\qquad (sJ^\mu )_{,\mu }=0 \,,\end{aligned}$$
(14)
$$\begin{aligned} \beta _A:&\qquad J^\mu \alpha ^A_{,\mu }=0 \,,\end{aligned}$$
(15)
$$\begin{aligned} \alpha ^A:&\qquad (\beta _AJ^\mu )_{,\mu }=0 \,, \end{aligned}$$
(16)

where \(\mu \) signifies the chemical potential and defined as \(\mu = \dfrac{\partial {}\rho }{\partial {}n} = \dfrac{\rho + P}{n}\). Because of the interaction term the fluid dynamical equations (Eqs. (11) and  (12)) become generalized. Equation (13) shows that the particle number of the system is conserved, whereas, Eq. (14) indicates that there is no flow of entropy across the fluid elements. Hence the system is in an adiabatic process. Together they constrain the system and we can write

$$\begin{aligned} \nabla _{\mu }(n U^{\mu }) = 0 \quad \textrm{and} \quad \nabla _{\mu }(n s U^{\mu }) = 0, \end{aligned}$$
(17)

where \( \nabla _{\mu }\) is covariant derivative. Suppose in the action the interaction term does not exists i.e., \( f=0 \). Hence from the first law of thermodynamics one can define the temperature of the fluid as \( T = \dfrac{1}{n} \left( \dfrac{\partial {}\rho }{\partial {}s }\right) \) [34]. In that case Eq. (12) gives \( U^\mu \theta _{,\mu } = \dfrac{1}{n} \left( \dfrac{\partial {}\rho }{\partial {}s }\right) = T \), hence, \(\theta \) behaves as the potential for temperature. Using Eqs. (11),  (12) and  (15) one can show that the \( F = \mu - Ts = \varphi _{,\mu } U^\mu \), where F is the chemical free energy and \(\varphi \) signifies the potential of chemical free energy.Footnote 1 Because of the presence of interaction in the system i.e., \( f \ne 0 \), Eq. (12) gives

$$\begin{aligned} U^\mu \theta _{,\mu } = \dfrac{1}{n} \left( \dfrac{\partial {}\rho }{\partial {}s } - \dfrac{\partial {}f }{\partial {}s} n U^{\mu } \partial {}_{\mu }\phi \right) = T + T_\textrm{int}. \end{aligned}$$
(18)

Using Eqs. (11),  (12) and  (15) the potential for chemical free energy \( \varphi \) gets generalized

$$\begin{aligned} \varphi _{,\mu } U^\mu= & {} \mu - sT + U^{\alpha } \partial {}_{\alpha } \phi \left( s \dfrac{\partial {}f}{\partial {}s} - n \dfrac{\partial {}f}{\partial {}n} -f \right) \nonumber \\= & {} F + F_\textrm{int}. \end{aligned}$$
(19)

The above equations specify how the non-minimal field-fluid coupling affects the basic fluid equations.

Varying the action in Eq. (1) with respect to \(\phi \) produces the modified field equation. Due to the presence of non-minimal interaction term the field equation becomes:

$$\begin{aligned}{} & {} \bigg [- \dfrac{\partial {}{\mathcal {L}}}{\partial {}\phi } + \dfrac{\partial {}f}{\partial {}\phi } \, n U^{\alpha }\partial {}_{\alpha }\phi \nonumber \\{} & {} \quad - \partial {}_{\mu } \left( \dfrac{\partial {}{\mathcal {L}}}{\partial {}X} \partial {}^\mu \phi \right) + \partial {}_{\mu }\left( \dfrac{\partial {}f}{\partial {}X}(\partial {}^\mu \phi ) \right) (n U^{\alpha } \partial {}_{\alpha } \phi ) \nonumber \\{} & {} \quad + f_{,X} (\partial {}^{\mu }\phi ) \partial {}_{\mu }(nU^\alpha \partial {}_{\alpha }\phi ) - \partial {}_{\mu }(f) n U^{\mu } \bigg ] = 0. \end{aligned}$$
(20)

This equation can be expressed in terms of covariant derivative as:

$$\begin{aligned}{} & {} \bigg [ - \dfrac{\partial {}{\mathcal {L}}}{\partial {}\phi } - \nabla _{\mu } \left( \dfrac{\partial {}{\mathcal {L}}}{\partial {}X} \partial {}^\mu \phi \right) + \dfrac{\partial {}f}{\partial {}\phi } \, n U^{\alpha }\partial {}_{\alpha }\phi \nonumber \\{} & {} \quad + \nabla _{\mu }\left( \dfrac{\partial {}f}{\partial {}X}(\partial {}^\mu \phi ) \right) (n U^{\alpha } \partial {}_{\alpha }\phi ) \nonumber \\{} & {} \quad + f_{,X}(\partial {}^\mu \phi ) \nabla _{\mu }(nU^{\alpha }\partial {}_{\alpha }\phi ) \bigr \}- ( \nabla _{\mu }f) (nU^{\mu }) \bigg ] = 0 . \end{aligned}$$
(21)

Before we end this section we want to specify that in this paper heat exchange between various fluid parcels will be assumed to be zero. The whole fluid is assumed to be having the same temperature. The constraints in Eq. (17) specify that \({\dot{s}}=0\). Henceforth we will always assume that the entropy density of the fluid is constant.

3 Conservation of energy–momentum tensor in presence of non-minimal coupling

Due to the presence of the non-minimal coupling term, the individual energy–momentum tensors, in general, are not conserved, but the total energy–momentum tensor is conserved, hence \( \nabla _{\mu }T^{{\mu \nu }}=0 \). To prove the conservation of total energy–momentum tensor we regroup the various energy–momentum tensors and define two new energy–momentum tensors as: \( T^{'{\mu \nu }} = T_{(\phi )}^{{\mu \nu }} + n \dfrac{\partial {}f}{\partial {}X} (\partial {}^\mu \phi \partial {}^\nu \phi ) U^{\alpha }\partial {}_{\alpha }\phi \) and \({\tilde{T}}^{\mu \nu } = T_{(M)}^{{\mu \nu }} + T_{(\textrm{int})}^{{\mu \nu }} - n \dfrac{\partial {}f}{\partial {}X} (\partial {}^\mu \phi \partial {}^\nu \phi ) U^{\alpha }\partial {}_{\alpha }\phi \). Using these redefined tensors we can now write the Einstein field equation as:

$$\begin{aligned} G^{{\mu \nu }} =\kappa ^2 (T'^{\mu \nu }+{\tilde{T}}^{\mu \nu }). \end{aligned}$$
(22)

Now we can show that \( \nabla _{\mu }T^{'{\mu \nu }} =Q^{\nu } \) and \( \nabla _{\mu } \tilde{T}^{{\mu \nu }} = - Q^{\nu } \). To show that we will first calculate \(\nabla _{\mu }T^{'{\mu \nu }}\). Using Eq. (21) we get

$$\begin{aligned} \nabla _{\mu }T^{'{\mu \nu }}= & {} \bigg [f_{,n}( \nabla _{\mu }n) + f_{,X}( \nabla _{\mu }X) \bigg ] (nU^{\mu } )\partial {}^{\nu }\phi \nonumber \\{} & {} - f_{,X}( \nabla ^{\nu }X) (nU^{\alpha } \partial {}_{\alpha }\phi ) = Q^{\nu }. \end{aligned}$$
(23)

Now to show \( \nabla _{\mu } \tilde{T}^{{\mu \nu }} = - Q^{\nu } \), we express the covariant derivative into parallel and perpendicular components of the fluid flow as:

$$\begin{aligned} \nabla _{\mu } \tilde{T}^{{\mu \nu }} = h^{\nu }_{\lambda } \nabla _{\mu } \tilde{T}^{\mu \lambda } - U^\nu U_\lambda \nabla _{\mu } \tilde{T}^{\mu \lambda }, \end{aligned}$$
(24)

where \(h_{{\mu \nu }} = g_{{\mu \nu }} + U_{\mu }U_{\nu }\). Using \( \nabla _{\mu }(nU^{\mu }) = 0, \; {\dot{s}} = 0\), \( U_{\nu } \nabla _{\mu }U^{\nu } = 0 = U^{\nu } \nabla _{\mu }U_{\nu }\) and \( U^{\nu }U_{\nu } =-1 \), one can show that the last term of Eq. (24) can be written as

$$\begin{aligned} U_\nu \nabla _{\mu } \tilde{T}^{\,{\mu \nu }} = n^2\dfrac{\partial {}f}{\partial {}n} U^\alpha \partial {}_{\alpha }\phi \nabla _{\mu }U^\mu . \end{aligned}$$
(25)

One can also rewrite the first term of Eq. (24) in a different way. The manipulations involving the first term on the right hand side (of Eq. (26)) are explicitly shown in Appendix A. Using the result obtained in the appendix we can write

$$\begin{aligned} \begin{aligned} \nabla _{\lambda } \tilde{T}^{\beta \lambda }&= h^{\beta }_{\nu } \nabla _{\lambda } \tilde{T}^{\nu \lambda } - U^{\beta }U_{\nu } \nabla _{\lambda } \tilde{T}^{\nu \lambda }\\&= g^{\beta \mu }h_{{\mu \nu }} \nabla _{\lambda } \tilde{T}^{\nu \lambda } - U^\beta n^2\dfrac{\partial {}f}{\partial {}n} U^\alpha \partial {}_{\alpha }\phi \nabla _{\nu }U^\nu \\&= g^{\beta \mu } \big [ \nabla _\lambda (nf_{,n} U^\nu \partial {}_\nu \phi )n U^\lambda U_\mu \\&\quad + n U^\lambda ( \nabla _\lambda U_\mu )nf_n U^\nu \partial {}_\nu \phi - \bigl \{ f_{,n} \nabla _\lambda n\\&\quad + f_{,X} \nabla _{\lambda } X \bigr \} nU^\lambda \partial {}_\mu \phi +n \nabla _\mu (nf_{,n} U^\nu \partial {}_\nu \phi ) \\&\quad + \left\{ f_{,n} \nabla _\mu n + f_{,X} \nabla _\mu X \right\} nU^\lambda \partial {}_\lambda \phi \big ] \\&\quad - g^{\beta \mu } n^2f_{,n} U^\alpha \partial {}_{\alpha }\phi (U^\lambda \nabla _{\lambda }U_\mu ) \\&\quad - g^{\beta \mu } h^\lambda _\mu \nabla _\lambda \left( n^2f_n U^\alpha \partial {}_{\alpha }\phi \right) \\&\quad - U^\beta n^2\dfrac{\partial {}f}{\partial {}n} U^\alpha \partial {}_{\alpha }\phi \nabla _{\nu }U^\nu . \end{aligned} \nonumber \\ \end{aligned}$$
(26)

After some rearrangement one can show that

$$\begin{aligned} \nabla _{\lambda } \tilde{T}^{\beta \lambda }= & {} Q^{\beta }=- \left( f_{,n} ( \nabla _\lambda n ) \right. \left. + f_{,X} ( \nabla _{\lambda } X) \right) nU^\lambda \partial {}^\beta \phi \nonumber \\{} & {} + f_{,X} \nabla ^\beta X (nU^\lambda )\partial {}_\lambda \phi .\nonumber \\ \end{aligned}$$
(27)

Comparing Eq. (23) and Eq. (27) it can be seen that the total energy–momentum tensor is conserved as:

$$\begin{aligned} \nabla _\mu (T'^{\mu \nu }+ \tilde{T}^{\mu \nu }) = 0. \end{aligned}$$
(28)

4 The dynamical equations

Using the above general framework of coupled field-fluid scenarios we investigate this system in the context of a spatially flat FLRW universe. The Friedmann equations get modified due to the presence of an interacting scenario. The Friedmann equation and other equations specifying the fluid, like the equations specifying the conservation of particle number density and entropy per particle, can be written as:

$$\begin{aligned} 3H^2= & {} \kappa ^2 \left( \rho +\rho _{\phi }+\rho _\textrm{int}\right) = \kappa ^2 \rho _\textrm{tot}, \end{aligned}$$
(29)
$$\begin{aligned} 2{\dot{H}}+3H^2= & {} -\kappa ^2 \left( P+P_{\phi }+P_\textrm{int}\right) = -\kappa ^2 P_\textrm{tot}. \end{aligned}$$
(30)
$$\begin{aligned} {\dot{n}}+3Hn= & {} 0,\quad \textrm{and} \quad {\dot{s}}=0. \end{aligned}$$
(31)

Here \(H={\dot{a}}/a\) is the Hubble parameter. We have used the form of k-essence Lagrangian, \({{\mathcal {L}}}=-V(\phi )F(X)\) where F(X) is a non-canonical kinetic function of X and potential \(V(\phi )\) depends on the scalar field \(\phi \). Hence the energy density and pressure of the k-essence sector can be written as

$$\begin{aligned} \rho _{\phi }= & {} V(\phi )(2XF_{,X}-F),\quad P_{\phi }=V(\phi )F(X). \end{aligned}$$
(32)

The equation of state (EOS), \(\omega _\phi \), and adiabatic sound speed \(c_{s(\phi )}^2\) in the k-essence sector are:

$$\begin{aligned} \omega _\phi = \frac{P_{\phi }}{\rho _{\phi }},\quad c_{s(\phi )}^2 = \frac{dP_\phi /dX}{d\rho _\phi /dX}. \end{aligned}$$
(33)

In the absence of an interacting term the above EOS determines the nature of the system at late times. The sound speed limit put additional constrained on the system as \( 0 \le c_{s(\phi )}^2 \le 1 \). For negative sound speed the system becomes unstable.

The modified scalar field equation, Eq. (21), can be expressed as:

$$\begin{aligned}{} & {} -{\mathcal {L}}_{,\phi } + {\mathcal {L}}_{,\phi X}{\dot{\phi }}^2 + {\mathcal {L}}_{,XX}{\dot{\phi }}^2 \ddot{\phi } + 3H{\mathcal {L}}_{,X} {\dot{\phi }} \nonumber \\{} & {} \quad + {\mathcal {L}}_{,X}\ddot{\phi } + \bigg [3nH f_{,nX}{\dot{\phi }} - f_{,\phi X} {\dot{\phi }}^2\nonumber \\{} & {} \quad - f_{,XX} {\dot{\phi }}^2 \ddot{\phi } - f_{,X} \ddot{\phi } \bigg ]\, n {\dot{\phi }} - n f_{,X} {\dot{\phi }} \ddot{\phi } \nonumber \\{} & {} \quad + 3n^2 H f_{,n} - n f_{,X} {\dot{X}} = 0. \end{aligned}$$
(34)

where subscript XX or \(\phi X\) implies a second-order partial derivative with respect to X or \(\phi \) and X. The \( Q^{\nu } \) in this background metric will become \( Q^{0} = n^2f_{,n} 3\,H {\dot{\phi }} \). Therefore the field-fluid exchanges energy density via this coupling. In the next section, we will set up the dynamical approach to study the coupling scenario.

5 Dynamical stability of non-minimally coupled system

In this section we will construct the dimensionless dynamical variables to study the dynamics of the interacting system. The basic dynamical variables are \( x,y,z,\sigma \) which are defined as

$$\begin{aligned} x = {\dot{\phi }}, \quad y = \frac{\kappa \sqrt{V(\phi )}}{\sqrt{3}H}, \quad \sigma = \frac{\kappa \sqrt{\rho }}{\sqrt{3}H}, \quad z = \frac{\kappa ^2 nf}{3H^2} \end{aligned}$$
(35)

The other variables A, B, C, D, E can be expressed in terms of the basic variables as:

$$\begin{aligned} \quad B= & {} \frac{nf_{,\phi }\kappa ^2}{H^3}, \quad E= \dfrac{n\kappa ^2}{H^3} \frac{\partial ^2 f}{\partial \phi \partial X},\quad D = \dfrac{ n\kappa ^2 }{3H^2}f_{,X}, \quad \nonumber \\ A= & {} -\dfrac{n^2\kappa ^2}{H^2} \frac{\partial f}{\partial n}, \quad C= \dfrac{\kappa ^2 P_\textrm{int}}{3H^2} \end{aligned}$$
(36)

Due to the presence of a constraint equation only x, y, z are used to study the dynamical system. The constrained equation of the system can be written by using Eq. (29) as:

$$\begin{aligned} \sigma ^2 = 1- y^2\left( 3x^4/4-x^2/2\right) - x^3 D. \end{aligned}$$
(37)

where \(\sigma \) denotes the fluid variable and it satisfies the inequality \( 0 \le \sigma ^2 \le 1 \). The other Friedmann equation can be expressed as

$$\begin{aligned} \dfrac{{\dot{H}}}{H^2} = -\frac{3}{2} [\omega \sigma ^2 + y^2 F+C+1]. \end{aligned}$$
(38)

The dynamical equations can be expressed compactly if we represent the dimensionless energy density of the scalar field as \(\Omega _{\phi }\) and interaction energy density as \( \Omega _\textrm{int} \)

$$\begin{aligned} \Omega _{\phi } = \frac{\kappa ^2 \rho _{\phi }}{3H^2} = y^2 (x^2 F_{,X} - F), \quad \Omega _\textrm{int} = \frac{\kappa ^2 \rho _\textrm{int}}{3H^2} = D x^3.\nonumber \\ \end{aligned}$$
(39)

Using the variables defined above the autonomous equations of the system are:

$$\begin{aligned} x^{\prime }= \frac{{\dot{x}}}{H}= & {} \dfrac{\lambda F- \lambda F_{,X}x^2-3 x y^2F_{,X}-x^2A_{,X}-\frac{E}{3} x^3 - A}{[y^2F_{,XX}+ D_{,X} x]x^2 + [y^2 F_{,X} + 3 D x]}\nonumber \\ \end{aligned}$$
(40a)
$$\begin{aligned} y^{\prime } = \frac{{\dot{y}}}{H}= & {} \dfrac{\lambda x }{2y} + \dfrac{3}{2}y \left( \omega \sigma ^2 + y^2 F+C+1 \right) \nonumber \\ \end{aligned}$$
(40b)
$$\begin{aligned} z^{\prime } = \frac{{\dot{z}}}{H}= & {} 3z\left( \omega \sigma ^2 + y^2 F+C \right) + A + \frac{Bx}{3} + D x x^{\prime }\nonumber \\ \end{aligned}$$
(40c)

The system possesses 3-D autonomous equations for any general interaction. Here prime signifies the derivative of the dynamical variables x, y, z with respect to Hdt. In these equations \(\lambda \) is the dimensionless parameter which signifies the derivative of k-essence sector potential \(\lambda = \dfrac{\kappa ^2 V_{,\phi }}{3\,H^3} \). Corresponds to these autonomous equations the critical points \( (x_{0}, y_{0}, z_{0}) \) can be found by solving \(x^\prime =0\), \(y^\prime =0\) and \(z^\prime =0\). To understand the behavior of the system near these critical points we linearize the autonomous equation around these critical points using the Taylor series expansion up to the first order and hence find the Jacobian matrix. Based on the eigenvalues of the matrix the nature of fixed points can be established. The \( 3 \times 3 \) matrix posses three eigenvalues and if all the eigenvalues have negative (positive) real parts then the fixed point is asymptotically stable (unstable). If any of these eigenvalues which have a positive real part and the rest have a negative real part; the fixed points become saddle. If any of the eigenvalues turns out to be zero, the linear stability theory fails, and one needs to employ the centre manifold theory to further understand the nature of those critical points.

To further study the autonomous equations we need to introduce the form of interaction and potential. Hence we will use a power-law type potential of k-essence which has been studied extensively [22, 25,26,27] as:

$$\begin{aligned} V(\phi ) = \dfrac{\delta }{\kappa ^2\phi ^2}. \end{aligned}$$
(41)
Table 1 Models parameters of the interacting k-essence field-fluid
Full size table

Here \(\delta \) is a dimensionless parameter. This type of inverse square law scalar field-dependent potential is widely used to probe the late time accelerating universe in cosmological models employing k-essence theory. We have also introduced two different interaction models in Table 1. The form of f as tabulated are some of the simplest choices one can make in the present context. Based on these forms of interaction we will establish the dynamics in the next section. To proceed further we have chosen the form of the kinetic part of the k-essence field as:

$$\begin{aligned} F(X) = -X+ X^2. \end{aligned}$$
(42)

This form of F is also widely used in the literature.

5.1 Study of model I

To study the field-fluid coupling of the system we have chosen the form of interaction \( f= \dfrac{\rho }{n} (\phi /\kappa )^m X^{\eta } \), where \( \rho \) signifies the matter-energy density, n is the fluid number density and \(\phi \), X depend on the field and its time derivatives. Here m and \(\eta \) are the model parameters and can take any value based on the system dynamics. The total energy density and pressure for this system can be expressed as:

$$\begin{aligned} \rho _\textrm{tot}= & {} \rho + \rho _{\phi } +\rho _\textrm{int} =\rho + V(\phi ) (2XF_{,X}-F) +n\frac{\partial {}f}{\partial {}X}{\dot{\phi }}^3, \nonumber \\ \end{aligned}$$
(43)
$$\begin{aligned} P_\textrm{tot}= & {} P + P_{\phi } + P_\textrm{int} = P + V(\phi )F -\left( n^2 \frac{\partial {}f}{\partial {}n} \right) {\dot{\phi }}.\nonumber \\ \end{aligned}$$
(44)

Constraint equation for this model can be written as:

$$\begin{aligned} \sigma ^2 = 1- y^2\left( \frac{3}{4} x^4-\frac{1}{2} x^2\right) - 2x\eta z. \end{aligned}$$
(45)

Using the above expressions we can define the grand EOS and adiabatic sound speed as:

$$\begin{aligned} \omega _\textrm{tot}= & {} \frac{P_\textrm{tot}}{\rho _\textrm{tot}} = \omega \sigma ^2 + y^2\left( \frac{1}{4} x^4-\frac{1}{2} x^2\right) +C , \end{aligned}$$
(46)
$$\begin{aligned} c_s^2= & {} \dfrac{dP_\textrm{tot}/dX}{d\rho _\textrm{tot}/dX} = \dfrac{y^2 F_{,X} - 2 \omega z(\eta + 1/2)/x}{y^2(F_{,X} + 2XF_{,XX}) + \frac{4\eta (\eta -1)z }{x} + \frac{6 \eta z}{x}}. \nonumber \\ \end{aligned}$$
(47)

5.2 Dynamics of the model

We will study this model in the presence of both radiation (\( \omega = 1/3 \)) and matter (\( \omega = 0 \)). The 3-D autonomous system is complicated and contains three model-dependent parameters \(\delta , m\) and \(\eta \). Hence, finding the critical points in model-dependent parameters is a rather complicated task. To resolve this issue, we will use the following technique. We have checked that our system contains eight critical points, when we deal with dust and six critical points, when we deal with radiation. Some of these critical points does not explicitly depend upon m and \(\eta \). The others may depend upon all the three parameters. We figure out the critical points which are independent of specific values of m and \(\eta \) in a particular way and later find out the other critical points in a more conventional approach.

To find out the m and \(\eta \) independent critical points we have expressed \(\omega _\textrm{tot}\), at a critical point, in terms of \( \gamma \), which is an arbitrary constant. We have replaced (the functional form of the kinetic term of k-essence scalar field) F, at the critical point, in terms of other variables present in the expression of Eq. (49). We have then

$$\begin{aligned} \omega _\textrm{tot}= & {} -1 + \gamma , \end{aligned}$$
(48)
$$\begin{aligned} F= & {} \dfrac{-1 + \gamma - \omega \sigma ^2 - C}{y^2}. \end{aligned}$$
(49)

Using these relations, at the critical points, we find out some of the critical points which do not explicitly depend on the values of m and \(\eta \).

5.2.1 Dynamics in presence of dust

By replacing F in Eqs. (40a)–(40c) by the above value we get one set of the critical point \( P_{\mp } \) of the system which is solely dependent on \( \gamma \) and \(\delta \). The general form of the (m and \(\eta \) independent) critical points in the matter-dominated background have been tabulated in the top row of Table 2. Using this set of critical point \( P_{\mp } \) in the expression of the EOS, in Eq. (46), we find the relation between \(\gamma \) and \(\delta \) shown in Fig. 1. There exist two solutions of \(\gamma \) as \(\gamma _{1}\) and \(\gamma _{2}\) which characterize the non-phantom and phantom solutions, respectively, for one set of \(\delta \) values. Inserting the value of \( \gamma _{1} \) and \(\gamma _{2}\) in the critical point expression \( P_{\mp } \) produces the four actual \(\delta \) dependent critical points listed in Table  2, which are labeled as \( P_{1\mp } \) and \( P_{2\mp } \).

Table 2 Critical points in the matter dominated background for general m and \(\eta \) for model I
Full size table
Fig. 1
figure 1

Variation of \(\gamma \) against \( \delta \) for model I

Full size image

From Fig. 1 we see that for \(0<\delta <1\), grand EOS parameter shows a non-accelerating behavior corresponding to critical point \( P_{2\mp } \), while \( P_{1\mp } \) shows phantom behavior. As \(\delta \) increases, value of grand EOS parameter corresponding to \( P_{2\mp } \) is going towards an accelerating EOS (\(\gamma _{1} \rightarrow 0\)) parameter while \( P_{1\mp } \) also converges towards \( (\gamma _{2} \rightarrow 0)\) and produces less phantom behavior. Hence, by choosing a larger value of \(\delta \) one can find both accelerating and phantom solutions. Using this technique, we have found four critical points which are independent of m and \(\eta \) for a large value of \(\delta \). As the system has eight critical points we have to find four more critical points. It turns out that the other four critical points, called \( Q_{\mp } \) and \( R_{\mp } \), have to be found out from the equations \(x^\prime =0\), \(y^\prime =0\) and \(z^\prime =0\) by. This critical points may depend on m and \(\eta \). We have evaluated all the other critical points by choosing the benchmark value of \(\delta = 20\). For \(\delta = 20\) the other sets of critical points \( Q_{\mp } \) and \( R_{\mp } \) are shown in Table 3, in which only \( Q_{\mp } \) are both m and \(\eta \) dependent. The other critical points \( R_{\mp } \) are evaluated numerically. For \( Q_{\mp } \), we evaluated the grand EOS parameter, which turns out to be only m dependent. Consequently this critical point can be set to any phase of the universe based on the value of m. We have chosen to set it in the radiation phase and m turns out to be \( m = 1/2 \). To further constraint the model parameter \( \eta \) corresponding to \( m=1/2 \) we have used the sound limit, which puts a bound on \(\eta \), and get \(\eta \le {31}/{20}\). We have listed the all critical points corresponding to the case where \(\omega =0\) in Table 3 for \( \delta =20, m= 1/2 \) and \(\eta = 1\) case.

We have verified that none of the critical points corresponding to \( w=0 \) exists at infinity. To understand the nature of the system, we have portrayed the 3-D phase space and evolution of system parameters with respect to the logarithmic function of the FLRW scale factor. In presence of dust two sets of critical points \(P_{1\mp }\) and \(P_{2\mp }\) are stable attractors. Near one point the universe shows accelerating phantom behavior with negative sound speed whereas near the other set we have accelerating behavior as \( -1< \omega _\textrm{tot}<-1/3 \) and the sound speed is positive around the second point. The dynamics of the critical points can be understood from the 3-D phase space plot in Fig. 2a in which trajectories corresponding to blue color show that they are attracted towards a critical point while red-colored trajectories signify that they are repelled from the critical points.

In the present case the phase space is not compact as \(-\infty< x < \infty \), \(-\infty< y < \infty \) and \(-\infty< z < \infty \). To represent the phase space in a bounded region we have used a transformation of phase space variables. These transformations connect the conventional phase space variables (x, y, z) to the new set of variables (X, Y, Z) where the new variables are defined as:

$$\begin{aligned} X=\tan ^{-1}x,\,\,\,\,Y=\tan ^{-1}y,\,\,\,\,Z=\tan ^{-1}z. \end{aligned}$$
(50)

These new variables map the unbounded phase space to a bounded region as

$$\begin{aligned} -\frac{\pi }{2}< X< \frac{\pi }{2},\,\,\,\,-\frac{\pi }{2}< Y< \frac{\pi }{2},\,\,\,\,-\frac{\pi }{2}< Z < \frac{\pi }{2}. \end{aligned}$$

We use these new variables to redefine the phase space in both the cases, in presence of dust or radiation, in Model I.

Table 3 The critical points and their nature in presence of dust for \(\delta = 20, \) for general m and \(\eta \) (the top two rows), and for \( m= 1/2, \eta = 1\) (for the last two rows)
Full size table

We see that the trajectories nearby \( Q_{\mp } \) evolve towards \(P_{2\mp }\), showing that \( Q_{\mp } \) signifies saddle type critical points. We can identify the points \( Q_{\mp } \) correspond to the radiation phase of the universe. Hence all the nearby trajectories are globally attracted towards \( P_{2\mp } \). The other trajectories nearby the points \( R_{\mp } \), which signify the matter phase, are repelled from that fixed points and eventually are attracted towards \( P_{2 \mp } \). Few trajectories are attracted towards \( P_{1\mp } \), making the point a stable attractor. Hence we see that for the dust solutions the field-fluid scenario can successfully depict all four phases of the universe.

Fig. 2
figure 2

3D phase space and evolution plot in presence of dust for model I

Full size image

We have evaluated the evolutionary dynamics of system parameters such as grand EOS (\(\omega _\textrm{tot}\)), k-essence energy density \( \Omega _{\phi } \), fluid density \( \sigma ^2 \), sound speed \( c_{s}^2 \) and interaction density \( \Omega _\textrm{int} \) and have shown their development in Fig. 2b. From the figure, we analyze that as the system enters the very early epoch the coupling between field-fluid becomes very strong and k-essence density and interaction density start with a high magnitude. Since k-essence energy density always has to be positive \( \Omega _{\phi } \ge 0 \), and system follows the constrained relation Eq. (37), hence \( \Omega _\textrm{int} \) becomes negative. Due to the high energy density of k-essence field, the grand EOS of the system tends to be stiff. The sound speed is also greater high initially. As the system evolves, the total EOS of the system becomes 1/3, and the k-essence energy density dominates over the fluid energy density \( \sigma ^2 \) and \( \Omega _\textrm{int} \). This phase is representative of the radiation dominated phase. The fluid energy density, \( \sigma ^2 \), then starts increasing and \( \Omega _{\phi } \) decreases and becomes subdominant in the matter phase; at this phase, the EOS becomes zero. The k-essence energy density is not zero during the matter phase and hence keeps increasing and becomes dominant at the late-time phase. With this phase transition the grand EOS also decreases towards a negative value and gets saturated at \( \approx -0.6 \). The sound speed also decreases and becomes \( \approx 10^{-3} \). Hence the system transforms from stiff matter to the late-time cosmic acceleration when the hydrodynamic fluid represents dust.

Table 4 Critical points, in presence of radiation, for general m and \( \eta \)
Full size table
Table 5 Model I critical points, in presence of radiation, for \(\delta = 20, \) and general m and \(\eta \), and stability for \( m= -2/3, \eta = 1\)
Full size table

5.2.2 Dynamics in presence of radiation

In this subsection we will explore the interaction between k-essence scalar field and radiation fluid. Using the same technique, as discussed above after Eq. (49), we find the critical points \( P_{\mp } \) listed in Table 4. We get the same \( \gamma -\delta \) relation shown in Fig. 1 for these critical points. Hence choosing \( \delta = 20 \), we get three sets of critical points, which are tabulated in Table 5 for any value of m and \(\eta \). Out of these three sets of critical points, \( P_{1\mp } \) and \( P_{2\mp } \) are similar to the previous case. Only critical points \( Q_{\mp } \) are both m and \(\eta \) dependent. At these points the grand EOS parameter depends on m and \(\eta \). For \( m=-2/3 \) the points signify the matter phase. The sound speed limit puts further constraint on \(\eta \) as \( \eta \le -0.19 \) or \( -0.03< \eta < 1.7 \). Using all these mentioned constraints we have tabulated the critical points and their corresponding properties in Table 5 for \( \delta = 20, m= -2/3\) and \(\eta = 1\).

From the phase space plot in Fig. 3a (where again we have used the variables X, Y, Z as defined in Eq. (50)) we notice that both \(P_{1\mp }\) and \(P_{2\mp }\) are stable attractors. The trajectories nearby \( Q_{\mp } \) evolve towards the \( P_{2\mp } \), and makes the \( P_{2\mp } \) into a global attractor near which the universe shows accelerating expansion. This figure depicts that the universe may start at some point and gets attracted towards the matter phase \( Q_{\mp } \) and then eventually end up towards the accelerating phase.

In Fig. 3b the dynamics of the system becomes highly non-linear. We have plotted the regions where the grand EOS parameter of the system enters the radiation phase. During this phase, fluid density (\( \sigma ^2 \)) dominates and interaction energy density remains subdominant. As the k-essence energy density evolves, the fluid density decreases, and the value of grand EOS parameters becomes zero. The change in phase is very steep, and at this phase, k-essence energy density becomes greater than fluid density. This signify that the field itself behaves like the dark-matter candidate. As k-essence dominates in the late-time phase the system enters the accelerating phase. During the whole evolution period sound speed lies between 0 to 1.

Fig. 3
figure 3

3D phase space and evolution plot in presence of radiation for model I

Full size image

5.3 Study of model II

The chosen form of the interaction that corresponds to this model is \( f= \dfrac{F}{n \kappa ^2 \phi ^2} \), where \( F= X^2-X\) is the usual kinetic term of the k-essence field. In this model the interaction term does not contain the energy density of the fluid \( \rho \); rather it only depends on one of the fluid parameters, the number density n. Because of this the dimensionality of the phase space reduces from three to two. Therefore the interacting system can be studied using only two dynamical variables x and y. The model parameters have been evaluated in Table 1. The constrain equation can be written as:

$$\begin{aligned} \begin{aligned} \sigma ^2&= 1 - y^2 x^2 \left[ (x^2 -1) \left( \dfrac{3}{4} - \dfrac{x}{\delta } \right) + \dfrac{1}{4} \right] . \end{aligned} \end{aligned}$$
(51)

The equation of the state and generalized adiabatic sound speed of the system can be defined as:

$$\begin{aligned} \omega _\textrm{tot}= & {} \omega \sigma ^2 + y^2 F\left( 1 + \dfrac{1}{\delta } x \right) , \quad \nonumber \\ c_{s}^2= & {} \dfrac{F_{,X} + x F_{,X}/\delta + F/(x \delta )}{(F_{,X} + 2 X F_{,XX}) + F_{,XX} x^3/\delta + 3 x F_{,X} /\delta } . \end{aligned}$$
(52)

Using these quantities we study the dynamical development of the system.

5.4 Dynamics of model II

In the present case the equation of state and the sound speed are not invariant under \( x \rightarrow -x \) transformation. Sound speed is solely a function of dynamical variable x. The critical points of the system can be found by solving the autonomous equations \( x' = 0, y' = 0\) in Eqs. (40a) and  (40b) for some values of \(\delta \). The system turns out to be complicated and hence it is difficult to calculate analytically the critical points in terms of the model parameters. Hence we have adopted the numerical technique to find the critical points. We choose \(\delta \) in such a way that at a late time it produces enough acceleration. Corresponding to the specific choice of \(\delta \) we will find the dynamics of the system. We have found no critical points at infinity in the present case.

5.4.1 Dynamics in the presence of dust

We have found six critical points in the present case, when the fluid has \(\omega =0\). The critical points are shown in Table 6. Out of these six critical points, \( P_{1}, \) and \( P_{2} \) are saddle points and the corresponding grand EOS parameters near these points are approximately zero; hence these critical points signify the matter-dominated phase of the universe. At these phase points the sound speed is \(\approx 0.2\). The other four critical points of this scenario are stable. To understand the nature of the critical points we have depicted the 2-D phase space, shown in Fig. 4a. The phase space is divided into red, green, and yellow regions, which signify the system’s physical characteristics. The red region signifies the phantom phase, the yellow region corresponds to the accelerating phase and the blue region shows the region where the sound speed limit is obeyed. The rest of the region shows a non-accelerating phase. The nearby trajectories of points \( P_{1} \) and \( P_{2} \) are repelled from these and attracted towards the points \( P_{6} \) and \( P_{5} \). A red phase space trajectory has been plotted by solving the autonomous system, which connects \( P_{1} \) to \( P_{6} \) and \( P_{2} \) to \( P_{5} \) by heteroclinic orbits. The evolution of the trajectory from \( P_{1} \) and \( P_{2} \) signifies that the universe evolves from the matter dominated phase and enters into the accelerating phase. It turns out that \( P_{5} \) and \( P_{6} \) are globally attractors near which the sound speed remains small. The other fixed points \( P_{3} \) and \( P_{4} \) correspond to the phantom phase with a negative sound speed. The nearby vector points are attracted towards these points; hence the phantom phase of the system is stable.

Table 6 The critical points for model II, in the presence of dust, for \(\delta = 45.11\)
Full size table
Fig. 4
figure 4

In phase space the blue region shows the sound speed limit \( 0<c_{s}^2<1 \), red region shows the phantom \( \omega _\textrm{tot}< -1 \) and yellow region signifies \( - 1\le \omega _\textrm{tot} \le -1/3 \)

Full size image

The dynamic evolution of the system is shown in Fig. 4b. It is seen that at an early epoch the grand EOS is \( \approx 0.25 \). At this phase both k-essence energy density and fluid energy density \( \sigma ^2 \) are subdominant while \(\Omega _\textrm{int}\) dominates. The early phase mimics a radiation dominated phase. As the universe evolves both k-essence energy density and \(\sigma ^2\) increase and grand EOS decreases. When \(\sigma ^2\) increases the grand EOS becomes zero. At this point both interaction density and \(\Omega _{\phi }\) are subdominant. It is seen that as \(\Omega _{\phi }\) starts increasing the grand EOS decreases and saturates to \( \approx -0.7 \). At the late time both fluid density and interaction density become subdominant, and the sound speed of the universe becomes close to zero. In presence of dust, the coupling term explains the late-time scenario of the universe and it also gives rise to some interesting behavior in the early phase.

5.4.2 Dynamics in presence of radiation

In the present case we have chosen the model parameter \( \delta = 20.63 \). For this benchmark point we found a total of five fixed points out of which one of the fixed points \( P_{1} \) corresponds to radiation dominated phase. We have listed all the fixed points and their nature in Table 7. The 2-D phase space is displayed in Fig. 5a. We see that \( P_{2} \) and \( P_{3} \) are fixed points near which we get accelerated expansion with positive sound speed and nearby trajectories are attracted towards them. The fixed point corresponding to the radiation phase is \( P_{1} \), some of the trajectories are attracted towards it and finally they get repelled. The evolution around this point is shown by a red trajectory which connects with \( P_{2} \) and completes a heteroclinic orbit. This shows that the radiation phase evolves into an accelerating phase. There are some nearby vector points of \( P_{1} \), which shows that they are repelled from this point and go towards \( P_{2} \) from the blue-white region, which may explain why the universe goes from radiation to matter phase and then finally enters into the accelerating phase. This point will become clear when we analyze the system dynamics. There exists no non-accelerating fixed point corresponding to \( P_{3} \), which signifies the non-symmetric nature of the autonomous system. Here \( P_{2} \) becomes the global attractor. Like the previous case, the phantom fixed points are stable with nearly zero sound speed. The dynamics as shown in Fig. 5b has similar behavior compared to our previous case (dynamics in presence of dust) in late-time phase. However, in the early epoch, the fluid density dominated with EOS 1/3. The interaction density increases for some time but then it decreases gradually.

Table 7 The critical points for model II, in the presence of radiation, for \(\delta = 20.63\)
Full size table
Fig. 5
figure 5

In phase space the blue region shows the sound speed limit \( 0<c_{s}^2<1 \), red region shows the phantom \( \omega _\textrm{tot}< -1 \) and yellow region signifies \( - 1\le \omega _\textrm{tot} \le -1/3 \)

Full size image

In the present case the k-essence density, \( \Omega _{\phi } \), is monotonically increasing and dominates at the late time. As the universe evolves from the early phase, the grand EOS parameter starts decreasing and crosses the 0 line, specifying a brief matter phase of the universe. This transition happens due to the dynamic characteristics of the k-essence field. Hence k- essence is solely responsible for the non-accelerated matter phase. During this transition the interaction density is non-zero but highly sub-dominated. Hence radiation background fluid coupled non-minimally with k-essence can explain the early phase by the non-zero interaction between the two sectors, but at the late-time k-essence explains the accelerating scenario of the universe.

6 Conclusion

This paper explores the cosmological dynamics of a non-minimally coupled k-essence scalar field and a relativistic fluid using a variational approach. The variational approach is a valuable technique to handle the theories involving fields and fluids in a covariant way. The main reason for choosing the variational approach is that it can produce the dynamical equations of cosmological development from a very simple premise. This approach can also explain more radical type of interactions between a scalar field and matter fluid. We have studied this field-fluid coupled system where the fluid component is radiation or dust. The coupling term in the level of the action has been chosen as \(f(n,s,\phi ,X)J^{\mu }\partial _{\mu } \phi \), in which the interaction term f depends on both field and fluid parameters. We have explicitly coupled the space-time derivative of the field \(\phi \) with the fluid 4-velocity \( U^{\mu } \) which gives rise to a highly non-linear dynamical evolution. The non-linear effects are amplified in this model because for non-zero interaction energy density, any change in the field value immediately affects the fluid sector via the derivative coupling term. If the field energy density is rapidly changing in some phase then the fluid properties will also rapidly change in that phase. On the other hand if the field value does not change with time this kind of interaction vanishes. Due to the presence of the specific derivative dependent interaction term several thermodynamic quantities such as system’s temperature \( T_\textrm{tot} \) and chemical-free energy \( F_\textrm{tot} \) have been generalized. We have established that the total energy momentum of the field-fluid coupled system is conserved. The other matter quantities like the equation of state (EOS) and sound speed (\( c_{s}^2 \)) have been generalized properly.

The resolution of the tension related to the mismatch of the values of \(H_0\) in the high and low-redshift measurements remains a difficult problem. Previous authors have tried several phenomenological approaches to address this problem using the concept of interacting dark energy scenario in which the coupling is introduced, by hand, at the level of the continuity equations (for energy densities). As these interactions are not introduced at the level of the action these theories may give rise to instabilities at the perturbation level [53]. In the present work, we have constructed a coupled dark sector, where the coupling is introduced at the Lagrangian level and consequently all the physical parameters can be consistently defined in the background and perturbation levels. From recent data analysis based on Planck and Hubble Space Telescope results, some theorists do think that a minimally coupled dark sector may not effectively reduce the \(H_0\) tension. A recent study on interacting dark matter with dark energy showed that non-minimal coupling in the dark sector can significantly reduce the \(H_0\) tension [54]. Following the success of interacting dark sector, we have tried to build a generalized model of derivative-dependent coupling in the field-fluid sector which relies on robust mathematical foundations. The physical significance of the presented models can be testified in the future when our predictions are matched with observational results.

We have constructed a dynamical system to harness the dynamics of the coupled field-fluid system. The primary variables have been selected as x, y, \(\sigma \) and z. Due to the presence of a constraint the dynamical system has three dimensional phase space. To compactly write the dynamical equations we have defined some other parameters. The other parameters are expressed in terms of the variables x, y, \(\sigma \) and z. Hence in general we have constructed the 3-D autonomous system which encodes all the information regarding the system. To explore the dynamics we have constructed two models, Model I and Model II using the inverse square potential \(V(\phi ) \propto 1/\phi ^2\). In all the cases the the \(k-\)essence kinetic term \(F(X)=X^2-X\). To work out the dynamics we have used the simplest models of interaction although these simplest models do contain multiple constant parameters. There can be other models of interaction which may contain more constant parameters.

In model I we studied the system in presence of dust and radiation separately. It was seen that the critical points had various dependencies on the model parameters \(\delta \), m and \(\eta \). Some of the critical points are only dependent on \(\delta \). Those critical points were found out in a heuristic manner. Our method showed that some critical points lie in the phantom domain and others lie in the non-phantom region. We used the grand EOS and the sound speed limit to select the other model parameters. In presence of dust we have found the critical points that correspond to different phases of the universe. These critical points specify various phases of the universe, from radiation domination to accelerating and phantom phases. Due to the presence of the derivative dependent interaction term, the interaction energy density and scalar field density plays important role in the early phase of the universe, whereas scalar field density dominates over other components in the late phase. In presence of radiation, due to the interaction, the early phase becomes highly non-linear, and hence most physical quantities increase abruptly. As the system evolves and enters the radiation phase the fluid density dominates. The transition from the radiation to matter phase is abrupt and the system slowly enters into the late-time cosmic acceleration phase.

In model II the system reduces from a 3-D dynamical system into a 2-D system governed by x and y only. When the fluid present corresponds to dust we have selected the model parameter \(\delta = 45.11\), for which we found a set of critical points which can describe the evolution of the universe from the matter dominated phase to the phantom phase of the universe. We found a heteroclinic orbit in phase space which shows how the matter dominated phase evolves and is finally attracted towards a stable accelerating point. The phantom points are also stable, but we have not found any critical points corresponding to the early phase which connects to them. Hence these phantom fixed points become isolated. In the early phase the interaction density dominates and the system is in a radiation dominated phase. While at the matter phase, fluid density becomes dominant, and the field density dominates the late-time phase.

In presence of radiation the phase space has five critical points, out of which only one critical point shows the radiation epoch. This point evolves towards the accelerating point and makes a heteroclinic orbit. Other critical points become isolated stable points. We see from the evolution plot that the fluid density becomes dominant at an early epoch which signifies the radiation phase, and at some point, interaction density also increases and then decreases. So in the presence of radiation the behavior of the system becomes a little interesting in the early epoch but then the system enters an accelerating phase after a brief matter epoch. Hence from both models it becomes vivid that the late-time dynamics is dominated by the k-essence scalar field, but as we move towards the early phases, the non-linear effects of interaction play a more significant role. We see that the field-fluid model collectively explains different epochs of the universe, from radiation to accelerating and phantom epochs. The problem of cosmic coincidence is taken care of by the dynamics of the field-fluid coupling.

In the previous study with algebraic coupling in the dark sector [43] the authors showed that one can attain an accelerating phase of the universe in the late phase of cosmic evolution. In the present work, we have worked with the derivative coupling in the dark sector and obtained results which can be used for the late cosmic acceleration phase. In such a circumstance a comparison between the previous work [43] and the present work can be useful. First of all, we specify that both the works have somethings in common, namely the form of F(X) and the form of the potential function \(V(\phi )\) remains the same. There is another similarity between these two works, both aim to produce a stable, accelerating expansion phase in the future. We will show that the last point is obtained in two different ways in these two attempts. Except for these similarities, there are important differences at the level of the models. The forms of the functions, f, which accompany the couplings are functionally different in the two works and the parameters in \(V(\phi )\) are also chosen differently so that the different models yield proper cosmic evolution. The number and nature of the various critical points of the systems are affected by the nature of the different types of couplings used in the two works. The qualitative difference between the nature of the critical points primarily distinguishes the present paper from the previous one. The previous study focussed mainly on the dark sector whereas the present work includes non-minimal coupling of the dark energy sector with radiation. One must note that in the previous work, the late phase of the universe permits the non-minimal coupling in the dark sector whereas in the present case as the scalar field value saturates near the stable critical point, specifying the late phase of the universe, the non-minimal coupling ceases to play any significant role. In the final phase of the universe, dark energy and dark matter are two separate entities in the present models as the interaction is directly proportional to the rate of change of the scalar field. This is another important difference between these two attempts related to different couplings in the dark sector. There are many other differences between the two models based on two different couplings, some of these points are related to formal theoretical issues while others are related to the nature of cosmological dynamics. We end this discussion with one prominent difference in the level of cosmological dynamics which can have observational effects. In model II of the previous work, the evolution of the effective (or total) EOS and \(c_s^2\) show jumps before the final stable phase is obtained. These jumps specify a phase transition period. In the present work, the same parameters smoothly enter the late phase and no jumps or kinks are noticed. This fact may have interesting cosmological effects.

In this paper we have introduced a new k-essence field derivative dependent field-fluid coupling. In some of the previous works the derivative independent field-fluid coupling was studied. Those previous works involved quintessence scalar field. In the future one may attempt to work with more complicated field-fluid couplings in cosmology. We expect that our present project will play an important role in those future attempts.