# Quantum-Enhanced Machine Learning for Covid-19 and Anderson Insulator
Predictions

Paul-Aymeric McRae  Michael Hilke  Department of Physics, McGill University,
Montréal, Canada H3A 2T8  Mila - Quebec AI Institute, Montréal, Québec, Canada
INTRIQ - Quebec Institute for Quantum Information

###### Abstract

Quantum Machine Learning (QML) algorithms to solve classifications problems
have been made available thanks to recent advancements in quantum computation.
While the number of qubits are still relatively small, they have been used for
“quantum enhancement” of machine learning. An important question is related to
the efficacy of such protocols. We evaluate this efficacy using common
baseline data sets, in addition to recent coronavirus spread data as well as
the quantum metal-insulator transition in three dimensions. For the
computation, we used the 16 qubit IBM quantum computer. We find that the
“quantum enhancement” is not generic and fails for more complex machine
learning tasks.

###### keywords:

Quantum Machine Learning , COVID-19 , Metal-Insulator , Quantum Computing

##  1 Introduction

In the last decade or so, quantum computing has seen huge leaps in innovation
and public awareness [1]. It has now become possible to run programs on real
quantum computers [2] or to realistically simulate them on a classical machine
[3]. Theoretically a quantum computer can offer great speed-up [4] or better
results [5] as compared to algorithms performed on a classical binary system.
Given these advances in the field it became interesting to combine it with
another rapidly-expanding topic in computer science: machine learning, [6, 7,
8] leading to the aptly-named domain of Quantum Machine Learning (QML). Thanks
to IBM’s QISKit (Quantum Information Software Kit) [2], QML has become
tangible and accessible [9], with interesting results.

It is useful to define and classify different forms of QML. An effective
representation is to consider the usual machine learning (ML) framework, where
we have an input, I, a computation, C and an output O. In conventional ML all
three (I, C and O) are performed classically. In QML, any of I, C, or O or any
combination thereof can be quantum. Arguably, the simplest approach is to take
a quantum system, for example an electron in a random potential (Anderson
localization) [10], which we label quI (quantum input) and use classical
machine learning techniques to improve the classification boundary of
localized and extended states [11]. Here we push the scheme one step further,
by using a partial quantum C (or quC). In this classification scheme, both the
input and the computation is partially quantum as we describe below.

In contrast, it is possible to use a purely classical I, for example the
Covid-19 infection rate in the US, but use a hybrid quantum C. In some cases,
this approach was shown to lead to an improvement in the classification
outcome, which was dubbed, “quantum enhancement” [12]. We use this scheme to
explore the possible “quantum enhancement” of different data samples, such as
the more complex Covid-19 data mentioned above. In all cases, we consider a
classical O, so that we only consider the cases I-quC-O and quI-quC-O as our
QML framework.

###  1.1 Support Vector Machines

The flavour of machine learning used here is a Support Vector Machine (SVM).
SVMs are used to classify data into different categories (usually two) [13,
14, 15]. They are a type of supervised learning [14, 15]; that is, the branch
of machine learning in which a already-classified data set is used as an input
on which the machine trains [15]. Upon successful training, the SVM can then
predict the category in which a new data point belongs [15, 16, 17].

Support Vector Machines work by taking the training data set and virtually
plotting the data in $d$-dimensional space where $d$ is the number of
parameters (or features) of each data point. It then attempts to find the
$(d-1)$-dimensional hyperplane which separates the data into two classes (in
the case of data with two parameters this is a line) [18]. Specifically, it
aims to find the hyperplane with the largest margin between points of either
class [14, 18]. Points lying on this margin are known as support vectors [18].

In many cases it is impossible to draw a separating hyperplane between classes
of data, especially when the boundary is highly nonlinear [18]. In such a case
where the SVM is unable to find a suitable separating hyperplane, a technique
known as a “kernel trick” is used, which involves using a non-linear function
(often called a Feature Map [19]) to project the data into a higher dimension
where a separating hyperplane may be found [20].

In this work, all classical SVM operations have been performed using the SVM
Python module developed by scikit-learn, whose documentation can be found in
[14].

Figure 1: Two-parameter ($x_{1}$, $x_{2}$) data set with two categories (clubs
and diamonds) along with possible separating hyperplanes. The green line does
not separate the data into two sets. The orange and purple lines do separate
the data but do not do so maximally. The blue line is the maximally separating
hyperplane with support vectors in blue circles.

###  1.2 Quantum Support Vector Machines

The quantum version of the SVM is best described as “quantum-assisted” or
“quantum-enhanced” [19] in the sense that the algorithm is largely classical
with certain operations performed by a quantum processor (real or simulated).
The Quantum SVM (or QSVM) made available by IBM QISKit functions the same way
as a classical SVM but the feature mapping [19] and some operations of the
hyperplane calculation [15] are handled by the quantum processor.

###  1.3 Variational Quantum Classifier

An important quantum alternative to the QSVM is the Variational Quantum
Classifier (VQC). This belongs to a class known as Variational Quantum
Algorithms, which are hybrid classical-quantum algorithms which iteratively
use the measurement of a parametrized quantum circuit as the input of a
classical process. The classical algorithm’s output is then passed as the
parameter of the quantum circuit, and the loop then repeats [21, 22]. The VQC
encodes a data input as a quantum state, which is then processed and the
quantum probability is then measured, which serves as the base inference for
the binary classification. By refining this inference, the VQC can make
predictions based on data.

(a)

(b)

Figure 2: Sample circuits from (a) QSVM and (b) VQC runs on a 5-qubit system.
Wrap-around for visualization.

###  1.4 IBM Quantum Experience and QISKit

The term “quantum computing” can refer to several different types of machines
[1]. Here we refer specifically to those publicly available via the IBM
Quantum Experience. These operate on the model of quantum circuits [23], a
system which can be defined as follows: a classical computer stores
information in bits which can hold a value of either 0 or 1, and can perform
operations on these bits using a combination of logic gates (specifically AND,
OR, NOT) [24]. By analogy a quantum circuit uses qubits which can store
superpositions of $\ket{0}$ and $\ket{1}$ states and operates on them using
quantum gates [25], which in quantum mechanical terms can be described by
operators on two-level quantum systems and product states of such systems
[26].

IBM has developed cloud-based access to real quantum circuit-based computers
as part of the IBM Quantum Experience [27] for which software can be written
using QISKit [1]. Several backends are available, with 5 [28] or 16 qubits
[29] or a simulator of over 30 qubits [30]. Their real computers are built
using transmons [31], a type of superconducting charge qubit [32, 33], with
algorithms designed to reduce errors that arise due to the fragility of
quantum information [34, 35].

###  1.5 Sources of Data

For our work, we use publicly available data sets provided by UC Irvine as
well as computed data. The first UCI data set is the “Breast Cancer Wisconsin
(Diagnostic) Data Set”, which contains 569 instances of tumor data each
containing 32 attributes which are used to determine whether the tumor is
benign or malignant. The second UCI set is known as “Wine Data Set” and
contains 178 instances of data sourced from wine, each containing 13
attributes which are used to classify which of three cultivars the wine was
sourced from [36]. These two data sets are often used as baselines for QML
trials, including in [37] as well as in IBM’s own QISKit tutorials [38].

Our generated data falls into three classes: ad-hoc data, Anderson data and
Coronavirus data. The ad-hoc data describes all data generated on a 20 by 20
grid, some of which use Python’s built-in pseudorandom functions. The
coronavirus data used was compiled by The New York Times and publicly released
on GitHub [39]. It contains the cumulative number of cases in the United
states on a county and state level. For this, the county-level data is used.
Finally the Anderson data is based on Anderson’s seminal result, where he
showed that certain materials can undergo a quantum phase transition between
conductor and an insulator based on the strength of disorder in the system
[10]. This can be numerically simulated [40] by varying energy and disorder to
get a two-dimensional data picture. The data is generated using an iterative
Green’s function technique to obtain the transmission coefficient as a
function of disorder strength [41].

##  2 Methods

###  2.1 Classification

As SVMs are designed to classify data into one of two sets (in most cases)
based on a training input, all data used needed some parameter that delineates
the set to which a data point would belong. The training input is defined by
taking the entirety of the classified data set and using 67% to fill the
training set. The remaining 33% is used as a “testing set” to verify the
accuracy of the SVM. The classifier used per data type is described as
follows.

For both the “Wine” and “Breast Cancer” data sets, QISKit provides functions
to procure the data such that it is already classified, and each point can be
reduced in dimensionality using the PCA algorithm. In our experiments, we
primarily reduced the feature dimension to two, though certain experiments use
three- or four-dimensional feature data.

In the case of the ad-hoc data the classifier is defined along with the data,
with the goal to create a particular “shape”. As mentioned, all data was
generated on a $20\times 20$ grid with increment size 1. The four different
generated data sets are presented in Figure 3.

(a)

(b)

(c)

(d)

Figure 3: Ad-Hoc data. (a) Perfectly linearly separable data; (b) Alternating
category between adjacent points; (c) Randomly assigned points; (d) “Circle”
of one category surrounded by the other

The coronavirus example uses data points consisting of a US county’s
population and its total number of cases, per day since the first recorded
case in the country. The classifier was whether a county is deemed “safe” or
“unsafe” based on whether the famous “curve” (i.e. the number of new cases per
day) was sufficiently flat or in fact decreasing. This was done by observing
the curve locally by taking the change in the number of cases over three days
preceding May 8th (an arbitrary choice). The derivative of cases per day over
the three days was then averaged. If this average was below a certain
threshold, i.e. the curve was flat or decreasing, the county was deemed safe.
Otherwise, it is classified as unsafe. The ML experiments were run using
thresholds of 1, 3, and 5 new cases per day on average. Only counties which
had a nonzero number of cases over the May $5^{th}$ to $8^{th}$ interval were
included in the data set. Likewise counties which changed the geographical
area of reported data over this time period were excluded.

Figure 4: Coronavirus data. Yellow counties are unsafe (curve not flattened),
teal counties are safe (curve flattened) and purple countries have no cases or
insufficient data. We define “Safe” as $\leq 5$ new cases per day on average.

The Anderson data was generated numerically by considering an $L\times L\times
L$ cube between two semi-infinite leads. The conductivity of the material is
computed by the transmission coefficient via the Landauer formula as a
function of Energy ($E$), disorder ($V_{a}$), cube side length ($L$) for
different disorder configurations. The transmission coefficient is obtained
from the Green’s function by tracing over all channels [41]. The full Green’s
function is computed using a layer by layer iterative computation of the
partial Green’s function [42]. The conductance (which we denote $G$) was
calculated for varying values of the parameters ($E$, $L$, $V_{a}$ and
disorder configurations). We used 100 different disorder configurations for
the averaging and a maximum of $L=14$. The dimensionality of this matrix was
reduced by log-averaging the value of $G$ for all disorder configurations for
fixed $E$, $V_{a}$, and $L$, leading to $\bar{G}$. The following finite-size
scaling method was used to differentiate metallic and insulating states [43].
A least-squares linear fit to approximate $\bar{G}$ as a function of $L$,
whose slope (or partial derivative $\frac{\partial\bar{G}}{\partial L}$) was
recorded for each $E$ and $V_{a}$. The classification was then made based on
this value: a positive $\frac{\partial\bar{G}}{\partial L}$ means the system
is conductive, while a negative number indicates the system is an insulator.
Values within $1\times 10^{-6}$ of $0$ were ignored.

(a)

(b)

Figure 5: Anderson data: Plot (a) shows the raw generated data where the color
scale corresponds to $\mbox{sgn}(\partial_{L}G)|\log|\partial_{L}G||$ and plot
(b) is the classified data: $\mbox{sgn}(\partial_{L}G)$ (red corresponds to a
conductor, blue to an insulator, and white was excluded with values
$|\partial_{L}G)|<10^{-6}$ ).

###  2.2 Randomization

To achieve an apparently random state, two principles were kept in mind:
reproducibility and invertibility. As such, “randomization” was achieved using
custom made Pseudo-Random Number Generators (PRNG), which generate numbers
that appear random but whose behaviour can actually be predicted [44].
Specifically given the same input and the same defining parameter (dubbed the
“seed”) a PRNG will always return the same value [44]. Here, each of the two
parameters for each data point was run through the same PNRG, though different
problems used different generators.

We can define a PRNG as follows:

$rand(x)=1+(a\cdot x+b)\mod m$ (1)

Where $b$ is an arbitrary integer, and $a$ and $m$ are prime numbers, with
$a<m$ and $m$ chosen such that it is greater than all the values of the data,
in order to guarantee invertibility. The inverse can be found as:

$rand^{-1}(y)=a^{-1}\cdot(y-1-b)\mod m$ (2)

Where $a^{-1}$ is the modular inverse of $a$ (that is, the value such that
$a\cdot a^{-1}\mod m\equiv 1$. This value can be found by iterating through
all integers between $1$ and $m$ until a number satisfying this condition is
found. If $a$ and $m$ are coprime (which they must be since they are both
prime) then this inverse must exist and is unique [45].

###  2.3 Data Reduction

During the experimentation, we developed a technique to reduce the size of the
data set. Notably this had the benefit of removing outliers, which we observed
increased the accuracy. As well, with a smaller data set the training time
would be reduced. To achieve this reduction it was important to keep data
point more relevant to the training (that is, points closer to the data
boundary where classification becomes more difficult). This was achieved by
using a classical linear SVM to draw the hyperplane that best separates the
data into two classes, then selecting data points that fall within some
distance of this hyperplane. The chosen distance was tuned manually for best
results.

###  2.4 Machine Learning Procedures

Classical machine learning experiments were done using the scikit-learn Python
package [14] and specifically the SVM module. Most experiments were done “out
of the box” by using a linear kernel, but occasionally the machine was tweaked
by tuning the hyperparameters, specifically by changing the kernel to a
polynomial or RBF. In some trials, the data was preprocessed using a standard
scaler, which rescales the data by subtracting the mean of the data from the
sample and dividing the result by the standard deviation. This drastically
decreases the training time, especially on large data sets. As a rule, on data
which was randomized, unscaled data would not be tested as the algorithm would
not converge in reasonable time.

All quantum machine learning was done using IBM QISKit’s QSVM module [46]. The
feature map used was Second-order Pauli-Z evolution circuit [47] using two
repeated circuits and linear entanglement. QML experiments were performed on
two backends: the QASM simulator [3] and the 16-qubit computer in Melbourne
[29].

##  3 Results and Analysis

###  3.1 Volatility of Results

Havenstein et. al. [37] show that a VQC can outperform a classical SVM on both
the Wine and Breast Cancer data sets, a result we were unable to reproduce.
This can be accounted by the volatility of running similar, but not identical
trials. For instance we notice that accuracy differences can arise from minor
changes in the way the data is partitioned into training/test data, including
the seed (called state in the code) used to pseudo-randomly shuffle the data
as well as the size of the sets. We report these difference in Table 1.

Table 1: Accuracy of different runs using slightly different train/test
partitions for the QASM simulator (first table) and the IBMQ 16-Qubit
Melbourne Backend (second table).

Data Partitioning | Accuracy (%)
---|---
Using 33-37 test/train split, random state 42 | 75.0
Using 30-70 test/train split, random state 42 | 72.22
Using 33-67 test/train split, random state 12 | 75.0
Using 30-70 test/train split, random state 12 | 66.67

Data Partitioning | Accuracy (%)
---|---
Using 33-37 test/train split, random state 42 | 75.0
Using 30-70 test/train split, random state 42 | 85.0
Using 33-67 test/train split, random state 12 | 66.67
Using 30-70 test/train split, random state 12 | 72.22

We also note that on successive trials run on either a real quantum backend or
the statevector simuluator, accuracy would fluctuate, which makes the outputs
of these machines stochastic and non-reproducible.

###  3.2 Time Complexity and Runtime

While experimenting with QML we noted very long training and prediction times
when using large data sets. Naturally this raised the question of how runtimes
scale with the size of the input. We observed that the duration of a QSVM run
on the QASM Simulator scaled as $O(n^{2})$ where $n$ is the size of the
training data (Fig. 6). We similarly observed linear $O(n)$ time using a VQC
with a QASM simulator (Fig. 6) as well as an exponential $O(2^{d})$ time
dependence on the feature dimension $d$ (see Figure 7). The simulator was used
rather than real backend since the time to run experiments on them depends on
the number of jobs in the queue which is a factor beyond our control. We
demonstrate the time dependence in the figures using least-squares fits.

As presented in an overview by [48], the theorized time complexity of the QSVM
algorithm is $\tilde{O}(log(nd)\kappa^{2}/\epsilon^{3})$ where $\kappa$ is a
quantity dependent on the kernel (specifically, it is the condition number of
a matrix which must be inverted), $\epsilon$ is the accuracy, and we use the
notation $\tilde{O}(h(n))=O(h(n)log(h(n))$ (i.e. suppressing slowly varying
multiplicative terms). This notably differs from our observation, which is
likely to arise from two factors: $(1)$ the fact that a simulator was used and
not a pure quantum circuit and $(2)$ the difference between implementation of
an algorithm and theoretical performance. Expanding on this second point,
consider the analogous example of binary search, which can be completed in
$O(log(n))$ time, but only if the list is sorted. If not, the list must be
sorted first, and this overhead will cost at least $O(nlog(n))$ time.

By exploring the implementation of the QSVM we can see the $O(n^{2})$ comes
from computation of the kernel matrix; that is, the matrix of the kernel
function computed between all pairs of data points. Since for $n$ datapoints
there are $n^{2}$ such pairs, the origin of the quadratic computational
complexity is obvious. The VQC’s linear dependence on $n$ is harder to
explain, though a heuristic explanation is the absence of any quadratic-time
functions and the fact that the algorithm likely makes at least one pass over
all the data. We were, however, unable to explain the exponential dependence
on feature dimension.

Figure 6: Computation time for quantum classification algorithms versus number
of training data points, along with least-squares fits. Figure 7: Computation
time for VQC versus feature dimension. We have both the exponential evolution
(which is captured by big-O, which drops slowly increasing terms) and the
closer fit: an exponential plus a linear function.

###  3.3 Baselines Revisited

A quick search through literature or online documentation will reveal that the
Wine and Breast Cancer data sets are commonly used to test QML algorithms [49,
37, 50]. However, these are often run without comparisons to a classical
equivalent, such as the case of the QISKit tutorials [38] or the thesis by
Mendozza [50]. Others, such as Havenstein et. al. [37] make explicit
comparisons, though we were not able to reproduce their exact results for
reasons stated in Subsection 3.1.

Here we present our findings for our experimentation, with varying parameters
such as feature dimension, number of data points, and backend used. In no
experiments using either the Breast Cancer or Wine data sets were we able to
observe a better accuracy of a quantum method as compared to classical
counterparts. This would seem to indicate that claims of a quantum advantage
using these baselines may be overstated, and that a new baseline should be
explored which more clearly presents a quantum advantage.

Figure 8: Accuracy of various classification methods versus number of training
samples from the Breast Cancer data set. The classical methods (Linear and RBF
kernel SVMs) consistently outperform the QML.

###  3.4 Results on Coronavirus and Anderson

While using a benchmark data set to reveal details on the performance of a
machine learning algorithm, the most important question when exploring a novel
technique is how well does it work in practice. For this we used two data sets
from “real” sources (Coronavirus and Anderson data) as well as several
manufactured data distributions (which we refer to as Ad-Hoc). In this
subsection we present our conclusions from this data. Our raw results can be
found in B.

Here again we note that in nearly every case, a form of classical
classification outperforms a quantum counterpart. In some rare instances,
quantum narrowly displays better accuracy than classical algorithms, generally
when the setup has a degree of randomness involved. On these data sets we
observe the larger prevalence of a potentially catastrophic behaviour: some of
the quantum methods achieve an artificially high accuracy by solely predicting
one class. This generally occurs when the data disparity is skewed towards one
class. However, unlike the classical methods implemented in scikit-learn,
QISKit classification algorithms do not allow the same degree of
hyperparameter tuning, which makes this behaviour potentially unavoidable.

##  4 Conclusion

While quantum computation is a field of rapidly evolving interest, the
practicality of quantum classification algorithms lags behind that of simple
classical methods. Lower accuracy and longer times, caused by the high
complexity of simulation and the relatively limited real backends, have caused
QML to be a novelty more than a great leap in innovation. Nevertheless, with
refinement of existing algorithms, including more variability in the tuning of
hyperparameters, QML may in the future emerge as a strong contender in the
field of machine learning.

## References

  * National Academies of Sciences, Engineering, and Medicine [2019] National Academies of Sciences, Engineering, and Medicine .  Quantum Computing: Progress and Prospects.  Washington, DC: The National Academies Press; 2019.  ISBN 978-0-309-47969-1.  URL: https://www.nap.edu/catalog/25196/quantum-computing-progress-and-prospects. doi:10.17226/25196.
  * IBM Research Editorial Staff [2019] IBM Research Editorial Staff .  QISKit – Quantum Information Software Kit for Quantum Computation.  IBM Research Blog; 2019.  URL: https://www.ibm.com/blogs/research/2018/02/qiskit-index/.
  * Qiskit Development Team [2020a] Qiskit Development Team .  Qasmsimulator.  Qiskit 0.19.6 documentation; 2020a.  URL: https://qiskit.org/documentation/stubs/qiskit.providers.aer.QasmSimulator.html.
  * Ronnow et al. [2014] Ronnow TF, Wang Z, Job J, Boixo S, Isakov SV, Wecker D, et al.  Defining and detecting quantum speedup.  Science 2014;345(6195):420–4.  URL: https://doi.org/10.1126/science.1252319. doi:10.1126/science.1252319.
  * Montanaro [2016] Montanaro A.  Quantum algorithms: an overview.  npj Quantum Information 2016;2(1).  URL: https://doi.org/10.1038/npjqi.2015.23. doi:10.1038/npjqi.2015.23.
  * Dunjko and Briegel [2018] Dunjko V, Briegel HJ.  Machine learning & artificial intelligence in the quantum domain: a review of recent progress.  Reports on Progress in Physics 2018;81(7):074001.
  * Arunachalam and de Wolf [2017] Arunachalam S, de Wolf R.  Guest column.  ACM SIGACT News 2017;48(2):41–67.  URL: https://doi.org/10.1145/3106700.3106710. doi:10.1145/3106700.3106710.
  * Ciliberto et al. [2018] Ciliberto C, Herbster M, Ialongo AD, Pontil M, Rocchetto A, Severini S, et al.  Quantum machine learning: a classical perspective.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 2018;474(2209):20170551.  URL: https://doi.org/10.1098/rspa.2017.0551. doi:10.1098/rspa.2017.0551.
  * Havlicek et al. [2017] Havlicek V, Temme K, Córcoles A, Liu P, Chen R, Pistoia M, et al.  Qiskit aqua: Experiment with classification problem with quantum-enhanced support vector machines.  2017\.
  * Anderson [1958] Anderson PW.  Absence of diffusion in certain random lattices.  Phys Rev 1958;109:1492–505.  URL: https://link.aps.org/doi/10.1103/PhysRev.109.1492. doi:10.1103/PhysRev.109.1492.
  * Ohtsuki and Ohtsuki [2017] Ohtsuki T, Ohtsuki T.  Deep learning the quantum phase transitions in random electron systems: Applications to three dimensions.  Journal of the Physical Society of Japan 2017;86(4):044708.
  * Dunjko et al. [2016] Dunjko V, Taylor JM, Briegel HJ.  Quantum-enhanced machine learning.  Physical review letters 2016;117(13):130501.
  * Cortes and Vapnik [1995] Cortes C, Vapnik V.  Support-vector networks.  Machine Learning 1995;20(3):273–97.  URL: https://doi.org/10.1007/bf00994018. doi:10.1007/bf00994018.
  * Pedregosa et al. [2011] Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 2011;12:2825–30.
  * Rebentrost et al. [2014] Rebentrost P, Mohseni M, Lloyd S.  Quantum support vector machine for big data classification.  Physical Review Letters 2014;113(13).  URL: https://doi.org/10.1103/physrevlett.113.130503. doi:10.1103/physrevlett.113.130503.
  * Mohri [2012] Mohri M.  Foundations of Machine Learning (Adaptive Computation and Machine Learning).  The MIT Press; 2012.  ISBN 026201825X.  URL: https://www.xarg.org/ref/a/026201825X/.
  * Russell [2009] Russell S.  Artificial Intelligence: A Modern Approach (3rd Edition).  Pearson; 2009.  ISBN 0136042597.  URL: https://www.xarg.org/ref/a/0136042597/.
  * Berwick [2003] Berwick R.  An idiot’s guide to support vector machines (SVMs).  2003\.
  * Havlíček et al. [2019] Havlíček V, Córcoles AD, Temme K, Harrow AW, Kandala A, Chow JM, et al.  Supervised learning with quantum-enhanced feature spaces.  Nature 2019;567(7747):209–12.  URL: https://doi.org/10.1038/s41586-019-0980-2. doi:10.1038/s41586-019-0980-2.
  * Boser et al. [1992] Boser BE, Guyon IM, Vapnik VN.  A training algorithm for optimal margin classifiers.  In: Proceedings of the fifth annual workshop on Computational learning theory - COLT 92. ACM Press; 1992,URL: https://doi.org/10.1145/130385.130401. doi:10.1145/130385.130401.
  * LaRose et al. [2019] LaRose R, Tikku A, O’Neel-Judy É, Cincio L, Coles PJ.  Variational quantum state diagonalization.  npj Quantum Information 2019;5(1).  URL: https://doi.org/10.1038/s41534-019-0167-6. doi:10.1038/s41534-019-0167-6.
  * Saldanha et al. [2020] Saldanha D, Agastya A, G P, and MRP.  A review of supervised variational quantum classifiers.  International Journal of Engineering Research and 2020;V9(04).  URL: https://doi.org/10.17577/ijertv9is040485. doi:10.17577/ijertv9is040485.
  * Qiskit Development Team [2020b] Qiskit Development Team .  Quantum circuits.  Qiskit 0.19.6 documentation; 2020b.  URL: https://qiskit.org/documentation/apidoc/circuit.html.
  * Patterson and Hennessy [2013] Patterson DA, Hennessy JL.  Computer Organization and Design: The Hardware/Software Interface.  Fifth ed.; Morgan Kaufmann; 2013\.  ISBN 0124077269.
  * Asfaw et al. [2020] Asfaw A, Bello L, Ben-Haim Y, Bravyi S, Capelluto L, Vazquez AC, et al.  Learn Quantum Computation Using Qiskit.  2020\.  URL: http://community.qiskit.org/textbook.
  * Nielsen and Chuang [2011] Nielsen MA, Chuang IL.  Quantum computation and quantum information.  10th anniversary ed.; Cambridge University Press; 2011.  ISBN 9781107002173.
  * ibm [2016] IBM Makes Quantum Computing Available on IBM Cloud to Accelerate Innovation.  2016\.  URL: https://www-03.ibm.com/press/us/en/pressrelease/49661.wss.
  * IBM Q Team [2019a] IBM Q Team .  IBM Q 5 Yorktown backend specification V1.3.0.  2019a.  URL: https://github.com/Qiskit/ibmq-device-information/tree/master/backends/yorktown/V1.
  * IBM Q Team [2019b] IBM Q Team .  IBM Q 16 Melbourne backend specification V1.3.0.  2019b.  URL: https://github.com/Qiskit/ibmq-device-information/blob/master/backends/melbourne/V1/.
  * IBM Research Editorial Staff [2018] IBM Research Editorial Staff .  An Open High-Performance Simulator for Quantum Circuits.  IBM Research Blog; 2018.  URL: https://www.ibm.com/blogs/research/2018/05/quantum-circuits/.
  * Gambetta et al. [2017] Gambetta JM, Chow JM, Steffen M.  Building logical qubits in a superconducting quantum computing system.  npj Quantum Information 2017;3(1):2.  URL: https://doi.org/10.1038/s41534-016-0004-0. doi:10.1038/s41534-016-0004-0.
  * Kjaergaard et al. [2020] Kjaergaard M, Schwartz ME, Braumüller J, Krantz P, Wang JIJ, Gustavsson S, et al.  Superconducting qubits: Current state of play.  Annual Review of Condensed Matter Physics 2020;11(1):369–395.  URL: http://dx.doi.org/10.1146/annurev-conmatphys-031119-050605. doi:10.1146/annurev-conmatphys-031119-050605.
  * Koch et al. [2007] Koch J, Yu TM, Gambetta J, Houck AA, Schuster DI, Majer J, et al.  Charge-insensitive qubit design derived from the cooper pair box.  Physical Review A 2007;76(4).  URL: http://dx.doi.org/10.1103/PhysRevA.76.042319. doi:10.1103/physreva.76.042319.
  * Jerry Chow [2014] Jerry Chow .  Dealing with errors in quantum computing.  IBM Research Blog; 2014.  URL: https://www.ibm.com/blogs/research/2014/06/dealing-with-errors-in-quantum-computing/.
  * Chow et al. [2014] Chow JM, Gambetta JM, Magesan E, Abraham DW, Cross AW, Johnson BR, et al.  Implementing a strand of a scalable fault-tolerant quantum computing fabric.  Nature Communications 2014;5(1):4015.  URL: https://doi.org/10.1038/ncomms5015. doi:10.1038/ncomms5015.
  * Dua and Graff [2017] Dua D, Graff C.  UCI machine learning repository.  2017\.  URL: http://archive.ics.uci.edu/ml.
  * Havenstein et al. [2018] Havenstein C, Thomas D, Chandrasekaran S.  Comparisons of performance between quantum and classical machine learning.  SMU Data Science Review 2018;1(4).  URL: https://scholar.smu.edu/cgi/viewcontent.cgi?article=1047&context=datasciencereview.
  * Abraham et al. [2019] Abraham H, AduOffei , Agarwal R, Akhalwaya IY, Aleksandrowicz G, Alexander T, et al.  Qiskit: An open-source framework for quantum computing.  2019\.  doi:10.5281/zenodo.2562110.
  * Smith et al. [2020] Smith M, Yourish K, Almukhtar S, Collins K, Ivory D, Harmon A.  Coronavirus (Covid-19) Data in the United States.  GitHub; 2020.
  * Croy et al. [2006] Croy A, Römer RA, Schreiber M.  Localization of electronic states in amorphous materials: Recursive green’s function method and the metal-insulator transition at e $\neq$ 0\.  In: Hoffmann KH, Meyer A, editors. Parallel Algorithms and Cluster Computing. Berlin, Heidelberg: Springer Berlin Heidelberg.  ISBN 978-3-540-33541-2; 2006, p. 203–26.
  * Datta [2000] Datta S.  Nanoscale device modeling: the green’s function method.  Superlattices and microstructures 2000;28(4):253–78.
  * MacKinnon [1985] MacKinnon A.  The calculation of transport properties and density of states of disordered solids.  Zeitschrift für Physik B Condensed Matter 1985;59(4):385–90.
  * MacKinnon and Kramer [1981] MacKinnon A, Kramer B.  One-parameter scaling of localization length and conductance in disordered systems.  Physical Review Letters 1981;47(21):1546.
  * Vadhan [2012] Vadhan SP.  Pseudorandomness.  Foundations and Trends in Theoretical Computer Science 2012;7(1–3):1–336.  URL: http://dx.doi.org/10.1561/0400000010. doi:10.1561/0400000010.
  * Shoup [2005] Shoup V.  A computational introduction to number theory and algebra.  Cambridge: Cambridge University Press; 2005.  ISBN 978-0521851541.
  * Qiskit Development Team [2020c] Qiskit Development Team .  Qsvm.  Qiskit 0.19.6 documentation; 2020c.  URL: https://qiskit.org/documentation/stubs/qiskit.aqua.algorithms.QSVM.html.
  * Qiskit Development Team [2020d] Qiskit Development Team .  Zzfeaturemap.  Qiskit 0.19.6 documentation; 2020d.  URL: https://qiskit.org/documentation/stubs/qiskit.circuit.library.ZZFeatureMap.html.
  * Duan et al. [2020] Duan B, Yuan J, Yu CH, Huang J, Hsieh CY.  A survey on hhl algorithm: From theory to application in quantum machine learning.  Physics Letters A 2020;384(24):126595.  URL: http://www.sciencedirect.com/science/article/pii/S037596012030462X. doi:https://doi.org/10.1016/j.physleta.2020.126595.
  * Sharma [2020] Sharma S.  Qeml (quantum enhanced machine learning): Using quantum computing to enhance ml classifiers and feature spaces.  2020\.  arXiv:2002.10453.
  * Mendozza et al. [2019] Mendozza R, Acampora G, Vitiello A.  A quantum implementation of support vector machine algorithm on ibm qx processors.  Ph.D. thesis; University of Naples Federico II Department of Physics; 2019\.  URL: http://www.fisica.unina.it/documents/12375590/13725484/2862_MendozzaR_16-10-2019.pdf/2fe12929-a8e7-46b7-bae8-b3f06d4c1f8f.

##  Appendix A RBF kernel classical SVM

An RBF kernel classical SVM without data preprocessing can be used to locate
the data boundary between the insulator/conductor phase transition in the
Anderson problem to 96.67% accuracy. Even when the data was reduced and the
boundary became less visibly obvious, this method of classification was able
to obtain 79.7% accuracy. We can conclude that this is a good metric for
classifying the properties of a a disordered solid.

We also observed that classically, using an RBF kernel, we were able to
predict whether a county in the USA was deemed safe to 90.99% accuracy. A
linear kernel yielded a slightly smaller accuracy at 90.56%. As a verification
of the quality of this process, we used all data collected from the 8th of May
as a training set and tested the machine on several other dates, with both RBF
and linear kernels, as reported in Table 2. We observe that although RBF is
incrementally better on present-day data, the predictive capacity of the
linear kernel surpasses it. We also note the linear kernel’s decaying
prediction accuracy, which we can likely attribute to the data drifting
further from its initial state.

Table 2: Prediction accuracy of a machine trained using May 8th COVID data. | May 9 | May 18 | May 28 | June 8
---|---|---|---|---
Linear | 91.55% | 91.07% | 90.20% | 88.34%
RBF | 81.01% | 82.61% | 81.84% | 82.72%

##  Appendix B Data Tables

Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy |  Only one Class Predicted?
---|---|---|---|---
SVM, Linear Kernel | Local Processor |  0 |  85.0 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  85.0 |  False
SVM, RBF Kernel | Local Processor |  0 |  85.0 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  85.0 |  False
QSVM | qasm_simulator |  11 |  75.0 |  False
VQC | qasm_simulator |  48 |  70.0 |  False
QSVM | ibmq_16_melbourne |  1966 |  75.0 |  False
QSVM | statevector_simulator |  0 |  75.0 |  False
VQC | statevector_simulator |  31 |  75.0 |  False
Table 3: Results of trials on breast cancer data set with 60 total samples. Trials run on ibmq_16_melbourne and statevector_simulator are the results from one trial and were not necessarily reproducible. Feature Dimension | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
2 | SVM, Linear Kernel, scaled | Local Processor |  0 |  100.0 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  100.0 |  False
QSVM | qasm_simulator |  6 |  100.0 |  False
VQC | qasm_simulator |  32 |  71.43 |  False
3 | SVM, Linear Kernel, scaled | Local Processor |  0 |  100.0 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  92.86 |  False
QSVM | qasm_simulator |  14 |  64.29 |  False
VQC | qasm_simulator |  58 |  42.86 |  False
4 | SVM, Linear Kernel | Local Processor |  0 |  92.86 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  92.86 |  False
SVM, RBF Kernel | Local Processor |  0 |  100.0 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  100.0 |  False
VQC | qasm_simulator |  128 |  50.0 |  False
Table 4: Results of trials of the Wine data set, with increasing feature dimension. Feature dimensions higher than 3 are not supported for QSVMs. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel | Local Processor |  0 |  100.0 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  100.0 |  False
SVM, RBF Kernel | Local Processor |  0 |  100.0 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  100.0 |  False
QSVM | qasm_simulator |  445 |  43.18 |  False
VQC | qasm_simulator |  315 |  54.55 |  False
All Data, Randomized | SVM, Linear Kernel | Local Processor |  16 |  50.76 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  59.85 |  False
SVM, RBF Kernel | Local Processor |  0 |  50.0 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  46.21 |  False
QSVM | qasm_simulator |  444 |  50.0 |  False
VQC | qasm_simulator |  314 |  54.55 |  False
Table 5: Results for trials using perfectly linearly separable Ad-Hoc data. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel | Local Processor |  0 |  50.0 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  50.0 |  False
SVM, RBF Kernel | Local Processor |  0 |  20.45 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  53.03 |  False
QSVM | qasm_simulator |  433 |  48.48 |  False
VQC | qasm_simulator |  301 |  52.27 |  False
All Data, Randomized | SVM, Linear Kernel | Local Processor |  17 |  47.73 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  45.45 |  False
SVM, RBF Kernel | Local Processor |  0 |  50.0 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  42.42 |  False
QSVM | qasm_simulator |  423 |  48.48 |  False
VQC | qasm_simulator |  300 |  62.12 |  False
Table 6: Results for trials using alternating grid Ad-Hoc data. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel | Local Processor |  0 |  52.63 |  True
SVM, Linear Kernel, scaled | Local Processor |  0 |  52.63 |  True
SVM, RBF Kernel | Local Processor |  0 |  52.63 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  46.62 |  False
QSVM | qasm_simulator |  424 |  49.62 |  False
VQC | qasm_simulator |  299 |  51.13 |  False
All Data, Randomized | SVM, Linear Kernel | Local Processor |  33 |  54.14 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  52.63 |  True
SVM, RBF Kernel | Local Processor |  0 |  52.63 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  50.38 |  False
QSVM | qasm_simulator |  422 |  50.38 |  False
VQC | qasm_simulator |  299 |  49.62 |  False
Table 7: Results for trials using random grid Ad-Hoc data. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel | Local Processor |  0 |  72.93 |  True
SVM, Linear Kernel, scaled | Local Processor |  0 |  72.93 |  True
SVM, RBF Kernel | Local Processor |  0 |  97.74 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  90.98 |  False
QSVM | qasm_simulator |  423 |  72.93 |  True
VQC | qasm_simulator |  299 |  57.89 |  False
All Data, Randomized | SVM, Linear Kernel | Local Processor |  42 |  66.17 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  72.93 |  True
SVM, RBF Kernel | Local Processor |  0 |  72.93 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  72.93 |  True
QSVM | qasm_simulator |  422 |  72.93 |  True
VQC | qasm_simulator |  300 |  54.89 |  False
Table 8: Results for trials using “circle” Ad-Hoc data. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel | Local Processor |  0 |  95.07 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  95.51 |  False
SVM, RBF Kernel | Local Processor |  0 |  96.67 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  96.0 |  False
QSVM | qasm_simulator |  14009 |  56.52 |  False
VQC | qasm_simulator |  1722 |  49.42 |  False
All Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  49.28 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  93.0 |  False
QSVM | qasm_simulator |  14042 |  54.64 |  False
VQC | qasm_simulator |  1707 |  54.2 |  False
Reduced Data | SVM, Linear Kernel | Local Processor |  0 |  76.69 |  False
SVM, Linear Kernel, scaled | Local Processor |  0 |  78.2 |  False
SVM, RBF Kernel | Local Processor |  0 |  79.7 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  77.0 |  False
QSVM | qasm_simulator |  493 |  54.14 |  False
VQC | qasm_simulator |  322 |  60.9 |  False
Reduced Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  50.38 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  60.0 |  False
QSVM | qasm_simulator |  485 |  60.15 |  False
VQC | qasm_simulator |  323 |  59.4 |  False
QSVM | ibmq_16_melbourne |  117525 |  57.14 |  False
QSVM | ibmq_16_melbourne |  (X) |  69.17 |  False
Table 9: Results of the trials using the Anderson data set. For the reduced data, we picked points within $\delta=2$ of the hyperplane learned using a linear SVM. Included are two separate runs on the reduced, randomized data with a real 16-qubit computer. The second run did not have time measured. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel, scaled | Local Processor |  0 |  90.56 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  90.99 |  False
QSVM | qasm_simulator |  28505 |  81.02 |  True
VQC | qasm_simulator |  2825 |  61.29 |  False
QSVM | ibmq_16_melbourne |  (X) |  81.02 |  True
All Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  81.02 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  90.35 |  False
QSVM | qasm_simulator |  28562 |  81.02 |  True
VQC | qasm_simulator |  2845 |  65.64 |  False
Reduced Data | SVM, Linear Kernel, scaled | Local Processor |  0 |  92.14 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  91.92 |  False
QSVM | qasm_simulator |  27485 |  84.61 |  True
VQC | qasm_simulator |  2709 |  61.13 |  False
Reduced Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  84.61 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  90.92 |  False
QSVM | qasm_simulator |  28533 |  54.04 |  False
VQC | qasm_simulator |  2826 |  64.78 |  False
Table 10: Results of Coronavirus trials where the cutoff for safety is 5 new cases per day. For the reduced data, we picked points within $\delta=2000$ of the hyperplane learned using a linear SVM. Time was not measured for the 16-qubit run. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel, scaled | Local Processor |  0 |  87.92 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  87.92 |  False
QSVM | qasm_simulator |  27450 |  74.89 |  True
VQC | qasm_simulator |  3003 |  58.9 |  False
All Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  74.89 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  86.33 |  False
QSVM | qasm_simulator |  28466 |  74.89 |  True
VQC | qasm_simulator |  3000 |  61.12 |  False
Reduced Data | SVM, Linear Kernel, scaled | Local Processor |  0 |  89.54 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  90.66 |  False
QSVM | qasm_simulator |  26579 |  78.53 |  True
VQC | qasm_simulator |  2869 |  59.96 |  False
Reduced Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  78.53 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  86.21 |  False
QSVM | qasm_simulator |  27890 |  78.53 |  True
VQC | qasm_simulator |  2849 |  64.07 |  False
Table 11: Results of Coronavirus trials where the cutoff for safety is 3 new cases per day. For the reduced data, we picked points within $\delta=2000$ of the hyperplane learned using a linear SVM. Data | Algorithm | Backend |  Time (Seconds, rounded) |  Accuracy (%) |  Only one Class Predicted?
---|---|---|---|---|---
All Data | SVM, Linear Kernel, scaled | Local Processor |  0 |  80.51 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  79.77 |  False
QSVM | qasm_simulator |  29125 |  60.06 |  True
VQC | qasm_simulator |  3005 |  52.86 |  False
All Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  60.06 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  78.5 |  False
QSVM | qasm_simulator |  30751 |  60.06 |  True
VQC | qasm_simulator |  3014 |  54.77 |  False
Reduced Data | SVM, Linear Kernel, scaled | Local Processor |  0 |  83.2 |  False
SVM, RBF Kernel, scaled | Local Processor |  0 |  84.66 |  False
QSVM | qasm_simulator |  26656 |  63.49 |  True
VQC | qasm_simulator |  2855 |  52.86 |  False
Reduced Data, Randomized | SVM, Linear Kernel, scaled | Local Processor |  0 |  63.49 |  True
SVM, RBF Kernel, scaled | Local Processor |  0 |  76.48 |  False
QSVM | qasm_simulator |  29351 |  63.49 |  True
VQC | qasm_simulator |  2775 |  57.11 |  False
Table 12: Results of Coronavirus trials where the cutoff for safety is 1 new
case per day. For the reduced data, we picked points within $\delta=2000$ of
the hyperplane learned using a linear SVM.


11institutetext: School of Aerospace Engineering, Tsinghua Universuty,
Beijing, 100084, China
11email<EMAIL_ADDRESS>
# Centaurs potentially in retrograde co-orbit resonance with Saturn

Miao Li 11 Yukun Huang 11 Shengping Gong<EMAIL_ADDRESS>
(Received XXX; accepted YYY)

###### Abstract

Aims. 2015 BZ509 is the first asteroid confirmed to be in retrograde co-orbit
resonance (or 1/-1 resonance) with the giant planets in the solar system.
While Saturn is the only giant planet whose Trojans are not discovered until
now, we identify some small bodies among Centaurs and Damocloids that are
potentially in 1/-1 resonance with Saturn in the present study.

Methods. We integrate numerically the motion of the 1000 clones (include the
nominal orbit) of each Centaur whose orbit has a semi-major axis between 9.3
au and 9.8 au and an inclination $i>90\degr$. To confirm and evaluate the 1/-1
resonant configurations mentioned above, we introduce an useful one degree
integrable approximation for planar 1/-1 resonance.

Results. We identify four candidates potentially in 1/-1 resonance with
Saturn. The capture in this particular resonant state during the 40000 yr
integration timespan is very common for 2006 RJ2 (906/1000 clones), 2006 BZ8
(878/1000 clones), and 2017 SV13 (998/1000 clones), and it is less likely for
2012 YE8 (426/1000 clones). According to our statistical results, 2006 RJ2 is
the best candidate to be currently in a 1/-1 mean motion resonance with
Saturn, and 2017 SV13 is another important potential candidate. Moreover, 2012
YE8 and 2006 BZ8 are also Centaurs of interest but their current and long-term
1/-1 resonant state with Saturn is less likely. The proportions of the clones
captured in the relative long-term stable co-orbit resonance (over 10000 yr)
are also given. The motions of the 2006 RJ2, 2015 BZ509, and 2006 BZ8 in the
solar system are just around the ideal equilibrium points of the 1/-1
resonance given by the planar semi-analytical model.

Conclusions. Small bodies in retrograde co-orbit resonance with giant planets
may be more common than previously expected. Identification of these potential
mysterious minor bodies encourages the search for such objects on a larger
scale in our solar system. The findings of this paper are also useful for
understanding the origin and dynamical evolution of the Centaurs and
Damocloids on retrograde orbits.

###### Key Words.:

celestial mechanics – minor planets, asteroids – planets and satellites:
dynamical evolution and stability – comets

##  1 Introduction

All of the major planets and the significant portions of the asteroids have
orbits around Sun with inclination $i<90\degr$ when viewed from above the
north ecliptic pole of our solar system (Wiegert et al., 2017). But small
bodies of the solar system on retrograde orbits do exist, such as irregular
satellites which have eccentric, highly inclined or even retrograde orbits
around their giant planets (Jewitt & Sheppard, 2003). Although the origins of
these mysterious retrograde objects are not clear until now, they are
generally believed to have been captured from heliocentric orbit during the
final phase of planetary accretion (Pollack et al., 1979; Heppenheimer &
Porco, 1977). Halley-type comets (HTCs) provide us some good candidates to
study the dynamics of retrograde objects, while their origin has long been a
matter of debate until now (Emel’yanenko & Bailey, 1998; Fernandez & Gallardo,
1999; Wiegert & Tremaine, 1999; Nurmi et al., 2002; Fouchard et al., 2014;
Olsson-Steel, 1988; Levison et al., 2001, 2006; Fernandez et al., 2016). Minor
bodies in retrograde orbits are rare in our solar system. Among the 753782
small bodies found so far, only 93 are in retrograde motion111Minor Planet
Center; https://www.minorplanetcenter.net/iau/MPCORB.html, retrieved 18
January 2018.. In recent years, the research about the heliocentric retrograde
orbital motion gradually multiplies but lots of details remain uncertain until
now. Morais & Namouni (2013a) identified some asteroids trapped into
retrograde resonance for thousands of years with Jupiter and Saturn.
Greenstreet et al. (2012) suggested that the retrograde orbits inside the
Jupiter’s heliocentric path may be produced through the gravitational
interaction of Jupiter and Saturn. Surprisingly, after the discovery of the
first retrograde extrasolar planet, which is known as HAT-P-7b (Narita et al.,
2009; Winn et al., 2009), Triaud et al. (2010) identified some other
exoplanets on retrograde orbits. Interestingly, Chen et al. (2016) identified
a new retrograde trans-Neptunian object (TNO), nicknamed Niku, which is
confirmed in 7/9 mean motion resonance with Neptune by Morais & Namouni
(2017b). These findings further stimulate the researchers’ interest in
retrograde motion. Understanding the dynamics of the retrograde minor bodies
is significant for the study of the mechanisms and origins of the major and
minor planets on retrograde orbits.

More than 7000 Trojans have been discovered until now, which are locked in
stable prograde co-orbit motion with Jupiter. Some Trojans of other planets
including Earth, Mars, Uranus and Neptune222Until now, 1 Earth Trojan, 9 Mars
Trojans, 1 Uranus Trojan and 17 Neptune Trojans have been discovered. Data
taken from website of JPL Small-Body Database Search Engine, retrieved 5 May
2018. have also been observed in recent years such as the first Trojan of
Earth, 2010 TK7 (Connors et al., 2011), but no Saturn Trojans have never been
identified (Hou et al., 2014). Interestingly, a retrograde co-orbit motion was
not found until recently. In recent years, some theoretical work on these
mysterious orbits has verified that the co-orbit retrograde motion can be
stable. Morais & Namouni (2013b) developed a semi-analytical model for
retrograde co-orbit motion through the averaged unexpanded disturbing
function. They have confirmed later that there is a great potential for the
retrograde orbits trapped in co-orbit resonance (Namouni & Morais, 2015). The
stability and properties of the spacial coorbital motion have been studied by
Morais & Namouni (2016) from numerical simulations. The authors confirmed that
a stable co-orbit motion can exist regardless of the orbital inclination,
which, of course, includes the retrograde case.

The recent discovery of the first asteroid trapped in retrograde co-orbit
resonance with Jupiter, 2015 BZ509, has renewed interest in the searching for
such mysterious objects (Wiegert et al., 2017). 2015 BZ509 is the first
asteroid in retrograde co-orbit resonance with a planet of our solar system,
and it has long-term stability. Wiegert et al. (2017) suggested that this
asteroid can keep this stable resonant state about a million years. The
retrograde co-orbit resonance with Jupiter protects it from disruptive close
encounters with the other planets. Namouni & Morais (2017) have found that
2015 BZ509 lies almost exactly at the capture efficiency peak by numerical
simulations. Indeed, Morais & Namouni (2013a) have identified that 2006 BZ8
could enter the co-orbit retrograde resonance with Saturn in the future. From
numerical simulations, Namouni & Morais (2015) showed that asteroids on
retrograde motion are more easily captured in mean motion resonances than
those on prograde motion. In this paper, we identify the dynamics of the
Centaurs which have co-orbit retrograde orbits around Saturn by numerical
simulations, trying to find some new small bodies in co-orbit retrograde
resonance with the giant planets of our solar system.

The article is structured as follows: in section 2, we introduce the main
properties of the retrograde mean motion resonance, and the co-orbit
retrograde resonance is discussed in detail. Section 3 gives a full
description of the numerical model used in our study. In section 4, we show
the general statistical results of our numerical simulations and analyze the
complicated resonant states we observed. In section 5, we analyze the phase-
space portrait of the co-orbit retrograde resonance and present a comparison
between the results obtained from the semi-analytical model and the numerical
simulations. Our conclusions and some discussion are given in the last
section.

##  2 Retrograde co-orbit resonance

A retrograde mean motion resonance between a small body and a planet (denoted
by superscript symbol ′) occurs when $pn-qn\arcmin\approx 0$, where $p$ and
$q$ are both positive integers, and $n$ and $n\arcmin$ are the mean motion
frequencies of the small body and the planet, respectively. We call such
resonance a retrograde $p/q$ resonance or $p/-q$ resonance.

To determine whether these minor bodies are in retrograde co-orbit resonances
with Saturn, we have to first get the correct expression of the particular
resonant angles defined by the Fourier expansion of the disturbing function
(Murray & Dermott, 1999). For retrograde mean motion resonance, Morais &
Namouni (2013b) presented a transformation of the relevant resonant term for
the corresponding retrograde resonance to keep with the conventional
definition of osculating orbital elements in prograde motion. The expression
of the longitude of the pericenter becomes $\varpi=\omega-\Omega$, while the
mean longitude of the asteroid becomes $\lambda=M+\omega-\Omega$, where
$\omega$, $\Omega$ and $M$ are the argument of pericenter, the longitude of
the ascending node and the mean anomaly, respectively.

Morais & Namouni (2013b, a, 2016) reported that the resonant angles for the
$p/-q$ resonance have the following forms:

$\phi=q\lambda-p\lambda\arcmin-(p+q-2k)\varpi+2k\Omega,$ (1)

where $p$, $q$, and $k$ are positive integers. $p+q\geq 2k$ must be satisfied,
and the lowest order terms are proportional to $\cos^{2k}(i/2)e^{p+q-2k}$,
where $i$ and $e$ are respectively the inclination and the eccentricity of the
small body. The resonant term corresponding to $k=0$ will be the dominant one
for the nearly coplanar retrograde motion because of $\cos(i/2)\approx 0$. In
particular, for retrograde co-orbit resonance, the lowest order resonant terms
are proportional to

$\begin{split}&e^{2}\cos(\phi),k=0,\\\ &\cos^{2}(i/2)\cos(\phi_{z}),k=1,\\\
\end{split}$ (2)

where $\phi=\lambda-\lambda\arcmin-2\varpi$ and
$\phi_{z}=\lambda-\lambda\arcmin+2\Omega$ are two resonant angles of the 1/-1
resonance. When we study the Centaurs in 1/-1 resonance with Saturn, we also
consider the possibility of the Centaurs in 2/-5 resonance with Jupiter as the
positions of Jupiter and Saturn are near the 5/2 resonance.

##  3 Descriptions about the numerical model

To find all of the small bodies in the 1/-1 resonance with Saturn, we use the
JPL Small-Body Database Search Engine333JPL Small-Body Database Search Engine;
https://ssd.jpl.nasa.gov/sbdb_query.cgi to search for small bodies with
retrograde orbits (inclination $i>90\degr$) and semi-major axes ranging from
9.3 au to 9.8 au (semi-major axis of Saturn $a_{\rm Saturn}$ = 9.54 au). There
are four Centaurs meet our limitations, which are 2006 RJ2, 2006 BZ8, 2012 YE8
and 2017 SV13. We studied Centaurs on retrograde orbits, which are usually
believed strongly chaotic because of their giant-planet crossing orbit. The
discovered Centaurs mainly have two kinds of dynamical evolution (Bailey &
Malhotra, 2009; Volk & Malhotra, 2013): random walk and resonance sticking. As
these selected small bodies among Damocloids (Jupiter Tisserand’s parameter
$T<2$), Jewitt (2005) suggested that the most obvious origin would be Oort
Cloud or Halley-type comet.

We numerically integrate the motions of these Centaurs in the solar system
using the MERCURY (Chambers, 1999) package with an accuracy parameter
$10^{-12}$. The Burlisch-Stoer integrator is employed while considering all
the gravity of the eight planets throughout the integration. As listed in
Table 1 and Table 2, the uncertainty parameters $U$444Uncertainty parameter is
defined by Minor Planet Center, on a scale 0-9, 9 being very uncertain. of
2006 RJ2, 2006 BZ8, 2012 YE8, 2017 SV13 according to Minor Planet Center are
5, 0, 5 and 4, respectively. To exclude the influence of uncertainties in the
orbit determination, we generated 1000 clone orbits (include the nominal one)
for the each selected Centaur by standard deviations (1-$\sigma$
uncertainties), which are corresponding to Gaussian distribution around the
nominal orbit. We integrate the nominal and clone orbits over a 40000 yr
period from -10000 yr to 30000 yr around present-day.

Table 1: Nominal orbital elements at JED2458000.5 of 2006 RJ2. The 1-$\sigma$ uncertainty and the Jupiter Tisserand’s parameter are also given. Data taken from website of JPL Small-Body Database Search Engine. The uncertainty parameter is given by the MPC. Full Name |  2006 RJ2
---|---
Epoch (JED) |  2458000.5
Semi-major axis, $a$(au) |  $9.676917955\pm 9.4898\times 10^{-2}$
Eccentricity, $e$ |  $0.7623217118\pm 2.2762\times 10^{-3}$
Inclination, $i$(°) |  $164.5713262\pm 5.9369\times 10^{-3}$
Longitude of the ascending node, $\Omega$(°) |  $191.70837\pm 2.9035\times 10^{-4}$
Argument of perihelion, $\omega$(°) |  $161.421009\pm 4.943\times 10^{-2}$
Mean anomaly, $M$(°) |  $133.2269958\pm 1.9675$
Jupiter Tisserand’s parameter, $T_{J}$ |  $-1.1674$
Uncertainty parameter, $U$ |  5

##  4 Crucial statistical results of numerical simulations

To evaluate this particular resonant state of each candidate, we focus on the
fractions of clones trapped in 1/-1 resonance, currently in 1/-1 resonance and
trapped in 1/-1 resonance over 10000 yr with Saturn, respectively. The crucial
statistical results of numerical simulations include the properties of each
candidate in 1/-1 mean motion resonance with Saturn are given in Table 3. The
following section shows detailed descriptions of dynamical states of the four
Centaurs under study.

2006 RJ2 was discovered on 14 September 2006 and referred to as Centaur. As
shown in Table 1, the Centaur has an orbit with large inclination and
eccentricity. As shown in Figure 1, the nominal orbit of 2006 RJ2 is currently
in a long-term stable 1/-1 resonance with Saturn. In fact, the fourth panel of
Figure 1 shows that libration of the resonant angle $\phi$ is of large
amplitude and is centered on 0°, which confirms that the nominal orbit of 2006
RJ2 is in 1/-1 resonance with Saturn. The oscillation of the resonant angles
for most of the clones of 2006 RJ2 also behaves in a similar way. As shown in
Table 3, 906 clones are trapped in retrograde co-orbit resonance and 752
clones among them keep this resonant state over 10000 yr. Above all, a
majority of (805/1000) clones are shown to be in the co-orbit resonance at
present time. Therefore, these results seem to suggest that 2006 RJ2 is very
likely a retrograde co-orbital Centaur of Saturn, and that finding is not
accidentally due to uncertainties in the orbit determination. From this, it is
probable that 2006 RJ2 is another minor body apart from 2015 BZ509 in
retrograde co-orbit resonance with giant planets. The last panel of Figure 1
shows that the resonance protects the Centaur from close encounters with the
giant planets. In fact, the minimum distance in a close encounter with the
giant planets is always larger than 2 Hill’s radii during the whole
integration timespan. The brief capture in the 2/-5 resonance with Jupiter
after the Centaur exits the 1/-1 resonance is also possible (363/1000 clones).
The path of the Centaur 2006 RJ2 in the Saturn’s rotating frame is shown in
Figure 2. When viewed from above the plane of Saturn’s orbit (Figure 2a), the
motion of 2006 RJ2 has a typical type of trisectrix curve defined by Morais &
Namouni (2017a). Morais & Namouni (2016) have confirmed by the numerical
investigations that this trisectrix motion of 2006 RJ2 can be long-term stable
if it’s influenced only by Saturn.

Figure 1: Properties of the nominal orbit of 2006 RJ2. From top to lower
panel: semimajor axis $a$ and location of exact 1/-1 resonance with Saturn;
eccentricity $e$; inclination $i$; angle for nearly coplanar motion involving
1/-1 resonance with Saturn $\phi$; 2/-5 resonant angle with Jupiter
$\phi_{2}$; argument of pericenter $\omega$; the minimum distance in close
encounters with Jupiter and Saturn within three Hill’s radii of the planet is
represented by red plus and green square, respectively. Figure 2: The path of
the Centaur 2006 RJ2 in the frame that rotates with Saturn. $X$, $Y$, and $Z$
are the cartesian axes of the Saturn’s rotating frame. The $X$-axis points
from Sun to Saturn and $Z$-axis is vertical to Saturn’s orbital plane. Saturn
is presented as the green filled circle. Panel a): the path of 2006 RJ2 in
Saturn’s rotating frame viewed from the $Z$-axis direction. The black dotted
curve is the Saturn’s orbit in the inertial reference frame. Panel b): the
motion of 2006 RJ2 viewed from the $Y$-axis direction. The blue curve
represents the motion of 2006 RJ2 above the Saturn’s orbital plane while the
red curve means the opposite.

Now, we focus on the other three potential retrograde co-orbital Centaurs of
Saturn, 2006 BZ8, 2012 YE8 and 2017 SV13, which have retrograde orbits around
Saturn and may be captured into the 1/-1 resonance within the whole
integration timespan. Table 2 shows their orbital elements and the 1-$\sigma$
uncertainties at JED2458000.5. 2006 BZ8 has been confirmed to be in the 2/-5
resonance with Jupiter (libration around 180°)(Morais & Namouni, 2013a).
Moreover, as listed in Table 3, our numerical statistical results show that
this Centaur has low probability of currently in 1/-1 resonant state with
Saturn (9/1000), which is not mentioned in Morais & Namouni (2013a). Moreover,
as shown in Figure 3, we choose one of the best clones captured in 1/-1
resonant state. Our statistical results indicate that the capture in 1/-1
resonance with Saturn (libration around 0°) within 40000 yr integration
timespan is possible (878/1000 clones), and 38 clones among them are captured
in 1/-1 resonant state over 10000 yr. Morais & Namouni (2013a) suggested that
2006 BZ8 is in Lidov-Kozai secular resonance, we find that 2006 BZ8 is locked
in Lidov-Kozai mechanism (the argument of pericenter $\omega$ librates around
0°) at present and will exit at around 120000 yr${}^{\ref{foot1}}$. The main
dynamical properties of the Lidov-Kozai mechanism for retrograde orbits
($i>90\degr$) outside the mean motion resonances are identical with the
symmetrical prograde orbits with inclination $180-i$. Kozai (1962) studied the
secular evolution of a test partical ($a$ = 3 au) in the main asteroid belt
while considering the effect of Jupiter only. Thomas & Morbidelli (1996) used
the semi-analytical method to study the Lidov-Kozai dynamics around the main
belt ($a$ = 3 au) and the Kuiper belt ($a$ = 45 au) by considering the effects
of the four giant planets. Gronchi & Milani (1999) introduced an example ($a$
= 7 au) to analyze the Lidov-Kozai dynamics on Centaurs. They pointed out that
orbits with moderate value of maximum orbital inclination ($i_{\rm max}$ =
25°), the argument of pericenter $\omega$ librates around 0°or 180°. For even
larger $i_{\rm max}$, there will be new libration regions where
$\omega=90\degr$ or 270°and large eccentricity. But there are not systematic
studies about the role of Lidov-Kozai mechanism on minor bodies near giant
planets until now. We will analyze in detail the effects of Lidov-Kozai
dynamics on retrograde Centaurs near giant planets in our future work. It’s
interesting to note that in the restricted three-body model (Sun-Saturn-2006
BZ8), the Centaur 2006 BZ8 is in 1/-1 resonance with Saturn very well over
10000 yr. It has the similar result in the Sun-Jupiter-2006 BZ8 model where
the Centaur is in 2/-5 resonance with Jupiter. We will try to analyze this
interesting phenomenon in our future work.

2012 YE8 (Figure 4) is currently near the 2/-5 resonance with Jupiter and the
1/-1 resonance with Saturn at the same time. Close approaches to Saturn and
Jupiter may lead the Centaur to enter or exit the resonant states. Since its
orbit is far from the planar case, we consider two possible 1/-1 resonant
angles ($\phi$ and $\phi_{z}$). In our numerical integrations for 2012 YE8, as
listed in Table 3, 140/1000 clones are currently in 1/-1 resonance with Saturn
(libration around 0°), 426 clones are trapped in this resonant state, and 100
clones among them have relative long-term stability. Meanwhile, 268/1000
clones are captured in the 2/-5 resonance with Jupiter (libration around
180°), most often occur between $-10000$ yr and $-5000$ yr and from 15000 yr
to 25000 yr. The argument of pericenter $\omega$ always close to 0°over the
integration timespan, which means a Lidov-Kozai secular resonance555Longer
integrations have been done to confirm these results..

2017 SV13 was just discovered on 17 September 2017 and referred to as Centaur.
Different from other three Centaurs, 2017 SV13 has a near-vertical orbit. As
shown in Figure 5, 2017 SV13 is another suitable candidate in 1/-1 resonance
with Saturn. Almost all of clones (998/1000 clones) are captured in the 1/-1
resonance with Saturn (libration around 0°), and the resonant state of 626
clones among them last beyond 10000 yr. Meanwhile, 434 clones are shown to be
in retrograde co-orbit resonance at present, and clones are always expected in
1/-1 resonant state between 5000 yr and 15000 yr timespan. These results make
2017 SV13 to be another good candidate for Saturn’s retrograde co-orbital
Centaur. After leaving the 1/-1 resonance around 15000 yr, 595/1000 clones
will be captured in 2/-5 resonance with Jupiter (libration around 180°). The
motion is not in Lidov-Kozai secular resonance presently but will be briefly
captured in the Lidov-Kozai mechanism (the argument of pericenter $\omega$
librates around $\pm$90°).${}^{\ref{foot1}}$

Table 2: Nominal orbital elements at JED2458000.5 of 2006 BZ8, 2012 YE8, 2017 SV13. The 1-$\sigma$ uncertainties and the Jupiter Tisserand’s parameter are also given. Data taken from website of JPL Small-Body Database Search Engine. The uncertainty parameters are given by the MPC. Full Name |  2006 BZ8
---|---
Epoch (JED) |  2458000.5
Semi-major axis, $a$(au) |  $9.608985507\pm 6.812\times 10^{-5}$
Eccentricity, $e$ |  $0.8024503954\pm 1.3873\times 10^{-6}$
Inclination, $i$(°) |  $165.318425\pm 1.8448\times 10^{-5}$
Longitude of the ascending node, $\Omega$(°) |  $183.7606180\pm 4.5765\times 10^{-5}$
Argument of perihelion, $\omega$(°) |  $82.46140331\pm 6.4218\times 10^{-5}$
Mean anomaly, $M$(°) |  $134.9664552\pm 1.4364\times 10^{-3}$
Jupiter Tisserand’s parameter, $T_{J}$ |  $-1.027$
Uncertainty parameter, $U$ |  0
Full Name |  2012 YE8
Epoch (JED) |  2458000.5
Semi-major axis, $a$(au) |  $9.38168449\pm 7.1464\times 10^{-2}$
Eccentricity, $e$ |  $0.5914212664\pm 3.037\times 10^{-3}$
Inclination, $i$(°) |  $136.1071898\pm 2.2571\times 10^{-2}$
Longitude of the ascending node, $\Omega$(°) |  $135.3614075\pm 2.6817\times 10^{-2}$
Argument of perihelion, $\omega$(°) |  $35.05260406\pm 3.5294\times 10^{-2}$
Mean anomaly, $M$(°) |  $134.9664552\pm 5.6912\times 10^{-1}$
Jupiter Tisserand’s parameter, $T_{J}$ |  $-1.006$
Uncertainty parameter, $U$ |  5
Full Name |  2017 SV13
Epoch (JED) |  2458000.5
Semi-major axis, $a$(au) |  $9.686282362\pm 3.4318\times 10^{-2}$
Eccentricity, $e$ |  $0.7928247257\pm 7.1768\times 10^{-4}$
Inclination, $i$(°) |  $113.2479445\pm 6.9472\times 10^{-3}$
Longitude of the ascending node, $\Omega$(°) |  $11.70939482\pm 7.4949\times 10^{-4}$
Argument of perihelion, $\omega$(°) |  $343.2582508\pm 1.0814\times 10^{-2}$
Mean anomaly, $M$(°) |  $0.6024677028\pm 4.0766\times 10^{-3}$
Jupiter Tisserand’s parameter, $T_{J}$ |  $-0.119$
Uncertainty parameter, $U$ |  4
Figure 3: A cloned orbit of 2006 BZ8, which is chosen as one of the best
clones captured in 1/-1 resonance. The basic meaning of this figure is the
same as Figure 1. Figure 4: A cloned orbit of 2012 YE8. Since its orbit is far
from the planar case, another 1/-1 resonant angle with Saturn $\phi_{z}$ is
included. The other basic meaning of this figure is the same as Figure 1.
Figure 5: A cloned orbit of 2017 SV13. Since its orbit is far from the planar
case, another 1/-1 resonant angle with Saturn $\phi_{z}$ is included. In the
last panel, the minimum distance in close encounters with Jupiter and Uranus
within three Hill’s radii of the planet is represented by red plus and cyan
cross, respectively. The other basic meaning of this figure is the same as
Figure 1. Table 3: Numerical statistics about the fraction of clones for each
of the four Centaurs under study. The fraction is the ratio of the total
clones in 1/-1 resonance to 1000. To evaluate the resonant state of each
Centaur, we focus on the fractions of clones trapped in resonance (top row),
currently in resonance (middle row), and trapped in resonance over 10000 yr
(bottom row), respectively

.  Centaurs 2006 RJ2 2006 BZ8 2012 YE8 2017 SV13 trapped in 1/-1 resonance
906/1000 878/1000 426/1000 998/1000 currently in 1/-1 resonance 805/1000
9/1000 140/1000 434/1000 captured in 1/-1 resonance over 10000 yr 752/1000
38/1000 100/1000 626/1000

##  5 Phase-space portrait of 1/-1 mean motion resonance

To confirm that the Centaurs may be captured in 1/-1 resonance with Saturn,
this section is devoted to the construction of the theory of the retrograde
mean motion resonance, and the phase-space portrait of 1/-1 mean motion
resonance specifically. Follow the method introduced by Huang et al. (2018) to
analyze the dynamics of 1/-1 resonance, firstly, we introduce the following
retrograde Poincare variables:

$\begin{split}&\Lambda=L,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\lambda=M+\omega-\Omega,\\\ &P=L-G,\ \ \ \ \ \ \ \ \ \ p=-\omega+\Omega,\\\
&Q=G+H,\ \ \ \ \ \ \ \ \ q=\Omega,\\\ \end{split}$ (3)

where $L=\sqrt{a}$, $G=L\sqrt{1-e^{2}}$, $H=G$$\cos$$i$. The canonical
variables $L,M,G,\omega,H,\Omega$ are usually called the Delaunay variables.
These new retrograde Poincare variables are different from these given by
Murray & Dermott (1999). But it is easy to check that the following contact
transformation relationship holds:

$\Lambda\lambda+Pp+Qq=LM+G\omega+H\Omega,$ (4)

So the new retrograde Poincare variables is a set of canonical action-angle
variables. The 1/-1 mean motion resonant Hamiltonian for the circular
restricted three-body model is given by Morbidelli (2002)

$\mathcal{H}=\mathcal{H}_{0}(\Lambda,\Lambda\arcmin)+\epsilon\mathcal{H}_{1}(\Lambda,P,Q,p,q,\lambda-\lambda\arcmin),$
(5)

where

$\mathcal{H}_{0}=-\tfrac{1}{2\Lambda^{2}}+n\arcmin\Lambda\arcmin,$ (6)

is the main term, and $\mathcal{H}_{1}$ is the disturbing term. To study the
dynamics of Hamiltonian (5), we introduce a set of canonical action-angle
variables as follows:

$\begin{split}&S=P,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\sigma=\tfrac{\phi}{2},\\\ &S_{z}=Q,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \sigma_{z}=\tfrac{\phi_{z}}{2},\\\ &N=-2\Lambda+P+Q,\ \ \ \ \ \ \ \
\ \upsilon=-\tfrac{\lambda-\lambda\arcmin}{2},\\\
&\tilde{\Lambda}\arcmin=\Lambda\arcmin+2\Lambda,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \tilde{\lambda}\arcmin=\lambda\arcmin,\\\ \end{split}$ (7)

where $\sigma$ and $\sigma_{z}$ are the critical angles of the retrograde mean
motion resonance. If we restrict our retrograde motion to planar case
$i=180\degr$, the second resonant angle $\phi_{z}$ is neglectable and then
$S_{z}$ or $Q$ is always equal to 0. Besides, it is obvious that $N$ is
constant of motion for $\mathcal{H}$, whose expression is now

$N=-2\Lambda+P+Q=-\sqrt{a}(1+\sqrt{1-e^{2}}),$ (8)

which define a curve in the $a-e$ plane. Thus the Hamiltonian of the planar
case only depends on one critical angle $\sigma$. The new integrable
Hamiltonian reads:

$\mathcal{H}_{\rm
planar}=\mathcal{H}_{0}(N,S)+\epsilon\mathcal{H}_{1}(N,S,2\sigma),$ (9)

Figure 6: The phase-space portrait of the 1/-1 mean motion resonance with
Jupiter associated with the orbital elements of the mysterious asteroid 2015
BZ509, presented in the $e-\phi$ polar space. The variance of color and
contours represent level curves of the $\mathcal{H}_{\rm planar}$ (equation
9). The red scatter diagram present the results of the numerical integration
in the N-body numerical model, while the yellow dashed line denotes the
collision curve with the perturber.

As for the disturbing term $\epsilon\mathcal{H}_{1}$, we numerically evaluate
its value by averaging over all fast angles (Morbidelli, 2002). Then we can
calculate the value of equation (9) on $e-\phi$ plane. We apply this planar
1/-1 resonance semi-analytical model to the first asteroid in 1/-1 resonance
with Jupiter, 2015 BZ509666Orbital elements of 2015 BZ509 taken from JPL
Small-Body Database Search Engine; https://ssd.jpl.nasa.gov/sbdb_query.cgi,
retrieved on January 2018., which has a nearly coplanar orbit. We can predict
its motion in numerical model very well (Figure 6).

Figure 7: The phase-space portrait of the 1/-1 mean motion resonance with
Saturn associated with the orbital elements of the Centaur 2006 RJ2, presented
in the $e-\phi$ polar space. The basic meaning of this figure is the same with
Figure 6. But we confine the plot range of the eccentricity to (0.7,0.78) for
better visual presentation. Figure 8: The phase-space portrait of the 1/-1
mean motion resonance with Saturn associated with the orbital elements of the
Centaur 2006 BZ8. The basic meaning of this figure is the same with Figure 6.
But we confine the plot range of the eccentricity to (0.62,0.75) for better
visual presentation.

As shown in Tables 1 and 2, the inclinations of 2006 RJ2 and 2006 BZ8 are both
close to 180°. We can apply this planar semi-analytical model to these two
Centaurs to evaluate their dynamical motions. The $e-\phi$ diagram for 2006
RJ2 is presented in Figure 7. We notice that the motion of 2006 RJ2 in the
numerical model is predicted pretty accurately by the planar semi-analytical
model. The motion of the Centaur is just around the ideal equilibrium point of
the 1/-1 resonance. The curves given by the numerical integration correspond
to the large-amplitude libration as defined by Murray & Dermott (1999). The
phase-space portrait of 2006 BZ8 after capturing in 1/-1 resonance also
presented in Figure 8. The motion in numerical model is not predicted pretty
accurately as 2006 RJ2 due to the non-negligible impact of the Jupiter. But it
also belongs to the scale of large-amplitude libration around the ideal
equilibrium point. As shown in Figures 6-8, the collision curve slices the
entire phase space into different region. This work further verify the 1/-1
resonant state of 2006 RJ2 and 2006 BZ8. We will research on the theory of the
retrograde mean motion resonance that can apply to three-dimensional cases in
future work.

##  6 Conclusions

Morais & Namouni (2016) predicted stable co-orbit resonances exist at all
inclinations, but an asteroid in stable retrograde co-orbit resonance with
planet has not yet been discovered until the first identified asteroid in 1/-1
resonance with Jupiter, 2015 BZ509 (Wiegert et al., 2017). In this paper, we
identify four potential retrograde co-orbital Centaurs of Saturn. The
statistical results of 1000 clones (include nominal orbit) of each candidate
by numerical simulations are given in Table 3. According to our statistical
results, 2006 RJ2 is the best candidate to be currently in a 1/-1 mean motion
resonance with Saturn, as another minor body found in retrograde co-orbit
resonance with the planet after 2015 BZ509. 2017 SV13 is another important
potential candidate whose clones are are always captured in this particular
resonant state during the whole integration timespan. Moreover, 2012 YE8 and
2006 BZ8 are also Centaurs of interest but their current and long-term 1/-1
resonant state with Saturn is less likely.

Our findings suggest that small bodies in retrograde co-orbit resonance with
giant planets may be more common in our solar system than previously expected.
The existence of these mysterious small bodies may give us a kind of unique
perspective to study the co-orbit motion in our solar system. As for the small
bodies on planet-crossing orbits such as Centaurs, Fernandez et al. (2016);
Bailey & Malhotra (2009); Namouni & Morais (2015) suggested that those locked
in stable resonances have longer lifetimes than others. That is to say, the
Centaurs in retrograde co-orbit resonance we identified here, 2006 RJ2, 2006
BZ8, 2012 YE8, 2017 SV13 may belong to the most dynamically stable objects
among Centaurs. This is significant for the research of the Centaurs’ origin
and dynamical evolution.

The dynamics of the planar retrograde co-orbit resonance is studied through
the averaged Hamiltonian. We introduce an useful integrable approximation for
planar 1/-1 resonance and confirm the 1/-1 resonant states by the comparison
of the phase-space portrait between the semi-analytical model and the
numerical results. The planar 1/-1 resonance semi-analytical model can predict
the motion of 2015 BZ509 (its inclination is close to 180°) quite well. The
real motions of the 2006 RJ2 and 2006 BZ8 (after capturing in 1/-1 resonance)
in the solar system are just around the ideal equilibrium points of the 1/-1
resonance given by the planar semi-analytical model. These results further
verify our findings given by the numerical simulations. Namouni & Morais
(2017) confirmed that the efficiency of co-orbit capture at large retrograde
inclination is an intrinsic feature of co-orbital resonance. Therefore, this
planar semi-analytical model can provide a good approximation to the real
motion of most small bodies in 1/-1 resonance with the planets.

Centaurs are small bodies of the solar system on the long-term unstable orbits
with a semi-major axis between those of the outer planets. As a result, these
retrograde Centaurs wander radially between the outer planets have lots of
”chance” to be captured in co-orbit resonance because of the frequent close
encounters with giants planets.

The conclusions of our present study encourage the search for such minor
bodies on a larger scale to analyze the origins and evolution of these
mysterious objects. Our results imply that retrograde mean motion resonance
plays a crucial role in the dynamical evolution of the Centaurs on retrograde
orbits. In this regard, the findings of this paper are useful for
understanding the origin and dynamical evolution of the Centaurs and
Damocloids on retrograde orbits.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grant No.11772167)

## References

  * Bailey & Malhotra (2009) Bailey, B. L. & Malhotra, R. 2009, Icarus, 203, 155
  * Chambers (1999) Chambers, J. E. 1999, Monthly Notices of the Royal Astronomical Society, 304, 793
  * Chen et al. (2016) Chen, Y.-T., Lin, H. W., Holman, M. J., et al. 2016, The Astrophysical Journal Letters, 827, L24
  * Connors et al. (2011) Connors, M., Wiegert, P., & Veillet, C. 2011, Nature, 475, 481, identifier: BFnature10233
  * Emel’yanenko & Bailey (1998) Emel’yanenko, V. V. & Bailey, M. E. 1998, Monthly Notices of the Royal Astronomical Society, 298, 212
  * Fernandez & Gallardo (1999) Fernandez, J. A. & Gallardo, T. 1999, in IAU Colloq. 173: Evolution and Source Regions of Asteroids and Comets, ed. J. Svoren, E. M. Pittich, & H. Rickman, 327
  * Fernandez et al. (2016) Fernandez, J. A., Gallardo, T., & Young, J. D. 2016, Monthly Notices of the Royal Astronomical Society, 461, 3075
  * Fouchard et al. (2014) Fouchard, M., Rickman, H., Froeschlé, C., & Valsecchi, G. B. 2014, Icarus, 231, 99
  * Greenstreet et al. (2012) Greenstreet, S., Gladman, B., Ngo, H., Granvik, M., & Larson, S. 2012, The Astrophysical Journal Letters, 749, L39
  * Gronchi & Milani (1999) Gronchi, G. F. & Milani, A. 1999, A&A, 341, 928
  * Heppenheimer & Porco (1977) Heppenheimer & Porco. 1977, Icarus, 30, 385
  * Hou et al. (2014) Hou, X. Y., Scheeres, D. J., & Liu, L. 2014, Monthly Notices of the Royal Astronomical Society, 437, 1420
  * Huang et al. (2018) Huang, Y., Li, M., Li, J., & Gong, S. 2018, The Astronomical Journal, 155, 262
  * Jewitt (2005) Jewitt, D. 2005, The Astronomical Journal, 129, 530
  * Jewitt & Sheppard (2003) Jewitt, D. C. & Sheppard, S. S. 2003, Nature, 423, 261
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Levison et al. (2001) Levison, H. F., Dones, L., & Duncan, M. J. 2001, The Astronomical Journal, 121, 2253
  * Levison et al. (2006) Levison, H. F., Duncan, M. J., Dones, L., & Gladman, B. J. 2006, Icarus, 184, 619
  * Morais & Namouni (2017a) Morais, H. & Namouni, F. 2017a, Nature, 543, 635
  * Morais & Namouni (2013a) Morais, M. H. M. & Namouni, F. 2013a, Monthly Notices of the Royal Astronomical Society: Letters, 436, L30
  * Morais & Namouni (2013b) Morais, M. H. M. & Namouni, F. 2013b, Celestial Mechanics and Dynamical Astronomy, 117, 405, identifier: 9519
  * Morais & Namouni (2016) Morais, M. H. M. & Namouni, F. 2016, Celestial Mechanics and Dynamical Astronomy, 125, 91
  * Morais & Namouni (2017b) Morais, M. H. M. & Namouni, F. 2017b, Monthly Notices of the Royal Astronomical Society: Letters, 472, L1
  * Morbidelli (2002) Morbidelli, A. 2002, Modern celestial mechanics : aspects of solar system dynamics
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics (Cambridge: Cambridge Univ. Press)
  * Namouni & Morais (2017) Namouni, F. & Morais, H. 2017, Computational and Applied Mathematics
  * Namouni & Morais (2015) Namouni, F. & Morais, M. H. M. 2015, Monthly Notices of the Royal Astronomical Society, 446, 1998
  * Narita et al. (2009) Narita, N., Sato, B., Hirano, T., & Tamura, M. 2009, Publications of the Astronomical Society of Japan, 61, L35
  * Nurmi et al. (2002) Nurmi, P., Valtonen, M. J., Zheng, J. Q., & Rickman, H. 2002, Monthly Notices of the Royal Astronomical Society, 333, 835
  * Olsson-Steel (1988) Olsson-Steel, D. 1988, Monthly Notices of the Royal Astronomical Society, 234, 389
  * Pollack et al. (1979) Pollack, J., A.Burns, J., & E.Tauber, M. 1979, Icarus, 37, 587
  * Thomas & Morbidelli (1996) Thomas, F. & Morbidelli, A. 1996, Celestial Mechanics and Dynamical Astronomy, 64, 209
  * Triaud et al. (2010) Triaud, A. H. M. J., Collier Cameron, A., Queloz, D., et al. 2010, Astronomy & Astrophysics, 524, A25, identifier: aa14525-10
  * Volk & Malhotra (2013) Volk, K. & Malhotra, R. 2013, Icarus, 224, 66
  * Wiegert et al. (2017) Wiegert, P., Connors, M., & Veillet, C. 2017, Nature, 543, 687, identifier: nature22029
  * Wiegert & Tremaine (1999) Wiegert, P. & Tremaine, S. 1999, Icarus, 137, 84
  * Winn et al. (2009) Winn, J. N., Johnson, J. A., Albrecht, S., et al. 2009, The Astrophysical Journal, 703, L99