# Numerical Simulations of Tidal Deformation and Resulting Light Curves of
Small Bodies: Material Constraints of 99942 Apophis and 1I/‘Oumuamua

Aster G. Taylor Dept. of Astronomy and Astrophysics, University of Chicago,
5640 S Ellis Ave, Chicago, IL 60637  Darryl Z. Seligman Dept. of Astronomy and
Carl Sagan Institute, Cornell University, 122 Sciences Drive, Ithaca, NY,
14853, USA  Douglas R. MacAyeal Dept. of Geophysical Sciences, University of
Chicago, 5734 S Ellis Ave, Chicago, IL 60637  Olivier R. Hainaut European
Southern Observatory, Karl-Schwarzschild-Strasse 2, Garching bei München,
D-85748, Germany  Karen J. Meech Institute for Astronomy, University of
Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

###### Abstract

In this paper, we present an open source software (SAMUS) which simulates
constant-density, constant-viscosity liquid bodies subject to tidal forces for
a range of assumed viscosites and sizes. This software solves the Navier-
Stokes equations on a finite-element mesh, incorporating the centrifugal,
Coriolis, self-gravitational, and tidal forces. The primary functionality is
to simulate the deformation of minor bodies under the influence of tidal
forces. It may therefore be used to constrain the composition and physical
structure of bodies experiencing significant tidal forces, such as 99942
Apophis and 1I/‘Oumuamua. We demonstrate that SAMUS will be useful to
constrain the material properties of Apophis during its near-Earth flyby in
2029. Depending on the material properties, Apophis may experience an area
change of up to 0.5%, with similar effects on the photometric brightness. We
also apply SAMUS to constrain the material dynamic viscosity of 1I/‘Oumuamua,
the first interstellar object discovered traversing the inner Solar System.
‘Oumuamua experienced a close approach to the Sun at perihelion ($q\simeq
0.25$ au) during which there were significant tidal forces that may have
caused deformation of the body. This deformation could have lead to observable
changes in the photometric light curve based on the material properties. The
application of SAMUS to produce synthetic observations which incorporate tidal
deformation effects demonstrate that no deformation — an infinite dynamic
viscosity — best reproduces the photometric data. While these results indicate
that ‘Oumuamua did not experience significant tidal deformation, a
sophisticated model incorporating non-principal axis rotation is necessary to
conclusively analyze both ‘Oumuamua and Apophis.

Interstellar Objects (52); Comets (280); Hydrodynamics (1963)

##  1 Introduction

Tidal gravitational forces are one of the primary drivers of the evolution of
the objects within the Solar System. For example, the geosynchronous orbit of
the Moon was explained by Darwin (1879, 1880) as a natural consequence of
tidal dissipation. The theory of spin evolution of satellites due to tidal
forces was improved vastly in the following centuries (Kaula, 1964; MacDonald,
1964). The discovery of the 3:2 spin orbit resonance of Mercury (Pettengill &
Dyce, 1965) led to significant theoretical advances on long-term tidal
dissipation (Peale & Gold, 1965; Liu & O’Keefe, 1965; Colombo, 1965; Goldreich
& Peale, 1966, 1968; Bagheri et al., 2022).

The theory of tides has led to multiple predictions that were later
corroborated by data. An archetypal example of this was the prediction of
vulcanism on Io by Peale et al. (1979) which was later verified by the Voyager
1 spacecraft (Smith et al., 1979). Similarly, Wisdom et al. (1984) predicted
that the Saturnian satellite Hyperion was chaotically tumbling, which was
confirmed by Voyager 2 (Black et al., 1995). For a more recent example, Fuller
et al. (2016) predicted the rapid outward migration of the Saturnian satellite
Titan due to sustained excitation of inertial waves via tidal forces, which
was then measured by the Cassini spacecraft (Lainey et al., 2020).

Tidal deformation can also affect small bodies like comets and asteroids, and
in the most dramatic cases, lead to catastrophic disruption events. A
historical example of this was the tidal disruption of the comet Shoemaker-
Levy 9 into Jupiter’s atmosphere (Shoemaker et al., 1993; Weaver et al., 1995;
Lellouch et al., 1995; Noll et al., 1995). Chains of craters on the surfaces
of Callisto and Ganymede have also been explained by similarly catastrophic
tidal events (Schenk et al., 1996). In milder cases tidal deformation can
result in satellite formation via material stripping (Richardson et al.,
1998), which has been invoked to explain the near ubiquity of doublet craters
on Solar System bodies (Melosh & Stansberry, 1991; Bottke & Melosh, 1996b, a;
Melosh et al., 1996; Cook et al., 2003). In this paper, we consider the
mildest case, where tidal forces lead to body deformation, with application to
the near-Earth object 99942 Apophis and the interstellar object 1I/‘Oumuamua
(MPEC, 2017; Williams et al., 2017). For recent reviews on interstellar
objects, see Jewitt & Seligman (2022) and Moro-Martín (2022).

Apophis was discovered in 2004 and subsequently identified as a potential
Earth impact threat. While follow-up radar observations have eliminated any
chance of impact within a century, Apophis will pass within $\sim$6 Earth
radii in 2029 (Brozović et al., 2018), providing an excellent opportunity for
observation and characterization. The OSIRIS-REx mission, after returning from
the asteroid 101955 Bennu, will be renamed OSIRIS-APEX and dispatched to
encounter and observe Apophis (Nolan et al., 2021). Radar observations
(Brozović et al., 2018) revealed that Apophis has a bi-lobed shape.
Observations also revealed that Apophis exhibits non-principal-axis rotation
(Brozović et al., 2018; Pravec et al., 2014; Lee et al., 2022), and measured a
non-zero Yarkovsky effect (Pérez-Hernández & Benet, 2022).

The potential effects of Apophis’ near-Earth flyby on the geophysical,
photometric, and rotational states have been investigated. Benson et al.
(2023) demonstrated that the close encounter will induce sufficient tidal
torques to modify the rotational state. The tidal forces will also be
sufficient to produce resurfacing events which could modify the photometric
properties (Yu et al., 2014; Kim et al., 2023). The magnitude and nature of
these effects are potentially effective probes of the geophysical properties
of Apophis, including density distribution (Dinsmore & de Wit, 2022), rigidity
and dissipation (Hirabayashi, 2022), and seismic response (DeMartini et al.,
2019). We investigate the material properties for which tidal forces may
induce a detectable shape deformation in the object, enabling SAMUS to
constrain the material properties of Apophis. Similar effects are investigated
for I1/‘Oumuamua.

After the October 2017 discovery of ‘Oumuamua, there was an immediate
acquisition of ground- and space-based observations of the rapidly fading
object. These observations produced a high-quality composite light curve
spanning approximately a month (29.3 days) and a spatial segment of $l\simeq
0.13$ au. In total there were 818 observations, reported by Meech et al.
(2017), Bolin et al. (2017), Bannister et al. (2017), Drahus et al. (2017),
Fraser et al. (2018), Jewitt et al. (2017), Knight et al. (2017), and Belton
et al. (2018). These observations were collectively reported in Belton et al.
(2018).

‘Oumuamua was interpreted to be highly elongated, with an aspect ratio
estimated to be $>$3:1 (Bolin et al., 2017; Knight et al., 2017), $>$5:1
(Bannister et al., 2017; Fraser et al., 2018; Jewitt et al., 2017), and up to
10:1 (Meech et al., 2017). Frequency analysis of the light curve showed a
maximum at a period of $p\simeq 4.3$ hours (Belton et al., 2018). This was
interpreted to be half of the rotational period — corresponding to a
revolution of 180∘. ‘Oumuamua’s variations in absolute magnitude of $H\simeq
22.5\pm 1.3$ over its rotation led to the conclusion that it was exhibiting
complex, non-principal axis rotation (Drahus et al., 2017; Meech et al., 2017;
Fraser et al., 2018).

This analysis was further refined by Mashchenko (2019) who demonstrated via
full light curve modeling that a near-symmetric oblate ellipsoid with
dimensions of 115:111:19$\sim$6:6:1 meters provided a best fit geometry for
the light curve data. This size estimate assumes a geometric albedo of $A=0.1$
and would change with a different albedo, but with appropriately scaled
dimensions. While a prolate ellipsoid with dimensions of 342:42:42 meters is
also allowable, the torques required to replicate the motion are highly tuned,
so the prolate geometry is disfavored.

Deep imaging revealed a notable lack of cometary activity, classifying
‘Oumuamua as an asteroidal body and restricting possible dust outputs (upper
limits range from $\sim 2\cdot 10^{-4}$ kg s-1 (Jewitt et al., 2017) to
$1.7\cdot 10^{-3}$ kg s-1 (Meech et al., 2017)). Additionally, while outbound
at 2 au, there was a significant non-detection of the object with the Spitzer
Space Telescope. This non-detection placed limits on the production of CO and
CO2 (Trilling et al., 2018).

Astrometric positional data revealed that the trajectory was inconsistent with
pure Keplerian motion (Micheli et al., 2018). The addition of a radially
outward non-gravitational acceleration of the form $a=4.92\cdot
10^{-4}\,(r/1\,{\rm au})^{-2}\,\boldsymbol{\hat{r}}$ cm s-2 provides a greatly
improved, $30\sigma$ fit to the trajectory.111An acceleration with a form of
$r^{-1}$ is nearly as good of a fit. Micheli et al. (2018) proposed a comet-
like outgassing as an explanation for this acceleration, ruling out radiation
pressure, the Yarkovsky effect, magnetic forces, and others which would
require extreme physical properties.

The restrictions on the coma and micron-scale dust presence in the vicinity
are quite stringent. Therefore, theories for the provenance of the object
positing cometary outgassing as the source of the acceleration require
additional complexity to avoid violating the Spitzer or photometric
observations.222Although the original Spitzer estimates had a computational
error, see Seligman et al. (2021). Even these revised CO limits, however, are
prohibitive to the $1/r^{2}$ fit. H2O ice was initially proposed as an
outgassing accelerant because it is the most common volatile in Solar System
comets (Rickman, 2010; A’Hearn et al., 2012; Ootsubo et al., 2012; Cochran et
al., 2015; Biver & Bockelée-Morvan, 2016; Bockelée-Morvan & Biver, 2017) and
its presence is not in tension with the Spitzer non-detection. However, the
relatively high enthalpy of sublimation (51 kJ mol-1) of H2O implies that
water sublimation would require more energy input than ‘Oumuamua received from
solar radiation (Sekanina, 2019).

Attempting to unify these constraints, Seligman & Laughlin (2020) argued that
only hypervolatiles could serve as the accelerant for ‘Oumuamua. They found
that only molecular hydrogen (H2), neon, molecular nitrogen (N2), and argon
were allowable accelerants for an oblate spheroid — although CO was also shown
to be energetically feasible. Those authors also investigated the feasibility
of hydrogen ice as the bulk constituent — originally hypothesized by
Füglistaler & Pfenniger (2018) — as it requires the lowest active surface
fraction to be explanatory. In this hypothesis, ‘Oumuamua would have formed in
a failed prestellar core in a Giant Molecular Cloud. This model naturally
explains the extreme shape (via continuous H2 ablation), the low excess
velocity speed, and young age (Mamajek, 2017; Gaidos et al., 2017; Feng &
Jones, 2018; Almeida-Fernandes & Rocha-Pinto, 2018; Hallatt & Wiegert, 2020;
Hsieh et al., 2021). However, there are theoretical barriers to the formation
of macroscopic bodies composed of solid hydrogen, such as the frigid
temperatures required for formation and rapid evaporation in the interstellar
medium (Hoang & Loeb, 2020; Phan et al., 2021; Levine et al., 2021; Levine &
Laughlin, 2021). Although the low condensation temperature of molecular
hydrogen ($<$10K) poses difficulties for its formation, Levine & Laughlin
(2021) demonstrated that adiabatic expansion pockets were a plausible
formation environment for such an object. Jackson & Desch (2021) instead
suggested that ‘Oumuamua was composed of molecular nitrogen (N2) ice, while
Desch & Jackson (2021) proposed that impacts on extrasolar Pluto analogues
would provide a plausible source for objects like ‘Oumuamua. However, Levine
et al. (2021) demonstrated that the necessary mass density for this formation
to be plausible is unreasonably high. Seligman et al. (2021) found that a
typographical mistake in Trilling et al. (2018) led to the reported outgassing
limits of CO to be underestimated by two orders of magnitude. When this error
was corrected for, those authors showed that a body characterized by a modest
covering fraction of CO exhibiting sporadic activity could explain both the
acceleration and Spitzer observations.

Bergner & Seligman (2023) demonstrated that the crystallization of amorphous
water ice would produce sufficient radiolytically produced and entrapped H2 to
provide the observed nongravitational acceleration. This crystallization would
occur in the absence of sublimation of the overall ice matrix, thereby
explaining the lack of dust coma observed. Further more, Chesley et al.
(2016), Farnocchia et al. (2022) and Seligman et al. (2023) reported
statistically significant non-gravitational accelerations on seven near-Earth
objects (NEOs) that did not display visible activity, similar to ‘Oumuamua.

Motivated by these advances, we reconsider the hypothesis that ‘Oumuamua was
an icy object. Flekkøy et al. (2019) considered the effects of tidal forces in
the context of a dust aggregate and demonstrated a surprising stability
against tidal stresses. We build upon this work to consider shearing effects
due to tidal forces in the context of an icy cometary body. At perihelion,
‘Oumuamua passed within $r_{H}\simeq 0.256$ au of the Sun, subjecting it to
non-trivial tidal stresses — although attempts to pre-discover the object here
resulted in non-detections (Hui & Knight, 2019). The tidal stresses, combined
with the large-magnitude centrifugal force produced by ‘Oumuamua’s high aspect
ratio and rapid rotation, produce significant shearing stress on the body. In
this paper, we investigate the effects of these shearing forces and provide
constraints on the size and dynamic viscosity of ‘Oumuamua. The latter could
potentially provide a constraint on the possible material composition of
‘Oumuamua.

##  2 Numerical Simulations of Tidal Deformation

Many small asteroids are loosely-bound granular bodies. Characterization of
plastic deformation of these bodies under rotational (Hirabayashi, 2014, 2015;
Hirabayashi & Scheeres, 2019) and tidal stress (Kim et al., 2021) has been the
subject of much research, and numerical models similar to those developed and
presented in this paper are often applied to these objects (Hirabayashi &
Scheeres, 2019; Sánchez et al., 2017). Numerical calculations incorporating
finite-element and soft-sphere discrete approximations are routinely used to
analyze the tensile strength of these objects with and without cohesive
forces. In this paper, we instead model these objects as fluid masses, which
greatly simplifies the analysis. While this approximation is a simplification,
it is appropriate for both semi-crystalline solids such as ices and — in
certain circumstances — for granular materials subjected to cohesive and
friction forces. This simplification reduces the problem to only two degrees
of freedom — the dynamic viscosity $\mu$ and the density $\rho$ — which fully
describe the material properties and state of ‘Oumuamua. This simplification
enables us to model the deformation explicitly over time, instead of relying
on the ‘deformation modes’ identified in Hirabayashi & Scheeres (2019).

In this section, we present a generalized software, Simulator of Asteroid
Malformation Under Stress (SAMUS), which simulates the deformation of
constant-density constant-viscosity liquid-body ellipsoids under forcing
pressures, which can be applied to constrain the dynamic viscosity and size of
minor objects. SAMUS incorporates tidal, centrifugal, Coriolis, and self-
gravitational forces for minor bodies, and allows for customized trajectory,
principal axes, rotational period, density, dynamic viscosity, rotational
axis, and simulation cutoffs. It is accessible on PyPi and at Taylor (2023),
and can be installed via pip.

SAMUS currently implements several simplifying assumptions due to
computational practicalities. SAMUS uses a “fixed-axis” NPA rotation and a
constant rotational period; tumbling and 3-dimensional rotation are not
currently included. However, this is not stringently enforced in the model.
The rotation axis is used to compute the non-inertial forces and to calculate
the tidal force in the body frame. However, SAMUS does allow the orientation
of the body to be shifted by the tidal forces, allowing for a more complex and
accurate evolution. However, it is well-known that tidal forces can induce
angular momentum axis drift, and so the incorporation of that effect into
SAMUS will be incorporated in future versions of this software — currently,
the rotational axis and magnitude are kept fixed, and the angular momentum
evolves with the changing shape.

SAMUS also assumes that the simulated object is not subjected to ablation, and
that shape changes only arise from specified shearing forces. This package
also assumes that the simulated body is homogeneous, both in density and in
dynamic viscosity. Given the temperature and material dependence of both of
these properties, this implies a constant temperature and material over the
body — at least across the scale of the cell size. Crystallization of
amorphous materials would induce (minor) changes in the viscosity of the
material and change the temperature, but this effect is currently not included
in these simulations.

###  2.1 Numerical Calculations

The SAMUS simulation software is written in Python 3.8.10 (Van Rossum & Drake,
2009), and is primarily based on the FEniCS (Alnaes et al., 2015; Logg et al.,
2012a), UFL (Alnaes et al., 2014), and DOLFIN (Logg & Wells, 2010; Logg et
al., 2012b) packages. It also has dependencies on NumPy (Harris et al., 2020),
SciPy (Virtanen et al., 2020), pandas (pandas development team, 2020;
McKinney, 2010), quaternion (Boyle et al., 2022), and MPI for Python (Dalcín
et al., 2005, 2008, 2011; Dalcín & Fang, 2021). All of these must be installed
in the user’s distribution. SAMUS is primarily structured as a Python class,
and solves the weak formulation of the partial differential Navier-Stokes
equations over a finite-element mesh.

The domain used by SAMUS is an $\mathcal{S}_{3}$ (3-ball) domain (created by
Gmsh (Geuzaine & Remacle, 2020) and loaded into DOLFIN), distorted into an
ellipsoid with principal axes $a,b,c$. After reading in the provided
trajectory data, body parameters, and simulation parameters, SAMUS uses an
Euler finite-difference approximation to iteratively solve the Navier-Stokes
equations. The mesh is advectively updated at each time step to simulate the
tidal deformation, where the computed fluid velocity is used to find the
displacement vector.

FEniCS is used to solve the weak formulation of the incompressible Navier-
Stokes equations with Dirichlet boundary conditions. In the following
equations, $\boldsymbol{u}=\partial\boldsymbol{r}/\partial t$ is the velocity,
$\boldsymbol{r}$ is the position in the co-moving, non-inertial frame, $\rho$
is the density, $\mu$ is the dynamic viscosity, $p$ is the pressure, and
$\boldsymbol{\Omega}$ is the angular velocity vector. $V$ and $Q$ are function
spaces over $\mathds{R}^{3}$ and $\mathds{R}$ respectively,333Defined as
continuous Galerkin domains. with $\boldsymbol{u},\boldsymbol{v}\in V$ and
$p,q\in Q$, and $\boldsymbol{v},q$ are test functions. The differential
‘$\text{d}x$’ represents a volume integral over the body domain.

The strong momentum equation is

$\rho\frac{\partial\boldsymbol{u}}{\partial
t}+\rho(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}-\nabla\cdot\boldsymbol{\sigma}(\boldsymbol{u},p)=\boldsymbol{f}\,.$
(1)

SAMUS uses the weak form, which is

$\displaystyle\rho\frac{\partial\boldsymbol{u}}{\partial
t}\cdot\boldsymbol{v}\
\text{d}x+\mu\nabla\boldsymbol{u}\cdot\nabla\boldsymbol{v}\ \text{d}x$ (2)
$\displaystyle+$
$\displaystyle\rho(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}\cdot\boldsymbol{v}\
\text{d}x-p\nabla\cdot\boldsymbol{v}\ \text{d}x$ $\displaystyle=$
$\displaystyle\boldsymbol{F}_{\text{tidal}}\cdot\boldsymbol{v}\
\text{d}x+\boldsymbol{g}\cdot\boldsymbol{v}\ \text{d}x$ $\displaystyle-$
$\displaystyle\rho(\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\boldsymbol{r}))\cdot\boldsymbol{v}\
\text{d}x$ $\displaystyle-$ $\displaystyle
2\rho(\boldsymbol{\Omega}\times\boldsymbol{u})\cdot\boldsymbol{v}\ \text{d}x\
\ \ \forall\boldsymbol{v}\in V.$

The mass continuity equation, in its strong form, is

$\nabla\cdot\boldsymbol{u}=0,$ (3)

and SAMUS again uses the weak form,

$q\nabla\cdot\boldsymbol{u}\ \text{d}x=0\ \ \ \forall q\in Q.$ (4)

The derivation of the weak form from the strong form of the Navier-Stokes
equations is given in Quarteroni (2014). The right-hand side of the momentum
equation represents the forcing, each term of which is defined in Table 1.

FORCING TERMS
---
Tidal force: | $\boldsymbol{F}_{\text{tidal}}$
Self-gravitational force: | $\boldsymbol{g}$
Centrifugal force: | $-\rho(\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\boldsymbol{r}))$
Coriolis force: |  $-2\rho(\boldsymbol{\Omega}\times\boldsymbol{u})$.
Table 1: The forcing terms in the Navier-Stokes equations.

In SAMUS, the acceleration $\partial\boldsymbol{u}/\partial t$ is estimated
with an Euler finite-difference method. For a given time step indexed by $i$,
the acceleration is approximated as $\partial\boldsymbol{u}/\partial
t\simeq(\boldsymbol{u}_{i}-\boldsymbol{u}_{i-1})/\Delta t$. The time step
$\Delta t$ is adaptively modified to ensure that the Courant–Friedrichs–Lewy
(CFL) condition ($|\boldsymbol{u}|\Delta t/\Delta x<C_{\text{max}}$) is met
(Courant et al., 1967). $C_{\text{max}}$ can be user-defined and is set to a
default of 1, which is the standard limit. We performed extensive stability
and convergence tests, which are available in the SAMUS package (and not
described in this paper).

SAMUS additionally uses FEniCS to rapidly solve the weak form of the self-
gravitational force as given by Gauss,

$c\nabla\cdot\boldsymbol{g}\ \text{d}x=-4\pi G\rho\ c\ \text{d}x\ \ \ \forall
c\in Q\,,$ (5)

where $c$ is a scalar test function in $Q$. This method is relatively rapid
and allows for efficient computation of self-gravity for even highly distorted
bodies.

SAMUS produces a `csv` file containing timestamps, the maximum dimension of
the body on each axis, and the moment of inertia $I$ at each step, which is
calculated using

$I=\int_{S}\rho\,\,\bigg{(}\,\frac{\|\boldsymbol{r}\times\boldsymbol{\Omega}\|^{2}}{\|\boldsymbol{\Omega}\|^{2}}\
\bigg{)}\,\text{d}x\,.$ (6)

SAMUS is further capable of incorporating a broad class of user-defined
functions in these outputs.

###  2.2 Trajectory Jump Method

SAMUS uses a “trajectory jump” method for efficiency, reducing the number of
computations necessary over the trajectory. SAMUS first computes the time-
averaged deformation over a (user-defined) number of rotational periods. It
then performs a linear extrapolation of this average and steps forward in
simulation time until either (i) the CFL condition is violated or (ii) the
heliocentric distance changes by 1% (this threshold is similarly user-
defined). This method assumes that the rotational period of the body is
significantly shorter than the timescale within which the object moves through
its trajectory significantly. The quality of this first-order linear
approximation of the distortion was validated with convergence tests, using
halved tolerances (available in the package). However, for bodies with slower
rotation, further testing should be performed to confirm the validity of this
methodology.

###  2.3 Tidal Force Computation

In SAMUS, quaternion is used to rotate points to the stationary frame, a
necessary step to compute the the continuum tidal force over
$\mathcal{S}_{3}$. Quaternions allow for rapid computation of rotation by an
arbitrary angle about an arbitrary axis, without gimbal locking (Kuipers,
2007). For a given point $\boldsymbol{r}$ and quaternion $\boldsymbol{q}$, the
rotated point is
$\boldsymbol{r}^{\prime}=\boldsymbol{q}\boldsymbol{r}\boldsymbol{q}^{*}$. A
rotation by an angle $\theta$ about an axis
$\boldsymbol{\hat{\Omega}}=(\Omega_{x},\Omega_{y},\Omega_{z})$ is described by
$\boldsymbol{q}=\cos{(\theta/2)}+\Omega_{x}\sin{(\theta/2)}\boldsymbol{i}+\Omega_{y}\sin{(\theta/2)}\boldsymbol{j}+\Omega_{z}\sin{(\theta/2)}\boldsymbol{k}$
and conjugate
$\boldsymbol{q^{*}}=\cos{(\theta/2)}-\Omega_{x}\sin{(\theta/2)}\boldsymbol{i}-\Omega_{y}\sin{(\theta/2)}\boldsymbol{j}-\Omega_{z}\sin{(\theta/2)}\boldsymbol{k}$.
Here, $\boldsymbol{i}$, $\boldsymbol{j}$, and $\boldsymbol{k}$ are the
imaginary unit quaternions.

The position quaternion
$\boldsymbol{r}^{\prime}=x^{\prime}\boldsymbol{i}+y^{\prime}\boldsymbol{j}+z^{\prime}\boldsymbol{k}$
in the co-rotating frame is given by
$\boldsymbol{r}=\boldsymbol{q^{*}}\boldsymbol{x^{\prime}}\boldsymbol{q}$ in
the non-rotating frame. The tidal force at each point on the body is then
computed by SAMUS using

$\boldsymbol{F}_{\text{tidal}}(\boldsymbol{r}^{\prime})=-\boldsymbol{\hat{x}}GM_{\sun}\rho\left(\frac{1}{(r_{H}-(\boldsymbol{q^{*}}\boldsymbol{r}^{\prime}\boldsymbol{q})_{x})^{2}}-\frac{1}{r_{H}^{2}}\right)\,,$
(7)

for an object with constant density $\rho$.

A pseudocode describing the algorithm used in SAMUS is given in Appendix B.

##  3 Analytic Approximation of Tidal Deformation

In this section, we present analytic estimates to approximate the tidal
deformation. For a point mass at some distance $R$ from an object of mass $M$,
the magnitude of the acceleration due to gravity is

$a_{g}=\frac{GM}{R^{2}}\,.$ (8)

In Equation 8, G is the gravitational constant. The difference in the
acceleration across a distance $\Delta r$ is the tidal acceleration $a_{\rm
tidal}$, which is given by

$a_{\rm tidal}=GM\,\bigg{|}\frac{1}{(R+\Delta
r)^{2}}-\frac{1}{R^{2}}\bigg{|}\,.$ (9)

This simplifies to

$a_{\rm tidal}=GM\bigg{|}\frac{\Delta r(\Delta r+2R)}{R^{2}(R+\Delta
r)^{2}}\bigg{|}\,.$ (10)

Assuming that $\Delta r\ll R$ and Taylor expanding about $\Delta r=0$ yields,

$a_{\rm tidal}\simeq 2GM\,\bigg{|}\frac{\Delta r}{R^{3}}\bigg{|}\,.$ (11)

Now let us consider the scale of deformation experienced by an object
subjected to an external acceleration $a$. The force experienced by a
continuum object of density $\rho$ is $a\rho$. The magnitude of the tidal
force per unit volume, $F_{\rm tidal}$, is

$F_{\rm tidal}\simeq 2GM\rho\,\bigg{|}\frac{\Delta r}{R^{3}}\bigg{|}\,.$ (12)

For the dynamic shear viscosity $\mu$, we use the definition that
$\sigma=2\mu\dot{\epsilon}$, where $\sigma$ is the stress and $\dot{\epsilon}$
is the time derivative of the strain. For a domain of length $L$, the strain
$\epsilon$ is defined as the fractional change in the length of the domain,
and its time derivative is simply $\dot{\epsilon}=\dot{L}(t)/L_{0}$, where
$L_{0}$ is the initial domain length. On the other hand, the stress is the
force per unit area, so we will write $\sigma=F_{\rm tidal}L_{0}$ with $L_{0}$
the characteristic length.

Figure 1:  Distance at which an object of 100 m is deformed by 30%, versus
dynamic viscosity (in poise) and interaction timescale (in hours). Results are
presented for cases in which the Sun, the Earth, and Jupiter are the primary
orbital bodies.

Therefore, we write

$\underbrace{\frac{2GM\,\rho\Delta
r}{R^{3}}}_{\sigma}=\underbrace{\frac{2\mu}{L_{0}}\frac{dL(t)}{dt}}_{2\mu\dot{\epsilon}}\,.$
(13)

In this equation, the left-hand side is $\sigma$, the force times a
characteristic length, and the right-hand side is $2\mu\dot{\epsilon}$. This
approximate equation can be solved, using Equation 12 and setting $\Delta
r=L(t)$, since we wish to find the deformation of the entire object:

$L(t)=L_{0}\,\exp{\bigg{[}\frac{GML_{0}^{2}\rho}{\mu R^{3}}t\bigg{]}}\,.$ (14)

This equation allows us to compute the expected deformation for a given force
and over a given timescale. This equation also exhibits the expected
dependence — increasing with time, density, and initial length and decreasing
with the viscosity and the orbital distance.

Equation 14 allows for analysis of tidal deformation and material cohesion in
general small bodies. In Figure 1, we present heatmaps showing the distance at
which an object will deform by 30%, which we assume to be significant. The
value of 30% deformation is only an estimate, and a greater understanding of
small-body cohesion is necessary to provide stricter constraints on structural
collapse. We present these distances for objects of 100 m in radius, in orbit
around the Sun, the Earth, and Jupiter for a range of dynamic viscosity and
interaction timescales. These results may provide constraints on the cohesive
stability of Sun-grazing comets and on tidally destructed comets such as
Shoemaker-Levy 9. These results can also be used to constrain the cohesive
properties of small bodies in the solar system in future research.

##  4 Simplified Light Curve Model

In this section, we present a model to generate synthetic light curves of
minor bodies with the simplifying assumption of a “fixed-axis” NPA rotation.
To validate our methodology, we adopt idealized rotational states and orbital
geometries in Appendix A, which we compare against the more complex derived
values. We then describe a methodology to fit these synthetic data to
photometric light curves to obtain parameters for a characteristic “fixed-
axis” NPA rotation for minor bodies.

###  4.1 Light Curve Model

In our light curve model, we incorporate the change in the phase angle due to
the minor body’s astrometric progression in its orbit, while restricting the
rotation to follow a single “fixed-axis” NPA rotation. This is a stringent
simplification, as this sort of rotation is non-physical. However, this model
allows us to semi-accurately reproduce observed light curves and incorporate
the simulated tidal deformation for comparison, balancing physical accuracy
and practical restrictions.

In order to model the light curve, we use Equation 10 from Muinonen & Lumme
(2015), which we refer to as ‘ML15’ for the remainder of this paper. We assume
that the relevant objects have a diffuse Lommel-Seeliger scattering surface,
which represents a closely-packed particulate medium with weak multiple
scattering. Integrating the Lommel-Seeliger scattering function over the
exposed-and-visible surface gives an expression for the brightness at any
orientation.

We assume a single-scattering albedo $A$ and an isometric single-scattering
phase function $P(\alpha)=1$. This incorporates modulation based on the body
orientation, but assumes no additional modulation from the scattering
function. We define the phase angle $\alpha$ as the interior Sun-object-Earth
angle, and note that
$\cos\alpha=\boldsymbol{\hat{e}}_{\sun}\cdot\boldsymbol{\hat{e}}_{\earth}$,
where $\boldsymbol{\hat{e}}_{\sun}$ and $\boldsymbol{\hat{e}}_{\earth}$ are
unit vectors in the direction of the Sun and Earth respectively. We also
define the matrix $\boldsymbol{C}$ as

$\boldsymbol{C}\equiv\begin{pmatrix}a^{-2}&0&0\\\ 0&b^{-2}&0\\\
0&0&c^{-2}\end{pmatrix}\,,$ (15)

and the parameters $T_{\sun}$ and $T_{\earth}$ to be

$\displaystyle T_{\sun}\equiv$
$\displaystyle\sqrt{\boldsymbol{\hat{e}}_{\sun}^{T}C\boldsymbol{\hat{e}}_{\sun}}$
(16) $\displaystyle T_{\earth}\equiv$
$\displaystyle\sqrt{\boldsymbol{\hat{e}}_{\earth}^{T}C\boldsymbol{\hat{e}}_{\earth}}\,.$

We further define the parameter $T$ as

$T\equiv\sqrt{T_{\sun}^{2}+T_{\earth}^{2}+2T_{\sun}T_{\earth}\cos\alpha^{\prime}}\,,$
(17)

the angle $\alpha^{\prime}$ as

$\displaystyle\cos\alpha^{\prime}=$
$\displaystyle\frac{\boldsymbol{\hat{e}}_{\sun}^{T}C\boldsymbol{\hat{e}}_{\earth}}{T_{\sun}T_{\earth}}$
(18) $\displaystyle\sin\alpha^{\prime}=$
$\displaystyle\sqrt{1-\cos^{2}\alpha^{\prime}}\,,$

and the angle $\lambda^{\prime}$ as

$\displaystyle\cos\lambda^{\prime}=$
$\displaystyle\frac{T_{\sun}+T_{\earth}\cos\alpha^{\prime}}{T}$ (19)
$\displaystyle\sin\lambda^{\prime}=$
$\displaystyle\frac{T_{\sun}\sin\alpha^{\prime}}{T}\,.$

Then the disk-integrated brightness $L_{\rm ML15}$ is

$\displaystyle L_{{\rm ML15}}(\alpha)=$ $\displaystyle\frac{1}{8}\pi
F_{0}AP(\alpha)abc\frac{T_{\sun}T_{\earth}}{T}$ (20)
$\displaystyle\bigg{(}\cos(\lambda^{\prime}-\alpha^{\prime})+\cos\lambda^{\prime}+\sin\lambda^{\prime}\sin(\lambda^{\prime}-\alpha^{\prime})$
$\displaystyle\ln{\big{[}\cot(\frac{1}{2}\lambda^{\prime})\cot(\frac{1}{2}(\alpha^{\prime}-\lambda^{\prime}))\big{]}}\bigg{)},$

where $\pi F_{0}$ is the incident flux density. The absolute magnitude is
given by

$\displaystyle H=\Delta V-$ $\displaystyle
2.5\log\bigg{(}abc\frac{T_{\sun}T_{\earth}}{T}$ (21)
$\displaystyle\big{(}\cos(\lambda^{\prime}-\alpha^{\prime})+\cos\lambda^{\prime}+\sin\lambda^{\prime}\sin(\lambda^{\prime}-\alpha^{\prime})$
$\displaystyle\ln{\big{[}\cot(\frac{1}{2}\lambda^{\prime})\cot(\frac{1}{2}(\alpha^{\prime}-\lambda^{\prime}))\big{]}}\big{)}\bigg{)},$

where $\Delta V$ is a constant to absorb the flux and magnitude conversion
(Mashchenko, 2019).

Figure 2: Synthetic light curves of a body rotating with an arbitrary period
for a range of phase angle (A) and aspect ratio (B).

###  4.2 Phase Angle and Aspect Ratio Dependence

In this subsection, we investigate the dependence of the ML15 model on phase
angle and aspect ratio. This will be applicable to ‘Oumuamua in Section 4.3.
We compute synthetic light curves for a single (arbitrary) period with varying
phase angle and aspect ratio, and set the parameters ($\theta$, $\phi$,
$\psi$, $\beta_{0}$, and $\Delta V$) to the optimized values given in Section
4.3. However, we verified that the phase angle– and aspect ratio– dependent
behavior of the light curve does not sensitively depend on the other
parameters. This analysis both validates the ML15 model and explains
components of the synthetic light curves in Figure 13.

We present 40 synthetic light curves with phase angles ($\alpha$) uniformly
distributed in [$0$,$\pi/2$) and with an aspect ratio of 6:6:1 in Figure 2A.
For low phase angle, there are two distinct peaks at $\beta=\pi/2$ and
$\beta=3\pi/2$, caused by the half-period cycles of illumination from the
larger and smaller cross-sections. As $\alpha\rightarrow\pi/2$, the light
curve becomes approximately sinusoidal, as at $\alpha=\pi/2$ only a single
face of the body is observable, with brightness variation due to changing
exposure.

We also show 40 synthetic light curves for aspect ratios ranging from 1 to 20
and with $\alpha=\pi/2$ in Figure 2B. Here, the magnitude variation scales
with increasing aspect ratio and is zero for a sphere. This effect is simply
due to the increasing cross-sectional area contrast for larger aspect ratios.

###  4.3 Obtaining a Characteristic Rotation Axis

In this subsection, we describe the use of the ML15 model to find a
characteristic rotational state for an arbitrary body under the assumption of
a “fixed-axis” NPA rotation. To find this rotational state, we fit a synthetic
light curve generated with Equation 21 to photometric data (assumed to be
corrected for light travel time, helio- and geo-centric distance, and Solar
magnitude). We must carefully note that a “fixed-axis” NPA rotation is non-
physical, and does not exist in nature. As such, this model cannot be used to
describe the physical rotation of an object, but only to provide a well-
fitting light curve which can be compared to a deformed model.

We use the SciPy package’s scipy.optimize.curve_fit for the optimization,
which uses a non-linear least squares algorithm to minimize $\chi^{2}$, where

$\displaystyle\chi^{2}\equiv$
$\displaystyle\sum_{i}\left(\frac{(y_{i}-\mu_{i})^{2}}{\sigma_{y_{i}}^{2}}\right)\,.$
(22)

Here, $y_{i}$ and $\sigma_{y_{i}}$ denote photometric measurements and
associated errors, while $\mu_{i}$ denotes corresponding synthetic values.

The parameters that we optimize define the rotation axis and the Earth-
pointing axis. The Sun-pointing direction $\boldsymbol{\hat{e}}_{\sun}$ is
fixed to be along the $\boldsymbol{\hat{x}}$ direction, reducing symmetric
degeneracy. For a given phase angle $\alpha$, the direction of the observer
$\boldsymbol{\hat{e}}_{\earth}$ is therefore constrained to a cone centered on
the x-axis such that
$\boldsymbol{\hat{e}}_{\sun}\cdot\boldsymbol{\hat{e}}_{\earth}=\boldsymbol{\hat{x}}\cdot\boldsymbol{\hat{e}}_{\earth}=\cos\alpha$.
We then define a new variable $\theta$ such that the direction of the
observation is
$\boldsymbol{\hat{e}}_{\earth}=\cos\alpha\,\boldsymbol{\hat{x}}+\sin\alpha\cos\theta\,\boldsymbol{\hat{y}}+\sin\alpha\sin\theta\,\boldsymbol{\hat{z}}$,
with $\theta=0$ if $\boldsymbol{\hat{e}}_{\earth}$ lies in the x-y plane.

We additionally parameterize the rotation axis $\boldsymbol{\hat{\Omega}}$ by
spherical coordinates about the x-axis. We use a polar angle $\phi$ such that
$\boldsymbol{\hat{\Omega}}\cdot\boldsymbol{\hat{x}}=\cos\phi$ and an azimuthal
angle $\psi$ measured from the y-axis and restricted to the y-z plane. With
these definitions,
$\boldsymbol{\hat{\Omega}}=\cos\phi\,\boldsymbol{\hat{x}}+\sin\phi\cos\psi\,\boldsymbol{\hat{y}}+\sin\phi\sin\psi\,\boldsymbol{\hat{z}}$.
We restrict the optimization to $\theta\in[0,2\pi)$, $\phi\in[0,\pi/2]$, and
$\psi\in[0,2\pi)$ to reflect both the modular domain and symmetry. We optimize
two additional variables from Equation 21 — the initial attitude about the
rotation axis $\beta_{0}$, and the constant $\Delta V$, which parameterizes
the flux and albedo. We also optimize the rotational period $p$, which we
restrict to sufficiently encompass all previously measured values. This period
maps the time of observation $t$ to the rotation attitude $\beta$ via
$\beta=2\pi\cdot(t\%p)/p$, a continuous linear mapping from
$(t\in[0,\infty))\rightarrow(\beta\in[0,2\pi))$. Validation of the search
method via fitting to simple rotators is presented in Section 4.4.

The parameter space has many minima because it is highly degenerate and
interdependent — (i) there are multiple well-fitting “fixed-axis” rotations
due to tumbling in these objects, (ii) there is often weak dependence of the
light curve on parameter sets, and (iii) there are distinct parameter sets
which produce relatively similar light curves due to the highly symmetric
equations. Therefore, the optimal parameters depend strongly on the initial
conditions. To account for this, we perform a grid-search optimization over
this data set, using a set of grid-spaced initial values for the parameters
($p$, $\theta$, $\phi$, $\psi$, $\beta_{0}$, and $\Delta V$). For each point
in the grid, curve_fit is used to obtain a parameter set that describes a
local minimum in the parameter space. The final optimized parameters are the
parameter set with the lowest $\chi^{2}$ out of all local minima found with
the initial parameter values.

Note that this model does not incorporate the effects of a changing period,
tumbling, or torques, limiting its applicability for objects with complex
rotational states. However, this model is sufficiently accurate for use in the
SAMUS software, due to the limited effects of the rotational state on the
magnitude of the tidal deformation.

Figure 3: Optimized fit (red line) for a random parameter choice. The true
(blue line) and noise-added data (grey points) are also plotted for
comparison.

Figure 4: The distribution of $\chi^{2}$ for $10^{4}$ synthetic light curves.
The histogram marginalizes over all selected parameters, which were randomly
chosen.

###  4.4 Validation with Simulated Simple Rotators

In this subsection, we validate our use of scipy.optimize.curve_fit to
optimize the rotation axis by sitting a series of simple rotators. We generate
random combinations of parameters $p$, $\theta$, $\phi$, $\psi$, $\beta_{0}$,
and $\Delta V$ which are used to generate synthetic light curves for a 1-day
period (using the ML15 model). We add Gaussian random errors to the synthetic
data, with standard deviations three times the mean of the errors in the
‘Oumuamua photometric data. We also re-scaled the error by the log of the
square root of the magnitude, to imitate photon noise.

We fit the ML15 model to the randomized data using random initial conditions
drawn from Gaussian distributions about the true values. These initial
conditions are unrealistically close to the true values compared to what we
can achieve when optimizing real data. However, the grid-search optimization
performed in Section 4.3 overcomes this drawback and produces high-quality
fits.

In Figure 4, we show the true and noise added data and the optimal fit for a
randomly drawn parameter set. This fit is extremely high quality, and the true
value, optimized value, and statistical difference of the parameters are given
in Table 2. Testing a variety of random parameters provides similar results
(not shown).

Parameters
---
Parameter | True Value | Fit Value | Stat. Diff.
$p$ | 7.2952 | 7.2954 | 0.414
$\theta$ | 5.7832 | 5.7832 | 0.022
$\phi$ | 0.5935 | 0.5916 | 1.418
$\psi$ | 0.2803 | 0.2797 | 1.043
$\beta_{0}$ | 0.3523 | 0.3522 | 0.252
$\Delta V$ | 33.3230 | 33.3265 | 1.402
Table 2: True and derived parameters for the synthetic light curve shown in
Figure 4.

We also perform this test for $10^{4}$ randomly sampled values, and plot their
accuracy versus the estimated values in Figure 4. The majority of fits to
synthetic data are very high-quality with $\chi^{2}\simeq 10^{-2}$.
Importantly, we verified that the distribution of $\chi^{2}$ shown in this
figure is independent of the value of the randomized parameters and of the
parameter itself (not shown for the sake of brevity). Clusters close to
$\chi^{2}\simeq 10^{-2}$ are nearly perfect fits, similar to those presented
in Figure 4. This indicates that this optimization method is effective for any
underlying set of parameters.

##  5 Applicability of SAMUS to 99942 Apophis

As discussed in the introduction, the near-Earth object (NEO) 99942 Apophis
will make a close approach to the Earth in 2029 (Giorgini et al., 2008;
Farnocchia et al., 2013). Nongravitational perturbations from the Yarkovsky
effect could produce an impact event in the 2068 close approach with a
probability $>10^{-6}$ (Farnocchia et al., 2013). The 2029 close approach will
expose the object to significant tidal forces that most likely will not
produce catastrophic disruption but could produce local failures (Scheeres et
al., 2005; Yu et al., 2014). The strength of such forces makes Apophis an
excellent candidate for further investigation with SAMUS. While large-scale
deformation is likely negligible (Yu et al., 2014), the magnitude (or absence)
of deformation in Apophis will provide information on its material properties.
Notably, in 2029, caution must be taken to correct for predicted variations in
the light curve due to tidal resurfacing (Yu et al., 2014; Binzel et al.,
2010; Kim et al., 2023) and alterations in Apophis’ spin state (Scheeres et
al., 2005; Benson et al., 2023). In this section, we discuss the potential
effects of this close approach on Apophis’ shape.

Figure 5: Order-of-magnitude area change of Apophis for a range of dynamic
viscosities, in m2.

At its closest approach, Apophis will be within 38,000 km of the Earth
(Farnocchia et al., 2013). Using Equation 11, we find that Apophis will
experience a maximum tidal acceleration of $4.14\times 10^{-4}$ cm s-2, using
$M=5.97\cdot 10^{27}$ g and $\Delta r=170$ m (Lee et al., 2022).

On closest approach, Apophis will spend approximately 4 hours closer than
50,000 km (https://ssd.jpl.nasa.gov/horizons.cgi). Using the tidal
acceleration computed above, we can use Equation 14 to compute the expected
deformation. The deformation for a range of viscosities is given by Figure 5
assuming a density of $\rho=1.5$ g cm-3, and a radius of $L_{0}\simeq 170$ m
(Lee et al., 2022). While high viscosities show little deformation in the
object, detectable changes are likely possible for a range of low viscosities.
For example, a dynamic viscosity of $\mu\simeq 10^{7}$ poise yields a
deformation of $L\sim 1$ m. This will increase the surface area of Apophis by
approximately $0.5\%$, assuming an ellipsoidal cross section with the vertical
axes. Therefore, low dynamic viscosities could create potentially detectable
effects on the photometric magnitude. However, larger viscosities will not
produce detectable signals.

SAMUS will then be a useful tool to constrain the material properties of
Apophis, using data from the closest approach. This analysis will be
especially relevant in combination with the OSIRIS-APEX mission, which will
encounter Apophis during the asteroid’s 2029 flyby (Nolan et al., 2021). The
in situ measurements, in combination with remote observations, will be useful
for the analysis of Apophis’ composition. The use of SAMUS would be especially
powerful in combination with the analysis of Hirabayashi (2022), who
demonstrated that changes in the rotational state are potential probes of
Apophis’ material properties.

However, Apophis’ rotation state is currently in the process of non-principal
axis (NPA) rotation and is slowly tumbling (Pravec et al., 2014; Benson et
al., 2023). Due to the difficulties of interpreting data from this sort of
light curve, we must develop a sophisticated NPA rotation model for
incorporation into SAMUS. This rotation model will also be invaluable in
accounting for changes in the rotation state due to gravitationally-induced
torques (Scheeres et al., 2005; Benson et al., 2023).

Figure 6: Synthetic light curve using the ML15 model (red line) alongside
photometric observations of ‘Oumuamua (grey points, data from Belton et al.
(2018)). The period, initial rotation state, and average magnitude are
optimized for the first 6 nights.  Figure 7: Flowchart showing the structure
of the simulation runs used to model ‘Oumuamua. Figure 8: Simulated change of
‘Oumuamua’s moment of inertia from 2017 May 4 to 2018 January 16. Red X’s
indicate simulations that were halted because of disintegration or numerical
divergence. Empty squares indicate parameters for simulations that were not
run, as convergence was achieved for lower dynamic viscosity values.

##  6 ‘Oumuamua Simulation

In this section, we use SAMUS to constrain the size and dynamic viscosity of
‘Oumuamua, demonstrating the effectiveness of this software for investigating
objects which are potential candidates for tidal deformation. We first use the
methodology described in Section 4 to obtain a characteristic “fixed-axis” NPA
rotation for ‘Oumuamua, which we then incorporate into a SAMUS model for the
simulation of tidal deformation. It is important to note that the rotation
model found for ‘Oumuamua is non-physical. Properly modeling NPA rotation is
necessary for a full understanding of the case of ‘Oumuamua, which we leave
for future research. In these models, we adopt an aspect ratio of 6:6:1,
following Mashchenko (2019). We also make the approximation that the granular
constituents of ‘Oumuamua are sufficiently small, such that the bulk aggregate
behaves as a fluid.

###  6.1 ‘Oumuamua Axis Fitting

In order to find the characteristic rotation axis, we use the methodology
described above to fit a “fixed-axis” NPA rotation to the photometric data
collected in Belton et al. (2018). These data are corrected for light travel
time, helio- and geo-centric distance, Solar magnitude, and filter color
(Belton et al., 2018), and the phase angle data are taken from JPL’s Horizons
database (https://ssd.jpl.nasa.gov/horizons.cgi). The optimal parameters (from
a grid search over 3,000 points) are given in Table 3. Although this model
does not accurately capture the physics of ‘Oumuamua’s rotation, it provides a
basic light curve for comparison to modeled deformation.

OPTIMAL PARAMETERS
---
Parameter: | Value:
$p$ [h] | $7.3975\pm 1.11\cdot 10^{-3}$
$\theta$ | $4.6995\pm 1.69\cdot 10^{-2}$
$\phi$ | $1.4193\pm 1.32\cdot 10^{-3}$
$\psi$ | $1.2545\pm 3.77\cdot 10^{-3}$
$\beta_{0}$ | $4.3950\pm 9.82\cdot 10^{-3}$
$\Delta V$ | $32.1560\pm 1.26\cdot 10^{-2}$
Table 3: The optimal parameters for the ML15 model in comparison to the
photometric data. Note that the errors here are not for the grid-search
optimization but for the single, maximal optimization, as the grid-search
uncertainty is unknown, and obtaining the uncertainty through Markov Chain
Monte Carlo methods is computationally infeasible for this problem.

These parameters correspond to a rotation axis of
$\boldsymbol{\hat{\Omega}}=0.1509\boldsymbol{\hat{x}}+0.3075\boldsymbol{\hat{y}}+0.9395\boldsymbol{\hat{z}}$.
This is relatively consistent with physical expectations, which suggests that
the rotation axis should approach the smallest principal axis over time due to
relative rotational energies.

We show the synthetic light curve for these optimal parameters in Figure 6.
The fit is qualitatively accurate for the October nights, but matches the
November data poorly. This effect may be a result of a secular change in the
spin period. However, proper evaluation of this effect requires the
computation of outgassing-induced rotational dynamics, which are addressed in
an accompanying paper. While the lack of tumbling in the model causes the fit
to be non-exact, even in October, this does not have a significant effect on
the tidal deformation effects addressed in this paper.

###  6.2 ‘Oumuamua Tidal Deformation

We adopt an aspect ratio of 6:6:1 (Mashchenko, 2019) and a bulk density of
$\rho=0.5$ g cm-3, which is typical for Solar System comets (Britt et al.,
2006). The viscosity is a parameter in these experiments, but is assumed to be
homogeneous within ‘Oumuamua. We run simulations with an initial primary axis
of $a=35,55,75,95,115$, and 150 meters in diameter, using the characteristic
“fixed-axis” NPA rotation described in Section 6.1. The rotation axis is
$\boldsymbol{\hat{\Omega}}=0.1509\boldsymbol{\hat{x}}+0.3075\boldsymbol{\hat{y}}+0.9395\boldsymbol{\hat{z}}$,
and we assume a constant rotational period of $p=7.3975$ hours. While these
are held constant in each simulation, if the simulated body undergoes
significant deformation, conservation of angular momentum will produce a
change in the spin period. The simulations are run from 2017 May 4 to 2018
January 16, using trajectory data obtained from the JPL Horizons database.

We also ran simulations with rotations about each principal axis, and verified
that this does not significantly affect the results, validating our use of the
characteristic rotation axis and our conclusions. These simulations are not
shown, but their results are available on GitHub.

While the dynamic viscosity depends on the temperature, neither are fully
known for ‘Oumuamua. Therefore, we simplify the problem by performing these
numerical experiments for a variety of effective viscosities, ignoring the
temperature dependence. We initialize each simulation with a dynamic viscosity
of $\mu=10^{5}$ g cm-1 s1, approximately equal to that of peanut butter under
high pressure (Citerne et al., 2001). The simulation is halted if the body is
distorted to more than twice its initial size, as the cometary materials
should disintegrate when subjected to such significant force (although the
shear tolerances of cometary materials are unknown). The simulation is reset
and the viscosity increased by an order of magnitude if (i) the solvers for
the Navier-Stokes equations fail to converge, indicating a divergence beyond
physical conditions, (ii) the simulation was halted due to a non-physical size
increase, or (iii) the moment of inertia changes by more than 1% over the
path. These simulations are run with a trajectory jump tolerance of 1% (see
Section 2.2), and $C_{\text{max}}=1$. We compute 10 time steps over each
rotational period, which fully samples the force over the rotation. For
clarity, a flowchart describing this structure is shown in Figure 7, and the
code used to create these simulations is available on GitHub.

Figure 9: Evolution of ‘Oumuamua’s spin period due to modulation in the moment
of inertia. The period is initialized at $p=7.3906$, ensuring that $p=7.3975$
hours at discovery (vertical dashed line) for the optimal fit. Rows represent
different initial sizes.  Figure 10: Evolution of ‘Oumuamua’s aspect ratios
due to tidal deformation. The primary aspect ratio (left column) is the ratio
of the largest to the smallest axis, while the secondary (right column) is the
ratio of the intermediate to the smallest axis. At $t=0$, both aspect ratios
are 6. The time of discovery is plotted as a vertical dashed line. Figure 11:
Synthetic light curves of ‘Oumuamua which incorporate simulated changes in
period and aspect ratio from tidal deformation. Rows correspond to the initial
sizes and colors correspond to the dynamic viscosities of the body in the
simulations. The time of discovery is indicated with a vertical dashed line.

In Figure 8, we present an overview of all of the simulations. The moment of
inertia is calculated as $I(t)=\int\rho\boldsymbol{r}(t)^{2}\text{d}x$, and we
show $I(t_{\text{max}})/I(0)$. These changes in the moment of inertia will
also affect the spin period.

To quantify this effect, we assume an idealized scenario in which the spin
period only changes in response to tidal deformation and subsequent changes in
the moment of inertia. This is not realistic, because outgassing torques
should dominate the rotational state. However, we are only considering the
tidal deformation of this object in the absence of outgassing. We assume that
the angular momentum is conserved such that $L(0)=L(t)$, where
$L(t)=I(t)\omega(t)$ and $\omega(t)$ is the angular frequency. We additionally
set $p(0)=7.3906$ hours, chosen such that the object has $p=7.3975$ hours at
the time of detection, for the simulation which produced an optimal fit to the
photometric data. The evolution of the spin period is shown in Figure 9 and
the aspect ratios in Figure 10. For high viscosities, the aspect ratios and
the spin period are approximately constant, as the viscous forces are much
stronger than the tidal forces. The rotational period increases in every case,
since the rotation about the z-axis forces the tidal deformation to increase
the moment of inertia. Although these simulations are highly idealized and do
not incorporate outgassing torques, it is clear that tidal deformation can
change the rotational state and light curve amplitude, based on the viscosity.

$\chi^{2}$ of SIMULATED LIGHT CURVES vs DATA
---
Size [m] | Dynamic Viscosity [g cm-1 s-1]
| $10^{5}$ | $10^{6}$ | $10^{7}$ | $10^{8}$ | $10^{9}$ | $\infty$
35 | - | 5.038e+04 | 3.882e+04 | 3.833e+04 | - | 3.831e+04
55 | - | 1.078e+05 | 3.981e+04 | 3.840e+04 | - | 3.835e+04
75 | - | - | 4.175e+04 | 3.851e+04 | 3.830e+04 | 3.829e+04
95 | - | - | 4.539e+04 | 3.867e+04 | 3.832e+04 | 3.830e+04
115 | - | - | 5.217e+04 | 3.886e+04 | 3.833e+04 | 3.831e+04
150 | - | - | 8.213e+04 | 3.935e+04 | 3.837e+04 | 3.833e+04
Table 4:  $\chi^{2}$ values for each synthetic light curve fit to the
photometric data. ’-’s indicate simulations which were numerically unstable.

Overall, objects with larger sizes undergo significantly more deformation and
require larger viscosities to maintain stability. A dynamic viscosity of
$\mu=10^{7}$ g cm-1 s-1 is sufficient for all tested sizes to maintain
physical conditions over the trajectory. For initial sizes of 35 and 55
meters, $\mu=10^{8}$ g cm-1 s-1 lead to only small-scale changes in the body
over its trajectory. For the remaining sizes (75, 95, 115, and 150 meters),
$\mu=10^{9}$ g cm-1 s-1 similarly allows only minor changes.

###  6.3 Synthetic Light Curves Incorporating Tidal Effects

In this section, we present synthetic light curves for ‘Oumuamua over 5 1-day
periods which incorporate the effects of tidal deformation. These 5 1-day
periods begin on 2017 May 4, 2017 August 9, 2017 September 15 (perihelion),
2017 October 24 (detection), and 2018 January 29 (final simulated point). As
before, trajectory and phase angle data are obtained from JPL’s Horizons
database.

The light curves are computed with Equation 21 (the ML15 model) assuming
optimal “fixed-axis” NPA rotation parameters (Section 4.3). The ML15 model
incorporates a Lommel-Seeliger scattering surface for an ellipsoidal body,
phase angle effects, and an arbitrary orientation and rotation. We compute the
period evolution using the time evolution of the aspect ratio and moment of
inertia (Figures 9 and 10) computed in Section 6. This period is used to
compute the orientation, $\beta$, as $\beta=2\pi\,(t\%p)/p$, with ‘%’
representing the modulo. The ML15 model also incorporates the instantaneous
aspect ratio into the light curve, using simulated data provided by SAMUS. The
final synthetic light curves that incorporate amplitude and period modulations
are presented in Figure 11. For comparison, we also show light curves with
constant period and aspect ratio, which correspond to the dynamic viscosity
$\mu\rightarrow\infty$.

Each model is initialized on 2017 May 4 with period $p=7.3975$ and aspect
ratio 6:6:1. The parameters for the observation vector and the “fixed-axis”
NPA rotation are the optimal results from Section 6.1. As these simulations
use initial sizes and rotation periods which were obtained from data measured
in October 2017, these light curves do not match up well with available
photometric data, but instead provide qualitative examples of an evolving
light curve due to tidal effects.

These light curves are corrected for brightness variation due to ‘Oumuamua’s
helio- and geo-centric distances. Therefore, amplitude variations are due to
modulations of the aspect ratio and/or phase angle. Larger objects have
brighter average magnitudes because the albedo is constant in all simulations.
There are curious features in the light curves in August, which are due to
changes in phase angle (Section 4.2).

Aside from the effects of phase angle and period modulations, there is little
variation in the shape of these light curves. Additionally, there is a notable
decrease in amplitude post-perihelion. While the aspect ratio does increase
during this time period, this feature in the light curve is entirely due to
the phase angle. This effect is most obvious for the low-viscosity
simulations. This implies that the evolving aspect ratio does not produce an
observable signature in the amplitude for these simulation parameters. While
tidal deformation can significantly alter the spin period, its effect on the
amplitude of the light curve is not detectable for ‘Oumuamua. Period and
amplitude changes in light curves of future interstellar objects may result
from tidal deformation, and could be used to constrain the dynamic viscosity,
and potentially the material composition.

Figure 12: Best fit $\chi^{2}$ values between synthetic and photometric data.
Empty spaces indicate parameter choices for simulations that were not run.
Fits for both October and November are both computed and the $\chi^{2}$ values
are added together. Figure 13: Synthetic light curves for ‘Oumuamua which
incorporate simulated aspect ratio changes and constant period. Colored lines
indicate the dynamic viscosity and photometric data are shown in grey points.
These curves use optimized values of $\beta_{0}$ and $\Delta V$. Note that the
rotation state is entirely fixed in these curves, with no evolution in the
rotation period or axis. This allows for comparison to the photometric data
and the original composite light curve, since the evolution of the rotation
will require a more complex and physically accurate model.

We also compare the synthetic light curves with the photometric data, using
the optimal period of $p=7.3975$ hours. As in Section 4, we optimize the
initial rotation state $\beta_{0}$ and constant $\Delta V$ with
scipy.optimize.curve_fit, although we keep the remaining parameters constant.
By keeping the rotation state constant, we are able to compare these light
curves to the optimal model computed in Section 6.1, despite the non-physical
rotation model. The November data is optimized independently, assuming
$p_{\text{Nov}}=7.1910$, which is found by a separate minimization. In order
to evaluate the validity of each fit, the $\chi^{2}$ values for both months
are added together, which we present for each simulation in Table 4 and Figure
12. The resulting synthetic light curves along with the photometric data are
shown in Figure 13 for October nights.

Figure 14: The best-fit synthetic light curve (red line) and photometric data
(grey points). This model incorporates the simulated changes in aspect ratio.
Here $\mu=10^{9}$ g cm-1 s-1 and the body is initially 75 meters.

A dynamic viscosity of $10^{9}$ g cm-1 s-1 and initial size of 75:75:12.5
meters produce optimal fits (discounting $\mu=\infty$), and this best-fit
synthetic light curve is plotted along with the photometric data in Figure 14.
It is worth noting that higher viscosity always produces a better fit for
every initial size, therefore, no deformation provides the optimal explanation
for ‘Oumuamua’s light curve.

##  7 Discussion

In this paper, we presented a novel open-source software (SAMUS) which models
the deformation of ellipsoidal minor bodies under tidal stresses. This code is
capable of constraining the material viscosity of such objects, as this
deformation can produce observable changes in the rotation state and the
photometric light curve. As a demonstration of the functionality of SAMUS, we
applied it to the interstellar object 1I/‘Oumuamua and investigated the
material viscosity. Our simulations show that tidal deformation would not
cause measurable amplitude variations for ‘Oumuamua. It is possible that tidal
deformation could be detected in other objects, depending on their bulk
properties and albedo. The non-detection of tidal deformation for ‘Oumuamua
indicates that tidal forces were not a significant factor operating within its
body, likely due to the rigidity/viscosity of this object.

We assume in these simulations and in SAMUS that the subject (‘Oumuamua) has a
constant dynamic viscosity. However, ‘Oumuamua was tumbling and receiving
solar irradiation across the surface. Fitzsimmons et al. (2018) modelled the
radial thermal profile of ‘Oumuamua, assuming a thermal conductivity of
$k=100$ erg s-1 cm-1 K-1, a heat capacity of $C=5.5\cdot 10^{6}$ erg g-1, and
a density of $\rho=1$ g cm-3 (see Figure 4 of that publication). Those authors
found that although the surface of the body reached a maximum of $T\simeq
600$K at perihelion, subsurface layers at depths $\geq$30 cm never reached a
temperature greater than $T\simeq 150$K. Therefore, it is feasible that there
was a strong temperature gradient immediately below the surface, which was not
present in the cometary core. While the assumption of constant dynamic
viscosity is less applicable for the surface half-meter, the majority of the
body has a minimal temperature gradient and likely a nearly-constant dynamic
viscosity.

We used SAMUS to constrain the effective dynamic viscosity and size of
‘Oumuamua. For the optimal 115:111:19 meter size found by Mashchenko (2019),
the dynamic viscosity must be $\mu\geq 10^{9}$ g cm-1 s-1, roughly equivalent
to the viscosity of bitumen pitch ($\mu\simeq 3\cdot 10^{9}$ g cm-1 s-1)
(Edgeworth et al., 1984).444This experiment measured pitch over temperatures
ranging from 9∘C to 30∘C. Allowable dynamic viscosities are provided for
different initial sizes in Table 5. ‘Stability’ here describes cases with
$<1\%$ change in the moment of inertia over the trajectory. These are
extremely high viscosities — for comparison, terrestrial fluids such as water
or olive oil have viscosities of $\mu\simeq 10^{-2}$ and $\mu\simeq 1$ g cm-1
s-1, respectively (Rumble, 2022; Fellows, 2022). However, the allowable
viscosity range is compatible with a variety of terrestrial materials,
including water ice ($\mu\simeq 10^{13}$ g cm-1 s-1) (Fowler, 1997) and the
terrestrial mantle ($\mu\simeq 10^{22}$ g cm-1 s-1), although the temperature
dependence of dynamic viscosity makes these direct comparisons difficult.

Unfortunately, these constraints on material composition are weakened by the
difficulty of measuring the viscosity of proposed materials, especially at the
relevant temperatures. While water ice has a known dynamic viscosity of
$\mu\simeq 10^{13}$ g cm-1 s-1 (over long-term high-latitude terrestrial
temperatures) (Fowler, 1997), the viscosities of exotic solids such as CO, H2,
and N2 ice are far more difficult to measure. Vilella & Deschamps (2017)
analyzed the dynamics of N2 ice glaciers in the Sputnik Planitia of Pluto, and
found that the observed polygonal structure is consistent with viscosities of
$10^{15}$–$10^{17}$ g cm-1 s-1, although these values are not well-constrained
and are highly uncertain. Due to experimental barriers, however, no
measurements are available for the viscosity of H2 or CO ice. Additionally,
the temperature dependence of viscosity renders the extension of such
measurements to vacuum difficult, as measurements at such temperatures are not
generally available. Finally, the structural composition of and material
mixing within the body further complicate the viscosity measurements. Despite
these practical difficulties, constraining the dynamic viscosity provides a
new methodology to potentially test proposed structures and compositions for
‘Oumuamua, future interlopers of this type, and Solar System minor bodies.

MINIMAL STABLE VISCOSITY
---
Primary Axis [m] | Viscosity [g cm-1 s-1]
35 | $10^{8}$
55 | $10^{8}$
75 | $10^{9}$
115 | $10^{9}$
150 | $10^{9}$
Table 5: Stable dynamic viscosities for initial sizes of ‘Oumuamua.

The second-largest dynamic viscosity remains relatively stable in all cases,
with small changes in moment of inertia slightly larger than the $<1\%$
simulation cutoff. These viscosities allow for changing period and amplitude,
while preventing nonphysical divergence of the body. It should be noted,
however, that the viscosity of $10^{9}$ g cm-1 s-1 is orders of magnitude
lower than the $10^{13}$–g cm-1 s-1 viscosity of water ice. Therefore, this
constraint provides little differentiation between various solid ices.

It is also worth noting that changes in the moment of inertia can arise from
sources other than tidal deformation, including ice crystallization, nucleus
size changes, and dense mantle formation (Watanabe, 1992). Like in tidal
deformation, these changes can cause the rotation state to detectably evolve.
However, many of the effects will cause a decrease in the moment of inertia,
in contrast to the increase that tidal deformation generally causes.

We have also demonstrated that for these models, an initial size of 75 meters
and a dynamic viscosity of $\mu\gg 10^{9}$ g cm-1 s-1 (no deformation)
provides an optimal match to the photometric data, although different object
sizes fit nearly equally well, and are relatively unconstrained by this model.
These results indicate that tidal deformation likely played little-to-no role
in the rotation state or the shape of ‘Oumuamua along its Solar System
trajectory. It is possible that the change in period detected by Flekkøy et
al. (2019) is due to non-tidal torques. While this result is not strongly
dependent on the rotation (see GitHub for simulations around the principal
axes), a fully physical NPA rotation model would be needed to confirm this
result.

‘Oumuamua left many unanswered questions as it exited the Solar System, and
despite intense scrutiny, there is still no general consensus regarding the
provenance of the object. The discovery implies a spatial number density of
similar objects of order $n_{o}\sim 1-2\times 10^{-1}\,$au-3 (Trilling et al.,
2017; Laughlin & Batygin, 2017; Jewitt et al., 2017; Moro-Martín, 2018; Zwart
et al., 2018; Do et al., 2018; Moro-Martín, 2019). Detection and
characterization of future interstellar objects offer the most promising
avenue for resolving these questions.

The forthcoming Rubin Observatory Legacy Survey of Space and Time (LSST)
(Jones et al., 2009; Ivezić et al., 2019) will effectively detect such
transient objects (Solontoi et al., 2011; Vereš & Chesley, 2017a, b; Jones et
al., 2018). The survey should detect $\sim$1 ‘Oumuamua-like interstellar
object every year (Moro-Martín et al., 2009; Engelhardt et al., 2017; Cook et
al., 2016; Trilling et al., 2017; Seligman & Laughlin, 2018; Hoover et al.,
2022; Marčeta, 2023). In addition, the forthcoming NEO Surveyor (Mainzer et
al., 2015) may also detect interstellar objects, and could offer information
about outgassing sources via its infrared capabilities. Space based in-situ
measurements of an interstellar object would provide valuable information
regarding the composition and bulk properties (Hein et al., 2017; Seligman &
Laughlin, 2018; Meech et al., 2019; Castillo-Rogez et al., 2019; Jones & ESA
Comet Interceptor Team, 2019; Hibberd et al., 2020; Donitz et al., 2021; Pau
Sánchez et al., 2021; Meech et al., 2021; Hibberd et al., 2022; Moore et al.,
2021a, b). Additionally, future observations of the amplitude variations (if
any) of interstellar objects may be used to constrain their viscosities — and
may be able to constrain composition — via the techniques developed in this
paper.

In the future, similar analyses of tidal deformation with SAMUS may be useful
for other small bodies. For objects with high rotation rates and closer solar
approaches, tidal and centrifugal forces will have a larger effect. An example
is 3200 Phaethon which also exhibits unexplained activity (Jewitt & Li, 2010;
Jewitt et al., 2013; Li & Jewitt, 2013; Hui & Li, 2017) and is a target for
the Japan Aerospace Exploration Agency (JAXA) DESTINY+ mission in situ mission
(Arai et al., 2021). For objects like 3200 Phaethon, tidal deformation may be
detectable in photometric data, and would enable more strenuous constraints on
the bulk physical properties of those objects. The recent detection of ‘dark
comets’ (Chesley et al., 2016; Farnocchia et al., 2022; Seligman et al., 2023)
and their still-unknown provenance provides another class of objects which may
benefit from the application of SAMUS and other techniques developed in this
paper.

The 2029 near-Earth flyby of 99942 Apophis is also an excellent candidate for
the application of SAMUS. We demonstrate that the significant tidal forces
experienced by Apophis may lead to deformation in its shape, which could be
detectable either by photometric observations, radar observations, or the
OSIRIS-APEX mission. Using SAMUS, the magnitude of such deformation would
provide constraints on the internal material properties of Apophis.

However, it is critical to note that the rotation model used in SAMUS is non-
physical “fixed-axis” NPA rotation. Because of the tumbling rotation states of
‘Oumuamua and Apophis, a fully physical model will be necessary to fully
address the evolution of ‘Oumuamua and to apply SAMUS to Apophis in the coming
years.

The implementation of additional factors into SAMUS is also a worthwhile
subject for future work. Particularly relevant is the addition of a shifting
and strictly-conserved angular momentum axis, under the influence of external
torques and tumbling. This addition will enable more accurate simulations of
the deformation and rotation of such objects and thereby more accurate light
curve fitting and analysis. The addition outgassing torques and ablation to
SAMUS, and a physical rotation model will further enable effective analysis of
rotation states and material properties of minor bodies.

##  8 Acknowledgements

We thank Faith Vilas for useful advice regarding the scientific content and
structure of this manuscript. We thank the two anonymous reviewers for
insightful and helpful suggestions which greatly strengthened the scientific
content of this manuscript. This research utilized the University of Chicago’s
Research Computing Center for numerical calculations. We thank Adina Feinstein
for useful conversations and suggestions. We thank David Jewitt for useful
feedback on the originally submitted version of the manuscript. DZS
acknowledges financial support from the National Science Foundation Grant No.
AST-17152, NASA Grant No. 80NSSC19K0444 and NASA Contract NNX17AL71A from the
NASA Goddard Space Flight Center.

## References

  * A’Hearn et al. (2012) A’Hearn, M. F., Feaga, L. M., Keller, H. U., et al. 2012, ApJ, 758, 29, doi: 10.1088/0004-637X/758/1/29
  * Almeida-Fernandes & Rocha-Pinto (2018) Almeida-Fernandes, F., & Rocha-Pinto, H. J. 2018, MNRAS, 480, 4903, doi: 10.1093/mnras/sty2202
  * Alnaes et al. (2014) Alnaes, M. S., Logg, A., Olgaard, K. B., Rognes, M. E., & Wells, G. N. 2014, ACM Transactions on Mathematical Software, 40, doi: 10.1145/2566630
  * Alnaes et al. (2015) Alnaes, M. S., Blechta, J., Hake, J., et al. 2015, Archive of Numerical Software, 3, doi: 10.11588/ans.2015.100.20553
  * Arai et al. (2021) Arai, T., Yoshida, F., Kobayashi, M., et al. 2021, in 52nd Lunar and Planetary Science Conference, Lunar and Planetary Science Conference, 1896
  * Bagheri et al. (2022) Bagheri, A., Efroimsky, M., Castillo-Rogez, J., et al. 2022, Advances in Geophysics, 63, 231, doi: 10.1016/bs.agph.2022.07.004
  * Bannister et al. (2017) Bannister, M. T., Schwamb, M. E., Fraser, W. C., et al. 2017, The Astrophysical Journal, 851, doi: 10.3847/2041-8213/aaa07c
  * Belton et al. (2018) Belton, M. J., Hainaut, O. R., Meech, K. J., et al. 2018, The Astrophysical Journal, 856, doi: 10.3847/2041-8213/aab370
  * Benson et al. (2023) Benson, C. J., Scheeres, D. J., Brozović, M., et al. 2023, Icarus, 390, 115324, doi: 10.1016/j.icarus.2022.115324
  * Bergner & Seligman (2023) Bergner, J. B., & Seligman, D. Z. 2023, Nature, 615, 610, doi: 10.1038/s41586-022-05687-w
  * Binzel et al. (2010) Binzel, R. P., Morbidelli, A., Merouane, S., et al. 2010, Nature, 463, 331
  * Biver & Bockelée-Morvan (2016) Biver, N., & Bockelée-Morvan, D. 2016, IAU Focus Meeting, 29A, 228, doi: 10.1017/S1743921316002945
  * Black et al. (1995) Black, G. J., Nicholson, P. D., & Thomas, P. C. 1995, Icarus, 117, 149, doi: 10.1006/icar.1995.1148
  * Bockelée-Morvan & Biver (2017) Bockelée-Morvan, D., & Biver, N. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160252, doi: 10.1098/rsta.2016.0252
  * Bolin et al. (2017) Bolin, B. T., Weaver, H. A., Fernandez, Y. R., et al. 2017, The Astrophysical Journal, 852, doi: 10.3847/2041-8213/aaa0c9
  * Bottke & Melosh (1996a) Bottke, William F., J., & Melosh, H. J. 1996a, Icarus, 124, 372, doi: 10.1006/icar.1996.0215
  * Bottke & Melosh (1996b) Bottke, W. F., & Melosh, H. J. 1996b, Nature, 381, 51, doi: 10.1038/381051a0
  * Boyle et al. (2022) Boyle, M., Long, J., Ling, M., et al. 2022, moble/quaternion: Release v2022.4.2, v2022.4.2, Zenodo, doi: 10.5281/zenodo.6499564
  * Britt et al. (2006) Britt, D. T., Consolmagno, G. J., & Merline, W. J. 2006, in 37th Annual Lunar and Planetary Science Conference, ed. S. Mackwell & E. Stansbery, Lunar and Planetary Science Conference, 2214
  * Brozović et al. (2018) Brozović, M., Benner, L. A., McMichael, J. G., et al. 2018, Icarus, 300, 115, doi: https://doi.org/10.1016/j.icarus.2017.08.032
  * Castillo-Rogez et al. (2019) Castillo-Rogez, J., Landau, D., Chung, S.-J., & Meech, K. 2019, in Approach to exploring interstellar objects and long-period comets
  * Chesley et al. (2016) Chesley, S. R., Farnocchia, D., Pravec, P., & Vokrouhlický, D. 2016, in Asteroids: New Observations, New Models, ed. S. R. Chesley, A. Morbidelli, R. Jedicke, & D. Farnocchia, Vol. 318, 250–258, doi: 10.1017/S1743921315008790
  * Citerne et al. (2001) Citerne, G. P., Carreau, P. J., & Moan, M. 2001, Rheol. Acta, 40, 86, doi: 10.1007/s003970000120
  * Cochran et al. (2015) Cochran, A. L., Levasseur-Regourd, A.-C., Cordiner, M., et al. 2015, Space Sci. Rev., 197, 9, doi: 10.1007/s11214-015-0183-6
  * Colombo (1965) Colombo, G. 1965, Nature, 208, 575, doi: 10.1038/208575a0
  * Cook et al. (2003) Cook, C. M., Melosh, H. J., & Bottke, W. F. 2003, Icarus, 165, 90, doi: 10.1016/S0019-1035(03)00177-5
  * Cook et al. (2016) Cook, N. V., Ragozzine, D., Granvik, M., & Stephens, D. C. 2016, ApJ, 825, 51, doi: 10.3847/0004-637X/825/1/51
  * Courant et al. (1967) Courant, R., Friedrichs, K., & Lewy, H. 1967, IBM Journal of Research and Development, 11, 215, doi: 10.1147/rd.112.0215
  * Dalcín & Fang (2021) Dalcín, L., & Fang, Y.-L. L. 2021, Computing in Science & Engineering, 23, 47, doi: 10.1109/MCSE.2021.3083216
  * Dalcín et al. (2005) Dalcín, L., Paz, R., & Storti, M. 2005, Journal of Parallel and Distributed Computing, 65, 1108, doi: https://doi.org/10.1016/j.jpdc.2005.03.010
  * Dalcín et al. (2008) Dalcín, L., Paz, R., Storti, M., & D’Elía, J. 2008, Journal of Parallel and Distributed Computing, 68, 655, doi: https://doi.org/10.1016/j.jpdc.2007.09.005
  * Dalcín et al. (2011) Dalcín, L. D., Paz, R. R., Kler, P. A., & Cosimo, A. 2011, Advances in Water Resources, 34, 1124, doi: https://doi.org/10.1016/j.advwatres.2011.04.013
  * Darwin (1879) Darwin, G. H. 1879, The Observatory, 3, 79
  * Darwin (1880) —. 1880, Nature, 21, 235, doi: 10.1038/021235a0
  * DeMartini et al. (2019) DeMartini, J. V., Richardson, D. C., Barnouin, O. S., et al. 2019, Icarus, 328, 93, doi: https://doi.org/10.1016/j.icarus.2019.03.015
  * Desch & Jackson (2021) Desch, S. J., & Jackson, A. P. 2021, Journal of Geophysical Research: Planets, 126, e2020JE006807, doi: https://doi.org/10.1029/2020JE006807
  * Dinsmore & de Wit (2022) Dinsmore, J. T., & de Wit, J. 2022, Monthly Notices of the Royal Astronomical Society, 520, 3459, doi: 10.1093/mnras/stac2866
  * Do et al. (2018) Do, A., Tucker, M. A., & Tonry, J. 2018, ApJ, 855, L10, doi: 10.3847/2041-8213/aaae67
  * Donitz et al. (2021) Donitz, B. P., Castillo-Rogez, J. C., & Matousek, S. E. 2021, in 2021 IEEE Aerospace Conference (50100), 1–12, doi: 10.1109/AERO50100.2021.9438142
  * Drahus et al. (2017) Drahus, M., Guzik, P., Waniak, W., et al. 2017, Nature Astronomy, 2, 407–412, doi: 10.1038/s41550-018-0440-1
  * Edgeworth et al. (1984) Edgeworth, R., Dalton, B. J., & Parnell, T. 1984, European Journal of Physics, 5, 198, doi: 10.1088/0143-0807/5/4/003
  * Engelhardt et al. (2017) Engelhardt, T., Jedicke, R., Vereš, P., et al. 2017, AJ, 153, 133, doi: 10.3847/1538-3881/aa5c8a
  * Farnocchia et al. (2013) Farnocchia, D., Chesley, S. R., Chodas, P. W., et al. 2013, Icarus, 224, 192, doi: 10.1016/j.icarus.2013.02.020
  * Farnocchia et al. (2022) Farnocchia, D., Seligman, D. Z., Granvik, M., et al. 2022, arXiv e-prints, arXiv:2212.08135, doi: 10.48550/arXiv.2212.08135
  * Fellows (2022) Fellows, P. 2022, Food Processing Technology: Principles and practice (Woodhead Publishing)
  * Feng & Jones (2018) Feng, F., & Jones, H. R. A. 2018, ApJ, 852, L27, doi: 10.3847/2041-8213/aaa404
  * Fitzsimmons et al. (2018) Fitzsimmons, A., Snodgrass, C., Rozitis, B., et al. 2018, Nature Astronomy, 2, 133, doi: 10.1038/s41550-017-0361-4
  * Flekkøy et al. (2019) Flekkøy, E. G., Luu, J., & Toussaint, R. 2019, The Astrophysical Journal, 885, doi: 10.3847/2041-8213/ab4f78
  * Fowler (1997) Fowler, A. C. 1997, The Mathematics of Models for Climatology and Environment, 301–336, doi: 10.1007/978-3-642-60603-8_9
  * Fraser et al. (2018) Fraser, W. C., Pravec, P., Fitzsimmons, A., et al. 2018, Nature Astronomy, 2, 383–386, doi: 10.1038/s41550-018-0398-z
  * Füglistaler & Pfenniger (2018) Füglistaler, A., & Pfenniger, D. 2018, A&A, 613, A64
  * Fuller et al. (2016) Fuller, J., Luan, J., & Quataert, E. 2016, MNRAS, 458, 3867, doi: 10.1093/mnras/stw609
  * Gaidos et al. (2017) Gaidos, E., Williams, J., & Kraus, A. 2017, RNAAS, 1, 13, doi: 10.3847/2515-5172/aa9851
  * Geuzaine & Remacle (2020) Geuzaine, C., & Remacle, J.-F. 2020, Gmsh, 4.6.0.  http://http://gmsh.info/
  * Giorgini et al. (2008) Giorgini, J. D., Benner, L. A. M., Ostro, S. J., Nolan, M. C., & Busch, M. W. 2008, Icarus, 193, 1, doi: 10.1016/j.icarus.2007.09.012
  * Goldreich & Peale (1966) Goldreich, P., & Peale, S. 1966, AJ, 71, 425, doi: 10.1086/109947
  * Goldreich & Peale (1968) Goldreich, P., & Peale, S. J. 1968, ARA&A, 6, 287, doi: 10.1146/annurev.aa.06.090168.001443
  * Hallatt & Wiegert (2020) Hallatt, T., & Wiegert, P. 2020, The Astronomical Journal, 159, 147, doi: 10.3847/1538-3881/ab7336
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hein et al. (2017) Hein, A. M., Perakis, N., Eubanks, T. M., et al. 2017, arXiv e-prints, arXiv:1711.03155.  https://arxiv.org/abs/1711.03155
  * Hibberd et al. (2022) Hibberd, A., Hein, A., Eubanks, M., & Kennedy, Robert, I. 2022, arXiv e-prints, arXiv:2201.04240.  https://arxiv.org/abs/2201.04240
  * Hibberd et al. (2020) Hibberd, A., Hein, A. M., & Eubanks, T. M. 2020, Acta Astronautica, 170, 136, doi: 10.1016/j.actaastro.2020.01.018
  * Hirabayashi (2014) Hirabayashi, M. 2014, Icarus, 236, 178, doi: https://doi.org/10.1016/j.icarus.2014.02.024
  * Hirabayashi (2015) —. 2015, Monthly Notices of the Royal Astronomical Society, 454, 2249, doi: 10.1093/mnras/stv2017
  * Hirabayashi (2022) Hirabayashi, M. 2022, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 54, AAS/Division for Planetary Sciences Meeting Abstracts, 115.06
  * Hirabayashi & Scheeres (2019) Hirabayashi, M., & Scheeres, D. J. 2019, Icarus, 317, 354, doi: https://doi.org/10.1016/j.icarus.2018.08.003
  * Hoang & Loeb (2020) Hoang, T., & Loeb, A. 2020, The Astrophysical Journal Letters, 899, L23, doi: 10.3847/2041-8213/abab0c
  * Hoover et al. (2022) Hoover, D. J., Seligman, D. Z., & Payne, M. J. 2022, , 3, 71, doi: 10.3847/PSJ/ac58fe
  * Hsieh et al. (2021) Hsieh, C.-H., Laughlin, G., & Arce, H. G. 2021, ApJ, 917, 20, doi: 10.3847/1538-4357/ac0729
  * Hui & Knight (2019) Hui, M.-T., & Knight, M. M. 2019, The Astronomical Journal, 158, 256, doi: 10.3847/1538-3881/ab50b8
  * Hui & Li (2017) Hui, M.-T., & Li, J. 2017, AJ, 153, 23, doi: 10.3847/1538-3881/153/1/23
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jackson & Desch (2021) Jackson, A. P., & Desch, S. J. 2021, Journal of Geophysical Research: Planets, 126, doi: 10.1029/2020je006706
  * Jewitt & Li (2010) Jewitt, D., & Li, J. 2010, AJ, 140, 1519, doi: 10.1088/0004-6256/140/5/1519
  * Jewitt et al. (2013) Jewitt, D., Li, J., & Agarwal, J. 2013, ApJ, 771, L36, doi: 10.1088/2041-8205/771/2/L36
  * Jewitt et al. (2017) Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, The Astrophysical Journal, 850, doi: 10.3847/2041-8213/aa9b2f
  * Jewitt & Seligman (2022) Jewitt, D., & Seligman, D. Z. 2022, arXiv e-prints, arXiv:2209.08182.  https://arxiv.org/abs/2209.08182
  * Jones & ESA Comet Interceptor Team (2019) Jones, G., & ESA Comet Interceptor Team. 2019, Comet Interceptor A Mission to a Dynamically New Solar System Object.  http://www.cometinterceptor.space/uploads/1/2/3/7/123778284/comet_interceptor_executive_summary.pdf
  * Jones et al. (2009) Jones, R. L., Chesley, S. R., Connolly, A. J., et al. 2009, Earth Moon and Planets, 105, 101, doi: 10.1007/s11038-009-9305-z
  * Jones et al. (2018) Jones, R. L., Slater, C. T., Moeyens, J., et al. 2018, Icarus, 303, 181, doi: 10.1016/j.icarus.2017.11.033
  * Kaula (1964) Kaula, W. M. 1964, Reviews of Geophysics and Space Physics, 2, 661, doi: 10.1029/RG002i004p00661
  * Kim et al. (2023) Kim, Y., DeMartini, J. V., Richardson, D. C., & Hirabayashi, M. 2023, MNRAS, 520, 3405, doi: 10.1093/mnras/stad351
  * Kim et al. (2021) Kim, Y., Hirabayashi, M., Binzel, R. P., et al. 2021, Icarus, 358, 114205, doi: https://doi.org/10.1016/j.icarus.2020.114205
  * Knight et al. (2017) Knight, M. M., Protopapa, S., Kelley, M. S., et al. 2017, The Astrophysical Journal, 851, doi: 10.3847/2041-8213/aa9d81
  * Kuipers (2007) Kuipers, J. B. 2007, Quaternions and rotation sequences: A Primer with applications to orbits, aerospace, and virtual reality (Princeton University Press)
  * Lainey et al. (2020) Lainey, V., Casajus, L. G., Fuller, J., et al. 2020, Nature Astronomy, 4, 1053, doi: 10.1038/s41550-020-1120-5
  * Laughlin & Batygin (2017) Laughlin, G., & Batygin, K. 2017, Research Notes of the American Astronomical Society, 1, 43, doi: 10.3847/2515-5172/aaa02b
  * Lee et al. (2022) Lee, H. J., Kim, M. J., Marciniak, A., et al. 2022, A&A, 661, L3, doi: 10.1051/0004-6361/202243442
  * Lellouch et al. (1995) Lellouch, E., Paubert, G., Moreno, R., et al. 1995, Nature, 373, 592, doi: 10.1038/373592a0
  * Levine et al. (2021) Levine, W. G., Cabot, S. H., Seligman, D., & Laughlin, G. 2021, The Astrophysical Journal, 922, 39, doi: 10.3847/1538-4357/ac1fe6
  * Levine & Laughlin (2021) Levine, W. G., & Laughlin, G. 2021, The Astrophysical Journal, 912, 3, doi: 10.3847/1538-4357/abec85
  * Li & Jewitt (2013) Li, J., & Jewitt, D. 2013, AJ, 145, 154, doi: 10.1088/0004-6256/145/6/154
  * Liu & O’Keefe (1965) Liu, H.-S., & O’Keefe, J. A. 1965, Science, 150, 1717, doi: 10.1126/science.150.3704.1717
  * Logg et al. (2012a) Logg, A., Mardal, K., & Wells, G. N. 2012a, Automated Solution of Differential Equations by the Finite Element Method (Springer), doi: 10.1007/978-3-642-23099-8
  * Logg & Wells (2010) Logg, A., & Wells, G. N. 2010, ACM Transactions on Mathematical Software, 37, 1–28, doi: 10.1145/1731022.1731030
  * Logg et al. (2012b) Logg, A., Wells, G. N., & Hake, J. 2012b, in Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, ed. A. Logg, K.-A. Mardal, & G. Wells (Berlin, Heidelberg: Springer Berlin Heidelberg), 173–225, doi: 10.1007/978-3-642-23099-8_10
  * MacDonald (1964) MacDonald, G. J. F. 1964, Science, 145, 881, doi: 10.1126/science.145.3635.881
  * Mainzer et al. (2015) Mainzer, A., Grav, T., Bauer, J., et al. 2015, AJ, 149, 172, doi: 10.1088/0004-6256/149/5/172
  * Mamajek (2017) Mamajek, E. 2017, Research Notes of the AAS, 1, 21, doi: 10.3847/2515-5172/aa9bdc
  * Marčeta (2023) Marčeta, D. 2023, Astronomy and Computing, 42, 100690, doi: 10.1016/j.ascom.2023.100690
  * Mashchenko (2019) Mashchenko, S. 2019, Monthly Notices of the Royal Astronomical Society, 489, 3003, doi: 10.1093/mnras/stz2380
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Meech et al. (2019) Meech, K., Castillo-Rogez, J., Hainaut, O., Lazio, J., & Raymond, S. 2019, BAAS, 51, 552
  * Meech et al. (2021) Meech, K., Castillo-Rogez, J., Bufanda, E., et al. 2021, in Bulletin of the American Astronomical Society, Vol. 53, 282, doi: 10.3847/25c2cfeb.ea404475
  * Meech et al. (2017) Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378–381, doi: 10.1038/nature25020
  * Melosh et al. (1996) Melosh, H. J., Ingram, J., & Bottke, W. F. 1996, in Lunar and Planetary Science Conference, Vol. 27, Lunar and Planetary Science Conference, 863
  * Melosh & Stansberry (1991) Melosh, H. J., & Stansberry, J. A. 1991, Icarus, 94, 171, doi: 10.1016/0019-1035(91)90148-M
  * Micheli et al. (2018) Micheli, M., Farnocchia, D., Meech, K. J., et al. 2018, Nature, 559, 223, doi: 10.1038/s41586-018-0254-4
  * Moore et al. (2021a) Moore, K., Castillo-Rogez, J., Meech, K. J., et al. 2021a, in Bulletin of the American Astronomical Society, Vol. 53, 481, doi: 10.3847/25c2cfeb.1d58e5af
  * Moore et al. (2021b) Moore, K., Courville, S., Ferguson, S., et al. 2021b, Planet. Space Sci., 197, 105137, doi: 10.1016/j.pss.2020.105137
  * Moro-Martín (2018) Moro-Martín, A. 2018, ApJ, 866, 131, doi: 10.3847/1538-4357/aadf34
  * Moro-Martín (2019) —. 2019, AJ, 157, 86, doi: 10.3847/1538-3881/aafda6
  * Moro-Martín (2022) Moro-Martín, A. 2022, arXiv e-prints, arXiv:2205.04277.  https://arxiv.org/abs/2205.04277
  * Moro-Martín et al. (2009) Moro-Martín, A., Turner, E. L., & Loeb, A. 2009, ApJ, 704, 733, doi: 10.1088/0004-637X/704/1/733
  * MPEC (2017) MPEC. 2017, The International Astronomical Union Minor Planet Center, IAU Minor Planet Center.  https://minorplanetcenter.net//mpec/K17/K17UI1.html
  * Muinonen & Lumme (2015) Muinonen, K., & Lumme, K. 2015, Astronomy & Astrophysics, 584, doi: 10.1051/0004-6361/201526456
  * Nolan et al. (2021) Nolan, M., DellaGiustina, D., Polit, A., et al. 2021, in AGU Fall Meeting Abstracts, Vol. 2021, NH13A–01
  * Noll et al. (1995) Noll, K. S., McGrath, M. A., Trafton, L. M., et al. 1995, Science, 267, 1307, doi: 10.1126/science.7871428
  * Ootsubo et al. (2012) Ootsubo, T., Kawakita, H., Hamada, S., et al. 2012, ApJ, 752, 15, doi: 10.1088/0004-637X/752/1/15
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Pau Sánchez et al. (2021) Pau Sánchez, J., Morante, D., Hermosin, P., et al. 2021, arXiv e-prints, arXiv:2107.12999.  https://arxiv.org/abs/2107.12999
  * Peale et al. (1979) Peale, S. J., Cassen, P., & Reynolds, R. T. 1979, Science, 203, 892, doi: 10.1126/science.203.4383.892
  * Peale & Gold (1965) Peale, S. J., & Gold, T. 1965, Nature, 206, 1240, doi: 10.1038/2061240b0
  * Pérez-Hernández & Benet (2022) Pérez-Hernández, J. A., & Benet, L. 2022, Communications Earth and Environment, 3, 10, doi: 10.1038/s43247-021-00337-x
  * Pettengill & Dyce (1965) Pettengill, G. H., & Dyce, R. B. 1965, Nature, 206, 1240, doi: 10.1038/2061240a0
  * Phan et al. (2021) Phan, V. H. M., Hoang, T., & Loeb, A. 2021, Erosion of Icy Interstellar Objects by Cosmic Rays and Implications for ‘Oumuamua, arXiv, doi: 10.48550/ARXIV.2109.04494
  * Pravec et al. (2014) Pravec, P., Scheirich, P., Ďurech, J., et al. 2014, Icarus, 233, 48, doi: 10.1016/j.icarus.2014.01.026
  * Quarteroni (2014) Quarteroni, A. 2014, Navier-Stokes equations (Milano: Springer Milan), 429–482, doi: 10.1007/978-88-470-5522-3_16
  * Richardson et al. (1998) Richardson, D. C., Bottke, W. F., & Love, S. G. 1998, Icarus, 134, 47, doi: 10.1006/icar.1998.5954
  * Rickman (2010) Rickman, H. 2010, in Lecture Notes in Physics, Berlin Springer Verlag, ed. J. Souchay & R. Dvorak, Vol. 790 (Springer Berlin Heidelberg), 341–399, doi: 10.1007/978-3-642-04458-8_7
  * Rumble (2022) Rumble, J. R. 2022, CRC Handbook of Chemistry and Physics, 103rd edn. (CRC Press)
  * Sánchez et al. (2017) Sánchez, P., Scheeres, D., Hirabayashi, M., & Tardivel, S. 2017, EPJ Web Conf., 140, 14004, doi: 10.1051/epjconf/201714014004
  * Scheeres et al. (2005) Scheeres, D., Benner, L., Ostro, S., et al. 2005, Icarus, 178, 281, doi: https://doi.org/10.1016/j.icarus.2005.06.002
  * Schenk et al. (1996) Schenk, P. M., Asphaug, E., McKinnon, W. B., Melosh, H. J., & Weissman, P. R. 1996, Icarus, 121, 249, doi: 10.1006/icar.1996.0084
  * Sekanina (2019) Sekanina, Z. 2019, Outgassing As Trigger of 1I/‘Oumuamua’s Nongravitational Acceleration: Could This Hypothesis Work at All?, arXiv, doi: 10.48550/ARXIV.1905.00935
  * Seligman & Laughlin (2018) Seligman, D., & Laughlin, G. 2018, AJ, 155, 217, doi: 10.3847/1538-3881/aabd37
  * Seligman & Laughlin (2020) Seligman, D., & Laughlin, G. 2020, The Astrophysical Journal, 896, doi: 10.3847/2041-8213/ab963f
  * Seligman et al. (2021) Seligman, D. Z., Levine, W. G., Cabot, S. H., Laughlin, G., & Meech, K. 2021, The Astrophysical Journal, 920, 28, doi: 10.3847/1538-4357/ac1594
  * Seligman et al. (2023) Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, Accepted for publication in Planet. Sci. J., arXiv:2212.08115.  https://arxiv.org/abs/2212.08115
  * Shoemaker et al. (1993) Shoemaker, C. S., Shoemaker, E. M., Levy, D. H., et al. 1993, IAU Circ., 5725, 1
  * Smith et al. (1979) Smith, B. A., Soderblom, L. A., Johnson, T. V., et al. 1979, Science, 204, 951, doi: 10.1126/science.204.4396.951
  * Solontoi et al. (2011) Solontoi, M., Ivezić, Ž., & Jones, L. 2011, in American Astronomical Society Meeting Abstracts# 217, Vol. 217, 252–11
  * Taylor (2023) Taylor, A. 2023, Simulator of Asteroid Malformation Under Stress, v1.0.0, Zenodo, doi: 10.5281/zenodo.7537861
  * Trilling et al. (2017) Trilling, D. E., Robinson, T., Roegge, A., et al. 2017, ApJ, 850, L38, doi: 10.3847/2041-8213/aa9989
  * Trilling et al. (2018) Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, The Astronomical Journal, 156, 261, doi: 10.3847/1538-3881/aae88f
  * Van Rossum & Drake (2009) Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)
  * Vereš & Chesley (2017a) Vereš, P., & Chesley, S. R. 2017a, AJ, 154, 13, doi: 10.3847/1538-3881/aa73d0
  * Vereš & Chesley (2017b) —. 2017b, AJ, 154, 12, doi: 10.3847/1538-3881/aa73d1
  * Vilella & Deschamps (2017) Vilella, K., & Deschamps, F. 2017, Journal of Geophysical Research: Planets, 122, 1056–1076, doi: 10.1002/2016je005215
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Watanabe (1992) Watanabe, J.-I. 1992, PASJ, 44, 163
  * Weaver et al. (1995) Weaver, H. A., A’Hearn, M. F., Arpigny, C., et al. 1995, Science, 267, 1282, doi: 10.1126/science.7871424
  * Williams et al. (2017) Williams, G. V., Sato, H., Sarneczky, K., et al. 2017, Central Bureau Electronic Telegrams, 4450, 1
  * Wisdom et al. (1984) Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137, doi: 10.1016/0019-1035(84)90032-0
  * Yu et al. (2014) Yu, Y., Richardson, D. C., Michel, P., Schwartz, S. R., & Ballouz, R.-L. 2014, Icarus, 242, 82, doi: 10.1016/j.icarus.2014.07.027
  * Zwart et al. (2018) Zwart, P., S., Torres, S., Pelupessy, I., Bédorf, J., & Cai, M. X. 2018, MNRAS, 479, L17, doi: 10.1093/mnrasl/sly088

##  Appendix A Light Curve Model Validation

Figure 15: Diagram of the cross-sectional ellipse. Unprimed coordinates are
fixed in the reference frame of the Solar System, while primed coordinates are
fixed in the reference frame of the object. The angle $\beta$ defines the
rotation state of the body with respect to the Sun.

In this subsection, we validate the use of our ML15 model (described in
Section 4) via comparison to a simple principal-axis rotation that the complex
model can be easily reduced to. To accomplish this, we produce a synthetic
photometric time series of ‘Oumuamua assuming principal-axis rotation, a
simplified viewing geometry, and no modulation from phase angle. We position
the Sun, the Earth, and ‘Oumuamua along the x-axis in conjunction, roughly
representative of the astrometric orbital arrangement when ‘Oumuamua was
discovered. The non-rotating unprimed (x, y, z) axes move with ‘Oumuamua along
its trajectory, while the primed coordinate axes (x’, y’, z’) align with the
principal axes and co-rotate with the body. We orient the body such that its
semi-major principal axes are a:b:c$\sim$6:6:1 (Mashchenko, 2019), and we
assume that ‘Oumuamua rotates solely about the y-axis. Rotation about the
x-axis produces a flat light curve, since the Sun-pointing projected area is
unchanged. Due to ‘Oumuamua’s (assumed) x-y symmetry, rotation about the
z-axis similarly produces no change in the projected area and a flat light
curve, leaving the y-axis as the only non-trivial principal axis rotation.

We derive the projected area of an ellipsoid with an arbitrary rotation angle
(with respect to the x-axis) $\beta$. It is known that any cross-section of an
ellipsoid is an ellipse, although possibly with different aspect ratio and
orientation. With rotation about the y-axis, the principal axes of the ellipse
are along the Cartesian axes, and we solve for the lengths of these
2-dimensional principal axes. The semi-major axis in the y-direction is $b$,
and so we now solve for the semi-major axis in the x-z plane which is
orthogonal to the rotation axis. In this plane, the cross-section has an
aspect ratio of 6:1, with $a=6$ in the x’-direction and $c=1$ in the
z’-direction (Figure 15). The observation vector, defined to be parallel to
$\boldsymbol{\hat{x}}$, is rotated about the y-axis such that the rotated axis
is
$\cos\beta\,\boldsymbol{\hat{x}^{\prime}}+\sin\beta\,\boldsymbol{\hat{z}^{\prime}}$,
in the primed frame. We then find the location of the points
$(x^{\prime}_{0},z^{\prime}_{0})$ where the tangent to the ellipse is parallel
to the observation vector, whose projection along the z-axis is this semi-
major axis.

The cross-sectional ellipse is defined by the function
$f_{\text{2d}}(x^{\prime},z^{\prime})$, where

$f_{\text{2d}}(x^{\prime},z^{\prime})=\frac{x^{\prime
2}}{a^{2}}+\frac{z^{\prime 2}}{c^{2}}-1.$ (A1)

The tangent $\partial z^{\prime}/\partial x^{\prime}$ is the derivative with
respect to x of Equation A1, given by

$\frac{\partial z^{\prime}}{\partial
x^{\prime}}(x)=\mp\frac{cx^{\prime}}{a^{2}\sqrt{1-x^{\prime 2}/a^{2}}}\,.$
(A2)

The $\mp$ reflects the degeneracy of the ellipse along slices in the x’ axis.
The slope of the observation vector in the x’-z’ plane is simply
$\tan(\beta)$, so by setting $\partial z^{\prime}/\partial
x^{\prime}=\tan(\beta)$, and simplifying, we obtain

$x^{\prime}_{0}=\pm\frac{a^{2}\tan\beta}{\sqrt{c^{2}+a^{2}\tan^{2}(\beta)}}\,.$
(A3)

Now by substituting Equation A3 into $z^{\prime}=\pm c\sqrt{1-x^{\prime
2}/a^{2}}$ (derived from Equation A1), the points
$(x^{\prime}_{0},z^{\prime}_{0})$ are given by

$(x^{\prime}_{0},z^{\prime}_{0})=(\pm\frac{a^{2}\tan\beta}{\sqrt{c^{2}+a^{2}\tan^{2}(\beta)}},\mp\frac{c^{2}}{\sqrt{c^{2}+a^{2}\tan^{2}(\beta)}}).$
(A4)

The distance between these points and the observation line which passes
through the origin is the second semi-major axis of the projected ellipse. For
a point $(x^{\prime}_{0},\ z^{\prime}_{0})$, the distance is the projection
onto the z-axis, so $d=|\sin\beta\,x^{\prime}_{0}+\cos\beta\,z^{\prime}_{0}|$.
Substituting Equation A4, the distance is

$\displaystyle
d=\frac{|a^{2}\sin\beta\tan\beta+c^{2}\cos\beta|}{\sqrt{c^{2}+a^{2}\tan^{2}(\beta)}}\,.$
(A5)

Since the area of the ellipse is $\pi bd$, the brightness $L_{\text{fix}}$
(with simplification) is given by

$L_{\text{PA}}=\pi b\sqrt{(a\sin\beta)^{2}+(c\cos\beta)^{2}}.$ (A6)

We verified (not shown) that this simple analytic model roughly reproduces the
photometric data given in Belton et al. (2018).

To compare and validate the ML15 model versus this one, we computed the error
between the ML15 model and the principal-axis model described in this section.
The error between these light curves is at machine precision. Similarly, we
constructed an additional, numeric model allowing for “fixed-axis” NPA
rotation, while maintaining the fixed and simplified astrometric arrangement.
The error between this model and ML15 is of $\mathcal{O}(10^{-5})$. These
comparisons validate the use of the ML15 model in this paper.

##  Appendix B Pseudocode

In this subsection we present a pseudocode which demonstrates the most basic
functionality of SAMUS. As a class method, SAMUS is capable of running in a
modular format with complex function call patterns. In steps 1-8, we
initialize the simulation. In steps 9.a-9.d, we compute the forcing functions
described in Section 2.3, Table 1, and Equation 5. In steps 9.e-9.i, we update
the model for a single time step. In steps 10 and 11, we implement the
trajectory jump method described in Section 2.2, and in step 12, we save and
output simulation products. A script to replicate this example is available in
the SAMUS package.

  1. SAMUS Pseudocode

  2. 1.

Set parameters.

  3. 2.

Read in 3-ball mesh and reshape to ellipsoid.

  4. 3.

Read in trajectory data.

  5. 4.

Create functions and function spaces for all variables.

  6. 5.

Create a UFL form for the Navier-Stokes equations.

  7. 6.

Create a FEniCS non-linear solver for the Navier-Stokes equations.

  8. 7.

Create a FEniCS non-linear solver for the Gaussian gravity formulation.

  9. 8.

Set t=0.

  10. 9.

WHILE

`Number of cycles < number in loop:`

    1. (a)

Update tidal force, computed using Eqn 7.

    2. (b)

Update gravitational force, computed using Eqn 5.

    3. (c)

Update Coriolis force, computed using the relevant term from Table 1.

    4. (d)

Update centrifugal force, computed using the relevant term from Table 1.

    5. (e)

Solve Navier-Stokes equations. If this diverges, STOP.

    6. (f)

Move mesh with velocity.

    7. (g)

Check if deformation crosses threshold. If it does, STOP.

    8. (h)

Compute moment of inertia using Equation 6.

    9. (i)

Update time step.

  11. 10.

Average velocities over the rotation cycles.

  12. 11.

WHILE

Change in the heliocentric distance is less than tolerance:

    1. (a)

Check to ensure that CFL<C${}_{\texttt{max}}$.

    2. (b)

Move mesh using average velocity.

    3. (c)

Update time step.

  13. IF:

The number of steps is less than the upper limit, then repeat Steps 8-11.

  14. 12.

Save the mesh, the functions, and the moments of inertia over the path.

This pseudocode structure is also shown in a flowchart available in the SAMUS
documentation (Taylor, 2023).