\sc Scalar field deformations of Lambda-CDM cosmology
{\sc Artur Alho
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This paper treats nonrelativistic matter and a scalar field $\phi$ with a monotonically decreasing potential minimally coupled to gravity in flat Friedmann-Lema\^{itre-Robertson-Walker cosmology. The field equations are reformulated as a three-dimensional dynamical system on an extended compact state space, complemented with cosmographic diagrams. A dynamical systems analysis provides global dynamical results describing possible asymptotic behavior. It is shown that one should impose global and asymptotic bounds on $\lambda=-V^{-1\,dV/d\phi$ to obtain viable cosmological models that continuously deform $\Lambda$CDM cosmology. In particular we introduce a regularized inverse power-law potential as a simple specific example. %Revisiting a select set of familiar potentials, but in a generic global %dynamical systems setting, suggests that one should impose global and %asymptotic bounds on $\lambda=-V^{-1\,dV/d\phi$ %, and not just demanding a shallow tail of a %scalar field potential $V(\phi)$, as is the case for e.g. an inverse %power-law potential, %in order to obtain an attracting separatrix surface that describes the %evolution of all solutions with $\Omega_m$ initially close to one, set %by, e.g., inflation. %An attracting separatrix surface is to be contrasted with %individually attracting solutions, such as tracker solutions that occur for %potentials with unbounded $\lambda$. %Furthermore, letting $\lambda \rightarrow 0$ for all $\phi$ and %$\Omega_m\rightarrow 1$ initially yields $\Lambda$CDM cosmology and thus it %is possible to continuously modify $\Lambda$CDM dynamics by choosing %potentials that continuously deform the $\Lambda$CDM attracting separatrix %urface, thereby also providing a useful testbed for $\Lambda$CDM cosmology; %specific examples of potentials and the associated global dynamics are given. %To facilitate %comparisons with $\Lambda$CDM cosmology and to complement the dynamical %systems pictures, we provide cosmographic diagrams. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DeformOct21_2015.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  \filename{Deform_Oct08_2015_resubmit.tex},
%  Updated cu:  2015-10-08 aa: 2015-04-28
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[12pt]{article}

\usepackage{amssymb,amsfonts,amsthm}%,psfig}
%\usepackage[draft]{graphicx}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{bm}             %% boldmath-command: \bm{}
\usepackage[mathscr]{eucal}

\usepackage{cancel}
%\usepackage{outbraces}

\usepackage{psfrag}

\usepackage{multirow}

%\usepackage{showlabels}

\usepackage{wasysym}

\usepackage{color}

%%\documentclass{article}
%\usepackage{amssymb,amsfonts,amsthm,graphicx,psfig}
%\usepackage{amssymb,amsfonts,amsthm}
%\usepackage[draft]{graphicx}
%\usepackage{graphicx}
%\usepackage{subfigure}
%%\usepackage{amssymb}
%%\usepackage{amsthm}
%%\usepackage{amsmath}
%%\usepackage{bm}             %% boldmath-command: \bm{}
%\usepackage[mathscr]{eucal}
%\usepackage{cancel}
%\usepackage{oubraces}
%\usepackage{psfrag}
%\usepackage{wasysym}


%\usepackage{fancyheadings}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\topskip 0mm \topmargin 0mm \headheight 0mm \textwidth 15cm
%%\textwidth 16cm
%\textheight 23.5cm \evensidemargin 0mm \oddsidemargin 0mm
\pagestyle{headings}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------- DIN A4: 21cm times 29.7cm ------------%

\topmargin-1.9cm %+1in %% top edge to top of (a box containing) the head
%\topmargin-0.9cm %+1in %% top edge to top of (a box containing) the head
%\headheight            %% height of (a box containing) the head
\headsep1.2cm             %% vertical space between head-box and body-box
%\topskip               %% space between top of body-box and bottom of first line

\textheight22.5cm
\textwidth15cm

\marginparwidth2.3cm    %% width of marginal notes
\marginparsep0.5cm      %% distance between margin of text and marginal note

\parskip0.2cm     %% space between paragraphs
\parindent0cm     %% indentation
%\baselineskip     %% vertical distance between lines

\footskip1.5cm     %% distance bottom last line to bottom of foot
%\footheight      %% height of foot-box

\oddsidemargin0.5cm %+1in  %% distance left edge to text on rh-pages
\evensidemargin-0.6cm %+1in  %% on lh-pages

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\parindent0cm     %% indentation
%\parskip2ex
%\topskip 0mm \topmargin 0mm \headheight 0mm
%\textwidth 15cm
%\textwidth 16cm
%\textheight 23.5cm \evensidemargin 0mm \oddsidemargin 0mm
%\pagestyle{headings}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\mnote}[1]{\marginpar{\raggedright\tiny\em
$\!\!\!\!\!\!\,\bullet$ #1}}       %Marginal Note

\newtheorem*{conjecture}{Conjecture}

\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\def\hsp5{\hspace{5mm}}
\newcommand{\marginalnote}[1]{\marginpar{\begin{footnotesize}{\it #1}\end{footnotesize}}}
\newcommand{\textfrac}[2]{{\textstyle \frac{#1}{#2}}}
\newcommand{\sfrac}[2]{{\textstyle{#1\over#2}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ----------------------------------------------------------------
\title{\sc Scalar field deformations of Lambda-CDM cosmology}
% ----------------------------------------------------------------

\begin{document}

\author{{\sc Artur Alho}\thanks{Electronic address:
{\tt aalho@math.ist.utl.pt}} \\[1ex]
Center for Mathematical Analysis, Geometry and Dynamical Systems,\\
Instituto Superior T\'ecnico, Universidade de Lisboa, \\
Av. Rovisco Pais, 1049-001 Lisboa, Portugal
\and \\
{\sc Claes Uggla}\thanks{Electronic address:
{\tt claes.uggla@kau.se}} \\[1ex]
Department of Physics, \\
Karlstad University, S-651 88 Karlstad, Sweden \\[2ex] }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\date{\normalsize{October 04, 2015}}
\date{}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This paper treats nonrelativistic matter and a scalar field $\phi$ with a
monotonically decreasing potential minimally coupled to gravity in flat
Friedmann-Lema\^{i}tre-Robertson-Walker cosmology. The field equations are
reformulated as a three-dimensional dynamical system on an extended compact
state space, complemented with cosmographic diagrams. A dynamical systems
analysis provides global dynamical results describing possible asymptotic
behavior. It is shown that one should impose \emph{global and asymptotic}
bounds on $\lambda=-V^{-1}\,dV/d\phi$ to obtain viable cosmological models
that continuously deform $\Lambda$CDM cosmology. In particular we introduce a
regularized inverse power-law potential as a simple specific example.
%Revisiting a select set of familiar potentials, but in a generic global
%dynamical systems setting, suggests that one should impose \emph{global and
%asymptotic} bounds on $\lambda=-V^{-1}\,dV/d\phi$
%, and not just demanding a shallow tail of a
%scalar field potential $V(\phi)$, as is the case for e.g. an inverse
%power-law potential,
%in order to obtain an \emph{attracting separatrix surface} that describes the
%evolution of \emph{all} solutions with $\Omega_m$ initially close to one, set
%by, e.g., inflation.
%An attracting separatrix \emph{surface} is to be contrasted with
%individually attracting solutions, such as tracker solutions that occur for
%potentials with unbounded $\lambda$.
%Furthermore, letting $\lambda \rightarrow 0$ for all $\phi$ and
%$\Omega_m\rightarrow 1$ initially yields $\Lambda$CDM cosmology and thus it
%is possible to continuously modify $\Lambda$CDM dynamics by choosing
%potentials that continuously deform the $\Lambda$CDM attracting separatrix
%urface, thereby also providing a useful testbed for $\Lambda$CDM cosmology;
%specific examples of potentials and the associated global dynamics are given.
%To facilitate
%comparisons with $\Lambda$CDM cosmology and to complement the dynamical
%systems pictures, we provide cosmographic diagrams.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\centerline{\bigskip\noindent PACS numbers: 04.20.-q, 98.80.-k, 98.80.Bp,
%98.80.Jk}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Recent observations, such as those of the cosmic background radiation, show
that $\Lambda$ cold dark matter ($\Lambda$CDM) cosmology is remarkably
consistent with observations, although some tensions remain,
see~\cite{plaI15,plaXIII15} and references therein. On the other hand, the
unknown origins and nature of dark matter and dark energy remain mysteries
which generate a steady flow of models and theories. Nonetheless, current
observations should be taken seriously, and viable cosmological models should
therefore presumably only deviate from $\Lambda$CDM cosmology marginally.
Moreover, although successful, the $\Lambda$CDM model still needs to be
observationally tested, which suggests that it should be continuously
deformed.
%or, preferably, a theory should admit a continuous set of models that each
%admit an open set of solutions that represents a continuous deformation of
%the $\Lambda$CDM models, thereby providing a natural testing ground for these
%models.
In addition, for a variety of reasons, it is of interest to consider a
dynamical dark energy content described by a field theoretical model instead
of a cosmological constant $\Lambda$ or a dark energy fluid. This leads to
the question: What phenomenological conditions must be imposed on an
effective classical field description of dark energy to continuously deform
$\Lambda$CDM cosmology for an open set of observationally viable solutions?
%What phenomenological conditions must be imposed on an effective classical
%field description of dark energy in order to obtain a continuous set of
%models where each model results in an open set of solutions that continuously
%deform $\Lambda$CDM cosmology?

Although observational compatibility arguably comes first, there also exist
other issues that are of interest such as various fine-tuning problems and if
there exist observationally viable models that emerge from, or at least have
ties to, a fundamental theory. The rather modest purpose of this paper,
however, is to
%
\begin{itemize}
\item[(i)] describe general features and determine asymptotic behavior of
    flat Friedmann-Lema\^{i}tre-Robertson-Walker (FLRW) models with a
    perfect fluid and a scalar field $\phi$, minimally coupled to
    gravity, with a positive monotonic potential $V(\phi)$,
\item[(ii)] shed some light on the above ``observationally viable field
    theoretical $\Lambda$CDM-deformation issue'' by considering some
    simple models of the above type with the equation of state for the
    perfect fluid specialized to dust, representing nonrelativistic
    matter (mainly baryons and nonbaryonic dark matter).
\end{itemize}
%
Since the literature about minimally coupled scalar fields and a perfect
fluid in a flat FLRW spacetime geometry is vast, we will start with only a
few general sample references: the influential
papers~\cite{peerat88,ratpee88,steetal99} and the recent review~\cite{tsu13},
and references therein.

It is of interest for discussions and future purposes to first consider a
perfect fluid with a general barotropic equation of state $p_m =
(\gamma_m-1)\rho_m$, where $p_m$ and $\rho_m \geq 0$ are the pressure and
energy density, respectively, and $\gamma_m=\gamma_m(\rho_m)$, together with
a scalar field $\phi$, minimally coupled to gravity, with a self-interaction
potential $V(\phi)$. Thus, the Einstein and matter field
equations for the present flat FLRW models can be written as (see,
e.g.,~\cite{waiell97,col03,wei08})
%
\begin{subequations}\label{Hphieq}
\begin{align}
3H^2 &= \frac12 \dot{\phi}^2 + V(\phi) + \rho_m,\label{Gauss1}\\
\dot{H} &= -\frac{1}{2}\left(\dot{\phi}^2 + \gamma_m\rho_m\right),\label{Ray1}\\
0 &= \ddot{\phi} + 3H\dot{\phi} + V_{\phi},\label{KG}\\
\dot{\rho}_m &= - 3H\gamma_m\rho_m ,\label{rhoeq}
\end{align}
\end{subequations}
%
where $V_\phi = dV/d\phi$ and where overdots denote derivatives with respect
to cosmic time $t$. Units are such that $8\pi G = 1 = c$, where $G$ is the
gravitational constant and $c$ the speed of light in vacuum. It is assumed
throughout the paper that the Universe is expanding, i.e., $H>0$.

Equation~\eqref{KG} can be heuristically interpreted as an equation for a
particle of unit mass with a one-dimensional coordinate $\phi$, moving in a
potential $V(\phi)$ with a friction force $-3H\dot{\phi}$.
Equation~\eqref{Gauss1} shows that $H$ can be expressed in terms of $\phi$,
$\dot{\phi}$, and $\rho_m$. Introducing $\dot{\phi}$ as an independent
variable therefore leads to a three-dimensional dynamical system for
$\phi,\dot{\phi},\rho_m$. Alternatively one can solve Eq.~\eqref{rhoeq}
to obtain $\rho_m=\rho_m(a)$ and consider the variables $\phi$, $\dot{\phi},
a$ as the dependent variables, where $a$ obeys the equation $\dot{a} = aH$,
where $H$ can be expressed in terms of $\phi,\dot{\phi},a$ by means
of~\eqref{Gauss1}.

It is of some interest to introduce effective equation of state parameters,
defined by
%
\begin{equation}
\gamma_\phi = \frac{\rho_\phi + p_\phi}{\rho_\phi} =
\frac{\dot{\phi}^2}{\frac12\dot{\phi}^2 + V(\phi)},\qquad
\gamma_{\mathrm{tot}} = \frac{\rho_{\mathrm{tot}} + p_{\mathrm{tot}}}{\rho_{\mathrm{tot}}},
\end{equation}
%
or, alternatively, $w_* = p_*/\rho_* = \gamma_* - 1$ ($*=\phi,\,
\mathrm{tot}$), where
%
\begin{subequations}
\begin{alignat}{2}
\rho_\phi &= \frac12\dot{\phi}^2 + V(\phi),&\qquad p_\phi &=
\frac12\dot{\phi}^2 - V(\phi),\\
\rho_{\mathrm{tot}} &= \rho_\phi + \rho_m, &\qquad p_{\mathrm{tot}} &= p_\phi
+ p_m.
\end{alignat}
\end{subequations}
%

Adding a cosmological constant to a perfect fluid with a barotropic equation
of state can be regarded as a problem with a single fluid with a changed
barotropic equation of state.
%, e.g., in the case of $\Lambda$CDM cosmology
%$\gamma_{\mathrm{tot}}$ is given by $\gamma_{\mathrm{tot}} = \rho_m/(\rho_m +
%\Lambda)$, which yields $\gamma_{\mathrm{tot}} =
%\Omega_{m0}\exp(-3\tau)/(\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0})$.
Replacing a cosmological constant in a general barotropic fluid description
with a scalar field therefore adds the dynamical degrees of freedom, $\phi$
and $\dot{\phi}$, to $\rho_m$ (or $a$), which implies a complication, since
this leads to a larger set of solutions with different histories. In the
context of $\Lambda$CDM deformations this problem is somewhat alleviated by
the following assumption: To obtain an evolutionary history such that primordial 
nucleosynthesis gives observationally compatible light
nuclear abundances, we will assume that there exists a mechanism such as
inflation that sets initial conditions for the present models at high
redshifts so that the matter content greatly dominates over the scalar field
(dark energy) content in the early Universe.
%We stress that our main goal is to consider as simple classes of models as
%possible in order to clearly display general features that lead to open sets
%of $\Lambda$CDM deformed families of solutions.


For simplicity we will from now on consider a linear equation of state with a
constant $\gamma_m$ restricted to the range $0<\gamma_m<2$, which avoids
bifurcations at $\gamma_m=0$, corresponding to a cosmological constant, and
$\gamma_m=2$, which yields a stiff perfect fluid with the speed of sound
equal to the speed of light; $\gamma_m=1$ corresponds to a fluid without
pressure, which will be the focus when describing and discussing $\Lambda$CDM
deformations, while $\gamma_m=4/3$ corresponds to a radiation fluid. It
follows from Eq.~\eqref{rhoeq} that
%
\begin{equation}\label{rhosol}
\rho_m=\rho_0(a/a_0)^{-3\gamma_m}.
\end{equation}
%

Furthermore, most of the scalar field potentials that have been 
considered as possible alternatives to a cosmological constant as dark energy
have strictly monotonic scalar field potentials, as exemplified by, e.g., the popular 
inverse power-law potential.  
%Furthermore, since we are primarily interested in
%continuous scalar field deformations of $\Lambda$CDM cosmology this requires
%that the potential $V(\phi)$ is such that it does not give rise to solutions
%with different qualitative character, which occurs if the potential has
%extrema and/or inflection points. 
Apart from the model with $V=\Lambda = \mathrm{const.}$, we 
will therefore assume that $V$ is defined, differentiable,
positive, and strictly monotone for $\phi \in (\phi_-,\phi_+)$. Without loss
of generality, we assume that $V$ is monotonically decreasing. We also assume
that the potential is such that $\phi_+=\infty$. If $\phi_-$ is finite or
$-\infty$ depends on how fast $V(\phi)$ increases when $\phi\rightarrow
\phi_-$, e.g., for an inverse power-law potential, $V \propto
\phi^{-\alpha}$, $\phi_-=0$, while an exponential potential, $V \propto
\exp(-\lambda \phi)$, leads to $\phi_-=-\infty$.

A constant or monotonically strictly decreasing potential makes it convenient
to use a particular type of dynamical systems formulation that brings
%
\begin{equation}
\lambda(\phi) = - \frac{V_\phi}{V}
\end{equation}
%
into focus, where $\lambda$ is zero for a constant potential and positive for
$\phi \in (\phi_-,\infty)$ when $V(\phi)$ is strictly monotonically
decreasing. Within this context it is natural to regard $\lambda$ as more
``fundamental'' than $V$ itself, where $V$ is obtained from $\lambda(\phi)$ via
%
\begin{equation}\label{Vfromlambda}
V = V_0\exp\left(-\int_{\phi_0}^\phi \lambda(\tilde{\phi}) d\tilde{\phi}\right).
\end{equation}
%
In the special case where $\lambda$ and $\lim_{\phi \rightarrow
\phi_-}\lambda$ are bounded it follows that the potential can be bounded by
an exponential and therefore $\phi_- = -\infty$. If $\lambda$ is unbounded,
as in the case of an inverse power-law potential, $\phi_-$ is finite.

%Here, however, we will only attempt to shed some light on the above `field
%theoretical $\Lambda$CDM-deformation issue' with some simple quintessence
%examples by means of a dynamical systems approach, which also will situate
%some previous work in a hopefully useful new context.

Finally we note that the system~\eqref{Hphieq} has unbounded variables and 
right-hand sides that blow up. It is therefore not suitable for 
a global or asymptotic analysis of
the solution space and its properties. Below we will therefore make a change of
variables in order to obtain a dynamical system on a compact state space that
allows these issues to be addressed, as well as giving global illustrative
pictures of the entire solution spaces of the models under
consideration.\footnote{For a few examples of work dealing with dynamical
systems methods in cosmology,
see~\cite{waiell97,col03,giamir10,ugg13,leofad14,tam14,alhugg15,alhetal15},
and references therein.} 

The outline of the paper is as follows. In the next section we (i) formulate
the dynamical systems approach that is used throughout the paper to deal with
a fluid and a minimally coupled scalar field; (ii) describe the state space
features; (iii) point out some global properties and determine asymptotic
behavior for models with a perfect fluid and a scalar field such that
$0<\lambda(\phi) <\infty$, $\phi \in (\phi_-,\infty)$, and
%
\begin{equation}
\lim_{\phi \rightarrow \infty} \lambda = \lambda_+,\qquad \lim_{\phi \rightarrow \phi_-}\lambda = \lambda_-,
\end{equation}
%
where $\lambda_+$ is finite and $\lambda_-$ is finite or $\infty$.

In Sec.~\ref{sec:frozenlambda} we (i) situate $\Lambda$CDM cosmology,
which occurs for $\lambda=0$, in the three-dimensional state space, which
leads to the concept of an attracting separatrix surface rather than
individual ``attractor solutions''; (ii) consider the $\lambda=\mathrm{const.}$
(i.e., an exponential potential) and dust models in the present context and
give a specific example that shows that some of these models are completely
solvable; (iii) discuss the symmetries of the above models and how the
associated structures might be somewhat misleading for the general picture.
%This is partly due to that the symmetries make it possible to obtain
%projected two-dimensional representations of the dynamics of these models,
%which is usually what is done, while it is necessary to consider them in the
%generic three-dimensional state space in order to situate them in a generic
%context.

In Sec.~\ref{sec:dynlambda} we turn from the previous ``frozen $\lambda$''
models to some ``dynamical $\lambda$'' models. We consider (i) the inverse
power-law potentials, characterized by $\lambda= \alpha/\phi$, in a global
dynamical systems setting, tying local results to the global features
discussed in Sec.~\ref{sec:dynsys}; in particular we point out that
observational viability requires fine-tuning of $\alpha$ to small values,
since $\lambda \rightarrow \infty$ when $\phi \rightarrow 0$ produces a
``memory'' that affects the entire evolution of, e.g., ``tracker'' (attractor)
solutions. For this reason we therefore (ii) consider the arguably simplest
possible ``$\lambda$-regularization'' of the inverse power-law potential,
namely, $\lambda = \alpha/\sqrt{C^2 + \phi^2}$, which, in contrast to the
inverse power-law case, gives rise to continuous deformations of the
attracting $\Lambda$CDM separatrix surface, and which for $C \gtrsim \alpha$ 
yields much less stringent observational compatibility
conditions on $\alpha$. Finally, the section is concluded with (iii) a
general discussion about observational viability conditions that $\lambda$
should satisfy for \emph{any} potential that might arise from some fundamental
theory.

Throughout we use dynamical systems pictures which clearly display general
key features, and, in order to compare with $\Lambda$CDM cosmology, these
pictures are complemented by diagrams that indicate the evolutionary history
of physically important quantities. Finally, we outline how it is possible to
treat two fluids, e.g. dust and radiation, and a scalar field in
the Appendix.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Dynamical systems approach}\label{sec:dynsys}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-----------------------------------------------------------------
\subsection{Dynamical systems formulation}
%-----------------------------------------------------------------

The presently used dynamical systems formulation is based on the following
dependent variables,
%
\begin{equation}\label{vardef1}
(x,\Omega_V,\Omega_m,Z) =
\left(\frac{\dot{\phi}}{\sqrt{6}H}, \frac{V}{3H^2}, \frac{\rho_m}{3H^2},Z(\phi)\right).
\end{equation}
%
These variables are by no means suitable for all scalar field potentials,
e.g., for potentials with a zero minimum such as monomial potentials it is
advisable to replace the scalar field variable $Z$ with an $H$-based variable
that takes into account a varying averaged oscillatory time scale at late
times, as done in~\cite{alhetal15}. Nevertheless, the above variables are
useful for positive monotonic potentials. %positive,
%monotone, differentiable potential $V(\phi)$ on an interval $\phi \in
%(\phi_-,\phi_+)$, where $\phi_+$ ($\phi_-$) may or may not be equal to
%$\infty$ ($-\infty$), and where $\lim_{\phi\rightarrow \phi_+}V(\phi)$ has a
%non-negative limit, which we henceforth assume.
An optimal choice of the variable $Z(\phi)$ depends on the potential that is
studied. However, to obtain a suitable global dynamical systems formulation
$Z$ should always be chosen to be a globally invertible monotone function in
$\phi$, defined on a \emph{bounded} interval $Z\in (Z_-,Z_+)$, $Z_\pm =
Z(\phi_\pm)$; without loss of generality, we choose $Z$ to be monotonically
increasing in $\phi$ so that $dZ/d\phi>0$ for $Z\in (Z_-,Z_+)$.

Apart from the dependent variables we also need to choose a new time
variable. In order to obtain a regular dynamical system this choice depends
on if $\lambda$ is bounded or not. In the case $\lambda$ is bounded, as
exemplified by, e.g., exponentially bounded potentials, it is convenient to use
${\tau}$, defined by $dt=H^{-1}d{\tau}$, where ${\tau}= \ln(a/a_0)$ sometimes
is referred to as $N$, the number of $e$-folds from a reference time $t_0$
[$a_0 = a(t_0)$], which for the present time leads to $\tau = -\ln(1+z)$,
where $z$ is the redshift. Thus, the field equations can be
written as the following coupled system:
%
\begin{subequations}\label{dynsys0}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}\lambda(Z)\,\Omega_V,\label{xeq}\\
\Omega_m^\prime &= 3\left[2x^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\label{Omeq}\\
Z^\prime &= \sqrt{6}\frac{dZ}{d\phi}\,x,\label{Zeq0}
\end{align}
\end{subequations}
%
where ${}^\prime$ denotes derivatives with respect to $\tau$. With some slight
abuse of notation we have written $\lambda(Z) = \lambda(\phi(Z)) =
-V_\phi/V$. Furthermore, the Gauss constraint~\eqref{Gauss1} is used to
globally express the Hubble-normalized scalar field potential energy,
$\Omega_V$, in terms of the state space variables:
%
\begin{equation}\label{Gauss2}
\Omega_V = 1 - x^2 - \Omega_m,
\end{equation}
%
while the deceleration parameter $q$, defined via ${H}^\prime = - (1+q)H$, is
given by
%
\begin{equation}\label{q}
q = - 1 + 3x^2 + \frac32\gamma_m\Omega_m.
\end{equation}
%

The above system is suitable when $\lambda$ is bounded, but not when
$\lim_{\phi \rightarrow \phi_-}\lambda = \infty$, as exemplified by, e.g.,
inverse power-law potentials. In order to obtain a regular dynamical system
when $\lim_{\phi \rightarrow \phi_-}\lambda = \infty$ we therefore choose a
new time variable $\bar{\tau}$, defined by $d\tau = g(Z)d\bar{\tau}$, where
$g(Z)$ is a suitable positive bounded function of $Z$ such that
$\lim_{Z\rightarrow Z_-}g(Z)=0$,
%
\begin{equation}
\lim_{Z\rightarrow Z_-}g_\lambda(Z)  = g_{\lambda_-} =\mathrm{const}. < \infty, \quad \text{where} \quad
g_\lambda(Z) := g(Z)\lambda(Z),
\end{equation}
%
and $\lim_{Z\rightarrow Z_+}g(Z)=1$, where, e.g., $g(Z) =1/(1 + \lambda(Z))$ is
a suitable choice for inverse power-law potentials. This results in the
system
%
\begin{subequations}\label{dynsysZ}
\begin{align}
\frac{dx}{d\bar{\tau}} &= -g(Z)(2-q)x + \sqrt{\frac{3}{2}}g_\lambda(Z)\Omega_V,\label{xeqg}\\
\frac{d\Omega_m}{d\bar{\tau}} &= 3g(Z)\left[2x^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\label{Omeqg}\\
\frac{dZ}{d\bar{\tau}} &= \sqrt{6}g(Z)\frac{dZ}{d\phi}\,x.\label{Zeq0g}
\end{align}
\end{subequations}
%
Note that~\eqref{dynsys0} is obtained by setting $g(Z)=1$, hence
$\bar{\tau}=\tau$, and so the above system formally covers both cases.


%-----------------------------------------------------------------
\subsection{State space structures}
%-----------------------------------------------------------------

The above assumptions lead to the following domain for the state vector
$(x,\Omega_m,Z)$ for the \emph{interior state space} ${\bf S}$ for models
with a perfect fluid and a scalar field:
%
\begin{equation}
\Omega_m>0, \qquad \Omega_V = 1- x^2 - \Omega_m > 0, \qquad Z_-<Z<Z_+.
\end{equation}
%
Furthermore, it follows from~\eqref{dynsys0}, and the following auxiliary
equation for $\Omega_V$,
%
\begin{equation}\label{OmV}
\Omega_V^\prime = 2\left(1 + q - \sqrt{\frac32}\lambda(Z)\,
x\right)\Omega_V,
\end{equation}
%
that it is possible to extend the state space and include the following
boundaries:
%
\begin{itemize}
\item The boundary $\Omega_m=0$ (i.e. $\rho_m=0$) yields the (pure)
    \emph{scalar field boundary subset} [see Eq.~\eqref{Omeq}], and thus
    %
    \begin{equation}\Omega_\phi = x^2 + \Omega_V =1,
    \end{equation}
    %
    with an interior state space ${\bf S}_\phi$ given by $\Omega_V = 1 -
    x^2 > 0$, $Z_-<Z<Z_+$.
\item The boundary $\Omega_V = 1 - x^2 - \Omega_m=0$ corresponds to a
    model with a perfect fluid and a massless scalar field; we refer to
    this subset as the \emph{massless scalar field boundary subset}. Note
    that this boundary also contains the pure massless scalar field
    boundary subset $\Omega_m=0$, $x=\pm 1$, and the \emph{perfect fluid
    boundary subset} $\Omega_m=1$ ($\Omega_\phi=0$, i.e.,
    $\Omega_V=0=x$), which is described by a line of fixed points,
    $\mathrm{FL}_Z$, with constant $Z$.
\item In addition we assume that the potential and $Z$ are such that the state space ${\bf S}$
can be extended to not only include the boundaries $\Omega_m=0$,
$\Omega_V=0$, but also the boundaries $Z=Z_-$ and $Z=Z_+$, which furthermore
are assumed to constitute invariant boundary subsets.
\end{itemize}
%

%With the exception of the inverse power-law potential, we assume that
%$\lambda$ is \emph{bounded} when $\phi\in (\phi_-,\phi_+)$, $\phi \rightarrow
%\phi_\pm$. In addition
%We , so that 
Thus the \emph{extended compact state spaces} $\bar{\bf S}$ and $\bar{\bf S}_\phi$ are characterized by
%
\begin{subequations}
\begin{align}
\Omega_m &\geq 0, &\quad \Omega_V &= 1 - x^2 - \Omega_m \geq 0, &\quad &Z_-\leq Z \leq Z_+ ,\\
\Omega_V &= 1- x^2 \geq 0, &\quad & Z_-\leq Z \leq Z_+, & &
\end{align}
\end{subequations}
%
respectively, and the dynamical systems on $\bar{\bf S}$ and $\bar{\bf
S}_\phi$ are of differentiability class $C^1$ or higher. We therefore assume
that the potential is such that it is possible to find a scalar field
variable $Z$ and, if needed, a function $g(Z)$ so that
%
\begin{itemize}
\item[(i)] $g_\lambda(Z)$ \emph{and} $g(Z)dZ/d\phi$ are bounded and
    differentiable on the extended interval $Z\in [Z_-,Z_+]$,
\item[(ii]) $g(Z)dZ/d\phi|_{Z=Z_\pm}=0$, i.e., $Z=Z_-$ and $Z=Z_+$ are
    invariant boundary subsets.
\end{itemize}
%

Although the above restrictive assumptions require that the potential
$V(\phi)$ is of differentiability class $C^2$, they cover a vast class of
potentials $V(\phi)$. Moreover, the extended state space treatment of the
admissible class of potentials can be used to provide bounds on asymptotic
behavior of potentials that do not admit $Z$-extended state spaces, and thus
they provide a natural starting point for very general situations. It should
also be noted that for many qualitative dynamical aspects, as will be seen,
there is no need to state any more information about $Z(\phi)$ or $V(\phi)$,
although quantitative treatments of course need specific functions and thus
we will give specific examples when dealing with continuous $\Lambda$CDM
deformations.

To physically compare solutions that are obtained in the dynamical systems
approach with $\Lambda$CDM cosmology, we will complement the state space
description with diagrams involving physically important quantities.
%To physically understand the dynamics given by the state space formulation we
%complement it by considering physically interesting quantities.
%Once one knows the dynamics it has to be physically interpreted. This can be
%done by plotting physically interesting quantities.
Apart from the physically important state space variable $\Omega_m$ we have
given the deceleration parameter $q$ in terms of state space variables in
Eq.~\eqref{q}, and below we will also express $H$ in state space variables.
In addition it is of interest to give the jerk parameter
$j$~\cite{blaetal04,rapetal07}, since even though the observational
feasibility of $j$ is questionable, it is still useful in order to
theoretically compare models with $\Lambda$CDM cosmology for which $j=1$; for
constant $\gamma_m$ the jerk parameter is given by
%
\begin{equation}
j = 1 + 9x^2 - \frac92\gamma_m(1-\gamma_m)\Omega_m - 3\sqrt{6}\lambda(Z)\,\Omega_V\,x.
\end{equation}
%

Next we turn to structures that are helpful when it comes to determining global
and asymptotic state space features.

%-----------------------------------------------------------------
\subsection{Nonlocal features}
%-----------------------------------------------------------------

%-----------------------------------------------------------------
\subsubsection*{The monotonicity of $H$}
%-----------------------------------------------------------------

The Hubble parameter $H$, or more conveniently $H^2$, can be expressed in the
state space variables in ${\bf S}$ and ${\bf S}_\phi$ and is given by
%
\begin{equation}\label{H2}
H^2 = \frac{V(\phi(Z))}{3\Omega_V} = \frac{V(\phi(Z))}{3(1 - x^2 - \Omega_m)},
\end{equation}
%
as follows from the definition of $\Omega_V$. This has far-reaching
consequences for the global dynamics due to
%
\begin{equation}
H^\prime = - (1+q)H \quad \Rightarrow  \quad (H^2)^\prime = - 2(1+q)H^2,
\end{equation}
%
since Eqs.~\eqref{Gauss2} and~\eqref{q} lead to
%
\begin{equation}\label{eqqineq}
1+q=3x^2 + \frac32\gamma_m\Omega_m, \qquad
2-q=3\Omega_V + \frac32(2-\gamma_m)\Omega_m,
\end{equation}
%
which together with $0<\gamma_m<2$ yield
%
\begin{equation}
-1\leq q\leq 2.
\end{equation}
%
As will be seen, the above inequalities for $q$, induced by energy and causality conditions,
play a central role for the dynamics, which strongly suggests that the field
equations should be expressed in terms of $q$, as done in~\eqref{dynsys0}.

The above relations imply that $q=2$ only occurs on the pure massless scalar
field subsets $x=\pm 1$, while $q=-1$ requires $\Omega_m=x=0$, $\Omega_V=1$.
Since $1+q>0$ on ${\bf S}$, $H^2$ is a monotonically decreasing function on
${\bf S}$. This is also true for ${\bf S}_\phi$, as can be seen as follows.
On ${\bf S}_\phi$ we have the inequalities $-1\leq q \leq 2$. As a
consequence $H^2$ is a monotonically decreasing function also on this subset,
except when $1+q=0$. However, when $1+q=0$ is not an invariant subset, then
$1+q=0$ just corresponds to an inflection point in the graph of $H^2$, and
$H^2$ is thereby still monotonically decreasing.\footnote{For an elaboration
of this in the case of monomial potentials, see~\cite{alhetal15}.} It is only
when $1+q=0$ is an invariant subset that $H^2$ stops decreasing and for which
$1+q=0$ can be an asymptotic state. For this to happen requires that $x=0$,
since $1+q=3x^2$ on ${\bf S}_\phi$, also means that $x=0$ must be
an invariant subset. Since $x^\prime = \sqrt{\frac{3}{2}}\lambda(Z)$ when
$\Omega_m=0=x$, $\Omega_V=1$, this leads to the condition that
$\lambda(Z)=0$, which by our assumptions about monotonically decreasing
potentials can only happen on the $Z_\pm$ boundaries. Hence $H^2$ is also
monotonically decreasing on ${\bf S}_\phi$. This implies that there can be no
fixed points or recurring orbits (solution trajectories) in ${\bf S}$ and
${\bf S}_\phi$, and due to~\eqref{H2} all solution trajectories in ${\bf S}$
and ${\bf S}_\phi$ originate and end at $Z_\pm$ and $\Omega_V=0$, which shows
that it is essential to extend the state space to include these boundaries 
in order to fully describe the dynamics of the present models .

%-----------------------------------------------------------------
\subsubsection*{The $Z_\pm$ boundaries}
%-----------------------------------------------------------------

Let us begin by considering the $Z_\pm$ boundaries for the case where
$\lim_{Z\rightarrow Z_\pm}\lambda = \lambda_\pm$ is finite. Dropping the
subscript $\pm$ on $\lambda_\pm$ leads to the system
%
\begin{subequations}\label{dynsysred}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}\lambda\,\Omega_V,\label{xeqp}\\
\Omega_m^\prime &= 3\left[2x^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\label{Omeqp}
\end{align}
\end{subequations}
%
on $Z_\pm$ where we recall that $\Omega_V = 1 - x^2 - \Omega_m$ and  $q = - 1
+ 3x^2 + \frac32\gamma_m\Omega_m$. As will be elaborated on in the next
section, these equations describe the reduced equations of an exponential
potential $V \propto \exp(-\lambda\phi)$ when $\lambda \neq 0$ and those of a
constant potential when $\lambda=0$. The above system has an extended state
space $\bar{\bf S}_\mathrm{red}$ given by
%
\begin{equation}
\Omega_m\geq 0, \qquad \Omega_V = 1 - x^2 - \Omega_m \geq 0,
\end{equation}
%
and admits a number of fixed points given in Table~\ref{tab:fixedpoint}.
%
\begin{table}[ht!]
\begin{center}
\begin{tabular}{|c|c|c|c|}\hline
% &  Fixed point values & &  \\ \hline
$\mathrm{FL}$ & $x = 0$, $\Omega_m = 1$ & $q = (3\gamma_m-2)/2$ & - \\ \hline
$\mathrm{EM}$ & $x = \sqrt{3/2}\gamma_m\lambda^{-1}$, $\Omega_m = 1 - 3\gamma_m\lambda^{-2}$; $\lambda >\sqrt{3\gamma_m}$ & $q =(3\gamma_m-2)/2$
& $\gamma_\phi=\gamma_m$ \\ \hline
$\mathrm{M}^\pm$ & $x = \pm 1$, $\Omega_m =0$ & $q = 2$ & $\gamma_\phi=2$ \\ \hline
$\mathrm{PL}$ & $x = \lambda/\sqrt{6}$, $\Omega_m = 0$; $\lambda <\sqrt{6}$ & $q = -1 + \lambda^2/2$ & $\gamma_\phi = \lambda^2/3$ \\ \hline
$\mathrm{dS}$ & $x = \Omega_m = 0$; $\lambda=0$ & $q = -1$ & $\gamma_\phi =0$ \\ \hline
\end{tabular}
\end{center}
\caption{Table depicting the type of fixed points that can occur for a FLRW
model with a monotonically decreasing scalar field potential with a finite
$\lambda$ and a perfect fluid with a linear equation of state parameter
$0<\gamma_m<2$.} \label{tab:fixedpoint}
\end{table}
%

We have here introduced a notation for the fixed points where the kernel
$\mathrm{M}$ stands for a massless scalar field state; $\mathrm{FL}$ stands
for a Friedmann-Lema\^{i}tre perfect fluid state, while $\mathrm{dS}$ stands
for a de Sitter state (for interpretation of various types of de Sitter
states in scalar field cosmology, see~\cite{alhetal15}); the kernel
$\mathrm{PL}$ stands for power law (inflation, when $\lambda <\sqrt{2}$, since
this leads to an accelerating state with $q<0$), while $\mathrm{EM}$ stands
for a scaling solution with an exponential potential and matter in the form of a
perfect fluid with a linear equation of state. In addition a superscript
describes the value of $x$, while in the full three-dimensional treatment we
also add a subscript that describes the value for $Z$.

A one-parameter set of orbits in $\bar{\bf S}_\mathrm{red}$ 
originates from each of the sources $\mathrm{M}^+$ and $\mathrm{M}^-$ when
${\lambda}<\sqrt{6}$, but when ${\lambda}>\sqrt{6}$ there are no solutions
that originate from $\mathrm{M}^+$ into ${\bf S}_\mathrm{red}$,
while one solution originates from $\mathrm{FL}$. Toward the future
$\mathrm{PL}$ is a sink when ${\lambda} < \sqrt{3\gamma_m}$
%(and the future attractor on $\bar{\bf S}_\mathrm{red}$
while $\mathrm{EM}$ is future stable when ${\lambda} > \sqrt{3\gamma_m}$.
%(and the future attractor on $\bar{\bf S}_\mathrm{red}$).
Thus bifurcations take place when ${\lambda}=\sqrt{6}$, which is when
$\mathrm{PL}$ leaves the physical state space via $\mathrm{M}^+$, and when
${\lambda}=\sqrt{3\gamma_m}$, which is when $\mathrm{EM}$ enters the state
space via $\mathrm{PL}$. In addition there is a bifurcation when $\lambda=0$,
which is when $\mathrm{PL}$ becomes the future sink $\mathrm{dS}$. The
solution space structures for representative examples of these different
cases are given in Fig.~\ref{fig:redrep}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[Solution space for $\lambda=0$ and $\gamma_m=1$.]{\label{fig:SS_CC2D}
\includegraphics[width=0.45\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_CC_2D.pdf}} \quad
\subfigure[Solution space for ${\lambda}=\sqrt{\frac{3}{2}}$ and $\gamma_m=1$.]{\label{fig:SSCC2D}
\includegraphics[width=0.45\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP2D.pdf}}\\
\subfigure[Solution space for ${\lambda}=2$ and $\gamma_m=1$.]{\label{fig:SSCCnum}
\includegraphics[width=0.45\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP_Lambda2_2D.pdf}} \qquad
\subfigure[Solution space for ${\lambda}=3$ and $\gamma_m=1$.]{\label{fig:SSCCLamda2}
\includegraphics[width=0.45\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP_Lambda3_2D.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{Representative solution space structures for $\bar{\bf
S}_\mathrm{red}$.}
\label{fig:redrep}
\end{figure}
%

The case when $\lambda \rightarrow \infty$ when $\phi \rightarrow \phi_-$
yields the following system on $Z_-$:
%
\begin{subequations}
\begin{align}
\frac{dx}{d\bar{\tau}} &= \sqrt{\frac{3}{2}}g_{\lambda_-}\,\Omega_V,\label{xm}\\
\frac{d\Omega_m}{d\bar{\tau}} &= 0,\label{Omm}
\end{align}
\end{subequations}
%
where we recall that $g_{\lambda_-} = \lim_{Z\rightarrow Z_-}
(g(Z)\lambda(Z))>0$. It follows that $\Omega_m = \mathrm{const}.$ while the
subset $\Omega_V=1-x^2-\Omega_m=0$ forms a line of fixed points
$\mathrm{M}^x$, and that $x$ is monotonically increasing in ${\bf
S}_\mathrm{red}$. Thus $\mathrm{M}^x$ with $x<0$ is the source while
$\mathrm{M}^x$ with $x>0$ is the sink of all orbits in ${\bf
S}_\mathrm{red}$; finally, note that $\mathrm{M}^x$ with $x=0$ is the fixed
point $\mathrm{FL}$.

%-----------------------------------------------------------------
\subsubsection*{The $\Omega_V=0$ boundary}
%-----------------------------------------------------------------

We first note that the equation for $x$ on $\Omega_V=0$ is invariant under
the transformation $x\rightarrow - x$, which enables one to use
$\Omega_\mathrm{stiff}=x^2$ as a variable instead of $x$, and as a result
the dynamics of the massless scalar field explicitly takes the same form as
that of a stiff fluid, $p=\rho$. This leads to the equations
%
\begin{subequations}\label{dynsysOmV}
\begin{align}
\Omega_m^\prime &= 3(2-\gamma_m)(1 - \Omega_m)\Omega_m,\label{OmeqV}\\
Z^\prime &= \sqrt{6}\,\frac{dZ}{d\phi}\,x = \pm\sqrt{6}\,\frac{dZ}{d\phi}\,\sqrt{1-\Omega_m},\label{ZeqV}
\end{align}
\end{subequations}
%
where we have used the Gauss constraint to write $\Omega_\mathrm{stiff}= 1 -
\Omega_m$ and $x=\pm \sqrt{\Omega_\mathrm{stiff}}$. Since the equation for
$\Omega_m$ decouples from that of $Z$, it follows that all orbits (solution
trajectories) in the two-dimensional state space look the same as the orbits
of the one-dimensional problem associated with $\Omega_m$, when projected
onto $\Omega_m$. As follows from~\eqref{OmeqV}, $\Omega_m$ monotonically
increases from $0$ to $1$. Thus the past asymptotic state of all orbits on
$\Omega_V=0$ resides on the pure massless scalar field subset $\Omega_m=0$,
$x^2=1$ and thereby $q=2$; thus either $x=+1$, which from~\eqref{ZeqV} leads to
$Z=Z_-$ and thus a fixed point $\mathrm{M}^+_{Z_-}$, since $dZ/d\phi>0$, or
$x=-1$, which yields $Z=Z_+$ and a fixed point $\mathrm{M}^-_{Z_+}$. It also
follows that the future asymptotic state of all orbits on $\Omega_V=0$
resides on $\Omega_m=1$, $x=0$, $Z=\mathrm{const}.$, which thereby form a
line of fixed points $\mathrm{FL}_Z$. To determine what happens with $Z$
toward the future on $\Omega_V=0$ we note that it is possible to solve the
equation for $\Omega_m$ explicitly, but the solution can also be obtained by
using the fact that $\rho_\mathrm{stiff} \propto \exp(-6\tau)$, $\rho_m \propto
\exp(-3\tau)$, $3H^2 = \rho_\mathrm{stiff} + \rho_m$, which gives
%
\begin{subequations}
\begin{align}
\frac{H}{H_0} &= \left(\Omega_{\mathrm{stiff}0}\exp(-6\tau)+ \Omega_{m0}\exp(-3\gamma_m\tau)\right)^{1/2},\\
x &= \pm \Omega_{\mathrm{stiff}0}^{1/2}\,\exp(-3\tau)/(H/H_0), \\
\Omega_m &= \Omega_{m0}\exp(-3\gamma_m\tau)/(H/H_0)^2.
\end{align}
\end{subequations}
%
To obtain the solution for $Z$ requires a specification of $Z$. However, it
is easier to obtain $Z$ by integrating the equation for $\phi$,
%
\begin{equation}
\phi^\prime = \sqrt{6}x = \pm(6\Omega_{\mathrm{stiff}0})^{1/2}\,\exp(-3\tau)/(H/H_0),
\end{equation}
%
and then inserting the result in the chosen definition of $Z$. As follows
from the above equation, $\phi \rightarrow \mathrm{const.}$ when $\tau
\rightarrow +\infty$, and therefore also $Z\rightarrow \mathrm{const.}$ in
this limit; i.e., each point on the line $\mathrm{FL}_Z$ attracts one
solution with positive and one with negative initial $x$ on the $\Omega_V=0$
boundary. For a representative depiction of the $\Omega_V=0$ boundary, see
Fig.~\ref{fig:V0Dust}.

%-----------------------------------------------------------------
\subsubsection*{Role of the $\Omega_V=0$ boundary}
%-----------------------------------------------------------------

To study the neighborhood of $\Omega_V=0$ we note that
%
\begin{equation}
\Omega_V^{-1}\Omega_V^\prime|_{\Omega_V=0} = 2\left(1 + q - \sqrt{\frac32}\lambda(Z)\, x\right),
\end{equation}
%
and hence
%
\begin{equation}\label{MFL}
\Omega_V^{-1}\Omega_V^\prime|_{\mathrm{M}^\pm_{Z_\mp}} = \sqrt{6}\left(\sqrt{6} \mp \lambda(Z_\mp)\right), \qquad
\Omega_V^{-1}\Omega_V^\prime|_{\mathrm{FL}_Z} = 3\gamma_m.
\end{equation}
%
As a consequence it follows that $\mathrm{M}^-_{Z_+}$ is a source for orbits
in ${\bf S}$ and ${\bf S}_\phi$ as is $\mathrm{M}^+_{Z_-}$ when $\lambda(Z_-)
< \sqrt{6}$, but if $\lambda(Z_-) > \sqrt{6}$, then no orbits enter ${\bf S}$
or ${\bf S}_\phi$ from $\mathrm{M}^+_{Z_-}$;\footnote{Solutions originating
from $\mathrm{M}^-_{Z_+}$ corresponds to ``scalar field particles'' that come
from $\infty$ moving toward decreasing $\phi$, while if $\lambda(Z_-) <
\sqrt{6}$ the potential is shallow enough so that scalar field particles also 
can come from $-\infty$ moving toward increasing $\phi$. Solutions originating
from $\mathrm{FL}_Z$ correspond to scalar field particles that initially lie
still at some initial value $\phi$ and then roll down the potential
$V(\phi)$.} a one-parameter set of orbits, one from each point on
$\mathrm{FL}_Z$, enters ${\bf S}$, where $\mathrm{FL}_Z$ acts as a
``transversal saddle line,'' forming a two-dimensional invariant subset in
${\bf S}$. By demanding initial data for which $\Omega_m$ is close to 1,
this implies the important feature that for \emph{all} the models we consider
such solutions can be approximated by this (at least initially)
\emph{attracting invariant} surface \emph{subset} of codimension 1 with
respect to the state space dimension. We will later contrast this feature
with discussions about \emph{single} ``attracting'' solutions that appear for
various models in the literature.

%-----------------------------------------------------------------
\subsubsection*{Monotonicity of negative $x$}
%-----------------------------------------------------------------

Further restrictions on the global dynamics of the present models come from
eqs.~\eqref{xeq} and~\eqref{Zeq0}. It follows from~\eqref{xeq} that $x$ is
monotonically increasing when $-1<x<0$, since then $2-q>0$. This in turn
implies that if an orbit resides or partly resides in the region $-1<x<0$,
then the orbit originates from $x=-1$, which due to~\eqref{Zeq0} and
$dZ/d\phi>0$ implies that the past asymptotic state must be
$\mathrm{M}^-_{Z_+}$. Furthermore, the monotonicity of $x$ when $-1<x<0$
implies that all orbits in ${\bf S}$ and ${\bf S}_\phi$ must end at $x\geq
0$.

The equation for $x$ can be written on the following form:
%
\begin{equation}
x^\prime = -\frac32\left(2-\gamma_m\right)x\,\Omega_m - 3\left(x - \frac{\lambda(Z)}{\sqrt{6}}\right)\Omega_V,
\end{equation}
%
which shows that $x$ is monotonically decreasing if $1>x >
\lambda(Z)/\sqrt{6}$. In combination with the above result for $x$ this
implies that all orbits in ${\bf S}$ and ${\bf S}_\phi$ must end at $0 \leq x
\leq \lambda(Z)/\sqrt{6}$. However, in order for an orbit to end at $x=0$
requires that $x=0$ be an invariant subset. Since
%
\begin{equation}
x^\prime|_{x=0} = \sqrt{\frac32}\lambda(Z)(1-\Omega_m),
\end{equation}
%
and since $\Omega_m=1$ is the $\mathrm{FL}_Z$ subset, it follows that this
can only be the case if $\lambda(Z)=0$, which can only happen on the $Z_\pm$
boundaries for a monotonically strictly decreasing potential. The above
equation also shows that $x$ is monotonically increasing when $x=0$ when
$Z_-<Z<Z_+$. Furthermore, since $dZ/d\phi>0$ it follows from~\eqref{Zeq0}
that $Z$ is monotonically increasing (decreasing) in $\tau$, or $\bar{\tau}$,
when $x>0$ ($x<0$).

If $x$ is future asymptotically positive it follows that $Z\rightarrow Z_+$.
However, to exclude or show that any of the orbits from $\mathrm{M}^-_{Z_+}$
end up at $\mathrm{dS}_{Z_-}$ if $\lambda_- =0$ requires more information
about the potential. If $\lim_{\phi \rightarrow \phi_-} V(\phi) =
\mathrm{const}.$, then there are solutions that end at $Z_-$, but if $V
\rightarrow \infty$ when $\phi \rightarrow \phi_-$, there are not. This can
be understood by considering~\eqref{KG}. A solution with negative $x$ can be
viewed as a particle moving in a direction where the potential is monotonically 
increasing. At the same time it is losing energy due to the friction force 
$-3H\dot{\phi}$ until $\dot{\phi}=0$,
%($V\neq 0$ implies $H\neq 0$
%due to the Gauss constraint~\eqref{Gauss1},),
but then $\ddot{\phi} = -V_\phi$, where $V_\phi$ is strictly negative.
%, also
%in the limit $\phi_-$ when $V \rightarrow \infty$, and 
As a consequence the particle changes direction and $\dot{\phi}$ and thereby also $x$ become
positive, and therefore $Z\rightarrow Z_+$ toward the future.

Thus all orbits in ${\bf S}$ and ${\bf S}_\phi$ end at the invariant boundary
subset $Z_+$ when the potential is monotonically decreasing and $V
\rightarrow \infty$ when $\phi \rightarrow \phi_-$, which we for simplicity
from now on assume (apart from the case of a constant potential, which is
treated separately, all our explicit examples obey these conditions).

%-----------------------------------------------------------------
\subsubsection*{The role of the $Z_-$ subset}
%-----------------------------------------------------------------

The four dynamically qualitatively different cases on $Z_-$ when $0\leq
\lambda_- < \infty$ lead to different interior dynamics:
%
\begin{itemize}
\item[(i)] When $\lambda_-<\sqrt{3\gamma_m}$, a two-parameter
    (one-parameter) set of orbits originates from the source
    $\mathrm{M}_{Z_-}^+$ into ${\bf S}$ (${\bf S}_\phi$), while a single
    orbit originates from $\mathrm{PL}_{Z_-}$ into ${\bf S}_\phi$.
\item[(ii)] When $\sqrt{3\gamma_m} < \lambda_- < \sqrt{6}$, again a
    two-parameter (one-parameter) set of orbits originates from the
    source $\mathrm{M}_{Z_-}^+$ into ${\bf S}$ (${\bf S}_\phi$), but in
    this case there is also a one-parameter set (a single orbit) that
    enters ${\bf S}$ (${\bf S}_\phi$) from $\mathrm{PL}_{Z_-}$ as well as
    a single orbit from $\mathrm{EM}_{Z_-}$.
\item[(iii)] When $\lambda_->\sqrt{6}$, a single orbit enters ${\bf S}$
    from $\mathrm{EM}_{Z_-}$.
\item[(iv)] When $\lambda_-=0$, a two-parameter (one-parameter) set of orbits originates
from the source $\mathrm{M}_{Z_-}^+$ into ${\bf S}$ (${\bf S}_\phi$), while a
single orbit originates from $\mathrm{dS}_{Z_-}$ into ${\bf S}_\phi$.
\end{itemize}
%

When $\lambda_- \rightarrow \infty$ at $Z_-$, there is just one orbit that
enters ${\bf S}$ from $\mathrm{FL}_{Z_-}$, the so-called tracker or attractor
solution. The exclusion of orbits into ${\bf S}$ and ${\bf S}_\phi$ from the
fixed points $\mathrm{M}^x_{Z_-}$ with $x\neq 0$ can be understood as
follows: There are no solutions coming from $\mathrm{M}^x_{Z_-}$ when $x<0$
since $Z$ is decreasing toward the future when $x<0$. When $x>0$, it follows
that $\Omega_V$ is decreasing in the vicinity of $\mathrm{M}^x_{Z_-}$ and
hence there are no solutions coming from this part of $\mathrm{M}^x_{Z_-}$
into ${\bf S}$ or ${\bf S}_\phi$ either. Hence $\mathrm{M}^x_{Z_-}$ with
$x=0$, i.e., $\mathrm{FL}_{Z_-}$ is the only fixed point on
$\mathrm{M}^x_{Z_-}$ which can give rise to a solution into ${\bf S}$, and it
does give rise to the tracker solution, which can be established by center
manifold theory or by consideration of nearby solutions and continuity, since
this fixed point acts as a kind of saddle in the full state space.

%-----------------------------------------------------------------
\subsubsection*{The role of the $Z_+$ subset}
%-----------------------------------------------------------------

The only solutions that originate from the $Z_+$ subset into ${\bf S}$ (${\bf
S}_\phi$) are those from the source $\mathrm{M}_{Z_+}^-$, from which a
two-parameter (one-parameter) set of orbits originates into ${\bf S}$ (${\bf
S}_\phi$).

For the present monotonically decreasing potentials, all orbits in ${\bf S}$
and ${\bf S}_\phi$ end at $Z_+$ as follows when $0 \leq\lambda_+<\infty$:
%
\begin{itemize}
\item[(i)] When $\lambda_+<\sqrt{3\gamma_m}$, a two-parameter
    (one-parameter) set of orbits ends at the global sink
    $\mathrm{PL}_{Z_+}$ from ${\bf S}$ (${\bf S}_\phi$).
\item[(ii)] When $\sqrt{3\gamma_m} < \lambda_+ < \sqrt{6}$,  a
    two-parameter set of orbits in ${\bf S}$ ends at the sink
    $\mathrm{EM}_{Z_+}$ while a one-parameter set of orbits in ${\bf
    S}_\phi$ end at $\mathrm{PL}_{Z_+}$.
\item[(iii)] When $\lambda_+>\sqrt{6}$, a two-parameter set of orbits in
    ${\bf S}$ ends at the sink $\mathrm{EM}_{Z_+}$ while a one-parameter
    set of orbits in ${\bf S}_\phi$ end at $\mathrm{M}_{Z_+}^+$.
\item[(iv)] When $\lambda_+=0$, a two-parameter (one-parameter) set of orbits 
    ends at the global sink $\mathrm{dS}_{Z_+}$ from ${\bf S}$ (${\bf S}_\phi$).
\end{itemize}
%

%-----------------------------------------------------------------
\subsubsection*{An additional useful quantity}
%-----------------------------------------------------------------

Finally, consider
%
\begin{equation}
\xi = \frac{\Omega_m}{\Omega_V},
\end{equation}
%
which obeys
%
\begin{equation}
\xi^\prime = -\sqrt{6}\left(\sqrt{\frac32}\,\gamma_m -
\lambda(Z)\,x\right)\xi,
\end{equation}
%
as follows from~\eqref{Omeq} and~\eqref{OmV}. When $\sqrt{\frac32}\,\gamma_m
- \lambda(Z)\,x>0$, %, which e.g. is the case when $x$ is negative or zero or
%when $\lambda(Z) < \sqrt{\frac32}\gamma_m$,
$\xi$ is monotonically decreasing. For globally bounded $\lambda$ such that
$\lambda(Z) < \sqrt{\frac32}\gamma_m$ it follows that $\xi$ is strictly
monotonically decreasing and hence $\xi \rightarrow 0$ and thereby
$\Omega_m\rightarrow 0$ toward the future; i.e., the future dynamics resides
on the scalar field subset. Toward the past $\xi \rightarrow \infty$ and
hence $\Omega_V\rightarrow 0$ and as a consequence it follows that the past
dynamics asymptotically resides on the $\Omega_V=0$ subset, which of course
is already implied by our previous discussion in a somewhat more complicated
manner, which, however, includes more details and general situations.

We will now turn to explicit examples of potentials and a fluid with a dust
equation of state, $\gamma_m=1$. Throughout we will choose scalar field
variables so that $Z_-=0$ and $Z_+ = 1$, so that $\bar{\bf S}$ is given by
%
\begin{equation}
\Omega_m\geq 0, \qquad \Omega_V = 1 -x^2 - \Omega_m \geq 0, \qquad 0\leq Z \leq 1.
\end{equation}
%
It should be pointed out that this type of choice of scalar field variable is
by no means an optimal choice for all scalar field potentials that can be dealt with
by means of a scalar field variable $Z(\phi)$, as illustrated by the treatment of the
generalized Chaplygin gas in~\cite{ugg13}, where $-1\leq Z \leq 1$, but
choices such that $0\leq Z \leq 1$ are ``good enough'' for our present
purposes. Since $Z_-=0$ and $Z_+=1$ we will introduce $0$ and $1$ as
subscripts instead of $Z_\pm$ to distinguish fixed points on the two
boundaries (in addition we have the subscript $Z$ for the line of fixed
points $\mathrm{FL}_Z$ and, for a constant potential, the line of fixed
points $\mathrm{dS}_Z$).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Frozen $\lambda$CDM dynamics}\label{sec:frozenlambda}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Before considering the constant and exponential potentials, let us first
recall some facts about $\Lambda$CDM cosmology. In their simplest form, the
$\Lambda$CDM models are solutions of Einstein's equations with (i) a
spatially isotropic and homogeneous flat FLRW geometry and (ii) a perfect fluid
that has negligible pressure, i.e., dust, representing nonrelativistic
matter, and a cosmological constant $\Lambda$. It is useful to describe
solutions by the following quantities:
%
\begin{subequations}
\begin{align}
H &= \dot{a}/a,\\
q &= -H^{-2}\left(\frac{\ddot{a}}{a}\right) = -\left(1 + \frac{H^\prime}{H}\right),\\
j &= H^{-3}\left(\frac{\dddot{a}}{a}\right) = -q^\prime + q + 2q^2,\\
\Omega_m &= \frac{\rho_m}{3H^2},\\
\Omega_\Lambda &= \frac{\Lambda}{3H^2}.
\end{align}
\end{subequations}
%
Here the cosmographic (or cosmokinetic) quantities $a(t), H(t), q(t), j(t)$
are the cosmological scale factor and the time-dependent Hubble,
%(i.e., $H$ is the time
%dependent Hubble parameter; also, recall that $H = \theta/3$, where $\theta$
%is the expansion of the symmetry surfaces and thereby also the fluid),
deceleration, and jerk parameters, respectively~\cite{blaetal04,rapetal07},
which take the following form for the $\Lambda$CDM models:
%
\begin{subequations}
\begin{align}
\frac{H}{H_0} &= \sqrt{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}},\\
q &= -1 + \frac32\Omega_m = -1 + \frac32\left(\frac{\Omega_{m0}\exp(-3\tau)}{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}}\right),\label{qLCDM}\\
j &= 1,\\
\Omega_m &= \frac{\Omega_{m0}\exp(-3\tau)}{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}}, \\
\Omega_\Lambda &= \frac{\Omega_{\Lambda0}}{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}},
\end{align}
\end{subequations}
%
where $\Omega_m + \Omega_\Lambda = 1$ and where $H_0$, $\Omega_{m0}$, and
$\Omega_{\Lambda0}$ are the values of $H$, $\Omega_{m}$, and
$\Omega_{\Lambda}$ at the reference time $t_0$ and hence
$\tau=0$.\footnote{Due to noisy data, the current observational feasibility
of the jerk parameter $j$ is questionable. Nevertheless, the fact that $j=1$
for the $\Lambda$CDM models makes it useful as a theoretical indicator of
deviations from $\Lambda$CDM cosmology, and it is for this reason we give
$j$. Moreover, as will be seen, the jerk parameter seems to be the parameter
that is most sensitive when it comes to deviations from $\Lambda$CDM
cosmology; it is unfortunate that it is also the most difficult quantity to
observationally measure of the presently discussed ones.} Throughout we will
use the following present values:
%
\begin{equation}
\Omega_{m0} = 0.3, \qquad \Omega_{\Lambda0} = 0.7.
\end{equation}
%
It follows from~\eqref{qLCDM} that $q_0=-0.55$.

Figure~\ref{fig:LCDMgraphs} depicts the evolution of $H/H_0$, $q$, and
$\Omega_m$ in terms of the redshift $z$ for the $\Lambda$CDM models with the
above present values.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$H(z)$-diagram for the $\Lambda$CDM models.]{\label{fig:hLCDM}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_LCDM.pdf}} \,
\subfigure[$q(z)$-diagram for the $\Lambda$CDM models.]{\label{fig:qLCDM}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_LCDM.pdf}}
\subfigure[$\Omega_m(z)$-diagram for the $\Lambda$CDM models.]{\label{fig:OmegaLCDM}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{OmegaM_LCDM.pdf}}\vspace{-0.5cm}
\end{center}
\caption{Depiction of the evolutionary history of $\Lambda$CDM cosmology for
$\Omega_{m0} = 0.3$, $\Omega_{\Lambda0} = 0.7$ in terms of the redshift $z$.}
\label{fig:LCDMgraphs}
\end{figure}
%
As can be seen, $\Omega_m$ monotonically decreases from 1 (corresponding to
the big bang limit $\tau \rightarrow -\infty$) to zero (corresponding to the
infinite future limit  $\tau \rightarrow + \infty$, described by a de Sitter
state with $\Omega_\Lambda=1$). Alternatively this also follows by
considering the differential equation for $\Omega_m$, which is given by (see,
e.g., p. 62 in~\cite{waiell97})
%
\begin{equation}\label{OmprimeLCDM}
\Omega_m^\prime = - 3(1-\Omega_m)\Omega_m.
\end{equation}
%

%-----------------------------------------------------------------
\subsection{$\Lambda$CDM dynamics as $\lambda=0$-CDM dynamics}
%-----------------------------------------------------------------

Setting $p_m=0$ and $V=\Lambda$ yields a model that can be interpreted as that
of dust, a massless scalar field, and a cosmological constant $\Lambda$. The
associated model does not solve any of the issues one might find unattractive
with the $\Lambda$CDM model, but it does show some aspects that are important
in a broader context. Since $\lambda=0$, one obtains a reduced coupled
two-dimensional system for $x$ and $\Omega_m$, while the decoupled equation
for $Z$ can be integrated once the evolution for $x$ and $\Omega_m$ is found.
However, instead of just considering the essential reduced extended
two-dimensional state space $\bar{\bf S}_\mathrm{red}$, we need to consider
the full three-dimensional extended state space $\bar{\bf S}$ to illustrate
how this model is connected with more general models for which $\lambda =
\lambda(Z)$. Note, however, that the decoupling of the equation for $Z$ from
the coupled system of $x$ and $\Omega_m$ results in that all solutions with the
same initial data of $x$ and $\Omega_m$ yield the same curves when projected
onto the $x-\Omega_m$-plane, irrespective of the initial value of $Z$. To
proceed requires that $Z$ be specified. There are many possible choices that
lead to an analytic dynamical system on a compact state space that cover the
domain $\phi_\pm = \pm \infty$, but here, due to future purposes, we choose
%
\begin{equation}\label{Zdef}
Z=\frac{1}{1+\exp(-{\lambda}\phi)},
\end{equation}
%
which is monotonically increasing in $\phi$ and where ${\lambda}$ is an
arbitrary positive constant; the scalar field variable $Z$ thereby has an
extended range $Z\in [0,1]$. This leads to the following dynamical system for
the present models:
%
\begin{subequations}\label{dynsys}
\begin{align}
x^\prime &= -(2-q)x,\label{xeq1}\\
\Omega_m^\prime &= 3\left[2x^2 - (1 - \Omega_m)\right]\Omega_m,\\
Z^\prime &= \sqrt{6}{\lambda} Z(1-Z)x,\label{Zeq1}
\end{align}
\end{subequations}
%
where $2-q=3(1-x^2) - \frac32\Omega_m = 3(1-x^2-\Omega_m) + \frac32\Omega_m$.
The reduced coupled system for $x$ and $\Omega_m$ is invariant under the
transformation $x \rightarrow -x$, which makes it possible to identify this
problem with that of a stiff fluid, characterized by $\Omega_\mathrm{stiff} =
x^2$, a cosmological constant, and dust. In a similar manner as for the
$\Omega_V=0$ subset, this leads to
%
\begin{subequations}
\begin{align}
\frac{H}{H_0} &= \left(\Omega_{\mathrm{stiff}0}\exp(-6\tau) + \Omega_{\Lambda 0} + \Omega_{m0}\exp(-3\tau)\right)^{1/2},\\
x &= \pm \sqrt{\Omega_{\mathrm{stiff}0}}\exp(-3\tau)/(H/H_0), \\
\Omega_m &= \Omega_{m0}\exp(-3\tau)/(H/H_0)^2.
%, \\
%\frac{x^2}{\Omega_m} &\propto \exp(-3\tau).
\end{align}
\end{subequations}
%
As for the $\Omega_V=0$ subset it is easiest to obtain $Z$ by integrating the
equation for $\phi$, given by
%
\begin{equation}
\phi^\prime = \sqrt{6}x = \pm(6\Omega_{\mathrm{stiff}0})^{1/2}\,\exp(-3\tau)/(H/H_0),
\end{equation}
%
and then inserting the solution into~\eqref{Zdef}.

We have previously dealt with the global dynamics of monotonically decreasing
potentials. We here give a complete description for the constant potential. 
The past and future states of all orbits follow from the exact solution, or from 
the monotone function $H^2$ in combination with the dynamical structure at the
boundaries and the local properties of the fixed points. The result is that a
two-parameter set of orbits in ${\bf S}$ and a one-parameter set in ${\bf
S}_\phi$ originate from each of the hyperbolic sources $\mathrm{M}_0^+$ and
$\mathrm{M}_1^-$ while a one-parameter set of solutions originates from the
transversally hyperbolic line of fixed points $\mathrm{FL}_{Z}$ into ${\bf
S}$; orbits in ${\bf S}$ and ${\bf S}_\phi$ end at the transversally
hyperbolic line of fixed points $\mathrm{dS}_{Z}$, which thereby constitute
the future attractor. The orbits that come from $\mathrm{M}_0^+$
($\mathrm{M}_1^-$) correspond to solutions with increasing (decreasing)
$\phi$ from $\phi \rightarrow - \infty$ ($\phi \rightarrow +\infty$),
starting from a massless state [heuristically this can be understood from
Eq.~\eqref{KG} where the friction force $-3H\dot{\phi}$ generates energy
toward the past, making the energy content of $\Lambda$ (and $\rho$) 
negligible compared to the kinetic energy of the scalar field], while
solutions that originate from $\mathrm{FL}_{Z}$ reside on the invariant
subset $x=0$.

Since the $x=0$ subset corresponds to solutions with different constant
$\phi$, and thereby also a constant $\rho_\phi = \Lambda$, these models are
identical to the $\Lambda$CDM models. We therefore denote the $x=0$ invariant
subset surface as the $\Lambda$CDM surface or subset. Furthermore, this
surface divides the state space into two disconnected domains with $x>0$ and
$x<0$, respectively, and the $\Lambda$CDM subset thereby acts as a separatrix
surface. Note that it is a necessary condition that the potential is a
constant for $x=0$ to be an invariant subset in ${\bf S}$ and thus no other
potentials yield a $\Lambda$CDM model.

Due to the regularity of the dynamical system and that each fixed point on
$\mathrm{FL}_{Z}$ acts as a transversal saddle point (exemplifying the
previous general discussion concerning the $\Omega_V=0$ subset and the
vicinity of $\Omega_V=0$ and $\mathrm{FL}_Z$), and since $\mathrm{dS}_{Z}$ is
the future attractor, it follows that the $\Lambda$CDM separatrix surface
$x=0$ is
%attracting nearby solutions,
%indeed, it follows from the dynamics that it is a temporally globally
%attracting surface, i.e., all solutions that are near the surface
%\emph{anywhere} are decreasing their distance to the surface; we therefore
%refer to the $\Lambda$CDM separatrix surface $x=0$ as
an ``attracting'' surface or invariant subset (indeed, $x^2$ is monotonically
decreasing). As a consequence, demanding that $\Omega_m$ is close to 1
initially leads to solutions with almost $\Lambda$CDM behavior. The solution
space for this set of models is depicted in Fig.~\ref{fig:CC} where we have
chosen ${\lambda} = \sqrt{3/2}$.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[The massless scalar field boundary $\Omega_V=0$.]{\label{fig:V0Dust}
\includegraphics[width=0.44 \textwidth, trim = 0cm 2.5cm 0cm 0cm]{V0_Boundary_Final.pdf}}\qquad
\subfigure[The scalar field boundary $\Omega_m=0$.]{\label{fig:M0CC}
\includegraphics[width=0.44\textwidth, trim = 0cm 2.5cm 0cm 0cm]{M0_Boundary_CC.pdf}}
\subfigure[The $(x,\Omega_{m},Z)$ state space for a constant potential and dust.]{\label{fig:SSCC}
\includegraphics[width=0.5\textwidth, trim = 0cm 2.5cm 2cm 0cm]{SS_CC.pdf}} \vspace{0.0cm}
%\subfigure[Projected state space picture onto the $(x,\Omega_{m})$-plane for a constant potential and dust.]{\label{fig:SSCCTop}
%\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{SS_CC_TopView.pdf}} \vspace{-0.5cm}
\end{center}
\caption{Depiction of the solution space of a fluid without pressure and a scalar field
with a constant potential $V(\phi)=\Lambda>0$, where 
$Z=(1+\exp(-\sqrt{3/2}\phi))^{-1}$, illustrating the $x=0$ $\Lambda$CDM
separatrix surface subset and its attracting nature.}
\label{fig:CC}
\end{figure}
%
Note that the structure on the $\Omega_V=0$ subset is qualitatively the same
for all the potentials we consider, as discussed in the previous section.

The present models in effect do not really have a \emph{dynamical} dark
energy since the dark energy content ``freezes'' toward an ``unnatural'' value
$\Lambda>0$ toward the future for models with $x\neq 0$, and they therefore
do not solve any coincidence or energy scale problems one might have with the
$\Lambda$CDM models. Nevertheless, since models with $\Omega_m$ close to 1
as an initial condition basically have the same evolution as the $\Lambda$CDM
models, these models with their \emph{attracting} $\Lambda$CDM
\emph{separatrix surface subset}, in an evolutionary history sense, are the
best a scalar field and dust model can accomplish when it comes to mimicking
$\Lambda$CDM cosmology, and they therefore, from a purely time development
perspective, set the standards for any truly dynamical dark energy model.
Next we turn to the well-known case of an exponential potential (see,
e.g.,~\cite{copetal98,Nunes2000}), but in the present three-dimensional context.


%-----------------------------------------------------------------
\subsection{Constant $\lambda$ $\lambda$CDM dynamics}
%-----------------------------------------------------------------

An exponential potential
%
\begin{equation}
V=V_0\exp(-{\lambda}\phi), \qquad {\lambda}>0,
\end{equation}
%
and a dust matter equation of state results in the dynamical system
%
\begin{subequations}\label{dynsysexp}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}{\lambda}(1 - x^2 - \Omega_m),\label{Sigexp}\\
\Omega_m^\prime &= 3\left[2x^2 - (1 - \Omega_m)\right]\Omega_m,\\
Z^\prime &= \sqrt{6}{\lambda} Z(1-Z)x,
\end{align}
\end{subequations}
%
where we have used the definition~\eqref{Zdef} for $Z$ (also, recall that
$2-q=3(1-x^2) - \frac32\Omega_m$). Once again there is a decoupling between
the scalar field variable $Z$ and a coupled system for $x$ and $\Omega_m$ on
a reduced state space ${\bf S}_\mathrm{red}$, parametrized by ${\lambda}$.
As for a constant potential, it therefore follows that solutions with the
same initial $x$ and $\Omega_m$, but with different initial $Z$, have the
same trajectories when projected onto the $x-\Omega_m$-plane.

%The dynamical system~\eqref{dynsysexp} admits a number of fixed points on
%$\bar{\bf S}$, all located on the boundary, as follows from the monotonicity
%of $H^2$; they are given in Table~\ref{tab:fixedpointexp} together with
%values for $q$ and, when $\Omega_m\neq 1$, $\gamma_\phi = 2x^2/(1-\Omega_m)$.
%
%\begin{table}[ht!]
%\begin{center}
%\begin{tabular}{|c|c|c|}\hline
% &  Fixed point values & \\ \hline
%$\mathrm{M}_0^\pm$ & $x = \pm 1$, $\Omega_m = \Omega_V = Z =0$ & $q = 2$, $\gamma_\phi=2$ \\ \hline
%$\mathrm{M}_1^\pm$ & $x = \pm 1$, $\Omega_m = \Omega_V = 0$, $Z=1$ & $q = 2$, $\gamma_\phi=2$ \\ \hline
%$\mathrm{PL}_{0,1}$ & $x = \bar{\lambda}/\sqrt{6}$, $\Omega_m = 0$, $Z = 0,1$; $\bar{\lambda} <\sqrt{6}$ & $q = -1 + \bar{\lambda}^2/2$, $\gamma_\phi = \bar{\lambda}^2/3$ \\ \hline
%$\mathrm{EM}_{0,1}$ & $x = \sqrt{3/2}\bar{\lambda}^{-1}$, $\Omega_m = 1 - 3\bar{\lambda}^{-2}$, $Z = 0,1$; $\bar{\lambda} >\sqrt{3}$ & $q =1/2$, $\gamma_\phi=1$ \\ \hline
%$\mathrm{FL}_Z$ & $x = \Omega_V= 0$, $\Omega_m = 1$, $Z=\mathrm{const.}$ & $q = 1/2$ \\ \hline
%%$\mathrm{dS}_{0,1}$ & $x = \Omega_m = 0$, $\Omega_V = 1$, $Z = 0,1$ & $q = -1$ \\ \hline
%\end{tabular}
%\end{center}
%\caption{Table depicting the fixed points in the full three-dimensional state
%space of a model with dust and an exponential scalar field potential $V=
%V_0\exp(-\bar{\lambda}\phi)$.} \label{tab:fixedpointexp}
%\end{table}
%
%
%\begin{table}[ht!]
%\begin{center}
%\begin{tabular}{|c|c|c|c|}\hline
%% &  Fixed point values & &  \\ \hline
%$\mathrm{FL}_*$ & $x = 0$, $\Omega_m = 1$ & $q = (3\gamma_m-2)/2$ & - \\ \hline
%$\mathrm{EM}_*$ & $x = \sqrt{3/2}\gamma_m\lambda^{-1}$, $\Omega_m = 1 - 3\gamma_m\lambda^{-2}$; $\lambda >\sqrt{3\gamma_m}$ & $q =(3\gamma_m-2)/2$
%& $\gamma_\phi=\gamma_m$ \\ \hline
%$\mathrm{M}_*^\pm$ & $x = \pm 1$, $\Omega_m =0$ & $q = 2$ & $\gamma_\phi=2$ \\ \hline
%$\mathrm{PL}_*$ & $x = \lambda/\sqrt{6}$, $\Omega_m = 0$; $\lambda <\sqrt{6}$ & $q = -1 + \lambda^2/2$ & $\gamma_\phi = \lambda^2/3$ \\ \hline
%$\mathrm{dS}_*$ & $x = \Omega_m = 0$; $\lambda=0$ & $q = -1$ & $\gamma_\phi =0$ \\ \hline
%\end{tabular}
%\end{center}
%\caption{Table depicting the fixed points that can occur for a FLRW model
%with a monotonically decreasing scalar field potential with a finite
%$\lambda$ and a perfect fluid with a linear equation of state parameter
%$0<\gamma_m<2$.} \label{tab:fixedpointexp}
%\end{table}
%

The qualitative structure of the solution space of $\bar{\bf S}$ is entirely
determined by the general monotonic features and the structures on the
boundary subsets of ${\bf S}$ and their neighborhoods described earlier. Thus
there are two-parameter sets of orbits in ${\bf S}$ that originate from the
sources $\mathrm{M}_0^+$ and $\mathrm{M}_1^-$ when ${\lambda}<\sqrt{6}$, but
only from $\mathrm{M}_1^-$ when ${\lambda}>\sqrt{6}$, while a one-parameter
set of solutions originates from the transversally hyperbolic line of fixed
points $\mathrm{FL}_{Z}$. Toward the future $\mathrm{PL}_1$ is a sink when
${\lambda} < \sqrt{3}$ (and the future attractor on ${\bf S}$ and ${\bf
S}_\phi$) while $\mathrm{EM}_1$ is future stable when ${\lambda} > \sqrt{3}$
(and the future attractor on ${\bf S}$).
%The present case
%can therefore only mimic $\Lambda$CDM dynamics for quite small
%$\bar{\lambda}$.
As previously shown, it is a general feature that there is a one-parameter
set of solutions that originate from the line $\mathrm{FL}_Z$, forming an
attracting invariant subset, but it is only when ${\lambda}<\sqrt{3}$ that this
invariant subset forms an ``attracting \emph{separatrix} surface,'' which in
the limit ${\lambda}\rightarrow 0$ leads to the previous $\Lambda$CDM
surface.

%Since $x$ is monotonically increasing when $-1<x<0$, it follows that the
%separatrix surface is located at $x\geq 0$; note that the same holds true for
%the attracting separatrix surfaces discussed in the future examples.

It is of historical interest to point out that it was the case with
${\lambda} > \sqrt{6}$ with $\mathrm{EM}_1$ as the future attractor (although
in the two-dimensional projected context) that was the original role model
for concepts such as attractor (the fixed point $\mathrm{EM}_1$ \emph{is} the
formal future attractor) and tracker solutions in scalar field cosmology
[recall that $\gamma_\phi = \gamma_m$ for general $\gamma_m$ (see,
Table~\ref{tab:fixedpoint}); i.e., the effective scalar field equation of
state ``tracks'' the matter equation of state]. As we will see, this is rather
misleading. Moreover, from current observational constraints these models are
no longer viable cosmological models since $q=1/2$ for the future attractor
$\mathrm{EM}_1$. Since models with small ${\lambda}$ are arguably more
relevant for understanding more general and observationally competitive
scalar field models, we focus on this class of models next.
%, while models with sufficiently small $\lambda$ are,
%although, as will be argued below, they are unlikely candidates for solving
%any problems one might have with $\Lambda$CDM cosmology. Nevertheless, from a
%general dynamical systems perspective they are of interest, and we therefore
%focus on this class next.

It is worthwhile to point out that since this does not seem to be generally
known, there are models with a perfect fluid and a scalar field with an
exponential potential that admits simple explicit solutions, as shown by
Uggla \emph{et al.} in~\cite{uggetal95} (where many other solvable scalar
field and modified gravity models can be found as well). One such example is
for dust when ${\lambda}=\sqrt{3/2}$. Using the methods
in~\cite{uggetal95}, we see that $x$ and $\Omega_m$ can be found in terms of
cosmological time $t$:
%
\begin{equation}\label{expsol}
x = \frac{2t^3-c}{4t^3+6t+c}, \qquad
\Omega_{m}=\frac{1-x^2}{1+t^2}=\frac{12t\left(t^3+3t+c\right)}{\left(4t^3+6t+c\right)^2},
\end{equation}
%
which thus explicitly gives the trajectories in the $x-\Omega_m$-plane in
parametrized form, where the constant $c$ characterizes the different
solutions, with $c=0$ yielding the separatrix subset, i.e., the heteroclinic
orbits from $\mathrm{FL}_Z$ to $\mathrm{PL}_1$ in ${\bf S}$ (a heteroclinic
orbit is a solution trajectory that connects two distinct fixed points).
%Although of course not necessary, we will specialize the potential several
%times so that this particular example occurs as a boundary subset when
%illustrating various properties of the solution space of different
%potentials.
The solution space structure is depicted in Fig.~\ref{fig:EXP1} [see also the
reduced state space in Fig.~\ref{fig:SSCC2D}].
%The
%interior solution space for models with ${\lambda}<\sqrt{3}$ for the
%representative value ${\lambda}=\sqrt{3/2}$ is depicted in
%Fig.~\ref{fig:SSEXP}, with focus on the `attractive separatrix surface.'
%In addition
%figures~\ref{fig:SSEXPTop1}, \ref{fig:SSEXPTop2}, \ref{fig:SSEXPTop3} show
%the projected dynamics onto the reduced $x-\Omega_m$ state space for
%${\lambda}=\sqrt{3/2}$, ${\lambda}=2$, and ${\lambda}=3$, respectively,
%representing the three classes of exponential potential models.
%
%\begin{figure}[ht!]
%\begin{center}
%\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=\sqrt{\frac{3}{2}}$]{\label{fig:Z0Dust_1}
%\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_Final.pdf}} \hspace{0.2cm}
%\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=0.2$]{\label{fig:Z0Dust_2}
%\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_v2_Final.pdf}}
%\end{center}
%\caption{Solutions on the invariant $Z=0$ boundary subset of $\bar{\bf S}$.
%In Fig. (a) the space-dot line corresponds to $a=10$, the dotted line to $a=1$, the dash-dotted line to $a=0$,
%the dashed line to $a=-1/2$ and the space-dotted line to $a=-2$. Fig. (b) represent numerical solutions.}
%\label{fig:exp}
%\end{figure}
%
%\begin{figure}[ht!]
%\begin{center}
%\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=\sqrt{\frac{3}{2}}$]{\label{fig:V0Dust}
%\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{V0_Boundary_Final.pdf}} \qquad
%\end{center}
%\caption{Solutions on the invariant $\Omega_V=0$ boundary subset of $\bar{\bf S}$.
%}
%\label{fig:V0}
%\end{figure}
%
\begin{figure}[ht!]
\begin{center}
\subfigure[The scalar field boundary $\Omega_m=0$; ${\lambda}=\sqrt{3/2}$]{\label{fig:M0EXP}
\includegraphics[width=0.44 \textwidth, trim = 0cm 2.5cm 0cm 0cm]{M0_Boundary_EXP.pdf}} \,
\subfigure[The $(x,\Omega_{m},Z)$ state space for dust and ${\lambda}=\sqrt{3/2}$.]{\label{fig:SSEXP}
\includegraphics[width=0.5\textwidth, trim = 0cm 2.5cm 2cm 0cm]{SS_EXP.pdf}} \vspace{-0.5cm}
\end{center}
\caption{State space structures for dust and $V=V_0\exp(-{\lambda}\phi)$ with
${\lambda}=\sqrt{3/2}$.
}
\label{fig:EXP1}
\end{figure}
%
%
%\begin{figure}[ht!]
%\begin{center}
%\subfigure[The $(x,\Omega_{m},Z)$ state space for dust and ${\lambda}=\sqrt{\frac{3}{2}}$.]{\label{fig:SSEXP}
%\includegraphics[width=0.5\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP.pdf}} \quad
%\subfigure[Projection of the $(x,\Omega_{m},Z)$ state space onto $(x,\Omega_{m})$ for dust and ${\lambda}=\sqrt{\frac{3}{2}}$.]{\label{fig:SSEXPTop1}
%\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP_TopView.pdf}}\\
%\subfigure[Projection of the $(x,\Omega_{m},Z)$ state space onto $(x,\Omega_{m})$ for dust and ${\lambda}=2$.]{\label{fig:SSEXPTop2}
%\includegraphics[width=0.44\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP_TopView_Lambda2.pdf}} \qquad
%\subfigure[Projection of the $(x,\Omega_{m},Z)$ state space onto $(x,\Omega_{m})$ for dust and ${\lambda}=3$.]{\label{fig:SSEXPTop3}
%\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXP_TopView_Lambda3.pdf}}
%\vspace{-0.5cm}
%\end{center}
%\caption{State space structure for $\bar{\bf S}$ for dust and
%$V=V_0\exp(-{\lambda}\phi)$.}
%\label{fig:EXP2}
%\end{figure}
%

Let us now consider the interior solutions in Fig.~\ref{fig:SSEXP} with
${\lambda}=\sqrt{3/2}$ in an $\Omega_m-q$ diagram and in several
\emph{cosmographic diagrams} (see Fig.~\ref{fig:EXP2c}), together with the
$\Lambda$CDM model with $\Omega_{m0}=0.3$ and $\Omega_{\Lambda0}=0.7$. To
have a feeling for the evolution in terms of redshift one can compare with
the graphs in Fig.~\ref{fig:LCDMgraphs}, which gives a good estimate for the
history also for the other solutions when expressed in the redshift $z$, as
long as they are close to the $\Lambda$CDM model.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\Omega_{m}-q$ diagram]{\label{fig:Oqec}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_M_EXP.pdf}} 
\subfigure[$H-q$ diagram]{\label{fig:Hqec}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_EXP.pdf}} \\
\subfigure[$H-j$ diagram]{\label{fig:Hjec}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_j_EXP.pdf}}\,\,
\subfigure[$q-j$ diagram]{\label{fig:qjec}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_j_EXP.pdf}}
\end{center}
\vspace{-0.5cm}
\caption{Diagrams for the observable $\Omega_m$ and the cosmographic parameters $H$, $q$, $j$ (the Hubble,
deceleration, and jerk parameters) for dust and a potential $V = V_0 \exp(-\sqrt{3/2}\,\phi)$
for the solutions (with the same notation) in Fig.~\ref{fig:SSEXP}.
The thick line represents $\Lambda$CDM cosmology with $\Omega_{m0}=0.3$, $\Omega_{\Lambda 0}=0.7$.
The horizontal lines with $\Omega_{m}=0.3$, $H/H_0=1$, and $q=-0.55$ represent the present situation, 
i.e., $z=0$.}
\label{fig:EXP2c}
\end{figure}
%
%
%\begin{figure}[ht!]
%\begin{center}
%\includegraphics[width=0.4\textwidth]{h_q_j_EXP.pdf} \vspace{-1.5cm}
%\end{center}
%\caption{$H-q-j$; $\lambda_0=\sqrt{\frac{3}{2}}$.}
%\label{fig:EXPHqj}
%\end{figure}
%
As can be seen, these models are deviating rather significantly from the
$\Lambda$CDM model, and therefore a considerably smaller value of ${\lambda}$
is needed to obtain an evolution that is close to $\Lambda$CDM dynamics, even
when $\Omega_m$ is close to 1 initially. Furthermore, in a certain sense,
discussed next, the present class of models is as special as the scalar
field models with a cosmological constant, and they are therefore unlikely
candidates for solving any of the issues one might have with $\Lambda$CDM
cosmology.
%, even though models with very small
%$\bar{\lambda}$, and with $\Omega_m$ close to one initially, are going to
%have an evolutionary history very similar to that of the $\Lambda$CDM models.


%-----------------------------------------------------------------
\subsection{Symmetries of frozen $\lambda$CDM dynamics}
%-----------------------------------------------------------------

It is only the cases of a constant or an exponential potential that
effectively reduce the problem from three to two dimensions by decoupling the
equation for the scalar field variable. The underlying reason for this is
that these models admit scaling symmetries.

In the case of an exponential potential, a translation of the scalar field
leads to a scaling of the potential, which in combination with an appropriate
scaling of $t$ yields a scaling of $\rho_\phi$; furthermore, a scaling of the
spatial coordinates leads to a scaling of the scale factor $a$ and thereby
also of $\rho_m$ in the case of a linear equation of state with $\gamma_m\neq
0$, which results in a scaling symmetry of Einstein's field equations. This
in turn gives rise to a one-parameter set of equivalent solutions, being the
reason for the decoupling property, and it is also the underlying reason for
why there exist ``scaling solutions,'' which more appropriately should be
referred to as homothetic self-similar solutions since the spacetime geometry
of these solutions admits a homothetic Killing vector field. For models with
both a perfect fluid with a linear equation of state and a scalar field with
an exponential potential, a scaling solution can only exist if $\rho_m$ and
$\rho_\phi$ behave in the same manner, which requires $\gamma_m=\gamma_\phi$.
The existence of a fixed point $\mathrm{EM}$ on the reduced state space
$x-\Omega_m$ therefore necessarily depends on this feature, but it turns out
to not always be possible since $\mathrm{EM}$ only exists when ${\lambda} >
\sqrt{3\gamma_m}$; furthermore, the \emph{global} attractive property of
$\mathrm{EM}$, which does \emph{not} follow from a fixed point analysis, follows
from a monotone function that can be \emph{derived} from another kind of
symmetry of the field equations that is associated with coordinate scalings
(this is a general mechanism, discussed for anisotropic models and perfect
fluids in~\cite{heiugg10}). For a constant potential scalings of the
coordinates instead bring a solution with a given $V=\Lambda$ to a solution
with a different $\Lambda$, and because $\Lambda$ carries dimension
(length$^{-2}$), there are no solutions admitting a homothetic symmetry when
$\Lambda \neq 0$; i.e., this scaling symmetry is not quite the same as that
in the exponential case.

The existence of symmetries is what makes these models tractable, and therefore
also popular. However, the symmetries endow these models with special
properties, which suggests that perhaps it is not a good idea to consider
these models as role models, since other models do not admit similar
symmetries and properties, or at least that some careful considerations may
be needed. Moreover, originally, before one realized that the deceleration
parameter $q$ is evolving, it was thought that scale invariance and thereby a
constant $q$ was quite desirable, since this would solve fine-tuning problems
such as coincidence and energy scale problems. On the other hand, the
symmetries also make the models quite inflexible. Due to the changing
evolutionary history of $q$, where $q$ is currently negative, the scaling
symmetries must be broken, while still leading to models that in some way at
least alleviate various fine-tuning problems; for these reasons, and from a
phenomenological perspective, scale-invariant models are no longer the favorite
candidates as dark energy models.

%Next we consider the simplest possible generalization by breaking the scaling
%symmetry by having a potential that consists of an exponential and a constant
%term. Although these models are unlikely to serve as useful dynamical dark
%energy models, they do exhibit some features that are characteristic of more
%general classes of models, and it is for this reason we consider them next.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Dynamical $\lambda$CDM dynamics}\label{sec:dynlambda}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-----------------------------------------------------------------
\subsection{Inverse power-law potentials}
%-----------------------------------------------------------------

In the influential papers~\cite{peerat88,ratpee88} Peebles and Ratra argued
that dark energy can be phenomenologically modeled by a minimally coupled
scalar field with a potential with a shallow tail, which results in the dark field 
energy density decreasing more slowly than the matter energy density to
its ``natural'' value -- zero. In particular they considered an inverse
power-law potential,
%
\begin{equation}
V=\frac{V_0}{\phi^{\alpha}}, \qquad \phi >0, \quad \alpha > 0 \quad \Rightarrow \quad \lambda = \frac{\alpha}{\phi}.
\end{equation}
%
In addition they assumed a matter dominated flat FLRW universe after
inflation, which resulted in a set of approximate flat FLRW equations for
which they found a particular solution, which was shown to be linearly stable
within the approximate context. The associated solution of the exact FLRW
equations has subsequently been referred to as an attractor or tracker
solution~\cite{steetal99}, and has resulted in numerous papers.
In~\cite{tam14} a thorough local dynamical systems analysis of the present
models was performed, based on previous work,
e.g.,~\cite{ngetal01,ure12,gon14} (for additional references see,
e.g.,~\cite{tsu13,tam14}). Here we will use a slight variation of this
approach and use the following bounded scalar field variable:
%
\begin{equation}
Z=\frac{1}{1+\lambda} = \frac{\phi}{\phi + \alpha} \quad \Rightarrow \quad \lambda=\frac{1-Z}{Z},
\quad \phi = \alpha\left(\frac{Z}{1-Z}\right),
\end{equation}
%
where $Z$ thereby is monotonically increasing in $\phi$ (in~\cite{tam14}
$1-Z$ was used as the scalar field variable).

The present case has an unbounded $\lambda$, and therefore we use the
dynamical system~\eqref{dynsysZ}, where we choose $g(Z) = Z$, and thus
$g_{\lambda_-}=1$ and
%
\begin{equation}
\frac{d{\tau}}{d\bar{\tau}} = Z, \qquad \frac{dt}{d\bar{\tau}} = H^{-1}Z.
\end{equation}
%
and the following dynamical system for the state vector $(x,\Omega_m,Z)$:
%
\begin{subequations}
\begin{align}
\frac{dx}{d\bar{\tau}} &=-(2-q)xZ+\sqrt{\frac{3}{2}}\left(1-Z\right)\left(1-x^2-\Omega_m\right), \\
\frac{d\Omega_{m}}{d\bar{\tau}} &= 3\left[2x^2-(1-\Omega_{m})\right]\Omega_m Z, \\
\frac{dZ}{d\bar{\tau}} &=\frac{\sqrt{6}}{\alpha}\left(1-Z\right)^2 Z x,
 \end{align}
\end{subequations}
%
where $2-q=3(1-x^2) - \frac32 \Omega_m$.

Due to the discussion about global dynamics in Sec.~\ref{sec:dynsys}, it
follows that the conclusions obtained from the local fixed point results
in~\cite{tam14} correspond to global asymptotic features. Hence a
two-parameter set of solutions originates from $\mathrm{M}^-_1$, and a
one-parameter set from $\mathrm{FL}_Z$, where the tracker solution originates
from $\mathrm{FL}_0$. Toward the future all orbits on ${\bf S}$ and ${\bf
S}_\phi$ end at the global future attractor $\mathrm{dS}_1$.\footnote{Since
$\lim_{\phi\rightarrow +\infty}V(\phi)= 0$ it follows that the future
asymptotic state $\mathrm{dS}_1$ now corresponds to the Minkowski spacetime
instead of the de Sitter spacetime. However, because the matter energy
density goes to zero faster than the scalar field energy density it still
follows that $q\rightarrow -1$; the asymptotic evolution thus resembles that
of the de Sitter spacetime, since asymptotically the spacetime approaches
that of a Minkowski spacetime in a (local) foliation where $q\rightarrow -1$.
A similar statement holds for the models in the next subsection.} The
solution structure for $\alpha=6$, and especially that of the subset that
originates from $\mathrm{FL}_Z$, is depicted in Fig.~\ref{fig:RP_alpha6}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[The $(x,\Omega_{m},Z)$ state space.]{\label{fig:SSRP_alpha6}
\includegraphics[width=0.45\textwidth, trim = 1cm 2cm 2cm 0cm]{SS_RP_6.pdf}}
\subfigure[Projection of $(x,\Omega_{m},Z)$ state space onto  $(x,\Omega_{m})$.]{\label{fig:SSRPTop_alpha6}
\includegraphics[width=0.45\textwidth, trim = 0.5cm 2cm 1cm 0cm]{SS_RP6_TopView.pdf}}
\end{center}
\vspace{-0.5cm}
\caption{Solution structure for the state space for models with dust and a scalar field
with inverse power-law potential $V=V_0\phi^{-\alpha}$ with $\alpha=6$.}
\label{fig:RP_alpha6}
\end{figure}
%

It is of some interest to complement the dynamical systems picture with a
heuristic description where we regard the scalar field as a particle moving
in a potential, while losing energy due to the friction force
$-3H\dot{\phi}$. Due to the existence of a potential wall that is so steep
that $\lambda \rightarrow \infty$ when $\phi\rightarrow 0$, all ``scalar field
particles'' (solutions) either come from $\phi \rightarrow -\infty$ and bounce
against the potential wall (the two-parameter set of solutions from
$\mathrm{M}^-_1$) or from being initially still at some $\phi$ (no initial
kinetic scalar field energy, i.e., $x=0$) and then rolling down the potential
in the positive $\phi$-direction; the tracker solution corresponds to the
solution that starts by being initially still with an infinite scalar field
energy associated with the limit $\phi \rightarrow 0$ (note that this result
also follows from the asymptotic explicit results concerning this solution
given originally in~\cite{peerat88,ratpee88}).

All solutions therefore eventually move in an increasingly shallow potential,
slowed by ``friction,'' and hence with decreasing scalar field energy,
although the matter energy decreases faster, thus leading to a future
asymptotic (quasi) de Sitter state at $\phi \rightarrow +\infty$ (i.e.,
$\mathrm{dS}_1$).

Loosely speaking, the tracker solution is connected with a bifurcation
associated with going from a potential with finite $\lambda$ to infinite
$\lambda$, transforming the attracting focus of the exponential case with
large $\lambda$ to a center. The orbits on $Z=0$ are characterized by
$\Omega_m = \mathrm{const}$. They therefore correspond to the periodic
orbits found in the monomial case at late times, given in different variables
in~\cite{alhetal15}, but in the present case they are interrupted by being
cut in half by the line of fixed points, $\mathrm{M}^x_0$, which is a formal
consequence of the new time variable. The physical reason for this feature is
that in contrast to the monomial case there are no scalar field oscillations
at late times for inverse power-law potentials.

As can be seen from Fig.~\ref{fig:RP_alpha6}, the tracker solution forms part
of the boundary of the $\mathrm{FL}_Z$ saddle subset, which therefore is not a
separatrix surface. Because of its initial linear stability property in the Peebles-Ratra 
formulation~\cite{peerat88,ratpee88}, which was seen as alleviating the
coincidence problem, nearby solutions are attracted to it.
However, there exists an open set of solutions that are further away from it
that are not, even if $\Omega_m$ is close to 1 initially (those with fairly
large initial $Z$). This is to be expected since the Peebles-Ratra
stability analysis is only local. It should also be pointed out that there is
no {\it a priori} reason for believing that the local attracting property
holds globally in time.\footnote{See, e.g.~\cite{ugg13} for an example of an
attracting solution that is locally but not globally attracting.} However,
all solutions near the tracker solution, indeed all solutions that are near
or on the subset that originates from $\mathrm{FL}_Z$, are eventually
shadowing the center manifold of the fixed point $\mathrm{dS}_1$ on ${\bf
S}_\phi$\footnote{Incidentally, it is the center manifold solution on ${\bf
S}_\phi$ that the slow-roll approximation, given in, e.g.,~\cite{wei08},
describes approximately; see~\cite{alhetal15,alhugg15} for a similar
situation at early times instead of late times.} and it is the center
manifold solution on ${\bf S}_\phi$ that is an ``attractor solution'' rather than the
tracker solution, since any solution close to this solution is attracted to
it because it is a center manifold, but note that it is $\mathrm{dS}_1$
that is the formal future attractor.

As seen from Fig.~\ref{fig:RP_alpha6}, the tracker solution and the solutions
that are near to it on the subset that originates from $\mathrm{FL}_Z$ are
those solutions on this subset that diverge from $\Lambda$CDM dynamics the
most (in this sense they are the observationally worst solutions in this
subset of solutions, which is also seen by plotting the solutions in diagrams
for the observables; see Fig.~\ref{fig:Cosm_RP_alpha6}) 
and quite substantially when $\alpha=6$. Let us
therefore consider the dynamics for models with lower $\alpha$.

%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\Omega_{m}-q$ diagram.]{\label{fig:q_M_RP_alpha6}
\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{q_M_RP6.pdf}}
\subfigure[$h-q$ diagram.]{\label{fig:h_q_RP_alpha6}
\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_RP6.pdf}}
\subfigure[$h-j$ diagram.]{\label{fig:h_j_RP_alpha6}
\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_j_RP6.pdf}}
\subfigure[$q-j$ diagram.]{\label{fig:j_q_RP_alpha6}
\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{q_j_RP6.pdf}}
%\subfigure[$h-q-j$.]{\label{fig:h_q_j_RP_alpha6}
%\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_j_RP6.pdf}}
\end{center}
\caption{Observables for dust and an inverse power-law potential $V=V_0\phi^{-\alpha}$ with $\alpha=6$. 
The thick solid line represents $\Lambda$CDM cosmology; the horizontal lines represent now, i.e., zero redshift; 
and the other lines correspond to the solutions in Fig.~\ref{fig:RP_alpha6}, where the grey line corresponds 
to the tracker solution. 
}
\label{fig:Cosm_RP_alpha6}
\end{figure}
%

The solution structure for $\alpha=1/2$, and especially that of the subset
that originates from $\mathrm{FL}_Z$, is depicted in
Fig.~\ref{fig:RP_alpha1_2}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$(x,\Omega_{m},Z)$ state space.]{\label{fig:SSRP_alpha1_2}
\includegraphics[width=0.45\textwidth, trim = 1cm 2cm 2cm 0cm]{SS_RP_1_2.pdf}}
\subfigure[Projection  of $(x,\Omega_{m},Z)$ state space onto  $(x,\Omega_{m})$.]{\label{fig:SSRPTop_alpha1_2}
\includegraphics[width=0.45\textwidth, trim = 0.5cm 2cm 1cm 0cm]{SS_RP1_2_TopView.pdf}}
\end{center}\vspace{-0.5cm}
\caption{The solution space of dust and an inverse power-law potential $V=V_0\phi^{-\alpha}$ with $\alpha=1/2$.
}
\label{fig:RP_alpha1_2}
\end{figure}
%
The topological structure is the same as for $\alpha=6$, but the subset that
originates from $\mathrm{FL}_Z$ is now located closer to the $x=0$-plane and
as a consequence the dynamics of solutions with $\Omega_m$ close to 1
initially more closely resembles that of the $\Lambda$CDM models. These
results are therefore in line with the conclusion in~\cite{tsu13}, that it is only
for quite small $\alpha$ that the tracker solution is observationally viable,
which restricts the observational relevance of these models. The deviations
from $\Lambda$CDM cosmology can be explicitly seen in the diagrams for the
observables in Fig.~\ref{fig:Cosm_RP_alpha1_2}.
%As for
%$\alpha=6$, the tracker solution and the nearby solutions are those with the
%worst evolution in this set of solutions.
%, in the sense that they all deviate
%substantially from $\Lambda$CDM evolution, in spite of a quite small value of
%$\alpha$.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\Omega_{m}-q$ diagram]{\label{fig:q_M_RP_alpha1_2}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_M_RP1_2.pdf}}
\subfigure[$h-q$ diagram]{\label{fig:h_q_RP_alpha1_2}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_RP1_2.pdf}}
\subfigure[$h-j$ diagram]{\label{fig:h_j_RP_alpha1_2}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_j_RP1_2.pdf}}\,\,
\subfigure[$q-j$ diagram]{\label{fig:j_q_RP_alpha1_2}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_j_RP1_2.pdf}}
%\subfigure[$h-q-j$.]{\label{fig:h_q_j_RP_alpha1_2}
%\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_j_RP1_2.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{Observables for dust and an inverse power-law potential $V=V_0\phi^{-\alpha}$ with $\alpha=1/2$. 
The thick solid line represents $\Lambda$CDM cosmology; the horizontal lines represent now, i.e., zero redshift;
and the other lines correspond to the solutions in Fig.~\ref{fig:RP_alpha1_2}, where the grey line corresponds 
to the tracker solution.}
\label{fig:Cosm_RP_alpha1_2}
\end{figure}
%

Note that Figs.~\ref{fig:RP_alpha6} and~\ref{fig:RP_alpha1_2} also show that it is not only the
tracker solution that attracts nearby solutions, but the whole
$\mathrm{FL}_Z$ saddle subset. The underlying reasons for this are that (i)
the whole line of fixed points $\mathrm{FL}_Z$ acts as a line of saddle
points, and (ii) all solutions in this subset are attracted to the
center manifold of $\mathrm{dS}_1$ at late times before asymptotically
approaching the global future attractor $\mathrm{dS}_1$.

Finally, it is worth noting that the underlying reason for why the early
behavior of the tracker solution exhibits a ``tracking property'' is because
the present case can be viewed as a bifurcation where
$\mathrm{EM}_0$ merges with $\mathrm{FL}_0$ in the limit $\lambda \rightarrow
\infty$. The property that $\gamma_\phi=\gamma_m=1$ for $\mathrm{EM}$,
irrespective of $\lambda$ (as long as $\mathrm{EM}$ exists until it merges
with $\mathrm{PL}$), can subsequently be exploited by comparing the reduced
two-dimensional dynamics of a sequence of reduced exponential scalar field
state spaces, as done in~\cite{tam14,ngetal01,ure12}, to yield an approximate
heuristic description of the behavior of the tracking solution. Nevertheless,
it should be clear that the dynamics of an inverse power-law potential model
is globally very different from that of a given exponential potential, which
is what can be expected from a model that does not exhibit scaling
symmetries.

%, although lately its observational relevance has been
%questioned~\cite{tsu13}. That this is the case should really not come as a
%surprise if one considers the situation from a dynamical systems perspective.

The observational viability difficulties of the dust and inverse power-law
potential models, including those of tracker solutions, arise because
$\lambda\rightarrow \infty$ when $\phi\rightarrow \infty$, since a large
$\lambda$ makes the dynamics deviate from $\Lambda$CDM dynamics. This
motivates us to consider other types of scalar field potentials for which
$\lambda$ is globally and asymptotically regularized. For this reason, we
therefore phenomenologically regularize $\lambda$ of the inverse power-law
potential and use~\eqref{Vfromlambda} to derive $V$, in arguably the
simplest possible way, so that $\lambda$ becomes finite and one obtains
continuous deformations of $\Lambda$CDM cosmology.

%-----------------------------------------------------------------
\subsection{$\lambda$-regularization of inverse power-law potentials}
%-----------------------------------------------------------------

To obtain a type of potential that has a regular $\lambda$ and still behaves
like an inverse power-law potential for large $\phi$, we generalize the
inverse power-law potential by choosing
%
\begin{equation}
\lambda = \frac{\alpha}{\sqrt{C^2 + \phi^2}} = \frac{\lambda_\mathrm{max}}{\sqrt{1 + (\phi/C)^2}}, \quad C>0.
\end{equation}
%
Here $\lambda_\mathrm{max} = \alpha/C$ is the maximum value of
$\lambda(\phi)$ while $C$ describes the peak width of the present
$\lambda(\phi)$. Using~\eqref{Vfromlambda}, this ``regularized'' inverse
power-law $\lambda(\phi)$ leads to the potential
%
\begin{equation}
V = V_0(\phi +\sqrt{C^2+\phi^2})^{-\alpha}.
\end{equation}
%
It follows that $\lim_{\phi \rightarrow \pm\infty}(\phi^{\pm\alpha} V) =
\mathrm{const}$. Note that this class of potentials contains the inverse
power-law potential by setting $C=0$ (restricting dynamics to $\phi>0$) and a
cosmological constant by setting $\alpha=0$. As an example, we give the
potential and $\lambda$ for $\alpha=6$ and $C=10$, depicted in
Fig.~\ref{fig:Va10}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[Potential diagram]{\label{fig:Va6C10}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{V_phi_Alpha6_C10.pdf}}
\subfigure[$\lambda$ diagram]{\label{fig:La6C10}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{Lambda_phi_Alpha6_C10.pdf}}
%\subfigure[$h-q-j$.]{\label{fig:h_q_j_RP_alpha1_2}
%\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_j_RP1_2.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{The potential $V = V_0(\phi +\sqrt{C^2+\phi^2})^{-\alpha}$ and associated $\lambda= \alpha/\sqrt{C^2 + \phi^2}$
for $\alpha=6$ and $C=10$.}
\label{fig:Va10}
\end{figure}
%

The above is (so far) a purely phenomenologically motivated potential, but so
was the original motivation for the inverse power-law
potential~\cite{peerat88,ratpee88}. The point here is to give a simple
specific example of continuous $\Lambda$CDM scalar field deformations,
illustrating features that \emph{any} observationally viable model arising from some
fundamental theory must have.

To continue we choose
%
\begin{equation}
Z= \frac{z}{1+z}, \qquad z = \phi + \sqrt{C^2 + \phi^2}, \qquad \phi = \frac{z^2-C^2}{2z} 
= \frac{Z^2-C^2(1-Z)^2}{2Z(1-Z)}
\end{equation}
%
as the scalar field variable (here $z$ should not be confused with the redshift). 
This leads to
%
\begin{equation}
\lambda = \frac{2\alpha Z(1-Z)}{Z^2+C^2(1-Z)^2}
\end{equation}
%
and the regular dynamical system
%
\begin{subequations}\label{dynsysreg}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}\lambda(Z)(1 - x^2 - \Omega_m),\label{xeqreg}\\
\Omega_m^\prime &= 3\left[2x^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\label{Omeqreg}\\
Z^\prime &= \sqrt{6} \alpha^{-1}\lambda Z(1-Z)\,x = \sqrt{6}\frac{2Z^2(1-Z)^2}{Z^2+C^2(1-Z)^2}\,x\label{Zeqreg}
\end{align}
\end{subequations}
%
%where
%
%\begin{equation}
%\frac{dZ}{d\phi} = \alpha^{-1}\lambda Z(1-Z) = \frac{2Z^2(1-Z)^2}{Z^2+C^2(1-Z)^2}
%\end{equation}
%
where $2-q=3(1-x^2) - \frac32 \Omega_m$. The qualitative properties of this
system follow from the general discussion in Sec.~\ref{sec:dynsys}. We
give an explicit representation of the solution space for the case $C=10$ and
$\alpha=6$ in Fig.~\ref{fig:sol_a6_C10}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$(x,\Omega_{m},Z)$ state space.]{\label{fig:SS_RegAlpha6}
\includegraphics[width=0.45\textwidth, trim = 1cm 2.2cm 2cm 0cm]{GenINVPL_Alpha6_C10.pdf}}
\subfigure[Projection  of $(x,\Omega_{m},Z)$ state space onto  $(x,\Omega_{m})$.]{\label{fig:SS_RegAlpha6_Top}
\includegraphics[width=0.45\textwidth, trim = 0.5cm 2cm 1cm 0cm]{GenINVPL_Alpha6_C10_Topview.pdf}}
%\subfigure[$h-q-j$.]{\label{fig:h_q_j_RP_alpha1_2}
%\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_j_RP1_2.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{Depiction of the solution space of flat FLRW models with dust and a scalar field with a potential
$V = V_0(\phi +\sqrt{C^2+\phi^2})^{-\alpha}$ and associated $\lambda= \alpha/\sqrt{C^2 + \phi^2}$
for $\alpha=6$ and $C=10$.}
\label{fig:sol_a6_C10}
\end{figure}
%
The associated diagrams for physical observables are given in
Fig.~\ref{fig:reglamdacosmkin}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\Omega_{m}-q$ diagram]{\label{fig:q_M_RegAlpha6}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_M_GenInv_Alpha6_C10.pdf}}
\subfigure[$h-q$ diagram]{\label{fig:h_q_RegAlpha6}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_GenInv_Alpha6_C10.pdf}}
\subfigure[$h-j$ diagram]{\label{fig:h_j_RegAlpha6}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{h_j_GenInv_Alpha6_C10.pdf}}\,\,
\subfigure[$q-j$ diagram]{\label{fig:j_q_RegAlpha6}
\includegraphics[width=0.48\textwidth, trim = 0cm 0cm 0cm 0cm]{q_j_GenINV_Alpha6_C10.pdf}}
%\subfigure[$h-q-j$.]{\label{fig:h_q_j_RP_alpha1_2}
%\includegraphics[width=0.45\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_j_RP1_2.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{Observables for flat FLRW models with dust and a scalar field with the
potential $V = V_0(\phi +\sqrt{C^2+\phi^2})^{-\alpha}$ and associated $\lambda= \alpha/\sqrt{C^2 + \phi^2}$
for $\alpha=6$ and $C=10$. The thick solid line represents $\Lambda$CDM cosmology; 
the horizontal lines represent now, i.e., zero redshift;
and the other lines correspond to the solutions in Fig.~\ref{fig:sol_a6_C10}.}
\label{fig:reglamdacosmkin}
\end{figure}
%
As can be seen the $\lambda$-regularized models give much better
results than the corresponding models with an inverse power-law potential with 
the same $\alpha$, when $C$ is chosen so that it is of the same order of 
magnitude as $\alpha$, and even better if $C$ is larger (cf. 
Figs.~\ref{fig:Cosm_RP_alpha6} and~\ref{fig:reglamdacosmkin}, which both have 
$\alpha=6$). The results basically linearly improve with increasing $C$ for a
fixed $\alpha$, which should not come as a surprise since
$\lambda_\mathrm{max} = \alpha/C$; for $\alpha =6$ and $C=100$, and therefore
$\lambda_\mathrm{max}=0.06$, the dynamics for models with $\Omega_m$ close to
1 initially is almost indistinguishable from $\Lambda$CDM dynamics. Toward
smaller $C$, the models increasingly yield similar results as the inverse
power-law potential, which again is to be expected.


%-----------------------------------------------------------------
\subsection{Observational $\lambda$ conditions for continuous $\Lambda$CDM deformation}
%-----------------------------------------------------------------

The structure of the flat FLRW models with dust and inverse power-law
potential strongly suggests that it does not suffice to just have a
monotonically decreasing potential with a shallow potential tail, since even
if $\lambda\rightarrow 0$ when $\phi \rightarrow +\infty$ a large
$\lambda(\phi)$ for some $\phi$
%for models with $\Omega_m$ close to one,
generally results in considerable deviations from $\Lambda$CDM cosmology.
There are significant differences in the solution structure if $\lambda$
is globally and asymptotically bounded or not.

This suggests that observational viability implies \emph{global and
asymptotic boundedness conditions} on $\lambda$, 
so that the whole subset that originates from $\mathrm{FL}_Z$ into
${\bf S}$ only deviates moderately from the $x=0$ $\Lambda$CDM surface. It
should be pointed out that this not only produces an observationally viable
set of models, but also alleviates the coincidence and possibly other
fine-tuning problems since \emph{all} orbits with initial values of
$\Omega_m$ close to 1 behave similarly because the invariant subset
that originates from $\mathrm{FL}_Z$ then constitutes an attracting separatrix
surface, eliminating the need for any measure.\footnote{Because open sets of 
solutions do not follow the tracker solution closely,
it was argued in~\cite{tam14} that solutions that were close
to the tracker solution were in some sense generic due to the center manifold
structures. However, such structures can be eliminated or changed by a change
of variables, while the fact that there are open sets of solutions that
behave quite differently than the tracker solution in terms of physical
observables cannot. It is therefore a question of measures if trackerlike
behavior is generic or not, and a consensus about measures seems rather
unlikely. On the other hand, if center manifold structures reflect behavior
in physical observables, they might very well play an important role in
cosmology.}



\subsection*{Acknowledgments}
A. A. is supported by the projects CERN/FP/123609/2011, EXCL/MAT-GEO/0222/2012,
and CAMGSD, Instituto Superior T{\'e}cnico by FCT/Portugal through
UID/MAT/04459/2013, and the FCT Grant No. SFRH/BPD/85194/2012. Furthermore, A. A.
also thanks the Department of Engineering and Physics at Karlstad University,
Sweden, for kind hospitality.


\begin{appendix}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{A scalar field and two fluids}\label{app:twofluid}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this work we have considered a matter content that consists of a fluid
without pressure. However, the early Universe after inflation is radiation
dominated and therefore somewhat more sophisticated models have both
radiation and dust as matter content. In a dynamical systems context this
situation can be treated as follows. Consider two noninteracting fluids with
energy densities $\rho_1$ and $\rho_2$, respectively. It is then convenient
to introduce
%
\begin{equation}\label{chi}
\rho_m = \rho_1 + \rho_2,\qquad \chi = \frac{\rho_1}{\rho_m},
\end{equation}
%
as the two perfect fluid variables instead of $\rho_1$ and
$\rho_2$.\footnote{For alternative ways of handling several fluids, see,
e.g.,~\cite{colwai92,heietal05}.} Assume further that the fluids obey linear
equations of state characterized by the constant equation of state parameters
$\gamma_1$ and $\gamma_2$, where we without loss of generality set
$\gamma_2>\gamma_1$ (if $\gamma_2=\gamma_1$ we treat the problem as a single
fluid). The system~\eqref{dynsys0} holds for fluids with general barotropic
equations of state, which includes the present case. However, having two
fluids instead of one with linear equations of state makes it convenient to
augment the system~\eqref{dynsys0} with the additional variable $\chi$, which
leads to
%
\begin{subequations}\label{dynsys2}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}\lambda(Z)\,\Omega_V,\label{xeq2}\\
\Omega_m^\prime &= 3\left[2x^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\label{Omeq2}\\
Z^\prime &= \sqrt{6}\frac{dZ}{d\phi}\,x,\label{Zeq2}\\
\chi^\prime &= 3(\gamma_2-\gamma_1)\chi(1-\chi),\label{chi2}
\end{align}
\end{subequations}
%
where, as before, $\Omega_V = 1 - x^2 - \Omega_m$ and $\Omega_m=
\rho_m/3H^2$. The deceleration parameter is given by
%
\begin{equation}\label{qchi}
q = - 1 + 3x^2 + \frac32\gamma_m\Omega_m = - 1 + 3x^2 + \frac32\left(\gamma_1\chi + \gamma_2(1-\chi)\right)\Omega_m.
\end{equation}
%
since $\gamma_m = \frac{\rho_1 + \rho_2 + p_1+p_2}{\rho_1+\rho_2} =
\gamma_1\chi + \gamma_2(1-\chi)$. Due to the similarity with the single fluid
case we choose to include the same boundaries, but with the addition of
$\chi$ we also consider $\chi \in [0,1]$, where $\chi=0$ ($\chi=1$)
corresponds to the single fluid case for $\rho_2$ ($\rho_1$).

Since $\gamma_2>\gamma_1$ it follows that $\chi$ is strictly monotonically
increasing when $0<\chi<1$. Furthermore, $\chi$ can be solved in terms of the
scale factor $a$ by using the conservation equation $d\rho_i/d\ln a = -
3\gamma_i\rho_i$, which yields $\rho_i \propto a^{-3\gamma_i}$ and therefore
%
\begin{equation}\label{chia}
\chi = \frac{1}{1 + k(a/a_0)^{\gamma_1-\gamma_2}} = \frac{1}{1 + k\exp[(\gamma_1-\gamma_2)\tau]},
\end{equation}
%
where $k= \rho_{20}/\rho_{10} = \Omega_{20}/\Omega_{10}$. It follows that
$\chi \rightarrow 0$ ($\chi \rightarrow 1$) when $a\rightarrow 0$, i.e.,
$\tau \rightarrow - \infty$ ($a\rightarrow \infty$, i.e. $\tau \rightarrow +
\infty$). Thus $\chi =0$ in the asymptotic past and $\chi=1$ in the
asymptotic future, which tells us that the fluid with the most soft (stiff)
equation of state dominates asymptotically over the other fluid at late
(early) times. As a consequence of the above features it follows that
asymptotic future (past) dynamics reside on the $\chi =1$ ($\chi =0$) subset.

Because it is possible to solve the conservation equations for each of
the two fluids in terms of the scale factor, this gives rise to a constant of
motion, which also reflects that the problem is really a
three-dimensional one. To obtain this conserved quantity we first note that
%
\begin{equation}
\Omega_1 = \chi\Omega_m, \qquad \Omega_2 = (1-\chi)\Omega_m.
\end{equation}
%
Then, since $\rho_1 \propto a^{-3\gamma_1}$ and $\rho_2 \propto
a^{-3\gamma_2}$, it follows that
%
\begin{equation}
\mathrm{const.} = \frac{\rho_1^{\gamma_2}}{\rho_2^{\gamma_1}} \propto
\frac{\Omega_1^{\gamma_2}}{\Omega_2^{\gamma_1}}(H^2)^{\gamma_2-\gamma_1}
= \chi^{\gamma_2}(1-\chi)^{-\gamma_1}\Omega_m^{\gamma_2-\gamma_1}(H^2)^{\gamma_2-\gamma_1},
\end{equation}
%
where $H^2$ is given in terms of the state space variables in Eq.~\eqref{H2}.

%Thus, since $\gamma_2 > \gamma_1$, if $\chi \rightarrow 0$ or $\Omega_m
%\rightarrow 0$ then $T \rightarrow 0$, and vice versa, or if $\chi
%\rightarrow 1$ then $T\rightarrow 1$, and vice versa.
%This, of course,
%assumes that $\rho_1\rho_2>0$ and say nothing about the special cases with
%only one or no fluids.

The matter subset $\Omega_m=1$, i.e., $x=0$, $\Omega_V=0$, now corresponds to
the two-dimensional system
%
\begin{subequations}\label{dynsysOm1}
\begin{align}
Z^\prime &= 0,\label{ZeqOm1}\\
\chi^\prime &= 3(\gamma_2-\gamma_1)\chi(1-\chi),\label{chiOm1}
\end{align}
\end{subequations}
%
which has the solution $Z=\mathrm{const}.$, $\chi = 1/(1 +
k\exp[(\gamma_1-\gamma_2)\tau])$. The analysis of the single fluid case is
directly translatable to the present context, and the line of fixed points
$\mathrm{FL}_Z$ is now replaced with the above two-dimensional subset. It
follows that the inclusion of radiation changes little of the previous
discussion as regards qualitative dynamical features.

It should perhaps be pointed out that instead of $\chi$ and $\Omega_m $ it is
possible to use $\Omega_1=\rho_1/3H^2$ and $\Omega_2=\rho_2/3H^2$  as
variables, where
%
\begin{equation}
\frac{d\Omega_i}{d\tau} = [2q - (3\gamma_i-2)]\Omega_i,\qquad q = - 1 + 3x^2 +
\frac32(\gamma_1\Omega_1 + \gamma_2\Omega_2),
\end{equation}
%
and where $\Omega_V = 1 - x^2 - \Omega_1 - \Omega_2$.

\end{appendix}


\begin{thebibliography}{99}

\bibitem{plaI15} R.~Adam {\it et al}. (Planck Collaboration), 
\newblock Planck Collaboration I, Planck 2015 results. I. Overview of products and results,
\newblock arXiv:1502.01582.

\bibitem{plaXIII15} P.~A.~R.~Ade {\it et al}. (Planck Collaboration) 
\newblock Planck 2015 results. XIII. Cosmological parameters,
\newblock  arXiv:1502.01589.

\bibitem{peerat88} P.~J.~E.~Peebles and B.~Ratra.
\newblock Cosmology with a time-variable cosmological ``constant''.
\newblock Astrophys.~J. {\bf 325}, L17 (1988).%L17-L20
\newblock DOI: 10.1086/185100

\bibitem{ratpee88} B.~Ratra and P.~J.~E.~Peebles.
\newblock Cosmological consequences of a rolling homogeneous scalar field.
\newblock Phys.~Rev. D{\bf 37}, 3406 (1988).
\newblock DOI: 10.1103/PhysRevD.37.3406

\bibitem{steetal99} P.~J.~Steinhardt, L.~Wang, and I.~Zlatev.
\newblock Cosmological tracking solutions.
\newblock Phys.~Rev. D{\bf 59}, 123504 (1999).
\newblock DOI: 10.1103/PhysRevD.59.123504

\bibitem{tsu13} S.~Tsujikawa.
\newblock Quintessence: A review.
\newblock Classical Quantum Gravity {\bf 30}, 214003 (2013).
\newblock DOI: 10.1088/0264-9381/30/21/214003

\bibitem{waiell97} J.~Wainwright and G.~F.~R.~Ellis.
\newblock {\em Dynamical Systems in Cosmology}.
\newblock {C}ambridge {U}niversity {P}ress, Cambridge, (1997).

\bibitem{col03} A.~A.~Coley.
\newblock {\em Dynamical Systems and Cosmology}.
\newblock Kluwer Academic Publishers, Dordrecht, (2003).

\bibitem{wei08} S.~Weinberg.
\newblock {\em Cosmology}.
\newblock Oxford University Press, Oxford, (2008).

\bibitem{giamir10} R.~Giambo and J.~Miritzis.
\newblock Energy exchange for homogeneous and isotropic universes with a scalar field
coupled to matter.
\newblock Classical Quantum Gravity  {\bf 27}, 095003 (2010).
\newblock DOI: 10.1088/0264-9381/27/9/095003.

\bibitem{ugg13} C.~Uggla.
\newblock Global cosmological dynamics for the scalar field representation of the modified Chaplygin gas.
\newblock Phys. Rev. D{\bf 88}, 064040 (2013).
\newblock DOI: 10.1103/PhysRevD.88.064040.

\bibitem{leofad14} G.~Leon and C.~R.~Fadragas.
\newblock Cosmological Dynamical Systems: And their Applications.
\newblock LAP LAMBERT Academic Publishing (2012).
\newblock arXiv:1412.5701

\bibitem{tam14} N.~Tamanini.
\newblock Dynamical systems in dark energy models.
\newblock PhD Thesis, University College, London, (2014).

\bibitem{alhugg15} A.~Alho and C. Uggla.
\newblock Global dynamics and inflationary center manifold and slow-roll approximants.
\newblock J. Math. Phys. {\bf 56}, 012502 (2015).% -
\newblock DOI: 10.1063/1.4906081.

\bibitem{alhetal15} A.~Alho, J.~Hell and C.~Uggla.
\newblock Global dynamics and asymptotics for monomial scalar field potentials and perfect fluids.
\newblock Classical Quantum Gravity {\bf 32}, (14:145005) (2015)
\newblock DOI: 10.1088/0264-9381/32/14/145005

\bibitem{blaetal04} R.~D.~Blandford, M.~Amin, E.~A.~Baltz, K.~Mandel, and P.~J.~Marshall.
\newblock Cosmokinetics.
\newblock Observing Dark Energy, ASP Conference Series, {\bf 339} 17 (2005) 
\newblock DOI: 2005ASPC..339...27B
          
\bibitem{rapetal07} D.~Rapetti, S.~W.~Allen, M.~A.~Amin, and R.~D.~Blandford.
\newblock A kinematical approach to dark energy studies.
\newblock Mon.~Not.R~.~Astron.~Soc. {\bf 375}, 1510 (2007)%1510-1520
\newblock Doi: 10.1111/j.1365-2966.2006.11419.x.

\bibitem{copetal98} E.~J.~Copeland, A.~R.~Liddle, and D.~Wands.
\newblock Exponential potentials and cosmological scaling solutions.
\newblock Phys. Rev. D{\bf 57}, 4686 (1998).
\newblock DOI: 10.1103/PhysRevD.57.4686

\bibitem{Nunes2000} A.~Nunes and J.~P.~Mimoso.
\newblock On the potentials yielding cosmological scaling solutions.
\newblock Phys. Lett. B{\bf 488}, 423 (2000)
\newblock DOI: 10.1016/S0370-2693(00)00919-9

\bibitem{uggetal95} C.~Uggla, R.~T.~Jantzen, and K.~Rosquist.
\newblock Exact hypersurface-homogeneous solutions in cosmology and astrophysics.
\newblock Phys. Rev. D{\bf 51} 5522 (1995). %5522-5557
\newblock DOI: 10.1103/PhysRevD.51.5522

\bibitem{heiugg10} J.~M.~Heinzle and C.~Uggla.
\newblock Monotonic functions in Bianchi models: why they exist and how to find them.
\newblock Classical Quantum Gravity {\bf 27}, 015009 (2010).
\newblock DOI: 10.1088/0264-9381/27/1/015009

\bibitem{ngetal01} S.~C.~C.~Ng, N.~J.~Nunes, and F.~Rosati.
\newblock Applications of scalar attractor solutions to cosmology.
\newblock Phys.~Rev. D{\bf 64}, 083510 (2001).
\newblock DOI:10.1103/PhysRevD.64.083510.

\bibitem{ure12} L.~A.~Urena-Lopez.
\newblock Unified description of the dynamics of quintessential scalar fields.
\newblock J. Cosmol. Astropart. Phys. 03 (2012) 035.
\newblock DOI: 10.1088/1475-7516/2012/03/035.

\bibitem{gon14} Y.~Gong.
\newblock The general property of dynamical quintessence field.
\newblock Phys.~Lett. B{\bf 342} (2014) 342.% 342-349
\newblock DOI:10.1016/j.physletb.2014.03.013.

\bibitem{colwai92} A.~A.~Coley and J.~Wainwright
\newblock Qualitative analysis of two-fluid Bianchi cosmologies.
\newblock Classical Quantum Gravity {\bf 9} 651 (1992).
\newblock DOI: 10.1088/0264-9381/9/3/008

\bibitem{heietal05} J.~M.~Heinzle, N.~R{\"o}hr, and C.~Uggla.
\newblock Matter and dynamics in closed cosmologies.
\newblock Phys.~Rev. D{\bf 71} 083506 (2005).
\newblock DOI: 10.1103/PhysRevD.71.083506.



\end{thebibliography}



\end{document}





\bibitem{peerat03} P.~J.~E.~Peebles and B.~Ratra.
\newblock The Cosmological Constant and Dark Energy.
\newblock Rev.~Mod.~Phys. {\bf 75} 559 (2003).



\bibitem{zlaetal99} I.~Zlatev, L.~Wang, and P.~J.~Steinhardt.
\newblock Quintissence, Cosmic Coincidence, and the Cosmological Constant.
\newblock Phys.~Rev.~Lett. {\bf 82}, 896 (1999).


\bibitem{colwai92} A.~A.~Coley and J.~Wainwright
\newblock Qualitative analysis of two-fluid Bianchi cosmologies.
\newblock Class.~Quantum Grav. {\bf 9} 651 (1992).

\bibitem{heietal05} J.~M.~Heinzle, N.~R{\"o}hr, and C.~Uggla.
\newblock Matter and dynamics in closed cosmologies.
\newblock Phys.~Rev. {\bf D71} 083506 (2005).
\newblock DOI: 10.1103/PhysRevD.71.083506.









\bibitem{car81} J.~Carr.
\newblock {\em Applications of center manifold theory}.
\newblock Springer Verlag, New York, 1981.


\bibitem{ratqui92} B.~Ratra and A.~Quillen.
\newblock Gravitational Lensing Effects in a Time-Variable Cosmological `Constant'
Cosmology.
\newblock Mon.~Not.~Roy.~Astr.~Soc. {\bf 259}, 738 (1992).

\bibitem{bramar00} P.~Brax and J.~Martin.
\newblock The Robustness of Quintissence.
\newblock Phys.~Rev. {\bf D 61}, 103502 (2000).

\bibitem{braetal00} P.~Brax, J.~Martin, and A.~Riazuelo.
\newblock Exhaustive study of cosmic microwave background anisotropies in quintessential
scenarios.
\newblock Phys.~Rev. {bf D 62}, 103505 (2000).
\newblock DOI:10.1103/PhysRevD.62.103505.

\bibitem{uza99} J.-P.~Uzan.
\newblock Cosmological scaling solutions of nonminimally coupled scalar
fields.
\newblock Phys.~Rev. {\bf D 59}, 123510 (1999).
\newblock  DOI:10.1103/PhysRevD.59.123510.

\bibitem{lidshe98} A.~B.~Liddle and R.~J.~Scherrer.
\newblock Classification of scalar field potentials with cosmological scaling solutions.
\newblock Phys.~Rev. {\bf D 59}, 023509 (1998).
\newblock DOI:10.1103/PhysRevD.59.023509.

\bibitem{chi10} T.~Chiba.
\newblock Equation of state of tracker fields.
\newblock Phys.~Rev. {\bf D 81}, 023515 (2010).
\newblock DOI:10.1103/PhysRevD.81.023515.






********************

The stability of the fixed points is as follows:

On the $Z=0$ boundary the fixed points M$^{0}_{\pm}$ are sources, PL is a
sink, and FL$^{0}$ a saddle, see eg. \cite{}.

On the $Z=1$ boundary the fixed points M$^{1}_{\pm}$ are sources, dS is a
sink, and FL$^1$ a saddle.

Since on $\Sigma_{\dagger}>0$ ($\Sigma_{\dagger}<0$), Z is monotonically
increasing (decreasing), it follows that in the full state space M$^{0}_{+}$
and M$^{1}_{-}$ are sources, M$^{0}_{-}$ and M$^{1}_{+}$ are saddles, and Es
is a saddle.

We will assume that $\lambda_{0}^2<3\gamma_m$, which results in that the
fixed point $\mathrm{PL}$ is outside the physical state space and that
$\mathrm{Es}$ is stable \emph{on} the $Z=0$ subset, since the assumption
$\gamma_m<2$ implies $\lambda^2_0<6$. Moreover the power-law solution
$\mathrm{PL}$, is accelerating if and only if $\lambda_{0}^2<2$.


Besides the line of fixed points FL$_{Z}$ which is transversally hyperbolic
and the dS fixed point which is non-hyperbolic having a zero eigenvalue
corresponding to a 1-dimensional stable center  manifold, all the remaining
fixed points are hyperbolic.

Since they lie at the intersection of invariant boundaries subsets, their
stability readily follows from the analysis of the flow on these boundaries.



We will assume that $\bar{\lambda}^2<3\gamma_m$, which results in that the
fixed point $\mathrm{Ei}$ is outside the physical state space and that
$\mathrm{Es}$ is stable \emph{on} the $Z=0$ subset, which also leads to that
$q < - (3\gamma_m-2)/2$, i.e., this fixed point corresponds to an
accelerating state if $\gamma_m>2/3$, which is the case for e.g. a
pressureless fluid and a radiation equation of state. The solution space on
this subset is depicted in Fig.~\ref{fig:Z0}.

\begin{figure}[ht!]
\begin{center}
\subfigure[$\gamma_{m}=1$ and $\lambda_{-}=1$]{\label{fig:Z0Dust}
\includegraphics[width=0.5\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_Final.pdf}} \qquad
\end{center}
\caption{Solutions on the invariant $Z=0$ boundary subset of $\bar{\bf S}$.
}
\label{fig:Z0}
\end{figure}

%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=\sqrt{\frac{3}{2}}$]{\label{fig:Z0Dust_1}
\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_Final.pdf}} \hspace{0.2cm}
\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=0.2$]{\label{fig:Z0Dust_2}
\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_v2_Final.pdf}}
\end{center}
\caption{Solutions on the invariant $Z=0$ boundary subset of $\bar{\bf S}$.
In Fig. (a) the space-dot line corresponds to $a=10$, the dotted line to $a=1$, the dash-dotted line to $a=0$,
the dashed line to $a=-1/2$ and the space-dotted line to $a=-2$. Fig. (b) represent numerical solutions.}
\label{fig:Z0}
\end{figure}
%





Note that since
%
\begin{equation}
H^\prime = - (1+q)H,
\end{equation}
%
it follows that $H$ is monotonically decreasing, except when $q=-1$, which is
only possible on the scalar field boundary. This will have significant
consequences for the global dynamics of the models we consider and will be
further discussed below. Here we note that
%
\begin{equation}
H^2 = \frac{V(\phi(Z))}{3\Omega_V} = \frac{V(\phi(Z))}{3(1 - x^2 - \Omega_m)}, \qquad (H^2)^\prime = - 2(1+q)H^2.
\end{equation}
%
Finally, note that $Z$ is monotonically increasing (decreasing) when $x>0$
($x<0$).

%
\begin{subequations}
\begin{align}
\gamma_\phi &= \frac{\rho_\phi + p_\phi}{\rho_\phi} =
\frac{\dot{\phi}^2}{\sfrac12\dot{\phi}^2 + V(\phi)},\\
\gamma_{\mathrm{tot}} &= \frac{\rho_{\mathrm{tot}} + p_{\mathrm{tot}}}{\rho_{\mathrm{tot}}} =
\gamma_\phi\Omega_\phi + \gamma_m\Omega_m,\label{gammatot}\\
\Omega_\phi &= \frac{\kappa\rho_\phi}{3H^2},\\
\Omega_m &= \frac{\kappa\rho_m}{3H^2},
\end{align}
\end{subequations}
%


To make contact with $\Lambda$CDM cosmology we first consider a fluid without
pressure and a cosmological constant. Defining $\Omega_m = \rho_m/3H^2$ and
$\Omega_\Lambda = \Lambda/3H^2$ leads to that the Friedmann equation, given
by $3H^2 = \rho_m + \Lambda$ (obtainable from~\eqref{Gauss1}), yields
$\Omega_m+\Omega_\Lambda=1$, while~\eqref{rhosol} yields
$\rho_m=\rho_0(a/a_0)^{-3} = \rho_0\exp(-3\tau)$, where we have defined
$\tau$ according to $\tau = \ln(a/a_0)$. This, together with $3H^2 = \rho_m +
\Lambda$ and the definitions, gives
%
\begin{equation}
\Omega_m = \frac{\Omega_{m0}\exp(-3\tau)}{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}}, \qquad
\Omega_\Lambda = \frac{\Omega_{\Lambda0}}{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}},
\end{equation}
%
where $\Omega_{m0}$ and $\Omega_{\Lambda0}$ are the values of $\Omega_{m}$
and $\Omega_{\Lambda}$ at a reference time $\tau=0$, which, if chosen to be
now, yields $\tau = -\ln(1+z)$ where $z$ is the redshift. Thus $\Omega_m$ is
monotonically decreasing from 1 in the limit $\tau \rightarrow -\infty$,
corresponding to the big bang, to zero when $\tau \rightarrow + \infty$,
i.e., the infinite future which is described by a de Sitter state with
$\Omega_\Lambda=1$. This can also be seen by considering the differential
equation for $\Omega_m$, which is given by (see e.g. p 62 in~\cite{waiell97})
%
\begin{equation}
\Omega_m^\prime = - 3(1-\Omega_m)\Omega_m,
\end{equation}
%
where we in this paper use a $^\prime$ for the $\tau$-derivative.

%
\begin{subequations}
\begin{align}
q &= -H^{-2}\left(\frac{\ddot{a}}{a}\right) = -\left(1 + \frac{H^\prime}{H}\right),\\
j &= H^{-3}\left(\frac{\dddot{a}}{a}\right) = -q^\prime + q + 2q^2,
\end{align}
\end{subequations}
%
which in the present $\Lambda$CDM case are given by
%
\begin{equation}
q = -1 + \frac32\Omega_m = -1 + \frac32\left(\frac{\Omega_{m0}\exp(-3\tau)}{\Omega_{m0}\exp(-3\tau) + \Omega_{\Lambda0}}\right),
\qquad j=1.
\end{equation}
%






%
\begin{subequations}
\begin{align}
\gamma_\phi &= \frac{\rho_\phi + p_\phi}{\rho_\phi} =
\frac{\dot{\phi}^2}{\sfrac12\dot{\phi}^2 + V(\phi)},\\
\gamma_{\mathrm{tot}} &= \frac{\rho_{\mathrm{tot}} + p_{\mathrm{tot}}}{\rho_{\mathrm{tot}}} =
\gamma_\phi\Omega_\phi + \gamma_m\Omega_m,\label{gammatot}\\
\Omega_\phi &= \frac{\kappa\rho_\phi}{3H^2},\\
\Omega_m &= \frac{\kappa\rho_m}{3H^2},
\end{align}
\end{subequations}
%
where
%
\begin{subequations}
\begin{alignat}{2}
\rho_{\mathrm{tot}} &= \rho_\phi + \rho_m, &\qquad
p_{\mathrm{tot}} &= p_\phi + p_m,\\
\rho_\phi &= \sfrac12\dot{\phi}^2 + V(\phi),&\qquad p_\phi &=
\sfrac12\dot{\phi}^2 - V(\phi).
\end{alignat}
\end{subequations}
%
and where we note that $\Omega_{\mathrm{tot}} =
\kappa\rho_{\mathrm{tot}}/3H^2 = 1$ due to eq.~\eqref{Gauss}. As a
consequence of eqs.~\eqref{qdef},~\eqref{Gauss}, and~\eqref{gammatot},
%
\begin{equation}
1+q = \frac32\gamma_{\mathrm{tot}} = \frac32\left(\gamma_\phi\Omega_\phi + \gamma_m\Omega_m\right).
\end{equation}
%

Expressed in Hubble-normalized quantities, the field equations yield the
following evolution equations
%
\begin{subequations}
\begin{align}
\dot{a} &= aH,\\
\dot{\phi} &= \sqrt{\frac{6}{\kappa}}H\Sigma_\dagger,\\
\dot{\Sigma}_\dagger &= -(2-q)H\Sigma_\dagger - \sqrt{\frac{\kappa}{6}}\frac{V_{,\phi}}{H} =
-H\left[(2-q)\Sigma_\dagger - \sqrt{\frac{3}{2}}\lambda\Omega_V\right],\\
\dot{\Omega}_V &= 2(1+q)H\Omega_V - \sqrt{\frac{6}{\kappa}}\left(\frac{\kappa V_{,\phi}}{3H^2}\right)H\Sigma_\dagger
= H\left[2(1+q) - \sqrt{6}\lambda\Sigma_\dagger\right]\Omega_V ,\\
\dot{\Omega}_m &= 3H(\gamma_\phi - \gamma_m)\Omega_\phi\Omega_m = 3H(2\Sigma_\dagger^2 - \gamma_m\Omega_\phi)\Omega_m,
\end{align}
\end{subequations}
%
where
%
\begin{equation}
\lambda = - \frac{V_{,\phi}}{\sqrt{\kappa}V},
\end{equation}
%
while the Gauss constraint takes the form
%
\begin{equation}
1 = \Omega_\phi + \Omega_m = \Sigma_\dagger^2 + \Omega_V + \Omega_m.
\end{equation}
%

Let us now use $\tau$ as time variable where
%
\begin{equation}
\frac{dt}{d\tau} = H^{-1}.
\end{equation}
%
Then
%
\begin{subequations}
\begin{align}
\frac{d\ln a}{d\tau} &= 1,\\
\frac{d\phi}{d\tau} &= \sqrt{6}\Sigma_\dagger,\\
\frac{d\Sigma_\dagger}{d\tau} &= -(2-q)\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda\Omega_V,\\
\frac{d\Omega_V}{d\tau} &= 2\left[1+q - \sqrt{\frac{3}{2}}\lambda\Sigma_\dagger\right]\Omega_V ,\\
\frac{d\Omega_m}{d\tau} &= 3(\gamma_\phi - \gamma_m)\Omega_\phi\Omega_m = 3(2\Sigma_\dagger^2 - \gamma_m\Omega_\phi)\Omega_m,\\
\frac{d\Omega_\phi}{d\tau} &= -3(\gamma_\phi - \gamma_m)\Omega_\phi\Omega_m = -3(2\Sigma_\dagger^2 - \gamma_m\Omega_\phi)\Omega_m,\\
\frac{dw_\phi}{d\tau} &= -(1 - w_\phi)\left[3(1 + w_\phi) - \sqrt{6}\lambda\Sigma_\dagger\right],\\
\frac{d\gamma_\phi}{d\tau} &= -(2 - \gamma_\phi)\left[3\gamma_\phi - \sqrt{6}\lambda\Sigma_\dagger\right],\\
\frac{d\lambda}{d\tau} &= -\sqrt{6}\Sigma_\dagger(\Gamma - 1)\lambda^2,
\end{align}
\end{subequations}
%
where
%
\begin{equation}
\Gamma = \frac{V\,V_{,\phi\phi}}{V_{,\phi}^2}.
\end{equation}
%

For multiple fluids
%
\begin{subequations}
\begin{align}
\frac{d\ln a}{d\tau} &= 1,\\
\frac{d\phi}{d\tau} &= \sqrt{6}\Sigma_\dagger,\\
\frac{d\Sigma_\dagger}{d\tau} &= -(2-q)\Sigma_\dagger + \sqrt{\sfrac{3}{2}}\lambda\Omega_V,\\
\frac{d\Omega_V}{d\tau} &= \left[2(1+q) - \sqrt{6}\lambda\Sigma_\dagger\right]\Omega_V ,\\
\frac{d\Omega_{mi}}{d\tau} &= \left[2q - (3\gamma_i-2)\right]\Omega_{mi}.
\end{align}
\end{subequations}
%


\end{document}

Current cosmology is an observationally and theoretically very active field
dominated by "standard cosmology," which is based on two ingredients:
spatially isotropic, homogeneous and flat Friedmann-Robertson-Walker (FLRW)
background models in general relativity (GR), and perturbations thereof. In
particular, $\Lambda$CDM and inflationary cosmology take key positions. The
standard scenario is remarkably successful, although some anomalies and
tensions between different data sets exist, but the origins and physical
nature of dark matter and energy, which according to standard cosmology
constitute about 95\% of the matter content in the observable universe, are
mysteries that are among the most intriguing ones in present day physics. The
primary purpose of this part of the research programme is to shed light on
why standard cosmology and the inflationary scenario are so successful, and
to find new nonlinear effects. Firstly, by situating standard cosmology in a
broader context as regards background models. The aim here is to critically
assess various assumptions and find new dynamical mechanisms that shed light
on the success of standard cosmology. This is to be accomplished by means of
novel ideas concerning: a) the role of physical scales in a hierarchy of
increasingly general background models, b) how to use physical scales in
order to formulate the associated equations so that it is possible to perform
an analysis that reveals the possibilities and consequences of background
cosmology, as well as leading to new analytic approximation schemes, e.g.
complementing the slow-roll approximation, which has and is playing a
significant role in inflationary cosmology. Secondly, by investigating
perturbations, especially higher order perturbations. The purpose with this
is to study the origins and consequences of nonlinearities on the cosmic
microwave background (CMB) and large scale structures, e.g., clarifying the
role and quantitative contributions of linear primordial and higher order
perturbations to the CMB B-mode polarization. In addition this part of the
programme aims at linking relativistic theory with Newtonian gravity, e.g.
setting initial conditions for structure formation at nonlinear order that
are consistent with relativistic cosmology.

Observations of the CMB, large scale structures, distant galaxies and Type Ia
supernovae, light element abundances, and, e.g., various astrophysical data
sets, yield increasingly restrictive demands on cosmological models, see e.g.
[1] and references therein. Nevertheless, mysteries abound. What are the
origins and nature of dark matter, dark energy, and inflation? How does
inflation end (what is the physics of reheating)? Why do we live in an era
when the dark matter and dark energy ratio is about 1:3 (the coincidence
problem)? Why is there such an enormous discrepancy between estimates of the
vacuum energy based on particle physics, dimensional analysis, and
cosmological observations (the cosmological constant problem)? These issues
have generated a multitude of theories of gravity, different couplings
between geometry and matter fields, and a variety of matter models and
interactions, see e.g. [2] and references therein. To minimize model
dependence, one can use observations to describe the evolutionary history of
the universe and construct e.g. a suitable matter equation of state, or a
scalar field potential, that yields this evolution [3]. However, such
reconstructed histories are hampered by derivatives of noisy data.
Comparisons with observations thereby lead to "priors," which eventually
result in more specific cosmological models. But what types of models fit
observational data - only $\Lambda$CDM cosmology? Or are there dynamical
mechanisms that produce classes of observationally competitive approximations
of $\Lambda$CDM cosmology, thus providing a useful testbed for $\Lambda$CDM
cosmology?

Below I will describe how I have recently started to tackle the above issues
for background FLRW models. This is partly accomplished by exploring how
different physical scales affect dynamical possibilities and properties and
how physical scales can be used to produce new powerful dynamical systems
formulations of the field equations for increasingly wide sets of models.
This novel approach for understanding background dynamics reveals that some
of the current attempts to resolve the above issues misinterpret presently
available results, and, more importantly, some essential features have
apparently been missed, as discussed below. Furthermore, the new formulations
make it possible to introduce a number of mathematical tools that have not
been used previously in cosmology, e.g., new approximation methods, which
e.g. complement and improve the slow-roll approximation.

The mystery of dark matter and dark energy has resulted in a plethora of
gravity theories, matter sources and interactions, which are observationally
tested as regards the evolution of FLRW models and the structures that arise
from perturbations thereof. Since it is the background FLRW models that
constitute the starting point for this, it is important to have a good
understanding of the various possibilities that exist in FLRW cosmology. One
of the first attempts to dynamically model dark energy was by a minimally
coupled scalar field (quintessence) in GR with an exponential potential,
while dark and normal matter were described by a perfect fluid with a linear
equation of state. Due to the scale-invariant properties of this type of
model, it was possible and fruitful to "Hubble-normalize" the problem by
appropriately scaling the variables with the Hubble variable (i.e. the
logarithmic time derivative of the cosmological scale factor) in order to
obtain a reduced dimensionless regular dynamical system. When certain matter
conditions are satisfied, this dynamical system has a future attractor
described by a fixed point, which corresponds to a scale-invariant solution
called an attractor solution. The matter and dark energy form a fixed ratio
for this solution, which has resulted in ideas associated with
scaling/tracker solutions and tracking behaviour [2]. In my opinion, this is
misleading. By necessity scale invariance must be broken if field theoretic
models are going to resemble and generalize $\Lambda$CDM cosmology in order
to be observationally competitive, since the existence of a cosmological
constant introduces a scale. Breaking scale invariance by a
non-scale-invariant scalar field potential, or by an effective potential
arising from some modified gravity theory, or by other types of fields,
necessarily leads to that one has to use other variables than
Hubble-normalized ones in order to obtain useful dynamical systems
formulations that reveal the solution space structure and thereby the
properties of a given model; furthermore, it also inevitably results in that
concepts such as attractor and tracker solutions must change, if appropriate
at all.

To provide a simple first warm up example of how one can use the physical
structure of a problem to reformulate the field equations into a regular
dynamical system on a global bounded state space, thus making it possible to
fully understand a background model, I considered the scalar field model for
a modified Chaplygin gas [8]. This problem is obtained by reconstructing a
minimally coupled scalar field potential from the perfect fluid description
of a modified Chaplygin gas. Since the scalar field problem has one more
degree of freedom than the original perfect fluid problem, the original
perfect fluid solution appears as a particular solution in the scalar field
state space. This poses a problem for scalar field reconstructions from
observational data that gives a description of the evolutionary history of
the universe (which the modified Chaplygin perfect fluid aims to do). The
scalar field model is only observationally relevant if the particular
"reconstruction solution" or history is an "attractor solution" (i.e. other
solutions are in some sense attracted to it during their evolution). Previous
work had considered the evolution along a single dimensional variable and
shown that it decayed along the modified Chaplygin perfect fluid solution. On
the other hand, the global dynamical systems analysis revealed that there is
an open set of scalar field solutions that never is close to the modified
Chaplygin solution and an open set that is during intermediate and late
times; the latter set, however, when expressed in physical dimensionless
quantities, is only attracted along part of the "attractor" Chaplygin
solution, not the whole solution. The underlying reason that solutions
"follow it" turns out to be due to another quite general feature: an
"attractor" solution is a one-dimensional unstable manifold of a fixed point
saddle on a boundary subset in an appropriately formulated dynamical system
that respects the physical scales of the problem, and it ends at a global
mathematically well-defined future attractor, which all solutions in the
physical state space approach. This illustrates that although fine-tuning
problems often are discussed in cosmology, they nevertheless have been
misunderstood and underestimated; new insights are needed if such problems
are to be solved. As a second illustrative example I and some collaborators
considered perfect fluids and a minimally coupled scalar field with a
monomial potential [9]. It should be pointed out that the recent analysis of
the CMB casts some doubt on the physical relevance of these models as
inflationary models, but this is not the point. We used these models as
examples of how one can use inherent physical scales to reformulate a problem
as a global regular dynamical system on a compact state space (next I will
discuss how similar ideas yield new observationally interesting results). The
novel formulation made it possible to import local and global dynamical
systems techniques, and resulted in a visual global description and complete
understanding of the solution space and asymptotic behaviour. At late times
we modified mathematical averaging techniques in order to prove statements
about how fluid and scalar field properties affect "asymptotic source
dominance" and asymptotic self-similarity breaking, where the latter is an
important physical feature with connections to other fields of physics. We
also showed that the "attractor solution" of these models is connected with
the slow-roll approximation, which turns out to be associated with the
tangency condition for a one-dimensional unstable centre manifold of a de
Sitter fixed point on a boundary subset that describes the dynamics at early
times. In fact this solution is only attracting nearby solutions for some of
the parameter values of this class of models; for other values it is only
attracting solutions for a brief time. We also showed how one can use methods
from dynamical systems theory to find approximations that complement the
slow-roll approximation, and noticed that the accuracy of the slow-roll
approximation quite unpredictably depends on the shape of the potential, and
thus this approximation should be used with care. As stated, this is just an
example; similar ideas to those introduced in [9] can be modified and
employed in quite general situations, revealing quite unexpected results that
shed new light on the possibilities to model e.g. dark energy. In an
influential paper [12], Ratra and Peebles argued that dark energy can be
phenomeno-logically modelled by a minimally coupled scalar field with a
shallow potential that leads to that the dark field energy density decreases
more slowly than the matter energy density. In particular they considered an
inverse power-law potential and assumed a matter dominated flat FLRW universe
after inflation, which resulted in a set of approximate flat FLRW equations
for which they found a particular solution, which was shown to be linearly
stable within the approximate context. The associated solution of the exact
FLRW equations has subsequently been referred to as an attractor or tracker
solution [13], and has resulted in numerous papers, although lately its
observational relevance has been questioned [2]. Recently I used the physical
scales associated with this problem to recast the field equations to a
regular dynamical system on a three-dimensional compact state space, which
has yielded a complete understanding of the solution space and its
properties. In particular, it was shown that the tracker solution originates
from a fixed point that corresponds to a flat FLRW model with matter only and
that it ends at a de Sitter fixed point on a boundary subset, which
constitutes the future global attractor of all solutions with a scalar field
and a perfect fluid. Moreover, before the future asymptotic state the tracker
solution shadows a special solution on the scalar field boundary subset that
corresponds to a one-dimensional stable center manifold of the de Sitter
fixed point. This latter solution is analogous to the attractor solution in
inflationary models, and its tangent space is described by the slow-roll
approximation. These features, together with the approximate linear stability
at early times, are the underlying reasons "over- and undershooting"
solutions roughly behave as the tracker solution. However, it turns out that
these solutions are not the most observationally interesting ones these
models admit.

The analysis of the new dynamical system reveals that there exists an orbit
(solution trajectory) on a boundary subset that is associated with
intermediate and late time dynamics that turns out to exactly describe
$\Lambda$CDM cosmology! It follows that the models admit an open set of
solutions that shadow this solution trajectory (which in addition is stable),
after having shadowed an early pure perfect fluid evolution, which in turn
yields a set of solutions that continuously deforms the $\Lambda$CDM models,
thereby being the most observationally interesting ones! It should be pointed
out that in contrast to the tracker solution, these solutions do not come
close to the above mentioned centre manifold; estimates that can be found in
the literature that use the slow-roll approximation, for e.g. the tracker
solution, are therefore not applicable for this new class of solutions. The
$\Lambda$CDM shadowing solutions have a continuous spectrum of much less dark
energy initially than the tracker and the over- and undershooting solutions,
and apparently they have been missed because this problem has not been
addressed adequately before. Moreover, in contrast to producing field
theoretic models from a given evolutionary history, which leads to
fine-tuning issues, the new type of dynamical systems approach I have
introduced suggests how one can determine properties that any field theoretic
models must have in order to yield an open set of solutions that continuously
deform $\Lambda$CDM cosmology, thus also somewhat alleviating fine-tuning
problems. It still, however, remains to be explained how inflation yields a
very small dark energy content after a reheating epoch, which presumably
involves unknown interactions between the inflaton and other matter fields,
but  the new developments will still give some restrictions concerning this
issue.



****************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Dynamical systems analysis}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-----------------------------------------------------------------
\subsection{Some global aspects}
%-----------------------------------------------------------------

We first note that by using the Gauss constraint~\eqref{Gauss2} we obtain the
following relation
%
\begin{equation}
2-q = 3\Omega_V + \frac32(2-\gamma_m)\Omega_m,
\end{equation}
%
and since we have assumed $0<\gamma_m<2$, $\rho_m\geq 0$, $V\geq 0$, it
follows from this equation and~\eqref{q} that
%
\begin{equation}\label{qineq}
-1 \leq q \leq 2.
\end{equation}
%
Furthermore, $q=-1$ implies that $\Sigma_\dagger = \Omega_m = 0, \Omega_V=1$,
which amounts to a de Sitter state, which corresponds to a fixed point
discussed below, while $q=2$ implies that $\Sigma_\dagger = \pm 1, \Omega_m =
\Omega_V = 0$, which corresponds to two one-dimensional boundary subsets
obtained by the intersection of the $\Omega_m=0$ and $Omega_V=0$ boundary
subsets, yielding a pure massless scalar field state.

We then note that since
%
\begin{equation}
H^\prime = - (1+q)H,
\end{equation}
%
it follows that $H$ is monotonically decreasing, except when $q=-1$; however,
this is only possible on the scalar field boundary, and then this only yields
an inflection point in the graph for $H$ when $q=-1$ is not an invariant
subset, which is only the case for the so-called de Sitter fixed point
discussed below on the boundary subset $Z=1$. By means of this monotonicity
property and~\eqref{Gauss1} it follows that for the present models $H$ is
monotonically decreasing on ${\bf S}$, with $H \rightarrow \infty$ when $\tau
\rightarrow - \infty$ and $H \rightarrow 0$ when $\tau \rightarrow + \infty$.
The monotonicity of $H$ induces a monotone function on ${\bf S}$ due to that
%
\begin{equation}
H^2 = \frac{V(\phi(Z))}{3\Omega_V} = \frac{V(\phi(Z))}{3(1 - \Sigma_\dagger^2 - \Omega_m)}, \qquad (H^2)^\prime = - 2(1+q)H^2.
\end{equation}
%
As a consequence there can be no fixed points, recurring orbits, etc, in
${\bf S}$ and fixed points and the past and future attractors must reside on
the boundary of ${\bf S}$.

Due to that $H\rightarrow \infty$ toward the past it follows that either
$V(\phi(Z)) \rightarrow \infty$, which due to the assumed monotonicity of $V$
corresponds to $\phi \rightarrow - \infty$ and hence $Z\rightarrow 0$, or
that $\Omega_V = 1 - \Sigma_\dagger^2 - \Omega_m$, and therefore the past
attractor for all interior orbits in ${\bf S}$ resides on either $Z=0$ or
$\Omega_V=0$. Furthermore, since $H\rightarrow 0$ toward the future it
follows that $V(\phi(Z)) \rightarrow 0$, which due to the assumed
monotonicity of $V$ corresponds to $\phi \rightarrow + \infty$ and hence
$Z\rightarrow 1$, and therefore the future attractor for all interior orbits
in ${\bf S}$ resides on $Z=1$.

Moreover, since we assume that $\lambda > 0$ and $\lim_{\phi\rightarrow
\infty} \lambda = 0$, it follows from~\eqref{qineq} and~\eqref{xeq} that if
$\Sigma_\dagger <0$ then $\Sigma_\dagger$ is monotonically increasing until
$\Sigma_\dagger \geq 0$. Furthermore, it follows from eq.~\eqref{Zeq0} that
$Z$ is monotonically increasing if $\Sigma_\dagger>0$.

%-----------------------------------------------------------------
\subsection{Invariant subsets}
%-----------------------------------------------------------------

%-----------------------------------------------------------------
\subsubsection*{Fixed points}
%-----------------------------------------------------------------

The dynamical system~\eqref{dynsys} admits a number of fixed points on
$\bar{\bf S}$, all located on the boundary, which are given in
Table~\ref{tab:fixedpoint}.
%
\begin{center}
\begin{tabular}{|c|c|c|}\hline\label{tab:fixedpoint}
 &  Fixed point & \\ \hline
$\mathrm{M}_0^\pm$ & $\Sigma_\dagger = \pm 1$, $\Omega_m = \Omega_V = Z =0$ & $q = 2$ \\ \hline
$\mathrm{M}_1^\pm$ & $\Sigma_\dagger = \pm 1$, $\Omega_m = \Omega_V = 0$, $Z=1$ & $q = 2$ \\ \hline
$\mathrm{Es}$ & $\Sigma_\dagger = \lambda_{0}/\sqrt{6}$, $\Omega_m = Z = 0$ & $q = -1 + \lambda_{0}^2/2$ \\ \hline
$\mathrm{PL}$ & $\Sigma_\dagger = \sqrt{3}\gamma_m/\sqrt{2}\lambda_{0}$, $\Omega_m = 1 - 3\gamma_m/\lambda_{0}^2$, $Z = 0$ & $q =(3\gamma_m-2)/2$ \\ \hline
$\mathrm{FL}_Z$ & $\Sigma_\dagger = \Omega_V= 0$, $\Omega_m = 1$, $Z=Z_0=\mathrm{const.}$ & $q = (3\gamma_m-2)/2$ \\ \hline
$\mathrm{dS}$ & $\Sigma_\dagger = \Omega_m = 0$, $\Omega_V = Z = 1$ & $q = -1$ \\ \hline
\end{tabular}
\end{center}
%
We have here introduced a notation for the fixed points where the kernel is
an abbreviation for the solutions they correspond to in terms of their state
space variables while a subscript suggests a value for $Z$ while a
superscript designates if $\Sigma_\dagger = 1$ or $\Sigma_\dagger = -1$; thus
the kernel $\mathrm{M}$ stands for a massless scalar field state and no
perfect fluid; $\mathrm{PL}$ for the matter scaling solution arising for
exponential potentials; $\mathrm{PL}$ for the power-law solution of
exponential potentials and no matter; $\mathrm{FL}$ stands for a
Friedmann-Lema\^{i}tre perfect fluid state, described by a line of fixed
points, while $\mathrm{dS}$ stands for a de Sitter state.

We will assume that $\lambda_{0}^2<3\gamma_m$, which results in that the
fixed point $\mathrm{PL}$ is outside the physical state space and that
$\mathrm{Es}$ is stable \emph{on} the $Z=0$ subset, since the assumption
$\gamma_m<2$ implies $\lambda^2_0<6$. Moreover the power-law solution
$\mathrm{PL}$, is accelerating if and only if $\lambda_{0}^2<2$.


Besides the line of fixed points FL$_{Z}$ which is transversally hyperbolic
and the dS fixed point which is non-hyperbolic having a zero eigenvalue
corresponding to a 1-dimensional stable center  manifold, all the remaining
fixed points are hyperbolic.

Since they lie at the intersection of invariant boundaries subsets, their
stability readily follows from the analysis of the flow on these boundaries.

%-----------------------------------------------------------------
\subsubsection*{The $Z=0$ subset}
%-----------------------------------------------------------------

At the exponential potential boundary subset $Z=0$  we have $\lambda(0) =
\lambda_{0}$ and therefore the dynamical system
%
\begin{subequations}\label{dynsys0}
\begin{align}
\Sigma_\dagger^\prime &= -(2-q)\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda_{0}(1 - \Sigma_\dagger^2 - \Omega_m),\label{Sigexp}\\
\Omega_m^\prime &= 3\left[2\Sigma_\dagger^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,
\end{align}
\end{subequations}
%
where we recall that $2-q=3(1-\Sigma_\dagger^2) - \frac32 \gamma_m\Omega_m$.

>From the previous stated assumptions, it follows that on this boundary,
$\mathrm{M}^{\pm}_{0}$ are sources, $\mathrm{FL}_{0}$ is a saddle, and
$\mathrm{PL}$ is a sink.

For $\lambda_{0}=\sqrt{\frac{3}{2}}$ and $\gamma_m=1$ the solutions are
explicitly given in terms of proper synchronous time by
\begin{equation}
\Sigma_{\dagger}=\frac{2t^3-a}{4t^3+6t+a}\qquad,\qquad
 \Omega_{m}=\frac{1-\Sigma^{2}_{\dagger}}{1+t^2}=\frac{12t\left(t^3+3t+a\right)}{\left(4t^3+6t+a\right)^2}
\end{equation}



The solution space on this subset is depicted in Fig.~\ref{fig:Z0}.

\begin{figure}[ht!]
\begin{center}
\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=\sqrt{\frac{3}{2}}$]{\label{fig:Z0Dust_1}
\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_Final.pdf}} \qquad
\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=0.2$]{\label{fig:Z0Dust_2}
\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_v2_Final.pdf}}
\end{center}
\caption{Solutions on the invariant $Z=0$ boundary subset of $\bar{\bf S}$.
In Fig. (a) the space-dot line corresponds to $a=10$, the dotted line to $a=1$, the dash-dotted line to $a=0$,
the dashed line to $a=-1/2$ and the space-dotted line to $a=-2$. Fig. (b) represent numerical solutions.}
\label{fig:Z0}
\end{figure}

%-----------------------------------------------------------------
\subsubsection*{The $\Omega_m=0$ subset}
%-----------------------------------------------------------------

The pure scalar field subset yields the following dynamical system:
%
\begin{subequations}\label{dynsyssc}
\begin{align}
\Sigma_\dagger^\prime &= \left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(Z)\right)(1-\Sigma_\dagger^2),\label{Sigeqsc}\\
Z^\prime &= \sqrt{6}\lambda_{0} Z(1-Z)\Sigma_\dagger.\label{Zeqsc}
\end{align}
\end{subequations}
%


In order to have a picture of the scalar field boundary a particular example
of a potential must be choosen.


\begin{itemize}
 \item {\bf Center manifold and slow-roll approximation}
\end{itemize}


A center manifold analysis of $\mathrm{dS}$ applied to the
system~\eqref{dynsys1} reveals that the center manifold resides on the scalar
field boundary $\Omega_m=0$ (for a pedagogical introduction to center
manifolds, see~\cite{car81}, and for recent applications in reminiscent
cosmological situations, see~\cite{alhugg14,alhugg14b}). Because of this we
will restrict the analysis to this boundary subset and
eq.~\eqref{dynsysscalar}. Without loss of generality we locate the fixed
point dS at the origin, by introducing the variable
$$\tilde{Z}=1-Z$$

% First we Taylor expand $\lambda(Z)$ around $Z=1$
% \begin{equation}
%  \lambda(Z)=\lambda(1)+\sum_{k=1}\frac{1}{k!}\frac{d^{k}\lambda}{dZ^{k}}|_{Z=1}(Z-1)^k = 0 +\bar{\lambda}_1 (Z-1)+...
% \end{equation}
%
%

%Introducing $\tilde{Z}=1-Z$ the system reads
and the system becomes
\begin{subequations}
\begin{align}
\Sigma_\dagger^\prime &= \left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(\tilde{Z})\right)(1-\Sigma_\dagger^2),\\
\tilde{Z}^\prime &= -\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger.
\end{align}
\end{subequations}
We now Taylor expand $\lambda(\tilde{Z})$ around $\tilde{Z}=0$
\begin{equation}
\lambda(\tilde{Z})=
\lambda(0)+\sum_{k=1}\frac{1}{k!}\frac{d^{k}\lambda}{d\tilde{Z}^{k}}|_{\tilde{Z}=0} \tilde{Z}^k =
0 +\tilde{\lambda}_1 \tilde{Z}+\frac{1}{2}\tilde{\lambda}_{2}\tilde{Z}^2 +\cdots
\end{equation}
A linearization yields
%
\begin{subequations}
\begin{align}
E^s &= \{(\tilde{Z},\Sigma_{\dagger})|\,\tilde{Z}=0\},\\
E^c &= \{(\tilde{Z},\Sigma_{\dagger})|\, \frac{\sqrt{6}}{\tilde{\lambda}_{1}}\Sigma_{\dagger}-\tilde{Z} = 0\},
\end{align}
\end{subequations}

To investigate the center manifold $W^c$ associated with the tangent space
$E^c$, we adapt the variables to this tangent space and replace
$\Sigma_\dagger$ with
%
\begin{equation}
 v = \frac{\sqrt{6}}{\tilde{\lambda}_1} \Sigma_{\dagger} - \tilde{Z},
\end{equation}
%
where $\mathrm{dS}$ is located $(\tilde{Z},v) = (0,0)$. This leads to the
transformed system
%
\begin{subequations}
\begin{align}
\frac{d\tilde{Z}}{d\tau} &=-\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger \\
                     %    &=-\tilde{\lambda}_{1}\lambda_{0} (1-\tilde{Z})\tilde{Z}(\tilde{Z}+v), \label{tildeZeq}\\
\frac{dv}{d\tau} &=\frac{\sqrt{6}}{\tilde{\lambda}_1}\left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(Z)\right)(1-\Sigma_\dagger^2)+\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger, \\
               %  &=\left[-3v+\frac{3}{\tilde{\lambda}_1}\left(\frac{1}{2}\tilde{\lambda}_2\tilde{Z}^2+\cdots\right)\right]\left(1-\frac{\tilde{\lambda}^{2}_{1}}{6}\left(\tilde{Z}+v\right)^2\right)
                %   +\tilde{\lambda}_{1}\lambda_{0} (1-\tilde{Z})\tilde{Z}(\tilde{Z}+v).\label{veq}
\end{align}
\end{subequations}
%
where $\Sigma_{\dagger}=\frac{\tilde{\lambda}_1}{\sqrt{6}}(\tilde{Z}+ v)$ The
linearization of eq.~\eqref{veq} yields $\frac{d v}{d\tau} = - 3 v$, while
eq.~\eqref{tildeZeq} only has higher order terms. The center manifold $W^c$
can be obtained as the graph $ v =
 h(\tilde{Z})$ near $(\tilde{Z}, v) = (0,0)$ (i.e., use
$\tilde{Z}$ as an independent variable), where $ h(0) =0$ (fixed point
condition) and $\frac{d h}{d\tilde{Z}}(0) =0$ (tangency condition). Inserting
this relationship into the above equations and using $\tilde{Z}$ as the
independent variable leads to

\begin{equation}\label{tildehZ}
-\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger\left(1+\frac{d h}{d\tilde{Z}}\right)
-\frac{\sqrt{6}}{\tilde{\lambda}_1}\left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(\tilde{Z})\right)(1-\Sigma_\dagger^2) = 0
\end{equation}

where $\Sigma_{\dagger}=\frac{\tilde{\lambda}_1}{\sqrt{6}}(\tilde{Z}+
h(\tilde{Z}))$. Solving this equation amounts to finding the center manifold
solution, which presumably is not possible. Instead we replace $h(\tilde{Z})$
with a truncated power series in $\tilde{Z}$. Solving the ensuing algebraic
equations order by order yields

 \begin{equation}\label{CMZ_Expansion}
h(\tilde{Z}) = \frac{1}{6}\frac{2\tilde{\lambda}^{2}_{1}\lambda_{0}+
3\tilde{\lambda}_{2}}{\tilde{\lambda}_1} \tilde{Z}^2 + \dots %\frac{\alpha^2}{108}\tilde{H}^4 + \dots
\quad \Rightarrow \quad
\Sigma_{\dagger} \approx \frac{\tilde{\lambda}_1}{\sqrt{6}}\tilde{Z}
\left(1 + \frac{1}{6}\frac{2\tilde{\lambda}^{2}_{1}\lambda_{0}+3\tilde{\lambda}_{2}}{\tilde{\lambda}_1} \tilde{Z} + \dots\right).
\end{equation}
%
Therefore, to leading order, the center manifold yields
%
\begin{equation}
\frac{d\tilde{Z}}{d\tau} = -\lambda_{0}\tilde{\lambda}_{1}\tilde{Z}^2 + \dots \quad \Rightarrow \quad
\tilde{Z} = \frac{1}{\lambda_{0}\tilde{\lambda}_{1}\tau},\quad\text{as}\quad \tau\rightarrow +\infty,
\end{equation}


The $[1/1]$ Pad\'e is then

\begin{equation}
[1/1]_{\Sigma_{\dagger}}=\frac{\sqrt{6}\tilde{\lambda}^{2}_{1}
\tilde{Z}}{6\tilde{\lambda}_1-(2\tilde{\lambda}^{2}_{1}\lambda_{0}+3\tilde{\lambda}_2)\tilde{Z}}
\end{equation}


%
% which via the definition of $\tilde{Z}$~\eqref{taudef} gives
% %
% \begin{equation}
% H =% \left(\sqrt{\frac{3}{\alpha c^2}}\right)^{\frac{\alpha}{2}}\tau^{-\frac{\alpha}{4}},\qquad\text{as}\quad \tau\rightarrow +\infty.
% \end{equation}
% %
% Since $dt/d\tau = dt/d\ln{a} = H^{-1}$, it follows that
% %
% \begin{equation}\label{attractorresults}
% \ln{a} \propto %t^{\frac{1}{1+\frac{\alpha}{4}}},
% \qquad
% H\propto
% %t^{-\frac{\alpha}{4\left(1+\frac{\alpha}{4}\right)}},
% \qquad\text{as}\quad t\rightarrow +\infty.
% \end{equation}





%-----------------------------------------------------------------
\subsubsection*{The $\Omega_V=0$ subset}
%-----------------------------------------------------------------

In the case of the massless scalar field subset, which describes a massless
scalar field and a perfect fluid with a linear equation of state, one can
solve the Gauss constraint $1-\Sigma_\dagger^2 - \Omega_m$ globally for
$\Omega_m$, i.e., $\Omega_m = 1-\Sigma_\dagger^2$, which yields the dynamical
system

%
\begin{subequations}\label{dynsysm}
\begin{align}
\Sigma_\dagger^\prime &= -\frac32(2 - \gamma_m)(1-\Sigma_\dagger^2)\Sigma_\dagger,\label{Sigeqm}\\
Z^\prime &= \sqrt{6}\lambda_{0} Z(1-Z)\Sigma_\dagger,\label{Zeqm}
\end{align}
\end{subequations}
%
In this case the equation for $\Sigma_\dagger$ decouples and can be solved
separately. Moreover, it is invariant under the transformation
$\Sigma_\dagger \rightarrow -\Sigma_\dagger$, which allows one to identify
this problem with that of a (non-interacting) fluid and a stiff fluid,
characterized by $\Omega_\mathrm{stiff} = \Sigma_\dagger^2$, with a known
simple solution obtained via $\rho_\mathrm{stiff} \propto a^{-6}$ and the
Gauss constraint~\eqref{Gauss1} and the definitions:
%
\begin{equation}
\Omega_\mathrm{stiff} = \left(1+r\exp[3(2-\gamma_m)\tau]\right)^{-1}\quad \Rightarrow \quad
\Sigma_\dagger = \pm\left(1+r\exp[3(2-\gamma_m)\tau]\right)^{-1/2},
\end{equation}
%
where $r$ is an integration constant describing the ratio of the perfect
fluid and stiff fluid energy densities at a reference time $\tau=0$. Note
also that the special case $\Sigma_\dagger =0$ corresponds to the line of
fixed points $\mathrm{FL}_Z$ with $\Omega_m=1$. It follows that
$\Sigma_\dagger \rightarrow \pm 1$ ($\Sigma_\dagger \rightarrow 0$) when
$\tau \rightarrow -\infty$ ($\tau \rightarrow +\infty$).

If one wants the full solution on this boundary subset, it is easier to solve
the equation $\phi^\prime = \sqrt{6}\Sigma_\dagger$, which leads to that
$\phi - \phi_0 \propto \mathrm{arccoth}(|\Sigma_\dagger|^{-1/2})$ and hence
$\phi \rightarrow \phi_0$ when $\tau \rightarrow \infty$ and thus $Z
\rightarrow Z_0$, in agreement with the transversally hyperbolic nature of
$\mathrm{FL}_Z$.

On this subset $\mathrm{M}^{0}_{+}$, $\mathrm{M}^{1}_{-}$ are sources, and
$\mathrm{M}^{0}_{-}$ and $\mathrm{M}^{1}_{+}$ are saddles.

The solution space on this subset is therefore easily found and depicted in
Fig.~\ref{fig:V0}.


\begin{figure}[ht!]
\begin{center}
\subfigure[$\gamma_{m}=1$ and $\lambda_{0}=\sqrt{\frac{3}{2}}$]{\label{fig:V0Dust}
\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{V0_Boundary_Final.pdf}} \qquad
\end{center}
\caption{Solutions on the invariant $\Omega_V=0$ boundary subset of $\bar{\bf S}$.
}
\label{fig:V0}
\end{figure}









%-----------------------------------------------------------------
\subsubsection*{The $Z=1$ subset}
%-----------------------------------------------------------------


%
\begin{subequations}\label{dynsys}
\begin{align}
\Sigma_\dagger^\prime &= -(2-q)\Sigma_\dagger,\label{Sigeq1}\\
\Omega_m^\prime &= 3\left[2\Sigma_\dagger^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\\
\end{align}
\end{subequations}
%
where $2-q=3(1-\Sigma_\dagger^2) - \frac32 \gamma_m\Omega_m$. This system is
invariant under the transformation $\Sigma_\dagger \rightarrow
-\Sigma_\dagger$, which allows one to identify this problem with that of a
stiff fluid, characterized by $\Omega_\mathrm{stiff} = \Sigma_\dagger^2$, a
cosmological constant and a perfect fluid, which leads to that
%
\begin{equation}
3H^2 = \rho_{s0}\exp(-6\tau) + \Lambda + \rho_{m0}\exp(-3(2-\gamma_m)\tau),
\end{equation}
%
which yields
%
\begin{subequations}
\begin{align}
\Sigma_\dagger &= \pm \sqrt{\rho_{s0}}\exp(-6\tau)/\sqrt{3}H, \\
\Omega_m &= \rho_{s0}\exp(-3(2-\gamma_m)\tau)/3H^2, \\
\frac{\Sigma_\dagger^2}{\Omega_m} &\propto \exp(3(2-\gamma_m)\tau).
\end{align}
\end{subequations}
%

On this boundary the fixed points $\mathrm{M}^{1}_{\pm}$ are sources,
$\mathrm{dS}$ is a sink, and $\mathrm{FL}^1$ a saddle.

The solution space on this subset is depicted in Fig.~\ref{fig:Z1}


\begin{figure}[ht!]
\begin{center}
\subfigure[$\gamma_{m}=1$]{\label{fig:Z1Dust}
\includegraphics[width=0.45\textwidth, trim = 0cm 3cm 0cm 0cm]{Z1_Boundary_Final.pdf}} \qquad
\end{center}
\caption{Solutions on the invariant $Z=1$ boundary subset of $\bar{\bf S}$.
}
\label{fig:Z1}
\end{figure}


\section{Global analysis}

Since on $\Sigma_{\dagger}>0$ ($\Sigma_{\dagger}<0$), Z is monotonically
increasing (decreasing), it follows that in the full state space
$\mathrm{M}^{0}_{+}$ and $\mathrm{M}^{1}_{-}$ are sources,
$\mathrm{M}^{0}_{-}$ and $\mathrm{M}^{1}_{+}$ are saddles, and $\mathrm{PL}$
is a saddle.

Therefore interior orbits either originate from $\mathrm{M}^{0}_{+}$,
$\mathrm{M}^{1}_{-}$ (2-parameter set of solutions) or from $\mathrm{FL}_{Z}$
(one-parameter set of solutions).

Moreover since at $\Sigma_{\dagger}=0$, $\Sigma^{\prime}_{\dagger}>0$, the
flow can only cross the the positive $\Sigma$ direction, this is
$\Sigma_{\dagger}$ acts as a semipermeable membrane. From this it then
follows that $\mathrm{dS}$ is a sink in the full state space. ...

\section{Picard analysis of $\mathrm{FL}_{Z}$}

First we Taylor expand $\lambda(Z)$ around $Z_0$
\begin{equation}
 \lambda(Z)=\lambda(Z_0)+\sum_{k=1}\frac{1}{k!}\frac{d^{k}\lambda}{dZ^{k}}|_{Z=Z_0}(Z-Z_0)^k = \lambda_0 +\lambda_1 (Z-Z_0)+...
\end{equation}

A linear stability analysis of the line of fixed points FL$_{Z}$ shows that
it is transversally hyperbolic, with eigenvalues
$-3\left(1-\frac{\gamma_m}{2}\right), 3\gamma_m, 0$. The eigenvector
associated with the unstable direction is

% and corresponding eigenvectors given by
% \begin{equation*}
%  w_{1}=-3\left(1-\frac{\gamma_m}{2}\right)\qquad,\qquad
%  v^{(1)}_{1}=\begin{pmatrix}
%        -\sqrt{\frac{3}{2}}\left(1-\frac{\gamma_m}{2}\right) \\
%        0 \\
%        \lambda_{0}Z_{0}\left(1-Z_0\right)
%        \end{pmatrix}
%\end{equation*}

\begin{equation*}
% w_{2}=3\gamma_m \qquad,\qquad
 \mathbf{u}^{(1)}=\left(\sqrt{\frac{3}{2}}\gamma_m, -3\frac{\gamma_m \left(1+\frac{\gamma_m}{2}\right)}{\lambda_{0}},
 \lambda_{0}Z_{0}\left(1-Z_0\right)\right)
\end{equation*}

% \begin{equation*}
%  w_{3}=0 \qquad,\qquad
%  v^{(1)}_{3}=\begin{pmatrix}
%              0 \\
%              0 \\
%              1
%              \end{pmatrix}
% \end{equation*}

Expanding to higher order

\begin{equation}
\mathbf{u}%=(\Sigma_{\dagger},\Omega_m, Z)
=(0,1,Z_0)+
\left(\sqrt{\frac{3}{2}}\gamma_m, -\frac{3\gamma_m \left(1+\frac{\gamma_m}{2}\right)}{\lambda_{0}},\lambda_{0}Z_{0}\left(1-Z_0\right)\right)T+
(c_1, c_2, c_3)T^2+\dots
\end{equation}
where for convenience we introduced $T=e^{3\gamma_m \tau}$. Inserting the
above expansion in the full system, we obtain
\begin{equation}
 c_1=\sqrt{\frac{3}{2}}\gamma_m \left(1+\frac{\gamma_m}{2}\right)\frac{A}{B}
 %\left(\lambda_1 \lambda_{0} Z_0 (1-Z_0)-3\gamma_{m}(1+\gamma_m)-\lambda^{2}_{0}\right)
\end{equation}
\begin{equation}
 c_2=3\frac{\gamma_m}{\lambda^{2}_{0}}\left(3\gamma_m \left(1+\frac{\gamma_m}{2}\right)^2+\lambda^{2}_{0}\right)
\end{equation}
\begin{equation}
 c_3 =-\lambda_0 C \left(Z0-\frac{1}{2}\right) +\frac{(1+\frac{\gamma_m}{2})}{2}\frac{AC}{B}
\end{equation}
where
\begin{equation}
 A=\left(\lambda_1 C-3\gamma_m (1+\gamma_m)-\lambda^{2}_{0}\right)
\end{equation}
\begin{equation}
 B=\lambda_0 \left(1+\frac{3}{2}\gamma_m\right)
\end{equation}
\begin{equation}
 C=\lambda_0 Z_0 (1-Z_0)
\end{equation}


For the $[1/1]$ Pad\'e we obtain
\begin{equation}
 [1/1]_{\mathbf{u}}=\frac{\left(0,\quad 1,\quad Z_0\right)
 +\left(\sqrt{\frac{3}{2}}\gamma_m ,\quad \frac{\lambda_{0}}{\left(1+\frac{\gamma_m}{2}\right)},\quad C+\lambda_0 Z_0 \left(Z_0-\frac{1}{2}\right)-\frac{Z_0 A}{2B}\left(1+\frac{\gamma_m}{2}\right)\right) T}
 {(1,\quad 1,\quad 1)+\left(-\frac{A}{B}\left(1+\frac{\gamma_m}{2}\right), \quad \left(\frac{3\gamma_m \left(1+\frac{\gamma_m}{2}\right)^2+\lambda^{2}_{0}}{\lambda_0 \left(1+\frac{\gamma_m}{2}\right)}\right),\quad \lambda_0 \left(Z_0-\frac{1}{2}\right)-\frac{A}{2B}\left(1+\frac{\gamma_m}{2}\right)\right) T}
\end{equation}


\section{Examples}


\subsection{Constant Potential}

\begin{equation}
 V(\phi)=V_0>0
\end{equation}
then
\begin{equation}
 \lambda(\phi)=0
\end{equation}

\begin{figure}[ht!]
\begin{center}
\subfigure[Scalar field boundary $\Omega_m=0$]{\label{fig:M0CC}
\includegraphics[width=0.4\textwidth, trim = 0cm 3cm 0cm 0cm]{M0_Boundary_CC.pdf}} \qquad
\subfigure[$(\Sigma_\dagger,\Omega_{m},Z)$ state space.]{\label{fig:SSCC}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_CC.pdf}}
\subfigure[Projection  of $(\Sigma_\dagger,\Omega_{m},Z)$ state space into  $(\Sigma_\dagger,\Omega_{m})$.]{\label{fig:SSCCTop}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_CC_TopView.pdf}}
\end{center}
\caption{$V(\phi)=V_0>0$.
}
\label{fig:CC}
\end{figure}


\subsection{Exponential Potential}

\begin{equation}
 V(\phi)=V_0 e^{-\lambda_0 \phi}
\end{equation}
then
\begin{equation}
 \lambda(\phi)=\lambda_0
\end{equation}

\begin{figure}[ht!]
\begin{center}
\subfigure[Scalar field boundary $\Omega_m=0$]{\label{fig:M0EXP}
\includegraphics[width=0.4\textwidth, trim = 0cm 3cm 0cm 0cm]{M0_Boundary_EXP.pdf}} \qquad
\subfigure[$(\Sigma_\dagger,\Omega_{m},Z)$ state space.]{\label{fig:SSEXP}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_EXP.pdf}}
\subfigure[Projection  of $(\Sigma_\dagger,\Omega_{m},Z)$ state space into  $(\Sigma_\dagger,\Omega_{m})$.]{\label{fig:SSEXPTop}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_EXP_TopView.pdf}}
\end{center}
\caption{$\lambda_0=\sqrt{\frac{3}{2}}$.
}
\label{fig:EXP}
\end{figure}






\subsection{Example 3}

\begin{equation}
 V(\phi)=V_0\left(A e^{-\lambda_0 \phi}+1\right).
\end{equation}
Then
\begin{equation}
 \lambda(\tilde{Z})=\lambda_0\tilde{Z}
\end{equation}
in linear in $\tilde{Z}$, so that
\begin{equation}
 \tilde{\lambda}_{1}=\lambda_0=-\lambda_{1}
\end{equation}
and all high order $\tilde{\lambda}_{k}$ are identically zero.

In this case the center manifold of $\mathrm{dS}$ fixed point on the scalar
field $\Omega_m=0$ boundary is
\begin{equation}
 h(\tilde{Z})=\frac{\lambda^{2}_{0}}{3}\tilde{Z}^2
\end{equation}
\begin{equation}
 \Sigma_{\dagger}(\tilde{Z})=\frac{\lambda_{0}}{\sqrt{6}}\tilde{Z}\left(1+\frac{\lambda^{2}_{0}}{3}\tilde{Z}+\dots\right)
\end{equation}
and the corresponding $[1/1]$ Pad\'e reads
\begin{equation}
 [1,1]_{\Sigma_{\dagger}}=\sqrt{\frac{3}{2}}\frac{\lambda_{0}\tilde{Z}}{3-\lambda^{2}_{0}\tilde{Z}}.
\end{equation}


\begin{figure}[ht!]
\begin{center}
\subfigure[$V(\phi)$]{\label{fig:POTEXPCC}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Potential_EXPCC.pdf}} \qquad
\subfigure[$\lambda(\phi)$.]{\label{fig:LAMBDAEXPCC}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Lambda_EXPCC.pdf}}
\subfigure[Scalar field boundary $\Omega_m=0$]{\label{fig:M0EXPCC}
\includegraphics[width=0.4\textwidth, trim = 0cm 3cm 0cm 0cm]{M0_Boundary_EXPCC.pdf}} \qquad
\subfigure[$(\Sigma_\dagger,\Omega_{m},Z)$ state space.]{\label{fig:SSEXPCC}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_EXPCC.pdf}}
\subfigure[Projection  of $(\Sigma_\dagger,\Omega_{m},Z)$ state space into  $(\Sigma_\dagger,\Omega_{m})$.]{\label{fig:SSEXPCCTop}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_EXPCC_TopView.pdf}}
\end{center}
\caption{$\lambda_0=\sqrt{\frac{3}{2}}$ and $A=1/2$.
}
\label{fig:EXPCC}
\end{figure}





\subsection{Example 4}

\begin{equation}
 V(\phi)=V_0 \left(K_1 e^{-\lambda_0 \phi}+\frac{1}{\left(K_2 +\phi^2\right)^{\frac{\alpha}{2}}}\right)
\end{equation}
which as $\phi\rightarrow+\infty$ the inverse power law term dominates, the
global variable $Z$ is not suitable for a center manifold analysis since
$\tilde{\lambda}_1=0$. Instead we shall use a local variable
\begin{equation}
 Y=\frac{\alpha}{\phi}
\end{equation}
Then
\begin{subequations}
\begin{align}
\Sigma_\dagger^\prime &= \left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(Z)\right)(1-\Sigma_\dagger^2),\label{Sigsb}\\
 Y^{\prime} &=-\frac{\sqrt{6}}{\alpha}Y^2 \Sigma_{\dagger}.\label{Yeqsc}
\end{align}
\end{subequations}

We now Taylor expand $\lambda(Y)$ around $Y=0$
 \begin{equation}
 \lambda(Y)=\lambda(0)+\sum_{k=1}\frac{1}{k!}\frac{d^{k}\lambda}{dY^{k}}|_{Y=0} Y^k = 0 +\bar{\lambda}_1 Y+\frac{1}{2}\bar{\lambda}_{2}Y^2 +\cdots
 \end{equation}
A linearization yields
%
\begin{subequations}
\begin{align}
E^s &= \{(Y,\Sigma_{\dagger})|\,Y=0\},\\
E^c &= \{(Y,\Sigma_{\dagger})|\, \frac{\sqrt{6}}{\bar{\lambda}_{1}}\Sigma_{\dagger}-Y = 0\},
\end{align}
\end{subequations}
%
% To investigate the center manifold $W^c$ associated with the tangent space
% $E^c$, we adapt the variables to this tangent space and replace
% $\Sigma_\dagger$ with
% %
% \begin{equation}
%  v = \frac{\sqrt{6}}{\tilde{\lambda}_1} \Sigma_{\dagger} - \tilde{Z},
% \end{equation}
% %
% where $\mathrm{dS}$ is located $(\tilde{Z},v) = (0,0)$. This leads
% to the transformed system
% %
% \begin{subequations}
% \begin{align}
% \frac{d\tilde{Z}}{d\tau} &=-\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger \\
%                      %    &=-\tilde{\lambda}_{1}\lambda_{0} (1-\tilde{Z})\tilde{Z}(\tilde{Z}+v), \label{tildeZeq}\\
% \frac{dv}{d\tau} &=\frac{\sqrt{6}}{\tilde{\lambda}_1}\left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(Z)\right)(1-\Sigma_\dagger^2)+\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger, \\
%                %  &=\left[-3v+\frac{3}{\tilde{\lambda}_1}\left(\frac{1}{2}\tilde{\lambda}_2\tilde{Z}^2+\cdots\right)\right]\left(1-\frac{\tilde{\lambda}^{2}_{1}}{6}\left(\tilde{Z}+v\right)^2\right)
%                 %   +\tilde{\lambda}_{1}\lambda_{0} (1-\tilde{Z})\tilde{Z}(\tilde{Z}+v).\label{veq}
% \end{align}
% \end{subequations}
% %
% where $\Sigma_{\dagger}=\frac{\tilde{\lambda}_1}{\sqrt{6}}(\tilde{Z}+ v)$ The linearization of eq.~\eqref{veq} yields $\frac{d v}{d\tau}
% = - 3 v$, while eq.~\eqref{tildeZeq} only has higher order terms. The
% center manifold $W^c$ can be obtained as the graph $ v =
%  h(\tilde{Z})$ near $(\tilde{Z}, v) = (0,0)$ (i.e., use
% $\tilde{Z}$ as an independent variable), where $ h(0) =0$ (fixed point
% condition) and $\frac{d h}{d\tilde{Z}}(0) =0$ (tangency condition).
% Inserting this relationship into the above equations and using $\tilde{Z}$ as
% the independent variable leads to
%
%  \begin{equation}\label{tildehZ}
% -\sqrt{6}\lambda_{0} (1-\tilde{Z})\tilde{Z}\Sigma_\dagger\left(1+\frac{d h}{d\tilde{Z}}\right)
% -\frac{\sqrt{6}}{\tilde{\lambda}_1}\left(-3\Sigma_\dagger + \sqrt{\frac{3}{2}}\lambda(\tilde{Z})\right)(1-\Sigma_\dagger^2) = 0,
%  \end{equation}
%
%  where $\Sigma_{\dagger}=\frac{\tilde{\lambda}_1}{\sqrt{6}}(\tilde{Z}+ h(\tilde{Z}))$. Solving
%  this equation amounts to finding the center manifold solution, which
%  presumably is not possible. Instead we replace $h(\tilde{Z})$ with a
%  truncated power series in $\tilde{Z}$. Solving the
%  ensuing algebraic equations order by order yields
%
%  \begin{equation}\label{CMZ_Expansion}
% h(\tilde{Z}) = \frac{1}{6}\frac{2\tilde{\lambda}^{2}_{1}\lambda_{0}+3\tilde{\lambda}_{2}}{\tilde{\lambda}_1} \tilde{Z}^2 + \dots %\frac{\alpha^2}{108}\tilde{H}^4 + \dots
% \quad \Rightarrow \quad
% \Sigma_{\dagger} \approx \frac{\tilde{\lambda}_1}{\sqrt{6}}\tilde{Z}\left(1 + \frac{1}{6}\frac{2\tilde{\lambda}^{2}_{1}\lambda_{0}+3\tilde{\lambda}_{2}}{\tilde{\lambda}_1} \tilde{Z} + \dots\right).
% \end{equation}
% %
% Therefore, to leading order, the center manifold yields
% %
% \begin{equation}
% \frac{d\tilde{Z}}{d\tau} = -\lambda_{0}\tilde{\lambda}_{1}\tilde{Z}^2 + \dots \quad \Rightarrow \quad
% \tilde{Z} = \frac{1}{\lambda_{0}\tilde{\lambda}_{1}\tau},\quad\text{as}\quad \tau\rightarrow +\infty,
% \end{equation}
%
%
% The $[1/1]$ Pad\'e is then
%
% \begin{equation}
%  [1/1]_{\Sigma_{\dagger}}=\frac{\sqrt{6}\tilde{\lambda}^{2}_{1} \tilde{Z}}{6\tilde{\lambda}_1-(2\tilde{\lambda}^{2}_{1}\lambda_{0}+3\tilde{\lambda}_2)\tilde{Z}}
% \end{equation}
%


\begin{figure}[ht!]
\begin{center}
\subfigure[$V(\phi)$]{\label{fig:POTEXPINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Potential_EXPINV.pdf}} \qquad
\subfigure[$\lambda(\phi)$.]{\label{fig:LAMBDAEXPINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Lambda_EXPINV.pdf}}
\subfigure[Scalar field boundary $\Omega_m=0$]{\label{fig:M0EXPINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{M0_Boundary_EXPINV.pdf}} \qquad
\subfigure[$(\Sigma_\dagger,\Omega_{m},Z)$ state space.]{\label{fig:SSEXPINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_EXPINV.pdf}}
\subfigure[Projection  of $(\Sigma_\dagger,\Omega_{m},Z)$ state space into  $(\Sigma_\dagger,\Omega_{m})$.]{\label{fig:SSEXPINVTop}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_EXPINV_TopView.pdf}}
\end{center}
\caption{$\lambda_0=\sqrt{\frac{3}{2}}$, $\alpha=2$, and $K_1=1/100$, $K_2=100$.
}
\label{fig:EXPINV}
\end{figure}











\subsection{Example 5}

\begin{equation}
 V(\phi)=V_0 \left(K_1\left(\left(\phi-K_2\right)^2+K_3\right)^{\frac{\beta}{2}}\frac{e^{\lambda_0\phi}}{1+e^{\lambda_0\phi}}+\frac{1}{\left(K_4+\phi^2\right)^{\frac{\alpha}{2}}}\right)
\end{equation}



\begin{figure}[ht!]
\begin{center}
\subfigure[$V(\phi)$]{\label{fig:POTMONINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Potential_MONINV.pdf}} \qquad
\subfigure[$\lambda(\phi)$.]{\label{fig:LAMBDAMONINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Lambda_MONINV.pdf}}
\subfigure[Scalar field boundary $\Omega_m=0$]{\label{fig:M0MONINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 3cm 0cm 0cm]{M0_Boundary_MONINV.pdf}} \qquad
\subfigure[$(\Sigma_\dagger,\Omega_{m},Z)$ state space.]{\label{fig:SSMONINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_MONINV.pdf}}
\subfigure[Projection  of $(\Sigma_\dagger,\Omega_{m},Z)$ state space into  $(\Sigma_\dagger,\Omega_{m})$.]{\label{fig:SSMONINVTop}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 2cm 0cm]{SS_MONINV_TopView.pdf}}
\end{center}
\caption{$\lambda_0=\sqrt{\frac{3}{2}}$, $\alpha=2$, $\beta=1$, and $K_1=1/100$, $K_2=20$, $K_3=20$, $K_4=100$.
}
\label{fig:MONINV}
\end{figure}



%-----------------------------------------------------------------
\subsection{A potential with an exponential and constant term}
%-----------------------------------------------------------------

As before the matter content is dust, i.e., $\gamma_m=1$, while
%
\begin{equation}
V= V_0(\exp(-\bar{\lambda}\phi) + 1), \qquad \lambda =
\bar{\lambda}\left(\frac{\exp(-\bar{\lambda}\phi)}{1 + \exp(-\bar{\lambda}\phi)} \right)= \bar{\lambda}(1-Z),
\end{equation}
%
where we have made a translation in $\phi$ to set the constant potential part
to $V_0$, see Fig.~\ref{fig:lambdaexpL} for $\lambda(\phi)$, and where we
have used the definition~\eqref{Zdef} for $Z$.
%
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=0.4\textwidth]{Lambda_EXPCC.pdf} \vspace{-0.5cm}
\end{center}
\caption{$\lambda(\phi)$ for $V = V_0(\exp(-\bar{\lambda}\phi) + 1)$, exemplified by $\bar{\lambda}=\sqrt{\frac32}$.}
\label{fig:lambdaexpL}
\end{figure}
%

This leads to the dynamical system
%
\begin{subequations}\label{dynsysexpLambda}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}\bar{\lambda}(1-Z)(1 - x^2 - \Omega_m),\label{Sigexp}\\
\Omega_m^\prime &= 3\left[2x^2 - (1 - \Omega_m)\right]\Omega_m,\\
Z^\prime &= \sqrt{6}\bar{\lambda} Z(1-Z)x,
\end{align}
\end{subequations}
%
where again $2-q=3(1-x^2) - \frac32\Omega_m$.

The global dynamics is easily understood from the existence of the monotone
function $H^2$ in combination with the local properties of the fixed points
and the invariant boundary structure, where the $Z=0$ subset is described by
the reduced $x-\Omega_m$ dynamics for a potential with an exponential term
only, while the $Z=1$ subset is described by the reduced system for a
constant potential.

Due to the bifurcations on the `exponential' $Z=0$ subset we have three
cases: $0<\bar{\lambda}<\sqrt{3}$, $\sqrt{3}<\bar{\lambda}<\sqrt{6}$,
$\sqrt{6}<\bar{\lambda}$. As in the exponential potential case there is a
two-parameter set of solutions originating from $\mathrm{M}_0^+$ and
$\mathrm{M}_1^-$ when $\bar{\lambda} < \sqrt{6}$, while $\mathrm{M}_1^-$ is
the only source when $\bar{\lambda}>\sqrt{6}$, and a one-parameter set from
$\mathrm{FL}_Z$ into ${\bf S}$. All interior solutions, as well as the
solutions on the scalar field subset $\Omega_m=0$ end at the global attractor
$\mathrm{dS}_1$, as follows from the monotone function $H^2$, in combination
with the boundary structure and the local properties of the fixed points.

In the first case with $0<\bar{\lambda}<\sqrt{3}$ the attracting invariant
saddle surface that originates from $\mathrm{FL}_Z$ is a separatrix surface,
see Fig.~\ref{fig:SS_EXPCC1}; in the second case with
$\sqrt{3}<\bar{\lambda}<\sqrt{6}$ the saddle subset is bounded by the
trajectory from $\mathrm{EM}_0$ to $\mathrm{dS}_1$, which in also the
boundary of the subset that originates from the saddle $\mathrm{PL}_0$, which
together form a separatrix surface, see Fig.~\ref{fig:SS_EXPCC2}; in the
third case with $\sqrt{6}<\bar{\lambda}$ there is no separatrix surface, see
Fig.~\ref{fig:SS_EXPCC3}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\bar{\lambda}=\sqrt{3/2}$]{\label{fig:SS_EXPCC1}
\includegraphics[width=0.46\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXPCC.pdf}} \,
\subfigure[$\bar{\lambda}=2$]{\label{fig:SS_EXPCC2}
\includegraphics[width=0.46\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXPCC_lambda2.pdf}} \\ \vspace{-0.2cm}
\subfigure[$\bar{\lambda}=3$]{\label{fig:SS_EXPCC3}
\includegraphics[width=0.46\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXPCC_lambda3.pdf}} \vspace{-0.5cm} \,
\subfigure[Projection  of $(x,\Omega_{m},Z)$ state space onto  $(x,\Omega_{m})$ for $\bar{\lambda}=3$.]{\label{fig:SSEXPCCTop_lambda3}
\includegraphics[width=0.46\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_EXPCC_lambda3_TopView.pdf}}
\end{center}
\vspace{-0.1cm}
\caption{Representative state space pictures for models with dust and a scalar field potential
$V=V_0(\exp(-\bar{\lambda}\,\phi) + 1)$, where the three cases $0<\bar{\lambda}<\sqrt{3}$,
$\sqrt{3}<\bar{\lambda}<\sqrt{6}$, $\sqrt{6}<\bar{\lambda}$ have been represented by the values $\bar{\lambda}= \sqrt{3/2}$, $\bar{\lambda}=2$ and
$\bar{\lambda}=3$, showing the different character of the invariant subset that originates from the line of
fixed points $\mathrm{FL}_Z$.}
\label{fig:EXPCC}
\end{figure}
%
Note that the heteroclinic orbit from $\mathrm{EM}_0$ to $\mathrm{dS}_1$
forms a part of the boundary of the invariant subset that originates from
$\mathrm{FL}_Z$ when $\bar{\lambda} > \sqrt{6}$; moreover, when $\lambda$ is
large $\mathrm{EM}_0$ is a focus, which leads to that solutions with initial
data close to $\mathrm{FL}_0$ spiral around it. The projection of the
dynamics given in Fig.~\ref{fig:SSEXPCCTop_lambda3} shows that the
heteroclinic orbit $\mathrm{EM}_0 \rightarrow \mathrm{dS}_1$, and thereby
also solutions with initial small $Z$ and an $\Omega_m$ close to one, are
quite far from the $x=0$-plane, and as a consequence their evolution is also
substantially different from that of $\Lambda$CDM cosmology. To obtain models
of the present type with an observationally competitive history requires a
quite small $\bar{\lambda}$ (for such models a measure of deviations from
$\Lambda$CDM dynamics is obtained by considering the location of the
$\mathrm{PL}_0$ fixed point, since this point then gives the maximum
deviation of the attracting separatrix surface from the `$\Lambda$CDM
$x=0$-plane'). Finally, note that the behavior at late times of all the
solutions on the attracting separatrix surface, and nearby solutions, is
dominated by the one-dimensional stable center manifold of $\mathrm{dS}_1$ on
${\bf S}_\phi$, which therefore plays an important role for behavior at late
times; a situation that is typical for the remaining examples.

%Next we consider the figures for the observables in~\ref{fig:obsEXPCC}, which
%again shows that increasing shadowing of the $\Lambda$CDM case occurs for
%initial values with large $Z$, while decreasing $Z$ yields an increasing
%deformation of the $\Lambda$CDM evolution. Nevertheless, it can be seen that
%models with a potential that consists of an exponential and a constant lead
%to a separaterix surface that describes solutions that continuously deform
%the $\Lambda$CDM evolution, and by requiring initial data for which
%$\Omega_m$ is close to one leads to an open set of solutions that exhibits
%this property.
%
%\begin{figure}[ht!]
%\begin{center}
%\subfigure[$\Omega_m-q$]{\label{fig:}
%\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{q_M_EXPCC.pdf}} \qquad
%\subfigure[$H-q$]{\label{fig:}
%\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{h_q_EXPCC.pdf}} \\
%\subfigure[$H-j$]{\label{fig:}
%\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{h_j_EXPCC.pdf}}\qquad
%\subfigure[$q-j$]{\label{fig:}
%\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{j_q_EXPCC.pdf}} \vspace{-0.5cm}
%\end{center}
%\caption{Observables for $V=V_0(\exp(-\sqrt{3/2}\,\phi) + 1)$.}
%\label{fig:obsEXPCC}
%\end{figure}
%

Letting $\bar{\lambda} \rightarrow \infty$ for the present models takes us
naturally to the next very popular model, that of dust with an inverse
power-law scalar field potential, which, however, was only recently studied
in the context of a global dynamical system~\cite{tam14}. Here we will add
some global and heuristic aspects to that analysis, as well as situating
these models in the context of almost $\Lambda$CDM cosmology.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Models with global bounds on $\lambda$}\label{sec:bound}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The present class of models admit cases for which $\lambda(\phi)$ is
non-monotonic, indeed, the models with a monotonically decreasing potential
for which $\lim_{\phi\rightarrow \pm\infty} \lambda =0$ by definition have a
non-monotone $\lambda$. A $\lambda$-based scalar field variable $Z$ such as
that in the inverse power-law example is thereby not appropriate. An optimal
choice of $Z$ depends on the type of potential, especially its asymptotic
properties. Below we will consider two examples of potentials, one with
$\lim_{\phi\rightarrow -\infty} \lambda = \bar{\lambda}
>0$ and one with $\lim_{\phi\rightarrow -\infty} \lambda = 0$, which will
exemplify general features of these types of models. In both cases we will
use
%
\begin{equation}\label{Zdef2}
Z=\frac{1}{1+\exp(-\bar{\lambda}\phi)},
\end{equation}
%
as the $Z$ variable. This is by no means an optimal choice of scalar field
variable for certain analytical consideration for all models with a globally
bounded $\lambda$. One can even question if there in general exists a
globally and analytically most convenient variable since the two asymptotic
regimes, $\phi\rightarrow + \infty$ and $\phi\rightarrow -\infty$, and the
intermediate behavior, in general require a different type of variable.
Nevertheless, for the present models and our main purpose, which is to obtain
a picture of continuous and global $\Lambda$CDM deformations for quite
general potentials, the above choice of $Z$ is `good enough.'


%-----------------------------------------------------------------
\subsection{Models with $\lim_{\phi\rightarrow -\infty} \lambda = \bar{\lambda}
>0$}
%-----------------------------------------------------------------

To illustrate the properties of these models we give a simple example with a
potential of the type
%
\begin{equation}
V(\phi)=V_0 \left(K_1 e^{-\bar{\lambda} \phi}+\frac{1}{\left(K_2 +\phi^2\right)^{\frac{\alpha}{2}}}\right),
\end{equation}
%
where the positive constants $\bar{\lambda}$, $K_1$, $K_2$, and $\alpha$ are
to be chosen so that the potential is monotonically decreasing. In this case
the potential is asymptotically described by the exponential term when $\phi
\rightarrow - \infty$ and by an inverse power law behavior $\phi^{-\alpha}$
when $\phi \rightarrow + \infty$. A specific example is obtained by setting
$\bar{\lambda}=\sqrt{\frac{3}{2}}$, $K_1=1/100$, $K_2=100$, $\alpha=2$, where
the values have been partly chosen for numerical convenience, see
Fig.~\ref{fig:EXPINVPot}; note that $\lambda$ is not monotonically
decreasing.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$V(\phi)$]{\label{fig:POTEXPINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Potential_EXPINV_LogScale.pdf}} \qquad
\subfigure[$\lambda(\phi)$.]{\label{fig:LAMBDAEXPINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Lambda_EXPINV.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{Features of the potential
$V(\phi)=V_0 \left(K_1 e^{-\bar{\lambda} \phi}+ (K_2 +\phi^2)^{-\alpha/2}\right)$,
with $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $\alpha=2$, and $K_1=1/100$, $K_2=100$.}
\label{fig:EXPINVPot}
\end{figure}
%

These models give rise to the solution space depicted in
Fig.~\ref{fig:EXPINVstate}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[The scalar field boundary $\Omega_m=0$]{\label{fig:M0EXPINV}
\includegraphics[width=0.42\textwidth, trim = 0cm 2.5cm 0cm 0cm]{M0_Boundary_EXPINV.pdf}} \quad
\subfigure[The $(x,\Omega_{m},Z)$ state space.]{\label{fig:SSEXPINV}
\includegraphics[width=0.52\textwidth, trim = 0cm 2.5cm 0cm 0cm]{SS_EXPINV.pdf}} \vspace{-0.3cm}
%\\
%\subfigure[Projection of the $(x,\Omega_{m},Z)$ state space onto the $(x,\Omega_{m})$-plane.]{\label{fig:SSEXPINVTop}
%\includegraphics[width=0.50\textwidth, trim = 0cm 2.5cm 0cm 0cm]{SS_EXPINV_TopView.pdf}}   \vspace{-0.5cm}
\end{center}
\caption{State space structures, including the deformed $\Lambda$CDM separatrix surface with solutions
originating from $\mathrm{FL}_Z$, for dust and a potential
$V(\phi)=V_0 \left(K_1 e^{-\bar{\lambda} \phi}+ (K_2 +\phi^2)^{-\alpha/2}\right)$
with $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $\alpha=2$, and $K_1=1/100$, $K_2=100$.}
\label{fig:EXPINVstate}
\end{figure}
%
The subset $Z=0$ have the same equations as those on the reduced $x-\Omega_m$
state space for the exponential potential, which in the present case have the
solutions given in eq.~\eqref{expsol}, while the subset $Z=1$ is described by
the reduced solvable equations for a constant potential. As can be seen the
structure of the deformed $\Lambda$CDM separatrix surface is qualitatively
similar to that for a potential consisting of an exponential and a non-zero
constant. The difference is more in the physical situation. For the present
models the dark energy is decreasing toward its `natural' value zero,
although at a slower rate than the energy density  of the matter content.
%The dynamical similarities with the case of a potential consisting of an
%exponential and a non-zero constant lead to qualitatively similar pictures
In order to physically compare the solutions with the $\Lambda$CDM models,
they are depicted in the $q-\Omega_m$ and cosmographic diagrams in
Fig.~\ref{fig:obsEXPCC}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\Omega_m-q$ diagram]{\label{fig:q_M_EXPINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{q_M_EXPINV.pdf}} \qquad
\subfigure[$H-q$ diagram]{\label{fig:h_q_EXPINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{h_q_EXPINV.pdf}} \\
\subfigure[$H-j$ diagram]{\label{fig:h_j_EXPINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{h_j_EXPINV.pdf}}\qquad
\subfigure[$q-j$ diagram]{\label{fig:j_q_EXPINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{j_q_EXPINV.pdf}} \vspace{-0.5cm}
\end{center}
\caption{Observables for models with dust and $V(\phi)=V_0 \left(K_1 e^{-\bar{\lambda} \phi}+ (K_2 +\phi^2)^{-\alpha/2}\right)$
with $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $\alpha=2$, and $K_1=1/100$, $K_2=100$. The full line represents $\Lambda$CDM cosmology,
the horizontal line represents now, i.e, zero redshift, while the other lines correspond to the solution in Fig.~\ref{fig:EXPINVstate}.}
\label{fig:obsEXPCC}
\end{figure}
%
Most models with $\Omega_m$ close to one initially quite closely resemble
those of $\Lambda$CDM cosmology; it is only those with a very small initial
value of $Z$ that deviate rather significantly. Since the fixed point
$\mathrm{PL}_0$ gives a measure of maximum deviation from the $\Lambda$CDM
$x=0$ surface, it follows that this situation can be improved by decreasing
$\bar{\lambda}$, and the limit $\bar{\lambda}\rightarrow 0$ takes us to the
next set of models.

%-----------------------------------------------------------------
\subsection{Models with $\lim_{\phi\rightarrow -\infty} \lambda =0$}
%-----------------------------------------------------------------

To illustrate the properties of these models we consider the following
example:
%
\begin{equation}
V(\phi)= V_0 \left[K_1\left(\left(\phi-K_2\right)^2+K_3\right)^{\frac{\beta}{2}}
\left(\frac{\exp(-\bar{\lambda}\phi)}{1+\exp(-\bar{\lambda}\phi)}\right)+\frac{1}{\left(K_4+\phi^2\right)^{\alpha/2}}\right],
\end{equation}
%
where the positive constants $\bar{\lambda}$, $K_1$, $K_2$, $K_3$, $K_4$,
$\alpha$, and $\beta$ are to be chosen so that the potential is monotonically
decreasing. The potential is dominated by a power-law $\phi^\beta$ when $\phi
\rightarrow - \infty$ and by an inverse power-law behavior $\phi^{-\alpha}$
when $\phi \rightarrow + \infty$. An example of such a potential is obtained
by setting $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $K_1=1/200$, $K_2=1$, $K_3=1$,
$K_4=50$, $\alpha=2$, and $\beta=1$, where the constant have been partly
chosen for numerical convenience, and is illustrated in
Fig.~\ref{fig:MONINVpot}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$V(\phi)$]{\label{fig:Potential_MONINV_v2}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Potential_MONINV_v2.pdf}} \qquad
\subfigure[$\lambda(\phi)$.]{\label{fig:LAMBDAMONINV}
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0cm]{Lambda_MONINV_v2.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{Properties of the scalar field potential
$V(\phi)= V_0 \left(K_1\left(\left(\phi-K_2\right)^2+K_3\right)^{\frac{\beta}{2}}
\left(\frac{\exp(-\bar{\lambda}\phi)}{1+\exp(-\bar{\lambda}\phi)}\right)+\frac{1}{\left(K_4+\phi^2\right)^{\alpha/2}}\right)$
with $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $K_1=1/200$, $K_2=1$, $K_3=1$,
$K_4=50$, $\alpha=2$, and $\beta=1$}
\label{fig:MONINVpot}
\end{figure}
%
These models give rise to the state space features depicted in
Fig.~\ref{fig:MONINVstate}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[Scalar field boundary $\Omega_m=0$]{\label{fig:M0_MONINV_v2}
\includegraphics[width=0.4\textwidth, trim = 0cm 2cm 0cm 0cm]{M0_MONINV_v2.pdf}} \quad
\subfigure[The $(x,\Omega_{m},Z)$ state space.]{\label{fig:SSMONINV_v2}
\includegraphics[width=0.45\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_MONINV_v2.pdf}}\\
\subfigure[Projection of the $(x,\Omega_{m},Z)$ state space onto the $(x,\Omega_{m})$-plane.]{\label{fig:SSMONINVTop}
\includegraphics[width=0.45\textwidth, trim = 0cm 2cm 0cm 0cm]{SS_MONINV_TopView_v2.pdf}}
\vspace{-0.5cm}
\end{center}
\caption{State space features for models with dust and a scalar field with the potential
$V(\phi)= V_0 \left(K_1\left(\left(\phi-K_2\right)^2+K_3\right)^{\frac{\beta}{2}}
\left(\frac{\exp(-\bar{\lambda}\phi)}{1+\exp(-\bar{\lambda}\phi)}\right)+\frac{1}{\left(K_4+\phi^2\right)^{\alpha/2}}\right)$ with
$\lambda_0=\sqrt{\frac{3}{2}}$, $\alpha=2$, $\beta=1$, and $K_1=1/200$, $K_2=1$, $K_3=1$,
$K_4=50$.}
\label{fig:MONINVstate}
\end{figure}
%

In this case both the $Z=0$ and $Z=1$ subsets have equations that are the
same as those of the reduced system for a constant potential. The deformed
$\Lambda$CDM separatrix surface is bounded by two $\Lambda$CDM orbits at
$Z=0$ and $Z=1$, the line of fixed points $\mathrm{FL}_Z$ and the
heteroclinic separatrix orbit on ${\bf S}_\phi$ that originates from
$\mathrm{dS}_0$ and ends at the future attractor $\mathrm{dS}_1$.
% (note the
%`attracting' property of the center manifold of $\mathrm{dS}_1$ on ${\bf
%S}_\phi$). The deformation of the $\Lambda$CDM models thereby take place
%continuously when $0<Z<1$ for models with this kind of deformed separatrix
%surface, where models with $Z\rightarrow 0\, \text{or}\,1$ initially exhibit
%$\Lambda$CDM dynamics.
The properties of these solutions in terms of physical observables are
illustrated in Figs.~\ref{fig:obsMONINV} and~\ref{fig:obsMONINV2}.
%
\begin{figure}[ht!]
\begin{center}
\subfigure[$\Omega_m-q$ diagram]{\label{fig:q_M_MONINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{q_M_MONINV_v2.pdf}} \qquad
\subfigure[$H-q$ diagram]{\label{fig:h_q_MONINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{h_q_MONINV_v2.pdf}} \\
\subfigure[$H-j$ diagram]{\label{fig:h_j_MONINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{h_j_MONINV_v2.pdf}}\qquad
\subfigure[$q-j$ diagram]{\label{fig:j_q_MONINV}
\includegraphics[width=0.45\textwidth, trim = 0cm 0.8cm 0cm 0cm]{j_q_MONINV_v2.pdf}} \vspace{-0.5cm}
\end{center}
\caption{Observables for solutions with dust and a scalar field potential $V(\phi)= V_0 \left(K_1\left(\left(\phi-K_2\right)^2+K_3\right)^{\frac{\beta}{2}}
\left(\frac{\exp(-\bar{\lambda}\phi)}{1+\exp(-\bar{\lambda}\phi)}\right)+\frac{1}{\left(K_4+\phi^2\right)^{\alpha/2}}\right)$
with $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $\alpha=2$, $\beta=1$, and $K_1=1/100$, $K_2=20$, $K_3=20$, $K_4=100$, as matter content.}
\label{fig:obsMONINV}
\end{figure}
%
%
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=0.5\textwidth]{h_q_j_MONINV_v2.pdf} \vspace{-0.5cm}
\end{center}
\caption{Cosmographic diagram for $V(\phi)= V_0 \left(K_1\left(\left(\phi-K_2\right)^2+K_3\right)^{\frac{\beta}{2}}
\left(\frac{\exp(-\bar{\lambda}\phi)}{1+\exp(-\bar{\lambda}\phi)}\right)+\frac{1}{\left(K_4+\phi^2\right)^{\alpha/2}}\right)$
with $\bar{\lambda}=\sqrt{\frac{3}{2}}$, $\alpha=2$, $\beta=1$, and $K_1=1/100$, $K_2=20$, $K_3=20$, $K_4=100$.}
\label{fig:obsMONINV2}
\end{figure}
%

Due to the globally small $\lambda(Z)$, the present models have a temporal
evolution that is very similar to that of $\Lambda$CDM cosmology when
$\Omega_m$ is initially close to one. They illustrate that it is possible to
deform the $x=0$ $\Lambda$CDM surface continuously, thus obtaining models
that are observationally indistinguishable from $\Lambda$CDM cosmology,
\emph{and} to explore the observational bounds, thus forming a natural
phenomenological scalar field testing ground for $\Lambda$CDM cosmology.

We will deal with the case $V=\Lambda$, and hence $\lambda =0$, in the next
section, so let us now focus on the case where $V$ is strictly monotonically
decreasing and $\lambda(Z)$ is hence positive, except in the limits $Z_\pm$
where $\lambda$ might be zero.

%-----------------------------------------------------------------
\subsubsection*{Asymptotic past dynamics}
%-----------------------------------------------------------------

Due to that $\rho_m \propto \exp(-3\gamma_m\tau) \propto a^{-3\gamma_m}$ it
follows that $H^2\rightarrow \infty$ toward the past and hence that either
(i) $V\rightarrow \infty$ (and thereby $Z\rightarrow Z_-$, since the
potential is monotonically decreasing) or (ii) $\Omega_V = 1 - x^2 -
\Omega_m\rightarrow 0$, toward the past. However, case (i), which assumes a
potential for which $V$ becomes unbounded in the direction of decreasing
$\phi$, implies a special subcase of case (ii), i.e., all orbits in ${\bf S}$
originate from the $\Omega_V=0$ subset, as can be seen as follows. It follows
from~\eqref{Zeq0} that $x$ must have been negative (and a priori possibly
asymptotically zero) toward the past (this can also be understood by noticing
that a `scalar field particle' must move toward decreasing $\phi$ in order
for $V \rightarrow \infty$ for a monotonically decreasing potential).
Equation~\eqref{xeq} shows that $x$ is monotonically increasing when
$-1<x<0$, since then $2-q>0$ (recall that $q=2$ only when $x^2=1$), which
implies that if an orbit resides or partly resides in the region $-1<x<0$
then the orbit originates from $x=-1$, and hence the past asymptotic state is
a massless scalar field state described by $x=-1$, $\Omega_m=\Omega_V = 0$.
Since all orbits in ${\bf S}$ originates from the massless scalar field
boundary subset $\Omega_V=0$, this suggests that we take a closer look at
this subset and its neighborhood.

%-----------------------------------------------------------------
\subsubsection*{The massless scalar field boundary and its neighborhood}
%-----------------------------------------------------------------

We first note that the equations on $\Omega_V=0$ are invariant under the
transformation $x\rightarrow - x$, and therefore one can use
$\Omega_\mathrm{stiff}=x^2$ as a variable instead of $x$, which explicitly
leads to that the dynamics of the massless scalar field take the same form as
that of a stiff fluid, $p=\rho$. This leads to the equations
%
\begin{subequations}\label{dynsysOmV}
\begin{align}
\Omega_m^\prime &= 3(2-\gamma_m)(1 - \Omega_m)\Omega_m,\label{OmeqV}\\
Z^\prime &= \pm\sqrt{6}\,\frac{dZ}{d\phi}\,\sqrt{1-\Omega_m},\label{ZeqV}
\end{align}
\end{subequations}
%
where we have used the Gauss constraint to write $\Omega_\mathrm{stiff}= 1 -
\Omega_m$. Since the equation for $\Omega_m$ decouple from that of $Z$, it
follows that all orbits in the two-dimensional state space look the same as
the orbits of the one-dimensional problem associated with $\Omega_m$, when
projected onto $\Omega_m$. As follows from~\eqref{OmeqV}, $\Omega_m$
monotonically increases from $0$ to $1$. Thus the past asymptotic state of
all orbits on $\Omega_V=0$ resides on the pure massless scalar field subset
$\Omega_m=0$, $x^2=1$ and thereby $q=2$; either on $x=+1$ which leads to
$Z=Z_-$ and thus on a fixed point $\mathrm{M}^+_{Z_-}$, or $x=-1$ which
yields $Z=Z_+$ and a fixed point $\mathrm{M}^-_{Z_+}$. It also follows that
the future asymptotic state of all orbits on $\Omega_V=0$ resides on
$\Omega_m=1$ and $x=0$, and thereby $q=\frac12$ when $\gamma_m=1$. To
determine what happens with $Z$ toward the future on $\Omega_V=0$ we note
that it is possible to solve the equation for $\Omega_m$ explicitly, but the
solution can also be obtained by using that $\rho_\mathrm{stiff} \propto
\exp(-6\tau)$, $\rho_m \propto \exp(-3\tau)$, $3H^2 = \rho_\mathrm{stiff} +
\rho_m$, which gives
%
\begin{subequations}
\begin{align}
\frac{H}{H_0} &= \left(\Omega_{\mathrm{stiff}0}\exp(-6\tau)+ \Omega_{m0}\exp(-3\tau)\right)^{1/2},\\
x &= \pm \Omega_{\mathrm{stiff}0}^{1/2}\,\exp(-3\tau)/(H/H_0), \\
\Omega_m &= \Omega_{m0}\exp(-3\tau)/(H/H_0)^2.
%,\\
%\frac{x^2}{\Omega_m} &\propto \exp(-3\tau).
\end{align}
\end{subequations}
%
To obtain the solution for $Z$ requires a specification of $Z$. However, it
is easier to obtain $Z$ by integrating the equation for $\phi$,
%
\begin{equation}
\phi^\prime = \sqrt{6}x = \pm(6\Omega_{\mathrm{stiff}0})^{1/2}\,\exp(-3\tau)/(H/H_0),
\end{equation}
%
and then inserting the result in the chosen definition of $Z$. As follows
from the above equation, $\phi \rightarrow \mathrm{const.}$ when $\tau
\rightarrow \infty$, and therefore also $Z\rightarrow \mathrm{const.}$ in
this limit, i.e., in the full state space the $\Omega_m=1$ subset on
$\Omega_V=0$ corresponds to a line of fixed points $\mathrm{FL}_Z$ where each
point on the line attracts one solution with positive (negative) $x$ on the
$\Omega_V=0$ boundary.

To study the neighborhood of $\Omega_V=0$ we note that
%
\begin{equation}
\Omega_V^{-1}\Omega_V^\prime|_{\Omega_V=0} = 2\left(1 + q - \sqrt{\frac32}\lambda(Z)\, x\right),
\end{equation}
%
and hence
%
\begin{equation}
\Omega_V^{-1}\Omega_V^\prime|_{\mathrm{M}^\pm_{Z_\mp}} = \sqrt{6}\left(\sqrt{6} \mp \lambda(Z_\mp)\right), \qquad
\Omega_V^{-1}\Omega_V^\prime|_{\mathrm{FL}_Z} = 3.
\end{equation}
%
As a consequence it follows that $\mathrm{M}^-_{Z_+}$ is a source for orbits
in ${\bf S}$ as is $\mathrm{M}^+_{Z_-}$ when $\lambda(Z_-) < \sqrt{6}$, but
if $\lambda(Z_-) > \sqrt{6}$ then no orbits enter ${\bf S}$ from this fixed
point;\footnote{Solutions originating from $\mathrm{M}^-_{Z_+}$ corresponds
to `scalar field particles' that come from $\phi_+$ moving toward decreasing
$\phi$, while if $\lambda(Z_-) < \sqrt{6}$ the potential is shallow enough so
that scalar field particles can come from $\phi_-$ moving toward increasing
$\phi$. Solutions originating from $\mathrm{FL}_Z$ correspond to scalar field
particles that initially lie still at some initial value $\phi$ in the
potential $V(\phi)$.} a one-parameter set of orbits, one from each value of
$Z$, thus forming an invariant subset in ${\cal S}$, enters ${\bf S}$ from
$\mathrm{FL}_Z$, which acts as a `transversal saddle.' In combination with
that we demand initial data for which $\Omega_m$ is close to one, this
implies the important feature that for \emph{all} the models we consider such
solutions can be approximated by this (at least initially) `\emph{attracting
invariant} surface \emph{subset}' of co-dimension one with respect to the
state space dimension. We will later contrast this feature with discussions
about \emph{single} `attracting' solutions that appear for various models in
the literature.

%-----------------------------------------------------------------
\subsubsection*{Restrictions on asymptotic future dynamics}
%-----------------------------------------------------------------

Let us now return to the implication of the monotonicity of $H$, but for
future asymptotic dynamics. It follows from~\eqref{Hmon} that either (i)
$H=0$ or (ii) $1+q=0$ (and if so, $\dot{\phi}=V=\rho_m=0$, and thereby
$Z=Z_+$) toward the future, or both. In case (i) it follows that
$H\rightarrow 0$ and hence, due to the Gauss constraint, $V(Z)\rightarrow 0$,
which, since $V$ is monotonically decreasing leads to $Z\rightarrow Z_+$,
i.e. case (i) can only happen for potentials for which $\lim_{\phi
\rightarrow \infty}V=0$. Case (ii) automatically brings the scalar field
subset ${\bf S}_\phi$ into focus, since that it is only on this subset
$1+q=0$ due to that eqs.~\eqref{eqqineq} and~\eqref{Gauss2} yield
%
\begin{equation}
1+q=0 \quad \Rightarrow \quad \Omega_m=0=x, \quad \Omega_V=1,
\end{equation}
%
which describes an asymptotic (quasi) de Sitter state.

On ${\bf S}_\phi$ we have the inequalities $-1\leq q \leq 2$. As a
consequence $H^2$ is a monotonically decreasing function also on this subset,
except when $1+q=0$. However, when $1+q=0$ is not an invariant subset, then
$1+q=0$ just corresponds to an inflection point in the graph of $H^2$, and
$H^2$ is thereby still monotonically decreasing.\footnote{For an elaboration
of this in the case of monomial potentials, see~\cite{alhetal15}.} It is only
when $1+q=0$ is an invariant subset $H^2$ stops decreasing and for which
$1+q=0$ can be an asymptotic state. For this to happen requires that $x=0$,
since $1+q=3x^2$ on ${\bf S}_\phi$, which results in that also $x=0$ must be
an invariant subset. Since $x^\prime = \sqrt{\frac{3}{2}}\lambda(Z)$ when
$\Omega_m=0=x$, $\Omega_V=1$, this leads to the condition that $\lambda(Z)=0$
asymptotically toward the future, which in turn implies that $Z=Z_-$ or
$Z=Z_+$ asymptotically toward the future. Due to the monotonicity of $H^2$ it
follows that the first case is only possible if $\phi_-=-\infty$ and if
$\lim_{\phi\rightarrow -\infty}V = \mathrm{const}. < \infty$, while the
monotonicity of $H^2$ yields that $Z=Z_+$ asymptotically toward the future if
the potential is unbounded; in this paper we will only consider unbounded
monotonically decaying potentials (with the exception of when $V=\Lambda$,
which is treated separately) and hence $Z=Z_+$ asymptotically toward the
future. Note that this statement holds for both ${\bf S}$ and ${\bf S}_\phi$.

In case (i) $\lim_{Z\rightarrow Z_+}\lambda(Z)$ may or may not have a limit,
but for simplicity we will assume that $\lim_{Z\rightarrow Z_+}\lambda(Z) =
\lambda_+ \geq 0$ and in both case (i) and case (ii) we assume that
$\lambda(Z)$ is differentiable on $\bar{\bf S}$. Since the future attractor
must reside on the subset $Z_+$ it is determined by the dynamics on this
subset, which is given by the dynamical system:
%
\begin{subequations}\label{dynsysp}
\begin{align}
x^\prime &= -(2-q)x + \sqrt{\frac{3}{2}}\lambda_+\,\Omega_V,\label{xeqp}\\
\Omega_m^\prime &= 3\left[2x^2 - \gamma_m(1 - \Omega_m)\right]\Omega_m,\label{Omeqp}
\end{align}
\end{subequations}
%
where we recall that
%
$$
\Omega_V = 1 - x^2 - \Omega_m, \qquad q = - 1 + 3x^2 + \frac32\gamma_m\Omega_m.
$$
%



Further restrictions on the global dynamics of the present models come from
eqs.~\eqref{xeq} and~\eqref{Zeq0}. It follows from~\eqref{xeq} that $x$ is
monotonically increasing when $-1<x<0$, since then $2-q>0$. This in turn
implies that if an orbit resides or partly resides in the region $-1<x<0$
then the orbit originates at $x=-1$, and hence the past asymptotic state is a
massless scalar field state described by $x=-1$, $\Omega_m=\Omega_V = 0$.
Note the central position of the energy and causality induced inequalities
for $q$, which strongly suggests that the field equation should be expressed
in terms of $q$, as done in~\eqref{dynsys0}. In addition, the equation for
$x$ can be written on the following form
%
\begin{equation}
x^\prime = -\frac32\left(2-\gamma_m\right)x\,\Omega_m - 3\left(x - \frac{\lambda(Z)}{\sqrt{6}}\right)\Omega_V,
\end{equation}
%
which shows that $x$ is monotonically decreasing if $1>x \geq
\lambda(Z)/\sqrt{6}>0$. Moreover, since $Z$ is monotonically increasing in
$\phi$, it follows from~\eqref{Zeq0} that $Z$ is monotonically increasing
(decreasing) in $\tau$ when $x>0$ ($x<0$). If $x$ is future asymptotically
positive it follows that $Z\rightarrow Z_+$. Taking this in combination with
the above implies that asymptotically toward the future $0 \leq x\leq
\lim_{Z\rightarrow Z_+}\lambda(Z)/\sqrt{6}$, where $x\neq 0$ is only possible
toward the future if $\lim_{Z\rightarrow Z_+}\lambda(Z) \neq 0$. If on the
other hand $\lim_{Z\rightarrow Z_+}\lambda(Z) =0$, then $x=0$ toward the
future, and due to~\eqref{Omeq} this implies that $\Omega_m\rightarrow 0$ and
thus the future state is a (quasi) de Sitter state with $x=0=\Omega_m$,
$\Omega_V=1$.


%Since $0 \leq x\leq \lim_{Z\rightarrow Z_+}\lambda(Z)/\sqrt{6}$ toward the
%future it follows that asymptotically toward the








Let us now turn to the issue at hand, `dynamical almost $\Lambda$CDM
cosmology.' To do so we begin with some familiar introductory dynamical
systems examples, which allow us to introduce some ideas and conventions used
later on, and which motivate some general conditions on the scalar field
potential that lead to an open set of solutions that continuously deform
$\Lambda$CDM cosmology.
%Apart from some general considerations in the
%discussion, we throughout assume that the matter equation of state is that of
%dust, i.e., $\gamma_m=1$.

These monotonic features severely restrict the global dynamics and is
reflected in the local analysis of fixed points in subsequent sections.
Finally, it is possible to extract additional information about the global
and asymptotic dynamics when the potential is further restricted, especially
if e.g. the monotonicity of $H^2$ is combined with an analysis of the
dynamics on the boundaries, including a local fixed point analysis.


Although this is not always most convenient, see e.g.~\cite{ugg13}, in this
paper we will choose scalar field variables so that $Z_-=0$ and $Z_+ = 1$ so
that $\bar{\bf S}$ is given by
%
\begin{equation}
\Omega_m\geq 0, \qquad \Omega_\phi = x^2 + \Omega_V = 1 - \Omega_m \geq 0, \qquad \Omega_V \geq 0, \qquad 0\leq Z \leq 1.
\end{equation}
%



In this paper we have investigated flat FLRW models with a minimally coupled
scalar field and dust by means of global dynamical systems on an extended
compact three-dimensional state space. This framework has been used in order
to obtain models that exhibit an open set of solutions that almost behave as
$\Lambda$CDM cosmology, assuming that there is a mechanism that sets initial
conditions so that $\Omega_m$ is close to one initially. It was found that it
does not suffice to have a shallow tail of a monotonically decreasing scalar
field potential, but that it is essential to have global as well as
asymptotic bounds on $\lambda=-V_\phi/V$. In particular, if $\lim_{\phi
\rightarrow \pm \infty} \lambda < \sqrt{3}$ there exists an `attracting
separatrix surface,' which severely restricts the dynamics, indeed the
dynamics on this invariant subset approximately describes the dynamics of
\emph{all} solutions with $\Omega_m$ close to one initially. This whole
invariant set of solutions of co-dimension one therefore take the role of the
single tracker solution, which is associated with the limit $\lambda
\rightarrow \infty$ of asymptotically very steep potentials, such as the
inverse power-law potential in the limit $\phi \rightarrow 0$. Furthermore,
the present models give quite similar results to $\Lambda$CDM cosmology, even
if $\lambda \rightarrow 0$ much more slowly when $\phi \rightarrow \infty$
than a reasonably observationally viable inverse power-law potential, i.e.,
global and asymptotic bounds make the models more insensitive to how the
asymptotic limit $\lambda \rightarrow 0$ is approached when $\phi \rightarrow
\infty$. In addition, since the surface $x=0$ precisely describes
$\Lambda$CDM dynamics when $V=\mathrm{const}.$, it is possible to
continuously deform the separatrix surface $x=0$ by gradually changing the
potential. We stress that the examples we have used to illustrate the above
in this paper are just that, examples, and it is easy to produce other
explicit examples. In forthcoming work we are going to use dynamical systems
methods to obtain results for general potentials, and we are also going to
explore physical consequences more thoroughly.

1505.06903.extracted_bib
{plaI15 R.~Adam {et al. (Planck Collaboration), \newblock Planck Collaboration I, Planck 2015 results. I. Overview of products and results, \newblock arXiv:1502.01582.
1505.06903.extracted_bib
{plaXIII15 P.~A.~R.~Ade {et al. (Planck Collaboration) \newblock Planck 2015 results. XIII. Cosmological parameters, \newblock arXiv:1502.01589.
1505.06903.extracted_bib
{peerat88 P.~J.~E.~Peebles and B.~Ratra. \newblock Cosmology with a time-variable cosmological ``constant''. \newblock Astrophys.~J. {325, L17 (1988).%L17-L20 \newblock DOI: 10.1086/185100
1505.06903.extracted_bib
{ratpee88 B.~Ratra and P.~J.~E.~Peebles. \newblock Cosmological consequences of a rolling homogeneous scalar field. \newblock Phys.~Rev. D{37, 3406 (1988). \newblock DOI: 10.1103/PhysRevD.37.3406
1505.06903.extracted_bib
{steetal99 P.~J.~Steinhardt, L.~Wang, and I.~Zlatev. \newblock Cosmological tracking solutions. \newblock Phys.~Rev. D{59, 123504 (1999). \newblock DOI: 10.1103/PhysRevD.59.123504
1505.06903.extracted_bib
{tsu13 S.~Tsujikawa. \newblock Quintessence: A review. \newblock Classical Quantum Gravity {30, 214003 (2013). \newblock DOI: 10.1088/0264-9381/30/21/214003
1505.06903.extracted_bib
{waiell97 J.~Wainwright and G.~F.~R.~Ellis. \newblock {\em Dynamical Systems in Cosmology. \newblock {Cambridge {University {Press, Cambridge, (1997).
1505.06903.extracted_bib
{col03 A.~A.~Coley. \newblock {\em Dynamical Systems and Cosmology. \newblock Kluwer Academic Publishers, Dordrecht, (2003).
1505.06903.extracted_bib
{wei08 S.~Weinberg. \newblock {\em Cosmology. \newblock Oxford University Press, Oxford, (2008).
1505.06903.extracted_bib
{giamir10 R.~Giambo and J.~Miritzis. \newblock Energy exchange for homogeneous and isotropic universes with a scalar field coupled to matter. \newblock Classical Quantum Gravity {27, 095003 (2010). \newblock DOI: 10.1088/0264-9381/27/9/095003.
1505.06903.extracted_bib
{ugg13 C.~Uggla. \newblock Global cosmological dynamics for the scalar field representation of the modified Chaplygin gas. \newblock Phys. Rev. D{88, 064040 (2013). \newblock DOI: 10.1103/PhysRevD.88.064040.
1505.06903.extracted_bib
{leofad14 G.~Leon and C.~R.~Fadragas. \newblock Cosmological Dynamical Systems: And their Applications. \newblock LAP LAMBERT Academic Publishing (2012). \newblock arXiv:1412.5701
1505.06903.extracted_bib
{tam14 N.~Tamanini. \newblock Dynamical systems in dark energy models. \newblock PhD Thesis, University College, London, (2014).
1505.06903.extracted_bib
{alhugg15 A.~Alho and C. Uggla. \newblock Global dynamics and inflationary center manifold and slow-roll approximants. \newblock J. Math. Phys. {56, 012502 (2015).% - \newblock DOI: 10.1063/1.4906081.
1505.06903.extracted_bib
{alhetal15 A.~Alho, J.~Hell and C.~Uggla. \newblock Global dynamics and asymptotics for monomial scalar field potentials and perfect fluids. \newblock Classical Quantum Gravity {32, (14:145005) (2015) \newblock DOI: 10.1088/0264-9381/32/14/145005
1505.06903.extracted_bib
{blaetal04 R.~D.~Blandford, M.~Amin, E.~A.~Baltz, K.~Mandel, and P.~J.~Marshall. \newblock Cosmokinetics. \newblock Observing Dark Energy, ASP Conference Series, {339 17 (2005) \newblock DOI: 2005ASPC..339...27B
1505.06903.extracted_bib
{rapetal07 D.~Rapetti, S.~W.~Allen, M.~A.~Amin, and R.~D.~Blandford. \newblock A kinematical approach to dark energy studies. \newblock Mon.~Not.R~.~Astron.~Soc. {375, 1510 (2007)%1510-1520 \newblock Doi: 10.1111/j.1365-2966.2006.11419.x.
1505.06903.extracted_bib
{copetal98 E.~J.~Copeland, A.~R.~Liddle, and D.~Wands. \newblock Exponential potentials and cosmological scaling solutions. \newblock Phys. Rev. D{57, 4686 (1998). \newblock DOI: 10.1103/PhysRevD.57.4686
1505.06903.extracted_bib
{Nunes2000 A.~Nunes and J.~P.~Mimoso. \newblock On the potentials yielding cosmological scaling solutions. \newblock Phys. Lett. B{488, 423 (2000) \newblock DOI: 10.1016/S0370-2693(00)00919-9
1505.06903.extracted_bib
{uggetal95 C.~Uggla, R.~T.~Jantzen, and K.~Rosquist. \newblock Exact hypersurface-homogeneous solutions in cosmology and astrophysics. \newblock Phys. Rev. D{51 5522 (1995). %5522-5557 \newblock DOI: 10.1103/PhysRevD.51.5522
1505.06903.extracted_bib
{heiugg10 J.~M.~Heinzle and C.~Uggla. \newblock Monotonic functions in Bianchi models: why they exist and how to find them. \newblock Classical Quantum Gravity {27, 015009 (2010). \newblock DOI: 10.1088/0264-9381/27/1/015009
1505.06903.extracted_bib
{ngetal01 S.~C.~C.~Ng, N.~J.~Nunes, and F.~Rosati. \newblock Applications of scalar attractor solutions to cosmology. \newblock Phys.~Rev. D{64, 083510 (2001). \newblock DOI:10.1103/PhysRevD.64.083510.
1505.06903.extracted_bib
{ure12 L.~A.~Urena-Lopez. \newblock Unified description of the dynamics of quintessential scalar fields. \newblock J. Cosmol. Astropart. Phys. 03 (2012) 035. \newblock DOI: 10.1088/1475-7516/2012/03/035.
1505.06903.extracted_bib
{gon14 Y.~Gong. \newblock The general property of dynamical quintessence field. \newblock Phys.~Lett. B{342 (2014) 342.% 342-349 \newblock DOI:10.1016/j.physletb.2014.03.013.
1505.06903.extracted_bib
{colwai92 A.~A.~Coley and J.~Wainwright \newblock Qualitative analysis of two-fluid Bianchi cosmologies. \newblock Classical Quantum Gravity {9 651 (1992). \newblock DOI: 10.1088/0264-9381/9/3/008
1505.06903.extracted_bib
{heietal05 J.~M.~Heinzle, N.~R{\"ohr, and C.~Uggla. \newblock Matter and dynamics in closed cosmologies. \newblock Phys.~Rev. D{71 083506 (2005). \newblock DOI: 10.1103/PhysRevD.71.083506.
1505.06903.extracted_bib
{peerat03 P.~J.~E.~Peebles and B.~Ratra. \newblock The Cosmological Constant and Dark Energy. \newblock Rev.~Mod.~Phys. {75 559 (2003).
1505.06903.extracted_bib
{zlaetal99 I.~Zlatev, L.~Wang, and P.~J.~Steinhardt. \newblock Quintissence, Cosmic Coincidence, and the Cosmological Constant. \newblock Phys.~Rev.~Lett. {82, 896 (1999).
1505.06903.extracted_bib
{colwai92 A.~A.~Coley and J.~Wainwright \newblock Qualitative analysis of two-fluid Bianchi cosmologies. \newblock Class.~Quantum Grav. {9 651 (1992).
1505.06903.extracted_bib
{heietal05 J.~M.~Heinzle, N.~R{\"ohr, and C.~Uggla. \newblock Matter and dynamics in closed cosmologies. \newblock Phys.~Rev. {D71 083506 (2005). \newblock DOI: 10.1103/PhysRevD.71.083506.
1505.06903.extracted_bib
{car81 J.~Carr. \newblock {\em Applications of center manifold theory. \newblock Springer Verlag, New York, 1981.
1505.06903.extracted_bib
{ratqui92 B.~Ratra and A.~Quillen. \newblock Gravitational Lensing Effects in a Time-Variable Cosmological `Constant' Cosmology. \newblock Mon.~Not.~Roy.~Astr.~Soc. {259, 738 (1992).
1505.06903.extracted_bib
{bramar00 P.~Brax and J.~Martin. \newblock The Robustness of Quintissence. \newblock Phys.~Rev. {D 61, 103502 (2000).
1505.06903.extracted_bib
{braetal00 P.~Brax, J.~Martin, and A.~Riazuelo. \newblock Exhaustive study of cosmic microwave background anisotropies in quintessential scenarios. \newblock Phys.~Rev. {bf D 62, 103505 (2000). \newblock DOI:10.1103/PhysRevD.62.103505.
1505.06903.extracted_bib
{uza99 J.-P.~Uzan. \newblock Cosmological scaling solutions of nonminimally coupled scalar fields. \newblock Phys.~Rev. {D 59, 123510 (1999). \newblock DOI:10.1103/PhysRevD.59.123510.
1505.06903.extracted_bib
{lidshe98 A.~B.~Liddle and R.~J.~Scherrer. \newblock Classification of scalar field potentials with cosmological scaling solutions. \newblock Phys.~Rev. {D 59, 023509 (1998). \newblock DOI:10.1103/PhysRevD.59.023509.
1505.06903.extracted_bib
{chi10 T.~Chiba. \newblock Equation of state of tracker fields. \newblock Phys.~Rev. {D 81, 023515 (2010). \newblock DOI:10.1103/PhysRevD.81.023515. ******************** The stability of the fixed points is as follows: On the $Z=0$ boundary the fixed points M$^{0_{\pm$ are sources, PL is a sink, and FL$^{0$ a saddle, see eg. \cite{. On the $Z=1$ boundary the fixed points M$^{1_{\pm$ are sources, dS is a sink, and FL$^1$ a saddle. Since on $\Sigma_{\dagger>0$ ($\Sigma_{\dagger<0$), Z is monotonically increasing (decreasing), it follows that in the full state space M$^{0_{+$ and M$^{1_{-$ are sources, M$^{0_{-$ and M$^{1_{+$ are saddles, and Es is a saddle. We will assume that $\lambda_{0^2<3\gamma_m$, which results in that the fixed point $\mathrm{PL$ is outside the physical state space and that $\mathrm{Es$ is stable on the $Z=0$ subset, since the assumption $\gamma_m<2$ implies $\lambda^2_0<6$. Moreover the power-law solution $\mathrm{PL$, is accelerating if and only if $\lambda_{0^2<2$. Besides the line of fixed points FL$_{Z$ which is transversally hyperbolic and the dS fixed point which is non-hyperbolic having a zero eigenvalue corresponding to a 1-dimensional stable center manifold, all the remaining fixed points are hyperbolic. Since they lie at the intersection of invariant boundaries subsets, their stability readily follows from the analysis of the flow on these boundaries. We will assume that $\bar{\lambda^2<3\gamma_m$, which results in that the fixed point $\mathrm{Ei$ is outside the physical state space and that $\mathrm{Es$ is stable on the $Z=0$ subset, which also leads to that $q < - (3\gamma_m-2)/2$, i.e., this fixed point corresponds to an accelerating state if $\gamma_m>2/3$, which is the case for e.g. a pressureless fluid and a radiation equation of state. The solution space on this subset is depicted in Fig.~\ref{fig:Z0. \begin{figure[ht!] \begin{center \subfigure[$\gamma_{m=1$ and $\lambda_{-=1$]{\label{fig:Z0Dust \includegraphics[width=0.5\textwidth, trim = 0cm 3cm 0cm 0cm]{Z0_Boundary_Final.pdf \qquad