Anisotropic cosmological dynamics in $f(T)$ gravity in the presence of a perfect fluid

TL;DR

This paper investigates how anisotropic cosmological models in $f(T)$ gravity evolve with different matter contents, revealing new singularity behaviors and multiple evolution scenarios through numerical analysis.

Contribution

It introduces a detailed analysis of anisotropic $f(T)$ gravity with perfect fluid, identifying new singularity types and evolution scenarios not present in general relativity.

Findings

Matter content influences the nature of cosmological singularities.

New anisotropic and isotropic regimes can replace Kasner solutions.

Six distinct evolution scenarios are numerically identified.

Abstract

We consider the cosmological evolution of a flat anisotropic Universe in $f(T)$ gravity in the presence of a perfect fluid. It is shown that the matter content of the Universe has a significant impact of the nature of a cosmological singularity in the model studied. Depending on the parameters of the $f(T)$ function and the equation of state of the perfect fluid in question the well-known Kasner regime of general relativity can be replaced by a new anisotropic solution, or by an isotropic regime, or the cosmological singularity changes its nature to a non-standard one with a finite values of Hubble parameters. Six possible scenarios of the cosmological evolution for the model studied have been found numerically.

Figures

Equations

$$S=\frac{1}{2K}\int ef(T)d^{4}x+S_m,$$

$$e_{\mu }^A=\mathrm{diag}(1,a(t),b(t),c(t)),$$

$$T=-2(H_a{H}_b+H_a{H}_c+H_b{H}_c),$$

$$\dot{T}=-2({\dot{H}}_a(H_b+H_c)+{\dot{H}}_b(H_a+H_c)+{\dot{H}}_c(H_a+H_b)),$$

$$f(T)-2T{f}_T=2K\rho ,$$

$$\dot{T}{f}_{TT}(H_b+H_c)+\frac{f}{2}+f_T\left({\dot{H}}_b+{\dot{H}}_c+{(H_b)}^2+{(H_c)}^2+2H_b{H}_c+H_a{H}_b+H_a{H}_c\right)=-Kw\rho ,$$

$$\dot{T}{f}_{TT}(H_a+H_c)+\frac{f}{2}+f_T\left({\dot{H}}_a+{\dot{H}}_c+{(H_a)}^2+{(H_c)}^2+2H_a{H}_c+H_a{H}_b+H_b{H}_c\right)=-Kw\rho ,$$

$$\dot{T}{f}_{TT}(H_a+H_b)+\frac{f}{2}+f_T\left({\dot{H}}_a+{\dot{H}}_b+{(H_a)}^2+{(H_b)}^2+2H_a{H}_b+H_a{H}_c+H_b{H}_c\right)=-Kw\rho ,$$

$$\dot{\rho }+(H_a+H_b+H_c)(1+w)\rho =0.$$

$$H_a\dot{T}{f}_{TT}+f_T\left({\dot{H}}_a+{H_a}^2-H_b{H}_c\right)=-\frac{1}{2}K(w+1)\rho .$$

$$H_b\dot{T}{f}_{TT}+f_T\left({\dot{H}}_b+{H_b}^2-H_a{H}_c\right)=-\frac{1}{2}K(w+1)\rho ,$$

$$H_c\dot{T}{f}_{TT}+f_T\left({\dot{H}}_c+{H_c}^2-H_a{H}_b\right)=-\frac{1}{2}K(w+1)\rho .$$

$$(H_a+H_b+H_c)\dot{T}{f}_{TT}+f_T\left(({\dot{H}}_a+{\dot{H}}_b+{\dot{H}}_c)+{(H_a+H_b+H_c)}^2+\frac{3}{2}T\right)=-\frac{3}{2}K(w+1)\rho .$$

$$(H_b-H_a)\dot{T}{f}_{TT}+f_T\left({\dot{H}}_b-{\dot{H}}_a+{H_b}^2-{H_a}^2+H_c(H_b-H_a)\right)=0,$$

$$(H_c-H_b)\dot{T}{f}_{TT}+f_T\left({\dot{H}}_c-{\dot{H}}_b+{H_c}^2-{H_b}^2+H_a(H_c-H_b)\right)=0,$$

$$(H_a-H_c)\dot{T}{f}_{TT}+f_T\left({\dot{H}}_a-{\dot{H}}_c+{H_a}^2-{H_c}^2+H_b(H_a-H_c)\right)=0.$$

$$2f_T\left(({\dot{H}}_b-{\dot{H}}_a)(H_c-H_b)+({\dot{H}}_c-{\dot{H}}_b)(H_a-H_b)\right)=0.$$

$$2f_T\left(({\dot{H}}_c-{\dot{H}}_b)(H_a-H_c)+({\dot{H}}_a-{\dot{H}}_c)(H_b-H_c)\right)=0,$$

$$2f_T\left(({\dot{H}}_a-{\dot{H}}_c)(H_b-H_a)+({\dot{H}}_b-{\dot{H}}_a)(H_c-H_a)\right)=0.$$

$${\dot{H}}_c(H_a-H_b)={\dot{H}}_a(H_c-H_b)+{\dot{H}}_b(H_a-H_c).$$

$$\begin{array}[]{l}\frac{\dot{H}_{b}-\dot{H}_{a}}{H_{b}-H_{a}}=\frac{\dot{H}_{c}-\dot{H}_{b}}{H_{c}-H_{b}}=\frac{\dot{H}_{a}-\dot{H}_{c}}{H_{a}-H_{c}}~{}~{}\Rightarrow\\ \Rightarrow~{}~{}\frac{d\Big{(}\ln(H_{a}-H_{b})\Big{)}}{dt}=\frac{d\Big{(}\ln(H_{c}-H_{b})\Big{)}}{dt}=\frac{d\Big{(}\ln(H_{a}-H_{c})\Big{)}}{dt}.\end{array}$$

$$H_b-H_a=C_1(H_c-H_b),$$

$$H_c-H_b=C_2(H_a-H_c),$$

$$H_b-H_a=C_3(H_a-H_c),$$

$$H_c=-\frac{1}{C_1}{H}_a+\frac{1+C_1}{C_1}{H}_b.$$

$$\begin{array}[]{l}T+f_{0}T^{N}-2T-2f_{0}NT^{N}=2K\rho~{}~{}\Rightarrow\\ \Rightarrow~{}~{}-T+f_{0}(1-2N)T^{N}=2K\rho,\end{array}$$

$$H_a\dot{T}{f}_{0}N(N-1)T^{N-2}+(1+f_{0}N{T}^{N-1})\left({\dot{H}}_a+{H_a}^2-H_b{H}_c\right)=-\frac{1}{2}K(w+1)\rho ,$$

$$H_b\dot{T}{f}_{0}N(N-1)T^{N-2}+(1+f_{0}N{T}^{N-1})\left({\dot{H}}_b+{H_b}^2-H_a{H}_c\right)=-\frac{1}{2}K(w+1)\rho ,$$

$$H_c\dot{T}{f}_{0}N(N-1)T^{N-2}+(1+f_{0}N{T}^{N-1})\left({\dot{H}}_c+{H_c}^2-H_a{H}_b\right)=-\frac{1}{2}K(w+1)\rho .$$

$$\begin{array}[]{l}(H_{a}+H_{b}+H_{c})\dot{T}N(N-1)T^{N-2}+\\ +(1+f_{0}NT^{N-1})\left((\dot{H}_{a}+\dot{H}_{b}+\dot{H}_{c})+{(H_{a}+H_{b}+H_{c})}^{2}+\frac{3}{2}T\right)=-\frac{3}{2}K(w+1)\rho.\end{array}$$

Full text

1 Introduction

Recently a new class of modified gravity theories has attracted a great deal of attention. Its story goes back to the beginning of the 20th when Einstein reformulated his General Relativity (GR) in therms of torsion instead of curvature [1] and the Weitzenböck connection [2] instead of the Levi-Civita one. The resulting theory, known as the teleparallel equivalent of general relativity (TEGR) [3], has the same equations of motion as in GR, but it has a different mathematical structure. This theory has been forgotten for a long time and only a decade ago it attracted new attention. As for TEGR itself, despite complete equivalence to GR on the level of the equations of motion, its different mathematical background leads to different global properties, such as conservation quantities, the topic which is currently the subject of intense investigations [3, 4]. What has also become understood only after the theory’s revival is that modifying TEGR the way analogous to the famous $f(R)$ modifications of GR leads to the theory (known as $f(T)$ gravity; see the review in [5]) which is different from $f(R)$ already at the level of equations of motion [6, 7, 8]. The function $f(T)$ is assumed to be smooth; other possible restrictions from thermodynamical considerations have been discussed in Ref. [9]. The equations of motion in $f(T)$ theory is of the second order, to be compared with the fourth-order equations in $f(R)$. However, $f(T)$ gravity had some still unresolved conceptional problems connected with the lack of local Lorentz covariance [10, 11] and an unknown number of dynamical degrees of freedom (see [12] and references therein).

For the equations of motion this leads to a situation when non-pathological equations require a special (so called “proper” or “good”) tetrad to be used [4, 13, 14, 15], or equivalently, a non-trivial spin connection to be calculated for a given tetrad [16]. Otherwise, the non-symmetric part of the equations of motion requires explicitly the condition $\frac{d^{2}f(T)}{dT^{2}}=0$, which means $f(T)\propto T$ and we return to unmodified TEGR.

Before this problem can be understood in detail, one of the possible ways to study $f(T)$ gravity is finding solutions for the cases where the “proper” tetrad is known. For example, it is well known that for a flat Friedmann-Lemaître-Robertson-Walker (FLRW) cosmology with a scale factor $a$ the diagonal cartesian tetrad $(1,a,a,a)$ is a proper tetrad. This means that the general problem does not affect studies of FLRW cosmology in $f(T)$ gravity [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Moreover, it was shown that the diagonal tetrad $(1,a,b,c)$ is a proper tetrad for a flat anisotropic Universe (here $a$, $b$, $c$ are three different scale factors of the Bianchi I metric), which allows us to study corresponding anisotropic cosmological dynamics [5, 28, 29, 30, 31, 32, 33, 34, 35, 36].

In our recent paper [36] we have considered Bianchi I cosmological dynamics of a vacuum Universe in $f(T)=T+f_{0}T^{N}$ gravity. We have found that, depending on the branch of solutions ($T=0$ or $T\neq 0$), the Kasner solution [37] is either an exact (for $T=0$ branch) or asymptotic solution (for $T\neq 0$ branch) in the high-energy regime where the corrections to Einstein gravity dominate. However, intuitively we would expect that deviations from the solutions of Einstein gravity (the Kasner solutions belong to this class) should increase as the higher order terms in the action of the theory start to dominate. What we found is the opposite behavior — deviations starts to dominate in low-energy dynamics (where a new de Sitter solution appears), while the high-energy regime is the same for TEGR and $f(T)$ gravity. Such a property leads to the suggestion that the model we have studied is oversimplified. Indeed, our considerations have been restricted by vacuum models, which are very different from the Universe we live in. In GR, vacuum approximation near a cosmological singularity is known to be relevant — any matter content of the Universe (apart from a stiff fluid) has negligible influence on the cosmological dynamics, and the nature of the cosmological singularity is determined solely by the vacuum solutions. However, it is not a priory evident whether this property is still valid in modified gravity. The goal of the present paper is to study cosmological evolution in a flat anisotropic Universe in the presence of a perfect fluid and to determine the influence of the matter content of the Universe upon its dynamics.

2 The equations of motion

We consider cosmological models with the action of $f(T)$ theory with matter

[TABLE]

where $e=\text{det}(e^{A}_{\mu})=\sqrt{-g}$ is the determinant consisting of the tetrad components $e^{A}_{\mu}$,   $f(T)$ is a general differentiable function of the torsion scalar $T$,   $S_{m}$ is the matter action and $K=8\pi G$. Units $\hbar=c=1$ will be used.

The following diagonal tetrad is chosen:

[TABLE]

which relates to the Bianchi I metric   $\mathrm{d}s^{2}=\mathrm{d}t^{2}-a^{2}(t)\mathrm{d}x^{2}-b^{2}(t)\mathrm{d}y^{2}-c^{2}(t)\mathrm{d}z^{2}$,   where   $a(t)$,   $b(t)$,   $c(t)$   are scale factors. The torsion scalar for the chosen tetrad (2) is

[TABLE]

where $H_{a}\equiv\frac{\dot{a}}{a}$,   $H_{b}\equiv\frac{\dot{b}}{b}$,   $H_{c}\equiv\frac{\dot{c}}{c}$ are anisotropic Hubble parameters, and a dot denotes the derivative with respect to time. Equation (3) reduces to $T=-6H^{2}$ in the isotropic case $a(t)=b(t)=c(t)$,   $H_{a}=H_{b}=H_{c}=H$. Then the time derivative of torsion scalar has the form

[TABLE]

We derive the equations of motion varying the action (1) with respect to the chosen tetrad (2) (see, for example, [28])

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho$ is the energy density of matter, $p$ the pressure of matter, $p=w\rho$ the matter equation of state; $w\in[-1;1]$ is a constant. Here we denote $f_{T}=\frac{df(T)}{dT}$,    $f_{TT}=\frac{d^{2}f(T)}{dT^{2}}$.

The continuity equation for matter is

[TABLE]

We subtract (6) from the sum of (7) and (8), substitute the constraint $K\rho=\frac{f}{2}-Tf_{T}$ and (3) to the obtained equation and find

[TABLE]

The following two equations are obtained analogously:

[TABLE]

[TABLE]

The sum of these three equations is

[TABLE]

3 The linear relation between $H_{a}$, $H_{b}$, $H_{c}$

The system under investigation contains four variables — three Hubble parameters and the matter energy density. We have also one constraint equation (5), so we could expect that the number of independent variables is equal to three. However, the special very symmetric nature of equations of motion induces one more relation between the Hubble parameters.

To show this we subtract in pairs the equations of the system (10)-(12) and find

[TABLE]

[TABLE]

[TABLE]

Now we add Eq. (14) multiplied by $(H_{c}-H_{b})$ to Eq. (15) multiplied by $(H_{a}-H_{b})$

[TABLE]

Analogously to (17) we obtain

[TABLE]

[TABLE]

The following expression is found from these system for $f_{T}\neq 0$:

[TABLE]

For $f_{T}\neq 0$, $H_{a}\neq H_{b}\neq H_{c}$ the system (17)-(19) gives us

[TABLE]

Solving these differential equations we find

[TABLE]

[TABLE]

[TABLE]

where $C_{3}=C_{1}C_{2}$ is obtained after the substitution (23) to (22). The sum (23) and (24) is $H_{c}-H_{a}=(C_{2}+C_{1}C_{2})(H_{a}-H_{c})$. Then $C_{2}+C_{1}C_{2}=-1~{}~{}\Rightarrow~{}~{}C_{2}=-\frac{1}{1+C_{1}}$ and $C_{3}=-\frac{C_{1}}{1+C_{1}}$. Therefore, $H_{c}$ is the linear combination of $H_{a}$ and $H_{b}$ for the assumptions $f_{T}\neq 0$ and $H_{a}\neq H_{b}\neq H_{c}$:

[TABLE]

4 Dynamical system for the model $f(T)=T+f_{0}T^{N}$

In what follows we consider cosmological models with the Lagrangian density function $f(T)=T+f_{0}T^{N}$, where $f_{0}$,   $N>0$ are parameters. Then the constraint (5) and others field equations (10)-(12) have the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We find the sum of Eqs. (27)-(29) to be

[TABLE]

Due to the linear relation (25) we can expect that two variables are enough for the corresponding dynamical system. For such variables we choose the torsion scalar and the sum of the three Hubble parameters; therefore, we shall use Eqs. (26) and (30) in the present section. Separate behavior of Hubble parameters will be considered later in the following sections.

New expansion-normalized variables are introduced as follows:

[TABLE]

where we denote $X=H_{a}+H_{b}+H_{c}$.

The variable $y$ depends on $T^{N-1}$, therefore we can express $T^{N-1}$ and $f_{T}$ through $y$:

[TABLE]

then

[TABLE]

[TABLE]

and

[TABLE]

It follows from the definition of variable $y$ that

1). if $y\to 0$ then $T^{N-1}\to const=\frac{1}{f_{0}(1-2N)}$, $N\neq\frac{1}{2}$,

2). if $y\to-1$ then $T\to 0$,

3). if $y\to(1-2N)/N$ then $T\to\pm\infty$.

Dividing the constraint equation (26) by $X^{2}f_{T}$ we find

[TABLE]

Therefore, the variable $r$ is not independent and can be excluded. Now we divide (13) by $X^{2}f_{T}$ and obtain

[TABLE]

Taking the time derivative of (26) we have

[TABLE]

Using (32), (33), (36) and (38) we obtain the expression for the quantities $\frac{\dot{T}}{X^{3}}$ and $\frac{\dot{T}}{TX}$:

[TABLE]

and

[TABLE]

Taking into account (35) we have

[TABLE]

We differentiate the new variables $x$, $y$ with respect to   $\ln(abc)$,   where d\Big{(}\ln(abc)\Big{)}=dt\left(\frac{\dot{a}}{a}+\frac{\dot{b}}{b}+\frac{\dot{c}}{c}\right)=Xdt, and we find the system of first-order differential equations

[TABLE]

[TABLE]

Using Eqs. (39), (40), (41) we finally obtain

[TABLE]

[TABLE]

5 Stationary points

We solve the system (44), (45), find stationary points and calculate the corresponding eigenvalues for each point in order to obtain their type of stability in the linear approach.

1. $x=0$, $y=(1-2N)/N$.

Eigenvalues are

[TABLE]

This stationary point is either an unstable node for $N>1$, $w\in\left(-1;\frac{1}{2N-1}\right)$ or a saddle point for $N>1$, $w\in\left(\frac{1}{2N-1};1\right]$. For $w=\frac{1}{2N-1}$, $N>1$, we find $\lambda_{1}=0$, $\lambda_{2}>0$ and this fixed point is unstable.

2. $x=-2/3$, $y=(1-2N)/N$.

The eigenvalues are calculated to be

[TABLE]

We see from the sings of the eigenvalues that this stationary point is either an unstable node for $N>1,w\in\left(\frac{1}{2N-1};1\right]$ or a saddle for $N>1,w\in\left(-1;\frac{1}{2N-1}\right)$. For $w=\frac{1}{2N-1}$, $N>1$, this point is unstable as $\lambda_{2}>0$.

3. $x=-2/3$, $y=-1$.

We find the eigenvalues

[TABLE]

It is a stable node for $N>1$, $w\in(-1;1)$.

4. $x=-2/3$, $y=0$.

The eigenvalues are obtained:

[TABLE]

which indicate that it is a stable node for $w\in(-1;1]$.

Using stationary point coordinates we calculate $\frac{\dot{X}}{X^{2}}=0$, $\frac{\dot{T}}{TX}=0$ and find $X(t)$, $T(t)$, $\rho(t)$,

[TABLE]

[TABLE]

[TABLE]

For $y=0$ we find ${T_{0}}^{N-1}=\frac{1}{f_{0}(1-2N)}$. As in this solution the sum of Hubble parameters $X$ and the torsion scalar $T$ are constants, then all Hubble parameters are equal, $H_{a}=H_{b}=H_{c}=H_{0}=\pm\sqrt{-\frac{T_{0}}{6}}$. This is the de Sitter solution.

5. In partial case of $w=w_{cr}=1/(2N-1)$ the stationary line exists:

$x\in(-\infty;+\infty)$, $y=(1-2N)/N$.

The eigenvalues for this stationary line are

[TABLE]

As $\lambda_{2}>0$ for $N>1$, this stationary line is unstable.

6 The asymptotic power-law solutions

Since most of known exact or asymptotic solutions in homogeneous cosmology have a power-law form, we check the existence of the asymptotic power-law solution of the form

[TABLE]

Then the torsion scalar $T$, its time derivative $\dot{T}$, the energy density of matter $\rho$ are

[TABLE]

[TABLE]

[TABLE]

We substitute the solution (54) into the initial system (26)-(29) for two asymptotic limits: 1). $t\to t_{0}$,  2). $t\to+\infty$. We shall consider only the case of $N>1$ in all calculations below.

1). For $t\to t_{0}$ the torsion scalar $|T|\to\infty$ and $|T^{N}|\gg|T|$. Then the term $-T$ can be neglected in comparison with $f_{0}(1-2N)T^{N}$ in the constraint (26)

[TABLE]

and using Eq. (57) we find

[TABLE]

We equate power indices and multipliers in Eq. (59) and obtain

[TABLE]

[TABLE]

Substituting the power-law solution (54) into Eq. (30) we obtain

[TABLE]

Then we rewrite the substitution of the solution (54) into the system (27)-(29) and subtract in pairs the obtained equations,

[TABLE]

[TABLE]

[TABLE]

Therefore, we have two asymptotic power-law solutions for $t\to t_{0}$:

a). Isotropic asymptotic regime with $p_{1}=p_{2}=p_{3}=p$,

[TABLE]

[TABLE]

where we have used (60), (61).

b). Anisotropic asymptotic regime with $p_{1}+p_{2}+p_{3}=2N-1$,

[TABLE]

[TABLE]

We find from (62)

[TABLE]

Here $p_{1}p_{2}+p_{1}p_{3}+p_{2}p_{3}\neq 0$.

The found isotropic solution a). corresponds to the fixed point 2 and the anisotropic one b). corresponds to the stationary line 5.

2). For $t\to+\infty$ we find that $T\to 0$, $|T^{N}|\ll|T|$. Neglecting $f_{0}(1-2N)T^{N}$ in (26) we obtain

[TABLE]

Taking into account (57) we have

[TABLE]

Equating power indices and multipliers in (72) we obtain

[TABLE]

[TABLE]

The substitution (54) into the system (27)-(29) and the subtraction in pairs of the resulting equations lead to

[TABLE]

[TABLE]

[TABLE]

Thus, we have the 3-d power-law solution in a form which exists for $t\to+\infty$

c). Isotropic asymptotic regime with $p_{1}=p_{2}=p_{3}=p$,

[TABLE]

where we have used (73),

[TABLE]

This isotropic solution c). corresponds to the fixed point 3.

7 More complicated asymptotic power-law solution

In the previous section we have identified the solutions in three fixed points from our list, one more (the de Sitter solution) has been identified in the Sect. 3. However, the nature of a solution in the point 1 is still to be determined. As $X=0$ at this point, it should correspond to some high-energy regime.

To reveal the meaning of the point 1 we try to check a more complicated form for the Hubble parameters:

[TABLE]

[TABLE]

[TABLE]

Here $p_{a1}$, $p_{a2}$, $p_{b1}$, $p_{b2}$, $p_{c1}$, $p_{c2}$, $\alpha$ are constants, which satisfy the following relations:

[TABLE]

[TABLE]

[TABLE]

The assumptions (83) and (84) lead to $p_{a1}\neq p_{b1}\neq p_{c1}$. Equation (85) shows that the first terms in (80)-(82) dominate for $t\to t_{0}$ and the Hubble parameters increase as ${(t-t_{0})}^{-1}$:

[TABLE]

[TABLE]

[TABLE]

The torsion scalar $T$ in this solution is

[TABLE]

Taking into account (84), (85) we find from (89) for $t\to t_{0}$

[TABLE]

Therefore, despite neglecting the second terms in the Hubble parameters $H_{a}$, $H_{b}$, $H_{c}$ for $t\to t_{0}$, these terms are important for the torsion scalar $T$.

The time derivative of torsion scalar is

[TABLE]

As $|T^{N}|\gg|T|$ for $T\to\pm\infty$, we find from the constraint equation (26)

[TABLE]

Substituting (80)-(82) to the continuity equation we have

[TABLE]

Then we substitute the expressions for $\rho$ and $T$ (93), (90) to (92) and find

[TABLE]

Equating the power-law indices and coefficients in (94) we obtain

[TABLE]

[TABLE]

We substitute complex power-law solution to the system (27)-(29), then subtract in pairs the obtained equations and find

[TABLE]

[TABLE]

[TABLE]

Here we keep only dominating terms with ${(t-t_{0})}^{-\alpha(N-1)+2}$. As our initial assumptions lead to $p_{a1}\neq p_{b1}\neq p_{c1}$, the factor $-1+p_{a1}+p_{b1}+p_{c1}-\alpha(N-1)$ in the system (97)-(99) should be equal to zero, which allows us to find $\alpha$:

[TABLE]

Taking into account (95) we have

[TABLE]

where $\alpha$ satisfies the initial assumption (85) $-1<\alpha-1<1$   for   $-1<w<\frac{1}{2N-1}$.

Therefore, the complex power-law solution exists with the following coefficients:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $p_{a1}p_{b1}+p_{a1}p_{c1}+p_{b1}p_{c1}=0$. The found solution corresponds to the stationary point 1.

8 Numerical examples of cosmological scenarios

Fixed points analysis can get information only about local properties of the dynamics. In the present section we consider global properties for models $f(T)=T+f_{0}T^{N}$ with integer $N$. Since they are different for even and odd $N$, we consider these cases separately, choosing $N=2$ and $N=3$ for numerical studies. First of all, the constraint equation imposes limitations on the possible values of $T$. In the case of $f(T)=T+f_{0}T^{2}$ the constraint equation (26) is a quadratic one for the torsion scalar $T$. It has two branches of solutions: $T_{1}$ and $T_{2}$. We find the intervals of allowed values of $T$ from the condition $\rho\geqslant 0$ (see Fig. 1, where the function $\rho(T)=\frac{-T+(1-2N)f_{0}T^{N}}{2K}$ is plotted for $N=2$):

1). for $f_{0}>0$ the torsion scalar $T\in\left[\frac{1}{-3f_{0}};0\right]$, where $\frac{1}{-3f_{0}}\leqslant T_{2}<\frac{1}{-6f_{0}}$ and $\frac{1}{-6f_{0}}\leqslant T_{1}\leqslant 0$,

2). for $f_{0}<0$ the torsion scalar T\in(-\infty;0]\cup\Big{[}\frac{1}{-3f_{0}};+\infty\Big{)}, where $T_{1}\leqslant 0$ and $T_{2}\geqslant\frac{1}{-3f_{0}}>0$.

In particular, we can see that the absolute value of the torsion scalar is restricted from above for $f_{0}>0$. This means that the points 1, 2 and the line 5 cannot be realized. Numerical studies (see below) indicate that instead of an analog of the Big Bang singularity the Universe starts its evolution from a non-standard singularity where $T$ is finite (and $H_{a}\neq H_{b}\neq H_{c}$ are finite) while $\dot{T}$ diverges. Note that for the isotropic case the same conclusion (large enough $H$ being unreachable if $f_{0}>0$) has been already drawn in [24].

We investigate numerically the model $f(T)=T+f_{0}T^{2}$ using the system (113)-(115) (see Appendix), which is obtained from the initial equations of motion. The following types of scenarios are realized.

For $f_{0}>0$ two variants of evolution are possible:

• Ia.

along the branch $T_{1}$ from the non-standard singularity with $T=\frac{1}{-6f_{0}}$ to the isotropic solution of the fixed point 3 (see the left plot in Fig. 2);

• Ib.

along the branch $T_{2}$ from the non-standard singularity with $T=\frac{1}{-6f_{0}}$ to de Sitter solution of the fixed point 4 (see the right graph in Fig. 2).

Note that the torsion scalar $T$ changes in a very restricted zone during a cosmological evolution along the Ib case — the initial absolute value of the torsion scalar |T_{in}|=|T_{0}|=\big{|}\frac{1}{-6f_{0}}\big{|} at the singularity ($T_{in}$ corresponds to the maximum of the function $\rho(T)$) is only two times less than its final absolute value |T_{fin}|=|T_{0}|=\big{|}\frac{1}{-3f_{0}}\big{|} at the de Sitter point. It is easy to see that for a general power-law function $f(T)=T+f_{0}T^{N}$ the ratio of final |T_{fin}|=|T_{0}|={\big{|}\frac{1}{f_{0}(1-2N)}\big{|}}^{\frac{1}{N-1}} and initial |T_{in}|={\big{|}\frac{1}{f_{0}N(1-2N)}\big{|}}^{\frac{1}{N-1}} values of the torsion scalar is equal to $N^{\frac{1}{N-1}}$. This means that this de Sitter point cannot be related to our Universe.

In the $f_{0}<0$ case points 1, 2 and the line 5 becomes accessible. We have found, however, one more regime which is not covered by the fixed points listed in the Sect. 3. Namely, this is a non-standard singular regime when $T$ approaches the point where $f_{T}=\frac{df(T)}{dT}=0$. Our numerical results indicate that for $f_{0}<0$ four scenarios exist:

• IIa.

along the branch $T_{1}$ from one of three solutions (corresponding to the stationary points 1, 2 or the line 5 depending on $w$) with infinite $T$ to the isotropic solution of the fixed point 3 (see Fig. 3);

• IIb.

along the branch $T_{2}$ with the start and the finish at the non-standard singularity with $T_{cr}=\frac{1}{-2f_{0}}$, where $T\to T_{cr}-0$ (see the left picture in Fig. 4);

• IIc.

along the branch $T_{2}$ from the anisotropic solution of the fixed point 1 for $w<\frac{1}{2N-1}$ to the non-standard singularity with $T_{cr}=\frac{1}{-2f_{0}}$, where $T\to T_{cr}+0$ (see the right plot in Fig. 4);

• IId.

along the branch $T_{2}$ with the start and the finish at the non-standard singularity with $T_{cr}=\frac{1}{-2f_{0}}$, where $T\to T_{cr}+0$ (see Fig. 5).

The scenarios similar to Ia, Ib, IIa are possible for an isotropic Universe and have been described in [24]. The scenarios IIb, IIc, IId are intrinsically anisotropic and have no isotropic analogs. It should be noted that the quantity $X=H_{a}+H_{b}+H_{c}=\frac{\dot{V}}{V}$, where $V\equiv a(t)b(t)c(t)$ is the volume factor. When $X=0$ the cosmological evolution passes through the turnaround (if $\dot{X}<0$) or the bounce (if $\dot{X}>0$). The turnaround exists in the case IIb and the bounce occurs in the scenario IId.

For the $N=3$ case the function $\rho(T)=\frac{-T+(1-2N)f_{0}T^{N}}{2K}$ is shown in Fig. 6. In this case the condition $\rho\geqslant 0$ gives us allowed values of $T$:

1). for $f_{0}>0$ the torsion scalar $T\in[-\infty;0]$,

2). for $f_{0}<0$ the torsion scalar T\in\Big{(}-\frac{1}{\sqrt{-5f_{0}}};0\Big{]}\cup\Big{[}\frac{1}{\sqrt{-5f_{0}}};+\infty\Big{)}.

From Fig. 6 we see that only an analog of the scenario IIa is possible for $f_{0}>0$. On the other hand, analogs of all remaining regimes found above for $N=2$ are possible for $N=3$ with $f_{0}<0$. We have found all of them in our numerical studies.

9 Conclusions

We have considered a cosmological evolution of a flat anisotropic Universe in $f(T)=T+f_{0}T^{N}$ gravity, with integer $N>1$ in the presence of a perfect fluid. Combining analytical and numerical methods we have identified possible asymptotic regimes. It is interesting that the Kasner solution does not appear in our list. We know that Kasner solution is an asymptotic regime (in the high-energy limit) for a vacuum cosmological models in $f(T)$ gravity. In GR it is an exact solution for vacuum models and an asymptotic solution for models with perfect fluid (apart from the stiff fluid). This means that the matter content of the Universe in $f(T)$ gravity, in contrast to GR, plays an important role in anisotropic dynamics near a cosmological singularity. For a power-law functions $f(T)=T+f_{0}T^{N}$ there exists a critical equation of state parameter for the matter $w_{cr}=1/(2N-1)$ such that for $w>w_{cr}$ the initial cosmological singularity appears to be isotropic. For $w<w_{cr}$ changes in comparison with GR are less drastic — the initial singularity is still anisotropic, leading terms in the time dependence of Hubble parameters are still power-law type, though the sum of the power indices is in general bigger than $1$. Moreover, for a particular sign of $f_{0}$ (positive for even $N$ and negative for odd $N$) the nature of the singularity changes completely from the standard one (where Hubble parameters as well as the matter energy density tend to infinity) to the non-standard one (where Hubble parameters and the matter energy density remain finite, and their time derivatives diverge). It is interesting that the de Sitter point exists only in such a situation. As a result, the torsion scalar $T$ can change at most two times during a cosmological evolution ending at de Sitter attractor. All these results indicate also that the famous GR situation described by the phrase “matter does not matter”, which describes the anisotropic cosmological dynamics near a cosmological singularity, does not occur in $f(T)$ gravity for any equation of state parameter $w$ bigger than $-1$.

Acknowledgements

M. A. S. and A. V. T. are supported by RSF Grant №16-12-10401 and A. V. T. is supported by the Russian Government Program of Competitive Growth of Kazan Federal University.

10 Appendix: The system for numerical integration

Using Eq. (20) for $\dot{H}_{c}$ we rewrite Eq. (4) for $\dot{T}$:

[TABLE]

where $H_{a}\neq H_{b}$, and we substitute (20) and (106) to Eqs. (10), (11),

[TABLE]

[TABLE]

Collecting terms with $\dot{H}_{a}$ and $\dot{H}_{b}$ we find

[TABLE]

[TABLE]

In order to express $\dot{H}_{a}$ from this system we multiply Eq. (109) to

\Big{(}f_{T}(H_{a}-H_{b})-4H_{b}f_{TT}({H_{a}}^{2}-H_{b}H_{c})\Big{)}, Eq. (110) by $4H_{a}f_{TT}({H_{a}}^{2}-H_{b}H_{c})$ and then we sum the obtained results,

[TABLE]

It is transformed to

[TABLE]

Finally, we have

[TABLE]

The expression for $\dot{H}_{b}$ is derived similarly to $\dot{H}_{a}$. It has the form

[TABLE]

Substituting Eqs. (113), (114) to Eq. (20) we find

[TABLE]

The system of equations (113)-(115), is used for the numerical integration, where $\rho=\frac{1}{2K}(f(T)-2Tf_{T})$, the function $f(T)$, its derivatives $f_{T}$, $f_{TT}$ and $T=-2(H_{a}H_{b}+H_{a}H_{c}+H_{b}H_{c})$ should be substituted.