# Non-principal axis rotation in binary asteroid systems and how it weakens
the BYORP effect

Alice C. Quillen<EMAIL_ADDRESS>Anthony LaBarca
<EMAIL_ADDRESS>YuanYuan Chen<EMAIL_ADDRESS>Department of Physics
and Astronomy, University of Rochester, Rochester, NY 14627, USA  Key
Laboratory of Planetary Sciences, Purple Mountain Observatory, Chinese Academy
of Sciences, Nanjing 210008, China

###### Abstract

Using viscoelastic mass/spring model simulations, we explore tidal evolution
and migration of compact binary asteroid systems. We find that after the
secondary is captured into a spin-synchronous state, non-principal axis
rotation in the secondary can be long-lived. The secondary’s long axis can
remain approximately aligned along the vector connecting secondary to primary
while the secondary rocks back and forth about its long axis. Inward orbital
semi-major axis migration can also resonantly excite non-principal axis
rotation. By estimating solar radiation forces on triangular surface meshes,
we show that the magnitude of the BYORP effect induced torque is sensitive to
the secondary’s spin state. Non-principal axis rotation within the 1:1 spin-
orbit resonance can reduce the BYORP torque or cause frequent reversals in its
direction.

##  1 Introduction

The field of research on binary asteroids currently has two dynamical
conundrums. Light curve based estimates of viscoelastic behavior in single
asteroids based on tumbling damping give elastic modulus times energy
dissipation parameter, $\mu Q\sim 10^{11}$ Pa (Harris, 1994; Pravec et al.,
2014), exceeding those estimated in binary asteroids (Goldreich and Sari,
2009; Taylor and Margot, 2011; Jacobson and Scheeres, 2011; Scheirich et al.,
2015; Nimmo and Matsuyama, 2019) by 1 to 4 orders of magnitude. Drift rates
estimated due to radiation forces, known as binary-YORP or BYORP, are so rapid
that binary NEA asteroid lifetimes should be short, $\lesssim 10^{5}$ years
(Ćuk and Burns, 2005). This short lifetime makes it difficult to account for
the abundance of binary NEA asteroids, about one sixth of NEAs, (Margot et
al., 2002; Pravec et al., 2006) within the NEA asteroid orbital lifetime,
$\sim 10^{7}$ years (Gladman et al., 2000).

With asteroid binary light curves, Pravec et al. (2016) find that synchronous
secondary rotation predominates in close binaries and in binaries with near
circular orbits. Only binary members that occupy a spin-orbit resonance, such
as the synchronous 1:1 state, are expected to have non-zero cumulative BYORP
effects. Otherwise the time-averaged torques are expected to cancel out (Ćuk
and Burns, 2005). The rapid predicted BYORP drift rate, of order a few cm/yr
(Ćuk and Burns, 2005), spurred Ćuk and Nesvorný (2010) to explore scenarios
where the tidal torque is neglected. In contrast, Jacobson and Scheeres (2011)
proposed that BYORP induced orbital drift could be halted by tidal torque with
sufficient tidal dissipation. Orbital eccentricity could cause an elongated
secondary to tumble, delaying BYORP drift-induced evolution and extending the
lifetime of NEA binaries (Naidu and Margot, 2015).

Possible ways to account for the spin synchronous rotation and the inferred
tidal spin-down of the secondary are earlier formation in the Main Belt,
formation at a slower spin rate, formation at a smaller orbital semi-major
axis where the spin down time is shorter, or if the secondary has lower $\mu
Q\sim 10^{9}$ Pa (Taylor and Margot, 2011). Rubble pile NEAs may have low $\mu
Q$ values due to stress on contact points arising from internal voids
(Goldreich and Sari, 2009), a dissipative regolith surface layer $\mu Q\sim
10^{8}$ Pa (Nimmo and Matsuyama, 2019) or because the rheology is better
described with a frequency dependent model (Efroimsky, 2015). These studies
have neglected the possible role of some processes that could affect the spin
down rate, such as mass shedding from the primary (Walsh et al., 2008;
Hirabayashi et al., 2015; Scheeres, 2015) that could give a population of
objects that might impact or accumulate on the secondary (e.g., Davis and
Scheeres 2020), ongoing mass loss, recently detected on asteroid (101955)
Bennu (Lauretta et al., 2019), and mass movements as inferred from asteroid
surface morphology (Jawin et al., 2020).

With multiple epoch photometric observations spanning 20 years, Scheirich et
al. (2015, 2021) have directly measured binary orbital drift rates in 3 binary
asteroid systems. In the case of (88710) 2001 SL9, the BYORP effect is the
only known physical mechanism that can cause the measured $\sim-4$ cm/yr
inward drift of its mutual binary orbit. For NEA (66391) 1999 KW4, known as
Moshup, Scheirich et al. (2021) measured a drift rate of 1.2 cm/s which is
lower than that predicted for the BYORP effect, $\sim 8.5$ cm/s, using the
secondary shape model by Ostro et al. (2006).

Apollo-class Near-Earth Asteroid (65803) Didymos is the target of the
international collaboration known as AIDA (abbreviation for Asteroid Impact &
Deflection Assessment) that supports the development and data interpretation
of the NASA DART (Cheng et al., 2016) and ESA Hera (Michel and et al, 2021)
space missions. AIDA associated missions will provide a unique opportunity to
study a binary asteroid system. The secondary, called Dimorphos, of the
Didymos system is the target of the Double Asteroid Redirection Test (DART)
mission. The DART spacecraft will impact the binary’s secondary fall 2022 to
study the momentum transfer caused by a kinetic impact. The Light Italian
Cubesat for Imaging of Asteroids (LICIA cube; Dotto et al. 2021), carried with
the DART spacecraft, will take in situ observations of the impact. Detailed
physical surface and interior characterization as well as images of the
surface morphology of the two bodies in the Didymos system by the Hera mission
in early 2027 will probe how deformation occurs in a microgravity environment,
giving constraints on binary asteroid internal structure, formation and
evolution mechanisms.

###  1.1 Outline

Inspired by the opportunity represented by DART and associated missions, we
examine the spin evolution of the Didymos binary system. In section 1.2 we
review estimates for the tidal spin down time of both bodies in the system. In
section 2 we explore rotation states and evolution of a binary asteroid system
using our mass/spring model simulation code (Quillen et al., 2019a, b, 2020).
Our simulations are similar to the soft sphere simulations described by Agrusa
et al. (2020), except internal dissipation occurs due to damping in the
springs. Because we are interested in long term evolution, our simulations are
carried out for more rotation periods than those by Agrusa et al. (2020). In
section 2.2 we use mass-spring model viscoelastic simulations to examine the
process of secondary tidal spin down and tidal lock. In section 2.3 we mimic
subsequent BYORP induced migration in the binary orbital semi-major axis
$a_{B}$ through forced migration and examine its affect on the spin dynamics.
In section 3 using triangular surface models we explore how obliquity and non-
principal axis (NPA) rotation affects estimates of the BYORP induced drift
rate. A summary and discussion follows in section 4.

###  1.2 Properties of the Didymos binary system

In Table 1, we list properties of the Didymos binary system. We refer to the
primary as Didymos and the secondary as Dimorphos. Many of the quantities in
this table were measured or recently reviewed by Naidu et al. (2020).

Table 1: Didymos and Dimorphos properties

Parameter Symbol Value Comments/Reference Binary orbital period $P_{B}$
11.9217 $\pm$ 0.0002 hour Scheirich and Pravec (2009) Binary orbital semi-
major axis $a_{B}$ 1.19 $\pm$ 0.03 km Naidu et al. (2020) Binary orbit
eccentricity $e_{B}$ $<0.03$ Scheirich and Pravec (2009) Didymos spin rotation
period $P_{p}$ 2.2600 $\pm$ 0.0001 hours Pravec et al. (2006) Didymos body
diameter $2R_{p}$ 780 $\pm$ 30 m Naidu et al. (2020) Diameter or radius ratio
$R_{s}/R_{p}$ 0.21 $\pm$ 0.01 Scheirich and Pravec (2009) Semi-major
axis/radius primary $a_{B}/R_{p}$ 3.05 Dimorphos diameter $2R_{s}$ 164 $\pm$
18 m From $R_{p}$ and $R_{s}/R_{p}$ measurements Dimorphos spin rotation
period $P_{s}$ 11.92 hours Spin synchronous state assumed Total mass
$M_{p}+M_{s}$ $5.37\pm 0.44\times 10^{11}$ kg Kepler’s third law Mass ratio
$q=M_{s}/M_{p}$ 0.0093 From ratio of radii Mass of Didymos $M_{p}$ $5.32\times
10^{11}$ From mass ratio Mass of Dimorphos $M_{s}$ $4.95\times 10^{9}$ From
mass ratio Mean density of Didymos $\rho_{p}$ $2.14\ {\rm g/cm}^{3}$ From mass
and radius Mean density of Dimorphos $\rho_{s}$ $\rho_{s}=\rho_{p}$ Assumed
Period ratio (orbit/spin primary) $P_{B}/P_{p}$ 5.28 Didymos body axis ratios
$a_{p}:b_{p}:c_{p}$ 1:0.98:0.95 Based on Naidu et al. (2020)
Radii $R_{p},R_{s}$ for the primary and secondary are assumed to be those of
the volume equivalent sphere. We refer to the primary as Didymos and the
secondary as Dimorphos. Dydimos’ body axis ratios are based on the Dynamically
Equivalent Equal Volume Ellipsoid which is an uniform density ellipsoid with
the same volume and moments of inertia as the shape model by Naidu et al.
(2020) and taken from their Table 3. Light curve observations by Pravec et al.
(2006) and radar data by Naidu et al. (2020) are consistent with the secondary
being in the spin synchronous state (tidally locked). The assumption of equal
primary and secondary densities is supported by radar observations of NEA 2000
DP107 (Margot et al., 2002).

Table 2: Didymos and Dimorphos estimated properties energy density | $e_{g,p}$ | 830 Pa
---|---|---
energy density | $e_{g,s}$ | 37 Pa
gravity timescale | $t_{\rm grav,p}$ | 1282 s
breakup spin | $\omega_{\rm breakup,p}$ | $7.8\times 10^{-4}{\rm s}^{-1}$

###  1.3 Tidal spin down rates for primary and secondary

To estimate the tidal torque, we consider a uniform density spherical body
that is part of binary system. The body is initially not tidally locked, so it
is not in the spin synchronous state which is also called the 1:1 spin-orbit
resonance. The spin-down rate $\dot{\omega}$ and a spin down time $t_{\rm
despin}$ can be estimated from the tidally generated torque $\tau_{\rm tidal}$
with $\dot{\omega}_{\rm tidal}=\tau_{\rm tidal}/I$ and $t_{\rm
despin}\sim\dot{\omega}_{\rm tidal}/\omega_{init}$ (e.g., Peale 1977; Gladman
et al. 1996) where the body’s moment of inertia, $I=\frac{2}{5}MR^{2}$, the
body’s initial spin is $\omega_{\rm init}$ and the body’s mass and radius are
$M$ and $R$.

The secular part of the semi-diurnal ($l=2$) term in the Fourier expansion of
the perturbing potential from point mass $M_{pert}$ (the other mass in the
binary), gives a tidally induced torque (e.g., Kaula 1964; Goldreich 1963;
Efroimsky and Makarov 2013)

$\tau_{tidal}=\frac{3}{2}\frac{GM_{pert}^{2}}{a_{B}}\left(\frac{R}{a_{B}}\right)^{5}k_{2}\sin\epsilon_{2}$
(1)

where orbital semi-major axis is $a_{B}$. Here $k_{2}\sin\epsilon_{2}$ depends
on frequency and is known as a quality function. It is often approximated as
$k_{2}/Q$ with energy dissipation or quality factor $Q$ describing the
dissipation in the body with mass $M$, and the Love number $k_{2}$
characterizes the deformation of the body with mass $M$. See Efroimsky (2015)
for discussion on the accuracy of this approximation.

For an incompressible homogeneous spherical elastic body, the Love number

$k_{2}\sim 0.038\frac{e_{g}}{\mu}$ (2)

where $\mu$ is the body’s rigidity (equal to the shear modulus if the body is
incompressible) and

$e_{g}\equiv\frac{GM^{2}}{R^{4}}$ (3)

is a measure of the body’s gravitational energy density, (following Equation 6
by Quillen et al. 2017 and based on Murray and Dermott 1999; Burns and
Safronov 1973). Note that at constant density, $M\propto R^{3}$, the energy
density $e_{g}\propto R^{2}$ and Love number $k_{2}\propto R^{2}$ which gives
a size dependent Love number similar to (but not identical to) the weaker
$k_{2}\propto R$ proposed by Goldreich and Sari (2009).

Using the torque of equation 1, and equation 2, we estimate a spin down rate

$\dot{\omega}_{tidal}\sim 0.14\left(\frac{GM^{2}/R^{4}}{\mu
Q}\right)\left(\frac{R}{a_{B}}\right)^{6}\left(\frac{M_{pert}}{M}\right)\frac{GM_{pert}}{R^{3}}.$
(4)

A particle resting on the surface can leave due to centrifugal acceleration if
the spin $\omega\gtrsim\omega_{\rm breakup}$ with

$\omega_{\rm breakup}\equiv\sqrt{\frac{GM}{R^{3}}}=\sqrt{\frac{4\pi
G\rho}{3}},$ (5)

where $\rho$ is the body’s mean density. If the initial spin is equal to the
that of centrifugal breakup $\omega_{\rm init}=\omega_{\rm breakup}$ then

$t_{\rm despin}\approx
P_{orb}\frac{2}{15\pi}\left(\frac{M}{M_{pert}}\right)^{\frac{3}{2}}\left(\frac{a_{B}}{R}\right)^{\frac{9}{2}}\frac{Q}{k_{2}}\sqrt{\frac{M_{pert}+M}{M_{pert}}}$
(6)

where $P_{orb}={2\pi}\sqrt{\frac{a_{B}^{3}}{G(M_{pert}+M)}}$ is the binary
orbital period.

Using equations 2 and 3 the spin down time is

$\frac{t_{\rm despin}}{P_{B}}\approx
1\left(\frac{M}{M_{pert}}\right)^{\frac{3}{2}}\left(\frac{a_{B}}{R}\right)^{\frac{9}{2}}\left(\frac{\mu
Q}{GM^{2}/R^{4}}\right)\sqrt{\frac{M_{pert}+M}{M_{pert}}}.$ (7)

This equation can be used for either a spinning primary or a spinning
secondary. In the case of a spinning primary

$\displaystyle\frac{t_{\rm despin,p}}{P_{B}}$
$\displaystyle\approx\left(\frac{M_{p}}{M_{s}}\right)^{\frac{3}{2}}\left(\frac{a_{B}}{R_{p}}\right)^{\frac{9}{2}}\left(\frac{\mu_{p}Q_{p}}{GM_{p}^{2}/R_{p}^{4}}\right)\sqrt{\frac{M_{p}+M_{s}}{M_{s}}}$
(8)

where subscript $p$ refers to quantities for the primary, subscript $s$ refers
to quantities for the secondary and subscript $B$ refers to quantities
describing the binary mutual orbit. In the case of a spinning secondary

$\displaystyle\frac{t_{\rm despin,s}}{P_{B}}$
$\displaystyle\approx\left(\frac{M_{s}}{M_{p}}\right)^{\frac{3}{2}}\left(\frac{a_{B}}{R_{s}}\right)^{\frac{9}{2}}\left(\frac{\mu_{s}Q_{s}}{GM_{s}^{2}/R_{s}^{4}}\right)\sqrt{\frac{M_{p}+M_{s}}{M_{p}}}.$
(9)

Wobble damping timescales estimated from asteroid light curves give an
estimate for asteroid viscoelastic material properties $\mu Q\sim 10^{11}$ Pa
(Harris, 1994; Pravec et al., 2014). Using this value for primary and
secondary we estimate

$\displaystyle t_{\rm despin,p}$ $\displaystyle\approx 2.4\times 10^{11}\ {\rm
yr}\left(\frac{\mu_{p}Q_{p}}{10^{11}\
{\rm{Pa}}}\right)\left(\frac{a_{B}}{1.19\ \mathrm{km}}\right)^{3}$
$\displaystyle\times\left(\frac{q}{0.0093}\right)^{-2}\left(\frac{R_{p}}{390\
{\rm m}}\right)^{-5}\left(\frac{\rho_{p}}{2.1\ {\rm g\
cm}^{-3}}\right)^{-\frac{3}{2}}$ (10) $\displaystyle t_{\rm despin,s}$
$\displaystyle\approx 5.8\times 10^{8}\ {\rm
yr}\left(\frac{\mu_{s}Q_{s}}{10^{11}\
{\rm{Pa}}}\right)\left(\frac{a_{B}}{1.19\ \mathrm{km}}\right)^{3}$
$\displaystyle\times\left(\frac{q}{0.0093}\right)^{2}\left(\frac{R_{s}}{82\
{\rm m}}\right)^{-5}\left(\frac{\rho_{s}}{2.1\ {\rm g\
cm}^{-3}}\right)^{-\frac{3}{2}}$ (11)

using quantities measured for Didymos and Dimorphos listed in Table 1. Here
$q=M_{s}/M_{p}$ is the secondary to primary mass ratio.

The long primary spin-down time in Equation 1.3 implies that tidal dissipation
does not affect the primary’s spin. All observed near-Earth and small main-
belt asteroid binaries with secondary to primary diameter ratio lower than 0.6
have quickly spinning primaries with periods of a few hours (Pravec and
Harris, 2007; Pravec et al., 2016). On an asteroid, the reflection, absorption
and emission of solar radiative energy produces a torque that changes the
rotation rate and obliquity of a small body. The Yarkovsky-O’Keefe-
Radzievskii-Paddack (YORP) effect can account for the primary high spin rates
(Rubincam, 2000; Scheeres, 2015; Hirabayashi et al., 2015). YORP could also
slow down the primary’s spin, though at slow spin rates, the binary system may
become unstable and disrupt, leading to fewer binaries with slowly spinning
primaries. For an illustration of an instability causing binary disruption,
see Figure 14 by Davis and Scheeres (2020).

The timescale in equation 11 for spin down time of the secondary is longer
than the estimated NEA lifetime and is in conflict with the large number of
binary asteroids with synchronous secondaries unless the binary asteroids were
born in the main belt before they became NEAs or the tidal drift rate is
underestimated by equation 11 (Taylor and Margot, 2011), as we discussed in
section 1.

Constraints on asteroid internal properties based on tidal evolution of the
orbital semi-major axis are only relevant if the BYORP induced drift rate is
low (e.g., Taylor and Margot 2011) or if the tidal dissipation is sufficiently
strong that the tidal torque can counter the BYORP torque (Jacobson and
Scheeres, 2011). The median value for small main belt binaries (primary radius
$\sim 3$ km), estimated with a collisional lifetime of $\Delta t=10^{9}$ yr,
is $\mu_{p}Q_{p}\sim 10^{11}$ Pa whereas that for small NEA binaries (with
primary radius less than 1 km), estimated with a lifetime of $\Delta t=10$ Myr
is $\mu_{p}Q_{p}\sim 5\times 10^{8}$ Pa (Taylor and Margot, 2011). Assuming an
equilibrium state between BYORP and tidal torque, Jacobson and Scheeres (2011)
estimate $BQ_{p}/k_{2p}\sim 10^{3}$ giving $Q_{p}/k_{2p}\sim 10^{6}$. With
BYORP coefficient $B\sim 10^{-3}$ this gives $Q_{p}\mu_{p}\sim 10^{8}$ Pa,
similar to that estimated for NEA asteroid binaries by Taylor and Margot
(2011). (Here we used equation 2 for the Love number, and estimated
$0.038G\left(4\pi/3\rho_{p}R_{p}\right)^{2}\sim 10^{2}$ Pa with values
$R_{p}=1$ km and $\rho_{p}=2$ g/cm3.)

As emphasized by Wisdom (1987), a non-round body must tumble before entering
the spin-synchronous state. Numerical simulations show that the time for
obliquity and libration and non-principal axis rotation to damp can
substantially prolong the time it takes to reach a spin synchronous state with
low obliquity, low free libration and undergoing principal axis rotation
(Naidu and Margot, 2015; Quillen et al., 2020). Hence the long spin down time
estimated for Dimorphos in equation 11 suggests that Dimorphos could have
spent a long period of time in a complex spin state that is nearly synchronous
but has non-principal axis (NPA) rotation, obliquity oscillations and high
amplitude libration.

Epochs of chaotic tumbling could prolong the lifetime of a binary asteroid as
BYORP induced drift in the tumbling state is expected to be negligible. Naidu
and Margot (2015) pointed out that the libration amplitude can be difficult to
detect from a light curve. With libration angle equal to the angle between the
secondary’s long axis and the vector between primary and secondary centers of
mass, Pravec et al. (2016) showed that libration angles greater than about
$20^{\circ}$ can be detected through analysis of light-curve measurements.
Pravec et al. (2016) attributed small offsets between the secondary rotational
light curve and models to irregular secondary shape. While secondary light
curves are sensitive to the orientation of the secondary’s long axis, they may
be less sensitive to the orientation of the secondary’s short axis.

Figure 1: A snap-shot of the mass/spring model simulation of a binary asteroid
with output shown in Figure 2. The green lines show springs and the gray
spheres show mass nodes. Both primary and secondary shapes are resolved with
mass nodes.

##  2 Mass/spring model simulations of an asteroid binary

To explore the spin evolution of the two non-spherical bodies in an asteroid
binary, we use the mass-spring model simulation code (Quillen et al., 2016a;
Frouard et al., 2016; Quillen et al., 2016b, 2017, 2019a, 2019b, 2020), that
is built on the modular N-body code rebound (Rein and Liu, 2012). In our code,
a viscoelastic solid is approximated as a collection of mass nodes that are
connected by a network of springs. Springs between mass nodes are damped and
the spring network approximates the behavior of a Kelvin-Voigt viscoelastic
solid.

In our previous work we only resolved a single spinning body with masses and
springs, though the body could be perturbed by point masses (e.g., Quillen et
al. 2017). Here we resolve both primary and secondary in the asteroid binary
with masses and springs. We do not include the effect of planets or tidal
perturbations associated with the heliocentric orbit. As we are interested in
the behavior of long integrations (as when we studied the obliquity of Pluto
and Charon’s satellites; Quillen et al. 2017) rather than details of internal
tidal dissipation (Quillen et al., 2019a) we only use a dozen mass nodes to
resolve the secondary and about a hundred mass nodes to resolve the primary.
The bodies are not spherical. Because we resolve the bodies with individual
mass nodes, we do not need to expand their gravitational perturbations in
multipole terms. See Figure 1 for a snap-shot of one of our binary
simulations.

The mass particles or nodes in both resolved spinning bodies are subjected to
three types of forces: the gravitational forces acting on every pair of mass
nodes in the simulation, and the elastic and damping spring forces acting only
between sufficiently close particle pairs. As all forces are applied equal and
oppositely to pairs of point particles and along the direction of the vector
connecting the pair, momentum and angular momentum conservation are assured.

We work with mass in units of $M_{p}$, the mass of the primary in the asteroid
binary, distances in units of volumetric radius, $R_{\rm vol,p}$, the radius
of a spherical body with the same volume as the primary, and time in units of
$t_{\rm grav}$

$\displaystyle t_{\rm grav}$ $\displaystyle\equiv\sqrt{\frac{R_{\rm
vol,p}^{3}}{GM_{p}}}=\sqrt{\frac{3}{4\pi G\rho_{p}}}.$ (12)

This timescale at the density of the primary for the Didymos system is listed
in Table 2.

In our simulations, the primary and secondary are approximately the same
density. Spin is given in units of $\omega_{\rm breakup}=t_{\rm grav}^{-1}$
(see equation 5) and computed using the mean density of the primary which is
$3/(4\pi)$ in our numerical units. This frequency is also listed in Table 2
for the estimated density of Didymos.

Initial node distribution and spring network for each body is chosen with the
tri-axial ellipsoid random spring model described by Frouard et al. (2016);
Quillen et al. (2016b). The confining surface used to generate the particles
obeys $\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$ with
$a,b,c$ equal to half the lengths of the desired principal body axes.
Particles cannot be closer than a minimum distance $d_{I}$ from each other.
Since we generate the resolved bodies with only a few particles, axis ratios
are more accurately computed using the eigenvalues of the moment of inertia
matrix after we generate the particle distribution. For principal moments of
inertia $A\geq B\geq C$, we compute the body axis ratios as
$\frac{b}{a}=\sqrt{\frac{A+C-B}{A+B-C}}$ and
$\frac{c}{a}=\sqrt{\frac{B+C-A}{A+B-C}}$.

Once the particle positions have been generated, every pair of particles
within distance $d_{Spr}$ of each other are connected with a single spring.
The parameter $d_{Spr}$ is the maximum rest length of any spring in the
network. Springs have a spring constant $k_{Spr}$, giving elastic behavior and
a damping rate parameter $\gamma_{Spr}$ giving viscous behavior (Frouard et
al., 2016). Our simulated material is compressible, so energy damping arises
from both deviatoric and volumetric stresses. Springs are created at the
beginning of the simulation and do not grow or fail during the simulation, so
there is no plastic deformation. For numerical stability, the simulation time-
step must be chosen so that it is shorter than the time it takes elastic waves
to travel between nodes (Frouard et al., 2016; Quillen et al., 2016b). All
mass nodes in the primary have the same mass and all springs in the primary
have the same spring constant and damping parameter. All mass nodes in the
secondary have the same mass and all springs in the secondary have the same
spring constant. However, the mass nodes in the primary are not the same mass
as those in the secondary. Similarly the spring constants of the springs in
the primary differ from those in the secondary.

###  2.1 Simulation output

At the beginning of the simulation we store the positions of the mass nodes
for both primary and secondary. The positions are in a coordinate system with
$x$ axis aligned with the long principal body axis and $z$ axis aligned with
the short principal body axis. The position of node $j$ from the center of
mass of the body in which it resides is ${\bf
r}^{0}_{j}=(r^{0}_{j,0},r^{0}_{j,1},r^{0}_{j,2})$ where $r^{0}_{j,0}$ is the
$x$ coordinate. The position of node $j$ from the center of mass of the body
in which it resides at a later time during the simulation is ${\bf
r}_{j}=(r_{j,0},r_{j,1},r_{j,2})$. At each simulation output we compute the
covariance matrix $S_{\alpha\beta}=\sum_{j}r^{0}_{j,\alpha}r_{j,\beta}$ using
the current node positions and the initial node positions. The covariance
matrix is stored so that we can later compute a rotation describing the body
orientation at each simulation output using the Kabsch algorithm (Kabsch,
1976; Berthold and Horn, 1986; Coutsias et al., 2004). At each simulation
output the current body orientation is given by a rotation that we describe
with a unit quaternion. The quaternion rotates the initial body node positions
to give the node positions at the simulation output time. This quaternion is
equivalent to the unit eigenvector associated with the maximum eigenvalue of a
4x4 matrix derived from the covariance matrix $S_{\alpha\beta}$ (Berthold and
Horn, 1986; Coutsias et al., 2004). From the covariance matrix at each
simulation output we compute a quaternion at each simulation output time. We
use this quaternion to compute the orientation of the body’s long and short
axes by rotating the $x$ and $z$ axes at each simulation output time.

At each simulation output time we compute the primary and secondary’s moment
of inertia matrix from the positions of their mass nodes with respect to the
body’s center of mass. The spin angular momenta for each body, ${\bf
L}_{s},{\bf L}_{p}$, are computed from the positions and velocities of each
body’s mass nodes. The spin vectors for each body,
${\boldsymbol{\omega}}_{s},{\boldsymbol{\omega}}_{p}$ are computed from the
body spin angular momentum vectors and the inverse of their moment of inertia
matrices, e.g., ${\boldsymbol{\omega}}_{s}={\bf I}_{s}^{-1}{\bf L}_{s}$.

We denote the directions of the secondary’s long and short body principal axes
as $\hat{\bf i}_{s}$ and $\hat{\bf k}_{s}$, respectively. The intermediate
secondary body principal axis is $\hat{\bf j}_{s}$. The inclination of the
secondary’s long principal axis is the angle between $\hat{\bf i}_{s}$ and the
binary orbital plane. The secondary’s non-principal angle $\theta_{\rm NPA,s}$
is the angle between the secondary’s short principal axis $\hat{\bf k}_{s}$
and the direction of the secondary’s spin angular momentum, ${\bf L}_{s}$. The
secondary’s short axis tilt is the angle between the secondary’s short axis
and the orbit normal, ${\rm arccos}(\hat{\bf k}_{s}\cdot\hat{\bf l}_{B})$. The
orbit normal is close to the $z$ axis during the simulation. The secondary’s
precession angle $\theta_{ls}$ is the angle on the $xy$ plane of the
secondary’s spin angular momentum vector, $\theta_{ls}={\rm
arctan2}(L_{sy},L_{sx}$) where ${\bf L}_{s}=(L_{sx},L_{sy},L_{sz})$. The
obliquities for primary and secondary $\epsilon_{s},\epsilon_{p}$ are the
angle between body spin angular momentum and orbital normal $\hat{\bf l}_{B}$.
The secondary libration angle $\phi_{lib,s}$ is that between the secondary’s
long principal body axis, projected onto the orbital plane, and the vector
connecting primary and secondary centers of mass.

The binary orbit normal $\hat{\bf l}_{B}$ direction is computed using the body
center of mass positions and velocities. The orbital semi-major axis $a_{B}$,
eccentricity $e_{B}$ and inclination $i_{B}$ are computed from the positions
and velocities of the primary and secondary center of mass and assuming
Keplerian motion (we assume point masses so the extended mass distributions
are not taken into account). The orbital inclination is zero if the binary
orbit lies in the $xy$ plane. We measure the angle of the secondary in its
orbit, $\theta_{B}$, from the vector between the two body center of mass
positions, $\bf r$, projected into the $xy$ plane. As orbital inclination and
eccentricity are low, $\theta_{B}={\rm arctan2}(r_{y},r_{x})$ is approximately
equal to the binary’s mean longitude, $\lambda_{B}$. Using an epicyclic
approximation we measure an epicylic angle $\theta_{kB}$ describing radial
oscillations of the orbit. The epicyclic angle is computed from the vector
between secondary and primary, the radius $r=|{\bf r}|$, the mean value of the
radius, $\bar{r}$, the orbital mean motion $n_{B}$, and the radial velocity
component, $v_{r}=\frac{dr}{dt}$. The mean radius $\bar{r}$ and mean motion
are estimated by smoothing the time series of $r$ and $n_{B}$ with a Savinsky-
Golay filter that has a window length of a few orbits. The epicyclic angle is
then $\theta_{kB}={\rm arctan2}((r-\bar{r})n_{B},v_{r})$. The angles
$\theta_{B},\theta_{kB}$ are used instead of orbital elements because a
Keplerian approximation is poor due to the non-spherical shapes of the primary
and secondary (e.g., Renner and Sicardy 2006; Agrusa et al. 2020).

####  2.1.1 Multipole coefficients

We describe our convention for calculating gravitational multipole
coefficients for the primary and secondary. Following Kaula (1966); Tricarico
et al. (2021), the gravitational potential outside a massive non-spherical
body, of mass $M$, where it satisfies Laplace’s equation, can be expanded in
terms of multipoles,

$\displaystyle V(r,\theta,\phi)$
$\displaystyle=-\frac{GM}{r}\Bigg{[}1+\sum_{l=2}^{\infty}\sum_{m=0}^{l}\left(\frac{R_{vol}}{r}\right)^{l}P_{lm}(\cos\theta)\times$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\left(C_{lm}\cos(m\phi)+S_{lm}\sin(m\phi)\right)\Bigg{]}$ (13)

where $r$ is radius, $\phi\in[0,2\pi]$ is longitude and $\theta\in[0,\pi]$ is
colatitude. The center of mass is at the origin and this is why the $l=1$
terms are absent from the expansion. The functions $P_{lm}()$ are the
associated Legendre polynomials which can be computed using Rodrigues’ formula

$P_{lm}(x)=\frac{(1-x^{2})^{m/2}}{2^{l}l!}\frac{d^{l+m}}{dx^{l+m}}(x^{2}-1)^{l}.$
(14)

The dimensionless multipole coefficients, $C_{lm},S_{lm}$, depend on a
characteristic length which we set to be that of the volume equivalent sphere
$R_{vol}$. Each mass node in our simulations is specified with an index $i$,
has mass $m_{i}$ and spherical coordinates $(r_{i},\theta_{i},\phi_{i})$. The
multipole coefficients of the body are

$\displaystyle C_{lm}$
$\displaystyle=\frac{1}{M}\frac{(l-m)!}{(l+m)!}(2-\delta_{m0})\sum_{i}m_{i}\left(\frac{r_{i}}{R_{vol}}\right)^{l}P_{lm}(\cos\theta_{i})\cos(m\phi_{i})$
$\displaystyle S_{lm}$
$\displaystyle=\frac{1}{M}\frac{(l-m)!}{(l+m)!}2\sum_{i}m_{i}\left(\frac{r_{i}}{R_{vol}}\right)^{l}P_{lm}(\cos\theta_{i})\sin(m\phi_{i})$
(15)

where total mass $M=\sum_{i}m_{i}$. The $l=0$ coefficients are called zonal
coefficients, $J_{l}=-C_{l0}$, and the $l\neq 0$ coefficients are called
tesseral coefficients.

Table 3: Viscoelastic Mass/Spring Model Simulation Parameters for Tidal spin
down

Description | Symbol | Quantity |
---|---|---|---
Total integration time | $t_{max}$ | $5\times 10^{5}=15400\ P_{B}$ |
Time step | $dt$ | 0.002 |
Initial orbital semi-major axis | $a_{B}$ | 3.0 |
Description | Symbols | For Primary | For Secondary
Mass | $M_{p},M_{s}$ | 1 | 0.01
Radius equiv. volume sphere | $R_{vol}$ | 1 | 0.22
Initial spin | $\omega_{init}$ | 0.8 | 0.5
Minimum interparticle distances | $d_{I}$ | 0.27 | 0.14
Maximum spring length/$d_{I}$ | $d_{Spr}/d_{I}$ | 2.4 | 2.3
Spring constants | $k_{Spr}$ | 0.20 | 0.03
Spring damping parameters | $\gamma_{Spr}$ | 1 | 30
Body axis ratios | $a:b:c$ | 1.0:0.958:0.881 | 1.0:0.773:0.698
Number of particles | $N$ | 165 | 18
Number of springs | $N_{Spr}$ | 1950 | 113

Axis ratios are those derived from the principal axis moments of inertia and
that are computed from the generated particle positions. Time is given in
units of $t_{\rm grav}$. Distances are in units of $R_{p}$, the primary’s
volume equivalent radius. Spin rates are in units of the break up spin rate
$\omega_{\rm breakup}$ or equivalently $t_{\rm grav}^{-1}$. Multipole
coefficients for primary and secondary are listed in Table 4.

###  2.2 Simulation of secondary tidal spin down

We carry out a simulation with simulation parameters listed in Table 3. The
secondary is initially spinning quickly and is not in a spin-orbit resonance.
Secondary to primary mass ratio, body axis ratios, primary spin, and binary
orbital semi-major axis are intended to be similar to those of the Didymos
system. The numbers of mass nodes and springs in each body and ellipsoid axis
ratios measured from the body’s moments of inertia are listed in Table 3. For
this simulation we set the damping parameters for the springs to be at
relatively high values, following Quillen et al. (2017, 2020), so that we can
see the process of tidal spin down due to dissipation.

Both bodies were initially set to be rotating about their shortest principal
axis and at zero obliquity. Here obliquity is the angle between body spin
angular momentum and binary orbit normal vector direction $\hat{\bf l}_{B}$.
The binary orbit is initially in the $xy$ plane with orbit normal $\hat{\bf
l}_{B}=\hat{\bf z}$ in the $z$ direction. We set the spin of the primary to be
$\omega_{p}=0.8\ \omega_{\rm breakup}=0.8\ t_{\rm grav}^{-1}$, below the near
1 value that Didymos currently has, to ensure that self-gravity can hold the
primary together while keeping the springs under slight compression. As
springs do not dissolve in our simulation, a rapidly spinning body near
rotational breakup would unphysicaly put the springs under tension. The
secondary is initially at $\omega_{s}=0.5\ \omega_{\rm breakup}$.

The binary orbit was begun at low orbital eccentricity. The initial orbit
eccentricity $e_{B}$ is not exactly zero because we compute it without taking
into account the non-spherical mass distributions of primary or secondary (for
discussion on this issue see Agrusa et al. 2020).

Figure 2: Simulation of a tidally evolving asteroid binary. Both primary and
secondary are extended non-spherical bodies. The parameters for this
simulation are listed in Table 3 and a simulation snap-shot is shown in Figure
1. The top panel shows primary and secondary spin rates in units of $t_{\rm
grav}^{-1}$. The second panel from the top shows primary and secondary
obliquity (which is the angle between spin angular momentum and orbit normal).
The third panel from top shows the secondary libration angle, $\phi_{lib,s}$.
In the middle panels, we show the inclination of the secondary’s long body
axis with respect to the orbital plane, the secondary’s non-principal axis
angle $\theta_{\rm NPA,s}$, the tilt of the secondary’s short axis with
respect to the orbit normal ${\rm acos}(\hat{\bf k}\cdot\hat{\bf l}_{B})$, a
resonant angle associated with precession and a resonant angle associated with
libration. The second from bottom panel shows the binary orbit semi-major axis
in units of the primary’s radius, $a_{B}/R_{p}$. The bottom panel shows binary
orbital eccentricity $e_{B}$ and inclination $i_{B}$. Angles are in degrees
except for the binary orbit inclination in the bottom panel which is in
radians. Time is in units of binary’s initial orbital period. This simulation
illustrates that the spin-down rate of the secondary as it approaches the spin
synchronous state (shown in the top panel at early times) is faster than the
obliquity and non-principal axis rotation damping rate at later times.

In Figure 2, we plot primary and secondary spin in units of $t_{grav}^{-1}$,
obliquity, secondary libration angle, orbital semi-major axis, eccentricity
and inclination as a function of time. The horizontal axis shows time in units
of the initial binary orbital period, $P_{B}$. All angles are in degrees
except for the orbital inclination $i_{B}$ which is in radians. In this figure
we also plot the secondary’s long axis inclination angle, non-principal axis
angle $\theta_{\rm NPA,s}$, short axis tilt angle, and resonant angles
associated with precession and libration.

In the top panel of Figure 2, we plot both primary and secondary spins. Dotted
lines show the locations of spin-orbit resonances where the spin is equal to
the binary orbit mean motion $\omega=n_{B}/j$ divided by integer $j$. The
secondary enters the 1:1 spin-orbit resonance at about $t\approx 1000\ P_{B}$
at which time the libration angle $\phi_{lib,s}$ begins to librate about 0.
Simultaneously, the secondary’s obliquity jumps to about $\epsilon_{s}\sim
40^{\circ}$ and non-principal axis rotation is excited. The non-principal axis
rotation is reflected in the large range exhibited by the short axis tilt
angle. A smaller jump in obliquity occurs at about $t\sim 400\ P_{B}$ where
the 2:1 spin-orbit resonance (where $\omega_{s}\approx 2n_{B}$) is crossed.
Because tidal dissipation also is present in the primary, the primary spin
slowly decreases during the simulation and this is accompanied by a slow
increase in the orbital semi-major axis.

####  2.2.1 Non-principal axis rotation and the time it takes for it to decay

The simulation shown in Figure 2 illustrates that the spin-down rate of the
secondary (shown in the top panel) is faster than the obliquity and non-
principal axis damping rate. Past a time of about 1000 orbit periods, the
secondary spin rate on average is consistent with being in the 1:1 spin-orbit
resonance, with secondary spin rate equal to the binary orbit mean motion
$\omega_{s}\sim n_{B}$. Long after ($10^{4}$ orbit periods after) the
secondary enters the 1:1 spin-orbit resonance, its obliquity and non-principal
axis rotation angle $\theta_{\rm NPA,s}$ remain high and vary chaotically.
Long lasting obliquity and non-principal axis angle variations were seen in
some simulations of tidally evolving Phobos and Deimos (Quillen et al., 2020),
though not with as high values of $\theta_{NPA}$. We confirm the findings of
Naidu and Margot (2015) who similarly emphasized that chaotic non-principal
axis rotation of the secondary could be long-lived in binary asteroid systems,
though they discussed tumbling rather than NPA rotation within the 1:1 spin
orbit resonance.

Figure 3: An illustration of the non-principal axis rotation state of the
secondary for the simulation shown in Figure 2 after the secondary has entered
the 1:1 spin-orbit resonance. The secondary’s long axis (shown with direction
$\hat{\bf i}_{s}$) remains aligned along the direction to the primary and this
maintains a low libration angle $\phi_{lib,s}$. However, the secondary’s short
axis, with direction $\hat{\bf k}_{s}$, can rotate as shown with red arrow.
The secondary can slowly rotate about its long axis without leaving the 1:1
spin-orbit resonance. The short axis tilt angle is the angle between $\hat{\bf
k}_{s}$ and the binary orbit normal $\hat{\bf l}_{B}$. The short axis can rock
back and forth, as shown by the red arrow, due to non-principal axis rotation.
Replacing the secondary with an airplane that has nose pointing to the
primary, long axis libration is equivalent to yaw, short axis tilt is
equivalent to roll and variations in the long axis inclination are equivalent
to pitch.

The nature of the secondary’s non-principal axis rotation state in the
simulation shown in Figure 2 is illustrated in Figure 3. The secondary’s long
body axis remains aligned with the direction to the primary giving it a low
libration angle that oscillates about zero. Meanwhile the short axis of the
secondary can rotate. The body can rotate about its long axis without
significantly increasing the libration angle. The secondary lies in the 1:1
spin-orbit resonance but also exhibits significant non-principal axis
rotation. The light curve of such a secondary might appear to reside in a low
obliquity and tidally locked principal axis rotation state. It might be
difficult to differentiate between a principal axis rotation state and a non-
principal axis rotation state if the libration angle $\phi_{lib,s}$ remains
low in both states.

In a frame synchronous with the orbit and replacing the secondary with a plane
that has nose pointing to the primary, long axis libration is equivalent to
yaw, short axis tilt is equivalent to roll and variations in the long axis
inclination are equivalent to pitch. Independently and concurrent with our
study, Agrusa et al. (2021) have also noticed that within the 1:1 spin-orbit
resonance, an elongated secondary can rotate about its long axis. This may be
related to a phenomenon called the “Barrel instability” (Ćuk et al., 2020).

Tidal spin down in our simulation takes place in 1000 orbital periods,
corresponding to only 1.4 years at the binary orbital period of Didymos. This
is about 8 orders of magnitude shorter than estimated for Dimorphos’ spin down
via equation 9. The ratio of the tidal spin down time to the binary orbital
period of 1000 is long enough that the system could be varying adiabatically.
However, the tidal dissipation rate is fast enough that the simulation may not
exhibit phenomena associated with weaker resonances that only would be noticed
at lower drift rates (Quillen, 2006). Below we compare the non-principal axis
rotation damping time to that of spin down. Even if we are not well into the
adiabatic regime, the two phenomena may exhibit similar sensitivity to the
tidal dissipation rate. Multiple simulations of the same body (with axis
ratios similar to those of our secondary) that were begun at slightly
different initial conditions showed a wider variety of phenomena in slowly
drifting simulations than in more rapidly drifting ones (Quillen et al.,
2020). Future work could check that longer more adiabatically drifting
simulations show similar phenomena as we find here and could quantitatively
explore the diversity of possible secondary NPA rotation exhibited within the
1:1 resonance.

The study of Phobos and Deimos, which compared tumbling decay at different
tidal dissipation rates (Quillen et al., 2020), found that non-principal axis
rotation is sometimes prolonged at lower tidal drift rates, however, the ratio
of time spent within the 1:1 spin orbit resonance at high levels of NPA
rotation to that required for spin down was typically much shorter than we see
in our binary asteroid simulation. The simulated Phobos and Deimos have
similar body axis ratios as the secondary simulated here, so body shape does
not account for the difference in length of time spent at high levels of NPA
rotation within the 1:1 resonance. The Phobos and Deimos simulations shown in
Figures 4 and 5 by Quillen et al. (2020) have a larger ratio of orbital semi-
major axis to moon radius ($a/R_{moon}\sim 48$) than the equivalent ratio for
the binary asteroid simulations presented here ($a_{B}/R_{s}\sim 14$, which is
similar to the actual value). These Phobos and Deimos simulations were done
with Mars modeled by a point mass. However, our simulated asteroid secondary
orbits a non-round primary with ratio of orbital major axis to primary radius
$a_{B}/R_{p}\sim 3$ which is similar to the actual value for the Didydmos
system. The ratios $a_{B}/R_{s}$ and $a_{B}/R_{p}$ are low enough in binary
asteroid systems that the non-spherical body shapes could affect the spin
dynamics (Agrusa et al., 2020). Higher order multipole components in the
gravitational potential are effectively stronger in the binary asteroid than
the Phobos/Deimos setting due to the non-spherical shapes of both primary and
secondary and their proximity. Mars has zonal coefficient $J_{2}\approx
8.75\times 10^{-4}$ (Genova et al., 2016), whereas asteroid (101955) Bennu has
$J_{2}\approx 0.01926$ (Chesley et al., 2020). The multipole coefficients of
the primary in a binary asteroid system are likely to be more than an order
larger than those of Mars.

Because the multipole terms in the gravitational field may be important we
have computed them using equations 15 for primary and secondary in the
simulation shown in Figures 2 and 4 and list them for $l=2$ to 5 in Table 4.
The coefficients are computed with origin at the body’s center of mass, with
$x$-axis oriented along the longest body principal axis and $z$-axis oriented
along the shortest principal axis. This orientation gives
$C_{21}=S_{21}\approx 0$. Zonal coefficients $J_{l}=-C_{l0}$.

Table 4: Multipole coefficients of primary and secondary

primary $C_{20}=$-4.6333E-02 $C_{21}=$1.2341E-08 $C_{22}=$5.1758E-03
$S_{21}=$-1.3732E-08 $S_{22}=$-1.2697E-08 $C_{30}=$8.6079E-03
$C_{31}=$2.0366E-03 $C_{32}=$3.3490E-04 $C_{33}=$4.1048E-04
$S_{31}=$-3.3649E-03 $S_{32}=$-4.1383E-04 $S_{33}=$5.3149E-04
$C_{40}=$2.3331E-03 $C_{41}=$-2.7991E-04 $C_{42}=$-4.5127E-04
$C_{43}=$-5.1480E-05 $C_{44}=$-2.2386E-05 $S_{41}=$-8.2225E-04
$S_{42}=$-2.9561E-04 $S_{43}=$1.2788E-04 $S_{44}=$-1.1213E-05
$C_{50}=$5.8134E-03 $C_{51}=$-9.6107E-04 $C_{52}=$-6.0360E-05
$C_{53}=$5.0140E-05 $C_{54}=$2.4168E-05 $C_{55}=$6.6303E-06
$S_{51}=$1.3849e-03 $S_{52}=$4.9828E-06 $S_{53}=$-4.1337E-05
$S_{54}=$-1.4656E-05 $S_{55}=$3.8767E-06 secondary $C_{20}=$-1.2240E-01
$C_{21}=$-1.8134E-07 $C_{22}=$3.9505E-02 $S_{21}=$-1.1728E-07
$S_{22}=$9.9006E-08 $C_{30}=$1.6585E-02 $C_{31}=$1.6529E-02
$C_{32}=$6.9445E-03 $C_{33}=$9.8851E-05 $S_{31}=$6.5180E-03
$S_{32}=$2.2701E-03 $S_{33}=$-1.1486E-03 $C_{40}=$2.1781E-02
$C_{41}=$-9.5380E-03 $C_{42}=$-3.3054E-03 $C_{43}=$-1.6999E-03
$C_{44}=$2.9416E-04 $S_{41}=$-7.2263E-03 $S_{42}=$-1.0916E-03
$S_{43}=$1.7710E-03 $S_{44}=$1.5515E-04 $C_{50}=$1.8595E-02
$C_{51}=$-6.2283E-03 $C_{52}=$-1.0997E-03 $C_{53}=$6.9053E-04
$C_{54}=$2.8126E-04 $C_{55}=$2.8868E-05 $S_{51}=$6.0582E-03
$S_{52}=$2.0111E-04 $S_{53}=$1.5493E-04 $S_{54}=$7.4877E-05
$S_{55}=$7.5874E-05

The primary in our simulations has multipole coefficients that have similar
magnitude to those of asteroid Bennu (in Table 5 by Chesley et al. 2020). The
multipole coefficients of the secondary are larger than those of the primary
because it is more elongated and because it is resolved in the simulation with
fewer mass nodes. Are the multipole coefficients sufficiently strong that they
could affect the spin dynamics of the secondary in our simulations? The
primary is spinning quickly. Perturbations on the secondary associated with
the multipoles of the primary would have frequencies that exceed the
secondary’s spin. However, via a higher order expansion, perturbations with
frequencies that are sums or differences of other perturbation frequencies (as
are commonly found in second order expansions) might be similar to the
rotation rate and so might account for the chaotic behavior. We tentatively
attribute the long lived non-principal axis rotation seen in our simulations
to the strength of perturbations associated with multipole potential terms
that are present because both primary and secondary are non-spherical.

Figure 2 shows that the time to damp the non-principal axis rotation and
obliquity in our simulation is longer than the time it takes for the secondary
to spin down and enter the spin synchronous state. Fitting time dependent
exponentially decreasing functions to smoothed secondary spin and obliquity
arrays, we find that the damping time is about 5 times longer for obliquity
after reaching the 1:1 spin-orbit resonance than for spin prior to reaching
the resonance.

As discussed in section 1.3, the time for the secondary to spin down could be
similar to the lifetime of NEA asteroids. With a value considered typical of
main belt asteroids, the tidal spin down time for Dimorphos exceeds the NEA
lifetime. If binary asteroid secondaries are predominantly in principal axis
rotation states, then much lower values of $\mu Q$ would be required to
account for their synchronous rotation. This follows because it takes longer
for tidal torque to reduce the level of NPA rotation after reaching the 1:1
spin orbit resonance. In our simulation, we have compared the decay time of
the non-principal axis rotation to the tidal spin down time. Assuming that the
tidal dissipation rate in the simulations is sufficiently adiabatic that the
ratio of the two time scales is insensitive to the tidal dissipation rate, we
infer that it takes even longer for NPA rotation to decay than to enter the
1:1 spin orbit resonance. The even longer time require for non-principal axis
rotation damping exacerbates the difficulty reconciling the estimates for $\mu
Q$ in binary and single asteroids.

What if the asteroid binary secondaries are commonly within complex rotation
states with non-principal axis rotation? They could still reside in the 1:1
spin-orbit resonance while exhibiting non-principal axis rotation, as seen in
our simulation shown in Figure 2. If binary asteroid secondaries are commonly
found in non-principal axis rotation states then their tidal dissipation rate
cannot be extremely high. Constraints on secondary rotation state may improve
upon tidal based estimates for their material properties.

In Figure 2, in the panel labeled ‘precession’ (fourth from bottom), we show a
resonant angle $\theta_{ls}-\theta_{B}$ which is the secondary spin angular
momentum direction projected onto the orbital plane as seen in the frame
rotating with the secondary in its orbit. At $t\sim 9500\ P_{B}$ this resonant
angle is mostly near zero and $180^{\circ}$ and at the same time there is a
jump in the secondary’s non-principal axis angle $\theta_{\rm NPA,s}$. The
distribution of $\theta_{ls}-\theta_{B}$ suggests that the short axis tilts
back and forth with an oscillation period that is twice the binary orbit
period. The secondary’s spin angular momentum ${\bf L}_{s}$ is often aligned
with the direction between secondary and primary. At the same time the
secondary’s short axis, when projected onto the $xy$-plane, spends time at
angle $\pm 90^{\circ}$. The secondary’s short axis periodically rocks back and
forth about the orbit normal because of rotation about its long axis. The jump
in the non-principal axis angle $\theta_{\rm NPA,s}$ suggests that it is a
resonant state, so it could last a really long time (10 to 100 times longer
than the spin down time, as did resonant non-principal axis rotation states
seen in simulations of Deimos; Quillen et al. 2020). In future, the
theoretical framework developed by Boué and Laskar (2009); Fahnestock and
Scheeres (2008); Agrusa et al. (2021) for binary systems might be applied to
better identify and create a model for this NPA rotation resonance. If sub-
resonances are identified within the 1:1 spin-orbit resonances, then perhaps a
theory accounting for the high level of NPA rotation exhibited by our
simulations could be developed.

We explored simulations similar to that shown in Figure 2 at larger orbital
semi-major axis with a larger secondary to primary mass ratio (similar to that
of the binary system Binary Near-Earth Asteroid (66391) 1999 KW4 now known as
Moshup) and saw similar behavior to that described here. We also explored
similar simulations with different axis ratios for the secondary (within 0.1
of those used here) we and also saw similar behavior.

Figure 4: A binary asteroid simulation that exhibits forced inward migration.
The initial conditions for this simulation are given by the last output of the
simulation shown in Figure 2. The panels of this figure are similar to those
of Figure 2. Time is given in orbital periods from the beginning of the
simulation. This simulation illustrates that secondary obliquity and non-
principal axis rotation can be excited during inward orbital migration. At
$t\sim 1070\ P_{B}$ the obliquity is excited, the long axis is tilted out of
the orbital plane and the resonant angle $\theta_{ls}-\theta_{B}$ ceases
spending time at 0 and $180^{\circ}$. At $t\sim 1300\ P_{B}$ a libration
resonance is excited with resonant angle $2\phi_{lib}+\theta_{kB}$. The
orbital eccentricity increases and the binary disrupts at $t\sim 1600\ P_{B}$.
This figure shows that obliquity and non-principal axis rotation can be
excited in the secondary via induced migration.

###  2.3 Simulations of inward migration

YORP and BYORP effect torques cause slow evolution so they can be considered
adiabatic. BYORP effect torque could cause the binary orbital semi-major axis
to drift either outward or inward (Ćuk and Burns, 2005). We explore whether
obliquity and non-principal axis rotation can be excited during orbital
migration. We carry out a simulation that forces the binary orbital semi-major
axis to shrink with a drift rate $\dot{a}_{B}=-10^{-5}R_{p}/t_{g}$. Migration
is forced by perturbing the center of mass velocities using the procedure
described by Beaugé et al. (2006) and which we used previously to force
outward migration of Charon when studying the obliquity evolution of Pluto and
Charon’s minor satellites (Quillen et al., 2017). Our migration simulation
uses the last output of the tidal spin down simulation (shown in Figure 2) as
its initial conditions. Thus its parameters are the same as those listed in
Table 3 except that we set the spring damping constants to 0.01 so as to
reduce tidal dissipation. The secondary begins in the 1:1 spin synchronous
state but with some NPA rotation. The migration simulation was run until the
binary disrupted at $t\sim 1600\ P_{B}$. The simulation is shown in Figure 4
which has panels similar to those of Figure 2. We restricted the time shown in
the Figure to make it easier to see variations in angles and other quantities
prior to the short period of instability that led to its breakup.

Figure 4 shows two obliquity jumps, one at $t\sim 1070\ P_{B}$ and the other
at $t\sim 1300\ P_{B}$. The first jump in obliquity occurs when the precession
angle $\theta_{ls}-\theta_{B}$ ceases to be near zero and $180^{\circ}$ and
spends more time near $90^{\circ}$ and $270^{\circ}$. At the same time the
secondary’s long axis is tilted out of the orbital plane giving oscillations
in the long axis inclination angle. The body’s long axis tilts above and below
the orbital plane with an oscillation period that is equal to 2 binary orbit
rotation periods. In a frame rotating with the secondary in its orbit, the
secondary shifts from rotating about its long axis (and tilting its short axis
w.r.t to the orbit normal) to rotating about its intermediate axis and giving
variations in the long axis inclination.

The second jump in obliquity at $t\sim 1300P_{B}$ is associated with an
increase in the libration angle and orbital eccentricity. This second jump is
also associated with a transition in the behavior or the resonant angle
$2\phi_{lib,s}+\theta_{kB}$. Prior to the second transition, this angle
oscillates and afterwards it librates about $270^{\circ}$. This angle
corresponds to the secondary rotating twice in the frame of orbital rotation
about its short axis per epicyclic period. The resonance is strong enough to
increase the orbital eccentricity until the binary disrupts. We suspect the
disruption mechanism is the same as shown in Figure 14 by Davis and Scheeres
(2020).

Because the primary spin was unaffected during the two obliquity transitions
in the migration simulation, we infer that spin/spin resonant coupling
(Seligman and Batygin, 2021) between the two bodies was not a contributing
factor. Rotation of a torque free triaxial body can be described in terms of
precession and nutation (Boué and Laskar, 2009; Fahnestock and Scheeres, 2008;
Agrusa et al., 2021). We tentatively associate the first obliquity transition
with nutation excitation. The precession and nutation frequencies can be
calculated using analytical or semi-analytical theoretical methods (Boué and
Laskar, 2009; Fahnestock and Scheeres, 2008; Agrusa et al., 2021) and the
location of resonances associated with them could in future be predicted.

The migration simulation shown in Figure 4 illustrates that if the system
migrates inward via BYORP, non-principal axis rotation can be excited due to
spin-orbit resonances that involve obliquity, libration, and non-principal
axis rotation. Thus even if tidal dissipation is strong enough to damp the
secondary’s obliquity and non-principal axis rotation, inward migration could
re-excite the secondary’s rotation state.

We also explored a simulation that forced the primary to spin down, as might
occur due to the YORP effect. The primary spin drift rate used was
$\dot{\omega}_{p}=-2\times 10^{5}t_{g}^{-2}$ and as in our migration
simulation, initial conditions were the end state of the tidal spin down
simulation shown in Figure 2. The spin state of the secondary was unaffected
until the primary spin reached the spin/spin resonance $\omega_{p}\sim
2\omega_{s}$ which caused the binary to disrupt. Only when the primary spin
was low, did we see evidence for spin-spin resonant coupling. In this setting
spin/spin resonant coupling (Seligman and Batygin, 2021) could be important.

We explored a simulation similar to that shown in Figure 4 except the
secondary migrated outward (instead of inward) from $a_{B}/R_{p}\approx 2$ to
4 at a drift rate $\dot{a}_{B}=5.0\times 10^{-6}$. This simulation showed a
similar obliquity jump at $a_{B}/R_{p}\sim 2.1$ as the first transition in
Figure 4. This implies that the resonance that excited the long axis
inclination is sensitive to the binary period. Otherwise this simulation was
dull. After that transition, we did not see obliquity or non-principal axis
excitation. Perhaps simulations that migrate outward more slowly
(adiabatically) would reveal weaker resonant behavior.

##  3 Implications for BYORP drift

In section 2.2 our tidal spin down simulation showed that secondary obliquity
and non-principal axis excitation can be long-lived within the 1:1 spin-orbit
resonance. In section 2.3 our migration simulation showed that obliquity, long
axis inclination and libration can be excited if the binary semi-major axis
drifts inward due to BYORP. In this section we explore effect of non-principal
axis rotation on the BYORP drift rate for a secondary that is within the 1:1
spin-orbit resonance.

Prior calculations of the BYORP affect (Ćuk and Burns, 2005; McMahon and
Scheeres, 2010; Steinberg and Sari, 2011; Scheirich et al., 2021) assume that
the secondary is tidally locked (in the 1:1 spin synchronous state), undergoes
principal axis rotation, is at zero obliquity and lacks free libration. To
estimate the role of non-principal axis rotation in affecting the BYORP
induced drift rate we compute the BYORP torque on the binary orbit using
different secondary shape models that are comprised of triangular surface
meshes.

The shapes of the secondaries in most binary systems are not well
characterized. An exception is with the Binary Near-Earth Asteroid (66391)
1999 KW4, with primary known as Moshup and secondary known as Squannit.
Squannit’s shape model was derived by Ostro et al. (2006) and a BYORP
generated drift rate was computed with this shape model by Scheirich et al.
(2021). We use the same Squannit shape model to compare our BYORP drift rate
estimates to that computed by Scheirich et al. (2021).

To estimate the torque from BYORP we neglect surface thermal inertia,
following Rubincam (2000), so that thermal radiation is re-emitted with no
time lag. The reflected and thermally radiated components are assumed to be
Lambertian so they are emitted with flux parallel to the local surface normal.
We ignore heat conduction and neglect shadows on the day-lit side. We ignore
variations of albedo or emissivity on the surface. Secondary intersections of
re-radiated energy are ignored.

The asteroid surface is described with a closed triangular oriented vertex-
face mesh, similar to the calculations by Scheeres (2007); McMahon and
Scheeres (2010); Steinberg and Sari (2011). We index each triangular facet
with an integer $i$. The $i$-th facet has area $S_{i}$, and unit normal
$\hat{\bf n}_{i}$. The radiation force on the $i$-th triangular facet is
computed as

${\bf F}_{i}=\begin{cases}-\frac{2}{3}\frac{F_{\odot}}{c}{S_{i}}(\hat{\bf
n}_{i}\cdot\hat{\bf s}_{\odot})\hat{\bf n}_{i}&\mathrm{if\ }\hat{\bf
n}_{i}\cdot\hat{\bf s}_{\odot}>0\\\ 0&\mathrm{otherwise}\end{cases},$ (16)

where $F_{\odot}$ is the solar radiation flux and $c$ is the speed of light.
The direction of the Sun is given by unit vector $\hat{\bf s}_{\odot}$. The
factor of 2/3 is consistent with Lambertian reflectance and emission. Facets
on the day side have $\hat{\bf n}_{i}\cdot\hat{\bf s}_{\odot}>0$ and we assume
that only these contribute to the radiation force. In terms of the radial
distance $r_{h}$ from the binary to the Sun,
$F_{\odot}\equiv\frac{P_{\odot}}{r_{h}^{2}}$ where radiation pressure constant

$P_{\odot}\equiv\frac{L_{\odot}}{4\pi c}=10^{17}\ {\rm kg\ m\ s}^{-2}$ (17)

and where $L_{\odot}$ is the solar luminosity.

The torque affecting the binary orbit from a single facet is

${\boldsymbol{\tau}}_{i,B}=\begin{cases}-\frac{2}{3}\frac{F_{\odot}}{c}{S_{i}}(\hat{\bf
n}_{i}\cdot\hat{\bf s}_{\odot})({\bf a}_{B}\times\hat{\bf n}_{i})&\mbox{if
}\hat{\bf n}_{i}\cdot\hat{\bf s}_{\odot}>0\\\ 0&\mbox{otherwise}\end{cases},$
(18)

where ${\bf a}_{B}$ is the secondary’s radial vector from the binary center of
mass. We neglect the distance of the facet centroids from the secondary center
of mass when computing the radiative torques.

The torque affecting the binary orbit is the sum of the torques from each
facet. This torque is averaged over different binary positions in its orbit
around the Sun, affecting the direction of ${\bf s}_{\odot}$, and over the
different positions in the binary mutual orbit and spin rotation of the
secondary,

${\boldsymbol{\tau}}_{BY}=\frac{1}{T}\int_{0}^{T}dt\
\sum_{i}{\boldsymbol{\tau}}_{i,B}.$ (19)

To represent this average, time $T$ significantly exceeds binary and
heliocentric orbit periods. In practice we sum over different orbital phases
to compute the average. The average can also be done by integrating over the
mean anomaly $M_{h}$ of the binary in its heliocentric orbit and the mean
anomaly of the binary orbit $M_{B}$

${\boldsymbol{\tau}}_{BY}=\int_{0}^{2\pi}\int_{0}^{2\pi}dM_{h}\ dM_{B}\
\sum_{i}{\boldsymbol{\tau}}_{i,B}.$ (20)

If $\hat{\bf l}_{B}$ is the binary orbit normal then the mutual semi-major
axis drift rate is

$\dot{a}_{B}=\frac{2{\boldsymbol{\tau}}_{BY}\cdot\hat{\bf
l}_{B}}{M_{s}n_{B}a_{B}}.$ (21)

BYORP is often described in terms of a dimensionless number (e.g., McMahon and
Scheeres 2010). By summing over different orientations and face model facets,
we calculate the dimensionless number

$B_{BY}\equiv({\boldsymbol{\tau}}_{BY}\cdot\hat{\bf
l}_{B})\frac{c}{F_{\odot}R_{s}^{2}a_{B}}.$ (22)

This is equivalent to letting $a_{B}=R_{s}=F_{\odot}/c=1$. For a binary in a
non-circular heliocentric orbit, averaging of the solar flux would give
$\bar{F}_{\odot}=\frac{P_{\odot}}{a_{h}^{2}\sqrt{1-e_{h}^{2}}}$ where $a_{h}$
and $e_{h}$ are the binary orbit’s heliocentric semi-major axis and
eccentricity. With our dimensionless number, $B_{BY}$, equations 21 and 22
give

$\displaystyle\dot{a}_{B}$
$\displaystyle=\frac{2P_{\odot}}{a_{h}^{2}\sqrt{1-e_{h}^{2}}}\frac{R_{s}^{2}}{M_{s}n_{B}}g_{BY}.$
(23)

We compare our expression to equation 3 by Scheirich et al. (2021)

$\displaystyle\dot{a}_{B}$
$\displaystyle=\frac{2P_{\odot}}{a_{h}^{2}\sqrt{1-e_{h}^{2}}}\frac{R_{mean}^{2}}{M_{s}n_{B}}B.$
(24)

Comparison between this equation and equation 23 shows that our dimensionless
$B_{BY}$ parameter of equation 22 is equivalent to the dimensionless BYORP $B$
parameter introduced by McMahon and Scheeres (2010). Our parameter $B_{BY}=\pi
f_{BY}$ where $f_{BY}$ is the dimensionless parameter computed by Steinberg
and Sari (2011).

###  3.1 Shape models

To explore how BYORP depends on surface shape we use 3 shape models that are
triangular surface meshes. The three shape models are an icosphere, a coarse
version of the Squannit shape model by Ostro et al. (2006) and a perturbed
version of the same Squannit shape model. With the Squannit shape model by
Ostro et al. (2006), Scheirich et al. (2021) computed a BYORP B-coefficient of
$B=2.082\times 10^{-2}$ for the Moshup/Squannit binary system. To test our
BYORP B-coefficient computation, we compare our computed B-coefficient to that
computed by Scheirich et al. (2021).

The icosphere is an icosahedron that is subdivided twice. Each triangular face
of the icosahedron is divided into 4 triangular faces and then each of them is
divided into 4 triangular faces giving a total of 320 triangular faces. The
vertex positions are then corrected so they have distance 1 from the origin.

The Squannit shape model by Ostro et al. (2006) contains 2292 faces is shown
with $xy$ and $yz$ plane views in the top two panels of Figure 5. The red
bounding squares in this Figure have a size of 3 in units of the radius of a
volume-equivalent sphere. To speed up our BYORP calculation we compute the
BYORP drift rate for a similar but coarser triangular surface mesh with fewer
faces. We reduce the face number by merging short edges that are less than a
specified length. We adjusted the specified minimum edge length to give a
triangular surface mesh with about 300 faces. The coarser version of the
Squannit shape model with 302 faces is shown in the middle panels of 5. We
checked that the results of our BYORP B-coefficient calculations are not
sensitive to the number of faces used in our model (coefficients were less
than a few percent different than those computed with a surface mesh that has
twice the number of faces).

To create a rougher version of the Squannit shape model we perturbed each
vertex of the coarse Squannit 302 face shape model by adding random values to
the $x,y$ and $z$ coordinates of each vertex. The random values are taken from
a uniform distribution spanning $[-\Delta,\Delta]$ where $\Delta=0.05$ in
units of the radius of the volume-equivalent sphere. The perturbed 302 face
Squannit shape model is shown in bottom two panels of Figure 5. The
perturbations to the vertices make it rougher than the original and coarse
Squannit shape models.

Figure 5: Triangular surface mesh shape model for Squannit by Ostro et al.
(2006) is shown in the top two panels. The $xy$ plane is viewed on the left
and the $xz$ plane is viewed on the right. This shape model has 2292 faces.
The middle two panels shows a coarser version of the same shape model with 302
faces. The bottom two panels show a model that is similar to the coarse one,
but each vertex position has been perturbed in a random direction to increase
the surface roughness. The red bounding squares have a length of 3 in units of
the radius of the volume equivalent sphere.  Figure 6: BYORP B-coefficients
computed using two Squannit secondary shape models shown in the bottom four
panels of Figure 5, as a function of the maximum tilt angle of the body’s
short axis from the orbit normal. The body tilt is randomly chosen at each
binary orbit orientation angle used in the torque averaging computation. The
black circle is the B-coefficient computed by Scheirich et al. (2021) based on
the Squannit model by Ostro et al. (2006). The green circles are computed for
an icosphere. The red squares show the coarse Squannit shape model with 302
faces. The brown crosses show a perturbed or rough Squannit shape model, also
with 302 faces. The BYORP B-coefficient reverses sign when the short angle
tilts by more than about $60^{\circ}$.

###  3.2 BYORP B-coefficients computed with randomly chosen short axis tilt
angles

Using the three shape models described in section 3.1, we numerically compute
BYORP B-coefficients in two settings. In our first set of calculations, we
compute B-coefficients assuming a circular binary orbit and choosing
orientations for the secondary at different positions in the orbit. In our
second set of calculations (section 3.3), we use the secondary orientations
measured directly from our simulations to compute the B-coefficients.

In our mass-spring model simulations described in section 2, non-principal
axis rotation of the secondary was present at low libration angle and was
characterized by rotation about the secondary’s long axis. The secondary’s
long axis remains aligned with the direction to the primary, corresponding to
$\phi_{lib,s}\approx 0$, however the short axis can tilt away from the orbit
normal. We described this rotation as the short axis tilt and it is measured
from the angle between the body short axis and the orbit normal. We assume
that the secondary angular rotation rate $\omega_{s}=n_{B}$ is equal to the
binary orbit mean motion. This condition places the secondary in the 1:1 spin
synchronous resonance. We orient the secondary so that its long axis remains
in the orbital plane and points directly toward the primary so that the
libration angle remains at zero. At each binary orbit orientation, the
secondary’s short axis orientation is randomly chosen using a uniform
distribution so that the angle between the short axis and the orbit normal
ranges from 0 to a maximum of $\theta_{\rm maxtilt}$. Specifically, after
rotating the secondary so that its long axis is pointed toward the primary,
the secondary is rotated about its long axis by an angle chosen from a uniform
distribution within $\\{-\theta_{\rm maxtilt},\theta_{\rm maxtilt}\\}$.

We assume a zero inclination circular orbit for the binary (about the Sun) and
a circular mutual binary orbit. To numerically calculate the BYORP
B-coefficient we first average over the Sun’s possible directions using 36
equally spaced solar angles. Then we average over the binary orbit using 144
equally spaced binary orbit orientations. In each of these 144 possible
orientations the short axis tilt is randomly chosen.

The BYORP B-coefficients computed using the three shape models shown in Figure
5 and with randomly chosen short axis tilt angles are shown in Figure 6 as a
function $\theta_{\rm maxtilt}$. The green circles show that the B-coefficient
computed for the icosphere is zero, as expected. The coarse Squannit face
model, shown with red squares, gives a BYORP B-coefficient of $B_{BY}=0.02$
for $\theta_{\rm maxtilt}=0$ and this agrees with the value computed by
Scheirich et al. (2021) for Squannit that is shown with a black circle. For
both Squannit and rough Squannit shape models, the BYORP B-coefficient drops
to zero at a maximum tilt angle of about $\theta_{\rm maxtilt}\sim
60^{\circ}$. This gives a quantitative estimate for how much non-principal
axis rotation within the 1:1 spin-orbit resonance is required to reduce BYORP
drift.

The rough or perturbed Squannit shape model gives lower values of the BYORP
B-coefficient than the Squannit shape model. This is opposite to the expected
trend. Usually rougher surface models have higher YORP and BYORP coefficients
(Steinberg and Sari, 2011). By computing BYORP B-coefficients for perturbed
ellipsoids we confirmed that rougher surfaces tend to have larger
B-coefficients. The specific shape of the Squannit shape model must be
responsible for its high B-coefficient. The ratio of B-coefficients for
Squannit and perturbed Squannit shape models we computed at $\theta_{\rm
maxtilt}=0$ is lower than 2. In the Moshup/Squannit system, Scheirich et al.
(2021) measured a drift rate about 8 times lower than predicted. A rougher
surface for Squannit does not give a large enough difference in the
B-coefficient to fully account for its low measured B-coefficient (assuming
principal axis rotation). In contrast, the BYORP B-coefficient can be
negligible if there is non-principal axis rotation with $\theta_{\rm
maxtillt}\sim 60^{\circ}$.

Figure 7: BYORP-B coefficient is computed using secondary and binary
orientations from the tidal spin-down simulation shown in Figure 2 and is
shown in the top panel as a function of time in orbital periods. The
B-coefficients were computed using our three shape models (see Figure 5). The
red squares are the coarse Squannit shape model, the brown pluses are the
perturbed (or rougher) Squannit shape model. The green circles show an
icosphere shape model. The icosphere has a B-coefficient of zero, as expected.
The secondary’s non-principal axis angle $\theta_{\rm NPA,s}$ is shown in the
bottom panel. The BYORP B-coefficients for the Squannit shape models
frequently reverse in sign. Only after the NPA angle has dropped below
$45^{\circ}$ near the end of the simulation does the BYORP B-coefficient
remain positive. The rougher Squannit shape model has B-coefficients similar
to but lower amplitude than the coarse Squannit shape model. Both Squannit
shape modes have BYORP B-coefficients lower than predicted for a tidally
locked body undergoing principal axis rotation, at zero obliquity and with no
free libration. The conventionally predicted B-coefficients are 0.02 for the
coarse Squannit model and 0.014 for the rough Squannit model. This figure
shows that long-lived non-principal axis rotation within the 1:1 spin-orbit
resonance can cause a reduction in the amplitude of the BYORP B-coefficient
and reversals in the direction of BYORP effect drift.  Figure 8: Similar to
Figure 7 except the BYORP B coefficient is computed using the migration
simulation shown in Figure 4. As the non-principal axis rotation angle drops,
the BYORP B-coefficients for the Squannit shape models increase.

###  3.3 BYORP B-coefficients computed using the mass-spring model numerical
simulations

We find the secondary orientation at each simulation output in our mass-spring
model simulations. We use the vector between secondary and primary, binary
orbit normal and body orientation at each simulation output to compute the
BYORP B-coefficient for three shape models, two of which are shown in Figure
5.

The secondary orientation (computed with a quaternion from a simulation output
as described in section 2.1) was used to rotate the shape model vertices. The
rotated shape model vertices, face normals and face areas were then used to
compute the radiative torque on each facet using binary orientation vector and
orbit normal from the concurrent simulation output. At each simulation output,
18 solar positions were used to compute an average BYORP torque. 100
consecutive simulation outputs, spanning about 6 binary orbit periods, were
then combined for an average of the BYORP induced torque at a particular
moment during the simulation. We note that the axis ratios of the Squannit
shape model ($b_{s}/a_{s}=0.762$, $c_{s}/a_{s}=0.576$; Ostro et al. 2006) are
similar to but not exactly the same as those of the secondary in the
simulation (0.773, 0.698).

To check our numerical routines we ran a short simulation of a secondary in
the spin synchronous state, zero obliquity and in a principal axis rotation
rate. We measured the BYORP B-coefficient in this short simulation. The
measured B-coefficient values for each shape model were consistent with those
computed at $\theta_{\rm maxtilt}=0$, and shown in Figure 6, for the same
shape models, as expected. The B-coefficient computed for the coarse Squannit
shape model again agreed with that computed by Scheirich et al. (2021),
$B_{BY}\approx 0.02$.

Figure 7 shows the BYORP B-coefficients computed for our three shape models
from the tidal spin-down simulation shown in Figure 2 as a function of time
and Figure 8 similarly shows the BYORP B-coefficients but computed for the
migration simulation shown in Figure 4. In these figures, green circles refer
to the icosphere shape model, red squares refer to the 302 face coarse
Squannit shape model and brown pluses refer to the rougher Squannit shape
model. In both figures the top panels show the B-coefficients and the bottom
panels show the non-principal axis angle $\theta_{\rm NPA,s}$.

Figure 7 shows that the BYORP B-coefficients are lower than those computed for
principal axis rotation state ($B_{BY}\approx 0.020$ for the Squannit shape
model and $B_{BY}\approx 0.014$ for the perturbed or rough Squannit shape
model) throughout the tidal-spin down simulation. The sign of the
B-coefficient frequently reverses which implies that the cumulative
B-coefficient averages to zero. Not until the end of the simulation do the
reversals cease and the BYORP B-coefficient remains above zero. At the end of
the simulation the average B-coefficient is about 4 times lower than that in a
principal axis rotation state. The simulation ran for 15000 orbital periods
but the 1:1 spin-orbit resonance was reached after only about 1000 orbital
periods. The BYORP torque was low for a duration an order of magnitude longer
than the spin down time. We infer that long-lived non-principal axis rotation
within the 1:1 spin-orbit resonance may prolong the lifetime of the binary by
reducing BYORP effect orbital drift.

In Figure 7 the BYORP B-coefficient for the Squannit shape models showed
reversals even for non-principal axis angles $\theta_{\rm NPA,s}$ as low as
$45^{\circ}$. This is smaller than the maximum short axis tilt angle that gave
a B-coefficient of zero in Figure 6. In the computations described in section
3.2, we assumed a uniform distribution for the short axis tilt and neglected
variations in libration and long axis inclination. This may not be a good
approximation for the orientation angle distributions exhibited by the
secondary in the numerical simulations. A comparison between the red squares
and brown pluses in Figure 7 shows that the rough Squannit model has BYORP
B-coefficients similar to but lower in amplitude than that of the coarse
Squannit model, as we found previously in section 3.2.

The migration simulation shown in Figures 4 and 8 has initial condition
starting where the tidal-spin down simulation ended. In Figure 8, as the NPA
angle decreases, the amplitude of the BYORP B-coefficient slowly increases
until it nears the value expected for a tidally locked body undergoing
principal axis rotation (0.02 for the Squannit shape model). This simulation
exhibited jumps in obliquity and long axis inclination angle. However the
BYORP B-coefficients slowly increase throughout the simulation (until right
before the binary is disrupted) and more closely follow the NPA angle than the
obliquity or long axis tilt variations. We suspect that the time averaged
BYORP B-coefficient is more sensitive to the NPA angle variations than to the
libration amplitude, or the long axis tilt and obliquity variations.

We have presented an estimate of the BYORP torque as a function of time using
two spin-dynamics simulations and two shape models. Future studies can improve
upon and extend this study by carrying out longer duration simulations (more
orbits) with more slower tidal evolution or migration and with better resolved
secondary shapes and internal structure. Future BYORP B-coefficients computed
using spin dynamics simulations could explore how the sensitivity of the BYORP
torque to the extent of non-principal axis rotation depends on body shape,
internal structure (e.g., homogeneity) and surface properties such as albedo,
reflectively, emissivity and thermal inertia.

##  4 Summary and Discussion

Using our viscoelastic mass/spring model N-body code, we have carried out long
timescale ($10^{4}$ binary orbital period or about 15 year long) simulations
of a tidally evolving asteroid binary system similar to the NEA (65803)
Didymos binary system. We begin a simulation with secondary outside of tidal
lock, undergoing principal axis rotation and at low obliquity. Due to attitude
instability when crossing into the 1:1 spin-orbit resonance, the secondary’s
obliquity is excited and it enters a non-principal axis rotation state
(Wisdom, 1987). We find that the obliquity and non-principal axis rotation
damping time exceeds the time to enter the spin synchronous state by about an
order of magnitude, so secondary obliquity and non-principal axis rotation can
be long-lived in NEA binary secondaries, perhaps even lasting millions of
years.

While in the 1:1 spin-orbit resonance and in a non-principal axis rotation
state, the secondary’s long axis librates about the direction of the primary
but the secondary’s short axis can exhibit wide excursions away from the orbit
normal. In such a state, it may be difficult to identify non-principal axis
rotation from a light curve, and yet, the secondary could be correctly
identified as residing in the 1:1 spin-orbit resonance. We propose that some
binary asteroids that have been classified as having synchronous secondaries
from their light curves may reside in the 1:1 spin-orbit resonance but may
simultaneously reside in a complex non-principal axis rotation state.

Our numerical simulations confirm the work by (Naidu and Margot, 2015) who
emphasized that non-principal axis rotation may be slowly damped. Long-lived
non-principal axis rotation was also previously seen in some simulations of a
tidally evolving Deimos (Quillen et al., 2020).

If binary secondaries are observed to undergo principal axis rotation, have
near-zero obliquity and low amplitudes of free libration, then they would
require even higher levels of tidal dissipation than previously estimated to
reach their tidally locked states. Lower required values of $\mu Q$ (shear
modulus times tidal dissipation parameter) would exacerbate the tension
between viscoelastic material properties estimated for single and binary
asteroids. In contrast, if the secondaries are observed to be in non-principal
axis rotation states, then higher values of $\mu Q$ would be possible and the
discrepancy between estimated single and binary asteroids material properties
may be ameliorated.

To mimic the effect of BYORP induced migration, we carried out a numerical
simulation of a binary system initially in a low obliquity, low principal axis
rotation and spin-synchronous rotation state. Spin orbit resonances, crossed
during migration, can excite obliquity, libration and non-principal axis
rotation. Even if tidal dissipation is strong enough to damp the secondary’s
obliquity and non-principal axis rotation, BYORP induced migration could re-
excite the secondary’s rotation state. We saw resonant excitation of obliquity
and non-principal axis rotation but only during inward migration. Impacts on
the secondary or surface movements on it might also re-excite its rotation
state.

Previous BYORP calculations have necessarily assumed zero secondary obliquity,
principal axis rotation and low libration amplitude within the 1:1 spin-orbit
resonance. An expectation is that when the secondary tumbles, the cumulative
BYORP torque will average to zero (Ćuk and Burns, 2005; Naidu and Margot,
2015). However, our simulations show an intermediate rotation state where the
secondary remains within the 1:1 spin-orbit resonance and simultaneously
exhibits chaotic non-principal axis variations. Using a shape model (Ostro et
al., 2006) for Squannit, the secondary of the Moshup system (NEA (66391) 1999
KW4), we estimate the BYORP B-coefficient. We find that the BYORP torque is
dependent upon the extent of minor axis tilt within the 1:1 spin synchronous
state. There is a range of possible BYORP torque values that is dependent upon
the extent of non-principal axis rotation.

We use secondary body orientations measured in our mass-spring model
simulations and Squannit shape models to compute the BYORP B-coefficient. The
B-coefficient is reduced compared to that computed for the tidally locked,
zero obliquity and principal axis rotation state. The B-coefficient shows
reversals in sign for an extended period of time after the spin-synchronous
state is entered. Only when the non-principal axis angle drops below about
$45^{\circ}$ does the B-coefficient and associated BYORP effect torque remain
unidirectional. The lifetime of the binary would be extended due to non-
principal axis rotation. However, we find that non-principal axis rotation
excitation due to inward migration may not be strong enough to entirely
prevent migration by the BYORP effect. Once the secondary non-principal axis
angle drops below about $45^{\circ}$, and unless non-principal axis rotation
in the secondary is re-excited by impacts or mass-movements, the BYORP induced
orbital migration rate can increase until it is as rapid as estimated
previously for the principal axis rotation state (Ćuk and Burns, 2005).

Non-principal axis rotation of Squannit, the secondary in the Moshup system,
could help account for its low BYORP measured drift rate, which was about 7
times lower than that predicted from the secondary shape model (Scheirich et
al., 2021). To a lesser extent, a shape model with a rougher surface would
also give a lower B-coefficient. As discussed by Scheirich et al. (2021), the
discrepancy between predicted and measured BYORP drift rate may simply be due
to uncertainties in the secondary’s topography.

The DART and Hera missions to Didymos, on-going photometric surveys and radar
studies can confirm or rule out the possibility that binary asteroid
secondaries are in complex rotation states. If binary secondaries are found in
complex rotation states, then tidal dissipation need not be as high as that
required to balance the conventionally estimated BYORP torque, as previously
proposed by Jacobson and Scheeres (2011). As the rotation state affects the
BYORP induced orbital drift rate, binary population evolution models (Jacobson
et al., 2016) could be modified to take into account the role and evolution of
the secondary rotation state. Long duration spin dynamics simulations could be
modified to take into account YORP and BYORP effect induced torques. Meteroid
impacts and mass movements can perturb the spin state and body shape, and
modify surface properties such as albedo, reflectively, emissivity and thermal
inertia, which affect the radiative torque. As a consequence, binary asteroid
population evolution models could take into account these additional
processes.

Acknowledgements

This material is based upon work supported in part by NASA grants
80NSSC21K0143 and 80NSSC17K0771, and National Science Foundation Grant No.
PHY-1757062.

We thank Patrick Michel for correspondence that helped us improve the
manuscript.

A repository for code developed and used in this paper is at
https://github.com/aquillen/asteroid_binary.

## References

  * Agrusa et al. (2021) Agrusa, H. F., Gkolias, I., Tsiganus, K., Richardson, D. C., Meyer, A. J., Scheeres, D. J., Ćuk, M., Jacobson, S. A., Michel, P., Karatekin, Ö.., Cheng, A. F., Hirabayashi, M., Zhang, Y., Fahnestock, E. G., Davis, A. B., 2021\. The excited spin state of Dimorphos resulting from the DART impact. Icarus.
  * Agrusa et al. (2020) Agrusa, H. F., Richardson, D. C., Davis, A. B., Fahnestock, E., Hirabayashi, M., Chabot, N. L., Cheng, A. F., Rivkin, A. S., Michel, P., the DART Dynamics Working Group, 2020. A benchmarking and sensitivity study of the full two-body gravitational dynamics of the DART mission target, binary asteroid 65803 Didymos. Icarus 349, 113849.
  * Beaugé et al. (2006) Beaugé, C., Michtchenko, T. A., Ferraz-Mello, S., 2006. Planetary migration and extrasolar planets in the 2/1 mean-motion resonance. Monthly Notices of the Royal Astronomical Society 365, 1160–1170.
  * Berthold and Horn (1986) Berthold, K., Horn, P., 1986. Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society of America 4, 629–642.
  * Boué and Laskar (2009) Boué, G., Laskar, J., 2009. Spin axis evolution of two interacting bodies. Icarus 201, 750–767.
  * Burns and Safronov (1973) Burns, J. A., Safronov, V. S., 1973. Asteroid nutation angles. Monthly Notices Royal Astronomical Society 165, 403–411.
  * Cheng et al. (2016) Cheng, A., Michel, P., Jutzi, M., Rivkin, A., Stickle, A., Barnouin, O., Ernst, C., Atchison, J., Pravec, P., Richardson, D., 2016. Asteroid impact & deflection assessment mission: Kinetic impactor. Planetary and Space Science 121, 27–35.
URL https://www.sciencedirect.com/science/article/pii/S0032063315300337

  * Chesley et al. (2020) Chesley, S. R., French, A. S., Davis, A. B., Jacobson, R. A., Brozović, M., Farnocchia, D., Selznick, S., Liounis, A. J., Hergenrother, C. W., Moreau, M. C., Pelgrift, J., Lessac-Chenen, E., Molaro, J. L., Park, R. S., Rozitis, B., Scheeres, D. J., Takahashi, Y., Vokrouhlický, D., Wolner, C. W. V., Adam, C., Bos, B. J., Christensen, E. J., Emery, J. P., Leonard, J. M., McMahon, J. W., Nolan, M. C., Shelly, F. C., Lauretta, D. S., aug 2020\. Trajectory estimation for particles observed in the vicinity of (101955) bennu. Journal of Geophysical Research: Planets 125 (9).
URL https://doi.org/10.1029%2F2019je006363

  * Coutsias et al. (2004) Coutsias, E. A., Seok, C., Dill., K. A., 2004. Using quaternions to calculate RMSD. Journal of Computational Chemistry 25, 1849–1857.
  * Ćuk and Burns (2005) Ćuk, M., Burns, J. A., 2005. Effects of thermal radiation on the dynamics of binary NEAs. Icarus 176, 418–431.
  * Ćuk et al. (2020) Ćuk, M., Jacobson, S., Walsh, K., 2020. ”Barrel Instability” for elongated secondaries in binary asteroids. In: AAS/Division of Dynamical Astronomy Meeting. p. 103.04.
  * Ćuk and Nesvorný (2010) Ćuk, M., Nesvorný, D., 2010. Orbital evolution of small binary asteroids. Icarus 207, 732–743.
  * Davis and Scheeres (2020) Davis, A. B., Scheeres, D. J., 2020. High-fidelity modeling of rotationally fissioned asteroids. The Planetary Science Journal 1, 25–42.
  * Dotto et al. (2021) Dotto, E., Corte, V. D., Amoroso, M., Bertini, I., Brucato, J. R., Capannolo, A., Cotugno, B., Cremonese, G., Tana, V. D., Gai, I., Ieva, S., Impresario, G., Ivanovski, S., Lavagna, M., Lucchetti, A., Epifani, E. M., Meneghin, A., Miglioretti, F., Modenini, D., Pajola, M., Palumbo, P., Perna, D., s. Pirrotta, Poggiali, G., Rossi, A., Simioni, E., Simonetti, S., Tortora, P., Zannoni, M., Zanotti, G., Zinzi, A., Cheng, A., Rivkin, A., Adams, E., Reynolds, E., Fretz, K., 2021. LICIACube - the light Italian Cubesat for imaging of asteroids in support of the NASA DART mission towards asteroid (65803) Didymos. Planetary and Space Science 199 (7590), 105185.
  * Efroimsky (2015) Efroimsky, M., 2015. Tidal evolution of asteroidal binaries. Ruled by viscosity. Ignorant of rigidity. Astronomical Journal 150, 98–110.
  * Efroimsky and Makarov (2013) Efroimsky, M., Makarov, V. V., Feb 2013. Tidal friction and tidal lagging. Applicability limitations of a popular formula for the tidal torque. Astrophysical Journal 764 (1), 26.
  * Fahnestock and Scheeres (2008) Fahnestock, E., Scheeres, D., 2008. Simulation and analysis of the dynamics of binary near-earth asteroid (66391) 1999 KW4. Icarus 194, 410–435.
  * Frouard et al. (2016) Frouard, J., Quillen, A. C., Efroimsky, M., Giannella, D., 2016. Numerical simulation of tidal evolution of a viscoelastic body modelled with a mass-spring network. Monthly Notices of the Royal Astronomical Society 458, 2890–2901.
  * Genova et al. (2016) Genova, A., Goossens, S., Lemoine, F. G., Mazarico, E., Neumann, G. A., Smith, D. E., Zuber, M. T., 2016. Seasonal and static gravity field of Mars from MGS, Mars Odyssey and MRO radio science. Icarus 272, 228–245.
  * Gladman et al. (2000) Gladman, B., Michel, P., Froeschlé, C., 2000. The near-Earth object population. Icarus 146, 176–189.
  * Gladman et al. (1996) Gladman, B., Quinn, D. D., Nicholson, P., Rand, P., 1996. Synchronous locking of tidally evolving satellites. Icarus 122, 166–192.
  * Goldreich (1963) Goldreich, P., 1963. On the eccentricty of satellite orbits in the solar system. Monthly Notices of the Royal Astronomical Society 126, 257–268.
  * Goldreich and Sari (2009) Goldreich, P., Sari, R., 2009. Tidal evolution of rubble piles. Astrophysical Journal 691, 54–60.
  * Harris (1994) Harris, A. W., 1994. Tumbling asteroids. Icarus 107, 209–211.
  * Hirabayashi et al. (2015) Hirabayashi, M., Sánchez, D. P., Scheeres, D. J., Jul 2015. Internal Structure of Asteroids Having Surface Shedding Due to Rotational Instability. Astrophysical Journal 808 (1), 63.
  * Jacobson and Scheeres (2011) Jacobson, S., Scheeres, D., 2011. Long-term stable equilibria for synchronous binary asteroids. Astrophysical Journal Letters 736, L19–L23.
  * Jacobson et al. (2016) Jacobson, S. A., Marzari, F., Rossi, A., Scheeres, D. J., 2016. Matching asteroid population characteristics with a model constructed from the YORP-induced rotational fission hypothesis. Icarus 277, 381–394.
  * Jawin et al. (2020) Jawin, E. R., Walsh, K. J., Barnouin, O. S., McCoy, T. J., Ballouz, R.-L., DellaGiustina, D. N., Connolly, H. C., Marshall, J., Beddingfield, C., Nolan, M. C., Molaro, J. L., Bennett, C. A., Scheeres, D. J., Daly, M. G., Asad, M. A., Daly, R. T., Bierhaus, E. B., Susorney, H. C. M., Kaplan, H. H., Enos, H. L., Lauretta, D. S., sep 2020. Global patterns of recent mass movement on asteroid (101955) Bennu. Journal of Geophysical Research: Planets 125 (9).
URL https://doi.org/10.1029%2F2020je006475

  * Kabsch (1976) Kabsch, W., 1976. A solution for the best rotation to relate two sets of vectors. Acta Crystallographica A32, 922–923.
  * Kaula (1964) Kaula, M., 1964. Tidal dissipation by solid friction and the resulting orbital evolution. Reviews of Geophysics 2, 661–684.
  * Kaula (1966) Kaula, W. M., 1966. Theory of Satellite Geodesy, Applications of Satellites to Geodesy. Blaisdell Publishing Company, Waltham Mass.
  * Lauretta et al. (2019) Lauretta, D. S., Hergenrother, C. W., Chesley, S. R., Leonard, J. M., Pelgrift, J. Y., Adam, C. D., Asad, M. A., Antreasian, P. G., Ballouz, R.-L., Becker, K. J., Bennett, C. A., Bos, B. J., Bottke, W. F., Brozović, M., Campins, H., Connolly, H. C., Daly, M. G., Davis, A. B., de León, J., DellaGiustina, D. N., d’Aubigny, C. Y. D., Dworkin, J. P., Emery, J. P., Farnocchia, D., Glavin, D. P., Golish, D. R., Hartzell, C. M., Jacobson, R. A., Jawin, E. R., Jenniskens, P., Kidd, J. N., Lessac-Chenen, E. J., Li, J.-Y., Libourel, G., Licandro, J., Liounis, A. J., Maleszewski, C. K., Manzoni, C., May, B., McCarthy, L. K., McMahon, J. W., Michel, P., Molaro, J. L., Moreau, M. C., Nelson, D. S., Owen, W. M., Rizk, B., Roper, H. L., Rozitis, B., Sahr, E. M., Scheeres, D. J., Seabrook, J. A., Selznick, S. H., Takahashi, Y., Thuillet, F., Tricarico, P., Vokrouhlický, D., Wolner, C. W. V., dec 2019. Episodes of particle ejection from the surface of the active asteroid (101955) Bennu. Science 366 (6470), eaay3544.
URL https://doi.org/10.1126%2Fscience.aay3544

  * Margot et al. (2002) Margot, J. L., Nolan, M. C., Benner, L. A. M., Ostro, S. J., Jurgens, R. F., Giorgini, J. D., Slade, M. A., Campbell, D. B., 2002. Binary asteroids in the near-Earth object population. Science 296 (5572), 1445–1448.
  * McMahon and Scheeres (2010) McMahon, J., Scheeres, D., 2010. Detailed prediction for the BYORP effect on binary near-Earth asteroid (66391) 1999 KW4 and implications for the binary population. Icarus 209, 494–509.
  * Michel and et al (2021) Michel, P., et al, 2021. The ESA Hera mission: Detailed characterisation of the DART impact outcome and of the asteroid binary (65803) Didymos, in preparation.
  * Murray and Dermott (1999) Murray, C. D., Dermott, S. F., 1999. Solar System Dynamics. Cambridge University Press.
  * Naidu et al. (2020) Naidu, S., Benner, L., Brozovic, M., Nolan, M., Ostro, S., Margot, J., Giorgini, J., Hirabayashi, T., Scheeres, D., Pravec, P., Scheirich, P., Magri, C., Jao, J., 2020. Radar observations and a physical model of binary near-earth asteroid 65803 Didymos, target of the DART mission. Icarus 348, 113777.
URL https://www.sciencedirect.com/science/article/pii/S0019103520301640

  * Naidu and Margot (2015) Naidu, S. P., Margot, J.-L., 2015. Near-Earth asteroid satellite spins under spin-orbit coupling. The Astronomical Journal 149 (2), 80.
  * Nimmo and Matsuyama (2019) Nimmo, F., Matsuyama, I., 2019. Tidal dissipation in rubble-pile asteroids. Icarus 321, 715–721.
  * Ostro et al. (2006) Ostro, S., Margot, J.-L., Benner, L., Giorgini, J., Scheeres, D., Fahnestock, E., Broschart, S., Bellerose, J., Nolan, M., Magri, C., Pravec, P., Scheirich, P., Rose, R., Jurgens, R., DeJong, E., Suzuki, S., 2006. Radar imaging of binary near-Earth asteroid (66391) 1999 KW4. Science 314 (5803), 1276–1280.
  * Peale (1977) Peale, S., 1977. Rotation histories of the natural satellites. In: Burns, J. A. (Ed.), IAU Colloq. 28: Planetary Satellites. University of Arizona, Tucson, pp. 87–112.
  * Pravec and Harris (2007) Pravec, P., Harris, A. W., 2007. Binary asteroid population. 1. angular momentum content. Icarus 190 (1), 250–259.
  * Pravec et al. (2014) Pravec, P., Scheirich, P., Ďurech, J., Pollock, J., Kušnirák, P., Hornoch, K., Galád, A., Vokrouhlický, D., Harris, A., Jehin, E., Manfroid, J., Opitom, C., Gillon, M., Colas, F., Oey, J., Vraštil, J., Reichart, D., Ivarsen, K., Haislip, J., LaCluyze, A., 2014. The tumbling spin state of (99942) Apophis. Icarus 233, 48 – 60.
URL http://www.sciencedirect.com/science/article/pii/S0019103514000578

  * Pravec et al. (2016) Pravec, P., Scheirich, P., Kušnirák, P., Hornoch, K., Galád, A., Naidu, S., Pray, D., Világi, J., Gajdoš, Š., Kornoš, L., Krugly, Y., Cooney, W., Gross, J., Terrell, D., Gaftonyuk, N., Pollock, J., Husárik, M., Chiorny, V., Stephens, R., Durkee, R., Reddy, V., Dyvig, R., Vraštil, J., Žižka, J., Mottola, S., Hellmich, S., Oey, J., Benishek, V., Kryszczyńska, A., Higgins, D., Ries, J., Marchis, F., Baek, M., Macomber, B., Inasaridze, R., Kvaratskhelia, O., Ayvazian, V., Rumyantsev, V., Masi, G., Colas, F., Lecacheux, J., Montaigut, R., Leroy, A., Brown, P., Krzeminski, Z., Molotov, I., Reichart, D., Haislip, J., LaCluyze, A., 2016. Binary asteroid population. 3. Secondary rotations and elongations. Icarus 267, 267–295.
URL https://www.sciencedirect.com/science/article/pii/S0019103515005722

  * Pravec et al. (2006) Pravec, P., Scheirich, P., Kušnirák, P., Šarounová, L., Mottola, S., Hahn, G., Brown, P., Esquerdo, G., Kaiser, N., Krzeminski, Z., Pray, D., Warner, B., Harris, A., Nolan, M., Howell, E., Benner, L., Margot, J.-L., Galád, A., Holliday, W., Hicks, M., Krugly, Y., Tholen, D., Whiteley, R., Marchis, F., DeGraff, D., Grauer, A., Larson, S., Velichko, F., Cooney, W., Stephens, R., Zhu, J., Kirsch, K., Dyvig, R., Snyder, L., Reddy, V., Moore, S., Gajdoš, Š., Világi, J., Masi, G., Higgins, D., Funkhouser, G., Knight, B., Slivan, S., Behrend, R., Grenon, M., Burki, G., Roy, R., Demeautis, C., Matter, D., Waelchli, N., Revaz, Y., Klotz, A., Rieugné, M., Thierry, P., Cotrez, V., Brunetto, L., Kober, G., 2006. Photometric survey of binary near-earth asteroids. Icarus 181 (1), 63–93.
URL https://www.sciencedirect.com/science/article/pii/S0019103505003891

  * Quillen (2006) Quillen, A. C., 2006. Reducing the probability of capture into resonance. Monthly Notices of the Royal Astronomical Society 365 (4), 1367–1382.
  * Quillen et al. (2016a) Quillen, A. C., Giannella, D., Shaw, J. G., Ebinger, C., 2016a. Crustal failure on icy moons from a strong tidal encounter. Icarus 275, 267–280.
  * Quillen et al. (2016b) Quillen, A. C., Kueter-Young, A., Frouard, J., Ragozzine, D., 2016b. Tidal spin down rates of homogeneous triaxial viscoelastic bodies. Monthly Notices of the Royal Astronomical Society 463, 1543–1553.
  * Quillen et al. (2020) Quillen, A. C., Lane, M., Nakajima, M., Wright, E., 2020. Excitation of tumbling in Phobos and Deimos. Icarus 340, 113641.
  * Quillen et al. (2019a) Quillen, A. C., Martini, L., Nakajima, M., 2019a. Near/far side asymmetry in the tidally heated Moon. Icarus 329, 182–196, arXiv:1810.10676.
  * Quillen et al. (2017) Quillen, A. C., Nichols-Fleming, F., Chen, Y.-Y., Noyelles, B., 2017. Obliquity evolution of the minor satellites of Pluto and Charon. Icarus 293, 94–113.
  * Quillen et al. (2019b) Quillen, A. C., Wagner, K. J., Sanchez, P., 2019b. Simulations of wobble damping in viscoelastic rotators. Monthly Notices of the Royal Astronomical Society 485 (1), 725–738.
  * Rein and Liu (2012) Rein, H., Liu, S.-F., 2012. REBOUND: An open-source multi-purpose N-body code for collisional dynamics. Astronomy and Astrophysics 537, A128.
  * Renner and Sicardy (2006) Renner, S., Sicardy, B., 2006. Use of the geometric elements in numerical simulations. Celestial Mechanics and Dynamical Astronomy, 94, 237–248.
  * Rubincam (2000) Rubincam, D. P., 2000. Radiative spin‐up and spin‐down of small asteroids. Icarus 148 (1), 2–11.
  * Scheeres (2007) Scheeres, D. J., 2007. The dynamical evolution of uniformly rotating asteroids subject to YORP. Icarus 188, 430–450.
  * Scheeres (2015) Scheeres, D. J., 2015. Landslides and mass shedding on spinning spheroidal asteroids. Icarus 247, 1–17.
  * Scheirich and Pravec (2009) Scheirich, P., Pravec, P., 2009. Modeling of lightcurves of binary asteroids. Icarus 200, 531–547.
  * Scheirich et al. (2015) Scheirich, P., Pravec, P., Jacobson, S., Ďurech, J., Kušnirák, P., Hornoch, K., Mottola, S., Mommert, M., Hellmich, S., Pray, D., Polishook, D., Krugly, Y., Inasaridze, R., Kvaratskhelia, O., Ayvazian, V., Slyusarev, I., Pittichová, J., Jehin, E., Manfroid, J., Gillon, M., Galád, A., Pollock, J., Licandro, J., Alí-Lagoa, V., Brinsfield, J., Molotov, I., 2015\. The binary near-Earth asteroid (175706) 1996 FG3 — an observational constraint on its orbital evolution. Icarus 245, 56–63.
URL https://www.sciencedirect.com/science/article/pii/S0019103514004850

  * Scheirich et al. (2021) Scheirich, P., Pravec, P., Kušnirák, P., Hornoch, K., McMahon, J., Scheeres, D., Čapek, D., Pray, D., Kučáková, H., Galád, A., Vraštil, J., Krugly, Y. N., Moskovitz, N., Avner, L., Skiff, B., McMillan, R., Larsen, J., Brucker, M., Tubbiolo, A., Cooney, W., Gross, J., Terrell, D., Burkhonov, O., Ergashev, K., Ehgamberdiev, S., Fatka, P., Durkee, R., Schunova, E. L., Inasaridze, R. Y., Ayvazian, V., Kapanadze, G., Gaftonyuk, N., Sanchez, J., Reddy, V., McGraw, L., Kelley, M., Molotov, I., 2021\. A satellite orbit drift in binary near-earth asteroids (66391) 1999 KW4 and (88710) 2001 SL9 — indication of the BYORP effect. Icarus 360, 114321\.
URL https://www.sciencedirect.com/science/article/pii/S0019103521000208

  * Seligman and Batygin (2021) Seligman, D., Batygin, K., 2021. The onset of chaos in permanently deformed binaries from spin-orbit and spin-spin coupling. The Astrophysical Journal 913 (1), 31.
  * Steinberg and Sari (2011) Steinberg, E., Sari, R., 2011. Binary YORP and evolution of binary asteroids. Astronomical Journal 141, 55–65.
  * Taylor and Margot (2011) Taylor, P. A., Margot, J.-L., 2011. Binary asteroid systems: Tidal end states and estimates of material properties. Icarus 212, 662–676.
  * Tricarico et al. (2021) Tricarico, P., Scheeres, D., French, A., McMahon, J., Brack, D., Leonard, J., Antreasian, P., Chesley, S., Farnocchia, D., Takahashi, Y., Mazarico, E., Rowlands, D., Highsmith, D., Getzandanner, K., Moreau, M., Johnson, C., Philpott, L., Bierhaus, E., Walsh, K., Barnouin, O., Palmer, E., Weirich, J., Gaskell, R., Daly, M., Seabrook, J., Nolan, M., Lauretta, D., dec 2021. Internal rubble properties of asteroid (101955) bennu. Icarus 370, 114665.
URL https://doi.org/10.1016%2Fj.icarus.2021.114665

  * Walsh et al. (2008) Walsh, K. J., Richardson, D. C., Michel, P., 2008. Rotational breakup as the origin of small binary asteroids. Nature Letters 454 (7), 188–191.
  * Wisdom (1987) Wisdom, J., 1987. Rotational dynamics of irregularly shaped natural satellites. Astronomical Journal 94 (5), 1350–1360.