Static Spherically Symmetric Einstein-aether models II: Integrability and the Modified Tolman-Oppenheimer-Volkoff approach

TL;DR

This paper explores the existence of analytic solutions in Einstein-ther theory for static spherically symmetric spacetimes, analyzing the modified TOV equations and their physical implications for perfect fluid and scalar field models.

Contribution

It provides a detailed dynamical system analysis and identifies conditions under which analytic solutions exist in Einstein-ther theory, highlighting modifications to the TOV equations.

Findings

Analytic solutions exist with Laurent expansion for perfect fluids with linear EoS.

The TOV equations are significantly modified in Einstein-ther theory.

Equilibrium points of fluid and scalar field models are characterized and analyzed.

Abstract

We investigate the existence of analytic solutions for the field equations in the Einstein-\ae ther theory for a static spherically symmetric spacetime and provide a detailed dynamical system analysis of the field equations. In particular, we investigate if the gravitational field equations in the Einstein-\ae ther model in the static spherically symmetric spacetime possesses the Painlev\`e property, so that an analytic explicit integration can be performed. We find that analytic solutions can be presented in terms of Laurent expansion only when the matter source consists of a perfect fluid with linear equation of state (EoS) $\mu =\mu _{0}+\left( \texttt{h} -1\right) p,~\texttt{h} >1$. In order to study the field equations we apply the Tolman-Oppenheimer-Volkoff (TOV) approach and other approaches. We find that the relativistic TOV equations are drastically modified in Einstein-\ae ther theory, and we explore the physical implications of this modification. We study perfect fluid models with a scalar field with an exponential potential. We discuss all of the equilibrium points and discuss their physical properties.

Figures

Tables

Table 1. Table 1: Equilibrium points and the corresponding eigenvalues of the dynamical system ( 57 ).

Labels	$(S_1,S_2,Y)$	Eigenvalues	Stability
$L_1$	$(S_1^{\ast },1-\text{b}S_1^{\ast }{}_{}{}^2,0)$	$0,-S_1^{\ast },-S_1^{\ast }$	Normally hyperbolic.
			Stable for $S_1^{\ast }>0$.
			Unstable for .
$L_2$	$(0,1,Y_0)$	$0,-Y_0,2Y_0$	Saddle.
$P_1$	$(\frac{2}{\text{h}+2},1-\frac{8\text{b}}{{\left(\text{h}+2\right)}^2},\frac{1}{\text{h}})$	$\frac{2}{\text{h}+2},-\frac{\text{h}+2+\sqrt{64\text{b}-7{\left(\text{h}+2\right)}^2}}{2\text{h}\left(\text{h}+2\right)},-\frac{\text{h}+2-\sqrt{64\text{b}-7{\left(\text{h}+2\right)}^2}}{2\text{h}\left(\text{h}+2\right)}$	Saddle.
$P_2$	$(\frac{1}{\sqrt{\text{b}}},0,\frac{1}{\text{h}})$	$\frac{2}{\text{h}},\frac{-\text{h}+4\sqrt{\text{b}}-2}{\sqrt{\text{b}}\text{h}},\frac{1}{\sqrt{\text{b}}}$	Source for
			$\text{h}>1,\text{b}>\frac{1}{16}{\left(\text{h}+2\right)}^2$,
			saddle otherwise.
$P_3$	$(-\frac{1}{\sqrt{\text{b}}},0,\frac{1}{\text{h}})$	$\frac{2}{\text{h}},-\frac{1}{\sqrt{\text{b}}},\frac{\text{h}+4\sqrt{\text{b}}+2}{\sqrt{\text{b}}\text{h}}$	Saddle.
$P_4$	$(\frac{\text{h}+2}{4\text{b}},0,\frac{1}{\text{h}})$	$\frac{\text{h}+2}{4\text{b}},\frac{{\left(\text{h}+2\right)}^2-8\text{b}}{4\text{b}\text{h}},\frac{{\left(\text{h}+2\right)}^2-16\text{b}}{8\text{b}\text{h}}$	Sink for
			.

Table 2. Table 2: Equilibrium points and the corresponding eigenvalues of the dynamical system defined by equations( 75 ). We use the notation Δ = 36 − 7 ​ h 2 − 28 ​ h Δ 36 7 superscript h 2 28 h \Delta=\sqrt{36-7\texttt{h}^{2}-28\texttt{h}} .

Labels	$(U,V,Y)$	Eigenvalues	Stability
$L_1$	$(1,0,Y_0)$	$0,{\left(1-Y_0\right)}^3,{\left(1-Y_0\right)}^3$	Source.
$L_2$	$(\frac{3}{4},0,Y_0)$	$0,\frac{1}{2}{\left(1-Y_0\right)}^3,-\frac{3}{4}{\left(1-Y_0\right)}^3$	Saddle.
$L_3$	$(U_0,0,1)$	$0,0,0$	Nonhyperbolic.
$L_4$	$(0,0,Y_0)$	$0,3{\left(1-Y_0\right)}^3,-{\left(1-Y_0\right)}^3$	Saddle.
$L_5$	$(0,V_0,1)$	$0,0,0$	Nonhyperbolic.
$L_6$	$(U_0,1,0)$	$0,-1+U_0,-1+U_0$	Sink.
$M_1$	$(0,\frac{1}{2},0)$, $(0,0)$	$\frac{3}{2},-\frac{1}{2},\frac{1}{2}$	Saddle.
$M_2$	$(0,1,0)$	$0,-1,-1$	Sink.
$M_3$	$(1,1,1)$	$1-\text{h},\text{h}-1,\text{h}-1$	Saddle.
$M_5$	$(1,1,0)$	$0,0,0$	Nonhyperbolic.
$M_6$	$(0,\frac{\text{h}-1}{\text{h}+1},\frac{1}{\text{h}})$	$-\frac{{\left(\text{h}-1\right)}^3}{{\text{h}}^2\left(\text{h}+1\right)},\frac{{\left(\text{h}-1\right)}^3}{{\text{h}}^2\left(\text{h}+1\right)},\frac{2{\left(\text{h}-1\right)}^3}{{\text{h}}^3\left(\text{h}+1\right)}$	Saddle.
$M_7$	$(1,\frac{\text{h}-1}{\text{h}+1},\frac{1}{\text{h}})$, $(\frac{2}{\sqrt{5}},\frac{1}{\sqrt{5}})$	$\frac{{\left(\text{h}-1\right)}^2}{{\text{h}}^2\left(\text{h}+1\right)},\frac{{\left(\text{h}-1\right)}^2}{{\text{h}}^2\left(\text{h}+1\right)},-\frac{2{\left(\text{h}-1\right)}^3}{{\text{h}}^3\left(\text{h}+1\right)}$	Saddle.
$M_8$	$(\frac{1}{2},\frac{2{\left(\text{h}-1\right)}^2}{3{\text{h}}^2-2},\frac{1}{\text{h}})$	$\frac{{\left(\text{h}-1\right)}^3}{\text{h}\left(3{\text{h}}^2-2\right)},-\frac{{\left(\text{h}-1\right)}^3\left(\text{h}+2-\Delta \right)}{4{\text{h}}^2\left(3{\text{h}}^2-2\right)},-\frac{{\left(\text{h}-1\right)}^3\left(\text{h}+2+\Delta \right)}{4{\text{h}}^2\left(3{\text{h}}^2-2\right)}$	Saddle.

Table 3. Table 3: Equilibrium points and the corresponding eigenvalues and stability conditions of the dynamical system ( 81 ).

Labels	$(U,V,Y)$	Eigenvalues	Stability
$L_1$	$(1,0,Y_0)$	$0,-\frac{{\left(Y_0-1\right)}^2\left(\left(3\text{b}-2\right)Y_0-1\right)}{\text{b}},-\frac{{\left(Y_0-1\right)}^2\left(\left(3\text{b}-2\right)Y_0-1\right)}{\text{b}}$	Attractor for
			, or
			.
			Source for
			, or
			.
$L_1\left(1/\text{h}\right)$	$(1,0,\frac{1}{\text{h}})$	$0,-\frac{{\left(\text{h}-1\right)}^2\left(3\text{b}-\text{h}-2\right)}{{\text{b}\text{h}}^3},-\frac{{\left(\text{h}-1\right)}^2\left(3\text{b}-\text{h}-2\right)}{{\text{b}\text{h}}^3}$.	Nonhyperbolic.
			2D unstable manifold for
			.
			2D stable manifold for
			,
			or $\text{h}>1,\text{b}>\frac{\text{h}+2}{3}$.
$L_2$	$(\frac{3}{4},0,Y_0)$ $\text{b}=1$	$0,-\frac{1}{2}{\left(Y_0-1\right)}^3,\frac{3}{4}{\left(Y_0-1\right)}^3$	Saddle.
$L_2\left(0\right)$	$(\frac{3\text{b}}{1+3\text{b}},0,0)$	$0,-\frac{3}{3\text{b}+1},\frac{2}{3\text{b}+1}$	Saddle.
$L_2\left(1\right)$	$(\frac{3\text{b}}{1+3\text{b}},0,1)$	$0,0,0$	Nonhyperbolic.
$L_3$	$(U_0,0,1)$	$0,0,0$	Nonhyperbolic.
$L_4$	$(0,0,Y_0)$	$0,-3{\left(Y_0-1\right)}^3,{\left(Y_0-1\right)}^3$	Nonhyperbolic for $Y_0=0$,
			saddle otherwise.
$L_5\left(1\right)$	$(0,1,1)$	$0,0,0$	Nonhyperbolic.
$L_6$	$(U_0,1,0)$	$0,-\frac{2\left(\sqrt{\text{b}}+1\right)\left(U_0-1\right)}{\text{b}-1},-\frac{2\left(\sqrt{\text{b}}+1\right)\left(U_0-1\right)}{\text{b}-1}.$	Nonhyperbolic.
			2D stable manifold for
			, or
			.
			2D unstable manifold for
			.
$M_1$	$(0,\frac{1}{2},0)$	$\frac{3}{2},-\frac{1}{2},\frac{1}{2}$	Saddle.
$M_2$	$(0,1,0)$	$\frac{2\left(1+\sqrt{\text{b}}\right)}{\text{b}-1},0,\frac{2\left(1+\sqrt{\text{b}}\right)}{\text{b}-1}$	Nonhyperbolic.
			2D stable manifold for
			, or .
$M_3$	$(1,1,1)$	$0,0,0$	Nonhyperbolic.
$M_5$	$(1,1,0)$	$0,0,0$	Saddle.
$M_6$	$(0,\frac{\text{h}-1}{\text{h}+1},\frac{1}{\text{h}})$	$0,0,\frac{4{\left(\text{h}-1\right)}^3}{{\text{h}}^3\left(\text{h}+1\right)}$	Nonhyperbolic.
			1D unstable manifold.
$M_7$	$(1,\frac{\text{h}-1}{\text{h}+1},\frac{1}{\text{h}})$	$0,0,\frac{2{\left(\text{h}-1\right)}^2\left(4\text{b}-\text{h}-2\right)}{{\text{b}\text{h}}^3\left(\text{h}+1\right)}$.	Nonhyperbolic.
			1D unstable manifold for
			, or
			.
			1D stable manifold for
			, or
			$\text{b}>\frac{3}{4},\text{h}>4\text{b}-2$.
$M_9$	$(1,1,Y_0)$	$0,\frac{2Y_0\left(\left(2-4\text{b}\right)Y_0+1\right)}{\text{b}},-\frac{2Y_0\left(\left(2-4\text{b}\right)Y_0+1\right)}{\text{b}}$	Saddle.
$M_{10}$	$(1,V_0,\frac{1-V_0}{1+V_0})$	$0,0,\frac{2{\left(Y_0-1\right)}^2{Y}_0\left(\left(4\text{b}-2\right)Y_0-1\right)}{\text{b}\left(Y_0+1\right)}$	Nonhyperbolic.
$M_{11}$	$(1,V_0,\frac{1}{2\left(2\text{b}-1\right)})$, $\text{b}>\frac{3}{4}$	$0,0,0$	Nonhyperbolic.
$M_{12}$	$(U_0,\frac{\text{h}-1}{\text{h}+1},\frac{1}{\text{h}})$	$0,0,-\frac{2{\left(\text{h}-1\right)}^2\left(U_0\left(2\text{b}\left(\text{h}-3\right)+\text{h}+2\right)-2\text{b}\left(\text{h}-1\right)\right)}{{\text{b}\text{h}}^3\left(\text{h}+1\right)}$	Nonhyperbolic.

Table 4. Table 4: Equilibrium points and the corresponding eigenvalues and stability conditions of the dynamical system ( 93 ).

Labels	$(S_1,S_2,Y)$	Eigenvalues	Stability
$L_1$	$({\Sigma }_c,1-\text{b}{\Sigma }_c^2,0)$ for ${\Sigma }_c\ne 1$	$0,-\frac{{\Sigma }_c}{1+n},-{\Sigma }_c$	Nonhyperbolic
$L_2$	$(0,1,y_c)$	$0,-y_c,2y_c$	saddle
$P_1$	$(\frac{2}{3},\frac{1}{9}\left(9-8\text{b}\right),1)$	$\begin{array}{c}\frac{2}{3\left(n+1\right)},\\ -\frac{1}{6}\left(3-\sqrt{64\text{b}-63}\right),\\ -\frac{1}{6}\left(3+\sqrt{64\text{b}-63}\right)\end{array}$	$\begin{array}{c}\text{Nonhyperbolic for}\text{b}=\frac{9}{8}\\ \text{saddle otherwise}\end{array}$
$P_2$	$(\frac{1}{\sqrt{\text{b}}},0,1)$	$\frac{1}{\sqrt{\text{b}}\left(n+1\right)},2,\frac{4\sqrt{\text{b}}-3}{\sqrt{\text{b}}}$	$\begin{array}{c}\text{nonhyperbolic for}\text{b}=\frac{9}{16}\\ \text{unstable for}n>0,\text{b}>\frac{9}{16}\\ \text{saddle, otherwise}\end{array}$
$P_3$	$(-\frac{1}{\sqrt{\text{b}}},0,1)$	$-\frac{1}{\sqrt{\text{b}}\left(n+1\right)},2,\frac{4\sqrt{\text{b}}+3}{\sqrt{\text{b}}}$	saddle
$P_4$	$(\frac{3}{4\text{b}},0,1)$	$\frac{3}{4\text{b}\left(n+1\right)},-\frac{8\text{b}-9}{4\text{b}},-\frac{16\text{b}-9}{8\text{b}}$

Table 5. Table 5: Equilibrium points and the corresponding eigenvalues of the dynamical system ( 109 ). We use the notations Δ = n ​ ( 5 ​ n + 4 ) − n 2 ​ n Δ 𝑛 5 𝑛 4 𝑛 2 𝑛 \Delta=\frac{\sqrt{n(5n+4)}-n}{2n} , n 0 = − ( 3 − 4 ​ b ) 2 16 ​ b 2 − 28 ​ b + 11 subscript 𝑛 0 superscript 3 4 b 2 16 superscript b 2 28 b 11 n_{0}=-\frac{(3-4\texttt{b})^{2}}{16\texttt{b}^{2}-28\texttt{b}+11} . λ c subscript 𝜆 𝑐 \lambda_{c} is a function of n 𝑛 n and b .

Labels	$(U,V,Y)$	Eigenvalues	Stability
$L_1$	$(1,0,Y_0)$	$0,-\frac{{\left(Y_0-1\right)}^2\left(\left(3\text{b}-2\right)Y_0-1\right)}{\text{b}},-\frac{{\left(Y_0-1\right)}^2\left(\left(3\text{b}-2\right)Y_0-1\right)}{\text{b}}$	Attractor for
			, or
			.
			Source for
			, or
			.
$L_2$	$(\frac{3\text{b}}{3\text{b}+1},0,Y_0)$	$0,\frac{\left(Y_0-1\right)\left(\left(9\text{b}-7\right)Y_0-2\right)}{3\text{b}+1},\frac{3\left(Y_0-1\right)\left(3\text{b}+\left(9\left(\text{b}-2\right)\text{b}+5\right)Y_0+1\right)}{{\left(3\text{b}+1\right)}^2}$	Nonhyperbolic
			2D stable manifold for
			, or
			.
			2D unstable manifold for
$L_3$	$(U_0,0,1)$	$0,0,0$	Nonhyperbolic.
$L_4$	$(0,0,Y_0)$	$0,-{\left(Y_0-1\right)}^2,3{\left(Y_0-1\right)}^2$	Saddle.
$L_5$	$(0,V_0,1)$	$0,0,0$	Nonhyperbolic.
$L_6$	$(U_0,1,0)$, $n=0$	$0,-\frac{2\left(\sqrt{\text{b}}+1\right)\left(U_0-1\right)}{\text{b}-1},-\frac{2\left(\sqrt{\text{b}}+1\right)\left(U_0-1\right)}{\text{b}-1}.$	Nonhyperbolic.
			2D stable manifold for
			, or
			.
			2D unstable manifold for
			.
$P_1$	$(0,\frac{n+1}{n+2},0)$	$-\frac{1}{n+2},\frac{1}{n+2},-\frac{n-3}{n+2}$	Saddle.
$P_2$	$(0,1,0)$	$-\frac{1}{n+1},-\frac{1}{n+1},-\frac{n}{n+1}$	Sink.
$P_3$	$(1,1,1)$	$-\frac{1}{n+1},\frac{1}{n+1},\frac{n}{n+1}$	Saddle.
$P_4$	$(\frac{\text{b}\left(n-3\right)}{\text{b}\left(n-3\right)+n-1},\frac{2\left(n+1\right)}{3n+1},0)$	$\frac{2-2n}{\left(3n+1\right)\left(\text{b}\left(n-3\right)+n-1\right)}$,
		$-\frac{\left(n-1\right)\left(n+\sqrt{1-n\left(7n-22\right)}-5\right)}{2\left(3n+1\right)\left(\text{b}\left(n-3\right)+n-1\right)}$,
		$\frac{\left(n-1\right)\left(-n+\sqrt{1-n\left(7n-22\right)}+5\right)}{2\left(3n+1\right)\left(\text{b}\left(n-3\right)+n-1\right)}$	Stable spiral for $\text{b}\ge 0,n>5$.
$P_5$	$(1,1,0)$	$0,0,0$	Nonhyperbolic.
$P_6$	$(0,\frac{1}{2\Delta +1},\frac{\Delta }{\Delta +1})$	$0,0,\frac{4\Delta }{{\left(\Delta +1\right)}^2\left(2\Delta +1\right)}$	Nonhyperbolic.
			1D unstable manifold.
$P_7$	$(1,\frac{1}{2\Delta +1},\frac{\Delta }{\Delta +1})$	$0,0,\frac{2\Delta \left(4\text{b}\Delta -3\Delta -1\right)}{\text{b}{\left(\Delta +1\right)}^2\left(2\Delta +1\right)}$.	Nonhyperbolic.
			1D unstable manifold for
			, or
			, or
			$\text{b}\ge \frac{1}{8}\left(7+\sqrt{5}\right),n>0$
			1D stable manifold for
			, or
			.
$P_8$	$(\frac{1}{2},\frac{2}{{\left(\Delta +1\right)}^2},1-\frac{4}{\Delta \left(\Delta +2\right)+3})$	$0,0,0$	Nonhyperbolic.
$P_9$	$(1,1,Y_0)$	$0,\frac{2Y_0\left(\left(2-4\text{b}\right)Y_0+1\right)}{\text{b}},-\frac{2Y_0\left(\left(2-4\text{b}\right)Y_0+1\right)}{\text{b}}$	Saddle.
$P_{10}\left(0\right)$	$(0,V_0,\frac{1-V_0}{1+V_0})$	$0,0,-\frac{8\left(V_0-1\right)V_0^{}{}_{}{}^2}{{\left(V_0+1\right)}^2}$	Nonhyperbolic.
			1D unstable manifold.
$P_{10}\left(1\right)$	$(1,V_0,\frac{1-V_0}{1+V_0})$	$0,0,\frac{2\left(V_0-1\right)V_0\left(4\text{b}\left(V_0-1\right)-V_0+3\right)}{\text{b}{\left(V_0+1\right)}^2}$	Nonhyperbolic.
$P_{11}$	$(1,V_0,\frac{1}{2\left(2\text{b}-1\right)})$, $\text{b}>\frac{3}{4}$	$0,0,0$	Nonhyperbolic.
$P_{12}$	$(U_0,\frac{1}{2\Delta +1},\frac{\Delta }{\Delta +1})$	$0,0,{\lambda }_c$	Nonhyperbolic.

Table 6. Table 6: Stability of the equilibrium points of the dynamical system ( 131 ).

Labels	Stability
$L_1$	Nonhyperbolic. 2D unstable (resp. stable) manifold for (resp. $S_1^{\ast }>0$).
$L_2^{\epsilon }$	$L_2^{+}$ (resp. $L_2^{-}$ ) nonhyperbolic. 3D unstable (resp. stable) manifold for .
$L_3^{\epsilon }$	$L_3^{+}$ (resp. $L_3^{-}$) is nonhyperbolic. 3D unstable (resp. stable) manifold for $Y_0>0,k>-2\sqrt{2}$.
$L_4^{\epsilon }$	saddle
$L_5^{\epsilon }$	$L_5^{+}$ (resp. $L_5^{-}$) is nonhyperbolic with a 4D unstable (resp. stable) manifold for
	$k\in ℝ,X_{\varphi }^{\ast }=0,\text{h}>1,\text{b}>\frac{1}{16}{\left(\text{h}+2\right)}^2$, or
	, or
	, or
	, or
	, or
	, or
	.
$L_6^{\epsilon }$	saddle.
$P_1^{\epsilon }$	nonhyperbolic.
$P_2^{\epsilon }$	saddle.
$P_3^{\epsilon }$	$P_3^{+}$ (resp. $P_3^{-}$) is nonhyperbolic with a 4D unstable (resp. stable) manifold for $\text{h}>1,\text{b}>\frac{1}{16}{\left(\text{h}+2\right)}^2$
$P_4^{\epsilon }$	saddle.
$P_5^{\epsilon }$	$P_5^{+}$ (resp. $P_5^{-}$) is a sink (resp. a source) for . They are saddle otherwise.
$P_6^{\epsilon }$	saddle.
$P_7^{\epsilon }$	$P_7^{+}$ (resp. $P_7^{-}$) is a source (resp. sink) for
	, or
	, or
	, or
	.
$P_8^{\epsilon }$	saddle.
$P_9^{\epsilon }$	For some choices of parameters, say $\text{b}=1,\text{h}=1$, or $\text{b}=1,\text{h}=2$, they are saddle.

Table 7. Table 7: Stability of the equilibrium points of the dynamical system ( 131 ) (cont.).

Labels	Stability
$P_{10}^{\epsilon }$	$P_{10}^{+}$ (resp. $P_{10}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for
	, or
	, or
	, or
	.
	$P_{10}^{+}$ (resp. $P_{10}^{-}$) nonhyperbolic with a 3D stable (resp. unstable) manifold for
	, or
	, or
	, or
	, or
	, or
	, or
	, or
	, or
	.
$P_{11}^{\epsilon }$	$P_{11}^{+}$ (resp. $P_{11}^{-}$) is nonhyperbolic with a 3D stable (resp. unstable) manifold for
	, or
	, or
	, or
	, or
	, or
	.
	$P_{11}^{+}$ (resp. $P_{11}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for
	, or
	, or
	, or
	, or
	.

Table 8. Table 8: Stability of the equilibrium points of the dynamical system ( 141 ).

Labels	Stability
$L_1^{\epsilon }$	Nonhyperbolic with a 2D unstable manifold (resp. stable) for (resp. $S_1^{\ast }>0$).
$L_2^{\epsilon }$	$L_2^{+}$ (resp. $L_2^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for .
$L_3^{\epsilon }$	$L_3^{+}$ (resp. $L_3^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for $Y_0>0,k>-2\sqrt{2}$.
$L_4^{\epsilon }$	saddle.
$L_5^{\epsilon }$	$L_5^{+}$ (resp. $L_5^{-}$) is nonhyperbolic, with a 4D unstable (resp. stable) manifold for
	$k\in ℝ,X_{\varphi }^{\ast }=0,\text{b}>\frac{9}{16},n>0$, or
	, or
	, or
	, or
	.
$L_6^{\epsilon }$	saddle.
$P_1^{\epsilon }$	Nonhyperbolic.
$P_2^{\epsilon }$	saddle.
$P_3^{\epsilon }$	saddle.
$P_4^{\epsilon }$	$P_4^{+}$ (resp. $P_4^{-}$) is nonhyperbolic with a 4D unstable (resp. stable) manifold for $n>0,\text{b}>\frac{9}{16}$.
$P_5^{\epsilon }$	$P_5^{+}$ (resp. $P_5^{-}$) is a sink (resp. a source) for or a saddle otherwise.
$P_6^{\epsilon }$	saddle.
$P_7^{\epsilon }$	$P_7^{+}$ (resp. $P_7^{-}$) is a source (resp. a sink) for
	$k\le -2\sqrt{2},\text{b}>\frac{1}{k^2},n>0$, or
	, or
	, or
	$k\ge 2\sqrt{2},\text{b}>\frac{1}{k^2},n>0$.
	It is a saddle otherwise.
$P_8^{\epsilon }$	saddle.
$P_9^{\epsilon }$	For $\text{b}=1$ they are saddle.

Table 9. Table 9: Stability of the equilibrium points of the dynamical system ( 141 ) (cont.).

Labels	Stability
$P_{10}^{\epsilon }$	$P_{10}^{+}$ (resp. $P_{10}^{-}$) has a 3D stable (resp. unstable) manifold for
	, or
	, or
	, or
	, or
	, or
	, or
	, or
	, or
	, or
	.
	$P_{10}^{+}$ (resp. $P_{10}^{-}$) has a 3D unstable (resp. stable) manifold for
	, or
	, or
	$k\ge 2\sqrt{2},n\ge 0,\text{b}>\frac{9k^2-8}{16k^2}$.
$P_{11}^{\epsilon }$	$P_{11}^{+}$ (resp. $P_{11}^{-}$) has a 3D stable (resp. unstable) manifold for
	$k\le -2\sqrt{2},n\ge 0,\text{b}>\frac{9k^2-8}{16k^2}$, or
	, or
	.
	$P_{11}^{+}$ (resp. $P_{11}^{-}$) has a 3D unstable (resp. stable) manifold for
	, or
	, or
	, or
	, or
	, or
	.

Equations

$$y^{\prime \prime }+3y{y}^{\prime }+y^3=0.$$

$$y_{0}p\left(p-1\right){\left(r-r_0\right)}^{p-2}+3{\left(y_0\right)}^{2}p{\left(r-r_0\right)}^{2p-1}+a_0{\left(r-r_0\right)}^{3p}=0,$$

$$y\left(r\right)={\left(r-r_0\right)}^{-1}+y_1+y_2{\left(r-r_0\right)}^1+y_3{\left(r-r_0\right)}^2+...$$

$$y_2=-y_1^2,y_3=y_1^3,...,y_N={\left(-1\right)}^{N+1}{\left(y_1\right)}^N.$$

$$p^{\prime }(\rho )=-\frac{\left(p{\rho }^3+2m\right)(p+\mu )}{2\rho (\rho -2m)},$$

$$m^{\prime }(\rho )=\frac{\mu {\rho }^2}{2},$$

$${\psi }^{\prime }(\rho )=\frac{\left(p{\rho }^3+2m\right)}{2\rho (\rho -2m)},$$

$$\mathcal{M}=\frac{m}{\rho },D=\frac{\mu {\rho }^2}{2},$$

$$\frac{dD}{d\tau }=D\left(\frac{n+5}{2}+\frac{(n+1)D}{2nx}+\frac{n+1}{4x}\right),$$

$$\frac{dx}{d\tau }=-x+D-\frac{1}{2}.$$

$$u=\frac{D}{x},v=\frac{1}{x},$$

$$\frac{du}{d\tau }=u\left(\frac{(1-n)u}{2n}+\frac{(n+3)v}{4}+\frac{n+7}{2}\right),$$

$$\frac{dv}{d\tau }=v\left(-u+\frac{1}{2}v+1\right).$$

$$p=a(\Gamma -1){\mu }^{\Gamma },\Gamma =1+\frac{1}{n},$$

$$r^2=\frac{2}{k}\times \frac{\mathcal{M}}{D},m^2=\frac{2}{k}\times \frac{{\mathcal{M}}^3}{D},$$

$$R_{\text{comp}}=\frac{\sqrt{k}r}{\sqrt{1+k{r}^2\left(1+{\mathcal{M}}^2\right)}},M_{\text{comp}}=\frac{\sqrt{k}\mathcal{M}r}{\sqrt{1+k{r}^2\left(1+{\mathcal{M}}^2\right)}},$$

$$R_{\text{comp}}=\frac{R/R_{\text{ref}}}{1+R/R_{\text{ref}}},M_{\text{comp}}=\frac{M/M_{\text{ref}}}{1+M/M_{\text{ref}}},$$

$${\Gamma }_N(p):=\frac{\mu }{p}\frac{dp}{d\mu },{\rm Y}(p):={\Gamma }_N(p{)}^{-1},\sigma (p):=\frac{p}{\mu },\Omega (p)=\frac{\omega (p)}{1+\omega (p)}$$

$$\underset{p\to 0}{\lim }{\rm Y}(p)={{\rm Y}}_0:=\frac{n_0}{1+n_0},\underset{p\to 0}{\lim }{\rm Y}(p)=1,\underset{p\to \infty }{\lim }\sigma (p)={\gamma }_1-1,$$

$$\underset{p\to 0}{\lim }\frac{d{\rm Y}}{d\ln p}=\underset{p\to \infty }{\lim }\frac{d{\rm Y}}{d\ln p}=0.$$

$$S=S_{GR}+S_u+S_m,$$

$$S_u=\int d^{4}x\sqrt{-g}\left(-K^{ab}{}_{cd}{\nabla }_a{u}^c{\nabla }_b{u}^d+\Lambda \left(u^c{u}_c+1\right)\right)$$

$$K^{ab}{}_{cd}\equiv c_1{g}^{ab}{g}_{cd}+c_2{\delta }_c^a{\delta }_d^b+c_3{\delta }_d^a{\delta }_c^b+c_4{u}^a{u}^b{g}_{cd},$$

$$G_{ab}=T_{ab}^{\text{ae}}+T_{ab}^{\text{mat}}$$

$$T_{ab}^{\text{ae}}$$

$$-2c_4{\dot{u}}_a{\dot{u}}_b+2\Lambda u_a{u}_b+g_{ab}{\mathcal{L}}_u,{J^a}_m=-{K^{ab}}_{mn}{\nabla }_b{u}^n,{\dot{u}}_a=u^b{\nabla }_b{u}_a,$$

$$\Lambda u_b$$

$$u^a{u}_a$$

$$\Lambda =-u^b{\nabla }_a{J}_b^a-c_4{\dot{u}}_a{\dot{u}}^a,$$

$$0=h^{bc}{\nabla }_a{J}_b^a+c_4{h}^{bc}{\dot{u}}_a{\nabla }_b{u}^a,$$

$$T_{ab}^{\text{mat}}\equiv \mu u_a{u}_b+p{h}_{ab},$$

$$d{s}^2=-N^2(r)d{t}^2+\frac{d{r}^2}{r^2}+K^{-1}(r)(d{\vartheta }^2+{\sin }^2\vartheta d{\phi }^2),$$

$$r\frac{dx}{dr}=\frac{\mu +3p}{2\text{b}}+2(\text{b}-1)y^2+3xy+K,$$

Full text

1. 1 Introduction

2. 1.1 Spherically symmetric models

3. 1.2 Background

4. 1.3 The Ablowitz, Ramani and Segur algorithm

5. 1.4 Tolman-Oppenheimer-Volkoff equation for relativistic star models

6. 1.5 Integrability and modified Tolman-Oppenheimer-Volkoff formulation

7. 2 Einstein-æther Gravity

8. 3 Static spherically symmetric spacetime with a perfect fluid

9. 3.1 Model 1: Linear equation of state

10. 3.1.1 Singularity analysis and integrability

11. 3.1.2 Equilibrium points in the finite region of the phase space

12. 3.1.3 Dynamical systems analysis based on the Newtonian homology invariants

13. 3.1.4 Physical discussion

14. 3.2 Model 2: Polytropic equation of state

15. 3.2.1 Singularity analysis and integrability

16. 3.2.2 Equilibrium points in the finite region of the phase space

17. 3.2.3 Dynamical systems analysis based on the Newtonian homology invariants

18. 3.2.4 Physical discussion

19. 4 Stationary comoving æther with perfect fluid and scalar field in static metric

20. 4.1 Model 3: Perfect fluid with linear equation of state and a scalar field with an exponential potential

21. 4.1.1 Singularity analysis and integrability

22. 4.1.2 Equilibrium points in the finite region of the phase space

23. 4.1.3 Physical discussion

24. 4.2 Model 4: Perfect fluid with polytropic of state and a scalar field with an exponential potential

25. 4.2.1 Singularity analysis and integrability

26. 4.2.2 Equilibrium points in the finite region of the phase space

27. 4.2.3 Physical discussion

28. 5 Conclusions

29. A Equilibrium points for Model 1

30. A.1 Equilibrium points of Model 1 at the finite part of the phase space $(S_{1},S_{2},Y)$

31. A.2 Equilibrium points of Model 1 at infinity.

32. A.3 Equilibrium points of Model 1 using homology invariant variables $(U,V,Y)$ in the case of General Relativity ($\texttt{b}=1$).

33. A.4 Equilibrium points of Model 1 using homology invariant variables $(U,V,Y)$ for the Einstein- æther modification ($\texttt{b}\neq 1$).

34. B Equilibrium points for Model 2

35. B.1 Equilibrium points of Model 2 at the finite part of the phase space $(S_{1},S_{2},Y)$

36. B.2 Equilibrium points of Model 2 at infinity

37. B.3 Equilibrium points of Model 2 using homology invariant variables $(U,V,Y)$ in the case of General Relativity ($\texttt{b}=1$)

38. B.4 Equilibrium points of Model 2 using homology invariant variables $(U,V,Y)$ for the Einstein- æther modification ($\texttt{b}\neq 1$)

39. C Equilibrium points for Model 3

40. D Equilibrium points for Model 4

1 Introduction

Einstein-æther theory [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] is an effective field theory preserving locality and covariance, and consists of General Relativity (GR) dynamically coupled to a time-like unit vector field (the æther), in which the local spacetime structure is determined by both the dynamical æther vector field and the geometrical metric tensor field [7]. Since a preferred frame is specified by the æther at each spacetime point, Lorentz invariance is broken spontaneously, more precisely it breaks the invariance under boosts. We note that every hypersurface-orthogonal Einstein-æther solution is a solution of the IR limit of (extended) Hořava gravity [13]; the relationship between Einstein-æther theory and Hořava-Lifschitz (HL) gravity is further clarified in [14].

Exact solutions and the qualitative analysis of Einstein-æther cosmological models have been presented by a number of authors, [15, 16, 17, 18, 19], with an emphasis on the impact of Lorentz violation on the inflationary scenario [6, 4, 5]. More recent developments of the HL theory were reviewed in [20], including universal horizons, black holes and their thermodynamics, gauge-gravity duality, and the possible quantization of the theory. HL cosmological scenarios were tested against new observational constraints in [21] using updated Cosmic Microwave Background, Type Ia supernovae, and Baryon Acoustic Oscillations cosmological data. It is known that cosmologically viable extended Einstein-æther theories are compatible with cosmological data [22].

1.1 Spherically symmetric models

The Einstein-æther field equations (FE) are derived from the Einstein-æther action [1, 3]. In such a model there will be additional terms in the FE due to the effects of the curvature of the underlying spherically symmetric geometry and from the stress tensor, $T_{ab}^{\ae}$, for the æther field (which depends on a number of dimensionless parameters $c_{i}$). In addition, because a FLRW limit is not required by the analysis, the æther field can be “tilted” with respect to the (perfect fluid) CMB rest frame, which provides additional terms in the matter stress tensor, $T_{ab}^{mat}$ (characterized by a “hyperbolic tilt angle”, $\alpha$, measuring the relative boost). It is anticipated that this tilt will decay to the future [23]. For example, (and consistent with the local linear stability analysis presented in [4]), $\alpha\rightarrow 0$ at late times in a tilted (æther) Bianchi I cosmological model in which the initial hyperbolic tilt angle $\alpha$ and its time derivative are sufficiently small [6].

The Einstein-æther parameters $c_{i}$ are dimensionless constants. In the case of spherical symmetry the æther is hypersurface orthogonal, so that the twist vanishes and $c_{4}$ can be set to zero without any loss of generality [7]. The spherically symmetric solutions of Einstein-æther theory are a subset of the solutions of Hořava gravity. In general the converse is not true, except when the solutions have a regular center [7]. The $c_{i}$ parameters also determine the effective Newtonian gravitational constant $G$ in the model, and consequently $8\pi G$ can be set to unity by an appropriate renormalization. The spherically symmetric model can thus be characterized by two non-trivial constant parameters. In GR $c_{i}=0$. In the qualitative analysis that follows we shall assume values for the non-GR parameters $c_{i}$ compatible with all current constraints [7, 24]. In the case of a static vacuum there is a 3-parameter family of spherically symmetric solutions. In particular, when the æther is parallel to the Killing vector, which is referred as a “static æther”, a static vacuum solution is explicitly known [25].

Spherically symmetric static and stationary Einstein-æther solutions in which there is an additional (to GR) radial tilting æther mode are of physical importance, and a number of black hole and time-independent solutions are known [7, 24, 25, 26]. In general, the dynamics of the cosmological scale factor and perturbations in non-rotating neutron star and black hole solutions are very similar to those in GR. Although a fully nonlinear positive energy result has been demonstrated for spherically symmetric solutions at an instant of time symmetry [9], a thorough study of the fully nonlinear solutions has not yet been carried out.

1.2 Background

Recently spherically symmetric Einstein-æther cosmological models with a perfect fluid matter source have been studied [27], including FRW models [28], Kantowski-Sachs models [29, 19], and spatially homogenous metrics [30]. In all cases the matter source was assumed to include a scalar field which is coupled to the expansion of the æther field through a generalized exponential potential.

In order to perform a dynamical systems analysis [31] it is useful to introduce suitable normalized variables [32, 33, 34] to simplify the FEs, which also facilities their numerical study. In particular, we derive the equilibrium points of the algebraic-differential system by introducing proper normalized variables [18]. For every point the stability and the physical variables can be given.

The static spherically symmetric model with perfect fluid, first introduced in Section 6.1 of [18], has been investigated in [35] using appropriate dynamical variables inspired by [36]. In the reference [35] two of us have presented asymptotic expansions for all of the equilibrium points in the finite region of the phase space. We have shown that the Minkowski spacetime can be given in explicit spherically symmetric form [36] irrespective of the æther parameters. We have shown that we can have non-regular self-similar perfect fluid solutions (like those in [37, 38, 39]), self-similar plane-symmetric perfect fluid models, and Kasner’s plane-symmetric vacuum solutions [40]. We discussed the existence of new solutions related with naked singularities or with horizons. The line elements have been presented in explicit form. We also discussed the dynamics at infinity and presented some numerical results supporting our analytical results. Scalar field configurations in static spherically symmetric metrics were also included. The physical consequences of the models were discussed in [36, 41, 42], most of which are astrophysical rather than cosmological ones. Relativistic polytropic EoS were examined in [43] in Hořava gravity and in Einstein-æther theories with anisotropic fluid. Assuming viable analytic and exact interior solutions, there was found the equation of state (EoS) by a polynomial fit, not accurate enough to describe the relation between the density and the radial and tangential pressure, so that a correction of the relativistic polytropic EoS was needed to give a satisfactory fit.

The present research is the follow-up of reference [35]. Here we continue the investigation of static spherically symmetric Einstein-æther models. Particularly, we investigate the integrability of the FE and we present the Modified Tolman-Oppenheimer-Volkoff approach.

1.3 The Ablowitz, Ramani and Segur algorithm

The modern treatment of singularity analysis is summarized in the ARS algorithm (from the initials of Ablowitz, Ramani and Segur). There are three main steps to the algorithm, they are: (i) to determine the movable singularity, which should be a single pole, (ii) to determine the resonances which consist of the position of the integration constants in the Laurent expansion and (iii) perform the consistency test.

In order to demonstrate the application of the singularity analysis we consider the well-known Painlevè-Ince equation [44]:

[TABLE]

We substitute $y\left(r\right)=y_{0}\left(r-r_{0}\right)^{p}$ in the latter equation and we find the polynomial equation

[TABLE]

for which balance occurs if $p=-1$ and consequence $y_{0}=1$ or $y_{0}=2$. The requirement that there is a balance in the exponents is necessary in order for the function $y\left(r\right)=y_{0}\left(r-r_{0}\right)^{p}$ describes the behavior of the solution near the singularity. As we can see there are two coefficient constants $y_{0}$ which mean that there are two possible different solutions.

For the second step of the ARS algorithm we take $y\left(r\right)=y_{0}\left(r-r_{0}\right)^{-1}+\varepsilon\left(r-r_{0}\right)^{-1+s}$

and we linearize around $\varepsilon=0$. We equate the remaining leading-order terms to zero and we find the values of the resonances, where for $y_{0}=1$ they are $s_{1}=-1$ and $s_{2}=1,$ while for $y_{0}=2$ the resonances are $s_{1}=-1$ and $s_{2}=-2$. Because the resonances are rational numbers, we shall say that the Painlevè-Ince equation passes the second step  of the ARS algorithm.

The final step is to write the solution in terms of a Laurent expansion, and for $y_{0}=1$, the solution is given by the following partial Laurent expansion

[TABLE]

where $y_{1}$ is arbitrary, and denotes an integration constant of the solution; the position of the singularity $r_{0}$ is considered another integration constant. By replacing the latter Laurent expansion in the Painlevè-Ince equation, we obtain the coefficients

[TABLE]

Hence, we conclude that the Painlevè-Ince possesses the Painlevè property and passes the ARS algorithm.

The study of whether a nonlinear dynamical system posses the integrability property is essential in all physical theories. Indeed, by knowing that a given dynamical system is integrable we know that numerical trajectories/solutions describe real solutions of the dynamical system, while small changes in the initial conditions do not affect the main evolution of the system.

Usually we refer to a dynamical system as integrable when there exists a sufficient number of invariant functions which can be used to reduce the system into an algebraic equation and when it is thus feasible to write the explicit solution in terms of closed-form functions. However, in the context of the singularity analysis the concept of integrability is strict; that is, the solution of the dynamical system must be analytic apart from isolated movable pole-like singularities (or branch point singularities) so that the complex plane of the independent variable can be divided into sections and in each section the solution is analytic. Because of this latter criteria, the free parameters of the theory are constrained in order for the dynamical system to possess the Painlevè property. This is essentially only a mathematical criterion and does not necessarily provide physical properties of the system.

For other values of the free parameters, where the dynamical system does not pass the singularity analysis, other methods for the determination of the conservation laws should be used in order to analyze the integrability of the system.

1.4 Tolman-Oppenheimer-Volkoff equation for relativistic star models

The Tolman-Oppenheimer-Volkoff approach for a relativistic star model is based on the equations (in units where $8\pi G=1$):

[TABLE]

where, $\mu(\rho)$, $p(\rho)$, $m(\rho)$ are the density, the pressure, and the mass enclosed in a sphere of radius $\rho$. Equation (5a) is the well-known Tolman-Oppenheimer-Volkoff equation. $\psi(\rho)$ is the gravitational potential. The last equation is decoupled, and we obtain a reduced 2D system of differential equations provided an EoS $\mu=\mu(p)$ is given.

Now we discuss some particular models. Following the reference [45], a model with the isothermal equation of state $\mu=np$ was studied. By utilizing the new variables

[TABLE]

where $\mathcal{M}$ denotes the Misner-Sharp mass [46], and introducing $x=\mathcal{M}-\frac{1}{2}$, and the new radial variable $\tau=\ln\rho$ the field equations can be cast as [45]:

[TABLE]

Defining

[TABLE]

the Lotka-Volterra model is obtained

[TABLE]

In [45] it was proven that the above system is Liouville-integrable if and only if $n\in\{-1,-3,-5,-6\}$, where by Liouville-integrable is it understood that a system of polynomial differential equations have first integrals given by elementary functions or integrals of elementary functions (that is, functions expressed in terms of combinations of exponential functions, trigonometric functions, logarithmic functions or polynomial functions). In the special cases $n=-1$ (cosmological de Sitter solution) and $n=-3$ , the resulting system can be expressed as Abel equations of second order that are exactly solvable.

Another important EoS model leads to ideal polytropes stars [47]:

[TABLE]

where $\Gamma$ is the so-called polytropic index. Is it well-known that the value $n=5$ is a bifurcation value that separates two regimes [47]:

$5<n<\infty$: corresponding to “gas spheres” of infinite extension;

• 2.

$0<n<5$: corresponding to finite stars with boundary.

Therefore, it is convenient to consider variables that can distinguish between configurations with infinite extensions from that of finite extension and, similarly, that distinguish between finite and infinite masses.

For fluid spheres of mass $m$ and radius $r$ we have the universal scaling ($m$-$r$-relation) [48]:

[TABLE]

where $k$ is a constant depending on the EoS, where it is presumed that $r$ and $m$ are finite. This motivates the following definition of compact variables:

[TABLE]

which are both finite when $r$ and $m$ are finite, and $R_{\text{comp}}\rightarrow 1$ as $r\rightarrow\infty$, and $M_{\text{comp}}\rightarrow 1$ as $m\rightarrow\infty$. Furthermore, the circle $R_{\text{comp}}^{2}+M_{\text{comp}}^{2}=1$ corresponds to any solution with $r\rightarrow\infty$, and $m\rightarrow\infty$. These variables are closely related (and have an analogous physical interpretation) to the variables used in [48]:

[TABLE]

where $R_{\text{ref}}$ and $M_{\text{ref}}$ are typical values of radii and mass. That is $R_{\text{comp}}\rightarrow 1$ as $r\rightarrow\infty$, and $M_{\text{comp}}\rightarrow 1$ as $m\rightarrow\infty$. The only difference is that in the second formulation the case of infinite mass and infinite radius corresponds to the point $(1,1)$ in the $M_{\text{comp}}$ vs $R_{\text{comp}}$ diagram, instead of the whole or part of the circle $R_{\text{comp}}^{2}+M_{\text{comp}}^{2}=1$ in our formulation (but they share the same qualitative features).

In reference [48]some relevant mass-radius theorems for relativistic spherically symmetric static perfect fluid models were proved. In particular, a general class of asymptotically polytropic EoS, where the relevant parameters $n_{0}$ and $\Gamma_{N}$ are defined as follows, was considered:

[TABLE]

where $p$ and $\mu$ are the density and the pressure of the fluid and $\omega(p)$, the so-called pressure variable, is a dimensionless continuous function for $p\in[0,\infty)$, strictly monotonically increasing and sufficiently smooth on $(0,\infty)$ with $\omega(0)=0$ and $\lim_{p\rightarrow\infty}\omega(p)=\infty$ [48]. Was also considered an EoS such that $\Upsilon(p)$ is a continuous function for $p\in[0,\infty)$, with continuous first order derivative for $p\in(0,\infty)$, and such that

[TABLE]

Here $\Upsilon_{0}$, and $n_{0}$ are constants that satisfy $0\leq\Upsilon_{0}<1\quad(0\leq n_{0}<\infty)$. Additionally, it is assumed

[TABLE]

The relevant Mass-radius theorem for relativistic spherically symmetric static perfect fluid models states (Theorems: 5.1, 5.2, 5.3, 5.4 of [48]):

All regular solutions have infinite masses and infinite radii if $\Gamma_{N}\leq\frac{6}{5}$ and $\sigma<1$.

• 2.

All regular solutions with $\Omega_{c}\in(0,\Omega_{\text{max}})$ have finite mass and radii if $\Gamma_{N}>\frac{6}{5}$ and $\sigma<\frac{1}{6}(-4+\sqrt{16+2(5\Gamma_{N}-6)})$.

• 3.

All regular and non-regular perfect fluid solutions have finite radii and masses if $n_{0}\leq 3$ (i.e., $\Gamma_{0}\geq\frac{4}{3}$).

• 4.

For $0<n_{0}\leq 3$, and for sufficient high central pressures (i.e., when $\Omega\rightarrow 1$), the mass-radius $M$-$R$ diagram exhibits a spiral structure, with $M$ and $R$ given by the formulas (29) and (30) of [48].

As a consequence, for relativistic static spherically symmetric perfect fluid models with exact polytropic EoS, all regular solutions possess finite $R$ and $M$ for $n_{0}\equiv n\leq 3$ ($\Gamma\geq\frac{4}{3}$). All regular solutions with small central pressures have finite extent when $3<n_{0}<5$ ($\frac{6}{5}<\Gamma<\frac{4}{3}$). In the case of asymptotic polytropic EoS (in the sense discussed above), the result is that all regular solutions with $n_{0}\lesssim 3.366$ possess finite $R$ and $M$ [48].

1.5 Integrability and modified Tolman-Oppenheimer-Volkoff formulation

In order to prove that there actually exist solutions, the integrability of the system should be established [49, 50, 51, 52, 53, 54, 55]. We first study the integrability of static spherically symmetric spacetimes in Einstein- æther gravity with a matter source of the form of a perfect fluid with EoS $\mu=\mu_{0}+\left(\texttt{h}-1\right)p,~{}\texttt{h}>1$. The governing equations for such a spacetime form an algebraic-differential system of first order which can be written as system of second-order differential equations with fewer independent variables by utilizing the ARS (Ablowitz, Ramani and Segur) algorithm [56]. We first show that the system passes the singularity test, which implies that they are integrable. Additionally, we consider a polytropic EoS $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$, and we extend the analysis by adding a scalar field with an exponential potential.

Specifically, we investigate if the gravitational FE in the Einstein-æther model in the static spherically symmetric spacetime posses the Painlevè property (so that an analytic explicit integration can be performed). In order to perform such an analysis we apply the classical treatment for the singularity analysis which is summarized in the ARS algorithm. As far as the dynamical system with only the perfect fluid is concerned, we show that the FE are integrable and we write the analytic solution in terms of mixed Laurent expansions. We show that if the scalar field is not present, and assuming that the perfect fluid has the EoS $\mu=\mu_{0}+\left(\texttt{h}-1\right)p,~{}\texttt{h}>1$, the system passes the singularity test, which implies that it is integrable (in a similar way to what has been done for Kantowski-Sachs Einstein-æther theories in [19]). That is, in the case of the static spherically symmetric spacetime with a perfect fluid there always exist a negative resonance which means that the Laurent expansion is expressed in a Left and Right Painlevé Series because we integrate over an annulus around the singularity which has two borders [57]. On the other hand, in the presence of a scalar field, or when the EoS of the perfect fluid is polytropic $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$, the FE do not pass the Painlevè test. Therefore, they are not integrable.

It is worth mentioning a few technicalities. First, we have to choose coordinates that can make the calculations simpler. Furthermore, since we know that the full system is, in general, non-integrable (including simpler models with realistic matter for neutron star solutions or for vacuum solutions in the context of modified gravity) we have to use techniques that do not involve actually solving the FE. The local semi-tetrad splitting [58] allows us to recast the FE into an autonomous system of covariantly defined quantities [59, 60]. This enables us to get ride of all the coordinate singularities that may appear because of badly defined coordinates. In addition, the autonomous systems can be simplified further by incorporating Killing symmetries of the spacetimes. For studying the integrability of the FE we combine the ARS method [19], Lie/Noether symmetries [61, 62] and Cartan symmetries [63]. In this research we use the ARS method.

One special application of interest is to use dynamical system tools to determine the conditions under which stable stars can form. By using the Tolman-Oppenheimer-Volkoff (TOV) approach [37, 38, 39], the relativistic TOV equations are drastically modified in Einstein-æther theory, and we can explore the physical implications of this modification. This approach to studying static stars uses the mass and pressure as dependent variables of the Schwarzschild radial coordinate. The resulting equations are however not regular at the center of the star, and a regularizing method for analyzing the solution space is required. In the relativistic and Newtonian case the resulting equations of the regularizing procedure are polynomial, and therefore they are regular. However, for the non-relativistic case the main difficulty is that the resulting system is rational. However, we are able to find some analytical results that are confirmed by numerical integrations of 2D and 3D dynamical system in compact variables, to obtain a global picture of the solution space for a linear EoS, that can visualized in a geometrical way. This study can be extended to a wide class of EoS. For example, we study the case of a polytropic EoS: $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$ in the static Einstein-æther theory, by means of the so-called $(S_{1},S_{2},Y)$-formulation given by the model (93). For $\texttt{b}=1$ the model of [41] is recovered.

The paper is organized as follows. In Section 2 we give the basic definitions of Einstein-æther gravity. In Section 3 we study static spherically symmetric spacetime with a perfect fluid with a linear EoS (Subsection 3.1) and with a polytropic EoS (Subsection 3.2). For each choice of the EoS for the perfect fluid we provide i) a singularity analysis and integrability analysis of the FE; ii) a local stability analysis of the equilibrium points in the finite region of the phase space, using the $S_{1},S_{2},Y$ formulation, that takes advantage of the local semi-tetrad splitting, iii) a local stability analysis of the equilibrium points at infinity (when $\texttt{b}<0$) using the dominant quantities to create compact variables; and finally iv) we present a dynamical systems analysis based on the Newtonian homology invariants, $(U,V,Y)$-formulation, that makes use of the TOV formulation of the FE. As we see, are obtained as a particular case the results for relativistic and Newtonian stellar models. For the TOV formulation of the Einstein-æther modification, the main difficulty is that although the non-regular nature of the $m,p$ equations at the center is overcome using Newtonian homology-like invariants, the resulting equations are rational (includes radicals) and are more difficult to analyze as compared with the local semi-tetrad splitting formulation. Finally, we investigate in Section 4 a stationary comoving æther with a perfect fluid with both linear and polytropic EoS and a scalar field with exponential potential in a static metric. We find that the models do not have the Painlevè property. Therefore, they are not integrable. We investigate the stability of the equilibrium points using the local semi-tetrad splitting formulation due to their advantages over the Newtonian homology-like invariants formulation in the Subsections 4.1.2 (for the linear EoS) and 4.2.2 (for polytropic EoS). Finally in Section 5 we discuss our results and draw our conclusions. In A and in B we present a detailed stability analysis of the equilibrium points for perfect fluid models (with linear or polytropic EoS, respectively) using the four approaches, i) - iv) previously described. In C and in D we present a detailed stability analysis of the equilibrium points when the matter content consists on a perfect fluid and a scalar field with an exponential potential.

2 Einstein-æther Gravity

The Einstein-æther Action Integral is given by [7, 8]

[TABLE]

where $S_{GR}=\int d^{4}x\sqrt{-g}\left(\frac{R}{2}\right)$ is the Einstein-Hilbert term, $S_{m}$ is the term which corresponds to the matter source and

[TABLE]

corresponds to the æther field. $\Lambda$ is a Lagrange multiplier enforcing the time-like constraint on the æther [9], for which we have introduced the coupling [7]

[TABLE]

which depends upon four dimensionless coefficients $c_{i}$. Finally, $u^{a}$ is the normalized time-like vector (observer) in which $u^{a}u_{a}=-1$.

Variation with respect to the metric tensor in (16) provides the gravitational FE

[TABLE]

in which $G_{ab}$ is the Einstein tensor, $T_{ab}^{\text{mat}}$ corresponds to $S_{m}$, and $T_{ab}^{\ae}$ is the æther tensor [64]

[TABLE]

and variation with respect to the vector field $u^{a}$ and the Lagrange multiplier gives us

[TABLE]

[TABLE]

and the compatibility conditions

[TABLE]

where $h_{ab}=(g_{ab}+u_{a}u_{b})$ is the projective tensor where $h^{ab}u_{b}=0$. We shall use the equation (22) as a definition for the Lagrange multiplier, whereas the second equation (23) leads to a set of restrictions that the æther vector must satisfy.

The energy momentum tensor of a matter source of the form of a perfect fluid (which includes a scalar field) with energy density $\mu$, and pressure $p$, in the 1+3 decomposition with respect to $u^{a}$, is

[TABLE]

where $\mu$ is the density and $p$ is the pressure of the perfect fluid, respectively.

For simplicity, in the following we introduce the constants redefinition $c_{\theta}=c_{2}+(c_{1}+c_{3})/3,c_{\sigma}=c_{1}+c_{3},\ c_{\omega}=c_{1}-c_{3},c_{a}=c_{4}-c_{1},\texttt{b}=c_{a}+1$.

3 Static spherically symmetric spacetime with a perfect fluid

We consider a static spherically symmetric spacetime with line element

[TABLE]

where we have fixed the spatial gauge to have $e_{1}{}^{1}(r)\equiv r$, and $x\equiv\frac{1}{2}r\frac{d\ln(K)}{dr},~{}y\equiv r\frac{d\ln(N)}{dr}$.

The FE are [18]:

[TABLE]

Furthermore, there exists the constraint equation

[TABLE]

From (26c) and (26d) we have that $y=-\frac{rp^{\prime}}{\mu+p}~{}\ ,~{}x=\frac{1}{2}\frac{rK^{\prime}}{K}$, where the prime ′ denotes a derivative with respect $r$. By substituting into (26a), and (26b) we find a system of two second-order ordinary differential equations,

[TABLE]

while (27) becomes

[TABLE]

In order to apply Tolman-Oppenheimer-Volkoff approach (recall that we use units where $8\pi G=1$) we use a coordinate change to write the line element as

[TABLE]

and define

[TABLE]

where $\psi$ is the gravitational potential, $m(\rho)$ denotes the mass up to the radius $\rho$, with range $0\leq\rho<\infty$, which in the relativistic case is defined by

[TABLE]

such that (in the relativistic case)

[TABLE]

Now we will calculate $m^{\prime}(\rho)$, $p^{\prime}(\rho)$ and $\psi^{\prime}(\rho)$ in the Einstein-æther theory and compare with the relativistic TOV equations. By definition

[TABLE]

Equation (34a) leads to the definition

[TABLE]

Taking derivative with respect to $\rho$ in both sides of (35) and using Eqs. (26a) and (31) we obtain

[TABLE]

On the other hand, the identity

[TABLE]

can be written as

[TABLE]

Equations (36), (38) have two solutions for $p^{\prime}(\rho)$ and $m^{\prime}(\rho)$. We consider the branch

[TABLE]

Taking the limit $\texttt{b}\rightarrow 1$ in the Eqs. (39) we obtain the proper general relativistic equations (5).

We have to consider the reality conditions:

[TABLE]

or

[TABLE]

For the analysis we have to specify the equation of state, at least on implicit form $p=p(s),\mu=\mu(s)$. One can be interested in determine conditions on which stable stars can form. The set of equations (39), determine the star’s structure and the geometry in the static spherically symmetric Einstein-æther theory for a perfect fluid. Equation (39a) is a modification of the well-known Tolman-Oppenheimer-Volkoff equation for relativistic star models.

3.1 Model 1: Linear equation of state

For the equation of state for the perfect fluid we consider that $\mu(r)=\mu_{0}+\left(\texttt{h}-1\right)p\left(r\right)$, with $\mu_{0}\geq 0,\texttt{h}>1$.

3.1.1 Singularity analysis and integrability

In this section we use the line element (25). Then, by introducing the radial rescaling $r=e^{\ell}$, such that $\ell\rightarrow-\infty$ as $r\rightarrow 0$ and $\ell\rightarrow\infty$ as $r\rightarrow\infty$, we obtain the equations

[TABLE]

while (27) becomes

[TABLE]

where ∗ means derivative with respect to $\ell$.

Next, we apply the ARS algorithm and we find that the dominant terms are111We have taken the position ${\ell}_{0}$ of the singularity to be at zero. $p_{d}=k{\ell}^{-2}~{},~{}K=K_{0}{\ell}^{-2},~{}$where

[TABLE]

For the resonances we have that

[TABLE]

where

[TABLE]

Furthermore, as for $\texttt{h}>1$ the resonance $s_{-}$ is negative, the Laurent expansion is a Left-Right Painlevé Series. Importantly, for the singularity analysis to succeed the resonances must be rational numbers.

Consider the special case $\texttt{h}=4,~{}$with $\texttt{b}=\frac{11}{2}$, which gives the resonances $s_{1}=-1~{},~{}s_{2}=3~{},~{}s_{+}=\frac{1}{2}~{},~{}s_{-}=-2.~{}$ Hence the solution is given by

[TABLE]

[TABLE]

For $\mu_{0}=0$, after we replace the solution (45), (46) in (28a) we find that the coefficient constants are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $p_{1},K_{6}$ are the two constants of integration. The third constant of integration is the position of the singularity while the fourth constant of integration appears in the next coefficients. Furthermore, from (41) we have the constraint for the constants of integration, $z_{6}\simeq 154.84\left(p_{1}\right)^{6}$.

For $\mu_{0}\neq 0$, we find the coefficients (47), (48),

[TABLE]

[TABLE]

[TABLE]

where again $p_{1},K_{6}$ are two constants of integration.

From  this analysis we performed the consistency test and we showed that the FE (40a)-(41) pass the singularity analysis and form an integrable system for values of the free parameters for which the resonances are rational numbers.

3.1.2 Equilibrium points in the finite region of the phase space

In this section we consider the line element (25), with the definition ($N=e^{\psi}$). For the dynamical systems formulation we introduce the phase space variables

[TABLE]

Assuming that the energy density $\mu$ and the pressure $p$ are both non-negative we obtain $0\leq Y\leq 1$, where we have attached the boundaries $Y=0$ and $Y=1$. On the other hand, $S_{2}\geq 0$ by definition.

The relation between the gravitational potential $\Omega$, and the matter field is given by

[TABLE]

Hence,

[TABLE]

where $\alpha$ is a freely specifiable constant corresponding to the freedom of scaling the time coordinate in the line element. This in turn reflects the freedom in specifying the value of the gravitational potential at some particular value of $r$. Matching an interior solution with the exterior Schwarzschild solution at the value where the pressure vanishes, however, fixes this constant. Using the variables $S_{1},S_{2},Y$ we obtain the dynamical system

[TABLE]

where we have defined the new radial coordinate $\lambda$ by

[TABLE]

We have the additional equation

[TABLE]

which implies that $P=0$ is an invariant subset. Thus, since $S_{2}\geq 0$, we have from (54b) that $P\geq 0$ for $0\leq Y\leq 1$. Hence, system (57) defines a flow on the (invariant subset) phase space:

[TABLE]

Finally, from the definitions of $x,y$ and $\theta$ we obtain the auxiliary equations

[TABLE]

The first equation implies that the sign of $\theta$ is invariant. Thus, we can choose $\theta>0$ and then, the direction of the radial variable $\lambda$ is well defined. The last two equations lead

[TABLE]

which relates the variable $S_{1}$ and the logarithmic derivative of the gravitational potential [41]. The third equation implies that the function $r$, defined by

[TABLE]

is monotone for $y\neq 0,r\neq 0$, and $S_{1}\neq 1$. Additionally, the function $Z=\frac{Y}{1-Y},\quad\frac{dZ}{d\lambda}=-\frac{S_{1}}{n+1}Z$ is monotonic for $S_{1}\neq 0,Y\neq 0$. Together with other monotone functions defined on the boundary subsets, this implies that all of the attractors for $\texttt{b}\geq 0$ (when the phase space is compact) are equilibrium points on the boundary sets $Y=0,Y=1$, $S_{2}=0$ and $P:=1-\texttt{b}S_{1}^{2}-S_{2}=0$ which are invariant for the flow of (57).

In the table 1 the equilibrium points and the corresponding eigenvalues of the dynamical system (57) are presented (see A.1 for details; the behavior at infinity is discussed in A.2).

We have the $m-r$-relation [48]

[TABLE]

where $\mathcal{M}(r)$ denotes the Misner-Sharp mass [46]. Using the equations (3.1.2), we can define the compact variables

[TABLE]

which depends only of the phase-space variables $(S_{1},S_{2},Y)$. Evaluating numerically the expressions $R_{\text{comp}},M_{\text{comp}}$ at the orbits of the system (57), we can see whether the resulting model leads to finite radius and finite mass (see discussion in Section 1.4, that follows from reference [48]).

In table 1 we present the equilibrium points and the corresponding eigenvalues of the dynamical system (57).

In the Figure (1(a)) are shown some orbits of the system (57) on the plane $S_{2}=0$ for $\texttt{b}=1,\texttt{h}=1$. In the Figure (1(b)) are shown some orbits of the system (57) on the plane $S_{2}=0$ for $\texttt{b}=1.25,\texttt{h}=1$. In the Figures 1(c) and 1(d) is presented how the compact variables $R_{\text{comp}},M_{\text{comp}}$, defined by (96) varies along the orbits of (57) for different choices of the parameters.

3.1.3 Dynamical systems analysis based on the Newtonian homology invariants

In this section we use dynamical system tools to determine conditions under which stable stars can form. We start with the set of equations (39), which determine the star’s structure and the geometry in the static spherically symmetric Einstein-æther theory for a perfect fluid. Equation (39a) is a modification of the well-known Tolman-Oppenheimer-Volkoff equation for relativistic star models. The idea is to construct a 3D dynamical system in compact variables and obtain a global picture of the solution space for a linear equation of state, that can be visualized in a geometrical way. This study can be extended to a wide class of equations of state. We therefore introduce the $(U,V,Y)$- formulation as follows: we define

[TABLE]

that are obtained by a compactification of the so- called Newtonian homology invariants [65]

[TABLE]

where $\mu$ and $p$ are the density and pressure of the background fluid. The homology invariants are defined in [65] by

[TABLE]

The relation between the variables $\{U,V,Y\}$ and $\{Q,S,C\}$ used in paper [35] are:

[TABLE]

Furthermore, the relation between the variables ${U,V}$ and $S_{1},S_{2}$ are

[TABLE]

To use the dynamical system tools, we have to define a new independent variable by

[TABLE]

Under the conditions $\mu\geq 0,p\geq 0$, it follows that the variables $U,V,Y$ are compact. The new variable ${\lambda}$ has the same range as $\ln(r)$, i.e., ${\lambda}\rightarrow-\infty$ as $r\rightarrow 0$ and ${\lambda}\rightarrow+\infty$ as $r\rightarrow+\infty$.

Assumptions of the positive energy density, pressure and mass, with the condition $\texttt{h}\geq 1$, leads to the inequalities:

[TABLE]

Since we are dealing with relativistic models the condition $2m/\rho<1$ has to be satisfied, so that we have

[TABLE]

As a first step before to studying the full system we consider the case $\texttt{b}=1$, which corresponds to the equations of GR (5). We compare our results with the results found in [41, 48].

General Relativity case ($\texttt{b}=1$)

Using the equations (5), and considering the time derivative (72), we obtain the dynamical system

[TABLE]

We have included the boundary subsets $U=0,U=1$, together with the “static surface” $(1+Y)V-(1-Y)=0$. From (75c) follows that the variable $Y$ is monotonically decreasing and all the orbits ends at the Newtonian boundary subset $Y=0$.

The equilibrium points and curves of equilibrium points of (75) are summarized in table 2 (see details in the A.3).

We have the useful $m-r$-relation

[TABLE]

where $\mathcal{M}(r)$ defines the Misner-Sharp mass [46]. Using the equations (3.1.3), we can define the compact variables

[TABLE]

which depends only of the phase-space variables $(U,V,Y)$, and evaluated along the orbits of the system (75), are useful to analyze whether the resulting model leads to finite radius and finite mass.

In the table 2 is summarized the information of the equilibrium points and the corresponding eigenvalues of the dynamical system defined by equations (75). In the Figure 2 it is presented how the compact variables $R_{\text{comp}},M_{\text{comp}}$, defined by (77), vary along the orbits of (75), and the attractor solutions has finite mass and finite radius. Additionally, are presented some orbits of the system (75).

Low pressure regime (The Newtonian subset $Y=0$)

We now consider the low pressure regime, also referred as the “Newtonian subset” $Y=0$, which is of great importance for relativistic stars. The stability analysis below is restricted to the dynamics onto the invariant set $Y=0$, and the stability along the $Y$-axis is not taken into account.

The evolution equations are given by

[TABLE]

defined on the phase plane

[TABLE]

The equilibrium (curves of equilibrium) points in this invariant set are:

$L_{1}(0):=(U,V)=\left(1,0\right)$, which is an endpoint of the line $L_{1}$ with $Y_{0}=0$ described in the previous section. The eigenvalues are $1,1$, such that it is a source. 2. 2.

$L_{2}(0):=(U,V)=\left(\frac{3}{4},0\right)$, which is an endpoint of the line $L_{2}$ with $Y_{0}=0$ described in the previous section. The eigenvalues are $-\frac{3}{4},\frac{1}{2}$, such that it is a saddle. 3. 3.

$L_{4}(0):=(U,V)=\left(0,0\right)$, which is an endpoint of the line $L_{4}$ with $Y_{0}=0$ described in the previous section. The eigenvalues are $3,-1$, such that it is a saddle. 4. 4.

$L_{6}:=(U,V)=(U_{0},1),U_{0}\in[0,1]$. The eigenvalues are $U_{0}-1,0$. It it a sink. 5. 5.

$M_{1}:=(U,V)=\left(0,\frac{1}{2}\right)$. The eigenvalues are $\frac{3}{2},\frac{1}{2}$, such that it is a source. 6. 6.

$M_{2}:=(U,V)=\left(0,1\right)$. The eigenvalues are $-1,0$, such that it is a non-hyperbolic equilibrium point (it is contained in the attractor line $L_{6}$ previously described). 7. 7.

$M_{5}:=(U,V)=\left(1,1\right)$. The eigenvalues are $0,0$. Hence, it is a non-hyperbolic equilibrium point (it is contained in the attractor line $L_{6}$ previously described).

In the Figure 3 are presented some streamlines of the system (78). The solid (red) line denotes the line $L_{6}$, which is a local attractor. The straight vertical line connecting $L_{2}(0)$ with the line $L_{6}$ is the so-called Tolman orbit (associated to the positive eigenvalue of $L_{2}(0)$, the so-called Tolman point).

High pressure regime ($Y=\texttt{h}^{-1}$)

In the limit $\texttt{b}=1,Y=\texttt{h}^{-1}$, the equations becomes

[TABLE]

The equilibrium points in this regime are

$L_{1}(1/\texttt{h})$: $(U,V)=(1,0)$, with eigenvalues $\frac{(\texttt{h}-1)^{3}}{\texttt{h}^{3}},\frac{(\texttt{h}-1)^{3}}{\texttt{h}^{3}}$. Local source. 2. 2.

$L_{2}(1/\texttt{h})$: $(U,V)=\left(\frac{3}{4},0\right)$, with eigenvalues $\frac{(\texttt{h}-1)^{3}}{2\texttt{h}^{3}},-\frac{3(\texttt{h}-1)^{3}}{4\texttt{h}^{3}}$. Saddle. 3. 3.

$L_{4}(1/\texttt{h})$: $(U,V)=(0,0)$, with eigenvalues $\frac{3(\texttt{h}-1)^{3}}{\texttt{h}^{3}},-\frac{(\texttt{h}-1)^{3}}{\texttt{h}^{3}}$. Saddle. 4. 4.

$M_{6}$: $(U,V)=\left(0,\frac{\texttt{h}-1}{\texttt{h}+1}\right)$, with eigenvalues $\frac{2(\texttt{h}-1)^{3}}{\texttt{h}^{3}(\texttt{h}+1)},-\frac{(\texttt{h}-1)^{3}}{\texttt{h}^{2}(\texttt{h}+1)}$. Exists for $\texttt{h}>1$. Saddle. 5. 5.

$M_{7}$: $(U,V)=\left(1,\frac{\texttt{h}-1}{\texttt{h}+1}\right)$. The eigenvalues are $-\frac{2(\texttt{h}-1)^{3}}{\texttt{h}^{3}(\texttt{h}+1)},\frac{(\texttt{h}-1)^{2}}{\texttt{h}^{2}(\texttt{h}+1)}$. Exists for $\texttt{h}>1$. Saddle. 6. 6.

$M_{8}$: $(U,V)=\left(\frac{1}{2},\frac{2(\texttt{h}-1)^{2}}{3\texttt{h}^{2}-2}\right)$, with eigenvalues

$-\frac{(\texttt{h}-1)^{3}\left(\texttt{h}+\sqrt{36-7\texttt{h}(\texttt{h}+4)}+2\right)}{4\texttt{h}^{2}\left(3\texttt{h}^{2}-2\right)},-\frac{(\texttt{h}-1)^{3}\left(\texttt{h}-\sqrt{36-7\texttt{h}(\texttt{h}+4)}+2\right)}{4\texttt{h}^{2}\left(3\texttt{h}^{2}-2\right)}$. Exists for $\texttt{h}>1$. It is an stable node for $1<\texttt{h}\leq\frac{2}{7}\left(4\sqrt{7}-7\right)$. Stable spiral for $\texttt{h}>\frac{2}{7}\left(4\sqrt{7}-7\right)$. 7. 7.

$N_{1}$: $(U,V)=(0,1)$, with eigenvalues $-\frac{(\texttt{h}-1)^{2}(\texttt{h}+6)}{\texttt{h}^{3}},-\frac{2(\texttt{h}-1)^{2}}{\texttt{h}^{3}}$. This point is over the static surface. It can be a local attractor on the invariant set $Y=1/\texttt{h}$. 8. 8.

$N_{2}$: $(U,V)=(1,1)$, with eigenvalues $\frac{2(\texttt{h}-1)^{2}}{\texttt{h}^{3}},-\frac{(\texttt{h}-2)(\texttt{h}-1)}{\texttt{h}^{3}}$. This point is over the static surface. It is a local source for $1<\texttt{h}<2$. It is a saddle for $\texttt{h}>2$. 9. 9.

$N_{3}$: $(U,V)=\left(\frac{(\texttt{h}-1)(\texttt{h}+6)}{\texttt{h}^{2}+6\texttt{h}-8},1\right)$, with eigenvalues $\frac{(\texttt{h}-2)(\texttt{h}-1)^{2}(\texttt{h}+6)}{\texttt{h}^{3}(\texttt{h}(\texttt{h}+6)-8)},\frac{2(\texttt{h}-1)^{2}(\texttt{h}(\texttt{h}+4)-4)}{\texttt{h}^{3}(\texttt{h}(\texttt{h}+6)-8)}$. Exists for $\texttt{h}\geq 2$. This point is over the static surface. It is a local source for $\texttt{h}>2$.

In the Figure 4 are presented some streamlines of the system (80) which corresponds to the high pressure regime. In this regime, the solution $M_{8}$, that has finite mass and radius, is a local attractor. The line $V=\frac{\texttt{h}-1}{\texttt{h}+1}$ acts as a separatrix. The orbits above this line are attracted by $N_{1}$.

Einstein-æther modification

In this case the dynamical system is given by

[TABLE]

where

[TABLE]

Observe that $\lim_{\texttt{b}\rightarrow 1}G=-1$. In the limit $\texttt{b}\rightarrow 1$, we obtain an equivalent dynamical system to (75).

The dynamical system (81) is defined on the bounded phase space

[TABLE]

From the equation (81c) we determine the invariant sets $U=1,Y=1,Y=\frac{1}{\texttt{h}}$. In the table 3 are presented the equilibrium points and the corresponding eigenvalues and stability conditions of the dynamical system (81). See the details of the analysis in A.4. In the analysis of the solution space we identify the low pressure regime (so-called, the Newtonian subset), which corresponds to the invariant set $Y=0$.

Low pressure regime (The Newtonian subset)

We now consider the “Newtonian subset” $Y=0$, where the evolution equations are given by

[TABLE]

defined on the phase plane

[TABLE]

The equilibrium (curves of equilibrium) points in the invariant set $Y=0$ are:

$M_{1}:=(U,V)=\left(0,\frac{1}{2}\right)$. The eigenvalues are $\frac{3}{2},\frac{1}{2}$, such that it is a source. 2. 2.

$M_{2}:=(U,V)=\left(0,1\right)$. The eigenvalues are $-1,0$, such that it is a non-hyperbolic equilibrium point (it is contained in the attractor line $L_{6}$ described below). 3. 3.

$M_{5}:=(U,V)=\left(1,1\right)$. The eigenvalues are $0,0$. Hence, it is a non-hyperbolic equilibrium point (it is contained in the attractor line $L_{6}$ described below). 4. 4.

$L_{1}(0):=(U,V)=\left(1,0\right)$. The eigenvalues are $\frac{1}{\texttt{b}},\frac{1}{\texttt{b}}$, such that it is a source for $\texttt{b}>0$ and a sink for $\texttt{b}<0$. 5. 5.

$L_{2}(0):=(U,V)=\left(\frac{3\texttt{b}}{1+3\texttt{b}},0\right)$, which is an endpoint of the line $L_{2}$ described in the previous section. It exists for $\texttt{b}\geq 0$. The eigenvalues are $-\frac{3}{3\texttt{b}+1},\frac{2}{3\texttt{b}+1}$, such that it is a saddle. 6. 6.

$L_{4}(0):=(U,V)=\left(0,0\right)$, which is an endpoint of the line $L_{4}$ described in the previous section. The eigenvalues are $3,-1$, such that it is a saddle. 7. 7.

The line $L_{6}:=(U,V)=(U_{0},1),U_{0}\in[0,1]$. The eigenvalues are $U_{0}-1,0$. This line is normally hyperbolic and due to the non-zero eigenvalue is negative it is a sink.

We see that the main difference with respect the case previously discussed is in that $L_{1}(0)$ can be an attractor for $\texttt{b}<0$, see figure 5(a). In the Figure 5(b) it is presented the case for $\texttt{b}=0.9>0$, essentially the same behavior as for the model (78).

3.1.4 Physical discussion

By using the $(S_{1},S_{2},Y)$ formulation, we have deduced the system (57). In this model the line of equilibrium points $L_{1}:(S_{1},S_{2},Y)=\left({S_{1}^{*}},1-\texttt{b}{S_{1}^{*}}^{2},0\right)$ is normally hyperbolic. The sector of $L_{1}$ with $S_{1}^{*}>0$ is stable, whereas, the sector of $L_{1}$ with $S_{1}^{*}<0$ is unstable. Another important equilibrium point is $P_{4}:\left(\frac{\texttt{h}+2}{4\texttt{b}},0,\frac{1}{\texttt{h}}\right)$; a sink for $\texttt{b}<0,\texttt{h}>1$.

For the low pressure regime governed by the equations (78), the line $L_{6}$ is a local attractor. The Tolman orbit is the one that connects the Tolman point $L_{2}(0)$ with the line $L_{6}$, and it is associated with the positive eigenvalue of $L_{2}(0)$. When $\texttt{b}<0$ (Einstein-æther modification), $L_{1}(0)$ can be an attractor, whereas for $\texttt{b}>0$, we essentially have the same behavior as for the model (78).

In the high pressure regime and for $\texttt{h}>1$ and $\texttt{b}=1$, it is obtained the system of differential equations (80), for which the attractor is $M_{8}$: $(U,V)=\left(\frac{1}{2},\frac{2(\texttt{h}-1)^{2}}{3\texttt{h}^{2}-2}\right)$. It is stable node for $1<\texttt{h}\leq\frac{2}{7}\left(4\sqrt{7}-7\right)$; Stable spiral for $\texttt{h}>\frac{2}{7}\left(4\sqrt{7}-7\right)$. This point represents physical solutions that have finite mass and radius.

Implementing the procedure described in A.2 (which applies for $\texttt{b}<0$), we conclude that there are not physically relevant equilibrium points at infinity, apart from the source point $I_{4}$ that has $S_{1}\rightarrow\infty$, $S_{2}$ finite and $Y=\frac{1}{\texttt{h}}$. This solution corresponds to $y$ finite and $x\rightarrow y$. That is, $\frac{dp}{d\rho}=\frac{1}{\mu+p},\frac{p}{p+\mu}=\frac{1}{\texttt{h}},$ where $\rho$ is the radial coordinate. Then $\frac{dp}{d\rho}=\frac{1}{\texttt{h}p}$ implies that $p(\rho)=\frac{\sqrt{\texttt{h}p_{c}^{2}+2\rho}}{\sqrt{\texttt{h}}},$ where $p_{c}$ is the central pressure (at $\rho=0$). As $\rho\rightarrow\infty$, the pressure increases as the solutions move away from this point.

3.2 Model 2: Polytropic equation of state

We use the line element

[TABLE]

where $\psi(\lambda)$ is the gravitational potential, $r(\lambda)$ is the usual Schwarzschild radial parameter, and $N(\lambda)$ is a dimensionless (under scalar transformations) freely specifiable function, that defines the spatial radial coordinate $\lambda$. The choice $N=1$ corresponds to isotropic coordinates and the function $N$ can viewed as a relative gauge function with respect to the isotropic gauge [41].

The FE are [18]:

[TABLE]

where for any function $f(\lambda)$, we define ${\mathbf{e}_{1}}(f)\equiv\frac{f^{\prime}(\lambda)}{N(\lambda)r(\lambda)}$, and

$x=-\frac{{\mathbf{e}_{1}}(r)}{r},\quad y={\mathbf{e}_{1}}(\psi),\quad K=\frac{1}{r^{2}}$. Furthermore there exists the constraint equation

[TABLE]

Choosing a new radial coordinate $\tau$ such that

[TABLE]

the system (87) becomes

[TABLE]

with restriction (88). The above equations transforms to eqs. (A.3) in [41] for the choice $\texttt{b}=1$, under the identifications $\theta=y-x,\sigma=y,B^{2}=K$. We consider a polytropic equation of state $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$.

3.2.1 Singularity analysis and integrability

The FE does not pass the Painlevè test, and therefore they are not integrable.

3.2.2 Equilibrium points in the finite region of the phase space

For the dynamical system’s formulation we introduce the variables

[TABLE]

Assuming that the energy density $\mu$ and the pressure $p$ are both non-negative we obtain $0\leq Y\leq 1$, where we have attached the boundaries $Y=0$ and $Y=1$. On the other hand, $S_{2}\geq 0$ by definition.

Now we define the independent dimensionless variable $\lambda$ by the choice $N^{2}r^{2}=Y^{2}/\theta^{2}$, that is

[TABLE]

and we obtain the dynamical system

[TABLE]

where we used a polytropic equation of state $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$, where $q$ is a nonnegative constant. For $\texttt{b}=1$ it is recovered the model of [41].

Using the same arguments as in Section 3.1.2, we have that $Y=0,Y=1,S_{2}=0$, and $1-\texttt{b}S_{1}^{2}-S_{2}=0$, are invariant sets for the flow of (93), and we can argue that system (93) defines a flow on the (invariant subset) phase space

[TABLE]

We have the useful the $m-r$ relation

[TABLE]

Using the equations (3.2.2), we can define the compact variables, which depends only of the phase-space variables $(S_{1},S_{2},Y)$,

[TABLE]

Evaluating numerically the expressions $R_{\text{comp}},M_{\text{comp}}$ at the orbits of the system (93), we can see whether the resulting model leads to finite radius and finite mass. Furthermore, imposing regularity at the center, that is, as $\lambda\rightarrow-\infty$, and are extrapolated the conditions for relativistic stars as given by the Buchdahl inequalities [66, 67], which in units where $8\pi G=1$ are expressed as

[TABLE]

where $Y_{c}$ is related with the central pressure of the fluid by

$Y_{c}=\frac{1}{p_{c}^{-\frac{1}{n+1}}q^{-\frac{n}{n+1}}+1},\quad\left[\frac{Y_{c}}{1-Y_{c}}\right]^{1+n}=p_{c}q^{n}$.

In the Table 4 it is summarized the information of the equilibrium points and the corresponding eigenvalues and stability conditions of the dynamical system (93) (see details in B.1; for the analysis at infinity see B.2).

In the Figure 6 are presented some streamlines of the system (93) in the invariant set $Y=0$. The red (thick) line denotes the singular curve $L_{1}$. The left branch $S_{1}<0$ is unstable and the right branch $S_{1}>0$ is stable. In the Figure (7) are presented some streamlines of the system (93) in the invariant set $Y=1$. The red (thick) line denotes the singular curve $L_{1}$. In figure (7) (a), and (b), the endpoint of line $L_{2}$ is a saddle. $P_{1}$ is the local attractor, $P_{2}$ and $P_{3}$ are local sources and $P_{4}$ is a saddle. In figure (7) (c) with $\texttt{b}=\frac{9}{8}$, $P_{1}$ and $P_{4}$ coincides and they are local attractors. In (d), $P_{4}$ is a local attractor (it is an hyperbolic saddle in the full phase space and there is a 1-parameter set of orbits that enters the interior of the phase space from it in the direction of the eigenvector associated to the positive eigenvalue for $\texttt{b}>0$). Recall for $\texttt{b}<0$ the attractor is $P_{4}$. In Figure 8 are displayed some streamlines of the system (93) in the invariant set $S_{2}=0$. This set corresponds to the plane-symmetric boundary. Finally, in the Figure 9 are presented some streamlines of the system (93) in the invariant set $S_{2}=1-\texttt{b}S_{1}^{2}$ (surface of pressure zero) projected on the plane $(Y,S_{1})$. The line $L_{1}$ is denoted by a vertical dashed red line. The line $L_{2}$ is denoted by an horizontal red solid line. In the Figure 10(c) it is presented how the compact variables $R_{\text{comp}},M_{\text{comp}}$, defined by (96) varies along the orbits of (93). In the Figure 10(d) it is presented how the compact variables $R_{\text{comp}},M_{\text{comp}}$, defined by (96) varies along the orbits of (93). The last two plots shows that the typical behavior for this choice of parameters is to obtain regions of infinite radius and finite mass, that is, $(R_{\text{comp}},M_{\text{comp}})=(1,0)$ is a local attractor in the $(R_{\text{comp}},M_{\text{comp}})$ diagram.

3.2.3 Dynamical systems analysis based on the Newtonian homology invariants

As in section 3.1.3 we can implement the $(U,V,Y)$- formulation introducing the variables

[TABLE]

where $\mathcal{M}$ is the Misner-Sharp mass defined by $\mathcal{M}(r)=\frac{m}{r}$, $r(\lambda)$ is the usual Schwarzschild radial parameter introduced in (86), $\mu$ and $p$ are the density and pressure of the scalar field with a polytropic equation of state $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0,q\geq 0$. Additionally, we introduce the radial derivative

[TABLE]

The relation between the variables ${U,V}$ and $S_{1},S_{2}$ are

[TABLE]

We will see that there are several equilibrium points in this formulation that cannot be studied with the $\{S_{1},S_{2},Y\}$-formulation.

Equations for General Relativity ($\texttt{b}=1$)

In this example we obtain the dynamical system

[TABLE]

defined in the phase space

[TABLE]

The static surface

[TABLE]

is an invariant set, due to the function $\mathcal{S}$ defined by

[TABLE]

satisfies

[TABLE]

Additionally, any of the sets $U=0,U=1,V=0,V=1,Y=0,Y=1$, defines an invariant set.

The equilibrium points of the system (101) are summarized in A.3 (they were found in the reference [41], we use here the same notation). We distinguish two asymptotic regimes: the low pressure regime corresponding to the Newtonian subset, and the high pressure regime.

Low pressure regime (the Newtonian subset)

We now consider the “Newtonian subset” $Y=0$. In this subset the dynamical system equations are reduced to

[TABLE]

The equilibrium points in the Newtonian subset are

$L_{1}(0):(U,V)=(1,0)$, with eigenvalues $1,1$. Source. 2. 2.

$L_{2}(0):(U,V)=\left(\frac{3}{4},0\right)$, with eigenvalues $-\frac{3}{4},\frac{1}{2}$. Saddle. 3. 3.

$L_{4}(0):(U,V)=(0,0)$, with eigenvalues $3,-1$. Saddle. 4. 4.

$L_{6}:(U,V)=(U_{0},1),U_{0}\in[0,1]$, with eigenvalues $0,-(1-U_{0})$. Exists for $n=0$. This line is normally hyperbolic and due to the non-zero eigenvalue is negative it is a sink. 5. 5.

$P_{1}:(U,V)=\left(0,\frac{n+1}{n+2}\right)$, with eigenvalues $\frac{1}{n+2},-\frac{n-3}{n+2}$. Source for $0\leq n<3$, nonhyperbolic for $n=3$, saddle when $n>3$. 6. 6.

$P_{2}:(U,V)=(0,1)$, with eigenvalues $-\frac{1}{n+1},-\frac{n}{n+1}$. Nonhyperbolic for $n=0$, sink for $n>0$. 7. 7.

$P_{4}:(U,V)=\left(\frac{n-3}{2(n-2)},\frac{2(n+1)}{3n+1}\right)$, with eigenvalues

$-\frac{(n-1)\left(n-5+\sqrt{1-7n^{2}+22n}\right)}{4(n-2)(3n+1)},-\frac{(n-1)\left(n-5-\sqrt{1-7n^{2}+22n-1}\right)}{4(n-2)(3n+1)}$. Exists for $n>3$. Unstable node for $3<n\leq\frac{1}{7}\left(11+8\sqrt{2}\right)\approx 3.18767$. Unstable spiral for $\frac{1}{7}\left(11+8\sqrt{2}\right)<n<5$. It is a center for $n=5$. Stable spiral for $n>5$. 8. 8.

$P_{5}:(U,V)=(1,1)$ with eigenvalues $0,0$. Nonhyperbolic.

We deduce the useful $m-r$-relation

[TABLE]

Using the equations (3.2.3) , we can define the compact variables

[TABLE]

which depends only of the phase-space variables $(U,V,Y)$. Evaluating numerically the expressions $R_{\text{comp}},M_{\text{comp}}$ at the orbits of the system (101), we can see whether the resulting model leads to finite radius and finite mass.

In the Figure 11 are presented some orbits of the system (101) for the choices $n=4,n=5,n=6$ for the relativistic case ($\texttt{b}=1$).

In the Figure 12 it is depicted the typical behavior of the streamlines of the system (106), for the choices (a) $n=0$, (b) $0<n\leq 3$, (c) $3<n\leq\frac{1}{7}\left(11+8\sqrt{2}\right)\approx 3.18767$, (d) $\frac{1}{7}\left(11+8\sqrt{2}\right)<n<5$, (e) $n=5$, and (f) $n>5$ in accordance with results summarized in the list above.

Einstein-æther modification

The dynamical system becomes

[TABLE]

where

[TABLE]

defined in the bounded phase space

[TABLE]

The details of the stability of the relevant equilibrium point of the dynamical system (109) are presented in Table 5. Evaluating numerically the expressions $R_{\text{comp}},M_{\text{comp}}$ defined by (108) at the orbits of the system (109), we can see whether the resulting model leads to finite radius and finite mass.

The system (109) admits the relevant invariant sets $Y=0$ (the low-pressure Newtonian subset) and the invariant set $Y=1$ (high-pressure subset). Now we discuss about the dynamics on these invariant sets (see full details in B.4).

Low pressure regime (The Newtonian subset)

In the invariant set $Y=0$ the dynamics is governed by the equations

[TABLE]

In the especial case $\texttt{b}=1$ are recovered the equations (106), and therefore we have the know results for stellar and relativistic dynamics investigated in the reference [41]. The equilibrium points that lies on the invariant set $Y=0$ are:

$L_{1}(0)$: $(U,V)=(1,0)$. The eigenvalues are $\frac{1}{\texttt{b}},\frac{1}{\texttt{b}}$. It is a source for $\texttt{b}\geq 0$ and a sink for $\texttt{b}<0$. Therefore, the dynamics is richer as compared with the previous case where it is always a source. 2. 2.

$L_{2}(0)$: $(U,V)=\left(\frac{3\texttt{b}}{3\texttt{b}+1},0\right)$. The eigenvalues are $-\frac{3}{3\texttt{b}+1},\frac{2}{3\texttt{b}+1}$. It is a saddle. 3. 3.

$L_{4}(0)$: $(U,V)=(0,0)$. The eigenvalues are $-1,3$. It is a saddle. 4. 4.

$L_{6}$: $(U,V)=(U_{0},1),U_{0}\in[0,1]$, with eigenvalues $0,-(1-U_{0})$. Exists for $n=0$. This line is normally hyperbolic and due to the non-zero eigenvalue is negative it is a sink. 5. 5.

$P_{1}$: $(U,V)=\left(0,\frac{n+1}{n+2}\right)$. The eigenvalues are $\frac{1}{n+2},-\frac{n-3}{n+2}$. Source for $0\leq n<3$, nonhyperbolic for $n=3$, saddle when $n>3$. This equilibrium point has the same properties as for the Newtonian subset. 6. 6.

$P_{2}$: $(U,V))=(0,1)$. The eigenvalues are $-\frac{n}{n+1},-\frac{1}{n+1}$. Nonhyperbolic for $n=0$, sink for $n>0$. 7. 7.

$P_{4}$: $(U,V)=\left(\frac{\texttt{b}(n-3)}{\texttt{b}(n-3)+n-1},\frac{2(n+1)}{3n+1}\right)$. Exists for $\texttt{b}<0,1\leq n\leq 3$, or $\texttt{b}=0,n>1$, or $\texttt{b}>0,n=1$, or $\texttt{b}>0,n\geq 3$, with eigenvalues

$-\frac{(n-1)\left(n+\sqrt{1-n(7n-22)}-5\right)}{2(3n+1)(\texttt{b}(n-3)+n-1)},\frac{(n-1)\left(-n+\sqrt{1-n(7n-22)}+5\right)}{2(3n+1)(\texttt{b}(n-3)+n-1)}$. As compared with the Newtonian case, the dynamical features of this point are more intricate (starting that it not only exists when $n>3$, but also that it is a two-parametric solution). It is an unstable node for $\texttt{b}\geq 0,3<n\leq\frac{1}{7}\left(11+8\sqrt{2}\right)\approx 3.18767$. Unstable spiral for $\texttt{b}\geq 0,\frac{1}{7}\left(11+8\sqrt{2}\right)<n<5$. Stable spiral for $\texttt{b}\geq 0,n>5$. Center for $\texttt{b}\neq 0,n=1$, or $n=3$. Saddle for $\texttt{b}\leq 0,1<n<3$. 8. 8.

$P_{5}$: $(U,V)=(1,1)$, with eigenvalues $0,0$. Nonhyperbolic.

In the Figure 13 are presented some streamlines of the system (112) for the choices $n=0,n=3,n=3.18$. In the Figure 14 are presented some streamlines of the system (112) for the choices $n=4,n=5,n=6$.

High pressure regime

In the invariant set $Y=1$ the dynamics is given by the equations

[TABLE]

The equilibrium points in the invariant set $Y=1$ are

$L_{3}$: $(U,V)=(U_{0},0),\quad 0\leq U_{0}\leq 1$, with eigenvalues $0,0$. Nonhyperbolic. 2. 2.

$L_{5}$: $(U,V)=(0,V_{0}),\quad 0\leq V_{0}\leq 1$, with eigenvalues $0,0$. Nonhyperbolic. 3. 3.

$P_{3}$: $(U,V)=(1,1)$, with eigenvalues $-\frac{2(4\texttt{b}-3)}{\texttt{b}},\frac{2(4\texttt{b}-3)}{\texttt{b}}$. Nonhyperbolic for $\texttt{b}=\frac{3}{4}$. It is a saddle otherwise.

In this example the system is integrable. The orbit passing through $(U_{0},V_{0})$ is

[TABLE]

3.2.4 Physical discussion

Case of Einstein- æther modification

Using the $(S_{1},S_{2},Y)$ formulation given by the system (93) we have derived the following results. For $\texttt{b}\neq 0$, one branch of $L_{1}$ is a local source and the other branch is a local sink. The line $P_{2}$ is a local source. The self-similar plane symmetric solution $P_{4}$ is unstable (saddle). The orbits coming from this point are associated with a negative mass singularity. The orbits in the upper half $S_{1}>0$ are attracted by either $L_{1}$ or $L_{2}$. $P_{2}$ is a local source. In the lower half $S_{1}<0$, the line $L_{1}$ is the local source and $L_{2}$ acts as the local sink. Each point located on $L_{2}$ corresponds to the flat Minkowski solution written in spherically symmetric form. An orbit, associated with the eigenvalue $2y_{c}$, with $y_{c}>0$, enters the interior of the phase space from each point of $L_{2}$ parametrized by $y_{c}=p_{c}/(p_{c}+\mu_{c})=q\mu_{c}^{1/n}/(1+q\mu_{c}^{1/n})$, where $p_{c}$ and $\mu_{c}$ are the central pressure and central density, respectively.

In GR, the point $(S_{1},S_{2},Y)=\left(\frac{2}{3},\frac{1}{9},1\right)$, corresponds to the so called self-similar Tolman solution discussed in [36]. In the Einstein-æther theory this solution is generalized to a 1-parameter set of solutions represented by the equilibrium points $P_{1}:(S_{1},S_{2},Y)=\left(\frac{2}{3},\frac{1}{9}(9-8\texttt{b}),1\right)$. The single orbit that enters the interior of the phase space from $P_{1}$, associated with the positive eigenvalue of $P_{1}$ is called the Tolman orbit. The solutions having $Y=1$ corresponds to the limiting situation $\frac{p}{\mu}\rightarrow\infty.$ In the case of a linear equation of state $p=(\gamma-1)\mu$, this is equivalent to $\gamma\rightarrow\infty$.

Concerning the singular curve $L_{1}$, the left branch $S_{1}<0$ is unstable and the right branch $S_{1}>0$ is stable. When an orbit ends on $L_{1}$, we can determine the total mass and total radius of the corresponding solution. Evaluating at $L_{1}$, where $S_{2}=1-\texttt{b}\Sigma_{c}^{2}$, we have $M/R=\Sigma_{c}(1-(1+\texttt{b})\Sigma_{c}/2)/(1-\texttt{b}\Sigma_{c}^{2})$. Thus, the solutions ending on $L_{1}$ when $\Sigma_{c}>0$ have finite radii and masses (say for $0<\Sigma_{c}(1-(1+\texttt{b})\Sigma_{c}/2)/(1-\texttt{b}\Sigma_{c}^{2})<\text{const.}$). The orbits ending on $L_{1}$ when $\Sigma_{c}=0$ or $\Sigma_{c}=2/(1+\texttt{b})$, describe solutions with infinite radii and finite masses, or masses that approach infinity slower than $r$.

Regarding the behavior of solutions at infinity, we refer to the discussion in B.2. In particular, for $\texttt{b}<0$ there are not physically relevant equilibrium points at infinity, apart from the point $I_{4}$, that has $S_{1}\rightarrow\infty$, $S_{2}$ finite and $Y=1$, which is a source (it exists only for $\texttt{b}<0$). This solution corresponds to $y$ finite and $x\rightarrow y$. That is, $\frac{dp}{d\rho}=\frac{1}{\mu+p},\frac{p}{p+\mu}=1,$ where $\rho$ is the radial coordinate. Then $\frac{dp}{d\rho}$ implies $\frac{1}{p}\implies p(\rho)=\sqrt{p_{c}^{2}+2\rho},$ where $p_{c}$ is the central pressure (at $\rho=0$). As $\rho\rightarrow\infty$ the pressure is increasing as the solutions move away from this point.

Case of General Relativity

The system (101) was fully studied in [41].

For $n\geq 0$, there exists regular Newtonian solutions, and when $0\leq n<5$ it is possible to obtain solutions that represents star models with finite radii and mass, while for $n\geq 5$ infinite regular models are obtained (for $n=5$ there are solutions with infinite radius but finite mass).

The choice $n=0$ corresponds to relativistic incompressible fluid models. The set $U=\frac{3}{4}$ contains all regular solutions. Since $L_{6}$ is the attractor line in this case, the solution (Tolman orbit) is the straight line that starts on $L_{2}$ and ends on the line $L_{6}$ at $U=\frac{3}{4},V=1$.

For the choice $0<n\leq 3$, the point $P_{2}$ is a hyperbolic sink. On the other hand, $(V(-1+Y)Y(1-Y+U(-1+2Y)))/(1+n)<0$ for any interior point of the phase space, which means that $Y$ is monotonically decreasing. Therefore, any orbit approaches the set $Y=0$. This result, combined with the existence of a monotonic function in the subset $Y=0$, implies that $P_{2}$ is a global sink, which means that the radii and masses of a regular models are finite for $0<n\leq 3$.

For $3<n<5$, the point $P_{2}$ is the only sink, but not a global one. There is a 1-parameter set of orbits that ends at $P_{1}$, acting as a separatrix surface in the state space. The corresponding solutions have finite masses but infinite radii. The boundary of the separatrix surface associated with $P_{1}$ is described by an orbit in the invariant set $U=0$; an orbit that connects $P_{1}$ with $P_{4}$ in the invariant set $Y=0$; and a relativistic orbit ending on $P_{4}$ (corresponding to a solution with infinite mass and infinite radius). If there is a regular separatrix surface intersecting the separatrix associated to $P_{1}$ or its boundary orbit at an interior point, the solutions will either end on $P_{1}$ or $P_{4}$. That is, the solutions do not necessarily end at $P_{2}$ (see the numerical elaborations performed in [41]).

When $n=5$, the separatrix surface that ends at the equilibrium point $P_{1}$ completely encloses the regular subset of orbits, which prevents any regular orbit ending at $P_{2}$. The Newtonian regular orbit ends on $P_{1}$, while all general relativistic regular orbits tend to $P_{4}$ (which is a center for $n=5$) and to the closed solutions surrounding $P_{4}$.

For $n>5$, the separatrix surface associated with $P_{1}$ encloses all regular solutions (including Newtonian orbit), and all regular solutions end on the hyperbolic sink $P_{4}$. Therefore, Newtonian and relativistic polytropic models have infinite radii and mass.

4 Stationary comoving æther with perfect fluid and scalar field in

static metric

We use now the metric

[TABLE]

where we defined $x=\mathbf{e}_{1}\ln e_{2}{}^{2},y=\mathbf{e}_{1}\ln N$, and the differential operator $\mathbf{e}_{1}=e_{1}{}^{1}\partial_{r}$. The equations for the variables $x,y,p,\phi,K$ are:

[TABLE]

[TABLE]

Taking the differential operator $\mathbf{e}_{1}(...)$ in both sides of (117), using the equations (116) to substitute the spatial derivatives, and using again the restriction (117) solved for $K,$ we obtain an identity. Thus, (116) (Gauss constraint) is a first integral of the system. The æther constraint is identically satisfied.

In order to apply Tolman-Oppenheimer-Volkoff (TOV) approach (we use units where $8\pi G=1$) we use the metric

[TABLE]

where $m(\rho)$ denotes the mass up to the radius $\rho$, with range $0\leq\rho<\infty$. Using the expressions (34b) and (35), and the equations (116a), (116d) and (117), we obtain the equations

[TABLE]

Equations (119) have two solutions for $p^{\prime}(\rho),m^{\prime}(\rho),\phi^{\prime\prime}(\rho)$. We consider the branch

[TABLE]

where $\Xi\equiv\sqrt{2\left(-2\texttt{b}m-(\texttt{b}-1)\rho^{3}(p-W(\phi))+\rho\right)-(\texttt{b}-1)\rho^{2}(\rho-2m)\phi^{\prime}(\rho)^{2}}$.

Taking the limit $\texttt{b}\rightarrow 1$ in Eqs. (120), we obtain the proper relativistic equations:

[TABLE]

where $\phi$ is the scalar field and $\psi$ is the gravitational potential. The above TOV formulation is general. For the analysis we have to specify the equation of state and the potential of the scalar field.

4.1 Model 3: Perfect fluid with linear equation of state and a scalar field with an exponential potential

Let us now show how (116) can be used to obtain exact solutions and in addition, be used in the dynamical systems approach for investigating the structure of the whole solution space. We use the linear equation of state

[TABLE]

where the constants $\mu_{0}$ and h satisfy $\mu_{0}\geq 0,\texttt{h}\geq 1.$ The case $\texttt{h}=1$ corresponds to an incompressible fluid with constant energy density, while the case $\mu_{0}=0$ describes a scale-invariant EoS.

4.1.1 Singularity analysis and integrability

Without loss of generality we assume the differential operator in the form $\mathbf{e}_{1}=\partial_{{\ell}}$, where we have introduced the radial rescaling $r=e^{\ell}$, such that $\ell\rightarrow-\infty$ as $r\rightarrow 0$ and $\ell\rightarrow\infty$ as $r\rightarrow\infty$. Furthermore with the use of (116c) and (116e) we rewrite the system (116a)-(116e) as that of three second-order ordinary differential equations in terms of the variables, $p\left(\ell\right)$, $K\left(\ell\right)$ and $\phi\left(\ell\right)$, where for the perfect fluid we assume the EoS (122) and assume that the potential for the scalar field is $W\left(\phi\right)=W_{0}e^{-k\phi}$. Recall that equations (116c) and (116e) provide

[TABLE]

Consider first that $\mu_{0}=0$, and we substitute $p\left({\ell}\right)=p_{0}{\ell}^{p_{1}},K\left({\ell}\right)=K_{0}{\ell}^{p_{2}},~{}\phi=\phi_{0}{\ell}^{p_{3}}$ in the system of second order equations. We search for the dominant behavior and we find that $p_{1}=p_{2}=p_{3}=-2$, while when the constants $p_{0},K_{0}$ and $\phi_{0}$ are

[TABLE]

and

[TABLE]

The dominant behavior is also a solution of the FE.

We follow the ARS algorithm which we described above and we find the six resonances to be

[TABLE]

[TABLE]

[TABLE]

We conclude that the FE pass the singularity test.

On the other hand, for $\mu_{0}\neq 0$, the dominant terms are not a specific solution of the FE. Furthermore, the ARS algorithm fails and the system does not pass the singularity test. Now that we have studied the integrability of the field FE, we continue with the phase-space analysis.

4.1.2 Equilibrium points in the finite region of the phase space

In this section we introduce a new dynamical system formulation best suited to the case when a scalar field is included in the matter content. We investigate here both the general relativistic case ($\texttt{b}=1$) and the Einstein-æther modification using the equations (116). The advantage of this formulation is that the resulting system is polynomial, and we do not have rational functions as in the TOV formulation. Therefore, we define the normalized variables

[TABLE]

and the time derivative

[TABLE]

to obtain the dynamical system

[TABLE]

Additionally,

[TABLE]

Therefore, the system defines a flow on the invariant set

[TABLE]

The system is invariant under the change of coordinates

$(X_{\phi},Y_{\phi},S_{1})\rightarrow(-X_{\phi},-Y_{\phi},-S_{1})$ with the simultaneous reversal in the independent variable ${\lambda}\rightarrow-{\lambda}$. In relation to the phase space dynamics this implies that for two points related by this symmetry, say $P^{+}$ and $P^{-}$, each has the opposite dynamical behavior to the other; that is, if the equilibrium point $P^{+}$ is an attractor for a given set of parameters, then $P^{-}$ is a source under the same choice.

The details of the stability of the equilibrium points (curves/surfaces of equilibrium points) of the dynamical system (131) are presented in the C. These results are summarized in Tables 6 and 7.

We have the useful $m-r$-relation

[TABLE]

Using the equations (4.1.2), we can define the compact variables, which depends only of the phase-space variables $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)$,

[TABLE]

For determine whether the resulting model leads to finite radius and finite mass one can evaluate $R_{\text{comp}},M_{\text{comp}}$ along the orbits of the system (131).

Now, we present some numerical integrations of the system (131). In Figure 15 we present some 2D projections of the orbits of the dynamical system (131) for $\texttt{b}=1,\texttt{h}=1,k=-3$. In figure 16 we present some 3D projections of the orbits of the system (113) for the same choice of parameters and the same selection of the initial conditions as in Figure 15. In this figures it is illustrated that the attractor as $\lambda\rightarrow\infty$, is $P_{7}^{-}$, that for this choice of parameters is is asymptotically the Minkowski solution (as discussed before). The past attractor as $\lambda\rightarrow-\infty$, is $P_{7}^{+}$.

4.1.3 Physical discussion

According to our analytical results, the equilibrium points with 5D unstable or stable manifolds are $P_{5}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(0,\varepsilon,\frac{\texttt{h}+2}{4\texttt{b}}\varepsilon,0,\frac{1}{\texttt{h}}\right)$, and $P_{7}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(\frac{k\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{2\sqrt{2}\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,0,\frac{1}{\texttt{h}}\right)$. $P_{5}^{+}$ (resp. $P_{5}^{-}$) is a sink (resp. a source) for $\texttt{b}<0,\texttt{h}>1$. They are saddle otherwise. $P_{7}^{+}$ (resp. $P_{7}^{-}$) is a source (resp. sink) for $k\leq-2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},1<\texttt{h}<\texttt{b}k^{2}$, or $-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or $2\sqrt{\frac{2}{3}}<k<2\sqrt{2},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or $k\geq 2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},1<\texttt{h}<\texttt{b}k^{2}$. Following the discussion in [35], $P_{5}^{\pm}$ satisfies $(Q,S,C,A_{\phi},A_{W})=\left(\pm 1,\pm\frac{\texttt{h}+2}{4\texttt{b}},0,0,0\right)$. The physical solutions have the asymptotic metric:

$ds^{2}=-\bar{N}_{0}^{2}e^{\pm\frac{(\texttt{h}+2)\bar{\lambda}}{2\texttt{b}}}dt^{2}+\frac{(\texttt{h}+2)^{2}e^{\pm\frac{\texttt{h}(\texttt{h}+2)\bar{\lambda}}{4\texttt{b}}}}{16c_{1}^{2}\texttt{b}^{2}}d\bar{\lambda}^{2}+{\bar{K}_{0}}^{-1}e^{\pm 2\left(1-\frac{\texttt{h}+2}{4\texttt{b}}\right)\bar{\lambda}}(d\vartheta^{2}+\sin^{2}\vartheta d\varphi^{2})$.

Introducing $\rho_{\pm}=e^{\pm\left(1-\frac{\texttt{h}+2}{4\texttt{b}}\right)\bar{\lambda}}$, the line element becomes

$-\bar{N}_{0}^{2}\rho_{\pm}^{\frac{2(\texttt{h}+2)}{4\texttt{b}-\texttt{h}-2}}dt^{2}+\frac{(\texttt{h}+2)^{2}}{2c_{1}^{2}(4\texttt{b}-\texttt{h}-2)^{2}}\rho_{\pm}^{-\frac{16\texttt{b}-\texttt{h}^{2}-6\texttt{h}+24}{2(4\texttt{b}-\texttt{h}-2)}}d\rho_{\pm}^{2}+\frac{\rho_{\pm}^{2}}{\bar{K}_{0}}(d\vartheta^{2}+\sin^{2}\vartheta d\varphi^{2})$.

For $\texttt{b}=1$, these solutions are the analogues of $P_{5}^{\pm}$ investigated in [36], which correspond to self-similar plane-symmetric perfect fluid models. For $\texttt{h}\geq 1,16\texttt{b}\geq(\texttt{h}+2)(\texttt{h}+4)$, the exponent of the component $tt$ is positive and the exponent of the $\rho_{\pm}\rho_{\pm}$ component is negative. Thus the singularity has a horizon at $\rho_{\pm}=0^{+}$ .

The points $P_{7}^{\varepsilon}$ have $(Q,S,C,A_{\phi},A_{W})=\left(1,\frac{\varepsilon}{2\texttt{b}},0,\frac{k\varepsilon}{2\sqrt{2}},\sqrt{\frac{-8\texttt{b}+\texttt{b}k^{2}+2}{8\texttt{b}}}\right)$. That is, $P_{7}^{\pm}$ is exactly the equilibrium point $P_{18}^{\pm}$ studied in [35]. According to the results in [35], the asymptotic metric is given by $N=\bar{N}_{0}e^{\pm\frac{\bar{\lambda}}{2\texttt{b}}},K=\bar{K}_{0}e^{\mp\frac{(2\texttt{b}-1)\bar{\lambda}}{\texttt{b}}},y=\bar{y}_{0}e^{\mp\frac{1}{4}k^{2}\bar{\lambda}}$. That is, the metric can be written as:

$ds^{2}=-{\bar{N}_{0}}^{2}e^{\pm\frac{\bar{\lambda}}{\texttt{b}}}dt^{2}+\frac{1}{4\bar{y}_{0}^{2}\texttt{b}^{2}}e^{\pm\frac{k^{2}\bar{\lambda}}{2}}d\bar{\lambda}^{2}+{\bar{K}_{0}}^{-1}e^{\pm\frac{(2\texttt{b}-1)\bar{\lambda}}{\texttt{b}}}\left(d\vartheta^{2}+\sin^{2}\vartheta d\varphi^{2}\right)$.

The sink is $P_{7}^{-}$ with asymptotic line element:

$ds^{2}=-{\bar{N}_{0}}^{2}e^{-\frac{\bar{\lambda}}{\texttt{b}}}dt^{2}+\frac{1}{4\bar{y}_{0}^{2}\texttt{b}^{2}}e^{-\frac{k^{2}\bar{\lambda}}{2}}d\bar{\lambda}^{2}+{\bar{K}_{0}}^{-1}e^{-\frac{(2\texttt{b}-1)\bar{\lambda}}{\texttt{b}}}\left(d\vartheta^{2}+\sin^{2}\vartheta d\varphi^{2}\right)$,

with some conditions on the parameters. Imposing the additional condition $\texttt{b}>\frac{1}{2}$ (i.e., when $k\leq-2\sqrt{2},\texttt{b}>\frac{1}{2},1<\texttt{h}<\texttt{b}k^{2}$, or $-2\sqrt{2}<k<-2,\frac{1}{2}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or $2<k<2\sqrt{2},\frac{1}{2}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or $k\geq 2\sqrt{2},\texttt{b}>\frac{1}{2},1<\texttt{h}<\texttt{b}k^{2}$) the solution becomes the Minkowski solution as $\lambda\rightarrow+\infty$. However, in all the cases $k\leq-2,\frac{1}{k^{2}}<\texttt{b}<\frac{1}{2},1<\texttt{h}<\texttt{b}k^{2}$, or $-2<k<-2\sqrt{\frac{2}{3}},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or $2\sqrt{\frac{2}{3}}<k<2,\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or $k\geq 2,\frac{1}{k^{2}}<\texttt{b}<\frac{1}{2},1<\texttt{h}<\texttt{b}k^{2}$ (all have $\texttt{b}<\frac{1}{2}$), the metric component $e^{-\frac{(2\texttt{b}-1)\bar{\lambda}}{\texttt{b}}}\rightarrow\infty$ as $\lambda\rightarrow\infty$. Therefore, the solution becomes singular (in two directions the length is stretched, and in one direction it becomes large; like a cigar singularity). This case does not appear in GR where $\texttt{b}=1$.

4.2 Model 4: Perfect fluid with polytropic of state and a scalar field with an exponential potential

Using the parametrization of Section (3.2), that is, introducing the line element

[TABLE]

the new radial coordinate $\tau$ such that

[TABLE]

and the differential operator $\bf{e}_{1}=\partial_{\tau}$, the equations for the variables $x,y,p,\phi,K$ become:

[TABLE]

We consider the polytropic equation of state $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$, and assume the scalar field potential to be $W\left(\phi\right)=W_{0}e^{-k\phi}$.

4.2.1 Singularity analysis and integrability

The FE do not pass the Painlevè test, and therefore they are not integrable.

4.2.2 Equilibrium points in the finite region of the phase space

Introducing the variables

[TABLE]

and using the independent dimensionless variable $\lambda$ defined by the choice $N^{2}r^{2}=Y^{2}/(\theta^{2}+W)$, that is

[TABLE]

we obtain the dynamical system

[TABLE]

Additionally,

[TABLE]

Therefore, the system defines a flow on the invariant set

[TABLE]

The system is invariant under the change of coordinates $(X_{\phi},Y_{\phi},S_{1})\rightarrow(-X_{\phi},-Y_{\phi},-S_{1})$ with the simultaneous reversal in the independent variable ${\lambda}\rightarrow-{\lambda}$. In relation to the phase space dynamics this implies that for two points related by this symmetry, say $P^{+}$ and $P^{-}$, each has the opposite dynamical behavior to the other; that is, if the equilibrium point $P^{+}$ is an attractor for a choice of parameters, then $P^{-}$ is a source under the same choice. The details of the stability of the equilibrium points (curves/surfaces of equilibrium points) of the dynamical system (141) are presented in D. These results are summarized in Tables 8 and 9.

We have the useful $m-r$ relation:

[TABLE]

Using these equation we define the compact variables

[TABLE]

which depend only of the phase-space variables $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)$. Evaluating numerically the expressions $R_{\text{comp}},M_{\text{comp}}$ on the orbits of the system (141), we can determine whether the resulting model leads to finite radius and finite mass. To complete this section we make some numerical integrations of the system (141). In the figure 17 are presented some 2D projections of the orbits of the system (141) for $\texttt{b}=1,n=6,k=-3$. In the figure 18 are presented some 3D projections of the orbits of the system (141) for the same choice of parameters.

4.2.3 Physical discussion

According to our analytical results, the points with the 5D unstable or stable manifolds are $P_{5}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}0,\varepsilon,\frac{3\varepsilon}{4\texttt{b}},0,1\Big{)}, exists for $\texttt{b}<0$, or $\texttt{b}\geq\frac{9}{16}$, and $P_{7}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\frac{k\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{2\sqrt{2}\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,0,1\Big{)}. $P_{5}^{+}$ (resp. $P_{5}^{-}$) is a sink (resp. a source) for $\texttt{b}<0,n>0$ or a saddle otherwise. $P_{7}^{+}$ (resp. $P_{7}^{-}$) is a source (resp. a sink) for $k\leq-2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},n>0$, or $-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},n>0$, or $2\sqrt{\frac{2}{3}}<k<2\sqrt{2},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},n>0$, or $k\geq 2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},n>0$. In the simulations presented in Figs. 17, 18 the attractor as $\lambda\rightarrow+\infty$ is thereof the equilibrium point $P_{7}^{-}$, and as $\lambda\rightarrow-\infty$ the past attractor is $P_{7}^{+}$. Comparing the models 3 and 4 it can be argued that the dynamics of the scalar field is qualitatively very similar for both linear and polytropic EoS.

5 Conclusions

In this paper we have investigated the FE in the Einstein-æther model in a static spherically symmetric spacetime. We have studied if the gravitational FE posses the Painlevè property, and consequently whether an analytic explicit integration can be performed for the FE. We have applied the classical treatment for the singularity analysis as summarized in the ARS algorithm (see subsection 1.3). As far as the dynamical system with only the perfect fluid is concerned, we have explicitly shown in subsection 3.1.1 that there always exist a negative resonance which means that the Laurent expansion is expressed in a Left and Right Painlevé Series, because we have integrated over annulus around the singularity which has two borders [57]. Therefore, when the scalar field is not present, and assuming the perfect fluid has the EoS $\mu=\mu_{0}+\left(\texttt{h}-1\right)p,~{}\texttt{h}>1$, the system passes the singularity test which means that the FE are integrable (in a similar way to what we have done for Kantowski-Sachs Einstein-æther theories in [19]). The very specific case when $\mu_{0}=0,\texttt{h}=n+1,\texttt{b}=1$, is exactly the example discussed in section 1.4 and studied in [45], that it is Liouville-integrable if and only if $n\in\{-1,-3,-5,-6\}$. Therefore, our system is integrable, but not Liouville-integrable due to the hypothesis $\texttt{h}>1$ implies $n\notin\{-1,-3,-5,-6\}$. On the other hand, in the presence of the scalar field, as shown earlier in subsections 4.1.1 and 4.2.1, the FE do not generally pass the Painlevè test, except in the special case of a linear EoS with parameter $\mu_{0}=0$.

Furthermore, in the section 3 we have applied the Tolman-Oppenheimer-Volkoff (TOV) approach [37, 38, 39] and found that the usual relativistic TOV equations are then drastically modified in Einstein-æther theory for a perfect fluid to (39), where $\texttt{b}=c_{a}+1$, and $c_{a}$ is the æther parameter. This modification has several physical implication that we have discussed earlier in subsections 3.1.4, and 3.2.4. Then we have constructed a 3D dynamical system in compact variables and obtained a global picture of the entire solution space for different choices of the EoS that can be visualized in a geometrical way. For higher dimensional systems we can still obtain pictorial information by numerical integrations and the use of projections like in sections 4.1.3 and 4.2.3. The results obtained have been discussed physically in each section where they are derived (see discussions in sections 3.1.4, 3.2.4, 4.1.3 and 4.2.3), and we have thus obtained an appropriate description of the universe both on local and larger scales for the models under investigation.

By considering a static spherically symmetric spacetime with a perfect fluid with a linear equation of state (Model 1, Section 3.1) we have deduced the following qualitative physical results (this is brief summary: the details are given in subsection 3.1.4):

• (i)

We derived the system (57). In this model the line of equilibrium points $L_{1}:(S_{1},S_{2},Y)=\left({S_{1}^{*}},1-\texttt{b}{S_{1}^{*}}^{2},0\right)$ is stable for $S_{1}^{*}>0$ or unstable for $S_{1}^{*}<0$. The point is $P_{4}:\left(\frac{\texttt{h}+2}{4\texttt{b}},0,\frac{1}{\texttt{h}}\right)$ is a sink for $\texttt{b}<0,\texttt{h}>1$.

• (ii)

For the low pressure regime governed by the equations (78), the line $L_{6}:(U,V,Y)=(U_{0},1,0)$ is a local attractor. The orbit connecting the so-called Tolman point $L_{2}(0):(U,V,Y)=\left(\frac{3}{4},0,0\right)$ with the line $L_{6}$ represents the so-called Tolman orbit (associated with the positive eigenvalue of $L_{2}(0)$). When $\texttt{b}<0$ (Einstein- æther modification) $L_{1}(0):(U,V,Y)=\left(1,0,0\right)$ can be an attractor, whereas for $\texttt{b}>0$, we have the same qualitatively behavior as for the model (78).

• (iii)

In the high pressure regime for $\texttt{h}>1$, and for GR ($\texttt{b}=1$), governed by the system (80), the attractor satisfies: $(U,V,Y)=\left(\frac{1}{2},\frac{2(\texttt{h}-1)^{2}}{3\texttt{h}^{2}-2},\frac{1}{\texttt{h}}\right)$. It is a stable node for $1<\texttt{h}\leq\frac{2}{7}\left(4\sqrt{7}-7\right)$ or a stable spiral for $\texttt{h}>\frac{2}{7}\left(4\sqrt{7}-7\right)$. This point represents solutions that have finite mass and radius.

• (iv)

Implementing the procedure described in the A.2 we conclude that there are not physical relevant equilibrium points at infinity, apart of the point $I_{4}$ that has $S_{1}\rightarrow\infty$, $S_{2}$ finite and $Y=\frac{1}{\texttt{h}}$, which is a source (it exists only for $\texttt{b}<0$). This solution corresponds to $y$ finite and $x\rightarrow y$. That is, $\frac{dp}{d\rho}=\frac{1}{\mu+p},\frac{p}{p+\mu}=\frac{1}{\texttt{h}},$ where $\rho$ is the radial coordinate. Then $\frac{dp}{d\rho}=\frac{1}{\texttt{h}p}\implies p(\rho)=\frac{\sqrt{\texttt{h}p_{c}^{2}+2\rho}}{\sqrt{\texttt{h}}},$ where $p_{c}$ is the central pressure (at $\rho=0$). As $\rho\rightarrow\infty$ the pressure is increasing as the solutions move away from this point.

We have also studied the case of a polytropic equation of state $p=q\mu^{\Gamma},\Gamma=1+\frac{1}{n},n>0$ in the static Eistein-æther theory, by means of the $(\Sigma,S_{2},Y)$-formulation given by the the model (93) (Model 2, Section 3.2). The main results are summarized as follows (this is brief summary: the details are given in subsection 3.2.4):

• (i)

Using the $(S_{1},S_{2},Y)$ formulation we deduce the system (93), and when $\texttt{b}\neq 0$ one branch of $L_{1}$ is a local source and the other branch is a local sink. The line $P_{2}$ is a local source. The self-similar plane symmetric solution $P_{4}$ is a saddle (the orbits coming from this point are associated with a negative mass singularity).

• (ii)

The orbits on the upper half $S_{1}>0$ are attracted by either $L_{1}$ or $L_{2}$. When $S_{1}<0$, the line $L_{1}$ is the local source and $L_{2}$ acts as the local sink. Each point located on $L_{2}$ corresponds to the flat Minkowski solution written in spherically symmetric form. An orbit parametrized by $y_{c}=p_{c}/(p_{c}+\mu_{c})=q\mu_{c}^{1/n}/(1+q\mu_{c}^{1/n})$, where $p_{c}$ and $\mu_{c}$ are the central pressure and central density, enters the interior of the phase space from each point on $L_{2}$.

• (iii)

In GR, the point $(S_{1},S_{2},Y)=\left(\frac{2}{3},\frac{1}{9},1\right)$, corresponds to the so called self-similar Tolman solution discussed in [36]. In the Einstein-æther theory this solution is generalized to a 1-parameter set of solutions given by the equilibrium points $P_{1}:(S_{1},S_{2},Y)=\left(\frac{2}{3},\frac{1}{9}(9-8\texttt{b}),1\right)$. The single orbit that enters the interior of the phase space from $P_{1}$, associated to the positive eigenvalue of $P_{1}$ is called the Tolman orbit. The solutions having $Y=1$ corresponds to the limiting situation $\frac{p}{\mu}\rightarrow\infty.$ In the case of a linear equation of state $p=(\gamma-1)\mu$, this is equivalent to $\gamma\rightarrow\infty$.

• (iv)

Regarding the singular curve $L_{1}$, the left branch $S_{1}<0$ is unstable and the right branch $S_{1}>0$ is stable. Evaluating at $L_{1}$ we have $M/R=\Sigma_{c}(1-(1+\texttt{b})\Sigma_{c}/2)/(1-\texttt{b}\Sigma_{c}^{2})$. Thus, the solutions ending on $L_{1}$ when $\Sigma_{c}>0$ have finite radii and masses. The orbits ending on $L_{1}$ when $\Sigma_{c}=0$ or $\Sigma_{c}=2/(1+\texttt{b})$ describe solutions with infinite radii and finite masses, or masses that approach infinity slower than $r$.

• (v)

Regarding the behavior of solutions at infinity, we refer to the discussion in B.2. In particular, for $\texttt{b}<0$, there are not physical relevant equilibrium points at infinity, apart of the point $I_{4}$, that has $S_{1}\rightarrow\infty$, $S_{2}$ finite and $Y=1$, which is a source. This solution corresponds to $y$ finite and $x\rightarrow y$. Therefore, we obtain $p(\rho)=\sqrt{p_{c}^{2}+2\rho},$ where $p_{c}$ is the central pressure (at $\rho=0$), and the pressure is increasing as $\rho\rightarrow\infty$.

The system (101) corresponds to GR, and it was comprehensively studied in reference [41]. Summarizing the results of [41] (see earlier subsection 3.2.4 for details), for $n\geq 0$ there exists regular Newtonian solutions, and when $0\leq n<5$ it is possible to obtain solutions that represents star models with finite radii and mass, while for $n\geq 5$ infinite regular models are obtained. In more detail,

• (i)

$n=0$ corresponds to relativistic incompressible fluid models, and the set $U=\frac{3}{4}$ contains all regular solutions,

• (ii)

for $0<n\leq 3$, the point $P_{2}$ is a hyperbolic sink, and due to the existence of a monotonic function in the subset $Y=0$, $P_{2}$ is consequently a global sink, so that the radii and masses of a regular models are finite for $0<n\leq 3$,

• (iii)

for $3<n<5$, the point $P_{2}$ is the only sink, but not a global one, and so, the solutions do not end necessarily at $P_{2}$,

• (iv)

when $n=5$, a separatrix surface completely encloses the regular subset of orbits, and prevents any regular orbit ending at $P_{2}$ (the Newtonian regular orbit ends on $P_{1}$, while all general relativistic regular orbits tends to $P_{4}$ and to the closed solutions surrounding $P_{4}$),

• (v)

for $n>5$, a separatrix surface encloses all regular solutions (including Newtonian orbits), and all regular solutions end at the hyperbolic sink $P_{4}$, and so, Newtonian and relativistic polytropic models have infinite radii and mass.

Finally, we have also naturally extended of the previous analysis in the GR static case [41], by studying a stationary comoving æther with a perfect fluid with a linear EoS (Model 3, Section 4.1), and with a polytropic EoS (Model 4, Section 4.2), in addition to a scalar field with an exponential potential.

First we show that the singularity test is not passed, except in the special case of a linear EoS with parameter $\mu_{0}=0$ (subsections 4.1.1 and 4.2.1). We then performed a phase space analysis and presented all of the equilibrium points in the cases of both a linear EoS (see subsection 4.1.2) and a polytropic EoS (see subsection 4.2.2), and described the corresponding exact solutions and discussed their physical properties in the appendix. We note that the sources and sinks in the linear EoS are the pairs $P_{5}^{\varepsilon}$ and $P_{7}^{\varepsilon}$ for appropriate values of the model parameters (and the corresponding pairs for the polytropic EoS case studied subsequently), which all have $S_{1}=\frac{K}{\theta^{2}+W}=0$. In particular in the case of $P_{5}^{\varepsilon}$ (where $C=0$ and $Q=\pm 1$) the asymptotic metric, corresponding a a self-similar plane symmetric model containing a singularity with a non-trivial horizon, was explicitly displayed; and the case of $P_{7}^{\varepsilon}$, that can be asymptotically a Minkowski solution or a singularity, depending on the parameters.

Acknowledgments

G. L. was funded by Comisión Nacional de Investigación Científica y Tecnológica (CONICYT) through FONDECYT Iniciación grant no. 11180126 and by Vicerrectoría de Investigación y Desarrollo Tecnológico at Universidad Católica del Norte. A. C. was supported by NSERC of Canada and partially supported by FONDECYT Iniciación grant no. 11180126. The authors would like to thank Prof. PGL Leach and to J. Latta for useful comments and suggestions.

Appendix A Equilibrium points for Model 1

A.1 Equilibrium points of Model 1 at the finite part of the phase space $(S_{1},S_{2},Y)$

The equilibrium points and the corresponding eigenvalues of the dynamical system (57) are the following:

$L_{1}$: $(S_{1},S_{2},Y)=\left({S_{1}^{*}},1-\texttt{b}{S_{1}^{*}}^{2},0\right)$. Exists for $\texttt{b}\leq 0,S_{1}^{*}\in\mathbb{R}$, or $\texttt{b}>0,-\sqrt{\frac{1}{\texttt{b}}}\leq S_{1}^{*}\leq\sqrt{\frac{1}{\texttt{b}}}$. The eigenvalues are $0,-{S_{1}^{*}},-{S_{1}^{*}}$. Nonhyperbolic with a 2D stable (resp. unstable) manifold for $S_{1}^{*}>0$ (resp. $S_{1}^{*}<0$). 2. 2.

$L_{2}$: $(S_{1},S_{2},Y)=(0,1,Y_{0})$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are $0,-Y_{0},2Y_{0}$. It is a saddle. 3. 3.

$P_{1}$: $(S_{1},S_{2},Y)=\left(\frac{2}{\texttt{h}+2},1-\frac{8\texttt{b}}{(\texttt{h}+2)^{2}},\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{h}\geq 1,0\leq\texttt{b}\leq\frac{1}{8}(\texttt{h}+2)^{2}$. The eigenvalues are $\frac{2}{\texttt{h}+2},-\frac{\texttt{h}+2+\sqrt{64\texttt{b}-7(\texttt{h}+2)^{2}}}{2\texttt{h}(\texttt{h}+2)},-\frac{\texttt{h}+2-\sqrt{64\texttt{b}-7(\texttt{h}+2)^{2}}}{2\texttt{h}(\texttt{h}+2)}$. It is a saddle. 4. 4.

$P_{2}$: $(S_{1},S_{2},Y)=\left(\frac{1}{\sqrt{\texttt{b}}},0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}>0$. The eigenvalues are $\frac{2}{\texttt{h}},\frac{-\texttt{h}+4\sqrt{\texttt{b}}-2}{\sqrt{\texttt{b}}\texttt{h}},\frac{1}{\sqrt{\texttt{b}}}$. It is a source for $\texttt{h}>1,\texttt{b}>\frac{1}{16}\left(\texttt{h}+2\right)^{2}$, saddle otherwise. 5. 5.

$P_{3}$: $\left(-\frac{1}{\sqrt{\texttt{b}}},0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}>0$. The eigenvalues are $\frac{2}{\texttt{h}},-\frac{1}{\sqrt{\texttt{b}}},\frac{\texttt{h}+4\sqrt{\texttt{b}}+2}{\sqrt{\texttt{b}}\texttt{h}}$. It is a saddle. 6. 6.

$P_{4}$: $(S_{1},S_{2},Y)=\left(\frac{\texttt{h}+2}{4\texttt{b}},0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}<0,\texttt{h}\geq 1$ or $\texttt{h}\geq 1,\texttt{b}\geq\frac{1}{16}(\texttt{h}+2)^{2}$. The eigenvalues are

$\frac{\texttt{h}+2}{4\texttt{b}},\frac{(\texttt{h}+2)^{2}-8\texttt{b}}{4\texttt{b}\texttt{h}},\frac{(\texttt{h}+2)^{2}-16\texttt{b}}{8\texttt{b}\texttt{h}}$. It is nonhyperbolic for $\texttt{b}=\frac{1}{16}\left(\texttt{h}+2\right)^{2}$. Is is a sink for $\texttt{b}<0,\texttt{h}>1$.

A.2 Equilibrium points of Model 1 at infinity.

In this section we cover the case $\texttt{b}<0$, which gives an unbounded phase space in which $S_{1},S_{2}$ can take infinite values. Due to in $S_{2}\leq 1-\texttt{b}S_{1}^{2},\texttt{b}<0$, $1-\texttt{b}S_{1}^{2}$ is the dominant quantity, and we define the normalized variables and radial derivative:

[TABLE]

to obtain the dynamical system

[TABLE]

defined in the compact phase space

[TABLE]

Now we discuss about the equilibrium points at infinity of the system (57) in the above compact coordinates:

1. $I_{1}$

: $(\Pi_{1},\Pi_{2},Y)=(-1,0,0)$, with eigenvalues $-1,-1,1$. Therefore, it is a saddle. 2. $I_{2}$

: $(\Pi_{1},\Pi_{2},Y)=(1,0,0)$, with eigenvalues $-1,1,1$. Therefore, it is a saddle. 3. $I_{3}$

: $(\Pi_{1},\Pi_{2},Y)=\left(-1,0,\frac{1}{\texttt{h}}\right)$, with eigenvalues $-\frac{\texttt{h}+2}{\texttt{h}},-1,\frac{\infty\sqrt{\left|\texttt{b}\right|}}{\texttt{h}}$. Therefore, it is a saddle. 4. $I_{4}$

: $(\Pi_{1},\Pi_{2},Y)=\left(1,0,\frac{1}{\texttt{h}}\right)$, with eigenvalues $1,\frac{\texttt{h}+2}{\texttt{h}},\frac{\infty\sqrt{\left|\texttt{b}\right|}}{\texttt{h}}$. Therefore, it is a source. 5. $I_{5}$

: $(\Pi_{1},\Pi_{2},Y)=\left(-1,1,\frac{1}{\texttt{h}}\right)$, with eigenvalues $-1,\frac{\texttt{h}+2}{\texttt{h}},\frac{\infty\sqrt{\left|\texttt{b}\right|}}{\texttt{h}}$. Therefore, it is a saddle. 6. $I_{6}$

: $(\Pi_{1},\Pi_{2},Y)=\left(1,1,\frac{1}{\texttt{h}}\right)$, with eigenvalues $-\frac{\texttt{h}+2}{\texttt{h}},1,\frac{\infty\sqrt{\left|\texttt{b}\right|}}{\texttt{h}}$. Therefore, it is a saddle.

A.3 Equilibrium points of Model 1 using homology invariant variables $(U,V,Y)$ in the case of General Relativity ($\texttt{b}=1$).

The equilibrium points and curves of equilibrium points of (75) are summarized below.

$L_{1}$, exists for $0\leq Y_{0}\leq\frac{1}{\texttt{h}}$ (analogous to line $L_{1}$ characterized in table II of [41]). It is a source. 2. 2.

$L_{2}$, exists for $0\leq Y_{0}\leq\frac{1}{\texttt{h}}$ (analogous to line $L_{2}$ characterized in table II of [41]). it is a saddle. 3. 3.

$L_{3}$, exists for $0\leq U_{0}\leq 1,\texttt{h}=1$ (analogous to line $L_{3}$ studied in table II of [41]). It is non-hyperbolic. 4. 4.

$L_{4}$, exists for $0\leq Y_{0}\leq\frac{1}{\texttt{h}}$ (analogous to line $L_{4}$ characterized in table II [41]). It is a saddle. 5. 5.

The system (75) admits the line of equilibrium points $L_{5}$ (which is the analogue of $L_{5}$ characterized in table II of [41]). It is non-hyperbolic. From the assumption $0\leq V\leq(1-Y)/(1+Y)$ the static solutions are forced to be away the line $L_{5}$. 6. 6.

$L_{6}$, exists for $0\leq U_{0}\leq 1$ (analogous to line $L_{6}$ studied in [41]). The center manifold of this line is tangent to the eigenvector $\left(1,0,0\right)^{T}$, such that the line is normally hyperbolic. Due to the negativity of the non-zero eigenvalues it is a sink. 7. 7.

$M_{1}$, always exists. It is the analogous to the point $P_{1}$ presented in table II of [41] by choosing $n=0$. It is a saddle. 8. 8.

$M_{2}$, always exists. This point is the analogous to the point $P_{2}$ showed in table II of [41]. It is a special point of the line $L_{6}$ when $U_{0}=0$. It is a sink. 9. 9.

$M_{3}$ is the analogous to the point $P_{3}$ showed in table II of [41]. It is a saddle. This point is of little interest for the static solutions. 10. 10.

$M_{5}$, always exists. It is the analogous to the point $P_{5}$ showed in table II of [41] (it is contained in the attractor line $L_{6}$ described below). It is nonhyperbolic. 11. 11.

The points $M_{6},M_{7}$, and $M_{8}$ only exists for linear equations of state as we have considered in our model. They are always saddle.

There is no equilibrium point of (75) analogous to $P_{4}$, presented in table II of [41].

Finally, there are three equilibrium points with $V=1,Y=\frac{1}{\texttt{h}}$, which are of little interest for the static solutions, because they are not over the “static surface”. Strictly speaking since the last inequality defining the phase space (75d) is not fulfilled at these points they do not exist. Therefore, we omit any discussion of them.

A.4 Equilibrium points of Model 1 using homology invariant variables $(U,V,Y)$ for the Einstein- æther modification ($\texttt{b}\neq 1$).

The equilibrium points/lines of equilibrium points of the system (81), and the corresponding eigenvalues, are enumerated as follows (see Table 3):

$L_{1}$: The line of equilibrium points $(U,V,Y)=\left(1,0,Y_{0}\right)$, with eigenvalues

$0,-\frac{(Y_{0}-1)^{2}((3\texttt{b}-2)Y_{0}-1)}{\texttt{b}},-\frac{(Y_{0}-1)^{2}((3\texttt{b}-2)Y_{0}-1)}{\texttt{b}}$. This a line that it is an attractor for $\texttt{h}>1,\texttt{b}<0,0<Y_{0}<\frac{1}{\texttt{h}}$, or $\texttt{h}>1,\texttt{b}>\frac{\texttt{h}+2}{3},\frac{1}{3\texttt{b}-2}<Y_{0}<\frac{1}{\texttt{h}}$. It is a source for $\texttt{h}>1,0<\texttt{b}\leq\frac{\texttt{h}+2}{3},0<Y_{0}<\frac{1}{\texttt{h}}$, or $\texttt{h}>1,\texttt{b}>\frac{\texttt{h}+2}{3},0<Y_{0}<\frac{1}{3\texttt{b}-2}$. 2. 2.

$L_{1}(1/\texttt{h})$: The equilibrium point $(U,V,Y)=\left(1,0,\frac{1}{\texttt{h}}\right)$, with eigenvalues

$0,-\frac{(\texttt{h}-1)^{2}(3\texttt{b}-\texttt{h}-2)}{\texttt{b}\texttt{h}^{3}},-\frac{(\texttt{h}-1)^{2}(3\texttt{b}-\texttt{h}-2)}{\texttt{b}\texttt{h}^{3}}$. Nonhyperbolic with a 2D unstable manifold for $\texttt{h}>1,0<\texttt{b}<\frac{\texttt{h}+2}{3}$. Nonhyperbolic with a 2D stable manifold for $\texttt{h}>1,\texttt{b}<0$, or $\texttt{h}>1,\texttt{b}>\frac{\texttt{h}+2}{3}$. 3. 3.

$L_{2}$: The equilibrium point $(U,V,Y)=\left(\frac{3}{4},0,Y_{0}\right)$. Exists only for $\texttt{b}=1$. The eigenvalues are $0,-\frac{1}{2}(Y_{0}-1)^{3},\frac{3}{4}(Y_{0}-1)^{3}$. Nonhyperbolic and it behaves as a saddle. 4. 4.

$L_{2}(0)$: The equilibrium point $(U,V,Y)=\left(\frac{3\texttt{b}}{1+3\texttt{b}},0,0\right)$. The eigenvalues are $0,-\frac{3}{3\texttt{b}+1},\frac{2}{3\texttt{b}+1}$. Nonhyperbolic and it behaves as a saddle. 5. 5.

$L_{2}(1)$: The equilibrium point $(U,V,Y)=\left(\frac{3\texttt{b}}{1+3\texttt{b}},0,1\right)$. The eigenvalues are $0,0,0$. It is nonhyperbolic. 6. 6.

$L_{3}$: The line of equilibrium points $(U,V,Y)=\left(U_{0},0,1\right)$. The eigenvalues are $0,0,0$. It is nonhyperbolic. 7. 7.

$L_{4}$: The lines of equilibrium point $(U,V,Y)=\left(0,0,Y_{0}\right)$, with eigenvalues $0,-3(Y_{0}-1)^{3},(Y_{0}-1)^{3}$. Nonhyperbolic for $Y_{0}=0$, saddle otherwise. 8. 8.

$L_{5}(1)$: The equilibrium point $(U,V,Y)=\left(0,1,1\right)$, with eigenvalues $0,0,0$. It is nonhyperbolic. 9. 9.

$L_{6}$: The line of equilibrium points $(U,V,Y)=\left(U_{0},1,0\right)$. The eigenvalues are

$0,-\frac{2\left(\sqrt{\texttt{b}}+1\right)(U_{0}-1)}{\texttt{b}-1},-\frac{2\left(\sqrt{\texttt{b}}+1\right)(U_{0}-1)}{\texttt{b}-1}.$ This line of equilibrium points is nonhyperbolic with a 2D stable manifold for $0<U_{0}<1,0\leq\texttt{b}<1$ (two negative eigenvalues) or $\texttt{b}<0,0<U_{0}<1$ (two complex conjugated eigenvalues with negative real parts). This line of equilibrium points is nonhyperbolic with a 2D unstable manifold for $0<U_{0}<1,\texttt{b}>1$. 10. 10.

$M_{1}$: The equilibrium point $(U,V,Y)=\left(0,\frac{1}{2},0\right)$, with eigenvalues $\frac{3}{2},\frac{1}{2},-\frac{1}{2}$. Saddle. 11. 11.

$M_{2}$: The equilibrium point $(U,V,Y)=\left(0,1,0\right)$, with eigenvalues

$\frac{2\left(1+\sqrt{\texttt{b}}\right)}{\texttt{b}-1},0,\frac{2\left(1+\sqrt{\texttt{b}}\right)}{\texttt{b}-1}$. Nonhyperbolic with a 2D stable manifold for $0<\texttt{b}<1$ (two negative real eigenvalues) or $\texttt{b}<0$ (two complex conjugated eigenvalues). 12. 12.

$M_{3}$: The equilibrium point $(U,V,Y)=\left(1,1,1\right)$, with eigenvalues $0,0,0$. It is nonhyperbolic. 13. 13.

$M_{5}$: The equilibrium point $(U,V,Y)=\left(1,1,0\right)$. The eigenvalues at the equilibrium points are $0,0,0$. Therefore, it is nonhyperbolic. However, taking $Y=\varepsilon\ll 1$, instead of exactly $Y=0$, we obtain the eigenvalues $0,-\frac{2\varepsilon}{\texttt{b}},\frac{2\varepsilon}{\texttt{b}}$, that is, the equilibrium point is indeed a saddle point. 14. 14.

$M_{6}$: The equilibrium point $(U,V,Y)=\left(0,\frac{\texttt{h}-1}{\texttt{h}+1},\frac{1}{\texttt{h}}\right)$, with eigenvalues

$0,0,\frac{4(\texttt{h}-1)^{2}}{\texttt{h}^{2}(\texttt{h}+1)}$. Non hyperbolic with a 1D unstable manifold. 15. 15.

$M_{7}$: The equilibrium point $(U,V,Y)=\left(1,\frac{\texttt{h}-1}{\texttt{h}+1},\frac{1}{\texttt{h}}\right)$, with eigenvalues

$0,0,\frac{2(\texttt{h}-1)(4\texttt{b}-\texttt{h}-2)}{\texttt{b}\texttt{h}^{2}(\texttt{h}+1)}$. Non hyperbolic with a 1D unstable manifold for $\texttt{b}<0,\texttt{h}>1$, or $\texttt{b}>\frac{3}{4},1<\texttt{h}<4\texttt{b}-2$. Nonhyperbolic with a 1D stable manifold for $0<\texttt{b}\leq\frac{3}{4},\texttt{h}>1$, or $\texttt{b}>\frac{3}{4},\texttt{h}>4\texttt{b}-2$.

Furthermore, there are four equilibrium points with $V=1,Y=\frac{1}{\texttt{h}}$, which are of little interest for the static solutions, because they do not satisfy the inequalities $0\leq V\leq\frac{1-Y}{1+Y},\texttt{h}\geq 1.$ We omit the discussion of them.

Finally, we have found new lines of/equilibrium points in comparison with the previous case (75) (see table 2), like:

$M_{9}$: The line of equilibrium points $(U,V,Y)=\left(1,1,Y_{0}\right)$. The eigenvalues are

$0,\frac{2Y_{0}((2-4\texttt{b})Y_{0}+1)}{\texttt{b}},-\frac{2Y_{0}((2-4\texttt{b})Y_{0}+1)}{\texttt{b}}$. This line is nonhyperbolic and it behaves as a saddle. 2. 2.

$M_{10}$: The line of equilibrium points $(U,V,Y)=\left(1,V_{0},\frac{1-V_{0}}{1+V_{0}}\right)$. This line of equilibrium points has the eigenvalues $\frac{2(V_{0}-1)V_{0}(4\texttt{b}(V_{0}-1)-V_{0}+3)}{\texttt{b}(V_{0}+1)^{2}},0,0$. It is nonhyperbolic. 3. 3.

$M_{11}$: The line of equilibrium points $(U,V,Y)=\left(1,V_{0},\frac{1}{4\texttt{b}-2}\right)$. The eigenvalues are $0,0,0$. It is nonhyperbolic. 4. 4.

$M_{12}$: The line of equilibrium points $(U,V,Y)=\left(U_{0},\frac{\texttt{h}-1}{\texttt{h}+1},\frac{1}{\texttt{h}}\right)$. The eigenvalues are

$0,0,-\frac{2(\texttt{h}-1)^{2}(U_{0}(2\texttt{b}(\texttt{h}-3)+\texttt{h}+2)-2\texttt{b}(\texttt{h}-1))}{\texttt{b}\texttt{h}^{3}(\texttt{h}+1)}$. It is nonhyperbolic.

Appendix B Equilibrium points for Model 2

B.1 Equilibrium points of Model 2 at the finite part of the phase space

$(S_{1},S_{2},Y)$

The dynamical system (93) has the (lines of) equilibrium points at the finite region:

Line of equilibrium points: $L_{1}:=\left(\Sigma_{c},1-\texttt{b}\Sigma_{c}^{2},0\right)$. Exists for $\texttt{b}\leq 0,\Sigma_{c}\in\mathbb{R}$, or $\texttt{b}>0,-\sqrt{\frac{1}{\texttt{b}}}\leq\Sigma_{c}\leq\sqrt{\frac{1}{\texttt{b}}}$. Eigenvalues: $0,-\frac{\Sigma_{c}}{1+n},-\Sigma_{c}$. Stable 2D manifold for $n>0,\texttt{b}\leq 0,\Sigma_{c}>0$, or $n>0,\texttt{b}>0,0<\Sigma_{c}\leq\sqrt{\frac{1}{\texttt{b}}}$. Unstable 2D manifold for $n>0,\texttt{b}\leq 0,\Sigma_{c}<0$, or $n>0,\texttt{b}>0,-\sqrt{\frac{1}{\texttt{b}}}\leq\Sigma_{c}<0$. This line contains as particular case the line $L_{1}$ characterized in Table I of [41] when $\texttt{b}=1$. 2. 2.

Line of equilibrium points: $L_{2}:=\left(0,1,y_{c}\right)$. Eigenvalues: $0,-y_{c},2y_{c}$. Saddle. This line coincides with the line $L_{2}$ characterized in Table I of [41]. 3. 3.

$P_{1}:=\left(\frac{2}{3},\frac{1}{9}(9-8\texttt{b}),1\right)$. Exists for $n>0,0\leq\texttt{b}\leq\frac{9}{8}$.

Eigenvalues: $\frac{2}{3(n+1)},\frac{1}{6}\left(-3-\sqrt{64\texttt{b}-63}\right),\frac{1}{6}\left(-3+\sqrt{64\texttt{b}-63}\right)$. Nonhyperbolic for $\texttt{b}=\frac{9}{8}$. Saddle (2D stable manifold for $n>0,0\leq\texttt{b}<\frac{9}{8}$). This point reduces to $P_{1}$ characterized in Table I of [41] when $\texttt{b}=1$. 4. 4.

$P_{2}:=\left(\frac{1}{\sqrt{\texttt{b}}},0,1\right)$. Exists for $\texttt{b}>0,n>0$. Eigenvalues: $\frac{1}{\sqrt{\texttt{b}}(n+1)},2,\frac{4\sqrt{\texttt{b}}-3}{\sqrt{\texttt{b}}}$. Nonhyperbolic for $\texttt{b}=\frac{9}{16}$. Unstable for $n>0,\texttt{b}>\frac{9}{16}$. Saddle, otherwise. This point reduces to $P_{2}$ characterized in Table I of [41] when $\texttt{b}=1$. 5. 5.

$P_{3}:=\left(-\frac{1}{\sqrt{\texttt{b}}},0,1\right)$. Exists for $\texttt{b}>0,n>0$. Eigenvalues: $-\frac{1}{\sqrt{\texttt{b}}(n+1)},2,\frac{4\sqrt{\texttt{b}}+3}{\sqrt{\texttt{b}}}$. Saddle. This point reduces to $P_{3}$ characterized in Table I of [41] when $\texttt{b}=1$. 6. 6.

$P_{4}:=\left(\frac{3}{4\texttt{b}},0,1\right)$. Exists for $\texttt{b}<0,n>0$, or $\texttt{b}\geq\frac{9}{16},n>0$. Eigenvalues:

$\frac{3}{4\texttt{b}(n+1)},-\frac{8\texttt{b}-9}{4\texttt{b}},-\frac{16\texttt{b}-9}{8\texttt{b}}$. Nonhyperbolic for $\texttt{b}=\frac{9}{8}$, or $\texttt{b}=\frac{9}{16}$. Stable for $n>0,\texttt{b}<0$. Saddle, otherwise. This point reduces to $P_{4}$ characterized in Table I of [41] when $\texttt{b}=1$.

B.2 Equilibrium points of Model 2 at infinity

We discussed before the phase space is bounded for $\texttt{b}>0$, such that in this section we cover the case $\texttt{b}<0$ which gives an unbounded phase space. $Y\in[0,1]$ is bounded but $S_{1},S_{2}$ can take infinite values, so that we may define a coordinate system in which the coordinates consist of $Y$, and the bounded variables $\Pi_{1},\Pi_{2}$, and the radial derivative $\bar{\lambda}$ defined in (146) to obtain the dynamical system

[TABLE]

defined in the compact phase space

[TABLE]

Now we discuss about the equilibrium points at infinity of the system (93) in the above compact coordinates:

1. $I_{1}$

: $(\Pi_{1},\Pi_{2},Y)=(-1,0,0)$, with eigenvalues $\frac{1}{n+1},-1,-1$. It is a saddle. 2. $I_{2}$

: $(\Pi_{1},\Pi_{2},Y)=(1,0,0)$, with eigenvalues $-\frac{1}{n+1},1,1$. It is a saddle. 3. $I_{3}$

: $(\Pi_{1},\Pi_{2},Y)=\left(-1,0,1\right)$, with eigenvalues $-\frac{1}{n+1},-3,\infty\sqrt{\left|\texttt{b}\right|}$. It is a saddle. 4. $I_{4}$

: $(\Pi_{1},\Pi_{2},Y)=\left(1,0,1\right)$, with eigenvalues $\frac{1}{n+1},3,\infty\sqrt{\left|\texttt{b}\right|}$. It is a source. 5. $I_{5}$

: $(\Pi_{1},\Pi_{2},Y)=\left(-1,1,1\right)$, with eigenvalues $-\frac{1}{n+1},3,\infty\sqrt{\left|\texttt{b}\right|}$. It is a saddle. 6. $I_{6}$

: $(\Pi_{1},\Pi_{2},Y)=\left(1,1,1\right)$, with eigenvalues $\frac{1}{n+1},-3,\infty\sqrt{\left|\texttt{b}\right|}$. It is a saddle.

B.3 Equilibrium points of Model 2 using homology invariant variables $(U,V,Y)$ in the case of General Relativity ($\texttt{b}=1$)

The equilibrium points of the system (101) are (they were found in the reference [41], we use here the same notation):

$L_{1}:(U,V,Y)=(1,0,Y_{0})$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are $0,(Y_{0}-1)^{2},(Y_{0}-1)^{2}$. It is a source. 2. 2.

$L_{2}:(U,V,Y)=\left(\frac{3}{4},0,Y_{0}\right)$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are $0,-\frac{3}{4}(Y_{0}-1)^{2},\frac{1}{2}(Y_{0}-1)^{2}$. It is nonhyperbolic, behaves as saddle. 3. 3.

$L_{3}:(U,V,Y)=(U_{0},0,1)$. Exists for $0\leq U_{0}\leq 1$. The eigenvalues are $0,0,0$. It is onhyperbolic. 4. 4.

$L_{4}:(U,V,Y)=(0,0,Y_{0})$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are $0,-(Y_{0}-1)^{2},3(Y_{0}-1)^{2}$. It is nonhyperbolic, behaves as saddle. 5. 5.

$L_{5}:(U,V,Y)=(0,V_{0},1)$. Exists for $-1\leq V_{0}\leq 1$. The eigenvalues are $0,0,0$. It is nonhyperbolic. 6. 6.

$L_{6}:(U,V,Y)=(U_{0},1,0)$. Exists for $n=0$. The eigenvalues are $0,-(1-U_{0}),-(1-U_{0})$. This line is normally hyperbolic and due to the non-zero eigenvalue is negative it is a sink. 7. 7.

$P_{1}:(U,V,Y)=\left(0,\frac{n+1}{n+2},0\right)$. Exists for $n\geq 0$. The eigenvalues are

$-\frac{1}{n+2},\frac{1}{n+2},-\frac{n-3}{n+2}$. Saddle with a 2D unstable manifold for $0<n<3$. Saddle with a 2D stable manifold for $n>3$. 8. 8.

$P_{2}:(U,V,Y)=(0,1,0)$. Always exists. The eigenvalues are $-\frac{1}{n+1},-\frac{1}{n+1},-\frac{n}{n+1}$. It is a sink. 9. 9.

$P_{3}:(U,V,Y)=(1,1,1)$. Always exists. The eigenvalues are $-\frac{1}{n+1},\frac{1}{n+1},\frac{n}{n+1}$. It is a saddle with a 2D unstable manifold. 10. 10.

$P_{4}:(U,V,Y)=\left(\frac{n-3}{2(n-2)},\frac{2(n+1)}{3n+1},0\right)$. Exists for $n>3$. The eigenvalues are

$-\frac{n-1}{(n-2)(3n+1)},-\frac{(n-1)\left(n-5+\sqrt{1-n(7n-22)}\right)}{4(n-2)(3n+1)},-\frac{(n-1)\left(n-5-\sqrt{1-n(7n-22)}\right)}{4(n-2)(3n+1)}$. Stable for $n>5$ (two complex eigenvalues with negative real parts). 11. 11.

$P_{5}:(U,V,Y)=(1,1,0)$. Always exists. The eigenvalues are $0,0,0$. It is nonhyperbolic.

B.4 Equilibrium points of Model 2 using homology invariant variables $(U,V,Y)$ for the Einstein- æther modification ($\texttt{b}\neq 1$)

The equilibrium points of the system (109), and the corresponding eigenvalues, are the following:

$L_{1}$: $(U,V,Y)=(1,0,Y_{0})$, with eigenvalues

$0,-\frac{(Y_{0}-1)^{2}((3\texttt{b}-2)Y_{0}-1)}{\texttt{b}},-\frac{(Y_{0}-1)^{2}((3\texttt{b}-2)Y_{0}-1)}{\texttt{b}}$. It is an attractor for $\texttt{h}>1,\texttt{b}<0,0<Y_{0}<\frac{1}{\texttt{h}}$, or $\texttt{h}>1,\texttt{b}>\frac{\texttt{h}+2}{3},\frac{1}{3\texttt{b}-2}<Y_{0}<\frac{1}{\texttt{h}}$. It is a source for $\texttt{h}>1,0<\texttt{b}\leq\frac{\texttt{h}+2}{3},0<Y_{0}<\frac{1}{\texttt{h}}$, or $\texttt{h}>1,\texttt{b}>\frac{\texttt{h}+2}{3},0<Y_{0}<\frac{1}{3\texttt{b}-2}$. 2. 2.

$L_{2}$: $(U,V,Y)=\left(\frac{3\texttt{b}}{3\texttt{b}+1},0,Y_{0}\right)$, with eigenvalues

$0,\frac{({Y_{0}}-1)((9\texttt{b}-7){Y_{0}}-2)}{3\texttt{b}+1},\frac{3({Y_{0}}-1)(3\texttt{b}+(9(\texttt{b}-2)\texttt{b}+5){Y_{0}}+1)}{(3\texttt{b}+1)^{2}}$. It is nonhyperbolic with a 2D stable manifold for $\texttt{b}<-\frac{1}{3},\frac{-3\texttt{b}-1}{9\texttt{b}^{2}-18\texttt{b}+5}<{Y_{0}}<1$, or $\texttt{b}>1,\frac{2}{9\texttt{b}-7}<{Y_{0}}<1$. It has a 2D unstable manifold for $\frac{2}{3}<\texttt{b}<1,\frac{-3\texttt{b}-1}{9\texttt{b}^{2}-18\texttt{b}+5}<{Y_{0}}<1$. 3. 3.

$L_{3}$: $(U,V,Y)=(U_{0},0,1)$, with eigenvalues $0,0,0$. It is nonhyperbolic. 4. 4.

$L_{4}$: $(U,V,Y)=(0,0,Y_{0})$, with eigenvalues $0,-(Y_{0}-1)^{2},3(Y_{0}-1)^{2}$. It is a saddle. 5. 5.

$L_{5}$: $(U,V,Y)=(0,V_{0},1)$, with eigenvalues $0,0,0$. It is nonhyperbolic. 6. 6.

$L_{6}$: $(U,V,Y)=(U_{0},1,0)$. Exists for $n=0$. The eigenvalues are

$0,-\frac{2\left(\sqrt{\texttt{b}}+1\right)(U_{0}-1)}{\texttt{b}-1},-\frac{2\left(\sqrt{\texttt{b}}+1\right)(U_{0}-1)}{\texttt{b}-1}$. It is nonhyperbolic. There is a 2D stable manifold for $0<U_{0}<1,0\leq\texttt{b}<1$, or $\texttt{b}<0,0<U_{0}<1$. There is a 2D unstable manifold for $0<U_{0}<1,\texttt{b}>1$. 7. 7.

$P_{1}$: $(U,V,Y)=\left(0,\frac{n+1}{n+2},0\right)$. The eigenvalues are $-\frac{1}{n+2},\frac{1}{n+2},-\frac{n-3}{n+2}$ . It is a saddle. 8. 8.

$P_{2}$: $(U,V,Y)=(0,1,0)$. The eigenvalues are $-\frac{1}{n+1},-\frac{1}{n+1},-\frac{n}{n+1}$. It is a sink. 9. 9.

$P_{3}$: $(U,V,Y)=(1,1,1)$, with eigenvalues $-\frac{1}{n+1},\frac{1}{n+1},\frac{n}{n+1}$. It is a saddle. 10. 10.

$P_{4}$: $(U,V,Y)=\left(\frac{\texttt{b}(n-3)}{\texttt{b}(n-3)+n-1},\frac{2(n+1)}{3n+1},0\right)$, with eigenvalues

$\frac{2-2n}{(3n+1)(\texttt{b}(n-3)+n-1)},-\frac{(n-1)\left(n+\sqrt{1-n(7n-22)}-5\right)}{2(3n+1)(\texttt{b}(n-3)+n-1)},\frac{(n-1)\left(-n+\sqrt{1-n(7n-22)}+5\right)}{2(3n+1)(\texttt{b}(n-3)+n-1)}$. It is a stable spiral for $\texttt{b}\geq 0,n>5$. 11. 11.

$P_{5}$: $(U,V,Y)=(1,1,0)$. The eigenvalues are $0,0,0$. It is nonhyperbolic. 12. 12.

$P_{6}$: $(U,V,Y)=\left(0,\frac{1}{2\Delta+1},\frac{\Delta}{\Delta+1}\right)$, where $\Delta=\frac{\sqrt{n(5n+4)}-n}{2n}$. The eigenvalues are $0,0,\frac{4\Delta}{(\Delta+1)^{2}(2\Delta+1)}$. It is nonhyperbolic with a 1D unstable manifold. 13. 13.

$P_{7}$: $(U,V,Y)=\left(1,\frac{1}{2\Delta+1},\frac{\Delta}{\Delta+1}\right)$, where $\Delta=\frac{\sqrt{n(5n+4)}-n}{2n}$. The eigenvalues are $0,0,\frac{2\Delta(4\texttt{b}\Delta-3\Delta-1)}{\texttt{b}(\Delta+1)^{2}(2\Delta+1)}$. It is nonhyperbolic. There is a 1D unstable manifold for $\texttt{b}<0,n>0$, or $\frac{3}{4}<\texttt{b}<\frac{1}{8}\left(7+\sqrt{5}\right),0<n<n_{0}$, or $\texttt{b}\geq\frac{1}{8}\left(7+\sqrt{5}\right),n>0$. There is a 1D stable manifold for $0<\texttt{b}\leq\frac{3}{4},n>0$, or $\frac{3}{4}<\texttt{b}<\frac{1}{8}\left(7+\sqrt{5}\right),n>n_{0}$, where $n_{0}=-\frac{(3-4\texttt{b})^{2}}{16\texttt{b}^{2}-28\texttt{b}+11}$. 14. 14.

$P_{8}$: $(U,V,Y)=\left(\frac{1}{2},\frac{2}{(\Delta+1)^{2}},1-\frac{4}{\Delta(\Delta+2)+3}\right)$, with eigenvalues $0,0,0$. It is nonhyperbolic. 15. 15.

$P_{9}$: $(U,V,Y)=\left(1,1,Y_{0}\right)$, with eigenvalues $0,\frac{2Y_{0}((2-4\texttt{b})Y_{0}+1)}{\texttt{b}},-\frac{2Y_{0}((2-4\texttt{b})Y_{0}+1)}{\texttt{b}}$. It is a saddle. 16. 16.

$P_{10}(0)$: $(U,V,Y)=\left(0,V_{0},\frac{1-V_{0}}{1+V_{0}}\right)$, with eigenvalues $0,0,-\frac{8({V_{0}}-1){V_{0}}^{2}}{({V_{0}}+1)^{2}}$. It is nonhyperbolic with a 1D unstable manifold. 17. 17.

$P_{10}(1)$: $(U,V,Y)=\left(1,V_{0},\frac{1-V_{0}}{1+V_{0}}\right)$, with eigenvalues $0,0,\frac{2({V_{0}}-1){V_{0}}(4\texttt{b}({V_{0}}-1)-{V_{0}}+3)}{\texttt{b}({V_{0}}+1)^{2}}$. It is nonhyperbolic. 18. 18.

$P_{11}$: $(U,V,Y)=\left(1,V_{0},\frac{1}{2(2\texttt{b}-1)}\right)$. Exists for $\texttt{b}>\frac{3}{4}$. The eigenvalues are $0,0,0$. It is nonhyperbolic. 19. 19.

$P_{12}$: $(U,V,Y)=\left(U_{0},\frac{1}{2\Delta+1},\frac{\Delta}{\Delta+1}\right)$, where $\Delta=\frac{\sqrt{n(5n+4)}-n}{2n}$, with eigenvalues $0,0,\lambda_{c}$, where $\lambda_{c}$ is a function of $n$ and b. It is nonhyperbolic.

Appendix C Equilibrium points for Model 3

The relation between the new variables (4.1.2) and the set of variables used in [35], say

$e_{1}{}^{1}\frac{d\bar{\lambda}}{dr}=\sqrt{\frac{\mu_{0}}{\texttt{h}}+\theta^{2}},Q=\frac{\theta}{\sqrt{\frac{\mu_{0}}{\texttt{h}}+\theta^{2}}},S=\frac{\sigma}{\sqrt{\frac{\mu_{0}}{\texttt{h}}+\theta^{2}}},C=\frac{\texttt{h}K}{\mu_{0}+\texttt{h}\theta^{2}},A_{\phi}=\frac{\Phi}{\sqrt{2}\sqrt{\frac{\mu_{0}}{\texttt{h}}+\theta^{2}}},A_{W}=\frac{\sqrt{W_{0}}e^{-\frac{k}{2}\phi}}{\sqrt{\frac{\mu_{0}}{\texttt{h}}+\theta^{2}}}$, is

[TABLE]

The details of the stability of the relevant equilibrium point (curves/surfaces of equilibrium point) of the dynamical system (131) are as follows.

$L_{1}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(X_{\phi}^{*},Y_{\phi}^{*},S_{1}^{*},1-\texttt{b}{S_{1}^{*}}^{2}-{X_{\phi}^{*}}^{2},0\right)$. Always exists. The eigenvalues are $0,0,0,-S_{1}^{*},-S_{1}^{*}$. It is nonhyperbolic with a 2D unstable manifold (resp. stable) for $S_{1}^{*}<0$ (resp. $S_{1}^{*}>0$). 2. 2.

$L_{2}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(\varepsilon,\varepsilon,0,0,Y_{0})$, where $\varepsilon=\pm 1$ is the sign of $\theta$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are $0,0,2Y_{0}\varepsilon,4Y_{0}\varepsilon,\left(4-\sqrt{2}k\right)Y_{0}\varepsilon$. $L_{2}^{+}$ (resp. $L_{2}^{-}$ ) is nonhyperbolic with a 3D unstable (resp. stable) manifold for $Y_{0}>0,k<2\sqrt{2}$. Otherwise the unstable (resp. stable) manifold is 2D. 3. 3.

$L_{3}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(-\varepsilon,\varepsilon,0,0,Y_{0})$. Always exists. The eigenvalues are

$0,0,4Y_{0}\varepsilon,2Y_{0}\varepsilon,\left(\sqrt{2}k+4\right)Y_{0}\varepsilon$. $L_{3}^{+}$ (resp. $L_{3}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for $Y_{0}>0,k>-2\sqrt{2}$. Otherwise the unstable (resp. stable) manifold is 2D. 4. 4.

$L_{4}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(0,\varepsilon,0,1,Y_{0})$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are

$0,-Y_{0}\varepsilon,-Y_{0}\varepsilon,2Y_{0}\varepsilon,2Y_{0}\varepsilon$. They are nonhyperbolic and behaves as saddle. 5. 5.

$L_{5}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(\varepsilon X_{\phi}^{*},\varepsilon,\varepsilon\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}},0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}<0,{X_{\phi}^{*}}^{2}\geq 1$ or $\texttt{b}>0,{X_{\phi}^{*}}^{2}\leq 1$. The eigenvalues are

$0,\frac{2}{\texttt{h}}\varepsilon,\frac{-\sqrt{2}k{X_{\phi}^{*}}-\frac{2\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}+4}{\texttt{h}}\varepsilon,\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}\varepsilon,\frac{4\sqrt{\texttt{b}}-(\texttt{h}+2)\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}\texttt{h}}\varepsilon$. $L_{5}^{+}$ (resp. $L_{5}^{-}$) is nonhyperbolic with a 4D unstable (resp. stable) manifold for

1. (a)

$k\in\mathbb{R},{X_{\phi}^{*}}=0,\texttt{h}>1,\texttt{b}>\frac{1}{16}(\texttt{h}+2)^{2}$, or 2. (b)

$-1<{X_{\phi}^{*}}<0,k\geq\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},\texttt{h}>1,\texttt{b}>-\frac{1}{16}(\texttt{h}+2)^{2}({X_{\phi}^{*}}^{2}-1)$, or 3. (c)

$0<{X_{\phi}^{*}}<1,k\leq\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},\texttt{h}>1,\texttt{b}>-\frac{1}{16}(\texttt{h}+2)^{2}({X_{\phi}^{*}}^{2}-1)$, or 4. (d)

$-1<{X_{\phi}^{*}}<0,\frac{2\sqrt{2}}{{X_{\phi}^{*}}}<k<\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},\texttt{h}>-\frac{2\sqrt{2}k{X_{\phi}^{*}}}{\sqrt{2}k{X_{\phi}^{*}}-4},\texttt{b}>-\frac{1}{16}(\texttt{h}+2)^{2}({X_{\phi}^{*}}^{2}-1)$, or 5. (e)

$0<{X_{\phi}^{*}}<1,\frac{2\sqrt{2}}{3{X_{\phi}^{*}}}<k<\frac{2\sqrt{2}}{{X_{\phi}^{*}}},\texttt{h}\geq-\frac{2\sqrt{2}k{X_{\phi}^{*}}}{\sqrt{2}k{X_{\phi}^{*}}-4},\texttt{b}>-\frac{1}{16}(\texttt{h}+2)^{2}({X_{\phi}^{*}}^{2}-1)$, or 6. (f)

$-1<{X_{\phi}^{*}}<0,\frac{2\sqrt{2}}{{X_{\phi}^{*}}}<k<\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},1<\texttt{h}\leq-\frac{2\sqrt{2}k{X_{\phi}^{*}}}{\sqrt{2}k{X_{\phi}^{*}}-4},\texttt{b}>\frac{2({X_{\phi}^{*}}^{2}-1)}{-k^{2}{X_{\phi}^{*}}^{2}+4\sqrt{2}k{X_{\phi}^{*}}-8}$, or 7. (g)

$0<{X_{\phi}^{*}}<1,\frac{2\sqrt{2}}{3{X_{\phi}^{*}}}<k<\frac{2\sqrt{2}}{{X_{\phi}^{*}}},1<\texttt{h}<-\frac{2\sqrt{2}k{X_{\phi}^{*}}}{\sqrt{2}k{X_{\phi}^{*}}-4},\texttt{b}>\frac{2({X_{\phi}^{*}}^{2}-1)}{-k^{2}{X_{\phi}^{*}}^{2}+4\sqrt{2}k{X_{\phi}^{*}}-8}$. 6. 6.

$L_{6}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(\varepsilon X_{\phi}^{*},\varepsilon,-\varepsilon\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}},0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}<0,{X_{\phi}^{*}}^{2}\geq 1$ or $\texttt{b}>0,{X_{\phi}^{*}}^{2}\leq 1$. The eigenvalues are

$0,\frac{2}{\texttt{h}}\varepsilon,\frac{-\sqrt{2}k{X_{\phi}^{*}}+\frac{2\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}+4}{\texttt{h}}\varepsilon,-\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}\varepsilon,\frac{4\sqrt{\texttt{b}}+(\texttt{h}+2)\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}\texttt{h}}\varepsilon$. Nonhyperbolic, behaves as saddle. 7. 7.

$P_{1}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(\varepsilon,\frac{2\sqrt{2}}{k}\varepsilon,0,0,0\right)$. Exists for $k\leq-2\sqrt{2}$ or $k\geq 2\sqrt{2}$. The eigenvalues are $0,0,0,0,0$. They are nonhyperbolic. 8. 8.

$P_{2}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(0,\varepsilon,0,1,\frac{1}{\texttt{h}}\right)$. Always exists. The eigenvalues are

$0,-\frac{1}{\texttt{h}}\varepsilon,-\frac{1}{\texttt{h}}\varepsilon,\frac{2}{\texttt{h}}\varepsilon,\frac{2}{\texttt{h}}\varepsilon$. They are nonhyperbolic. Behaves as a saddle. 9. 9.

$P_{3}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(0,\varepsilon,\frac{1}{\sqrt{\texttt{b}}}\varepsilon,0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}>0$. The eigenvalues are

$0,\frac{4-\frac{2}{\sqrt{\texttt{b}}}}{\texttt{h}}\varepsilon,\frac{2}{\texttt{h}}\varepsilon,\frac{4\sqrt{\texttt{b}}-\texttt{h}-2}{\sqrt{\texttt{b}}\texttt{h}}\varepsilon,\frac{1}{\sqrt{\texttt{b}}}\varepsilon$. $P_{3}^{+}$ (resp. $P_{3}^{-}$) is nonhyperbolic with a 4D unstable (resp. stable) manifold for $\texttt{h}>1,\texttt{b}>\frac{1}{16}(\texttt{h}+2)^{2}$. 10. 10.

$P_{4}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(0,\varepsilon,-\frac{1}{\sqrt{\texttt{b}}}\varepsilon,0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}>0$. The eigenvalues are

$0,\frac{2\left(\frac{1}{\sqrt{\texttt{b}}}+2\right)}{\texttt{h}}\varepsilon,\frac{2}{\texttt{h}}\varepsilon,-\frac{1}{\sqrt{\texttt{b}}}\varepsilon,\frac{4\sqrt{\texttt{b}}+\texttt{h}+2}{\sqrt{\texttt{b}}\texttt{h}}\varepsilon$. They are nonhyperbolic and behaves as saddle. 11. 11.

$P_{5}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(0,\varepsilon,\frac{\texttt{h}+2}{4\texttt{b}}\varepsilon,0,\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{b}<0,\texttt{h}\geq 1$ or $\texttt{h}\geq 1,\texttt{b}\geq\frac{1}{16}(\texttt{h}+2)^{2}$. The eigenvalues are $\frac{\texttt{h}+2}{4\texttt{b}}\varepsilon,\frac{\texttt{h}+2}{4\texttt{b}}\varepsilon,\frac{(\texttt{h}+2)^{2}-8\texttt{b}}{4\texttt{b}\texttt{h}}\varepsilon,\frac{(\texttt{h}+2)^{2}-16\texttt{b}}{8\texttt{b}\texttt{h}}\varepsilon,\frac{(\texttt{h}+2)^{2}-16\texttt{b}}{8\texttt{b}\texttt{h}}\varepsilon$. $P_{5}^{+}$ (resp. $P_{5}^{-}$) is a sink (resp. a source) for $\texttt{b}<0,\texttt{h}>1$. They are saddle otherwise. 12. 12.

$P_{6}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(0,\varepsilon,\frac{2}{\texttt{h}+2}\varepsilon,1-\frac{8\texttt{b}}{(\texttt{h}+2)^{2}},\frac{1}{\texttt{h}}\right)$. Exists for $\texttt{h}\geq 1,0\leq\texttt{b}\leq\frac{1}{8}(\texttt{h}+2)^{2}$. The eigenvalues are

$\frac{2}{\texttt{h}+2}\varepsilon,\frac{2}{\texttt{h}+2}\varepsilon,-\frac{1}{\texttt{h}}\varepsilon,-\frac{\sqrt{64\texttt{b}-7(\texttt{h}+2)^{2}}+\texttt{h}+2}{2\texttt{h}(\texttt{h}+2)}\varepsilon,\frac{\sqrt{64\texttt{b}-7(\texttt{h}+2)^{2}}-\texttt{h}-2}{2\texttt{h}(\texttt{h}+2)}\varepsilon$. They are saddle. 13. 13.

$P_{7}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\left(\frac{k\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{2\sqrt{2}\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,0,\frac{1}{\texttt{h}}\right)$. Exists for

1. (a)

$k\leq-2\sqrt{2},\texttt{b}\geq 0$, or 2. (b)

$-2\sqrt{2}<k<2\sqrt{2},0\leq\texttt{b}\leq-\frac{4}{2k^{2}-16}$, or 3. (c)

$k\geq 2\sqrt{2},\texttt{b}\geq 0$.

The eigenvalues are

$\frac{\sqrt{2}\left(\texttt{b}\left(k^{2}-4\right)+2\right)}{\sqrt{\texttt{b}}\texttt{h}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\texttt{b}\left(k^{2}-8\right)+2}{\sqrt{2}\sqrt{\texttt{b}}\texttt{h}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\texttt{b}\left(k^{2}-8\right)+2}{\sqrt{2}\sqrt{\texttt{b}}\texttt{h}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}\left(\texttt{b}k^{2}-\texttt{h}\right)}{\sqrt{\texttt{b}}\texttt{h}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon$.

$P_{7}^{+}$ (resp. $P_{7}^{-}$) is a source (resp. sink) for

1. (a)

$k\leq-2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},1<\texttt{h}<\texttt{b}k^{2}$, or 2. (b)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or 3. (c)

$2\sqrt{\frac{2}{3}}<k<2\sqrt{2},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<\texttt{b}k^{2}$, or 4. (d)

$k\geq 2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},1<\texttt{h}<\texttt{b}k^{2}$. 14. 14.

$P_{8}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=$

$\left(\frac{k\texttt{b}\varepsilon}{\sqrt{\frac{\texttt{b}k^{2}}{2}+1}\sqrt{\left(k^{2}-2\right)\texttt{b}+2}},\frac{\varepsilon\sqrt{\texttt{b}k^{2}+2}}{\sqrt{\left(k^{2}-2\right)\texttt{b}+2}},\frac{2\varepsilon}{\sqrt{\texttt{b}k^{2}+2}\sqrt{\left(k^{2}-2\right)\texttt{b}+2}},\frac{\left(k^{2}-4\right)\texttt{b}+2}{\left(k^{2}-2\right)\texttt{b}+2},\frac{1}{\texttt{h}}\right)$. Exists for

1. (a)

$k<0,-\frac{2}{k^{2}}\leq\texttt{b}\leq 0$, or 2. (b)

$k=0,\texttt{b}\leq 0$, or 3. (c)

$k>0,-\frac{2}{k^{2}}\leq\texttt{b}\leq 0$.

The eigenvalues are $-\frac{\sqrt{\texttt{b}k^{2}+2}}{\texttt{h}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon,\frac{-\sqrt{\texttt{b}k^{2}+2}-\sqrt{\texttt{b}\left(-\left(7k^{2}-32\right)\right)-14}}{2\texttt{h}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon$,

$\frac{\sqrt{\texttt{b}\left(-\left(7k^{2}-32\right)\right)-14}-\sqrt{\texttt{b}k^{2}+2}}{2\texttt{h}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon,\frac{2\texttt{b}k^{2}-2\texttt{h}}{\texttt{h}\sqrt{\texttt{b}k^{2}+2}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon,\frac{2}{\sqrt{\texttt{b}k^{2}+2}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon$. They are saddle. 15. 15.

$P_{9}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=$

${\left(\varepsilon\frac{k\texttt{h}(2+\texttt{h})}{\sqrt{32k^{4}\texttt{b}^{2}-16\texttt{h}^{3}+k^{2}\texttt{h}\left(16\texttt{b}(-2+\texttt{h})+(2+\texttt{h})^{3}\right)}},\right.}{\varepsilon\frac{2\sqrt{2}\left(2k^{2}\texttt{b}+\texttt{h}^{2}\right)}{\sqrt{32k^{4}\texttt{b}^{2}-16\texttt{h}^{3}+k^{2}\texttt{h}\left(16\texttt{b}(-2+\texttt{h})+(2+\texttt{h})^{3}\right)}},}\\ {\left.\varepsilon\frac{\sqrt{2}k^{2}(2+\texttt{h})}{\sqrt{32k^{4}\texttt{b}^{2}-16\texttt{h}^{3}+k^{2}\texttt{h}\left(16\texttt{b}(-2+\texttt{h})+(2+\texttt{h})^{3}\right)}},0,\frac{1}{\texttt{h}}\right)}$. Exists for

1. (a)

$\texttt{h}\geq 1,k>2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 2. (b)

$\texttt{h}\geq 1,k=-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 3. (c)

$\texttt{h}\geq 1,k=2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 4. (d)

$\texttt{h}\geq 1,k<-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 5. (e)

$\texttt{h}\geq 1,k=0,\texttt{b}<\frac{1}{4}$, or 6. (f)

$\texttt{h}\geq 1,k=-2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 7. (g)

$\texttt{h}\geq 1,k=2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 8. (h)

$\texttt{h}\geq 1,k>2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 9. (i)

$\texttt{h}\geq 1,k<-2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 10. (j)

$\texttt{h}\geq 1,k>2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 11. (k)

$\texttt{h}\geq 1,k<-2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 12. (l)

$\texttt{h}\geq 1,k=-2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 13. (m)

$\texttt{h}\geq 1,k=2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 14. (n)

$\texttt{h}\geq 1,0<k<2\sqrt{2},\texttt{b}<-\frac{2}{k^{2}}$, or 15. (o)

$\texttt{h}\geq 1,-2\sqrt{2}<k<0,\texttt{b}<-\frac{2}{k^{2}}$, or 16. (p)

$\texttt{h}\geq 1,0<k<2\sqrt{2},\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 17. (q)

$\texttt{h}\geq 1,-2\sqrt{2}<k<0,\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 18. (r)

$\texttt{h}\geq 1,0<k<2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 19. (s)

$\texttt{h}\geq 1,-2\sqrt{2}<k<0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$.

The eigenvalues are

${0,0,\varepsilon\frac{2\left(2k\texttt{b}+\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}-4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)\texttt{h}},}\\ {\varepsilon\frac{\left(k+2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}-4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},}\\ {\varepsilon\frac{\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}-4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}\left(k(-2+8\texttt{b}-\texttt{h})-2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\texttt{h}\right)}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)\texttt{h}}}$.

For some choices of parameters, say $\texttt{b}=1,\texttt{h}=1$, or $\texttt{b}=1,\texttt{h}=2$, they are saddle. 16. 16.

$P_{10}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=$${\left(\varepsilon\frac{\sqrt{4+2\texttt{b}\left(-8+k\left(k+4k\texttt{b}-4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{2+k^{2}\texttt{b}},\right.}\\ {\varepsilon\frac{\left(2k\texttt{b}+\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}-4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},}\\ {\left.\varepsilon\frac{\left(k+2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}-4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},0,\frac{1}{\texttt{h}}\right)}.

Exists for

1. (a)

$\texttt{h}\geq 1,k>2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 2. (b)

$\texttt{h}\geq 1,k=-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 3. (c)

$\texttt{h}\geq 1,k=2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 4. (d)

$\texttt{h}\geq 1,k<-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 5. (e)

$\texttt{h}\geq 1,k=0,\texttt{b}<\frac{1}{4}$, or 6. (f)

$\texttt{h}\geq 1,k=-2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 7. (g)

$\texttt{h}\geq 1,k=2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 8. (h)

$\texttt{h}\geq 1,k>2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 9. (i)

$\texttt{h}\geq 1,k<-2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 10. (j)

$\texttt{h}\geq 1,k>2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 11. (k)

$\texttt{h}\geq 1,k<-2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 12. (l)

$\texttt{h}\geq 1,k=-2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 13. (m)

$\texttt{h}\geq 1,k=2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 14. (n)

$\texttt{h}\geq 1,0<k<2\sqrt{2},\texttt{b}<-\frac{2}{k^{2}}$, or 15. (o)

$\texttt{h}\geq 1,-2\sqrt{2}<k<0,\texttt{b}<-\frac{2}{k^{2}}$, or 16. (p)

$\texttt{h}\geq 1,0<k<2\sqrt{2},\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 17. (q)

$\texttt{h}\geq 1,-2\sqrt{2}<k<0,\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 18. (r)

$\texttt{h}\geq 1,0<k<2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 19. (s)

$\texttt{h}\geq 1,-2\sqrt{2}<k<0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$.

The eigenvalues are

${0,0,\varepsilon\frac{2\left(2k\texttt{b}-\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)\texttt{h}},}\\ {\varepsilon\frac{\left(k-2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},}\\ {\varepsilon\frac{\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}\left(k(-2+8\texttt{b}-\texttt{h})+2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\texttt{h}\right)}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)\texttt{h}}}$.

$P_{10}^{+}$ (resp. $P_{10}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for

1. (a)

$-2\sqrt{2}<k<0,\texttt{b}<-\frac{2}{k^{2}},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 2. (b)

$2\sqrt{\frac{2}{3}}<k<2\sqrt{2},\frac{9k^{2}-8}{16k^{2}}<\texttt{b}\leq-\frac{2}{k^{2}-8},1<\texttt{h}<-\frac{2k^{2}}{k^{2}-8}-4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}$, or 3. (c)

$k=2\sqrt{2},\texttt{b}>\frac{1}{2},1<\texttt{h}<4\texttt{b}-1$, or 4. (d)

$k>2\sqrt{2},\texttt{b}>\frac{9k^{2}-8}{16k^{2}},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$.

$P_{10}^{+}$ (resp. $P_{10}^{-}$) nonhyperbolic with a 3D stable (resp. unstable) manifold for

1. (a)

$k=-2\sqrt{\frac{2}{5}},-\frac{5}{4}<\texttt{b}<\frac{1}{4},1<\texttt{h}<\frac{1}{2}\sqrt{5-16\texttt{b}}+\frac{1}{2}$, or 2. (b)

$k=-2\sqrt{\frac{2}{3}},\frac{1}{4}<\texttt{b}<\frac{3}{8},1<\texttt{h}<\sqrt{3-8\texttt{b}}+1$, or 3. (c)

$k=-2\sqrt{\frac{2}{3}},-\frac{3}{4}<\texttt{b}<\frac{1}{4},1<\texttt{h}<\sqrt{3-8\texttt{b}}+1$, or 4. (d)

$-2\sqrt{\frac{2}{3}}<k<-2\sqrt{\frac{2}{5}},\frac{1}{4}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 5. (e)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},\texttt{b}=-\frac{2}{k^{2}-8},1<\texttt{h}<8\texttt{b}-2$, or 6. (f)

$-2\sqrt{\frac{2}{5}}<k<0,-\frac{2}{k^{2}}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 7. (g)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},\frac{1}{4}<\texttt{b}<-\frac{2}{k^{2}-8},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 8. (h)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 9. (i)

$-2\sqrt{\frac{2}{3}}<k<-2\sqrt{\frac{2}{5}},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$. 17. 17.

$P_{11}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=$${\left(\varepsilon\frac{\sqrt{4+2\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{2+k^{2}\texttt{b}},\right.}\\ {\varepsilon\frac{\left(2k\texttt{b}-\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},}\\ {\left.\varepsilon\frac{\left(k-2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},0,\frac{1}{\texttt{h}}\right)}. Exists for

1. (a)

$\texttt{h}=1,k<-2\sqrt{\frac{2}{3}},\texttt{b}=\frac{9k^{2}-8}{16k^{2}}$, or 2. (b)

$\texttt{h}=1,k>2\sqrt{\frac{2}{3}},\texttt{b}=\frac{9k^{2}-8}{16k^{2}}$, or 3. (c)

$\texttt{h}=1,0<k\leq 2\sqrt{\frac{2}{3}},\texttt{b}\leq\frac{9k^{2}-8}{16k^{2}}$, or 4. (d)

$\texttt{h}=1,k<-2\sqrt{\frac{2}{3}},\texttt{b}\leq\frac{1}{k^{2}}$, or 5. (e)

$\texttt{h}=1,k>2\sqrt{\frac{2}{3}},\texttt{b}\leq\frac{1}{k^{2}}$, or 6. (f)

$\texttt{h}=1,-2\sqrt{\frac{2}{3}}\leq k<0,\texttt{b}\leq\frac{9k^{2}-8}{16k^{2}}$, or 7. (g)

$\texttt{h}>1,k>2\sqrt{2}\sqrt{\frac{\texttt{h}}{\texttt{h}+2}},\texttt{b}=\frac{-8\texttt{h}^{2}+\texttt{h}^{2}k^{2}+4\texttt{h}k^{2}+4k^{2}}{16k^{2}}$, or 8. (h)

$\texttt{h}>1,k<-2\sqrt{2}\sqrt{\frac{\texttt{h}}{\texttt{h}+2}},\texttt{b}=\frac{-8\texttt{h}^{2}+\texttt{h}^{2}k^{2}+4\texttt{h}k^{2}+4k^{2}}{16k^{2}}$, or 9. (i)

$\texttt{h}>1,0<k<2\sqrt{2}\sqrt{\frac{\texttt{h}}{\texttt{h}+2}},\texttt{b}\leq\frac{-8\texttt{h}^{2}+\texttt{h}^{2}k^{2}+4\texttt{h}k^{2}+4k^{2}}{16k^{2}}$, or 10. (j)

$\texttt{h}>1,-2\sqrt{2}\sqrt{\frac{\texttt{h}}{\texttt{h}+2}}<k<0,\texttt{b}\leq\frac{-8\texttt{h}^{2}+\texttt{h}^{2}k^{2}+4\texttt{h}k^{2}+4k^{2}}{16k^{2}}$, or 11. (k)

$\texttt{h}>1,k\geq 2\sqrt{2}\sqrt{\frac{\texttt{h}}{\texttt{h}+2}},\texttt{b}\leq\frac{\texttt{h}}{k^{2}}$, or 12. (l)

$\texttt{h}>1,k\leq-2\sqrt{2}\sqrt{\frac{\texttt{h}}{\texttt{h}+2}},\texttt{b}\leq\frac{\texttt{h}}{k^{2}}$.

The eigenvalues are

${0,0,\frac{2\left(2k\texttt{b}-\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)\texttt{h}},}\\ {\frac{\left(k-2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)},}\\ {\frac{\sqrt{2+\texttt{b}\left(-8+k\left(k+4k\texttt{b}+4\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\right)\right)}\left(k(-2+8\texttt{b}-\texttt{h})+2\sqrt{2+\left(-8+k^{2}\right)\texttt{b}}\texttt{h}\right)}{(-1+4\texttt{b})\left(2+k^{2}\texttt{b}\right)\texttt{h}}}$.

$P_{11}^{+}$ (resp. $P_{11}^{-}$) is nonhyperbolic with a 3D stable (resp. unstable) manifold for

1. (a)

$\texttt{b}>\frac{1}{2},k=-2\sqrt{2},1<\texttt{h}<4\texttt{b}-1$, or 2. (b)

$\texttt{b}\geq\frac{9}{16},k<-2\sqrt{2},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 3. (c)

$\texttt{b}<-\frac{1}{4},\sqrt{2}\sqrt{-\frac{1}{\texttt{b}}}<k<2\sqrt{2},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 4. (d)

$\frac{1}{2}<\texttt{b}<\frac{9}{16},-2\sqrt{2}\sqrt{-\frac{1}{16\texttt{b}-9}}<k<-2\sqrt{2},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 5. (e)

$\texttt{b}>\frac{1}{2},-2\sqrt{2}<k\leq-\sqrt{2}\sqrt{\frac{4\texttt{b}-1}{\texttt{b}}},1<\texttt{h}<-\frac{2k^{2}}{k^{2}-8}-4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}$, or 6. (f)

$\frac{3}{8}<\texttt{b}\leq\frac{1}{2},-2\sqrt{2}\sqrt{-\frac{1}{16\texttt{b}-9}}<k\leq-\sqrt{2}\sqrt{\frac{4\texttt{b}-1}{\texttt{b}}},1<\texttt{h}<-\frac{2k^{2}}{k^{2}-8}-4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}$.

$P_{11}^{+}$ (resp. $P_{11}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for

1. (a)

$\texttt{b}>\frac{3}{8},k=\sqrt{2}\sqrt{\frac{4\texttt{b}-1}{\texttt{b}}},1<\texttt{h}<8\texttt{b}-2$, or 2. (b)

$\texttt{b}>\frac{3}{8},\sqrt{2}\sqrt{\frac{4\texttt{b}-1}{\texttt{b}}}<k<2\sqrt{2},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 3. (c)

$-\frac{1}{4}<\texttt{b}<\frac{1}{4},2\sqrt{2}\sqrt{-\frac{1}{16\texttt{b}-9}}<k<2\sqrt{2},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 4. (d)

$\frac{1}{4}<\texttt{b}\leq\frac{3}{8},2\sqrt{2}\sqrt{-\frac{1}{16\texttt{b}-9}}<k<2\sqrt{2},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$, or 5. (e)

$\texttt{b}\leq-\frac{1}{4},2\sqrt{2}\sqrt{-\frac{1}{16\texttt{b}-9}}<k<\sqrt{2}\sqrt{-\frac{1}{\texttt{b}}},1<\texttt{h}<4\sqrt{\frac{\texttt{b}k^{4}-8\texttt{b}k^{2}+2k^{2}}{\left(k^{2}-8\right)^{2}}}-\frac{2k^{2}}{k^{2}-8}$.

Appendix D Equilibrium points for Model 4

The details of the stability of the relevant equilibrium point (curves/surfaces of equilibrium point) of the dynamical system (141) are as follows.

$L_{1}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}X_{\phi}^{*},Y_{\phi}^{*},S_{1}^{*},1-\texttt{b}{S_{1}^{*}}^{2}-{X_{\phi}^{*}}^{2},0\Big{)}. Always exists. The eigenvalues are $0,0,0,-\frac{{S_{1}^{*}}}{n+1},-{S_{1}^{*}}$. Nonhyperbolic with a 2D unstable manifold (resp. stable) for $S_{1}^{*}<0$ (resp. $S_{1}^{*}>0$). 2. 2.

$L_{2}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(\varepsilon,\varepsilon,0,0,Y_{0})$. Exists for $1\leq Y_{0}\leq 1$. The eigenvalues are

$0,0,2Y_{0}\varepsilon,4Y_{0}\varepsilon,\left(4-\sqrt{2}k\right)Y_{0}\varepsilon$. $L_{2}^{+}$ (resp. $L_{2}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for $Y_{0}>0,k<2\sqrt{2}$. Otherwise the unstable (resp. stable) manifold is 2D. 3. 3.

$L_{3}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(-\varepsilon,\varepsilon,0,0,Y_{0})$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are

$0,0,2Y_{0}\varepsilon,4Y_{0}\varepsilon,\left(\sqrt{2}k+4\right)Y_{0}\varepsilon$. $L_{3}^{+}$ (resp. $L_{3}^{-}$) is nonhyperbolic with a 3D unstable (resp. stable) manifold for $Y_{0}>0,k>-2\sqrt{2}$. Otherwise the unstable (resp. stable) manifold is 2D. 4. 4.

$L_{4}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(0,\varepsilon,0,1,Y_{0})$. Exists for $-1\leq Y_{0}\leq 1$. The eigenvalues are

$0,-Y_{0}\varepsilon,-Y_{0}\varepsilon,2Y_{0}\varepsilon,2Y_{0}\varepsilon$. They are nonhyperbolic and behaves as saddle. 5. 5.

$L_{5}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\varepsilon{X_{\phi}^{*}},\varepsilon,\varepsilon\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}},0,1\Big{)}. Exists for $\texttt{b}<0,{X_{\phi}^{*}}^{2}\geq 1$ or $\texttt{b}>0,{X_{\phi}^{*}}^{2}\leq 1$. The eigenvalues are

$0,\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}(n+1)}\varepsilon,2\varepsilon,\left(-\sqrt{2}k{X_{\phi}^{*}}-\frac{2\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}+4\right)\varepsilon,\left(4-\frac{3\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}\right)\varepsilon$. $L_{5}^{+}$ (resp. $L_{5}^{-}$) is nonhyperbolic, with a 4D unstable (resp. stable) manifold for

1. (a)

$k\in\mathbb{R},{X_{\phi}^{*}}=0,\texttt{b}>\frac{9}{16},n>0$, or 2. (b)

$-1<{X_{\phi}^{*}}<0,k>\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},\texttt{b}>\frac{9}{16}\left(1-{X_{\phi}^{*}}^{2}\right),n>0$, or 3. (c)

$0<{X_{\phi}^{*}}<1,k\leq\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},\texttt{b}>\frac{9}{16}\left(1-{X_{\phi}^{*}}^{2}\right),n>0$, or 4. (d)

$-1<{X_{\phi}^{*}}<0,\frac{2\sqrt{2}}{{X_{\phi}^{*}}}<k\leq\frac{2\sqrt{2}}{3{X_{\phi}^{*}}},\texttt{b}>\frac{2-2{X_{\phi}^{*}}^{2}}{k^{2}{X_{\phi}^{*}}^{2}-4\sqrt{2}k{X_{\phi}^{*}}+8},n>0$, or 5. (e)

$0<{X_{\phi}^{*}}<1,\frac{2\sqrt{2}}{3{X_{\phi}^{*}}}<k<\frac{2\sqrt{2}}{{X_{\phi}^{*}}},\texttt{b}>\frac{2-2{X_{\phi}^{*}}^{2}}{k^{2}{X_{\phi}^{*}}^{2}-4\sqrt{2}k{X_{\phi}^{*}}+8},n>0$. 6. 6.

$L_{6}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\varepsilon{X_{\phi}^{*}},\varepsilon,-\varepsilon\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}},0,1\Big{)}. Exists for $\texttt{b}<0,{X_{\phi}^{*}}^{2}\geq 1$ or $\texttt{b}>0,{X_{\phi}^{*}}^{2}\leq 1$. The eigenvalues are

$0,-\frac{\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}(n+1)}\varepsilon,2\varepsilon,\left(-\sqrt{2}k{X_{\phi}^{*}}+\frac{2\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}+4\right)\varepsilon,\left(\frac{3\sqrt{1-{X_{\phi}^{*}}^{2}}}{\sqrt{\texttt{b}}}+4\right)\varepsilon$. They are nonhyperbolic and behaves as saddle. 7. 7.

$P_{1}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\varepsilon,\frac{2\sqrt{2}}{k}\varepsilon,0,0,0\Big{)}. Exists for $k\leq-2\sqrt{2}$, or $k\geq 2\sqrt{2}$. The eigenvalues are $0,0,0,0,0$. They are nonhyperbolic. 8. 8.

$P_{2}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}0,\varepsilon,\frac{2}{3}\varepsilon,1-\frac{8\texttt{b}}{9},1\Big{)}. Exists for $0\leq\texttt{b}\leq\frac{9}{8}$. The eigenvalues are

$\frac{2}{3(n+1)}\varepsilon,-\varepsilon,\frac{2}{3}\varepsilon,\frac{1}{6}\left(-\sqrt{64\texttt{b}-63}-3\right)\varepsilon,\frac{1}{6}\left(\sqrt{64\texttt{b}-63}-3\right)\varepsilon$. They are saddle. 9. 9.

$P_{3}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=(0,\varepsilon,0,1,1)$. The eigenvalues are $0,-\varepsilon,-\varepsilon,2\varepsilon,2\varepsilon$. They are nonhyperbolic and behaves as saddle. 10. 10.

$P_{4}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}0,\varepsilon,\frac{\varepsilon}{\sqrt{\texttt{b}}},0,1\Big{)}. Exists for $\texttt{b}>0$. The eigenvalues are

$0,\frac{1}{\sqrt{\texttt{b}}(n+1)}\varepsilon,\frac{2\left(2\sqrt{\texttt{b}}-1\right)}{\sqrt{\texttt{b}}}\varepsilon,\frac{4\sqrt{\texttt{b}}-3}{\sqrt{\texttt{b}}}\varepsilon,2\varepsilon$. $P_{4}^{+}$ (resp. $P_{4}^{-}$) is nonhyperbolic with a 4D unstable (resp. stable) manifold for $n>0,\texttt{b}>\frac{9}{16}$. 11. 11.

$P_{5}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}0,\varepsilon,\frac{3\varepsilon}{4\texttt{b}},0,1\Big{)}. Exists for $\texttt{b}<0$, or $\texttt{b}\geq\frac{9}{16}$. The eigenvalues are

$\frac{3}{4\texttt{b}(n+1)}\varepsilon,\frac{3}{4\texttt{b}}\varepsilon,-\frac{8\texttt{b}-9}{4\texttt{b}}\varepsilon,-\frac{16\texttt{b}-9}{8\texttt{b}}\varepsilon,-\frac{16\texttt{b}-9}{8\texttt{b}}\varepsilon$. $P_{5}^{+}$ (resp. $P_{5}^{-}$) is a sink (resp. a source) for $\texttt{b}<0,n>0$ or a saddle otherwise. 12. 12.

$P_{6}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}0,\varepsilon,-\frac{\varepsilon}{\sqrt{\texttt{b}}},0,1\Big{)}. Exists for $\texttt{b}>0$. The eigenvalues are

$0,-\frac{1}{\sqrt{\texttt{b}}(n+1)}\varepsilon,\frac{2\left(2\sqrt{\texttt{b}}+1\right)}{\sqrt{\texttt{b}}}\varepsilon,\frac{4\sqrt{\texttt{b}}+3}{\sqrt{\texttt{b}}}\varepsilon,2\varepsilon$. They are nonhyperbolic and behaves as saddle. 13. 13.

$P_{7}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\frac{k\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{2\sqrt{2}\sqrt{\texttt{b}}}{\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,0,1\Big{)}. Exists for

1. (a)

$k\leq-2\sqrt{2},\texttt{b}\geq 0$, or 2. (b)

$-2\sqrt{2}<k<2\sqrt{2},0\leq\texttt{b}\leq-\frac{4}{2k^{2}-16}$, or 3. (c)

$k\geq 2\sqrt{2},\texttt{b}\geq 0$.

The eigenvalues are

$\frac{\sqrt{2}\left(\texttt{b}\left(k^{2}-4\right)+2\right)}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}}{\sqrt{\texttt{b}}(n+1)\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\sqrt{2}\left(\texttt{b}k^{2}-1\right)}{\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\texttt{b}\left(k^{2}-8\right)+2}{\sqrt{2}\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon,\frac{\texttt{b}\left(k^{2}-8\right)+2}{\sqrt{2}\sqrt{\texttt{b}}\sqrt{\texttt{b}k^{2}+2}}\varepsilon$.

$P_{7}^{+}$ (resp. $P_{7}^{-}$) is a source (resp. a sink) for

1. (a)

$k\leq-2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},n>0$, or 2. (b)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},n>0$, or 3. (c)

$2\sqrt{\frac{2}{3}}<k<2\sqrt{2},\frac{1}{k^{2}}<\texttt{b}<-\frac{2}{k^{2}-8},n>0$, or 4. (d)

$k\geq 2\sqrt{2},\texttt{b}>\frac{1}{k^{2}},n>0$.

It is a saddle otherwise. 14. 14.

$P_{8}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=$

$\left(\frac{k\texttt{b}\varepsilon}{\sqrt{\left(\frac{\texttt{b}k^{2}}{2}+1\right)\left(\left(k^{2}-2\right)\texttt{b}+2\right)}},\varepsilon\sqrt{\frac{{\texttt{b}k^{2}+2}}{{\left(k^{2}-2\right)\texttt{b}+2}}},\frac{2\varepsilon}{\sqrt{\left(\texttt{b}k^{2}+2\right)\left(\left(k^{2}-2\right)\texttt{b}+2\right)}},\frac{\left(k^{2}-4\right)\texttt{b}+2}{\left(k^{2}-2\right)\texttt{b}+2},1\right)$. Exists for

1. (a)

$k<0,-\frac{2}{k^{2}}\leq\texttt{b}\leq 0$, or 2. (b)

$k=0,\texttt{b}\leq 0$, or 3. (c)

$k>0,-\frac{2}{k^{2}}\leq\texttt{b}\leq 0$.

The eigenvalues are

$\frac{2\texttt{b}k^{2}-2}{\sqrt{\texttt{b}k^{2}+2}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon,-\frac{\sqrt{\texttt{b}k^{2}+2}}{\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon,\frac{2}{(n+1)\sqrt{\texttt{b}k^{2}+2}\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon$,

$\frac{-\sqrt{\texttt{b}k^{2}+2}-\sqrt{\texttt{b}\left(-\left(7k^{2}-32\right)\right)-14}}{2\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon,\frac{\sqrt{\texttt{b}\left(-\left(7k^{2}-32\right)\right)-14}-\sqrt{\texttt{b}k^{2}+2}}{2\sqrt{\texttt{b}\left(k^{2}-2\right)+2}}\varepsilon$. They are saddle. 15. 15.

$P_{9}^{\varepsilon}$: $(X_{\phi},Y_{\phi},S_{1},S_{2},Y)=$

\Big{(}\frac{3k\varepsilon}{\sqrt{k^{2}\left(16\texttt{b}\left(2k^{2}\texttt{b}-1\right)+27\right)-16}},\frac{2\left(2\texttt{b}k^{2}+1\right)\varepsilon}{\sqrt{16\texttt{b}^{2}k^{4}+\left(\frac{27}{2}-8\texttt{b}\right)k^{2}-8}},\frac{3k^{2}\varepsilon}{\sqrt{16\texttt{b}^{2}k^{4}+\left(\frac{27}{2}-8\texttt{b}\right)k^{2}-8}},0,1\Big{)}. Exists for

1. (a)

$k<-2\sqrt{\frac{2}{3}},\texttt{b}=\frac{9k^{2}-8}{16k^{2}}$, or 2. (b)

$k>2\sqrt{\frac{2}{3}},\texttt{b}=\frac{9k^{2}-8}{16k^{2}}$, or 3. (c)

$k<-2\sqrt{\frac{2}{3}},\texttt{b}\leq\frac{1}{k^{2}}$, or 4. (d)

$k>2\sqrt{\frac{2}{3}},\texttt{b}\leq\frac{1}{k^{2}}$, or 5. (e)

$0<k\leq 2\sqrt{\frac{2}{3}},\texttt{b}\leq\frac{9k^{2}-8}{16k^{2}}$, or 6. (f)

$-2\sqrt{\frac{2}{3}}\leq k<0,\texttt{b}\leq\frac{9k^{2}-8}{16k^{2}}$.

For $\texttt{b}=1$ the eigenvalues are

$\frac{3k^{2}\varepsilon}{(n+1)\sqrt{16k^{4}+\frac{11k^{2}}{2}-8}},\frac{-7k^{2}-8}{\sqrt{64k^{4}+22k^{2}-32}}\varepsilon,\frac{-7k^{2}-i\sqrt{7k^{2}+8}\sqrt{17k^{2}-32}-8}{2\sqrt{64k^{4}+22k^{2}-32}}\varepsilon$,

$\frac{-7k^{2}+i\sqrt{7k^{2}+8}\sqrt{17k^{2}-32}-8}{2\sqrt{64k^{4}+22k^{2}-32}}\varepsilon,\frac{k^{2}-4}{\sqrt{16k^{4}+\frac{11k^{2}}{2}-8}}\varepsilon$. For $\texttt{b}=1$ they are saddle. 16. 16.

$P_{10}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\varepsilon\frac{\sqrt{2\texttt{b}\left(k\left(4\texttt{b}k+k-4\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)-8\right)+4}}{\texttt{b}k^{2}+2},

$\varepsilon\frac{\left(2k\texttt{b}+\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)\sqrt{\texttt{b}\left(k\left(4\texttt{b}k+k-4\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)}$,

\varepsilon\frac{\left(k+2\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)\sqrt{\texttt{b}\left(k\left(4\texttt{b}k+k-4\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)},0,1\Big{)}. Exists for

1. (a)

$k=0,\texttt{b}<\frac{1}{4}$, or 2. (b)

$k=2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 3. (c)

$k=2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 4. (d)

$k=2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 5. (e)

$k=-2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 6. (f)

$k=-2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 7. (g)

$k=-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 8. (h)

$0<k<2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 9. (i)

$0<k<2\sqrt{2},\texttt{b}<-\frac{2}{k^{2}}$, or 10. (j)

$0<k<2\sqrt{2},\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 11. (k)

$-2\sqrt{2}<k<0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 12. (l)

$-2\sqrt{2}<k<0,\texttt{b}<-\frac{2}{k^{2}}$, or 13. (m)

$-2\sqrt{2}<k<0,\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 14. (n)

$k<-2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 15. (o)

$k<-2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 16. (p)

$k<-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 17. (q)

$k>2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 18. (r)

$k>2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 19. (s)

$k>2\sqrt{2},\texttt{b}>\frac{1}{4}$.

The eigenvalues are

$0,0,\frac{2\left(\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+2\texttt{b}k\right)\sqrt{\texttt{b}\left(k\left(-4\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+4\texttt{b}k+k\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)}\varepsilon$,

$\frac{\left(2\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+k\right)\sqrt{\texttt{b}\left(k\left(-4\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+4\texttt{b}k+k\right)-8\right)+2}}{(4\texttt{b}-1)(n+1)\left(\texttt{b}k^{2}+2\right)}\varepsilon$,

$\frac{\left((8\texttt{b}-3)k-2\sqrt{\texttt{b}\left(k^{2}-8\right)+2}\right)\sqrt{\texttt{b}\left(k\left(-4\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+4\texttt{b}k+k\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)}\varepsilon$.

$P_{10}^{+}$ (resp. $P_{10}^{-}$) has a 3D stable (resp. unstable) manifold for

1. (a)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},n\geq 0,\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 2. (b)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 3. (c)

$k=-2\sqrt{\frac{2}{3}},n\geq 0,-\frac{3}{4}<\texttt{b}<\frac{1}{4}$, or 4. (d)

$k=-2\sqrt{\frac{2}{3}},n\geq 0,\frac{1}{4}<\texttt{b}<\frac{3}{8}$, or 5. (e)

$-2\sqrt{\frac{2}{3}}<k<-2\sqrt{\frac{2}{5}},n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 6. (f)

$-2\sqrt{\frac{2}{3}}<k<-2\sqrt{\frac{2}{5}},n\geq 0,\frac{1}{4}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}}$, or 7. (g)

$k=-2\sqrt{\frac{2}{5}},n\geq 0,-\frac{5}{4}<\texttt{b}<\frac{1}{4}$, or 8. (h)

$-2\sqrt{\frac{2}{5}}<k<-\frac{2\sqrt{2}}{3},n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}}$, or 9. (i)

$k=-\frac{2\sqrt{2}}{3},n\geq 0,-\frac{9}{4}<\texttt{b}<0$, or 10. (j)

$-\frac{2\sqrt{2}}{3}<k<0,n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}}$.

$P_{10}^{+}$ (resp. $P_{10}^{-}$) has a 3D unstable (resp. stable) manifold for

1. (a)

$-2\sqrt{2}<k<0,n\geq 0,\texttt{b}<-\frac{2}{k^{2}}$, or 2. (b)

$2\sqrt{\frac{2}{3}}<k<2\sqrt{2},n\geq 0,\frac{9k^{2}-8}{16k^{2}}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 3. (c)

$k\geq 2\sqrt{2},n\geq 0,\texttt{b}>\frac{9k^{2}-8}{16k^{2}}$. 17. 17.

$P_{11}^{\varepsilon}$: (X_{\phi},Y_{\phi},S_{1},S_{2},Y)=\Big{(}\varepsilon\frac{\sqrt{2\texttt{b}\left(k\left(4\texttt{b}k+k+4\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)-8\right)+4}}{\texttt{b}k^{2}+2},

$\varepsilon\frac{\left(2k\texttt{b}-\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)\sqrt{\texttt{b}\left(k\left(4\texttt{b}k+k+4\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)}$,

\varepsilon\frac{\left(k-2\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)\sqrt{\texttt{b}\left(k\left(4\texttt{b}k+k+4\sqrt{\left(k^{2}-8\right)\texttt{b}+2}\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)},0,1\Big{)}. Exists for

1. (a)

$k=0,\texttt{b}<\frac{1}{4}$, or 2. (b)

$k=2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 3. (c)

$k=2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 4. (d)

$k=2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 5. (e)

$k=-2\sqrt{2},-\frac{1}{4}<\texttt{b}<\frac{1}{4}$, or 6. (f)

$k=-2\sqrt{2},\texttt{b}<-\frac{1}{4}$, or 7. (g)

$k=-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 8. (h)

$0<k<2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 9. (i)

$0<k<2\sqrt{2},\texttt{b}<-\frac{2}{k^{2}}$, or 10. (j)

$0<k<2\sqrt{2},\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 11. (k)

$-2\sqrt{2}<k<0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 12. (l)

$-2\sqrt{2}<k<0,\texttt{b}<-\frac{2}{k^{2}}$, or 13. (m)

$-2\sqrt{2}<k<0,\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 14. (n)

$k<-2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 15. (o)

$k<-2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 16. (p)

$k<-2\sqrt{2},\texttt{b}>\frac{1}{4}$, or 17. (q)

$k>2\sqrt{2},-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 18. (r)

$k>2\sqrt{2},-\frac{2}{k^{2}-8}\leq\texttt{b}<-\frac{2}{k^{2}}$, or 19. (s)

$k>2\sqrt{2},\texttt{b}>\frac{1}{4}$.

The eigenvalues are

$0,0,\frac{\left(k-2\sqrt{\texttt{b}\left(k^{2}-8\right)+2}\right)\sqrt{\texttt{b}\left(k\left(4\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+4\texttt{b}k+k\right)-8\right)+2}}{(4\texttt{b}-1)(n+1)\left(\texttt{b}k^{2}+2\right)}\varepsilon$,

$\frac{2\left(2\texttt{b}k-\sqrt{\texttt{b}\left(k^{2}-8\right)+2}\right)\sqrt{\texttt{b}\left(k\left(4\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+4\texttt{b}k+k\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)}\varepsilon$,

$\frac{\left(2\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+(8\texttt{b}-3)k\right)\sqrt{\texttt{b}\left(k\left(4\sqrt{\texttt{b}\left(k^{2}-8\right)+2}+4\texttt{b}k+k\right)-8\right)+2}}{(4\texttt{b}-1)\left(\texttt{b}k^{2}+2\right)}\varepsilon$.

$P_{11}^{+}$ (resp. $P_{11}^{-}$) has a 3D stable (resp. unstable) manifold for

1. (a)

$k\leq-2\sqrt{2},n\geq 0,\texttt{b}>\frac{9k^{2}-8}{16k^{2}}$, or 2. (b)

$-2\sqrt{2}<k<-2\sqrt{\frac{2}{3}},n\geq 0,\frac{9k^{2}-8}{16k^{2}}<\texttt{b}\leq-\frac{2}{k^{2}-8}$, or 3. (c)

$0<k<2\sqrt{2},n\geq 0,\texttt{b}<-\frac{2}{k^{2}}$.

$P_{11}^{+}$ (resp. $P_{11}^{-}$) has a 3D unstable (resp. stable) manifold for

1. (a)

$k=\frac{2\sqrt{2}}{3},n\geq 0,-\frac{9}{4}<\texttt{b}<0$, or 2. (b)

$2\sqrt{\frac{2}{5}}<k<2\sqrt{2},n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{1}{4}$, or 3. (c)

$0<k<\frac{2\sqrt{2}}{3},n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}}$, or 4. (d)

$\frac{2\sqrt{2}}{3}<k\leq 2\sqrt{\frac{2}{5}},n\geq 0,-\frac{2}{k^{2}}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}}$, or 5. (e)

$2\sqrt{\frac{2}{5}}<k\leq 2\sqrt{\frac{2}{3}},n\geq 0,\frac{1}{4}<\texttt{b}<\frac{9k^{2}-8}{16k^{2}}$, or 6. (f)

$2\sqrt{\frac{2}{3}}<k<2\sqrt{2},n\geq 0,\frac{1}{4}<\texttt{b}\leq-\frac{2}{k^{2}-8}$.