We list the guarantees of the sketch that we will use in the error analysis,
and point to where they come from in Section 6.2. The guarantees of Theorem
6.5 individually fail with probability $O(\delta)$, so we will rescale to say
that they all hold with probability $\geqslant 1-\delta$.

$\displaystyle\MT@delim@RIGHT\norm@nostar{[AS]_{*,j}}^{2}$
$\displaystyle\leqslant 2\lVert A\rVert_{F}^{2}/s\text{ for all }j\in[s]$ by
Definition 6.4 ($AS$ col bound)
$\displaystyle\MT@delim@RIGHT\norm@nostar{S^{\dagger}b}^{2}$
$\displaystyle\leqslant 2\MT@delim@RIGHT\norm@star{b}^{2}$ by Definition 6.4
($\|S^{\dagger}b\|$ bound)
$\displaystyle\MT@delim@RIGHT\norm@nostar{AA^{\dagger}-AS(AS)^{\dagger}}$
$\displaystyle\leqslant\mu\varepsilon$ by Theorem 6.5 ($AA^{\dagger}$ AMP)
$\displaystyle\MT@delim@RIGHT\norm@nostar{Ab-ASS^{\dagger}b}$
$\displaystyle\leqslant\mu\varepsilon$ by Theorem 6.5 ($Ab$ AMP)
$\displaystyle\lVert TAS\rVert_{F}^{2}$ $\displaystyle=\lVert
AS\rVert_{F}^{2}\leqslant 2\lVert A\rVert_{F}^{2}$ by Definition 6.4
($\|TAS\|_{F}$ bound)
$\displaystyle\MT@delim@RIGHT\norm@nostar{(AS)^{\dagger}(AS)-(TAS)^{\dagger}TAS}$
$\displaystyle\leqslant\mu\varepsilon$ by Theorem 6.5 ($(AS)^{\dagger}AS$ AMP)

A direct consequence of Eq. $AA^{\dagger}$ AMP and Eq. $(AS)^{\dagger}AS$ AMP
is that $AS$ and $TAS$ have bounded spectral norms.

$\displaystyle\MT@delim@RIGHT\norm@nostar{AS}^{2}$
$\displaystyle=\MT@delim@RIGHT\norm@nostar{AS(AS)^{\dagger}}\leqslant\MT@delim@RIGHT\norm@nostar{AA^{\dagger}}+\mu\varepsilon=\MT@delim@RIGHT\norm@star{A}^{2}+\mu\varepsilon$
$\displaystyle\MT@delim@RIGHT\norm@nostar{TAS}^{2}$
$\displaystyle=\MT@delim@RIGHT\norm@nostar{(TAS)^{\dagger}(TAS)}\leqslant\MT@delim@RIGHT\norm@nostar{(AS)^{\dagger}(AS)}+\mu\varepsilon=\MT@delim@RIGHT\norm@star{AS}^{2}+\mu\varepsilon\leqslant\MT@delim@RIGHT\norm@star{A}^{2}+2\mu\varepsilon$
($\|TAS\|$ bound)

######  Remark 10.5.

In the streaming model, this procedure can be done in a constant number of
passes. Further, only need sampling and query access to $b$ to give a
$\log\frac{1}{\delta}$ failure probability: if we only have query access to
$b$, then sketch of $S$ can use the probability distribution corresponding to
$\frac{\|A_{*,i}\|^{2}}{\lVert A\rVert_{F}^{2}}$ suffices to give the
necessary guarantee with probability $\geqslant 0.99$.

#### One Clenshaw iteration.

We are trying to perform the odd Clenshaw recurrence defined in Section 8.2.
The matrix analogue of this is

$\displaystyle u_{k}$
$\displaystyle=2(2AA^{\dagger}-I)u_{k+1}-u_{k+2}+2a_{2k+1}Ab,$ (Odd Matrix
Clenshaw) $\displaystyle u$ $\displaystyle=\tfrac{1}{2}(u_{0}-u_{1}).$

We now show how to compute the next iterate $u_{k}$ given the previous two
iterates in the form $u_{k+1}=(SA)^{\dagger}v_{k+1}$ and
$u_{k+2}=(SA)^{\dagger}v_{k+2}$, along with $b$, for all $k\geqslant 0$. The
analysis begins by showing that $u_{k}$ is
$\varepsilon\MT@delim@RIGHT\norm@star{v_{k+1}}$-close to the output of an
exact, zero-error iteration. Next, we show that
$\MT@delim@RIGHT\norm@star{v_{k}}$ is $O(d)$-close to its zero-error iteration
value. Finally, we show that these errors don’t accumulate too much towards
the final outcome.

First, we note that the runtime of each iteration is $O(r)$, since this is the
cost of producing the sketches $B^{(k)}$ and $B_{\dagger}^{(k)}$, and actually
computing the iteration costs $O(r+s+t)=O(r)$.

Starting from Eq. Odd Matrix Clenshaw, we take the following series of
approximations by applying intermediate sketches:

$\begin{split}4AA^{\dagger}u_{k+1}-\,&2u_{k+1}-u_{k+2}+a_{2k+1}Ab\\\
&\approx_{1}4AS(AS)^{\dagger}u_{k+1}-2u_{k+1}-u_{k+2}+a_{2k+1}Ab\\\
&=AS\MT@delim@RIGHT\Paren@star{4(AS)^{\dagger}(AS)v_{k+1}-2v_{k+1}-v_{k+2}}+a_{2k+1}Ab\\\
&\approx_{2}AS\MT@delim@RIGHT\Paren@star{4(TAS)^{\dagger}(TAS)v_{k+1}-2v_{k+1}-v_{k+2}}+a_{2k+1}Ab\\\
&\approx_{3}AS\MT@delim@RIGHT\Paren@star{4(TAS)^{\dagger}(TAS)v_{k+1}-2v_{k+1}-v_{k+2}+a_{2k+1}S^{\dagger}b}\\\
&\approx_{4}AS(4(TAS)^{\dagger}B^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+a_{2k+1}S^{\dagger}b)\\\
&\approx_{5}AS(4B_{\dagger}^{(k)}B^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+a_{2k+1}S^{\dagger}b).\end{split}$

We note that sketching without the best suffices to get a
$d\|A\|_{F}^{4}/(\mu\varepsilon)^{4}$ runtime, since we can now compute the
expression
$4(TAS)^{\dagger}(TAS)v_{k+1}-2v_{k+1}-v_{k+2}+a_{2k+1}S^{\dagger}b$ and call
it $v_{k}$. We improve on this using the best sketch from Section 6.1.

As for approximation error, let
$\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\varepsilon_{4}$ and
$\varepsilon_{5}$ be the errors introduced in the corresponding approximation
above. Using the previously established bounds on $S$ and $T$,

$\displaystyle\varepsilon_{1}$ $\displaystyle\leqslant
4\|AA^{\dagger}-AS(AS)^{\dagger}\|\|u_{k+1}\|\leqslant
4\mu\varepsilon\|u_{k+1}\|$ by Eq. $AA^{\dagger}$ AMP
$\displaystyle\varepsilon_{2}$ $\displaystyle\leqslant
4\|AS\|\|((AS)(AS)^{\dagger}-(TAS)(TAS)^{\dagger})v_{k+1}\|\leqslant
5\mu\varepsilon\|v_{k+1}\|$ by Eq. $(AS)^{\dagger}AS$ AMP
$\displaystyle\varepsilon_{3}$
$\displaystyle\leqslant\MT@delim@RIGHT\abs@star{a_{2k+1}}\|Ab-
ASS^{\dagger}b\|\leqslant\MT@delim@RIGHT\abs@star{a_{2k+1}}\mu\varepsilon$ by
Eq. $Ab$ AMP

The bounds on $\varepsilon_{4}$ and $\varepsilon_{5}$ follow from the bounds
in Section 6.1 applied to $TAS$. With probability $\geqslant 1-1/(100d)$, the
following hold:

$\displaystyle\varepsilon_{4}$ $\displaystyle\leqslant
4\|AS(TAS)^{\dagger}(TAS-B^{(k)})v_{k+1}\|$ $\displaystyle\leqslant
4\sqrt{100d/r}\lVert AS(TAS)^{\dagger}\rVert_{F}\lVert
TAS\rVert_{F}\MT@delim@RIGHT\norm@star{v_{k+1}}$ by Corollary 6.3
$\displaystyle\leqslant
d^{-3/2}\mu\varepsilon\MT@delim@RIGHT\norm@star{v_{k+1}}$ Similarly,
$\displaystyle\varepsilon_{5}$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@nostar{4AS((TAS)^{\dagger}-B_{\dagger}^{(k)})B^{(k)}v_{k+1}}$
$\displaystyle\leqslant 4\sqrt{100d/r}\lVert
AS\rVert_{F}\MT@delim@RIGHT\norm@nostar{B^{(k)}v_{k+1}}$ by Corollary 6.3
$\displaystyle\leqslant 4\sqrt{100d/r}\lVert
AS\rVert_{F}\MT@delim@RIGHT\parens@star{\MT@delim@RIGHT\norm@nostar{TASv_{k+1}}+\sqrt{100d/r}\lVert
I_{t}\rVert_{F}\lVert TAS\rVert_{F}\MT@delim@RIGHT\norm@star{v_{k+1}}}$ by
Corollary 6.3 $\displaystyle\leqslant
d^{-3/2}\mu\varepsilon\MT@delim<EMAIL_ADDRESS>by taking
$r=\Theta(d^{4}(s+t))=\Theta(\frac{d^{4}\lVert
A\rVert_{F}^{4}}{(\mu\varepsilon)^{2}}\log\frac{n}{\delta})$.

Call the cumulative error for this round $\varepsilon^{(k)}$.

$\displaystyle\varepsilon^{(k)}:=\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{3}+\varepsilon_{4}+\varepsilon_{5}\lesssim\mu\varepsilon\Big{(}\|v_{k}\|+\MT@delim@RIGHT\abs@star{a_{2k+1}}\Big{)}.$

#### Error accumulation across iterations.

A simple argument shows that the error of the full algorithm is the sum of the
errors for each round. In other words, after completing the iteration, we have
a vector $u$ such that

$\displaystyle\MT@delim@RIGHT\norm@star{u-p(A)b}\leqslant\sum_{k=0}^{d}\varepsilon^{(k)}\lesssim\mu\varepsilon\sum_{k=0}^{d}(\|v_{k}\|+\MT@delim@RIGHT\abs@star{a_{2k+1}})\leqslant\varepsilon+\sum_{k=0}^{d}\mu\varepsilon\|v_{k}\|.$
(48)

So, it suffices to bound the $v_{k}$’s. The recursions defining them is

$\displaystyle v_{k}$
$\displaystyle=4R_{\dagger}^{(k)}R^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+a_{2k+1}S^{\dagger}b$
$\displaystyle=2(2(TAS)^{\dagger}(TAS)-I)v_{k+1}-v_{k+2}+a_{2k+1}S^{\dagger}b+4(B_{\dagger}^{(k)}B^{(k)}-(TAS)^{\dagger}(TAS))v_{k+1}$
This solves to $\displaystyle v_{k}$
$\displaystyle=\sum_{i=k}^{d}U_{i-k}(2(TAS)^{\dagger}(TAS)-I)\MT@delim@RIGHT\parens@nostar[\Big{]}{a_{2i+1}S^{\dagger}b+4(B_{\dagger}^{(i)}B^{(i)}-(TAS)^{\dagger}(TAS))v_{i+1}}.$

We bound the second moment (averaging over the randomness of the $B^{(k)}$’s),
using that the $B^{(k)}$, $B^{(k)}_{\dagger}$ and $B^{(k^{\prime})}$ are all
drawn independently. First, notice that the means are bounded.

$\bar{v}_{k}=\operatorname*{\mathbf{E}}_{[k,d]}\left[v_{k}\right]=\sum_{i=k}^{d}U_{i-k}(2(TAS)^{\dagger}(TAS)-I)a_{2k+1}Sb.$

Using sub-multiplicativity of the operator norm,

$\displaystyle\MT@delim@RIGHT\norm@star{\bar{v}_{k}}$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@nostar[\Big{]}{\sum_{i=k}^{d}U_{i-k}(2(TAS)^{\dagger}(TAS)-I)a_{2k+1}}\MT@delim@RIGHT\norm@star{S^{\dagger}b}$
$\displaystyle=\MT@delim@RIGHT\norm@nostar[\Big{]}{\sum_{i=k}^{d}U_{i-k}(T_{2}(x))a_{2k+1}}_{\operatorname{Spec}(TAS)}\MT@delim@RIGHT\norm@star{S^{\dagger}b}$
$\displaystyle\leqslant
e\Big{\lVert}\sum_{i=k}^{d}U_{i-k}(2x^{2}-1)a_{2k+1}\Big{\rVert}_{\mathrm{sup}}\MT@delim@RIGHT\norm@star{S^{\dagger}b}$
by Lemma 5.4 and Eq. $\|TAS\|$ bound $\displaystyle\leqslant
2e\lVert\sum_{i=k}^{d}U_{i-k}(2x^{2}-1)a_{2k+1}\rVert_{\mathrm{sup}}$
$\displaystyle\leqslant 2e\frac{d-k+1}{\mu d^{2}}$ by Item 10.2(c)

We now compute the second moment of $v_{k}$.

$\displaystyle\operatorname*{\mathbf{E}}_{[k,d]}\left[\MT@delim@RIGHT\norm@star{v_{k}}^{2}\right]$
$\displaystyle=\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\operatorname*{\mathbf{E}}_{[k,d]}\left[\MT@delim@RIGHT\Norm@star{v_{k}-\bar{v}_{k}}^{2}\right]$
$\displaystyle=\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\operatorname*{\mathbf{E}}_{[k,d]}\left[\MT@delim@RIGHT\Norm@star{\sum_{i=k}^{d}U_{i-k}(2(TAS)^{\dagger}(TAS)-I)4(B_{\dagger}^{(i)}B^{(i)}-(TAS)^{\dagger}(TAS))v_{i+1}}^{2}\right]$
We use independence of the $B^{(i)}$’s to note that the expectation of cross
terms is $0$, and therefore,
$\displaystyle=\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\operatorname*{\mathbf{E}}_{[k,d]}\left[\sum_{i=k}^{d}\MT@delim@RIGHT\Norm@star{U_{i-k}(2(TAS)^{\dagger}(TAS)-I)4(B_{\dagger}^{(i)}B^{(i)}-(TAS)^{\dagger}(TAS))v_{i+1}}^{2}\right]$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\sum_{i=k}^{d}16(i-k+1)^{2}\operatorname*{\mathbf{E}}_{[i,d]}\left[\MT@delim@RIGHT\Norm@star{(B_{\dagger}^{(i)}B^{(i)}-(TAS)^{\dagger}(TAS))v_{i+1}}^{2}\right]$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\sum_{i=k}^{d}\frac{16(i-k+1)^{2}}{d^{4}}\operatorname*{\mathbf{E}}_{[i+1,d]}\left[\MT@delim@RIGHT\norm@star{v_{i+1}}^{2}\right],$

where the last line follows from the computation

$\displaystyle\operatorname*{\mathbf{E}}_{k}\left[\MT@delim@RIGHT\norm@star{(B_{\dagger}^{(k)}B^{(k)}-(TAS)^{\dagger}(TAS))v_{k+1}}^{2}\right]$
$\displaystyle=\operatorname*{\mathbf{E}}_{k}\left[\MT@delim@RIGHT\norm@star{B_{\dagger}^{(k)}B^{(k)}v_{k+1}}^{2}\right]-\MT@delim@RIGHT\Norm@star{(TAS)^{\dagger}(TAS)v_{k+1}}^{2}$
$\displaystyle=\operatorname*{\mathbf{E}}_{k}\left[(B^{(k)}v_{k+1})^{\dagger}(B_{\dagger}^{(k)})^{\dagger}B_{\dagger}^{(k)}(B^{(k)}v_{k+1})\right]-\MT@delim@RIGHT\Norm@star{(TAS)^{\dagger}(TAS)v_{k+1}}^{2}$
$\displaystyle\leqslant\operatorname*{\mathbf{E}}_{k}\left[(B^{(k)}v_{k+1})^{\dagger}((TAS)(TAS)^{\dagger}+\frac{1}{r}\operatorname{Tr}(I_{s})\lVert
TAS\rVert_{F}^{2}I_{t})(B^{(k)}v_{k+1})\right]-\MT@delim@RIGHT\Norm@star{(TAS)^{\dagger}(TAS)v_{k+1}}^{2}$
by Eq. 11
$\displaystyle\leqslant\operatorname*{\mathbf{E}}_{k}\Bigl{[}v_{k+1}^{\dagger}(TAS)^{\dagger}((TAS)(TAS)^{\dagger}+\frac{s}{r}\lVert
TAS\rVert_{F}^{2}I_{t})(TAS)v_{k+1}$
$\displaystyle\qquad+v_{k+1}^{\dagger}\frac{1}{r}\operatorname{Tr}((TAS)(TAS)^{\dagger}+\frac{s}{r}\lVert
TAS\rVert_{F}^{2}I_{t})\lVert
TAS\rVert_{F}^{2}v_{k+1}\Bigr{]}-\MT@delim@RIGHT\Norm@star{(TAS)^{\dagger}(TAS)v_{k+1}}^{2}$
by Eq. 11 $\displaystyle=\frac{s}{r}\lVert
TAS\rVert_{F}^{2}\MT@delim@RIGHT\norm@star{TASv_{k+1}}^{2}+\frac{1}{r}\lVert
TAS\rVert_{F}^{4}\MT@delim@RIGHT\norm@star{v_{k+1}}^{2}+\frac{st}{r^{2}}\lVert
TAS\rVert_{F}^{4}\MT@delim@RIGHT\norm@star{v_{k+1}}^{2}$
$\displaystyle\lesssim\MT@delim@RIGHT\parens@star{\frac{s\lVert
A\rVert_{F}^{2}}{r}+\frac{\lVert A\rVert_{F}^{4}}{r}+\frac{st\lVert
A\rVert_{F}^{4}}{r^{2}}}\MT@delim@RIGHT\norm@star{v_{k+1}}^{2}$ by Eqs.
$\|TAS\|_{F}$ bound and $\|TAS\|$ bound $\displaystyle=\frac{\lVert
A\rVert_{F}^{2}}{r}\MT@delim@RIGHT\parens@nostar[\Big{]}{s+\lVert
A\rVert_{F}^{2}+\frac{t\lVert
A\rVert_{F}^{2}}{r}}\MT@delim@RIGHT\norm@star{v_{k+1}}^{2}\leqslant\frac{1}{d^{4}}\MT@delim@RIGHT\norm@star{v_{k+1}}^{2}$

Now, we can bound the second moments of all $\MT@delim@RIGHT\norm@star{v_{k}}$
by bounding the recurrence. We define the following recurrence $c_{k}$ to
satisfy
$\operatorname*{\mathbf{E}}_{[k,d]}[\MT@delim@RIGHT\norm@star{v_{k}}^{2}]\leqslant
c_{k}$, and then inductively show that
$c_{k}\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\frac{(d-k+1)^{2}}{\mu^{2}d^{5}}\leqslant
30\frac{(d-k+1)^{2}}{\mu^{2}d^{4}}$:

$\displaystyle c_{k}$
$\displaystyle=\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\sum_{i=k}^{d}\frac{(i-k+1)^{2}}{d^{4}}c_{i+1}$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\sum_{i=k}^{d}\frac{(i-k+1)^{2}}{d^{4}}\frac{30(d-i+1)^{2}}{\mu^{2}d^{4}}$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\frac{30}{\mu^{2}d^{8}}\sum_{i=k}^{d}(i-k+1)^{2}(d-i+1)^{2}$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\frac{(d-k+1)^{5}}{\mu^{2}d^{8}}$
$\displaystyle\leqslant\MT@delim@RIGHT\norm@star{\bar{v}_{k}}^{2}+\frac{(d-k+1)^{2}}{\mu^{2}d^{5}}$

We have shown that
$\operatorname*{\mathbf{E}}[\MT@delim@RIGHT\norm@star{v_{k}-\bar{v}_{k}}^{2}]\leqslant\frac{1}{d}\MT@delim@RIGHT\parens@nostar{\frac{(d-k+1)}{\mu
d^{2}}}^{2}$. By Markov’s inequality, with probability $\geqslant 0.999$, we
have that for all $k$,
$\MT@delim@RIGHT\norm@star{v_{k}}\lesssim\frac{d-k+1}{\mu d^{2}}$. Returning
to the final error bound Eq. 48,

$\displaystyle\sum_{k=0}^{d}\varepsilon^{(k)}\lesssim\varepsilon+\sum_{k=0}^{d}\mu\varepsilon\|v_{k}\|\lesssim\varepsilon+\sum_{k=0}^{d}\varepsilon(d-k+1)/d^{2}\lesssim\varepsilon.$
(49)

This concludes the error analysis.

#### Output description properties.

After the iteration concludes, we can compute $u$ by computing
$x=\tfrac{1}{2}S(v_{0}-v_{1})$ in linear $O(s)$ time. Then,
$u=\tfrac{1}{2}(u_{0}-u_{1})=\tfrac{1}{2}AS(v_{0}-v_{1})=Ax$. Note that though
$x\in\mathbbm{C}^{n}$, its sparsity is at most the sparsity of $x$, which is
bounded by $s$.

Further, using the prior bounds on $v_{0}$ and $v_{1}$, we have that

$\displaystyle\sum_{j=1}^{n}\MT@delim@RIGHT\norm@star{A_{*,j}}^{2}\MT@delim@RIGHT\abs@star{x_{i}}^{2}$
$\displaystyle=\sum_{j=1}^{s}\MT@delim@RIGHT\norm@star{[SA]_{*,j}}^{2}\MT@delim@RIGHT\abs@star{\tfrac{1}{2}(v_{0}-v_{1})_{i}}^{2}$
$\displaystyle\leqslant\sum_{j=1}^{s}\tfrac{2}{s}\lVert
A\rVert_{F}^{2}\MT@delim@RIGHT\abs@star{\tfrac{1}{2}(v_{0}-v_{1})_{i}}^{2}$
$\displaystyle\leqslant\tfrac{2}{s}\lVert
A\rVert_{F}^{2}\MT@delim@RIGHT\norm@star{\tfrac{1}{2}(v_{0}-v_{1})}^{2}$
$\displaystyle\leqslant\tfrac{2}{s}\lVert
A\rVert_{F}^{2}(\MT@delim@RIGHT\norm@star{v_{0}}^{2}+\MT@delim@RIGHT\norm@star{v_{1}}^{2})$
$\displaystyle\leqslant 2\lVert A\rVert_{F}^{2}/(\sqrt{s}\mu d)^{2}.$

###  10.2 Generalizing to even polynomials

We also obtain an analogous result for even polynomials. For the most part,
changes are superficial; the red text indicates differences from Theorem 10.2.
The major difference is that the representation of the output is
$A^{\dagger}x+\eta b$ instead of $Ax$, which results from a difference in the
dimension of the output, and the allowing for constant terms in $p$.

######  Theorem 10.6.

Suppose we are given sampling and query access to $A\in\mathbb{C}^{m\times n}$
and $b\in\mathbb{C}^{n}$ with
$\MT@delim@RIGHT\norm@star{A},\MT@delim@RIGHT\norm@star{b}\leqslant 1$, and
suppose we are given a $(2d)$-degree even polynomial, written in its Chebyshev
coefficients as

$\displaystyle{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}p(x)=\sum_{i=0}^{d}a_{2i}T_{2i}(x).}$

Then we can output a vector
$x\in\mathbbm{C}^{{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}m}}$
and $\eta\in\mathbbm{C}$ such that
$\|{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}A^{\dagger}x+\eta
b}-p(A)b\|\leqslant\varepsilon\MT@delim@RIGHT\norm@star{p}_{\operatorname{Spec}(A)}$
in time

$\displaystyle
O\MT@delim@RIGHT\parens@star{\min\MT@delim@RIGHT\braces@star{\operatorname{nnz}(A),\frac{\|A\|_{F}^{4}}{(\mu\varepsilon)^{4}}\log^{2}(n)}+\frac{d^{5}\|A\|_{F}^{4}}{(\mu\varepsilon)^{2}}\log(n)},$

assuming $\mu$ satisfies the following bounds:

  1. (a)

$\mu\varepsilon\leqslant\frac{1}{100d^{2}}$;

  2. (b)

$\mu\sum_{i=0}^{d}\MT@delim@RIGHT\abs@star{\tilde{a}_{2i}}\leqslant 1$;

  3. (c)

$\mu\lVert\sum_{i=k}^{d}4\tilde{a}_{2i+2}x\cdot
U_{i-k}(T_{2}(x))\rVert_{\mathrm{sup}}\leqslant\frac{d-k+1}{d^{2}}$ for all
$0\leqslant k\leqslant d$.

The output description has the additional properties that

$\displaystyle\sum_{i}\MT@delim@RIGHT\norm@star{A_{{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}i,*}}}^{2}\MT@delim@RIGHT\abs@star{x_{i}}^{2}\lesssim\frac{\varepsilon^{2}}{\log(n)d^{2}}\qquad\lVert
x\rVert_{0}\lesssim\frac{\lVert
A\rVert_{F}^{2}}{(\mu\varepsilon)^{2}}\log(n),$

so that by Corollary 5.11, we can

  1. (a)

Compute entries of $Ax$ in $\lVert x\rVert_{0}=\lVert
A\rVert_{F}/\varepsilon^{2}$ time;

  2. (b)

Sample $i\in[n]$ with probability
$\frac{\MT@delim@RIGHT\abs@star{y_{i}}^{2}}{\MT@delim@RIGHT\norm@star{y}^{2}}$
in $(\frac{\lVert A\rVert_{F}^{2}}{(\mu
d)^{2}}{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}+p(0)^{2}})\frac{\lVert
A\rVert_{F}^{2}\log(n)}{(\mu\varepsilon\MT@delim@RIGHT\norm@star{y})^{2}}$
time;

  3. (c)

Estimate $\MT@delim@RIGHT\norm@star{y}^{2}$ to $\nu$ relative error in
$(\frac{\lVert A\rVert_{F}^{2}}{(\mu
d)^{2}}{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}+p(0)^{2}})\frac{\lVert
A\rVert_{F}^{2}\log(n)}{(\nu\mu\varepsilon\MT@delim@RIGHT\norm@star{y})^{2}}$
time.

######  Proposition 10.7.

In Theorem 10.6, we can always take $1/\mu\lesssim
d^{2}{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}\log(d)}\MT@delim@RIGHT\norm@star{p}_{\operatorname{Spec}(A)}$.

######  Algorithm 10.8 (Even singular value transformation).

Input: A matrix $A\in\mathbb{C}^{m\times n}$, vector $b\in\mathbb{C}^{n}$,
$\varepsilon>0$, a degree-$2d$ polynomial $p$, and the corresponding
coefficients in the Chebyshev basis, i.e.
$p(x)=\sum_{i=0}^{d}a_{2i}T_{2i}(x)$. Preprocessing sketches: Let
$s,t=\Theta(\frac{\lVert
A\rVert_{F}^{2}}{(\mu\varepsilon)^{2}}\log(\frac{n}{\delta}))$. Compute all
$\tilde{a}_{2k}=a_{2k}-a_{2k+2}+\cdots\pm a_{2d}$. P1. If
$\operatorname{SQ}({\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}A})$
and $\operatorname{SQ}(b)$ are not given, compute data structures to simulate
them in $O(1)$ time; P2. Sample $S\in\mathbbm{C}^{s\times m}$ to be
${\textsc{AMP}}(A^{\dagger},A,\\{\frac{\|A_{*,i}\|^{2}}{\lVert
A\rVert_{F}^{2}}\\})$ (Theorem 6.5) with sketch size $s$; P3. Sample
$T\in\mathbbm{C}^{n\times t}$ to be
${\textsc{AMP}}(SA,A^{\dagger}S^{\dagger},\\{\tfrac{1}{2}(\frac{\MT@delim@RIGHT\norm@nostar{[SA]_{*,i}}^{2}}{\lVert
SA\rVert_{F}^{2}}+\frac{\MT@delim@RIGHT\abs@star{b_{i}}^{2}}{\MT@delim@RIGHT\norm@star{b}^{2}})\\})$
(Theorem 6.5) with sketch size $t$; P4. Compute a data structure that can
respond to
$\operatorname{SQ}({\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}SAT})$
queries in $O(1)$ time; Clenshaw iteration: Let $r=\Theta(d^{4}\lVert
A\rVert_{F}^{2}(s+t))=\Theta(\frac{d^{4}\lVert
A\rVert_{F}^{4}}{(\mu\varepsilon)^{2}}\log\frac{n}{\delta})$. Then starting
with $v_{d+1}=v_{d+2}=\vec{0}^{s}$ and going until $v_{0}$, I1. Let
$B^{(k)}={\textsc{best}}({\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}SAT})$
and
$B_{\dagger}^{(k)}={\textsc{best}}(({\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}SAT})^{\dagger})$
with parameter $r$; I2. Compute
$v_{k}=2(2{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}B^{(k)}B_{\dagger}^{(k)}}-I)v_{k+1}-v_{k+2}+{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}4\tilde{a}_{2k+2}B^{(k)}b}$.
Output: Output
$x=\tfrac{1}{2}{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}S^{\dagger}}(v_{0}-v_{1})$
and $\eta=\tilde{a}_{0}$ that satisfies
$\MT@delim@RIGHT\Norm@star{{\color[rgb]{0.828125,0.2109375,0.21484375}\definecolor[named]{pgfstrokecolor}{rgb}{0.828125,0.2109375,0.21484375}A^{\dagger}x+\eta
b}-p(A)b}\leqslant\varepsilon$.

Recall the odd and even recurrences defined in Sections 8.2 and Even Clenshaw.

$\displaystyle u_{k}$
$\displaystyle=2(2AA^{\dagger}-I)u_{k+1}-u_{k+2}+2a_{2k+1}Ab,$ (Odd Matrix
Clenshaw) $\displaystyle u$ $\displaystyle=\tfrac{1}{2}(u_{0}-u_{1}).$

The matrix analogue of the even recurrence is identical except for the final
term is $4\tilde{a}_{2k+2}A^{\dagger}Ab$ instead of $2a_{2k+1}Ab$.

$\displaystyle\tilde{a}_{2k}$ $\displaystyle\coloneqq
a_{2k}-a_{2k+2}+a_{2k+4}-\cdots\pm a_{2d}$ $\displaystyle u_{k}$
$\displaystyle=2(2A^{\dagger}A-1)u_{k+1}-u_{k+2}+4\tilde{a}_{2k+2}A^{\dagger}Ab,$
(Even Matrix Clenshaw) $\displaystyle u$
$\displaystyle=\tilde{a}_{0}b+\tfrac{1}{2}(u_{0}-u_{1}).$

So, a roughly identical analysis works upon making the appropriate changes.

#### One (even) Clenshaw iteration.

We maintain $u_{k}$ as $(SA)^{\dagger}v_{k}$. The error analysis proceeds by
bounding

$\displaystyle 4A^{\dagger}Au_{k+1}-\,$ $\displaystyle
2u_{k+1}-u_{k+2}+4\tilde{a}_{2k+2}A^{\dagger}Ab$
$\displaystyle\approx_{1}4(SA)^{\dagger}(SA)u_{k+1}-2u_{k+1}-u_{k+2}+4\tilde{a}_{2k+2}A^{\dagger}Ab$
$\displaystyle=(SA)^{\dagger}\MT@delim@RIGHT\Paren@star{4(SA)(SA)^{\dagger}v_{k+1}-2v_{k+1}-v_{k+2}}+4\tilde{a}_{2k+2}A^{\dagger}Ab$
$\displaystyle\approx_{2}(SA)^{\dagger}\MT@delim@RIGHT\Paren@star{4(SAT)(SAT)^{\dagger}v_{k+1}-2v_{k+1}-v_{k+2}}+4\tilde{a}_{2k+2}A^{\dagger}Ab$
$\displaystyle\approx_{3}(SA)^{\dagger}\MT@delim@RIGHT\Paren@star{4(SAT)(SAT)^{\dagger}v_{k+1}-2v_{k+1}-v_{k+2}+4\tilde{a}_{2k+2}SATT^{\dagger}b}$
$\displaystyle\approx_{4}(SA)^{\dagger}\MT@delim@RIGHT\Paren@star{4(SAT)B_{\dagger}^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+4\tilde{a}_{2k+2}SATT^{\dagger}b}$
$\displaystyle\approx_{5}(SA)^{\dagger}\MT@delim@RIGHT\Paren@star{4B^{(k)}B_{\dagger}^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+4\tilde{a}_{2k+2}SATT^{\dagger}b}$
$\displaystyle\approx_{6}(SA)^{\dagger}\MT@delim@RIGHT\Paren@star{4B^{(k)}B_{\dagger}^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+4\tilde{a}_{2k+2}B^{(k)}T^{\dagger}b}$

Using the same approaches as in the odd case, we can show that the error
incurred to compute $u_{k}$ is
$O(\mu\varepsilon(\MT@delim@RIGHT\norm@star{v_{k+1}}+\MT@delim@RIGHT\abs@star{\tilde{a}_{2k+2}}))$.

#### Error accumulation across iterations.

The corresponding recurrence for $v_{k}$ is

$\displaystyle v_{k}$
$\displaystyle=4B^{(k)}B_{\dagger}^{(k)}v_{k+1}-2v_{k+1}-v_{k+2}+4\tilde{a}_{2k+2}B^{(k)}T^{\dagger}b$
$\displaystyle=(4(SAT)(SAT)^{\dagger}-2I)v_{k+1}-v_{k+2}$ $\displaystyle\hskip
50.00008pt+\Big{[}4\tilde{a}_{2k+2}SATT^{\dagger}b+4(SAT(SAT)^{\dagger}-B^{(k)}B_{\dagger}^{(k)})v_{k+1}+(B^{(k)}-SAT)T^{\dagger}b\Big{]}$
This solves to $\displaystyle v_{k}$
$\displaystyle=\sum_{i=k}^{d}U_{i-k}(T_{2}(SAT))\Big{[}4\tilde{a}_{2k+2}SATT^{\dagger}b+4(SAT(SAT)^{\dagger}-B^{(k)}B_{\dagger}^{(k)})v_{k+1}+(B^{(k)}-SAT)T^{\dagger}b\Big{]}$

From here, the identical analysis applies. The bound from Item 10.6(c)
corresponds to $\MT@delim@RIGHT\norm@star{\operatorname*{\mathbf{E}}[v_{k}]}$
being bounded by $\frac{d-k+1}{\mu d^{2}}$.

#### Output description properties.

The main change is that the description is $A^{\dagger}x+\tilde{a}_{0}b$, with
an additional scalar $b$ term. We first use the prior analysis to get
$\operatorname{SQ}_{\phi}(A^{\dagger}x)$ for $\phi=\lVert
x\rVert_{0}\frac{\sum_{i}\MT@delim@RIGHT\norm@star{A_{i,*}}^{2}\MT@delim@RIGHT\abs@star{x_{i}}^{2}}{\MT@delim@RIGHT\norm@star{A^{\dagger}x}^{2}}=\frac{\lVert
A\rVert_{F}^{2}}{(\mu d)^{2}\MT@delim@RIGHT\norm@star{A^{\dagger}x}^{2}}$.
Then we use Lemma 5.8 to get $\operatorname{SQ}_{\varphi}(A^{\dagger}x+\eta
b)$ for $\phi=2(\frac{\lVert A\rVert_{F}^{2}}{(\mu
d)^{2}}+\eta^{2}\MT@delim@RIGHT\norm@star{b}^{2})/\MT@delim@RIGHT\norm@star{A^{\dagger}x+\eta
b}^{2}$. Finally, we note that $\tilde{a}_{0}=p(0)$.

###  10.3 Bounding iterates of singular value transformation

######  Theorem 10.9.

Let $p(x)$ be a odd degree $2d+1$ polynomial such that $\lVert
p\rVert_{\mathrm{sup}}\leqslant 1$. Then it suffices for Theorem 10.2 to take

$\displaystyle\mu=O(d^{2}\log(d)).$

###### Proof.

All there is to prove is that

$\displaystyle\Big{\lVert}\sum_{i=k}^{d}a_{i}U_{i-k}(T_{2}(x))\Big{\rVert}_{\mathrm{sup}}\leqslant(d-k+1)\log^{2}(d+1)$

First, we note that it suffices to prove the bound for
$\sum_{i=k}^{d}a_{i}U_{i-k}(y)$, since we get the bound by setting
$y=T_{2}(x)$. The strategy to proving this is the same as the one we discussed
before:

$\displaystyle\sum_{i}a_{i}U_{i-k}(T_{2}(x))$
$\displaystyle=\sum_{i}a_{i}\sum_{j\geqslant
0}T_{i-k-2j}(y)(1+\MT@delim@RIGHT\iver@star{i-k-2j\neq 0})$
$\displaystyle=\sum_{j\geqslant
0}\sum_{i}a_{i}T_{i-k-2j}(y)(1+\MT@delim@RIGHT\iver@star{i-k-2j\neq 0})$
$\displaystyle=\sum_{j\geqslant
0}\sum_{i}a_{i+k+2j}T_{i}(x)(1+\MT@delim@RIGHT\iver@star{i\neq 0})$
$\displaystyle=\sum_{i}T_{i}(x)(1+\MT@delim@RIGHT\iver@star{i\neq
0})\sum_{j\geqslant 0}a_{i+k+2j}$
$\displaystyle\leqslant\sum_{i}(1+\MT@delim@RIGHT\iver@star{i\neq
0})\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{j\geqslant 0}a_{i+k+2j}}$
$\displaystyle\leqslant\sum_{i}(1+\MT@delim@RIGHT\iver@star{i\neq
0})(32+8\ln^{2}(d+1))\MT@delim@RIGHT\iver@star{i\leqslant d-k}$
$\displaystyle\lesssim(d-k+1)\log^{2}(d+1)$

∎

######  Theorem 10.10.

Let $p(x)$ be a even degree-$(2d)$ polynomial such that $\lVert
p\rVert_{\mathrm{sup}}\leqslant 1$. Then it suffices for Theorem 10.6 to take

$\displaystyle\mu=O(d^{2}\log^{2}(d))$

###### Proof.

All there is to prove is that

$\displaystyle\Big{\lVert}\sum_{i=k}^{d}4\tilde{a}_{2i+2}xU_{i-k}(T_{2}(x))\Big{\rVert}_{\mathrm{sup}}=\Big{\lVert}\sum_{i=k}^{d}2\tilde{a}_{2i+2}U_{2(i-k)+1}(x)\Big{\rVert}_{\mathrm{sup}}$
$\displaystyle\leqslant(d-k+1)\log^{2}(d+1)$

The strategy to proving this is the same as the one we discussed before:

$\displaystyle\sum_{i}\tilde{a}_{2i+2}U_{2(i-t)+1}(x)$
$\displaystyle=\sum_{i}\sum_{k\geqslant
0}(-1)^{k}a_{2(i+k+1)}\sum_{j=0}^{i-t}T_{2j+1}(x)$
$\displaystyle=\sum_{i,j,k}(-1)^{k}a_{2(i+k+1)}T_{2j+1}(x)\MT@delim@RIGHT\iver@star{j,k\geqslant
0}\MT@delim@RIGHT\iver@star{j\leqslant i-t}$ $\displaystyle=\sum_{j\geqslant
0}T_{2j+1}\sum_{i,k}(-1)^{k}a_{2(i+k+1)}\MT@delim@RIGHT\iver@star{k\geqslant
0}\MT@delim@RIGHT\iver@star{j\leqslant i-t}$ $\displaystyle=\sum_{j\geqslant
0}T_{2j+1}\sum_{i,k,\ell}(-1)^{k}a_{2\ell}\MT@delim@RIGHT\iver@star{\ell=i+k+1}\MT@delim@RIGHT\iver@star{k\geqslant
0}\MT@delim@RIGHT\iver@star{j\leqslant i-t}$ $\displaystyle=\sum_{j\geqslant
0}T_{2j+1}\sum_{\ell}a_{2\ell}\sum_{k}(-1)^{k}\MT@delim@RIGHT\iver@star{k\geqslant
0}\MT@delim@RIGHT\iver@star{j\leqslant\ell-k-1-t}$
$\displaystyle=\sum_{j\geqslant
0}T_{2j+1}\sum_{\ell}a_{2\ell}\sum_{k}(-1)^{k}\MT@delim@RIGHT\iver@star{0\leqslant
k\leqslant\ell-j-1-t}$ $\displaystyle=\sum_{j\geqslant
0}T_{2j+1}\sum_{\ell}a_{2\ell}\MT@delim@RIGHT\iver@star{\ell-j-1-t\in
2\mathbb{Z}_{\geqslant 0}}$ $\displaystyle=\sum_{j\geqslant
0}T_{2j+1}\sum_{\ell\geqslant 0}a_{2(2\ell+j+1+t)}$
$\displaystyle\leqslant\sum_{j\geqslant
t+1}^{d}\MT@delim@RIGHT\abs@nostar[\Big{]}{\sum_{\ell\geqslant
0}a_{2(2\ell+j)}}$ $\displaystyle\leqslant(d-t)(32+8\ln^{2}(d+1))$
$\displaystyle\lesssim(d-k+1)\log^{2}(d+1)$

∎

##  11 Dequantizing algorithms

###  11.1 Recommendation systems

######  Lemma 11.1 (Polynomial approximations of the sign function (Lemma 25
of [GSLW19])).

For all $\delta>0$ and $\varepsilon\in(0,\frac{1}{2})$ there exists an
efficiently computable odd polynomial $p\in\mathbb{R}[x]$ with
$\deg(p)=O(\frac{1}{\delta}\log\frac{1}{\varepsilon})$

  * •

for all $x\in[-2,2]$, $\MT@delim@RIGHT\abs@star{p(x)}\leqslant 1$

  * •

for all $x\in[-2,-\delta]\cup[\delta,2]$,
$\MT@delim@RIGHT\abs@star{p(x)-\operatorname{sign}(x)}\leqslant\varepsilon$

######  Lemma 11.2.

For $\varsigma,\eta\in(0,\tfrac{1}{2})$ and $\sigma\in[-1,1]$, there exists an
even polynomial $p(x)\in\mathbb{R}[x]$ of degree
$O(\frac{1}{\eta\sigma}\log\frac{1}{\varsigma})$ such that $p(x)\in[0,1]$ for
all $x\in[-1,1]$ and

$\displaystyle p(x)\in\begin{cases}[1-\varsigma,1]&\text{for }1\leqslant
x\leqslant-(1+\eta)\sigma\\\ [0,\varsigma]&\text{for }-(1-\eta)\sigma\leqslant
x\leqslant(1-\eta)\sigma\\\ [1-\varsigma,1]&\text{for }(1+\eta)\sigma\leqslant
x\leqslant 1\end{cases}$

###### Proof.

Let $s(x)$ be the approximation in Lemma 11.1 with error paramters
$\delta\leftarrow\eta\sigma$ and $\varepsilon\leftarrow\frac{\varsigma}{2}$.
Then

$\displaystyle
p(x)\coloneqq(1-\varsigma)\Big{(}1+\frac{s(x-\sigma)+s(-x-\sigma)}{2}\Big{)}$

satisfies the desired parameters. ∎

######  Corollary 11.3 (De-quantizing Recommendation Systems).

The recommendation systems problem is to sample from
$[A_{\geqslant\sigma,\frac{1}{6}}]_{i,*}$. Using Theorem 10.6 on the
polynomial from Lemma 11.2, given $\operatorname{nnz}(A)$ preprocessing time
to make a data structure, we can compute an $x$ such that
$\MT@delim@RIGHT\norm@star{A^{\dagger}x-A_{\geqslant\sigma,\frac{1}{3}}^{\dagger}e_{i}}\leqslant\varepsilon$
in $\widetilde{O}\MT@delim@RIGHT\Paren@star{\frac{\lVert
A\rVert_{F}^{4}}{\sigma^{9}\varepsilon^{2}}}$ time, and sample from this $x$
in $\widetilde{O}\MT@delim@RIGHT\Paren@star{\frac{\lVert
A\rVert_{F}^{4}}{\sigma^{6}\varepsilon^{2}\MT@delim@RIGHT\norm@star{A^{\dagger}x}^{2}}}$
time.

###  11.2 Linear regression

######  Lemma 11.4 (Polynomial approximations of $1/x$, [GSLW19, Lemma 40],
following [CKS17]).

Let $\kappa>1$ and $0<\varepsilon\leqslant\delta\leqslant\tfrac{1}{2}$. There
is an odd polynomial $p(x)$ such that

$\displaystyle\MT@delim@RIGHT\abs@star{p(x)-\frac{\delta}{2x}}\leqslant\varepsilon$

is $\varepsilon$-close to $1/x$ on the domain
$[-1,1]\setminus(-\tfrac{1}{\kappa},\tfrac{1}{\kappa})$. Let
$J\coloneqq\MT@delim@RIGHT\ceil@star{\sqrt{b\log(4b/\varepsilon)}}$, then the
$O(\kappa\log\frac{\kappa}{\varepsilon})$-degree odd real polynomial

$\displaystyle g(x)\coloneqq
4\sum_{j=0}^{J}(-1)^{j}\MT@delim@RIGHT\bracks@star{\frac{\sum_{i=j+1}^{b}\binom{2b}{b+i}}{2^{2b}}}T_{2j+1}(x$

is $\varepsilon$-close to $f(x)$ on the interval $[-1,1]$, moreover
$\MT@delim@RIGHT\abs@star{g(x)}\leqslant
4J=O(\kappa\log(\frac{\kappa}{\varepsilon}))$ on this interval.

This is the polynomial that is applied to perform regression in the QSVT
framework, with slight modifications to make the polynomial suitably bounded.

######  Corollary 11.5.

If we have $\operatorname{SQ}(A)$ and $\operatorname{SQ}(b)$, we can apply
Theorem 10.2 on $g(x)/\kappa$ to error $\varepsilon/\kappa$, where $g(x)$ is
the polynomial in Lemma 11.4. Then after $O(\operatorname{nnz}(A))$ or
$O(\kappa^{12}\|A\|_{F}^{4}\varepsilon^{-4}\operatorname{polylog}(n\kappa/\varepsilon))$,
we compute This allows us to get a representation of a vector $y$ such that
$\MT@delim@RIGHT\norm@star{y-g(A)b}\leqslant\varepsilon\MT@delim@RIGHT\norm@star{b}$
in

$\displaystyle
O\MT@delim@RIGHT\parens@star{\frac{\kappa^{11}\|A\|_{F}^{4}}{\varepsilon^{2}}\operatorname{polylog}(n\kappa/\varepsilon)}$

time. We can also sample and query from the output in this amount of time.

###  11.3 Hamiltonian simulation

In this Subsection, we provide the corollary for de-quantizing Hamiltonian
simulation. We begin by recalling the following definition:

######  Definition 11.6 (Bessel functions of the first kind).

Given $i\in\mathbb{N}$,

$J_{i}(x)=\sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!\Gamma\MT@delim@RIGHT\Paren@star{i+m+1}}\MT@delim@RIGHT\Paren@star{\frac{x}{2}}^{2m+i},$

where $\Gamma$ is the gamma function.

We also need the following polynomial approximation to trigonometric
functions:

######  Lemma 11.7 (Polynomial approximation to trigonometric functions,
[GSLW19, Lemma 57]).

Given $t\in\mathbb{R}$ and $\varepsilon\in(0,1/e)$, let
$r=\Theta\MT@delim@RIGHT\Paren@star{t+\frac{\log(1/\varepsilon)}{\log\log(1/\varepsilon)}}$.
Then, the following degree $2r$ and $2r+1$ satisfy that for all $x\in[-1,1]$,

$\begin{split}\MT@delim@RIGHT\Norm@star{\cos(tx)-J_{0}(t)+2\sum_{i\in[1,r]}\MT@delim@RIGHT\Paren@star{-1}^{i}J_{2i}(t)T_{2i}(x)}_{\textrm{sup}}&\leqslant\varepsilon,\\\
\MT@delim@RIGHT\Norm@star{\sin(tx)-2\sum_{i\in[0,r]}\MT@delim@RIGHT\Paren@star{-1}^{i}J_{2i+1}(t)T_{2i+1}(x)}_{\textrm{sup}}&\leqslant\varepsilon,\end{split}$

where $J_{i}$ is the $i$-th Bessel function of the first kind.

######  Corollary 11.8 (Hamiltonian simulation).

Given a symmetric Hamiltonian $H\in\mathbb{C}^{n\times n}$ with
$\MT@delim@RIGHT\norm@star{H}\leqslant 1$ and a vector $b\in\mathbb{C}^{n}$,
after $O(\operatorname{nnz}(A))$ preprocessing, we can output a description of
a vector $v$ such that, with probability $\geqslant 0.9$,
$\MT@delim@RIGHT\abs@star{v-e^{iHt}b}\leqslant\varepsilon\MT@delim@RIGHT\abs@star{b}$
in $O(t^{9}\lVert
H\rVert_{F}^{4}\log(1/\varepsilon)^{9}/(\varepsilon^{2}\log\log(1/\varepsilon)^{9}))$.

###### Proof.

Let
$p_{\cos}=J_{0}(t)-2\sum_{i\in[1,r]}\MT@delim@RIGHT\Paren@star{-1}^{i}J_{2i}(t)T_{2i}(x)$
and let
$p_{\sin}=2\sum_{i\in[0,r]}\MT@delim@RIGHT\Paren@star{-1}^{i}J_{2i+1}(t)T_{2i+1}(x)$.

We apply Theorem 10.6 to get a description of $c$ such that
$\MT@delim@RIGHT\norm@star{c-p_{\cos}(tH)b}\leqslant\varepsilon\MT@delim@RIGHT\norm@star{b}$
and apply Theorem 10.2 to get a description of $s$ such that
$\MT@delim@RIGHT\norm@star{s-p_{\sin}(tH)b}\leqslant\varepsilon\MT@delim@RIGHT\norm@star{b}$.
Then

$\displaystyle
e^{iHt}b=\cos(Ht)b+i\sin(Ht)b\approx_{\varepsilon\MT@delim@RIGHT\norm@star{b}}p_{\cos}(Ht)b+ip_{\sin}(Ht)b\approx_{\varepsilon\MT@delim<EMAIL_ADDRESS>
This gives us a description $O(\varepsilon)$-close to $e^{iHt}$. Using
Corollary 5.11, we can get a sample from this output in
$O\MT@delim@RIGHT\Paren@star{t^{6}\lVert
A\rVert_{F}^{4}\log(n)\frac{1}{\varepsilon^{2}}\operatorname{polylog}(\frac{1}{\varepsilon})}$
time. ∎

## Acknowledgments

ET and AB thank Nick Trefethen, Simon Foucart, Alex Townsend, Sujit Rao, and
Victor Reis for helpful discussions. ET thanks t.f. for the support. AB is
supported by Ankur Moitra’s ONR grant. ET is supported by the NSF GRFP
(DGE-1762114). Work on this paper was initiated at the Simon Institute’s
“Probability, Geometry, and Computation in High Dimensions” program in 2020;
we thank the institute for its hospitality.

## References

  * [AA18] Scott Aaronson and Andris Ambainis.  Forrelation: A problem that optimally separates quantum from classical computing.  SIAM Journal on Computing, 47(3):982–1038, 2018.
  * [AHK06] Sanjeev Arora, Elad Hazan, and Satyen Kale.  A fast random sampling algorithm for sparsifying matrices.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 272–279. Springer, 2006.
  * [AM07] Dimitris Achlioptas and Frank McSherry.  Fast computation of low-rank matrix approximations.  Journal of the ACM (JACM), 54(2):9–es, 2007.
  * [AZL16] Zeyuan Allen-Zhu and Yuanzhi Li.  Lazysvd: Even faster svd decomposition yet without agonizing pain.  Advances in neural information processing systems, 29, 2016.
  * [BCW20] Ainesh Bakshi, Nadiia Chepurko, and David P Woodruff.  Robust and sample optimal algorithms for psd low rank approximation.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 506–516. IEEE, 2020.
  * [BCW22] Ainesh Bakshi, Kenneth L Clarkson, and David P Woodruff.  Low-rank approximation with $1/\varepsilon^{1/3}$ matrix-vector products.  pages 1130–1143, 2022.
  * [BKKS21] Vladimir Braverman, Robert Krauthgamer, Aditya R Krishnan, and Shay Sapir.  Near-optimal entrywise sampling of numerically sparse matrices.  In Conference on Learning Theory, pages 759–773. PMLR, 2021.
  * [BKL+19] Fernando G. S. L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu.  Quantum SDP solvers: Large speed-ups, optimality, and applications to quantum learning.  In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), pages 27:1–27:14, 2019.  arXiv: 1710.02581
  * [BKM22] Vladimir Braverman, Aditya Krishnan, and Christopher Musco.  Sublinear time spectral density estimation.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2022, page 1144–1157, New York, NY, USA, 2022. Association for Computing Machinery.
  * [BMN+21] Ryan Babbush, Jarrod R McClean, Michael Newman, Craig Gidney, Sergio Boixo, and Hartmut Neven.  Focus beyond quadratic speedups for error-corrected quantum advantage.  PRX Quantum, 2(1):010103, 2021.
  * [BW18] Ainesh Bakshi and David Woodruff.  Sublinear time low-rank approximation of distance matrices.  Advances in Neural Information Processing Systems, 31, 2018.
  * [BWP+17] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  Nature, 549(7671):195--202, 2017.
  * [CCH+22] Nadiia Chepurko, Kenneth Clarkson, Lior Horesh, Honghao Lin, and David Woodruff.  Quantum-inspired algorithms from randomized numerical linear algebra.  In International Conference on Machine Learning, pages 3879--3900. PMLR, 2022.
  * [CdW21] Yanlin Chen and Ronald de Wolf.  Quantum algorithms and lower bounds for linear regression with norm constraints, 2021.
  * [CEM+15] Michael B Cohen, Sam Elder, Cameron Musco, Christopher Musco, and Madalina Persu.  Dimensionality reduction for k-means clustering and low rank approximation.  In Proceedings of the forty-seventh annual ACM symposium on Theory of computing, pages 163--172, 2015.
  * [CGJ19] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation.  In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), pages 33:1--33:14, 2019.  arXiv: 1804.01973
  * [CGL+20a] Nai-Hui Chia, András Gilyén, Tongyang Li, Han-Hsuan Lin, Ewin Tang, and Chunhao Wang.  Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing - STOC 2020. ACM Press, 2020.
  * [CGL+20b] Nai-Hui Chia, András Gilyén, Han-Hsuan Lin, Seth Lloyd, Ewin Tang, and Chunhao Wang.  Quantum-inspired algorithms for solving low-rank linear equation systems with logarithmic dependence on the dimension.  In Proceedings of the 31st International Symposium on Algorithms and Computation (ISAAC), pages 47:1--47:17, 2020.
  * [CHI+18] Carlo Ciliberto, Mark Herbster, Alessandro Davide Ialongo, Massimiliano Pontil, Andrea Rocchetto, Simone Severini, and Leonard Wossnig.  Quantum machine learning: a classical perspective.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2209):20170551, 2018.
  * [CKS17] Andrew M. Childs, Robin Kothari, and Rolando D. Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920--1950, 2017.
  * [Cle55] C. W. Clenshaw.  A note on the summation of Chebyshev series.  Math. Tables Aids Comput., 9:118--120, 1955.
  * [CNW15] Michael B Cohen, Jelani Nelson, and David P Woodruff.  Optimal approximate matrix product in terms of stable rank.  arXiv preprint arXiv:1507.02268, 2015.
  * [DKM06] P. Drineas, R. Kannan, and M. Mahoney.  Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix.  SIAM Journal on Computing, 36(1):158--183, January 2006.
  * [DW20] Vedran Dunjko and Peter Wittek.  A non-review of quantum machine learning: trends and explorations.  Quantum Views, 4:32, 2020.
  * [DZ11] Petros Drineas and Anastasios Zouzias.  A note on element-wise matrix sparsification via a matrix-valued bernstein inequality.  Information Processing Letters, 111(8):385--389, 2011.
  * [FP68] L. Fox and I. B. Parker.  Chebyshev polynomials in numerical analysis.  Oxford University Press, London-New York-Toronto, Ont., 1968.
  * [GLG22] Sevag Gharibian and François Le Gall.  Dequantizing the quantum singular value transformation: Hardness and applications to quantum chemistry and the quantum pcp conjecture.  In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 19--32, 2022.
  * [GLM08] Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone.  Quantum random access memory.  Physical review letters, 100:160501, Apr 2008.
  * [GS18] Neha Gupta and Aaron Sidford.  Exploiting numerical sparsity for efficient learning: faster eigenvector computation and regression.  Advances in Neural Information Processing Systems, 31, 2018.
  * [GSLW19] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st ACM Symposium on the Theory of Computing (STOC), pages 193--204, 2019.  arXiv: 1806.01838
  * [GST22] András Gilyén, Zhao Song, and Ewin Tang.  An improved quantum-inspired algorithm for linear regression.  Quantum, 6:754, 2022.
  * [GT09] Alex Gittens and Joel A Tropp.  Error bounds for random matrix approximation schemes.  arXiv preprint arXiv:0911.4108, 2009.
  * [HHL09] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15):150502, 2009.  arXiv: 0811.3171
  * [IVWW19] Pitor Indyk, Ali Vakilian, Tal Wagner, and David P Woodruff.  Sample-optimal low-rank approximation of distance matrices.  In Conference on Learning Theory, pages 1723--1751. PMLR, 2019.
  * [JLGS20] Dhawal Jethwani, François Le Gall, and Sanjay K. Singh.  Quantum-inspired classical algorithms for singular value transformation.  In Proceedings of the 45th International Symposium on Mathematical Foundations of Computer Science (MFCS), pages 53:1--53:14, 2020.  arXiv: 1910.05699
  * [KD14] Abhisek Kundu and Petros Drineas.  A note on randomized element-wise matrix sparsification.  arXiv preprint arXiv:1404.0320, 2014.
  * [KDMI17] Abhisek Kundu, Petros Drineas, and Malik Magdon-Ismail.  Recovering pca and sparse pca via hybrid-(l1, l2) sparse sampling of data elements.  The Journal of Machine Learning Research, 18(1):2558--2591, 2017\.
  * [KLLP19] Iordanis Kerenidis, Jonas Landman, Alessandro Luongo, and Anupam Prakash.  q-means: A quantum algorithm for unsupervised machine learning.  Advances in Neural Information Processing Systems, 32, 2019.
  * [KP17] Iordanis Kerenidis and Anupam Prakash.  Quantum recommendation systems.  In Proceedings of the 8th Innovations in Theoretical Computer Science Conference (ITCS), pages 49:1--49:21, 2017.  arXiv: 1603.08675
  * [KP22] Iordanis Kerenidis and Anupam Prakash.  Quantum machine learning with subspace states.  arXiv preprint arXiv:2202.00054, 2022.
  * [Lan50] Cornelius Lanczos.  An iteration method for the solution of the eigenvalue problem of linear differential and integral operators.  J. Research Nat. Bur. Standards, 45:255--282, 1950.
  * [LC17] Guang Hao Low and Isaac L. Chuang.  Optimal Hamiltonian simulation by quantum signal processing.  Physical Review Letters, 118(1):010501, 2017.  arXiv: 1606.02685
  * [LGZ16] Seth Lloyd, Silvano Garnerone, and Paolo Zanardi.  Quantum algorithms for topological and geometric analysis of data.  Nature communications, 7(1):1--7, 2016.
  * [LMR13] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum algorithms for supervised and unsupervised machine learning.  arXiv, 2013.
  * [LMR14] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum principal component analysis.  Nature Physics, 10(9):631--633, July 2014.
  * [Mah11] Michael W. Mahoney.  Randomized algorithms for matrices and data.  Foundations and Trends® in Machine Learning, 3(2):123--224, 2011.
  * [MH02] John C Mason and David C Handscomb.  Chebyshev polynomials.  Chapman and Hall/CRC, 2002.
  * [MI11] Malik Magdon-Ismail.  Using a non-commutative bernstein bound to approximate some matrix algorithms in the spectral norm.  arXiv preprint arXiv:1103.5453, 2011.
  * [MM15] Cameron Musco and Christopher Musco.  Randomized block krylov methods for stronger and faster approximate singular value decomposition.  Advances in neural information processing systems, 28, 2015.
  * [MMS18] Cameron Musco, Christopher Musco, and Aaron Sidford.  Stability of the lanczos method for matrix function approximation.  In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1605--1624. Society for Industrial and Applied Mathematics, January 2018.
  * [MRTC21] John M. Martyn, Zane M. Rossi, Andrew K. Tan, and Isaac L. Chuang.  Grand unification of quantum algorithms.  Physical Review X, 2(4):040203, 2021.  arXiv: 2105.02859
  * [MW17] Cameron Musco and David P Woodruff.  Sublinear time low-rank approximation of positive semidefinite matrices.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 672--683. IEEE, 2017.
  * [MZ11] Avner Magen and Anastasios Zouzias.  Low rank matrix-valued chernoff bounds and approximate matrix multiplication.  In Proceedings of the twenty-second annual ACM-SIAM symposium on Discrete Algorithms, pages 1422--1436. SIAM, 2011.
  * [NN13] Jelani Nelson and Huy L. Nguyen.  OSNAP: Faster numerical linear algebra algorithms via sparser subspace embeddings.  In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science. IEEE, October 2013.
  * [Oli77] J. Oliver.  An error analysis of the modified Clenshaw method for evaluating Chebyshev and Fourier series.  J. Inst. Math. Appl., 20(3):379--391, 1977.
  * [Oli79] J. Oliver.  Rounding error propagation in polynomial evaluation schemes.  J. Comput. Appl. Math., 5(2):85--97, 1979.
  * [Pai76] C. C. Paige.  Error analysis of the lanczos algorithm for tridiagonalizing a symmetric matrix.  IMA Journal of Applied Mathematics, 18(3):341--349, 1976.
  * [RML14] Patrick Rebentrost, Masoud Mohseni, and Seth Lloyd.  Quantum support vector machine for big data classification.  Physical review letters, 113:130503, September 2014.
  * [RWC+20] Alessandro Rudi, Leonard Wossnig, Carlo Ciliberto, Andrea Rocchetto, Massimiliano Pontil, and Simone Severini.  Approximating hamiltonian dynamics with the nyström method.  Quantum, 4:234, 2020.
  * [Saa81] Yousef Saad.  Krylov subspace methods for solving large unsymmetric linear systems.  Mathematics of computation, 37(155):105--126, 1981.
  * [Sch41] A. C. Schaeffer.  Inequalities of A. Markoff and S. Bernstein for polynomials and related functions.  Bull. Amer. Math. Soc., 47:565--579, 1941.
  * [Sha04] Aleksei Shadrin.  Twelve proofs of the Markov inequality.  In Approximation theory: a volume dedicated to Borislav Bojanov, pages 233--298. Prof. M. Drinov Acad. Publ. House, Sofia, 2004.
  * [SM21] Changpeng Shao and Ashley Montanaro.  Faster quantum-inspired algorithms for solving linear systems.  arXiv, 2021.
  * [SV+14] Sushant Sachdeva, Nisheeth K Vishnoi, et al.  Faster algorithms via approximation theory.  Foundations and Trends® in Theoretical Computer Science, 9(2):125--210, 2014.
  * [Tan19] Ewin Tang.  A quantum-inspired classical algorithm for recommendation systems.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing - STOC 2019. ACM Press, 2019.
  * [Tan21] Ewin Tang.  Quantum principal component analysis only achieves an exponential speedup because of its state preparation assumptions.  Phys. Rev. Lett., 127:060503, Aug 2021.
  * [Tao13] T Tao.  Matrix identities as derivatives of determinant identities, 2013.  URL https://terrytao. wordpress. com/2013/01/13/matrix-identities-as-derivatives-of-determinant-identities, 2013\.
  * [TBI97] Lloyd N Trefethen and David Bau III.  Numerical linear algebra, volume 50.  Siam, 1997.
  * [Tre19] Lloyd N. Trefethen.  Approximation Theory and Approximation Practice, Extended Edition.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2019\.
  * [Tro12] Joel A Tropp.  User-friendly tail bounds for sums of random matrices.  Foundations of computational mathematics, 12(4):389--434, 2012.
  * [VdN11] Maarten Van den Nest.  Simulating quantum computers with probabilistic methods.  Quantum Information and Computation, 11(9&10):784--812, 2011.  arXiv: 0911.1624
  * [Vos91] Michael D Vose.  A linear algorithm for generating random numbers with a given distribution.  IEEE Transactions on software engineering, 17(9):972--975, 1991\.
  * [WKS14] Nathan Wiebe, Ashish Kapoor, and Krysta M Svore.  Quantum deep learning.  arXiv preprint arXiv:1412.3489, 2014.
  * [WKS16] Nathan Wiebe, Ashish Kapoor, and Krysta M. Svore.  Quantum perceptron models.  In Advances in Neural Information Processing Systems 29 (NIPS), pages 3999--4007, 2016.  arXiv: 1602.04799
  * [Woo14] David P. Woodruff.  Sketching as a tool for numerical linear algebra.  Foundations and Trends® in Theoretical Computer Science, 10, 2014\.
  * [ZFF19] Zhikuan Zhao, Jack K. Fitzsimons, and Joseph F. Fitzsimons.  Quantum-assisted Gaussian process regression.  Physical Review A, 99(5):052331, 2019.  arXiv: 1512.03929


11institutetext: Institut UTINAM, UMR 6213 CNRS - Univ. Bourgogne Franche-
Comté, OSU THETA, BP 1615, F-25010 Besançon Cedex, France
11email<EMAIL_ADDRESS>22institutetext: Aix Marseille Univ,
Institut Origines, CNRS, CNES, LAM, Marseille, France

# Volatile-rich comets ejected early on during Solar System formation

S.E. Anderson 11 J.-M. Petit 11 B. Noyelles 11 O. Mousis 22 P. Rousselot 11

(Accepted September 8, 2022)

###### Abstract

Context. Comet C/2016 R2 PanSTARRS (hereafter C/2016 R2) presents an unusually
high N2/CO abundance ratio, as well as a heavy depletion in H2O, making it the
only known comet of its kind. Understanding its dynamical history is therefore
of essential importance as it would allow us to gain a clearer understanding
of the evolution of planetesimal formation in our Solar System. Two studies
have independently estimated the possible origin of this comet from building
blocks formed in a peculiar region of the protoplanetary disk, near the ice
line of CO and N2.

Aims. We intend to investigate the fates of objects formed from the building
blocks in these regions. We hope to find a possible explanation for the lack
of C/2016 R2-like comets in our Solar System.

Methods. Using a numerical simulation of the early stages of Solar System
formation, we track the dynamics of these objects in the Jumping Neptune
scenario based on five different initial conditions for the protosolar disk.
We integrate the positions of 250 000 planetesimals over time in order to
analyze the evolution of their orbits and create a statistical profile of
their expected permanent orbit.

Results. We find that objects formed in the region of the CO- and N2\- ice
lines are highly likely to be sent towards the Oort Cloud or possibly ejected
from the Solar System altogether on a relatively short timescale. In all our
simulations, over 90% of clones formed in this region evolved into a
hyperbolic trajectory, and between 1% and 10% were potentially captured by the
Oort Cloud. The handful of comets that remained were either on long-period,
highly eccentric orbits like C/2016 R2, or absorbed into the Edgeworth–Kuiper
belt.

Conclusions. Comets formed ¡15 au were predominantly ejected early in the
formation timeline. As this is the formation zone likely to produce comets of
this composition, this process could explain the lack of similar comets
observed in the Solar System.

###### Key Words.:

Comets:general – Comets:individual:C/2016R2(PanSTARRS)– Protoplanetary Disks

##  1 Introduction

Comets are some of the most pristine bodies in the Solar System, having
remained relatively unchanged since their formation 4.6 billion years ago.
Cometary nuclei provide insights into the composition of the early
protoplanetary disk (PPD) through their isotopic abundance ratios. As their
composition reflects the physico-chemical conditions of the disk at the
location of their formation in the protosolar nebula (PSN), understanding
where each comet was formed reveals details as to the evolution of the Solar
System.

Decades of remote sensing of comets have revealed these objects to be water-
ice rich, with a typical carbon monoxide composition of CO/H2O = 4% (Bockelée-
Morvan & Biver, 2017), and depleted in N2 despite the abundance of this
molecule in the atmospheres and surfaces of the outer Solar System bodies,
such as Triton or Pluto (Cochran et al., 2000). However, radio observations of
the long-period comet C/2016 R2 (PanSTARRS) revealed that its composition is
unlike any comet observed before, with the spectrum dominated by bands of CO+.
This CO-rich comet is remarkably depleted in water, with a H2O/CO ratio of
only $\sim 0.32$ % (McKay et al., 2019) with an upper limit of H2O/CO ¡ 0.1
(Biver et al., 2018). Further, it has a peculiar abundance of N${}_{2}^{+}$,
with N2/CO estimated to be between 0.05 $\pm$ 0.01 (McKay et al., 2019), 0.06
$\pm$ 0.01 (Opitom et al., 2019), and 0.08 $\pm$ 0.02 (Biver et al., 2018),
which had never been seen in such high quantities in comets before. This
composition changes our perception of comet formation, as it was previously
understood that CO ice is unlikely freeze out without abundant water ice,
which has a higher binding energy than CO (Boogert et al., 2015). Most
volatile species would also be expected to deplete with each subsequent
passage of this comet within the inner Solar System. Understanding the
dynamical history of this comet is therefore of essential importance to
understanding the timeline of planetesimal formation in our Solar System.

Other potential N2-presenting candidates have been identified, such as C/1908
R1 (Morehouse), C/1961 R1 (Humason), C/1987 P1 (Bradfield), C/2001 Q4 (NEAT)
with N2/CO=0.027 (Feldman, 2015), and C/2002 VQ94 (LINEAR) for which
N2/CO=0.06 (Korsun et al., 2008). A few short-period comets also show an
increased N2/CO ratio, such as comet 29P/Schwassmann-Wachmann 1 with
N2/CO=0.013 (Ivanova et al., 2016), or comet 67P presenting N2/CO=0.0287, but
this result came from in situ measurements (Rubin et al., 2020). Some others
present moderately unusual water-poor compositions; for example, interstellar
comet 2I/Borisov, measured to have CO/H2O between 35% and 173% (Cordiner et
al., 2020; Bodewits et al., 2020), which is significantly higher than the
average cometary values for our Solar System, and could be explained by an
unusual formation environment beyond the CO snow line of its own system (Price
et al., 2021). Comet C/2009 P1 (Garradd) is another outlier with a CO
production rate of 63% of that of water, yet no N2 was detected (Feaga et al.,
2014). This simultaneously CO- and N2-rich and water-poor composition, along
with none of the usual neutrals seen in most cometary spectra, makes C/2016 R2
a unique and intriguing specimen, the only one of its kind to ever be
observed.

Such a small sample size makes it impossible to draw conclusions as to a
shared formation reservoir. The long-period comets share highly eccentric,
almost parabolic orbits —even hyperbolic in the case of C/1908 R1 (Morehouse)
and C/2001 Q4 (NEAT)—, while Comet 29P/Schwassmann-Wachmann 1 is likely a
captured Oort Cloud object (Neslušan et al., 2017). It is clear these objects
must have spent the majority of their lifetime at high heliocentric distance,
else they would have already lost their volatile content. Unfortunately,
attempts to trace back their dynamical history with any degree of certainty is
made impossible by the inherent chaotic nature of their motion due to frequent
close encounters with the gas giants, which strip comets of their dynamical
memory. As a result, there is no peculiar C/2016 R2-like orbit despite its
otherwise peculiar nature, and we cannot trace its dynamical history backwards
to a potential shared formation reservoir.

The origin of the unusual composition of C/2016 R2 is highly disputed. It may
be a fragment of a differentiated object as suggested by Biver et al. (2018),
similar to the CO-rich interstellar comet 2I/Borisov (Cordiner et al., 2020).
If CO is absent in the upper layers of an as-of-yet undiscovered
differentiated comet, as suggested by De Sanctis et al. (2001), then it is
possible that C/2016 R2 is a fragment of the core of such a comet. Desch &
Jackson (2021) theorize that 1I/’Oumuamua may be an N2 iceberg chipped off
from the surface of an ex-Pluto by an impact during a period of dynamical
instability, which could be applied to C/2016 R2. Another possibility is that
the particular composition of C/2016 R2 simply arises from where it formed in
the PSN: Perhaps this disk could evolve over time to create exotic
compositions at different disk locations in unique proportions, in ”special”
comet-forming annuli. Two studies independently estimated the possible origin
of this comet from building blocks formed in a peculiar region of the PSN,
near the ice lines of CO and N2. By evaluating the radial transport of
volatiles in the PSN, Mousis et al. (2021) found that the peculiar N2/CO ratio
of C/2016 R2 could be replicated by agglomeration from particles near the N2
and CO ice lines, within the 10-15 au region. Meanwhile, the CO/H2O ratio
would remain deeply depleted inward of the CO ice line, around the 8-11 au
region. Cold traps of hypervolatiles in the PSN in a small, specific region of
the disk could explain the peculiar composition of this latter comet.
Similarly, Price et al. (2021) model the effect of drifting solid material in
the PPD and find that the ideal location for the formation of CO-rich,
H2O-poor objects is beyond CO ice line. However, this would seem to indicate
that more CO-rich comets should exist than have previously been observed. The
N2/CO ratio was not a part of their study.

Here we explore the potential fates of comets formed from these building
blocks using a numerical simulation of early Solar System formation. By
examining the dynamical evolution of only the objects formed in a small exotic
pocket, or ”Sweet Spot,” of the PSN, which allows for peculiar-composition
comets to form, we hope to understand why so few are observed today. In
Section 2, we describe the model we use to simulate the early Solar System and
the dynamical evolution of these small bodies. In Section 3, we report on
these results and examine more closely the fates of all comets, then narrow
our interest to comets that would have formed inside the Sweet Spot. Finally,
in Section 4, we provide our conclusions as to what these fates will be.

##  2 Methods

We employ the Jumping Neptune scenario from Nesvorný (2015). We begin with
five planets: Jupiter, Saturn, and three ice giants of comparable mass, as
described by Deienno et al. (2017). This third ice giant, henceforth I1,
undergoes a series of encounters with Jupiter and Saturn which causes a
divergent jump in their semi-major axes before inducing a jump in Neptune’s
orbit as well. Finally, I1 is ejected onto a hyperbolic orbit, leaving the
remaining four planets near their present-day orbits. We tested several
alternative simulations, varying the multi-resonance configuration, the
distance from the last planet to the inner edge of the disk (1 or 2 au), the
mass of the disk (20 or 40 $M_{\Earth}$), and the inclination of the disk in
relation to the plane of the planets, with five different evolutions in each
case. We selected the simulations that best satisfy the criteria of similarity
with the Solar System today, consistent with the current orbital structure of
the trans-Neptunian population, in line with Deienno et al. (2017), which were
all from the 3:2, 3:2, 2:1, 3:2 multi-resonance configuration, with a disk of
40 $M_{\Earth}$. The initial multi-resonant configurations we choose for
Jupiter, Saturn, Uranus, Neptune, and I1, along with parameters for the
location and mass of disk, begin in a 3:2, 3:2, 2:1, 3:2 resonance, as Baguet
et al. (2019) find this is able to place a secular tilt resonance in the area
of the cold Edgeworth-Kuiper belt (between 39 and 48 au). This provides us
with five scenarios to explore, as defined in Table 1. All these
configurations require the existence of the fifth giant planet, with a mass
comparable to those of Uranus or Neptune, which is eventually ejected during
the instability. This ice planet would have formed within the volatile-rich
zone identified by Mousis et al. (2021). All the planetary evolution
simulations were run self-consistently with the five planets and a swarm of
1000 massive particles of the same mass, each 1/1000th of the mass of the
disk. The disk extends from its inner edge to 30 au.

Figure 1: Dynamical evolution of the planets and 1000 comet clones in Scenario 1 with eccentricity as a function of semi-major axis in log scale. Planets are represented in black, with Jupiter, Saturn, I1, Uranus, and Neptune, respectively, from left to right. Blue indicates comets formed between 4 and 8 au. Turquoise indicates comets formed between 8 and 12 au. Green indicates comets formed between 12 and 20 au. Yellow indicates comets formed between 20 and 50 au. After 100 Myr, the area around the giant planets is entirely cleared. Comets formed between 4 and 12 au are the first to be lost, and by 10 Myr almost none remain. By 100 Myr, the comets that remain in our simulation are almost entirely from the 20-50 au population. Figure 2: Final positions in log scale of the first 1000 clones from Scenario 1 color coded for the last recorded time before they are removed from the simulation (see main text). Violet clones were lost early in the simulation, while yellow clones indicate those that remain. The black line corresponds to an equal initial and final position: comets above would have moved away from their initial formation position, and comets below would have migrated inwards. The area near the orbits of Jupiter and Saturn clears quickly, sending the clones on highly eccentric orbits even before Neptune’s migration occurs. We see the relatively stable location of the current day Edgeworth–Kuiper Belt. Table 1: Initial conditions for the five scenarios explored in this study. The multi-resonance configuration is the 3:2, 3:2, 2:1, 3:2, with the outermost planet at 20.18 au from the Sun. In all cases studied here, we used a disk of 40 $M_{\Earth}$. (a) Distance of the inner bound of the disk (au), (b) inclination of the disk with respect to the invariable plane (∘), (c) node of the disk (∘), and (d) running number, i.e., number of generations used. | (a) | (b) | (c) | (d)
---|---|---|---|---
Scenario 1 | 1 | 1 | 100.2 | 08
Scenario 2 | 1 | 0 | 0 | 02
Scenario 3 | 2 | 0 | 0 | 04
Scenario 4 | 2 | 1 | 100.2 | 01
Scenario 5 | 2 | 1 | 100.2 | 05

The present-day Edgeworth–Kuiper belt extends from the orbit of Neptune at 30
au to approximately 50 au from the Sun. However, most of the small bodies of
the outer Solar System originated from the region between Jupiter and $\sim$30
au (Gomes, 2003; Tsiganis et al., 2005; Levison et al., 2008; Kaib & Quinn,
2008). With this in mind, we limit our simulations to planetesimals formed in
the 4 - 50 au range. This allows us to neglect the influence of the inner
planets, which, having small orbits, require more integration steps and longer
calculations on each of our clones. While the CO-rich comet-forming zone could
extend to 100 au (Price et al., 2021), the mass depletion of the classical
belt is already well explored.

We then run a modified SWIFT numerical integrator which uses a pre-recorded
evolution of the giant planets (Petit et al., 1999) and evolves our system
over 100 Myr. The previously calculated evolution of the planets is recorded
every 1000 yr or less and the positions of the planets are interpolated at
each time-step necessary for the integration of the motion of the test
particles. This ensures that each simulation for a given planetary evolution
will use exactly the same planetary evolution track, avoiding divergence due
to the intrinsic chaotic nature of planetary motion. Thus, our final planetary
system is sure to correctly reproduce the structure of the Solar System. The
major difference in planetary behavior between these scenarios is the moment
of ejection of I1. This occurs at 5 Myr, 6 Myr, 7 Myr, 8 Myr, and 13 Myr for
scenario 1, 2, 3, 4, and 5, respectively.

For each scenario, we run 50 sets of 1000 massless comet facsimiles or
”clones.” Each clone has randomly generated orbital elements setting them on
the same plane as the disk with varying semi-major axes between 4 au —to avoid
the inner Solar System— and 50 au. The clones are distributed with a number
density that varies as $r^{-1/2}$, or a surface density that varies as
$r^{-3/2}$. We therefore have a total of 250 000 clones for our five
scenarios. Our simulations count a clone as lost if it reaches beyond 10000 au
as we do not yet have the ability to estimate the effects of the Galactic
tidal forces. If a clone moves under 0.005 au from the Sun, or in collision
with a planet, it is also removed from the integration, as it is most likely
destroyed.

##  3 Results and Discussion

We examine the orbital elements of each clone, identified by its formation
location (initial semi-major axis). This is shown in Fig. 1. Within the first
1 Myr, 21% of all clones are lost from the simulation. This number rises to
nearly half (49%) after 10 Myr. By the end of the 100 Myr simulation run, we
have lost three quarters (76%) of our initial population. Only a quarter (24%)
of our clones remain. A snapshot of the first 1000 clones in our first
scenario is shown in Fig. 2 with their first and final positions. Each clone
is color coded for the moment it is lost, with earlier losses in purple and
those that remain in the end shown in yellow. The major loss of clones occurs
before $\sim$10 Myr: after this time, the area around the giant planets (¡ 15
au) is entirely cleared. It is important to note here that if we had used a
four-planet model, based on current planetary orbits, Saturn would play the
role of I1 and clear this region, leading to the same outcome.

We examine the percentage of clones lost in our simulations more closely for
each 1 au annulus from 4 au to 50 au, as shown in Fig. 3. We see that for
every 1 au annulus between 4 and 10 au, over 95% of the clones are lost before
the end of the 100 Myr in each scenario. This number dips to 90% around 12 au.
Then, between 12 and 20 au, each scenario still loses a minimum of 80% of
their clones within the simulation time. In comparison, annuli beyond 40 au
—the current location of the Classical Edgeworth–Kuiper Belt— only lose half
their clones, showing a zone that is relatively stable, containing objects
that do not move far from where they are formed. The behavior of these clones
is consistent between scenarios and independent of the moment of ejection of
I1.

Figure 3: Percentage of clones lost per formation location for each of the
five scenarios. The gray zone indicates the limitation of our simulation. The
blue zone indicates the N2/CO enrichment zone as predicted by Mousis et al.
(2021), while the overlaid green zone indicates the location of the ideal
CO/H2O enrichment zone.

Based on the ranges proposed by Price et al. (2021) and Mousis et al. (2021),
we examine different formation zones. Price et al. (2021) suggest a wide
range, arguing that the CO/H2O ice-enrichment zone is likely between 20 and
100 au, though these authors do not investigate a N2/CO ice-enrichment zone.
As the CO/H2O ice enrichment zone evolves over time, and without seeing how N2
would evolve in the simulations of Price et al. (2021), we cannot determine
where a specific C/2016 R2 formation zone could occur. In light of the fact
that both CO and N2 have similar sublimation temperatures, the two ice lines
should be near each other and make the 20-30 au annulus an area to explore.
Meanwhile, the results of Mousis et al. (2021) would indicate a narrow area,
as they find a CO/H2O ice-enrichment zone of $\sim$1-2 au wide, near 10 au.
Their N2/CO ice-enrichment zone is narrower still, seemingly less than 1 au.
The overlapping formation zone for a C/2016 R2-like comet would therefore be
incredibly narrow. We examine both a wide C/2016 R2 forming annulus between 8
and 20 au; a narrow one, only 8-10 au; and the narrowest one between 10 and 11
au, that is, the Sweet Spot. Interestingly, I1 is initialized and subsequently
ejected from this narrow region as well.

Figure 4: Final semi-major axes and eccentricities (left) and perihelion
distances (right) of all clones from all simulations remaining after 100 Myr.
Comets formed in the 8-11 au zone are shown in black. Any remaining object
with a perihelion $q<35$ au will likely be sent to the Oort Cloud by Neptune,
or will lose its hypervolatile majority ices to vacuum via insolation heating
Lisse et al. (2022).

The resulting statistics are shown in Table 2. On average, each simulation
loses 75% of its clones by 100 Myr, losing 90% of all clones formed in the
8-20 au range, 97% of all clones formed in the 8-10 au range, and $\sim 80\%$
in the 20-30 range. Consistently, each simulation loses 96% of all comets
formed between 10 and 11 au. If we examine the region of clones initialized
between 8 and 20 au, we find that half the clones are already lost by 5 Myr,
with two-thirds of clones being ejected after 15 Myr. If we narrow that region
further to 8-11 au, we find that 60% of clones formed in this region are
ejected in the first 1 Myr and 90% after 10 Myr. A handful of clones
($\sim$0.1%) are lost to collisions with the giant planets. Depending on the
chronology, these could help account for the delivery of the building blocks
of the Galilean and Saturnian satellites necessary for their formation (Ronnet
et al., 2018; Anderson et al., 2021). In each simulation, irrespective of the
scenario, only $\sim$1% of all remaining clones were from the initial 8-11 au
population; these are shown in black in Fig. 4. They will either find
themselves on highly eccentric orbits, be absorbed into the Edgeworth–Kuiper
belt, or join the scattered disk. These clones seem to be evenly distributed
within the population of remaining comets. The 10-11 au population makes up
only 0.4% of all surviving clones.

Table 2: Statistical loss outcomes of each of the scenarios after 100 Myr for each formation zone. | Total loss | 8-10 au | 10-11 au | 8-20 au | 20-30 au
---|---|---|---|---|---
S1 | 73$\%$ | 97$\%$ | 96$\%$ | 88$\%$ | 76$\%$
S2 | 79$\%$ | 97$\%$ | 96$\%$ | 90$\%$ | 80$\%$
S3 | 75$\%$ | 98$\%$ | 96$\%$ | 91$\%$ | 79$\%$
S4 | 77$\%$ | 97$\%$ | 96$\%$ | 90$\%$ | 79$\%$
S5 | 73$\%$ | 97$\%$ | 96$\%$ | 89$\%$ | 76$\%$

We must now estimate how many C/2016 R2-like comets could be captured by the
Oort Cloud, so as to then evolve dynamically over the next 4 Gyr and return to
visit the inner Solar System on C/2016 R2-like orbits. While it is tempting to
say that the cometesimals lost from our simulation were ejected from our Solar
System, a further investigation of the orbital elements at the moment they
were removed from the simulation is required in order to estimate their
capture rate by the Oort Cloud. This rate is poorly constrained as of yet as
this would depend greatly on the timeline of evolution coinciding with our
Sun’s ejection from its parent cluster. Further numerical simulations are
required in order to investigate the behavior of these comets beyond the 10000
au cutoff, although the effects of Galactic tides 4 Gyr ago are still unknown.
Nevertheless, we can make a safe estimate of which comets are bound to the
Solar System from the energy $z$ of the cometesimals at the moment they are
lost:

$z=i_{\alpha}\frac{GM_{\Sun}}{2a},$ (1)

where $i_{\alpha}$ is 1, 0, and -1 for e¿1, e=1, and e¡1, respectively,
$M_{\Sun}$ the mass of the Sun, and $G$ the gravitational constant. We
consider an object captured by the Oort Cloud if the final semi-major axis is
$a>10000$ au and its final energy is $z/GM_{\Sun}\leq 0.00005$ at the moment
it is lost from our simulation. Otherwise, we consider it has a truly
hyperbolic orbit and is seen as detached from the Solar System.

Under these conditions, 11% of all our clones have potentially reached the
Oort Cloud. The distribution of contributions to the Oort Cloud from each
formation zone can be seen in Fig. 5. Of the cometesimals formed between 8 and
11 au, 13% may have reached the Oort Cloud, which represents 12% of all the
cometesimals potentially captured. The 20-30 au formation zone contributes 12%
(making up 22% of the total number of those captured); the 30-40 au formation
zone contributes 9% (making up 14% of the total number of those captured); and
the 40-50 au formation zone contributes only 5% (making up 6% of the total
number of those captured). If this estimate is accurate, we should have far
more CO-rich comets in the Oort Cloud, between 10% and 20% of all long-period
comets we observe today. However, this is not the case.

As the timeline of evolution coinciding with our Sun’s ejection from its
parent cluster is still poorly constrained, we re-estimate our capture
conditions. Portegies Zwart et al. (2021) estimate that, so long as the Sun is
a cluster member, clones with an eccentricity of $e>0.98$ and a semi-major
axis of $a>2400$ au would be vulnerable to being stripped by the cluster
potential or by passing stars. This could not apply to the entire course of
our simulation as the Oort Cloud would be unable to form if this were the
case. If we align ourselves with Portegies Zwart et al. (2021) and consider
that our Sun is still within its parent cluster for the first 10 Myr of our
simulation, we consider a clone captured if it fulfills the first criteria
(final $a>10000$ au final $z/GM_{\Sun}\leq 0.00005$) along with a new
criterion, $e<0.98$ if this event occurs within the first 10 Myr of the
simulation. We have drastically different results, as seen in Fig. 5. Under
these conditions, of the cometesimals formed between 8 and 11 au, only 1% may
have reached the Oort Cloud, which represents 4% of all the cometesimals
potentially captured. The 20-30 au formation zone contributes 7% (making up
32% of the total number of those captured); the 30-40 au formation zone
contributes 8% (making up 33% of the total number of those captured); and the
40-50 au formation zone contributes only 5% (making up 15% of the total number
of those captured). These results are coherent with our current understanding
of the chronology of Oort Cloud formation (Portegies Zwart et al., 2021),
whereby it is estimated that the bulk or 70% of the Oort Cloud material
originates from the 15-40 au region, and are near what we observe today.

Figure 5: Probability of the fates of cometesimals from each formation zone at
the end of 100 Myr if the sun has left its parent cluster (upper pannel) and
if the Sun is still in its parent cluster for the first 10 Myr of our
simulation (lower pannel). In orange we see the cometesimals ejected from the
Solar System, in blue those captured by the Oort Cloud, and in green those
remaining in the simulation. If the sun has left its parent cluster before the
beginning of the simulation, each region ¡ 35 au would contribute $\sim$ 11%
of its clones to the Oort Cloud. If the Sun is still in its parent cluster for
the first 10 Myr of the simulation, only 1% of clones formed between 8 and 11
au are likely to have reached the Oort Cloud, representing 4% of all the
cometesimals potentially captured.

##  4 Conclusions

We find that the majority of objects formed between Saturn and the N2 ice line
are ejected early and rapidly in the simulation, meaning that even by the time
the Jumping Neptune scenario occurs, the clones are already lost. This could
explain the lack of N2-rich, CO-rich, and H2O-depleted comets: these were
formed in a very narrow region, and that region empties rapidly because of the
influence of giant planets. Portegies Zwart et al. (2021) call this procession
of comet ejections from the 5-11 au zone the ‘Conveyor Belt’, which aptly
describes the phenomenon we see here. Objects formed in this region would be
ejected early from the Solar System, in less than $\sim$10 Myr, and be
unlikely to join the Oort Cloud. Therefore, if N2-rich, CO-rich, and
H2O-depleted comets were to have formed under 11 au, $>$90% of this population
would have been ejected from the Solar System without having been captured by
the Oort Cloud. As the Price et al. (2021) formation zone for CO-enrichment is
tens of astronomical units in width, suggesting that nearly half of all
observable comets ( 40%) will be CO-rich, the Oort Cloud should be full of
comets of this type. As this is not the case, then we should rule out this
model in favor of that of Mousis et al. (2021). It should also be noted that
if this mechanism is indeed the one by which the comets formed and were
subsequently captured, then C/2016 R2-like comets may be some of the first
long-period comets to have formed, and the earliest to reach the Oort Cloud.
There, objects with a nucleus larger than 5 km could survive thousands of
orbits, and hypervolatile loss upon possible perihelion passage. This would
indicate that C/2016 R2 represents one of the first Oort Cloud objects, which
could provide a direct measurement of CO/N2/CH4 ratios in the PSN (Steckloff
et al., 2021; Davidsson, 2021; Prialnik, 2021; Lisse et al., 2022), and would
explain why so few N2-rich, CO-rich, and H2O-depleted comets have been
observed today.

We must also consider the possibility that many of the cometesimals remaining
may have lost their bulk hypervolatile species in the billions of years since
their formation, or even within the time frame of our simulation. Pure
hypervolatile ices are only stable on gigayear timescales beyond a
heliocentric distance of 100 au (Lisse et al., 2021). If this ejection period
were to take place at the same time as the sublimative period of the
Edgeworth–Kuiper Belt (Lisse et al., 2021; Steckloff et al., 2021), then they
would only have $\sim$20 Myr to be placed on a trajectory toward the Oort
Cloud before they lose their mostly hypervolatile ices to vacuum via
insolation heating (Lisse et al., 2022). When looking at both the sublimation
chronology and the Oort Cloud formation chronology together, we have a small
window of only $\sim$10 Myr in which an object could be ejected from the giant
planet region and inserted into the Oort Cloud. Our results are in line with
the hypothesis of Lisse et al. (2022), which states that interstellar object
2I/Borisov was ejected early from its parent system. However, these
chronologies (the Sun’s ejection from its parent cluster; the sublimative
period of the Edgeworth–Kuiper Belt; the Jumping Jupiter/Neptune scenario) are
still poorly constrained, even more so when considering how and when these
timelines overlap. The window for the formation of C/2016 R2 in the Conveyor
Belt region and subsequent capture by the Oort Cloud could be longer or
shorter than our 10 Myr estimate.

We understand that a more quantitative simulation should account for the
effect of Galactic tides and perturbations from passing stars from the Sun’s
birth cluster. However, such an endeavor is far beyond the scope of the
present work. The mass distribution of the Galaxy and the position of the Sun
in it 4 billion years ago is still unknown, as are the initial mass function
and orbital distribution of the Sun’s birth cluster. The dynamics of the Solar
System, in particular in its infancy, are inherently chaotic. Therefore, one
would need to run a large number of simulations, varying the Galactic
potential and the Sun’s birth cluster influence, in order to only get an
average statistics of possible events, given that it is unclear whether our
Solar system is generic or peculiar.

Another possible explanation for the composition of C/2016 R2 was presented by
Desch & Jackson (2021), who suggest that C/2016 R2 could be a fragment of a
differentiated KBO surface that was created from an impact during the period
of energetic impacts during the 2:1 Jupiter:Saturn resonance epoch. If the CO-
rich formation zone is further out, as suggested by Price et al. (2021), then
the objects formed there would have longer to form, impact, and to travel to
the Oort Cloud. Further studies of the chronology of the formation of these
objects and the time frame for the dynamic instability should be explored in
order to investigate the likelihood of these scenarios. A geo-chemical study
of how these objects form accompanied with a detailed isotopic and chemical
analysis of their current composition would also be beneficial. Understanding
whether a C/2016 R2-like object can form with its peculiar composition in situ
in the disk or the composition arises from the differentiation of a Pluto-like
object would shed light on which of these processes is more likely.

This also allows for the existence of possible exotic comets, with peculiar
enrichments stemming from unique composition pockets in the disk.
Hypothetically, the ice line of each species would create a small enrichment
zone, producing small bodies dominated by this species rather than H2O. By
examining the ice lines of the volatile molecules, we can estimate the
probabilities of finding comets with each composition.

###### Acknowledgements.

We acknowledge the region of Bourgogne-Franche-Comté for their funding of the
DIAZOTE project and this work. The project leading to this publication has
received funding from the Excellence Initiative of Aix-Marseille Université -
A*Midex, a French “Investissements d’Avenir programme” AMX-21-IET-018.

## References

  * Anderson et al. (2021) Anderson, S. E., Mousis, O., & Ronnet, T. 2021, The Planetary Science Journal, 2, 50
  * Baguet et al. (2019) Baguet, D., Morbidelli, A., & Petit, J. M. 2019, Icarus, 334, 99
  * Biver et al. (2018) Biver, N., Bockelée-Morvan, D., Paubert, G., et al. 2018, Astronomy & Astrophysics, 619, A127
  * Bockelée-Morvan & Biver (2017) Bockelée-Morvan, D. & Biver, N. 2017, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 375, 20160252
  * Bodewits et al. (2020) Bodewits, D., Noonan, J. W., Feldman, P. D., et al. 2020, Nature Astronomy, 4, 867
  * Boogert et al. (2015) Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015, Annual Review of Astronomy and Astrophysics, 53, 541
  * Cochran et al. (2000) Cochran, A. L., Cochran, W. D., & Barker, E. S. 2000, Icarus, 146, 583
  * Cordiner et al. (2020) Cordiner, M. A., Milam, S. N., Biver, N., et al. 2020, Nature Astronomy, 4, 861
  * Davidsson (2021) Davidsson, B. J. R. 2021, MNRAS, 505, 5654
  * De Sanctis et al. (2001) De Sanctis, M. C., Capria, M. T., & Coradini, A. 2001, AJ, 121, 2792
  * Deienno et al. (2017) Deienno, R., Morbidelli, A., Gomes, R. S., & Nesvorný, D. 2017, AJ, 153, 153
  * Desch & Jackson (2021) Desch, S. J. & Jackson, A. P. 2021, Journal of Geophysical Research (Planets), 126, e06807
  * Feaga et al. (2014) Feaga, L. M., A’Hearn, M. F., Farnham, T. L., et al. 2014, AJ, 147, 24
  * Feldman (2015) Feldman, P. D. 2015, ApJ, 812, 115
  * Gomes (2003) Gomes, R. S. 2003, Icarus, 161, 404
  * Ivanova et al. (2016) Ivanova, O. V., Luk‘yanyk, I. V., Kiselev, N. N., et al. 2016, Planet. Space Sci., 121, 10
  * Kaib & Quinn (2008) Kaib, N. A. & Quinn, T. 2008, Icarus, 197, 221
  * Korsun et al. (2008) Korsun, P. P., Ivanova, O. V., & Afanasiev, V. L. 2008, Icarus, 198, 465
  * Levison et al. (2008) Levison, H. F., Morbidelli, A., Vanlaerhoven, C., Gomes, R., & Tsiganis, K. 2008, Icarus, 196, 258
  * Lisse et al. (2022) Lisse, C. M., Gladstone, G. R., Steckloff, J. K., et al. 2022, in LPI Contributions, Vol. 2678, LPI Contributions, 2045
  * Lisse et al. (2021) Lisse, C. M., Young, L. A., Cruikshank, D. P., et al. 2021, Icarus, 356, 114072
  * McKay et al. (2019) McKay, A. J., DiSanti, M. A., Kelley, M. S. P., et al. 2019, AJ, 158, 128
  * Mousis et al. (2021) Mousis, O., Aguichine, A., Bouquet, A., et al. 2021, arXiv e-prints, arXiv:2103.01793
  * Neslušan et al. (2017) Neslušan, L., Tomko, D., & Ivanova, O. 2017, Contributions of the Astronomical Observatory Skalnate Pleso, 47, 7
  * Nesvorný (2015) Nesvorný, D. 2015, The Astronomical Journal, 150, 73
  * Opitom et al. (2019) Opitom, C., Hutsemékers, D., Jehin, E., et al. 2019, Astronomy & Astrophysics, 624, A64
  * Petit et al. (1999) Petit, J.-M., Morbidelli, A., & Valsecchi, G. B. 1999, Icarus, 141, 367
  * Portegies Zwart et al. (2021) Portegies Zwart, S., Torres, S., Cai, M. X., & Brown, A. G. A. 2021, A&A, 652, A144
  * Prialnik (2021) Prialnik, D. 2021, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 53, AAS/Division for Planetary Sciences Meeting Abstracts, 307.10
  * Price et al. (2021) Price, E. M., Cleeves, L. I., Bodewits, D., & Öberg, K. I. 2021, ApJ, 913, 9
  * Ronnet et al. (2018) Ronnet, T., Mousis, O., Vernazza, P., Lunine, J. I., & Crida, A. 2018, AJ, 155, 224
  * Rubin et al. (2020) Rubin, M., Engrand, C., Snodgrass, C., et al. 2020, Space Sci.Rev., 216, 102
  * Steckloff et al. (2021) Steckloff, J. K., Lisse, C. M., Safrit, T. K., et al. 2021, Icarus, 356, 113998
  * Tsiganis et al. (2005) Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459