text:
\section{Introduction}
Throughout this paper, all rings are commutative with nonzero identity and if $R$ is a subring of $T$, then the identity element of $R$ and $T$ is same. By an overring of $R$, we mean a subring of the total quotient ring of $R$ containing $R$. By a local ring, we mean a ring with unique maximal ideal. The symbol $\subseteq$ is used for inclusion, while $\subset$ is used for proper inclusion. We use $T(R)$ to denote the total quotient ring of $R$, $R'$  to denote the integral closure of a ring $R$ in $T(R)$, $Nil(R)$ to denote the set of nilpotent elements of $R$, and $Z(R)$ to denote the set of zero-divisors of $R$. If $I$ and $J$ are $R$-submodules of $T(R)$, then $(I : J) = \{x\in T(R) : xJ\subseteq I\}$. Badawi defined the divided prime ideal of $R$ in $\textup{\cite{bd}}$, as a prime ideal which is comparable to every ideal of $R$ and a ring $R$ is said to be divided if each prime ideal of $R$ is divided, see $\textup{\cite{bd}}$. In $\textup{\cite{badawi3}}$, Anderson and Badawi introduced the notion $\mathcal{H}$ to be the set of all rings $R$ such that $Nil(R)$ is a divided prime ideal and named these rings as $\phi$-rings. They used $\mathcal{H}_0$ to denote the subset of $\mathcal{H}$ such that $Nil(R) = Z(R)$. For further study on these rings, see $\textup{\cite{badawi3}}$, $\textup{\cite{badawi4}}$, $\textup{\cite{badawi}}$, $\textup{\cite{badawi8}}$, $\textup{\cite{badawi6}}$, $\textup{\cite{badawi9}}$, $\textup{\cite{badawi10}}$, $\textup{\cite{badawi7}}$, $\textup{\cite{badawi1}}$, $\textup{\cite{badawi11}}$, $\textup{\cite{badawi12}}$.

In this paper, we will focus on the ring extensions in class $\mathcal{H}$ having the same set of prime ideals which was extensively studied by Anderson and Dobbs for integral domains, see $\textup{\cite{rahul}}$. Naturally, one may think to study these ring extensions for some larger class. Note that if $R$ is an integral domain, then $R$ is in $\mathcal{H}$. This motivates us to study the ring extensions in class $\mathcal{H}$ having the same set of prime ideals. Note that if $R\in \mathcal{H}$, then there is a ring homomorphism $\phi$ from $T(R)$ to $R_{Nil(R)}$ given by $\phi(r/s) = r/s$, for all $r\in R$ and $s\in R\setminus Z(R)$. Moreover, the restriction of $\phi$ to $R$ is also a ring homomorphism given by $\phi(r) = r/1$, for all $r\in R$, see $\textup{\cite{badawi}}$. Note that if $R\in \mathcal{H}_0$, then $\phi(R) = R$.

In this paper, we discuss examples using this idealization theory. In $\textup{\cite{nagata}}$, Nagata defined a new class of rings, namely idealization of a module. If $R$ is a ring and $M$ is an $R$-module, then the idealization $R(+)M$ is the ring defined as follows: Its additive structure is that of the abelian group $R \oplus M$, and its multiplication is defined by $\left(r_1,m_1\right) \left(r_2,m_2\right) := \left(r_1 r_2,r_1 m_2 + r_2 m_1\right)$, for all $ r_1, r_2\in R$ and $m_1, m_2 \in M$. It will be convenient to view $R$ as a subring  of $R(+)M$ via the canonical injective ring homomorphism that sends $r$ to $\left( r,0\right)$.

\section{Results}

The ring extensions $R\subseteq T$ for which Spec($R$) = Spec($T$) are studied by Anderson and Dobbs in $\textup{\cite{rahul}}$. They discussed various type of domains $R$ for which there exist domains $T$ satisfying Spec($R$) = Spec($T$), and also given conditions on $R$ for which no such $T$ exists. In Lemma $3.1$, they proved that if $R\subseteq T$ is a ring extension with a common nonzero ideal $I$ which contains a non-zerodivisor of $R$, then $T$ is contained in the total quotient ring of $R$, and indeed $T\subseteq (I: I)$. The result does not seems to be correct. The proof goes wrong if for any nonzero $t \in T$ and for any non-zerodivisor $y \in I$ (of $R$), $ty = 0$. In fact, there is a class of ring extensions for which $ty = 0$ and hence thereby proving that $\textup{\cite[Lemma~3.1]{rahul}}$ is not correct. For example, consider $R = \mathbb{Z}$ and $T = \mathbb{Z}(+) \mathbb{Z}/2\mathbb{Z}$. Then $I = 2\mathbb{Z}$ is a common nonzero ideal with every nonzero element of $I$ is a non-zerodivisor of $R$. But $T\not\subseteq T(R)$. Note that this is true for all domains $R$ with nonzero maximal ideal $M$. In this case, take $T = R(+)R/M$. Then $M$ is a common nonzero ideal but $T\not\subseteq T(R)$. Note that if $ty \neq 0$, then the proof of $\textup{\cite[Lemma~3.1]{rahul}}$ will work. Thus, if we take $y$ to be a non-zerodivisor of $T$, then $\textup{\cite[Lemma~3.1]{rahul}}$ is correct with the same proof. We now state the modified result.

\begin{Lemma}\label{lem}
Let $R\subseteq T$ be a ring extension with a common nonzero ideal $I$. If $I$ contains a non-zerodivisor of $T$, then $T$ is contained in the total quotient ring of $R$, and indeed $T\subseteq (I: I)$.
\end{Lemma}

Also, in $\textup{\cite[Remark~3.4(a)]{rahul}}$, Anderson and Dobbs concluded the following:

Let $R\subset T$ be a ring extension such that Spec($R$) = Spec($T$) and $M$ be the unique maximal ideal of both $R$ and $T$. If $M$ contains a non-zerodivisor of $R$, then $T\subseteq (M: M)$.

They mentioned that the same proof of $\textup{\cite[Lemma~3.1]{rahul}}$ will work to prove the above result. As $\textup{\cite[Lemma~3.1]{rahul}}$ has been modified, some argument is needed to establish the above result. Note that $M$ contains a non-zerodivisor of $T$. Otherwise $M = Z(T)$ which gives $M = Z(R)$ as Spec($R$) = Spec($T$). This contradicts that $M$ contains a non-zerodivisor of $R$. Thus, the result follows by Lemma $\ref{lem}$.

Let $R\subset S$ be an extension of integral domains. Then $R$ is a said to be a pseudo-valuation subring (PV, for short) in $S$, see $\textup{\cite{ayache1}}$, if for each $x\in S\setminus R$ and non-unit $a$ of $R$, we have $x^{-1}a\in R$. Note that PV in a domain is a generalization of PVD. We now continue to investigate the properties of ring extensions having the same set of prime ideals initiated in $\textup{\cite{rahul}}$. Here is our first result.

\begin{Theorem}\label{r1}
Let $R\subseteq T\subset U$ be an extension of integral domains such that $Spec(R) = Spec(T)$. Then $R$ is a PV in $U$ if and only if $T$ is a PV in $U$.
\end{Theorem}

\begin{proof}
Assume that $R$ is a PV in $U$. Take $x\in U\setminus T$ and a non-unit $y$ in $T$. Then $y$ is also a non-unit in $R$ as Spec($R$) = Spec($T$). It follows that $x^{-1}y\in R$. Thus, $T$ is a PV in $U$. Conversely, suppose that $T$ is a PV in $U$. Take $x\in U\setminus R$ and $y$ a non-unit in $R$. Then $y$ is also a non-unit in $T$ as Spec($R$) = Spec($T$). If $x\notin T$, then $x^{-1}y\in T$ which is not a unit in $T$ as otherwise $x\in T$. It follows that $x^{-1}y\in R$ and we are done. On the other hand, if $x\in T$, then $x$ is a unit in $T$ as Spec($R$) = Spec($T$). It follows that $x^{-1}y$ is a non-unit in $T$ and so $x^{-1}y\in R$. Thus, $R$ is a PV in $U$.
\end{proof}

Note that the above proposition generalizes $\textup{\cite[Proposition~2.2]{rahul}}$ to PV in a domain. In Proposition $\ref{r3}$, we generalize the same to $\phi$-PVR as $\phi$-PVR generalizes PVD. Recall from $\textup{\cite{badawi}}$ that a ring $R \in \mathcal{H}$ is said to be a $\phi$-pseudo-valuation ring ($\phi$-PVR, for short), if each prime ideal $P$ of $R$ is $\phi$-strongly prime, that is, if $xy\in \phi(P)$ for $x, y\in R_{Nil(R)}$, then either $x\in \phi(P)$ or $y\in \phi(P)$. First, we need the following lemma.

\begin{Lemma}\label{r2}
Let $R$ be a ring in $\mathcal{H}$. If $T$ is a ring containing $R$ such that $Spec(R) = Spec(T)$, then $T\in \mathcal{H}$.
\end{Lemma}

\begin{proof}
Clearly, $Nil(T) = Nil(R)$ is a prime ideal of $T$. It remains to show that $Nil(R)$ is divided in $T$. Let $I$ be a proper ideal of $T$. Then $I$ is an ideal of $R$ as Spec($R$) = Spec($T$). Thus, $Nil(R)$ is divided in $T$.
\end{proof}

\begin{Proposition}\label{r3}
Let $R\subseteq T$ be a ring extension such that $R\in \mathcal{H}$ and $Spec(R) = Spec(T)$. Then $R$ is a $\phi$-PVR if and only if $T$ is a $\phi$-PVR.
\end{Proposition}

\begin{proof}
Note that $Nil(R)$ is a prime ideal of both $R$ and $T$. It follows that Spec($R/Nil(R)$) = Spec($T/Nil(R)$). Now, if $R$ is a $\phi$-PVR, then $R/Nil(R)$ is a PVD, by $\textup{\cite[Proposition~2.9]{badawi9}}$. It follows that $T/Nil(R)$ is a PVD, by $\textup{\cite[Proposition~2.2]{rahul}}$. Also, $T \in \mathcal{H}$, by Lemma $\ref{r2}$. Thus, by another appeal to $\textup{\cite[Proposition~2.9]{badawi9}}$, $T$ is a $\phi$-PVR. The converse is similar.
\end{proof}

A domain $R$ is called a going-down domain if $R\subseteq T$ satisfies the going-down property for each overring $T$ of $R$, see $\textup{\cite{dobbs}}$; and a ring $R$ is called a going-down ring if $R/P$ is a going-down domain for all prime ideals $P$ of $R$, see $\textup{\cite{dobbs1}}$. The next theorem generalizes $\textup{\cite[Corollary~A.4]{rahul}}$ and $\textup{\cite[Proposition~A.1]{rahul}}$ to rings in $\mathcal{H}$.

\begin{Theorem}\label{r20}
Let $R\in \mathcal{H}$. If $T$ is a ring containing $R$ properly such that $Spec(R) = Spec(T)$, then
\begin{enumerate}
\item[(i)] $R$ is not a going down ring.
\item[(ii)] $R$ is divided if and only if $T$ is divided.
\end{enumerate}
\end{Theorem}

\begin{proof}
To prove $(i)$, note that $Nil(R)$ is a prime ideal of both $R$ and $T$. It follows that Spec($R/Nil(R)$) = Spec($T/Nil(R)$). Now, if $R$ is a going down ring, then $R/Nil(R)$ is a going down domain, which is a contradiction by $\textup{\cite[Corollary~A.4]{rahul}}$.

For $(ii)$, let $R$ be divided and $P$ be any prime ideal of $T$. Then $P$ is a divided prime ideal of $R$. Let $I$ be a proper ideal of $T$. Then $I$ is a proper ideal of $R$ and so is comparable with $P$. Thus, $T$ is divided. Conversely, assume that $T$ is divided. Then $T/Nil(R)$ is divided and so $R/Nil(R)$ is divided, by $\textup{\cite[Proposition~A.1]{rahul}}$. Since $Nil(R)$ is divided, we conclude that $R$ is divided.
\end{proof}

A ring $R$ in $\mathcal{H}$ is said to be a $\phi$-chained ring, if for each $x\in R_{Nil(R)}\setminus \phi(R)$, we have $x^{-1}\in \phi(R)$, see $\textup{\cite{badawi6}}$. As a direct consequence of $\textup{\cite[Proposition~3.3]{badawi6}}$, $\textup{\cite[Corollary~3.4]{badawi6}}$, and $\textup{\cite[Proposition~3.8]{rahul}}$, we observe the generalization of $\textup{\cite[Proposition~2.5]{rahul}}$ to rings in $\mathcal{H}$.

\begin{Proposition}\label{r4}
Let $R\in \mathcal{H}$ be a local ring with maximal ideal $M$. If $M$ contains a non-zerodivisor, then the following are equivalent:
\begin{enumerate}
\item[(1)] $R$ is a $\phi$-PVR;
\item[(2)] $(M : M)$ is a $\phi$-chained overring of $R$ with maximal ideal $M$;
\item[(3)] $R$ has a $\phi$-chained overring $V$ such that $Spec(R) = Spec(V)$.
\end{enumerate}
\end{Proposition}

The next corollary is a companion to Theorem $\ref{r20}$ and Proposition $\ref{r4}$ which can be seen as a generalization of $\textup{\cite[Proposition~2.2]{hed}}$ to rings in $\mathcal{H}$.

\begin{Corollary}\label{r22}
Let $R\in \mathcal{H}$ be a local ring such that the maximal ideal contains a non-zerodivisor of $R$. Then $R$ is a $\phi$-chained ring if and only if $R$ is both a going down ring and a $\phi$-PVR.
\end{Corollary}

\begin{proof}
Let $R$ be a $\phi$-chained ring. Then it is a $\phi$-PVR trivially. Also, $R/Nil(R)$ is a valuation domain, by $\textup{\cite[Theorem~2.7]{badawi3}}$. It follows that $R/Nil(R)$ is a going down domain. Consequently, by $\textup{\cite[Proposition~2.1(a)]{dobbs1}}$, $R$ is going down.

Conversely,  assume that $R$ is both a going down ring and a $\phi$-PVR. Then by Proposition $\ref{r4}$, it follows that $R$ has a $\phi$-chained overring $V$ such that Spec($R$) = Spec($V$). Thus, $R = V$, by part $(i)$ of Theorem $\ref{r20}$.
\end{proof}

The next proposition was proved by Anderson and Dobbs for domains, see $\textup{\cite[Proposition~3.5(b)]{rahul}}$. Note that the same is true for rings also as the proof of $\textup{\cite[Proposition~3.5(b)]{rahul}}$ do not requires the ring to be a domain. Thus, the proof of next proposition follows mutatis mutandis from the proof of $\textup{\cite[Proposition~3.5(b)]{rahul}}$.

\begin{Proposition}\label{r5}
Let $R\subseteq T$ be rings such that $Spec(R) = Spec(T)$. Then $R_{P} = T_{T\setminus P}$ for each non-maximal prime ideal $P$ of $R$.
\end{Proposition}

The next proposition generalizes $\textup{\cite[Corollary~3.20(2)]{rahul}}$ to rings in $\mathcal{H}$.

\begin{Proposition}\label{r6}
Let $R\subseteq T$ be a ring extension such that $T\in \mathcal{H}$ and $R$ is local with maximal ideal $M$. Assume that $M$ is a finitely generated ideal of $R$. Then $Spec(R) = Spec(T)$ if and only if $M\in Spec(T)$.
\end{Proposition}

\begin{proof}
The ``only if'' assertion is trivial. For the converse part, assume that $M\in \mbox{Spec}(T)$. Then $Nil(T) = Nil(R)$. Also, $R/Nil(R)$ is local with finitely generated maximal ideal $M/Nil(R)$ and $T/Nil(R)$ is a domain. It follows that Spec($R/Nil(R)$) = Spec($T/Nil(R)$), by $\textup{\cite[Corollary~3.20]{rahul}}$. Thus, the result holds.
\end{proof}

\begin{Remark}\label{r7}
Let $R$ be a Pr\"ufer domain which is not a field. Then there is no ring $T$ which contains $R$ properly such that Spec($R$) = Spec($T$), see $\textup{\cite[Remark~3.17(b)]{rahul}}$. Note that this can be extended to $\phi$-Pr\"ufer rings also (A ring $R\in \mathcal{H}$ is said to be a $\phi$-Pr\"ufer ring, see $\textup{\cite{badawi3}}$, if $\phi(R)$ is a Pr\"ufer ring). More precisely, if $R$ is a $\phi$-Pr\"ufer ring such that $Nil(R)$ is not a maximal ideal of $R$ and $T$ is a ring containing $R$ such that Spec($R$) = Spec($T$), then $R = T$. To see this, first note that $Nil(R) = Nil(T)$. It follows that both $R/Nil(R)$ and $T/Nil(R)$ have the same set of prime ideals. As by $\textup{\cite[Theorem~2.6]{badawi3}}$, $R/Nil(R)$ is a Pr\"ufer domain,  we have $R/Nil(R)  = T/Nil(R)$ and so $R = T$.
\end{Remark}

In the next proposition, we investigate the overrings of $R \in \mathcal{H}_0$ having the same set of prime ideals. Moreover, it can be seen as a generalization of $\textup{\cite[Corollary~3.21]{rahul}}$.

\begin{Proposition}\label{r8}
Let $R\in \mathcal{H}_0$ be a local ring with maximal ideal $M$ such that $R'$ is a Pr\"ufer ring. Let $T$ be an overring of $R$. Then $Spec(R) = Spec(T)$ if and only if $M\in Spec(T)$.
\end{Proposition}

\begin{proof}
It is easy to see that both $R'$ and $T$ are in $\mathcal{H}_0$. Moreover, by $\textup{\cite{badawi}}$, $Nil(T) = Nil(R') = Nil(R)$ and $T(R/Nil(R)) = T(R)/Nil(R)$. It follows that $T/Nil(R)$ is an overring of $R/Nil(R)$. Moreover, by $\textup{\cite[Lemma~2.8]{badawi4}}$, $(R/Nil(R))' = R'/Nil(R)$. Now, by $\textup{\cite[Corollary~2.10]{badawi3}}$, $(R/Nil(R))'$ is a Pr\"ufer domain. Consequently, Spec($R/Nil(R)$) = Spec($T/Nil(R)$) if and only if $M/Nil(R)\in \mbox{Spec}(T/Nil(R))$, by $\textup{\cite[Corollary~3.21]{rahul}}$. Hence, the result follows.
\end{proof}

\begin{Remark}
In $\textup{\cite{rahul}}$, Anderson and Dobbs raised the following question: given a ring $R$, how does one can obtain all overrings $T$ of $R$ which satisfies Spec($R$) = Spec($T$). They answered the question in $\textup{\cite[Theorem~3.5]{rahul}}$ for domains. However, observe that the proof of $\textup{\cite[Theorem~3.5]{rahul}}$ works for rings also.
\end{Remark}

We now discuss an example that gives a class of ring extensions having the same set of prime ideals.

\begin{Example}\label{r11}
Let $B$ be a ring of the form $K + M$, where $K$ is a field and $M$ is a maximal ideal of $B$. Assume that $D$ is a subring of $K$, and
$A = D + M$. Set $U = A(+)T(B)$ and $V = B(+)T(B)$. Then Spec($U$) = Spec($V$) if and only if $V$ is local and $D$ is a field. To see this, note that if Spec($U$) = Spec($V$), then Spec($A$) = Spec($B$), by $\textup{\cite[Theorem~25.1(3)]{huckaba}}$. It follows that $B$ is local and $D$ is a field, by $\textup{\cite[Corollary~3.11]{rahul}}$. Consequently, $V$ is local, by another appeal to $\textup{\cite[Theorem~25.1(3)]{huckaba}}$. Conversely, assume that $V$ is local and $D$ is a field. Then $B$ is local, by $\textup{\cite[Theorem~25.1(3)]{huckaba}}$. It follows that Spec($A$) = Spec($B$), by $\textup{\cite[Corollary~3.11]{rahul}}$. Consequently, by another application of $\textup{\cite[Theorem~25.1(3)]{huckaba}}$, Spec($U$) = Spec($V$).
\end{Example}

The next proposition gives a class of rings $R$ in $\mathcal{H}_0$ which does not admit any ring $T$ properly containing $R$ and having the same set of prime ideals. For domains, it was observed in $\textup{\cite[pg.~372]{rahul}}$.

\begin{Proposition}\label{r12}
Let $R\in \mathcal{H}_0$ be an integrally closed local ring such that $R\neq T(R)$. If the maximal ideal $M$ of $R$ is finitely generated, then there is no ring $T$ properly containing $R$ such that $Spec(R) = Spec(T)$.
\end{Proposition}

\begin{proof}
If possible, suppose there exists a ring $T$ properly containing $R$ such that Spec($R$) = Spec($T$). Then we have $Nil(T) = Nil(R)$ and Spec($R/Nil(R)$) = Spec($T/Nil(R)$). Moreover, by $\textup{\cite{badawi}}$ and $\textup{\cite[Lemma~2.8]{badawi4}}$, we conclude that $T(R/Nil(R)) = T(R)/Nil(R)$ and $(R/Nil(R))' = R'/Nil(R)$, respectively. By hypothesis, it follows that $R/Nil(R)$ is an integrally closed local domain which is not a field. Moreover, the maximal ideal $M/Nil(R)$ of $R/Nil(R)$ is finitely generated. Consequently, we have a contradiction, by last paragraph of $\textup{\cite[pg.~372]{rahul}}$.
\end{proof}

The next example shows that the condition on $M$ of being finitely generated can not be removed.

\begin{Example}
Let $F\subset K$ be a purely transcendental field extension. Take a valuation domain $B = K + M$ with unique nonzero maximal ideal $M$. Set $R = A(+)T(B)$ and $T = B(+)T(B)$ where $A = F + M$. Then it is easy to see that $Z(R) = Nil(R) = Nil(T) = \{0\}(+)T(B)$. Let $(0, x)\in Nil(R)$ and $(y, z)\in R\setminus Nil(R)$. Then $y\neq 0$ and $(0, x) = (y, z)(0, x/y)$. It follows that $R\in \mathcal{H}_0$. Since $A$ is an integrally closed domain, $R$ is an integrally closed ring, by $\textup{\cite[Theorem~25.6]{huckaba}}$. Also, by $\textup{\cite[Corollary~25.5(3)]{huckaba}}$, $R\neq T(R)$. Since $M$ is not finitely generated, $M(+)T(B)$ is not finitely generated. Note that $M(+)T(B)$ is the unique maximal ideal of $R$, by $\textup{\cite[Theorem~25.1(3)]{huckaba}}$. However, by Example $\ref{r11}$, Spec($R$) = Spec($T$).
\end{Example}

Let $C(R)$ denotes the complete integral closure of $R$ in $T(R)$. In the next theorem, we show that the ring extensions in $\mathcal{H}_0$ having the same set of prime ideals will have the same complete integral closure. This generalizes $\textup{\cite[Proposition~3.15]{rahul}}$ for rings in $\mathcal{H}_0$.

\begin{Theorem}\label{r13}
Let $R \in \mathcal{H}_0$ be such that $R\neq T(R)$. If $T$ is a ring containing $R$ such that $Spec(R) = Spec(T)$, then $C(R) = C(T)$.
\end{Theorem}

\begin{proof}
It is easy to conclude that $Nil(R) = Nil(T)$, $Z(T) = Z(R)$ and Spec($R/Nil(R)$) = Spec($T/Nil(R)$). Also, $T(R/Nil(R)) = T(R)/Nil(R)$, by $\textup{\cite{badawi}}$, and so $R/Nil(R)$ is not a field. By $\textup{\cite[Proposition~3.15]{rahul}}$, it follows that $C(R/Nil(R)) = C(T/Nil(R))$. Now, by Lemma $\ref{r2}$, $T\in \mathcal{H}_0$. Consequently, we have $C(R)/Nil(R) = C(T)/Nil(R)$, by $\textup{\cite[Lemma~2.8]{badawi4}}$. Thus, $C(R) = C(T)$.
\end{proof}

The following corollary is an immediate consequence of Theorem $\ref{r13}$ when $R = C(R)$.

\begin{Corollary}\label{r14}
Let $R\in \mathcal{H}_0$ be a completely integrally closed ring such that $R\neq T(R)$. Then there is no ring $T$ properly containing $R$ such that $Spec(R) = Spec(T)$.
\end{Corollary}

In $\textup{\cite[Lemma~2.4]{rahul}}$, Anderson and Dobbs proved the next result for local domains which are not fields. Here, we generalize the same for rings. In fact, the proof follows mutatis mutandis from the proof of $\textup{\cite[Lemma~2.4]{rahul}}$. Note that an ideal $I$ of a ring $R$ is said to be divisorial if $I = (R : (R: I))$.

\begin{Theorem}\label{r17}
Let $R$ be a local ring with maximal ideal $M$. If $M$ contains a non-zerodivisor of $R$, then the following hold:
\begin{enumerate}
\item[(i)] $(R : M) = (M : M)$ if and only if $M$ is not a principal ideal of $R$. Moreover, if $M = Rr$ for some $r\in R$, then $(R : M) = Rr^{-1}$ and $(M : M) = R$.
\item[(ii)] If $M$ is not a divisorial ideal of $R$, then $(R : M) = (M : M)$.
\end{enumerate}
\end{Theorem}

The next proposition generalizes $\textup{\cite[Proposition~3.23]{rahul}}$.

\begin{Proposition}\label{r18}
Let $R$ be a local ring with maximal ideal $M$ such that $M$ contains a non-zerodivisor of $R$. Then $(M : M)\neq R$ if and only if $M$ is a non-principal divisorial ideal of $R$.
\end{Proposition}

\begin{proof}
The ``only if'' assertion follows mutatis mutandis from the proof of $\textup{\cite[Proposition~3.23]{rahul}}$ and ``if'' assertion follows from Theorem $\ref{r17}$.
\end{proof}

The proof of the next corollary follows directly from Proposition $\ref{r18}$ and Lemma $\ref{lem}$.

\begin{Corollary}\label{r19}
Let $R$ be a local ring such that the maximal ideal of $R$ contains a non-zerodivisor of $R$. Assume that the maximal ideal of $R$ is either principal or a non-divisorial ideal of $R$. Then there is no ring $T$ properly containing $R$ such that $Spec(R) = Spec(T)$.
\end{Corollary}

It is easy to see that $\textup{\cite[Corollary~3.30]{rahul}}$ is not true for fields as if $R$ is a field, then $R'$ is not of rank one. However, if we exempt $R$ to be a field, then $R'$ is a discrete (rank 1) valuation domain. In Proposition $\ref{r23}$, we generalize the same to rings in $\mathcal{H}_0$. First, we recall some definitions. An ideal of a ring $R$ is said to be a nonnil ideal if $I\not \subseteq Nil(R)$. A ring $R\in \mathcal{H}$ is called a nonnil-Noetherian ring if every nonnil ideal of $R$ is finitely generated, see $\textup{\cite{badawi10}}$. Recall from $\textup{\cite{badawi4}}$ that a ring $R\in \mathcal{H}$ is said to be a discrete $\phi$-chained ring if $R$ is a $\phi$-chained ring with at most one nonnil prime ideal and every nonnil ideal of $R$ is a principal ideal.

\begin{Proposition}\label{r23}
Let $R\in \mathcal{H}_0$. If $R$ is a nonnil-Noetherian $\phi$-PVR, then dim($R)\leq 1$ and $R'$ is a discrete $\phi$-chained ring.
\end{Proposition}

\begin{proof}
Note that $(R/Nil(R))' = R'/Nil(R)$, by $\textup{\cite[Lemma~2.8]{badawi4}}$. Also, by $\textup{\cite[Theorem~2.2]{badawi10}}$ and $\textup{\cite[Proposition~2.9]{badawi9}}$, $R/Nil(R)$ is a Noetherian PVD. It follows that $dim(R) \leq 1$ and $R'/Nil(R)$ is a discrete valuation domain, by $\textup{\cite[Corollary~3.30]{rahul}}$. Thus, the result holds, by $\textup{\cite[Lemma~2.9]{badawi4}}$.
\end{proof}

Recall from $\textup{\cite{m}}$ that a ring $R$ is said to be a weakly finite-conductor ring if $Ra\cap Rb$ is a finitely generated ideal for all $a, b\in R$.  We end this paper with the generalization of $\textup{\cite[Corollary~A.3]{rahul}}$ to rings in $\mathcal{H}_0$. Note that $\textup{\cite[Corollary~A.3]{rahul}}$ fails to hold when $R\subset T$ is a field extension such that $T$ is not a finitely generated $R$-module. Therefore, the exemption of $R$ to be a field is required in the statement of $\textup{\cite[Corollary~A.3]{rahul}}$.

\begin{Proposition}\label{r24}
Let $R\subset T$ be a ring extension such that $R\in \mathcal{H}_0$ is a weakly finite-conductor ring, $R\neq T(R)$, and $Spec(R) = Spec(T)$. Let $M$ be the maximal ideal of $R$. Then both $M$ and $T$ are finitely generated $R$-modules.
\end{Proposition}

\begin{proof}
Note that $R/Nil(R)$ is not a field as $T(R/Nil(R)) = T(R)/Nil(R)$, by $\textup{\cite{badawi}}$. Also, it is easy to see that $R/Nil(R)$ is a finite-conductor domain with maximal ideal $M/Nil(R)$. Thus, $M/Nil(R)$ and $T/Nil(R)$ are finitely generated $R/Nil(R)$-module, by $\textup{\cite[Corollary~A.3]{rahul}}$. Since $R\neq T(R)$, $M$ is a nonnil ideal of $R$. Now, the proof follows mutatis mutandis from the proof of $\textup{\cite[Lemma~2.4]{badawi3}}$.
\end{proof}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

   In the \emph{nucleated instability\/} (also called core
   instability) hypothesis of giant planet
   formation, a critical mass for static core  envelope
   protoplanets has been found. \citet{mizuno} determined
   the critical mass of the core to be about $12 \,M_\oplus$
   ($M_\oplus=5.975 \times 10^{27}\,\mathrm{g}$ is the Earth mass), which
   is independent of the outer boundary
   conditions and therefore independent of the location in the
   solar nebula. This critical value for the core mass corresponds
   closely to the cores of today's giant planets.

   Although no hydrodynamical study has been available many workers
   conjectured that a collapse or rapid contraction will ensue
   after accumulating the critical mass. The main motivation for
   this article
   is to investigate the stability of the static envelope at the
   critical mass. With this aim the local, linear stability of static
   radiative gas  spheres is investigated on the basis of Baker's
   (\citeyear{baker}) standard one-zone model.

   Phenomena similar to the ones described above for giant planet
   formation have been found in hydrodynamical models concerning
   star formation where protostellar cores explode
   (Tscharnuter \citeyear{tscharnuter}, Balluch \citeyear{balluch}),
   whereas earlier studies found quasi-steady collapse flows. The
   similarities in the (micro)physics, i.e., constitutive relations of
   protostellar cores and protogiant planets serve as a further
   motivation for this study.

\section{Baker's standard one-zone model}

   \begin{figure*}
   \centering

   \caption{Adiabatic exponent $\Gamma_1$.
               $\Gamma_1$ is plotted as a function of
               $\lg$ internal energy $\mathrm{[erg\,g^{-1}]}$ and $\lg$
               density $\mathrm{[g\,cm^{-3}]}$.}
              \label{FigGam}%
    \end{figure*}
   In this section the one-zone model of \citet{baker},
   originally used to study the Cephe{\"{\i}}d pulsation mechanism, will
   be briefly reviewed. The resulting stability criteria will be
   rewritten in terms of local state variables, local timescales and
   constitutive relations.

   \citet{baker} investigates the stability of thin layers in
   self-gravitating,
   spherical gas clouds with the following properties:
   \begin{itemize}
      \item hydrostatic equilibrium,
      \item thermal equilibrium,
      \item energy transport by grey radiation diffusion.
   \end{itemize}
   For the one-zone-model Baker obtains necessary conditions
   for dynamical, secular and vibrational (or pulsational)
   stability (Eqs.\ (34a,\,b,\,c) in Baker \citeyear{baker}). Using Baker's
   notation:
   \[
      \begin{array}{lp{0.8\linewidth}}
         M_{r}  & mass internal to the radius $r$     \\
         m               & mass of the zone                    \\
         r_0             & unperturbed zone radius             \\
         \rho_0          & unperturbed density in the zone     \\
         T_0             & unperturbed temperature in the zone \\
         L_{r0}          & unperturbed luminosity              \\
         E_{\mathrm{th}} & thermal energy of the zone
      \end{array}
   \]
\noindent
   and with the definitions of the \emph{local cooling time\/}
   (see Fig.~\ref{FigGam})
   \begin{equation}
      \tau_{\mathrm{co}} = \frac{E_{\mathrm{th}}}{L_{r0}} \,,
   \end{equation}
   and the \emph{local free-fall time}
   \begin{equation}
      \tau_{\mathrm{ff}} =
         \sqrt{ \frac{3 \pi}{32 G} \frac{4\pi r_0^3}{3 M_{\mathrm{r}}}
}\,,
   \end{equation}
   Baker's $K$ and $\sigma_0$ have the following form:
   \begin{eqnarray}
      \sigma_0 & = & \frac{\pi}{\sqrt{8}}
                     \frac{1}{ \tau_{\mathrm{ff}}} \\
      K        & = & \frac{\sqrt{32}}{\pi} \frac{1}{\delta}
                        \frac{ \tau_{\mathrm{ff}} }
                             { \tau_{\mathrm{co}} }\,;
   \end{eqnarray}
   where $ E_{\mathrm{th}} \approx m (P_0/{\rho_0})$ has been used and
   \begin{equation}
   \begin{array}{l}
      \delta = - \left(
                    \frac{ \partial \ln \rho }{ \partial \ln T }
                 \right)_P \\
      e=mc^2
   \end{array}
   \end{equation}
   is a thermodynamical quantity which is of order $1$ and equal to $1$
   for nonreacting mixtures of classical perfect gases. The physical
   meaning of $ \sigma_0 $ and $K$ is clearly visible in the equations
   above. $\sigma_0$ represents a frequency of the order one per
   free-fall time. $K$ is proportional to the ratio of the free-fall
   time and the cooling time. Substituting into Baker's criteria, using
   thermodynamic identities and definitions of thermodynamic quantities,
   \begin{displaymath}
      \Gamma_1      = \left( \frac{ \partial \ln P}{ \partial\ln \rho}
                           \right)_{S}    \, , \;
      \chi^{}_\rho  = \left( \frac{ \partial \ln P}{ \partial\ln \rho}
                           \right)_{T}    \, , \;
      \kappa^{}_{P} = \left( \frac{ \partial \ln \kappa}{ \partial\ln P}
                           \right)_{T}
   \end{displaymath}
   \begin{displaymath}
      \nabla_{\mathrm{ad}} = \left( \frac{ \partial \ln T}
                             { \partial\ln P} \right)_{S} \, , \;
      \chi^{}_T       = \left( \frac{ \partial \ln P}
                             { \partial\ln T} \right)_{\rho} \, , \;
      \kappa^{}_{T}   = \left( \frac{ \partial \ln \kappa}
                             { \partial\ln T} \right)_{T}
   \end{displaymath}
   one obtains, after some pages of algebra, the conditions for
   \emph{stability\/} given
   below:
   \begin{eqnarray}
      \frac{\pi^2}{8} \frac{1}{\tau_{\mathrm{ff}}^2}
                ( 3 \Gamma_1 - 4 )
         & > & 0 \label{ZSDynSta} \\
      \frac{\pi^2}{\tau_{\mathrm{co}}
                   \tau_{\mathrm{ff}}^2}
                   \Gamma_1 \nabla_{\mathrm{ad}}
                   \left[ \frac{ 1- 3/4 \chi^{}_\rho }{ \chi^{}_T }
                          ( \kappa^{}_T - 4 )
                        + \kappa^{}_P + 1
                   \right]
        & > & 0 \label{ZSSecSta} \\
     \frac{\pi^2}{4} \frac{3}{\tau_{ \mathrm{co} }
                              \tau_{ \mathrm{ff} }^2
                             }
         \Gamma_1^2 \, \nabla_{\mathrm{ad}} \left[
                                   4 \nabla_{\mathrm{ad}}
                                   - ( \nabla_{\mathrm{ad}} \kappa^{}_T
                                     + \kappa^{}_P
                                     )
                                   - \frac{4}{3 \Gamma_1}
                                \right]
        & > & 0   \label{ZSVibSta}
   \end{eqnarray}
   For a physical discussion of the stability criteria see \citet{baker} or \citet{cox}.

   We observe that these criteria for dynamical, secular and
   vibrational stability, respectively, can be factorized into
   \begin{enumerate}
      \item a factor containing local timescales only,
      \item a factor containing only constitutive relations and
         their derivatives.
   \end{enumerate}
   The first factors, depending on only timescales, are positive
   by definition. The signs of the left hand sides of the
   inequalities~(\ref{ZSDynSta}), (\ref{ZSSecSta}) and (\ref{ZSVibSta})
   therefore depend exclusively on the second factors containing
   the constitutive relations. Since they depend only
   on state variables, the stability criteria themselves are \emph{
   functions of the thermodynamic state in the local zone}. The
   one-zone stability can therefore be determined
   from a simple equation of state, given for example, as a function
   of density and
   temperature. Once the microphysics, i.e.\ the thermodynamics
   and opacities (see Table~\ref{KapSou}), are specified (in practice
   by specifying a chemical composition) the one-zone stability can
   be inferred if the thermodynamic state is specified.
   The zone -- or in
   other words the layer -- will be stable or unstable in
   whatever object it is imbedded as long as it satisfies the
   one-zone-model assumptions. Only the specific growth rates
   (depending upon the time scales) will be different for layers
   in different objects.

   \begin{table}
      \caption[]{Opacity sources.}
         \label{KapSou}
     $$
         \begin{array}{p{0.5\linewidth}l}
            \hline
            \noalign{\smallskip}
            Source      &  T / {[\mathrm{K}]} \\
            \noalign{\smallskip}
            \hline
            \noalign{\smallskip}
            Yorke 1979, Yorke 1980a & \leq 1700^{\mathrm{a}}     \\
            Kr\"ugel 1971           & 1700 \leq T \leq 5000 \\
            Cox \& Stewart 1969     & 5000 \leq             \\
            \noalign{\smallskip}
            \hline
         \end{array}
     $$
   \end{table}
   We will now write down the sign (and therefore stability)
   determining parts of the left-hand sides of the inequalities
   (\ref{ZSDynSta}), (\ref{ZSSecSta}) and (\ref{ZSVibSta}) and thereby
   obtain \emph{stability equations of state}.

   The sign determining part of inequality~(\ref{ZSDynSta}) is
   $3\Gamma_1 - 4$ and it reduces to the
   criterion for dynamical stability
   \begin{equation}
     \Gamma_1 > \frac{4}{3}\,\cdot
   \end{equation}
   Stability of the thermodynamical equilibrium demands
   \begin{equation}
      \chi^{}_\rho > 0, \;\;  c_v > 0\, ,
   \end{equation}
   and
   \begin{equation}
      \chi^{}_T > 0
   \end{equation}
   holds for a wide range of physical situations.
   With
   \begin{eqnarray}
      \Gamma_3 - 1 = \frac{P}{\rho T} \frac{\chi^{}_T}{c_v}&>&0\\
      \Gamma_1     = \chi_\rho^{} + \chi_T^{} (\Gamma_3 -1)&>&0\\
      \nabla_{\mathrm{ad}}  = \frac{\Gamma_3 - 1}{\Gamma_1}         &>&0
   \end{eqnarray}
   we find the sign determining terms in inequalities~(\ref{ZSSecSta})
   and (\ref{ZSVibSta}) respectively and obtain the following form
   of the criteria for dynamical, secular and vibrational
   \emph{stability}, respectively:
   \begin{eqnarray}
      3 \Gamma_1 - 4 =: S_{\mathrm{dyn}}      > & 0 & \label{DynSta}  \\
      \frac{ 1- 3/4 \chi^{}_\rho }{ \chi^{}_T } ( \kappa^{}_T - 4 )
         + \kappa^{}_P + 1 =: S_{\mathrm{sec}} > & 0 & \label{SecSta} \\
      4 \nabla_{\mathrm{ad}} - (\nabla_{\mathrm{ad}} \kappa^{}_T
                             + \kappa^{}_P)
                             - \frac{4}{3 \Gamma_1} =: S_{\mathrm{vib}}
                                      > & 0\,.& \label{VibSta}
   \end{eqnarray}
   The constitutive relations are to be evaluated for the
   unperturbed thermodynamic state (say $(\rho_0, T_0)$) of the zone.
   We see that the one-zone stability of the layer depends only on
   the constitutive relations $\Gamma_1$,
   $\nabla_{\mathrm{ad}}$, $\chi_T^{},\,\chi_\rho^{}$,
   $\kappa_P^{},\,\kappa_T^{}$.
   These depend only on the unperturbed
   thermodynamical state of the layer. Therefore the above relations
   define the one-zone-stability equations of state
   $S_{\mathrm{dyn}},\,S_{\mathrm{sec}}$
   and $S_{\mathrm{vib}}$. See Fig.~\ref{FigVibStab} for a picture of
   $S_{\mathrm{vib}}$. Regions of secular instability are
   listed in Table~1.

   \begin{figure}
   \centering

      \caption{Vibrational stability equation of state
               $S_{\mathrm{vib}}(\lg e, \lg \rho)$.
               $>0$ means vibrational stability.
              }
         \label{FigVibStab}
   \end{figure}

\section{Conclusions}

   \begin{enumerate}
      \item The conditions for the stability of static, radiative
         layers in gas spheres, as described by Baker's (\citeyear{baker})
         standard one-zone model, can be expressed as stability
         equations of state. These stability equations of state depend
         only on the local thermodynamic state of the layer.
      \item If the constitutive relations -- equations of state and
         Rosseland mean opacities -- are specified, the stability
         equations of state can be evaluated without specifying
         properties of the layer.
      \item For solar composition gas the $\kappa$-mechanism is
         working in the regions of the ice and dust features
         in the opacities, the $\mathrm{H}_2$ dissociation and the
         combined H, first He ionization zone, as
         indicated by vibrational instability. These regions
         of instability are much larger in extent and degree of
         instability than the second He ionization zone
         that drives the Cephe{\"\i}d pulsations.
   \end{enumerate}

\begin{acknowledgements}
      Part of this work was supported by the German
      \emph{Deut\-sche For\-schungs\-ge\-mein\-schaft, DFG\/} project
      number Ts~17/2--1.
\end{acknowledgements}

\section{Introduction}

HR~8799 (also called HD~218396 and HIP~114189) is a nearby \citep[$39.4 \pm 1.1$ pc;][]{2007A&A...474..653V}, young \citep[$30^{+20}_{-10}$~Myr assuming member of the Columba association;][]{2010Natur.468.1080M} A5V (up to F0V) star with mass $1.5 \pm 0.3\; M_{\odot}$ \citep{1999AJ....118.2993G, 2003AJ....126.2048G}. HR~8799 has four confirmed, directly imaged planets \citep{2008Sci...322.1348M, 2010Natur.468.1080M} and a debris disk \citep{2009ApJ...705..314S,2014ApJ...780...97M}. The planet masses are estimated to be $5$~(b) and $7$~(cde)~$M_{\text{Jup}}$ \citep{2010Natur.468.1080M}. The on-sky separations between planets and star are estimated to be 67.9~(b), 38.0~(c), 24.5~(d), and 14.5~(e) AU. \cite{2014ApJ...780...97M} estimated the inclination of the system with respect to the plane of the sky to be $26 \pm 3\ensuremath{{}^{\circ}}\xspace$ by observing the outer part of the debris disk and assuming coplanarity with the planetary orbits. They estimate the position angle of the inclination to be $64 \pm 3\ensuremath{{}^{\circ}}\xspace$ (G.~Kennedy, private communication). \figref{fig:obs_8799} shows all observations of the planets published at the time of writing.

\begin{figure}
\centering
\includegraphics[scale=1]{tikz/obs_8799.pdf}
\caption{Observations of the HR~8799 planets. The black dashed line corresponds to the position angle of the inclination as estimated by \cite{2014ApJ...780...97M} with the $\pm 1\sigma$ uncertainty shown by the grey region.}
\label{fig:obs_8799}
\end{figure}

The estimates of the planetary masses in the HR~8799 system depend on the age of the system through cooling models, indicating higher planet masses with higher age \citep{2003A&A...402..701B}. The age of the HR~8799 system is uncertain and has been estimated using different techniques to range between 30~Myr and 1~Gyr \citep[see][and references therein]{2010Natur.468.1080M,2008Sci...322.1348M}. The planet mass estimates reach the brown dwarf regime when an age of $\gtrsim 100$~Myr is assumed. In this paper we consider HR~8799~bcde as planets; the case of brown dwarfs is more closely investigated by \cite{2010ApJ...721L.199M}. We assume the age 30~Myr in this paper as it is estimated from the Columba association. An age of 30~Myr leads to multi-Jovian masses that, however, are still in the planetary regime. Whether it is true that HR~8799 is associated with the Columba association could be revealed by Gaia in the near future \citep{2001A&A...369..339P}. Other uncertainties in system properties include distance \citep[$39.4\pm1.1$\,pc][]{2007A&A...474..653V} and stellar mass \citep[$1.5\pm0.3\mathrm{\,M}_\odot$][]{1999AJ....118.2993G,2003AJ....126.2048G}.

The formation of the massive HR~8799~bcde planets can be explained by the recently developed models of pebble accretion \citep{2012A&A...544A..32L}, but the dynamics of the system has been a puzzle unless resonances are considered. The orbits of the HR~8799~bcde planets are difficult to constrain as the planets have been observed during a small fraction of their orbits (see \figref{fig:obs_8799}). In recent work by \cite{2015ApJ...803...31P}, however,  fits are made to constrain the orbital elements and a variety of solutions are found, favouring a slightly more eccentric orbit for planet d. It is reasonable to believe that the planet orbits are close to circular, as with high eccentricity the orbits are likely to cross. However, the HR~8799 planetary system appears to be unstable on a timescale much shorter than the estimated age of the system when integrating the orbits of the planets, assumed initially circular, seen pole-on and non-resonant \citep{2010ApJ...710.1408F, 2009A&A...503..247R, 2012ApJ...755...38S}. A possible and popular explanation to the survival of the HR~8799 planetary system is that the planets are locked deep in strong mean motion resonances that can lead to stability timescales up to 1~Gyr \citep[e.g. ][]{2010ApJ...710.1408F, 2009MNRAS.397L..16G, 2009A&A...503..247R, 2010ApJ...721L.199M, 2011ApJ...729..128C, 2011ApJ...741...55S, 2013A&A...549A..52E,  2014MNRAS.440.3140G}. These resonant configurations are characterised by small-amplitude librations of one or more resonant arguments corresponding to two-, three-, or four-planet mean motion resonances, thus ensuring that systems remain protected for extended periods from close encounters between the planets. \cite{2013MNRAS.430..320M} find stabilising effects from the outer debris disk, although this requires a very massive debris disk.

In this paper we show simulations of planetary systems that look like the HR~8799 system, are stable for longer than the estimated age, and are not stabilised by strong resonant lock. We create these systems by starting the planetary orbits slightly more widely separated than the observed on-sky separation assuming circular orbits. A more widely separated configuration can also look like the HR~8799 system as the orbits are not exactly circular. Larger separation between the planetary orbits results in a much longer stability timescale, which makes systems remain stable for longer than the estimated age of HR~8799. This solution is different from previously published explanations, as we do not need protective resonant lock or any other stabilising effects to create a long-lived system.

In \secref{sec:simulations} we describe how we simulate planetary systems. \secreftwo{sec:HR8799unstable}{sec:oursol} present first results consistent with earlier work and later long-lived solutions in agreement with the observed system. In \secref{sec:discussion} we introduce a concept of evolutionary phases in the lives of wide-orbit planetary systems. These phases could be used to better understand the past and future of a wide-orbit planetary system. We summarise our findings and conclusions in \secref{sec:conclusions}.

\section{Simulations}\label{sec:simulations}

\begin{table*}
\caption{Initial conditions of the simulations.}
\label{tab:simulations}
\centering
\begin{tabular}{l R{9mm} R{9mm} R{9mm} R{10mm} R{10mm} R{10mm} R{10mm} L{55mm}}
\toprule\midrule
Name & \multicolumn{3}{c}{Separation between orbits [$R_H$]} & \multicolumn{4}{c}{Initial semi-major axis [AU]} & Comment \\
& $\Delta _{\mathrm{bc}}$ & $\Delta _{\mathrm{cd}}$ & $\Delta _{\mathrm{de}}$ & $a_{\mathrm{b}}$ & $a_{\mathrm{c}}$ & $a_{\mathrm{d}}$ & $a_{\mathrm{e}}$ & \\
\midrule
\textbf{1} & 4.14 & 3.02 & 3.55 & 67.9 & 38.0 & 24.5 & 14.5 & Face-on assumption. (\secref{sec:sim1})\\ \\
\textbf{2} & 3.47 & 3.31 & 3.90 & 67.9 & 41.9 & 25.8 & 14.5 & Observed inclination, $i = 26\ensuremath{{}^{\circ}}\xspace$, position angle$ = 64\ensuremath{{}^{\circ}}\xspace$. (\secref{sec:sim2})\\ \\

\textbf{2$\alpha$} &&&&&&&& Same initial conditions as system~8 in simulation~2. Differ in output rates and thus timestep. (\secref{sec:sim2and4alpha}) \\ \\

\textbf{3} & 3.56 & 3.56 & 3.56 & 68.8 & 41.9 & 24.8 & 14.7 & Optimal inclination, $i = 25\ensuremath{{}^{\circ}}\xspace$, position angle$ = 42\ensuremath{{}^{\circ}}\xspace$. (\secref{sec:sim3}) \\ \\
\textbf{4} &&&&&&&& Equal initial orbit separation in $\Delta$. Planet~e fixed at 14.3~AU. (\secref{sec:sim4})\\
$\quad$a & 3.6 & 3.6 & 3.6 & 68.1 & 41.2 & 24.3 & 14.3 \\
$\quad$b & 3.65 & 3.65 & 3.65 & 69.7 & 41.9 & 24.5 & 14.3\\
$\quad$c & 3.7 & 3.7 & 3.7 & 71.2 & 42.5 & 24.7 & 14.3\\
$\quad$d & 3.75 & 3.75 & 3.75 & 72.9 & 43.2 & 24.9 & 14.3\\
$\quad$e & 3.8 & 3.8 & 3.8 & 74.6 & 43.9 & 25.0 & 14.3\\
$\quad$f & 3.85 & 3.85 & 3.85 & 76.3 & 44.6 & 25.2 & 14.3\\
$\quad$g & 3.9 & 3.9 & 3.9 & 78.1 & 45.3 & 25.4 & 14.3\\
$\quad$h & 3.95 & 3.95 & 3.95 & 79.9 & 46.0 & 25.6 & 14.3\\
$\quad$i & 4 & 4 & 4 & 81.8 & 46.7 & 25.8 & 14.3\\
$\quad$f\_long & 3.85 & 3.85 & 3.85 & 76.3 & 44.6 & 25.2 & 14.3 & Initialised in the same way as 4f, but not stopping after close encounter. Run time 100~Myr. (\secref{sec:discussion})\\ \\

\textbf{4$\alpha$} &&&&&&&& Same initial conditions as system~94 in simulation~4e. Differ in output rates and thus timestep. (\secref{sec:sim2and4alpha}) \\ \\

\textbf{5} &&&&&&&& Same initial conditions as simulation~4, but planet masses are 4~(b) and 6~(cde)~$M_{\mathrm{Jup}}$. (\secref{sec:sim5and6})\\
$\quad$a & 3.83 & 3.79 & 3.79 & 68.1 & 41.2 & 24.3 & 14.3 \\
$\quad$b & 3.88 & 3.84 & 3.84 & 69.7 & 41.9 & 24.5 & 14.3\\
$\quad$c & 3.93 & 3.90 & 3.90 & 71.2 & 42.5 & 24.7 & 14.3\\
$\quad$d & 3.99 & 3.95 & 3.95 & 72.9 & 43.2 & 24.9 & 14.3\\
$\quad$e & 4.04 & 4.00 & 4.00 & 74.6 & 43.9 & 25.0 & 14.3\\
$\quad$f & 4.09 & 4.05 & 4.05 & 76.3 & 44.6 & 25.2 & 14.3\\
$\quad$g & 4.14 & 4.11 & 4.11 & 78.1 & 45.3 & 25.4 & 14.3\\
$\quad$h & 4.20 & 4.16 & 4.16 & 79.9 & 46.0 & 25.6 & 14.3\\
$\quad$i & 4.25 & 4.21 & 4.21 & 81.8 & 46.7 & 25.8 & 14.3\\ \\

\textbf{6} &&&&&&&& Same initial conditions as simulation~4, but planet masses are 3~(b) and 5~(cde)~$M_{\mathrm{Jup}}$. (\secref{sec:sim5and6})\\
$\quad$a & 4.12 & 4.03 & 4.03 & 68.1 & 41.2 & 24.3 & 14.3 \\
$\quad$b & 4.18 & 4.08 & 4.08 & 69.7 & 41.9 & 24.5 & 14.3\\
$\quad$c & 4.24 & 4.14 & 4.14 & 71.2 & 42.5 & 24.7 & 14.3\\
$\quad$d & 4.29 & 4.20 & 4.20 & 72.9 & 43.2 & 24.9 & 14.3\\
$\quad$e & 4.35 & 4.25 & 4.25 & 74.6 & 43.9 & 25.0 & 14.3\\
$\quad$f & 4.41 & 4.31 & 4.31 & 76.3 & 44.6 & 25.2 & 14.3\\
$\quad$g & 4.46 & 4.36 & 4.36 & 78.1 & 45.3 & 25.4 & 14.3\\
$\quad$h & 4.52 & 4.42 & 4.42 & 79.9 & 46.0 & 25.6 & 14.3\\
$\quad$i & 4.58 & 4.47 & 4.47 & 81.8 & 46.7 & 25.8 & 14.3\\ \\
\bottomrule
\end{tabular}
\end{table*}

\subsection{Orbital parameters}\label{sec:orbital_parameters}

All planets in all our simulations are analogues of the HR~8799~bcde planets and they are initiated on circular, heliocentric orbits.
 The initial orbit inclination we find by drawing a random inclination between 0\ensuremath{{}^{\circ}}\xspace\ and 5\ensuremath{{}^{\circ}}\xspace\ and random longitude of ascending node between 0 and 360\ensuremath{{}^{\circ}}\xspace\ \citep[$\beta = 5\ensuremath{{}^{\circ}}\xspace$; see][]{2012ApJ...758...39J}. The initial planet position on the orbit (true anomaly) is randomised. We assign the mass $1.5\; M_{\odot}$ to the star and the planets are given the estimated values of the masses, $5\; M_{\text{Jup}}$ (b) and $7\; M_{\text{Jup}}$ (cde) if nothing else is stated (the masses are decreased in simulations~5 and 6).
 The planets orbit the barycentre of the system, which is not centred on the star as a result of the high planet
mass. This means that the planetary orbits have a low initial eccentricity as they are initiated on orbits circular around the star. The initial orbital eccentricities vary between different systems as the inclinations, longitude of ascending nodes, and true anomalies are randomised.

We vary the separation between the initial planet orbits (regarding the planets as test particles) in terms of the mutual Hill radius, which is defined as follows:
\begin{equation}\label{eq:RH}
R_H = \frac{r_{p1} + r_{p2}}{2}\left( \frac{M_{p1} + M_{p2}}{3M_{\star}} \right)^{1/3} ,
\end{equation}
where $r$ signifies star-planet separation in the orbital plane and $M$ mass \citep{1996Icar..119..261C}. The indices $p1$ and $p2$ refer to the two planets and $\star$ refers to the star. We then write the separation between the initial, circular orbits  in mutual Hill radii as
\begin{equation}\label{eq:Delta}
\Delta = (r_{p2} - r_{p1})/R_H ,
\end{equation}
where in this case $p2$ refers to the outer planet and $p1$ to the inner planet. For a given planetary system there is a critical $\Delta$ at which in 50\% of the cases the system leads to orbit crossing within a given time \citep{1996Icar..119..261C,2008ApJ...686..580C}, although for any $\Delta$ there is considerable scatter in the times in which particular realisations become unstable.

\subsection{MERCURY6 parameters}

We perform our simulations using the MERCURY6 package \citep{1999MNRAS.304..793C}. We use the Bulirsch-Stoer integrator and a Jacobian coordinate system in the output files. Ejection distance is 11~000~AU and fragmentation from collisions is not allowed. We have redefined a close encounter as when two planets are separated by less than one mutual Hill radius, defined as in \eqref{eq:RH}.

\subsection{Simulation parameters}

In all simulations apart from 4f\_long we stop the integration when a close encounter has occurred or 100~Myr is reached. In 4f\_long we continue running until 100~Myr to investigate the future of multi-planet systems in \secref{sec:discussion}. We use output rates of 10, 100 or 1000~years.  Every simulation contains 100 runs with the difference in initial placement on orbits and the inclinations of the orbits. The simulations we carried out are summarised in \tabref{tab:simulations} and described in detail in \secreftwo{sec:HR8799unstable}{sec:oursol}. We provide the initial conditions of our example system in \tabref{tab:ICs}, but we note that the system is so chaotic that a small change in time-step completely alters the lifetime of the system (see \secref{sec:sim2and4alpha}). A more successful method is to create statistically similar, long-lived systems as explained in this section and described in \secref{sec:oursol} (see also \tabref{sec:simulations}).

Throughout the paper we compare our simulations to the observations of \cite{2011ApJ...729..128C} from 8 October 2009, as these observations have low measurement uncertainties and all four planets are observed simultaneously.

\section{HR~8799 -- unstable as observed}\label{sec:HR8799unstable}

Below we describe in more detail the three simulations 1, 2, and 3 seen in \tabref{tab:simulations}. We use the term stability timescale as the timescale for which the systems in a simulation have not had any close encounters. This timescale we find is roughly the same as the time it takes for the systems to lose a planet.

In \figref{fig:nbrp_iroptsep_comb} the number of systems with 4, 3, 2, and 1 planets in simulations 1, 2, and 3 is shown as a function of time. The black line indicates the estimated age of HR~8799. The figure shows that none of the systems in simulations 1, 2, or 3 still have four planets when the age of HR~8799 is reached. The figure also shows that $\sim$80\% of the systems have two planets left and $\sim$20\% have one planet left after 100~Myr.

\begin{figure}
\centering
\includegraphics[width=\hsize]{tikz/nbrp_iroptsepfo.pdf}
\caption{Number of planets left in the systems as a function of time. Dotted line corresponds to the face-on assumption (simulation 1), solid line to the inclination inferred from observations (simulation 2), and dashed line to the largest separation achieved from inclination (simulation 3). See \tabref{tab:simulations} and \secref{sec:HR8799unstable} for descriptions about the simulations. The age estimate of HR~8799 ($30^{+20}_{-10}$~Myr) is indicated with a solid line with the grey region corresponding to the uncertainty.}
\label{fig:nbrp_iroptsep_comb}
\end{figure}

\subsection{Simulation 1: Face-on assumption}\label{sec:sim1}

For simplicity, we begin by assuming the HR~8799 system is seen entirely face-on, i.e. inclination $i = 0\ensuremath{{}^{\circ}}\xspace$ compared to the plane of the sky. We use the observed on-sky star-planet separations of 67.9 (b), 38.0 (c), 24.5 (d), and 14.5 (e) AU as initial star-planet separations. The initial separations between the orbits are then $\Delta = 4.14$ (b-c), 3.02 (c-d) and 3.55 (d-e)~$R_H$.  The dynamics of the face-on assumption has been tested earlier by \cite{2010ApJ...710.1408F} in the three-planet case and \cite{2010Natur.468.1080M} with all four planets. We confirm their results by finding a median time at which the planetary systems have the first close encounter of $0.10$~Myr and there is no planetary system without close encounters for longer than $9.4$~Myr. We refer to this simulation as simulation 1, see \tabref{tab:simulations} and dotted lines in \figref{fig:nbrp_iroptsep_comb}. Although none of our systems live for the age of HR~8799, we point out here the large range of lifetimes with the longest-lived system surviving 100 times longer than the median.

\subsection{Simulation 2: Inclined system according to observations}\label{sec:sim2}

\begin{figure}
\centering
\includegraphics[scale=1]{tikz/Delta_min}
\caption{Smallest orbit separation, $\Delta _{\min}$, between any planet pair changes with inclination and position angle. We indicate the inclination and position angle of the HR~8799 system measured by \cite{2014ApJ...780...97M} with a white square (simulation~2), which is $i = 26 \pm 3\ensuremath{{}^{\circ}}\xspace$ and position angle $64 \pm 3\ensuremath{{}^{\circ}}\xspace$. The largest $\Delta _{\min}$ occurs at $i = 25\ensuremath{{}^{\circ}}\xspace$ and position angle $42\ensuremath{{}^{\circ}}\xspace$, when $\Delta _{\min} = 3.56\; R_H$. This value is indicated with a filled, white circle (simulation~3). This plot is based on the observations from \cite{2008Sci...322.1348M,2010Natur.468.1080M}.}
\label{fig:Deltamin_cont}
\end{figure}

The disk of the HR~8799 system is observed to be inclined by $i = 26\ensuremath{{}^{\circ}}\xspace$ compared to the plane of the sky and  at a position angle of 64\ensuremath{{}^{\circ}}\xspace E of N \citep[$64 \pm 3\ensuremath{{}^{\circ}}\xspace$, best fit; G.~Kennedy, private communication; ][]{2014ApJ...780...97M}. Here we assume that the planets and disk are mutually coplanar. These inclination parameters we use also later in simulations 4, 5, and 6. The observed, on-sky separations between planets and star are then 67.9~(b), 41.9~(c), 25.8~(d), and 14.5~(e) AU, which means that planets c and d are located slightly further away from the star compared to the face-on case (simulation 1). The initial separations between the orbits are then $\Delta = 3.47$~(b-c), 3.31~(c-d), and 3.90~(d-e) $R_H$. We refer to this simulation as simulation~2 (see \tabref{tab:simulations}, solid lines in \figref{fig:nbrp_iroptsep_comb} and white square marker in \figref{fig:Deltamin_cont}).

The median time during which the systems have not had any close encounter is in this case $0.094$~Myr, while there are no systems that still have four planets after $9.39$~Myr. Compared to simulation~1 the stability timescale is similar and it remains too short to explain the current state of HR~8799.

\subsection{Simulation 3: Optimal separation from inclination}\label{sec:sim3}

Assume we did not have any system inclination estimate. Which inclination and position angle would then give the largest planetary orbit spacing and thus the longest stability timescale? In \figref{fig:Deltamin_cont} the minimum separation between any two planetary orbits in the system is given for sets of inclination and position angle assuming circular orbits centred on the observed mean \citep[using the observations of][]{2008Sci...322.1348M, 2010Natur.468.1080M}. The largest $\Delta _{\min}$ is found with an inclination of $i = 25\ensuremath{{}^{\circ}}\xspace$ and a position angle of 42\ensuremath{{}^{\circ}}\xspace . The separation between the orbits are then $\Delta = 3.56\; R_H$ for all planet pairs, which translates into star-planet separations of 68.8 (b), 41.9 (c), 24.8 (d), and 14.7 (e) AU. We use these optimal separation parameters in simulation 3 (see \tabref{tab:simulations}, dashed lines in \figref{fig:nbrp_iroptsep_comb} and white circle in \figref{fig:Deltamin_cont}). The median time the systems have not had any close encounter is $0.57$~Myr and all systems have had close encounters after $6.85$~Myr.
\figref{fig:nbrp_iroptsep_comb} shows that an increase of $0.25\; R_H$ in $\Delta _{\min}$ results in a factor of $\sim$2 longer stability timescale (compare medians of simulation 2 with simulation 3).

Preferred estimates of the HR~8799 system inclination from a dynamical point of view lie around $20-30\ensuremath{{}^{\circ}}\xspace$ \citep{2009A&A...503..247R, 2012ApJ...755...38S}, while astroseismological measurements favour inclinations higher than 40\ensuremath{{}^{\circ}}\xspace assuming that the planetary orbits are aligned with the stellar spin \citep{2011ApJ...728L..20W}.

None of simulations 1, 2 or 3 manage to explain the existence of the HR~8799 system as they simply fall apart before the age of the system is reached.

\section{A stable solution without resonant lock}\label{sec:oursol}

We find long-lived, simulated systems that look like HR~8799 during some time in their evolution. These systems are created with initially wider separation between the orbits than calculated from observations assuming circular orbits. This means that we assume the planets are not moving on completely circular orbits.

Our simulated planetary orbits have a low eccentricity of $e \lesssim 0.05$ from the beginning of the integration (for description on how this is set up see \secref{sec:orbital_parameters}). The eccentricity enables planets on more widely separated orbits to look like HR~8799 during parts of the orbits. A wider separation between the orbits in a system can make the system survive longer than the estimated age of HR~8799 as the stability timescale is strongly dependent on orbital separation \citep{2014prpl.conf..787D}.

In this section we present our simulations~4, 5, and 6 (see \tabref{tab:simulations}), which all contain long-lived systems that look like HR~8799. In \secref{sec:sim4} we describe simulation~4 and show an example of a long-lived system that seems to fit the observations of HR~8799. In \secref{sec:stat} we go into detail on how we determine if a simulated system fits the observations. \secref{sec:sim2and4alpha} presents a way of finding even more well-fitting systems. In \secref{sec:MMR} we analyse the simulations in terms of resonant behaviour. \secref{sec:sim5and6} describes simulations~5 and 6, which consider lower planet masses that are still within the estimated range.

\subsection{Simulation 4: Wider initial orbital separation}\label{sec:sim4}

We test a wider-orbit assumption by simulating systems with the following initial equal separations between the planet orbits: $\Delta = 3.6$, 3.65, 3.7, 3.75, 3.8, 3.85, 3.9, 3.95, and 4~$R_H$. We place planet~e at an initial star-planet separation of 14.3~AU, which then fixes the initial semi-major axes of the other planets. The planet masses are set to the nominal masses $M = $ 5~$M_{\rm Jup}$ (b) and 7~$M_{\rm Jup}$ (cde) \citep{2010Natur.468.1080M}. These simulations are called simulation 4a--i; see \tabref{tab:simulations}. Simulation~4 thus contains 900 runs compared to simulations~1, 2, and 3, which only contain 100 runs each.
It makes sense to talk about simulations 4a--i for the initialisation of the runs, but not in terms of the architecture of the outcome. The measured minimum orbit separation $\Delta _{\min}$ for each system ranges between 3.6 and 4~$R_H$ in almost each of simulations 4a--i, depending on the orbital eccentricities. Therefore we bin the simulations~4 with $\Delta _{\min}$, regardless of the initial set-up. We make nine bins with each 100~runs. These bins contain different ranges of orbit separations ($\Delta _{\min}$) and thus represent slightly different architectures of the planetary system.

\begin{figure}
\centering
\includegraphics[width=\hsize]{tikz/surv_simatoi}
\caption{Number of systems without close encounters shown with time. The different coloured lines correspond to bins of each 100~systems in simulation~4, sorted in measured minimum orbit separations ($\Delta _{\min}$). The legend indicates the bins and their systems' mean value of $\Delta _{\min}$. The black line and the grey zone show the age estimate with uncertainty of HR~8799.}
\label{fig:4ness_Delta}
\end{figure}

\figref{fig:4ness_Delta} shows the time during which the systems in the simulation~4 bins have not had any close encounter. The figure shows the trend of stability timescale with minimum orbit separation, $\Delta _{\min}$. Compared to the systems in \figref{fig:nbrp_iroptsep_comb} a substantial fraction of systems in \figref{fig:4ness_Delta} have not experienced any close encounter at the estimated age of HR~8799 ($\sim$30~Myr). \cite{2008ApJ...686..580C} find that a larger number of mutual Hill radii are needed to keep their planetary systems stable compared to what we find. For various reasons (e.g. our higher planet masses, lower initial eccentricities, and different stability criterion), we consider both results valid and our results not contradicting the results of \cite{2008ApJ...686..580C}. We see that at any given separation there is a wide range of stability times, a fact to which we return later.

\begin{figure*}
\centering
\begin{subfigure}[t]{.45\textwidth}
\centering
\includegraphics[scale=1]{tikz/sys8_ir_or10_ver9}
\caption{System 8 in simulation 2 shown together with the observations of the HR~8799 planetary system. The system seems to fit the observations well, but suffers a close encounter at an age of 1.6~Myr between planet d and e. In the figure we show the positions of the simulated planets every 10 years (black dots).}
\label{fig:sys8_sim2}
\end{subfigure}
\hspace{5mm}
\begin{subfigure}[t]{.45\textwidth}
\centering
\includegraphics[scale=1]{tikz/sys94_4e_or100_ver9}
\caption{System 94 in simulation~4e shown together with the observations of the HR~8799 planetary system. This system seems to fit the observations well until a close encounter occurs between planets d and e at an age of 45.4~Myr. We use this system as an example of a long-lived system without strong resonances that also look like HR~8799. In the figure we show the positions of the simulated planets every 100 years (black dots).}
\label{fig:sys94_sim4e}
\end{subfigure}
\caption{Comparison between a system in simulation~2 and a system in simulation~4e. Both systems seem to fit the observed positions of the planets, but only the system in simulation~4e is stable for longer than the estimated age of HR~8799. This is an on-sky projection of the simulations assuming a distance of 39.4~pc \citep{2007A&A...474..653V} and inclination 26\ensuremath{{}^{\circ}}\xspace\ \citep[position angle 64\ensuremath{{}^{\circ}}\xspace ,][]{2014ApJ...780...97M}. The axis around which the system is inclined is indicated with a dashed line. The position of the star is indicated with an asterisk. }
\label{fig:tyre}
\end{figure*}

Do the surviving systems in simulation~4 look like HR~8799? In \figref{fig:tyre} we compare a system from simulation~2 (left panel, system 8) with a system from simulation~4e (right panel, system 94). Both runs seem to fit the observations, but the system from simulation~4 survives over 45~Myr because of the slightly more widely separated planetary orbits. We want to stress that system~94, simulation~4e is not a special system, but seems normal compared to the other systems in simulation~4. As we explain in \secref{sec:stat} we find 12 systems fitting the observations in simulation~4. If we consider a distance to the HR~8799 system $1\sigma$ away from the mean of the estimate ($39.4+1.1$~pc $= 40.5$~pc) we almost double the number to 22 fitting systems in simulation~4. How we determine if a system fits the observations is discussed in detail in \secref{sec:stat}.

To demonstrate the similarity of our example system, system 94, to the observed HR~8799 system, we show a snapshot of the system at 37.6~Myr in \figref{fig:match}.

\begin{figure}
\centering
\includegraphics[scale=1]{tikz/match}
\caption{System 94, simulation 4e at 37.6~Myr (small black dots) compared to the observations of \cite{2011ApJ...729..128C} (coloured, larger dots). The axis around which the system is inclined is shown with a dashed line \citep[position angle 64\ensuremath{{}^{\circ}}\xspace\ and inclination 26\ensuremath{{}^{\circ}}\xspace ,][]{2014ApJ...780...97M}. The position of the star is indicated with an asterisk.}
\label{fig:match}
\end{figure}

\subsection{Determining whether a system fits the observations}\label{sec:stat}

In this section we describe in detail how we determine how each simulated system fits the observations. We have chosen to compare the simulated systems with the observations taken on 8 October 2009 by \cite{2011ApJ...729..128C} because the planets are observed simultaneously during this observation and the measurement uncertainties are comparably small.

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{Pictures/hist_94_4e_Currie11_all_d}
\caption{Comparison of simulations to observed planet positions for planet e, d, c, and b subsequently. The histograms show the star-planet distance in our system 94, simulation 4e, after having moved all output to the observed position angle of the planet, bootstrapped a data set, and added a normally distributed error to the simulation data. We show, in black, dark grey, and light grey lines, the observed on-sky star-planet distances from \cite{2011ApJ...729..128C} assuming the mean Earth-HR~8799 distance (39.4~pc) and adding 1 and 2$\sigma$ of this distance estimate \citep{2007A&A...474..653V}. This procedure is described in detail in \secref{sec:stat}.}
\label{fig:hist_stat}
\end{figure*}

We start by checking whether the system has a longer stability timescale than the estimated age of HR~8799
\citep[30~Myr,][]{2010Natur.468.1080M}. If that is the case, we proceed with checking whether the system looks like HR~8799. We check whether the architecture of the simulated system is similar to the architecture of HR~8799 in the following way:

\begin{enumerate}
\item For every output, move each planet along its orbit to the observed position angle. This is done because it is highly unlikely that we find a single output where all four planets happen to be located within the uncertainties of the observations. The space the observations take on compared to the orbit is about 0.04~AU/200~AU = $2 \times 10^{-4}$. The chance of then catching a point in time when all four planets are within the observations would be $(2 \times 10^{-4})^{3} = 8 \times 10^{-12}$ for a very well-fitting system (fixing for example b to always be at the observed position angle). To catch one output where all planets fit simultaneously would therefore require $\sim 10^{11}$ outputs  -- a number about three orders of magnitude higher than our output rate. As during the simulation the planets are at some point located at exactly the observed position angle, but we do not catch that moment, it is validated to move the planets along their orbits in this way to really measure how close they come to the observations. \\
\item Draw one million random star-planet separations from the distributions created in step number 1 for each planet. This is carried out to create a new distribution of star-planet separations with a method called bootstrapping with resampling. \\
\item To simulate the uncertainty of the observation, draw a normally distributed error to each star-planet separation in the bootstrapped distribution using the standard deviation of the observations. We show these distributions for the simulated planets in system 94, simulation~4e in \figref{fig:hist_stat}.\\
\item Measure for each planet the quantile for the observed star-planet separation in the created distribution. The quantile is the fraction of outputs residing to the left of the observed star-planet separation. The star-planet separations are indicated with black, dark grey, and light grey lines in \figref{fig:hist_stat} for Earth-HR~8799 distances of 39.4~pc, then $+1$ and $+2\sigma$ on the distance estimate.
\end{enumerate}
The quantiles are the measures for how well the observed planet position and the simulated planet positions fit. Quantiles close to 0.5 would indicate a very good fit. We assess systems with quantiles between $10^{-4}$ and $1 - 10^{-4}$ as systems fitting to the observations. A low limit for the quantiles means that we allow systems that rarely are found with the architecture of the HR~8799 system to also count as fitting systems. Our example system, system 94, simulation~4e actually has a lower number for the quantile for planet~b and does not count as a fitting system.

A visualisation of the assessment procedure can be seen in \figref{fig:hist_stat} where we show, for system 94, simulation~4e, the star-planet distributions plotted together with the star-planet distances inferred from the Earth-HR~8799 distances 39.4, 40.5, and 41.6~pc (corresponding to the mean of the estimate, then one and two standard deviations away from the mean). The quantiles for system 94, simulation~4e are $q_{\rm e} = 0.89$, $q_{\rm d} = 0.89$, $q_{\rm c} = 0.02$ and $q_{\rm b} = 8.7 \times 10^{-5}$ using the mean of the Earth-HR~8799 distance estimate. If we consider the distance one standard deviation away from the estimate (40.5~pc to HR~8799), however, we find the quantiles $q_{\rm e} = 0.98$, $q_{\rm d} = 0.99$, $q_{\rm c} = 0.16$ and $q_{\rm b} = 0.01$. At this distance we consider the system 94, simulation~4e to fit the observations well. To not miss systems such as system 94, simulation~4e in our analysis, we consider two values of the Earth-HR~8799 mean of the estimate distance (39.4~pc) and the distance one standard deviation away from the mean of the estimate (40.5~pc).

\begin{figure}
\centering
\includegraphics[scale=1]{tikz/grouped_frac8799surv_art11_logy_change}
\caption{We show how many of our systems in simulation~4 (top and middle panels) and simulation~5 (bottom panel) fit the observations. We binned the simulations in ranges of minimum orbit separations, $\Delta _{\min}$, 100 systems in each bin. The light blue shows the fraction of the systems in the bin that have a stability timescale longer than the estimated age of HR~8799. The dark blue curve shows the fraction of the systems in the bin that have architectures fitting to that of HR~8799. We determine whether a system fits through the procedure described in \secref{sec:stat}. The grey bars indicate the fraction of systems in each bin that fulfil both criteria. The top panel shows the results from simulation~4 using the mean of the Earth-HR~8799 distance estimate (39.4~pc), the middle panel shows the results from simulation~4 using a one standard deviation longer Earth-HR~8799 distance (40.5~pc), and the bottom panel shows the results from simulation~5, which has slightly lower planet masses, but is still within the estimated range.}
\label{fig:frac8799surv}
\end{figure}

We summarise our findings of simulation~4 in the top two panels of \figref{fig:frac8799surv}. Again we have binned the simulation~4 with measured minimum orbit separation, $\Delta _{\min}$, 100 systems in each bin. This makes in total 9 bins and 900~systems. We plot in light blue the fraction of the systems in each bin that have a stability timescale longer than the estimated age of HR~8799. In dark blue we show the fraction of systems in each bin that fit the observations of the HR~8799 system according to the above described method. The grey bars show the fraction of systems in each bin that fulfil both criteria and thus fit the HR~8799 system. We find for the nominal Earth-HR~8799 distance (39.4~pc), in the best case (where $\Delta _{\min} = 3.64-3.67$) 5 systems of 100 that fit the observations. In the entire range of orbit separations, we find 12 of 900 systems that fit the observations (see \figref{fig:frac8799surv}a). For the Earth-HR~8799 distance of 40.5~pc we find, in the best case (where $\Delta _{\min} = 3.74 - 3.79$) 8 systems of 100 fitting the observations and 22 of 900 when including all orbit separation ranges (see \figref{fig:frac8799surv}b).

In simulation~4 we create many varieties of architectures in the planetary systems that they would fit many ranges of inclinations and position angles. Therefore we are not concerned by the exact values of the inclination and position angle that we apply \citep[from][]{2014ApJ...780...97M}.

\subsection{Simulations 2$\alpha$ and 4$\alpha$: Verification}\label{sec:sim2and4alpha}

\begin{figure}
\centering
\includegraphics[scale=1]{tikz/copies_tce_94_8}
\caption{Comparison between simulations 4$\alpha$ (violet) and 2$\alpha$ (yellow). The lines show the number of systems in the simulation that have not experienced a close encounter, plotted with time. The slightly wider orbits in simulation~4$\alpha$ compared to simulation~2$\alpha$ makes the systems close encounter time shift about an order of magnitude. The black line and the grey zone show the age estimate and uncertainty of HR~8799.}
\label{fig:copsims}
\end{figure}

SSystem~94, simulation~4e (see e.g. \figref{fig:sys94_sim4e}) is not a special system, but there are many evolutions leading to systems that fit the observations of HR~8799. To show this, we simulate 100 copies of system~94, simulation~4e with the only difference being the simulation output rate, which changes the timestep slightly. This very small timestep variation changes the fate of the systems in$\sim 10^6$~years because they are chaotic. The systems have similar architecture, but with different evolution. We call them simulation~4$\alpha$.

Because our simulated systems are very chaotic it is likely that an integration made with our exact initial conditions would evolve differently from that which we present here. A small change in timestep alters the evolution on a timecale of $10^6$~yrs and our interesting systems have a stability timescale much longer than that.
Therefore, we note that our initial conditions in \tabref{tab:ICs} are likely to not reproduce our system 94, simulation 4e. Instead, we present our recipes on how to create similar systems.

We also create simulation~2$\alpha$ from system~8, simulation~2, using the same procedure as when creating simulation~4$\alpha$. \figref{fig:copsims} shows with time the number of systems in simulations~$2\alpha$ and $4\alpha$ that have not experienced any close encounter. The slopes of the curves show that the systems indeed have different evolutions despite their identical initial conditions. The location of the curves in the diagram shows that the stability timescale of simulation~$4\alpha$ is about one order of magnitude longer than the stability timescale of simulation~$2\alpha$. This leads to 26 systems in simulation~$4\alpha$ surviving longer than the estimated age of HR~8799 and 10 of these also matching the observations of HR~8799.

\subsection{Mean motion resonance analysis}\label{sec:MMR}

In this section we investigate whether resonances affect our systems in simulation~4. Stabilisation because the system is deep in a mean-motion resonance is a common explanation to the survival of the HR~8799 system \citep[see e.g.][]{2010ApJ...710.1408F, 2009MNRAS.397L..16G, 2009A&A...503..247R, 2010ApJ...721L.199M, 2011ApJ...729..128C, 2011ApJ...741...55S, 2013A&A...549A..52E, 2011ApJ...741...55S, 2014MNRAS.440.3140G}. In the following we show that we simulate systems that fit the HR~8799 system without being stabilised by strong resonant lock.

We consider resonant angles correlated to the Laplace resonances, 2:1, 4:2:1, and 8:4:2:1 \citep{1999ssd..book.....M}. We also consider angles related to the secular motion. All these angles are calculated using the mean longitude ($\lambda$) and longitude of periastron ($\varpi$), which are calculated from the longitude of ascending node ($\Omega$), argument of periastron ($\omega$), and mean anomaly ($M$) for each planet in the following way:
\begin{eqnarray}
\lambda &=& \Omega + \omega + M  \\
\varpi &=& \Omega + \omega
.\end{eqnarray}
The resonant angles we consider of the 2:1 resonance are
\begin{eqnarray}
\phi _{1, \; 2:1} =& 2\lambda _{\rm outer} - \lambda _{\rm inner} - \varpi _{\rm outer} \label{eq:phi1_21} \\
\phi _{2, \; 2:1} =& 2\lambda _{\rm outer} - \lambda _{\rm inner} - \varpi _{\rm inner} \label{eq:phi2_21}
,\end{eqnarray}
which for the three planet pairs (b-c, c-d and d-e) result in six angles to consider.

For the 4:2:1 resonance we consider the angle
\begin{equation}\label{eq:phi1_421}
\phi _{1, \; 4:2:1} = \lambda _{\rm inner} - 3\lambda _{\rm middle} + 2\lambda _{\rm outer}
,\end{equation}
which results in two angles as there are two sets of three-planet systems; b-c-d and c-d-e.

For the 8:4:2:1 resonance we take into account the four angles as follows:
\begin{eqnarray}
\phi _{1, \; 8:4:2:1} =& \lambda _{\rm e} - 2\lambda _{\rm d} - \lambda _{\rm c} + 2\lambda _{\rm b}  \\ \label{eq:phi1_8421}
\phi _{2, \; 8:4:2:1} =& 2\lambda _{\rm e} - 5\lambda _{\rm d} + \lambda _{\rm c} + 2\lambda _{\rm b} \\
\phi _{3, \; 8:4:2:1} =& \lambda _{\rm e} - \lambda _{\rm d} - 4\lambda _{\rm c} + 4\lambda _{\rm b} \\
\phi _{4, \; 8:4:2:1} =& \lambda _{\rm e} - 4\lambda _{\rm d} + 5\lambda _{\rm c} - 2\lambda _{\rm b}
,\end{eqnarray}
which we now label with the planet names.

We also look at the difference in longitude of periastron, which is an angle related to secular motion. It is defined as
\begin{equation}\label{eq:varphi_sec}
\varphi  = \varpi _{\rm inner} - \varpi _{\rm outer}
,\end{equation}
which leads to three angles from the three planet pairs (b-c, c-d and d-e).

\begin{figure*}
\centering
\includegraphics[scale=1]{tikz/res_94_4e}
\caption{Evolution of from top to bottom the resonant angles $\phi _{2, \; {\rm1c:2d}}$, $\phi _{1, \; {\rm 1b:2c:4d}}$, $\phi _{1, \; {\rm 1b:2c:4d:8e}}$ and $\varpi _{\rm c} - \varpi _{\rm d}$ for system 94, simulation~4e. The black dots show the output from our simulation, sampled every 100~years. In the top panel (evolution of the $\phi _{2, \; {\rm1c:2d}}$) the output clusters around 0~rad, which is a sign of weak resonant behaviour. The system alternates between libration and circulation in this resonant angle, which is a different behaviour than what is seen in systems affected by strong resonance in which all output would be narrowly confined about 0~rad. None of the bottom three panels show signs of resonant libration, unlike in the cases considered by \cite{2009MNRAS.397L..16G, 2014MNRAS.440.3140G}. The red line indicates the estimated age of HR~8799.}
\label{fig:res_angles_94_4e}
\end{figure*}

In \figref{fig:res_angles_94_4e} we show the evolution of four of the above angles in system~94, simulation~4e. We want to stress that system 94, simulation~4e is a system similar to many of our other systems in the resonant behaviour seen in \figref{fig:res_angles_94_4e} also. The evolution of the angles in our panels 2, 3, and 4 in \figref{fig:res_angles_94_4e} are also shown in the mean motion resonance analysis of \cite{2009MNRAS.397L..16G} (their Figs.~4, middle panel, and 6, top panel) and \cite{2014MNRAS.440.3140G} (their Fig.~1, bottom right panel).
The previously presented solutions \citep[see e.g.][]{2009MNRAS.397L..16G, 2010ApJ...710.1408F, 2014MNRAS.440.3140G} all show small-amplitude librations of one or more resonant angles. In contrast, our solution shows at most transitions between large-amplitude libration and circulation. This configuration does not guarantee the indefinite protection from close encounters that small-amplitude librations provide.

\begin{figure*}
\centering
\includegraphics[width=\hsize]{Pictures/pcolor_res_log_stabtime.eps}
\caption{Resonant behaviour of our simulation~4, binned with stability time ($x$-axis) and 100 systems per bin. Colour shows for different resonant angles ($y$-axis) how many systems in the bin are circulating (left panel), librating (right panel) or transitioning between a resonant and non-resonant state (middle panel). We consider the angles in Eqs.~\ref{eq:phi1_21}--\ref{eq:varphi_sec} (see \secref{sec:MMR}). Our stability times in simulation~4 stretch between 0.04~Myr up to stable systems after 100~Myr. We see a resonant behaviour in the 2:1 resonant angles where it is common that the systems transition between libration and circulation. However, there seems to be no dependence on stability time with resonances.}
\label{fig:pcolor_MMR}
\end{figure*}

To show the overall behaviour of our simulation~4 we bin the systems by stability time into 9 bins, each containing 100 systems. For each system we assess whether the resonant angles presented in Eqs.~\ref{eq:phi1_21}--\ref{eq:varphi_sec} are circulating, librating, or falling in and out of resonance.  We visualise the resonant behaviour of our simulation~4 in \figref{fig:pcolor_MMR}. \figref{fig:pcolor_MMR} shows no trend of longer stability time with librating angles (right panel) and no clear trend of stability time angles transitioning between librational and circulatory behaviour. We infer from this that our systems are not being stabilised by tight resonant lock. However, the separatrix-crossing behaviour we observe in some systems is characteristic of chaotic evolution driven by multiple resonances \citep{1979PhR....52..263C}. Indeed, this may offer an explanation as to why there is a long-lived tail of surviving systems with no consistent resonant protection: \cite{2010ApJ...712..819S} describe how some trajectories stick for extended periods at the edge of the chaotic region of phase space.

\subsection{Simulations 5 and 6: Lower planet masses}\label{sec:sim5and6}

Another way to increase the stability timescale of a planetary system is to decrease the planet masses. In terms of $\Delta$ (\eqreftwo{eq:RH}{eq:Delta}), lowering planet masses widens the orbit separations and thus increases the stability timescale.

We use the exact same initial conditions for simulation~4 and lower the planet masses to create simulation~5 and 6 (see \tabref{tab:simulations}). Simulation~5 has planet masses 4 (b) and 6 (cde) $M_{\rm Jup}$ and simulation~6 has planet masses 3 (b) and 5 (cde) $M_{\rm Jup}$. We want to point out that these planet masses are still comfortably within the estimated masses from cooling models \citep[$5\pm 2\;M_{\rm Jup}$ (b) and $7^{+3}_{-2}\;M_{\rm Jup}$ (cde);][]{2008Sci...322.1348M, 2010Natur.468.1080M}.

\begin{figure}
\centering
\includegraphics[scale=1]{tikz/sys24_5d}
\caption{System~24 in simulation 5d (using lower planet masses of 4~$M_{\rm Jup}$ for planet e and 6~$M_{\rm Jup}$ for planets b, c, and d). Black dots are our output data points assuming the inclination and position angle of 26\ensuremath{{}^{\circ}}\xspace and 64\ensuremath{{}^{\circ}}\xspace from \cite{2014ApJ...780...97M}. The coloured dots are the observations of the HR~8799 planetary system.
The system matches the data and does not experience any close encounter during our 100~Myr integration time.}
\label{fig:tyre_245d}
\end{figure}

\begin{figure*}
\centering
\includegraphics[width=\textwidth]{Pictures/hist_24_5d_Currie11_all_d.eps}
\caption{Same as \figref{fig:hist_stat} but for the system 24, simulation 5d (planet masses are lower than nominal, $M = $ 4~$M_{\rm Jup}$(e), and 6~$M_{\rm Jup}$ (bcd)).}
\label{fig:hist_245d}
\end{figure*}

An example of a system that fits from simulation~5 is system 24, simulation~5d seen in \figreftwo{fig:tyre_245d}{fig:hist_245d}. The system does not experience any close encounter during the 100~Myr integration time and we measure the quantiles to be $q_{\rm e}= 0.93$, $q_{\rm d} = 0.98$, $q_{\rm c} = 0.15$ and $q_{\rm b} =0.07$ at the nominal Earth-HR~8799 distance 39.4~pc. We consider this system to be a good fit to the observations of HR~8799. The system is not stabilised by low-amplitude librational behaviour.

Using the technique described in \secref{sec:stat} we find that for simulation~5, in the best case (where $\Delta _{\min} = 3.81 - 3.85$) 18 systems of 100 fit the observations. Looking at a broader range of orbit separations in total 82 out of 900 systems fit the observations (see \figref{fig:frac8799surv}c). Considering these 1~$M_{\rm Jup}$ lower planet masses in simulation~5 the architecture of HR~8799 is becoming common in our systems also after 30~Myr. The boost of fitting systems in simulation~5 is clearly seen in \figref{fig:frac8799surv} when comparing the bottom panel (simulation~5) with the top two panels (simulation~4).

We include simulation~6 to prove the point that lower planet masses implies much longer stability timescale. In simulation 6 we find 64 fitting systems of 100 in the best case (when $\Delta _{\min} = 3.98 - 4.02$) and 241 of 900 when considering all orbit separation ranges.

Considering a slightly larger Earth-HR~8799 distance increased the number of fitting systems with a factor of two \secref{sec:sim4}. Considering slightly lower planet masses boosts the number of fitting systems by a factor of ten.

\section{Evolutionary phases}\label{sec:discussion}

In this section we use simulation~4f\_long (see \tabref{tab:simulations}) to investigate the evolution of wide-orbit planetary systems. Simulation~4f\_long is created in the same way as simulation~4f, but the integrations are not stopped after a close encounter has occurred and therefore always continue up to 100~Myr.

We find three distinct phases of dynamical evolution in our wide-orbit planetary systems. We call the phases A, B, and C and they correspond to the time before the first close encounter between planets (A), the time between the first and last close encounter (B), and the time after the last close encounter (C). We define a close encounter between planets to be whenever two planets enter their mutual Hill radius, see \eqref{eq:RH}. We can use these evolutionary phases to predict what will happen and what has happened to the HR~8799 system and systems similar to it.

\begin{figure*}
\centering
\includegraphics[scale=1]{tikz/hist_ae_ver8}
\caption{Eccentricity and semi-major axes of the planets in each system of simulation~4f\_long at a randomly picked time for each system in each phase. From top to bottom we show the cases of 1-, 2-, 3-, and 4-planet systems and from left to right we show phase A, B, and C. All phase~A systems have low eccentricity and four planets. In phase~B, the systems have three or more planets. Typically, systems in phase~C have two planets,
one at $a\sim 10$~AU and one at $30 \lesssim a \lesssim 1000$~AU. Both phase~B and C show generally higher eccentricities.}
\label{fig:ae}
\end{figure*}

\figref{fig:ae} shows how the eccentricity and semi-major axes of the planetary orbits change between phase~A, B, and C. In the figure we picked one random time in phase~A, B, and C  in each system in simulation~4f\_long, respectively.

Phase~A is characterised by low-eccentricity orbits, not deviating much from the initial set-up. Phase~B and C show an even distribution of eccentricities ($0 < e < 1$) and planets pile-up around $a\sim 10$~AU. All systems in phase B have three planets or more, while the most common number of planets in phase~C is two planets. In the two-planet case of phase~C, one planet typically has a semi-major axis of $a\sim 10$~AU and the other $a\sim 30-1000$~AU. The timescale of phase~A is what we referred to as the stability timescale, which is the time it takes before the first close encounter occurs. This timescale is strongly dependent on orbital separation \citep{2014prpl.conf..787D}. The timescale of phase~B is independent of initial orbit separation and we measure it to vary between $10^4$ to $10^7$~years, peaking at $10^6$~years. These timescales for the clearing of wide-orbit planets are consistent with previous studies \citep{2009ApJ...696.1600V}. This timescale is typically shorter than the duration of phase~A and the lifetime of the observed systems (see \tabref{tab:obs_planets}). Therefore statistically it is unlikely to observe a system in phase B unless the star is very young.

\subsection{Application to observed wide-orbit planetary systems}\label{sec:obs_systems}

\begin{table*}
\caption{Observed wide-orbit planets around stars.}
\label{tab:obs_planets}
\centering
\begin{tabular}{lllcclcl}
\toprule\midrule
System 	& Age & Spectral type & Planets & $r_{\mathrm{proj}}$ & $M_p$ & Predicted phase & References\\
& [Myr] & & & [AU] & [$M_{\mathrm{Jup}}$] \\
\midrule

1RXS J1609 & 5
& M0$\pm 1$V
 & b & $\sim$330
 & $8^{+3}_{-2}$
 & B/C
 & 1, 2\\

AB Pic & $\sim$30 & K2V
& b
& $\sim$260 & $13-14$
& C
& 3\\

$\beta$ Pic & $12^{+8}_{-4}$
& A6V
& b
& $8-15$
& $9\pm 3$
& A/C
& 4, 5\\

Fomalhaut & $440$
& A3V
& b$^{*}$
& $\sim 119$ & $\lesssim 3$
& B/C
& 6, 7\\

GJ 504 & $160^{+350}_{-60}$ & G0V & b & 43.5 & $4.0^{+4.5}_{-1.0}$ & C
& 8\\

HD 106906 & $13 \pm 2$ & F5V & b & $\sim 650$ & $11 \pm 2$ & B/C
& 9, 10\\

HD 95086
& $10 -17$
& A8V
& b & $56.4 \pm 0.7$
& $5 \pm 2$
 & A/B/C & 11, 12\\

HR~8799 & $30^{+20}_{-10}$ & A5V to F0V & b & 68 & $7^{+3}_{-2}$ & A
& 13, 14, 15, 16\\
& & & c & 38 & $7^{+3}_{-2}$ & & 13, 14, 15\\
& & & d & 24 & $7^{+3}_{-2}$ & & 13, 14, 15\\
& & & e & 15 & $5 \pm 2$ & & 14, 15\\

$\kappa$ And & $30^{+20}_{-10}$
& B9IV
& b
& $55 \pm 2$
& $12.8^{+2.0}_{-1.0}$
& A/C
& 17, 18\\

\bottomrule
\end{tabular}
\tablefoot{1 -- \cite{2014ApJ...784...65B}, 2 -- \cite{2010ApJ...719..497L}, 3 -- \cite{2005A&A...438L..29C}, 4 -- \cite{2010Sci...329...57L}, 5 -- \cite{2001ApJ...562L..87Z}, 6 -- \cite{2012ApJ...754L..20M}, 7 -- \cite{2008Sci...322.1345K}, 8 -- \cite{2013ApJ...774...11K}, 9 -- \cite{2014ApJ...780L...4B}, 10 -- \cite{2012ApJ...746..154P}, 11 -- \cite{2013ApJ...772L..15R}, 12 -- \cite{2014A&A...565L...4G}, 13 -- \cite{2008Sci...322.1348M}, 14 -- \cite{2010Natur.468.1080M}, 15 -- \cite{2011ApJ...729..128C}, 16 -- \cite{1999AJ....118.2993G}, 17 -- \cite{2013ApJ...763L..32C}, 18 -- \cite{2011A&A...525A..71W}
\\
$^*$ -- Note news of Fomalhaut~b being dust cloud from collision with asteroid belt \citep{2011MNRAS.412.2137K,2014ApJ...786...70K}.
}
\end{table*}

In this section we predict the evolutionary phase of directly imaged, massive, wide-orbit planetary systems assuming their evolution is similar to that which we find for HR~8799 in simulation~4f\_long. \tabref{tab:obs_planets} gives some observed properties of the systems we include.

\paragraph{1RXS J1609}

1RXS~J1609~b is found at a projected distance of 330~AU and is probably scattered to this distance as a protoplanetary disk of this size is unlikely. Its age is not sufficient to distinguish between phase~B or C, but if there are no other planets present we can predict phase~C.

\paragraph{AB Pic}

AB~Pic~b is also found at a very large distance from the star ($\sim 260$~AU). It is highly likely the planet has been scattered there and also taking the age of $\sim 30$~Myr into
account, we predict phase~C as phase~B generally lasts for shorter time.

\paragraph{$\beta$ Pic}

$\beta$~Pic~b is separated from $\beta$~Pic by $8-15$~AU. The fact that no outer planets have been found indicates either that the system is a one-planet system in phase~C or that the outer planets are not massive enough to be visible and the system is in phase~A or C. However, the low eccentricity of the orbit \citep[$e \sim 0.06$,][]{2014PNAS..11112661M} favours phase~A.

\paragraph{Fomalhaut}

Fomalhaut~b is found at a projected distance of $\sim 119$~AU, which is far out and makes in situ formation less probable. The predicted eccentric orbit necessary to fit the observations would indicate phase~C considering the advanced age of the system ($\sim 440$~Myr).

\paragraph{GJ 504}

The recently imaged planet around the solar-type star GJ~504 \citep{2008Sci...322.1345K} is found at a projected star-planet distance of $43.5$~AU. GJ~504~b could be an outer planet in a phase~C system (consider the age of $\sim 160$~Myr), which means that there might be another planet on a closer orbit.

\paragraph{HD~106906}

HD~106906~b has an on-sky projected separation of $\sim 650$~AU, which makes it clear that the planet is scattered and the system is in phase~B or C as such a large protoplanetary disk is unlikely. It may be that there is also an inner planet in the HD~106906 system as the inner part cannot yet be resolved.

\paragraph{HD~95086}

HD~95086~b has a projected distance to the star of $56.4 \pm 0.7$~AU and could be an outer planet in a phase~B or C configuration, which then indicates possible inner companions. It could also be a phase~A system with an unresolved inner region because of the distance to the star \citep[$90.4\pm 3.4$,][]{2007A&A...474..653V}. The debris disk of HD~95086 is found to be similar to that of HR~8799 \citep{2015ApJ...799..146S}.

\paragraph{HR~8799}

The number of planets in HR~8799 and their low eccentricities make us classify the system as still being in phase~A.

\paragraph{$\kappa$ And}

The $\kappa$~And system seems like a higher mass version of the systems HD~95086 and GJ~504. The inner regions might hide a phase~A system whilst the outer regions might reveal an outer planet indicating phase~C. The system is old enough for us to predict phase~B to be unlikely.

\section{Summary and conclusions}\label{sec:conclusions}

The star HR~8799 is observed to be accompanied by four massive planets on wide orbits. Puzzlingly, the planetary system has been shown to be unstable in a timescale much shorter than its estimated age. A common explanation of its existence has been stabilisation by mean motion resonances \citep{2014MNRAS.440.3140G}.

We have simulated long-lived systems similar to HR~8799 without forcing them into resonance. The initial orbits of our planets are slightly more widely separated than circular orbits that intersect the observed star-planet on-sky separations. Slightly larger separations between the orbits result in a significantly longer stability timescale. The low eccentricity of the orbits can in some cases make a system have the observed architecture of HR~8799. We place the planets with random azimuth angle, assuming initial low eccentricities.

When using the estimated planetary masses (5, 7, 7 and 7~$M_{\rm Jup}$) together with the estimated Earth-HR~8799 distance we rarely find a fitting system; in the best case, $\Delta _{\min} = 3.64-3.67,$ we find 5 of 100 fitting systems, and in the broader range of orbit separations we find 12 out of 900 fitting systems. We almost double the number of fitting systems by considering the $1\sigma$ longer distance to HR~8799 of 40.5~pc; the best case now considers orbit separation range $\Delta _{\min} = 3.74-3.79$ and yields 8 of 100 fitting systems: in total we find 22 of 900 fitting systems. It is easier to find a fitting system by considering slightly lower planet masses (4, 6, 6, and 6~$M_{\rm Jup}$) that are still comfortably within the planet mass estimate. That results in the best case ($\Delta _{\min} = 3.98 - 4.02$) 18 of 100 systems fit the observations and in the broader range of orbit separations we find 82 out of 900 systems fitting. The above described numbers are visualised in \figref{fig:frac8799surv}. Examples of a system that fits the observations can be seen for the nominal planet masses in \figrefthree{fig:sys94_sim4e}{fig:match}{fig:hist_stat} and for the slightly lower planet masses in \figreftwo{fig:tyre_245d}{fig:hist_245d}.

We test our systems for mean motion resonances (see \secref{sec:MMR}) and find that most of them fall sporadically in and out of the 2:1 resonance between various planet pairs, however, we see no trend in long stability timescale in the systems with more resonant behaviour. This means that these resonances do not work in a stabilising way, unlike the resonances treated earlier in the literature for stability of the HR~8799 system.

We conclude that the HR~8799 system does not necessarily have to be caught in a strong resonance. We simulated many systems that at some point in their evolution look like HR~8799 and we want to stress that our example system (system 94, simulation 4e) is not a lucky coincidence. We present a different solution than before as our fitting systems are stable for longer than the lifetime of HR~8799 and do not show signs that they are deep in resonance.

A bonus from our simulations is a prediction of the future evolution of the HR~8799 system. We see typical architectures in different evolutionary phases of the system. Before close encounters between the planets, the orbits have low eccentricity ($e < 0.2$). After the first close encounter has occurred, the first planet ejection follows within $\sim 10^6$ years. The close encounter and ejection change the architecture of the system and the eccentricities are now more evenly distributed ($0 < e < 1$). In $\sim80$\% of the cases, the systems eventually lose another planet and end up with two planets on eccentric orbits: one with a semimajor axis of $a_{\text{inner}}\sim 10$~AU and the other at $a_{\text{outer}}\sim 30-1000$~AU  (see \figref{fig:ae}). This configuration is comparably more similar to other massive, wide-orbit planets, such as AB~Pic~b or HD~106906~b, than the current HR~8799~bcde.

Our HR~8799 analogues naturally evolve into architectures that resemble other wide-orbit planetary systems. Therefore we find it most plausible that the HR~8799 system is not of a different kind than other observed wide-orbit planetary systems, but simply in an earlier evolutionary phase. This phase seems to be relatively short for massive, wide-orbit planetary systems and it thus makes sense that only one system has been observed in this phase, that is HR~8799.

\begin{acknowledgements}
The authors thank the anonymous referees for their comments which improved the manuscript.
The authors thank Phil Uttley for expertise in statistics. YG acknowledges the support of a PhD studentship at the University of Amsterdam under the supervision of Selma E.\ de Mink. AJ and YG were partially funded by the European Research Council under ERC Starting Grant agreement 278675--PEBBLE2PLANET and by the Swedish Research Council (grant 2010–3710). MBD was supported by the Swedish Research Council (grant 2011-3991). RPC was supported by the Swedish Research Council (grants 2012-2254 and 2012-5807). Calculations presented in this paper were carried out at LUNARC using computer hardware funded in part by the Royal Fysiographic Society of Lund. AM was supported by the Swedish Research Council (grant 2011-3991) and KAW 2012.0150 from the Knut and Alice Wallenberg foundation.
\end{acknowledgements}

\bibliographystyle{aa-package/bibtex/aa}