21-785

# Gravity Estimation at Small Bodies via Optical Tracking of Hopping
Artificial Probes

Jacopo Villa Graduate Student, Aerospace Engineering Sciences, University of
Colorado Boulder, 3775 Discovery Drive, Boulder, CO 80303 Andrew French
Navigation Engineer, Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA Jay McMahon
Assistant Professor, Aerospace Engineering Sciences, University of Colorado
Boulder, 3775 Discovery Drive, Boulder, CO 80303 Daniel Scheeres Distinguished
Professor, Aerospace Engineering Sciences, University of Colorado Boulder,
3775 Discovery Drive, Boulder, CO 80303 and Benjamin Hockman Robotics
Technologist, Jet Propulsion Laboratory, California Institute of Technology,
4800 Oak Grove Drive, Pasadena, CA 91109, USA

###### Abstract

Despite numerous successful missions to small celestial bodies, the gravity
field of such targets has been poorly characterized so far. Gravity estimates
can be used to infer the internal structure and composition of small bodies
and, as such, have strong implications in the fields of planetary science,
planetary defense, and in-situ resource utilization. Current gravimetry
techniques at small bodies mostly rely on tracking the spacecraft orbital
motion, where the gravity observability is low. To date, only lower-degree and
order spherical harmonics of small-body gravity fields could be resolved. In
this paper, we evaluate gravimetry performance for a novel mission
architecture where artificial probes repeatedly hop across the surface of the
small body and perform low-altitude, suborbital arcs. Such probes are tracked
using optical measurements from the mothership’s onboard camera and orbit
determination is performed to estimate the probe trajectories, the small
body’s rotational kinematics, and the gravity field. The suborbital motion of
the probes provides dense observations at low altitude, where the gravity
signal is stronger. We assess the impact of observation parameters and mission
duration on gravity observability. Results suggest that the gravitational
spherical harmonics of a small body with the same mass as the asteroid Bennu,
can be observed at least up to degree 40 within months of observations.
Measurement precision and frequency are key to achieve high-performance
gravimetry.

##  1 Introduction

Our understanding of primitive bodies in the Solar System has grown
significantly over the last two decades, especially thanks to an ever
increasing number of exploration missions. Just recently, NASA’s OSIRIS-REx
and JAXA’s Hayabusa2 collected a wealth of scientific data around the Near-
Earth Asteroids Bennu and Ryugu, respectively, and successfully sampled
surface material to be studied on Earth. Several missions toward small bodies
are planned for the near future: NASA’s Janus, Psyche, and Lucy will explore
the exotic surroundings of binary asteroids, a nickel-iron small body, and
Main Belt’s as well as Trojan asteroids, respectively; AIDA, a NASA/ESA joint
mission, will test and assess asteroid deflection strategies in orbit for the
first time. Whilst past missions provided invaluable insights on the
geophysics of these worlds, many questions about small bodies remain
unanswered. One key feature requiring further investigation is the interior
structure and composition. Evidence shows that small body interiors can vary
dramatically, from mostly homogeneous to largely heterogeneous composition,
with under- or over-dense regions, including voids typical of rubble piles.
However, there is a lack of direct observations to confirm current geophysical
models[1]. Enhancing the knowledge of small body interiors could help to (1)
validate and refine current models of the solar system formation and
evolution, (2) assess the effectiveness of asteroid deflection techniques, and
(3) augment subsurface situational awareness for asteroid mining and resource
exploitation purposes[2].

Three main techniques have been used and proposed to map the interiors of
small bodies: radar tomography, seismic imaging, and gravimetry. Radar
tomography can characterize the 3D interior structure of a small body by
measuring the reflection and transmission of ground-penetrating radio signals.
This technique has been successfully employed to measure the dielectric
constant of comet 67P during the Rosetta mission and several other mission
concepts based on radar tomography have been proposed[3]. However, it requires
more complex mission architectures and instrumentation, and presents multiple
challenges for the global mapping of a small body. First, 3D reconstruction of
the deeper areas can be difficult, due to signal attenuation and unfavorable
observation geometry. Second, data inversion can be challenging and
unreliable, and it is recommended that this technique is informed by other
interior-mapping strategies, such as gravimetry. Third, certain materials
composing small bodies can be opaque to electromagnetic energy, undermining
the effectiveness of this method. Finally, overall performance strongly
depends on the target body’s dielectric properties as well as the selected
mission architecture and instrumentation, requiring mission-specific design
and pre-mission calibration[4, 5].

Seismic imaging measures the elastic waves traversing the body to infer its
interior properties, and has been widely used to map the internal structure of
the Earth. Several mission concepts have been proposed to apply the same
principle at small bodies, deploying seismometers and detonators on the
surface[6, 7]. However, a number of practical challenges arise, such as the
precise deployment and localization of surface probes and the unknown
dissipation properties of the target body. A simpler technique to infer crude
estimates of a small body’s interior properties consists in artificially
inducing seismic shaking with a kinetic impactor. Observing the displacement
of surface features caused by the impact event can reveal some of the body’s
internal properties, such as the bulk density. This technique was recently
adopted by the Hayabusa2 spacecraft at the asteroid Ryugu; however, no
relevant regolith displacements were detected, possibly due to the porosity
and plastic behavior of the asteroid’s surface, which dissipated seismic waves
caused by the impact event[8].

Finally, gravity measurements (or gravimetry) can be used to interpret the
interior properties of a small body, in particular, its mass distribution.
Gravimetry has a long history of successes in small body missions, and
consists in tracking the spacecraft motion (usually with ground-based
radiometric measurements and onboard optical measurements) in the vicinity of
the small body, co-estimating the spacecraft trajectory and the central body’s
gravity field. Since this technique is inherently part of the navigation
process, it does not necessarily require additional instruments and
operations. Ideally, a high-precision estimate of the gravity field would
require low-altitude observations (where the gravity signal is stronger),
large spatial distribution of measurements around the body and precise
measurements. This is especially important for low-mass target bodies, as the
gravity observability is lower. However, mission requirements usually
constrain spacecraft to fly far from the surface, where orbits are more
stable, but the gravity signal is weaker. As such, the low-signal gravity
anomalies remain undetected and, in turn, the heterogeneous mass distribution
largely undiscovered.

To date, gravity estimation based on spacecraft tracking only provided
estimates for the lower degrees and orders of the gravity field spherical
harmonics, which in turn only permits low-fidelity characterization of the
interior mass distribution. On the other hand, the OSIRIS-REx mission recently
pioneered a novel gravimetry approach: the mission team coincidentally
discovered natural ejecta flying in the vicinity of the asteroid Bennu and, by
tracking the particles motion, could estimate the asteroid’s gravity field to
about degree and order 8 [9]. For comparison, only the $J_{2}$ coefficient of
Bennu was recovered by tracking the spacecraft, at the end of the Orbital A
phase of that same mission[10]. This unprecedented success proved the
potential of measuring low-altitude, spatially-distributed objects to
drastically improve gravimetry performance at small bodies. However, tracking
natural ejecta comes with limitations: these natural events are not
controllable nor predictable, and they may not occur at all when visiting a
given small body. Furthermore, the ejecta are fragments of a low-albedo
celestial body and their size in the camera plane is sub-pixel. Hence, they
are associated with large astrometric uncertainty which, in the OSIRIS-REx
case, with a median of 1.05 pixels[9]. Additionally, such particles can hardly
be detected in front of the small body, since their photometric properties
match the one of the body surface. As such, a large fraction of low-altitude
measurements (those carrying the stronger gravity signal) cannot be detected
and used for gravimetry.

In this paper, we evaluate gravimetry performance for a novel mission concept
that leverages measurements of low-altitude, spatially-distributed objects for
gravimetry at small bodies. Instead of relying on serendipitous and
unpredictable observations of natural particles, we propose to track a swarm
of artificial probes designed to maximize gravimetry performance. We assume
that such probes are characterized by two key features: (i) they can
repeatedly hop from the surface of the small body, performing suborbital
motion, and (ii) they are visible both on the dark background and in front of
the small body surface, from the observer’s viewpoint. The probes motion is
observed by a mother spacecraft hovering in the vicinity of the small body.
This mission architecture was firstly presented as a NASA Innovative Advance
Concept (NIAC) study, named Gravity Poppers, by Hockman et al. and is proposed
as a low-budget approach for high-precision gravimetry at small
bodies111https://www.nasa.gov/directorates/spacetech/niac/2020_Phase_I_Phase_II/Gravity_Poppers/.
An infographic of this concept is shown in Figure 1. In its simplest design,
this concept only requires a camera aboard the mothership to observe the
probes motion. The probes could be actively emitting light with LEDs or
passively reflecting the sunlight. This paper is structured as follows: we
firstly define the problem statement and the simulation setup. Second, we
present the assumptions on the dynamical environment and measurements. Third,
we present orbit determination (OD) results and how they can inform the
mission design.

Figure 1: Infographic of the Gravity Poppers mission concept.

##  2 Problem Statement

The following mission architecture is considered: a set of artificial probes
performs suborbital motion around a small body, repeatedly hopping from its
surface once they come at rest. The probes motion is observed by a mother
spacecraft, hovering in the vicinity of the small body, which detects the
probes using the onboard camera. Using optical measurements, the gravity field
of the small body is estimated using an OD process. The objective of this work
is threefold: (1) quantify gravity estimation performance for a notional and
representative mission design, (2) analyze the sensitivity of gravimetry
performance to the most relevant design parameters, and (3) present order-of-
magnitude requirements for mission duration, measurement performance, and
mission architecture.

##  3 Simulation Setup

This study is carried out using numerical simulations. The simulated scenarios
are implemented using the Mission Analysis, Operations, and Navigation Toolkit
Environment (MONTE) developed by NASA Jet Propulsion Laboratory (JPL)[11].
MONTE provides a state-of-the-art toolkit for mission analysis, trajectory
design, flight path control, and navigation and is used as the prime orbit
determination software for missions flown by JPL. MONTE can be used as an
importable Python module and we integrate it in a broader Python code for
post-processing and data analysis.

Simulating the hopping motion of the probes around the small body calls for
some considerations on the probe dynamics. A probe can leave the surface and
enter a suborbital arc mainly because of two different events: either an
active hop initiated by the probe itself, or a natural bouncing event caused
by the probe’s residual kinetic energy. In the latter case, the probe’s
bouncing velocity is largely unpredictable, as it is affected by local
topography and surface absorption properties, which are usually heterogeneous
throughout the small body. Figure 2 shows the random nature of the hopping
motion around a small body. As such, modeling the deterministic links between
consecutive arcs proves challenging and would likely provide marginal gains
for gravity estimation. In addition, the trajectory arcs produced by natural
bounces may not be long enough to substantially enhance gravity knowledge.
Hence, for gravity estimation, we only consider the arcs caused by the probe’s
active hopping events, i.e., we neglect all bounces of the probes with the
surface. We assume that the arcs are uncorrelated with each other, and we
process them separately for OD.

Figure 2: Simulated trajectories of hopping probes around comet 67P Churyumov-
Gerasimenko, in the small body-fixed frame. The dynamics includes both
artificial hops (longer arcs) and natural bounces with the surface
(progressively shorter arcs). Velocity magnitude and direction are randomized
to simulate the effect of local surface properties.

###  3.1 Target Body

The population of small bodies in the solar system is highly diverse in many
regards, such as composition, orbit, size, bulk density, and mass. Since the
gravity signal is proportional to the small body’s mass, this is a key
parameter for gravity performance. Past missions have proven that higher-order
gravity terms of small bodies are unobservable when gravimetry is performed by
tracking an orbiting spacecraft. This is especially true for low-mass bodies.
In order to consider the worst-case scenario for gravimetry, we assume that
the simulated target body is a low-mass asteroid with the same mass, rotation
rate, and pole as the asteroid Bennu visited by the OSIRIS-REx
mission[12][13]. Table 2 captures the key parameters for the simulated small
body, whereas its non-uniform gravity model is described in the dynamics
section. We assume that the body orbits at 1 Astronomical Unit (AU) from the
Sun.

Table 1: Parameters of the simulated small body. Parameter |  Value |  Units
---|---|---
$\mu$ |  4.892 |  $\mathrm{m}^{3}/\mathrm{s}^{2}$
$\mu_{\mathrm{mascon}}$ |  0.0087 |  $\mathrm{m}^{3}/\mathrm{s}^{2}$
Mascon [$\phi$,$\lambda$] |  [45,45] |  deg
Radius |  246 |  m
Pole [RA,Dec,PM] |  [85.65,-60.17,54.77] |  deg
Rotation Period |  4.2 |  hours
Albedo |  0.04 |  –
Emissivity |  0.9 |  –

Table 2: Parameters of the notional probes. Parameter |  Value |  Units
---|---|---
Shape |  Sphere |  –
Radius |  3 |  cm
Mass |  0.2 |  kg
Albedo |  0 |  –
Emissivity |  0 |  –
Hop. $\Delta v$ Range |  [3,14] |  cm/s

##  4 Dynamical Environment

The ground-truth probe trajectories are simulated accounting for the main
forces acting in the small-body environment. We consider the central body’s
non-uniform gravity, N-body gravitational pulls, the solar radiation pressure
(SRP), as well as the thermal radiation pressure (TRP) and the albedo pressure
(AP) from the small body. Figure 3 shows the evolution of some representative
accelerations as a function of the probe’s altitude. In addition, stochastic
accelerations are included to simulate the presence of unmodeled forces. This
dynamics model is only used to simulate and estimate the probe’s motion, as we
assume that the observer, i.e. the mother spacecraft, is constantly hovering
at a predefined altitude with no forces acting on it.

###  4.1 Probe Design

We consider a notional probe whose parameters are reported in Table 2. We
consider the two probe designs proposed by Hockman et al. in the NIAC study:
(1) an active, 5-cm, cubic-shaped probe with strobing LEDs to enable detection
and (2) a passive, white-coated, spherical-shaped probe detected by passively
reflecting the sunlight. The cubic design presents multiple advantages, such
as reduced size (and hence reduced non-gravitational perturbations) and better
photometric properties, thanks to the use of low-power LEDs which make the
probe visible even when in the body’s umbra. Hence, this second design is
considered here for dynamics modeling. Since the probe is likely subject to
fast rotation during its motion, we consider its interaction with radiation
pressures as isotropic, and hence we approximate the cube with its
circumscribing sphere. Furthermore, we do not account for reflected and
infrared radiations from the probe, hence we set its albedo and emissivity to
zero. These assumptions are further discussed and justified in the radiation
pressure section.

###  4.2 Gravity

The gravity of major solar system bodies is considered. These are modeled as
point-mass gravity sources using the de430 ephemeris data set[14].

The non-uniform gravity of the target small body is simulated using an
artificial gravity field, since no estimates of higher degree and order for
small bodies are available to date. We use an exterior spherical harmonic
model, where the gravity potential can be written as the sum of the point-mass
potential $\mu/r$ and a spherical harmonic expansion. In spherical coordinates
$(r,\phi,\lambda)$ and with respect to a body-fixed reference frame, this can
be written as shown in Equation 1[15].

$\begin{split}U(r,\phi,\lambda)=&\frac{\mu}{r}+U^{\prime}\\\
U^{\prime}=&-\frac{\mu^{\prime}}{r}\sum_{l=1}^{\infty}\left(\frac{a_{e}}{r}\right)^{l}P_{l}(\mathrm{sin}(\phi))J_{l}+\\\
&+\frac{\mu^{*}}{r}\sum_{l=1}^{\infty}\sum_{m=1}^{l}\left(\frac{a_{e}}{r}\right)^{l}P_{lm}(\mathrm{sin}(\phi))[C_{lm}\mathrm{cos}(m\lambda)+S_{lm}\mathrm{sin}(m\lambda)]\end{split}$
(1)

where $\mu$ is the standard gravitational parameter of the central body, $r$
is the range from the central body’s center of mass to the spacecraft of
interest, $a_{e}$ and $\mu^{*}$ are the reference radius and the reference
gravitational parameter used to normalize the harmonic coefficients $C_{lm}$
and $S_{lm}$. $P_{lm}$ are the Legendre Associated Functions of degree $l$ and
order $m$. $J_{l}$, $C_{lm}$ and $S_{lm}$ are the normalized spherical
harmonic coefficients, where $J_{l}=-C_{l0}$ represents the so-called zonal
harmonics; $C_{lm}$ and $S_{lm}$ are referred to as sectoral harmonics when
$l\neq m$ and as tesseral harmonics when $l=m$. The spherical harmonic
expansion can be truncated to model the number of degree and orders of
interest.

We artificially generate the above harmonic coefficients. It has been shown
that the amplitude of such parameters decays with degrees by following a power
law, i.e., the so-called Kaula rule[16]. McMahon et al. forward-modeled a
modified Kaula rule for the asteroid Bennu, which is defined by Equation 2 and
whose parameters are reported in Table 3[17]. We generate the ground-truth
gravity field by randomly sampling the harmonic coefficients from the Gaussian
distribution associated with the Kaula rule for Bennu, up to degree and order
20. Figure 4 shows the resulting ground-truth gravity coefficients. As shown
in Equation 1, the number of coefficients for each degree increases with the
degree itself, hence simulating a large number of trajectories using a gravity
field larger than 20$\times$20 proved challenging and requires substantial
computational resources. This will be addressed in future work.

$\begin{split}C_{l0}=&-J_{l}\approx\frac{K_{\mathrm{zonal}}}{n^{\alpha_{\mathrm{zonal}}}},\\\
C_{lm}(m\geq 1)\;\mathrm{and}\;S_{lm}(m\geq
1)=&\sqrt{\frac{\sum_{m=1}^{n}(C_{lm}^{2}+S_{lm}^{2})}{2l}}\approx\frac{K_{\mathrm{rms}}}{n^{\alpha_{\mathrm{rms}}}}\end{split}$
(2)

Table 3: Mean and $1\sigma$ standard deviation of the Bennu’s modified Kaula rule parameters. | K | $\alpha$
---|---|---
Zonal | $0.084\pm 0.021$ | $2.08\pm 0.13$
RMS | $0.026\pm 0.005$ | $2.01\pm 0.12$

Figure 3: Sample gravitational and non-gravitational accelerations as a
function of altitude.
Figure 4: Simulated harmonic coefficients of the ground-truth gravity field.

It is known that the exterior spherical harmonic expansion only converges when
the position vector is located outside of the central body’s Brillouin
(circumscribing) sphere. Oftentimes, the shapes of small bodies present
concave regions, hence said condition for orbital motion is not always
satisfied. Several alternative gravity models have been proposed to model the
spacecraft motion within the body’s Brillouin sphere, e.g. the so-called
interior models[18]. To limit the scope of this study, we assume that the
small body considered has a perfectly spherical shape. As such, The motion of
the hopping probes around its surface can be fully represented with the
exterior spherical harmonics model. Future work will include interior gravity
models in the simulations.

In addition to the non-uniform gravity provided by the spherical harmonic
model, we simulate the presence of a surface mass concentration (or mascon) to
assess the observability of gravity anomalies introduced by local surface
feature, such as large boulders. We assume the mascon’s mass to be a rough
estimate of BenBen Saxum’s mass (Table 2), the largest boulder on the surface
of Bennu.

###  4.3 Radiation Pressure

The probe acceleration given by the solar radiation pressure is computed as

$\ddot{\mathbf{r}}_{SRP}=\frac{f\Psi}{mr_{\Sun}^{2}}\mathbf{F}_{N}$ (3)

where $f$ is the shadowing scaling factor which is 1 when the probe is outside
the shadow of the body, 0 when the probe is in the umbra of the body shadow,
and between 0 and 1 when the probe is in the penumbra of the body shadow;
$\Psi$ is the solar flux at 1 AU, in Newton; $m$ is the probe mass;
$\mathbf{r}_{\Sun}$ is the position vector from the probe to the Sun [];
$\mathbf{F}_{N}$ is the SRP geometry vector for a spherical body, computed as

$\mathbf{F}_{N}=-\pi
R^{2}\left(1+\frac{4}{3}\kappa_{d}\nu_{d}\right)\hat{r}_{\Sun}$ (4)

where $R$ is the sphere radius, $\kappa_{d}$ is the diffuse reflectivity
degradation factor and $\nu_{d}$ is the diffuse reflectivity factor. $\nu_{d}$
is defined as

$\nu_{d}=\frac{1}{3}(\gamma(1-\beta)+k(1-\gamma))$ (5)

where $\gamma$ is the fraction of incoming radiation that is reflected and
$\beta$ is the fraction of reflected radiation that is specular. Since we
assume $\gamma=0$ (i.e., the probe albedo is zero), we only consider the
effect of direct radiation pressure. Hence, $\ddot{\mathbf{r}}_{SRP}$ can be
rewritten as

$\ddot{\mathbf{r}}_{SRP}=f\Psi\frac{A}{m}\frac{\hat{r}}{r_{\Sun}^{2}}$ (6)

which is often called the cannonball model.

We use a standard cannonball model to account for SRP accelerations, which
assumes the body is a sphere[19]. This is usually a reasonable assumption even
for non-spherical bodies subject to fast tumbling, such as the notional cubic
probes, since the anisotropic SRP contributions even out over time. We assume
that the body has zero albedo and emissivity, i.e. the effect of the SRP only
accounts for the direct photon flux, and not for reflections nor re-emitted
infrared radiation. Note that this assumption mostly affects the ground-truth
trajectory and not gravity estimation results.

The small body is also a source of radiation pressure, with a component given
by the reflected sunlight (the albedo pressure, or AP), and one from the
thermal re-emission due to the body’s heating from the Sun (the thermal
radiation pressure, or TRP). In this case, the AP is given by:

$\displaystyle\ddot{\mathbf{r}}_{AP}=\iint_{S}\left(\frac{\Psi}{r_{\Sun}^{2}}\frac{A}{m}a(\mathbf{r}_{S})\,\mathrm{cos}(\Phi_{S})\frac{1}{\pi
r_{S}^{2}}\hat{r}_{\Sun}\right)dS$ (7)

where $S$ is the surface of the spherical cap visible by the spacecraft,
$a(\mathbf{r}_{S})$ is the surface albedo corresponding to the surface point
$\mathbf{r}_{S}$ and $\Phi_{S}$ is the angle that the sunlight direction makes
with the surface normal at $\mathbf{r}_{S}$. The same expression can be used
to compute the albedo pressure component given by the small body’s thermal
emission, by substituting the albedo $a(\mathbf{r}_{S})$ with the emissivity
function $e(\mathbf{r}_{S})$. Here, we assume that both the albedo and the
emissivity are constant parameters throughout the body surface (Table 2).

###  4.4 Stochastic Accelerations

In addition to the deterministic forces described above, stochastic
accelerations are included in the dynamics model to compensate for both
unmodeled forces and modeling errors in the dynamics. As discussed by Chesley
et al. for the OSIRIS-REx particle tracking case, there is a number of other
relevant forces affecting the motion in the vicinity of the small body, such
as the SRP and AP radiation reflected and re-emitted by the probes and the
Poynting-Robertson effect, due to stellar aberration affecting the incoming
SRP[9]. In this case study, mismodeling of the probe dynamics may be caused by
residual probe outgassing as well as changes in probe reflectance and
emissivity properties, due to repeated interactions with the surface regolith.

We model stochastic accelerations as an uncorrelated white random process with
a given mean and standard deviation, divided into consecutive batches. For
estimation, each batch can have a different acceleration mean, i.e.:

$\hat{\mathbf{a}}=\\{\hat{\mathbf{a}}_{1},...,\hat{\mathbf{a}}_{n}\\}\;\;\mathrm{for}\;\;\Delta
t=\\{\Delta t_{1},...,\Delta t_{n}\\}$ (8)

where $\hat{\mathrm{a}}$ are the estimated accelerations and $\Delta t$ are
the corresponding $n$ batch intervals. In the ground-truth dynamics model, we
assume that all batches are zero-mean. In the OSIRIS-REx scenario, the
unmodeled forces were always smaller than $1\times
10^{-11}\mathrm{km}/\mathrm{s}^{2}$[9]. As such, we set the stochastic
accelerations’ standard deviation as $\sigma_{\mathrm{st}}=1\times
10^{-11}\mathrm{km}/\mathrm{s}^{2}$ to make a conservative assumption on the
unmodeled forces. Such stochastic accelerations also account for errors in
modeling the SRP and AP, whose magnitude is always below said value (Figure
3). The batch interval is set to 10 minutes. Since the batch interval is a
tunable parameter for estimation, future work will focus on more rigorous
tuning strategies to optimize OD performance.

###  4.5 Hopping Events

As previously discussed, all arcs considered in this study are initiated by a
hopping event: the probe, initially steady on the terrain, activates the
hopping mechanism which manifests as a $\Delta v$ with respect to the small
body surface. Even if we assume that the hopping mechanism is calibrated to
provide a constant impulse, both the magnitude and direction of the hopping
velocity is randomized by the heterogeneous properties of the small body
surface. The hop velocity magnitude depends on the local absorption properties
of the surface regolith, which can vary from dusty to rocky regions. The
velocity direction is largely affected by the terrain topography, which can
feature slopes and obstacles and also varies across the surface.

To model these effects in our simulations, we randomly sample the hopping
velocity as follows. The velocity magnitude is sampled from a uniform
distribution:

$v_{\mathrm{hop}}=\mathcal{U}(3,14)\;\mathrm{cm/s}$ (9)

The distribution range is selected via trial and error to enable arcs with
good gravity observability (especially for the higher-degree and order terms).
The direction of the hopping velocity is defined using its right ascension
(RA) and declination (Dec) Euler angle, defined with respect to a local-
vertical local-horizontal frame. These angles are also sampled using uniform
distributions:

$\begin{cases}RA_{\mathrm{hop}}=2\pi u,\;\mathrm{where}\;u=\mathcal{U}(0,1)\\\
Dec_{\mathrm{hop}}=\mathrm{cos}^{-1}(\nu),\;\mathrm{where}\;\nu=\mathcal{U}(0,0.5)\\\
\end{cases}$ (10)

where the declination angle can span between 0 and 60 degrees. The
$\mathrm{cos}^{-1}$ formulation is used to avoid biasing the velocity. (If
spherical coordinates are directly sampled from a uniform distribution, the
resulting direction is effectively biased toward the pole(s).)

The initial position of the probes is also randomized to ensure that the
trajectory arcs spread across the entire body. The latitude $\phi_{0}$ and
longitude $\lambda_{0}$ of the probe’s initial position are randomly sampled
using the same approach previously described, whereas the radial coordinates
is set equal to the small body radius.

###  4.6 Arc Population

Using the above considerations, we simulate a baseline data set with 4000
randomized trajectory arcs around the target body. A sample of 100 arcs is
visualized in Figure 5, whereas some descriptive parameters of the resulting
data set are shown in Figure 6. The mean arc lifetime is about 1 hour and the
altitude of the apoapsis with respect to the surface is between tens and
hundreds of meters, for most cases. The idea behind such short, low-altitude
arcs is to maximize the time spent close to the surface, where the gravity
signal is stronger (Figure 3). Additionally, we note that most arcs feature
high eccentricities. This is arguably because suborbital arcs always present
subsurface periapsis, as discussed by McMahon et al.[19]

The small body rotation affects the probe arcs by changing the probe’s
velocity in the inertial frame. This effect is more relevant for fast-rotating
bodies and when the probes are located at low latitudes, as shown in Figure 6.
Furthermore, the hopping velocity is always biased toward the east direction.
This effect could be mitigated by actively controlling the hopping velocity
according to the latitude the probe is located at.

A very small fraction of the arc population can result in escape trajectories.
For this data set, about 0.3% of the orbit becomes hyperbolic and escape from
the small body. While the hopping velocity can be easily designed to avoid
direct escape from the body, it has been shown that escape events can also
occur after periapsis passage due to the perturbed dynamical environment [19].
This is more likely for higher-altitude arcs, where perturbations act for a
longer time. In this study, arcs evolving into hyperbolic orbits are not
considered. Future work will assess in detail escape statistics and related
design choices.

Figure 5: 3D visualization of 100 sampled probe arcs (colored trajectories)
around the small body (purple sphere) in the asteroid-centered inertial frame.
Figure 6: Histograms for the simulated arc parameters. Hyperbolic trajectories
are not reported. The red bar indicated the mean value for each bin.

##  5 Probe Measurements

Gravity estimation is performed by tracking the probe trajectories from an
observing mother spacecraft, using an onboard camera that provides relative-
angle measurements. These are the only data processed for OD, i.e., no ground-
based radiometric measurements are considered.

###  5.1 Mother Spacecraft

We assume that the mother spacecraft is hovering relative to the small body,
at a 4 km distance from the small body center and at a zero Sun phase angle
(i.e., the Sun is behind the observer). The altitude is chosen such that the
small body spans approximately 80% of the camera field-of-view width, which is
considered a good trade-off between probe resolution and width of the observed
scene.

The spacecraft position is modeled as a constant, inertial offset with respect
to the small body center. To compensate for the lack of the observer’s
dynamics model (e.g., uncertainties introduced by station-keeping and other
perturbations), we implement a position’s consider covariance which is
processed by the OD filter (as discussed in the filter setup section). Future
work will also tackle trajectory design for the mother spacecraft. The
spacecraft attitude is nadir-pointing, with a white-noise stochastic pointing
error associated to a 20-minute batch interval. The stochastic model is the
same as the one used for stochastic accelerations.

Under the assumption of a low-budget mission, we define the camera model based
on the JPL’s ASTERIA camera, a small yet capable instrument used to perform
astrophysical measurements from aboard a CubeSat[20]. The parameters used for
the camera model are reported in Table 4. The medium field of view enables
global coverage of the small body surface while observing from a safety
distance, whereas the high resolution ensures high-quality astrometry. We
assume a pinhole camera model, neglecting the effect of distortion and other
error sources[21], though periodic camera calibration can be used to model
distortions in practice.

Table 4: Camera model parameters. Parameter | Value | Units
---|---|---
Aperture diameter | 60.7 | mm
Focal length | 85 | mm
Detector dimensions | $2592\times 2192$ | pixels
Pixel size | $6.5\times 6.5$ |  $\mu$m
Detector Field of View | $11.2\times 9.6$ | deg
Measurement Uncertainty | 0.1 | pixels

###  5.2 Probe Visibility

One strength of using bright artificial probes is that they can be detected
not only on the dark background, but also when moving in front of the small
body surface. In this second case, the probes altitude is usually lower (e.g.
see Figure 5) and hence the gravity signal stronger. To assess the quality of
probe detections in front of the surface we consider its photometric and
astrometric properties. The probe’s photometry is given by comparing the
relative magnitude of the photon flux emitted by the probe with the one from
the small body surface. The astrometric uncertainty represents the detection
error given by the spatial distribution and discretization of the probe signal
in the camera plane. Hockman et al. presented a probe visibility analysis for
the Gravity Poppers scenario, considering the same mission scenario presented
in this work, and concluded that both passive probes and strobing probes can
be detected in front of the small body surface with an photometric SNR of 2 or
more. The astrometric error presented in said work is on the order of
hundredth of a pixel; however, it does not account for camera calibration and
other error sources. Hence, we conservatively assume an astrometric error of
0.1 pixel, which matches the measurement uncertainty in the OD filter (Table
4). Future work will further analyze how photometric and astrometric errors
affect the measurement uncertainty for probe detections.

###  5.3 Tracking Frequency

In addition to measurement precision, tracking frequency also plays an
instrumental role in gravity observability. As shown in Equation 1, both the
amplitude and the “wavelength” of the gravity harmonics decrease with
increasing degree. As such, observing higher-degree gravity harmonics requires
(1) taking measurements with high spatial resolution with respect to the
surface and (2) observing the motion at a relatively low altitude.

State-of-the-practice techniques for optical navigation require all navigation
images to be downloaded to the ground, in order to manually process surface
landmarks. Given the limitations in the spacecraft-to-ground link data rate,
there is an upper limit to the achievable imaging frequency, usually in the
order of one image every 10 minutes. As discussed in the results section, such
a tracking frequency greatly limits the amount of observable gravity terms,
even when using low-altitude detections. Since the probe photometric
properties make them easy to detect in front of the small body surface, we
claim that an autonomous detection algorithm (e.g. based on thresholding or
match filtering) could be implemented onboard. This routine could extract the
probe centroids and download them to the ground, in place of the entire image,
for OD and gravity estimation purposes. This process would dramatically reduce
the amount of data volume required, and hence would enable much higher
operational tracking frequencies.

With this reasoning, we assume that our notional spacecraft can sample one
image every 10 seconds. This parameter will also be refined in future work by
including more system-level considerations. For the arc population considered
in this scenario, the difference in spatial resolution between a 10-minute and
a 10-second sampling interval is massive, as shown by the ground tracks in
Figure 7. Importantly, the advantage of increasing the tracking frequency is
twofold: a higher spatial coverage of the gravitational potential and the
increased number of measurements for each arc. Since the apoapsis altitude of
most trajectory is very low, a 10-minute imaging interval would only provide a
handful of measurements for each arc, which limits the observability of high-
degree gravity terms.

Figure 7: Scatter plots of the observed probe positions in the small body-
fixed frame. Top: the entire latitude-longitude domain. Bottom: a magnified
image comparing the 10-minute tracking interval (black triangles) with the
10-second interval (colored dots). The altitude color bar is saturated at 200
meters.

##  6 Covariance Analysis

In order to quantify gravity estimation performance, we perform a covariance
analysis on the probe tracking scenario. To do so, we use the ground-truth
dynamics to simulate perfect measurements of the probes. Each estimated
parameter is associated with an a-priori covariance, and measurement
uncertainty is accounted.

###  6.1 Multi-Arc Filter

We use a multi-arc filtering approach to treat the uncorrelated trajectories
spanning over multiple time intervals[22]. This technique has been effectively
used in flight operations, e.g. to perform particles OD at the asteroid Bennu
and to estimate the Earth and Moon gravity field with unprecedented precision,
during the GRACE and GRAIL missions, respectively[23, 24].

The key principle of the multi-arc filter is that the estimated quantities can
be divided into local and global parameters, where the former are valid only
for a specific OD fit and the latter affect all OD solutions. In this study,
the probe states for each arc can be considered the local parameters, whereas
all the small body-related parameters (gravity field, pole, etc.) are global.
By framing the problem in this way, each arc’s OD can be performed separately.
Then, each OD solution can be used to recursively update the global parameters
estimates and related uncertainties. The estimates of the global parameters
improve and the uncertainties decrease as more arcs are processed. For a
covariance analysis, we are only interested in the parameters uncertainties
and not the estimated values since the measurement residuals are always zero.
In particular, we are interested in the covariance evolution of the gravity
field coefficients. The gravity covariance obtained by processing multiple
arcs is computed using the following fundamental property of information
matrices: if $N$ different information matrices are produced using $N$
different sets of observations, and the observation sets are independent, then
the total information matrix can be written as the sum of each information
matrix. Then, the estimated gravity covariance matrix $\hat{\mathbf{P}}_{G}$
can be written as

$\hat{\mathbf{P}}_{G}=\left(\mathbf{P}_{G,0}^{-1}+\sum_{j=1}^{N}\Lambda_{G}^{j}\right)^{-1}$
(11)

where $\Lambda_{G}^{j}$ is the gravity information matrix given by the $j$-th
set of observations (i.e. by the $j$-th arc) and $\mathbf{P}_{G,0}$ is the
a-priori gravity covariance matrix. Note that this expression is only valid if
the arcs are uncorrelated with each other. Moreover, it is important to avoid
that the a-priori information is counted for every single arc, as this would
bias the information gain toward higher, unrealistic values. The multi-arc
filtering process can be summarized as follows. After defining the a-priori
covariance for the global parameters, for every arc:

  1. 1.

random initial position and velocity, together with the a-priori arc
covariance, is defined (as described in the hopping event section);

  2. 2.

the ground-truth trajectory and related partials are propagated;

  3. 3.

perfect measurements are simulated, based on the ground-truth trajectory;

  4. 4.

the posterior covariance for the local and global parameters is computed;

  5. 5.

the gravity-field covariance is updated based on the information from the
given arc.

###  6.2 Filter Setup

The OD process is performed using a Square Root Information Filter (SRIF)[15].
Table 5 shows the parameters processed by the filter and related a-priori
uncertainties, whose values are assumed very large, compared to the nominal
value, to avoid double-counting the information in the multi-arc filter, as
previously discussed. In this way, the information gain can be considered
data-driven. The probe trajectories are dynamic parameters, whereas all the
small body-related quantities are assumed to be constant values (i.e. bias
parameters). The position uncertainty of the observing spacecraft is modeled
as a 5-meter consider covariance, which is in line with position uncertainty
values obtained during landmark-based navigation in proximity of the small
body, for flight scenarios[10]. The stochastic pointing uncertainty of 0.003
deg is in line with the OSIRIS-REx case. Finally, we assume that some subset
of the probes that are at rest on the surface can be used as body-fixed
landmarks. Due to their relative brightness with respect to the surface,
steady probes can provide high-quality surface-relative measurements to
enhance the gravity estimation. Here, we assume that 10 probes are used as
landmarks for any given arc and we randomly sample the latitude and longitude
of such landmarks in the equatorial band corresponding to
$\phi\in[-20^{o},+20^{o}]$.

Table 5: A-priori uncertainties for estimated parameters Name | Type | A-Priori $\sigma$ | Units
---|---|---|---
Probes pos. | Dynamic | 100 | km
Probes vel. | Dynamic | 100 | m/s
Pole RA | Bias | 100 | deg
Pole Dec | Bias | 100 | deg
Pole PM | Bias | 100 | deg
Pole $\omega$ | Bias | 100 | deg/day
$\mu$ | Bias | 1 | $\mathrm{km}^{3}/\mathrm{s}^{2}$
$\mu_{\mathrm{mascon}}$ | Bias | 1 | $\mathrm{m}^{3}/\mathrm{s}^{2}$
Normalized SH coeff. | Bias | 1 | –
Landmarks | Bias | 100 | m
S/C pos. | Consider | 5 | m
Stoch. Accel. | Stochastic | $\mathrm{1\cdot 10^{-11}}$ | $\mathrm{km/s^{2}}$
Stoch. Pointing | Stochastic | 0.003 | deg

##  7 Gravity Estimation Results

In this section, we present gravity estimation performance for the above
simulated scenario. The quality of gravity estimates depends on several design
parameters and on the duration of the observation campaign. We firstly present
gravimetry results obtained by processing all the 4000 simulated arcs. Second,
we discuss the evolution of gravity information over time and use a back-of-
the-envelope analysis to infer required mission durations. Third, we assess
the sensitivity of performance to key mission and arc parameters.

###  7.1 Gravity Field Observability

Gravimetry results obtained by processing all the simulated arcs are shown in
Figure 8 (left). The metric we use to assess performance is the signal-to-
noise ratio (SNR) of the estimated harmonic coefficients, which is computed as
the coefficient’s ground-truth magnitude over the posterior 1-$\sigma$ for
that same coefficient, provided by the OD solution. Note that the posterior
uncertainty is given by the sum of all information matrices from each
processed arc, as shown in Equation 11. We separate the SNR of the zonal
coefficients from the remaining (tesseral and sectoral) coefficients, as the
former are systematically larger (which can be deduced by the parameters in
Table 3). The tesseral and sectoral terms (here defined non-zonal) are
considered as a whole using their RMS value.

We compare our gravimetry results with those obtained during the OSIRIS-REx
mission, from both the particle tracking campaign and the standard spacecraft-
tracking OD during the Orbital-B mission phase[10]. Results suggest that,
using the proposed approach, the SNR for the zonal harmonic terms increases by
about two orders of magnitude with respect to the OSIRIS-REx’s particle
tracking campaign, which could estimate the gravity field up to about degree
10, and by three orders of magnitude when compared to OSIRIS-REx’s spacecraft
tracking, which only estimated the $J_{2}$ term. It is important to note that
these three cases differ in a number of ways, such as the amount and type of
processed measurements, hence the comparison can only be made qualitatively.
Finally, the mascon’s SNR is 107, proving that even surface features such as
large boulders, and potentially subsurface voids, can be observed with low
uncertainties.

Since the SNR of all simulated coefficients is always above the observability
threshold (i.e., SNR$>$1) for all 20 degrees and orders, we fit a model to our
SNR data to assess observability for higher-degree and order coefficients. The
numerator of the SNR model is given by the Kaula’s harmonic magnitudes
(Equation 2). For the denominator, we fit an exponential function $\alpha
e^{\beta x}$ to the posterior 1-$\sigma$ data. The estimated parameters are
reported in Table 6. This regression approach suggests that, under these
assumptions, the zonal harmonic coefficients are observable up to about degree
60, whereas the remaining coefficients are observable up to degree 40,
approximately. The residuals between the fitted model and the actual SNR
obtained via OD is shown in Figure 8 (right). The residual values are
normalized with respect to each degree’s estimated SNR (i.e. using the
harmonic value estimated by OD) and have mean of 0.16 and -0.13 for zonal and
non-zonal terms, respectively. The model appear to be valid to infer crude
estimates of the highest-degree observable harmonics; however, future work
will validate and refine this model by simulating higher-degree and order
gravity fields.

Table 6: Estimated parameters of the exponential model $\alpha e^{\beta x}$ for the posterior 1-$\sigma$ evolution. | $\alpha$ | $\beta$
---|---|---
Zonal | $2.0975\times 10^{-6}$ | $3.6937\times 10^{-2}$
RMS | $2.4817\times 10^{-6}$ | $4.2667\times 10^{-2}$

Figure 8: SNR of the estimated harmonic coefficients. Left: comparison between
the coefficients SNR estimated in this work (blue and orange) and by the
OSIRIS-REx mission operations using particle tracking (black) and spacecraft
tracking (purple). The dots are the SNR obtained by OD and the dashed lines
represent the estimated model for SNR evolution. Right: residuals between the
estimated model for SNR evolution and the actual SNR obtained by OD.

###  7.2 Mission Duration

To implement the proposed gravity estimation concept in an actual mission,
assessing requirements on the mission duration is paramount. Usually,
proximity operations around a small body require continuous ground-based
support, navigation, and planning, hence the probe observation time frame is
strictly linked to the overall mission cost and complexity. Since our OD
pipeline evaluates the probe arcs separately, the underlying mission duration
is not explicitly defined. However, we can use the arc population statistics
and a back-of-the-envelope approach to roughly quantify the mission duration.

We start by analyzing the gravity information gain as a function of the number
of arcs. We fit the SNR evolution with a power law model $\alpha x^{\beta}$,
as shown in Figure 9 (left). For example, the estimated power-law parameters
for $J_{2}$ are $\alpha=64.4709$ and $\beta=0.5590$. Note that there is a
transient state in SNR in the first hundreds of arcs, which is likely due to
the fact that the surface coverage with few arcs is sparse and the information
gain for each arc is more dependent upon the location of the arc itself. As
such, we fit the SNR evolution starting from the 2000-th arc onward. The model
residuals are lower than 2% of the SNR value (Figure 9, right). However,
future work will simulate more arcs to validate this model.

Figure 9: SNR evolution for some representative gravity harmonics. Left:
actual evolution compared to the fitted model. Right: residuals of the actual
values compared to the model.

The power-law model and its derivative, shown in Figure 9 (left) suggest that
there is a diminishing return in the information gain when the arc number
increases. We use this relationship to assess the required mission duration as
a function of the maximum observable zonal harmonic term, as shown in Figure
10. In other words, we assume a given gravity science objective (e.g.
observing the zonal harmonics up to degree 40) and estimate a mission duration
needed to achieve the goal. The mission duration $T_{\mathrm{mission}}$ is
evaluated as:

$T_{\mathrm{mission}}=\frac{\bar{T}_{\mathrm{arc}}N_{\mathrm{arc}}}{f_{\mathrm{vis}}n_{\mathrm{probe}}}$
(12)

where $\bar{T}_{\mathrm{arc}}=69$ minutes is the mean arc lifetime (Figure 6),
$f_{\mathrm{vis}}=0.5$ is a visibility factor assuming that only half of the
arcs are visible from the mother spacecraft, at any given time (due to the
body occlusion; note that this is a conservative assumption);
$n_{\mathrm{probe}}$ indicates the number of hopping probes utilized during
the mission (which does not account for probes potentially used as steady
surface landmarks). Results show that good estimates (i.e. SNR$>$10) up to
$J_{40}$ can be obtained within a period spanning from a few months to about
one year, depending on the number of probes used. On the other hand, the
mission time increases as a power law with the target degree to observe, hence
estimating higher-degree terms may not be done within reasonable time frames.
In addition, there is little advantage in increasing the number of hopping
probes above 20-30.

Figure 10: Mission duration as a function of the maximum observable zonal
harmonic coefficient J. The plot is parametrized with the number of hopping
probes.

###  7.3 Sensitivity To Measurements

The presented mission scenario is based on a number of assumptions on the
observation type and performance. In this section, we assess gravimetry
performance for cases where such assumptions are relaxed or new types of
observations are introduced. We simulate 1000 arcs for 5 different case
studies: (1) the baseline simulation setup (described in the previous
sections), (2) increasing the imaging rate to one image every 10 minutes, (3)
increasing the optical measurement uncertainty to 0.5 pixels, (4) including
ranging and Doppler measurements in the OD process, and (5) increasing the
camera pointing uncertainty. Case 2 is used to simulate performance when the
probes (or particles) are detected on the ground, and hence the imaging rate
is limited by the downlink data volume (e.g. OSIRIS-REx particle tracking).
Case 3 assumes that the photometric properties of the observed targets are
poorer than the artificial probes’, e.g. when tracking natural particles are
tracked. In case 4, we treat the artificial probes as Two-Way coherent
transponders and assume that a Deep-Space Atomic Clock (DSAC) is used aboard
the mother spacecraft[25]. We assume $\sigma_{\rho}=10\,\mathrm{cm}$ and
$\sigma_{\dot{\rho}}=0.1\,\mathrm{mm/s}$ for the range and range-rate
measurement uncertainty, respectively. We model such measurements as perfect
range and range rate. Note that, for this measurement strategy to be adopted,
a more advanced probe design is required to enable the spacecraft-to-probe
link.

For each case, we fit the harmonic SNR with the exponential model previously
introduced; results are shown in Figure 11. While all scenarios accurately
recover the low-degree harmonics, the maximum observable degree varies widely
among cases. As expected, introducing range and range-rate measurements
outperforms all other scenarios, observing the zonal harmonics up to about
degree 70, compared to degree 50 from the baseline. On the other hand,
reducing the imaging frequency and measurement uncertainty decreases the
maximum observed harmonic below 25 and 30, respectively. These conclusions
highlight the importance of autonomous onboard probe detection and good
photometric properties of the targets. Finally, gravimetry performance is not
particularly sensitive to the magnitude of the stochastic pointing errors, and
remain almost unchanged when the error is increased by an order of magnitude.
Future work will include more thorough sensitivity analyses that explore a
broader domain of the trade space for mission design.

Figure 11: Fitted SNR model for the estimated zonal harmonics, when different
observation parameters and when adding range and Doppler measurements.

###  7.4 Sensitivity to Probe Trajectories

Finally, we analyze the sensitivity of gravimetry performance to some
representative arc’s orbital parameters. We selected the arcs lifetime,
hopping velocity (both in the inertial and body-fixed frame), the altitude of
the orbital apoapsis, the eccentricity, and the initial latitude. The
information gain is represented for a given arc and a given harmonic degree,
and is computed as the RMS of all gravity SNR for that degree (accounting for
zonal and non-zonal coefficients with such a degree). Figure 12 presents
results for degrees 2 and 20. As expected, gravimetry for the lower degree
(which has a lower spatial resolution) is better for higher altitudes and
hopping velocities, whereas the opposite is true for higher-degree terms. We
also observe that the information gain variation for low-degree harmonics is
much larger for low-degree terms, compared to higher degrees. This might be
due to higher-altitude orbits being more diverse, as the orbital perturbations
act on a longer time of flight. However, excluding some outliers, orbital
parameters have little impact on gravity information gain for the vast
majority of the arc population. The hopping velocity can be calibrated, hence
it is arguably the most interesting parameter from a mission design
standpoint. However, there is little advantage in fine-tuning this parameter
to a specific value as the correlation with gravimetry performance is weak,
especially for the (most interesting) high-degree gravity terms.

Figure 12: Gravity information gain as a function of selected orbital
parameters, for degree 20 (top) and degree 2 (bottom). The red bar indicates
the mean for the given data bin.

##  8 Conclusions

In this study, we assess gravity estimation performance at small bodies for a
novel mission architecture. We reproduced a scenario similar to the OSIRIS-REx
mission, around a low-mass small body such as asteroid Bennu. We used
conservative assumptions on unmodeled forces and measurement uncertainty, but
assumed that autonomous onboard probe detection is available to enable high
tracking frequency. Results suggest that it is possible to recover the gravity
field of the low-mass body up to degree 40 within days or months, depending on
the target SNR, using 20 hopping probes. The mission time frame largely
depends on the gravity science objective. For the asteroid Bennu, the gravity
SNR outperforms spacecraft-tracking techniques by three orders of magnitude.
The gravity contribution from surface features, such as large boulders, can
also be precisely estimated. High measurement frequency and good probe
photometric properties are key to achieve high gravimetry performance within a
limited time frame, and spacecraft-to-probe radiometric measurements
substantially enhanced the gravity estimates. On the other hand, gravimetry is
not particularly sensitive to the orbital parameters of the hopping probes.

Future work will increase simulation realism, e.g. including the mother
spacecraft orbital dynamics, using non-spherical models for the target body
shape, and interior gravity models. Additionally, we will extend sensitivity
analyses to a broader trade space of mission design and measurement
strategies. Finally, we will implement a density estimation scenario to assess
how this gravimetry technique could be applied for specific planetary science
objectives.

##  9 Acknowledgements

A portion of this research was carried out at the Jet Propulsion Laboratory,
California Institute of Technology, under a contract with the National
Aeronautics and Space Administration. ©2021 California Institute of
Technology. Government sponsorship acknowledged.

## References

  * [1] D. Scheeres, D. Britt, B. Carry, and K. Holsapple, “Asteroid interiors and morphology,” Asteroids iv, Vol. 745766, 2015.
  * [2] P. Michel, F. E. DeMeo, and W. F. Bottke, Asteroids iv.  University of Arizona Press, 2015.
  * [3] V. Ciarletti, A. C. Levasseur-Regourd, J. Lasue, C. Statz, D. Plettemeier, A. Hérique, Y. Rogez, and W. Kofman, “CONSERT suggests a change in local properties of 67P/Churyumov-Gerasimenko’s nucleus at depth,” Astronomy & Astrophysics, Vol. 583, 2015, p. A40.
  * [4] L.-I. Sorsa, M. Takala, P. Bambach, J. Deller, E. Vilenius, and S. Pursiainen, “Bistatic full-wave radar tomography detects deep interior voids, cracks, and boulders in a rubble-pile asteroid model,” The Astrophysical Journal, Vol. 872, No. 1, 2019, p. 44.
  * [5] M. Haynes, A. K. Virkki, F. C. Venditti, D. Hickson, N. Pinilla-Alonso, J. Brisset, L. A. Benner, C. A. Raymond, J. Lazio, A. Freeman, et al., “Asteroids Inside Out: Radar Tomography,” 2020.
  * [6] D. Scheeres, E. Asphaug, J. Colwell, R. Dissly, P. Geissler, L. McFadden, V. Petr, R. Reinert, and H. Yano, “Asteroid surface science with pods,” Lunar and Planetary Science Conference, 2003, p. 1444.
  * [7] J. Plescia, O. Barnouin, D. Richardson, N. Schmerr, D. Lawrence, B. Denevi, C. Ernst, and H. Yu, “The asteroid probe experiment (apex) mission,” 2017.
  * [8] G. Nishiyama, T. Kawamura, N. Namiki, B. Fernando, K. Leng, K. Onodera, S. Sugita, T. Saiki, H. Imamura, Y. Takagi, et al., “Simulation of Seismic Wave Propagation on Asteroid Ryugu Induced by The Impact Experiment of The Hayabusa2 Mission: Limited Mass Transport by Low Yield Strength of Porous Regolith,” Journal of Geophysical Research: Planets, Vol. 126, No. 2, 2021, p. e2020JE006594.
  * [9] S. Chesley, A. French, A. Davis, R. Jacobson, M. Brozović, D. Farnocchia, S. Selznick, A. Liounis, C. Hergenrother, M. Moreau, et al., “Trajectory estimation for particles observed in the vicinity of (101955) Bennu,” Journal of Geophysical Research: Planets, Vol. 125, No. 9, 2020, p. e2019JE006363.
  * [10] J. M. Leonard, J. L. Geeraert, B. R. Page, A. S. French, P. G. Antreasian, C. D. Adam, D. R. Wibben, M. C. Moreau, and D. S. Lauretta, “OSIRIS-REx orbit determination performance during the navigation campaign,” 2019 AAS/AIAA Astrodynamics Specialist Conference, 2019, pp. 1–20.
  * [11] S. Evans, W. Taber, T. Drain, J. Smith, H.-C. Wu, M. Guevara, R. Sunseri, and J. Evans, “MONTE: The next generation of mission design and navigation software,” CEAS Space Journal, Vol. 10, No. 1, 2018, pp. 79–86.
  * [12] E. Mazarico, D. D. Rowlands, T. J. Sabaka, K. M. Getzandanner, D. P. Rubincam, J. B. Nicholas, and M. C. Moreau, “Recovery of Bennu’s orientation for the OSIRIS-REx mission: implications for the spin state accuracy and geolocation errors,” Journal of Geodesy, Vol. 91, No. 10, 2017, pp. 1141–1161.
  * [13] D. Scheeres, J. McMahon, A. French, D. Brack, S. Chesley, D. Farnocchia, Y. Takahashi, J. Leonard, J. Geeraert, B. Page, et al., “The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements,” Nature Astronomy, Vol. 3, No. 4, 2019, pp. 352–361.
  * [14] W. M. Folkner, J. G. Williams, D. H. Boggs, R. S. Park, and P. Kuchynka, “The planetary and lunar ephemerides DE430 and DE431,” Interplanetary Network Progress Report, Vol. 196, No. 1, 2014.
  * [15] B. Schutz, B. Tapley, and G. H. Born, Statistical orbit determination.  Elsevier, 2004.
  * [16] W. M. Kaula, Theory of satellite geodesy: applications of satellites to geodesy.  Courier Corporation, 2013.
  * [17] J. McMahon, D. Scheeres, S. Hesar, D. Farnocchia, S. Chesley, and D. Lauretta, “The OSIRIS-REx radio science experiment at Bennu,” Space Science Reviews, Vol. 214, No. 1, 2018, pp. 1–41.
  * [18] D. J. Scheeres, Orbital motion in strongly perturbed environments: applications to asteroid, comet and planetary satellite orbiters.  Springer, 2016.
  * [19] J. W. McMahon, D. J. Scheeres, S. R. Chesley, A. French, D. Brack, D. Farnocchia, Y. Takahashi, B. Rozitis, P. Tricarico, E. Mazarico, et al., “Dynamical evolution of simulated particles ejected from asteroid Bennu,” Journal of Geophysical Research: Planets, Vol. 125, No. 8, 2020, p. e2019JE006229.
  * [20] C. M. Pong and M. W. Smith, “Camera modeling, centroiding performance, and geometric camera calibration on ASTERIA,” 2019 IEEE Aerospace Conference, IEEE, 2019, pp. 1–17.
  * [21] W. M. Owen Jr, “Methods of optical navigation,” 2011.
  * [22] A. S. French, Precise Orbit Determination and Gravity Field Estimation During Small Body Missions.  PhD thesis, University of Colorado at Boulder, 2020.
  * [23] A. S. French, J. M. Leonard, J. L. Geeraert, B. R. Page, P. G. Antreasian, M. C. Moreau, J. W. McMahon, D. J. Scheeres, and D. S. Lauretta, “Multi-arc filtering during the navigation campaign of the osiris-rex mission,” AAS/AIAA Astrodynamics Specialist Conference, 2019, Univelt Inc., 2020, pp. 3051–3068.
  * [24] S. W. Asmar, A. S. Konopliv, M. M. Watkins, J. G. Williams, R. S. Park, G. Kruizinga, M. Paik, D.-N. Yuan, E. Fahnestock, D. Strekalov, et al., “The scientific measurement system of the Gravity Recovery and Interior Laboratory (GRAIL) mission,” GRAIL: Mapping the Moon’s Interior, 2013, pp. 25–55.
  * [25] T. A. Ely, E. A. Burt, J. D. Prestage, J. M. Seubert, and R. L. Tjoelker, “Using the deep space atomic clock for navigation and science,” IEEE transactions on ultrasonics, ferroelectrics, and frequency control, Vol. 65, No. 6, 2018, pp. 950–961.


# Energy Spectrum of a constrained Quantum Particle
and the Willmore Energy of the constraining Surface

Vicent Gimeno i Garcia<EMAIL_ADDRESS>Departament de Matemàtiques- IMAC,
Universitat Jaume I, Castelló, Spain.  Steen Markvorsen<EMAIL_ADDRESS>DTU
Compute, Technical University of Denmark, Lyngby, Denmark

###### Abstract

Geometric and topological bounds are obtained for the first energy level gap
of a particle constrained to move on a compact surface in 3-space. Moreover,
geometric properties are found which allows for stationary and uniformly
distributed wave functions to exist on the surface.

PACS number(s): 03.65.Ge, 68.65.2k

††preprint: APS/123-QED

##  I Introduction

During the first and second decades of this century the study of surface
physics has attracted attention from both theoretical and experimental
perspectives. Nanostructures, such as fullerenes, carbon-based low-dimensional
shapes, and nanotubes have been researched due to their potential
technological applications (see references [1, 2, 3, 4, 5, 6, 7, 8, 9]).

From a theoretical perspective, as developed e.g. in [10], the dynamics and
energy levels of a quantum particle on a surface can be described as a
particle in the curved space defined by the surface, and can thus be viewed
from a purely intrinsic perspective. However, when considering a particle that
is confined to an immersed surface in a locally Euclidean 3-space, an
extrinsic perspective can also be adopted.

In fact, according to [11] and [12], an extrinsically motivated effective
potential is present in the Schrödinger equation for a surface in 3-space.
That potential involves both the Gaussian curvature (i.e. an intrinsic
property of the surface) and the mean curvature (an extrinsic property of the
surface). This idea of inducing a combination of extrinsic and intrinsic
contributions into the governing equations for the quantum phenomena on a
surface can be further extended, for example by adding electromagnetic fields
or in the setting of Quantum Field Theory (as discussed in the following
references [13, 14, 15, 16, 17, 18, 19] for instance).

Da Costa in [12], and first Jensen-Koppe in [11], considered the classical
Schrödinger equation for a quantum scalar particle confined to a surface
$\Sigma$ in $\mathbb{R}^{3}$ and obtained an equation for the particle which
contains an effective potential expressed in terms of the mean curvature and
the Gaussian curvature of the surface. More precisely, the classical
Schrödinger equation was first considered in a tubular neighbourhood of the
the surface in $\mathbb{R}^{3}$:

$\Sigma_{\epsilon}=\\{x\in\mathbb{R}^{3}\,:\,{\rm
dist}(x,\Sigma)\leq\epsilon\\}.$

Assuming that the surface $\Sigma$ is compact, regular and embedded without
boundary the Schrödinger equation is – for sufficiently small $\epsilon$ –
well defined in the tubular 3-dimensional domain $\Sigma_{\epsilon}$, when
this domain is equipped with Dirichlet boundary conditions. After shrinking
$\epsilon\to 0$, da Costa then obtained the following Schrödinger equation for
the particle confined in the surface $\Sigma$,

$i\hbar\frac{\partial}{\partial
t}\Psi=\frac{\hbar^{2}}{2m}\left[-\Delta-\frac{1}{4}\left(k_{1}-k_{2}\right)^{2}\right]\Psi,$
(1)

where $k_{1}$ and $k_{2}$ denote the principal curvatures of $\Sigma$. An
effective potential $V_{\rm eff}$ appears in this way with the following
alternative geometric expression

$V_{\rm
eff}=-\frac{\hbar^{2}}{2m}\left(H^{2}-K\right)=-\frac{\hbar^{2}}{8m}S^{2},$

where $H=(k_{1}+k_{2})/2$ is the mean curvature, $K=k_{1}k_{2}$ is the
Gaussian curvature and $S=k_{1}-k_{2}$ is the skew curvature. The skew
curvature measures how far the surface is from being umbilical. See the
introduction in [20] for a well motivated description of the skew curvature
and [21] for the extension of the skew curvature for surfaces in the Lorentz-
Minkowski space and the interpretation of the skew curvature in terms of the
standard deviation of the normal curvature seen as a random variable.

Although equation (1) is commonly seen in the context of embedded, compact
regular surfaces, where it has direct and significant physical implications,
it is also an interesting and mathematically well-defined problem to examine
the development of wave functions in terms of equation (1) for compact regular
surfaces that are merely immersed – not necessarily embedded. Thus, in this
paper, we will consider general immersed surfaces with the Schrödinger
equation (1).

To find the energy spectrum of the particle, the eigenvalue problem for the
corresponding Hamiltonian has to be solved, i.e. we have to analyze the family
of stationary wave functions $\\{\Psi_{k}\\}_{k\in\mathbb{N}}$ and the energy
spectrum $\\{E_{k}\\}_{k\in\mathbb{N}}$ such that

$\frac{\hbar^{2}}{2m}\left[-\Delta-\frac{1}{4}S^{2}\right]\Psi_{k}=E_{k}\Psi_{k}.$
(2)

In this report we focus on the problem of obtaining upper bounds for the gap
in the energy spectrum (see inequalities (3) and (7)) using topological and
conformal invariants for the surfaces in question. The paper is organized as
follows: in Sect. II, the Willmore energy and the Euler characteristic will be
explained. These fundamental concepts will be the conformal and topological
invariants that we will use to obtain the upper bounds for the gap in the
spectrum of the energy of a particle confined in a given surface. In Sect.
III, we prove inequality (3), which gives an upper bound for the difference of
the two first bounded energy levels $E_{0}$ and $E_{1}$. Sect. IV is devoted
for the study of $E_{k}-E_{0}$ for $k>1$ and we obtain the qualitative
inequality (7). In Sect. V, we give a simple construction and description of
the family of surfaces of revolution with constant skew curvature. In Sect.
VI, we show how our techniques and results can be extended to non-compact
surfaces when they satisfy certain symmetry conditions. Sect. VII is devoted
to concluding remarks.

##  II Willmore energy and Euler Characteristic

The Willmore energy of an immersed surface is a conformal invariant (see [22])
and is defined as the integral of the square of the mean curvature over the
surface,

$\mathcal{W}(\Sigma)=\int_{\Sigma}H^{2}dA.$

It is known to be always greater than or equal to $4\pi$, with equality only
for the round spheres. Moreover, the Gauss-Bonnet Theorem relates the integral
of the Gaussian curvature over the surface to the Euler characteristic, which
is a topological invariant in the following way:

$\int_{\Sigma}K\,dA=2\pi\chi(\Sigma).$

The Euler characteristic $\chi(\Sigma)$ of an oriented and compact surface
without boundary is either equal to $2$ (for the sphere and for any surface
homeomorphic to the sphere) or equal to $2-2g$, where $g$ is the genus of the
surface.

##  III First Gap in the spectrum of energy and Willmore energy

In this section we provide an upper bound for the difference of the first and
second energies of a quantum particle in a two-dimensional surface with mean
curvature $H$, Gauss curvature $K$, and skew curvature $S$. The Hamiltonian
$\frac{\hbar^{2}}{2m}\left(-\Delta-V\right)$, with $V=H^{2}-K$, describes the
system and its spectrum of energies will be denoted by
$\\{E_{i}=\frac{\hbar^{2}}{2m}\lambda_{i}\\}_{i\in\mathbb{N}}$. Equation (2)
can be expressed as

$\left[-\Delta-V\right]\Psi=\frac{2mE}{\hbar^{2}}\Psi,\quad V=H^{2}-K.$

Let $\\{\lambda_{i}\\}_{i\in\mathbb{N}}$ be the spectrum of the operator
$-\Delta-V$, and hence, $E_{i}=\frac{\hbar^{2}}{2m}\lambda_{i}$. By Theorem
2.1 of [23],

$\displaystyle\lambda_{1}\leq$ $\displaystyle\frac{1}{{\rm
A}(\Sigma)}\int_{\Sigma}\left(2H^{2}-H^{2}+K\right)d{\rm A}$ $\displaystyle=$
$\displaystyle\frac{2\mathcal{W}(\Sigma)}{{\rm
A}(\Sigma)}-\frac{1}{4}\overline{S^{2}},$

where $\lambda_{1}$ is the second eigenvalue of the operator $-\Delta-V$,
${\rm A}(\Sigma)$ is the area of $\Sigma$, and
$\overline{S^{2}}:=\frac{\int_{\Sigma}S^{2}d{\rm A}}{{\rm A}(\Sigma)}$ is the
mean of the squared skew curvature of $\Sigma$. The first eigenvalue
$\lambda_{0}$ of the operator $-\Delta-V$ can also be estimated from below.
Let $u$ be an eigenfunction associated with $\lambda_{0}$. Then using the
Rayleigh quotient we get

$\displaystyle\lambda_{0}=$ $\displaystyle\frac{-\int_{\Sigma}u\Delta
udA-\int_{\Sigma}Vu^{2}dA}{\int_{\Sigma}u^{2}dA}$ $\displaystyle=$
$\displaystyle\frac{\int_{\Sigma}\|\nabla
u\|^{2}dA}{\int_{\Sigma}u^{2}dA}-\frac{\int_{\Sigma}Vu^{2}dA}{\int_{\Sigma}u^{2}dA}$
$\displaystyle\geq-\max_{\Sigma}V=-\frac{1}{4}\max_{\Sigma}S^{2}.$

Finally, since $E_{i}=\frac{\hbar^{2}}{2m}\lambda_{i}$, we can state the
following upper bound for the energy difference:

$\displaystyle E_{1}-E_{0}\leq$
$\displaystyle\frac{\hbar^{2}}{2m}\left(\frac{2\mathcal{W}(\Sigma)}{{\rm
A}(\Sigma)}+\frac{1}{4}\left(\max_{\Sigma}S^{2}-\overline{S^{2}}\right)\right).$
(3)

To achieve this upper bound, it is essential that each inequality employed in
its derivation is transformed into an equality. Specifically, a necessary
condition for achieving this upper bound is that the surface must be a surface
with constant skew curvature. Additionally, the upper bound given by equation
(3) is sharp (the equality is achieved for at least one surface). This is
because for a sphere with a radius of $R$ in $\mathbb{R}^{3}$, it is known
that $\mathcal{W}=4\pi$, ${\rm A}(S_{R})=4\pi R^{2}$, and $S=0$. Therefore,
our upper bound given in inequality (3) can be expressed as follows:

$0\leq E_{1}-E_{0}\leq\frac{\hbar^{2}}{2m}\cdot\frac{2}{R^{2}}$

However, it should be noted that in this particular case, since $H^{2}-K=0$,
the spectrum is equivalent to the Laplacian spectrum. Hence, according to
[24], we know that: $E_{1}-E_{0}=\frac{\hbar^{2}}{2m}\frac{2}{R^{2}}.$ As far
as we know it is still unknown if inequality (3) becomes an equality for other
surfaces aside from the round sphere. Nevertheless, as we mentioned before, a
prerequisite for achieving equality is that the surface must have constant
skew curvature.

We must notice that $\max_{\Sigma}S^{2}-\overline{S^{2}}\geq 0$ with equality
only for surfaces of constant skew curvature and moreover by using the Gauss-
Bonnet Theorem we have

$\frac{1}{4}\overline{S^{2}}=\frac{\int_{\Sigma}H^{2}-K}{{\rm
A}(\Sigma)}=\frac{\mathcal{W}(\Sigma)-2\pi\chi(\Sigma)}{{\rm A}(\Sigma)}.$

Therefore, for orientable surfaces:

$\frac{1}{4}\overline{S^{2}}=\frac{\mathcal{W}(\Sigma)-4\pi(1-g)}{{\rm
A}(\Sigma)}\cdot$

##  IV Further gaps in the spectrum

To provide an upper bound for the difference between consecutive eigenvalues
$E_{k-1}$ and $E_{0}$ we will analyze the number of non-positive eigenvalues
$\mathcal{N}(V_{\alpha})$ of the Schrödinger operator $-\Delta-V_{\alpha}$,
with $V_{\alpha}=H^{2}-K+\alpha$, with $\alpha$ being a real constant.

The number of non-positive eigenvalues of the Schrödinger operator $-\Delta-
V_{\alpha}$ have been bounded from below by

$\mathcal{N}(V_{\alpha})\geq c(g)\int_{\Sigma}V_{\alpha}dA$

in [25], where $c(g)$ depends only on the genus of $\Sigma$. By choosing
$\alpha$ such that the lower bound is $k$, _i.e._ ,

$c(g)\int_{\Sigma}V_{\alpha}dA=k,$ (4)

the first $k$-th eigenvalues $\\{\lambda_{i}(V_{\alpha})\\}_{i=0}^{k-1}$
associated with $-\Delta-V_{\alpha}$ are non-positive

$\lambda_{0}(V_{\alpha})\leq\lambda_{1}(V_{\alpha})\leq\cdots\leq\lambda_{k-1}(V_{\alpha})\leq
0.$

Therefore an upper bound of $\lambda_{k-1}-\lambda_{0}$ can now be derived and
expressed in terms of $\alpha$ and $\overline{S^{2}}$ because

$0\leq\lambda_{k-1}(V_{\alpha})-\lambda_{0}(V_{\alpha})\leq-\lambda_{0}(V_{\alpha})=-\lambda_{0}+\alpha.$
(5)

But applying equation (4)

$\displaystyle k=c(g)\int_{\Sigma}V_{\alpha}dA=$ $\displaystyle
c(g)\left(\int_{\Sigma}VdA+\alpha{\rm A}(\Sigma)\right)$ $\displaystyle=$
$\displaystyle c(g)\left(\frac{1}{4}\int_{\Sigma}S^{2}d{\rm A}+\alpha{\rm
A}(\Sigma)\right).$

Therefore, the parameter $\alpha$ can be obtained as

$\alpha=\frac{k}{{\rm A}(\Sigma)c(g)}-\frac{1}{4}\overline{S^{2}}.$ (6)

By using this expression for the value of $\alpha$, the upper bound of
inequality (5) can be written as

$\displaystyle\lambda_{k-1}-\lambda_{0}\leq$
$\displaystyle-\lambda_{0}+\frac{k}{{\rm
A}(\Sigma)c(g)}-\frac{1}{4}\overline{S^{2}}.$

But we can estimate $\lambda_{0}$ from below as before and get

$\displaystyle\lambda_{0}$ $\displaystyle\geq-\frac{1}{4}\max_{\Sigma}S^{2}.$

Then,

$E_{k-1}-E_{0}\leq\frac{\hbar^{2}}{2m}\left(\frac{k}{{\rm
A}(\Sigma)c(g)}+\frac{1}{4}\left(\max_{\Sigma}S^{2}-\overline{S^{2}}\right)\right).$
(7)

Although the theoretical constant $c(g)$ has not yet been determined, and its
exact value is unknown, it is not possible to determine the equality case.
However, we do know that when equality is achieved, as in the case of
inequality (3), the surface must be a surface of constant skew curvature.

Inequality (7) leads to the Weyl’s type asymptotic formula

$\liminf_{k\to\infty}\frac{E_{k}-E_{0}}{k}\leq\frac{\hbar^{2}}{2m}\frac{1}{c(g)A(\Sigma)}\cdot$

We have to remark here that the energy levels are continuous (non-discrete
energy levels) in the classical limit of $\hbar\to 0$ or in the high mass
limit of $m\to\infty$, and for any given compact surface
$\Sigma\subset\mathbb{R}^{3}$ if we consider the re-scaled family of surfaces
$\\{\lambda\Sigma\\}_{\lambda>0}$, then

$\lim_{\lambda\to\infty}E_{k}(\lambda\Sigma)-E_{0}(\lambda\Sigma)=0.$

In short, the classical limit is obtained at large mass or large scale – as
expected.

##  V Surfaces of constant skew curvature

The induced Schrödinger (1) equation leads to the study of surfaces of
constant skew curvature. From a physical point of view, a constant wave
function, $\Psi=c$, is a stationary solution of the Schrödinger equation if
there is a value of $E$ such that $-\frac{1}{4}S^{2}c=\frac{2mE}{\hbar^{2}}c$,
which means that the surface has constant skew curvature. For a compact
surface, if the particle is initially in a uniformly distributed wave
function, $\|\Psi\|^{2}=\frac{1}{{\rm A}(\Sigma)}$, it will remain so if it is
an eigenfunction of the Hamiltonian:

$\|\nabla\Theta\|^{2}+i\Delta\Theta=\frac{1}{4}S^{2}+\frac{2mE}{\hbar^{2}}.$

The necessary condition for this to occur is that the surface has constant
skew curvature. Because in a uniformly distributed wave function, we have a
real function $\Theta:\Sigma\to\mathbb{R}$ such that

$\Psi=\frac{e^{i\Theta}}{\sqrt{{\rm A}(\Sigma)}}\cdot$

Then

$\|\nabla\Theta\|^{2}+i\Delta\Theta=\frac{1}{4}S^{2}+\frac{2mE}{\hbar^{2}},$

which implies that $\Theta$ is a harmonic function ($\Delta\Theta=0$). But the
only harmonic functions in a compact surface is the constant function, and
thus the surface must have constant skew curvature. On the other hand, if the
compact surface has constant skew curvature, the geometrically induced
potential in this case is constant, leading to the spectrum of energies
coinciding with the pure Laplace spectrum, with the ground state being the
uniformly distributed wave function itself. To obtain a non-uniformly
distributed wave function, a certain amount of energy must be supplied, the
size of which is determined by our results. Indeed, from (3) and assuming
constant skew curvature $S=c$ on $\Sigma$ and $\Sigma$ to be an oriented
compact surface of genus $g$, and by using the Gauss-Bonnet theorem we have

$\frac{\hbar^{2}}{m}\frac{\mathcal{W}(\Sigma)}{{\rm
A}(\Sigma)}=\frac{\hbar^{2}}{m}\left(c^{2}+\frac{4\pi(1-g)}{A(\Sigma)}\right)\geq
E_{1}-E_{0}.$ (8)

The bounds in the gap spectrum of energy are improved when the surface has
constant skew curvature. The gap in the spectrum of energy is given by the
difference between the first two energy levels and can be bounded by the ratio
of the Willmore functional to the area of the surface. In this specific case,
with constant geometrically induced potential, since the energy spectrum can
be deduced from the Laplace spectrum, by using Theorem 1 of [26], we can
provide the following improvement of inequality (8):

$E_{1}-E_{0}\leq\frac{\hbar^{2}}{mA(\Sigma)}\cdot\min\left\\{4\pi(1+g),\mathcal{W}(\Sigma)\right\\}.$
(9)

In our main results for upper bounds in the gaps in the spectrum of energy,
inequalities (3) and (7), we always have the term

$\frac{1}{4}\left(\max_{\Sigma}S^{2}-\overline{S^{2}}\right),$

which vanishes exactly when the surface has constant skew curvature. What we
have stated then shows that surfaces of constant skew curvature play an
important role in understanding the behavior of a particle confined to a
surface. They are the only surfaces where stationary solutions with uniformly
distributed wave function are allowed and moreover, they are the unique
surfaces where equality could be attained for our upper bounds for the first
gap in the energy spectrum.

Surfaces of constant skew curvature have attracted much attention, see [27],
[28] and [20] for instance. In fact, R. López and A. Pámpano have recently
classified the surfaces of revolution with constant skew curvature, see [20].
One important result from this classification is that there are no other
regular and compact surfaces of revolution with constant skew curvature than
the round spheres.

This motivates the following pertinent open problem:

> Are the round spheres the only compact regular surfaces without boundary and
> with constant skew curvature that are immersed in Euclidean $3$-space?

This geometric problem is similar to the classical geometric problems of
constant curvature. For example, in 1853, J.H. Jellet demonstrated that if a
compact, star-shaped surface in 3-dimensional space, $\mathbb{R}^{3}$, has a
constant mean curvature, it is isometric to a standard round sphere.
Alexandrov later confirmed that any compact, embedded surface in the Euclidean
$3$-dimensional space with a non-zero constant mean curvature must also be a
sphere. The assumption on the embeddedness cannot be removed. Indeed, in 1984,
Henry C. Wente created the Wente torus, which is an immersion of a torus with
constant mean curvature in $\mathbb{R}^{3}$. In the same vein it is well known
that a compact simply connected surface with constant Gaussian curvature must
also be a round sphere.

Surfaces of revolution that have a constant, non-zero skew curvature
constitute a family of surfaces parametrized by a single parameter. As already
alluded to above, there are no compact regular surfaces of revolution with
constant non-zero skew curvature. However, in the non compact case regular
surfaces of revolution do appear, see [20].

For a brief illustration of this point, consider a surface of revolution with
the following parametrization for $t\geq 0$ and generator function $g(t)$:

$r(t,u)=(t\cdot\cos(u),t\cdot\sin(u),g(t))\quad.$ (10)

The absolute value of the difference between the two principal curvatures,
i.e. the skew curvature of the surface, is then:

$S(t)=\frac{-h^{3}(t)+t\cdot
h^{\prime}(t)-h(t)}{2\left(1+h^{2}(t)\right)^{3/2}}\quad,$ (11)

where $h(t)=g^{\prime}(t)$. The solution $h(t)=0$ to $S(t)=0$ gives horizontal
planes. The non-zero solutions to $-h^{3}(t)+t\cdot h^{\prime}(t)-h(t)=0$ are
the spheres:

$h(t)=\pm\frac{t}{\sqrt{\alpha-t^{2}}}\quad,$ (12)

where $\alpha$ is an integration constant.

For constant non-zero skew curvature we only need to consider $S(t)=1$ since
all other cases appear by scaling. Such surfaces (choosing $h(t)\geq 0$, which
corresponds to the red curves in the Figures below) satisfy

$\ln(t)+\frac{h(t)}{t\cdot\sqrt{1+h^{2}(t)}}=k\quad,$ (13)

where $k$ is a constant of integration. We solve for $h(t)$:

$h(t)=\frac{2(k-\ln(t))\cdot t}{\sqrt{1+(8k-4\ln(t))\cdot
t^{2}\cdot\ln(t)-4k^{2}\cdot t^{2}}}\quad.$ (14)

And thence, for each allowed $k$ we obtain $g(t)$ by (numerical) integration
of $h(t)$. For each value of $k$ there are positive $t$-intervals for which
the square root in the denominator is well defined in (14). There are at most
3 roots of the denominator, i.e. there are at most two $t$-intervals to
consider. There is one such interval (i.e one constant skew curvature surface)
for $k<k_{0}$ and two intervals (i.e. two constant skew curvature surfaces)
for $k>k_{0}$ where $k_{0}$ is the bifurcation value corresponding to the
vertical cylinder of radius $1/2$:

$k_{0}=1-\ln(2)\approx 0.307\quad.$ (15)

We integrate (numerically) the functions $h(t)$ for each value of $k$ and
obtain the corresponding generator functions $g(t)$ and thus the corresponding
surfaces of revolution. See figures 1, 2, 3, where the results are shown for
three values of $k>k_{0}$. The corresponding ’outer’ solutions – which do not
intersect the axis of revolution – are the most interesting here, since they
can be stacked periodically to give smooth non-compact surfaces with constant
skew curvature.

Figure 1:  Profile curves for surfaces of revolution with constant skew
curvature $1$ and $k=0.33>k_{0}$. The curves are obtained by periodically
stacking the (red) base curve solutions (with $h(t)\geq 0$) and their
symmetric (blue) companions (with $h(t)\leq 0$). Only the outer-most stacking
is smooth. The corresponding outer-most surface is also shown with two periods
only. That surface is compact in $\mathbb{R}^{3}/\mathbb{Z}$ as discussed
below.

Figure 2: Similar to figure 1, except here with $k=0.5>k_{0}$ .

Figure 3: Similar to figure 1, except here with $k=0.7>k_{0}$ .

##  VI Upper bounds for the spectral gap via geometric compactification

The techniques of this paper aim to derive upper bounds for the first gap in
the spectrum for compact surfaces immersed in the Euclidean space. However,
for non-compact surfaces with certain symmetries, some – or all – of our
results can still be adapted.

One specific scenario where our technique can be applied is for non-compact
surfaces of revolution with a profile curve that has both a translational
symmetry and a reflection symmetry. This is precisely the case for the profile
curves shown in figures 1, 2, and 3 and thence the associated surfaces of
revolution.

These surfaces are actually immersed surfaces in $\mathbb{R}^{3}$, obtained as
the image of an immersion
$\varphi:M=\mathbb{R}\times\mathbb{S}^{1}\to\mathbb{R}^{3}$. The surface is
said to have a translation symmetry if there exists a vertical vector
$\vec{A}\in\mathbb{R}^{3}$ such that the surface is invariant under
translation by $\vec{A}$. Additionally, it is said to have a reflection
symmetry around the $z$ axis if the surface is invariant under reflection
about this axis.

We say that $\Sigma=\varphi(M)$ admits a translation symmetry because there is
a vertical vector $\vec{A}\in\mathbb{R}^{3}$ such that the surface sigma is
invariant under the translation of $\vec{A}=(0,0,A)$,

$T_{\vec{A}}(\Sigma):=\\{x+\vec{A}\,:\,x\in\Sigma\\}=\Sigma.$

Moreover, we say that $\Sigma$ admits a reflection $R_{z}$ around the $z$ axis
because,

$R_{z}(\Sigma):=\\{(x_{1},x_{2},-x_{3}):\,(x_{1},x_{2},x_{3})\in\Sigma\\}=\Sigma.$

When the profile curve is arc-length parametrized, $M$ can be equipped with
global coordinates $s$ and $\theta$, where $s$ is the arc-length parameter and
$\theta$ is the corresponding rotation angle. The Hamiltonian operator
$\widetilde{H}=\frac{\hbar^{2}}{2m}\left[-\Delta-\frac{1}{4}\left(k_{1}-k_{2}\right)^{2}\right]$
is symmetric with respect to a translation $\tau_{a}(s,\theta)=(s+a,\theta)$
and a reflection parity change $P(s,\theta)=(-s,\theta)$. Therefore, there
exists unitary operators $\widetilde{\tau}_{a}$ and $\widetilde{P}$ such that

$[\widetilde{H},\widetilde{\tau}_{a}]=0,\quad[\widetilde{H},\widetilde{P}]=0.$

However, it is not possible to obtain a basis of simultaneous eigenstates for
$\widetilde{H}$, $\widetilde{\tau}_{a}$, and $\widetilde{P}$, because
$\widetilde{\tau}_{a}$ and $\widetilde{P}$ do not commute.

We are interested in the composition $T_{a}=\tau_{a}\circ P$, which has the
property that $T_{a}\circ T_{a}=I$. This allows us to obtain a simultaneous
basis of eigenstates for $\widetilde{H}$ and $\widetilde{T}_{a}$, with
eigenvalues of $\widetilde{T}_{a}$ being $\pm 1$. Therefore, any eigenfunction
of $\widetilde{H}$ can be decomposed into eigenfunctions of
$\widetilde{T}_{a}$ with eigenvalues $\pm 1$.

$\Psi_{E}(s,\theta)=A\Psi^{+}_{E}(s,\theta)+B\Psi^{-}_{E}(s,\theta),\quad
A,B\in\mathbb{C},$

where

$\widetilde{T}_{a}^{2}\Psi^{+}_{E}(s,\theta)=\Psi^{+}_{E}(s,\theta),\quad\widetilde{T}_{a}^{2}\Psi^{-}_{E}(s,\theta)=-\Psi^{-}_{E}(s,\theta).$

Then, when we look for the wave functions

$\langle
s,\theta|\widetilde{T}_{a}|E_{i},\epsilon_{i}\rangle=\epsilon_{i}\Psi_{E_{i},\epsilon_{i}}(s,\theta),$

But since $\widetilde{T}_{a}|s,\theta\rangle=|-s+a,\theta\rangle$ we can
conclude that

$\Psi_{E_{i},\epsilon_{i}}(-s+a,\theta)=\epsilon_{i}\Psi_{E_{i},\epsilon_{i}}(s,\theta).$
(16)

Taking $s=0$ we obtain

$\Psi_{E_{i},\epsilon_{i}}(a,\theta)=\epsilon_{i}\Psi_{E_{i},\epsilon_{i}}(0,\theta).$

Hence, in the case $\epsilon_{i}=1$,

$\Psi_{E_{i},1}(a,\theta)=\Psi_{E_{i},1}(0,\theta).$ (17)

In summary, the spectrum $\sigma(\Sigma)$ of $\Sigma$ can be divided into two
parts: one, denoted as $\sigma^{+}(\Sigma)$, is related to wave functions of
positive parity with respect to $\widetilde{T}_{a}$, and the other, denoted as
$\sigma^{-}(\Sigma)$, is related to wave functions of negative parity with
respect to $\widetilde{T}_{a}$. If $\Sigma$ has positive parity, the wave
function must satisfy the boundary value problem (17).

From a geometric perspective, the spaces $\mathbb{R}^{3}$ and $M$ can be
divided into equivalence classes using the relation of equivalence
$(x,y,z)\sim(x^{\prime},y^{\prime},z^{\prime})$ in $\mathbb{R}^{3}$ if and
only if there exist $k\in\mathbb{Z}$ such that

$(x,y,z)-(x^{\prime},y^{\prime},z^{\prime})=k\vec{A},$

and $(s,\theta)\sim(s^{\prime},\theta^{\prime})$ in $M$ if and only if, there
exist $k\in\mathbb{Z}$ such that

$T^{k}_{a}(s,\theta)=(s^{\prime},\theta^{\prime}).$

The quotient spaces will be denoted as $\mathbb{R}^{3}/\mathbb{Z}$ and
$M/\mathbb{Z}$, respectively. Since the translations act by isometries we can
use the following immersion $\beta:M/\mathbb{Z}\to\mathbb{R}^{3}/\mathbb{Z}$
by using the following commutative diagram:

${M}$${\mathbb{R}^{3}}$${M/\mathbb{Z}}$${\mathbb{R}^{3}/\mathbb{Z}}$$\scriptstyle{\varphi}$$\scriptstyle{\pi}$$\scriptstyle{\pi}$$\scriptstyle{\beta}$

where $\pi$ denotes the canonical map to the classes of equivalence,
$x\mapsto\pi(x):=[x]$. Because the spaces $M$ and $M/\mathbb{Z}$ and the
spaces $\mathbb{R}^{3}$ and $\mathbb{R}^{3}/\mathbb{Z}$ are locally isometric,
since $\varphi$ is an isometric immersion, $\beta$ is an isometric immersion
as well with the same values of $H$ and $K$.

The immersion $\beta:M/\mathbb{Z}\to\mathbb{R}^{3}/\mathbb{Z}$ can be used to
obtain the first gap in the spectrum of $\Sigma$ with positive parity. If
$M/\mathbb{Z}$ is smooth enough, this gap is bounded by the following:

$E_{1}-E_{0}\leq\frac{\hbar^{2}}{2m}\left(\frac{2\mathcal{W}(\beta(M/\mathbb{Z}))}{{\rm
A}(\beta(M/\mathbb{Z}))}+\frac{1}{4}\left(\max_{\beta(M/\mathbb{Z})}S^{2}-\overline{S^{2}}\right)\right)$

for $E_{1}$ and $E_{0}$ in $\sigma^{+}(\Sigma)$. If the surface $\Sigma$ is
constructed as the translation of a fundamental domain $L$ and has a
reflection symmetry, the first gap in the spectrum of positive parity is
bounded by:

$E_{1}-E_{0}\leq\frac{\hbar^{2}}{2m}\left(\frac{2\mathcal{W}(L)}{{\rm
A}(L)}+\frac{1}{4}\left(\max_{L}S^{2}-\overline{S^{2}}\right)\right).$

In the case of constant skew curvature, this becomes:

$E_{1}-E_{0}\leq\frac{\hbar^{2}}{2m}\frac{2\mathcal{W}(L)}{{\rm A}(L)}.$

Observe that this upper bounds proportional to the square mean curvature are
similar to the bounds obtained in [7], where truncated surfaces of revolution
are used with general boundary conditions (Dirichlet, Neumann, or Robin).

In the case shown by figures 1, 2 and 3, since $M/\mathbb{Z}$ is a torus with
constant skew curvature $S=1$, by using the Gauss-Bonnet theorem,

${\rm
A}(M/\mathbb{Z})=\int_{M/\mathbb{Z}}S^{2}dA=\int_{M/\mathbb{Z}}H^{2}dA=\mathcal{W}(M/\mathbb{Z}),$

and therefore

$E_{1}-E_{0}\leq\frac{\hbar^{2}}{m}.$

##  VII Conclusions

The main objective of this study was to investigate which geometric or
topological properties play a significant role in describing a particle
confined to a compact surface. With a compact surface possessing a bounded
effective and geometrically induced potential, the energy levels are discrete.
Determining the energy spectrum of a particle confined to a surface can be a
challenging task, even numerically. However, this study demonstrates that the
difference between the energy of the first excited state and the ground level
can be upper-bounded by two terms. One of these terms is related to the total
area and the Willmore energy, or some topological constant that depends only
on the genus of the surface. The other term measures how far the surface
deviates from being a surface of constant skew curvature.

This leads naturally to a study of the properties of surfaces with constant
skew curvature, specifically the surfaces of revolution of constant skew
curvature. They constitute a one-parameter family with a bifurcation point
between the compact and non-compact cases. Although the results in this report
are stated for the compact case, if a non-compact surface possesses enough
symmetries, the results can still be used via the shown compactification of
the surface to obtain bounds for a specific part of the original surface’s
spectrum.

## Acknowledgments

This Work was partially supported by the Research grant PID2020-115930GA-100
funded by MCIN/ AEI /10.13039/501100011033, and AICO/2021/252. The authors
would like to thank Luiz C. B. da Silva for comments and suggestions that have
improved the presentation of the present work. The authors would also like to
acknowledge the hospitality of the departments of mathematics at DTU and at
Universitat Jaume I, respectively.

## References

  * Zheng and Ke [2010] M. Zheng and C. Ke, Elastic deformation of carbon-nanotube nanorings, Small 6, 1647 – 1655 (2010).
  * Li  _et al._ [2020] F. Li, H. Li, M. Li, and Z. Kang, Mechanics of folding of nanorings, Mechanics of Materials 148, 10.1016/j.mechmat.2020.103493 (2020).
  * Ortix and Van Den Brink [2010] C. Ortix and J. Van Den Brink, Effect of curvature on the electronic structure and bound-state formation in rolled-up nanotubes, Physical Review B - Condensed Matter and Materials Physics 81, 10.1103/PhysRevB.81.165419 (2010).
  * Sadykov and Yudina [2020] N. R. Sadykov and N. V. Yudina, Geometrical potential in giant fullerenes, Technical Physics 65, 369 – 376 (2020).
  * Cheng  _et al._ [2021] R. Cheng, Y.-L. Wang, H. Zhao, C.-Z. Ye, G.-H. Liang, and H.-S. Zong, Bound states and energy shifts resulting from corrugations, Results in Physics 22, 10.1016/j.rinp.2021.103974 (2021).
  * Gravesen  _et al._ [2005a] J. Gravesen, M. Willatzen, and L. C. Lew Yan Voon, Quantum-mechanical particle confined to surfaces of revolution -truncated cone and elliptic torus case studies, Physica Scripta 72, 105 – 111 (2005a).
  * Willatzen  _et al._ [2010] M. Willatzen, J. Gravesen, and L. C. Lew Yan Voon, Analytic theory of curvature effects for wave problems with general boundary conditions, Physical Review A \- Atomic, Molecular, and Optical Physics 81 (2010).
  * Gravesen  _et al._ [2005b] J. Gravesen, M. Willatzen, and L. C. Lew Yan Voon, Schrödinger problems for surfaces of revolution-the finite cylinder as a test example, Journal of Mathematical Physics 46, 10.1063/1.1829376 (2005b).
  * Gravesen and Willatzen [2020] J. Gravesen and M. Willatzen, Quantum eigenstates of curved and varying cross-sectional waveguides, Applied Sciences (Switzerland) 10, 1 – 19 (2020).
  * DeWitt [1957] B. S. DeWitt, Dynamical theory in curved spaces. i. a review of the classical and quantum action principles, Reviews of Modern Physics 29, 377 – 397 (1957).
  * Jensen and Koppe [1971] H. Jensen and H. Koppe, Quantum mechanics with constraints, Annals of Physics 63, 586 (1971).
  * da Costa [1981] R. C. T. da Costa, Quantum mechanics of a constrained particle, Phys. Rev. A (3) 23, 1982 (1981).
  * Jensen and Dandoloff [2011] B. Jensen and R. Dandoloff, Remarks on quantum mechanics on surfaces, Physics Letters, Section A: General, Atomic and Solid State Physics 375, 448 – 451 (2011).
  * Ortix and Van Den Brink [2011] C. Ortix and J. Van Den Brink, Absence of anomalous couplings in the quantum theory of constrained electrically charged particles, Physical Review B - Condensed Matter and Materials Physics 83, 10.1103/PhysRevB.83.113406 (2011).
  * Matsutani [1993] S. Matsutani, Quantum field theory on curved low-dimensional space embedded in three-dimensional space, Physical Review A 47, 686 – 689 (1993).
  * Burgess and Jensen [1993] M. Burgess and B. Jensen, Fermions near two-dimensional surfaces, Physical Review A 48, 1861 – 1868 (1993).
  * Jensen and Dandoloff [2009] B. Jensen and R. Dandoloff, Quantum mechanics of a constrained electrically charged particle in the presence of electric currents, Physical Review A - Atomic, Molecular, and Optical Physics 80, 10.1103/PhysRevA.80.052109 (2009).
  * Homma  _et al._ [1990] T. Homma, T. Inamoto, and T. Miyazaki, Schrödinger equation for the nonrelativistic particle constrained on a hypersurface in a curved space, Physical Review D 42, 2049 – 2056 (1990).
  * Ferrari and Cuoghi [2008] G. Ferrari and G. Cuoghi, Schrödinger equation for a particle on a curved surface in an electric and magnetic field, Physical Review Letters 100, 10.1103/PhysRevLett.100.230403 (2008).
  * López and Pámpano [2020] R. López and A. Pámpano, Classification of rotational surfaces with constant skew curvature in 3-space forms, Journal of Mathematical Analysis and Applications 489, 10.1016/j.jmaa.2020.124195 (2020).
  * da Silva [2021] L. C. B da Silva, Surfaces of revolution with prescribed mean and skew curvatures in lorentz-minkowski space, Tohoku Mathematical Journal 73, 317 – 339 (2021).
  * Toda [2017] M. D. Toda,  _Willmore energy and willmore conjecture_ (2017) p. 1 – 142.
  * El Soufi and Ilias [2000] A. El Soufi and S. Ilias, Second eigenvalue of Schrödinger operators and mean curvature, Comm. Math. Phys. 208, 761 (2000).
  * Chavel [1984] I. Chavel,  _Eigenvalues in Riemannian geometry_ , Pure and Applied Mathematics, Vol. 115 (Academic Press Inc., Orlando, FL, 1984) pp. xiv+362, including a chapter by Burton Randol, With an appendix by Jozef Dodziuk.
  * Grigor’yan  _et al._ [2016] A. Grigor’yan, N. Nadirashvili, and Y. Sire, A lower bound for the number of negative eigenvalues of Schrödinger operators, J. Differential Geom. , 395 (2016).
  * Yang and Yau [1980] P. C. Yang and S. T. Yau, Eigenvalues of the Laplacian of compact Riemann surfaces and minimal submanifolds, Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) 7, 55 (1980).
  * da Silva  _et al._ [2017] L. C. B da Silva, C. C Bastos, and F. G Ribeiro, Quantum mechanics of a constrained particle and the problem of prescribed geometry-induced potential, Annals of Physics 379, 13 (2017).
  * Toda and Pigazzini [2017] M. Toda and A. Pigazzini, A note on the class of surfaces with constant skew curvatures, Journal of Geometry and Symmetry in Physics 46, 51 – 58 (2017).