20XX Vol. X No. XX, 000–000

11institutetext:  Key Laboratory of Planetary Sciences, Purple Mountain
Observatory, Chinese Academy of Sciences, Nanjing 210023, China;
<EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: School of Astronomy and Space Science, University of Science
and Technology of China, Hefei 230026, China
33institutetext: Center for Excellence in Comparative Planetology, Chinese
Academy of Sciences, Heifei 230026, China
Received 20XX Month Day; accepted 20XX Month Day

# Shape model and rotation acceleration of (1685) Toro and (85989) 1999 JD6
from optical observations

Jun Tian 1122 Hai-Bin Zhao 112233 Bin Li 1122

###### Abstract

The Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect is a net torque caused
by solar radiation directly reflected and thermally re-emitted from the
surface of small asteroids and is considered to be crucial in their dynamical
evolution. By long-term photometric observations of selected near-Earth
asteroids, it’s hoped to enlarge asteroid samples with a detected YORP effect
to facilitate the development of a theoretical framework. Archived light-curve
data are collected and photometric observations are made for (1685) Toro and
(85989) 1999 JD6, which enables measurement of their YORP effect by inverting
the light curve to fit observations from a convex shape model. For (1685)
Toro, a YORP acceleration $\upsilon=(3.2\pm 0.3)\times 10^{-9}\ \rm{rad\cdot
d^{-2}}(1\sigma\ error)$ is updated, which is consistent with previous YORP
detection based on different light-curve data; for (85989) 1999 JD6, it is
determined that the sidereal period is $7.667749\pm 0.000009$ h, the rotation
pole direction locates is at $\lambda=232\pm 2^{\circ},\ \beta=-59\pm
1^{\circ}$, the acceleration is detected to be $\upsilon=(2.4\pm 0.3)\times
10^{-8}\ \rm{rad\cdot d^{-2}}(1\sigma\ error)$ and in addition to obtaining an
excellent agreement between the observations and model. YORP should produce
both spin-up and spin-down cases. However, including (85989) 1999 JD6, the
$\rm{d}\omega/\rm{d}t$ values of eleven near-Earth asteroids are positive
totally, which suggests that there is either a bias in the sample of YORP
detections or a real feature needs to be explained.

###### keywords:

mechanisms: YORP effect – methods: data analysis – methods: observational –
techniques: photometric – minor planets, asteroids: individual: (85989) 1999
JD6 – minor planets, asteroids: individual: (1685) Toro

##  1 Introduction

The Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect is one of the
mechanisms of the long-term dynamical evolution of small asteroids in the
Solar system. It is a net torque caused by two principal mechanisms: the
anisotropic reflection of sunlight and thermal emission of an asteroid. The
YORP effect, also called non-gravitational effect, was introduced by Rubincam
(2000). Since then, YORP has been recognized as an important mechanism of
physical and dynamical evolution. YORP can change the rotation rate and spin-
axis obliquity and affects the distribution of rotation rates and obliquities,
especially for asteroids in the size range from $\sim$1 m to $\sim$40 km
(Bottke et al. 2006).

Asteroids can be accelerated to rotational fission (Jacobson et al. 2016),
mass shedding (Scheeres 2015), reshape (Cheng et al. 2021) and creat asteroid
binaries and pairs (Pravec et al. 2010; Margot et al. 2015). Asteroids also
can be decelerated to a tumbling state. Lowry and Taylor have directly
detected that the rotation period of asteroid (54509) YORP is decreasing
continuously through the light-curve data, which is consistent with the
prediction of the YORP theory. To date, the YORP effect has been detected on
ten near-Earth asteroids (NEAs): (54509) YORP, (1862) Apollo, (1620)
Geographos, (3103) Eger, (25143) Itokawa, (161989) Cacus, (101955) Bennu,
(68346) 2001 KZ66, (1685) Toro, and (10115) 1992 SK (Lowry et al. 2007;
Kaasalainen et al. 2007; Lowry et al. 2014; Taylor et al. 2007; Ďurech et al.
2008; Ďurec et al. 2008; Ďurech et al. 2012, 2018; Nolan et al. 2019;
Hergenrother et al. 2019; Zegmott et al. 2021; Ďurech et al. 2022). Crucially,
all of the detections are rotational accelerations. While that the YORP effect
can accelerate or decelerate the rotation rate of asteroids is not only a
prediction given by the YORP theory(Rubincam 2000), but also a conclusion
drawn through simulation experiments (Rossi et al. 2009; Golubov & Krugly
2012). Although the tangential-YORP (TYORP) theory – a component of the recoil
force parallel to the surface caused by re-emission of absorbed solar light
from centimeter- to decimeter-sized structures on the asteroid’s surface, has
been now proposed to explain the lack of decelerating rotation asteroids
(Golubov et al. 2014), in order to explain the real feature and further verify
the YORP theory, it is still necessary to continuously enlarge the asteroid
samples with a detected YORP effect.

NEAs tend to be more clearly affected by the YORP effect, as they tend to be
small and close proximity to the Sun. Therefore, the photometric data of NEAs
and the published light curves are being collected, which include, but are not
limited to, the Light Curve Database (LCDB) and the photometric data obtained
by China Near-Earth Object Survey Telescope (CNEOST), to screen for NEAs. In
order to search for NEAs whose rotations are decelerated by YORP effect,
special attention is also paid to those targets whose rotation periods are
between 6 and 9 hours, as they can be observed full rotation in one night
which is beneficial to obtain accurate rotation periods and maximize the
ability of the light-curve inversion model to detect YORP.

In this paper, we will present a shape model and spin-state analysis of two
NEAs, (1685) Toro and (85989) 1999 JD6. In Section 2, observing campaign of
(1685) Toro is described, and the spin state and the shape model are presented
as well as the approach to detect YORP-induced rotational acceleration. In
Section 3, the collected light-curve data, results of the sidereal period and
rotational pole described, and a shape model and the value of YORP rotational
acceleration by the YORP model are also given. Section 4 provides a discussion
of the results and their implications, and a summary of the main conclusion is
in Section 5.

##  2 (1685) Toro

(1685) Toro is an Apollo-type NEA discovered by C.A.Wirtanen at Mount Hamilton
on July 17,1948. It’s absolute magnitude is $H=14.48\pm 0.13$ mag, and the
assuming slope parameter $G=0.24\pm 0.11$ (Ďurech et al. 2022; Warner et al.
2009). A tentative YORP acceleration by Ďurech et al. (2018) was based on a
data set from apparitions from 1972 to 2016, then its YORP acceleration was
updated by Ďurech et al. (2022) through adding new photometric observations
for 2018, 2020, and 2021. On the basis of the photometric data from 1972 to
2016, light-curve data in 2020 is collected and continuous photometric
observations of (1685) Toro in February 2021 is also carried out. With all of
these light-curve datasets, its YORP value is further verified.

###  2.1 Optical light-curve

A telescope of diameter 80 cm in Yaoan Station of China is used to carry out
light curve observation of Toro. The telescope is equipped with a
4128$\times$4104 Balor 17F-12 sCMOS camera, providing a field of view (FOV) of
$1.6^{\circ}\times 1.6^{\circ}$ and a pixel scale of $1.4^{\prime\prime}\times
1.4^{\prime\prime}$. The Johnson R filter is chosen for the imaging
observations over eight nights from 11 to 18 February 2021. Each exposure
takes 60 s with cadence 22 s. The raw image data is reduced following standard
procedures, including bias, dark subtraction and flat-fielding correction.
Eight light curves are then obtained and are labeled with IDs 9-16 in Table 1.
New continuous photometric observations of Toro also including previously
published data at the Palmer Divide Station in 2020, the light curves are
labeled with IDs 1-8 in Table 1. It includes (1685) Toro’s distance from the
Sun($r$) and Earth ($\Delta$). The solar phase angle($\alpha$) at the middle
point of the observation interval is given, as well as the geocentric ecliptic
longitude ($\lambda_{0}$) and ecliptic latitude ($\beta_{0}$) of the asteroid,
the Universal Time (UT) “Date” at the beginning of the night, the apparent
peak-to-peak “Amplitude” and “Total” length of the light curve in Table 1. The
other processed light curves are obtained from the ALCDEF database (Warner &
Stephens 2020) and the Database of Asteroid Models from Inversion
Techniques111https://astro.troja.mff.cuni.cz/projects/damit/ (DAMIT) web page.

Table 1: Aspect data for new observations of (1685) Toro. Each light curve has
a numerical “ID” listed. The data sources are given in “Reference” : (1)Warner
& Stephens (2020); (2) this work. Observing facility key: 40-cm SC: 40 cm
Schmidt-telescope; 80-cm, Balor 17F-12: 80 cm telescope with the Balor 17F-12
sCMOS camera.

ID | UT Date | $r$ | $\Delta$ | $\alpha$ | $\lambda$ | $\beta$ | Amplitude | Total | Observing | Reference
---|---|---|---|---|---|---|---|---|---|---
| [yyyy-mm-dd] | [AU] | [AU] | [∘] | [∘] | [∘] | [mag] | [h] | facility |
1 | 2020-06-15 | 0.593 | 1.317 | 47.6 | 337.5 | 5.7 | 0.99 | 2.1 | 40-cm SC | 1
2 | 2020-06-17 | 0.572 | 1.304 | 48.2 | 339.3 | 6.2 | 0.96 | 2.2 | 40-cm SC | 1
3 | 2020-06-18 | 0.562 | 1.297 | 48.5 | 340.2 | 6.5 | 0.72 | 2.3 | 40-cm SC | 1
4 | 2020-06-19 | 0.552 | 1.290 | 48.9 | 341.2 | 6.9 | 0.75 | 2.2 | 40-cm SC | 1
5 | 2020-06-21 | 0.532 | 1.277 | 49.6 | 343.2 | 7.5 | 0.80 | 2.3 | 40-cm SC | 1
6 | 2020-06-22 | 0.523 | 1.270 | 50.0 | 344.3 | 7.8 | 0.95 | 2.3 | 40-cm SC | 1
7 | 2020-06-23 | 0.513 | 1.263 | 50.5 | 345.3 | 8.2 | 0.53 | 2.5 | 40-cm SC | 1
8 | 2020-06-27 | 0.477 | 1.236 | 52.3 | 349.8 | 9.7 | 0.79 | 2.5 | 40-cm SC | 1
9 | 2021-02-11 | 0.938 | 1.588 | 35.4 | 211.0 | -15.8 | 1.21 | 5.5 | 80-cm, Balor 17F-12 | 2
10 | 2021-02-12 | 0.933 | 1.593 | 35.0 | 211.0 | -16.0 | 1.00 | 4.8 | 80-cm, Balor 17F-12 | 2
11 | 2021-02-13 | 0.928 | 1.598 | 34.6 | 211.0 | -16.1 | 0.91 | 6.2 | 80-cm, Balor 17F-12 | 2
12 | 2021-02-14 | 0.923 | 1.603 | 34.2 | 210.9 | -16.3 | 1.39 | 5.8 | 80-cm, Balor 17F-12 | 2
13 | 2021-02-15 | 0.918 | 1.608 | 33.8 | 210.9 | -16.5 | 0.99 | 6.0 | 80-cm, Balor 17F-12 | 2
14 | 2021-02-16 | 0.913 | 1.613 | 33.4 | 210.8 | -16.6 | 1.02 | 5.8 | 80-cm, Balor 17F-12 | 2
15 | 2021-02-17 | 0.909 | 1.618 | 33.0 | 210.7 | -16.8 | 1.20 | 6.1 | 80-cm, Balor 17F-12 | 2
16 | 2021-02-18 | 0.904 | 1.623 | 32.5 | 210.7 | -16.9 | 1.08 | 6.3 | 80-cm, Balor 17F-12 | 2

###  2.2 Shape model and YORP rotational acceleration

The constant torque provided by the YORP effect produces a linear change in
the rotation rate. Therefore, the spin-state analysis and optimization require
investigation of the precise timing of the light curve. To achieve the best
match between the observed data and the model, a free parameter
$\upsilon\equiv d\omega/dt$ is added to the convex inversion in the YORP
model, which describes the change of rotation rate $\omega$ and is optimized
during the light-curve inversion together with the shape and spin parameters.
This approach is the same as what used by Ďurech et al. (2018). If a non-zero
$\upsilon$ provides a significantly better fit than $\upsilon=0$, it is
interpreted as detecting rotation acceleration or deceleration. The parameter
is a linear change of the rotation rate in time $\upsilon\equiv d\omega/dt$ or
a quadratic change of the rotation phase in time. The rotation phase in
radians $\varphi(t)$ can be expressed for any given time as:

$\varphi(t)=\varphi(T_{0})+\omega(t-T_{0})+\frac{1}{2}\upsilon(t-T_{0})^{2},$
(1)

where:

$t$: the time of observation(JD),

$T_{0}$ : the epoch from which the model is propagated,

$\varphi(T_{0})$ : initial rotation phase in radians,

$\omega$ : rotation rate in $\rm{rad\cdot d^{-1}}$; $\omega\equiv 2\pi/P$, $P$
is rotation period in days,

$\upsilon$ : the change of rotation rate in $\rm{rad\cdot d^{-2}}$;
$\upsilon\equiv d\omega/dt$ (the YORP strength).

The YORP model iteratively converges to best-fit parameters to minimize the
difference between the observed and modeled light curves. These spin-state
parameters published by Ďurech et al. (2018) are $P=10.19782\pm 0.00003$ h
(for JD 2441507.00000), $(\lambda,\beta)=(71\pm 10^{\circ},\ -69\pm
5^{\circ})$ $(3\sigma\ \rm{errors})$, and $\upsilon=3.0\times 10^{-9}\
\rm{rad\cdot d^{-2}}$. Our study uses these spin-state parameters as input
values for further optimization of the YORP model, which is reasonable and
robust as verified by independent scan and optimization calculations.
Consequently, the updated values ($1\sigma$ errors) are given as
follows:$P=10.197827\pm 0.000002$ h (for JD 2441507.0),
$(\lambda,\beta)=(60\pm 4^{\circ},-72\pm 2^{\circ})$, $\upsilon=(3.2\pm
0.3)\times 10^{-9}\ \rm{rad\cdot d^{-2}}$, which are further optimized by the
YORP model as input values. In the article uncertainties are estimated by a
bootstrap method. For (1685) Toro, 98 light curves are randomly selected from
all 115 light curves as input for optimizing the shape model and spin-state
parameters, then the above process is repeated for 8000 times. Uncertainly are
estimated from the distributions of there parameters given in Fig.1.

Figure 1: The uncertainties of the four parameters are calculated by the
bootstrap method. The four parameters are the spin axis direction in ecliptic
latitude $\beta$ and longitude $\lambda$, the sidereal period and the YORP
acceleration. The black dotted line is the Gaussian fitting curve of these
parameter distributions.

These spin parameters are also given by Ďurech et al. (2022) ($1\sigma$
error): $\upsilon=(3.3\pm 0.3)\times 10^{-9}\ \rm{rad\cdot d^{-2}}$ with
period $P=10.197826\pm 0.000002$ h (for JD 2441507.0), pole direction $(75\pm
3^{\circ},-69\pm 1^{\circ})$. Except for the pole’s ecliptic
longitude($\lambda$) , our results are in good agreement with them. It might
be a different light-curve data set used for the inversion model and different
weights of individual light curves. Nonetheless, they all give a well-matched
YORP rotational acceleration and sidereal period, and YORP rotational
acceleration indicates a weak but robust YORP strength. The convex shape model
is shown in Fig.2 and Fig.3 shows four light curves in which the difference
between the two models is the largest.

Figure 2: Shape model of asteroid (1685) Toro shown from equatorial level
(y-axis, left and x-axis, center , 90∘ apart) and pole-on (z-axis, right).The
model’s z-axis is aligned with the rotation axis and axis of maximum inertia.
Figure 3: Example light curves (black dots) of (1685) Toro shown with the
synthetic light curves produced by the best YORP model(red dashed curves) and
the best constant-period model(dotted black curves). The geometry of
observation is described by the aspect angle $\theta$, the solar aspect angle
$\theta_{0}$, and the solar phase angle $\alpha$. A full set of light curves
are provided in Fig. 11

To have a realistic and independent estimate of the uncertainty of YORP
acceleration, the same approach as Vokrouhlický et al. (2011) and Polishook
(2014) is adopted and $\chi^{2}$ for different fixed values of $\upsilon$ is
computed (all other parameters are optimized). Then the 1$\sigma$ uncertainty
interval of $\upsilon$ is defined so have been that $\chi^{2}$ increases by a
factor of $1+\sqrt{2/\nu}$, where $\nu$ is the number of degrees of freedom
($\nu\ \sim$ 5000 in case of Toro). The computed $\chi^{2}$ for different
fixed values of $\upsilon$ is shown in Fig. 4.

Figure 4: Dependence of the goodness of the fit measured by the $\chi^{2}$ on
the YORP parameter $\upsilon$ for asteroid (1685) Toro. The dashed curve is a
quadratic fit of the data points. The horizontal dotted line indicates a 2.0%
increase in the $\chi^{2}$ , which defines 1$\sigma$ uncertainty interval of
$\pm 0.3\times 10^{-8}\ \rm{rad\cdot d^{-2}}$ given the number of degress of
freedom.

Across the entire set of light-curve data, we identified a weak but reliable
YORP value. To detect the YORP effect of Toro, the observations in 1972 and
1988 are crucial. However, confirming the accuracy of the data is a challenge
due to the lack of independent observations of the same apparition in 1972 and
1988. Moreover, the present light-curve data do not provide adequate
constraints on the ecliptic longitude of the pole direction. Therefore,
photometric observation should always be needed in the future to further
confirm this YORP detection and refine the spin parameters.

##  3 (85989) 1999 JD6

(85989) 1999 JD6 is both an NEA of the Aten class and a potentially hazardous
asteroid (PHA) and was discovered on May 12, 1999 at Anderson Mesa by Lowell
Observatory Near-Earth-Object Search(LONEOS) project in Flagstaff (Tichy et
al. 1999). Photometric and radar observations indicate it is a contact-binary
asteroid belonging to Barbarian asteroids. The visible spectrum data showed
that it is a K-type asteroid (Binzel et al. 2001; de León et al. 2010), and
later it is identified as an L-type asteroid with the aid of the near-infrared
spectrum. Other physical characteristics have been determined, such as albedo
0.05-0.11 (Campins et al. 2009; Thomas et al. 2011; Reddy et al. 2012; Mainzer
et al. 2014; Nugent et al. 2016), maximum diameter 2 km (Marshall et al.
2015), rotation period $7.6638\pm 0.0001$ h (Szabó et al. 2001). Given the
foregoing, is the shape model and spin-state will be given below with light-
curve data.

###  3.1 Optical light-curve

The optical light-curve data for (85989) 1999 JD6 span twenty years. The first
light curves of this asteroid came from the four nights of July 2,3,5, and 6,
2000, and they were observed by a 59 cm Schmidt-telescope. Light curves were
published in VizieR222https://cdsarc.cds.unistra.fr/ (Szabó et al. 2001). In
addition , light-curve data also includes previously published photometry
obtained at the Palmer Divide Station (Warner 2014, 2015, 2018; Warner &
Stephens 2019, 2020). The processed light curves are retrieved from the
Asteroid Light-curve Data Exchange Format (ALCDEF) database (Warner et al.
2011a).

The observational circumstances of the data used in this paper are summarized
in Table 2. It includes (85989) 1999 JD6’s distance from the Sun($r$) and
Earth ($\Delta$). The solar phase angle($\alpha$) at the middle point of the
observation interval is given, as well as the geocentric ecliptic longitude
($\lambda_{0}$) and ecliptic latitude ($\beta_{0}$) of the asteroid, the
Universal Time (UT) “Date” at the beginning of the night, the apparent peak-
to-peak “Amplitude” and “Total” length of the light curve in Table 2.

Table 2: A list of optical light curves of asteroid (85989) 1999 JD6 used in
this study. Each light curve has a numerical “ID” listed. The data sources are
given in “Reference”: (3) Szabó et al. (2001); (4) Warner (2014); (5) Warner
(2015); (6) Warner (2018); (7) Warner & Stephens (2019); (8) Warner & Stephens
(2020). Observing facility key: 59-cm ST, 59 cm Schmidt-telescope, and PDS,
Palmer Divide Station(California, USA).

ID | UT Date | $r$ | $\Delta$ | $\alpha$ | $\lambda$ | $\beta$ | Amplitude | Total | Observing | Reference
---|---|---|---|---|---|---|---|---|---|---
| [yyyy-mm-dd] | [AU] | [AU] | [∘] | [∘] | [∘] | [mag] | [h] | facility |
1 | 2000-07-02 | 0.379 | 1.337 | 27.6 | 266.7 | 34.9 | 1.22 | 2.6 | 59-cm ST | 3
2 | 2000-07-03 | 0.378 | 1.336 | 27.8 | 266.2 | 34.9 | 1.22 | 2.2 | 59-cm ST | 3
3 | 2000-07-05 | 0.371 | 1.323 | 29.6 | 261.9 | 34.0 | 1.24 | 5.2 | 59-cm ST | 3
4 | 2000-07-06 | 0.370 | 1.321 | 29.8 | 261.4 | 33.9 | 1.24 | 3.6 | 59-cm ST | 3
5 | 2014-05-20 | 0.533 | 1.394 | 35.9 | 275.9 | 42.7 | 1.18 | 5.5 | PDS | 4
6 | 2014-05-21 | 0.530 | 1.397 | 35.4 | 274.9 | 42.9 | 1.15 | 5.6 | PDS | 4
7 | 2014-05-22 | 0.528 | 1.399 | 34.9 | 273.9 | 43.1 | 1.04 | 3.1 | PDS | 4
8 | 2015-06-07 | 0.575 | 1.336 | 45.1 | 323.1 | 22.0 | 0.92 | 3.2 | PDS | 5
9 | 2015-06-11 | 0.526 | 1.317 | 44.9 | 324.2 | 22.8 | 1.09 | 3.4 | PDS | 5
10 | 2015-06-12 | 0.514 | 1.312 | 44.9 | 324.5 | 22.9 | 1.07 | 3.5 | PDS | 5
11 | 2015-06-14 | 0.490 | 1.302 | 44.8 | 325.1 | 23.4 | 1.13 | 3.4 | PDS | 5
12 | 2015-06-15 | 0.477 | 1.297 | 44.8 | 325.4 | 23.6 | 1.18 | 3.1 | PDS | 5
13 | 2018-06-01 | 0.517 | 1.180 | 58.8 | 167.5 | 41.2 | 1.41 | 4.0 | PDS | 6
14 | 2018-06-02 | 0.527 | 1.188 | 58.2 | 168.9 | 40.6 | 1.10 | 4.4 | PDS | 6
15 | 2018-06-03 | 0.536 | 1.195 | 57.6 | 170.1 | 39.9 | 1.24 | 4.8 | PDS | 6
16 | 2018-06-04 | 0.546 | 1.202 | 57.0 | 171.3 | 39.3 | 1.25 | 4.4 | PDS | 6
17 | 2019-06-03 | 0.522 | 1.437 | 29.3 | 269.2 | 41.1 | 0.89 | 4.0 | PDS | 7
18 | 2019-06-04 | 0.520 | 1.438 | 28.9 | 268.0 | 41.1 | 0.96 | 4.1 | PDS | 7
19 | 2019-06-05 | 0.518 | 1.439 | 28.5 | 266.8 | 41.1 | 1.06 | 4.0 | PDS | 7
20 | 2019-06-06 | 0.516 | 1.439 | 28.2 | 265.6 | 41.1 | 1.13 | 4.5 | PDS | 7
21 | 2019-06-07 | 0.514 | 1.440 | 28.0 | 264.4 | 41.0 | 1.17 | 4.5 | PDS | 7
22 | 2020-06-18 | 0.445 | 1.231 | 51.6 | 337.8 | 20.0 | 0.85 | 2.7 | PDS | 8
23 | 2020-06-19 | 0.432 | 1.224 | 52.0 | 338.6 | 20.1 | 1.39 | 3.2 | PDS | 8
24 | 2020-06-20 | 0.420 | 1.217 | 52.3 | 339.3 | 20.2 | 1.12 | 2.5 | PDS | 8
25 | 2020-06-21 | 0.408 | 1.211 | 52.6 | 340.1 | 20.3 | 0.92 | 1.9 | PDS | 8
26 | 2020-06-22 | 0.395 | 1.204 | 53.0 | 341.0 | 20.5 | 1.28 | 3.3 | PDS | 8
27 | 2020-06-23 | 0.383 | 1.197 | 53.5 | 341.8 | 20.6 | 1.20 | 3.3 | PDS | 8

###  3.2 Sidereal period and pole search

The convex inversion program starts from initial values , which includes an
initial period, finds the local minimum of $\chi^{2}$, and obtains the
corresponding shape model and period solution. The method described in
Kaasalainen et al. (2001), a range of period values and six initial poles for
each trial period are scanned across, from which the lowest $\chi^{2}$ value
is selected. In consideration of the previously reported periods for (85989)
1999 JD6 (Szabó et al. 2001; Warner 2014, 2015, 2018; Polishook & Brosch 2008;
Warner & Stephens 2019, 2020) and saving computer calculation time, scanned
period value range from 7.0 to 8.0 h. The search result of the best period
indicates the initial rotational period is at 7.667714 h (Fig. 5). This value
is utilized in the next step as an optimal starting point for the subsequent
optimization process.

Figure 5:  This plot shows the best sidereal period (P = 7.667714 h) with the
local minimum of $\chi^{2}$. The found best sidereal period found is used as
the input value to search for the pole direction.

The next step in the shape modeling procedure is to search for the pole
direction with a rotational period of 7.667714 h as an input value and to
determine a best-fitting convex shape for 1999 JD6. The convex inversion
techniques described by Kaasalainen & Torppa (2001) and Kaasalainen et al.
(2001) is applied in the shape model. The first thing to do is to set up a
grid of pole positions covering the entire celestial sphere with a resolution
of $3^{\circ}\times 3^{\circ}$. The $\chi^{2}$ is calculated to fit the light
curves at each fixed pole position and the best sidereal period, and the
corresponding of the minimum of $\chi^{2}$ value is recorded. At this stage,
the initial epoch $T_{0}$ and the initial rotation phase $\varphi(T_{0})$, are
needed for the transformation between vectors $\mathbf{r}_{\rm{ast}}$ in the
asteroid co-rotating coordinate frame and vectors $\mathbf{r}_{\rm{ecl}}$ in
the ecliptic coordinate frame, and are fixed during the optimization. So
$T_{0}$ is set to be 2451728.000000, corresponding to the date of the first
light curve (2000-07-02), and $\varphi(T_{0})$ is set to be $0^{\circ}$. This
model assumes a constant rotation period and the results of the pole search
are shown in Fig.6. As can be seen, the pole coordinates are well constrained,
thanks to the large range of geocentric ecliptic latitude sampled by the
light-curve dataset.

Figure 6: The $\chi^{2}$ values map for all possible orientations of the pole
direction given in an Aitoﬀ projection of the sky in ecliptic coordinates.
These results are obtained by an inversion method using only light-curve data.
There is one local minima for pole position (lowest $\chi^{2}$ values are
indicated in dark blue), and the corresponding with pole direction is
$(\lambda,\beta)=(231^{\circ},-60^{\circ})$, where $\lambda$ and $\beta$ are
pole ecliptic longitude and latitude, respectively.

A model with a constant rotation period provides the best pole direction
$(\lambda,\beta)=(231^{\circ},-60^{\circ})$ and rotation period $P=7.667714$
h. However, this constant-period model is significantly worse than the YORP
model in fitting the light curves. In the following subsections, these spin
parameters are further optimized by the YORP model as input values.

###  3.3 Shape model and YORP rotational acceleration

Using the YORP model, the shape model and the spin-state parameters of 1999
JD6 are obtained with the whole light-curve data and fit very well, and their
uncertainties are estimated by a Monte Carlo (MC) method. However, there are
no photometric error data for four light curves in 2000. To realistically
estimate uncertainties of light-curve data, a Fourier series of maximum order
determined by an F-test is used to fit the light curves (Magnusson et al.
1996):

$F_{n}=\frac{(\chi^{2}_{0}-\chi^{2}_{n})/(\nu_{0}-\nu_{n})}{\chi^{2}_{n}/\nu_{n}},$
(2)

where $\chi^{2}_{n}$ is the chi-square for the fit of a Fourier series
truncate after order n, with $\nu_{n}=N-(2n+1)$ degrees of freedom, and N is
the number of each light curve point. $\chi^{2}_{0}$ is the chi-square for the
light curve, with $\nu_{0}=N$ degrees of freedom. The root-mean-square
residual is used as the uncertainty of individual light curve points. (85989)
1999 JD6’s light curves are weighted according to their precision.

8000 virtual light-curve data sets are created by randomly adding noise within
its corresponding photometric error range. Each new data set, is repeatedly
inverted and spin parameters are obtained. From the distribution of these
parameters in Fig.7, their uncertainties are estimated.

Figure 7: The uncertainties of the four parameters are calculated by the MC
method. The four parameters are the spin axis direction in ecliptic latitude
$\beta$ and longitude $\lambda$, the sidereal period and the YORP
acceleration. The black dotted line is the Gaussian fitting curve of these
parameter distributions.

The YORP model gives the best-fit values including a pole direction in
ecliptic coordinates $(232\pm 2^{\circ},-59\pm 1^{\circ})$, the sidereal
period $P=7.667749\pm 0.000009\ \rm{h}$ (for JD 2451728.0), and the YORP
rotational acceleration $\upsilon=(2.4\pm 0.3)\times 10^{-8}\ \rm{rad\cdot
d^{-2}}\ (1\sigma\ errors)$. The errors of these parameters are standard
deviation of 1$\sigma$ uncertainties. Fig.8 shows the convex shape model,
which is quite elongated and flat shape. The planar features are the result of
the procedure attempting to match the large amplitude of the light curves. The
agreement between synthetic light curves produced by this shape and real
observations is demonstrated in Fig.9.

Figure 8: Shape model of asteroid (85989) 1999 JD6 shown from equatorial level
(y-axis, left and x-axis, center , 90∘ apart) and pole-on (z-axis, right).The
model’s z-axis is aligned with the rotation axis and axis of maximum inertia.
Figure 9: Example light curves (black dots) of (85989) 1999 JD6 shown with the
synthetic light curves produced by the best YORP model(red dashed curves) and
the best constant-period model(dotted black curves). The geometry of
observation is described with the aspect angle $\theta$, the solar aspect
angle $\theta_{0}$, and the solar phase angle $\alpha$. A full set of light
curves are provied in Fig. 12.

To have an independent estimate, the uncertainty of the $\upsilon$ parameter
is also estimated by varying it around its best value and seeing the increase
in $\chi^{2}$. The same method is also used to compute $\chi^{2}$ for
different fixed values of $\upsilon$ as shown in Fig.10 for $\upsilon$ between
$1.0\times 10^{-8}$ and $3.0\times 10^{-8}\ \rm{rad\cdot d^{-2}}$. These
intervals are larger than those determined by the MC method.

Figure 10: Dependence of the goodness of the fit measured by the $\chi^{2}$ on
the YORP parameter $\upsilon$ for asteroid (85989) 1999 JD6. The dashed curve
is a quadratic fit of the data points. The horizontal dotted line indicates a
4.4% increase in the $\chi^{2}$ , which defines 1$\sigma$ uncertainty interval
of $\pm 0.5\times 10^{-8}\ \rm{rad\cdot d^{-2}}$ given the number of degrees
of freedom.

This paper largely eliminates any possible observational errors by weighting
the light-curve data and gives a YORP rotational acceleration. In upcoming
apparitions, more photometric observations will be required to ensure a robust
detection and decrease the uncertainty of the $\upsilon$ value.

##  4 Discussion

YORP effect is weak in magnitude, and the short-time-scale observation effect
is inconspicuous, making direct detection of YORP effect challenging. Using
light-curve data from a sufficiently extended observation span.

In the article, for (1685) Toro we used archived light curves from 1972 to
2021. From this data set, $\upsilon=3.2\times 10^{-9}\ \rm{rad/d^{2}}$ is
updated in a YORP value. The formal phase shift in $\Delta T\simeq 49$ years
produced by the YORP term is $\upsilon(\Delta T)^{2}/2\simeq 29^{\circ}$. For
this value of $\upsilon$, the $\chi^{2}$ drops by 11% with respect to
$\chi^{2}$ for $\upsilon=0$. The phase offset is $14^{\circ}$ when the light
curves in 1972 is omitted, and the difference between the YORP model and a
constant period model with $\upsilon=0$ is only 3% in $\chi^{2}$. A constant
period model is well fitted the light-curve data, and no YORP signal is
detected, except for the light curves in 1972(Dunlap et al. 1973) and
1988(Hoffmann & Geyer 1990) — only using the light-curve data from 1996 to
2021. This indicates that the YORP model’s direct detection of the YORP effect
is strongly dependent on the observation span.

For (85989) 1999 JD6, the YORP model measured a YORP rotational acceleration
of $(2.4\pm 0.3)\times 10^{-8}\ \rm{rad\cdot d^{-2}}$ only with the optical
light curves. If the four light curves from 2000 are exclude, the observations
span only 6 years, and the difference between the constant period and YORP
model is not significant:the two models fit the data essentially the same, and
the phase offset between the models for $2.4\times 10^{-8}\ \rm{rad\cdot
d^{-2}}$ is only $\sim\ 1^{\circ}$. However, there is also no reasonable
reason to suppose the four light curves from 2000 were all time-shifted in the
same way to mimic the YORP effect. Additional optical and radar data are
required to further confirm this tentative YORP rotational acceleration.

(85989) 1999 JD6 with a tentative YORP detection marks the eleventh direct
detection of YORP to date (all detections are reported in Table 3), and 1999
JD6 is the third contact-binary asteroid with YORP detections (besides Itokawa
and 2001 KZ66).

Table 3: Summary of all detections of YORP to date. Parameters listed in
table: the YORP parameter $\upsilon$, sidereal period (with uncertainty given
in parenthesis), secular rotation-period change($dP/dt$), pole
orientation($\lambda,\beta$), orbital obliquity, the total observation
span(Obs. yrs), and the diameters, whose value was taken from Reddy et al.
(2012) for 1999 JD6 .

Asteroid | $\upsilon$ | Peroid | $dP/dt$ | Pole | Obliquity | Diameters | Obs. yrs | Reference
---|---|---|---|---|---|---|---|---
| [$\times 10^{-8}\ \rm{rad\cdot d^{-2}}$] | [h] | [$\rm{ms\cdot yr^{-1}}$] | [∘] | [∘] | [km] |  |
YORP | $349\pm 30$ | 0.20283333(1) | -1.25 | (180, -85) | 174.3 | 0.113 | 2001-2005 | Lowry et al. (2007)
|  |  |  |  |  |  |  | Taylor et al. (2007)
Eger | $1.1\pm 0.5$ | 5.710156(7) | -3.1 | (226,-70) | 155.6 | 1.5 | 1987-2016 | Ďurech et al. (2018)
Apollo | $5.3\pm 1.3$ | 3.065447(3) | -4.3 | (50,-71) | 161.6 | 1.4 | 1980-2005 | Kaasalainen et al. (2007)
| $5.5\pm 1.2$ | 3.065448(3) | -4.5 | (48,-72) | 162.3 | 1.45 | 1980-2007 | Ďurech et al. (2008)
Cacus | $1.9\pm 0.3$ | 3.755067(2) | -2.3 | (254, -62) | 143.2 | 1.0 | 1978-2016 | Ďurech et al. (2018)
Geographos | $1.15\pm 0.15$ | 5.223336(2) | -2.7 | (58,-49) | 149.9 | 2.56 | 1969-2008 | Ďurec et al. (2008)
Itokawa | $3.54\pm 0.38$ | 12.132371(6) | -45.4 | (128.5,-89.7) | 178.4 | 0.33 | 2000-2007 | Lowry et al. (2014)
Bennu | $264\pm 105$ | 4.2960477(19) | -42.5 | (87,-65) | 161.0 | 0.492 | 1999-2012 | Nolan et al. (2019)
| $363\pm 52$ | 4.296007(2) | -58.4 |  |  |  | 1999-2018 | Hergenrother et al. (2019)
2001 KZ66 | $8.43\pm 0.69$ | 4.985997(42) | -18.3 | (170,-85) | 158.5 | 0.797 | 2010-2019 | Zegmott et al. (2021)
Toro | $0.33\pm 0.03$ | 10.197826(2) | -3.0 | (75,-69) | 160.1 | 3.5 | 1972-2021 | Ďurech et al. (2022)
| $0.32\pm 0.03$ | 10.197827(2) | -2.9 | (60,-72) | 166.7 | - | 1972-2021 | This work
1992 SK | $8.3\pm 0.6$ | 7.320232(10) | -38.8 | (94,-56) | 141.0 | 8.3 | 1999-2020 | Ďurech et al. (2022)
1999 JD6 | $2.4\pm 0.3$ | 7.667749(9) | -12.3 | (232, -59) | 165.3 | 1.6 | 2000-2020 | This work

The contact-binary asteroids are estimated to make up 15% to 35% of the NEA
population (Benner et al. 2006; Jacobson et al. 2016). There are three
possible formation mechanisms of the contact-binary asteroid 1) Two separate
bodies could slowly collided to form a bifurcated shape; 2) The two components
collapsed to form one body due to the BYORP effect shrinking their mutual
orbit (Scheeres 2007; Jacobson & Scheeres 2011); 3) A rubble-pile asteroid
with a weak core undergoes deformation as the rotation rate increases (Sánchez
& Scheeres 2018). For a small core relative to the asteroid size, and a weak-
enough cohesion in the core layer, the deformation is gradually changing with
a dent in the surface emerging first. As the asteroid is rotation increases,
it stretches and transitions to a shape resembling a contact binary before
breaking up into two different sized components. Due to the detection of the
YORP effect of 1999 JD6, the former two evolutionary mechanisms seem to
contradict this scenario. In addition, the overall surface of 1999 JD6 is
uniform with no significant differences in albedo or particle size (Kuroda et
al. 2021). It further reveals that the two components of the (85989) 1999 JD6
came from one body. To sum up, 1999 JD6 might be a rotationally deformed body
currently undergoing the process of rotational fission. The question is
whether the evolutionary mechanisms require a fast rotation rate due to the
YORP effect, whereas 1999 JD6 has a rotation period close to 8 h. Therefore,
it is possible that 1999 JD6’s bifurcated shape was formed during a previous
YORP cycle when it had a faster rotation period. The change in the rotational
state from a previous YORP cycle to the present is possibly caused by small
impacts (Scheeres 2018) and small-scale topographical changes (Statler 2009).
In addition, who the asteroids migrate from obliquity of 0∘(180∘) to obliquity
of 90∘ or from 90∘ to 0∘(180∘) (Golubov & Scheeres 2019), depends on the shape
and rotation state of the asteroids. The current spin-state of 1999 JD6 with a
high obliquity of 165.3∘ is close to an end state of YORP-induced obliquity
shift, but it is more possible that (85989) 1999 JD6 has left the tumbling
regime and is now in the process of migrating towards an obliquity of 90∘.

Under the assumption of zero-conductivity, the shapes of asteroids are roughly
classified into four types ($\rm{I}/\rm{II}/\rm{III}/\rm{IV}$) according to
their response to the YORP rotational moment (Vokrouhlický & Čapek 2002). The
behavior of both the spin and obliquity components of YORP for each type of
asteroid varies with obliquity. Among them, only the type $\rm{IV}$ asteroids
conform to the possible evolution trend of 1999 JD6, that is, the spin and
obliquely components of YORP are both positive for obliquities of $\sim$ 150∘
to 180∘.

It is interesting that for three contact-binary asteroids(Itokawa, 2001 KZ66
and 1999 JD6) for which YORP have been detected, a secular rotation-period
decrease value is larger than other asteroids(while the YORP strength is not
much different). This shows that it is more pronounced to change the rotation
rate of irregular or elongated asteroids. Among the eleven asteroids with YORP
effect, (54509) YORP has the largest YORP intensity and the fastest rotation
period, but its period only decreases by 1.25 $\rm{ms\cdot yr^{-1}}$. It
suggests that changing the rotation rate of already fast-spinning asteroids is
more complicated.

##  5 Conclusion

In this paper, 8 optical light curves are obtained from the monitoring of
(1685) Toro from 11 to 18 February, 2021. With these light-curve data and
published optical light curves during the 1972-2016, a robust but very weak
YORP signal is got, which is verified based on different photometric datasets.
These updated values are as follows:$P=10.197827\pm 0.000002\
\rm{h},(\lambda,\beta)=(60\pm 4^{\circ},-72\pm 2^{\circ})$, $\upsilon=(3.2\pm
0.3)\times 10^{-9}\ \rm{rad\cdot d^{-2}}$ ($1\sigma$ errors).

With published optical light curves, the shape model, spin-state parameters,
and a YORP detection for (85989) 1999 JD6 has been derived. These values are
as follows: a pole direction in ecliptic coordinates $(232\pm 2^{\circ},-59\pm
1^{\circ})$, the sidereal period $P=7.667749\pm 0.000009\ \rm{h}$ (for JD
2451728.0), and the YORP rotational acceleration $\upsilon=(2.4\pm 0.3)\times
10^{-8}\ \rm{rad\cdot d^{-2}}\ (1\sigma\ errors)$. 1999 JD6 is a contact-
binary asteroid. It is likely to form the deformation of a rubble-pile
asteroid with a weak-tensile-strength core due to YORP spin-up and before
that, it might have experienced a fast rotation YORP cycle. We plan to
periodically optical monitor 1999 JD6 in the future, and these additional
optical observations of 1999 JD6 could be used to refine the YORP detection. A
thermophysical analysis is also planned to determine the theoretical YORP
strength, which could lead to determination of the density in-homogeneity for
1999 JD6.

As mentioned above, all these asteroids are with a positive YORP value. While
these asteroids present large light-curve amplitudes and limit the morphology
and observation geometry of asteroids probed, it is dangerous to draw far-
reaching conclusions from a limited sample of only eleven objects. Therefore,
in the future, it is necessary to continue to enlarge the sample of YORP
detection, and begin to periodically monitor photometric observations on YORP
candidates such as 2100 Ra-Shalom (Ďurech et al. 2018). The main targets
include NEAs with a larger period (6$\sim$9 h) and inner main-belt asteroids
(MBAs). It could be crucial for constraining or reacquainting theoretical
concepts of YORP in planetary science.

###### Acknowledgements.

This work has been supported by the B-type Strategic Priority Program of the
Chinese Academy of Sciences (Grant No.XDB41010104), the National Natural
Science Foundation of China (Grant No. 11633009), the Space Debris and Near-
Earth Asteroid Defense Research Project (Grant Nos. KJSP2020020204,
KJSP2020020102), the Civil Aerospace Pre-research Project (Grant Nos. D020304,
D020302) and Minor Planet Foundation. We thank Professor Josef Ďurech for
providing the convexinv_YORP software package and helpful discussion. We also
thank Professor PING Yiding for providing observation data of the 0.8-m
telescope, and thank Dr. LI Xinran and Dr. CHEN Yuanyuan for their helpful
comments to improve the manuscript. This work makes use of the NASA/JPL
HORIZONS ephemeris-generating program. All image reduction and processing are
performed using the Image Reduction and Analysis Facility (IRAF) (Tody 1993).
The work makes use of the Asteroid Lightcurve Photometry Database repository
(Warner et al. 2011b) and the Database of Asteroid Models from Inversion
Techniques (Durech et al. 2010).

## Data Availability

The light curves and python programs utilized in this article are available on
CNEOST at http://157.0.0.68:32280/alc/, where they can be found and downloaded
by searching for the asteroid by name or number. The convex shape model
package (convexinv) is available on DAMIT at
https://astro.troja.mff.cuni.cz/projects/damit/.

## References

  * Benner et al. (2006) Benner, L. A. M., Nolan, M. C., Ostro, S. J., et al. 2006, Icarus, 182, 474
  * Binzel et al. (2001) Binzel, R. P., Harris, A. W., Bus, S. J., & Burbine, T. H. 2001, Icarus, 151, 139
  * Bottke et al. (2006) Bottke, W. F., Vokrouhlický, D., Rubincam, D. P., & Nesvorný, D. 2006, Annual Review of Earth and Planetary Sciences, 34, 157
  * Campins et al. (2009) Campins, H., Kelley, M. S., Fernández, Y., Licandro, J., & Hargrove, K. 2009, Earth Moon and Planets, 105, 159
  * Cheng et al. (2021) Cheng, B., Yu, Y., Asphaug, E., et al. 2021, Nature Astronomy, 5, 134
  * de León et al. (2010) de León, J., Licandro, J., Serra-Ricart, M., Pinilla-Alonso, N., & Campins, H. 2010, A&A, 517, A23
  * Dunlap et al. (1973) Dunlap, J. L., Gehrels, T., & Howes, M. L. 1973, AJ, 78, 491
  * Durech et al. (2010) Durech, J., Sidorin, V., & Kaasalainen, M. 2010, A&A, 513, A46
  * Golubov & Krugly (2012) Golubov, O., & Krugly, Y. N. 2012, ApJ, 752, L11
  * Golubov & Scheeres (2019) Golubov, O., & Scheeres, D. J. 2019, AJ, 157, 105
  * Golubov et al. (2014) Golubov, O., Scheeres, D. J., & Krugly, Y. N. 2014, The Astrophysical Journal, 794, 22
  * Hergenrother et al. (2019) Hergenrother, C. W., Maleszewski, C. K., Nolan, M. C., et al. 2019, Nature Communications, 10, 1291
  * Hoffmann & Geyer (1990) Hoffmann, M., & Geyer, E. H. 1990, Acta Astron., 40, 389
  * Jacobson et al. (2016) Jacobson, S. A., Marzari, F., Rossi, A., & Scheeres, D. J. 2016, Icarus, 277, 381
  * Jacobson & Scheeres (2011) Jacobson, S. A., & Scheeres, D. J. 2011, Icarus, 214, 161
  * Kaasalainen & Torppa (2001) Kaasalainen, M., & Torppa, J. 2001, Icarus, 153, 24
  * Kaasalainen et al. (2001) Kaasalainen, M., Torppa, J., & Muinonen, K. 2001, Icarus, 153, 37
  * Kaasalainen et al. (2007) Kaasalainen, M., Ďurech, J., Warner, B. D., Krugly, Y. N., & Gaftonyuk, N. M. 2007, Nature, 446, 420
  * Kuroda et al. (2021) Kuroda, D., Ishiguro, M., Naito, H., et al. 2021, A&A, 646, A51
  * Lowry et al. (2007) Lowry, S. C., Fitzsimmons, A., Pravec, P., et al. 2007, Science, 316, 272
  * Lowry et al. (2014) Lowry, S. C., Weissman, P. R., Duddy, S. R., et al. 2014, A&A, 562, A48
  * Magnusson et al. (1996) Magnusson, P., Dahlgren, M., Barucci, M. A., et al. 1996, Icarus, 123, 227
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Grav, T., et al. 2014, ApJ, 784, 110
  * Margot et al. (2015) Margot, J. L., Pravec, P., Taylor, P., Carry, B., & Jacobson, S. 2015, Asteroids IV (Univerity of Arizona Press)
  * Marshall et al. (2015) Marshall, S. E., Howell, E. S., Brozović, M., et al. 2015, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 47, AAS/Division for Planetary Sciences Meeting Abstracts #47, 204.09
  * Nolan et al. (2019) Nolan, M. C., Howell, E. S., Scheeres, D. J., et al. 2019, Geophysical Research Letters, 46, 1956
  * Nugent et al. (2016) Nugent, C. R., Mainzer, A., Bauer, J., et al. 2016, AJ, 152, 63
  * Polishook (2014) Polishook, D. 2014, Icarus, 241, 79
  * Polishook & Brosch (2008) Polishook, D., & Brosch, N. 2008, Icarus, 194, 111
  * Pravec et al. (2010) Pravec, P., Vokrouhlický, D., Polishook, D., et al. 2010, Nature, 466, 1085
  * Reddy et al. (2012) Reddy, V., Gaffey, M. J., Abell, P. A., & Hardersen, P. S. 2012, Icarus, 219, 382
  * Reddy et al. (2012) Reddy, V., Gaffey, M. J., Abell, P. A., & Hardersen, P. S. 2012, Icarus, 219, 382
  * Rossi et al. (2009) Rossi, A., Marzari, F., & Scheeres, D. J. 2009, Icarus, 202, 95
  * Rubincam (2000) Rubincam, D. P. 2000, Icarus, 148, 2
  * Sánchez & Scheeres (2018) Sánchez, P., & Scheeres, D. J. 2018, Planet. Space Sci., 157, 39
  * Scheeres (2015) Scheeres, D. 2015, Icarus, 247, 1
  * Scheeres (2018) Scheeres, D. 2018, Icarus, 304, 183, asteroids and Space Debris
  * Scheeres (2007) Scheeres, D. J. 2007, Icarus, 189, 370
  * Statler (2009) Statler, T. S. 2009, Icarus, 202, 502
  * Szabó et al. (2001) Szabó, G. M., Csák, B., Sárneczky, K., & Kiss, L. L. 2001, A&A, 375, 285
  * Taylor et al. (2007) Taylor, P. A., Margot, J.-L., Vokrouhlický, D., et al. 2007, Science, 316, 274
  * Thomas et al. (2011) Thomas, C. A., Trilling, D. E., Emery, J. P., et al. 2011, AJ, 142, 85
  * Tichy et al. (1999) Tichy, M., Moravec, Z., Lopez, A., et al. 1999, Minor Planet Electronic Circulars, 1999-J40
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Vokrouhlický et al. (2011) Vokrouhlický, D., Ďurech, J., Polishook, D., et al. 2011, AJ, 142, 159
  * Vokrouhlický & Čapek (2002) Vokrouhlický, D., & Čapek, D. 2002, Icarus, 159, 449
  * Warner (2014) Warner, B. D. 2014, Minor Planet Bulletin, 41, 213
  * Warner (2015) Warner, B. D. 2015, Minor Planet Bulletin, 42, 256
  * Warner (2018) Warner, B. D. 2018, Minor Planet Bulletin, 45, 366
  * Warner et al. (2009) Warner, B. D., Harris, A. W., & Pravec, P. 2009, Icarus, 202, 134
  * Warner & Stephens (2019) Warner, B. D., & Stephens, R. D. 2019, Minor Planet Bulletin, 46, 423
  * Warner & Stephens (2020) Warner, B. D., & Stephens, R. D. 2020, Minor Planet Bulletin, 47, 290
  * Warner et al. (2011a) Warner, B. D., Stephens, R. D., & Harris, A. W. 2011a, Minor Planet Bulletin, 38, 172
  * Warner et al. (2011b) Warner, B. D., Stephens, R. D., & Harris, A. W. 2011b, Minor Planet Bulletin, 38, 172
  * Zegmott et al. (2021) Zegmott, T. J., Lowry, S. C., Rożek, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 4914
  * Ďurec et al. (2008) Ďurec, J., Vokrouhlický, D., Kaasalainen, M., et al. 2008, A&A, 489, L25
  * Ďurech et al. (2008) Ďurech, J., Vokrouhlický, D., Kaasalainen, M., et al. 2008, A&A, 488, 345
  * Ďurech et al. (2012) Ďurech, J., Vokrouhlický, D., Baransky, A. R., et al. 2012, A&A, 547, A10
  * Ďurech et al. (2018) Ďurech, J., Vokrouhlický, D., Pravec, P., et al. 2018, A&A, 609, A86
  * Ďurech et al. (2022) Ďurech, J., Vokrouhlický, D., Pravec, P., et al. 2022, A&A, 657, A5

##  Appendix A additional figures

Figure 11: All light curves generated with the constant-period of convex-
inversion shape model (dotted black curves) and the YORP model (red dashed
curves) of asteroid (1685) Toro plotted over all available light-curve data
(black dots).

Figure 11: -continued.

Figure 11: -continued.

Figure 11:  -continued.

Figure 12:  All light curves generated with the constant-period of convex-
inversion shape model (red dots) and the YORP model (blue dots) of asteroid
(85989) 1999 JD6 plotted over all available light-curve data (black dots).
Light-curve details can be found in Table 2.


# Multimodal dynamics modeling for off-road autonomous vehicles

Jean-François Tremblay1, Travis Manderson1, Aurélio Noca1, Gregory Dudek1 and
David Meger1 1All authors are with the Mobile Robotics Laboratory (MRL),
McGill University, Montréal, Canada, H3A 2A7. Corresponding e-mail:
<EMAIL_ADDRESS>This work was supported by the Fonds de Recherche du Québec,
Nature et Technologie (FRQNT), Hydro-Québec and the NSERC Canadian Robotics
Network (NCRN).

###### Abstract

Dynamics modeling in outdoor and unstructured environments is difficult
because different elements in the environment interact with the robot in ways
that can be hard to predict. Leveraging multiple sensors to perceive maximal
information about the robot’s environment is thus crucial when building a
model to perform predictions about the robot’s dynamics with the goal of doing
motion planning. We design a model capable of long-horizon motion predictions,
leveraging vision, lidar and proprioception, which is robust to arbitrarily
missing modalities at test time. We demonstrate in simulation that our model
is able to leverage vision to predict traction changes. We then test our model
using a real-world challenging dataset of a robot navigating through a forest,
performing predictions in trajectories unseen during training. We try
different modality combinations at test time and show that, while our model
performs best when all modalities are present, it is still able to perform
better than the baseline even when receiving only raw vision input and no
proprioception, as well as when only receiving proprioception. Overall, our
study demonstrates the importance of leveraging multiple sensors when doing
dynamics modeling in outdoor conditions.

EKF

    extended Kalman filter
IMU

    inertial measurement unit
ANN

    artificial neural network
CNN

    convolutional neural network
RNN

    recurrent neural network
GRU

    gated recurrent unit
SLAM

    simultaneous localization and mapping
MBRL

    model-based reinforcement learning
RL

    reinforcement learning
RSSM

    recurrent state-space model
MVAE

    multimodal variational autoencoder
VAE

    variational autoencoder
RSSM

    recurrent state-space model
SLAC

    stochastic latent actor-critic
ICP

    iterative closest point
ELBO

    evidence lower bound
SAC

    soft actor-critic
RMSE

    root mean square error
MRSSM

    multimodal recurrent state-space model

##  I Introduction

A wheeled robot moving in the forest interacts with rocks, trees, ice, mud,
and a multitude of other elements, making its motions complex and difficult to
predict – a fact known to generations of humans who have previously been the
pilots of off-road vehicles. The prediction of such future motion is key to
control autonomous off-road vehicles. Understanding likely physical outcomes
can begin with knowing the proprioceptive and command information, which are
highly predictive on human-made roads, but exteroception is critical to
predict motions and to reason about roughness, ground traction or slipping,
and contact with foliage. Furthermore, each sensor in this environment should
be treated as unreliable, as they can fail from a variety of causes, from low
light to being covered with mud and soil.

This paper proposes a novel dynamics modeling approach suited to motion
prediction in the forest. We extend a recent multimodal, latent space approach
for dynamics modeling that is robust to missing modalities. Our method scales
well with the number of sensors present and is able to model the vehicle’s
dynamics from different sensor combinations, taking advantage of multiple
exteroceptive sensors when available. These predictions are propagated through
time building upon a recurrent neural network (RNN) based model [1] enabling
prediction over long horizons.

Our latent embedding allows us to build an abstract representation of the
state of the robot in a self-supervised manner, which automatically
incorporates information such as terrain classes and obstacles such as tree
branches. We learn without any manual terrain labeling, which would be
challenging in our scenario: while it is possible to label the few classes
present in urban driving (road, sidewalk, grass, etc) there are much wider
range of less clearly demarcated terrain elements present in the forest (Is
the robot on a rock? Which wheel is on the rock? Is the rock wet?) Rather, our
method uses _post hoc_ state-estimation from on-board vehicle lidar [2] as a
direct training signal, enabling future prediction for several seconds forward
in time, in new trajectories in the same forest region.

Multimodality enables our model to leverage all sensors present on the robot,
giving it access to maximal information about the state of the robot and its
surrounding environment. However, sensors can be prone to failure, especially
in natural and remote environments where replacements or repair are not
immediate Thus, we design our dynamics model to be able to perform predictions
without retraining and with minimal loss of accuracy in the case of missing
modalities. For example, a forest robot can benefit from the combination of
vision, lidar and inertial measurement units during the daytime and seamlessly
fall back to lidar and IMUs only during the nighttime with limited impact on
the performance. Additionally, our approach can accommodate sensors with
different rates where one or more sensors are effectively unavailable at a
given timestep.

We demonstrate the accuracy of our method, in terms of the final pose
predictions, on our real-world forest dataset, where we have 17 trajectories
of a robot traversing forest terrain. Our results show that our model is able
to benefit from all modalities when they are present, performing better than
our baseline, while still handling missing data when only receiving raw vision
or proprioception.

Overall, our contributions include:

  * •

The derivation of our multimodal recurrent state-space model (MRSSM), an
approach that combines multimodal learning with latent time-series prediction.

  * •

The derivation of a new multimodal variational autoencoder (MVAE) evidence
lower bound (ELBO) that maximizes the likelihood of held-out modalities and
takes into account all available training data.

  * •

The first use, to our knowledge, of MVAEs for challenging outdoor robotics
data.

  * •

Demonstration of missing-modality prediction robustness, again within this
challenging domain.

##  II Related work

We first review our application domain, autonomous navigation in outdoor
environments, and then move to the related algorithms for model-based
reinforcement learning and multimodal learning.

#### Navigation in difficult outdoor environments

Autonomous robots navigating in unstructured, outdoor environments face
adverse conditions in terms of complex, uneven terrain that is often slippery
and contains rubble that causes slippage. Additionally, vegetation often
intermittently occludes sensors such as lidar and cameras. A common strategy
for learning where to navigate in these environments is to use behavior
cloning [3] where a user provides demonstrations of traversing natural paths
[4, 5]. Similarly, end-to-end techniques have been used to effectively learn a
navigation policy by providing expert steering demonstrations for autonomous
driving [6] as well as difficult underwater environments [7, 8]. Another
common approach is to perform appearance-based segmentation, and train a
classifier to label the segmentation classes with a score corresponding to the
navigability difficulty [9, 10, 11, 12]. Often, learning the navigability
score is difficult, so some authors [13, 14] have used expert demonstrations
to learn the traversable classes.

####  MBRL

is an approach to reinforcement learning (RL) where we learn to predict the
results of our action on the state of the robot and the reward through
(self-)supervised learning. This is usually viewed as a regression problem,
which can be solved though Gaussian processes, for example [15], [16]. Neural
networks can also be used [17]. Until recently, these models did not work on
high-dimensional observations such as vision. Visual foresight [18, 19, 20,
21] was developed in a series of papers and used a video prediction model to
learn visual dynamics.

#### Latent space models

are one category of approaches for model-based reinforcement learning (MBRL).
They rely on mapping observations to a latent space, usually of lower
dimensionality than the observation space, and perform prediction in that
space. This dimensionality reduction helps sample efficiency and is learned,
usually with reconstruction losses using some encoder-decoder combination [1,
22, 23]. This abstract latent space, in our case, is useful for fusing
different observation modalities.

#### Multimodal learning

There are multiple facets contributing to the utility of multimodal learning
algorithms, the most important being: their scalability to a high number of
modalities, their robustness to missing modalities as well as their
performance on the task at hand achieved by taking advantage of all of the
modalities present. [24] performed the fusion of drone and ground-robot
imagery, in the context of robot dynamics modeling, by concatenating the
output of two convolutional encoders and feeding the result into a RNN. [25]
proposed Gated Recurrent Fusion Units to perform weighted fusion of multimodal
data in the context of self-driving vehicles. Their approach allows the model
to give more importance to different modalities through weighting, while our
approach does the same through uncertainty-aware fusion using a product of
experts. However, their approach does not allow for missing modalities. Using
a different philosophy altogether, [26] described an approach to perform
sensor selection for a reactive steering policy. Their approach did not
perform fusion but rather dynamically selected one sensor from four cameras
and one lidar to perform the prediction. However, we argue that leveraging all
available information rather than choosing one particular sensor is a better
choice in most cases. Closer to [27], detailed in subsection III-B, [28]
develop another approach to multimodal variational autoencoder, effectively
setting all missing modalities values to -2 randomly at training time. While
this approach is capable of representing data with many modalities, the
approach of [27] has shown the ideal combination of flexibility to missing
modalities and prediction accuracy, which we argue is crucial to deploy
multimodal models on field robots.

##  III Methods

We seek a predictive dynamics model capable of capturing the motions of a
robot over several seconds through environments with complex wheel-terrain-
obstacle interactions. High-dimensional sensing with images and/or lidar
allows interpretation of the environment, and we therefore select generative
time-series models that can interpret and predict physical robot states as
well as images. In particular, we turn to the so-called recurrent state-space
model (RSSM) introduced in [1], for its ability to accurately model
transitions in high-dimensional sensory spaces. A RSSM is a time-series latent
variable model, where observations $o_{t}$ at time $t$, do not directly
contain all relevant information, such as wheel contact and environment
structure. Therefore, our method will attempt to recover a latent state
representation $s_{t}$ sufficient to explain observations and, along with
actions $a_{t}$, transition behaviors.

To capture the multimodal sensor inputs required to sense a dense forest, we
utilize the multimodal variational autoencoder (MVAE) model [27]. This
approach scales well to several high-dimensional input modalities and
gracefully handles prediction when modalities are missing, such as will occur
each night for a robot using imagery from passive illumination. In our method,
the multimodal recurrent state-space model (MRSSM) combines these two recent
advances, which are compatible through a modification of the reconstruction
likelihood term in RSSMs training objective. We further extend the training-
time sampling procedure of multimodal variational autoencoder (MVAE), which
provides its strong performance with missing modalities.

###  III-A Recurrent state-space model

Our work builds on the RSSM introduced in PlaNet [1], which combines a latent
representation from a variational autoencoder (VAE) model with an RNN to
capture dynamics transitions. While several previous approaches combine a VAE
and an RNN (e.g., World Models [23] and [29]), PlaNet has shown an increase in
accuracy of prediction over longer prediction horizons by incorporating end-
to-end training.

The RSSM consists of the following components:

  1. 1.

a transition model: $p(s_{t}|s_{t-1},a_{t-1})$,

  2. 2.

an observation model: $p(o_{t}|s_{t})$,

  3. 3.

the variational state distribution: $q(s_{t}|s_{t-1},a_{t-1},o_{t})$.

The transition model is implemented as an hybrid deterministic/stochastic RNN,
as described by [1]. The observation model used to reconstruct images in the
previous work was a transposed convolutional neural network (CNN). We will
describe our extension to multimodal observations in Section III-C. Training
of each of the required models is performed in a unified procedure based on
variational inference. Specifically, the desired parameters are those that
maximize the likelihood of the observations gathered during training, under
the time-series model. While this objective is intractable, a lower bound can
be stated to guide optimization, as presented by [1]:

$\displaystyle\ln
p(o_{1:T}|a_{1:T-1})=\ln\int\prod_{t=1}^{T}p(s_{t}|s_{t-1},a_{t-1})p(o_{t}|s_{t})ds_{1:T}$
(1)
$\displaystyle\geq\sum_{t=1}^{T}\Big{(}\operatorname*{\mathbb{E}}_{q(s_{t}|o_{\leq
t},a_{<t})}\left[\ln p(o_{t}|s_{t})\right]$ $\displaystyle\hskip
10.0pt-\operatorname*{\mathbb{E}}_{\hskip 80.0ptq(s_{t-1}|o_{\leq
t-1},a_{<t-1})}\left[\mathcal{KL}\left[q(s_{t}|o_{\leq
t},a_{<t})||p(s_{t}|s_{t-1},a_{t-1})\right]\right]\Big{)}$ (2)
$\displaystyle:=\sum_{t=1}^{T}\text{ELBO}(o_{t}),$ (3)

where $s_{0}$ is sampled from our initial state distribution, a Gaussian with
mean 0 and variance 1. From observed sequences of observations and actions,
the parameters of each component model (neural network weights) can be
optimized, leading to effective predictions of the future at test time.

###  III-B Multimodal variational autoencoder

Now, we need some way to map observations of different modalities to a common
latent space. This is non-trivial, especially when modalities can be missing
or are not amenable to the same model architectures. Instead of having single
observation $o_{t}$ as in subsection III-A, we now have
$O_{t}=\\{o_{t}^{1},\dots,o_{t}^{N}\\}$, a set containing $N$ sensor readings
from the modalities present at time $t$. The modalities present can change,
for example due to sensor failure. Methods which concatenate the output of two
encoders [24] are inappropriate for this scenario.

MVAEs [27] provide a way to encode multimodal data, while being robust to any
combination of missing modalities. Each modality has one encoder and decoder,
which can have arbitrary architectures. Each encoder of the modalities present
at time $t$ separately outputs an estimate of the mean and variance for the
latent variable $s_{t}$. These estimates are then fused with what is called a
product of experts [30]:

$p(s_{t}|O_{t},s_{t-1},a_{t-1})=\frac{1}{Z}p(s_{t}|s_{t-1},a_{t-1})\prod_{o^{i}_{t}\in
O_{t}}q(s_{t}|o^{i}_{t})$ (4)

where $p(s_{t}|s_{t-1},a_{t-1})$ is a prior coming from our RSSM transition
model, called the prior expert, and $q(s_{t}|o^{i}_{t})$ is the distribution
outputted by the encoder for the modality $i$. This assumes the independence
of the different modalities conditional on the state. $Z$ is the partition
function, or normalization constant, which is intractable in a general product
of expert. However, when our variational distribution $q$ and our prior is
Gaussian, which is our case, we can compute the product’s mean and variance in
differentiable closed form [27].

Because the observations of different modalities are assumed to be
conditionally independent given the states, we have that:
$p(O_{t}|s_{t})=\prod_{o_{t}^{i}\in O_{t}}p(o_{t}^{i}|s_{t})$ and $\ln
p(O_{t}|s_{t})=\sum_{o_{t}^{i}\in O_{t}}\ln p(o_{t}^{i}|s_{t})$. In practice,
$\sum_{o_{t}^{i}\in O_{t}}\lambda_{i}\ln p(o_{t}^{i}|s_{t})$ is used in the
ELBO of the MVAE to weight the reconstruction error of modalities differently
[27]; this amounts to having the factorization
$p(O_{t}|s_{t})=\prod_{o_{t}^{i}\in O_{t}}p(o_{t}^{i}|s_{t})^{\lambda_{i}}/W$
where $W$ is a partition function which can be ignored when maximizing the
ELBO. Having low $\lambda_{i}$ leads to having a flat likelihood for modality
$i$, assigning higher likelihood to reconstructions with high reconstruction
error compared to high $\lambda_{i}$.

###  III-C Our approach: multimodal recurrent state space models

Figure 1: Diagram of how our model can do state estimation and perform
predictions. In the filtering step, the model takes the previous posterior and
action from time $t-1$ to generate a prior for time $t$. This prior is then
fused with experts for the modalities present at time $t$ to generate a
posterior. In the prediction step, one can take the posterior and action at
time $t$ and run it through the transition model to get the prior for time
$t+1$. To predict the future observations for time $t+1$, one can then sample
from this prior and feed it to the observation model of the modalities one
wants to predict. This prior can also be given to the transition model, along
with a planned action, to perform predictions farther in the future,

As our goal is to predict outcomes for several steps from multimodal
observations, we seek to combine RSSM and MVAE into a coherent approach.
Noting that PlaNet already utilizes a VAE to represent the observation model
of a single mode, $p(o_{t}|s_{t})$, the key for combining the methods is to
properly integrate the multimodal observation likelihood into the RSSM
variational lower bound training objective. The derivation is identical to
[1], just by substituting the observation likelihood $p(o_{t}|s_{t})$ for the
MVAE’s observation likelihood of $p(O_{t}|s_{t})$ Equation 4.

As suggested by [27], one needs to randomly remove modalities at training time
to be able to train the individual experts and be robust to missing
modalities. An extension of what [27] suggest to time series is randomly
selecting a subset of present modalities, $\tilde{I}$, and creating a new time
series, $\tilde{O}_{t}=\\{o^{i}_{t}|i\in\tilde{I},o^{i}_{t}\in O_{t}\\}$ and
maximizing ELBO($\tilde{O}_{t}$). However, this does not take advantage of the
fact that we have access to all information present at time $t$, $O_{t}$,
because the modalities were randomly sub-sampled as part of the training
procedure rather than being truly missing. It is interesting to train the
model to be able to predict modalities that are not present in
$\tilde{O}_{t}$. For example, in our dynamics modeling application, it could
be beneficial to train the vision expert to be able to predict velocities and
accelerations from vision data alone, something that is not explicitly learnt
by maximizing $\text{ELBO}(\tilde{O}_{t})$. We can thus reformulate the ELBO
for $O$ instead of $\tilde{O}$:

$\displaystyle\text{ELBO}(O_{t})=\operatorname*{\mathbb{E}}_{q(s_{t}|\tilde{O}_{\leq
t},a_{<t})}\left[\gamma_{t}\sum_{o^{i}_{t}\in O_{t}}\lambda_{i}\ln
p(o^{i}_{t}|s_{t})\right]$ $\displaystyle\hskip
10.0pt-\beta\operatorname*{\mathbb{E}}_{\hskip 80.0ptq(s_{t-1}|\tilde{O}_{\leq
t-1},a_{<t-1})}\left[\gamma_{t-1}\mathcal{KL}\left[q(s_{t}|O_{\leq
t},a_{<t})||p(s_{t}|s_{t-1},a_{t-1})\right]\right]$ (5) $\displaystyle\hskip
43.0pt:=\text{ELBO}(O_{t}|\tilde{O}_{t})$ (6)

where

$\gamma_{t}=\frac{q(s_{t}|O_{\leq t},a_{<t})}{q(s_{t}|\tilde{O}_{\leq
t},a_{<t})}$ (7)

This shows that by doing importance sampling, we can express the ELBO for
$O_{t}$ using the product of experts for $\tilde{O}_{t}$. We can get the
posterior $q(s_{t}|O_{\leq t},a_{<t})$ by always choosing all experts as our
first subset, and then using it when computing the $\mathcal{KL}$-divergence
term of the other sub-sampled $\tilde{O}_{t}$. It makes sense, in an intuitive
way, to use this posterior instead of the one coming from $\tilde{O}_{t}$.
Indeed, this $\mathcal{KL}$-divergence term is what enables our model to
predict future states; by using the posterior from $O_{t}$, we are training
our prior to be as close to the posterior containing the most information, no
matter the information that was used to generate this prior. Note that this
works for computing the ELBO of any other subset of $O_{t}$, but we restrict
our attention to $O_{t}$ only. When framed this way, training different
combinations of experts then simply amounts to choosing different importance
sampling distributions. In practice, due to the high dimensionality of our
latent state-space, we do not compute the importance sampling ratio for both
expectations in the ELBO and keep them to one, giving us a biased estimator of
the ELBO. The approach, which we call a multimodal recurrent state-space model
(MRSSM), is illustrated in Figure 1.

##  IV Experiments

To verify that our model uses vision in a meaningful way, we first conduct
controlled experiment in simulation. We then move to more extensive real-world
evaluation where we compare our new proposed ELBO of subsection III-C to [27]
and to two other baselines. The experiments consists of building a model
capable of predicting future velocities up to a chosen horizon $H$, given past
observations $O_{\leq t}$ and intended actions $a_{<t+H}$. Our code is based
on an existing PyTorch implementation of
PlaNet111https://github.com/Kaixhin/PlaNet and
MVAE222https://github.com/mhw32/multimodal-vae-public although it has been
through substantial changes.

###  IV-A Simulation experiments

Figure 2: Sample image from our off-road robot simulator, where we can see the
transition between terrain types, the gray rocky terrain having maximal
traction and the white icy terrain having very low traction. Furthermore,
small rocks visible here will affect the robot’s dynamics, and bigger trees
and rock add non-traversable obstacles.

We first validate our method using an off-road wheeled robot simulator based
on Unreal Engine [31], shown in Figure 2, which gives rich visuals and
realistic dynamics. We test in simulation because data is hard to get and we
can have full control of the terrains, as well as access to ground-truth. We
generated 4.38 hours of data in one map using a random exploration policy,
where the throttle and steering were sampled from a normal distribution
centered at half throttle and centered steering respectively. The dataset is
split into 684 trajectories, 10% of which we keep as a validation dataset. For
simplicity, only three modalities were recorded: linear and angular velocity
for proprioception, and a front-facing camera for exteroception. The map has
three terrain types with different friction coefficients, enabling us to
simulate situations where traction can change rapidly. Furthermore, there are
trees and rocks acting as rigid obstacles. The simulation enables us to build
a dataset where we know exactly when the robot’s traction changes. By using
this information only for evaluation purposes, we can see the contribution of
vision, where we expect it to help predict changes in traction, but not
provide much benefit on flat terrain with constant traction. The error metric
we use is the Euclidean norm of the translation error of the predicted final
pose, obtained by integrating the linear and angular velocity outputted by our
model over the prediction horizon. This metric allows us to evaluate both the
angular and linear velocity predictions and how their errors accumulate with
time. The results of this basic experiment aiming to assess the influence of
different modalities in varied scenario is visible in Figure 3, where we
isolate sequences where the robot transitions between different terrains. Our
results show that using all modalities results in a 40% decrease in
translation error during terrain transitions, and even during non-transitions,
our model using all modalities performs better than only proprioception
sensing. We demonstrated that our robot is capable of predicting traction
changes from visual information in an unsupervised manner with a controlled
simulation experiment, now we quantify the performance of our method using a
real-world dataset.

Figure 3: Median final pose error for 1s predictions, when the robot cross
different terrain types and for the whole test dataset. For scale, the robot
was going at an average speed of 5 m/s. We can see that when using all
modalities including vision, performance degrades minimally when the traction
situation changes. However, that is not the case when using only
proprioception has the model has no access to information on the changing
terrain. This shows that overall our model performs well when using only
proprioception, but it is able to benefit from vision to predict traction
changes.

###  IV-B Forest robot dynamics modeling

(a) Branches

(b) Downhill

Figure 4: Different events that affect the dynamics of the robot in our
dataset. Branches can cause slippage or block the robot. Going up or downhill
will also affect the dynamics through gravity. These images illustrate the
diversity of scenario our method is being tested on, and the real-world
conditions of the dataset used.

Here we are working with the Montmorency dataset [32], consisting of 17
trajectories with high quality simultaneous localization and mapping (SLAM)
based on off-line iterative closest point (ICP), allowing us to know the
velocity of the robot. These velocities are affected by the environment and
contain four sites with different terrain types, causing varied effects on the
dynamics. This includes going over branches, loose moss and steep hill
ascents/descents illustrated in Figure 4. This dataset consists of visual,
lidar and inertial data for a skid-steered robot with a maximal speed of 1m/s,
along with the target forward and turning velocity given by a human driven
operating the robot with a remote controller. Thus, the five observation
modalities with each their own expert are: linear velocity, angular velocity,
tri-axial accelerometer, lidar height map of a 4 m2 area around the robot and
the front camera view. Both the lidar height map and front camera are given to
our model with a $64\times 64$ resolution, with the three channels for the
lidar height map being the average $z$ coordinate of points falling in that
pixel, the standard deviation of the height of these points and the average
intensity of those points. We used a fully connected encoder and observation
model for the low-dimensional proprioceptive modalities, and a convolutional
encoder and observation model for the vision and lidar. The angular velocity
was represented as a quaternion for the relative rotation in the last 0.1s
(using the reconstruction loss proposed by [33]), which represents the
discrete time-step used in our model at which all modalities are sampled. One
baseline we look at is assuming that the robot is always at its target
velocity, called the Control baseline. For completeness, another baseline uses
an RSSM where the output of multiple encoders are concatenated; while this
approach allows the fusion of multiple modalities, it is not compatible with
missing modalities at test time.

We use the same pose prediction error metric as in subsection IV-A. In Figure
5, we can see that our model beats the Control baseline, even when doing
predictions up to three seconds into the future. To demonstrate the robustness
to missing modalities of our approach, we remove different combination of
modalities at test time, for the entire duration of the trajectory. We test
when all modalities are present, when only vision is present, and when only
proprioception is present (accelerometer and velocity measurements).
Furthermore, we test both the ELBO from [27], and the new proposed ELBO from
subsection III-C. We focus our evaluation on the site in a Maple grove, as
this is the only site where we have more than 30 minutes of driving, with a
total of 126 minutes in eight trajectories at this site. Validating on
trajectories on sites where there was only a few minutes of data present in
the training set did not lead to conclusive results. We perform validation in
a leave-one-out manner, reporting root mean square error (RMSE) statistics on
the eight validation trajectories.

Interestingly, when using the new ELBO of subsection III-C we could remove all
modalities, including velocity measurements, and only keep vision. Even in
this scenario, our model is able to keep track of the robot’s velocity (doing
a form of deep visual odometry) and perform predictions that were better than
the two baselines. This was not the case for the previous ELBO obtained from
[27], where modality ablation lead to a much larger increase, doubling the
median error. Overall, for three seconds ahead pose predictions, our model
leveraging all modalities with the ELBO of subsection III-C has the best
average median error of 11.3cm, compared to 23.6cm for the Control baseline.

Perhaps unsurprisingly, the RSSM which concatenates the output of multiple
encoders performed similarly to our approach using the ELBO from [27] when the
latter has access to all modalities. However, both of these approaches were
bested by our method using the new proposed ELBO, showing that there is a gain
in predictive power when learning to reconstruct missing modalities. There was
no modality ablation study for the concatenation baseline as it is unable to
perform predictions when modalities are missing. Altogether, our results
demonstrate that the combination of the proposed MRSSM and the ELBO of
subsection III-C allows for 1) lower RMSE of final pose predictions 2)
robustness to missing modalities at test time. Examples of predicted linear
velocity and images for our approach with the new ELBO are visible in Figure
6.

Figure 5: Error distribution for our baselines (left of dotted line) and our
proposed MRSSM under various modality ablations (right of dotted line) for the
Montmorency dataset, with a prediction horizon of 3s. In each distribution,
the first quartile, the median and the third quartile are drawn as dotted
lines, and the median is displayed on top. One can see that the distribution
of error, both in terms of median but also in terms of the tail of large
errors, is lower for the MRSSM with the new proposed ELBO. Furthermore, the
new proposed ELBO performs better when modalities are removed.

(a) Turning maneuver  (b) Uneven terrain  (c) Acceleration from stopped state

Figure 6: Sample of predicted linear velocities and images, over an horizon of
1s, or 10 timesteps of 0.1s. Note how the visual prediction on the right,
shown one timestep ahead of the current time, can become blurry for sequences
that are more difficult to predict. One can also see the visual diversity even
in the same testing site: (a) and (b) are in late fall with no leaves left on
trees while (c) is in mid-October with colored leaves and a low sun vastly
changing the visual appearance of the scene. This highlights the difficulty of
learning from visual data in outdoor scenarios.

##  V Conclusion

Driving off-road, across the wide spectrum of natural terrain requires
integrating all available information. In this paper, we have demonstrated
that modern approaches for time-series modeling and learning multimodal
representations are suitable for the problem of motion prediction in forests.
We derive a novel variational bound for learning time-series models using
multimodal data that requires no human labeling, but still captures important
environmental effects such as changing traction caused by different types. Our
quantitative results show that the motions of a four wheeled robot can be
predicted accurately as it traverses a wide variety of slopes and contacts
branches. Our future work will include closing the control loop around our
model, both by performing active exploration to collect informative data and
by implementing planning and control using the valuable predictive power of
our technique.

## References

  * [1] Danijar Hafner et al.  “Learning Latent Dynamics for Planning from Pixels”  In _ICML_ 97, Proceedings of Machine Learning Research  Long Beach, California, USA: PMLR, 2019, pp. 2555–2565
  * [2] F. Pomerleau et al.  “Long-term 3D map maintenance in dynamic environments”  In _ICRA_ , 2014, pp. 3712–3719
  * [3] Dean A Pomerleau  “Alvinn: An autonomous land vehicle in a neural network”  In _NeurIPS_ , 1989, pp. 305–313
  * [4] A. Giusti et al.  “A Machine Learning Approach to Visual Perception of Forest Trails for Mobile Robots”  In _IEEE RA-L_ 1.2, 2016, pp. 661–667
  * [5] N. Smolyanskiy, A. Kamenev, J. Smith and S. Birchfield  “Toward low-flying autonomous MAV trail navigation using deep neural networks for environmental awareness”  In _IROS_ , 2017, pp. 4241–4247
  * [6] Mariusz Bojarski et al.  “End to End Learning for Self-Driving Cars”  In _CoRR_ , 2016
  * [7] Travis Manderson, Juan Camilo Gamboa Higuera, Ran Cheng and Gregory Dudek  “Vision-based autonomous underwater swimming in dense coral for combined collision avoidance and target selection”  In _IROS_ , 2018, pp. 1885–1891  IEEE
  * [8] Travis Manderson et al.  “Vision-Based Goal-Conditioned Policies for Underwater Navigation in the Presence of Obstacles”  In _RSS_ , 2020
  * [9] Boris Sofman et al.  “Improving robot navigation through self-supervised online learning”  In _JFR_ 23.11‐12, 2006, pp. 1059–1075
  * [10] Brandon Rothrock et al.  “SPOC: Deep Learning-based Terrain Classification for Mars Rover Missions”  In _AIAA SPACE_ , 2016
  * [11] Luqi Wang et al.  “A Collaborative Aerial-Ground Robotic System for Fast Exploration”, 2018  arXiv:1806.02487 [cs.RO]
  * [12] Jeffrey Delmerico, Elias Mueggler, Julia Nitsch and Davide Scaramuzza  “Active autonomous aerial exploration for ground robot path planning”  In _IEEE RA-L_ 2.2, 2017, pp. 664–671
  * [13] Dan Barnes, Will Maddern and Ingmar Posner  “Find your own way: Weakly-supervised segmentation of path proposals for urban autonomy”  In _ICRA_ , 2017, pp. 203–210  arXiv:1610.01238v3
  * [14] L. Tang et al.  “From one to many: Unsupervised traversable area segmentation in off-road environment”  In _ROBIO_ , 2017, pp. 787–792
  * [15] Jonathan Ko and Dieter Fox  “GP-BayesFilters: Bayesian filtering using Gaussian process prediction and observation models”  In _Autonomous Robots_ 27.1 Springer, 2009, pp. 75–90
  * [16] M. P. Deisenroth, D. Fox and C. E. Rasmussen  “Gaussian Processes for Data-Efficient Learning in Robotics and Control”  In _IEEE T-PAMI_ 37.2, 2015, pp. 408–423
  * [17] J. C. Gamboa Higuera, D. Meger and G. Dudek  “Synthesizing Neural Network Controllers with Probabilistic Model-Based Reinforcement Learning”  In _IROS_ , 2018, pp. 2538–2544
  * [18] Frederik Ebert et al.  “Visual foresight: Model-based deep reinforcement learning for vision-based robotic control”  In _arXiv preprint arXiv:1812.00568_ , 2018
  * [19] C. Finn and S. Levine  “Deep visual foresight for planning robot motion”  In _ICRA_ , 2017, pp. 2786–2793
  * [20] Suraj Nair and Chelsea Finn  “Hierarchical Foresight: Self-Supervised Learning of Long-Horizon Tasks via Visual Subgoal Generation”  In _arXiv preprint arXiv:1909.05829_ , 2019  arXiv:1909.05829 [cs.LG]
  * [21] Chelsea Finn, Ian Goodfellow and Sergey Levine  “Unsupervised Learning for Physical Interaction through Video Prediction”  In _NeurIPS_ Curran Associates, Inc., 2016, pp. 64–72
  * [22] Carles Gelada et al.  “DeepMDP: Learning Continuous Latent Space Models for Representation Learning”  In _ICML_ 97, Proceedings of Machine Learning Research  Long Beach, California, USA: PMLR, 2019, pp. 2170–2179
  * [23] David Ha and Jürgen Schmidhuber  “World Models”  In _CoRR_ abs/1803.10122, 2018  arXiv:1803.10122
  * [24] Travis Manderson, Stefan Wapnick, David Meger and Gregory Dudek  “Learning to Drive Off Road on Smooth Terrain in Unstructured Environments Using an On-Board Camera and Sparse Aerial Images”  In _ICRA_ , 2020
  * [25] A. Narayanan, A. Siravuru and B. Dariush  “Gated Recurrent Fusion to Learn Driving Behavior from Temporal Multimodal Data”  In _IEEE RA-L_ 5.2, 2020, pp. 1287–1294
  * [26] Shihong Fang and Anna Choromanska  “Multi-modal Experts Network for Autonomous Driving”  In _ICRA_ , 2020
  * [27] Mike Wu and Noah Goodman  “Multimodal Generative Models for Scalable Weakly-Supervised Learning”  In _NeurIPS_ Curran Associates, Inc., 2018, pp. 5575–5585
  * [28] Martina Zambelli, Antoine Cully and Yiannis Demiris  “Multimodal representation models for prediction and control from partial information”  In _RAS_ 123, 2020, pp. 103312
  * [29] D. Park, Y. Hoshi and C. C. Kemp  “A Multimodal Anomaly Detector for Robot-Assisted Feeding Using an LSTM-Based Variational Autoencoder”  In _IEEE RA-L_ 3.3, 2018, pp. 1544–1551
  * [30] Geoffrey E Hinton  “Training products of experts by minimizing contrastive divergence”  In _Neural computation_ 14.8 MIT Press, 2002, pp. 1771–1800
  * [31] Epic Games  “Unreal Engine”, https://www.unrealengine.com, 2019
  * [32] Jean-François Tremblay et al.  “Automatic 3D Mapping for Tree Diameter Measurements in Inventory Operations” Special issue: Field and Service Robotics 2019  In _JFR_ , 2020  Wiley
  * [33] Alex Kendall and Roberto Cipolla  “Geometric Loss Functions for Camera Pose Regression With Deep Learning”  In _CVPR_ , 2017

  *[RNN]: recurrent neural network
  *[IMUs]: inertial measurement unit
  *[MRSSM]: multimodal recurrent state-space model
  *[MVAE]: multimodal variational autoencoder
  *[ELBO]: evidence lower bound
  *[MVAEs]: multimodal variational autoencoder
  *[RL]: reinforcement learning
  *[MBRL]: model-based reinforcement learning
  *[RSSM]: recurrent state-space model
  *[RSSMs]: recurrent state-space model
  *[VAE]: variational autoencoder
  *[CNN]: convolutional neural network
  *[SLAM]: simultaneous localization and mapping
  *[ICP]: iterative closest point
  *[RMSE]: root mean square error