$C_{r,s}\coloneqq(C_{k}^{(2s)})^{1/2}d^{s/2}e^{2d\gamma_{d}(\ell+1)}\gamma_{d}^{s}(\ell+1)^{s}/\sqrt{(2s)!}$.
∎

####  B.5 Proof of Theorem 4

###### Proof of Theorem 4.

The assumption on the growth of the derivatives of $k_{e}$ implies that
$C_{k}^{(2s)}\leq C_{k}^{2s}$ for the assumed constant $C_{k}$, and we employ
this bound throughout. Now, to more easily track the $s$-dependent constants
in the bound of Theorem 2 we employ a simpler upper bound

$\displaystyle C_{r,s}h^{s}\rho_{X_{n},\mathcal{X}}^{s}$
$\displaystyle=\frac{C_{k}^{s}d^{s/2}e^{2d\gamma_{d}(r+2s)}\gamma_{d}^{s}(r+2s)^{s}}{\sqrt{(2s)!}}h^{s}\rho_{X_{n},\mathcal{X}}^{s}\leq
e^{2d\gamma_{d}r}\left(\frac{C_{k}d^{1/2}e^{4d\gamma_{d}+1}h}{2s}\right)^{s}$
(31)

which holds whenever $\rho_{X_{n},\mathcal{X}}\leq 1/(\gamma_{d}(r+2s))$, and
where we used the elementary fact $1/\sqrt{(2s)!}\leq(e/(2s))^{s}$ to obtain
this simpler bound.

The idea of this proof is to pick a particular value of $s\in\mathbb{N}_{0}$
that will be $\rho_{X_{n},\mathcal{X}}$-dependent. For (31) to hold we require
$\rho_{X_{n},\mathcal{X}}\leq 1/(\gamma_{d}(r+2s))$, and the largest such $s$
for which this requirement is satisfied is

$\displaystyle
s^{\star}\coloneqq\left\lfloor-\frac{r}{2}+\frac{1}{2\gamma_{d}\rho_{X_{n},\mathcal{X}}}\right\rfloor.$
(32)

The assumption $\rho_{X_{n},\mathcal{X}}\leq 1/(2\gamma_{d}(r+1))$ implies
$s^{\star}\geq 1/(4\gamma_{d}\rho_{X_{n},\mathcal{X}})\geq 0$, so $s$ is well-
defined as an element of $\mathbb{N}_{0}$. Further, the assumption
$\rho_{X_{n},\mathcal{X}}\leq 1/(2d^{1/2}\gamma_{d}e^{4d\gamma_{d}+1})$
implies $t\mapsto[(\frac{1}{2}d^{1/2}e^{4d\gamma_{d}+1}h)/t]^{t}$ is
decreasing on $t\in[1/(4\gamma_{d}\rho_{X_{n},\mathcal{X}}),\infty)$ for any
$h\in(0,1]$. Thus we can replace $s$ by
$1/(4\gamma_{d}\rho_{X_{n},\mathcal{X}})$ in (31) to obtain the upper bound

$\displaystyle C_{r,s}h^{s}\rho_{X_{n},\mathcal{X}}^{s}$ $\displaystyle\leq
e^{2d\gamma_{d}r}\left(\frac{\frac{1}{2}C_{k}d^{1/2}e^{4d\gamma_{d}+1}h}{1/(4\gamma_{d}\rho_{X_{n},\mathcal{X}})}\right)^{1/(4\gamma_{d}\rho_{X_{n},\mathcal{X}})}=C_{n,r,s}h^{\frac{1}{4\gamma_{d}\rho_{X_{n},\mathcal{X}}}}$

where
$C_{n,r,s}\coloneqq(2C_{k}d^{1/2}\gamma_{d}e^{4d\gamma_{d}+1}\rho_{X_{n},\mathcal{X}})^{1/(4\gamma_{d}\rho_{X_{n},\mathcal{X}})}$
is a $h$-independent constant. The assumptions we made above on the box fill
distance are satisfied when
$\rho_{X_{n},\mathcal{X}}\leq\min\\{1/(2\gamma_{d}(r+1)),1/(2d^{1/2}\gamma_{d}e^{4d\gamma_{d}+1})\\}$.
∎

####  B.6 Verifying the Assumptions for Example 5

The aim of this appendix is to verify that the function
$\psi(x)=\sin(10x)+1_{x>0}x^{s+4}$ satisfies
$\psi(x)=c_{0}+c_{1}x+c_{2}x^{2}+c_{3}(x)x^{3}$ for some
$c_{0},c_{1},c_{2}\in\mathbb{R}$ and $c_{3}\in H^{s+1}(\mathcal{O})$ where
$\mathcal{O}=(-\delta,\delta)$, $\delta>0$, is an open neighbourhood of 0 and
we recall that $H^{s+1}(\mathcal{O})$ is the Sobolev space of $s+1$ times
weakly differentiable functions on $\mathcal{O}$; see Appendix A. It then
follows from Appendix A that $x\mapsto c_{3}(x)$ and $x\mapsto c_{3}(-x)$ are
elements of $\mathcal{H}_{k_{e}}(\mathcal{X})$ whenever $k_{e}$ is either the
Mátern or Wendland kernel with smoothness $s$ and $\mathcal{X}=[0,\delta)$,
since these kernels reproduce (up to an equivalent norm)
$H^{s+1}(\mathcal{X})$ in dimension $d=1$.

From linearity, it suffices to establish this fact separately for
$\psi_{1}(x)\coloneqq\sin(10x)$ and $\psi_{2}(x)\coloneqq 1_{x>0}x^{s+4}$. For
the first term, since the trigonometric functions are real-analytic, we have a
convergent power series $\psi_{1}(x)=\sum_{i=0}^{\infty}\tilde{c}_{i}x^{i}$
for $x\in[-\epsilon,\epsilon]$ for some $\epsilon>0$. Thus
$\psi_{1}(x)=\tilde{c}_{0}+\tilde{c}_{1}x+\tilde{c}_{2}x^{2}+c_{3}(x)x^{3}$
with $c_{3}(x)=\sum_{i=3}^{\infty}\tilde{c}_{i}x^{i-3}$, and our task is to
show that this latter series is convergent in a neighbourhood of $x=0$. To
this end, we restrict attention to $x\in[-\epsilon/2,\epsilon/2]$, and then
use the ratio test

$\frac{|\tilde{c}_{i}x^{i-3}|}{|\tilde{c}_{i}\epsilon^{i}|}=\frac{1}{|x|^{3}}\left|\frac{x}{\epsilon}\right|^{i}\leq\frac{1}{|x|^{3}}\frac{1}{2^{i}}\rightarrow
0\quad\text{as}\quad i\rightarrow\infty$

to deduce absolute convergence of the series, as required. For the second
term, we take $c_{0}=c_{1}=c_{2}=0$, so that $\psi_{2}(x)=c_{3}(x)x^{3}$ where
$c_{3}(x)=1_{x>0}x^{s+1}$, observing that $c_{3}$ has an $(s+1)$-order weak
derivative $c_{3}^{(s+1)}(x)$ taking the value $0$ on $x\leq 0$ and $(s+1)!$
on $x>0$, so that $c_{3}\in H^{s+1}(\mathcal{O})$.

####  B.7 Proof of Proposition 6 and Corollary 7

###### Proof of Proposition 6.

Our assumption implies $(e_{n})_{n\in\mathbb{N}}$ is Cauchy and thus this
sequence converges to a limit, denoted $e_{\infty}\in\mathbb{R}$. Let
$\varphi:[0,\infty)\rightarrow[0,1]$ be a smooth function with $\varphi(0)=0$,
$\varphi=1$ on $[1,\infty)$, and derivatives uniformly bounded. Let

$\displaystyle e(\mathbf{x})\coloneqq\begin{cases}e_{\infty}&\text{if
}\quad\mathbf{x}=\mathbf{0},\\\
e_{n+1}+(e_{n}-e_{n+1})\prod_{i=1}^{d}\varphi\left(\scaleto{\frac{x_{i}-\mathbf{x}_{n+1,i}}{\mathbf{x}_{n,i}-\mathbf{x}_{n+1,i}}}{20pt}\right)&\text{if
}\quad\mathbf{x}\in[\mathbf{0},\mathbf{x}_{n}]\setminus[\mathbf{0},\mathbf{x}_{n}],\\\
e_{1}&\text{if }\quad\mathbf{x}\notin[\mathbf{0},\mathbf{x}_{1}]\end{cases}$

so that $e(\mathbf{x}_{n})=e_{n}$ and $e(\mathbf{0})=e_{\infty}$. Now set
$f(\mathbf{x})=y_{\infty}+b(\mathbf{x})e(\mathbf{x})$. To establish
$|f|_{\mathcal{H}_{k}(\mathcal{X})}<\infty$ it suffices to establish that
$\lim_{\mathbf{x}\rightarrow\mathbf{0}}\partial^{\bm{\beta}}e(\mathbf{x})$
exists and is finite for all $|\bm{\beta}|=p$, since then we will have $e\in
C^{p}(\mathcal{X})\subset\mathcal{H}(k_{e})$. This is indeed the case, since
on $[\mathbf{0},\mathbf{x}_{n}]\setminus[\mathbf{0},\mathbf{x}_{n}]$ we have
that

$\partial^{\bm{\beta}}e(\mathbf{x})=(e_{n}-e_{n+1})\prod_{i=1}^{d}\varphi^{(\beta_{i})}\left(\frac{x_{i}-\mathbf{x}_{n+1,i}}{\mathbf{x}_{n,i}-\mathbf{x}_{n+1,i}}\right)\frac{1}{(\mathbf{x}_{n,i}-\mathbf{x}_{n+1,i})^{\beta_{i}}}.$

where the $\varphi^{(\beta_{i})}$ terms are uniformly bounded and the
remaining terms vanish since

$|e_{n}-e_{n+1}|\prod_{i=1}^{d}\frac{1}{(\mathbf{x}_{n,i}-\mathbf{x}_{n+1,i})^{\beta_{i}}}\leq\frac{|e_{n}-e_{n+1}|}{\min(\mathbf{x}_{n}-\mathbf{x}_{n+1})^{p}}$

where $p=\beta_{1}+\dots+\beta_{d}$, and this final bound was assumed to
vanish. ∎

###### Proof of Corollary 7.

The renormalised error is $e_{n}=C_{1}+C_{2}x_{n}^{p+1}+O(x_{n}^{p+2})$, so
that

$\frac{e_{n}-e_{n+1}}{(x_{n}-x_{n+1})^{p}}=C_{2}\frac{(x_{n}^{p+1}-x_{n+1}^{p+1})}{(x_{n}-x_{n+1})^{p}}+\frac{O(x_{n}^{p+2})+O(x_{n+1}^{p+2})}{(x_{n}-x_{n+1})^{p}}\rightarrow
0.$

Our assumptions on the sequence $(x_{n})_{n\in\mathbb{N}}$ imply that both
terms vanish in the $n\rightarrow\infty$ limit. ∎

####  B.8 Proof of Proposition 9

First we establish the correctness of the algebraic identity (16) given in the
main text. Again, it is useful to recall that $\mathcal{H}_{k}(\mathcal{X})$
is endowed with the semi-norm $|f|_{\mathcal{H}_{k}(\mathcal{X})}$ which can
be computed using either of the equivalent formulations $\|\mathbf{x}\mapsto
f(\mathbf{x})-f(\mathbf{0})\|_{\mathcal{H}_{k_{b}}(\mathcal{X})}$ or
$\|\mathbf{x}\mapsto(f(\mathbf{x})-f(\mathbf{0}))/b(\mathbf{x})\|_{\mathcal{H}_{k_{e}}(\mathcal{X})}$
(Paulsen and Raghupathi, 2016, Theorem 5.16); here the latter formulation will
be used.

######  Proposition 19.

The estimator
$\sigma_{n}^{2}[f]:=\frac{1}{n}|m_{n}[f]|_{\mathcal{H}_{k}(\mathcal{X})}^{2}$
has the explicit form

$\displaystyle\sigma_{n}^{2}[f]=\frac{1}{n}\left[f(X_{n})^{\top}\mathbf{K}_{b}^{-1}f(X_{n})-\frac{(\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}f(X_{n}))^{2}}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right]$

given in (16) of the main text.

###### Proof.

Let $\bar{m}_{n}[f]=m_{n}[f]-m_{n}[f](\mathbf{0})$. From direct calculation,

$\displaystyle|m_{n}[f]|_{\mathcal{H}_{k}(\mathcal{X})}^{2}$
$\displaystyle=\|\bar{m}_{n}[f]\|_{\mathcal{H}_{k_{b}}(\mathcal{X})}^{2}$
$\displaystyle=\left\|\mathbf{k}_{b}(\cdot)^{\top}\mathbf{K}_{b}^{-1}\left\\{f(X_{n})-\left(\frac{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}f(X_{n})}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right)\mathbf{1}\right\\}\right\|_{\mathcal{H}_{k_{b}}(\mathcal{X})}^{2}$
$\displaystyle=\left\\{f(X_{n})-\left(\frac{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}f(X_{n})}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right)\mathbf{1}\right\\}^{\top}\mathbf{K}_{b}^{-1}\left\\{f(X_{n})-\left(\frac{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}f(X_{n})}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right)\mathbf{1}\right\\}$
$\displaystyle=f(X_{n})^{\top}\mathbf{K}_{b}^{-1}f(X_{n})-\frac{(\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}f(X_{n}))^{2}}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}},$

which completes the proof. ∎

Next we present a general result, which forms the crux of the argument in our
proof of Proposition 9:

######  Proposition 20.

Assume that $f\in\mathcal{H}_{k}(\mathcal{X})$. Then

$\displaystyle\frac{|f(\mathbf{0})-m_{n}[f](\mathbf{0})|}{\sqrt{k_{n}[f](\mathbf{0},\mathbf{0})}}$
$\displaystyle\leq\frac{|f|_{\mathcal{H}_{k}(\mathcal{X})}}{\sigma_{n}[f]}.$

###### Proof.

From (27) in the the proof of Proposition 18, we have
$\mathrm{wce}(\mathsf{s})^{2}=\mathbf{w}^{\top}\mathbf{K}_{b}\mathbf{w}$ for
any algorithm $\mathsf{s}[f]=\mathbf{w}^{\top}f(X_{n})$ with
$\mathbf{1}^{\top}\mathbf{w}=1$. In particular, we can consider the algorithm
$s[f]=m_{n}[f](\mathbf{0})$, which has
$\mathbf{w}=\mathbf{K}_{b}^{-1}\mathbf{1}/(\mathbf{1}^{\top}\mathbf{K}_{n}^{-1}\mathbf{1})$,
to see that

$\displaystyle\mathrm{wce}(\mathsf{s})^{2}=\left(\frac{\mathbf{K}_{b}^{-1}\mathbf{1}}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right)^{\top}\mathbf{K}_{b}\left(\frac{\mathbf{K}_{b}^{-1}\mathbf{1}}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right)=\frac{1}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}=\frac{k_{n}[f](\mathbf{0},\mathbf{0})}{\sigma_{n}^{2}[f]}.$

This implies that

$\displaystyle\frac{|f(\mathbf{0})-m_{n}[f](\mathbf{0})|}{\sqrt{k_{n}[f](\mathbf{0},\mathbf{0})}}$
$\displaystyle\leq\frac{\mathrm{wce}(\mathsf{s})|f|_{\mathcal{H}_{k}(\mathcal{X})}}{\sqrt{k_{n}[f](\mathbf{0},\mathbf{0})}}\leq\frac{|f|_{\mathcal{H}_{k}(\mathcal{X})}}{\sigma_{n}[f]},$

from which the result is established. ∎

###### Proof of Proposition 9.

For $f\in\mathcal{H}_{k}(\mathcal{X})$ we have the decomposition
$f(\mathbf{x})=f(\mathbf{0})+\bar{f}(\mathbf{x})$ with
$\bar{f}\in\mathcal{H}_{k_{b}}(\mathcal{X})$, and
$|f|_{\mathcal{H}_{k}(\mathcal{X})}=\|\bar{f}\|_{\mathcal{H}_{k_{b}}(\mathcal{X})}$.
From Proposition 20 and
$\sigma_{n}^{2}[f]=\frac{1}{n}|m_{n}[f]|_{\mathcal{H}_{k}(\mathcal{X})}^{2}$,
it follows that

$\frac{\lvert
f(\mathbf{0})-m_{n}^{h}[f](\mathbf{0})\rvert}{\sqrt{k_{n}^{h}[f](\mathbf{0},\mathbf{0})}}\leq\frac{n\lvert
f\rvert_{\mathcal{H}_{k}(\mathcal{X})}}{\lvert
m_{n}^{h}[f]\rvert_{\mathcal{H}_{k}(\mathcal{X})}}=\frac{n\lVert\bar{f}\rVert_{\mathcal{H}_{k_{b}}(\mathcal{X})}}{\lVert\bar{m}_{n}^{h}[f]\rVert_{\mathcal{H}_{k_{b}}(\mathcal{X})}},$

where similarly we have decomposed
$m_{n}^{h}[f](\mathbf{x})=m_{n}^{h}[f](\mathbf{0})+\bar{m}_{n}^{h}[f](\mathbf{x})$.
Since $\bar{f}$ and $\bar{m}_{n}^{h}[f]$ are elements of
$\mathcal{H}_{k_{b}}(\mathcal{X})$ and
$k_{b}(\mathbf{x},\mathbf{x}^{\prime})=b(\mathbf{x})b(\mathbf{x}^{\prime})k_{e}(\mathbf{x},\mathbf{x}^{\prime})$,
we may write

$\bar{f}(\mathbf{x})=b(\mathbf{x})e(\mathbf{x})\quad\text{ and
}\quad\bar{m}_{n}^{h}[f](\mathbf{x})=b(\mathbf{x})\tilde{e}_{h}(\mathbf{x})$

for certain $e,\tilde{e}_{h}\in\mathcal{H}_{k_{e}}(\mathcal{X})$, with
$\|\bar{f}\|_{\mathcal{H}_{k_{b}}(\mathcal{X})}=\|e\|_{\mathcal{H}_{k_{e}}(\mathcal{X})}$
and
$\|\bar{m}_{n}^{h}[f]\|_{\mathcal{H}_{k_{b}}(\mathcal{X})}=\|\tilde{e}_{h}\|_{\mathcal{H}_{k_{e}}(\mathcal{X})}$.
Therefore

$\frac{\lvert
f(\mathbf{0})-m_{n}^{h}[f](\mathbf{0})\rvert}{\sqrt{k_{n}^{h}[f](\mathbf{0},\mathbf{0})}}\leq\frac{n\lVert
e\lVert_{\mathcal{H}_{k_{e}}(\mathcal{X})}}{\lVert\tilde{e}_{h}\rVert_{\mathcal{H}_{k_{e}}(\mathcal{X})}}.$
(33)

By construction, $m_{n}^{h}[f](h\mathbf{x}_{i})=f(h\mathbf{x}_{i})$ for
$i=1,\ldots,n$. Because
$f(\mathbf{x})=f(\mathbf{0})+b(\mathbf{x})e(\mathbf{x})$, the function
$\tilde{e}_{h}$ satisfies

$\tilde{e}_{h}(h\mathbf{x}_{i})=e(h\mathbf{x}_{i})+b(h\mathbf{x}_{i})^{-1}[f(\mathbf{0})-m_{n}^{h}[f](\mathbf{0})]$

for each $i=1,\ldots,n$. The assumptions of Theorem 2 with $s=1$ are
satisfied, and thus

$b(h\mathbf{x}_{i})^{-1}\lvert
f(\mathbf{0})-m_{n}^{h}[f](\mathbf{0})\rvert\leq
C_{r,1}h\rho_{X_{n},\mathcal{X}}\frac{\|b\|_{L^{\infty}(\mathcal{X}_{h})}}{b(h\mathbf{x}_{i})}|f|_{\mathcal{H}_{k}(\mathcal{X})}\to
0$

as $h\to 0$ since, $b$ being a polynomial, $\lim_{h\to
0}\|b\|_{L^{\infty}(\mathcal{X}_{h})}b(h\mathbf{x}_{i})^{-1}<\infty$. The
continuity of $e$, which follows from $e\in\mathcal{H}_{k_{e}}(\mathcal{X})$
and the continuity of $k_{e}$, then implies that
$\tilde{e}_{h}(h\mathbf{x}_{i})\to e(\mathbf{0})$ as $h\to 0$. Let
$c\in\mathbb{R}$. The function
$g(\mathbf{x}):=c\,k_{e}(\mathbf{x},\mathbf{x}^{\prime})k_{e}(\mathbf{x}^{\prime},\mathbf{x}^{\prime})^{-1}$
has minimal norm among all functions in $\mathcal{H}_{k_{e}}(\mathcal{X})$
that equal $c$ at $\mathbf{x}^{\prime}$ (Paulsen and Raghupathi, 2016,
Corollary 3.5). By the reproducing property $\langle
k_{e}(\cdot,\mathbf{x}_{1}),k_{e}(\cdot,\mathbf{x}_{2})\rangle_{\mathcal{H}_{k_{e}}(\mathcal{X})}=k_{e}(\mathbf{x}_{1},\mathbf{x}_{2})$
in the RKHS inner product, from which it follows that the norm of this
function is

$\lVert g\rVert_{\mathcal{H}_{k_{e}}(\mathcal{X})}=\langle
g,g\rangle_{\mathcal{H}_{k_{e}}(\mathcal{X})}^{1/2}=\lvert
c\rvert\,k_{e}(\mathbf{x}^{\prime},\mathbf{x}^{\prime})^{1/2}.$

Consequently, setting $\mathbf{x}^{\prime}=h\mathbf{x}_{1}$ and
$c=\tilde{e}_{h}(h\mathbf{x}_{1})$ gives

$\lVert\tilde{e}_{h}\rVert_{\mathcal{H}_{k_{e}}(\mathcal{X})}\geq\lvert\tilde{e}_{h}(h\mathbf{x}_{1})\rvert
k_{e}(h\mathbf{x}_{1},h\mathbf{x}_{1})^{-1/2}\to\lvert e(\mathbf{0})\rvert
k_{e}(\mathbf{0},\mathbf{0})^{-1/2}.$

Using this bound in (33) then shows that the ratio on the left-hand side is
bounded as $h\to 0$, since we have assumed $e(\mathbf{0})\neq 0$. ∎

####  B.9 Proof of Proposition 12

###### Proof of Proposition 12.

Recall that the estimator $\mathbf{r}_{n}^{h}[f]$ is any maximiser of

$\displaystyle\mathcal{L}_{n}^{h}(\mathbf{r})$
$\displaystyle=-\underbrace{\left\\{f(X_{n}^{h})^{\top}\mathbf{K}_{b_{\mathbf{r}},h}^{-1}f(X_{n}^{h})-\frac{(\mathbf{1}^{\top}\mathbf{K}_{b_{\mathbf{r}},h}^{-1}f(X_{n}^{h}))^{2}}{\mathbf{1}^{\top}\mathbf{K}_{b_{\mathbf{r}},h}^{-1}\mathbf{1}}\right\\}}_{=:Q_{n}^{h}(f,\mathbf{r})}-\log\det\mathbf{K}_{b_{\mathbf{r}},h}.$

Maximising $\mathcal{L}_{n}^{h}(\mathbf{r})$ is equivalent to maximising

$\displaystyle\mathcal{J}_{n}^{h}(\mathbf{r})=\mathcal{L}_{n}^{h}(\mathbf{r})-\mathcal{L}_{n}^{h}(\mathbf{r}_{0}).$

Let $\bar{f}(\mathbf{x}):=f(\mathbf{x})-f(\mathbf{0})$. By assumption,
$\bar{f}$ is an element of the RKHS of the covariance function
$k_{b_{\mathbf{r}_{0}}}$. According to Proposition 19,
$Q_{n}^{h}(f,\mathbf{r})=|m_{n}^{h}[f]|_{\mathcal{H}_{k}(\mathcal{X})}^{2}\geq
0$. Also from Proposition 19 and the minimal norm characterisation of the
interpolant, we obtain the $h$-independent bound
$Q_{n}^{h}(f,\mathbf{r}_{0})\leq\lVert\bar{f}\rVert_{\mathcal{H}_{k_{b_{\mathbf{r}_{0}}}}(\mathcal{X})}^{2}$.
Therefore

$\displaystyle\mathcal{J}_{n}^{h}(\mathbf{r})$
$\displaystyle=Q_{n}^{h}(f,\mathbf{r}_{0})+\log\det\mathbf{K}_{b_{\mathbf{r}_{0}},h}-Q_{n}^{h}(f,\mathbf{r})-\log\det\mathbf{K}_{b_{\mathbf{r}},h}$
$\displaystyle\leq\log\det\mathbf{K}_{b_{\mathbf{r}_{0}},h}-\log\det\mathbf{K}_{b_{\mathbf{r}},h}+\lVert\bar{f}\rVert_{\mathcal{H}_{k_{b_{\mathbf{r}_{0}}}}(\mathcal{X})}^{2}.$

Because the covariance matrix factorises as
$\mathbf{K}_{b_{\mathbf{r}},h}=\mathbf{B}_{\mathbf{r},h}\mathbf{K}_{e,h}\mathbf{B}_{\mathbf{r},h}$,
where $\mathbf{K}_{e,h}$ is the covariance matrix for $k_{e}$ at $X_{n}^{h}$
and $\mathbf{B}_{\mathbf{r},h}$ is a diagonal matrix with entries
$b_{\mathbf{r}}(h\mathbf{x}_{i})$, we get

$\displaystyle\mathcal{J}_{n}^{h}(\mathbf{r})\leq
2\sum_{i=1}^{n}\log\frac{b_{\mathbf{r}_{0}}(h\mathbf{x}_{i})}{b_{\mathbf{r}}(h\mathbf{x}_{i})}+\lVert\bar{f}\rVert_{\mathcal{H}_{k_{b_{\mathbf{r}_{0}}}}(\mathcal{X})}^{2}.$

Then, since the class of error bounds $b_{\mathbf{r}}$ is monotonically
parametrised, $\mathcal{J}_{n}^{h}(\mathbf{r})\to-\infty$ as $h\to 0$
uniformly over $\mathbf{r}\in[\mathbf{0},\mathbf{r}_{0}-\bm{\epsilon}]$ for
any $\bm{\epsilon}>\mathbf{0}$. Because
$\mathcal{J}_{n}^{h}(\mathbf{r}_{0})=0$, this establishes that $\liminf_{h\to
0}\mathbf{r}_{n}^{h}[f]\geq\mathbf{r}_{0}$, as claimed. ∎

####  B.10 Calculations for Multidimensional Output

To simplify the presentation we assume throughout that $\sigma=1$, since
$\sigma$ enters only as a multiplicative constant that can be propagated
through the calculations at the end. Let
$k_{\mathcal{X}}(\mathbf{x},\mathbf{x}^{\prime}):=k_{0}^{2}+b(\mathbf{x})b(\mathbf{x}^{\prime})k_{e}(\mathbf{x},\mathbf{x}^{\prime})$.
From the Kronecker decomposition
$\mathbf{K}=\mathbf{K}_{\mathcal{X}}\otimes\mathbf{K}_{\mathcal{T}}$ we have
that

$\displaystyle\mathbf{K}^{-1}=(\mathbf{K}_{\mathcal{X}}^{-1})\otimes(\mathbf{K}_{\mathcal{T}}^{-1})$

where from (24) we know that

$\displaystyle\mathbf{K}_{\mathcal{X}}^{-1}=\mathbf{K}_{b}^{-1}-\mathbf{K}_{b}^{-1}\mathbf{1}(k_{0}^{-2}+\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1})^{-1}\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}.$

Let $\mathbf{k}(\mathbf{x},\mathbf{t})$ be the column vector with entries
$k((\mathbf{x}_{i},\mathbf{t}_{i}),(\mathbf{x},\mathbf{t}))$, and analogously
define $\mathbf{k}_{\mathcal{X}}(\mathbf{x})$ and
$\mathbf{k}_{\mathcal{T}}(\mathbf{t})$ as the column vectors with respective
entries $k_{\mathcal{X}}(\mathbf{x}_{i},\mathbf{x})$ and
$k_{\mathcal{T}}(\mathbf{t}_{i},\mathbf{t})$. In this notation we have also
the Kronecker decomposition
$\mathbf{k}(\mathbf{x},\mathbf{t})=\mathbf{k}_{\mathcal{X}}(\mathbf{x})\otimes\mathbf{k}_{\mathcal{T}}(\mathbf{t})$
and the algebraic result that
$\mathbf{k}(\mathbf{x},\mathbf{t})^{\top}=(\mathbf{k}_{\mathcal{X}}(\mathbf{x})^{\top})\otimes(\mathbf{k}_{\mathcal{T}}(\mathbf{t})^{\top})$.
For the conditional mean function we therefore have that

$\displaystyle m_{n}[f](\mathbf{x},\mathbf{t})$
$\displaystyle=\mathbf{k}(\mathbf{x},\mathbf{t})^{\top}\mathbf{K}^{-1}f(X_{n})$
$\displaystyle=[(\mathbf{k}_{\mathcal{X}}(\mathbf{x})^{\top})\otimes(\mathbf{k}_{\mathcal{T}}(\mathbf{t}))^{\top}][(\mathbf{K}_{\mathcal{X}}^{-1})\otimes(\mathbf{K}_{\mathcal{T}}^{-1})]f(X_{n})$
$\displaystyle=[\mathbf{k}_{\mathcal{X}}(\mathbf{x})^{\top}\mathbf{K}_{\mathcal{X}}^{-1}]\otimes[\mathbf{k}_{\mathcal{T}}(\mathbf{t})^{\top}\mathbf{K}_{\mathcal{T}}^{-1}]f(X_{n})$

where in the final equality we have exploited the _mixed-product_ property of
the Kronecker product. The dependence on $k_{0}^{2}$ of this expression occurs
only in the term
$\mathbf{k}_{\mathcal{X}}(\mathbf{x})^{\top}\mathbf{K}_{\mathcal{X}}^{-1}$.
From the calculations in Section B.2, we can take the
$k_{0}^{2}\rightarrow\infty$ limit of the term
$\mathbf{k}_{\mathcal{X}}(\mathbf{x})^{\top}\mathbf{K}_{\mathcal{X}}^{-1}$, to
obtain that

$\displaystyle m_{n}[f](\mathbf{x},\mathbf{t})$
$\displaystyle=\left\\{\mathbf{k}_{b}(\mathbf{x})^{\top}\mathbf{K}_{b}^{-1}+[1-\mathbf{k}_{b}(\mathbf{x})^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}]\frac{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right\\}\otimes\left[k_{\mathcal{T}}(\mathbf{t})^{\top}\mathbf{K}_{\mathcal{T}}^{-1}\right]f(X_{n}).$

For the conditional covariance function, we have from a similar argument based
on the mixed-product property that

$\displaystyle
k_{n}[f]((\mathbf{x},\mathbf{t}),(\mathbf{x}^{\prime},\mathbf{t}^{\prime}))$
$\displaystyle=k((\mathbf{x},\mathbf{t}),(\mathbf{x}^{\prime},\mathbf{t}^{\prime}))-\mathbf{k}(\mathbf{x},\mathbf{t})^{\top}\mathbf{K}^{-1}\mathbf{k}(\mathbf{x}^{\prime},\mathbf{t}^{\prime})$
$\displaystyle=k_{\mathcal{X}}(\mathbf{x},\mathbf{x}^{\prime})k_{\mathcal{T}}(\mathbf{t},\mathbf{t}^{\prime})-\underbrace{[\mathbf{k}_{\mathcal{X}}(\mathbf{x})^{\top}\mathbf{K}_{\mathcal{X}}^{-1}\mathbf{k}_{\mathcal{X}}(\mathbf{x}^{\prime})]}_{(*)}[\mathbf{k}_{\mathcal{T}}(\mathbf{t})^{\top}\mathbf{K}_{\mathcal{T}}^{-1}\mathbf{k}_{\mathcal{T}}(\mathbf{t}^{\prime})].$
(34)

The term $(*)$ can be read off from (23) and (25):

$\displaystyle(*)=(k_{0}^{2}\mathbf{1}+\mathbf{k}_{b}(\mathbf{x}))^{\top}\left\\{\mathbf{K}_{b}^{-1}-\frac{\mathbf{K}_{b}^{-1}\mathbf{1}\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}}{k_{0}^{-2}+\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}\right\\}(k_{0}^{2}\mathbf{1}+\mathbf{k}_{b}(\mathbf{x}^{\prime}))$

The fractional term can again be treated using the Taylor expansion (26) for
$k_{0}^{-2}$ at 0, which upon expanding brackets yields

$\displaystyle(*)=k_{0}^{2}+\underbrace{\mathbf{k}_{b}(\mathbf{x})^{\top}\mathbf{K}_{b}^{-1}\mathbf{k}_{b}(\mathbf{x}^{\prime})-\frac{[\mathbf{k}_{b}(\mathbf{x})^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}-1][\mathbf{k}_{b}(\mathbf{x}^{\prime})^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}-1]^{\top}}{\mathbf{1}^{\top}\mathbf{K}_{b}^{-1}\mathbf{1}}+O(k_{0}^{-2})}_{(**)}.$

Substituting this expression into (34), and using the definition of
$k_{\mathcal{X}}$, we obtain that

$\displaystyle
k_{n}[f]((\mathbf{x},\mathbf{t}),(\mathbf{x}^{\prime},\mathbf{t}^{\prime}))$
$\displaystyle=(k_{0}^{2}+k_{b}(\mathbf{x},\mathbf{x}^{\prime}))k_{\mathcal{T}}(\mathbf{t},\mathbf{t}^{\prime})-(k_{0}^{2}+(**))[\mathbf{k}_{\mathcal{T}}(\mathbf{t})^{\top}\mathbf{K}_{\mathcal{T}}^{-1}\mathbf{k}_{\mathcal{T}}(\mathbf{t}^{\prime})]$
$\displaystyle=k_{0}^{2}\underbrace{[k_{\mathcal{T}}(\mathbf{t},\mathbf{t}^{\prime})-\mathbf{k}_{\mathcal{T}}(\mathbf{t})^{\top}\mathbf{K}_{\mathcal{T}}^{-1}\mathbf{k}_{\mathcal{T}}(\mathbf{t}^{\prime})]}_{(***)}$
$\displaystyle\qquad+k_{b}(\mathbf{x},\mathbf{x}^{\prime})k_{\mathcal{T}}(\mathbf{t},\mathbf{t}^{\prime})-(**)[\mathbf{k}_{\mathcal{T}}(\mathbf{t})^{\top}\mathbf{K}_{\mathcal{T}}^{-1}\mathbf{k}_{\mathcal{T}}(\mathbf{t}^{\prime})]$

where for $\mathbf{t},\mathbf{t}^{\prime}$ in the training set we have
$(***)=0$, which gives the desired result.

###  Appendix C Details for the Cardiac Model

This appendix summarises the most important aspects of the cardiac model that
we used, and is principally intended for researchers working in the area of
cardiac modelling who are interested in understanding the technical modelling
aspects of the case study that we report in Section 3.

##### Geometry

The model was based on a heart geometry obtained from a computer tomography
(CT) dataset of a single patient, used in Strocchi et al. (2020, 2023). The
original patient geometry was recorded using an average mesh resolution of
1.06 $\pm$ 0.16 mm; below we describe how our coarser and finer meshes were
constructed. Fibres in the atria and the ventricles were generated with
_universal atrial coordinates_ (Roney et al., 2019) and a rule-based method
developed by Bayer et al. (2012). Full details about the geometry are provided
in Strocchi et al. (2020, 2023).

##### Atrial and Ventricular Activation

For the cardiac dynamics, we simulated atrial and ventricular activation using
an _Eikonal model_ (Neic et al., 2017)

$\displaystyle\sqrt{\nabla
t_{a}(\mathbf{x})^{\top}\mathbf{V}(\mathbf{x})\nabla t_{a}(\mathbf{x})}$
$\displaystyle=1\qquad\mathbf{x}\in\Omega\;$ $\displaystyle t_{a}(\mathbf{x})$
$\displaystyle=t_{0}\qquad\mathbf{x}\in\Gamma$

where $t_{a}(\mathbf{x})$ is the local activation time in the active domain
$\Omega$, $\mathbf{V}$ is the tensor of the squared conduction velocities (CV)
in the fibres, sheet and normal directions, and $\Gamma$ is a subset of nodes
in the domain that get activated at time $t_{0}$.

##### Myocardium

The cardiac model that we used treats the myocardium as a transversely
isotropic conduction medium. The fibre CV in the atrial and ventricular
myocardium was set to 0.6 m/s and 0.9 m/s, consistent with experimental
measurements (Taggart et al., 2000; Hansson et al., 1998), while the
anisotropy ratio was set to 0.4 in all tissues (Strocchi et al., 2023). We
simulated fast endocardial conduction in the ventricles by defining a 1 mm
thick layer and assigning it with increased CV by five fold compared to normal
ventricular myocardium (Ono et al., 2009). Fast conduction in the Bachmann
bundle area in the atria was simulated by defining a region between the left
and the right atrium as in Strocchi et al. (2023), and by assigning it with CV
increased by 3.5 fold compared to normal atrial myocardium (Strocchi et al.,
2023). The atria and ventricles were electrically isolated by defining an
insulating layer to avoid non-physiological activation and to control the
atrio-ventricular delay. Atrial and ventricular activation were initiated at
the site of the right atrial and right ventricular pacing lead, respectively,
identified from the CT images. The atrio-ventricular delay was set to 150 ms.

##### Active Tension

Atrial and ventricular activation triggered a rise in local transmembrane
potential (Neic et al., 2017), which then triggered a rise in active tension.
Active tension was modelled with a phenomenological model developed by
Niederer et al. (2011):

$\displaystyle S_{a}(\mathbf{x},t)$
$\displaystyle=T_{\text{peak}}\phi(\lambda)\tanh^{2}\Bigg{(}\frac{t_{s}}{\tau_{r}}\Bigg{)}\tanh^{2}\Bigg{(}\frac{t_{\text{dur}}-t_{s}}{\tau_{d}}\Bigg{)},\qquad
0<t_{s}<t_{\text{dur}}\;,$ $\displaystyle\phi(\lambda)$
$\displaystyle=\tan(\text{ld}(\lambda-\lambda_{0}))\;,\qquad
t_{s}=t-t_{\text{a}}(\mathbf{x})-t_{\text{emd}}\;$

where $T_{\text{peak}}$, $\tau_{r}$, $\tau_{d}$, $t_{\text{dur}}$, ld,
$\lambda_{0}$ and $t_{\text{emd}}$ represent the reference tension, the rise
and the decay time, the twitch duration, the reference stretch and the
electro-mechanical delay. The following parameter values were used:

| $T_{\text{peak}}$ | $\tau_{r}$ | $\tau_{d}$ | $t_{\text{dur}}$ | ld | $\lambda_{0}$ | $t_{\text{emd}}$
---|---|---|---|---|---|---|---
| kPa | ms | ms | ms | - | - | ms
Atria | 80 | 100 | 50 | 450 | 0.7 | 6 | 20
Ventricles | 60 | 50 | 50 | 200 | 0.7 | 6 | 20

##### Passive Properties

We simulated passive properties of atrial and ventricular myocardium with a
transversely isotropic Guccione law

$\displaystyle\Psi(\mathbf{E})$
$\displaystyle=\frac{C}{2}\big{[}e^{Q}-1\big{]}\;,$ $\displaystyle Q$
$\displaystyle=b_{\text{ff}}E_{\text{ff}}^{2}+2b_{\text{fs}}(E_{\text{fs}}^{2}+E_{\text{fn}}^{2})+b_{\text{ss}}(E_{\text{ss}}^{2}+E_{\text{nn}}^{2}+2E_{\text{sn}}^{2})\;,$

where f, s, n in the Cauchy–Green strain tensor $\mathbf{E}$ represent the
strain in the local fibres, sheet and normal to sheet directions, and $C$,
$b_{\text{ff}}$, $b_{\text{fs}}$ $b_{\text{ss}}$ are the bulk stiffness, and
the stiffness in the fibre, cross-fibre and transverse directions. The
following parameter values were used:

| $C$ | $b_{\text{ff}}$ | $b_{\text{fs}}$ | $b_{\text{ss}}$
---|---|---|---|---
| kPa | - | - | -
Atria | 3 | 25 | 11 | 9
Ventricles | 4 | 25 | 11 | 9

Parameters in tissues other than the atria and the ventricles were represented
with a neo-Hookean law and the parameters were set according to Strocchi et
al. (2020, 2023). For all tissues, near-incompressibility was enforced with a
penalty method with a bulk modulus $\kappa=1000.0$ kPa. To constrain the
motion of the heart, we applied omni-directional springs on the superior vena
cava and the two right pulmonary veins, as well as a region around the apex
(Figure 5, left). Spring stiffness was set to 0.1 kPa/$\mu$m.

##### Adjusting the Spatial Resolution

We used meshtool, open-source software for mesh manipulation (Neic et al.,
2020), to refine or coarsen the mesh to a desired spatial resolution for the
experiments that we report. The spatial resolution $x_{1}$ that we used for
extrapolation in the main text is the _nominal_ resolution fed into meshtool,
but we note that the nominal resolution is typically not exactly achieved in
the re-meshing process; we therefore also report the effective average mesh
resolution:

Target dx [mm] |  Effective dx [mm] | # nodes | # elements
---|---|---|---
1.7 | 1.7 | 119,104 | 540,621
1.4 | 1.36 | 213,580 | 1,018,699
1.0 | 1.06 | 417,863 | 1,988,945
0.7 | 0.86 | 764,094 | 3,679,517
0.6 | 0.7 | 1,498,007 | 7,512,728
0.5 | 0.56 | 3,006,080 | 15,768,118
0.4 | 0.43 | 6,217,838 | 30,699,422

##### Computational Resources

All simulations were carried out using the _Cardiac Arrhythmia Research
Package_ (CARP) (Vigmond et al., 2003; Augustin et al., 2021) on ARCHER2, a UK
national super computing service (https://www.archer2.ac.uk/).


# Longterm Stability of Planetary Systems formed from a Transitional Disk

Rory Bowens<EMAIL_ADDRESS>Department of Astronomy & Astrophysics, The
Pennsylvania State University, State College, PA, USA  Center for Exoplanets
and Habitable Worlds, The Pennsylvania State University, State College, PA,
USA  Department of Astronomy, The University of Michigan, Ann Arbor, MI, USA
Andrew Shannon LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne
Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France
Department of Astronomy & Astrophysics, The Pennsylvania State University,
State College, PA, USA  Center for Exoplanets and Habitable Worlds, The
Pennsylvania State University, State College, PA, USA  Rebekah Dawson
Department of Astronomy & Astrophysics, The Pennsylvania State University,
State College, PA, USA  Center for Exoplanets and Habitable Worlds, The
Pennsylvania State University, State College, PA, USA  Jiayin Dong Flatiron
Research Fellow  Department of Astronomy & Astrophysics, The Pennsylvania
State University, State College, PA, USA  Center for Exoplanets and Habitable
Worlds, The Pennsylvania State University, State College, PA, USA  Center for
Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY
10010, USA

###### Abstract

Transitional disks are protoplanetary disks with large and deep central holes
in the gas, possibly carved by young planets. Dong, R., & Dawson, R. 2016,
ApJ, 825, 7 simulated systems with multiple giant planets that were capable of
carving and maintaining such gaps during the disk stage. Here we continue
their simulations by evolving the systems for 10 Gyr after disk dissipation
and compare the resulting system architecture to observed giant planet
properties, such as their orbital eccentricities and resonances. We find that
the simulated systems contain a disproportionately large number of circular
orbits compared to observed giant exoplanets. Large eccentricities are
generated in simulated systems that go unstable, but too few of our systems go
unstable, likely due to our demand that they remain stable during the gas disk
stage to maintain cavities. We also explore whether transitional disk inspired
initial conditions can account for the observed younger ages of 2:1 resonant
systems orbiting mature host stars. Many simulated planet pairs lock into a
2:1 resonance during the gas disk stage, but those that are disrupted tend to
be disrupted early, within the first 10 Myr. Our results suggest that systems
of giant planets capable of carving and maintaining transitional disks are not
the direct predecessors of observed giant planets, either because the
transitional disk cavities have a different origin or another process is
involved, such as convergent migration that pack planets close together at the
end of the transitional disk stage.

Transitional Disks — Protoplanets — Resonance

††software: Numpy

##  1 Introduction

Planets form in disks of gas and dust that surround young stars, known as
protoplanetary disks. Protoplanetary disks are sometimes observed with a deep
and wide gap in the dust distribution (e.g., Strom et al. 1989). These disks
are called transitional disks (Piétu et al., 2005). About $10\%$ of
protoplanetary disks are transitional disks (Luhman et al., 2010), but it
remains unclear whether this fraction indicates that $\sim 10\%$ of
protoplanetary disks spend most of their lives as transitional disks or that
most protoplanetary disks spend $\sim 10\%$ of their lives as transitional
disks (Owen, 2016). Numerous theories have been proposed to explain the
origins of the gaps (e.g., Dullemond & Dominik, 2005; Chiang & Murray-Clay,
2007; Krauss et al., 2007; Suzuki & Inutsuka, 2009; Vorobyov et al., 2015);
the leading theories are photoevaporation (Clarke et al., 2001; Alexander et
al., 2006a, b; Owen et al., 2011, 2012) and planetary sculpting (Calvet et
al., 2005; Dodson-Robinson & Salyk, 2011; Zhu et al., 2011, 2012; Dong et al.,
2015). Photoevaporative winds can deplete the inner disk when the
photoevaporative mass loss rate exceeds the accretion rate. Although early
photoevaporation models (e.g., Owen et al. 2011, 2012) produced
photoevaporative mass loss rates too low to be consistent with the high
observed accretion rates in transitional disks (Ercolano & Pascucci, 2017a),
recent models (e.g., Ercolano et al. 2021; Picogna et al. 2021) that more
accurately model the temperature structure of the disk can meet or exceed half
of observed accretion rates. Transitional disks with higher accretion rates
may be the consequence of stronger photoevaporative winds in carbon-depleted
disks (Ercolano et al., 2018; Wölfer et al., 2019). Alternatively, they may be
the result of magnetically driven supersonic accretion flows (Wang & Goodman
2017; but see Ercolano et al. 2018 for observational arguments against this
hypothesis). More in-depth reviews of the observational properties of
transitional disks and theory of their origins can be found in Espaillat et
al. (2014); Owen (2016); Ercolano & Pascucci (2017b); van der Marel (2017).

In this work, we focus on the question of whether gaps in transitional disks
form by planetary sculpting (Papaloizou & Lin, 1984; Marsh & Mahoney, 1993;
Paardekooper & Mellema, 2004). Planetary sculpting is a process in which
planets or brown dwarfs formed in the protoplanetary disk clear a gap (Dodson-
Robinson & Salyk 2011; Zhu et al. 2011, Rosenthal et al. 2020; see
Paardekooper et al. 2022 for a recent review). This theory is supported by the
discovery of protoplanets forming in the gap of the transitional disk PDS 70
(Keppler et al., 2018; Haffert et al., 2019; Isella et al., 2019; Muley et
al., 2019), possibly LkCa 15 (Kraus & Ireland 2012; Sallum et al. 2015; but
see also Currie et al. 2019), and a few other candidates (HD 100546, HD
142527, HD 169142; Quanz et al., 2013; Biller et al., 2012; Reggiani et al.,
2014). Further supporting this hypothesis, sub-structures in proto-planetary
disks have also been attributed to planetary sculpting (e.g., Huang et al.
2018; Long et al. 2018; Zhang et al. 2018; Choksi & Chiang 2021). Although the
sculpting planets creating the deep and wide gaps in transitional disks are
typically assumed to be Jovian – and we will make that assumption here – Fung
& Chiang (2017) showed that compact configurations of super-Earths can clear
inner cavities in disks with very low viscosity. See Ginzburg & Sari (2018)
and Garrido-Deutelmoser et al. (2022) for more on the properties of such gaps.

To account for the gaps seen in transitional disks, planetary systems must
remain stable over the observed disk timescale (Tamayo et al., 2015). To
assess which planetary system architectures are capable of producing
transitional disks, Dong & Dawson (2016) (hereafter DD16) performed $N$-body
simulations of such planetary systems, using three to six giant planets spaced
from 3 to 30 AU. The depletion of the gaps the planets produced was based on
different assumptions for the disk viscosity and scale height, and these
depletions dictated the spacings of the planets necessary to produce a
continuous gap. DD16 included an eccentricity damping force from the gas in
the gap. They found that a subset of planetary system configurations remained
stable over a typical one million year disk lifetime. Thus they concluded it
was plausible that the subset of planetary systems that contain Jovians
planets in the 3 to 30 AU range all appear as transitional disks during the
protoplanetary stage. Under this hypothesis, the phenomenon of transitional
disks is not pervasive, affecting all protoplanetary disks for $\sim$10% of
their lifetime, but instead is restricted to the $\sim$10% of disks that
happen to harbor giant planets.

Because DD16 conducted simulations only during the gas disk stage, we cannot
directly compare those systems to mature planetary systems observed around
main sequence field stars. After the protoplanetary disk dissipates in
$\lesssim 10^{7}$ years (Haisch et al., 2001; Pfalzner et al., 2014), the
systems may go unstable, as systems of massive planets in such compact
configurations often do (Chambers et al., 1996a; Smith & Lissauer, 2009;
Morrison & Kratter, 2016). However, the damping during the protoplanetary disk
state may allow these systems to find long-term stable compact dynamical
states (e.g., Melita & Woolfson 1996; Lee & Peale 2002; Dawson et al. 2016;
Morrison et al. 2020). The four tightly packed jovian planets around HR 8799
(Marois et al., 2008, 2010) may be in such a dynamical state (Fabrycky &
Murray-Clay, 2010; Goździewski & Migaszewski, 2014; Wang et al., 2018), and
thus perhaps an example of the post-gas evolution of such systems.

Here, we seek to understand the longer-term dynamical evolution of the giant
planet planetary systems capable of sculpting deep and wide gaps in
transitional disks. Simbulan et al. (2017) performed a case study of HL Tau
along these lines, finding ejection of one or more planets to be the most
common outcome. How this result can be generalized remains an open question.
Previous works simulating the longterm evolution of planetary systems compared
simulated systems’ eccentricity and semimajor axis distributions (e.g.,
Chatterjee et al., 2008; Jurić & Tremaine, 2008; Malmberg & Davies, 2009;
Petrovich et al., 2014; Carrera et al., 2019) to those of observed planets
(e.g., Mayor et al., 2011; Dawson & Johnson, 2012; Winn & Fabrycky, 2015; Xie
et al., 2016). However, those works typically began with ad hoc tightly packed
initial configurations of orbits to quickly induce instabilities. Here, we use
initial conditions physically motivated from DD16, as our interest is in how
well those hypothesized transitional-disk sculpting systems match observed
systems. We also explore the prevalence and evolution of mean motion
resonances in these simulated systems, which were common during the stage
simulated by DD16. Motivated by Koriski & Zucker (2011)’s finding that 2:1
mean motion resonance (MMR) systems are younger on average, we will assess
whether the systems that begin in MMR will break within observable timescales.

We describe our simulations in Section 2. We identify the presence and
behavior of orbital resonances Section 3. We assess longterm stability and
which characteristics affect it in Section 4. We investigate the planets’
eccentricities and compare to those of observed exoplanets in Section 5. We
present our conclusions in Section 6.

##  2 Simulations

We simulate the long-term evolution of the multi-planet systems that are
capable of opening the gaps observed in transitional disks. Most of our
simulations start where DD16’s left off, at the end of transitional disk
stage, with gas damping forces turned off, which is an approximation that the
gas disk is immediately removed. Their final planet masses, positions, and
velocities from the gas stage are our initial post-gas conditions. However,
since we do run some additional gas stage simulations, we describe the gas
stage simulations in detail in Appendix A.

For our post-gas simulations use the mercury6 Bulirsch-Stoer integrator with a
1000 AU ejection distance, a solar mass and solar radius central body, an
accuracy parameter of $10^{-12}$, and medium precision outputs every million
years. Our approximation of the instantaneous removal of the gas disk does not
induce problems for several reasons: 1) as we will show the majority of
systems remain stable after the removal, 2) the gas surface density was low to
begin with due to depletion in the gas, and 3) simulations show that
instantaneously removing a depleted gas disk does not significantly alter the
resonant dynamics (Morrison et al., 2020).

The configurations are summarized in Table 1. Configuration names are defined
as planet number followed by planet mass in Jupiters (with letters used for
configurations with the same planet numbers and masses). Being continuations
of DD16, the systems had between 3 and 6 equal mass planets, where the number
of planets required was dictated by the gap widths. We did not include
configurations that DD16 found to be unsuitable for creating proto-planetary
disk cavities because the amount of gas expected to be present in the cavity
would drive the planets too far apart via resonant repulsion to create
overlapping gaps (their 4-10acd, 5-5ac) or would be insufficient to stabilize
the configuration during the gas disk stage (their 6-2). For each
configuration, DD16 used various assumed gas surface densities since ALMA gas
observations and chemical disk modeling are unable to constrain gas depletion
in the inner few AU of a disk. For each configuration, we use the simulations
with the gas surface density consistent with the expected depletion in the
cavity (reported in Dong & Dawson 2016 based on Fung et al. 2014). For their
gas stage simulations, DD16 ran 10 random realizations for each configuration.
We supplement with an additional 10 gas stage realizations for each
configuration, that we then continue in the post-gas stage. See Appendix A for
more details about the gas stage simulations. In summary, we ran eighteen
realizations of twenty variations of the disk initial conditions for 10 Gyr,
for a total of 360 simulations. We found that 14 simulations had gone unstable
during the gas disk stage (i.e., lost a planet via ejections or collisions);
those simulations are considered to have instability events at 0 Gyr in all
future analysis. We also found that two simulations ended the 10 Gyr
simulations with zero planets due to central body collisions.

Table 1: Summary of 3 – 30 AU Systems Name | $\Sigma_{30}$ | $\Delta_{0}$ | GDI | 2:1 MMR % | 10 Gyr
---|---|---|---|---|---
|  g cm-2 |  ($R_{H}$) | Sims. | (Near %) | Stability %
3-5 | 0.001 | 6.2 | 0 | 15, (20) | 100
3-10a | 0.0001 | 5.2 | 0 | 40, (40) | 100
3-10b | 0.1 | 4.8 | 0 | 75, (25) | 100
3-10c | 0.0001 | 4.4 | 0 | 40, (25) | 95
3-10d | 0.1 | 4.4 | 0 | 75, (25) | 100
3-10e | 0.0001 | 4.0 | 1 | 0, (25) | 80
4-2 | 0.001 | 6.1 | 0 | 85, (15) | 95
4-5a | 0.001 | 4.7 | 0 | 85, (15) | 95
4-5b | 0.1 | 4.6 | 0 | 20, (80) | 100
4-5c | 0.1 | 4.6 | 0 | 20, (75) | 100
4-5d | 0.001 | 4.0 | 1 | 45, (45) | 20
4-10b | 0.01 | 3.9 | 2 | 40, (35) | 40
5-1 | 0.1 | 6.2 | 0 | 60, (40) | 90
5-2a | 0.1 | 4.6 | 0 | 95, (5) | 30
5-2b | 0.1 | 4.2 | 7 | 65, (10) | 30
5-5b | 0.01 | 3.6 | 1 | 70, (25) | 40
6-0.5 | 1 | 5.6 | 0 | 50, (40) | 5
6-1 | 1 | 4.4 | 2 | 45, (25) | 5

To better compare giant planets discovered by the radial velocity method –
which are commonly observed at $\sim$ 1–3 AU – we ran an additional set of
simulations that added an interior planet. We refer to the original set as
3–30 AU simulations and additional set as 1–30 AU simulations. The additional
planet is placed interior to the others by the average Hill spacing of
original configuration. First the average Hill spacing of the original
configuration is determined by averaging every pair’s separation in all ~20
simulations for the configuration, excluding those that went unstable during
the gas disk stage. Then the position for the innermost planet is determined
relative to the initial location of the 3 AU planet according to the average
Hill spacing. The newly created 1 AU planet is appended to the corresponding
original 3–30 AU simulation during the gas stage, in order to make the final
results more comparable. Similar to the 3–30 AU gas stage simulations
(Appendix A), these systems are simulated for 1 Myr with gas damping present
using the hybrid integrator (timestep of 3 days, accuracy parameter of
$10^{-12}$). However, they are then simulated for only 1 Gyr with gas damping
shut off using the Bulirsch-Stoer integrator. This reduced simulation time is
necessary for the efficient completion of the simulations.

##  3 Resonant behavior

Here we identify systems containing planets in two-body and three-body
resonances. We assess resonance only for adjacent planet pairs (or triplets in
the case of three body resonances) and by looking for librating resonant
angles. For example, for the 2:1 resonance, the resonant angles are
$\phi_{\rm{in}}$ and $\phi_{\rm{out}}$:

$\phi_{in}=2\lambda_{in}-\lambda_{out}-\varpi_{in}$ (1)

and

$\phi_{out}=2\lambda_{in}-\lambda_{out}-\varpi_{out}$ (2)

where $\lambda_{\rm{in}}$ and $\lambda_{\rm{out}}$ are the mean longitude of
the inner and outer planet, respectively, and $\varpi_{\rm{in}}$ and
$\varpi_{\rm{out}}$ are the longitude of pericenter of the inner and outer
planets, respectively. To ensure the resonant angles are well sampled, we run
simulations with more frequent output (every 10 years) for $10^{5}$ years,
using the same starting point as the standard simulations (i.e., immediately
after gas disk stage). The 3–30 AU systems are run with a Bulirsch-Stoer
integrator ($10^{-12}$ accuracy parameter) while the 1–30 AU are run with a
hybrid integrator (3 day timestep, $10^{-12}$ accuracy parameter).

The libration centers for the massive planets in our simulations are 0 or 180∘
(Dong & Dawson, 2016). We consider a planet pair resonant with libration about
0 if the resonant angle does not come within $\pm$0.5∘ of 180∘ during the
simulation (Fig. 1). Likewise, we consider a planet pair resonant with
libration about 180∘ if the resonant angle does not come within $\pm$0.5∘ of 0
during the simulation.

Furthermore, we label systems with at least one planet pair whose resonant
angle remained within $\pm$170∘ of 0 or 180 for 97.5% of the time as near
resonance. These angles linger near a specific value as they circulate (Fig.
2). If the resonant angles circulate without lingering for all the planet
pairs in a system, we classify the system as non-resonant. We inspect a subset
of 54 of the resonant angle plots by eye to ensure the classification was
reliable and that 0.5∘ worked well as a cut-off with our sampling frequency.
We also considered a stricter resonant angle range, requiring it to stay
$\pm$50∘ of 0∘ or 180 ∘ during the simulation. This did reduce the number of
resonant systems by 12% but had no significant impacts on other results.

Our resonance identification approach can fail for systems that go unstable
quickly. To solve this problem, we truncate the resonance data at 75% of the
instability timescale, defined as the timescale when a collision or ejection
occurs (e.g., for a system that went unstable after four hundred thousand
years, we use the resonant angle during the first three hundred thousand
years). We label all planets that collided or were ejected before 60,000 years
as having undergone "rapid instability." Less then ten systems had rapid
instability so this classification is only a minor part of the results. Some
systems that went unstable during the gas disk stage only had a single planet
remaining at the start of the gas-free stage. These systems are likewise given
a unique label apart from resonance or no resonance.

Figure 1:  The resonance angle for planet pairs in a simulation from Config
3-10b, starting in the post-gas stage. Planets 2 and 3 (the inner two of the
three; panel 3) happened to begin the gas disk stage in resonance and
maintained this configuration in the post-gas stage. In other cases, our
simulated planets get captured into resonance during the gas disk stage.
Figure 2:  The resonance angle among planet pairs in a simulation from Config
3-10a, starting in the post-gas stage. Each resonance angle circulates but we
consider planets 2 and 3 (the inner two of the three; panel 3) to be “near
resonance” because the angle lingers within $\pm 90^{\circ}$ of 0.The pair
reached this configuration part way through the earlier gas stage.

We find that it is more likely for the inner planets to be in 2:1 MMR than the
outer planets. Approximately 50% of the innermost planet pairs are in
resonance in all planet systems. For outermost pairs, 14% of those in the five
planet systems are in 2:1 MMR, 3% of the outer pairs in the six planet
systems, and none of those in the three or four planet systems. The higher
fraction of true resonance among inner pairs vs. outer may be because the
higher surface gas surface density and shorter orbital timescales in the inner
disk facilitate resonance capture during the gas disk stage.

In our four, five, and six planet systems, outer pairs are near 2:1 MMR
(“lingering” systems with at least one planet pair whose resonant angle
remained within $\pm$170∘ of 0 or 180 for 97.5% of the time), in approximately
70% pairs while inner pairs experience closer to 30% with near 2:1 MMR.
Intermediate planet pairs show intermediate rates of near 2:1 MMR. However, in
the three planet systems, 30% of the inner pairs are near resonance but none
of the outer pairs are near resonance.

We also examine the period ratios of 2:1 MMR systems (Fig. 3). Consistent with
DD16’s findings, planets with period ratios far outside of their nominal
resonance can have librating resonant angles. There are some systems that are
librating in the 2:1 that have a period ratio greater than 2.5, even up to 6
(i.e., their 2:1 resonant angle is librating). Libration of the 2:1 resonant
angle at such large period ratios is possible when the longitude of periapse
precesses quickly and is caused in our simulations by the eccentricity damping
during the gas disk stage. This phenomenon, which is more generally
established by a dissipative change to the eccentricity (which may or may not
be accompanied by a change in semi-major axis) is known as resonant repulsion
(e.g., Lithwick & Wu 2012). We find 47% of systems that lie within 10% of a
period ratio of 2 are in 2:1 MMR.

In later sections, we will focus on the 2:1 MMR because other two body
resonant angles rarely librate in our simulations (e.g., 3:2, 3:1, 4:3). Among
the 360 simulations, we find four systems contain one or more pairs that
librate in the 3:2 resonance (twenty-five near resonance), one system in the
3:1 resonance (one-hundred sixty-nine near-resonance), and zero systems in the
4:3 resonance. We also examined various three body resonances. First, we
define the three body librating angle equation:

$\phi_{3b/p,q}(1,2,3)=p\lambda_{1}-(p+q)\lambda_{2}+q\lambda_{3}$ (3)

In the above, $1$, $2$, and $3$ refer to the outer, middle, and inner planet,
respectively. We examined multiple potential resonances and found that only
9:6:4, 15:12:8, 4:2:1, and particularly 3:2:1 had a significant number of
three body resonance cases. 22% of systems displayed at least one type of
three body resonance (1% for three planet systems, 23% for four planet
systems, 45% for five planet systems, and 35% for six planet systems). In all
cases, there were less than 5 near resonance simulations.

Figure 3:  The ratio of the periods of adjacent planets for all 3–30 AU
simulated systems (blue line) and systems in 2:1 MMR (red-dotted/dashed), as
defined by libration of the resonant angle. Systems in 2:1 MMR do not all have
period ratios near 2. Many have ratios much larger than 2, even up to 6.
Similar trends were observed for the 1–30 AU simulations.

##  4 Stability

In approximately 31% of the simulations, an ejection or collision occurs over
10 Gyr. We define a system as stable if none of the member planets were
ejected or collided during either the gas or post-gas stage.

###  4.1 Impact of Planet Number and 2:1 MMR

Our data consists of 360 3–30 AU systems. Each set of 20 systems is defined by
a planet number and initial gas surface density profile (Table 1). In total
there are 6 three planet sets, 6 four planet sets, 4 five planet sets, and 2
six planet sets.

We find that systems with more planets experience more collisions and/or
ejections. Examining systems that remained stable during the gas disk stage,
approximately 13% of the outermost planets are lost (ejected or collided)
during the 10 Gyr gas-free phase in contrast to 20% of inner planets lost. The
difference is mostly due to more collisions among inner planets. In six planet
systems, approximately 33% of planets were ejected and 15% collided. These
values are significantly higher than the averages for the four planet systems
(9% ejected and 3% collided). We conclude that additional planets reduce a
system’s stability but that origin of the ejected planets is fairly uniform
across the semimajor axis range.

We plot the initial fraction of 3–30 AU systems in 2:1 MMR against the
stability fraction, color coded by planet count (Fig. 4). Recall that for a
given planet number, different sets have different planet masses and spacings.
The initial fraction of systems in resonance alone does not appear to impact
the stability but the initial planet count does. Moreover, the fraction of
systems in resonance is not strongly correlated with the initial planet count,
so it appears that planet count is independently driving the stability. The 3
planet systems remain stable ~90% of the time. The four and five planet
systems display a range of stability from low (~20% of a set stable) to high
(~95% of a set stable) values. Both sets of six planet systems have a low
stability fraction (5%).

Figure 4:  The fraction of the 3–30 AU simulations with at least one 2:1 MMR
in a configuration versus the overall stability of the configuration. Each
configuration contains 20 simulations. For readability, the three pairs of
identical sets (3-10b and 3-10d, 4-2 and 4-5a, 4-5b and 4-5c) have had their
2:1 MMR % (75%, 85%, and 20%, respectively) split slightly.

We define the stability timescale as the time until any ejection or collision
occurs within the system. We find some systematic differences in stability
timescale for systems with 2:1 MMR and without (see Figure 5). The timescales
at which systems went unstable differs slightly between the three MMR
categories (resonance, near resonance, and no resonance). However, the overall
rate of instabilities by 10 Gyrs between the three categories was similar: 73%
stable, 68% stable, and 67% stable for resonance, near resonance, and no
resonance; respectively. 10 systems could not have their resonance status
determined due to rapid instability and thus were not used in the above
assessments. From this data, we conclude that 2:1 MMR can improve the
stability of a system in the short term (Myr timescale) but that this
advantage is nullified with time. Based on only small differences in the final
stability rates, we conclude the presence of one or more 2:1 MMR pairs within
a system did not significantly impact the final number of planets that
collided or were ejected of 10 Gyr.

The above analysis of how resonances affect system stability looks at a system
level: if at least one pair is in resonance, we consider it a resonant system.
However, some trends may only be observable at the individual pair level. In
Figure 6, we report stability timescales in a fashion identical to Figure 5
but for adjacent pairs rather than systems. A pair is considered stable until
one member experiences an instability event. The pairs follow a similar trend
to the systems, though we find more significant evidence (compared to the
system-level analysis) that near 2:1 MMR pairs are less stable over the long
term. The higher instability of near systems may be the result of close
spacing without the stabilizing influence of resonance or chaotic behavior at
the resonance separatrix.

We also plot a calculation of resonance pairs based on period ratios rather
than angle libration (Figure 7). Identifying systems near integer period
ratios allows us to more directly compare to observational data where the
orbital parameters are typically not well-constrained enough to determine if
the resonance angle is librating. We determine if a period ratio is resonant
according to the criteria in Koriski & Zucker (2011):

$\delta=2\frac{|r-r_{c}|}{r+r_{c}}\leq 0.1$ (4)

where $r$ is the measured period ratio and $r_{c}$ the resonance period ratio
(in this case 2). The output timesteps are every Myr and once a system leaves
resonance we do not consider it capable of reentering. About 13% of our pairs
begin in resonance post-gas stage according to this criterion, compared to 24%
using the angle libration criterion. Most of these pairs are in five or six
planet configurations, with some sets having up to 40% of pairs near the 2:1.
Although the fraction of pairs near the 2:1 period ratio declines over time,
most disruptions occur in the first 50 Myr. Therefore, this behavior is not a
good explanation for the trend of younger ages for resonant pairs found by
Koriski & Zucker (2011), which requires a typical disruption timescale
$\sim$1–10 Gy and near 100% initial fraction (Dong & Dawson, 2016).

Figure 5:  A CDF comparing instability times based on whether the system
contains resonant, near-resonant, or no resonant planets. There are 185
resonant systems (of which 135 or 73% were stable), 111 near resonant systems
(of which 75 or 68% were stable), and 54 no resonant systems (of which 36 or
67% were stable). Non-resonant systems start below 100% because some
simulations went unstable during the gas disk stage and are marked as unstable
at 0 Gyr. Notably, there is little discrepancy between the final stability
rates, suggesting the presence of 2:1 MMR does not significantly impact
overall stability.  Figure 6:  A CDF comparing instability times for the three
categories of resonance, now plotting per pair rather than per system. There
are 262 resonant pairs (of which 67% were stable), 400 near resonant pairs (of
which 57% were stable), and 411 non- resonant pairs (of which 73% were
stable).  Figure 7:  A CDF showing the fraction of pairs with 2:1 MMR based on
the period ratio. About 13% of pairs begin in resonance according to this
metric.

###  4.2 Impact of Three-Body MMR

We consider the impact of three-body MMR on the stability of systems. These
classifications are once again carried out at the system level. The most
common type of three-body MMR in our systems is 3:2:1, occurring in 11% of
systems. We found that ~51% systems with at least one of the three-body
resonances present had a final stability rate, similar to the rate of 55% of
all systems with four or more planets (our three planet systems have a higher
stability rate but no three body resonances). The other three body MMRs have
small occurrence rates and thus limited statistics, precluding any conclusions
on their significance, if any.

###  4.3 Impact of Mutual Hill radii

We assess the impact of the initial post-gas spacing in mutual Hill radii and
its relationship to stability of the 10 Gyr simulations. We present the
starting (i.e., right after the gas disk phase) and final (i.e., after 10 Gyr)
mutual Hill radii separations for adjacent planet pairs in Figure 8. At the
start of the lifetimes all the planets have roughly equal separations of
mutual Hill radii. Some sets display average mutual Hill radii values that
differ from the majority of sets (see the blue diamond clump at 3,2 and the
purple square clump at 2,1). These pairs tend to remain stable at these wider
separations.

By comparing the two plots, we find that the planets furthest from the star
(i.e., planets on the right) experience the greatest change in mutual Hill
radii separation (especially amongst systems with four or more planets). Those
pairs with Hill radii above 10 (right-side of Figure 8) are exclusively
survivors of instability events (or pairs that began the post-gas stage with
such wide separations, as mentioned in the previous paragraph).

We find that if a system had an individual pair with a initial mutual Hill
radii values smaller than 3.5, it had a greater likelihood of instability.
However, these represented a small fraction of the total pairs. We conclude
that initial mutual Hill radius– within the narrow range encompassed in our
starting conditions – was not a primary driver of the stability rate for a
system, although its final value can be used to assess which systems
experienced instabilities.

Figure 8:  The mutual Hill radii separation of adjacent planet pairs for the
3–30 AU simulations after the gas disk simulations are complete (left) and
after 10 Gyr (right). X position within a column is slightly randomized for
readability. Planet 1 is the farthest from the star. Planets are numbered
upwards so a system with four planets has planets 1, 2, 3, and 4. After the 10
Gyr integration, many systems have far higher hill separations then they
initially had, owing to the tendency of instabilities to greatly alter a
system. If a planet is ejected, the number corresponding to a planet is
updated an new comparisons are made. That is, if planet 1 is ejected, planet 2
is renamed to planet 1 and a pair is then determined between the new planet 1
and planet 0.

###  4.4 1 AU Planets

The simulations by DD16 only included planets in a range from 3 to 30 AU (as
based upon the protoplanetary disk observations, it is unclear if the deep and
wide gaps extend to within 3 AU). For simulation set, we run 10 new
simulations with an additional planet near 1 AU (see Section 2 for details).
The 1 AU planets still fulfill the transitional disk criteria presented by
DD16: i.e., a planet massive enough to carve out a deep gap; packed closely
enough with other planets to create a continuous gap; and close enough to the
star to clear the disk at 1 AU. We run the simulations for 1 Myr the same gas
damping parameters as the corresponding 3–30 AU systems. Then we run the
simulations for 1 Gyr post-gas. The simulation timescale is shorter than for
the 3–30 AU systems to keep the run time feasible. We run 200 simulations in
total, each with 4 to 7 planets.

The increased number of planets reduces the stability of the system. Several 1
– 30 AU sets had stability rates similar to their 3–30 AU counterparts but a
majority of the 1–30 AU sets saw substantial decreases in stability compared
to their 3–30 AU counterparts. The reduction in stability is expected from
past studies. For example, Chambers et al. (1996b) found that adding planets
to a system led to shorter instability timescales. However, the impact was
more modest for systems with higher multiplicity: for example, they found that
a four planet system would go unstable significantly faster than a three
planet system, but a seven planet system would only go unstable a little
faster than a six planet system. We find that the 1–30 AU systems have shorter
instability timescales, particularly for sets that originally contained 3 or 4
planets. The 5 and 6 planet sets were already frequently unstable by 1 Gyr:
thus the additional instability from the extra planet has less impact on the
final stability fraction.

We present a CDF for pair-by-pair resonance status based on period ratio for
the 1 – 30 AU systems (Figure 9; in a fashion identical to the 3–30 AU systems
in Figure 7). The observed trends are mostly identical between the two types
of systems, barring a slightly lower initial resonance fraction for the 1–30
AU systems (approximately 10% compared to the prior 13%). The fraction of 2:1
period pairs at 9 Gyrs is likewise approximately 2% lower for the 1–30 AU
systems. Adding a planet at 1 AU is unable to replicate the trend of younger
ages for resonant pairs found by Koriski & Zucker (2011), which requires a
typical disruption timescale $\sim 1--10$ Gyr and near 100% initial fraction
(Dong & Dawson, 2016). Instead, our systems go unstable usually within $\sim
10--100$ Myr.

Figure 9:  A CDF showing the fraction of pairs with 2:1 MMR based on the
period ratio for 1–30 AU systems. About 10% of pairs begin in resonance
according to this metric, a slight decrease from the same analysis for 3–30 AU
systems (Figure 7). The two types of systems otherwise should identical
trends.

##  5 Eccentricity

The eccentricity distributions of our final systems is a relic of how the
planets evolved with time. In Figure 10, we plot the 1 Gyr eccentricity versus
semimajor axis distribution for the 3–30 AU systems (left) and 1–30 AU systems
(center). Furthermore, in Figure 11, we plot a CDF for the 1 Gyr eccentricity
of several categories. Of particular importance are the "unstable" categories
for both the 3–30 AU and 1–30 AU systems: the unstable category shows the CDF
of 1 Gyr eccentricities for surviving planets in systems that experienced
instability events. Since extending the 3–30 AU systems to 10 Gyr only only
decreased the fraction of circular orbits (<0.05 eccentricity) from 80% to
70%, we show the results at 1 Gyr to compare to the 1–30 AU systems, which
were only simulated for 1 Gyr (Section 4.4).

We studied the influence of ejections and collisions on the eccentricity of
surviving planets. In four and five planet systems (which had 80% and 50%
stability, respectively), the mean eccentricity of all survivors was 0.16
$\pm$ 0.25 while the eccentricity of survivors from unstable systems was 0.47
$\pm$ 0.23. As anticipated instability events increased the eccentricity of
the survivors. Further trends were found in the type of instability events.
Those systems which only experienced collisions had a final mean eccentricity
of 0.33 $\pm$ 0.35 while those survivors that only experienced ejections had a
final mean eccentricity of 0.50 $\pm$ 0.20. Although both types of instability
event correlated with higher eccentricity, systems that experienced ejections
saw a larger increase in the eccentricity of survivors.

The 3–30 AU systems and the 1–30 AU systems exhibit similar trends: nearly all
planets with an eccentricity greater than 0.2 are in systems that experienced
some instability event (primarily ejections, which occurred about 2.5 times
more frequently than collisions) during the simulation. Both sets of
simulations display a high concentration of low eccentricity planets: trending
below 0.1 within 30 AU and below 0.2 beyond 30 AU. Including the additional
planet near $\sim 1au$ slightly increased the typical eccentricities of
surviving planets, though the additional planet itself typically remains at
very low eccentricity if it survives.

We compare the eccentricities between our simulations and observed exoplanets
in both Figure 10 (right) and Figure 11. We use known Jovian mass planets (0.3
to 10 Jupiter M*$\sin(i)$) taken from the Exoplanet Archive on Aug. 4th, 2022.
For our comparison sample, we select from the database those planets with
eccentricities $>0$. Planets discovered by various methods including radial
velocity surveys and direct imaging are included in the comparison sample.
Although there can be some biases in the measured values (Pan et al., 2010),
observers can often measure eccentricity for giant planets detected by radial
velocity, which is how the vast majority of planets suitable for comparison
(i.e., a similar mass and semimajor axis range) to our sample are detected.

The observed planets show eccentricity values with a concentration towards low
eccentricity ($e<0.15$). However, the simulations show a much more extreme
concentration towards very low eccentricity ($e<0.05$). Unstable
configurations can reach high eccentricities (Figure 11), but our
configurations (even if we limit to certain sets) apparently do not produce
the correct mix of stable and unstable systems. Based on the large discrepancy
between the simulations (even when considering only unstable simulations or
only 5 and 6 planet systems) and the observed exoplanet eccentricities, we
conclude that the systems created in DD16, while capable of creating
transitional disks, can not reproduce evolved systems.

Figure 10:  Eccentricity versus semimajor axis for the 3–30 AU systems at 1
Gyr (left), the 1-30 AU systems at 1 Gyr (middle), and a selection of observed
exoplanets (right). Simulated systems are color coded by the initial number of
planets. Most simulated systems remain at low eccentricity $e\lesssim 0.1$ for
the entire evolution, although systems that began with more planets are more
likely to undergo instabilities and end up on elliptical orbits.  Figure 11:
A CDF comparing the final eccentricities of five different system pools. All
3-30 AU results are computed at 1 Gyr to make them directly comparable to the
1-30 AU results. Unstable systems are systems where at least one planet
experienced an instability event.

In contrast, other studies of planet-planet scattering that did not include a
gas disk stage and/or require stability during the gas disk stage (e.g.,
Chatterjee et al. 2008; Jurić & Tremaine 2008) produced eccentricity
distributions that matched the observations well. Our use of equal-mass
planets likely does not account for the difference, as such configurations are
just as likely to lead to elliptical orbits, albeit with a narrow distribution
(Ford et al., 2001). Instead, our requirement that configurations remain
stable during the gas disk stage to maintain a cavity (Dong & Dawson, 2016)
apparently ensures too much stability to excite eccentricities. Our
configurations do not achieve the “dynamically active” state identified by
Jurić & Tremaine (2008) that erases the memory of initial conditions and are
reminiscent of the “dynamically cold” configurations explored by Dawson et al.
(2016) that remain stable after the gas disk stage.

##  6 Conclusions

We simulated the long-term evolution of planetary systems capable of carving
out and maintaining a transitional disk during the gas stage, using initial
conditions from DD16. We subjected the planets to eccentricity damping during
the disk stage and simulated the systems another 10 Gyr post-disk.

Our main finding (Section 5) was that our systems tend to remain on stable,
circular orbits. The typically very low eccentricities are at odds with those
observed in real systems of giant exoplanets discovered via the radial
velocity method. For example, fraction of planets with eccentricities between
0 and 0.2 was approximately 25% higher in our simulated systems compared to
real systems. Among our simulated systems that experience an instability, the
eccentricities can reach the observed high values, but too few of our systems
go unstable. The stability of systems showed a dependence on multiplicity:
three and four planet systems were more stable (i.e., experienced no
collisions or ejections) compared to five and six planet systems by a
significant margin (86% stable versus 33% stable, respectively). This trend is
consistent with prior studies showing that higher planet multiplicity
decreases stability (Chambers et al., 1996b). However, even when considering
only high multiplicity systems and adding planets interior to those needed to
carve observed gaps, our systems were too stable and circular.

We also found that the presence of a 2:1 MMR in a system – which were commonly
established in our simulations during the gas disk stage – did not
significantly impact the overall likelihood of going unstable over 10 Gyr.
Koriski & Zucker (2011) found that the presence of a 2:1 period ratio for a
planetary pair indicates a system is younger on average. However, in contrast
to Koriski & Zucker (2011) who found that the typical lifetime of a 2:1 MMR is
near 4 Gyr, we found that half of our pairs broke their 2:1 period ratio by 10
Myr. We noted a slightly higher instability rate for the near 2:1 MMR systems
(those systems where the resonance angle still circled through all degrees but
with a preference for a certain value) when considering resonance on a pair-
by-pair basis (57% stable compared to approximately 70% stable for the
resonance or no resonance systems).

In future work, we could explore configurations of unequal mass planets,
though past work has shown this is unlikely to significantly boost the
resulting eccentricities (Ford et al., 2001). It may be possible to further
fine-tune our gas disk stage initial conditions (i.e., gas disk properties and
planet spacing) to produce a higher fraction of post-gas unstable systems.
However, more likely stability during the transitional disk stage precludes
wide-spread instabilities after the gas disk disappears. If so, transitional
disks may be typically caused by other processes besides giant planets, such
as photoevaporation (e.g., Picogna et al. 2021) or compact configurations of
super-Earths and mini-Neptunes in low viscosity disk (Fung & Chiang, 2017).
Another possibility is that giant planets undergo convergent widescale
migration at the very end of the transitional disk stage that packs them much
closer together. As proposed by van der Marel & Mulders (2021), this
explanation would also help account for the large size of transitional disks
cavities, compared to the peak in giant planet occurrence at smaller semi-
major axes. However, it would need to be explored what could trigger this just
in-time migration.

Computations for this research were performed on the Pennsylvania State
University’s Institute for Computational & Data Sciences Advanced
CyberInfrastructure (ICS-ACI). This content is solely the responsibility of
the authors and does not necessarily represent the views of the Institute for
CyberScience. The Center for Exoplanets and Habitable Worlds is supported by
the Pennsylvania State University and the Eberly College of Science. This
project was supported in part by NASA XRP NNX16AB50G and NASA XRP
80NSSC18K0355, the National Science Foundation under Grant No. NSF
PHY-1748958, and the Alfred P. Sloan Foundation’s Sloan Research Fellowship.
This research has made use of the NASA Exoplanet Archive, which is operated by
the California Institute of Technology, under contract with the National
Aeronautics and Space Administration under the Exoplanet Exploration Program.
Exoplanet Archive

## References

  * Alexander et al. (2006a) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2006a, MNRAS, 369, 216
  * Alexander et al. (2006b) —. 2006b, MNRAS, 369, 229
  * Biller et al. (2012) Biller, B., Lacour, S., Juhász, A., et al. 2012, ApJ, 753, L38
  * Calvet et al. (2005) Calvet, N., D’Alessio, P., Watson, D. M., et al. 2005, ApJ, 630, L185
  * Carrera et al. (2019) Carrera, D., Raymond, S. N., & Davies, M. B. 2019, A&A, 629, L7
  * Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793
  * Chambers et al. (1996a) Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996a, Icarus, 119, 261
  * Chambers et al. (1996b) —. 1996b, Icarus, 119, 261
  * Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580
  * Chiang & Murray-Clay (2007) Chiang, E., & Murray-Clay, R. 2007, Nature Physics, 3, 604
  * Choksi & Chiang (2021) Choksi, N., & Chiang, E. 2021, Monthly Notices of the Royal Astronomical Society, 510, 1657
  * Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485
  * Currie et al. (2019) Currie, T., Marois, C., Cieza, L., et al. 2019, ApJ, 877, L3
  * Dawson et al. (2015) Dawson, R. I., Chiang, E., & Lee, E. J. 2015, ArXiv:1506.06867, arXiv:1506.06867 [astro-ph.EP]
  * Dawson & Johnson (2012) Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122
  * Dawson et al. (2016) Dawson, R. I., Lee, E. J., & Chiang, E. 2016, ApJ, 822, 54
  * Dodson-Robinson & Salyk (2011) Dodson-Robinson, S. E., & Salyk, C. 2011, ApJ, 738, 131
  * Dong & Dawson (2016) Dong, R., & Dawson, R. 2016, ApJ, 825, 77
  * Dong et al. (2015) Dong, R., Zhu, Z., & Whitney, B. 2015, ApJ, 809, 93
  * Dullemond & Dominik (2005) Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971
  * Ercolano & Pascucci (2017a) Ercolano, B., & Pascucci, I. 2017a, Royal Society Open Science, 4, 170114
  * Ercolano & Pascucci (2017b) —. 2017b, Royal Society Open Science, 4, 170114
  * Ercolano et al. (2021) Ercolano, B., Picogna, G., Monsch, K., Drake, J. J., & Preibisch, T. 2021, MNRAS, 508, 1675
  * Ercolano et al. (2018) Ercolano, B., Weber, M. L., & Owen, J. E. 2018, MNRAS, 473, L64
  * Espaillat et al. (2014) Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 497
  * Fabrycky & Murray-Clay (2010) Fabrycky, D. C., & Murray-Clay, R. A. 2010, ApJ, 710, 1408
  * Ford & Chiang (2007) Ford, E. B., & Chiang, E. I. 2007, ApJ, 661, 602
  * Ford et al. (2001) Ford, E. B., Havlickova, M., & Rasio, F. A. 2001, Icarus, 150, 303
  * Fung & Chiang (2017) Fung, J., & Chiang, E. 2017, ApJ, 839, 100
  * Fung et al. (2014) Fung, J., Shi, J.-M., & Chiang, E. 2014, ApJ, 782, 88
  * Garrido-Deutelmoser et al. (2022) Garrido-Deutelmoser, J., Petrovich, C., Krapp, L., Kratter, K. M., & Dong, R. 2022, ApJ, 932, 41
  * Ginzburg & Sari (2018) Ginzburg, S., & Sari, R. 2018, MNRAS, 479, 1986
  * Goździewski & Migaszewski (2014) Goździewski, K., & Migaszewski, C. 2014, MNRAS, 440, 3140
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 329
  * Haisch et al. (2001) Haisch, Karl E., J., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42
  * Isella et al. (2019) Isella, A., Benisty, M., Teague, R., et al. 2019, arXiv e-prints, arXiv:1906.06308
  * Jurić & Tremaine (2008) Jurić, M., & Tremaine, S. 2008, ApJ, 686, 603
  * Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44
  * Kominami & Ida (2002) Kominami, J., & Ida, S. 2002, Icarus, 157, 43
  * Koriski & Zucker (2011) Koriski, S., & Zucker, S. 2011, ApJ, 741, L23
  * Kraus & Ireland (2012) Kraus, A. L., & Ireland, M. J. 2012, ApJ, 745, 5
  * Krauss et al. (2007) Krauss, O., Wurm, G., Mousis, O., et al. 2007, A&A, 462, 977
  * Lee & Peale (2002) Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596
  * Lithwick & Wu (2012) Lithwick, Y., & Wu, Y. 2012, ApJ, 756, L11
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17
  * Luhman et al. (2010) Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111
  * Malmberg & Davies (2009) Malmberg, D., & Davies, M. B. 2009, MNRAS, 394, L26
  * Marois et al. (2008) Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348
  * Marois et al. (2010) Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080
  * Marsh & Mahoney (1993) Marsh, K. A., & Mahoney, M. J. 1993, ApJ, 405, L71
  * Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv e-prints, arXiv:1109.2497
  * Melita & Woolfson (1996) Melita, M. D., & Woolfson, M. M. 1996, MNRAS, 280, 854
  * Morrison et al. (2020) Morrison, S. J., Dawson, R. I., & MacDonald, M. 2020, ApJ, 904, 157
  * Morrison & Kratter (2016) Morrison, S. J., & Kratter, K. M. 2016, ApJ, 823, 118
  * Muley et al. (2019) Muley, D., Fung, J., & van der Marel, N. 2019, ApJ, 879, L2
  * Owen (2016) Owen, J. E. 2016, PASA, 33, e005
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880
  * Owen et al. (2011) Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13
  * Paardekooper et al. (2022) Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, arXiv e-prints, arXiv:2203.09595
  * Paardekooper & Mellema (2004) Paardekooper, S. J., & Mellema, G. 2004, A&A, 425, L9
  * Pan et al. (2010) Pan, M., Zakamska, N. L., & Ford, E. B. 2010, in EAS Publications Series, Vol. 42, EAS Publications Series, ed. K. Gożdziewski, A. Niedzielski, & J. Schneider, 169
  * Papaloizou & Lin (1984) Papaloizou, J., & Lin, D. N. C. 1984, ApJ, 285, 818
  * Papaloizou & Larwood (2000) Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823
  * Petrovich et al. (2014) Petrovich, C., Tremaine, S., & Rafikov, R. 2014, ApJ, 786, 101
  * Pfalzner et al. (2014) Pfalzner, S., Steinhausen, M., & Menten, K. 2014, ApJ, 793, L34
  * Picogna et al. (2021) Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611
  * Piétu et al. (2005) Piétu, V., Guilloteau, S., & Dutrey, A. 2005, A&A, 443, 945
  * Quanz et al. (2013) Quanz, S. P., Amara, A., Meyer, M. R., et al. 2013, ApJ, 766, L1
  * Reggiani et al. (2014) Reggiani, M., Quanz, S. P., Meyer, M. R., et al. 2014, ApJ, 792, L23
  * Rein (2012) Rein, H. 2012, MNRAS, 422, 3611
  * Rosenthal et al. (2020) Rosenthal, M. M., Chiang, E. I., Ginzburg, S., & Murray-Clay, R. A. 2020, MNRAS, 498, 2054
  * Sallum et al. (2015) Sallum, S., Follette, K. B., Eisner, J. A., et al. 2015, Nature, 527, 342
  * Simbulan et al. (2017) Simbulan, C., Rein, H., Tamayo, D., Petrovich, C., & Murray, N. 2017, Monthly Notices of the Royal Astronomical Society, 469, 3337
  * Smith & Lissauer (2009) Smith, A. W., & Lissauer, J. J. 2009, Icarus, 201, 381
  * Strom et al. (1989) Strom, K. M., Strom, S. E., Edwards, S., Cabrit, S., & Skrutskie, M. F. 1989, AJ, 97, 1451
  * Suzuki & Inutsuka (2009) Suzuki, T. K., & Inutsuka, S.-I. 2009, in American Institute of Physics Conference Series, ed. T. Usuda, M. Tamura, & M. Ishii, Vol. 1158, 161
  * Tamayo et al. (2015) Tamayo, D., Triaud, A. H. M. J., Menou, K., & Rein, H. 2015, The Astrophysical Journal, 805, 100
  * van der Marel (2017) van der Marel, N. 2017, in Astrophysics and Space Science Library, ed. M. Pessah & O. Gressel, Vol. 445, 39
  * van der Marel & Mulders (2021) van der Marel, N., & Mulders, G. D. 2021, AJ, 162, 28
  * Vorobyov et al. (2015) Vorobyov, E. I., Lin, D. N. C., & Guedel, M. 2015, A&A, 573, A5
  * Wang et al. (2018) Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192
  * Wang & Goodman (2017) Wang, L., & Goodman, J. J. 2017, ApJ, 835, 59
  * Winn & Fabrycky (2015) Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409
  * Wölfer et al. (2019) Wölfer, L., Picogna, G., Ercolano, B., & van Dishoeck, E. F. 2019, MNRAS, 490, 5596
  * Xie et al. (2016) Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Science, 113, 11431
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47
  * Zhu et al. (2012) Zhu, Z., Nelson, R. P., Dong, R., Espaillat, C., & Hartmann, L. 2012, ApJ, 755, 6
  * Zhu et al. (2011) Zhu, Z., Nelson, R. P., Hartmann, L., Espaillat, C., & Calvet, N. 2011, ApJ, 729, 47

##  Appendix A Gas Stage Simulations

Some of the gas stage simulations are taking directly from DD16 but, as
described in Section 2, we run supplemental gas stage simulations as part of
this course. Gas stage simulations use the mercury6 Bulirsch-Stoer hybrid
integrator with a timestep of and Bulirsch-Stoer accuracy parameter of
$10^{-12}$) (Chambers, 1999). Each configuration features three to six equal
mass planets, with initial semimajor axes between 2 and 20 AU. The initial
semi-major axes were selected so that the planets just barely open a common
gap in the disk spanning 3–30 AU. Each configuration (Table 1) meets this
condition for an assumed disk viscosity, disk scale height, and planet mass
(0.5–10 Jupiter masses).

For each configuration, we run 20 random realizations; 10 were run as part of
DD16 and 10 are new to this paper. The semi-major was randomized from the
default values listed in Table 2 by about 5%, drawing from a normal
distribution centered on the default value. Initial eccentricities were 0.
Initial inclinations were random set to $\sim 0.01^{\circ}$ to avoid perfectly
coplanar planets. The initial mean longitude, longitude of periapse, and
longitude of ascending node of each planet are randomized between 0 and
$2\pi$.

The effects of the gas within the gap were implemented using prescriptions
from Papaloizou & Larwood (2000), Kominami & Ida (2002), Ford & Chiang (2007),
and Rein (2012) as user-defined forces in mercury6, following Dawson et al.
(2015) and DD16. The timescales for gas damping are:

$\displaystyle\tau=0.029\frac{{\rm g\ cm}^{-2}}{\Sigma_{30}}\left(\frac{a}{\rm
AU}\right)^{2}\frac{M_{\odot}}{M_{\rm p}}{\rm yr}\times$ $\displaystyle 1$
$\displaystyle,v<c_{s}$ $\displaystyle\left(\frac{v}{c_{s}}\right)^{3}$
$\displaystyle,v>c_{s},i<c_{s}/v_{\rm kep}$
$\displaystyle\left(\frac{v}{c_{s}}\right)^{4}$ $\displaystyle,i>c_{s}/v_{\rm
kep}$

where $v=\sqrt{e^{2}+i^{2}}v_{\rm kep}$ and $v_{\rm kep}$ is the Keplerian
velocity, and the sound speed $c_{s}=1.29{\rm km/s}\left(\frac{a}{\rm
AU}\right)^{-1/4}$. We impose $\dot{e}/e=-1/\tau$ and $\dot{i}/i=-2/\tau$
(Kominami & Ida, 2002). The value of $\Sigma_{30}$ for each configuration is
listed in Table 1. Following DD16, we do not impose migration ($\dot{a}$) or
precession ($\dot{\varpi}$, $\dot{\Omega}$) because DD16 found that these
effects were negligible except for fine-tuned values. Furthermore, the
direction and magnitude of migration is sensitive to uncertain disk
conditions; the magnitude is typically small in a depleted cavity. See Section
4.2.3 of DD16 for more detail.

Table 2: 3 – 30 AU Systems Default Semimajor Axis Name |  $P_{1}$ $a$ |  $P_{2}$ $a$ |  $P_{3}$ $a$ |  $P_{4}$ $a$ |  $P_{5}$ $a$ |  $P_{6}$ $a$
---|---|---|---|---|---|---
| (AU) | (AU) | (AU) | (AU) | (AU) | (AU)
3-5 | 2.6 | 7.1 | 19.4 | - | - | -
3-10a | 2.2 | 6.2 | 18.1 | - | - | -
3-10b | 2.6 | 6.9 | 18.5 | - | - | -
3-10c | 3.3 | 8.1 | 19.6 | - | - | -
3-10d | 3.3 | 8.0 | 19.0 | - | - | -
3-10e | 4.3 | 9.4 | 20.5 | - | - | -
4-2 | 2.8 | 5.7 | 11.5 | 23.1 | - | -
4-5a | 2.3 | 4.8 | 10.1 | 21.0 | - | -
4-5b | 2.4 | 4.9 | 10.8 | 20.7 | - | -
4-5c | 2.5 | 5.0 | 10.4 | 21.3 | - | -
4-5d | 3.4 | 6.3 | 11.8 | 22.1 | - | -
4-10b | 2.2 | 4.6 | 9.8 | 20.8 | - | -
5-1 | 2.5 | 4.3 | 7.6 | 13.2 | 23.1 | -
5-2a | 2.9 | 4.8 | 8.2 | 13.8 | 23.3 | -
5-2b | 3.7 | 5.8 | 9.3 | 14.9 | 23.8 | -
5-5b | 2.5 | 4.3 | 7.5 | 13.1 | 22.9 | -
6-0.5 | 3.5 | 5.1 | 7.6 | 11.3 | 16.7 | 24.8
6-1 | 3.6 | 5.3 | 7.8 | 11.4 | 16.8 | 24.8