1 Introduction

The last several decades have seen cosmology radically altered with unprecedented observational evidence, first with the discovery of a Universe that is accelerating [1, 2] due to some form of dark energy and more recently with the increasingly convincing Hubble tension [3,4,5,6,7]. This tension has arisen between local measurements of the Hubble constant \(H_0\) [8, 9] and predictions of this cosmological parameter from observations from the early Universe [10, 11] which require the use of the concordance model or \(\Lambda \)CDM. To a lesser extent, the tension also appears in other cosmological parameters related to large scale structure measurements [12,13,14], which has prompted various attempts in the community to resolve the possible issue using additional contributions from the matter sector, as well as renewed interest in physics beyond general relativity (GR) which are now becoming mainstream.

GR acts as the base gravitational model on which \(\Lambda \)CDM rests as a foundation. However, there exist many possible directions for modified gravity to work toward. Recently, there have been a plethora of possible observationally motivated theories in the literature [15,16,17,18] which are yet to show promise of behaving better than \(\Lambda \)CDM when confronted with observations. By and large, these are mostly built as correction terms to the Einstein–Hilbert action of GR [19, 20]. In these works, gravitational interactions continue to be described through the curvature associated with the Levi-Civita connection which is the sole source of curvature in GR [21]. However, this is not the only way to construct gravitational models. In recent years there has been growing interest in teleparallel gravity (TG) where torsion is considered as the base mode of interactions for the gravitational sector [22,23,24,25].

In TG, the Levi-Civita connection is replaced with the teleparallel connection [22, 26] which expresses geometric deformations through torsion rather than curvature. In this regime, all measures of curvature vanish identically such as the Ricci scalar , where as its regular curvature form does not vanish (over-circles represent quantities calculated with the Levi-Civita connection). By relating the connections together, a torsion scalar T can be produced which is equal to the regular Ricci scalar up to a boundary term, meaning that it will produce field equations that are dynamically equivalent to GR, also called the Teleparallel equivalent of General Relativity (TEGR). Thus, observations are indistinguishable between GR and TEGR. This boundary term is important because it embodies the fourth-order terms of the Ricci scalar, which are boundary terms in the Einstein–Hilbert action. Its only when generalizations such as \(f(\mathring{R})\) gravity are considered [20] do these terms impact the order of the field equations.

The division of the torsion scalar and boundary term means that a weaker Lovelock theorem is developed [27,28,29] where much more general theories that are second order can be produced such as f(T) gravity [25]. f(T) theories of gravity [30,31,32,33,34,35,36,37,38] emerged using the same rationale as \(f(\mathring{R})\) gravity [15, 19, 20] where the TEGR Lagrangian is generalized to an arbitrary function of the torsion scalar. These theories are generically second order and have shown promise in meeting some observational challenges in the literature [24, 39,40,41,42,43]. In this work, we explore the space of cosmological models that feature a scalar field. In curvature-based gravity, there has been extensive analyses of scalar couplings in this setting [44, 45] that have shown real promise in constructing viable models. An interesting collection of such models is that of Horndeski gravity [46] in which a single scalar field is allowed to couple arbitrarily to the Einstein–Hilbert action provided that it produces second order field equations [47,48,49]. In the TG regime, a teleparallel analogue of Horndeski gravity has been proposed in the literature [29] which has the added benefit that it allows for a much wider range of models that are consistent with recent measurements of the gravitational propagation speed [50], as well as hosting a vast array of models that are consistent with solar system tests [51]. These models also show a rich array of possible gravitational wave polarizations [52] which is a growing topic for gravitational wave astronomy. More recently, the teleparallel analogue of Horndeski gravity has been used to construct well motivated models using well-tempered cosmological methods in either Minkowski [53] or cosmological backgrounds [54]. Another interesting direction related to constructing models in this framework is that of using Noether symmetries to construct cosmology oriented models as was performed in Ref. [55].

In this study, we probe possible cosmological behaviours using a dynamical systems approach which can reveal the evolution of the individual models under consideration in the context of their potential to attain the known critical points of the Universe [56]. By taking a flat homogeneous and isotropic background solution, dynamical systems can be used to determine the number and nature of the critical points of the system which express whether these positions in the cosmic evolution are stable or not. In the literature, scalar couplings with torsion have shown interesting results such as Ref. [57] where a nontrivial scalar field is allowed to contribute while the kinetic term remains canonical. In the present study, we extend these efforts to allow for several models with a nontrivial kinetic term. Our aim is to assess whether this type of generalization can produce a dynamical system consistent with our expectations for the Universe. To that end, we first briefly review the literature on TG and the structure of the teleparallel analogue of Horndeski gravity in Sect. 2. The cosmological system is then introduced in Sect. 3, while the dynamical system analysis for the models is conducted in Sects. 4, 5. In these models we explore the impact of a power-law kinetic term coupled with the torsion scalar and another nonzero (for this background) scalar of the theory respectively. Finally, we summarize our results in Sect. 6 where we discuss how these results compare with the present literature.

2 Scalar-tensor teleparallel gravity

The curvature associated with GR and other models based on the Levi-Civita connection [16] is reformulated in TG through the teleparallel connection where torsion replaces curvature as the means by which gravity is expressed [22]. The origin of curvature in GR and related theories is not the metric, but rather the Levi-Civita connection (over-circles are used throughout to denote quantities determined using the Levi-Civita connection) which characterizes how any geometric deformation is to be characterized. On the other hand, TG characterizes gravitation as torsion through the teleparallel connection which is curvature-less and satisfies metricity [23, 58]. While the regular Riemann tensor does not vanish, its teleparallel analogue does, as do all quantitative measures of curvature, meaning that an entirely new approach to forming gravitational theory needs to be adopted (see reviews in Refs. [23,24,25]).

The most efficient path to forming teleparallel theories of gravity is through the tetrad (and its inverses ) which takes the place of the metric as the fundamental variable of theory through the relations

(1)

where Latin indices represent coordinates on the tangent space while Greek indices continue to represent indices on the general manifold [24]. While they do appear in GR, the direct use of tetrads in GR is largely suppressed [59]. Similar to the metric, the tetrads must satisfy orthogonality conditions which take of the form of

(2)

preserving internal consistency.

The teleparallel connection can be directly defined as [26, 60]

(3)

where contributions of the tetrad are complemented by which is a flat spin connection and is responsible for incorporating local Lorentz invariance into teleparallel theories. This arises due to the explicit appearance of the Lorentz indices and thus the Lorentz frames. As tetrads appear in GR, so too do spin connections but one important distinction is that they are not flat in the GR case [21]. The tetrad-spin connection pair represent the gravitational and local degrees of freedom in the equations of motion of TG. Now, analogous to the way that the Levi-Civita connection builds up to the Riemann tensor, the teleparallel connection directly leads to the torsion tensor [58]

(4)

where square brackets denote the anti-symmetry operator, and where acts as the field strength of gravity [23]. This tensor is covariant under both diffeomorphisms and local Lorentz transformations. By an appropriate combination of contractions of torsion tensors, a torsion scalar can be written down such that [22,23,24,25]

(5)

which is the result of a requirement that T be equivalent to the curvature scalar \(\mathring{R}\) (up to a boundary term). Similar to the way in which the curvature scalar is dependent only on the Levi-Civita connection, the torsion tensor is only dependent on the teleparallel connection.

Exchanging the Levi-Civita connection with the teleparallel connection means that measures of curvature identically vanish, such as \(R\equiv 0\) (where we emphasize that and ). In this context, we can write the following relation for the curvature and torsion scalars [39, 61]

$$\begin{aligned} R=\mathring{R} + T - B = 0. \end{aligned}$$
(6)

where is a total divergence term, where is the determinant of the tetrad. This guarantees that the GR and TEGR actions will generate identical field equations.

Taking a similar rationale as with many other extensions to GR, such as \(f(\mathring{R})\) gravity [20, 62], TEGR can be directly generalized to f(T) by raising the TEGR Lagrangian to an arbitrary function [30,31,32,33,34, 63]

$$\begin{aligned} \mathcal {S}_{\mathcal {F}(T)}^{} = \frac{1}{2\kappa ^2}\int \mathrm {d}^4 x\; e\left( -T + \mathcal {F}(T)\right) + \int \mathrm {d}^4 x\; e\mathcal {L}_{\text {m}}, \end{aligned}$$
(7)

where \(\kappa ^2=8\pi G\), and \(\mathcal {L}_{\text {m}}\) is the matter Lagrangian in the Jordan frame. The tetrad and spin connection are independent variables in TG and thus produce independent field, which correspond to the ten metrical field equations and the six local Lorentz degrees of freedom. However, the tetrad variation has an interesting property that when acted on by the anti-symmetric operator, it produces the spin connection field equations. Thus, all the equations of motion can be determined with this variation as [22]

$$\begin{aligned} W_{(\mu \nu )} = \kappa ^2 \Theta _{\mu \nu }, \quad \text {and} \quad W_{[\mu \nu ]} = 0, \end{aligned}$$
(8)

which also holds for any other teleparallel gravity theory, and where \({{\Theta }^{\nu }_\rho }\) is the regular energy-momentum tensor for matter. In this setting, the Weitzenböck gauge can defined as the tetrad choice in which the spin connection vanishes, i.e. where the anti-symmetric field equations are identically zero for the choice (1) of tetrad.

In the present work, we explore the dynamical systems of a particular class of scalar-tensor models which uniquely appear in teleparallel gravity. To form the broader class of scalar-tensor extensions in TG, we first consider the irreducibles pieces of the torsion tensor [58, 64]

(9)
(10)
(11)

which are respectively the axial, vector, and purely tensorial parts, and where \(\epsilon _{\mu \nu \lambda \rho }\) is the totally antisymmetric Levi-Civita tensor in four dimensions. Taking appropriate contractions leads to the scalar invariants [61]

(12)
(13)
(14)

These three scalars form the most general parity preserving scalars that are quadratic in contractions of the torsion tensor, and even reproduce the torsion scalar \(T\,{:=}\,\frac{3}{2}T_{\text {ax}}+\frac{2}{3}T_{\text {ten}}-\frac{2}{3}T{_{\text {vec}}}\). Recently, this has led to the proposal of a teleparallel analogue of Horndeski gravity [29, 50,51,52,53,54,55], also called Bahamonde–Dialektopoulos–Levi Said (BDLS) theory. As in curvature-based gravity, this is grounded on the Lovelock theorem [27, 28, 65], and leads to the linear torsion scalar contraction scalar invariants [29]

$$\begin{aligned} I_2 = v^{\mu } \phi _{;\mu }, \end{aligned}$$
(15)

where \(\phi \) is the scalar field, and while for the quadratic scenario, we find

$$\begin{aligned} J_{1}&=a^{\mu }a^{\nu }\phi _{;\mu }\phi _{;\nu }, \end{aligned}$$
(16)
$$\begin{aligned} J_{3}&=v_{\sigma }t^{\sigma \mu \nu }\phi _{;\mu }\phi _{;\nu }, \end{aligned}$$
(17)
(18)
(19)
(20)
(21)

where semicolons represent covariant derivatives with respect to the Levi-Civita connection. The Levi-Civita connection enters into the scalar field sector through the minimal coupling prescription of TG (See Ref. [29] for further details).

Naturally, the regular Horndeski terms from curvature-based gravity also appear in this framework [46]

$$\begin{aligned} \mathcal {L}_{2}&:= G_{2}(\phi ,X), \end{aligned}$$
(22)
$$\begin{aligned} \mathcal {L}_{3}&:= G_{3}(\phi ,X)\mathring{\Box }\phi , \end{aligned}$$
(23)
$$\begin{aligned} \mathcal {L}_{4}&:= G_{4}(\phi ,X)\left( -T+B\right) \nonumber \\&\quad +G_{4,X}(\phi ,X)\left( \left( \mathring{\Box }\phi \right) ^{2}-\phi _{;\mu \nu }\phi ^{;\mu \nu }\right) , \end{aligned}$$
(24)
(25)

where the kinetic term is defined as \(X\,{:=}\,-\frac{1}{2}\partial ^{\mu }\phi \partial _{\mu }\phi \). BDLS theory simply adds the further Lagrangian component [29]

$$\begin{aligned} \mathcal {L}_{\text {{ Tele}}}\,{:=}\, G_{\text {{ Tele}}} \left( \phi ,X,T,T_{\text {ax}},T_{\text {vec}},I_2,J_1,J_3,J_5,J_6,J_8,J_{10}\right) .\nonumber \\ \end{aligned}$$
(26)

This results in the BDLS action given by

$$\begin{aligned} \mathcal {S}_{\text {BDLS}}= & {} \frac{1}{2\kappa ^2}\int d^4 x\, e\mathcal {L}_{\text {{ Tele}}} \nonumber \\&+ \frac{1}{2\kappa ^2}\sum _{i=2}^{5} \int d^4 x\, e\mathcal {L}_i+ \int d^4x \, e\mathcal {L}_\mathrm{m}, \end{aligned}$$
(27)

with \(\mathring{G}_{\mu \nu }\) is the standard Einstein tensor. The curvature-based regular Horndeski theory is recovered for the limit where \(G_{\text {{ Tele}}} = 0\). BDLS theory is invariant under local Lorentz transformations and diffeomorphisms due to it being based on the torsion tensor. One minor difference with regular Horndeski theory is that calculations are now based on the tetrad and spin connection components rather than the metric tensor, but this will not impact the values for the \(\mathcal {L}_{2} - \mathcal {L}_{5}\) contributions. Another important point to highlight is that this version of the popular scalar-tensor theory provides a much more general framework on which to construct models. Indeed, this allows for a path to circumvent the strong constraints imposed in regular Horndeski gravity from gravitational wave observations [66].

3 Teleparallel scalar-tensor flat FLRW cosmology

We consider a flat isotropic and homogeneous background cosmology through the Friedmann–Lemaître–Robertson–Walker (FLRW) metric [21]

$$\begin{aligned} \mathrm{d}s^2 = - \mathrm{d}t^2 + a(t)^2(\mathrm{d}x^2 + \mathrm{d}y^2 + \mathrm{d}z^2), \end{aligned}$$
(28)

Where a(t) is the scale factor depends on cosmic time t. This can be described by the tetrad

(29)

which is consistent with the Weitzenböck gauge described in Sect. 2. We take the standard definition of the Hubble parameter \(H= \frac{\dot{a}}{a}\), where dots refer to derivatives with respect to cosmic time. We also consider the equation of state (EoS) for matter \(\omega _{m}=\frac{p_\mathrm{m}}{\rho _\mathrm{m}} = 0\) and radiation \(\omega _{r}=\frac{p_{r}}{\rho _\mathrm{r}} = 1/3\), which will both contribute to our representation of cosmology. In this work, we consider the class of models in which

$$\begin{aligned} G_2&= X - V(\phi ), \end{aligned}$$
(30)
$$\begin{aligned} G_3&= 0 = G_5, \end{aligned}$$
(31)
$$\begin{aligned} G_4&= 1/2\kappa ^2, \end{aligned}$$
(32)

where we take a generalization of a canonical scalar field together with a TEGR term. This then lets use probe different forms of the \(G_\mathrm{Tele}\) term in the action (27). As discussed, our aim is to probe the nature of power-law couplings with the kinetic term. To that end, we consider two models that embody nonvanishing terms for an FLRW background cosmology, which are

$$\begin{aligned} G_{\mathrm{Tele}_1}&= X^{\alpha } T, \end{aligned}$$
(33)
$$\begin{aligned} G_{\mathrm{Tele}_2}&= X^{\alpha } I_2, \end{aligned}$$
(34)

where the other terms effectively do not contribute to the Friedmann equations [29], and where

$$\begin{aligned} T&= 6H^2, \end{aligned}$$
(35)
$$\begin{aligned} I_2&= 3H\dot{\phi }. \end{aligned}$$
(36)

Thus, we can write the effective Friedmann equations as

$$\begin{aligned} \dfrac{3}{\kappa ^{2}}H^{2}&= \rho _\mathrm{m} + \rho _\mathrm{r} + X + V + 6H\dot{\phi } G_{\mathrm{Tele},I_2} \nonumber \\&\quad + 12H^2 G_{\mathrm{Tele},T} + 2X G_{\mathrm{Tele},X} - G_{\mathrm{Tele}},\end{aligned}$$
(37)
$$\begin{aligned} -\dfrac{2}{\kappa ^{2}}\dot{H}&= \rho _\mathrm{m} + \frac{4}{3}\rho _\mathrm{r} \nonumber \\&\quad + 2X + 3H\dot{\phi }G_{\mathrm{Tele},I_2} + 2XG_{\mathrm{Tele},X} - \frac{d}{dt}\nonumber \\&\quad \times \left( 4HG_{\mathrm{Tele},T} + \dot{\phi } G_{\mathrm{Tele},I_2}\right) , \end{aligned}$$
(38)

and where the scalar field equation is given by

$$\begin{aligned}&\frac{1}{a^3}\frac{d}{dt}\left[ a^3 \dot{\phi } \left( 1+ G_{\mathrm{Tele},X}\right) \right] = -V'(\phi ) - 9H^2 G_{\mathrm{Tele},I_2} \nonumber \\&\quad + G_{\mathrm{Tele},\phi } - 3\frac{d}{dt}\left( HG_{\mathrm{Tele},I_2}\right) . \end{aligned}$$
(39)

This effective fluid observes the energy-conservation equation

$$\begin{aligned} \dot{\rho }_{DE}+3H(\rho _\mathrm{DE}+p_\mathrm{DE}) = 0, \end{aligned}$$
(40)

but this description breaks down at perturbative level. Moreover, in what follows we will use the density parameters

$$\begin{aligned} \Omega _{m}=\frac{\kappa ^2\rho _m}{3H^2}, \quad \Omega _\mathrm{DE}=\frac{\kappa ^2\rho _\mathrm{DE}}{3H^2}, \quad \Omega _{r}=\frac{\kappa ^2\rho _r}{3H^2}. \end{aligned}$$
(41)

which satisfies the conservation relation

$$\begin{aligned} \Omega _m + \Omega _\mathrm{DE} + \Omega _r = 1. \end{aligned}$$
(42)

4 Model 1 – kinetic term coupled with torsion scalar

The action for Model 1 is given by a coupling between a power-law-like term and the torsion scalar represented by

$$\begin{aligned} S = \int d^4x e [X-V(\phi )-\frac{T}{2\kappa ^{2}}+X^{\alpha }T]+S_{m}+S{r}, \end{aligned}$$
(43)

where \(V(\phi )\) is the scalar potential, \(S_m\) represents the action of for matter and \(S_r\) describes the action for the radiation component. We shall perform the dynamical system analysis for the general case \(\alpha \), followed by two examples with \(\alpha =1\), \(\alpha =2\). Taking background FLRW cosmology (29) and the action above, we can obtain the Friedmann equations in Eqs. (44–45) whereas the Klein–Gordon equation can be obtained in Eq. (46), altogether giving

$$\begin{aligned}&\frac{T}{2\kappa ^{2}}-V(\phi )-X-X^{\alpha }T-2\alpha X^{\alpha }T = \rho _\mathrm{m}+\rho _\mathrm{r}, \end{aligned}$$
(44)
$$\begin{aligned}&-V(\phi )+\frac{T}{2\kappa ^{2}}+ X-X^{\alpha }T+\frac{2\dot{H}}{\kappa ^{2}}-4X^{\alpha }\dot{H}-8\alpha X^{\alpha }H\frac{\ddot{\phi }}{\dot{\phi }} = -p_{r}, \end{aligned}$$
(45)
$$\begin{aligned}&V^{'}(\phi )+3H\dot{\phi }+\frac{6X^{\alpha }T\alpha H}{\dot{\phi }}+\frac{4\alpha X^{\alpha }T \dot{H}}{\dot{\phi } H}\nonumber \\&\quad +\ddot{\phi }[1-\alpha X^{\alpha -1}T+\alpha ^{2}X^{\alpha -2}T] = 0. \end{aligned}$$
(46)

Now using the background expressions defined in Sect. 3, we can obtain the expression for energy density and pressure of effective dark energy as

$$\begin{aligned} \rho _\mathrm{DE}&= X^{\alpha }T+2\alpha X^{\alpha }T+V(\phi )+X, \end{aligned}$$
(47)
$$\begin{aligned} p_\mathrm{DE}&= -V(\phi )+X-X^{\alpha }T-4X^{\alpha }\dot{H}-\frac{8\alpha X^{\alpha } H \ddot{\phi }}{\dot{\phi }}. \end{aligned}$$
(48)

To study the phases of cosmic evolution, the autonomous dynamical system for the above set of cosmological expressions can be defined using following dimensionless variables.

$$\begin{aligned} x&=\dfrac{\kappa \dot{\phi }}{\sqrt{6}H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3}H},\quad u=2 X^\alpha \kappa ^2,\quad \rho = \frac{\kappa \sqrt{\rho _\mathrm{r}}}{\sqrt{3}H},\nonumber \\ \lambda&=\frac{-V^{'}(\phi )}{\kappa V(\phi )},\quad \Gamma =\dfrac{V(\phi )V^{''}(\phi )}{V^{'}(\phi )^{2}}, \end{aligned}$$
(49)

where the constraint equation for the dimensionless variables can be obtained as,

$$\begin{aligned} x^{2}+y^{2}+(1+2\alpha )u+\Omega _{m}+\rho ^{2}=1. \end{aligned}$$
(50)

Now the autonomous dynamical system can be defined by differentiating the dimensionless variables with respect to \(N=\ln a\) as,

$$\begin{aligned} \dfrac{dx}{dN}&= \frac{x \left( -x^2 \left( \rho ^2+3 (\alpha (2 \alpha +5)+1) u-3 y^2-3\right) +\sqrt{6} \lambda x y^2 (2 \alpha u+u-1)-3 x^4\right) }{2 (u-1) x^2-2 \alpha u (2 \alpha (u+1)+u-1)}\nonumber \\&\quad - \frac{\alpha u x \left( 2 \alpha \left( \rho ^2+3\right) +\rho ^2+(6 \alpha +3) u-3 (2 \alpha +1) y^2-3\right) }{2 (u-1) x^2-2 \alpha u (2 \alpha (u+1)+u-1)}, \end{aligned}$$
(51)
$$\begin{aligned} \dfrac{dy}{dN}&= \frac{-y \left( x^2 \left( \rho ^2+\left( 6 \alpha ^2+9 \alpha -3\right) u-3 y^2+3\right) -2 \sqrt{6} \alpha \lambda u x y^2+3 x^4\right) }{2 (u-1) x^2-2 \alpha u (2 \alpha (u+1)+u-1)} \nonumber \\&\quad + \frac{\alpha u (-y) \left( (6 \alpha +3) u-(2 \alpha -1) \left( -\rho ^2+3 y^2-3\right) \right) }{2 (u-1) x^2-2 \alpha u (2 \alpha (u+1)+u-1)}-y\sqrt{\frac{3}{2}} \lambda x, \end{aligned}$$
(52)
$$\begin{aligned} \dfrac{du}{dN}&= \frac{\alpha u \left( 2 \alpha u \left( \rho ^2+3 x^2-3 y^2\right) +(u-1) x \left( 6 x-\sqrt{6} \lambda y^2\right) \right) }{\alpha u (2 \alpha (u+1)+u-1)-(u-1) x^2}, \end{aligned}$$
(53)
$$\begin{aligned} \frac{d\rho }{dN}&= \frac{\rho \left( -x^2 \left( \rho ^2+6 \alpha ^2 u+9 \alpha u+u-3 y^2-1\right) +2 \sqrt{6} \alpha \lambda u x y^2-3 x^4\right) }{2 (u-1) x^2-2 \alpha u (2 \alpha (u+1)+u-1)} \nonumber \\&\quad + \frac{\alpha \rho u \left( 2 \alpha u+u+(2 \alpha -1) \left( -\rho ^2+3 y^2+1\right) \right) }{2 (u-1) x^2-2 \alpha u (2 \alpha (u+1)+u-1)}, \end{aligned}$$
(54)
$$\begin{aligned} \frac{d\lambda }{dN}&= \sqrt{6}(\Gamma -1)\lambda ^{2}x. \end{aligned}$$
(55)

Unless the parameters \(\Gamma \) is known, the dynamical systems presented in this work are not an autonomous systems. We will now on wards focus on the exponential potential \(V(\phi ) =V_{0}e^{-\lambda \kappa \phi }\) with \(\lambda \) is a dimensionless constant. This particular form of potential function leads to, \(\Gamma = 1\) and can give rise to an accelerated expansion of the universe. We obtain the critical points (or fixed points) for autonomous dynamical system presented in Eqs. (51–55) by imposing conditions \(\frac{dx}{dN}=0\), \(\frac{dy}{dN}=0\), \(\frac{du}{dN}=0\), \(\frac{d\rho }{dN}=0\). The critical points are titled with capital letters and presented in corresponding tables. To study the cosmological implications, the value of the deceleration parameter and value for the total EoS \(\omega _{tot}\) are also presented in the tables. From the Table 1 observation, it can be concluded that parameter \(\alpha \) contributes in the co-ordinates of critical points D and E and represents the de Sitter solution for the dynamical system presented in Eqs. (51–55). While analyzing the critical points for the case general \(\alpha \), there is a chance to get more number of critical points than for a particular value of \(\alpha \). The critical points L, M, F, G and A represent same cosmological implication. These critical points show deceleration parameter value \(q=\frac{1}{2}\) and \(\omega _{tot}=0\), hence explaining the cold dark matter-dominated era. Similarly the critical points B, C, and N represent the same phase of evolution with value of \(\omega _{tot} =1\), hence cannot describe current accelerated phase of evolution and behave as stiff matter. The critical points J and K are defined at \(\lambda =2\) and show value of \(\omega _{tot}=\frac{1}{3}\) hence describe the radiation dominated era. Since the critical points H and I represent deceleration parameter value, \(q=-1+\frac{\lambda ^2}{2}\), these critical points can describe current acceleration of the universe for any real value of \(\lambda \) and are compatible with current observational data.

The stability of critical points can be studied by obtaining eigenvalues of linear perturbation matrix at critical points. Depending upon the signature of eigenvalues one can classify the stability properties as, if all the eigenvalues possesses positive signature it is an unstable node; if all the eigenvalues possesses negative signature then it is a stable node; if one of the eigenvalue possesses positive signature and other possesses negative sign in this case it is saddle point and if the determinant of linear perturbation matrix is negative and the real parts of all the eigenvalues possesses negative signature then it is a stable spiral. The eigenvalues and stability conditions for dynamical system in Eqs. (51–55) are presented in Table 2. The existence of positive and negative eigenvalues for permutation matrix at the critical points A, B, C, J, K and N describe saddle point behaviour at these critical points hence these critical points are unstable. However value of deceleration parameter at these critical points clarify that these critical points can not explain accelerated expansion phase of the universe. Critical points F and G ensure stability in the range of parameter \(\alpha >0\) and \(\left( -2 \sqrt{\frac{6}{7}}\le \lambda<-\sqrt{3}\vee \sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\right) \) these critical points explain the standard matter dominated era. The critical points D and E shows stable behaviour and explain the de Sitter solution.

Table 1 Critical points for dynamical system corresponding to Model-I, for general \(\alpha \)
Full size table
Table 2 Eigenvalues and stability of eigenvalue at corresponding critical point
Full size table

For any real value of \(\lambda \), the critical points H and I represent dark energy dominated era and this point shows stable behaviour for the parameters obey the value \(\alpha >0\) and \((-\sqrt{3}<\lambda<0\vee 0<\lambda <\sqrt{3})\). During the study of stability conditions, for critical points F, G, H and I, we get condition on \(\alpha \), \(\alpha >0\). This implies stability of critical points can be assessed for positive value of \(\alpha \) for Model-I. All the eigenvalues at critical points L and M are less than zero for \(\lambda \in \mathbb {R}\wedge \lambda \ne 0\wedge \frac{7 \lambda ^2-24}{7 \lambda ^2}\le \chi <\frac{\lambda ^2-3}{\lambda ^2}\), hence show stable behaviour at these parametric values.

Fig. 1
figure 1

Phase portrait for above dynamical system, the upper left plot is for \(u=0\), \(\rho =0\), \(\lambda =\sqrt{\frac{2}{9}}\). The upper right plot having parameter values \(u=0\), \(\rho =0\), \(\zeta =\frac{1}{9}\), \(\tau =1\), \(\alpha =1.1\), lower left phase portrait is for \(u=0\), \(\rho =0\), \(\lambda =\sqrt{\frac{2}{9}}\), lower right phase portrait is for \(u=0\), \(\rho =0\), \(\delta =1\)

Full size image

We have analyse the phase space for all the critical points by fixing some parameters to a appropriate value. The phase space plots for dynamical system described in Eqs. (51–55) are presented in Fig. 1. From phase space diagram for Model-I, we can conclude that the phase space trajectories for critical points L, F, M, G, A, J, N and K are moving away from the critical point hence confirm saddle point behaviour. Critical points F and G show stability for \(-2 \sqrt{\frac{6}{7}}\le \lambda <-\sqrt{3}\) or \(\sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\) but we choose \(\lambda =\sqrt{\frac{2}{9}}\) hence these are showing saddle points nature in phase diagram. For the critical points L and M the stability conditions depend on \(\chi \) which represents u co-ordinate and the phase plots are analysed in xy-axis plane hence critical points L and M may show saddle point behaviour. From the phase diagram, it can observe that at critical points H, I, D and E trajectories are attracted towards the critical point hence describe attracting behaviour of these critical points, also these critical points can explain dark energy dominated universe. Although eigenvalues at B and C contains both positive and negative signature, from the phase portrait critical points B and C represent an unstable node leading to the positive eigenvalues only (due to consideration of \(u=0\), \(\rho =\) may the negative eigenvalues are not contributing in the phase space plot). We shall present below two examples of this model for \(\alpha =1\) and \(\alpha =2\).

4.1 Case A: \(\alpha =1\)

In this case we have consider \(\alpha =1\), in the action equation Eq. (43). The evolution equations can be obtained by limiting Eqs. (44–46) from general \(\alpha \) to \(\alpha =1\). To study cosmic evolution through dynamical system analysis approach, the set of dimensionless variables associated with the above set of cosmological equations can be defined as follow [57]

$$\begin{aligned} x= & {} \dfrac{\kappa \dot{\phi }}{\sqrt{6}H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3}H}, u=3\dot{\phi }^2 \kappa ^2,\quad \rho = \frac{\kappa \sqrt{\rho _\mathrm{r}}}{\sqrt{3}H},\nonumber \\ \lambda= & {} \frac{-V^{'}(\phi )}{\kappa V(\phi )},\quad \Gamma =\dfrac{V(\phi )V^{''}(\phi )}{V^{'}(\phi )^{2}}. \end{aligned}$$
(56)

In this case the dimensionless variables are selected such that they can linked with each other in following constraint equation form, note that in this expression u is considered as it is (without any scalar multiplier).

$$\begin{aligned} x^{2}+y^{2}+u+\Omega _{m}+\rho ^{2}=1 \end{aligned}$$
(57)

Using these dimensionless variables the cosmological expressions in this case can be written in terms of autonomous dynamical system as follow,

$$\begin{aligned} \dfrac{dx}{dN}= & {} \frac{3 x y^2 \left( \sqrt{6} \lambda (u-1) x+3 u+3 x^2\right) -3 x \left( u^2+u \left( \rho ^2+8 x^2+1\right) +x^2 \left( \rho ^2+3 x^2-3\right) \right) }{2 (u-3) x^2-2 u (u+1)}, \end{aligned}$$
(58)
$$\begin{aligned} \dfrac{dy}{dN}= & {} -y \left( \frac{-3 u^2-u \left( \rho ^2+12 x^2+3\right) +u y^2 \left( 2 \sqrt{6} \lambda x+3\right) -3 x^2 \left( \rho ^2+3 x^2-3 y^2+3\right) }{2 u (u+1)-2 (u-3) x^2}+\sqrt{\frac{3}{2}} \lambda x\right) , \end{aligned}$$
(59)
$$\begin{aligned} \dfrac{du}{dN}= & {} \frac{2 u \left( \rho ^2 u+(6 u-9) x^2\right) -u y^2 \left( \sqrt{6} \lambda (u-3) x+6 u\right) }{u (u+1)-(u-3) x^2}, \end{aligned}$$
(60)
$$\begin{aligned} \frac{d\rho }{dN}= & {} \frac{\rho \left( -u^2+u \left( \rho ^2+16 x^2-y^2 \left( 2 \sqrt{6} \lambda x+3\right) -1\right) +3 x^2 \left( \rho ^2+3 x^2-3 y^2-1\right) \right) }{2 u (u+1)-2 (u-3) x^2}, \end{aligned}$$
(61)
$$\begin{aligned} \dfrac{d\lambda }{dN}= & {} -\sqrt{6}(\Gamma -1)\lambda ^{2}x. \end{aligned}$$
(62)
Table 3 Critical points for dynamical system corresponding to Model-I, \(\alpha =1\)
Full size table

From Table 3 observations, we can conclude that, at critical points A, F, G the value of deceleration parameter is \(\frac{1}{2}\) with \(\omega _{tot}=0\). These critical points not represent the accelerating universe, but the cold dark matter-dominated era. The critical points B and C behave as stiff matter showing \(\omega _{tot}=1\). The critical points J and K represent the radiation-dominated solutions. The critical points D, E, H and I show the value of the deceleration parameter negative, these critical points can represent the accelerating behavior of the universe. The critical points D and E are the de Sitter solutions with the value of \(\omega _{tot}=-1\) and can be obtained only at particular value of \(\lambda =0\). At the critical points H and I deceleration parameter shows negative value for \(-\sqrt{2}\) < \(\lambda \) < \(\sqrt{2}\), these points explains dark energy dominated universe. To analyse the stability behavior of all these critical points the eigenvalues and stability conditions are presented in Table 4

From Table 4, we can conclude that, the critical points D, E shows stable behaviour and these critical points can attract the universe at late time. At the critical points H and I eigenvalues show stability at \(-\sqrt{3}<\lambda <0\) or \(0<\lambda <\sqrt{3}\) these points explain dark energy domination at late time. The critical points A to C are saddle points and hence unstable for all values of \(\lambda \). However these points cannot describe the current accelerated expansion of the universe. The radiation dominated representation belong to the critical points J and K, these points are also saddle points for any value of \(\lambda \) and hence unstable in nature. Although critical points F and G represents cold dark matter dominated universe these critical point obey stability at \(-2 \sqrt{\frac{6}{7}}\le \lambda <-\sqrt{3}\) or \(\sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\).

Table 4 Eigenvalues and stability of corresponding critical point
Full size table

From the Fig. 2 observation, it can be conclude that critical points H, I are attractors and can be analysed by setting \(\lambda =\sqrt{\frac{2}{9}}\) which belongs to the stability range of \(\lambda \) (that is \(-\sqrt{3}<\lambda<0\vee 0<\lambda <\sqrt{3}\) ). Since F and G show stability for \(-2 \sqrt{\frac{6}{7}}\le \lambda <-\sqrt{3}\) or \(\sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\) and here, we have analyse phase plots at \(\lambda =\sqrt{\frac{2}{9}}\) which does not belong to stability range of \(\lambda \), the critical points F and G show saddle point behaviour. The critical points A, J, K having eigenvalues with both positive and negative sign and hence these are saddle points. The particular parametric value \(\lambda =2\) helps us to analyse the stability at J and K also from phase space analysis it can be observe that phase space trajectories are moving away at these critical points and hence show the saddle point behaviour. The de Sitter solution is represented by critical points D and E, these critical points are exists only for parametric value \(\lambda =0\) and from the phase space analysis it can conclude that these points represent attracting solution. The critical points A to K can obtained in particular case \(\alpha =1\) similar to general \(\alpha \), but the critical points L, M and N discussed in general \(\alpha \) case are not contributing in the case \(\alpha =1\).

Fig. 2
figure 2

Phase portrait for dynamical system presented in Eqs. (58–62), the upper left plot is for the parametric value \(u=0\), \(\rho =0\), \(\tau =1\) and \(\lambda =\sqrt{\frac{2}{9}}\). The upper right plot is for the parameteric values \(u=0\), \(\rho =0\) and \(\zeta =\frac{1}{9}\). The lower left phase portrait is for \(u=0\), \(\rho =0\) and \(\lambda =\sqrt{\frac{2}{9}}\) and lower right phase portrait is for \(u=0\) and \(\rho =0\)

Full size image

4.2 Case B: \(\alpha =2\)

In this case we have analyse cosmological implications by using dynamical system approach for particular value \(\alpha =2\) in action Eq. (43). In this case, the set of dimensionless variables to obtain autonomous dynamical system can be defined as follow,

$$\begin{aligned} x= & {} \dfrac{\kappa \dot{\phi }}{\sqrt{6}H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3}H},\quad u=\frac{5}{2}\kappa ^2\dot{\phi }^{4},\quad \rho = \frac{\kappa \sqrt{\rho _\mathrm{r}}}{\sqrt{3}H},\nonumber \\ \lambda= & {} \frac{-V^{'}(\phi )}{\kappa V(\phi )},\quad \Gamma =\dfrac{V(\phi )V^{''}(\phi )}{V^{'}(\phi )^{2}}. \end{aligned}$$
(63)

One can observe that the dimensionless variables defined in the study of these scalar tensor models are not same in [57], but these types of variables are usually used to obtain viable critical points in cosmology.. The dimensionless variables defined in Eq. (63) also satisfy the constraint equation Eq. (57) and the dynamical system in this case can be defined as follow,

$$\begin{aligned} \frac{dx}{dN}&= \frac{x \left( -5 \rho ^2 \left( 2 u+x^2\right) +5 y^2 \left( \sqrt{6} \lambda (u-1) x+6 u+3 x^2\right) +3 (5-19 u) x^2-6 u (u+3)-15 x^4\right) }{2 (u-5) x^2-4 u (u+3)}, \end{aligned}$$
(64)
$$\begin{aligned} \frac{dy}{dN}&= -y \left( \frac{-6 u^2-3 u \left( 2 \rho ^2+13 x^2+6\right) +2 u y^2 \left( 2 \sqrt{6} \lambda x+9\right) -5 x^2 \left( \rho ^2+3 x^2-3 y^2+3\right) }{4 u (u+3)-2 (u-5) x^2}+\sqrt{\frac{3}{2}} \lambda x\right) , \end{aligned}$$
(65)
$$\begin{aligned} \frac{du}{dN}&= \frac{4 u \left( 2 \rho ^2 u+3 (3 u-5) x^2\right) -2 u y^2 \left( \sqrt{6} \lambda (u-5) x+12 u\right) }{2 u (u+3)-(u-5) x^2}, \end{aligned}$$
(66)
$$\begin{aligned} \frac{d\rho }{dN}&= \frac{\rho \left( -2 u^2+u \left( 6 \rho ^2+43 x^2-2 y^2 \left( 2 \sqrt{6} \lambda x+9\right) -6\right) +5 x^2 \left( \rho ^2+3 x^2-3 y^2-1\right) \right) }{4 u (u+3)-2 (u-5) x^2}, \end{aligned}$$
(67)
$$\begin{aligned} \frac{d\lambda }{dN}&= -\sqrt{6}(\Gamma -1)\lambda ^{2}x. \end{aligned}$$
(68)

To study cosmological implications, the study of dynamical system at critical points obtained from cosmological evolution equations is very important. Critical points for autonomous dynamical system presented in Eqs. (64–68) are presented in Table 5. From the table observations it can conclude that, although critical points have different co-ordinates than the critical points in the Table 3 (for \(\alpha =1\)) case, but the cosmological implications are almost similar in nature. When we observe critical points in Tables 3 and 5, we can easily see that the deceleration parameter (q) and \(\omega _{tot}\) for critical points with the same name are same. From Table 5 it can be clearly observe that critical points D, E, H and I can show deceleration parameter value in negative range and hence these critical points can deal with the dark energy dominated era. For critical points H and I, we get accelerating behaviour for \(-\sqrt{2}\) < \(\lambda \) < \(\sqrt{2}\) and critical points D and E are defined only for parametric value \(\lambda =0\) and represent de Sitter solution for the system. The other critical points do not gives negative value for deceleration parameter and hence defines non-accelerating phase of evolution. The critical points A, F and G represents cold dark matter dominated era with \(\omega _{tot}=0\). In this case (for \(\alpha =2\)) also we are getting critical points B and C representing stiff matter. The critical points J and K are defined for \(\lambda =2\) and deliver value for \(\omega _{tot}=\frac{1}{3}\), hence represent radiation-dominated era.

Table 5 Critical points for dynamical system corresponding to Model-I, \(\alpha =2\)
Full size table

The stability conditions for critical points corresponding to dynamical system in Eqs. (64–68) are presented in Table 6. The signature of eigenvalues confirms the stability of corresponding critical point. From the table observations, we can conclude that critical points A, B and C are unstable for all values of \(\lambda \) and are saddle points. At the critical points H and I eigenvalues show stability at \(-\sqrt{3}<\lambda <0\) or \(0<\lambda <\sqrt{3}\) and these points explain dark energy domination at late time. Critical points D, E show stable behaviour and in the further analysis it is noticed that they can attract the universe at late time. The critical points J and K represent radiation dominated and from signature of eigenvalues these points are saddle points for any value of \(\lambda \), hence unstable in nature. Critical points F and G represents cold dark matter dominated universe and these critical point obey stability at \(-2 \sqrt{\frac{6}{7}}\le \lambda <-\sqrt{3}\) or \(\sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\). From Tables 4 and 6 it can observe that, the stability conditions for the case (\(\alpha =1\)) and (\(\alpha =2\)) are showing similar nature to explain the evolution of the universe.

Table 6 Eigenvalues and stability of eigenvalue at corresponding critical point
Full size table

The phase space diagram is presented in Figs. 3 and 4 for the parametric values \(\lambda =\sqrt{\frac{2}{9}}\) which belongs to the stability range for \(\lambda \) for critical points H and I and other parameters \(\tau =1\), \(\zeta =\frac{1}{9}\) are chosen such that the phase space diagram explain the stability conditions for the corresponding critical points. Critical points D and E, represent de Sitter solution the phase space analysis confirms the attracting nature of these critical points. From the Fig. 2 observations, it can be conclude that due to different co-ordinates the critical point A is presented in Fig. 3 moves in positive X-axis than in Fig. 2, but it does not impact on it’s stability nature. Critical points F and G show stability for \(-2 \sqrt{\frac{6}{7}}\le \lambda <-\sqrt{3}\) or \(\sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\) but we choose \(\lambda =\sqrt{\frac{2}{9}}\) hence these are a saddle points. The phase space diagram allows us to conclude that critical points J and K are saddle points.

Fig. 3
figure 3

Phase portrait for dynamical system in Eqs. (64–68), the left plot is for \(u=0\), \(\rho =0\), \(\tau =1\), \(\lambda =\sqrt{\frac{2}{9}}\), the right plot having parametric values \(u=0\), \(\rho =0\), \(\zeta =\frac{1}{9}\)

Full size image
Fig. 4
figure 4

Phase portrait for dynamical system in Eqs. (64–68), left phase portrait is for \(u=0\), \(\rho =0\), \(\lambda =\sqrt{\frac{2}{9}}\) and the right phase portrait is for \(u=0\), \(\rho =0\)

Full size image

5 Model 2 – kinetic term coupled with \(I_2\)

The model second action equation consist of a coupling between the \(X^\alpha \) term with \(I_2\), the action for Model 2 is described as below with

$$\begin{aligned} S = \int d^4x e \left[ X-V(\phi )-\frac{T}{2 \kappa ^{2}}+X^{\alpha }I_{2} \right] +S_{m}+S{r}, \end{aligned}$$
(69)

In the action equation, we have \(I_2=3H\dot{\phi }\) and other notations are same as in the first model. The expression for sum of energy density for matter and radiation and negative of pressure for radiation can be obtained on varying of action equation for Model 2 with respect to the tetrad field presented in Eqs. (70–71) respectively. The motion equation in this case can be obtained by taking variation of action equation with respect to \(\phi \), presented in Eq. (72).

$$\begin{aligned}&\frac{T}{2\kappa ^{2}}-V(\phi )-X-X^{\alpha }I_2-2\alpha X^{\alpha }I_2 = \rho _\mathrm{m}+\rho _\mathrm{r}, \end{aligned}$$
(70)
$$\begin{aligned}&\qquad -V(\phi )+\frac{T}{2\kappa ^{2}}+ X+\frac{2\dot{H}}{\kappa ^{2}}-X^{\alpha } \ddot{\phi }-2\alpha X^{\alpha }\ddot{\phi } = -p_{r}, \end{aligned}$$
(71)
$$\begin{aligned}&\quad V^{'}(\phi )+3H\dot{\phi }+\frac{(3+6\alpha )}{2}X^\alpha T + (3+6\alpha )X^\alpha \dot{H}\nonumber \\&\qquad +\ddot{\phi } \left[ 1+\frac{\alpha X^{\alpha }6H}{\dot{\phi }}+\frac{\alpha ^{2} X^{\alpha }12H}{\dot{\phi }} \right] = 0. \end{aligned}$$
(72)

Using background expressions discussed in Sect. 3, we can obtain expression for energy density and pressure for effective dark energy are presented in Eqs. (73) and (74) respectively.

$$\begin{aligned} \rho _\mathrm{DE}&= V(\phi )+X+X^{\alpha }I_2+2\alpha X^{\alpha }I_2, \end{aligned}$$
(73)
$$\begin{aligned} p_\mathrm{DE}&= -V(\phi )+X-(1+2\alpha )X^{\alpha }\ddot{\phi }, \end{aligned}$$
(74)

To define autonomous dynamical system, the dimensionless variables defined in this case are as follow

$$\begin{aligned} x= & {} \dfrac{\kappa \dot{\phi }}{\sqrt{6}H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3}H}, \quad u=\frac{\kappa ^2 X^\alpha \dot{\phi } }{H},\nonumber \\ \rho= & {} \frac{\kappa \sqrt{\rho _\mathrm{r}}}{\sqrt{3}H}, \quad \lambda =\frac{-V^{'}(\phi )}{\kappa V(\phi )}, \quad \Gamma =\dfrac{V(\phi )V^{''}(\phi )}{V^{'}(\phi )^{2}}. \end{aligned}$$
(75)

These dimensionless variables also satisfy constraint equation Eq. (50), and dynamical system can be obtained as presented below

$$\begin{aligned} \dfrac{dx}{dN}&= \frac{x \left( (2 \alpha +1) u \left( 2 \alpha \rho ^2+6 \alpha +\rho ^2+6 \alpha x^2+9 x^2-y^2 \left( 6 \alpha +\sqrt{6} \lambda x+3\right) -3\right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2} \nonumber \\&\quad +\frac{x \left( 3 (2 \alpha u+u)^2+2 x \left( 3 x^3+x \left( \rho ^2-3 y^2-3\right) +\sqrt{6} \lambda y^2\right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}, \end{aligned}$$
(76)
$$\begin{aligned} \dfrac{dy}{dN}&= -\sqrt{\frac{3}{2}} \lambda x y+\frac{y \left( 3 (2 \alpha u+u)^2+2 x^2 \left( \rho ^2+3 x^2-3 y^2+3\right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}\nonumber \\&\quad +\frac{(2 \alpha +1) u y \left( 6 (\alpha +1) x^2-\sqrt{6} \lambda x y^2+2 \alpha \left( \rho ^2-3 y^2+3\right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}, \end{aligned}$$
(77)
$$\begin{aligned} \dfrac{du}{dN}&= \frac{u \left( 3 (2 \alpha u+u)^2+(2 \alpha +1) u \left( 4 \alpha \rho ^2+\rho ^2+12 \alpha x^2+9 x^2-y^2 \left( 12 \alpha +\sqrt{6} \lambda x+3\right) -3\right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}\nonumber \\&\quad + \frac{2 u x \left( 3 x^3+x \left( \rho ^2-3 \left( 4 \alpha +y^2+1\right) \right) +\sqrt{6} (2 \alpha +1) \lambda y^2\right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}, \end{aligned}$$
(78)
$$\begin{aligned} \frac{d\rho }{dN}&= \frac{\rho \left( (2 \alpha u+u)^2+2 x^2 \left( \rho ^2+3 x^2-3 y^2-1\right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}\nonumber \\&\quad +\frac{\rho u \left( (2 \alpha +1) \left( 6 (\alpha +1) x^2-\sqrt{6} \lambda x y^2+2 \alpha \left( \rho ^2-3 y^2-1\right) \right) \right) }{(2 \alpha +1) u (2 \alpha (u+2)+u)+4 x^2}, \end{aligned}$$
(79)
$$\begin{aligned} \frac{d\lambda }{dN}&= \sqrt{6}(\Gamma -1)\lambda ^{2}x. \end{aligned}$$
(80)
Table 7 Critical points for dynamical system corresponding to Model-II, for general \(\alpha \)
Full size table
Table 8 Eigenvalues and stability of eigenvalue at corresponding critical point
Full size table

The critical points for the dynamical system described in Eqs. (76–80) are presented in Table 7. From table observations it can conclude that, critical points J to O represent similar cosmological implications in terms of deceleration parameter \(q=1\) and \(\omega _{tot}=\frac{1}{3}\) describe radiation dominated phase of the universe. Amongst these critical points the critical points N and O are defined for \(\alpha =\frac{1}{2}\). The critical points F, G and P also describe the same cosmological implication and represent cold dark matter dominated era with \(\omega _{tot}=0\). \(\alpha \) play role in the co-ordinate representation of critical points D and E, these critical points describe de Sitter solution and defined for \(\lambda =0\). The critical points H and I represent value for deceleration parameter \(q=-1+\frac{\lambda ^2}{2}\), these critical points can explain dark energy dominated universe. Critical points A, B and C deliver same value for q and \(\omega _{tot}\) hence these critical points also represent similar phase of universe evolution and behaves as stiff matter.

The stability conditions of the critical points are discussed in Table 8. The critical points J to O are presenting radiation dominated universe, eigenvalues for the linear perturbation matrix at these critical points having at least one eigenvalue with positive signature hence showing unstable behaviour. The critical points F, G show eigenvalues in negative range for \(\alpha >0\wedge (-2 \sqrt{\frac{6}{7}}\le \lambda<-\sqrt{3}\vee \sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}})\) and critical point P is stable at \(\left( \lambda<0\wedge \sigma <\frac{\sqrt{\frac{3}{2}}}{\lambda }\right) \vee \left( \lambda>0\wedge \sigma >\frac{\sqrt{\frac{3}{2}}}{\lambda }\right) \), hence show stability within this range, also these points express \(\omega _{tot}=0\), hence addressing cold dark matter phase of the universe evolution. The critical points B and C represent stiff matter era and possessing at least one eigenvalue with positive signature hence are unstable in nature, where critical point A show stable behaviour in the parametric range same as critical point P. The critical points D, E , H and I represent dark energy dominated era, where D and E are non-hyperbolic critical points but are stable in nature and critical points H and I show their stability in the parametric range \(\alpha >0\wedge \left( -\sqrt{3}<\lambda<0\vee 0<\lambda <\sqrt{3}\right) \). The stability conditions of critical points F, G, H and I show \(\alpha >0\) condition which implies stability of critical points representing cold dark matter and dark energy dominated era can be obtained for positive value of \(\alpha \).

Fig. 5
figure 5

Phase portrait for above dynamical system Eqs. (76–80), the upper left plot is for \(u=0\), \(\rho =0\), \(\lambda =\sqrt{\frac{2}{9}}\) and the upper right plot having parameter values \(x=0\), \(\rho =0\),, \(\sigma =1\), \(\alpha =1.1\). The lower left phase portrait is for \(x=0\), \(\rho =0\), \(\alpha =1.1\), \(\varphi =\frac{1}{2}\), the lower right phase portrait is for \(u=0\), \(\rho =0\), \(\eta =1\), \(\epsilon =\frac{1}{3}\)

Full size image

The critical points are plotted in the phase diagram Fig. 5. These phase plots are plotted for the dynamical system presented in Eqs. (76–80). The lower right plot shows that the phase space trajectories are moving away from critical points L, M, N and O hence these points represent instability with saddle point behaviour. The critical points N and O are defined for \(\alpha =\frac{1}{2}\). The critical points D and E are de Sitter solutions and phase diagram clarify that these points behaves as attracting solutions. The upper left phase plot describe phase space trajectories behaviour of critical points H, I, B, C, F and G. We can observe that the critical points B and C which is unstable saddle point but showing unstable node behaviour which is leading to the positive eigenvalues at critical points B and C. Also point H and I showing attracting point behaviour, these critical points represent dark energy dominated era with stability as described in Table 8. Critical points F and G represent cold dark matter dominated era, we have plotted plots for \(\lambda =\sqrt{\frac{2}{9}}\) which is not in the stability range of F and G. The upper right diagram represent phase space trajectories at critical points K, J and P, since the phase space trajectories are moving away from these critical points, these critical points are showing saddle point behaviour hence unstable. Since the parametric value \(\sigma =1\) do not follow stability range of critical point A hence getting saddle point behaviour.

5.1 Case A: \(\alpha =1\)

We have consider \(\alpha =1\) in the action equation Eq. (69) to discuss role of \(\alpha \) in getting number of critical points and their participation in describing different phases of Universe evolution. In this case, the evolution equations can be obtained by limiting the general \(\alpha \) to \(\alpha =1\) from Eqs. (70–72) the set of dimensionless variables can be defined as follow

$$\begin{aligned} x= & {} \dfrac{\kappa \dot{\phi }}{\sqrt{6}H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3}H},\quad u=\frac{3\kappa ^2\dot{\phi }^3}{2H},\quad \rho = \frac{\kappa \sqrt{\rho _\mathrm{r}}}{\sqrt{3}H},\nonumber \\ \lambda= & {} \frac{-V^{'}(\phi )}{\kappa V(\phi )},\quad \Gamma =\dfrac{V(\phi )V^{''}(\phi )}{V^{'}(\phi )^{2}}. \end{aligned}$$
(81)

These dimensionless variables satisfy constraint equation Eq. (57), and the dynamical system in this case can be defined as follow,

$$\begin{aligned} \dfrac{dx}{dN}= & {} \frac{x \left( 3 u^2+3 \rho ^2 u-y^2 \left( \sqrt{6} \lambda (u-2) x+9 u+6 x^2\right) +15 u x^2+3 u+2 \rho ^2 x^2+6 x^4-6 x^2\right) }{u (u+4)+4 x^2}, \end{aligned}$$
(82)
$$\begin{aligned} \dfrac{dy}{dN}= & {} -y \left( \frac{-3 u^2-2 u \left( \rho ^2+6 x^2+3\right) +u y^2 \left( \sqrt{6} \lambda x+6\right) -2 x^2 \left( \rho ^2+3 x^2-3 y^2+3\right) }{u (u+4)+4 x^2}+\sqrt{\frac{3}{2}} \lambda x\right) , \end{aligned}$$
(83)
$$\begin{aligned} \dfrac{du}{dN}= & {} \frac{u \left( 3 u^2+5 \rho ^2 u-y^2 \left( \sqrt{6} \lambda (u-6) x+15 u+6 x^2\right) +21 u x^2-3 u+2 \rho ^2 x^2+6 x^4-30 x^2\right) }{u (u+4)+4 x^2}, \end{aligned}$$
(84)
$$\begin{aligned} \frac{d\rho }{dN}= & {} \frac{\rho \left( u^2+2 u \left( \rho ^2+6 x^2-1\right) -u y^2 \left( \sqrt{6} \lambda x+6\right) +2 x^2 \left( \rho ^2+3 x^2-3 y^2-1\right) \right) }{u (u+4)+4 x^2}, \end{aligned}$$
(85)
$$\begin{aligned} \dfrac{d\lambda }{dN}= & {} -\sqrt{6}(\Gamma -1)\lambda ^{2}x. \end{aligned}$$
(86)

The critical points with value of deceleration parameter and \(\omega _{tot}\) are presented in the Table 9. From the table observations we can conclude that we get less number of critical points than the general case. The critical points J and K show deceleration parameter value \(q=1\), hence describe radiation dominated era. The critical points F and G are showing similar cosmological implication in terms of value of deceleration parameter and value of \(\omega _{tot}\), representing cold dark matter dominated era. Amongst all the critical points, an accelerated phase of evolution can be described by the critical points D, E, H and I. The critical points D and E describe de Sitter solution with \(q=-1\) and critical points H and I deliver deceleration parameter value \(q=-1+\frac{\lambda ^2}{2}\) hence represent the dark energy dominated era. The critical points B and C gives \(q=2\), these points can describe stiff matter. The critical point A showing deceleration parameter value in positive range, hence can not describe the accelerated phase of the evolution.

Table 9 Critical points for the dynamical system corresponding to updated calculations for Model-I
Full size table

To analyse the stability conditions eigenvalues for all critical points are presented in Table 10. The critical points F, G show stability for the parametric values \(-2 \sqrt{\frac{6}{7}}\le \lambda<-\sqrt{3}\vee \sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\) and describe cold dark matter dominated era. The critical points H and I are shows stable behaviour for the parametric values \(-\sqrt{3}<\lambda<0\vee 0<\lambda <\sqrt{3}\) these critical points represent dark energy dominated era for any real values of \(\lambda \). The critical points J and K are defined at \(\lambda =2\) and presence of eigenvalues with both positive and negative signature leads to saddle point hence unstable. The critical points B and C is also showing saddle point behaviour and unstable in nature. In this case, we get critical point A with deceleration parameter value \(q=\frac{4}{5}\) with existence of positive eigenvalues at linear perturbation matrix hence unstable in nature.

Table 10 Eigenvalues and stability of eigenvalue at corresponding critical point
Full size table
Fig. 6
figure 6

Left phase portrait is plotted for the parametric values \(u=0\), \(\rho =0\) and \(\lambda =\sqrt{\frac{2}{9}}\) and right side plot parametric values are \(x=0\), \(\rho =0\), \(\sigma =\frac{1}{4}\) for right side Fig. 

Full size image

The phase space plots for dynamical system presented in Eqs (82–86) are described in Figs. 6 and 7. The phase space analysis concludes that critical points H and I shows attracting nature, the accelerating expansion of the universe can be described at these critical points. The phase trajectories at critical points B and C are moving away from the critical point hence showing unstable node behaviour leading to the positive eigenvalues. The phase space plot for critical points F and G describe saddle point behaviour and represent cold dark matter dominated era. The phase plots at critical points D and E show attracting behaviour and these critical points describe de Sitter solution. The critical point A is a saddle point and can be confirmed by observing phase trajectories at A. The plot in Fig. 7 also describe that phase space trajectories are moving away from critical points hence showing unstable. behaviour.

5.2 Case B: \(\alpha =2\)

In this case we have discussed dynamical system analysis for Model-II, \(\alpha =2\). The evolution expressions can be obtained by limiting Eqs. (70–72) from general \(\alpha \) to \(\alpha =2\). The dimensionless variables to obtain autonomous dynamical system can be defined as follow

$$\begin{aligned} x= & {} \dfrac{\kappa \dot{\phi }}{\sqrt{6}H},\quad y=\frac{\kappa \sqrt{V}}{\sqrt{3}H},\quad u=\dfrac{5\kappa ^{2}\dot{\phi }^{5}}{4H},\quad \rho = \frac{\kappa \sqrt{\rho _\mathrm{r}}}{\sqrt{3}H},\nonumber \\ \lambda= & {} \frac{-V^{'}(\phi )}{\kappa V(\phi )},\quad \Gamma =\dfrac{V(\phi )V^{''}(\phi )}{V^{'}(\phi )^{2}}. \end{aligned}$$
(87)

These dimensionless variables satisfy constraint equation presented in Eq. (57), the dynamical system in this case is as follow

$$\begin{aligned} \frac{dx}{dN}= & {} \frac{x \left( 3 u^2+5 \rho ^2 u-y^2 \left( \sqrt{6} \lambda (u-2) x+15 u+6 x^2\right) +21 u x^2+9 u+2 \rho ^2 x^2+6 x^4-6 x^2\right) }{u (u+8)+4 x^2}, \end{aligned}$$
(88)
$$\begin{aligned} \frac{dy}{dN}= & {} -y \left( \frac{-3 u^2-2 u \left( 2 \rho ^2+9 x^2+6\right) +u y^2 \left( \sqrt{6} \lambda x+12\right) -2 x^2 \left( \rho ^2+3 x^2-3 y^2+3\right) }{u (u+8)+4 x^2}+\sqrt{\frac{3}{2}} \lambda x\right) , \end{aligned}$$
(89)
$$\begin{aligned} \frac{du}{dN}= & {} \frac{u \left( 3 u^2+9 \rho ^2 u-y^2 \left( \sqrt{6} \lambda (u-10) x+27 u+6 x^2\right) +33 u x^2-3 u+2 \rho ^2 x^2+6 x^4-54 x^2\right) }{u (u+8)+4 x^2}, \end{aligned}$$
(90)
$$\begin{aligned} \frac{d\rho }{dN}= & {} \frac{\rho \left( u^2+2 u \left( 2 \rho ^2+9 x^2-2\right) -u y^2 \left( \sqrt{6} \lambda x+12\right) +2 x^2 \left( \rho ^2+3 x^2-3 y^2-1\right) \right) }{u (u+8)+4 x^2}, \end{aligned}$$
(91)
$$\begin{aligned} \frac{d\lambda }{dN}= & {} -\sqrt{6}(\Gamma -1)\lambda ^{2}x. \end{aligned}$$
(92)

The critical points with value of deceleration parameter and \(\omega _{tot}\) for dynamical system in Eqs. (88–92) are presented in Table 11. From table observations in Model-II for critical point A we get different positive deceleration parameter value for different values of \(\alpha \). For \(\alpha =2\) we are getting \(q=\frac{2}{3}\) and \(\omega _{tot}=\frac{1}{9}\). The critical points D, E represent the de Sitter solution to the system and defined for \(\lambda =0\). While critical points H and H deliver deceleration parameter value \(q=-1+\frac{\lambda ^2}{2}\) which may describe dark energy dominated era of the universe. Critical points F and G represent cold dark matter dominated era with \(\omega _{tot}=0\). The critical points J and K are defined for \(\lambda =2\) and describe radiation dominated era of the universe evolution. The critical points B and C behaves as stiff matter with \(\omega _{tot}=1\).

The stability conditions for this case \(\alpha =2\) are presented in Table 12. Table observations conclude that critical points H and I show stability for parametric condition \(-\sqrt{3}<\lambda<0\vee 0<\lambda <\sqrt{3}\) and can describe dark energy dominated era. The critical points F and G show stability at \(-2 \sqrt{\frac{6}{7}}\le \lambda<-\sqrt{3}\vee \sqrt{3}<\lambda \le 2 \sqrt{\frac{6}{7}}\) and describe cold dark matter dominated era. The critical point A, B and C with deceleration parameter in positive range are unstable since eigenvalues are with both positive and negative signature. The critical points D and E represent de Sitter solution and show stable behaviour. The critical points J and K are representing radiation dominated era at \(\lambda =2\) and are unstable in nature.

The phase space diagram for dynamical system Eqs. (88–92) are plotted in Fig. 8. The critical points D and E showing attracting stable de Sitter solution and represent dark energy dominated era. While at critical point A phase space trajectories are moving away hence unstable addressing saddle point behaviour. Similarly critical points J and K are showing unstable behaviour representing radiation dominated era. The upper right plot is presented for critical points B, C, H, I, F and G. The critical points B and C behaves as unstable node with respect to the existence of positive eigenvalues. The phase space trajectories are moving away from critical points F and G, these critical points represent unstable behaviour and critical points H and I display attracting behaviour of phase space trajectories and show consistent with current observational studies.

Fig. 7
figure 7

The phase portrait is plotted for the parametric values \(\lambda =2\), \(u=0\), \(\rho =0\)

Full size image

6 Conclusion

In this work we have explored the cosmological dynamics of dark energy through the prism of scalar-torsion gravity [36, 37] in the context of power-law couplings with the kinetic term. In particular, we have studied one nontrivial extension of the recently proposed teleparallel analogue of Horndeski gravity [29]. This new frameworks offers a pathway to circumvent the severe constraints on the speed of gravitational wave constraint. This is possible due to the lower order nature of TG where the curvature associated with the Levi-Civita connection is exchanged with the torsion from the teleparallel connection. In the present case, we are interested in extending the analysis of dynamical systems by looking into how power-law-like couplings with the nonvanishing terms for background FLRW cosmologies. The teleparallel analogue of Horndeski theory offers a broad spectrum of possible extensions to regular Horndeski gravity and thus the space of possible functional forms are vast. This work goes some way to elucidating the behaviour that such functions adhere to.

For our FLRW background cosmology, we explore these models in the presence of both radiation and cold dark matter through the density parameters in Eq. (41). The Friedmann equations in Eqs. (37, 38) and Klein–Gordon equation in Eq. (39) directly lead to a set of autonomous equations for each of the models under investigation. These are then used in each case to derive the critical points of the particular cosmologies from which we can expose the model behaviour using the dynamical analysis in the parameter phase space. This also opens the doorway to understanding the stability of the models in question.

In the first model for the action in Eq. (43), we utilize the dynamical variables defined in Eq. (49), which using the constraint in Eq. (50) together with the equations of motion, are then used to derive the system of autonomous equations given in Eqs. (51–55) which express the behaviour of the model in phase space. The critical points are then arrived at by imposing that each of these derivatives vanishes. These first order equations of motion of the dynamical variables are represented as derivatives with respect to \(N = \ln a\) which shows the behaviour of the system in a more direct way. The result of this analysis is shown in Table 1. In this table we show the values of the dynamical variables at which these critical points occur together with the value of the deceleration and EoS parameters which already show an indication of the cosmological behaviour at those points in the evolution of the cosmological model. Each of these critical points is then further analyzed for their stability in Table 2. In some circumstances, stability occurs for a smaller change of parameter values as described in the last column of the table. For transparency we also show the corresponding Eigenvalues at each critical point. In Fig. 1, we show the phase portraits of this model for four specific examples of representative parameter values. In these plots the nature of the critical points is further exposed through their impact on the evolution contours.

Table 11 Critical points for the dynamical system corresponding to updated calculations for Model-I
Full size table
Table 12 Eigenvalues and stability of eigenvalue at corresponding critical point
Full size table
Fig. 8
figure 8

Phase portrait upper left is plotted for the parametric values \(u=0\), \(\rho =0\) and \(\lambda =\sqrt{\frac{2}{9}}\) and for upper right side plot parametric values are \(x=0\), \(\rho =0\), \(\sigma =\frac{1}{4}\). The lower phase portrait is for the parametric values \(u=0\), \(\rho =0\), \(\lambda =2\)

Full size image

We further probe the behaviour of this model in the specific cases of the kinetic term being linear and quadratic which represent the first cases that a Taylor expansion would open to. We do this in Sects. 4.1 and 4.2 respectively. We also show the phase portraits for these cases in Figs. 2, 3, 4 for the two respective cases. Finally, we show the nuanced critical points for both cases in Tables 3, 5 and their respective stability in Tables 4, 6.

The second model we explore the coupling between a power-law kinetic term and \(I_2\) scalar written in Eq. (36). This scalar represents the only nonvanishing term that is linear in its contractions with the torsion tensor. In this case, we take the action Eq. (69) where TEGR and the canonical scalar field are complemented by this new coupling term together with the matter and radiation contributions. This leads directly to the Friedmann equations in Eqs (70, 71), and the Klein–Gordon equation in Eq. (72). Now, by defining the dynamical variables in Eq. (75) we can explore the dynamical system that the background cosmology represents through the system of linear autonomous differential equations given in Eqs. (76–80). By performing a similar critical point analysis as in first case, we find the critical points as listed in Table 7. These are then further analyzed for their stability nature in Table 8. As in the first case, we find a rich structure of critical points for the various model parameter values, which we show through the phase portrait in Fig. 5. As in the first model, we consider the cases where the power index takes linear and quadratic forms in the model action. For the linear case, the dynamical system is represented by Eqs. (82–86) which produce the critical points in Table 9 with natures shown in Table 10. This interesting scenario produces the phase portraits given in Figs. 6 and 7. On the other hand, the quadratic case is represented through the dynamical system given in Eqs. (88–92). The corresponding critical point analysis produces Table 11 which have stability conditions described in Table 12. The final phase portraits for this case are then shown in Fig. 8.

The action presented Eq. (27) produces the most general second order theory that contains only one scalar field in TG. Naturally, this will produce a wealth of dynamics for which reason we explore two prominent models and expose their dynamical systems in the ensuing sections. While the models produce a vast array of critical points in Tables 1, 7, these are not all realized in each possible evolution of the individual models. Another aspect to highlight is the impact of boundary conditions on these potential evolution histories, that is, some critical points will not be accessible for certain boundary. This is further highlighted through the phase portrait examples where some critical points are common to all cosmological histories while other only appear for some iterations of the system. Thus, as explored in the interesting review in Ref. [56], while some behaviours remain common to the model produced by one of the actions such as a de Sitter future critical point, others only partially appear in the broad range of possible cosmological evolution histories available.

If we glance at a study made in Ref. [68], we can easily compare the cosmological implications and stability conditions of critical points from Table-1 in Ref. [68] and Table 1 for the Model-I general \(\alpha \) case. From this comparison, we can note that we get four more critical points than in Ref. [68], this may be possible because of the model construction and background formalism in the teleparallel Horndeski theory. This comparison allows us to know more about minor differences; in our study we get four critical points which are able to describe the dark energy era that is critical points H, I D, E. From these, the two critical points, D and E, are the stable de-Sitter solution (the study for the Model-II, Table 7 also showing similar comparison conclusion with sixteen critical points), but in Ref. [68], we will likely deal with four critical points explaining the dark energy era with only one de-Sitter solution, critical point “l.” With this progress we also have the number of critical points describing the matter-dominated solution with \(\omega =0\) is more than in Ref. [68]. These novel critical points plays an important role to get more clarity in the description of the matter and the dark energy dominated era. Also, we can compare our analysis with the study of various dark energy models in the modified theory of gravity in Ref. [56].