text:
\begin{document}

\maketitle

\begin{abstract}
A frame representation is used to derive a first order quasi-linear
symmetric hyperbolic system for a scalar field minimally coupled to
gravity. This procedure is inspired by similar evolution equations
introduced by Friedrich to study the Einstein-Euler system.  The
resulting evolution system is used to show that, for some classes of
scalar field potentials, small nonlinear perturbations of an
expanding Friedmann-Robertson-Walker background with a
selfinteracting scalar field source, decay exponentially to zero or
converge to constant values. As an application of our results we
discuss the nonlinear stability of flat power-law solutions produced by an
exponential potential. It is found that for a certain range of the
parameter $p$, there is stability even in the absence of
accelerated expansion.

\end{abstract}

\section{Introduction}

An important problem of classical mathematical Cosmology concerns the
asymptotic stability of spatially homogeneous and isotropic
spacetimes. Within this class of spacetimes, those having a nonlinear
scalar field as the matter model have been extensively used to model
early and late times cosmological scenarios ---see
e.g. \cite{2004CQG_Rendall,2005CQG_Rendall,2006LNP_Rendall} and
references therein.  Recently, Ringstr\"{o}m
\cite{2008Ringstrom,2009Ringstrom} has proved that small perturbations
of the initial data of scalar field cosmological solutions to the {\it Einstein
Field Equations} (EFE) with accelerated expansion have maximal
globally hyperbolic developments that are future causally geodesically
complete. In particular, in \cite{2009Ringstrom} the case of power-law
solutions due to an exponential potential has been considered as an
application of the more general methods developed in
\cite{2008Ringstrom}. The nonlinear asymptotic stability of power-law
inflation solutions has also been discussed by Heinzle \& Rendall
\cite{Rendall2007}  using Kaluza-Klein reductions and the methods of
Anderson \cite{2005AHP_MTAnderson} ---the latter, in turn, inspired by
Friedrich's analysis of the stability of the de Sitter spacetime
\cite{Friedrich1986}.

\medskip
A natural way to analyse the stability of spacetimes is to ask whether
small perturbations of a given solution to the EFE asymptotically
decay to a background solution. Most of approaches to this question
have been limited to the use of linear or higher-order truncated
perturbation theory, and thus, they never take fully into account the
nonlinearity of the EFE ---see
e.g. \cite{CQG1990_Stewart,CQG1992_Bruni,BMT} and also
\cite{Rendall2009}. This type of analysis has been hampered by the
lack of a suitable formulation of the EFE for which the theory of
systems of first order hyperbolic partial differential equations can
be applied.

\medskip
In \cite{1996CQGra..13.1451F} Friedrich has introduced a frame
representation of the vacuum EFE. The evolution equations implied by
this alternative representation of the equations of General
Relativity constitute a \emph{first-order quasi-linear symmetric
hyperbolic  system (FOSH)}. In general, these systems are of the form
\begin{equation}
 \mathbf{A}^{0}(t,\bmath{x},\mathbf{u})\partial_{t}\mathbf{u}-\mathbf{A}^{j}(t,\bmath{x},\mathbf{u})\partial_{j}\mathbf{u}=\mathbf{B}(t,\bmath{x},\mathbf{u})\mathbf{u}
\label{FOSHsystem}
\end{equation}
where $\mathbf{u}=\mathbf{u}(\bmath{x},t)$ is a smooth vector-valued
function of dimension $s$ with domain in $\Sigma\times[0,T]$ with
$\Sigma$ a spacelike 3-dimensional manifold. Moreover, $\mathbf{A}^{0}$,
$\mathbf{A}^{j}$, $j=1,\,2,\,3$, and $\mathbf{B}$ denote smooth
$s\times s$ matrix valued-functions. The matrices  $\mathbf{A}^0$ and $\mathbf{A}^{j}$ are
symmetric. In addition, the matrix $\mathbf{A}^{0}$ is positive definite.  The
operators $\partial_{t}$ and $\partial_{j}$ stand, respectively, for
the partial derivatives with respect to the coordinates $t$ and
$x^{j}$.

\medskip
The construction for vacuum spacetimes given in
\cite{1996CQGra..13.1451F} has been extended in
\cite{PhysRevD.57.2317,Friedrich_Rendall} to the case of a
\emph{perfect fluid} using a Lagrangian description of the fluid flow
---see also \cite{2002TMNA_Bruhat_York,2002LNP...592...29C}. In both
the vacuum and the perfect fluid cases the introduction of a frame
formalism gives rise to extra gauge freedom.
This freedom is associated to the evolution of the spatial frame coefficients along
the flow of the time-like frame. If one fixes conveniently this gauge
(using, for example, the Fermi gauge) one obtains a hyperbolic reduction for
the evolution equations. As a consequence, given smooth initial data
satisfying the constraints, local in time existence and uniqueness of a
solution to the EFE can be established ---see
e.g. \cite{Friedrich_Rendall,lrr-1998-3} and also
\cite{2002TMNA_Bruhat_York,2002LNP...592...29C} for details.

\medskip

A natural way of performing a stability analysis is to consider a sequence of smooth initial data sets
$\mathbf{u}^{\epsilon}$ for the EFE satisfying the constraints
equations on a Cauchy hypersurface $\Sigma$. The sequence is assumed
to depend continuously on the parameter $\epsilon$ in such a way
that the limit $\epsilon\rightarrow0$ renders the reference
solution $\mathring{\mathbf{u}}$. In particular, one can write the full
solution to the EFE as the Ansatz
\begin{equation}
\label{Ansatz}
 \mathbf{u}^{\epsilon}=\mathring{\mathbf{u}}+\epsilon\breve{\mathbf{u}}^{\epsilon},
\end{equation}
where $\breve{\mathbf{u}}^{\epsilon}$ is a (nonlinear)
perturbation whose size is controlled by the parameter
$\epsilon$. Using the Ansatz \eqref{Ansatz} in equation
\eqref{FOSHsystem}, and writing
\begin{equation*}
\begin{aligned}
 \mathbf{B}(\mathring{\mathbf{u}}+\epsilon\breve{\mathbf{u}}^{\epsilon})&=\mathring{\mathbf{B}}(\mathring{\mathbf{u}})+\epsilon\breve{\mathbf{B}}(\breve{\mathbf{u}},\epsilon),\\
 \mathbf{A}^{\mu}(\mathring{\mathbf{u}}+\epsilon\breve{\mathbf{u}}^{\epsilon})&=\mathring{\mathbf{A}}^{\mu}(\mathring{\mathbf{u}})+\epsilon\breve{\mathbf{A}}^{\mu}(\breve{\mathbf{u}},\epsilon),\quad \mu=0,1,2,3
\end{aligned}
\end{equation*}
we are led to consider the following initial value problem for the
nonlinear perturbations:
\begin{equation}
\begin{aligned}
\left(\mathring{\mathbf{A}}^{0}(t,\bmath{x},\mathring{\mathbf{u}})+\epsilon \breve{\mathbf{A}}^{0}(t,\bmath{x},\breve{\mathbf{u}},\epsilon)\right)\partial_{t}\breve{\mathbf{u}}-\left(\mathring{\mathbf{A}}^{j}(t,\bmath{x},\mathring{\mathbf{u}})+\epsilon \breve{\mathbf{A}}^{j}(t,\bmath{x},\breve{\mathbf{u}},\epsilon)\right)\partial_{j}\breve{{\bf u}}&=\left(\mathring{{\bf B}}(t,\bmath{x},\mathring{{\bf u}})+\epsilon \breve{{\bf B}}(t,\bmath{x},\breve{{\bf u}},\epsilon)\right)\breve{{\bf u}}, \\
\breve{{\bf u}}(\bmath{x},0)&=\breve{{\bf u}}_{0}(\bmath{x}).
\end{aligned} \label{nonlinearequation}
\end{equation}
A particular approach to the existence and stability
of solutions to the Cauchy problem
\eqref{nonlinearequation}  has been
discussed in \cite{Acta98,1997funct.an..3003K,JMP2}. In this approach,
the stability of solutions follows from the existence of eigenvalues
for the nonprincipal part of the linearised system ($\epsilon=0$)
having a negative real part. Stability results for the case where the
coefficients of the linearized system are constant matrices are
well known \cite{Acta98}. In the case where the system is only
strongly hyperbolic, the inner product in $L^{2}$ has to be replaced
by the so-called $\mathcal{H}$-inner product ---see
\cite{1997funct.an..3003K}. A procedure to analyse stability in the
case of systems where $\mathring{{\bf B}}$ has vanishing eigenvalues
has been given in \cite{1997funct.an..3003K} ---see also
\cite{JMP1,JMP2}. These methods can be easily generalized to systems
of the type considered here where the matrices $\mathring{{\bf B}}$,
$\mathring{\mathbf{A}}^{\mu}$ are not constant but depend smoothly on
time ---see also \cite{PhysRevD.60.083507}.

\medskip
The approach described in the previous paragraph has been applied by
Reula \cite{PhysRevD.60.083507} to the Einstein-perfect fluid system of
\cite{PhysRevD.57.2317} with a cosmological constant, to prove the exponential decay of nonlinear
perturbations for a wide class of homentropic fluids in flat {\it Friedmann-Robertson-Walker} (FRW)
backgrounds. The analysis of \cite{PhysRevD.60.083507} gives, as only
restriction on the equation of state, that the sound speed is less than
one-third of the speed of light. An advantage of this approach is that
it avoids the problem of gauge-dependence in perturbation theory.

\medskip
In the present article we pursue a similar approach to the one used by
Reula in \cite{PhysRevD.60.083507} to analyse the nonlinear stability
of FRW spacetimes with a nonlinear scalar field. To this end, we
first construct a first order symmetric hyperbolic system for the EFE
with a scalar field as the matter source. This construction is
performed by splitting the wave equation for the scalar field into two
first order equations. In our analysis, the scalar field is used
to construct an adapted orthogonal frame.
Written in terms of this adapted frame, the energy-momentum tensor is
diagonal, independently of further gauge choices.
A similar construction has been considered
in the analysis of linear perturbations in
\cite{CQG1992_Bruni,CQG1988_Madsen}.

\medskip
An important difference between the system treated in the present
paper and the one in reference \cite{PhysRevD.60.083507} is that in our case the
matrix $\mathring{\mathbf{B}}$ has eigenvalues which are zero. These
vanishing eigenvalues appear in the block containing the perturbed
matter variables and expansion.
As it will be discussed in the sequel, this feature is related to the fact
that a subspace of the perturbations converge asymptotically to
constant values. Our analysis of the eigenvalues of the nonprincipal
part of the linearized system uses results from the so-called
\emph{Routh-Hurwitz problem}.  In particular, for an ever expanding
FRW-nonlinear scalar field background with spatial topology
$\mathbb{T}^{3}$,
i.e. $([0,\infty)\times\mathbb{T}^{3},\mathring{\bmath{g}},\mathring{\phi})$,
we show that, for scalar field potentials satisfying certain
inequalities, nonlinear perturbations exist and asymptotically
exponentially decay to zero or converge to constant values. In the
suitable limit, we recover the results of Ringstr\"om in
\cite{2008Ringstrom}. As an application of our results, we show that
for power-law solutions produced by an exponential potential with
\[
p>\tfrac{\sqrt{61}-1}{15} \approx 0.454
\]
and
flat spatial sections, sufficiently small nonlinear perturbations
decay exponentially to zero or converge to constant values. We have
confirmed this behaviour in some simple numerical computations of linear
perturbations. Our analysis shows that nonlinear stability of power-law
solutions with exponential potentials is possible even in the absence
of accelerated expansion ---a possibility that, hitherto, has not been
discussed before.

\medskip
Our main result is proved in the Sobolev norm
$H^{k}\left(\mathbb{T}^{3}\right)$ for $k\geq5$. As a consequence of
the use of results from the Routh-Hurwitz problem, the conditions we
obtain on the potential are, in a sense, sharp and cannot be improved
by the present methods. Given a particular background Cosmological
solution to the Einstein-nonlinear scalar field system ---calculated,
say, numerically--- one can verify explicitly whether the conditions
are satisfied or not.

\subsubsection*{Overview of the article}
The paper is organized as follows: in Section 2, we recall Friedrich's
frame formulation of the EFE. In Section 3, we discuss some relevant
properties of scalar field satisfying a nonlinear wave
equation. In Section 4 we discuss the conditions under which the Einstein-Friedrich-nonlinear
scalar field system is well-posed, in the sense that it forms a
symmetric hyperbolic system ---see  Theorem \ref{Theorem:wellposedness}.
Finally, in Section 5 we give the conditions for which there is convergence to
constant values or exponential decay of small nonlinear perturbations
on a FRW-nonlinear scalar field background. This is the main result
of the paper and we summarise it in Theorem \ref{Theorem:Main}.  We use units
such that $8\pi G=c=1$.

\section{Friedrich's frame formulation of Einstein Field Equations}

In this section we provide a brief introduction to Friedrich's frame
formulation of the Einstein field equations. The basic equation of
Friedrich's construction is the contracted Bianchi identity. From the
latter it is possible to deduce hyperbolic propagation equations for
the conformal Weyl tensor for a wide spectrum of gauge choices.

\subsection{Basic definitions}
Let $(\mathcal{M},\bmath{g})$ denote a 4-dimensional spacetime with
$\bmath{g}$ a Lorentzian metric of signature $(-,+,+,+)$. In order to
implement the frame formulation of the Einstein field equations, one
defines locally an orthonormal moving frame or $\it{tetrad}$ with
respect to the metric $\bmath{g}$ in an open neighbourhood
$\mathcal{U} \subset \mathcal{M}$.  The frame is a set
$\left\{\bmath{e}_{a}\right\}$ of linearly independent vector fields
in the tangent space $T_{p}\left(\mathcal{M}\right)$ at each
point $p\in\mathcal{U}$ such that
\begin{equation}
 \bmath{g}\left(\bmath{e}_{a},\bmath{e}_{b}\right)=\eta_{ab}, \quad a,
 \, b=0,\ldots 3,
\label{orthon}
\end{equation}
where $\eta_{ab}=\text{diag}\left(-1,1,1,1\right)$. The \emph{norm} of
a vector field, $\bmath{v}\in T_{p}\left(\mathcal{M}\right)$, in an orthonormal frame is defined as
\begin{equation*}
 |\bmath{v}|^{2}\equiv\bmath{g}\left(\bmath{v},\bmath{v}\right)=v^{a}v^{b}\eta_{ab}
\label{norm}
\end{equation*}
and in terms of a coordinate basis set
$\left\{\partial_{\alpha}\right\}$
we have $\bmath{e}_{a}=e_{a}^{\,\,\mu}\partial_{\mu}$. Condition \eqref{orthon} gives
\begin{equation*}
 \eta_{ab}=e_{a\,}^{\,\,\mu}e_{b}^{\,\,\nu}g_{\mu\nu}.
\end{equation*}
where $\mu,\nu=0,1,2,3$ denote spacetime tensor indices. The \emph{frame commutator} is written as
\begin{equation}
 [\bmath{e}_{a},\bmath{e}_{b}]=c^{c}_{\,\,ab}\bmath{e}_{c},
\label{commutatorframes}
\end{equation}
where $c^{c}_{\,\,ab}$ are the \emph{structure coefficients}.  The
\emph{dual basis} or \emph{coframe} is the set of linear forms
$\{\bmath{\theta}^{b}\}$ belonging to the dual space
$T^{*}_{p}\left(\mathcal{M}\right)$ at each point
$p\in\mathcal{U}$ defined by the pairing
$\langle\bmath{\theta}^{b},\bmath{e}_{a}\rangle=\delta^{b}_{a}$. In
terms of the dual basis we can write condition \eqref{orthon} as
\begin{equation*}
 \bmath{g}=-\left(\bmath{\theta}^{0}\right)^{2}+\sum^{3}_{a=1}\left(\bmath{\theta}^{a}\right)^{2}.
\end{equation*}
The \emph{spacetime (Levi-Civita) connection} in an orthonormal basis
is defined by
\begin{equation*}
\bmath{\nabla}_{a}\bmath{e}_{b}\equiv\gamma^{c}_{\,\,\,ba}\bmath{e}_{c},
\end{equation*}
where $\gamma^{c}_{\,\,ba}$ are the \emph{connections
coefficients}. The covariant derivative of a tensor in $\mathcal{M}$
can be written as
\begin{eqnarray*}
&& \bmath{\nabla}_{a}v_{q_{1}...q_{s}}{}^{p_{1}...p_{r}}=\bmath{e}_{a}\left(v_{q_{1}\cdots q_{s}}{}^{p_{1} \cdots p_{r}}\right)+\gamma_{\,\,\,fa}^{p_{1}}v_{q_{1} \cdots q_{s}}{}^{f \cdots p_{r}}+ \cdots \\
&& \hspace{4cm} \cdots+\gamma_{\,\,\,fa}^{p_{r}}v_{q_{1} \cdots q_{s}}{}^{p_{1} \cdots f}-\gamma_{\,\,\,q_{1}a}^{f}v_{f \cdots q_{s}}{}^{p_{1} \cdots p_{r}}-\cdots-\gamma_{\,\,\,q_{s}a}^{f}v_{q_{1} \cdots f}{}^{p_{1} \cdots p_{r}}.
\end{eqnarray*}
The torsion free and metric compatibility conditions imply, respectively, that
\begin{equation*}
 c^{c}_{\,\,ab}=\gamma^{c}_{\,\,ba}-\gamma^{c}_{\,\,ab}, \quad \gamma^{e}_{\,\,ba}\eta_{ec}-\gamma^{e}_{\,\,ca}\eta_{eb}=0.
\end{equation*}
The equations for the frame coefficients
$\left\{e_{a}^{\,\,\mu}\right\}$ are given by equation
\eqref{commutatorframes} in terms of the connection coefficients. In
turn, equations for the connection coefficients are obtained from the
Ricci identity
\begin{equation}
 R^{a}_{\,\,bcd}=\bmath{e}_{c}\left(\gamma^{a}_{\,\,bd}\right)-\bmath{e}_{d}\left(\gamma^{a}_{\,\,bc}\right)+\gamma^{a}_{\,\,fc}\gamma^{f}_{\,\,bd}-\gamma^{a}_{\,\,fd}\gamma^{f}_ {\,\,bc}-\gamma^{a}_{\,\,bf}\left(\gamma^{f}_{\,\,dc}-\gamma^{f}_{\,\,cd}\right).
\label{RicciIdentity}
\end{equation}
The Riemann tensor can be decomposed in terms of the conformal Weyl
tensor $\bmath{C}$ and the Schouten tensor $\bmath{S}$ as
\begin{equation}
 R^{a}{}_{bcd}=C^{a}{}_{bcd}+\delta^{a}{}_{[c}S_{d]b}-\eta_{b[c}S_{d]}{}^a.
\label{decompositionRiemann}
\end{equation}
For future use we introduce the \emph{Friedrich tensor} $\bmath{F}$ via
\begin{equation}
 F_{abcd}\equiv {C}_{abcd}-\eta_{a[c}S_{d]b},
\label{F}
\end{equation}
and its dual with respect to the last pair of indices
\begin{equation}
 {^{\star}}F_{abcd}={^{\star}}{C}_{abcd}+\frac{1}{2}S_{pb}\epsilon^{p}{}_{acd},
\label{Fdual}
\end{equation}
where $\epsilon_{abcd}$ is the usual Levi-Civita totally antisymmetric
symbol with $\epsilon_{0123}=1$. In terms of the Friedrich tensor, one
finds that the contracted Bianchi identities read
\begin{equation}
 \nabla_{a}F^{a}_{\,\,bcd}=0,\quad\nabla_{a}{^{\star}}F^{a}_{\,\,bcd}=0.
\label{BianchiIdentities}
\end{equation}

\subsection{Orthonormal decomposition of the field equations}

The field equations of Friedrich's frame formulation of the Einstein
equations are given by equations \eqref{commutatorframes},
\eqref{RicciIdentity} and \eqref{BianchiIdentities}, together with the
decomposition \eqref{decompositionRiemann}. The independent variables
of the system are
\begin{equation*}
 \left(e_{a}{}^\mu,\gamma^{a}{}_{bc},C^{a}{}_{bcd},S_{bc}\right).
\end{equation*}

\medskip
In what follows, we will shall decompose the field equations and
relevant tensors in terms of their parallel and orthogonal components
with respect to the time-like frame. We write
$\bmath{N}\equiv \bmath{e}_{0}$ and set
\begin{equation*}
\bmath{N}=N^{a}\bmath{e}_{a},\quad N^{a}=\delta_0{}^a,
\end{equation*}
where $N_{a}=-\delta_a{}^0$ in our signature. In terms of these
objects, tensor fields which are orthogonal to the timelike
frame-vector are defined by
\begin{equation*}
 T_{a_{1}\dots a_{p}\dots a_{q}}N^{a_{p}}=0,\quad p=1, 2,\ldots, q.
\end{equation*}

\medskip
Next, one defines
\begin{equation*}
 h_{ab} \equiv \eta_{ab}+N_{a}N_{b},
\end{equation*}
where the projector onto the orthogonal 3-subspaces satisfies
$h_{a}{}^{c}=\eta^{bc}h_{ab}$. The \emph{spatial covariant derivative} is
then given by
\begin{equation*}
 D_{a}T_{q_{1}\cdots q_{r}}=h_{a}{}^{b}h_{q_{1}}{}^{p_{1}}\cdots h_{q_{r}}{}^{p_{r}}\nabla_{b}T_{p_{1}\cdots p_{r}}.
\end{equation*}
In particular, one has that
\begin{equation*}
 D_{a}h_{bd}=0,\quad D_{a}\epsilon_{bcd}=0,
\end{equation*}
where $\epsilon_{bcd}$ is the spatial Levi-Civita symbol and the indices
run from $1$ to $3$.

\medskip
In order to further proceed with the geometric decomposition one
defines the \emph{acceleration vector} by
\begin{equation*}
 \bmath{a}\equiv \bmath{\nabla}_{0}\bmath{e}_{0}=\gamma^{p}{}_{00}\bmath{e}_{p}, \qquad p=1,2,3.
\end{equation*}
It follows then that $a^{p}=\gamma^{p}_{\,\,\,00}$ or equivalently,
$a_{p}=\gamma^{0}_{\,\,\,p0}$. We will also consider the so-called
\emph{Weingarten map} given by
\begin{equation*}
 \bmath{\chi}(\bmath{e}_{a})\equiv \bmath{\nabla}_{a}\bmath{e}_{0}=\gamma^{p}_{\,\,\,0a}\bmath{e}_{p},\quad a,p=1,2,3.
\end{equation*}
It follows then that $\chi_{a}{}^{p}=\gamma^{p}{}_{0a}$.  The tensor
$\chi_{ab}$ can be written in terms of its irreducible parts as
\begin{equation*}
 \chi_{ab}=\gamma^{0}{}_{ba}=(\chi^{ST})_{ab}+\frac{1}{3}\chi h_{ab}+(\chi^{A})_{ab},
\end{equation*}
where $(\chi^{ST})_{ab}$, $\chi$, $(\chi^{A})_{ab}$ denote,
respectively, its symmetric trace-free, trace and antisymmetric
parts. If the flow of $\bmath{e}_{0}$ is hypersurface orthogonal then
one has that $(\chi^{A})_{ab}=0$ and that
\begin{equation}
 \frac{1}{2}\pounds_{N}h_{ab}=\chi_{(ab)}=(\chi^{ST})_{ab}+\frac{1}{3}\chi h_{ab},
\label{LieDer3metric}
\end{equation}
where $\pounds_{N}$ denotes the Lie derivative along $\bmath{N}$ and
$\nabla_{a}N^{p}=-N_{a}a^{p}+\chi_{a}{}^{p}$.

\medskip
Finally, the 4-dimensional Levi-Civita symbol is also decomposed using
\begin{equation*}
 \epsilon_{abcd}=2\epsilon_{ab[c}N_{d]}-2N_{[a}\epsilon_{b]cd}.
\end{equation*}
Now, defining $\tilde{F}_{bcd}\equiv \nabla_{a}F^{a}{}_{bcd}$ it follows that
the first contracted Bianchi identity can be written as
\begin{equation}
 \tilde{F}_{bcd}=N_{b}\left[\tilde{F}^{\prime}_{0c0}N_{d}-\tilde{F}^{\prime}_{0d0}N_{c}\right]+2\tilde{F}^{\prime}_{b0[c}N_{d]}-N_{b}\tilde{F}^{\prime}_{0cd}+\tilde{F}^{\prime}_{bcd}=0,
\label{GeomDecBI}
\end{equation}
where contractions with $\bmath{N}$ are denoted by the index $0$, and
the prime ${}^\prime$ indicates that the remaining indices are
spatial. For example,
$\tilde{F}^{\prime}{}_{b0d}\equiv h_{b}{}^{q}N^{r}h_{d}{}^{s}\tilde{F}_{qrs}$.

\medskip
Given the vector $\bmath{N}$, the Weyl tensor is uniquely determined
through its \emph{electric} and \emph{magnetic} parts. These are
defined, respectively, by
\begin{equation*}
 E_{ab}\equiv h_{a}{}^{p}h_{b}{}^{q}N^{c}N^{d}C_{pqcd},\quad B_{bd}\equiv h_{b}{}^{p}h_{d}{}^{q}N^{a}N^{c}C^{\star}_{apcq}.
\end{equation*}
In terms of the latter, the Weyl tensor and its dual can be written as
\begin{eqnarray}
&& C_{abcd}=2\left[l_{a[c}E_{d]b}-l_{b[c}E_{d]a}\right]-2\left[N_{[c}B_{d]p}\epsilon^{p}{}_{ab}+N_{[a}B_{b]p}\epsilon^{p}{}_{cd}\right]
\label{DecompositionWeyl}\\
&& {^{\star}}C_{abcd}=2N_{[a}E_{b]p}\epsilon^{p}{}_{cd}-4E_{p[a}\epsilon_{b]}{}^{p}{}_{[c}N_{d]}-4N_{[a}B_{b][c}N_{d]}-B_{pq}\epsilon^{p}{}_{ab}\epsilon^{q}{}_{cd}
\label{DecompositionWeylDual}
\end{eqnarray}
where $l_{ab}\equiv h_{ab}+N_{a}N_{b}$.

\section{Nonlinear Scalar fields in the frame formalism}

In this section, we introduce a description on nonlinear scalar
fields which is particularly well adapted to the analysis to be
carried out in the present article.

\subsection{Basic equations}
 In what follows we will consider scalar fields $\phi\in C^{\infty}\left(\mathcal{M}\right)$
with a smooth positive selfinteracting potential $\mathcal{V}(\phi)$, and energy-momentum tensors of the form
\begin{equation*}
\bmath{T}=\bmath{\psi}\otimes\bmath{\psi}-\left(\frac{1}{2}
  |\bmath{\psi}|^{2} +\mathcal{V}(\phi)\right)\bmath{g},
\end{equation*}
where we have defined the 1-form
\begin{equation*}
\bmath{\psi}\equiv \bmath{\nabla}\phi
\end{equation*}
with $\bmath{\nabla}$ the spacetime connection. Accordingly, we define
\begin{equation}
 \psi_{a}\equiv \bmath{\psi}\left(\bmath{e}_{a}\right)=(\psi,\psi^{\prime}_{a}),
\label{definitionpsi}
\end{equation}
where we have written
\begin{equation}
 \psi\equiv \psi_{0}=\pounds_{N}\phi
\label{psi_{0}_frame}
\end{equation}
and
\begin{equation}
\psi^{\prime}_{a}\equiv h_{a}{}^b \psi_{b}=D_{a}\phi.
\end{equation}
The components of the energy-momentum tensor $\bmath{T}$ with respect
to the tetrad $\{\bmath{e}_a\}$ are given by
\begin{equation}
\label{stress}
 T_{ab}=\psi_{a}\psi_{b}-\left(\frac{1}{2} |\bmath{\psi} |^{2}+\mathcal{V}(\phi)\right)\eta_{ab},
\end{equation}
while its trace is
\[
T=-|\bmath{\psi} |^{2}-4\mathcal{V}(\phi).
\]
The Einstein field equations in the form
\[
\bmath{R}=\bmath{T}-\frac{1}{2}\bmath{g}\left(\text{Tr}_{\bmath{g}}\bmath{T}\right),
\]
then imply for the components of the Ricci tensor that
\begin{equation*}
 R_{ab}=\psi_{a}\psi_{b}+\mathcal{V}(\phi)\eta_{ab},
\end{equation*}
while the Ricci scalar is given by
\[
R=-T=|\bmath{\psi}|^{2}+4\mathcal{V}(\phi).
\]
From these expressions it follows that the components of the Schouten
tensor with respect to the frame $\{\bmath{e}_a\}$ are given by
\begin{equation*}
 S_{ab}=\psi_{a}\psi_{b}-\frac{1}{3}\left(\frac{1}{2}|\bmath{\psi}|^{2}-\mathcal{V}(\phi)\right)\eta_{ab}.
\end{equation*}

\subsection{Gauge considerations}
In order to construct and adapted frame to our particular problem,
we let $\bmath{\psi}\equiv\alpha\bmath{e}_{0}$. It follows that
\begin{equation}
 \quad\psi^{a}=\alpha\delta^{a}_{0},
\label{fixingpsi}
\end{equation}
so that
\[
\alpha=-\psi\quad\text{and}\quad D^{a}\phi=0.
\]
Accordingly,
\begin{equation}
 |\bmath{\psi}|^{2}=\bmath{g}\left(\bmath{\psi},\bmath{\psi}\right)=\alpha^{2}\eta_{00}=-\alpha^{2},
 \quad\alpha=\pm\sqrt{-|\bmath{\psi}|^{2}}.
\label{definitionalpha}
\end{equation}
If the vector $\bmath{\psi}$ is taken to be future oriented, then one must choose $\alpha$ to be positive. In terms of a coordinate basis the latter implies
\begin{equation}
 \psi^{\mu}=\alpha e_{0}{}^{\mu}=-\psi e_{0}{}^{\mu}, \quad  e_{0}{}^{\mu}=\frac{\nabla^{\mu}\phi}{\sqrt{-|\bmath{\psi}|^{2}}}.
\label{fixingpsicoordinates}
\end{equation}
and
\begin{equation}
D_{a}\phi=0, \quad  e{}^\prime_{a}{}^\mu\nabla_{\mu}\phi=0 .
\label{constraint_Da_phi}
\end{equation}
With this choice, we have that
\[
\psi_{a}=-\psi N_{a},
\]
 and therefore
\begin{eqnarray}
 &&T_{ab}=\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)N_{a}N_{b}+\left(\frac{1}{2}\psi^{2}-\mathcal{V}(\phi)\right)h_{ab},\\
 &&S_{ab}=\frac{1}{3}\left(\frac{5}{2}\psi^{2}-\mathcal{V}(\phi)\right)N_{a}N_{b}+\frac{1}{3}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)h_{ab}.
\label{Shoutenscalarfield}
\end{eqnarray}

\medskip
Using \eqref{definitionpsi} and \eqref{definitionalpha}, the equation
for the conservation of the energy-momentum tensor takes the form:
\begin{eqnarray}
 && \nabla^{a}T_{ab}=\nabla^{a}\left(\psi^{2}N_{a}N_{b}+\left(\frac{1}{2}\psi^{2}-\mathcal{V}(\phi)\right)\eta_{ab}\right) \\
&& \phantom{\nabla^{a}T_{ab}}=2\psi N_{b}N^{a}\left(\nabla_{a}\psi\right)+\psi^{2}\left(N_{b}\left(\nabla_{a}N^{a}\right)+N^{a}\left(\nabla_{a}N_{b}\right)\right)+\nabla_{b}\left(\frac{1}{2}\psi^{2}-\mathcal{V}(\phi)\right) \\
&&\phantom{\nabla^{a}T_{ab}} =\left(2\psi\pounds_{\text{\tiny{\bmath{N}}}}\psi+\psi^{2}\chi+\psi\frac{d\mathcal{V}}{d\phi}\right)N_{b}+\psi^{2}a_{b}+\psi\nabla_{b}\psi=0.
\label{div_EMtensor}
\end{eqnarray}
From the latter, projecting with respect to the timelike frame one obtains:
\begin{eqnarray}
&&N^{b}\left(\nabla^{a}T_{ab}\right)=0, \quad \pounds_{\text{\tiny{\bmath{N}}}}\psi+\chi\psi+\frac{d\mathcal{V}}{d\phi}=0,
\label{KGequation} \\
&& h_{c}^{\,\,\,b}\left(\nabla^{a}T_{ab}\right)=0, \,\,\,\quad\quad\quad D_{c}\psi+\psi a_{c}=0.
\label{momentumcons}
\end{eqnarray}
Moreover, using that $D_{a}\phi=0$ in the orthogonal subspaces to $\bmath{e}_{0}$, one obtains from equation \eqref{commutatorframes} that
\begin{equation*}
\left[\bmath{e}^{\prime}_{a},\bmath{e}^{\prime}_{b}\right]\phi=2\,(\chi^{A})_{ab}\psi=0,
\end{equation*}
which implies that
\begin{equation}
(\chi^{A})_{ab}=0.
\label{commutatorframes_scalarfield}
\end{equation}

\medskip
\noindent
\textbf{Remark.} Following Friedrich in \cite{Friedrich_Rendall}, one
could as well have defined
\begin{equation}
\nabla^{a}T_{ab}=q_{b}+qN_{b}, \quad J_{ab}=\nabla_{[a}q_{b]}.
\label{Alternative}
\end{equation}
Then, instead of using the  condition on the vanishing of the divergence of the
energy-momentum tensor, one could include the equations $q=0$ and
$q_{b}=0$ as a part of the equations determining the Einstein-scalar
field system in the frame representation. Once the gauge is fixed, the
first equation in \eqref{Alternative} appears in the reduced system of
evolution equations while the second part is regarded as a {\it zero
quantity} ---see equation (4.44) in \cite{Friedrich_Rendall}. It can
be shown that the zero quantities satisfy a system of \emph{subsidiary
evolution equations}. For this, it can be shown that the zero
quantities vanish if they are zero on the initial hypersurface. For
the quantity $q_{b}$ the relevant subsidiary equation is given in
equation (4.70) of \cite{Friedrich_Rendall}. We also notice that the
evolution for the acceleration can be computed from the tensor
$J_{ab}$.

\section{The Einstein-Friedrich-nonlinear scalar field system}

In this section, we derive a first order symmetric hyperbolic system
for the EFE coupled to a nonlinear scalar field. Making use of the
Bianchi identity and the energy-momentum tensor of equation
\eqref{stress}, we derive the propagation equations for the
\emph{electric} and \emph{magnetic} parts of the conformal Weyl
tensor. After fixing the gauge, we complete the reduced system of
evolution equations by deriving equations for the frame and the
connection coefficients. In the last part of this section we make some
remarks concerning the hyperbolicity of the system.

\subsection{Basic expressions}
We start by computing the various components for the Friedrich tensor
$\bmath{F}$. Using equations \eqref{DecompositionWeyl} and
(\ref{Shoutenscalarfield}), one finds that the definition \eqref{F} implies
\begin{equation}
 \begin{aligned}
  &F^{\prime}_{00c0}=0=-F^{\prime}_{000c},\quad F^{\prime}_{00cd}=0=-F^{\prime}_{00dc}, \\
  &F^{\prime}_{a00d}=-E_{ad}+\frac{1}{6}\left(\frac{5}{2}\psi^{2}-\mathcal{V}(\phi)\right)h_{ad}=-F^{\prime}_{a0d0}, \\
  &F^{\prime}_{ab0d}=B_{dp}\epsilon^{p}{}_{ab}=-F^{\prime}_{abd0}=-F^{\prime}_{ba0d}, \\
  &F^{\prime}_{0bcd}=B_{bp}\epsilon^{p}{}_{cd}=-F^{\prime}_{0bdc}=-F^{\prime}_{b0cd}, \\
  &F^{\prime}_{0b0d}=E_{bd}+\frac{1}{6}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)h_{bd}=-F^{\prime}_{0bd0,} \\
  &F^{\prime}_{abcd}=-2\left(h_{b[c}E_{d]a}-h_{a[c}E_{d]b}\right)-\frac{1}{6}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)\left(h_{ac}h_{db}-h_{ad}h_{cb}\right),
 \end{aligned}
\label{Fcomponents}
\end{equation}
with the nonvanishing traces
\begin{equation}
 \begin{aligned}
  & h^{ac}F^{\prime}_{a0c0}=\frac{1}{2}\left(\frac{5}{2}\psi^{2}-\mathcal{V}(\phi)\right), \\
  &h^{bd}F^{\prime}_{0b0d}=\frac{1}{2}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)=-h^{bd}F^{\prime}_{0bd0}, \\
  &h^{bd}F^{\prime}_{abcd}=E_{ac}-\frac{1}{3}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)h_{ac}=F^{\prime b}{}_{abc},\\
  &h^{ac}h^{bd}F^{\prime}_{abcd}=-\frac{1}{2}\psi^{2}-\mathcal{V}(\phi).
\end{aligned}
\label{Fcomponentstrace}
\end{equation}
Using the expression \eqref{Fdual} with equations
\eqref{DecompositionWeylDual} and \eqref{Shoutenscalarfield}, we get
the following components of the dual ${^{\star}}\bmath{F}$:
\begin{equation}
 \begin{aligned}
  &{^{\star}}F^{\prime}_{00c0}=-{^{\star}}F^{\prime}_{000c}=0,\quad {^{\star}}F^{\prime}_{00cd}=-{^{\star}}F^{\prime}_{00dc}=0,\\
  &{^{\star}}F^{\prime}_{a0c0}=-{^{\star}}F^{\prime}_{a00c}=B_{ac}, \\
  &{^{\star}}F^{\prime}_{abc0}=-2E_{p[b}\epsilon_{a]}{}^p{}_c-\frac{1}{6}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)\epsilon_{bac}=-{^{\star}}F^{\prime}_{ab0c}, \\
  &{^{\star}}F^{\prime}_{a0cd}=E_{ap}\epsilon^{p}{}_{cd}-\frac{1}{6}\left(\frac{5}{2}\psi^{2}-\mathcal{V}(\phi)\right)\epsilon_{acd}, \\
  &{^{\star}}F^{\prime}_{0bcd}=-E_{bp}\epsilon^{p}{}_{cd}-\frac{1}{6}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)\epsilon_{bcd}, \\
  &{^{\star}}F^{\prime}_{0b0d}=B_{bd}, \quad {^{\star}}F^{\prime}_{abcd}=-B_{pq}\epsilon^{p}{}_{ab}\epsilon^{q}{}_{cd}.
 \end{aligned}
\label{Fdualcomponents}
\end{equation}

\subsection{The Bianchi equations}

If one substitutes the expressions for the Friedrich tensor derived in
the previous section into the first Bianchi identities
$(\ref{GeomDecBI})$ one obtains the following relations for the
components of the zero quantity $\tilde{F}_{abc}$:
\begin{equation}
\begin{aligned}
 \tilde{F}^{\prime}_{0c0}=&-\pounds_{\text{\tiny{$\bmath{N}$}}}F^{\prime}_{00c0}+D^{q}F^{\prime}_{q0c0}+\chi_{c}^{\,\,\,s}F^{\prime}_{00s0}-\chi F^{\prime}_{00c0}-\chi^{qb}\left(F^{\prime}_{qbc0}+F^{\prime}_{q0cb}\right)+a^{b}\left(F^{\prime}_{0bc0}+F^{\prime}_{00cb}+F^{\prime}_{b0c0}\right),\\
 \tilde{F}^{\prime}_{0cd}=&-\pounds_{\text{\tiny{$\bmath{N}$}}}F^{\prime}_{00cd}+D^{q}F^{\prime}_{q0cd}+a^{b}F^{\prime}_{0bcd}+a^{q}F^{\prime}_{q0cd}-\chi^{qb}F^{\prime}_{qbcd}-\chi F^{\prime}_{00cd}-\chi^{q}_{\,\,\,c}F^{\prime}_{q00d}-\chi^{q}_{\,\,\,d}F^{\prime}_{q0c0}& \\
		 &+\chi_{c}^{\,\,\,s}F^{\prime}_{00sd}+\chi_{d}^{\,\,\,s}F^{\prime}_{00cs}+a_{c}F^{\prime}_{000d}+a_{d}F^{\prime}_{00c0},\\
  \tilde{F}^{\prime}_{b0d}=&-\pounds_{\text{\tiny{$\bmath{N}$}}}F^{\prime}_{0b0d}+D^{a}F^{\prime}_{ab0d}-\chi F^{\prime}_{0b0d}-\chi^{a}_{\,\,\,b}F^{\prime}_{a00d}-\chi^{ac}F^{\prime}_{abcd}+\chi_{b}^{\,\,\,s}F^{\prime}_{0s0d}+a_{b}F^{\prime}_{000d}+\chi_{d}^{\,\,\,s}F^{\prime}_{0b0s} \\
		   &+a^{q}F^{\prime}_{qb0d}+a^{c}F^{\prime}_{0bcd},\\
 F^{\prime}_{bcd}=&-\pounds_{\text{\tiny{$\bmath{N}$}}}F^{\prime}_{0bcd}+D^{a}F^{\prime}_{abcd}+a^{q}F^{\prime}_{qbcd}+\left(a_{b}F^{\prime}_{00cd}+a_{c}F^{\prime}_{0b0d}+a_{d}F^{\prime}_{0bc0}\right)-\chi F^{\prime}_{0bcd} \\
		  &+\chi_{b}^{\,\,\,q}F^{\prime}_{0qcd}+\chi_{c}^{\,\,\,q}F^{\prime}_{0bqd}+\chi_{d}^{\,\,\,q}F^{\prime}_{0bcq}-\chi^{q}_{\,\,\,b}F^{\prime}_{q0cd}-\chi^{q}_{\,\,\,c}F^{\prime}_{qb0d}-\chi^{q}_{\,\,\,d}F^{\prime}_{qbc0},
\end{aligned}
\label{Fevolution}
\end{equation}
where we have used the fact that $\bmath{F}$ is anti-symmetric in the
last two indices ---see e.g. \cite{PhysRevD.57.2317}. Similar relations
hold for the dual ${^{\star}}\tilde{\bmath{F}}$.

\medskip
\noindent
\textbf{Remark.} In reference \cite{Friedrich_Rendall} ---cfr. equation
(4.47)--- suitable zero quantities are defined by using the
decomposition in terms of irreducible components of
$\tilde{\bmath{F}}$.

\subsubsection{The evolution equation for the electric part of the Weyl tensor}

An evolution equation for the electric part of the Weyl tensor can be
obtained using the third equation of \eqref{Fevolution} together with
the expressions \eqref{Fcomponents} and \eqref{Fcomponentstrace}, and
then symmetrising with respect to the indices $(bd)$. One obtains the
equation
\begin{equation*}
\begin{aligned}
 \tilde{F}^{\prime}_{(b|0|d)}=&-\pounds_{\text{\tiny{\bmath{N}}}}E_{bd}-\frac{1}{6}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)\pounds_{\text{\tiny{\bmath{N}}}}h_{bd}-\frac{1}{6}h_{bd}\pounds_{\text{\tiny{\bmath{N}}}}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)+D_{a}B_{p(d}\epsilon_{b)}{}^{pa} \\
		      &+2a_{a}B_{p(b}\epsilon_{d)}{}^{pa}-2\chi E_{bd}+2\chi^{a}{}_{(b}E_{d)q}+3\chi_{(b}{}^{q}E_{d)q}-h_{db}\chi^{ac}E_{ac}-\frac{1}{3}\left(\psi^{2}-\mathcal{V}(\phi)\right)\chi_{(bd)}.
\end{aligned}
\end{equation*}
Similarly, using equation \eqref{LieDer3metric} we get
\begin{equation*}
 \begin{aligned}
  \tilde{F}^{\prime}_{(b|0|d)}=&-\pounds_{\text{\tiny{\bmath{N}}}}E_{bd}+D_{a}B_{p(d}\epsilon_{b)}{}^{pa}+2a_{a}B_{p(b}\epsilon_{d)}{}^{pa}-2\chi E_{bd}+2\chi^{a}{}_{(b}E_{d)q}+3\chi_{(b}{}^{q}E_{d)q}-h_{db}\chi^{ac}E_{ac} \\
		       &-\frac{1}{2}\psi^{2}\chi_{(bd)}-\frac{1}{6}h_{bd}\pounds_{\text{\tiny{\bmath{N}}}}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right).
\end{aligned}
\end{equation*}
The trace of the previous expression is given by
\begin{equation*}
 h^{rs}\tilde{F}^{\prime}_{(r|0|s)}=-\frac{1}{2}\psi^{2}\chi-\frac{1}{2}\pounds_{\text{\tiny{\bmath{N}}}}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right),
\end{equation*}
which is the evolution equation for the scalar field ---i.e. the equation
expressing the conservation of energy. The evolution equation
for the components of the tensor $E_{ab}$ is obtained by taking the
difference of the last two equations so that:
\begin{equation*}
\begin{aligned}
 \pounds_{\text{\tiny{\bmath{N}}}}E_{bd}-D_{a}B_{p(d}\epsilon_{b)}{}^{pa}=&2a_{a}B_{p(b}\epsilon_{d)}{}^{pa}-2\chi E_{bd}+2\chi^{q}{}_{(b}E_{d)q}+3\chi_{(b}{}^{q}E_{d)q}-h_{db}\chi^{ac}E_{ac} \\
										&-\frac{1}{2}\psi^{2}\left(\chi_{(bd)}-\frac{1}{3}\chi h_{bd}\right).
\end{aligned}
\end{equation*}
Finally using (\ref{commutatorframes_scalarfield}) we can write the last equation as
\begin{equation}
 \pounds_{\text{\tiny{\bmath{N}}}}E_{bd}-D_{a}B_{p(d}\epsilon_{b)}{}^{pa}=-\frac{1}{2}\psi^{2}\left(\chi^{ST}\right)_{bd}-\frac{1}{3}\chi E_{bd}+5\left(\chi^{ST}\right)^{q}_{(b}E_{d)q}+2a_{a}B_{p(b}\epsilon_{d)}{}^{pa}.
 \label{evolutionE} \\
\end{equation}

\subsubsection{The evolution equation for the magnetic part of the Weyl tensor}

An evolution equation for the magnetic part of the Weyl tensor can be
derived from the third equation of \eqref{Fevolution} using the expressions
\eqref{Fdualcomponents}. A computation yields
\begin{equation*}
 \begin{aligned}
  \tilde{F}_{b0d}=&-\pounds_{\text{\tiny{\bmath{N}}}}B_{bd}+D^{a}\left(2E_{p[b}\epsilon_{a]}{}^p{}_d +\frac{1}{6}\left(\frac{1}{2}\psi^{2}+\mathcal{V}(\phi)\right)\epsilon_{bad}\right)-\chi B_{bd}+\chi^{a}{}_{b}B_{ad}+2\chi_{(b}{}^{a}B_{d)a} \\
          &+2a_{q}B_{pb}\epsilon^{qp}{}_{d}+a^{q}E_{pq}\epsilon^{p}{}_{bd}+\chi^{qb}B_{pq}\epsilon^{p}{}_{ab}\epsilon^{q}{}_{cd}.
 \end{aligned}
\end{equation*}
Now, since $B_{bd}$ is a symmetric tensor, all the information about
its evolution is contained in the symmetrised expression
$\tilde{F}_{(b|0|d)}$. Consequently, symmetrising the previous
equation with respect to the spatial indices $(bd)$ we find
\begin{equation*}
\quad\pounds_{\text{\tiny{\bmath{N}}}}B_{bd}-D_{a}E_{p(b}\epsilon_{d)}{}^{ap}=-2a_{a}E_{p(b}\epsilon_{d)}{}^{pa}+2\chi_{(b}{}^{a}B_{d)a}+\chi^{a}{}_{(d}B_{b)a}-\chi B_{bd}+\chi_{ac}B_{pq}\epsilon^{pa}{}_{(b}\epsilon_{d)}{}^{qc},
\end{equation*}
so that using (\ref{commutatorframes_scalarfield}) we have
\begin{equation}
\partial_{t}B_{bd}-D_{a}E_{p(b}\epsilon_{d)}{}^{ap}=-\frac{1}{3}\chi B_{bd}+3\left(\chi^{ST}\right)^{q}{}_{(b}B_{d)q}+\left(\chi^{ST}\right)_{ac}B_{pq}\epsilon^{pa}{}_{(b}\epsilon_{d)}{}^{qc}-2a_{a}E_{p(b}\epsilon_{d)}{}^{pa}.
\label{evolutionB}
\end{equation}
The principal part of the evolution equations \eqref{evolutionE} and \eqref{evolutionB}
is given by the symmetric matrix
\begin{equation*}
\begin{pmatrix}
			    \bmath{e}_{0}& 0 & 0 & 0 & -\frac{1}{2}D_{1} & \frac{1}{2}D_{2}  \\
			    0 & \bmath{e}_{0} & 0 & \frac{1}{2}D_{1} & 0 & -\frac{1}{2}D_{3}  \\
		            0 & 0 & \bmath{e}_{0} & -\frac{1}{2}D_{2} & \frac{1}{2}D_{3} & 0  \\
       		            0 & \frac{1}{2}D_{1} & -\frac{1}{2}D_{2} & \bmath{e}_{0} & 0 & 0  \\
		            -\frac{1}{2}D_{1} & 0 & \frac{1}{2}D_{3} & 0 & \bmath{e}_{0} & 0  \\
         		    \frac{1}{2}D_{2} & -\frac{1}{2}D_{3} & 0 & 0 & 0 & \bmath{e}_{0}
		    \end{pmatrix}.
\end{equation*}

\subsection{The Lagrangian description and  Fermi transport}

In order to deduce the remaining evolution equations, we will adopt a
Lagrangian description. This point of view amounts to requiring the
timelike vector of the orthonormal frame to follow the matter flow
lines. Accordingly, we introduce coordinates $(t,\bmath{x})$ such that
\begin{equation}
 \bmath{e}_{0}=\partial_{t}, \quad e_{0}{}^{\mu}=\delta^{\mu}_{0}.
\label{lagrangiandescription}
\end{equation}
This particular choice is equivalent to setting
$\bmath{\theta}^{b}=\theta^{b}{}_{j}\mbox{d}x^{j}$ while at the same
time fixing the lapse function to one\footnote{See also
\cite{2002TMNA_Bruhat_York}, where a symmetric hyperbolic system was
obtained for the Einstein-Euler system. This construction holds for an
arbitrary Eulerian picture.}.  With this choice (since
$\pounds_{N}=\partial_{t}$), we have from equations
(\ref{psi_{0}_frame}), (\ref{KGequation}), and
(\ref{constraint_Da_phi}) that
\begin{eqnarray}
&& \partial_{t}\phi=\psi \equiv N^{a}\psi_{a}=-\alpha<0,
\label{defipsi} \\
&& \partial_{t}\psi=-\psi\chi-\frac{d\mathcal{V}}{d\phi},
\label{evolutionpsi} \\
&&e^{\prime\,\,0}_{a}\psi=-e^{\prime\,\,j}_{a}D_{j}\phi.
\label{D_j_phi}
\end{eqnarray}

\medskip
Now, the timelike coframe is given in terms of the natural cobasis
through the relation
\begin{equation*}
 \bmath{\theta}^{0}=\mbox{d}t+\beta_{j}\mbox{d}x^{j},
\end{equation*}
while the spatial frame vectors are found to be
\begin{equation}
 \bmath{e}^{\prime}_{a}=(\theta_{a}{}^{j})^{-1}\left(\partial_{j}-\beta_{j}\partial_{t}\right), \quad e^\prime_{a}{}^{0}=(\theta_{a}{}^{j})^{-1}\beta_{j}, \quad e^\prime_{a}{}^{j}=(\theta_{a}{}^{j})^{-1}.
\label{spatialframe_coeff}
\end{equation}
It then follows from equations \eqref{D_j_phi} and \eqref{spatialframe_coeff} that
\begin{equation}
 \beta_{j}=-\frac{1}{\psi}D_{j}\phi.
\label{Beta_j}
\end{equation}
Thus, since $\beta_{j}$ is nonzero, the surfaces of constant time are
not necessarily spacelike for the characteristic cone and this could
be a problem for the hyperbolicity of the system ---see
\cite{2002TMNA_Bruhat_York,2002LNP...592...29C}. \\
Finally, the remaining frame components are chosen to be Fermi propagated along $\bmath{e}_{0}$. That is, we require
\begin{equation*}
 \nabla_{0}\bmath{e}_{a}-\left(\bmath{g}\left(\bmath{e}_{a},\nabla_{0}\bmath{e}_{0}\right)\bmath{e}_{0}-\bmath{g}\left(\bmath{e}_{a},\bmath{e}_{0}\right)\nabla_{0}\bmath{e}_{0}\right)=0.
\end{equation*}
The latter condition implies
\begin{equation*}
 \gamma^{\prime a}{}_{b0}=0.
\end{equation*}

\subsection{Evolution equation for the frame coefficients}

As already mentioned, the evolution equations are obtained from the
relation \eqref{commutatorframes} to yield
\begin{equation*}
 \left[\bmath{e}_{0},\bmath{e}^{\prime}_{b}\right]=a_{b}\bmath{e}_{0}-\gamma^{\prime c}{}_{0b}\bmath{e}^{\prime}_{c},
\end{equation*}
where it has been used that $\gamma^{\prime c}{}_{b0}=0$ (Fermi
gauge). Therefore, the evolution equations for the remaining frame
coefficients read
\begin{eqnarray*}
 && \partial_{t}e^\prime_{b}{}^{i}=-\chi_{b}{}^{c}e^\prime_{c}{}^{i}, \\
 && \partial_{t}e^\prime_{b}{}^{0}=a_{b}-\chi_{b}{}^{c}e^\prime_{c}{}^{0},
\end{eqnarray*}
which, together with \eqref{spatialframe_coeff}, imply propagation
equations for the components of the metric in the local coordinate
system.  In particular, one has that
\begin{equation*}
 \partial_{t}\beta_{j}=\theta^{\prime b}{}_{j}a_{b}
\end{equation*}
with $\beta_{j}$ given by equation \eqref{Beta_j} ---see also equation
(6.2) in \cite{2002TMNA_Bruhat_York} for arbitrary lapse $U$. Finally,
since from \eqref{commutatorframes_scalarfield} we have $(\chi^A)_{ab}=0$,
the evolution equations for the frame coefficients can be written as
\begin{equation}
\begin{aligned}
 \partial_{t}e^\prime_{b}{}^{i}&=-\frac{1}{3}\chi e^\prime_{b}{}^{i}-\left(\chi^{ST}\right)_{b}{}^{c}e^\prime_{c}{}^{i} \\
 \partial_{t}e^\prime_{b}{}^{0}&=-\frac{1}{3}\chi e^\prime_{b}{}^{0}-\left(\chi^{ST}\right)_{b}{}^{c}e^\prime_{c}{}^{0}+a_{b}.
\end{aligned}
\label{evolutionspatialframes}
\end{equation}

\subsection{Evolution equations for the connection coefficients}

The equations for the connection coefficients are obtained from the
splitting of the Riemann tensor with respect to the frame
$\{\bmath{e}_a\}$. In general, we have that
\begin{equation}
\begin{aligned}
R^{\prime a}{}_{b0d}&=\bmath{e}_{0}\left(\gamma^{a}{}_{bd}\right)-D_{d}\gamma^{a}{}_{b0}-a_{d}\gamma^{a}{}_{b0}
                              -\left(\gamma^{p}{}_{d0}-\chi_{d}{}^{p}\right)\gamma^{\prime a}{}_{bp}-a_{b}\chi_{d}{}^{a}+a^{a}\chi_{db}, \\
R^{\prime 0}{}_{b0d}&=\bmath{e}_{0}\left(\chi_{db}\right)-D_{d}a_{b}-a_{b}a_{d}+\chi_{pb}\chi_{d}{}^{p}-\chi_{dp}\gamma^{p}{}_{b0}-\chi_{pb}\gamma^{p}{}_{d0}, \\
R^{\prime a}{}_{0cd}&=D_{c}\chi_{d}{}^{a}-D_{d}\chi_{c}{}^{a}-a^{a}\left(\chi_{cd}-\chi_{dc}\right), \\
R^{\prime a}{}_{bcd}&=\tilde{R}^{a}{}_{bcd}+\chi_{c}{}^{a}\chi_{db}-\chi_{d}{}^{a}\chi_{cb}-\gamma^{a}{}_{b0}\left(\chi_{cd}-\chi_{dc}\right),
\end{aligned}
\label{Riemann_connection}
\end{equation}
where $\tilde{R}^{a}{}_{bcd}$ denotes the Riemann tensor
constructed only with the spatial connection coefficients
$\gamma^{\prime c}{}_{ab}$.  The first two identities give
evolution equations once the Lagrangian gauge is introduced. The
remaining two equations are the \emph{quasi-constraints} for the
connection coefficients ---see
\cite{2002TMNA_Bruhat_York,2002LNP...592...29C}.  No equations for the
connection coefficient associated to the acceleration can be deduced
from these identities. In the sequel, it will be shown how evolution
equations for the acceleration can be obtained for our particular
problem.

\medskip
From equations \eqref{Riemann_connection}, we can also deduce two
important equations relating the Ricci tensor to the connection:
\begin{equation}
 \begin{aligned}
  R_{00}&=-\bmath{e}_{0}\left(\chi\right)+D_{p}a^{p}-\chi_{p}{}^{b}\chi_{b}{}^{p}+a_{p}a^{p} \\
  R^\prime_{0d}&=D_{c}\chi_{d}{}^{c}-D_{d}\chi-2\,a^{c}\left(\chi^{A}\right)_{cd}.
 \end{aligned}
\label{Ricci_connection}
\end{equation}

\medskip
The first identity in \eqref{Riemann_connection}, together with the
conditions for the Lagrangian and Fermi gauge provide the equation
\begin{equation*}
\partial_{t}\gamma^{a}{}_{bd}=-\gamma^{a}{}_{bp}\chi_{d}{}^{p}+2h^{ap}\chi_{d[p}a_{b]}+B_{dp}\epsilon^{pa}{}_{b},
\end{equation*}
describing the evolution of the spatial connection coefficients
$\gamma^{\prime c}{}_{ab}$. To obtain the last equation it has been
used that $R^{\prime a}{}_{b0d}=C^{\prime a}{}_{b0d}=B_{dp}\epsilon^{pa}{}_{b}$.
Moreover, since from \eqref{commutatorframes_scalarfield}
we have $(\chi^A)_{ab}=0$, we can write
\begin{equation}
\partial_{t}\gamma^{a}{}_{bd}=-\frac{1}{3}\chi\gamma^{a}{}_{bd}+\frac{2}{3}\chi h^{ap}h_{d[p}a_{b]}+B_{dp}\epsilon^{pa}{}_{b}-\left(\chi^{ST}\right)_{d}{}^{p}\gamma^{a}{}_{bp}+2h^{ap}\left(\chi^{ST}\right)_{d[p}a_{b]}.
\label{evolutionspatial_connectioncoefficients}
\end{equation}

The evolution equation for the part of the connection described by $\chi_{bd}$ is obtained from the second identity
in \eqref{Riemann_connection}. In order to do so, first, we will derive the evolution and the quasi-constraint equations
for the acceleration. The evolution equation for the acceleration can be obtained from
\begin{equation*}
\begin{aligned}
 \left[\bmath{e}_{0},\bmath{e}^{\prime}_{c}\right]\psi&=c^{0}{}_{0c}\bmath{e}_{0}(\psi)+c^{p}{}_{0c}\bmath{e}^{\prime}{}_{p}(\psi) \\
                                                      &=\gamma^{0}{}_{c0}(\partial_{t}\psi)+\left(\gamma^{p}{}_{c0}-\gamma^{p}{}_{0c}\right)(D_{p}\psi) \\
                                                      &=a_{c}\left(\partial_{t}\psi\right)-\chi_{c}{}^{p}\left(D_{p}\psi\right)  ,
\end{aligned}
\end{equation*}
where the Lagrangian and Fermi gauge have been employed. Now,
expanding the left hand side and making use of the evolution and the
quasi-constraint equation for the energy-momentum tensor of the scalar
field one has that
\begin{equation*}
 \partial_{t}a_{c}-D_{c}\chi=-\chi_{c}^{\,\,\,p}a_{p}+\left(\chi+\frac{2}{\psi}\frac{d\mathcal{V}}{d\phi}\right)a_{c}
\end{equation*}
so that using the second equation in \eqref{Ricci_connection} we arrive at
\begin{equation}
 \partial_{t}a_{c}-D_{p}\chi_{c}{}^{p}=-2\left(\chi^{A}\right)_{pc}a^{p}-\chi_{c}{}^{p}a_{p}+\left(\chi+\frac{2}{\psi}\frac{d\mathcal{V}}{d\phi}\right)a_{c}.
\label{evolution_a}
\end{equation}
Finally, since for our particular problem one has that $(\chi^A)_{ab}=0$, this last equation takes the form
\begin{equation}
\quad\partial_{t}a_{c}-D_{p}\chi_{c}{}^{p}=\left(\frac{2}{\psi}\frac{d\mathcal{V}}{d\phi}+\frac{2}{3}\chi\right)a_{c}-\left(\chi^{ST}\right)_{c}{}^{b}a_{b}.
\label{evolutionA}
\end{equation}
In the case of the quasi-constraint a computation yields
\begin{equation}
D_{c}a_{b}-D_{b}a_{c}=2\left(\chi+\frac{1}{\psi}\frac{d\mathcal{V}}{d\phi}\right)\left(\chi^{A}\right)_{cb}.
\end{equation}
Thus, making use of this equation in the second identity in
equation \eqref{Riemann_connection} and recalling the properties of the Fermi
gauge one finds that
\begin{equation}
\partial_{t}\chi_{db}-D_{b}a_{d}=E_{db}+\frac{1}{3}\left(\mathcal{V}(\phi)-\psi^{2}\right)h_{db}-\chi_{d}{}^{p}\chi_{pb}+2\left(\chi+\frac{1}{\psi}\frac{d\mathcal{V}}{d\phi}\right)\left(\chi^{A}\right)_{db}+a_{d}a_{b},
\label{evolution_chi}
\end{equation}
where it has been used
\[
R^{\prime 0}{}_{b0d}=C^{\prime
  0}{}_{b0d}+\frac{1}{2}\delta^{0{}}_{0}S_{bd}+\frac{1}{2}h_{bd}S_{0}{}^{0}=E_{bd}+\frac{1}{3}\left(\mathcal{V}(\phi)-\psi^{2}\right)h_{bd}.
\]
The principal part of the combined system of equations
\eqref{evolution_a} and \eqref{evolution_chi} is given by
\begin{equation*}
 \begin{pmatrix}
  \bmath{e}_{0} & -D_{1} & -D_{2} & -D_{3} \\
  -D_{1} & \bmath{e}_{0} & 0 & 0 \\
  -D_{2} & 0 & \bmath{e}_{0} & 0 \\
  -D_{3} & 0 & 0 & \bmath{e}_{0}
 \end{pmatrix}
 \begin{pmatrix}
  a_{d} \\
  \chi_{d}{}^{1} \\
  \chi_{d}{}^{2} \\
  \chi_{d}{}^{3}
 \end{pmatrix}.
\end{equation*}

Since $(\chi^A)_{ab}=0$, then $\chi_{bd}$ is symmetric and in this case, we have after symmetrisation that
\begin{equation}
2\partial_{t}\chi_{(bd)}-2D_{(b}a_{d)}
				    =\frac{2}{3}\left(\mathcal{V}(\phi)-\psi^{2}-\frac{1}{3}\chi^{2}\right)h_{bd}-\frac{4}{3}\chi\left(\chi^{ST}\right)_{bd}-2\left(\chi^{ST}\right)_{b}{}^{p}\left(\chi^{ST}\right)_{pd}+2a_{b}a_{d}+2E_{bd}.
\label{evolutionsecondfundform}
\end{equation}
Also, from the first equation in \eqref{Ricci_connection} it follows that
\begin{equation*}
 \partial_{t}\chi-D_{p}a^{p}=\mathcal{V}(\phi)-\psi^{2}-\frac{\chi^2}{3}-\left(\chi^{ST}\right)^{2}+a^{2},
\end{equation*}
where $(\chi^{ST})^{2}=(\chi^{ST})_{ab}(\chi^{ST})^{ab}$ and $a^{2}=a_{c}a^{c}$.

\subsection{Hyperbolicity considerations}

The system consisting of equations \eqref{defipsi},
\eqref{evolutionpsi}, \eqref{evolutionE}, \eqref{evolutionB},
\eqref{evolutionspatialframes},
\eqref{evolutionspatial_connectioncoefficients}, \eqref{evolutionA}
and \eqref{evolutionsecondfundform} can be written matricially as
\begin{equation}
 \mathbf{A}^{0}\partial_{t}\mathbf{u}-\mathbf{A}^{p}\bmath{e}^{\prime}_{p}(\mathbf{u})=\mathbf{B}(\mathbf{u})\mathbf{u}.
\label{framederiv_syst}
\end{equation}
As discussed in \cite{2002TMNA_Bruhat_York}, these systems are not
hyperbolic in the usual sense as, in general, the time lines are not
hypersurface orthogonal and the ``spatial'' frame vectors
$\bmath{e}^{\prime}_{a}$ have  components in the time direction ---see
\eqref{spatialframe_coeff}. Since the surfaces of constant time $t$
are not necessarily spacelike, this type of system is referred to as a
{\it quasi FOSH} system ---see  \cite{2002TMNA_Bruhat_York}. In terms
of the partial derivatives, equation \eqref{framederiv_syst} reads
\begin{equation}
 \tilde{\mathbf{A}}^{0}(\mathbf{u})\partial_{t}\mathbf{u}-\mathbf{A}^{j}(\mathbf{u})\partial_{j}\mathbf{u}=\mathbf{B}(\mathbf{u})\mathbf{u}
\label{partialderiv_syst}
\end{equation}
with
\begin{equation}
\tilde{\mathbf{A}}^{0}(\mathbf{u})\equiv
\mathbf{A}^{0}-\mathbf{A}^{p}e_{p}^{\prime}{}^{0}, \quad \mathbf{A}^{j}(\mathbf{u})\equiv\mathbf{A}^{p}e_{p}^{\prime}{}^{j}.
\label{A0_scalarfield}
\end{equation}

\medskip
In order to have a well posed initial value problem, the matrix $\tilde{\mathbf{A}}^{0}(\mathbf{u})$ must be positive definite.
This is the case as long the quadratic form
\begin{equation}
\sum_{b=1,2,3}\theta^{b}{}_{i}\theta^{b}{}_{j} -\beta_{i}\beta_{j}
\label{positive_quadratic_form}
\end{equation}
is positive definite -- see {\it Proposition 9} in
\cite{2002TMNA_Bruhat_York}.  In the next section we will consider a
reference solution admitting a foliation by homogeneous spacelike
hypersurfaces. As a consequence,  the linearisation of the system \eqref{nonlinearequation} is well posed
without the need to control the smallness of $\beta_{i}$. The
smallness of these terms is taken care by the perturbation fields ---see also \cite{PhysRevD.60.083507}.

\medskip
Written in terms of partial derivatives, our system of evolution equations reads
\begin{equation}
 \begin{aligned}
 \partial_{t}\phi&=\psi, \\
  \partial_{t}\psi&=-\psi\chi-\frac{d\mathcal{V}}{d\phi}, \\
 2\,\partial_{t}\chi_{(bd)}-2\,e_{(b}^{\prime}{}^{0}\partial_{t}a_{d)}-2\,e_{(b}^{\prime}{}^{j}\partial_{j}a_{d)}&=\frac{2}{3}\left(\mathcal{V}(\phi)-\psi^{2}-\frac{1}{3}\chi^{2}\right)h_{bd}-\left(\gamma^{p}{}_{bd}+\gamma^{p}{}_{db}\right)a_{p} \\
                                                                                &\quad-\frac{4}{3}\chi\left(\chi^{ST}\right)_{bd}-2\left(\chi^{ST}\right)_{b}{}^{p}\left(\chi^{ST}\right)_{pd}+2a_{b}a_{d}+2E_{bd},\\
\partial_{t}a_{c}-e_{p}^{\prime}{}^{0}\partial_{t}\chi_{c}{}^{p}-e_{p}^{\prime}{}^{j}\partial_{j}\chi_{c}{}^{p}&=\left(\frac{2}{\psi}\frac{d\mathcal{V}}{d\phi}+\frac{2}{3}\chi\right)a_{c}-\left(\chi^{ST}\right)_{c}{}^{b}a_{b} \\
                                                                                &\quad-\gamma^{q}{}_{cp}(\chi^{ST}){}^{p}{}_{q}-\gamma^{q}{}_{p}{}^p(\chi^{ST})_{cq}, \\
 \partial_{t}E_{bd}-\epsilon^{pa}{}_{(b|}e^{\prime}_{a}{}^0\partial_{t}B_{p|d)}-\epsilon^{pa}{}_{(b|}e^{\prime}_{a}{}^j\partial_{j}B_{p|d)}&=-\frac{1}{2}\psi^{2}\left(\chi^{ST}\right)_{bd}-\frac{1}{3}\chi E_{bd}+5\left(\chi^{ST}\right)^{q}_{(b}E_{d)q} \\
                                                                                &\quad+2a_{a}B_{p(b}\epsilon_{d)}{}^{pa}-\gamma^{q}{}_{pa}B_{q(d}\epsilon_{b)}{}^{pa}-\epsilon^{pa}{}_{(d}\gamma^{q}{}_{b)a}B_{pq}, \\
 \partial_{t}B_{bd}-\epsilon^{ap}{}_{(d}e_{|a}^{\prime}{}^{0}\partial_{t|}E_{b)p}-\epsilon^{ap}{}_{(d}e_{|a}^{\prime}{}^{j}\partial_{j|}E_{b)p}&=-\frac{1}{3}\chi B_{bd}+3\left(\chi^{ST}\right)^{q}_{(b}B_{d)q}+\left(\chi^{ST}\right)_{ac}B_{pq}\epsilon^{pa}{}_{(b}\epsilon_{d)}{}^{qc}\\
									   &\quad-2a_{a}E_{p(b}\epsilon_{d)}{}^{pa}-\gamma^{q}{}_{pa}E_{q(b}\epsilon_{d)}{}^{ap}-\epsilon^{ap}{}_{(b}\gamma^{q}{}_{d)a}E_{pq}, \\
 \partial_{t}e_{b}^{\prime}{}^{i}&=-\frac{1}{3}\chi e_{b}^{\prime}{}^{i}-\left(\chi^{ST}\right)_{b}{}^{c}e_{c}^{\prime}{}^{i},\\
 \partial_{t}e_{b}^{\prime}{}^{0}&=-\frac{1}{3}\chi e^\prime_{b}{}^{0}-\left(\chi^{ST}\right)_{b}{}^{c}e^\prime_{c}{}^{0}+a_{b}, \\
 \partial_{t}\gamma^{\prime a}{}_{bd}&=-\frac{1}{3}\chi\gamma^{\prime a}{}_{bd}+\frac{2}{3}\chi h^{ap}h_{d[p}a_{b]}+B_{dp}\epsilon^{pa}{}_{b}-\left(\chi^{ST}\right)_{d}{}^{p}\gamma^{\prime a}{}_{bp} \\
                                     &\quad+2h^{ap}\left(\chi^{ST}\right){}_{d[p}a_{b]}.
\end{aligned}
\label{EFsf}
\end{equation}

\medskip
This system has clearly the form given by equation \eqref{partialderiv_syst}. If one writes
\begin{equation}
 \mathbf{u}^{T}=\left(\phi,\,\,\psi,\,\,\mathbf{z}^{T},\,\,\mathbf{w}^{T},\,\,\mathbf{x}^{T},\,\,\mathbf{y}^{T}\right),
\label{u}
\end{equation}
where
\begin{equation*}
\begin{aligned}
 \mathbf{z}^{T}&=\left(\chi_{11},\,\,\chi_{22},\,\,\chi_{33},\,\,\left(\chi^{ST}\right)_{12},\,\,\left(\chi^{ST}\right)_{13},\,\,\left(\chi^{ST}\right)_{23},\,\,a_{1},\,\,a_{2},\,\,a_{3}\right), \\
 \mathbf{w}^{T}&=\left(E_{12},\,\,E_{13},\,\,E_{23},\,\,B_{12},\,\,B_{13},\,\,B_{23}\right), \\
 \mathbf{x}^{T}&=\left(e_{1}{}^{0},\,\,e_{2}{}^{0},\,\,e_{3}{}^{0},\,\,e_{1}{}^{1},\,\,e_{1}{}^{2},\,\,e_{1}{}^{3},\,\,e_{2}{}^{1},\,\,e_{2}{}^{2},\,\,e_{2}{}^{3},\,\,e_{3}{}^{1},\,\,e_{3}{}^{2},\,\,e_{3}{}^{3}\right),  \\
 \mathbf{y}^{T}&=\left(\gamma^{1}{}_{22},\,\,\gamma^{1}{}_{33},\,\,\gamma^{1}{}_{23},\,\,\gamma^{2}{}_{11},\,\,\gamma^{2}{}_{33},\,\,\gamma^{2}{}_{31},\,\,\gamma^{3}{}_{11},\,\,\gamma^{3}{}_{22},\,\,\gamma^{3}{}_{12}\right),
\end{aligned}
\end{equation*}
then the matrices given in equations \eqref{partialderiv_syst} and \eqref{A0_scalarfield} have the explicit form
\begin{equation}
 \tilde{\mathbf{A}}^{0}(\mathbf{u})=\begin{pmatrix}
			      \mathbf{I}_{2\times2} & 0 & 0 & 0\\
			      0 & \tilde{\mathbf{A}}^{0}_{9\times9} & 0 & 0 \\
                              0 & 0 &  \tilde{\mathbf{A}}^{0}_{6\times6} & 0 \\
                              0 & 0 & 0 & \mathbf{I}_{21\times21}
			    \end{pmatrix}, \quad \quad
\mathbf{A}^{j}(\mathbf{u})=\begin{pmatrix}
			      \mathbf{0}_{2\times2} & 0 & 0 & 0\\
			      0 & \mathbf{A}^{j}_{9\times9} & 0 & 0 \\
                              0 & 0 &  \mathbf{A}^{j}_{6\times6} & 0 \\
                              0 & 0 & 0 & \mathbf{0}_{21\times21}
			   \end{pmatrix},
\label{A_matrices}
\end{equation}
with

\begin{equation*}
\tilde{\mathbf{A}}^{0}_{9\times9}=\begin{pmatrix}
	             1 & 0 & 0 & 0 & 0 & 0 & -e_{1}^{\,\,\,0} & 0 & 0 \\
		     0 & 1 & 0 & 0 & 0 & 0 & 0 & -e_{2}^{\,\,\,0} & 0 \\
		     0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -e_{3}^{\,\,\,0} \\
		     0 & 0 & 0 & 2 & 0 & 0 & -e_{2}^{\,\,\,0} & -e_{1}^{\,\,\,0} & 0 \\
		     0 & 0 & 0 & 0 & 2 & 0 & -e_{3}^{\,\,\,0} & 0 & -e_{1}^{\,\,\,0} \\
		     0 & 0 & 0 & 0 & 0 & 2 & 0 & -e_{3}^{\,\,\,0} & -e_{2}^{\,\,\,0} \\
		     -e_{1}^{\,\,\,0} & 0 & 0 &  -e_{2}^{\,\,\,0} & -e_{3}^{\,\,\,0} & 0 &  1 & 0 & 0 \\
		     0 & -e_{2}^{\,\,\,0} & 0 &  -e_{1}^{\,\,\,0} & 0 & -e_{3}^{\,\,\,0} & 0 & 1 & 0 \\
		     0 & 0 & -e_{3}^{\,\,\,0} & 0 & -e_{1}^{\,\,\,0} & -e_{2}^{\,\,\,0} &  0 & 0 & 1 \\
		    \end{pmatrix},
\end{equation*}
\\
\begin{equation*}
\tilde{\mathbf{A}}^{0}_{6\times6}=\begin{pmatrix}
			    1 & 0 & 0 & 0 & -\frac{1}{2}e_{1}^{\,\,\,0} & \frac{1}{2}e_{2}^{\,\,\,0}  \\
			    0 & 1 & 0 & \frac{1}{2}e_{1}^{\,\,\,0} & 0 & -\frac{1}{2}e_{3}^{\,\,\,0}  \\
		            0 & 0 & 1 & -\frac{1}{2}e_{2}^{\,\,\,0} & \frac{1}{2}e_{3}^{\,\,\,0} & 0  \\
       		            0 & \frac{1}{2}e_{1}^{\,\,\,0} & -\frac{1}{2}e_{2}^{\,\,\,0} & 1 & 0 & 0  \\
		            -\frac{1}{2}e_{1}^{\,\,\,0} & 0 & \frac{1}{2}e_{3}^{\,\,\,0} & 0 & 1 & 0  \\
         		    \frac{1}{2}e_{2}^{\,\,\,0} & -\frac{1}{2}e_{3}^{\,\,\,0} & 0 & 0 & 0 & 1
		    \end{pmatrix},
\end{equation*}
\\
\begin{equation*}
\mathbf{A}^{j}_{9\times9}=\begin{pmatrix}
	             0 & 0 & 0 & 0 & 0 & 0 & e_{1}^{\,\,\,j} & 0 & 0 \\
		     0 & 0 & 0 & 0 & 0 & 0 & 0 & e_{2}^{\,\,\,j} & 0 \\
		     0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & e_{3}^{\,\,\,j} \\
		     0 & 0 & 0 & 0 & 0 & 0 & e_{2}^{\,\,\,j} & e_{1}^{\,\,\,j} & 0 \\
		     0 & 0 & 0 & 0 & 0 & 0 & e_{3}^{\,\,\,j} & 0 & e_{1}^{\,\,\,j} \\
		     0 & 0 & 0 & 0 & 0 & 0 & 0 & e_{3}^{\,\,\,j} & e_{2}^{\,\,\,j} \\
		     e_{1}^{\,\,\,j} & 0 & 0 &  e_{2}^{\,\,\,j} & e_{3}^{\,\,\,j} & 0 &  0 & 0 & 0 \\
		     0 & e_{2}^{\,\,\,j} & 0 &  e_{1}^{\,\,\,j} & 0 & e_{3}^{\,\,\,j} & 0 & 0 & 0 \\
		     0 & 0 & e_{3}^{\,\,\,j} & 0 & e_{1}^{\,\,\,j} & e_{2}^{\,\,\,j} &  0 & 0 & 0 \\
		    \end{pmatrix},
\end{equation*}
and
\begin{equation*}
\mathbf{A}^{j}_{6\times6}=\begin{pmatrix}
			    0& 0 & 0 & 0 & \frac{1}{2}e_{1}^{\,\,\,j} & -\frac{1}{2}e_{2}^{\,\,\,j}  \\
			    0 & 0 & 0 & -\frac{1}{2}e_{1}^{\,\,\,j} & 0 & \frac{1}{2}e_{3}^{\,\,\,j}  \\
		            0 & 0 & 0 & \frac{1}{2}e_{2}^{\,\,\,j} & -\frac{1}{2}e_{3}^{\,\,\,j} & 0  \\
       		            0 & -\frac{1}{2}e_{1}^{\,\,\,j} & \frac{1}{2}e_{2}^{\,\,\,j} & 0 & 0 & 0  \\
		            \frac{1}{2}e_{1}^{\,\,\,j} & 0 & -\frac{1}{2}e_{3}^{\,\,\,j} & 0 & 0 & 0  \\
         		    -\frac{1}{2}e_{2}^{\,\,\,j} & \frac{1}{2}e_{3}^{\,\,\,j} & 0 & 0 & 0 & 0
		    \end{pmatrix}.
\end{equation*}

\medskip
It can, therefore, be verified that $\tilde{\mathbf{A}}^{0}(\mathbf{u})$, $\mathbf{A}^{j}(\mathbf{u})$
are symmetric and that, furthermore, $\tilde{\mathbf{A}}^{0}(\mathbf{u})$ is positive definite as long
as \eqref{positive_quadratic_form} is satisfied.
We summarise the results of this section in the following theorem:
\begin{thm}
\label{Theorem:wellposedness}
The Einstein-Friedrich-nonlinear scalar field (EFsf) system
consisting of  the equations in \eqref{EFsf}
forms a quasi-linear first-order symmetric
hyperbolic (FOSH) system for the scalar field, its momentum-density,
the frame coefficients, the connection coefficients and the electric
and magnetic parts of the Weyl tensor, relatively to the slices of constant time $t$,
as long as the quadratic form
\begin{equation*}
\sum_{a=1,2,3}\theta^{a}{}_{i}\theta^{a}{}_{j}-\frac{D_{i}\phi}{\psi}\frac{D_{j}\phi}{\psi},
\end{equation*}
is positive definite.
\end{thm}

Using the standard theory of symmetric hyperbolic systems one can then
conclude the local existence in time and uniqueness of smooth solutions
for the evolution equations implied by  the Einstein-nonlinear scalar
field. In order to conclude the existence of solutions to the full
Einstein-scalar field equations one has to verify that the constraint
equations are satisfied during the evolution if these hold
initially. This will not be discussed here as the purpose of the present work  is to study the evolution
of the perturbations. In any case, it follows by general arguments that these are preserved during
evolution ---see e.g.  \cite{PhysRevD.60.083507,lrr-1998-3,Friedrich_Rendall}.

\medskip
\noindent
\textbf{Remark.} As part of the procedure to make the system
\eqref{EFsf} explicitly symmetric hyperbolic one has to divide the evolution equation
for the acceleration \eqref{evolutionA} by  $\psi^{2}$. This could
imply that the system is not well behaved when $\psi=0$. An inspection
shows that the potentially troublesome term is the one containing the first derivative of the
potential, which,  by virtue of the EFE must be zero in the asymptotic limit. Thus, the coefficient
\begin{equation}
 \lim_{t\rightarrow\infty}\frac{\mathcal{V}^{\prime}}{\psi}= const.
\label{AsymptoticCondition}
\end{equation}
must be required to exist. Moreover, we notice that our gauge
restricts the set of solutions, since it only admits scalar fields
$\phi$ which are monotonically decaying functions. Accordingly, the
stability of background solutions for which $\phi\rightarrow\infty$
cannot be studied by our method.

\section{Stability Analysis}

In this section we use the symmetric hyperbolic system derived in
last section to show that, for some classes of potentials, the
evolution of sufficiently small
nonlinear perturbations of a FRW-nonlinear scalar field background
prescribed on a Cauchy hypersurface with the topology of a 3-torus $\mathbb{T}^{3}$, either converge to
constant values or have an asymptotic exponential decay.

\subsection{The background solution}
As it is well known, the metric of a Friedman-Robertson-Walker (FRW)
spacetime ---i.e. an homogeneous and isotropic spacetime---  can be written as
\begin{equation*}
 \mbox{d}s^{2}=-\mbox{d}t^{2}+\left(\frac{a(t)}{\omega}\right)^{2}\delta_{ij}
 \mbox{d}x^{i}\mbox{d}x^{j},
\end{equation*}
where
\begin{equation*}
 \omega^{2}=1+\frac{k}{4}\delta_{ij}x^{i}x^{j}, \quad \quad\partial_{k}\omega^{2}=\frac{k}{2}x_{k},
\end{equation*}
and the constant $k=-1,0,1$ is the curvature of the spatial hypersurfaces.
Since the metric is conformal flat, it follows that
\[
\mathring{E}_{bd}=\mathring{B}_{bd}=0.
\]
Furthermore, if one foliates the spacetime with the surfaces of
constant $t$, one has that
\[
\left(\mathring{\chi}^{ST}\right)_{bd}=\left(\mathring{\chi}^{A}\right)_{bd}=0,
\quad \mathring{a}_{c}=0.
\]
For such background metrics
\[
\mathring{\chi}=3\dot{a}/a\equiv 3H,
\]
 where $\dot{}$ denotes differentiation with respect to time t, and $H$ is the so-called
\emph{Hubble function}.
In the case of a FRW cosmology the
Einstein-scalar field system reduces to the evolution equations
\begin{equation}
\begin{aligned}
 \frac{\mbox{d}\mathring{\phi}}{\mbox{d}t}&=\mathring{\psi}, \\
 \frac{\mbox{d}\mathring{\psi}}{\mbox{d}t}&=-3H\mathring{\psi}-\frac{\mbox{d}\mathcal{V}}{\mbox{d}\mathring{\phi}}, \\
 \frac{\mbox{d}H}{\mbox{d}t}&=-H^{2}-\frac{1}{3}\mathring{\psi}^{2}+\frac{1}{3}\mathcal{V}(\mathring{\phi}),
\end{aligned}
\label{evol_background}
\end{equation}
subject to the Friedmann-scalar field constraint equation
\begin{equation}
 {H}^{2}(t)-\frac{1}{6}\mathring{\psi}^{2}(t)-\frac{1}{3}\mathcal{V}(\mathring{\phi})=-\frac{k}{{a}^{2}} .
\label{Friedmannequation}
\end{equation}
Now, the gauge conditions for the frame are satisfied if one sets
\begin{equation*}
 \mathring{e}_{0}{}^{\mu}=\delta_{0}{}^\mu,\quad\mathring{e}_{b}{}^\mu={\left(\frac{\omega}{a}\right)}\delta_{b}{}^\mu,\quad b=1,2,3,
\end{equation*}
so that the spatial connection coefficients are given by
\begin{equation*}
 \mathring{\gamma}^{c}{}_{bd}=\frac{k}{4a^{2}}\left(h_{db}x^{c}-h_{d}{}^{c}x_{b}\right),\quad b,c,d=1,2,3,
\end{equation*}
with $x^{\mu}=(\omega/a)\delta^{\mu}{}_{c}x^{c}$. The remaining nonvanishing connection coefficients are
\begin{equation*}
 \mathring{\gamma}^{0}{}_{bd}=Hh_{bd} \quad , \quad \mathring{\gamma}^{b}{}_{0d}=Hh^{b}{}_{d},\quad b,d=1,2,3.
\end{equation*}

\subsection{Linearised evolution equations}

In this subsection, we derive the linearised system corresponding to
the nonlinear equations of Theorem \ref{Theorem:wellposedness}, for a
FRW background with a selfinteracting scalar field. In order to
perform the linearisation procedure we compute
\begin{equation*}
\left.\frac{d\breve{\mathbf{u}}^{\epsilon}}{d\epsilon}\right |_{\epsilon=0}
\end{equation*}
and drop all (nonlinear) terms of coupled perturbations.
In this way, we obtain the following system
\begin{equation}
 \begin{aligned}
    \partial_{t}\breve{\phi}&=\breve{\psi}, \\
    \partial_{t}\breve{\psi}&=-\left({\frac{\mbox{d}^{2}\mathcal{\mathring{V}}}{\mbox{d}\mathring{\phi}^{2}}}\right)\breve{\phi}-3H\breve{\psi}-\mathring{\psi}\breve{\chi},\\
   2\partial_{t}\breve{\chi}_{(bd)}-2\left(\frac{\omega}{{a}}\right)\delta_{(d}{}^{j}\partial_{j}\breve{a}_{b)}&=\frac{2}{3}\left(\left({\frac{\mbox{d}\mathcal{\mathring{V}}}{\mbox{d}\mathring{\phi}}}\right)\breve{\phi}-2\mathring{\psi}\breve{\psi}-2H\breve{\chi}\right)h_{bd}-4H\left(\breve{\chi}^{ST}\right)_{bd} \\
                                                                                &\quad-\left(\mathring{\gamma}^{p}{}_{bd}+\mathring{\gamma}^{p}{}_{db}\right)\breve{a}_{p}+2\breve{E}_{bd},\\
    \partial_{t}\breve{E}_{bd}-\left(\frac{\omega}{{a}}\right)\epsilon^{pa}{}_{(b|}\delta_{a}{}^{j}\partial_{j}\breve{B}_{p|d)}&=-\frac{1}{2}\mathring{\psi}^{2}(\breve{\chi}^{ST})_{bd}-H\breve{E}_{bd}-\mathring{\gamma}^{q}{}_{pa}\breve{B}_{q(d}\epsilon_{b)}{}^{pa}-\epsilon^{pa}{}_{(d|}\mathring{\gamma}^{q}{}_{|b)a}\breve{B}_{pq}, \\
    \partial_{t}\breve{B}_{bd}-\left(\frac{\omega}{{a}}\right)\epsilon_{(d|}{}^{ap}\delta_{a}{}^{j}\partial_{j}\breve{E}_{p|b)}&=-\mathring{\gamma}^{q}{}_{pa}\breve{E}_{q(b}\epsilon_{d)}{}^{ap}-\epsilon^{ap}{}_{(b}\mathring{\gamma}^{q}{}_{d)a}\breve{E}_{pq}-H\breve{B}_{bd} ,\\
    \partial_{t}\breve{a}_{c}-\left(\frac{\omega}{{a}}\right)\delta_{p}{}^{j}\partial_{j}\breve{\chi}_{c}{}^p&=-\mathring{\gamma}^{q}{}_{cp}(\breve{\chi}^{ST})_{q}{}^{p}-\mathring{\gamma}^{q}{}_p{}^p(\breve{\chi}^{ST})_{cq}+\left(2H+\frac{2}{\mathring{\psi}}{\frac{\mbox{d}\mathcal{\mathring{V}}}{\mbox{d}\mathring{\phi}}}\right)\breve{a}_{c}+\left(\frac{\mbox{d}H}{\mbox{d}t}\right)\breve{e}_{c}{}^0 ,\\
    \partial_{t}\breve{e}_{b}{}^{0}&=-H\breve{e}_{b}{}^{0}+\breve{a}_{b},\\
    \partial_{t}\breve{e}_{b}{}^{i}&=-\frac{1}{3}\left(\frac{\omega}{a}\right)\delta_{b}{}^{i}\breve{\chi}-\left(\frac{\omega}{a}\right)\delta_{c}{}^{i}\left(\breve{\chi}^{ST}\right)_{b}{}^{c}-H\breve{e}_{b}{}^{i},\\
    \partial_{t}\breve{\gamma}^{a}{}_{bd}&=-\frac{1}{3}\mathring{\gamma}^{a}{}_{bd}\breve{\chi}-\mathring{\gamma}^{a}{}_{bp}\left(\breve{\chi}^{ST}\right)_{d}{}^{p}+H\left(\delta^{a}_{d}\breve{a}_{b}-h_{bd}\breve{a}^{a}\right) \\
                                         &\quad+\breve{B}_{dp}\epsilon^{pa}{}_{b}-H\breve{\gamma}^{a}{}_{bd}.
\end{aligned}
\label{linearized_syst}
\end{equation}
\medskip
As a consequence, the linearised system has the following form
\begin{equation}
\mathring{\mathbf{A}}^{0}\partial_{t}\breve{\mathbf{u}}-\mathring{\mathbf{A}}^{j}(t,\bmath{x})\partial_{j}\breve{\mathbf{u}}=\mathring{\mathbf{B}}(t,\bmath{x})\breve{\mathbf{u}}.
\label{GeneralLinearisedSystem}
\end{equation}
If one defines $\breve{\mathbf{u}}$ in the same way as in equation \eqref{u} one
obtains the linearised matrix $\mathring{\mathbf{B}}(t,\bmath{x})$ given by
\begin{equation*}
 \begin{pmatrix}
 	\mathbf{B}^{(1)}_{2\times2} & \mathbf{B}^{(2)}_{2\times3} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 	\mathbf{B}^{(3)}_{3\times2} & \mathbf{B}^{(4)}_{3\times3} & 0 & -\frac{k}{4a\omega}\mathbf{B}^{(5)}_{3\times3} & 0 & 0 & 0 & 0 & 0\\
   	0 & 0 & -4H\mathbf{I}_{3\times3} & \frac{k}{4a\omega}\mathbf{B}^{(6)}_{3\times3} & 2\mathbf{I}_{3\times3} & 0 & 0 & 0 & 0\\
 	0 & 0 & -\frac{3k}{4a\omega}\mathbf{B}^{(7)}_{3\times3} &  \left(\frac{2}{\mathring{\psi}}\mathring{\mathcal{V}}^{\prime}+2H\right)\mathbf{I}_{3\times3} & 0 & 0 & \frac{dH}{dt}\mathbf{I}_{3\times3} & 0 & 0\\
        0 & 0 & -\frac{1}{2}\mathring{\psi}^{2}\mathbf{I}_{3\times3} & 0 &-H\mathbf{I}_{3\times3} & \frac{k}{8a\omega}\mathbf{B}^{(8)}_{3\times3} & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 & \frac{k}{8a\omega}\mathbf{B}^{(9)}_{3\times3} &-H\mathbf{I}_{3\times3} & 0 & 0 & 0 \\
	0 & 0 & 0 & \mathbf{I}_{3\times3} & 0 & 0 & -H\mathbf{I}_{3\times3} & 0 & 0 \\
        0 & -\frac{\omega}{3a}\mathbf{B}^{(10)}_{9\times3} & -\frac{\omega}{a}\mathbf{B}^{(11)}_{9\times3} & 0 & 0 & 0 & 0 & -H\mathbf{I}_{9\times9} & 0 \\
	0 & -\frac{k}{12a\omega}\mathbf{B}^{(12)}_{9\times3} & \frac{k}{4a\omega}\mathbf{B}^{(13)}_{9\times3} & -H\mathbf{B}^{(14)}_{9\times3} & 0 & \mathbf{B}^{(15)}_{9\times3} & 0 & 0 & -H\mathbf{I}_{9\times9}
\end{pmatrix}
\end{equation*}
where
\begin{equation*}
 \mathbf{B}^{(1)}_{2\times2}=\begin{pmatrix}
                     0 & 1 \\
		     -\mathring{\mathcal{V}}^{\prime\prime} & -3H
                    \end{pmatrix},\quad \mathbf{B}^{(2)}_{2\times3}=\begin{pmatrix}
                     0 & 0 & 0 \\
		     -\mathring{\psi} & -\mathring{\psi} & -\mathring{\psi}
                     \end{pmatrix},\quad
\mathbf{B}^{(3)}_{3\times2}=\begin{pmatrix}
                     \frac{1}{3}\mathring{\mathcal{V}}^{\prime} & -\frac{2}{3}\mathring{\psi} \\
		     \frac{1}{3}\mathring{\mathcal{V}}^{\prime} & -\frac{2}{3}\mathring{\psi} \\
		     \frac{1}{3}\mathring{\mathcal{V}}^{\prime} & -\frac{2}{3}\mathring{\psi}
                     \end{pmatrix},
\end{equation*}

\begin{equation*}
\mathbf{B}^{(4)}_{3\times3}=\begin{pmatrix}
                     -\frac{2}{3}H & -\frac{2}{3}H & -\frac{2}{3}H \\
		     -\frac{2}{3}H & -\frac{2}{3}H & -\frac{2}{3}H \\
		     -\frac{2}{3}H & -\frac{2}{3}H & -\frac{2}{3}H \\
                    \end{pmatrix},\quad
\mathbf{B}^{(5)}_{3\times3}=\begin{pmatrix}
							      0 & x^{j}\delta^{2}_{j} & x^{j}\delta^{3}_{j} \\
							      x^{j}\delta^{1}_{j} & 0 & x^{j}\delta^{3}_{j} \\
							      x^{j}\delta^{1}_{j} & x^{j}\delta^{2}_{j} & 0
	                                                   \end{pmatrix},\quad
\mathbf{B}^{(6)}_{3\times3}=

						\begin{pmatrix}
						      x_{j}\delta^{j}_{2} & x_{j}\delta^{j}_{1} & 0 \\
						      x_{j}\delta^{j}_{3} & 0 & x_{j}\delta^{j}_{1} \\
						      0 & x_{j}\delta^{j}_{3} & x_{j}\delta^{j}_{2}
			                        \end{pmatrix},
\end{equation*}

\begin{equation*}
			    \begin{pmatrix}
				x^{j}\delta^{2}_{j} & x^{j}\delta^{3}_{j} & 0 \\
                                x^{j}\delta^{1}_{j} & 0 & x^{j}\delta^{3}_{j} \\
				0 & x^{j}\delta^{1}_{j} & x^{j}\delta^{2}_{j}
			    \end{pmatrix},\quad
\mathbf{B}^{(8)}_{3\times3}=
					       \begin{pmatrix}
						  0 & -x^{j}\delta^{1}_{j} & x^{j}\delta^{2}_{j} \\
						  x^{j}\delta^{1}_{j} & 0 & -x^{j}\delta^{3}_{j} \\
						 -x^{j}\delta^{2}_{j} & x^{j}\delta^{3}_{j} & 0
			                       \end{pmatrix},\quad
\mathbf{B}^{(9)}_{3\times3}=
			    \begin{pmatrix}
				0 & x^{j}\delta^{1}_{j} & -x^{j}\delta^{2}_{j} \\
				-x^{j}\delta^{1}_{j} & 0 & x^{j}\delta^{3}_{j} \\
				x^{j}\delta^{2}_{j} & -x^{j}\delta^{3}_{j} & 0
			    \end{pmatrix},
\end{equation*}

\begin{equation*}
\mathbf{B}^{(10)}_{9\times3}=\begin{pmatrix}
                     1 & 1 & 1 \\
		     0 & 0 & 0 \\
		     0 & 0 & 0 \\
		     0 & 0 & 0 \\
		     1 & 1 & 1 \\
		     0 & 0 & 0 \\
		     0 & 0 & 0 \\
                     0 & 0 & 0 \\
                     1 & 1 & 1 \\
                    \end{pmatrix}, \quad
\mathbf{B}^{(11)}_{9\times3}=\begin{pmatrix}
                     0 & 0 & 0 \\
		     1 & 0 & 0 \\
		     0 & 1 & 0 \\
		     1 & 0 & 0 \\
		     0 & 0 & 0 \\
		     0 & 0 & 1 \\
		     0 & 1 & 0 \\
		     0 & 0 & 1 \\
		     0 & 0 & 0
                    \end{pmatrix}, \quad
\mathbf{B}^{(12)}_{9\times3}=\begin{pmatrix}
                     \delta^{1}_{j}x^{j} & \delta^{1}_{j}x^{j} & \delta^{1}_{j}x^{j} \\
                     \delta^{1}_{j}x^{j} & \delta^{1}_{j}x^{j} & \delta^{1}_{j}x^{j} \\
		     0 & 0 & 0 \\
		     \delta^{2}_{j}x^{j} & \delta^{2}_{j}x^{j} & \delta^{2}_{j}x^{j} \\
                     \delta^{2}_{j}x^{j} & \delta^{2}_{j}x^{j} & \delta^{2}_{j}x^{j} \\
		     0 & 0 & 0 \\
                     \delta^{3}_{j}x^{j} & \delta^{3}_{j}x^{j} & \delta^{3}_{j}x^{j} \\
		     \delta^{3}_{j}x^{j} & \delta^{3}_{j}x^{j} & \delta^{3}_{j}x^{j} \\
		     0 & 0 & 0
                    \end{pmatrix},
\end{equation*}

\begin{equation*}
 \mathbf{B}^{(13)}_{9\times3}=\begin{pmatrix}
                     \delta_{2}^{j}x_{j} & 0 & 0 \\
		     0 & \delta_{3}^{j}x_{j} & 0 \\
		     0 & \delta_{2}^{j}x_{j} & -\delta^{1}_{j}x^{j} \\
		     \delta_{1}^{j}x_{j} & 0 & 0 \\
		     0 & 0 & \delta_{3}^{j}x_{j} \\
		     \delta_{3}^{j}x_{j} & -\delta_{j}^{2}x^{j} & 0 \\
		     0 & \delta_{1}^{j}x_{j} & 0  \\
		     0 & 0 & \delta^{j}_{2}x_{j} \\
		     -\delta^{3}_{j}x^{j} & 0 & \delta^{j}_{1}x_{j}
                    \end{pmatrix},\quad
\mathbf{B}^{(14)}_{9\times3}=\begin{pmatrix}
                     1 & 0 & 0 \\
		     1 & 0 & 0 \\
		     0 & 0 & 0 \\
                     0 & 1 & 0 \\
		     0 & 1 & 0 \\
		     0 & 0 & 0 \\
		     0 & 0 & 1 \\
		     0 & 0 & 1 \\
                     0 & 0 & 0
		    \end{pmatrix},\quad
\mathbf{B}^{(15)}_{9\times3}=\begin{pmatrix}
                     0 & 0 & 1 \\
		     0 & 0 & -1 \\
		     0 & 0 & 0 \\
		     0 & -1 & 0 \\
		     0 & 1 & 0 \\
		     0 & 0 & 0 \\
                     1 & 0 & 0 \\
		     -1 & 0 & 0 \\
		     0 & 0 & 0
                    \end{pmatrix}.
\end{equation*}

This matrix $\mathring{\mathbf{B}}$ is made of blocks around the
diagonal and is close to being of lower triangular form. In fact, it
could easily be cast into that form
by `removing' the upper matrix $\mathring{\mathbf{B}}^{(5)}$. In the
more general case where $k\neq0$ one has
three main blocks.  A first block consists of the matrices $\mathring{\mathbf{B}}^{(1)}, \mathring{\mathbf{B}}^{(2)}, \mathring{\mathbf{B}}^{(3)}, \mathring{\mathbf{B}}^{(4)}$,
related to the perturbation variables $(\breve{\phi},\,\,
\breve{\psi},\,\, \breve{\chi}_{11},\,\, \breve{\chi}_{22},\,\,
\breve{\chi}_{33})$. A second block is related to the variables
$((\breve{\chi}^{ST})_{bd},\,\,\breve{a}_{c},\,\,\breve{E}_{bd},\,\,\breve{B}_{bd},\,\,\breve{e}^{\prime}_{b}{}^0)$
and a third one for the perturbation of the spatial connection and the remaining spatial
frame coefficients. This last block does not play an important role in the stability results.

\subsection{Results on the nonlinear stability of symmetric
  hyperbolic systems}

In this subsection we discuss the results of the theory of symmetric
hyperbolic systems required to discuss the nonlinear
stability of perturbations of a reference solution to \emph{the
Einstein-scalar field system with flat spatial sections} ---i.e. $k=0$. The analysis
carried out in the previous sections shows that the Cauchy problem for
these nonlinear perturbations takes the form
\begin{equation}
 \begin{aligned}
 &\left(\mathring{\mathbf{A}}^{0}-\epsilon \breve{{\bf A}}^{0}(t,\breve{{\bf u}},\epsilon)\right)\partial_{t}\breve{{\bf u}}-\left(\mathring{\mathbf{A}}^{j}(t)+\epsilon \breve{{\bf A}}^{j}(t,\breve{{\bf u}},\epsilon)\right)\partial_{j}\breve{{\bf u}}=\left(\mathring{{\bf B}}(t)+\epsilon \breve{{\bf B}}(t,\breve{{\bf u}},\epsilon)\right)\breve{{\bf u}}, \\
 &\breve{{\bf u}}(\bmath{x},0)=\breve{\mathbf{u}}_{0}(\bmath{x}),
\end{aligned}
\label{IVP_nonlinear_pert}
\end{equation}
where $\mathring{\mathbf{A}}^{0}$ is a constant diagonal matrix with
positive entries, while $\mathring{\mathbf{A}}^{j}$ and
$\mathring{\mathbf{B}}$ are matrices whose entries are smooth
functions of the time coordinate $t$.  The matrices
$\mathring{\mathbf{A}}^{j}$, $j=1,2,3$, are
symmetric. Furthermore, $\breve{{\bf A}}^{0}$, $\breve{{\bf A}}^{j}$ and $\breve{{\bf B}}$
vanish at least linearly with $\breve{\mathbf{u}}$. If $k=0$, the linearised system
\eqref{GeneralLinearisedSystem} reduces to
\[
\mathring{\mathbf{A}}^{0}\partial_{t}\breve{\mathbf{u}}-\mathring{\mathbf{A}}^{j}(t)\partial_{j}\breve{\mathbf{u}}=\mathring{\mathbf{B}}(t)\breve{\mathbf{u}}.
\]

\medskip
Our analysis will require the following definition:

\begin{definition}[Relaxed stability eigenvalue condition]
\label{Definition:RSEC}
A system of the form given by \eqref{IVP_nonlinear_pert} is said to
satisfy the relaxed stability eigenvalue condition if the following
conditions hold:
\begin{itemize}
\item[(i)] There is a constant $\delta_{0}$ with
$0<\delta_{0}\leq\delta(t)$ such that, the eigenvalues
$\lambda(t,\bmath{\xi})$ of the symbol
\[
\mbox{\emph{i}}\xi_{j}\mathring{\mathbf{A}}^{j}(t)+\mathring{\mathbf{B}}(t),
\quad \mbox{\emph{i}}=\sqrt{-1}.
\]
 satisfy
\begin{equation*}
 \text{\emph{Re}}\left(\lambda(t,\bmath{\xi})\right)\leq-\delta(t)
\end{equation*}
for all n-tuples of integers $\bmath{\xi}=\left(\xi_{1},...,
  \xi_{n}\right)\neq0$.
\item[(ii)] There is a constant $\delta_{0}$ with $0<\delta_{0}\leq\delta(t)$ such that, for $\bmath{\xi}=\bmath{0}$,
the eigenvalues $\lambda(t,\bmath{0})$ of $\mathring{{\bf B}}(t)$, satisfy either
\begin{equation*}
 \text{Re}\left(\lambda(t,\bmath{0})\right)\leq-\delta(t)\quad\text{or}\quad\lambda(t,\bmath{0})=0.
\end{equation*}
If the multiplicity of the vanishing eigenvalue is $m$, then there are $m$ linearly independent eigenvectors.
\item[(iii)] The Kernel of $\breve{{\bf B}}$ contains the Kernel of $\mathring{{\bf B}}$.
\end{itemize}
\end{definition}

\medskip
In order to state the result that allows us to ascertain the nonlinear stability for
the solutions of the Cauchy problem \eqref{IVP_nonlinear_pert}, we
introduce some notation. Let $H^{k}\left(\mathbb{T}^{n};\mathbb{R}^{s}\right)$ be the
space of all summable functions ${\bf
u}(\cdot,t):\mathbb{T}^{n}\rightarrow\mathbb{R}^{s}$ such that for
each multi-index $|\alpha|\leq k$, $D^{\alpha}{\bf u}(\bmath{x},t)$
exists in the weak sense and belongs to
$L^{2}\left(\mathbb{T}^{n}\right)$. The norm in $H^{k}\left(\mathbb{T}^{n};\mathbb{R}^{s}\right)$ is defined by
\begin{equation*}
\|{\bf u}(t)\|_{H^{k}(\mathbb{T}^{n})} \equiv \left(\sum_{|\alpha|=0}^{k}\int_{\mathbb{T}^{n}}|D^{\alpha}{\bf u}|^{2}\mbox{d}\bmath{x}\right)^{1/2},
\end{equation*}
 where
\[
D^{\alpha}=\frac{\partial^{|\alpha|}}{\partial
x^{\,\,\alpha_{1}}_{1}...\partial x^{\,\,\alpha_{n}}_{n}}.
\]
Given a perturbation $\breve{\mathbf{u}}(\bmath{x},t)$, we denote by
$\hat{\breve{\mathbf{u}}}(\bmath{\xi},t)$ the Fourier coefficient
corresponding to the frequency vector $\bmath{\xi}$ ---recall that
it is assumed that the topology of the spatial sections is taken to
be that of $\mathbb{T}^3$. In particular,
$\hat{\breve{\mathbf{u}}}(\bmath{0},t)$ denotes the coefficient corresponding
to $\bmath{\xi}=\bmath{0}$.

\medskip
The following theorem shows the exponential decay or convergence to
constant values of a solution to the nonlinear system
\eqref{IVP_nonlinear_pert}. Its proof follows closely the methods of
\cite{1997funct.an..3003K,JMP2} and will be omitted for the sake of
brevity.

\begin{thm}
\label{Theorem:stabilitytheory}
Let $\mathbf{P}$ denote the projector into the Kernel of
$\mathring{\mathbf{B}}(t)$. Furthermore,  let
$\breve{\mathbf{u}}^{(0)}\equiv \mathbf{P}\hat{\breve{u}}(\bmath{0},t)$
and $\breve{\mathbf{v}} \equiv \breve{\mathbf{u}}-\breve{\mathbf{u}}^{(0)}$.
If the symmetric hyperbolic system of the Cauchy problem \eqref{IVP_nonlinear_pert} satisfies the relaxed stability eigenvalue condition,
then the system is stable in the sense that, for sufficiently small $\epsilon$ there is $\epsilon_{0}\geq\epsilon>0$
such that
\begin{equation*}
\lim_{t\rightarrow\infty}\|\breve{\mathbf{v}}(t)\|_{H^{n+2}(\mathbb{T}^{n};\mathbb{R}^{s})}=0
\end{equation*}
and $\breve{\mathbf{u}}^{(0)}(t)\rightarrow const.$ as time goes to infinity.
\end{thm}

\medskip
\noindent
\textbf{Remark.} The use of Fourier methods in
Theorem \ref{Theorem:stabilitytheory} does not allow to disentagle
which fields (or combination thereof) decay to a constant. However, it
is possible to identify the blocks in $\mathring{\mathbf{B}}(t)$
responsible for this behaviour.

\medskip
In the remainder of the article, it will be analysed under which
conditions the symmetric hyperbolic evolution equations for
the perturbations of a reference solution to the Einstein-scalar field system
satisfy the requirements of the relaxed stability eigenvalue condition
(Definition \ref{Definition:RSEC}) so that Theorem
\ref{Theorem:stabilitytheory} can be used.

\subsection{Satisfying Assumption {\it (ii)} of Definition \ref{Definition:RSEC}}

In the case of a flat background ($k=0$, $\omega=1$), the linearized
matrices $\mathring{\mathbf{A}}$ and $\mathring{\mathbf{B}}$ are
functions of time only. Because of the block structure of the matrix, the
characteristic polynomial of $\mathring{\mathbf{B}}$ is the product of
the characteristic polynomials of the block matrices around the
diagonal. A computation gives the following polynomial:
\begin{eqnarray}
&&\left(\lambda+{H}\right)^{21}\times\left(\lambda^2-\left(2\frac{\mathring{\mathcal{V}}^{\prime}}{{\mathring{\psi}}}+H\right)\lambda-\frac{\mbox{d}H}{\mbox{d}t}-2H^{2}-2H\frac{\mathring{\mathcal{V}}^{\prime}}{{\mathring{\psi}}}\right)^{3} \times\left(\lambda^{2}+5{H}\lambda+\mathring{\psi}^{2}+4{H}^{2}\right)^{3}\nonumber
\\
&&\hspace{2cm}\times\bigg(\lambda^{5}+5{H}\lambda^{4}-\left(2\mathring{\psi}^{2}-6{H}^{2}-\mathring{\mathcal{V}}^{\prime\prime}\right)\lambda^{3}+\left(2{H}\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}\right)\lambda^{2}\bigg). \label{characteristicpolynomial}
\end{eqnarray}
The first term in the expression \eqref{characteristicpolynomial} is
the characteristic polynomial of the block related with the
perturbation variables
$(\breve{B}_{bd},\,\,\breve{e}^\prime_{b}{}^j,\,\,\gamma{}^{\prime
c}_{\,\,\,ab})$.  The second term is related to the variables
$(\breve{a}_{c},\,\,\breve{e}^{\prime}_{b}{}^0)$, and the third one to
$((\breve{\chi}^{ST})_{bd},\,\,\breve{E}_{bd})$. Finally, the last
term ---the fifth order polynomial--- arises from the block of
unknowns $(\breve{\phi},\,\,\breve{\psi},\,\, \breve{\chi}_{pp})$.

\medskip

In order to obtain conditions from the characteristic polynomial we will make use of results from
the {\it Li\'enard-Chipart theorem}. The latter gives necessary and sufficient
conditions for a polynomial with real coefficients to have roots
with negative real part ---see e.g. \cite{AnalPropPol2002}. This type
of polynomials are called \emph{Hurwitz polynomials}.

\begin{thm}[Li\'enard-Chipart]
Let
\[
f(z)= a_{0}z^{n}+a_{1}z^{n-1}+a_{2}z^{n-2}+\cdots+a_{n}, \quad a_{0}>0,
\]
be a polynomial with real coefficients. Then the following statements are equivalent:
\begin{itemize}
 \item[(i)] the polynomial is a Hurwitz polynomial;
 \item[(ii)] The coefficients of f are positive and $\delta_{2}>0,\,\delta_{4}>0,...,\delta_{n},\quad n$ even;
 \item[(iii)] The coefficients of f are positive and $\delta_{1}>0,\,\delta_{3}>0,...,\,\delta_{n},\quad n$ odd,
\end{itemize}
where the Hurwitz determinants are defined by
\begin{equation*}
 \delta_{0}\equiv1,\quad\delta_{l}\equiv\mbox{\emph{det}}\begin{pmatrix}
                                    a_{1}&a_{3}&a_{5}&\cdots&a_{2l-1} \\
				        1&a_{2}&a_{4}&\cdots&a_{2l-2} \\
					0&a_{1}&a_{3}&\cdots&a_{2l-3} \\
				   \vdots&\vdots&\vdots&\quad&\vdots  \\
					0&0&0&\cdots&a_{l}
                                   \end{pmatrix}\quad \quad l=1,\dots,n.
\end{equation*}
\end{thm}
It is easy to see that the first term in \eqref{characteristicpolynomial} requires
\begin{equation}
{H}(t)>0
\label{condition_H}
\end{equation}
which is simply the condition for an (ever) expanding background.
From the second factor in the characteristic polynomial \eqref{characteristicpolynomial},
one obtains the following conditions on $\mathring{\mathcal{V}}^{\prime}$:
\begin{equation}
\label{cond_1stderpotential}
\frac{\mathring{\mathcal{V}}^{\prime}}{\mathring{\psi}}<-\frac{H}{2},   \quad \quad 2H\frac{\mathring{\mathcal{V}}^{\prime}}{{\mathring{\psi}}}<-\frac{\mbox{d}H}{\mbox{d}t}-2H^{2} .
\end{equation}
Notice that by construction we have that $\psi<0$ ---see equation
\eqref{defipsi}. The first condition in \eqref{cond_1stderpotential}
implies that the first derivative of the scalar field potential is
positive. The second condition is related to
the evolution equation for the Hubble function \eqref{evol_background}.
Moreover, by virtue of condition \eqref{AsymptoticCondition} the ratio
$\mathring{\mathcal{V}}^{\prime}/\mathring{\psi}$ tends to a constant
which is bigger than the asymptotic value of $H/2$.

\medskip
The fifth order polynomial in \eqref{characteristicpolynomial} can be rewritten as
\begin{equation*}
 \lambda^2\left(\lambda^{3}+5{H}\lambda^{2}-\left(2\mathring{\psi}^{2}-6{H}^{2}-
 \mathring{\mathcal{V}}^{\prime\prime}\right)\lambda+\left(2{H}\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}\right)\right)
\end{equation*}
from where one sees that the matrix $\mathring{{\bf B}}(t)$ has a
double vanishing eigenvalue. A direct application of the
Li\'enard-Chipart Theorem to the remaining third
order polynomial renders the following two conditions:
\begin{equation}
\label{cond_2ndderpotential}
 \mathring{\mathcal{V}}^{\prime\prime}>-\mathring{\psi}\frac{\mathring{\mathcal{V}}^{\prime}}{2H}, \quad
                  \quad \mathring{\mathcal{V}}^{\prime\prime}-\mathring{\psi}\frac{\mathring{\mathcal{V}}^{\prime}}{3H}>10\left(-H^{2}+\frac{\mathring{\psi}^{2}}{3}\right).
\end{equation}

The first condition relates the second derivative of the potential to its first derivative, and implies that
$\mathring{\mathcal{V}}^{\prime\prime}$ has to be positive and
bounded. The second condition is sharper and also relates to the evolution of the Hubble function
\eqref{evol_background}.
Taking the limit $\psi\rightarrow0$, one recovers the necessary asymptotic conditions for
stability found in \cite{2008Ringstrom}.

\medskip
 \subsubsection{Nonflat backgrounds}
For the sake of conciseness, the case of  a reference solution with nonflat spatial sections will
not be  treated here. Nevertheless,  we point out some features of the
required  analysis. For this class of solutions the characteristic
polynomial of the matrix $\mathring{\mathbf{B}}$ is given by
\begin{eqnarray}
 && f(\lambda)\times\bigg(\lambda^{5}+5H\lambda^{4}-\left(2\mathring{\psi}^{2}-6H^{2}-\mathring{\mathcal{V}}^{\prime\prime}\right)\lambda^{3}+\left(2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}\right)\lambda^{2}\bigg)\times\left(\lambda+H\right)^{18},
\label{OpenBack_CharPol}
\end{eqnarray}
where $f(\lambda)$ is a polynomial of degree 15 in $\lambda$.
The coefficients of this polynomial are, unfortunately, large expressions
without an obvious structure. Accordingly, an
analysis of the roots by means of the Lienard-Chipart Theorem,
although in principle possible, becomes extremely cumbersome.
The study of perturbations of cosmological models with nonflat spatial
sections  also requires to take into account extra terms coming from
partial differentiation  when constructing the Sobolev norms. This is because the matrices $\mathring{\mathbf{A}}^{j}$
and $\mathring{\mathbf{B}}$ are now functions of the coordinates
$x^{i}$.

\medskip
Finally, it is mentioned that a discussion of the perturbations of
reference solutions with nonflat spatial sections requires suitable
adaptations of the relaxed stability eigenvalue condition (Definition
\ref{Definition:RSEC}) and of Theorem \ref{Theorem:stabilitytheory}.

\subsection{Satisfying Assumption {\it (i)} of Definition \ref{Definition:RSEC}}

The purpose of the present section is to show that Assumption (i) of
the  {\it stability eigenvalue condition} follows from Assumption
(ii). This turns out to be a consequence of a theorem due to Ortiz
in \cite{JMP2} ---see also \cite{JMP1}.

\begin{thm}
\label{Theorem:Ortiz}
Let $\mathbf{S}$ be a nonsingular matrix which takes the $s\times s$ matrix
$\mathring{\mathbf{B}}$ into block diagonal form.
If for all
unit frequency vectors $\bmath{\xi}$, the linear map
\begin{equation*}
 \mbox{\emph{i}}\mathbf{S}^{-1}\mathring{\mathbf{A}}^{0}\mathbf{S}+\mbox{\emph{i}}\mathring{\mathbf{A}}(t,{\bmath{\xi}}):\,\, Ker(\mathring{\mathbf{B}}(t))\rightarrow \mathbb{C}^{s}
\end{equation*}
with
\[
\mathring{\mathbf{A}}(t,{\bmath{\xi}}) \equiv -\mbox{\emph{i}}\mathbf{S}^{-1}\xi_{j}\mathring{\mathbf{A}}^{j}\mathbf{S},
\]
is injective,  then there are positive constants $\xi_{0}$,
$\delta_{0}$ and $\mu_{0}=\delta_{0}/c$,
such that all eigenvalues $\lambda(t,\bmath{\xi})$ of the symbol $\mbox{\emph{i}}\xi_{j}\mathring{\mathbf{A}}^{j}(t)+\mathring{\mathbf{B}}(t)$ satisfy
\begin{equation*}
 Re{\lambda(t,\bmath{\xi})}\leq\begin{cases}
				  -\mu_{0}|\bmath{\xi}|^{2}\,\,\,\text{or}\,\, -\delta_{0},\,\,\,&\text{if}\,\,\, |\bmath{\xi}|\leq\xi_{0} \\
				  -\delta_{0},       & \text{if } |\bmath{\xi}|\geq\xi_{0}
				\end{cases}.
\end{equation*}
\end{thm}

\medskip
\noindent
\textbf{Remark.} This theorem states that if the determinant of the
linear map is different from zero, then the eigenvalues of the
linearised symbol have negative real part for nonzero Fourier
frequency, decaying no faster than $|\bmath{\xi}|^{2}$ when
$|\bmath{\xi}|\rightarrow0$.

\medskip
In order to apply this theorem, we must first put $\mathring{\mathbf{B}}$ into block diagonal
form. To this end we need to find a transformation, $\mathbf{S}$, such that
\begin{equation}
\mathbf{S}^{-1}\mathring{\mathbf{B}}\mathbf{S}
                   =\begin{pmatrix}
                     \mathbf{0} & \mathbf{\mathbf{0}}  \\
		     \mathbf{0} & \mathring{\mathbf{B}}
                    \end{pmatrix}
\label{Bblockform}.
\end{equation}
To achieve this, it is sufficient to consider the block matrix
corresponding to the degree 5 polynomial which
contains the two zero eigenvalues. Let us define the $5\times5$ matrix $\mathring{\mathbf{M}}$ by
\begin{equation*}
\mathring{\mathbf{M}}=\begin{pmatrix}
                            \mathbf{B}^{(1)}_{2\times2} & \mathbf{B}^{(2)}_{2\times3} \\
 	                    \mathbf{B}^{(3)}_{3\times2} & \mathbf{B}^{(4)}_{3\times2}
                       \end{pmatrix}
                     =\begin{pmatrix}
                            0 & 1 & 0 & 0 & 0 \\
		            -\mathring{\mathcal{V}}^{\prime\prime} & -3H & -\mathring{\psi} & -\mathring{\psi} & -\mathring{\psi} \\
			    \frac{1}{3}\mathring{\mathcal{V}}^{\prime} & -\frac{2}{3}\mathring{\psi} & -\frac{2}{3}H & -\frac{2}{3}H & -\frac{2}{3}H \\
			    \frac{1}{3}\mathring{\mathcal{V}}^{\prime} & -\frac{2}{3}\mathring{\psi} & -\frac{2}{3}H & -\frac{2}{3}H & -\frac{2}{3}H \\
			    \frac{1}{3}\mathring{\mathcal{V}}^{\prime} & -\frac{2}{3}\mathring{\psi} & -\frac{2}{3}H & -\frac{2}{3}H & -\frac{2}{3}H
                      \end{pmatrix}.
\end{equation*} \\
Once $\mathring{\mathbf{M}}$ has put into block diagonal form, it
turns out that the condition \eqref{Bblockform} holds.
In order to find a suitable transformation, $\mathring{\mathbf{S}}$,
we first make a change of basis taking $\mathring{\mathbf{M}}$
into block triangular form. For this, we compute the eigenvectors
corresponding to the two vanishing eigenvalues.  From  the condition
\[
\mathring{\mathbf{M}}\mathbf{v}=\mathbf{0},
\]
it follows that the eigenvectors can be taken to be
\begin{equation*}
 \mathbf{v}_{1}=\begin{pmatrix}
             -\frac{3\mathring{\psi}}{\mathring{\mathcal{V}}^{\prime\prime}} \\
               0   \\
               1   \\
               1  \\
               1
            \end{pmatrix},\quad
\mathbf{v}_{2}=\begin{pmatrix}
	       \frac{6H}{\mathring{\mathcal{V}}^{\prime}} \\
               0   \\
               1  \\
               1  \\
               1
	    \end{pmatrix}.
\end{equation*}
With the help of these two eigenvectors we can easily transform
$\mathring{\mathbf{M}}$ into block triangular form by completing the
basis with three further linearly independent eigenvectors. For
instance, one can consider
\begin{equation*}
 \mathbf{v}_{3}=\begin{pmatrix}
               0 \\
               1   \\
               0   \\
               0  \\
               0
            \end{pmatrix}\quad,\quad
\mathbf{v}_{4}=\begin{pmatrix}
	       0 \\
               0   \\
               1  \\
               0  \\
               1
	    \end{pmatrix}\quad,\quad
\mathbf{v}_{5}=\begin{pmatrix}
	       0 \\
               0   \\
               0  \\
               1  \\
               1
	    \end{pmatrix}.
\end{equation*}
Using these basis vectors as the columns of a transformation matrix
one arrives to
\begin{equation*}
 \mathbf{P}\equiv \begin{pmatrix}
             -\frac{3\mathring{\psi}}{\mathring{\mathcal{V}}^{\prime\prime}} & \frac{6H}{\mathring{\mathcal{V}}^{\prime}} & 0 & 0 & 0 \\
              0 & 0 & 1 & 0 & 0 \\
              1 & 1 & 0 & 1 & 0 \\
              1 & 1 & 0 & 0 & 1 \\
              1 & 1 & 0 & 1 & 1
            \end{pmatrix}
\end{equation*}
with inverse given by
\begin{equation*}
\mathbf{P}^{-1}=\begin{pmatrix}
                 -\frac{1}{3}\frac{\mathring{\mathcal{V}}^{\prime\prime}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & 0 & \frac{2H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & \frac{2H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & -\frac{2H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                  \frac{1}{3}\frac{\mathring{\mathcal{V}}^{\prime\prime}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & 0 & \frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}  & \frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}  &  -\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                  0 & 1 & 0 & 0 & 0 \\
		  0 & 0 & 0 & -1 & 1 \\
		  0 & 0 & -1 & 0 & 1
                \end{pmatrix}.
\end{equation*}
As a result of these computations one finds that
\begin{equation*}
\begin{aligned}
 \mathbf{P}^{-1}\mathring{\mathbf{M}}\mathbf{P}&=\begin{pmatrix}
					  \mathbf{L}_{3\times3} & \mathbf{H}_{3\times2} \\
					  \mathbf{0}_{2\times3} & \mathbf{0}_{2\times2}
				      \end{pmatrix}
\end{aligned}
\end{equation*}
with
\begin{equation*}
\mathbf{L}_{3\times3}\equiv\begin{pmatrix}
                         -2H & 0 & -\frac{\mathring{\mathcal{V}}^{\prime\prime}}{3}\frac{4H\mathring{\psi}+\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                         -\mathring{\psi}\frac{\mathring{\mathcal{V}}^{\prime}}{\mathring{\mathcal{V}}^{\prime\prime}} & 0 & -\frac{\mathring{\mathcal{V}}^{\prime}}{3}\frac{2\mathring{\psi}^{2}-\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                         0 & -\frac{3}{\mathring{\mathcal{V}}^{\prime}}\left(2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}\right) & -3H
		      \end{pmatrix},\quad
\mathbf{H}_{3\times2}\equiv \begin{pmatrix}
                         -\frac{8}{3}\frac{H^{2}\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & -\frac{8}{3}\frac{H^{2}\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                         -\frac{4}{3}\frac{H\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & -\frac{4}{3}\frac{H\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                         -2\mathring{\psi} & -2\mathring{\psi}
                      \end{pmatrix}.
\end{equation*}

It can be readily verified that the matrix $\mathbf{L}$ is negative
definite as long as the conditions \eqref{cond_2ndderpotential} of the
previous section hold. Now, in order to proceed further we need another
transformation which converts the block upper triangular form into
block diagonal.  Choosing
\begin{equation*}
 \mathbf{R}\equiv \begin{pmatrix}
             \mathbf{I}_{3\times3} & \mathbf{Q}_{3\times2} \\
             \mathbf{0}_{2\times3} & \mathbf{I}_{2\times2}
            \end{pmatrix},\quad
 \mathbf{R}^{-1}=\begin{pmatrix}
             \mathbf{I}_{3\times3} & -\mathbf{Q}_{3\times2} \\
             \mathbf{0}_{2\times3} & \mathbf{I}_{2\times2}
            \end{pmatrix}
\end{equation*}
one finds that
\begin{equation*}
 \mathbf{R}^{-1}\left(\mathbf{P}^{-1}\mathring{\mathbf{M}}\mathbf{P}\right)\mathbf{R}=\begin{pmatrix}
                                      \mathbf{L} & \mathbf{L}\mathbf{Q}+\mathbf{H} \\
                                      \mathbf{0} & \mathbf{0}
                                     \end{pmatrix}.
\end{equation*}
Thus, in order to obtain the required block-diagonal form we just need to solve
\begin{equation*}
 \mathbf{L}\mathbf{Q}=-\mathbf{H}.
\end{equation*}
A computation shows that
\begin{equation*}
 \mathbf{Q}=\begin{pmatrix}
             -\frac{4}{3}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & -\frac{4}{3}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
             -\frac{2}{3}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & -\frac{2}{3}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2\mathring{H}\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
              0 & 0
            \end{pmatrix},
\end{equation*}
so that the whole transformation matrix $\mathbf{S}$ is given by
\begin{equation*}
 \mathbf{S} \equiv \begin{pmatrix}
              (\mathbf{P}\mathbf{R})_{5\times5} & \mathbf{0}_{5\times30} \\
              \mathbf{0}_{30\times5} & \mathbf{I}_{30\times30}
             \end{pmatrix}
\end{equation*}
where
\begin{equation*}
 \mathbf{P}\mathbf{R}=\begin{pmatrix}
				  -\frac{3\mathring{\psi}}{\mathring{\mathcal{V}}^{\prime\prime}} & \frac{6H}{\mathring{\mathcal{V}}^{\prime}} & 0 & 0 & 0 \\
				    0 & 0 & 1 & 0 & 0 \\
				    1 & 1 & 0 & \frac{1}{3} & -\frac{2}{3} \\
				    1 & 1 & 0 & -\frac{2}{3} & \frac{1}{3} \\
				    1 & 1 & 0 & \frac{1}{3} & \frac{1}{3}
		      \end{pmatrix}
\end{equation*}
with inverse given by
\begin{equation*}
 (\mathbf{PR})^{-1}=\begin{pmatrix}
                      -\frac{1}{3}\frac{\mathring{\mathcal{V}}^{\prime}\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & 0 & \frac{2}{3}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & \frac{2}{3}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & \frac{2}{3}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                      \frac{1}{3}\frac{\mathring{\mathcal{V}}^{\prime}\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & 0 & \frac{1}{3}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & \frac{1}{3}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} & \frac{1}{3}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}} \\
                      0 & 1 & 0 & 0 & 0 \\
                      0 & 0 & 0 & -1 & 1 \\
                      0 & 0 & -1 & 0 & 1
                     \end{pmatrix}.
\end{equation*}
Summarising: we have constructed the matrix $\mathbf{S}$ which
transforms the linearised matrix $\mathring{\mathbf{B}}$ into block
diagonal form.  In the following we  use this matrix to show that
the linear map arising in the formulation of Theorem \ref{Theorem:Ortiz} is injective.

\medskip
Applying the transformation $\mathbf{S}$, to the linearised matrices
of equation \eqref{A_matrices} ---see also equation
\eqref{linearized_syst}--- we get that
\begin{equation*}
\mathbf{S}^{-1}\xi_{j}\mathring{\mathbf{A}}^{j}\mathbf{S}=\begin{pmatrix}
						0 & 0 & \mathbf{A}^{(1)}_{5\times3} & 0 & 0 & 0 \\
						0 & 0 & \mathbf{A}^{(2)}_{3\times3} & 0  & 0 & 0 \\
						\mathbf{A}^{(3)}_{3\times5} & (\mathbf{A}^{(2)})^{T}_{3\times3} & 0 & 0 & 0 & 0 \\
                                                0 & 0 & 0 & 0 & \mathbf{A}^{(4)}_{3\times3} & 0 \\
                                                0 & 0 & 0 & (\mathbf{A}^{(4)})^{T}_{3\times3} & 0 & 0\\
                                                 0 & 0 & 0 & 0 & 0 & 0
						\end{pmatrix}
\end{equation*}
where
\begin{equation*}
 \mathbf{A}^{(1)}=\begin{pmatrix}
                           \frac{2\omega}{3a}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}\xi_{1} & \frac{2\omega}{3a}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}\xi_{2} & \frac{2\omega}{3a}\frac{H\mathring{\mathcal{V}}^{\prime\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}\xi_{3} \\
                           \frac{\omega}{3a}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}\xi_{1} & \frac{\omega}{3a}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}\xi_{2} & \frac{\omega}{3a}\frac{\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}{2H\mathring{\mathcal{V}}^{\prime\prime}+\mathring{\psi}\mathring{\mathcal{V}}^{\prime}}\xi_{3} \\
                            0 & 0 & 0 \\
                            0 & -\frac{\omega}{a}\xi_{2} & \frac{\omega}{a}\xi_{3} \\
                            -\frac{\omega}{a}\xi_{1} & 0 & \frac{\omega}{a}\xi_{3}
                          \end{pmatrix},\quad
\mathbf{A}^{(2)}=\begin{pmatrix}
                          \frac{\omega}{a}\xi_{2} & \frac{\omega}{a}\xi_{1} & 0 \\
                          \frac{\omega}{a}\xi_{3} & 0 & \frac{\omega}{a}\xi_{1} \\
                           0 & \frac{\omega}{a}\xi_{3} & \frac{\omega}{a}\xi_{2}
                         \end{pmatrix},
\end{equation*}
\begin{equation*}
 \mathbf{A}^{(3)}=\begin{pmatrix}
                            \frac{\omega}{a}\xi_{1} & \frac{\omega}{a}\xi_{1} & 0 & \frac{\omega}{3a}\xi_{1} & -\frac{2\omega}{3a}\xi_{1} \\
                            \frac{\omega}{a}\xi_{2} & \frac{\omega}{a}\xi_{2} & 0 & -\frac{2\omega}{3a}\xi_{2} &  \frac{\omega}{3a}\xi_{2} \\
                            \frac{\omega}{a}\xi_{3} & \frac{\omega}{a}\xi_{3} & 0 & \frac{\omega}{3a}\xi_{3} & \frac{\omega}{3a}\xi_{3}
                          \end{pmatrix},\quad
\mathbf{A}^{(4)}=\begin{pmatrix}
                          0 & \frac{\omega}{2a}\xi_{1} & -\frac{\omega}{2a}\xi_{2} \\
                          -\frac{\omega}{2a}\xi_{1} & 0 & \frac{\omega}{2a}\xi_{3} \\
                           \frac{\omega}{2a}\xi_{2} & -\frac{\omega}{2a}\xi_{3} & 0
                         \end{pmatrix}.
\end{equation*}
It is also noticed that the diagonality of $\mathring{\mathbf{A}}^{0}$
implies that it is invariant by $\mathbf{S}$.  It follows that the
linear map
$\mbox{i}\mathring{\mathbf{A}}^{0}+\mbox{i}\mathring{\mathbf{A}}(t,{\bmath{\xi}})$
has determinant
\begin{equation*}
\begin{aligned}
\frac{\mbox{i}}{8a^{10}}\left(\omega^{6}|\bmath{\xi}|^{6}+5a^{2}\omega^{4}|\bmath{\xi}|^{4}+8a^{4}\omega^{2}|\bmath{\xi}|^{2}+4a^{6}\right)
\left(\omega^{4}|\bmath{\xi}|^{4}+8a^{2}\omega^{2}|\bmath{\xi}|^{2}+16a^{4}\right),
\end{aligned}
\end{equation*}
where $|\bmath{\xi}|^{2}\equiv
\xi_{1}^{2}+\xi_{2}^{2}+\xi_{3}^{2}$. Therefore, for all $\xi$,
$a$ and $\omega$ the determinant is nonzero, and by Theorem
\ref{Theorem:Ortiz} one finds that Assumption \emph{(i)} of Definition
\ref{Definition:RSEC} holds.

\subsection{Verifying Assumption {\it (iii)} of Definition
  \ref{Definition:RSEC}}

In order to verify Assumption \emph{(iii)} of the {\it relaxed
stability eigenvalue condition}, we just need to consider the block
matrix of the nonlinear matrix, $\breve{\mathbf{B}}$, associated to
the perturbation variables $(\breve{\phi},\,\, \breve{\psi},\,\,
\breve{\chi}_{11},\,\, \breve{\chi}_{22},\,\, \breve{\chi}_{33})$. A
simple inspection of the linearisation process allow us to verify that
this block matrix has zeros in all its entries. Therefore Assumption
\emph{(iii)} follows directly.

\section{The main result}

We summarise the discussion of previous sections in the following theorem:

\begin{thm}
\label{Theorem:Main}
Let $\phi$ be a smooth real scalar field with a selfinteracting potential
$\mathcal{V}(\phi)$ in a flat expanding Friedmann-Robertson-Walker spacetime subject
to the nonlinear evolution equations
\begin{equation*}
 \frac{\mbox{\em d}\mathring{\phi}}{\mbox{\em d}t}=\mathring{\psi}
 \quad \text{and} \quad\frac{\mbox{\em d}\mathring{\psi}}{\mbox{\em
     d}t}=-3H\mathring{\psi}-{\frac{\mbox{\em
       d}\mathring{\mathcal{V}}}{\mbox{\em d}\mathring{\phi}}}. \\
\end{equation*}
Suppose there exists a constant $H_{0}$ such
that $0<H_{0}\leq H(t)$, for all $t\geq0$, and a potential
satisfying conditions
\begin{equation*}
  \frac{\mathring{\mathcal{V}}^{\prime}}{\mathring{\psi}}<-\frac{H}{2} \quad,\quad  2H\frac{\mathring{\mathcal{V}}^{\prime}}{{\mathring{\psi}}}<-\frac{\mbox{\em d}H}{\mbox{\em d}t}-2H^{2}  \quad, \quad  \mathring{\mathcal{V}}^{\prime\prime}>-\mathring{\psi}\frac{\mathring{\mathcal{V}}^{\prime}}{2H} \quad\text{and}
                  \quad \mathring{\mathcal{V}}^{\prime\prime}-\mathring{\psi}\frac{\mathring{\mathcal{V}}^{\prime}}{3H}>10\left(-H^{2}+\frac{\mathring{\psi}^{2}}{3}\right).
\end{equation*}
Then the FRW-nonlinear scalar field solution is stable in the sense that, given any initial
data for small nonlinear perturbations $\breve{\bmath{u}}_{0}$ whose norm
$\|\breve{\bmath{u}}\|_{H^{5}(\mathbb{T}^{3})}$ is finite, there exists $\epsilon_{0}>0$ such that for all
$0<\epsilon\leq\epsilon_{0}$ a solution to the nonlinear
perturbations exists for all times and decays exponentially to zero or
converges to constant values.
\end{thm}

\medskip
\noindent
\textbf{Remark.} Our main result can be understood as a \emph{cosmic no-hair theorem} for a wide class of scalar
field cosmologies. We note that our results can be useful when
considering numerical simulations of these models as
it gives conditions for exponential decay during the evolutions. It
could,  therefore, be used to test the consistency of the numerical
steps.

\subsection{Power law solutions with exponential potential}
As a direct application of Theorem \ref{Theorem:Main} we discuss the
so-called \emph{power-law} solutions to an exponential potential and show that this class of
solutions is covered by our analysis.

\medskip
 Consider the reference solution for a scalar field $\phi$ with exponential potential
\begin{equation*}
 \mathcal{V}(\phi)=\Lambda e^{\lambda\phi},
\end{equation*}
with $\Lambda$, $\lambda$ positive constants, given by the metric
\begin{equation*}
 \mbox{d}s^{2}=-\mbox{d}t^{2}+t^{2p}\delta_{ij}\mbox{d}x^{i}\mbox{d}x^{j}
\end{equation*}
with
\[
 \phi(t)=\frac{\sqrt{2p}}{2}\ln\left(\frac{p(3p-1)}{\Lambda
     t^{2}}\right), \quad \psi(t)=-\frac{\sqrt{2p}}{t}
\]
with $p>1/3$ and $\lambda^{2}p=2$. A direct computation shows that the
conditions in Theorem \ref{Theorem:Main} are satisfied if
\[
p> \tfrac{\sqrt{61}-1}{15} \approx 0.454.
\]
Thus, for this range of values one can conclude that \emph{sufficiently small
nonlinear perturbations to power-law solutions decay
exponentially to zero or converge to constant values}.

\medskip
\noindent
\textbf{Remark 1.} It is well
known that for $p>1$ (i.e. $\lambda<\sqrt{2}$), the exact solution
described in the previous lines has
accelerated expansion --- the so-called \emph{power law inflation}. Our main
result, Theorem \ref{Theorem:Main} shows that one has decay of the nonlinear
perturbations even in the absence of accelerated expansion if $p\in
(\tfrac{\sqrt{61}-1}{15},1]$. This is, to the best of our knowledge, a new result ---the
possibility of having decay in the marginal case $p=1$ has been
discussed in \cite{2009Ringstrom,Rendall2007}.

\medskip
\noindent
\textbf{Remark 2.} Simple numerical simulations of the evolution of
linearised perturbations supports the results of our analysis. A
detailed discussion of these is beyond the purposes of the present work.

\section*{ Acknowledgements}

AA and FM are supported by projects PTDC/MAT/108921/2008 and CERN/FP/116377/2010
and by CMAT, Univ. Minho, through FCT plurianual funding.
AA thanks the Relativity Group at the School of Mathematical Sciences,
Queen Mary, Univ. London, for their warm hospitality while most of
this work was done and FCT for grant SFRH/BD/48658/2008. FM thanks FCT
grant SFRH/BSAB/967/2010 and hospitality from CBPF, Rio de
Janeiro. JAVK was supported by an EPSRC Advanced Research Fellowship
and by a project research grant from the Leverhulme Trust  (F/07
476/AI). JAVK thanks the Centre of Mathematics of the University of
Minho for its hospitality.

config: arxiv