Dynamical system analysis of dark energy models in scalar coupled metric-torsion theories

TL;DR

This paper analyzes the phase space dynamics of dark energy models within scalar coupled metric-torsion theories, identifying stable solutions that can explain cosmic acceleration without direct interaction with matter.

Contribution

It introduces a detailed phase space analysis of non-minimal scalar-torsion couplings, identifying conditions for stable, accelerating solutions in both Einstein and Jordan frames.

Findings

Stable critical points correspond to accelerated expansion.

Evolution profiles match earlier exact solutions.

Physical viability depends on coupling parameter regions.

Abstract

We study the phase space dynamics of cosmological models in the theoretical formulations of non-minimal metric-torsion couplings with a scalar field, and investigate in particular the critical points which yield stable solutions exhibiting cosmic acceleration driven by the {\em dark energy}. The latter is defined in a way that it effectively has no direct interaction with the cosmological fluid, although in an equivalent scalar-tensor cosmological setup the scalar field interacts with the fluid (which we consider to be the pressureless dust). Determining the conditions for the existence of the stable critical points we check their physical viability, in both Einstein and Jordan frames. We also verify that in either of these frames, the evolution of the universe at the corresponding stable points matches with that given by the respective exact solutions we have found in an earlier work (arXiv: 1611.00654 [gr-qc]). We not only examine the regions of physical relevance for the trajectories in the phase space when the coupling parameter is varied, but also demonstrate the evolution profiles of the cosmological parameters of interest along fiducial trajectories in the effectively non-interacting scenarios, in both Einstein and Jordan frames.

Equations

$$\tilde{\mathcal{R}}:=\mathcal{R}-2\nabla {}_{\mu }{\mathcal{T}}^{\mu }-\frac{2}{3}{\mathcal{T}}_{\mu }{\mathcal{T}}^{\mu }+\frac{1}{24}{\mathcal{A}}_{\mu }{\mathcal{A}}^{\mu }+\frac{1}{2}{\mathcal{Q}}_{\mu \nu \sigma }{\mathcal{Q}}^{\mu \nu \sigma }.$$

$$\mathcal{S}=\int d^{4}x\sqrt{-g}[\frac{\beta {\varphi }^2}{2}\left(\mathcal{R}+4{\mathcal{T}}^{\mu }\frac{{\partial }_{\mu }\varphi }{\varphi }-\frac{2}{3}{\mathcal{T}}_{\mu }{\mathcal{T}}^{\mu }+\frac{1}{24}{\mathcal{A}}_{\mu }{\mathcal{A}}^{\mu }+\frac{1}{2}{\mathcal{Q}}_{\mu \nu \sigma }{\mathcal{Q}}^{\mu \nu \sigma }\right)$$

$$-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -V(\varphi )+{\mathcal{L}}^{(m)}],$$

$$\mathcal{S}=\int d^{4}x\sqrt{-g}\left[{\left(\frac{\varphi }{{\varphi }_{{}_0}}\right)}^2\frac{\mathcal{R}}{2{\kappa }^2}-\frac{\left(1-6\beta \right)}{2}{g}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -\frac{1}{2}{m}^2{\varphi }^2+{\mathcal{L}}^{(m)}\right],$$

$$e^{\tau }:={\left({}^{\varphi }{╱}_{{\varphi }_{{}_0}}\right)}^3,$$

$${\kappa }_{\text{eff}}(\tau )\equiv \sqrt{8\pi G_{\text{eff}}(\tau )}=\kappa e^{{}^{-\tau }{╱}_3}.$$

$$\mathcal{S}=\int d^{4}x\sqrt{-g}\left[e^{{}^{2\tau }{╱}_3}\left\{\frac{\mathcal{R}}{2{\kappa }^2}-\frac{\epsilon }{2}{g}^{\mu \nu }{\partial }_{\mu }\tau {\partial }_{\nu }\tau -\Lambda \right\}+{\mathcal{L}}^{(m)}\right],$$

$$\epsilon =\frac{\left(1-6\beta \right){\varphi }_{{}_0}^2}{9}=\frac{1-6\beta }{9{\kappa }^2\beta }\text{mbox}and\Lambda =\frac{1}{2}{m}^2{\varphi }_{{}_0}^2=\frac{m^2}{2{\kappa }^2\beta }.$$

$$\left|\mathcal{T}\right|:=\sqrt{-g^{\mu \nu }{\mathcal{T}}_{\mu }{\mathcal{T}}_{\nu }}=\sqrt{-g^{\mu \nu }{\partial }_{\mu }\tau {\partial }_{\nu }\tau }\text{mbox}and\left|\mathcal{A}\right|:=\sqrt{-g^{\mu \nu }{\mathcal{A}}_{\mu }{\mathcal{A}}_{\nu }}=4\kappa \sqrt{3\Lambda }.$$

$$\hat{\mathcal{S}}=\int d^{4}x\sqrt{-\hat{g}}\left[\frac{\hat{\mathcal{R}}}{2{\kappa }^2}-\frac{{\zeta }^2}{2}{\hat{g}}^{\mu \nu }{\partial }_{\mu }\tau {\partial }_{\nu }\tau -\Lambda e^{{}^{-2\tau }{╱}_3}+{\hat{\mathcal{L}}}^{(m)}\left(\hat{g},\tau \right)\right],$$

$$e^{{}^{\tau }{╱}_3}=\frac{1}{3}\int dt\left|\mathcal{T}\right|\text{mbox}and\Lambda =\frac{{\left|\mathcal{A}\right|}^2}{48{\kappa }^2},$$

$$H^2=\frac{{\kappa }^2}{3}\left[{\rho }^{(m)}+{\rho }^{(\tau )}\right]\text{mbox}and\dot{H}=-\frac{{\kappa }^2}{2}\left[{\rho }^{(m)}+{\rho }^{(\tau )}+p^{(\tau )}\right],$$

$${\rho }^{(\tau )}=\frac{{\zeta }^2}{2}{\dot{\tau }}^2+\Lambda e^{{}^{-2\tau }{╱}_3}\text{mbox}and{p}^{(\tau )}=\frac{{\zeta }^2}{2}{\dot{\tau }}^2-\Lambda e^{{}^{-2\tau }{╱}_3}.$$

$${\dot{\rho }}^{(m)}+3H{\rho }^{(m)}=-{\rho }^{(m)}\frac{\dot{\tau }}{3}\text{mbox}and{\dot{\rho }}^{(\tau )}+3H\left({\rho }^{(\tau )}+p^{(\tau )}\right)={\rho }^{(m)}\frac{\dot{\tau }}{3},$$

$${\rho }^{(m)}(t)=\frac{{\rho }_{{}_0}^{(m)}}{a^3(t)}{e}^{{}^{-\tau (t)}{╱}_3},$$

$$\rho :=\frac{3H^2}{{\kappa }^2}={\rho }_{\text{eff}}^{(m)}+{\rho }_{{}_X},$$

$${\rho }_{\text{eff}}^{(m)}=\frac{{\rho }_{{}_0}^{(m)}}{a^3}={\rho }^{(m)}{e}^{{}^{\tau }{╱}_3}\text{mbox}and{\rho }_{{}_X}={\rho }^{(\tau )}+\left(e^{{}^{-\tau }{╱}_3}-1\right){\rho }_{\text{eff}}^{(m)}.$$

$${\Omega }^{(m)}+{\Omega }^{(\tau )}={\Omega }_{\text{eff}}^{(m)}+{\Omega }_{\text{ssmallX}}=1,$$

$${\dot{\rho }}_{{}_X}+3H\left({\rho }_{{}_X}+p_{{}_X}\right)=0.$$

$$X:=\frac{\dot{\tau }}{3\sqrt{6}H}\text{mbox}andY:=\frac{\kappa \sqrt{\Lambda }{e}^{{}^{-\tau }{╱}_3}}{\sqrt{3}H},$$

$$\frac{dX}{dN}=\frac{3}{2}\left(X-\beta \sqrt{\frac{2}{3}}\right)\left(\frac{X^2}{\beta }-Y^2-1\right),$$

$$\frac{dY}{dN}=\frac{3Y}{2}\left(\frac{X^2}{\beta }-Y^2-2\sqrt{\frac{2}{3}}X+1\right),$$

$$\frac{X^2}{\beta }+Y^2-1+{\Omega }^{(m)}=0.$$

$$\tau (N)=3\sqrt{6}F(N),\text{mbox}whereF(N)\equiv {\int }_0^{N}X(\mathcal{N})d\mathcal{N}.$$

$${\Omega }_{\text{eff}}^{(m)}:=\frac{{\rho }_{\text{eff}}^{(m)}}{\rho }=\left(1-\frac{X^2}{\beta }-Y^2\right)e^{\sqrt{6}F}.$$

$$\mathsf{w}:=\frac{p}{\rho }=\frac{X^2}{\beta }-Y^2,$$

$${\mathsf{w}}_{{}_X}:=\frac{p_{{}_X}}{{\rho }_{{}_X}}=\frac{X^2-\beta Y^2}{\beta {\Omega }_{\text{ssmallX}}},$$

$$\frac{d}{dN}\begin{pmatrix}\delta X\\ \delta Y\end{pmatrix}=\,\mathcal{M}\begin{pmatrix}\delta X\\ \delta Y\end{pmatrix}\,,\qquad\mbox{where}\qquad\mathcal{M}\,\equiv\left[\begin{array}[]{cc}{}^{\displaystyle\partial\mathcal{F}}\!\!\!\diagup_{\!\!\displaystyle\partial X}&{}^{\displaystyle\partial\mathcal{F}}\!\!\!\diagup_{\!\!\displaystyle\partial Y}\\ {}^{\displaystyle\partial\mathcal{G}}\!\!\!\diagup_{\!\!\displaystyle\partial X}&{}^{\displaystyle\partial\mathcal{G}}\!\!\!\diagup_{\!\!\displaystyle\partial Y}\end{array}\right]_{\left(X_{c},Y_{c}\right)}\!\!\!\!.$$

$$\frac{X_c^2}{\beta }+Y_c^2=1.$$

$$H^2=\frac{{\kappa }^2\Lambda }{3-s}{a}^{-2s},\text{mbox}underthesubstitution:s=2\beta ,$$

$$X(0)=\sqrt{\frac{2}{3}}\beta =0.0082\text{mbox}andY(0)=\sqrt{1-\frac{2\beta }{3}-{\Omega }_{{}_0}^{(m)}}=0.8327.$$

Full text

1 Introduction

Dynamical stability is a major requirement for cosmological solutions representing dark energy (DE) that supposedly drives the late-time cosmic acceleration [1]. While the question as to how the DE evolves has been contemplated by a plethora of theoretical surmises and conjectures [2, 3, 4], observations have mostly been in favour of a non-dynamical DE, reminiscent of a cosmological constant $\Lambda$, at low to moderately high redshifts [5, 6, 7]. However, some scope is there to look for (albeit mild) deviations from the concordant $\Lambda$CDM model, comprising of $\Lambda$ and cold dark matter (CDM) as the dominant constituents of the universe [8, 9, 10]. In fact, the dynamical aspects of the DE are always worth examining, for a sufficiently longer span of evolution, tracing back from deep in the past, till extrapolating to high blueshifts in the future [3, 4, 11]. The theoretical motivation for this is obvious, in view of the well-known fine tuning and coincidence problems affecting the $\Lambda$CDM cosmology [8, 9].

Extensive searches for the dynamical DE, within the standard Friedmann-Robertson-Walker (FRW) framework, have mostly accounted for the scalar field candidates, such as quintessence, k-essence, tachyon, dilaton, chameleon, etc. [12, 13, 14, 15, 16], which have had many intriguing features [2, 3]. However, in recent years the focus has shifted to a purely geometric characterization of the DE in the so-called modified gravity theories [17] of e.g. the $f(\mathcal{R})$ type [18], where $\mathcal{R}$ is the Riemannian curvature scalar. Such theories can also be mapped to scalar-tensor theories [19, 20, 21, 22], and hence give rise to interacting (or unified) dark energy–matter scenarios [23] under conformal transformations. One’s perception though, of a ‘geometrical’ DE, is not limited to the formulations in the Riemannian space-time only. We may equally well look into the cosmologies emerging from the rather conventional extensions of General Relativity (GR), such as that formulated in the four-dimensional Riemann-Cartan ($U_{4}$) space-time with torsion — an antisymmetric tensor field that generalizes the Levi-Civita connections in GR [24, 25, 26, 27, 28, 29, 30]. Torsion is often considered as a geometric entity that provides a classical background for quantized spinning matter, and is therefore an inherent part of a fundamental (quantum gravitational) theory, such as string theory [28, 31, 32]. A completely antisymmetric torsion can have its source in the closed string massless Kalb-Ramond mode [33, 34], with interesting implications in cosmology and astrophysics [35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. Among other torsion scenarios of interest in the cosmological context, most notable are those based on the teleparallel $f(T)$ theories [45], extended gravity theories [46], Poincaré gauge theory of gravity [47, 48, 49], etc.

We in this paper turn our attention to the formalism of a metric-scalar-torsion (MST) theory developed in an earlier work (henceforth ‘paper I’) [50]. Such a theory deals with the $U_{4}$ Lagrangian non-minimally coupled to a scalar field $\phi$ (of presumably primordial origin), in a way that no uniqueness problem arises [28, 50, 51]. Now, in the standard cosmological framework, the torsion degrees of freedom get restricted by the FRW metric structure [30]. Also since $\phi$ acts as the source of the trace mode of torsion (via the corresponding equation of motion), we effectively have a scalar-tensor equivalent MST setup. The pseudo-trace mode of torsion can give rise to a mass term for $\phi$, via suitable augmentation of the effective action with say, some higher order torsion terms [50]. Considering further a pressureless dust-like cosmological matter, viable DE solutions in Einstein and Jordan frames have been worked out analytically in paper I [50], keeping the cosmological parameters within the corresponding error estimates for the $\Lambda$CDM model from recent observations. However, there remains the important question:

“would these (and possibly a few other) solutions, persist over time (i.e. stable), once subjected to fluctuations in the solution space (or the phase space)?”

Answering this requires an in-depth analysis of the MST-cosmological dynamics in both Einstein and Jordan frames. Our objective in this paper is to carry out such an analysis, by constructing (from cosmological equations in either cases) the autonomous system of equations in terms of suitable phase space variables.

For simplicity, we take into account only the two dominant components of the universe, viz. the dust and the scalar field $\phi$, whence the phase space is a two-dimensional (i.e. a phase plane). However, instead of working with $\phi$ and its mass $m$, throughout the analysis we resort to a torsion scalar $\tau\left(\sim 3\ln\phi\right)$ and a torsion constant $\Lambda\left(\sim m^{2}\right)$, so as to have a clear understanding of torsion’s effect on the dynamics. Whereas $\tau$ is equal to the time-integral of the norm of the torsion trace vector (as defined in the original Jordan frame), $\Lambda$ is given by the norm of the pseudo-trace vector of torsion (modulo some numerical factor).

We follow the standard methodology based on linear perturbation theory [2, 11] to determine the critical points (CPs) in the phase plane and their characteristic type and nature. This is essential, since each CP represents an equilibrium state of the system (the universe) in the asymptotic limit ($N=\ln a\rightarrow\infty$, where $a$ is the cosmological scale factor). Now, a given cosmological solution is considered stable if it transpires to the dynamical evolution of the universe at a stable CP. However, as is common in a plethora of contexts in the literature [2, 11, 52, 53, 54], there are instances of more than one CPs existing in a certain parametric domain, i.e. the range of values of a system parameter (e.g. a coupling parameter). The same is the situation we find here, in both Einstein and Jordan frames, for certain domains of our MST-coupling parameter $\beta$. This compels us to analyse the dynamical evolution at each individual CP and figure out the appropriate one(s) in the respective (Einstein or Jordan) frame. Additionally, we have to determine the parametric domain(s) in which a stable CP supports solution(s) that exhibit cosmic acceleration in the asymptotic limit. We do so in both Einstein and Jordan frames, and hence show that the corresponding exact solutions found analytically in paper are indeed stable. Numerically solving the autonomous equations, subject to appropriate initial conditions, we work out a host of relevant trajectories in the phase plane, in order to examine the overall dynamics of the system and its constituents leading up to the CPs. Also, for certain fiducial settings, we demonstrate the evolution of cosmological parameters of interest, such as the effective DE density and equation of state (EoS) parameters, $\Omega_{\text{\ssmall{\it X}}}$ and $\mathsf{w}_{\!{}_{X}}$, along the corresponding trajectories.

We carry out the dynamical analysis first in the Einstein frame, in which the cosmological equations are rather simple and have resemblance with those for quintessence. There are however two major differences. Firstly, the torsion scalar $\tau$ interacts with the (apriori dust-like) fluid, thus affecting the dynamics of both. As such, it is not possible to make a direct comparison of the MST-cosmological parameters with those estimated (from observations) for known models, such as $\Lambda$CDM. It is rather convenient to resort to the scenario in which the critical density of the universe, $\rho$, is decomposed into two effective non-interacting components, viz. the dust-like matter and a left-over, supposedly the DE [50]. In such a scenario, the physical relevance of the existent CPs in the phase plane is implicated by the eventual extinction of the matter sector, irrespective of the initial conditions. Secondly, it is not desirable to have the chosen phase space variables depending explicitly on the system parameter $\beta$. Otherwise their calibration would keep on changing with $\beta$, different domains of which are assigned for the existence of the CPs and(or) their physical relevance. Therefore, if instead of $\tau$ we choose to work with a redefined (quintessence-like) field $\varphi$ which absorbs $\beta$ in it, the physically admissible regions $\mathfrak{M}$ for the trajectories in the phase plane would get altered in shape and size [55].

Repeating the analysis in the Jordan frame is straightforward, but cumbersome because of an explicit $\tau$-dependence of the gravitational coupling factor $\kappa_{\text{\scriptsize eff}}^{2}\equiv 8\pi G_{\text{\scriptsize eff}}$, where $G_{\text{\scriptsize eff}}$ is the generalization of the Newton’s constant. There are some interesting consequences of this though, culminating from two legitimate standpoints. In principle, we may resort to one of the following: (a) a conventional scenario in which the critical density $\rho_{\!{}_{J}}$ varies with $\kappa_{\text{\scriptsize eff}}$ and is not conserved, although the matter density is conserved [11], and (b) an effective scenario in which the critical density $\rho$, defined as in a minimally coupled theory, is the sum of the densities of the dust and a left-over (supposedly the DE), which are individually conserved [50]. Now, since it is the same Jordan frame MST setup looked from different perspectives, the general outcomes of the dynamical analysis remain the same in both the scenarios, viz. the same number of CPs of the same type and nature. One difference is there though — in the conventional scenario, the physically admissible regions $\mathfrak{M}$ are confined within two similar curves (conic sections) in the phase plane, whereas only the outer curves are there in the effective scenario, for the same values of the effective Brans-Dicke parameter $\mathfrak{w}$ ($\propto\beta^{-1}$)111We choose to take $\mathfrak{w}$ as the Jordan frame system parameter, different domains of which ascribe to the existent and(or) the physical relevant CPs. The limiting values set on $\mathfrak{w}$ from extensive studies [56, 57], provide an independent credibility check of the cosmological solutions [50].. The cosmological dynamics in the effective scenario, although by-and-large similar to that in the Einstein frame, has an intriguing feature, viz. the existence of a stable CP that supports solutions which not only exhibit cosmic acceleration in the asymptotic limit, but also a super-accelerating or phantom regime in course of their evolution. One such stable solution is actually that found in paper I, for which the phantom barrier crossing takes place at an epoch in the near past, whereafter the phantom regime continues eternally [50].

This paper is organized as follows: starting with a general description of the MST formalism in §2, we write down the scalar-tensor equivalent actions in both Jordan and Einstein frames, in terms of the torsion scalar $\tau$. Considering first the Einstein frame MST-cosmological setup in §3, we proceed as follows: (i) work out (in §3.1) the equations of motion for the effectively non-interacting dust and DE sectors, (ii) construct from them (in §3.2) the autonomous system of equations and examine the domains of existence and(or) physical relevance of the CPs and also their type and nature, (iii) study (in §3.3) the dynamical evolution of the universe at each of the CPs and examine the stability of viable DE solutions such as the one found in paper I, and finally (iv) obtain the phase plane trajectories (in §3.4) by numerically solving the autonomous equations for appropriate sets of initial conditions, and hence illustrate the evolution of the DE density and EoS parameters along a fiducial trajectory that leads up to a stable CP. Almost the same chronology is maintained while repeating the dynamical analysis in §4 for the MST-cosmological setup in the Jordan frame. Some characteristic differences are there of course, compared to the Einstein frame analysis. Accordingly, a detailed account of, for e.g. the shape and size of the physically admissible region(s) $\mathfrak{M}$ with the variation in the system parameter $\mathfrak{w}$, is given. Also the phantom barrier crossing in the effective Jordan frame scenario is illustrated clearly, by working out the evolution along trajectories for fiducial settings corresponding to two different values of $\mathfrak{w}$. We conclude in §5 with a summary and a discussion on some implications and possible extensions.

We use the same notations and conventions as in paper I, viz. the metric signature throughout is $\left(-,+,+,+\right)$, units are chosen so that the speed of light $c=1$, and the determinant of the metric tensor $g_{\mu\nu}$ is denoted by $g$.

2 The general MST formalism in the cosmological setup

Let us discuss the basic tenets of the metric-scalar-torsion (MST) formalism, viz. that of scalar field couplings to the four-dimensional Riemann-Cartan ($U_{4}$) Lagrangian [50]. The $U_{4}$ space-time is characterized by an asymmetric affine connection: $\widetilde{\Gamma}^{\lambda}_{~{}~{}\mu\nu}(\neq\widetilde{\Gamma}^{\lambda}_{~{}~{}\nu\mu})$, which incorporates the third-rank torsion tensor defined as: $T^{\lambda}_{~{}~{}\mu\nu}:=\widetilde{\Gamma}^{\lambda}_{~{}\mu\nu}-\widetilde{\Gamma}^{\lambda}_{~{}\nu\mu}\left(=-\,T^{\lambda}_{~{}~{}\nu\mu}\right)$. Essentially, the Riemannian ($R_{4}$) covariant derivative $\nabla\!_{\mu}$ in GR (defined via the symmetric Levi-Cevita connections $\Gamma^{\lambda}_{~{}\mu\nu}$) is replaced with that ($\widetilde{\nabla}\!_{\mu}$) defined via $\widetilde{\Gamma}^{\lambda}_{~{}\mu\nu}$, preserving the metricity condition $\widetilde{\nabla}\!_{\mu}g_{\alpha\beta}=0$.

The torsion tensor can be decomposed into three irreducible modes, viz. the trace $\,\mathcal{T}_{\mu}:=T^{\nu}_{~{}\mu\nu}$, the pseudo-trace $\,\mathcal{A}^{\sigma}:=\epsilon^{\alpha\beta\gamma\sigma}\,T_{\alpha\beta\gamma}$ and the (pseudo-)tracefree tensorial residue $\,\mathcal{Q}_{\mu\nu\sigma}$, whence the $U_{4}$ analogue of the Ricci scalar curvature $\mathcal{R}$ is

[TABLE]

Accordingly, the free $U_{4}$ Lagrangian $\,L_{U_{4}}=\sqrt{-g}\,\widetilde{\mathcal{R}}$ has a purely algebraic dependence on torsion222Note that the $\sqrt{-g}\,\nabla\!_{\mu}\,\mathcal{T}^{\mu}$ term in $L_{U_{4}}$ is merely a total divergence (or, a boundary term).. Now, while coupling a scalar field $\phi$ to $L_{U_{4}}$, one encounters the well-known problem of non-uniqueness of the resulting action under the minimal coupling scheme ($\partial_{\mu}\rightarrow\widetilde{\nabla}\!_{\mu}$) [28, 58]. A simple (and convenient) way to avoid this is to assume a non-minimal term $\phi^{2}\,L_{U_{4}}$, so that upto a total divergence the action is [50]

[TABLE]

where $\beta$ is a dimensionless coupling constant, $V(\phi)$ is the scalar field potential, and $\mathcal{L}^{(m)}$ is the Lagrangian density for other matter fields in the theory. Eq. (2), dubbed as the ‘MST action’ [50], leads to the equation of motion $\,\mathcal{T}_{\mu}=\left({}^{\displaystyle 3\,}\!\!\!\diagup_{\!\!\displaystyle\!\phi}\right)\partial_{\mu}\phi$, which implies that the scalar field $\phi$ acts a source of the trace mode $\mathcal{T}_{\mu}$ of torsion. Moreover, in order to preserve the FRW metric structure in a standard cosmological setup, one requires the tensor mode $\mathcal{Q}_{\mu\nu\sigma}$ of torsion to vanish altogether, and the vector modes $\mathcal{T}_{\mu}$ and $\mathcal{A}^{\mu}$ to have only their temporal components existent [30, 59]. Also since the torsion field is generally taken to be massless [24, 32], one expects the scalar field source $\phi$, of its trace mode $\mathcal{T}_{\mu}$ to be massless as well. However, the pseudo-trace mode of torsion, $\mathcal{A}^{\mu}$, may effectively lead to a scalar field potential $V(\phi)={}^{\displaystyle m^{2}\phi^{2}}\!\!\!\diagup_{\!\!\displaystyle 2}$, where $m$ is a mass parameter for $\phi$. Such a possibility arises from a suitable augmentation of the MST action (2) with say, some higher order torsion term(s) such as $\,\phi^{2}\!\left(\mathcal{A}_{\mu}\mathcal{A}^{\mu}\right)^{2}$, whence one gets $\,\mathcal{A}_{\mu}\mathcal{A}^{\mu}=\,-\,{}^{\displaystyle 48\,m^{2}\!}\!\!\!\diagup_{\!\!\displaystyle\beta}$ via the corresponding equation of motion333Note that, a mass term for $\phi$ may also result from a norm-fixing constraint on $\mathcal{A}^{\mu}$ [50], similar to that in the vector-tensor gravity theories of Einstein-æther type [60, 61, 62], or in the mimetic gravity theories [63, 64, 65]. However, such an analogy has no specific physical motivation. [50]. The effective MST action then assumes the form of a scalar-tensor action in the Jordan frame [50]:

[TABLE]

where $\kappa=\sqrt{8\pi G}$ and $\phi_{{}_{0}}=\left(\kappa\sqrt{\beta}\right)^{-1}$ is the value444The parameter $\beta$ is of course taken to be positive definite, as otherwise the underlying quantum gravitational theory would be unbounded from below. of $\phi$ at the present epoch $t_{{}_{0}}$, such that the running gravitational coupling parameter $\,G_{\text{\scriptsize eff}}\sim\phi^{-2}$ has its present-day value $\,G_{\text{\scriptsize eff}}(t_{{}_{0}})=G$, the Newton’s constant.

Let us now define a dimensionless scalar field $\tau$ as:

[TABLE]

so that at $t=t_{{}_{0}}$, $\tau=0$, and the effective gravitational coupling factor is given by

[TABLE]

The Jordan frame action (2.3) is then expressed as

[TABLE]

where $\varepsilon$ and $\Lambda$ are two dimensionful constants, given by

[TABLE]

One may note that the field $\tau$ and the constant $\Lambda$ facilitates a clear understanding of the roles of the individual torsion modes in the MST-cosmological dynamics. If one makes interpretations in terms of torsion parameters chosen as the norms of the trace and pseudo-trace vector modes of torsion [50], then it is easy to see that

[TABLE]

So, we always have $\Lambda\propto\left|\mathcal{A}\right|^{2}$, and simply $\,\tau=\int dt\left|\mathcal{T}\right|\,$ in a cosmological space-time described by the spatially flat FRW metric, viz. $\,g_{\mu\nu}=$ diag$[-1,\,a(t),\,a(t),\,a(t)]$, where $t$ is the comoving time and $a(t)$ is the scale factor. Henceforth, we shall appropriately refer to $\tau$ as the ‘torsion scalar’ and $\Lambda$ as the ‘torsion constant’.

Now, as is usual in a scalar-tensor equivalent theory, the equations of motion are much simpler in the Einstein frame than in the Jordan frame [20, 21]. The Einstein frame MST-action can be obtained from Eq. (2.6), under the conformal transformation $\,g_{\mu\nu}\longrightarrow\widehat{g}_{\mu\nu}=e^{{}^{\scriptstyle{2\tau\!}}\!\!\diagup_{\!\!\scriptstyle{3}}}\,g_{\mu\nu}$ [50]:

[TABLE]

where $\zeta=(3\kappa\sqrt{\beta})^{-1}$ is a dimensionful constant, $\widehat{g}\equiv\text{det}\left(\widehat{g}_{\mu\nu}\right)$ is the Einstein frame metric determinant, $\widehat{\mathcal{R}}$ is the corresponding curvature scalar, and $\,\widehat{\mathcal{L}}^{(m)}\!\left(\widehat{g},\tau\right)=e^{{}^{\scriptstyle{-4\tau\!}}\!\!\diagup_{\!\!\scriptstyle{3}}}\,\mathcal{L}^{(m)}\!\left(g(\widehat{g},\tau)\right)$ is the corresponding matter Lagrangian density. Despite their mathematical equivalence, the Einstein and Jordan frames in general have different outcomes of physical measurements. The reason is obviously the gravitational coupling factor, which varies in one frame and not in the other. In fact, there is a longstanding debate as to which of these frames is actually of physical relevance [20, 21]. For completeness therefore, we shall subsequently carry out the dynamical analysis for our MST-cosmological formalism in both the Einstein and Jordan frames, taking one or the other to be physically relevant. We shall set up first the corresponding (Einstein or Jordan frame) cosmological equations, for the two system constituents, viz. the torsion scalar $\tau$ and the cosmological matter in the form of a pressureless (non-relativistic) dust. Defining suitable variables for the corresponding two-dimensional phase space (or the phase plane), we shall then construct the autonomous system of equations and look for stable solutions representing an effective DE evolution. Since the dust feels the effect of torsion555That is, the dust-like fluid and the torsion scalar $\tau$ have a mutual interaction, which results from either the conformal transformation or the varying gravitational coupling factor $\kappa_{\text{\scriptsize eff}}(\tau)$., we shall resort to an effectively non-interacting picture in the respective (Einstein or Jordan) frame.

3 Phase plane analysis in the Einstein frame

Let us consider, in this section, the Einstein frame to be suitable for physical observations, and drop for brevity the hats over all quantities defined in this frame. We shall however continue with the expressions (2.8) for the norms $\left|\mathcal{T}\right|$ and $\left|\mathcal{A}\right|$, as defined in the Jordan frame, in order to keep track of the individual terms of our original MST-action (2). We have therefore the relationships:

[TABLE]

in the Einstein frame, with $t$ being the corresponding comoving time coordinate.

3.1 Cosmological equations and the effective scenario

The Friedmann and Raychaudhuri equations, obtained from the action (2.9), are

[TABLE]

where $\,H:={}^{\displaystyle\dot{a}}\!\!\!\diagup_{\!\!\displaystyle a}$ is the Hubble parameter corresponding to the Einstein frame scale factor $a(t)$ (the overhead dot $\{\cdot\}\equiv{}^{\displaystyle d}\!\!\!\diagup_{\!\!\displaystyle dt}$), $\,\rho^{(m)}$ is the energy density of the fluid matter, whereas $\rho^{(\tau)}$ and $p^{(\tau)}$ are respectively the field energy density and pressure:

[TABLE]

The corresponding energy-momentum conservation relation, given by

[TABLE]

imply that the cosmological fluid does not retain its ‘dust’ interpretation in the Einstein frame, as its energy density depends explicitly on the torsion scalar $\tau$:

[TABLE]

where $\,\rho^{(m)}_{{}_{0}}=\rho^{(m)}\rvert_{t=t_{{}_{0}}}=\rho^{(m)}\rvert_{a=1}$ is the present-day value of $\rho^{(m)}$.

Now, it is easy to see that under a field redefinition $\zeta\tau(t)\equiv\varphi(t)$, the above Eqs. (3.2)–(3.5) correspond to those in an interacting system of an apriori pressureless cosmological fluid and a quintessence scalar field $\varphi$ with an exponential potential [12]. Such a correspondence is however misleading, since the dynamical analysis of the system in terms of the torsion scalar $\tau$ has a radical difference with that in terms of the redefined field $\varphi$ (see the discussion in the next subsection). In other words, working with $\tau$ not only pinpoints the dynamical effects of torsion on the evolution of the universe, but also leads to results different from those of the dynamical analysis for the standard scalar-tensor cosmologies in the Einstein frame (which are essentially the systems of interacting quintessence and cosmological matter) [11]. Moreover, Eq. (3.5) suggests that the expression for the critical (or total) density of the universe, viz. $\rho={}^{\displaystyle 3H^{2}\!\!}\!\!\!\diagup_{\!\!\displaystyle\!\kappa^{2}}=\rho^{(m)}+\rho^{(\tau)}$, is barely of any use when it comes to making a comparison with the parametric estimations of well-known models, such as $\Lambda$CDM, from physical observations. It is rather convenient to express [50]

[TABLE]

where $\,\rho^{(m)}_{\text{\footnotesize eff}}$ is an effective (dust-like) matter density and $\,\rho_{\!{}_{X}}$ is a surplus density (which we consider to be due to the DE). These are given respectively as

[TABLE]

The Friedmann equation can then be recast as

[TABLE]

where $\,\Omega^{(m)}={}^{\displaystyle\rho^{(m)}\!}\!\!\!\diagup_{\!\!\displaystyle\rho}$ and $\,\Omega^{(m)}_{\text{\footnotesize eff}}={}^{\displaystyle\rho^{(m)}_{\text{\footnotesize eff}}\!}\!\!\!\diagup_{\!\!\displaystyle\rho}$ are the actual and the effective matter density parameters respectively, $\,\Omega^{(\tau)}={}^{\displaystyle\rho^{(\tau)}\!}\!\!\!\diagup_{\!\!\displaystyle\rho}$ is the density parameter for the field $\tau$, whereas $\,\Omega_{\text{\ssmall{\it X}}}={}^{\displaystyle\rho_{\!{}_{X}}}\!\!\!\diagup_{\!\!\displaystyle\rho}$ is that for the DE. Identifying the DE pressure as $\,p_{\!{}_{X}}=p^{(\tau)}$, we also have the DE conservation relation (obtained using Eqs. (3.4)):

[TABLE]

Eqs. (3.7)–(3.9) govern the dynamical evolution of the system. One may in principle look for their outright solutions, for e.g. by guessing suitable solution ansatze, and then examine the physical viability of those solutions [50]. A rather general alternative is to construct an autonomous system of first order coupled differential equations, out of Eqs. (3.7)–(3.9), and look for such a system the plausible real roots and their stability against small fluctuations in the solution space (or the phase space).

3.2 Autonomous equations and the critical points

Defining the phase space variables as

[TABLE]

we obtain from the above cosmological equations, the autonomous equations:

[TABLE]

where $\,N(t)\equiv\ln a(t)$ is the number of e-foldings, and one also has the constraint

[TABLE]

Inverting the definition of the variable $X$ in Eq. (3.10), we express

[TABLE]

Eqs. (3.6), (3.7) and (3.13) imply that the effective matter density parameter is

[TABLE]

Moreover, the total pressure being $\,p=p_{\!{}_{X}}=p^{(\tau)}$, we have from Eqs. (3.6), (3.7) and (3.10) the total equation of state (EoS) parameter of the system given by

[TABLE]

i.e. the EoS parameter for the DE is

[TABLE]

since $\,\mathsf{w}=\Omega_{\text{\ssmall{\it X}}}\mathsf{w}_{\!{}_{X}}$, where $\,\Omega_{\text{\ssmall{\it X}}}=1-\Omega^{(m)}_{\text{\footnotesize eff}}$, by Eq. (3.8).

It is worth noting here that the autonomous system of equations (3.11)–(3.13) is symmetric under the interchange $Y\rightarrow-Y$, which means that we can restrict our analysis to the region $Y\geq 0$ of the $XY$ phase plane without loss of generality. Moreover, as mentioned in the previous subsection, the MST-cosmological equations (3.2) and (3.3) correspond to those for quintessence, under the redefinition $\,\varphi(t)\equiv\zeta\tau(t)$. Such a correspondence may not in general be reflected in the dynamical analysis though, when the (dimensionless) phase space variables are defined using $\varphi$. Actually, in comparison to the standard (and even the interacting) quintessence scenarios, the MST setup has the intriguing aspect of the coupling parameter $\beta$ playing a potentially active role in determining the viable cosmologies. So it is imperative to allow for a discrete alteration of the value of $\beta$ in the dynamical analysis. However, for simplicity we may keep the other parameter in the theory, viz. $\Lambda$, to remain fixed. Now in such a situation, while working with $\varphi$ instead of $\tau$, we cannot use the above definitions (3.10) of the phase space variables $X$ and $Y$, with a mere substitution of $\tau$ by ${}^{\displaystyle\varphi}\!\!\!\diagup_{\!\!\displaystyle\zeta}$ therein. The reason is that $\zeta$ being proportional to $\beta^{-{}^{1}\!/_{2}}$, such a substitution would mean $X$ and $Y$ explicitly dependent on $\beta$. Therefore their calibration would change when the value of $\beta$ is changed, thus giving rise to an ambiguity in the analysis. So, in terms of $\varphi$ one has to define altogether different phase space variables, under the demand that they need to be free from any explicit $\beta$-dependence666Their implicit dependence on $\beta$ is not a worry though, as that would not affect their calibration. Hence, we would have a different set of autonomous equations which may lead to a different dynamics of the same system, if we resort to the redefined field $\varphi$, instead of persisting with the torsion scalar $\tau$.

Now, the objective of analysing autonomous equations, say ${}^{\displaystyle dX}\!\!\!\diagup_{\!\!\displaystyle dN}=\mathcal{F}(X,Y)$ and ${}^{\displaystyle dY}\!\!\!\diagup_{\!\!\displaystyle dN}=\mathcal{G}(X,Y)$, is to determine the critical points, or the equilibrium solutions, and consequently examine the type and nature of such solutions, i.e. their stability in the two-dimensional phase space formed by $X$ and $Y$. Here, $\mathcal{F}$ and $\mathcal{G}$ are given functions of $X$ and $Y$ for a particular system, for e.g. the right hand sides of Eqs. (3.11) and (3.12). By definition, a critical point (CP) is assigned coordinates $\left(X_{c},Y_{c}\right)$ at which $\mathcal{F}$ and $\mathcal{G}$ vanish. Now, for small changes $\left(\delta X,\delta Y\right)$ about $\left(X_{c},Y_{c}\right)$, we have the following eigenvalue equation in the linear perturbation theory [2]:

[TABLE]

The eigenvalues $\mu_{1}$ and $\mu_{2}$ of the matrix $\mathcal{M}$, i.e. their type (real, imaginary or complex) and their sign, determine whether the solutions for $\delta X$ and $\delta Y$ have modes that are exponentially growing or decaying with $N$, and as such assert the type and the nature of the CPs. Specifically, there are the following cases:

1. (i)

Real $\mu_{1}$ and $\mu_{2}$: (a) If they are both negative, then $\left(\delta X,\delta Y\right)\rightarrow 0$ as $N\rightarrow\infty$, which means that all phase space trajectories in the vicinity would terminate at the CP $\left(X_{c},Y_{c}\right)$, i.e. the latter acts as a stable nodal sink, or an attractor. (b) If they are both positive, then the perturbations $\left(\delta X,\delta Y\right)$ build up over time (or $N$), taking the trajectories away from the CP $\left(X_{c},Y_{c}\right)$, which therefore acts as an unstable nodal source. (c) If they are of opposite sign, then the CP $\left(X_{c},Y_{c}\right)$ is a saddle point, i.e. it acts as an attractor along a particular direction (viz. the attractor axis) and as an unstable point along the direction normal to that. In general, the trajectories may tend towards the saddle point, but eventually get repelled away. Only for the trajectories that start from somewhere on the attractor axis, proceed along that, and manage to reach the saddle point, the latter acts as an attractor. Otherwise, the saddle point is unstable in nature. 2. (ii)

Imaginary $\mu_{1}$ and $\mu_{2}$: The trajectories would in general describe an ellipse, having at its center the CP $\left(X_{c},Y_{c}\right)$, which is stable and is called a center. 3. (iii)

Complex $\mu_{1}$ and $\mu_{2}$: Depending on whether their real parts are negative or positive, the trajectories would spiral towards or away from the CP $\left(X_{c},Y_{c}\right)$, rendering the latter to be a stable spiral sink or an unstable spiral source.

In the case of vanishing $\mu_{1}$ and $\mu_{2}$, the type and nature of a CP are indeterministic if we persist with the linear perturbation theory, which actually breaks down [2].

For our autonomous set of equations (3.11)–(3.13), the above methodology implies the existence of five distinct CPs (discounting for the multiplicity of course). The coordinates $\left(X_{c},Y_{c}\right)$ of these CPs and the domains of their existence and physical relevance (given by appropriate range of values of the parameter $\beta$) are listed in Table 3.2.

The eigenvalues $\mu_{1}$ and $\mu_{2}$ of the linear perturbation matrix $\mathcal{M}$ at each CP, and also the type and nature of the CPs are shown in Table 3.2.

The criterion for their physical relevance actually follows from Eq. (3.15) for the effective matter density parameter $\Omega^{(m)}_{\text{\footnotesize eff}}$ and the constraint (3.13). The presence of the exponential factor $e^{\sqrt{6}\,F}$ in Eq. (3.15), where $\,F(N)\equiv\int_{0}^{N}X(\mathcal{N})\,d\mathcal{N}$, implies that $\Omega^{(m)}_{\text{\footnotesize eff}}$ may keep on evolving with time (or $N$) even after the system reaches a CP. So the condition for a physically realistic matter density, viz. $\Omega^{(m)}_{\text{\footnotesize eff}}<1$, may get violated at some epoch, leaving the corresponding cosmological model unphysical at that CP. The exception(s) though is(are) the scenario(s) in which $\Omega^{(m)}_{\text{\footnotesize eff}}=0$ in the asymptotic limit. From Eq. (3.13) we find the corresponding (physically relevant) CP(s) to be on the circumference of an ellipse, whose center is at the origin of the phase plane:

[TABLE]

The region enclosed by this ellipse (let us call it $\mathfrak{M}$) is therefore the physically admissible region in the phase plane. The identification of this region should be emphasized as the key difference between the dynamical system analysis here for the Einstein frame MST-cosmology, and that for the standard (quintessence-type) scalar field models in which the energy densities due to the field and the cosmological fluid matter are individually conserved. See below the comparison in the next subsection.

3.3 Dynamical evolution of the universe at each critical point

For viable cosmologies, particularly from the perspective of the late-time cosmic acceleration (or of the DE), it is necessary to identify the supportive CP(s). As the criterion (3.19) of the physical relevance of the CPs follows from the argument that $\Omega^{(m)}_{\text{\footnotesize eff}}=0$ (i.e. $\Omega_{\text{\ssmall{\it X}}}=1$) asymptotically, the cosmological solutions at each of the CPs are given entirely by the energy component with the density parameter $\Omega_{\text{\ssmall{\it X}}}$. Although we are choosing to call it the ’DE’, it is yet to be verified whether this component actually complies with an accelerated expansion of the universe at late times (preceded by a regime of a decelerated one). Let us therefore look into the characteristics of the CPs and the dynamical aspects of the universe at each CP:

CP $E_{1}$: Exists for all values of the parameter $\beta$ and is situated at the intersection of the ellipse (3.19) and the negative $X$-axis. It acts an unstable nodal source, since the eigenvalues $\mu_{1}$ and $\mu_{2}$ of the linear perturbation matrix $\mathcal{M}$ are real and positive. The equilibrium state of the solutions, implicated by this CP and given by the DE with density parameter $\Omega_{\text{\ssmall{\it X}}}=1$, is of an extremely decelerated expansion of the universe, because the EoS parameter of the DE is $\mathsf{w}_{\!{}_{X}}=1$, i.e. the DE behaves like a steep fluid.

CP $E_{2}$: Exists for all values of $\beta$ and is at the intersection of the ellipse (3.19) and the positive $X$-axis. It acts as an unstable nodal source when $\beta<{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ and as a stable nodal sink when $\beta>{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$. The solutions at this CP, in either case, are given by a steep fluid-like DE, exhibiting an extremely decelerated expansion of the universe. For $\beta={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ however, the type and nature of this CP are indeterministic, as both the eigenvalues $\mu_{1}$ and $\mu_{2}$ of $\mathcal{M}$ vanish.

CP $E_{3}$: Exists for all values of $\beta$ and lies on the $X$-axis, but inside the ellipse (3.19) except for $\beta={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$. It acts as an unstable saddle point (of no physical relevance though) for $\beta\neq{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, whereas for $\beta={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ it coincides with $E_{2}$, whence its type and nature are indeterministic. .

CP $E_{4}$: Exists whenever $\beta\leq{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ and is constrained to lie on the circumference of the ellipse (3.19). Whereas for $\beta={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ it coincides with $E_{2}$ and $E_{3}$ (i.e. its type and nature are indeterministic), for $\beta<{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ it acts as a stable nodal sink and leads to a DE dominated accelerated expansion of the universe if further $\beta<{}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$ (as can be checked easily777From Eq. (3.16) we see that at this CP $E_{4}$, viz. $\left(\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta,\sqrt{1-{}^{\scriptstyle{2\beta}}\!\!\diagup_{\!\!\scriptstyle{3}}}\right)$, the total EoS parameter for the system is $\mathsf{w}=-1+{}^{\scriptstyle{4\beta}}\!\!\diagup_{\!\!\scriptstyle{3}}$. So, the acceleration condition $\mathsf{w}<-{}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{3}}$ implies $\beta<{}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$.). This is in fact the only CP (among the five) that can support solution(s) in presence of a non-vanishing potential term for the torsion scalar $\tau$, and hence the viable DE model(s).

CP $E_{5}$: Exists whenever $\beta\leq{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, but is neither on the $X$-axis nor on the ellipse (3.19), unless for $\beta={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, whence it coincides with $E_{2},E_{3}$ and $E_{4}$ (i.e. its type and nature are indeterministic). For $\beta<{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ it acts as an unstable saddle point (which is physically irrelevant as well).

Let us now make a comparison with the (saddle and stable) CPs found in the dynamical analysis for the standard scenario of (non-interacting) dust and quintessence scalar field $\varphi$ with an exponential potential $\sim e^{\kappa\lambda\varphi}$, where $\lambda$ is some numerical factor [2, 11]. Among the saddle point(s) that could exist, depending on the value of $\lambda$, there is one always at the origin of the phase plane. This point supports solutions which require $\varphi$ to get obliterated asymptotically, leaving the dust as the sole constituent of the universe. Among the stable point(s) that could exist, there is one that supports solutions which, depending on $\lambda$, are given in the asymptotic limit either entirely by a non-dynamic DE component (i.e. a cosmological constant $\Lambda$) or by a $\Lambda$CDM configuration with the densities of $\Lambda$ and the dust of the same order of magnitude. In contrast, the analysis of our Einstein frame MST-cosmological dynamics leads to, depending on the value of the parameter $\beta$, the saddle points $E_{3}$ ($\forall\beta\neq{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$) and $E_{5}$ ($\forall\beta<{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$), and the stable point $E_{2}$ ($\forall\beta>{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$) or $E_{4}$ ($\forall\beta<{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$). However, we additionally have the criterion (3.19) that permits only those solutions for which the cosmological matter gets obliterated asymptotically. None of the saddle points ($E_{3}$ and $E_{5}$) satisfies this criterion though. Nevertheless, as demonstrated below, the CP $E_{3}$ that lies in the region $\mathfrak{M}$ enclosed by the ellipse (3.19), has significance in funnelling physical trajectories towards the stable point $E_{4}$. The latter is of course the only CP which supports solutions exhibiting the cosmic acceleration in the asymptotic limit, and that too for $\beta$ restricted up to a maximum value $\beta_{\text{max}}={}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$. So the DE models in the Einstein frame MST setup are viable only for a fixed parametric range $0<\beta<{}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$. This also suggests (from Table 3.2) that only four CPs are of practical importance, viz. two unstable points $E_{1}$ and $E_{2}$ at the intersections of the $X$-axis and the ellipse (3.19), a saddle point $E_{3}$ on the $X$-axis and inside this ellipse, and a stable point $E_{4}$ on the circumference of the ellipse.

As to the status of the exact solution we have found in paper I by explicitly solving the Einstein frame MST-cosmological equations [50], note that:

(i) The small parametric bound $\beta_{\text{max}}\sim 10^{-2}$, for the viability of an almost $\Lambda$CDM-like DE model described by such a solution [50], is compatible with the rather loose bound $\beta_{\text{max}}={}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$ we have obtained here from the dynamical analysis.

(ii) The universe described by such a solution must transpire to the dynamics of the DE and the cosmological matter at the stable point $E_{4}$. This can be seen by working out that at $E_{4}$, i.e. at $\left(X_{c}=\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta,\,Y_{c}=\sqrt{1-{}^{\scriptstyle{2\beta}}\!\!\diagup_{\!\!\scriptstyle{3}}}\right)$, the torsion scalar is given by $e^{\tau}=a^{6\beta}$, whence the expression for the Hubble parameter:

[TABLE]

is precisely the same as the asymptotic (i.e. the $a\rightarrow\infty$ limiting) form of that we have had in paper I, while deriving the exact solution in the Einstein frame (see section 4.1 therein). The stability of such a solution is thus established.

Now, for a clear understanding of the qualitative aspects of the evolution of the universe leading up to the stable point $E_{4}$, let us refer back for convenience, to the original decomposition of the critical density into the densities of the (interacting) torsion scalar $\tau$ and the cosmological fluid in the Einstein frame. Note the following:

I. Any point on the $Y$-axis of the phase plane implicates a non-dynamic $\tau$ (i.e. $\dot{\tau}=0$, since $X=0$), whose contribution to the total energy density of the universe, via the potential $\Lambda e^{{}^{\scriptstyle{-2\tau}}\!\!\diagup_{\!\!\scriptstyle{3}}}$, has a fixed value. The all-important torsion mode is therefore the pseudo-trace $\mathcal{A}^{\mu}$, which is assumed to give rise to the potential, that acts a cosmological constant $\Lambda$ (with the corresponding EoS parameter $\mathsf{w}_{\!{}_{X}}=-1$). Hence the overall configuration for a point on the $Y$-axis is that of $\Lambda$CDM. Moreover, the constraint (3.13) implies that the further such a point is from the origin, the greater is the contribution of torsion to the energy content of the universe (subject always to the condition ${}^{\displaystyle X^{2}\!\!}\!\!\!\diagup_{\!\!\displaystyle\beta}+Y^{2}\leq 1$, however). The $\Lambda$CDM trajectory hence shows the evolution along the $Y$-axis, upto the point $(0,1)$, i.e. the apex of the elliptic boundary (3.19) of the physically admissible region $\mathfrak{M}$ in the phase plane.

II. A point anywhere except on the $Y$-axis of the phase plane represents a system configuration in which the torsion scalar $\tau$ is dynamical. The extent of such dynamics is determined by the magnitude of $X$, or equivalently by the contribution of the trace mode $\mathcal{T}_{\mu}$ of torsion to the total energy content of the universe. However, the viability of a DE evolution depends on how dominant is the potential $\Lambda e^{{}^{\scriptstyle{-2\tau}}\!\!\diagup_{\!\!\scriptstyle{3}}}$, and hence the torsion pseudo-trace $\mathcal{A}^{\mu}$, over the dynamical mode $\mathcal{T}_{\mu}$. In other words, the latter has to be quite subdued, which is commensurable with the smallness of the parameter $\beta$. In fact, it is easy to see that for $\beta\ll 1$, the stable point $E_{4}$ located at the boundary (3.19) of $\mathfrak{M}$ supports cosmologies which can be summed up as small deviations from $\Lambda$CDM. A DE model of such sort has been the one studied earlier [50], in which statistical bounds on $\beta$ are found by demanding that the value of, say, the Hubble constant $H_{{}_{0}}$ has to be within the $1\sigma$ error limits of that predicted for $\Lambda$CDM from physical observations. Although these bounds are important from the observational perspective, in order to see the overall qualitative nature of stable cosmologies represented by $E_{4}$ it suffices one to resort to the general (model-independent) upper limits — the coarse one, viz. $\beta_{\text{max}}={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, for the stable point $E_{4}$ to exist in the first place, and the tighter one, viz. $\beta_{\text{max}}={}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$, in order that a phase of accelerated expansion of the universe is supported by $E_{4}$.

3.4 Numerical solutions of the autonomous equations

The limitation of the equilibrium solutions of the autonomous equations (or the critical points in the phase plane) is that they do not provide any quantitative information as to what the state of a system has been prior to reaching them. To overcome this (atleast partially), we require to find particular solutions of the autonomous equations. Doing so analytically is however a formidable proposition for the coupled set (3.11)–(3.13). We therefore resort to numerical techniques for a given range of initial values of the variables $X$ and $Y$. Each set of numerical solutions $X(N)$ and $Y(N)$ traces out a trajectory representing the system’s evolution in the phase plane, right from the point of origin (in accord with the initial conditions) till the termination at one of the CPs. This also enables us to see the variation (with $N$) of any explicit function of $X$ or(and) $Y$, over the lifetime of a phase plane trajectory. Hence, in a cosmology resulting from the chosen initial conditions, we can in principle plot the quantities of interest, such as $\Omega_{\text{\ssmall{\it X}}}$, $\mathsf{w}_{\!{}_{X}}$ and $\mathsf{w}$, over a significantly large span of time (or $N$).

Figs. 1(a) – (e) show the trajectories the system follows in the $XY$ plane to reach a CP for different initial values and for the parametric settings $\beta=0.01,0.05,0.1,0.5$ and $1$. As mentioned in the last subsection, there are actually four CPs to be taken into account, viz. those in the physically admissible region $\mathfrak{M}$, or on its elliptic boundary (3.19). These CPs are superimposed for comparison in each of the Figs. 1(a)–(e).

Our main interest however, from the point of view of the accelerating cosmologies, is in the Figs. 1(a), (b) and (c), corresponding to the settings $\beta=0.01,0.05$ and $0.1$ respectively. Figs. 1(d) and (e), corresponding to $\beta=0.5$ and $1$ respectively, are only for the sake of completeness in the illustration. We see that the trajectories originating somewhere in the elliptically bounded region $\mathfrak{M}$ of the phase plane, excluding the $X$-axis, terminate at the stable point $E_{4}$. On the other hand, the trajectories which originate at a point on the $X$-axis, tend to terminate at the saddle point $E_{3}$. In fact, the trajectories originating anywhere in $\mathfrak{M}$ except the $X$-axis are deflected in a direction parallel to the $X$-axis towards $E_{3}$, which deflects them in a direction vertically above it, i.e. towards $E_{4}$ situated at $\left(X=\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta,\,Y=\sqrt{1-{}^{\scriptstyle{2\beta}}\!\!\diagup_{\!\!\scriptstyle{3}}}\right)$. In other words, the saddle point $E_{3}$ has the effect of funnelling trajectories towards the stable point $E_{4}$. Of course, the latter exhibits its own attractive nature as well. The $X$-axis (i.e. $Y=0$) is the stable axis for $E_{3}$ and the line $X=\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta$ is its unstable axis. As $\beta$ is increased from a small value (say, $0.01$), the area of the region $\mathfrak{M}$ increases, with the decrease in the eccentricity of its elliptic boundary. Accordingly, the saddle point $E_{3}$ and the stable point $E_{4}$ shift away from the $Y$-axis, and so do the unstable points $E_{1}$ and $E_{2}$. The shift continues till $b\rightarrow{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, whence the CPs $E_{2},E_{3}$ and $E_{4}$ tend to coincide with the CP $E_{5}$ (not shown in the Figs. 1(a) – (e)), which approaches the elliptic boundary of the region $\mathfrak{M}$ from outside. Further increase in the value of $\beta$ (beyond ${}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$) would take $E_{1},E_{2}$ and $E_{3}$ outside the physical realm, whereas $E_{4}$ and $E_{5}$ would cease to exist.

Let us now examine the evolution of some cosmological parameters of interest, and the torsion parameters, along a fiducial trajectory corresponding to a particular setting, say $\beta=0.01$. As to the initial conditions for this fiducial trajectory, we may conveniently set them at the present epoch ($N=0$), i.e. by appropriately choosing the values of $X(0)$ and $Y(0)$. One obvious choice is that in line with the exact solution obtained in paper I, whose stability we have already established in this paper. Such a solution corresponds to taking the ansatz $\,e^{\tau}=a^{6\beta}$, whence it follows from Eq. (3.10) that $X=$ constant $=\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta$. Moreover, since $\tau(0)=0$ (or $F(0)=0$ by the definition (3.14)), we have from Eq. (3.15) $\,\Omega^{(m)}_{{}_{0}}\equiv\Omega^{(m)}_{\text{\footnotesize eff}}(0)=1-{}^{\displaystyle X^{2}(0)}\!\!\!\diagup_{\!\!\displaystyle\beta}-Y^{2}(0)$. Considering now the fiducial values:

(i) $\,\Omega^{(m)}_{{}_{0}}=0.3$ (which is roughly the observational prediction [5, 6, 7] for most of the model-independent and model-dependent DE parametrizations), and

(ii) $\,\beta=0.01$ (which is the order of magnitude estimation [50] for the above ansatz, using the WMAP and Planck results [5, 7]),

we have the initial conditions

[TABLE]

Fig. 2(a) shows the evolution of the density parameters $\Omega^{(m)}_{\text{\footnotesize eff}}(N)$ and $\Omega_{\text{\ssmall{\it X}}}(N)$ over the fiducial trajectory (with the above initial conditions) for a fairly wide range of e-foldings $N\in\{-5,5\}$888considering that $N=-5$ corresponds to a redshift $z=e^{-N}-1=147.4$.. As expected (in view of the small value of $\beta$), these parameters vary with $N$ similar to their $\Lambda$CDM analogues. There are some subtleties however. Note that the curves in Fig. 2(a) are asymmetric about $N=0$. Actually, as we go from the present ($N=0$) to the future regime ($N>0$), $\Omega^{(m)}_{\text{\footnotesize eff}}$ and $\Omega_{\text{\ssmall{\it X}}}$ rapidly approach a near-saturation to the values [math] and $1$ respectively. That is, the DE tends to dominate completely over the dust-like matter even in the not-so-distant future, which is quite identical to the case in $\Lambda$CDM. On the other hand, as we go back in the past ($N<0$), $\Omega^{(m)}_{\text{\footnotesize eff}}$ and $\Omega_{\text{\ssmall{\it X}}}$ first approach each other rapidly, reach an equality point, then diverge with the same rapidity, attain extremum values, and finally approach each other once again (albeit very slowly) further back in the past. This is of course dissimilar to what happens for $\Lambda$CDM, and its root cause can be traced to the original interaction between the torsion scalar $\tau$ and the cosmological fluid. Although this interaction gets obscured in the effective picture, it leaves its imprint on the density profiles. In fact, the dissimilarity with $\Lambda$CDM is also evident from the evolution patterns of the EoS parameters, viz. $\mathsf{w}_{\!{}_{X}}(N)$ and $\mathsf{w}(N)$ corresponding to the dynamical DE and the system respectively, over the fiducial trajectory. These are shown in Fig. 2(b), for the same range of $N$, viz. $\{-5,5\}$.

As to the evolution of the torsion scalar $\tau$, and that of the torsion parameters $|\mathcal{T}|$ and $|\mathcal{A}|$, over the fiducial trajectory, first note that since $X=\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta$ is constant throughout, Eq. (3.14) implies $F(N)=\sqrt{{}^{\scriptstyle{2}}\!\!\diagup_{\!\!\scriptstyle{3}}}\beta N$. Therefore, $\tau(N)=6\beta N$, i.e. $\tau$ varies linearly with $N$, as shown in Fig. 3(a). Using Eqs. (3.1) and (3.10) we can now work out the functional forms of the parameters $|\mathcal{T}|$ and $|\mathcal{A}|$, and plot them after appropriate dimensional scaling. We choose to scale $|\mathcal{T}|$ with the Hubble parameter $H$, in order to have a direct measure of the effect of the trace mode of torsion on the cosmological evolution. However, $|\mathcal{A}|$ being a constant, it is imperative to consider its ratio with $|\mathcal{T}|$ and see how that evolves with $N$, for the fiducial setting. The expressions of ${}^{\displaystyle\left|\mathcal{T}\right|}\!\!\!\diagup_{\!\!\displaystyle H}$ and ${}^{\displaystyle\left|\mathcal{A}\right|}\!\!\!\diagup_{\!\!\displaystyle\left|\mathcal{T}\right|}$ are found to be

[TABLE]

Fig. 3(a) shows almost a constant variation of ${}^{\displaystyle\left|\mathcal{T}\right|}\!\!\!\diagup_{\!\!\displaystyle H}$ with $N$, for the entire range of $N\in\{-5,5\}$. This is expected, since the exponential factor in Eq. (3.22) hardly affects the evolution for the small value of $\beta$ ($=0.01$). Moreover, the magnitude of ${}^{\displaystyle\left|\mathcal{T}\right|}\!\!\!\diagup_{\!\!\displaystyle H}$ is small as well ($\sim 6\beta=0.06$), implying that the torsion trace has a minor contribution to the Hubble rate throughout the evolution (in the above range of $N$). The evolution of the ratio ${}^{\displaystyle\left|\mathcal{A}\right|}\!\!\!\diagup_{\!\!\displaystyle\left|\mathcal{T}\right|}$ over the fiducial trajectory is shown in 3(b). Deep in the past, this ratio has been very small, meaning a complete domination of $\left|\mathcal{T}\right|$ over $\left|\mathcal{A}\right|$. However, as $\left|\mathcal{T}\right|$ decreases with $N$, the ratio grows and in the near past the growth becomes very rapid. At the present epoch its value is ${}^{\displaystyle 2Y(0)}\!\!\!\diagup_{\!\!\displaystyle\beta}=167.3$. It continues to increase in the future until saturating to a value of about $200$.

Overall thus, in the MST scenario we have a DE evolution slightly deviated from $\Lambda$CDM, when the norm of the pseudo-trace of torsion is about two orders of magnitude stronger than that of its trace, in the near past and future.

4 Phase plane analysis in the Jordan frame

Recall that while dealing with the MST formalism in a standard cosmological setup in section 2, we first arrived at the scalar-tensor action (2.6) in the Jordan frame with a characteristic non-minimal coupling between the torsion scalar $\tau$ and gravity. We then went on to obtain the corresponding Einstein frame action (2.9) and carried out the phase plane analysis in the Einstein frame, because it is straightforward and in line precisely with the canonical formulation of GR. Having performed the stability analysis of the Einstein frame solutions in section 3, we now turn our attention to their Jordan frame counterparts in this section.

It should be pointed out here that the Jordan frame action (2.6) can always be cast in the standard Brans-Dicke (BD) form [20, 21, 50], once we identify the BD scalar field to be $\kappa_{\text{\scriptsize eff}}^{-2}(\tau)=\kappa^{-2}e^{{}^{\scriptstyle{2\tau\!}}\!\!\diagup_{\!\!\scriptstyle{3}}}$ and the effective BD parameter to be

[TABLE]

where $\varepsilon$ is the dimensionful parameter [cf. Eq. (2.7)], that appears in the kinetic term of the torsion scalar $\tau$ in Eq. (2.6). Throughout the analysis in this section, we shall categorize the results in terms of $\mathfrak{w}$, instead of $\beta$ or $\varepsilon$. The advantage is that the analysis could be restricted by the lower bound on $\mathfrak{w}$ ($\simeq 40$, to as large as $40000$) estimated from a plethora of independent studies [56, 57]. Of course, we also have the theoretical limitation $\mathfrak{w}>-{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, in accord with the presumption $\beta>0$.

4.1 Cosmological equations and the effective scenario

For the spatially flat FRW metric $\,g_{\mu\nu}=$ diag$[-1,\,a(t),\,a(t),\,a(t)]$, where $t$ is the Jordan frame comoving time coordinate, we simply have $\tau=\int dt\left|\mathcal{T}\right|$. The corresponding action (2.6) leads to Friedmann and Raychaudhuri equations:

[TABLE]

where $\,H:={}^{\displaystyle\dot{a}}\!/_{\!\displaystyle a}$ is the Jordan frame Hubble parameter, and

[TABLE]

are the corresponding total (or critical) energy density and pressure, $\,\rho^{(m)}$ being the energy density of the cosmological matter in the form of the pressureless dust. Unlike the situation in the Einstein frame, the dust density $\rho^{(m)}$ is conserved here, implying $\,\rho^{(m)}=\rho^{(m)}_{{}_{0}}a^{-3}$, where $\,\rho^{(m)}_{{}_{0}}=\rho^{(m)}\rvert_{a=1}$ is the present-day value of $\rho^{(m)}$. That is, the ‘dust’ has its usual interpretation in the Jordan frame. However, the critical density $\rho_{\!{}_{J}}$ is itself satisfies a somewhat non-standard conservation relation:

[TABLE]

Therefore, while making predictions for the Jordan frame observables, one requires to take account of the rate at which the torsion scalar $\tau$ changes (or, the gravitational coupling varies999Since the coupling factor $\kappa_{\text{\scriptsize eff}}(\tau)=\kappa\,e^{{}^{\scriptstyle{-\tau\!}}\!\!\diagup_{\!\!\scriptstyle{3}}}$, we have $\,{}^{\displaystyle\dot{\tau}}\!\!\!\diagup_{\!\!\displaystyle 3}=-\,{}^{\displaystyle d}\!\!\!\diagup_{\!\!\displaystyle dt}\left[\ln\kappa_{\text{\scriptsize eff}}\right]$.) with time. This is reminiscent of the weirdness of all Brans-Dicke type of cosmological models, especially when it comes to making statistical estimates of the changes due to such models on the values of cosmological parameters predicted by some known model, for e.g. $\Lambda$CDM. One way to avoid this, i.e. bypass the direct confrontation with the running ($\tau$-dependent) gravitational coupling, is to define and decompose the critical density similar to that in a minimally coupled theory [50]:

[TABLE]

where $\,\rho^{(m)}=\rho^{(m)}_{{}_{0}}a^{-3}$ is the dust density, $\rho_{\!{}_{X}}$ is a surplus density, which we shall consider as the DE density. We demand $\rho$ to satisfy the usual conservation relation:

[TABLE]

where $p$ is the total effective pressure. Under such a demand, we have the effective DE density and pressure, given respectively as

[TABLE]

satisfying the conservation relation

[TABLE]

It should however be pointed out here that the phase plane analysis in the Jordan frame can in principle be carried out with either of the sets $\left\{\rho_{\!{}_{J}},p_{\!{}_{J}},\rho^{(m)}\right\}$ and $\left\{\rho,\rho^{(m)},\rho_{\!{}_{X}},p_{\!{}_{X}}\right\}$. The stability criterion for the system, determined from the autonomous equations (to be constructed), is anyway not commensurate with the choice of the system constituents. Since the analysis for the original description of the system in terms of the critical density $\rho_{\!{}_{J}}$ [cf. Eq. (4.3)] would, in many respects, be similar to that for the standard scalar-tensor cosmologies in the Jordan frame [11], we shall proceed to carry out that first. Then in due course, we shall examine the characteristic changes in the effective scenario, described by the Eqs. (4.6)–(4.9), most importantly with the objective of finding a stable DE model.

4.2 Autonomous equations and the physically admissible region(s)

Defining the phase space variables as:

[TABLE]

we obtain from Eqs. (4.2)–(4.5) the autonomous equations

[TABLE]

Moreover, the conventional forms (4.3) and (4.4) of the Jordan frame critical density $\rho_{\!{}_{J}}$ and the total pressure $p_{\!{}_{J}}$, lead to the following expressions for the corresponding matter density parameter $\Omega^{(m)}_{J}$ and the total EoS parameter $\mathsf{w}_{\!{}_{J}}$:

[TABLE]

Similar to what we have seen in the Einstein frame, the autonomous system is symmetric under $Y\rightarrow-Y$. All points in the $XY$ phase plane do not represent physically relevant solutions as well. However, as opposed to Eq. (3.15) for the Einstein frame matter density parameter $\Omega^{(m)}_{\text{\footnotesize eff}}$, the above Eq. (4.13) for the corresponding parameter $\Omega^{(m)}_{J}$ in the Jordan frame does not mean that that matter density would keep on evolving with time (or $N$) even after the system reaches a critical point (CP). Consequently, in order to comply with the physical relevance condition $0\leq\Omega^{(m)}_{J}\leq 1$ in the Jordan frame, a CP is not required to lie on an ellipse that corresponds to vanishing matter density, as in the Einstein frame. Instead, for all points (critical or otherwise) in the phase plane, Eq. (4.13) defines the physically admissible region(s) $\mathfrak{M}$, bounded from both inside and outside by the curves

[TABLE]

These curves evidently have identical shape, typified by conic sections of all possible sorts (viz. circle, ellipse, parabola and hyperbola)101010Note that, when the curves are open (i.e. parabolic or hyperbolic), they are actually the ‘demarcation lines’ in the phase plane, rather than ‘boundaries’. Moreover, for hyperbolic demarcations the sets of disjoint branches would define two disjoint regions $\mathfrak{M}$, instead of a solitary one otherwise., depending on the value of the parameter $\mathfrak{w}$, as examined below:

Case I. $\,\mathfrak{w}>0$ : Eq. (4.15) can be recast as

[TABLE]

So the curves corresponding to $C=0$ and $1$ in general represent concentric ellipses, an inner and an outer one respectively, which enclose a solitary region $\mathfrak{M}$ accommodating the physically relevant trajectories in the phase plane. This region is symmetric about about the point $(\sigma,0)$. Moreover, the inner ellipse passes through the origin, i.e. one of its antipodes (or vertices), along either the minor axis or the major axis, is coincident with the origin. Consider further the following sub-cases:

(i) $\,\mathfrak{w}>1$ : We have $\,b_{x}^{2}<b_{y}^{2}$. So, $b_{x}$ and $b_{y}$ are respectively the semi-minor and semi-major axes of the elliptic boundaries.

These ellipses have the same eccentricity $e=\sqrt{1-{}^{\displaystyle b_{x}^{2}}\!\!\!\diagup_{\!\!\displaystyle b_{y}^{2}}}=\sqrt{1-{}^{\displaystyle 1}\!\!\!\diagup_{\!\!\displaystyle\mathfrak{w}}}$, but differ in the focal length $\,f=e\,b_{y}$ and the semi-latus rectum $\,\ell=\left(1-e^{2}\right)b_{y}$. Moreover, one vertex along the minor axis of the inner ellipse coincides with the origin. Refer to Figs. 4(a), (b) and (c) for exemplary settings $\mathfrak{w}=250$, $25$ and $5$ respectively, and with the same $X$-range, $Y$-range and axes calibration. As $\mathfrak{w}$ increases (from the unit value), we see the following: {itemlist}

The ellipses becomes more and more eccentric and when $\mathfrak{w}\rightarrow\infty$, $\,e\rightarrow 1$.

Both the ellipses shift towards the origin along the $X$-axis, and their (common) center located at $(\sigma,0)$, coincides with the origin in the limit $\mathfrak{w}\rightarrow\infty$.

The semi-latus recta of the ellipses decrease, and for large $\mathfrak{w}$, $\,\ell\rvert_{{}_{C=0}}\sim\mathfrak{w}^{-{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}}$ whereas $\ell\rvert_{{}_{C=1}}\sim\mathfrak{w}^{-1}$, i.e. the reduction is faster for the inner ellipse compared to the outer one. The semi-major and semi-minor axes of the ellipses also decrease, and in the limit $\mathfrak{w}\rightarrow\infty$, $\,b_{y}\rvert_{{}_{C=1}}\rightarrow 1$ whereas $b_{y}\rvert_{{}_{C=0}}\rightarrow 0$. Hence, for a very large $\mathfrak{w}$ the physically admissible region $\mathfrak{M}$ nearly becomes a very eccentric elliptical zone about the origin, since its inner boundary gets almost obliterated (see Fig. 4(a) as an example). Some weirdness is there though in the limit $\mathfrak{w}\rightarrow\infty$, since $e\rightarrow 1$, but none of the elliptic boundaries reduce to parabola. Instead, the region $\mathfrak{M}$ tends to become a line segment of unit length, coincident with the $Y$-axis and symmetric about the origin. Such a line segment in the phase plane is reminiscent of the $\Lambda$CDM model, which one expects while taking the parametric limit $\beta\rightarrow 0$ (or $\mathfrak{w}\rightarrow\infty$) in the MST-cosmological setup. In the other limit $\mathfrak{w}\rightarrow 1$, there is nothing unusual — $e\rightarrow 0$ and the ellipses reduce to circles.

(ii) $\,\mathfrak{w}=1$: We have $\,b_{x}^{2}=b_{y}^{2}={}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}+C$, and therefore both the boundaries of $\mathfrak{M}$ are concentric circles, with the radii $\sqrt{{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}}$ and $\sqrt{{}^{\scriptstyle{5}}\!\!\diagup_{\!\!\scriptstyle{2}}}$ respectively for the inner one (passing through the origin) and the outer one. The common center is at \big{(}\sqrt{{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}},0\big{)}, and the (annular) region $\mathfrak{M}$ has width $\sqrt{{}^{\scriptstyle{5}}\!\!\diagup_{\!\!\scriptstyle{2}}}-\sqrt{{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}}=0.3564$. See Fig. 4(d).

(iii) $\,0<\mathfrak{w}<1$ : We have $\,b_{x}^{2}>b_{y}^{2}$. So, $b_{x}$ and $b_{y}$ are respectively the semi-major and semi-minor axes for either of the bounding ellipses, having the same eccentricity $e=\sqrt{1-{}^{\displaystyle b_{y}^{2}}\!\!\!\diagup_{\!\!\displaystyle b_{x}^{2}}}=\sqrt{1-\mathfrak{w}}$, but different focal lengths $\,f=e\,b_{x}$ and semi-latus recta $\,\ell=\left(1-e^{2}\right)b_{x}$. Compared to the sub-case (i) above, the minor and major axes of the ellipses are interchanged (see Fig. 4(e) for the exemplary setting $\mathfrak{w}=0.5$). As $\mathfrak{w}$ decreases (from the unit value), the eccentricity $e$ decreases, and the common center of the ellipses shifts more and more away from the origin along the $X$-axis. So the region $\mathfrak{M}$ becomes more and more narrow (wide) along the common minor (major) axis. In the limit $\mathfrak{w}\rightarrow 0$, $\mathfrak{M}$ becomes a parabolic zone, discussed below.

Case II. $\,\mathfrak{w}=0$ : Eq. (4.15) reduces to the equation of a parabola

[TABLE]

The inner parabolic demarcation line has its vertex coincident with the origin, whereas the outer one has vertex shifted from the origin by an amount $\,\sigma={}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{\!\!\sqrt{6}}}$ along the negative $X$-axis. Both the parabolae have the same focal length $f=\sqrt{{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{8}}}$, and the region $\mathfrak{M}$ between them is a solitary zone that extends (with gradually diminishing width) to infinity along the $X$-direction (see Fig. 4(f)). Such a region in the phase plane is in fact quite unique since it depicts the bizarre scenario in which there is no kinetic term for the torsion scalar $\tau$ in the Jordan frame action (2.6), and yet $\tau$ is dynamical by virtue of its non-minimal coupling with gravity.

Case III. $\,\mathfrak{w}<0$ : Eq. (4.15) can now be recast as

[TABLE]

provided $\,|\mathfrak{w}|<{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2C}}$, which of course applies here for ($C=0,1$) by the presumption $\beta>0$. The curves represented by Eq. (4.18) are hyperbolae with major axis along the $X$-direction and the (common) center, i.e. the point of intersection of the asymptotes, at $(-\sigma,0)$. Trajectories of physical relevance are therefore contained in two disjoint open regions $\mathfrak{M}$, between the pairs of hyperbolic branches (corresponding to $C=0,1$), which are symmetric about $(-\sigma,0)$. The region on the right has the vertex of the inner hyperbolic demarcation line located at the origin, whereas the one on the left is entirely in the negative $X$ quadrant.

The inner and outer hyperbolae have the same eccentricity $e=\sqrt{1+{}^{\displaystyle b_{y}^{2}}\!\!\!\diagup_{\!\!\displaystyle b_{x}^{2}}}=\sqrt{1+|\mathfrak{w}|}$, but different focal lengths $\,f=e\,b_{x}$ and semi-latus recta $\,\ell=\left(e^{2}-1\right)b_{x}$. However, unlike the case I above, change in the values of $\mathfrak{w}$ do not make any major alteration in the characteristics of the curves, for e.g. their eccentricity, focal length, etc. The only point to note is that the curves become rectangular hyperbolae for $\mathfrak{w}=-1$. Figs. 5(a), (b) and (c) illustrate these curves and the regions they demarcate, for exemplary settings $\mathfrak{w}=-0.5$, $\mathfrak{w}=-1$ and $\mathfrak{w}\approx-1.5$ respectively. With the increase in $|\mathfrak{w}|$, the hyperbolae become more and more eccentric and their branches come closer together, i.e. their focal lengths decrease. Also, the (common) center of the hyperbolae move towards the origin along the $X$-axis. As $\,\mathfrak{w}\rightarrow-{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, the two hyperbolic regions $\mathfrak{M}$ tend to merge together. In fact, $\,\mathfrak{w}\rightarrow-{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ is a discontinuous limit, since the axes of hyperbolae get interchanged for $\mathfrak{w}\gtrless-{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$ (not shown in Fig. 5). Anyway, we exclude here the possibility $\mathfrak{w}<-{}^{\scriptstyle{3}}\!\!\diagup_{\!\!\scriptstyle{2}}$, because it implies $\beta<0$.

4.3 Critical points and the dynamical evolution of the universe

Following the procedure described in subsection 3.2, we find that corresponding to the Jordan frame autonomous Eqs. (4.11) and (4.12) there could be five distinct critical points (CPs). These CPs $(X_{c},Y_{c})$ are listed in Table 4.3, alongwith the domains of their physical relevance, given by the appropriate range of values of $\mathfrak{w}$.

Table 4.3 shows the eigenvalues ($\mu_{1},\mu_{2}$) of the linear perturbation matrix $\mathcal{M}$ at each CP, and the type and nature of the CPs they determine (see subsection 3.2).

Similar to what have seen in the Einstein frame, at most two of the five CPs (viz. $J_{2}$ and $J_{4}$ here) could be stable points, whereas one CP (viz. $J_{3}$) could be a saddle point. Particularly intriguing are the CPs $J_{3}$ and $J_{4}$, just as their Einstein frame counterparts $E_{3}$ and $E_{4}$ respectively, in pointing out the difference of the MST-cosmological dynamics over that of the non-interacting quintessence and dust models [2, 11]. Whereas the asymptotic forms of the solutions represented by $J_{3}$ imply co-existence of the dust and the field $\tau$, those represented by $J_{4}$ imply cosmic acceleration111111Verify using Eq. (4.14) that at the CP $J_{4}$, the total EoS parameter of the system is $\mathsf{w}_{\!{}_{J}}=-1$, i.e. the acceleration condition $\mathsf{w}_{\!{}_{J}}<-{}^{\displaystyle 1}\!\!\!\diagup_{\!\!\displaystyle 3}$ is always satisfied. resulting from a $\tau$-dominance (the equivalent of the Einstein frame DE-dominance). Also the evolution of the universe leading up to $J_{4}$ has qualitative similarity with that leading up to the stable point $E_{4}$ in the Einstein frame.

There are some differences though in the interpretation of the results of the Jordan frame and Einstein frame analyses. Note in particular that for $\mathfrak{w}\in(-1,0)$, the CPs $J_{2}$ and $J_{4}$ are both stable points of physical relevance, whereas in the Einstein frame we never had two stable points in any given parametric domain. Moreover, unlike its Einstein frame analogue $E_{2}$, the stable point $J_{2}$ supports solutions exhibiting cosmic acceleration121212Verify using Eq. (4.14) that at the CP $J_{2}$, $\mathsf{w}_{\!{}_{J}}=1+\left({}^{\displaystyle 2}\!\!\!\diagup_{\!\!\displaystyle\mathfrak{w}}\right)\left[1+\sqrt{1+{}^{\displaystyle 2\mathfrak{w}}\!\!\!\diagup_{\!\!\displaystyle 3}}\right]$, so that in the domain $\mathfrak{w}\in(-1,0)$ the acceleration condition $\mathsf{w}_{\!{}_{J}}<-{}^{\displaystyle 1}\!\!\!\diagup_{\!\!\displaystyle 3}$ is satisfied. at the asymptotic limit (as do the stable point $J_{4}$). So for $-1<\mathfrak{w}<0$ (whence there are two disjoint hyperbolic regions $\mathfrak{M}$) we require to determine first which among $J_{2}$ and $J_{4}$ is the appropriate stable CP at which the physical trajectories would terminate. Whereas at $J_{4}$ the total EoS parameter of the system $\mathsf{w}_{\!{}_{J}}=-1$ always, it is easy to check that at $J_{2}$ one has $\mathsf{w}_{\!{}_{J}}<-1$ for $-1<\mathfrak{w}<0$. Therefore, in this domain the trajectories terminating at $J_{2}$ would imply that in the asymptotic limit the universe is at a super-accelerating or phantom state which is so strong that even $\mathsf{w}_{\!{}_{J}}<-1$. A super-acceleration with a component EoS parameter, such as that for $\tau$, less than $-1$ is acceptable. But if the system’s EoS parameter $\mathsf{w}_{\!{}_{J}}<-1$, it may lead to physical instabilities (against the cosmological metric perturbations [66, 67]). So there is a strong reason to discard $J_{2}$. Besides, $J_{2}$ being always shifted from the origin (along the $X$-axis), the trajectories which represent a nearly $\Lambda$CDM evolution cannot terminate at $J_{2}$. The reason is that the $\Lambda$CDM configuration (recovered in the limit $\mathfrak{w}\rightarrow\infty$) is represented by a line segment of finite size along the $Y$-axis and symmetric about the origin.

In fact, within the physically admissible zone $\mathfrak{M}$, any point on the $Y$-axis implicates a configuration of a non-dynamical $\tau$, whose potential (originating from the torsion pseudo-trace $\mathcal{A}^{\mu}$) plays the role of a cosmological constant $\Lambda$. In other words, for any point on the $Y$-axis, $\mathcal{A}^{\mu}$ is the key torsion constituent whose contribution to the critical density of the universe has a fixed value, reminiscent of that due to $\Lambda$. On the other hand, the points which are not located on the $Y$-axis represent the system configurations for a dynamical $\tau$. The extent of such dynamics is determined by the magnitude of $X$, or the torsion trace mode $\mathcal{T}_{\mu}$, because in the Jordan frame we have the torsion trace parameter $\left|\mathcal{T}\right|=\dot{\tau}\sim HX$.

As to solving numerically the autonomous Eqs. (4.11)–(4.12), since we have the admissible region $\mathfrak{M}$ bounded from both inside and outside now, the choice of initial conditions for the phase plane trajectories are much restricted compared to that in the Einstein frame. Also, for simplicity in the numerical calculations and to emphasize on the viable MST-cosmologies not much deviated from $\Lambda$CDM, we resort to the situations corresponding to $\mathfrak{w}>0$ only, whence $\mathfrak{M}$ has closed (elliptic or circular) boundaries. Figs. 6 (a)–(e) show the evolutions of some select physical trajectories in the $XY$ phase plane, for the exemplary cases $\mathfrak{w}=0.5,1,5,25$ and $50$.

We see that the phase portraits are quite similar to those in the Einstein frame (see Figs. 1 (a)–(e)), except of course the presence of the inner boundary of $\mathfrak{M}$, and that the size of both the boundaries reduce with increasing $\mathfrak{w}$ (i.e. decreasing $\beta$) not only along $X$ but also along $Y$. The trajectories terminate at the stable point $J_{4}$ on the outer boundary. On the $X$-axis, and vertically beneath $J_{4}$, we have the saddle point $J_{3}$ which, although remains outside the admissible region $\mathfrak{M}$, has significance in funnelling the trajectories towards $J_{4}$. As $\mathfrak{w}$ increases, both $J_{3}$ and $J_{4}$ move towards the origin (with the shrinkage of $\mathfrak{M}$), and so do the unstable points $J_{1}$ and $J_{2}$ at the two intersections of the outer elliptic boundary and the $X$-axis. The other CP $J_{5}$ (not shown in Figs. 6 (a)–(e)) is beyond the outer boundary of $\mathfrak{M}$ for all values of $\mathfrak{w}>0$.

4.4 Dynamical evolution in the effective scenario

The task that remains now is to verify whether the dynamical evolution of the universe at a stable point is in accord with the exact solution we have found in the paper I by explicitly solving the Jordan frame MST-cosmological equations [50]. To do this, it is imperative to consider the effective scenario in the Jordan frame, in which we have critical density $\rho$ [cf. Eq. (4.6)] that satisfies the usual conservation relation (4.7). In fact, so far all our arguments have been based on the Jordan frame matter density parameter $\Omega^{(m)}_{J}$ or(and) the total EoS parameter $\mathsf{w}_{\!{}_{J}}$. The interpretations of the dynamical evolution of the universe in terms of an effective DE constituent is facilitated by the decomposition (4.6) of the (effective) critical density $\rho$ in the dust-like matter density $\rho^{(m)}$ and the DE density $\rho_{\!{}_{X}}$. Let us define the corresponding density parameters respectively as $\Omega^{(m)}:={}^{\displaystyle\rho^{(m)}\!\!}\!\!\!\diagup_{\!\!\displaystyle\rho}$ and $\Omega_{\text{\ssmall{\it X}}}:={}^{\displaystyle\rho_{\!{}_{X}}}\!\!\!\diagup_{\!\!\displaystyle\rho}$, and also the EoS parameters for the DE and the system respectively as $\mathsf{w}_{\!{}_{X}}:={}^{\displaystyle p_{\!{}_{X}}}\!\!\!\diagup_{\!\!\displaystyle\rho_{\!{}_{X}}}$ and $\mathsf{w}:={}^{\displaystyle p}\!\!\!\diagup_{\!\!\displaystyle\rho}$. Then from Eqs. (4.6), (4.8), (4.13) and (4.14) we have

[TABLE]

where

[TABLE]

What the exponential factor $e^{\sqrt{6}\,F}$ in Eq. (4.20) does is that it brings back the situation we have had in the Einstein frame (see subsection 3.2). That is, the matter density parameter $\Omega^{(m)}$ would keep on evolving with time (or $N$) even after the system reaches a CP. So the condition $\Omega^{(m)}<1$ would get violated eventually (at some later epoch), unless we set $\Omega^{(m)}=0$ in the asymptotic limit. Hence, the physical relevance of the phase plane trajectories requires the corresponding CPs to lie on the curve whose equation is

[TABLE]

This curve therefore encloses (or demarcates) the physically admissible region(s) $\mathfrak{M}$ in the phase. Note also that this curve is actually the outer boundary (or demarcation line(s)) of $\mathfrak{M}$ we have had while working in terms of the conventional matter density parameter $\Omega^{(m)}_{J}$ (in subsections and 4.2 and 4.3). We do not have any inner boundary of $\mathfrak{M}$ now, or we may say that the inner boundary is trivial (i.e. coincident with the origin $(0,0)$) as in the Einstein frame. Nevertheless, the situation here still has one difference with that in the Einstein frame. That is, with increasing (presumably positive) values of $\mathfrak{w}$ (i.e. with decreasing $\beta$), the size of $\mathfrak{M}$ shrinks along both $X$ and $Y$ (instead of along $X$ only, as in the Einstein frame). As to the phase portraits in Figs. 6 (a)–(e), since $\mathfrak{M}$ gets enlarged with the removal of the inner boundary, we have more choices for obtaining the trajectories numerically. Accordingly the saddle point $J_{3}$ (situated on the $X$-axis) now has more trajectories to deflect towards the stable point $J_{4}$ at the boundary (4.22).

Now, the cosmic acceleration being the main aspect of the exact Jordan frame solution found in paper I, one expects this solution to transpire to one of the stable points $J_{2}$ and $J_{4}$ in the asymptotic limit, since they are the only CPs that support accelerating cosmologies. By inspection we find that $J_{4}$ is the appropriate one, since at this CP $\left(X_{c}={}^{\displaystyle n}\!\!\!\diagup_{\!\!\displaystyle\!\!\!\sqrt{6}},\,Y_{c}=\sqrt{\left(1+{}^{\displaystyle n}\!\!\!\diagup_{\!\!\displaystyle 2}\right)\left(1+{}^{\displaystyle n}\!\!\!\diagup_{\!\!\displaystyle 3}\right)}\right)$, where $n=(1+\mathfrak{w})^{-1}$, the torsion scalar is given by $e^{\tau}=a^{{}^{\scriptstyle{3n\!}}\!\!\diagup_{\!\!\scriptstyle{\!2}}}$, whence the Hubble parameter:

[TABLE]

turns out to be precisely the asymptotic (i.e. the $a\rightarrow\infty$ limiting) form of that we have had in paper I, while deriving the exact solution in the Jordan frame (see section 5.1 therein). The stability of such a solution is thus established.

Let us finally look into the evolution of say, the EoS parameters $\mathsf{w}_{\!{}_{X}}(N)$ and $\mathsf{w}(N)$ over some fiducial trajectory, typically for two cases $\mathfrak{w}=5$ and $\mathfrak{w}=50$ (see Figs. 7 (a) and (b) respectively). As before, we take this trajectory (in either of these cases) to be the one in compliance with the initial conditions $X(0)$ and $Y(0)$, set with $\Omega^{(m)}_{{}_{0}}=0.3$ and keeping in mind the power-law ansatz $\,e^{\tau}=a^{{}^{\scriptstyle{3n\!}}\!\!\diagup_{\!\!\scriptstyle{\!2}}}$ for the stable solution in [50]. Eqs.(4.10) and (4.13) then give

[TABLE]

With $n=(\mathfrak{w}+1)^{-1}$, we therefore have $\left(X(0),Y(0)\right)=(0.068,0.918)$ for $\mathfrak{w}=5$, and $\left(X(0),Y(0)\right)=(0.008,0.846)$ for $\mathfrak{w}=50$.

The EoS parameter $\mathsf{w}$ for the system evolves in a similar way as in the Einstein frame, viz. remains almost fixed at a small positive value deep in the past, and then in the near past makes a steep transition to a value very close to $-1$, at which it saturates. At the present epoch ($N=0$) the universe is in the transition phase, which ends in the near future. The smaller the value of the parameter $\mathfrak{w}$, the steeper is the transition. The EoS parameter $\mathsf{w}_{\!{}_{X}}$ for the DE, on the other hand, also shows a transition from a positive value to a negative value, which is even more steep (compared to that of $\mathsf{w}$ for the same value of $\mathfrak{w}$). Deep in the past $\mathsf{w}_{\!{}_{X}}$ increases very slowly with $N$ from a small positive value, until reaching a maximum in the near past. Then it decreases rapidly and attains a minimum value less than $-1$, increases again and asymptotically tends to $-1$ (from below). The transition to the minimum value ends before the present epoch, whence the universe continues to be in the super-accelerating or phantom regime (characterized by $\mathsf{w}_{\!{}_{X}}<-1$)131313Note however that the total EoS parameter $\mathsf{w}$ of the system never goes below $-1$. Therefore, physical instabilities against metric perturbations [66, 67] can not arise.. Such an effective super-acceleration is actually the striking feature of the Jordan frame solution obtained in paper I. Furthermore, the smaller the value of $\mathfrak{w}$ the steeper is the transition of $\mathsf{w}_{\!{}_{X}}$ from its maximum value to its minimum value, and the greater are these values in magnitude (see the exaggerated portions shown in the insets of Figs. 7 (a) and (b)).

5 Conclusion

We have thus made a systematic analysis of the dynamical stability of cosmological DE solutions in the scalar-tensor equivalent formulation of a non-minimal metric-scalar coupling with torsion (viz. the MST coupling). The roles of the individual torsion modes, viz. the trace $\mathcal{T}_{\mu}$ and the pseudo-trace $\mathcal{A}_{\mu}$, on the cosmological dynamics are envisaged by working in terms of a suitably defined torsion scalar $\tau$ and a torsion constant $\Lambda$. Apart from $\tau$, we have taken for simplicity just a single component cosmological matter, viz. the pressureless ‘dust’, thus restricting our analysis to a two-dimensional phase space, or the phase plane. Not only the analysis applies to the exact solutions found in paper I, in Einstein and Jordan frames [50], but also to other solutions plausible for a cosmic acceleration in the MST-cosmological setup. Essentially, we have had the Einstein frame and Jordan frame cosmological equations cast in the form of the respective autonomous sets, whose equilibrium solutions (in the asymptotic limit) are the critical points (CPs). Identifying the type and nature of these CPs, by applying the linear perturbation theory, we have analysed the phase plane trajectories of physical relevance, and the solutions they support, satisfying the requirement of cosmic acceleration in the asymptotic limit.

In the Einstein frame, we have found as many as five CPs existing, including stable and saddle points. Nevertheless, the location of these CPs in the phase plane, and the cosmological dynamics leading up to them, are quite different from that in the case of a quintessence field (with an exponential potential) in presence of dust [2, 11]. The reason is two-fold: firstly, we have the torsion scalar $\tau$ interacting with the cosmological matter in the Einstein frame MST setup, and secondly, the phase variables need to be such that their calibrations are not affected by the MST coupling parameter $\beta$. For convenience, we have chosen to work in an effectively non-interacting scenario, in which the physical acceptability of cosmological solutions crucially depends on the eventual extinction of the (dust-like) matter. This implies that the CPs are of physical relevance only when they lie on an ellipse centred at the origin. Such an ellipse therefore forms the boundary of the physically admissible region $\mathfrak{M}$ for the phase plane trajectories, at a given value of $\beta$. With increasing $\beta$, the eccentricity of the ellipse decreases, and the area of $\mathfrak{M}$ increases, whence the CPs move further and further away from the origin. Of prime importance is the domain $0<\beta<{}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$, within which four CPs exist, including the stable point that supports solutions exhibiting cosmic acceleration asymptotically. An exact DE solution of such sort is that found in paper I. We have not only verified the stability of this solution, but also have numerically obtained the evolution profiles of the DE density and EoS parameters, $\Omega_{\text{\ssmall{\it X}}}$ and $\mathsf{w}_{\!{}_{X}}$, over a fiducial trajectory whose initial conditions are set in accord with such a solution. Nevertheless, our emphasis has mostly been on a much smaller value of $\beta$ ($\sim 10^{-1}$ or lesser), so that the DE evolution remains within the proximity of the concordance $\Lambda$CDM model.

In the Jordan frame, we have had the option to resort to either the conventional scenario in which the critical density of the universe is not conserved and depends explicitly on the running gravitational coupling factor, or the effective scenario in which the critical density is by definition conserved and decomposed into the densities due to the dust and a DE component, which are individually conserved. Although the effective scenario is more useful from the observational perspective, we have carried out the dynamical analysis first in the conventional scenario, which is in line with the standard scalar-tensor approaches in the literature [11]. In due course, we have shown that the general results of such an analysis do not differ in the effective scenario, except the broadening of the admissible regions $\mathfrak{M}$ for physical trajectories in the phase plane (see the full account given in §4.2 on the boundaries of such regions at different domains of the effective Brans-Dicke (BD) parameter $\mathfrak{w}\,(\sim\beta^{-1})$). We have found five distinct CPs existent, as in the Einstein frame. Examining the type and nature of these CPs and the solutions they support, we have once again identified the stable point at which the universe remains in an accelerating state of expansion asymptotically, without any physical instability against cosmological metric perturbations. In the effective scenario, we have not only verified that the exact Jordan frame solution found in paper I indeed transpires to this stable CP in the asymptotic limit, but also have worked out numerically the evolution profiles of the EoS parameters $\mathsf{w}$ and $\mathsf{w}_{\!{}_{X}}$, for the system and the DE respectively, over fiducial trajectories leading up to the stable CP. While the evolution of $\mathsf{w}$ is found to be similar to that in the Einstein frame, $\mathsf{w}_{\!{}_{X}}$ has shown the striking feature of crossing the phantom barrier at an epoch in the recent past.

On the whole, the main aspect of the MST-cosmological analysis is the indispensable role of the torsion pseudo-trace mode $\mathcal{A}^{\mu}$ in the stable DE configurations, in both Einstein and Jordan frames. The trace mode $\mathcal{T}_{\mu}$ (or its scalar field source $\phi$), although crucial in the MST formalism, only has a supplementary role in the DE evolution. This is evident from the fact that the effective measure of $\mathcal{A}^{\mu}$, viz. the torsion constant $\Lambda$, ensures the culmination of a phase trajectory at the stable point that supports the accelerated expansion of the universe. On the other hand, the torsion scalar $\tau$, derived from $\mathcal{T}_{\mu}$, is responsible for a dynamical evolution of the DE. The dynamics is weak enough as long as the coupling parameter $\beta$ is small. While in the Einstein frame we have smallness of $\beta$ acquiescing to the stringent bound $\beta<{}^{\scriptstyle{1}}\!\!\diagup_{\!\!\scriptstyle{2}}$, in the Jordan frame it implies a large BD parameter $\mathfrak{w}$, which is in accord with many independent studies [56, 57]. So it is a generic outcome of our analysis that stable DE models (including those found in paper I) [50] result from a feeble metric-scalar coupling with torsion (no matter how weak the latter is).

In principle, the dynamical analysis carried out here can be extended to that applicable to any scalar field DE model in which the scalar interacts with the dust-like matter. The effective scenario is best suited for this. One has to simply consider the (usual) critical density, viz. $\rho:={}^{\displaystyle 3H^{2}\!\!}\!\!\!\diagup_{\!\!\displaystyle\!\!\kappa^{2}}$, and its splitting into two non-interacting components — the dust and whatever is left over (let that be the DE), and carry out the analysis just as in this work. However, note that the results of the analysis would always be model specific, because the phase variables for any given model are required to be so chosen that they are not explicitly dependent on a system parameter, different domains of which have diverse implications in the analysis. In that sense, the results obtained in this paper have their uniqueness attributed to the analysis we have carried out by keeping the phase variables independent of our system parameter $\beta$ (or $\mathfrak{w}$) explicitly. Of actual importance is of course our definition of the phase variables in terms of the torsion scalar $\tau$ or the torsion constant $\Lambda$, rather than in terms of a redefined quintessence-like scalar field or its potential.

Let us end this paper with some open issues, such as: (i) what happens if instead of the torsion constant $\Lambda$, we consider some other contribution of the torsion pseudo-trace $\mathcal{A}^{\mu}$ (sourced by say, the Kalb-Ramond axion in string-inspired scenarios [33, 34])? (ii) what if we introduce an interaction of a generic form between torsion and the cosmological dust (or any other perfect fluid matter)? (iii) how about correlating the results of the dynamical analysis here for our MST-cosmological system, with that for the models of coupled quintessence, coupled tachyon etc. (as for e.g. in some recent works [68] and references therein)? (iv) how about extending the dynamical analysis in this paper to say, the Chaplygin gas or Chameleon cosmologies in the MST framework, or to generalize it for the modified gravity theories, such as $f(\mathcal{R})$ or mimetic gravity [18, 17, 64], in presence of torsion? and so on. Works addressing some of these are in progress [69, 70] and we hope to report them soon.

Acknowledgements

The work of ASB was supported by the Council of Scientific and Industrial Research (CSIR), Government of India. SS acknowledges the R & D Grant DRCH/R & D/2013-14/4155, Research Council, University of Delhi.