# Detecting Oceans on Exoplanets with Phase-Dependent Spectral Principal
Component Analysis

Dominick J. Ryan  Department of Astronomy and Planetary Science, Northern
Arizona University, Box 6010, Flagstaff, AZ 86011, USA  Habitability,
Atmospheres, and Biosignatures Laboratory, Northern Arizona University,
Flagstaff, AZ 86011, USA  Tyler D. Robinson Department of Astronomy and
Planetary Science, Northern Arizona University, Box 6010, Flagstaff, AZ 86011,
USA  Habitability, Atmospheres, and Biosignatures Laboratory, Northern Arizona
University, Flagstaff, AZ 86011, USA  NASA Astrobiology Institute’s Virtual
Planetary Laboratory, University of Washington, Box 351580, Seattle, WA 98195,
USA

###### Abstract

Stable surface liquid water is a key indicator of exoplanet habitability.
However, few approaches exist for directly detecting oceans on potentially
Earth-like exoplanets. In most cases, specular reflection of host starlight
from surface bodies of water — referred to as ocean glint — proves to be an
important aspect of liquids that can enable detection of habitable conditions.
Here, we propose that spectral principal component analysis (PCA) applied to
orbital phase-dependent observations of Earth-like exoplanets can provide a
straightforward means of detecting ocean glint and, thus, habitability. Using
high-fidelity, orbit-resolved spectral models of Earth, and for instrument
capabilities applicable to proposed exo-Earth direct imaging concept missions,
the extreme reddening effect of crescent-phase ocean glint is demonstrated as
the primary spectral component that explains phase-dependent variability for
orbital inclinations spanning 60–90°. At smaller orbital inclinations where
more-extreme crescent phases cannot be accessed, glint can still significantly
increase planetary brightness but reddening effects are less pronounced and,
as a result, glint is not plainly indicated by phase-dependent spectral PCA.
Using instrument models for future exoplanet direct imaging mission concepts,
we show that brightness enhancements due to glint could be detected across a
wide range of orbital inclinations with typical exposure times measured in
hours to weeks, depending on system distance and mission architecture. Thus,
brightness increases due to glint are potentially detectable for Earth-like
exoplanets for most system inclinations and phase-dependent spectral PCA could
indicate reddening due to glint for a subset of these inclinations.

Principal component analysis (1944) — Earth (planet) (439) — Habitable planets
(695) — Direct imaging (387)

††journal: PSJ

##  1 Introduction

The role of stable liquid water (i.e., habitability) has long been highlighted
as critical to the origin and evolution of life on Earth (e.g., Brack, 1993).
As a result, foundational works on the potential for habitable worlds outside
the Solar System have largely emphasized planets that may be able to maintain
surface oceans (Hart, 1978; Matsui & Abe, 1986; Abe & Matsui, 1988; Whitmire
et al., 1991; Kasting et al., 1993). Since these early works, myriad studies
have explored processes that impact exoplanet surface thermal conditions,
liquid water stability, and habitability (e.g., Wolf & Toon, 2010; von Paris
et al., 2010; Pierrehumbert & Gaidos, 2011; Goldblatt et al., 2013; Kopparapu
et al., 2013; Rugheimer et al., 2013; Wordsworth & Pierrehumbert, 2013;
Leconte et al., 2013; Shields et al., 2013; Ramirez & Kaltenegger, 2014; Way
et al., 2016; Meadows et al., 2018; Fauchez et al., 2018; Turbet et al.,
2018). For recent reviews related to exoplanet habitability, see Seager
(2013), Meadows & Barnes (2018), and Kopparapu et al. (2020).

Near- or far-future space telescopes that could potentially characterize
exoplanets for habitability indicators (such as NASA’s James Webb Space
Telescope as well as the Habitable Exoplanet Observatory [HabEx], Large
UltraViolet-Optical-InfraRed Surveyor [LUVOIR], Origins Space Telescope, or
Large Interferometer for Exoplanets [LIFE] concepts; Gardner et al., 2006;
Gaudi et al., 2018; Roberge & Moustakas, 2018; Battersby et al., 2018; Quanz
et al., 2018) will not spatially resolve their exoplanet targets. Thus, signs
of a potentially habitable surface environment must be gleaned from unresolved
observations, be these disk-integrated spectra from directly-imaged worlds or
transit and/or secondary eclipse data for transiting exoplanets. While
exoplanet spectra from any observational approach may indirectly indicate
habitable surface conditions, there is only one leading direct indicator of
surface liquid water — ocean glint.

Ocean glint is a mirror-like (i.e., specular) reflection of sunlight from a
liquid surface. Rooted in the Fresnel equations (Fresnel, 1821), glint is
expected to produce polarization and reflectance features in phase curves of
ocean-bearing exoplanets (McCullough, 2006). Critically, glint can cause
planetary reflectance to increase dramatically at crescent phases as the glint
spot becomes both more reflective as well as a larger portion of the
illuminated planetary disk (Williams & Gaidos, 2008). Observations and
realistic models of Earth’s phase-dependent brightness reveal that glint can
make our planet roughly twice as bright at crescent phases (as compared to a
non-glinting Earth), especially at red-visible and near-infrared wavelengths
(Robinson et al., 2010; Zugger et al., 2011; Robinson et al., 2014).
Polarization signatures due to glint from an Earth-like exoplanet would reach
a maximum when the planet is slightly less than half-illuminated, with peak
polarization fractions potentially approaching 20% (Stam, 2008; Zugger et al.,
2010). Importantly, polarization and/or near-infrared reflectance observations
could rule out a proposed false positive signature for glint where the high-
latitude (presumably more ice-covered and reflective) regions of a planet are
preferentially probed at crescent phases (Cowan et al., 2012), as discussed by
Robinson (2018).

Recently, Lustig-Yaeger et al. (2018) proposed the concept of “glint mapping”
for directly-imaged exoplanets, which takes advantage of a strong contribution
from specular reflection when an exo-Earth is viewed near crescent phase
(Cowan et al., 2009). Here, the “blinking” of the glint spot as the specular
point rotates across oceans and Lambert-reflecting continents can enable the
inference of longitudinally-resolved ocean maps from crescent-phase,
rotationally-resolved exo-Earth observations. Comparisons of reflectance maps
generated at quadrature phase versus crescent phase reveals the emergence of a
highly-reflective surface feature at the latter phase, strongly indicative of
surface oceans.

A challenge for glint mapping is that the approach requires reflected-light
phase curves of potential exo-Earths to be rotationally resolved at crescent
phase. Here, the planet is relatively quite faint and potentially near the
inner working angle of a high contrast imaging instrument. Such a requirement
then restricts the number of potential targets for glint mapping as well as
the scale of proposed future space telescopes that could apply the technique.
Thus, a glint detection method that does not require rotationally-resolved
lightcurves could complement the glint mapping technique and broaden the range
of exo-ocean detectability.

Robinson (2018) noted that, since the path that ocean glinted photons take
through the atmosphere is simple (direct trajectories with a single scattering
event at the surface), the glint spot has a very characteristic top-of-
atmosphere spectrum. More specifically, the glint spectrum is a solar spectrum
modulated by simple loses due to clearsky atmospheric gas absorption and
Rayleigh scattering opacity for straight-line trajectories towards and away
from the surface. Furthermore, as glint can be the dominant reflection source
at crescent phases for an exo-Earth, the simple glint spectrum should be
readily apparent in the disk-integrated spectrum of such a world near crescent
phases. At these crescent-phase geometries, long atmospheric pathlengths lead
to characteristic glint spectra that is heavily reddened.

The underlying idea that drives the work that follows is that spectral
principal component analysis applied to phase-resolved observations of an exo-
Earth could reveal crescent-phase glint reddening as a dominant spectral
contributor. Section 2 outlines the Earth spectral models, instrument
simulators, and numerical methods adopted in this work. Then, Section 3
demonstrates (1) a “proof of concept,” (2) an exploration of the impacts of
both orbital inclination (which impacts the range of accessible phase angles)
and number of target revisits on the detectability of a glint principal
component, and (3) results for glint detectability for future exoplanet direct
imaging mission concepts. Section 4 connects these new results to previous
studies, notes potential caveats, and highlights future avenues for
investigation. Finally, Section 5 outlines key conclusions.

##  2 Methods

The study at hand aims to understand how principal component analysis (PCA)
could be used to detect a glint component in phase-dependent spectra of an
ocean-bearing world. As is described in this section, we adopt high-fidelity,
phase- and time-dependent spectra of Earth as a stand-in for a hypothetical
exo-Earth. Observational quality (in terms of spectral resolving power and
accessible range of planetary phases) are degraded to be roughly consistent
with future exo-Earth direct imaging missions, and standard principal
component analysis techniques are applied to resulting datasets.

###  2.1 Spectral Models

Phase-dependent spectra of Earth come from the Virtual Planetary Laboratory
3-D spectral Earth model (Tinetti et al., 2006a, b; Robinson et al., 2011).
Here, spectral flux density at a given distance from Earth is computed as the
integral of the projected area-weighted intensity in the direction of the
observer,

$F_{\lambda,{\rm E}}\\!\left(\bm{\hat{o}},\bm{\hat{s}}\right)=\frac{R_{\rm
E}^{2}}{d^{2}}\int_{2\pi}\left(\bm{\hat{n}}\cdot\bm{\hat{o}}\right)I\\!\left(\bm{\hat{n}},\bm{\hat{o}},\bm{\hat{s}}\right)\mathrm{d}\omega\,$
(1)

where $F_{\lambda,{\rm E}}$ is Earth’s disk-integrated spectral flux density,
$R_{\rm E}$ is the radius of Earth, $d$ is the Earth-observer distance,
$I\\!\left(\bm{\hat{n}},\bm{\hat{o}},\bm{\hat{s}}\right)$ is the specific
intensity in the direction of the observer as a function of location on the
disk, $\bm{\hat{o}}$ is a unit vector in the direction of the observer,
$\bm{\hat{s}}$ is a unit vector in the direction of the Sun, $\bm{\hat{n}}$ is
a surface normal vector corresponding to $\mathrm{d}\omega$, and the solid
angle integral is over the observable hemisphere. The integral over the
planetary disk is performed numerically by pixelating the planet according to
the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) model (Górski
et al., 2005)111healpix.sourceforge.io. High-resolution, position-dependent
intensity spectra are computed using the one-dimensional, plane-parallel,
multiple-scattering Spectral Mapping Atmospheric Radiative Transfer (SMART)
model (developed by D. Crisp; Meadows & Crisp, 1996). Results from the VPL 3-D
Spectral Earth model have been validated against a wide range of satellite,
interplanetary spacecraft, and Earthshine observations of our planet (Robinson
et al., 2010, 2011, 2014; Schwieterman et al., 2015).

Analyses presented here are based on a suite of publicly-available VPL 3-D
Earth model spectra that capture Earth through an entire orbit from multiple
vantages (Robinson et al., 2010). In these datasets, high resolution spectra
are generated every four hours for the entire calendar year of 2008 (taken to
be concurrent with observations of Earth from NASA’s EPOXI mission; Livengood
et al., 2011; Robinson et al., 2011). The three-dimensional atmospheric state,
including water vapor mixing ratios and cloud coverage, are interpolated from
daily observations derived from Earth-observing satellite observations (Hall
et al., 1995; Riggs et al., 1999; Aumann et al., 2003; Beer et al., 2001;
Waters et al., 2006; Wolfe, 2006). Distinct spectral datasets were generated
for observing geometries that encompass orbits viewed edge-on (i.e.,
inclination, $i$, of 90°) through face-on (i.e., inclination of 0°). As
planets with orbits viewed at or near face-on do not achieve the crescent
phases required for significant glint contributions, only the datasets with
inclinations of 45°, 60°, 75°, and 90° are investigated here.

Phase-dependent spectral flux densities are converted into spectra of apparent
albedo (defined by normalizing the planetary flux to that from a perfectly
reflecting Lambert sphere observed at the same phase angle) with,

$A_{\rm app}=\frac{3}{2}\frac{F_{\lambda,{\rm E}}}{F_{\lambda,{\rm
S}}}\frac{\pi}{\sin\alpha+(\pi-\alpha)\cos\alpha}\ ,$ (2)

where $F_{\lambda,{\rm E}}$ is measured at the top-of-atmosphere (i.e., with
$d=R_{\rm E}$), $F_{\lambda,{\rm S}}$ is the solar spectral flux density at
normal incidence at the top of Earth’s atmosphere, and $\alpha$ is the star-
planet-observer (i.e., phase) angle. (Note that
$\cos\alpha=\bm{\hat{s}}\cdot\bm{\hat{o}}$.) As mentioned in Robinson &
Reinhard (2020), the factor of $3/2$ is due to the conversion between
geometric albedo and spherical albedo for a Lambertian reflector. Apparent
albedo spectra are useful when studying phase-dependent planetary spectra as
they help to identify phases where the planet preferentially scatters
radiation. In the case of Earth, apparent albedo values tend to increase
dramatically towards crescent phases, indicating strong forward scattering
from the planet (Qiu et al., 2003; Pallé et al., 2003). Note, however, that
the results below are tied most fundamentally to phase-dependent color
variations and raw brightness enhancements, so that the findings are
independent of the (wavelength-independent) conversion to apparent albedo.

Of course, high contrast imaging of exoplanets in reflected light does not
record the top-of-atmosphere planetary flux, but is instead sensitive to the
planetary flux observed at Earth or the planet-to-star flux ratio,

$\frac{F_{\rm p}}{F_{\rm s}}=A_{\rm g}\phi(\alpha)\frac{R_{\rm p}^{2}}{r^{2}}\
,$ (3)

where $A_{\rm g}$ is the wavelength-dependent planetary geometric (i.e., full-
phase) albedo, $\phi$ is the wavelength- and phase angle-dependent planetary
phase function, $R_{\rm p}$ is the planetary radius, and $r$ is the star-
planet distance. Given that the planet-to-star flux ratio is measured at the
star-Earth distance, and generalizing the Earth and Sun spectra in Equation 2,
the apparent albedo can be written in terms of the planet-to-star flux ratio
as,

$A_{\rm app}=\frac{3}{2}\frac{r^{2}}{R_{\rm p}^{2}}\frac{F_{\rm p}}{F_{\rm
s}}\frac{\pi}{\sin\alpha+(\pi-\alpha)\cos\alpha}\ .$ (4)

Thus, the apparent albedo can be inferred from exoplanet observations. Here,
though, the planetary radius may not be known, so that either the apparent
albedo scaled by $R_{\rm p}^{-2}$ or the planet-to-star flux ratio relative to
the Lambert phase function could be more immediately useful quantities for
understanding phase-dependent scattering from the planet.

###  2.2 Simulated Observations

The planetary phases that are observationally accessible for an exoplanet
depend both on the planetary orbital configuration and the details of the high
contrast imaging instrument (especially its inner working angle). For the
former, the phase angle at any point in the planetary orbit is given by,

$\cos\alpha=\sin\left(\theta+\omega_{\rm p}\right)\sin i\ ,$ (5)

where $\theta$ is the true anomaly, $\omega_{\rm p}$ is the argument of
periastron, and, as before, $i$ is the orbital inclination (Madhusudhan &
Burrows, 2012). Assuming an Earth-like orbit (effectively circular, 1 au from
a Sun-like host) for a hypothetical exo-Earth, and given the preceding
expression, an orbital inclination, a distance to the exoplanet system, and an
inner working angle, the analyses below assume that observations are equally
spaced (in time or, equivalently, true anomaly) from where the planet emerges
from the inner working angle (at gibbous phase) to where the planet disappears
into the inner working angle (at crescent phase).

For simplicity in this concept study, impacts of varying the instrument inner
working angle and planetary system distance are not explored. Instead, we
adopt an inner working angle of 100 mas, which is consistent with the near-
infrared performances of the Habitable Exoplanet Observatory (HabEx; proposed
to deliver an inner working angle of 110 mas at 1.8 $\upmu$m with its
starshade) and Large UltraViolet-Optical-InfraRed Explorer (LUVOIR; proposed
to deliver an inner working angle of 90 mas at 1.8 $\upmu$m with its 15-meter
architecture) concepts (Gaudi et al., 2018; Roberge & Moustakas, 2018). Red
and near-infrared wavelengths have already been identified as key ranges for
glint detection (Robinson et al., 2010; Zugger et al., 2011). Also for
simplicity, we fix the distance to the hypothetical exo-Earth system studied
below at 5 pc. While this is a somewhat optimistic assumption, the discussion
section highlights the circumstances where this distance could be increased.
Putting all of this together, and as an example, five equally-spaced
observations of this hypothetical exo-Earth in an edge-on orbit would sample
phase angles of 30°, 60°, 90°, 120°, and 150° (i.e., well-spanning gibbous to
crescent phases). Figure 1 demonstrates spectrally-degraded apparent albedo
spectra from the VPL Earth models across this range of phase angles, where
reddening and increasing apparent albedo towards crescent phase is indicative
of glint.

Figure 1:  Apparent albedo spectra from the VPL Earth model across a range of
phase angles and degraded to a resolving power of 140.

Simulated observations are created by degrading the wavelength coverage and
resolving power of high-resolution VPL Earth models to values consistent with
both HabEx and LUVOIR. Specifically, spectra are degraded to a resolving power
($\lambda$/$\Delta\lambda$) of 140 across the 0.4–1.8 $\upmu$m range. While
HabEx and LUVOIR may adopt coarser resolving powers for exo-Earth observations
in the near-infared (e.g., 40–70), glint features occur in the spectral
continuum and, thus, likely do not require even low-resolution spectroscopy
for identification. Finally, observations are taken as a 24-hour average of
VPL Earth model spectra at a given phase angle, thereby smoothing over any
rotational effects.

###  2.3 PCA Tool

Principal component analysis was performed using a straightforward, publicly
available PCA tool distributed as part of the sklearn.decomposition software
suite. Here, singular value decomposition is applied, which uses linear
algebraic techniques to decompose a dataset into its eigencomponents. The
number of eigencomponents derived for a spectral dataset is, at most, equal to
the number of independent observations in the set.

##  3 Results

The orbital inclination of a directly-imaged exo-Earth directly determines the
range of planetary phase angles that could be accessible to observation while
mission architecture and timing constraints strongly influence the number of
phase observations that could feasibly be obtained. As a proof of concept,
results below begin with an ideal scenario, where the exo-Earth orbit is
viewed edge on (meaning a wide range of phase angles could be sampled,
depending on the instrument inner working angle) and a large number of
revisits have occurred. Later subsections then explore sensitivity to viewing
geometry and revisits. Finally, results related to observability of glint
brightness enhancements (which are distinct from crescent-phase glint
reddening) with future exoplanet direct imaging missions are given.

###  3.1 Proof of Concept

Figure 2 demonstrates the functionality of glint detection via phase-resolved
spectral PCA. The planetary orbit was assumed to be edge-on ($i=90$°) and ten
observations were taken to span gibbous phase at the detector IWA through to
crescent phase at the IWA. Thus, this setup represents a near “best case”
scenario. To best demonstrate this “proof of concept,” the experiment was
repeated for an identical phase-dependent Earth dataset where ocean glint was
removed. Finally, a scaled spectrum of crescent-phase, cloud-free Earth
uniformly covered with ocean is shown to help indicate the spectral shape of
the raw glint signature. In the realistic Earth case, a glint-like spectrum is
identified as the main principle component, and this glint-like component is
not present for the non-glinting Earth case. Strikingly, the glint-like
component explains nearly 80% of the phase-dependent variance.

Figure 2: Results of PCA applied to phase-dependent Earth spectra. Orbital
inclination is edge-on ($i=90$°) and ten observations span gibbous to crescent
phase. Realistic Earth results (left) can be compared to a glint-free Earth
(right). The first three principal components are indicated, and a crescent-
phase, cloud-free ocean-covered Earth spectrum is used to indicate the stand-
alone glint signature (gray). Proportion of variance explained by components
one, two, and three are 0.79, 0.20, and 0.006, respectively.

###  3.2 Orbital Inclination Effects

Equation 5 indicates that the range of phase angles that could be accessed for
a given planet is reduced as the orbital inclination decreases away from a
edge-on geometry. Figures 3, 4, and 5 explore the impacts of apparent orbital
orientation, and correspond to inclinations of 75°, 60°, and 45°,
respectively. As in the previous section, ten phase-dependent observations are
assumed. For the adopted inclinations of 75°, 60°, and 45°, the maximum
possible phase angles are 165°, 150°, and 135°, respectively. However, except
for the $i=45$° case, all observing scenarios have their most extreme
accessible crescent phase limited by the detector IWA, so that these cases
sample phase angles between 30–150°. The cases with inclinations of either 75°
or 60° detect a glint contribution as a primary component, as is evidenced by
the disappearance of a heavily-reddened principal component in the simulations
that have glint removed. However, the $i=45$° scenario does not detect a
strong glint contribution to phase-dependent color variations, as is indicated
by the PCA returning roughly the same principal components when applied to
identical phase-dependent models where glint has been removed.

Figure 3: Same as Figure 2 except for an orbital inclination of $75$°. For the
case that includes glint (left), proportion of variance explained by
components one, two, and three are 0.70, 0.29, and 0.007, respectively.

Figure 4: Same as Figure 2 except for an orbital inclination of $60$°. For the
case that includes glint (left), proportion of variance explained by
components one, two, and three are 0.72, 0.28, and 0.008, respectively. Figure
5: Same as Figure 2 except for an inclination of 45° as applied to models that
include glint (solid) as well as identical models where glint is removed
(dashed). For the case that includes glint, proportion of variance explained
by components one, two, and three are 0.85, 0.15, and 0.015, respectively.

###  3.3 Revisit Number Effects

Revisit observations can consume substantial amounts of observatory time,
especially for space-based telescopes using an external occulter (Gaudi et
al., 2020). To better understand revisit requirements for glint detection,
Figures 6 and 7 perform phase-resolved spectral PCA when only two observations
are considered — one at gibbous phase and one at crescent phase — for orbital
inclinations spanning 90–60°. In these circumstances, only two spectral
principal components can be extracted. In all inclination cases, a glint-like
spectrum is identified as the first principal component.

Figure 6: Results of PCA applied to phase-dependent Earth spectra when only
two phase-resolved observations are adopted — one at gibbous phase and one at
crescent phase. Orbital inclination is edge-on ($i=90$°). Component spectra
(colors) can be compared to the stand-alone glint signature (gray). Proportion
of variance explained by components one and two are 0.84 and 0.16,
respectively.

Figure 7: Same as Figure 6 except for an inclination of 75° (left) and 60°
(right). For the 75° case, proportion of variance explained by components one
and two are 0.74 and 0.26, respectively. For the 60° case, proportion of
variance explained by components one and two are 0.64 and 0.36, respectively.

###  3.4 Detection Considerations

The results above indicate that phase-dependent spectral PCA could reveal
reddening due to glint, especially for orbital inclinations larger than 60°.
Uncovering such a reddened component is distinct from detecting a brightness
enhancement due to glint, tough. In practice, detection of glint will likely
be limited by exposure times required to reach a requisite signal-to-noise
ratio (SNR) at crescent phases. At increasingly larger phase angles,
decreasing illumination works against increasing apparent albedo and a
stronger glint signature. For less-inclined orbits, the limited range of
accessible planetary phase angles (through Equation 5) implies that the
planetary brightness can remain relatively large while glint contributions are
more limited. Figure 8 shows modeled apparent albedo spectra of Earth over a
range of crescent phases for cases that do and do not include glint.
Additionally, this figure shows the fractional increase in apparent albedo due
to glint over a broader range of phase angles, thereby quantifying the phase-
dependent glint signal. Note that the maximum phase angle achieved for orbits
with inclinations of 30°, 45°, and 60° are 120°, 135°, and 150°, respectively.

Figure 8: Left: Apparent albedo spectra of Earth with (solid) and without
(dashed) glint at select crescent phases (colors). Right: Wavelength-dependent
fractional brightness increase due to glint over a range of crescent phase
angles.

The brightness increases due to glint in Figure 8 can be used to inform
calculations of required exposure times for glint detection when considering
future direct imaging missions. (This is distinct from detecting crescent
phase glint reddening, with the detection of a brightness enhancement due to
glint being more analogous to a scenario where spectral retrieval could be
used to fit for a specular component to surface reflectance.) For example, at
a phase angle of 135°, a glinting Earth is roughly 50% brighter than a glint-
free Earth, averaged over the I, Y, and J photometric bands. Thus, to achieve
a confident (SNR of 5) glint detection at this phase angle, observations would
need to achieve a spectral SNR of at least 10 ($=5/0.5$) in the near-infrared.
Figure 9 then shows the wavelength-dependent exposure time required for both
HabEx (with its starshade) and the LUVOIR “A” architecture (with a 15-meter
primary and coronagraph instrument) to reach an SNR of 10 for an Earth-twin at
5 pc and a phase angle of 135°. Exozodical light at the level of 4.5
“exozodis” is assumed (Mennesson et al., 2019), and noise calculations were
performed using the instrument modeling techniques described in Robinson et
al. (2016).

Figure 9: Exposure time required to reach a spectrally-resolved SNR of 10 for
both the HabEx (brown) and LUVOIR-A (green) mission concepts for a crescent-
phase ($\alpha=$135°) Earth-twin exoplanet orbiting a Sun-like host at 5 pc.

Figures 8 and 9 indicate that Y band (at 1.02 $\upmu$m) is a leading spectral
region for glint detection, as this band demonstrates a strong glint
signature, overall lower integration times (due, in part, to a decrease in
spectral resolution in the near-infrared for both concepts), and reasonable
inner working angle performance. Highlighting this band, Figure 10 shows the
exposure time required to achieve a spectral SNR large enough to detect the
crescent-phase contribution due to glint at a SNR of 5 for both the HabEx and
LUVOIR-A concepts at different distances and for different system
inclinations. At smaller distances, the maximum achievable phase angle is
limited only by inclination and requisite exposure times increase with
distance primarily due to planet count rates scaling with the inverse square
of distance. Eventually, though, the telescope inner working angle begins to
limit the range of accessible phase angles. This introduces a discontinuity in
the curves, and causes the requisite exposure times to increase more strongly
with distance due to the combination of limited access to crescent phases
(where the glint signature is strongest) and the aforementioned distance-
squared effect.

Figure 10: Exposure time required to detect the Y band crescent-phase glint
signature at a SNR of 5 as a function of distance to the exo-Earth twin and
assuming a Sun-like host. Results for both the HabEx (brown) and LUVOIR-A
(green) concepts are shown, and curves for different orbital inclinations are
indicated by line style.

##  4 Discussion

The Fresnel equations and in situ observations on Earth indicate that ocean
reflectivity will increase by 1–2 orders of magnitude from normal incidence
(applicable at full phase) to glancing scattering angles (applicable at
crescent phases). While aerosols also typically demonstrate single-scattering
forward enhancements, realistic, multiple-scattering clouds do not maintain
such a strong asymmetry (Thomas & Stamnes, 1999). Thus, it may not be
surprising that spectral PCA applied to simulated phase-dependent Earth
observations plainly identifies ocean glint as the primary variability
component, especially in an edge-on orbital geometry (Figure 2).

In the “proof of concept” study (and all subsequent glint detection
demonstrations), the glint component does not identically resemble the raw
glint signature obtained from a cloud-free, ocean-covered Earth model. While
the resemblance is generally strong at red and near-infrared wavelengths,
similarity decreases at blue wavelengths. Here, the first principal component
is likely also capturing phase-dependent variations due to Rayleigh
scattering, which has weak forward (and backward) scattering lobes in its
phase function. While not an emphasis of this work, the second and third
spectral principal components identified in the realistic Earth cases for
Figures 2–4 can likely be attributed to phase-dependent variations in cloud
reflectivity (for the second principal component, and noting that water clouds
tend to increase in absorptivity towards longer near-infrared wavelengths) and
Lambertian reflection from land surfaces (for the third principal component).

The inclination study presented in Section 3.2 points to (at least) two key
considerations for any future efforts to detect ocean glint from Earth-like
exoplanets. First, sensitivity to the windows between water absorption bands
at red and near-infrared wavelengths (again, especially in the I, Y, and J
photometric bands) is critical to the spectral detection of glint. Second,
access to larger phase angles ensures observations at sufficiently extreme
crescent phases for the most straightforward detection of reddening due to
glint. Thus, smaller telescope inner working angles can enable the types of
low-separation observations needed for strong glint detection at crescent
phases. Promisingly, PCA results show that glint is the dominant component
explaining phase-dependent variability at orbital inclinations greater than
60°, which would encapsulate roughly half of all detected potential exo-Earths
in a blind direct imaging survey. For exo-Earths with inclinations between
60°–90°, crescent phase reddening effects from glint are so strong that only
two phase-resolved observations are required to reveal the glint component
(Figures 6 and 7).

Inclination study results also indicate the crescent phases at which glint
reddening becomes too weak to impact PCA results. A glint-reddened component
is strongly inferred for studies at inclinations at 90°–60°. For these
inclinations, the component indicative of glint explains more than 70% of the
phase-dependent color variance, implying even low-SNR observations could
reveal the glint signature. By contrast, PCA applied to phase-resolved spectra
of Earth at $i=45$° uncovers similar principal components even if glint is
removed from the underlying models (Figure 5), implying a glint-reddened
component is not inferred. Here, the largest accessible phase angle is 135°,
and the corresponding spectrum in Figure 8 shows that the continuum is not
particularly reddened at this phase angle. (Although the planet is overall
brighter due to glint at these continuum wavelengths.) Thus, phase-dependent
color variations due to glint over this more narrow range of phase angles are
grey and similar to those of clouds.

Detailed studies of the detectability of brightness increases due to glint
(Figure 10) indicate that a Y band crescent-phase glint enhancement could be
detected across a wide range of orbital inclinations and system distances. For
a typical HabEx target (at roughly 8 pc), the glint enhancement could be
detected for systems with inclinations above 60° in less than 200 hr. Glint
enhancements could be detected for exo-Earths on less-inclined orbits with
comparable integration times at distances below 5–7 pc. For the LUVOIR-A
concept, a glint enhancement could be detected for a target at 8 pc in 3–40 hr
for inclinations spanning 75°–30°, or, for a typical target at roughly 15 pc,
within 200 hr for inclinations above 45°.

Finally, ocean glint detection via phase-dependent spectral PCA offers an
approach to habitability characterization that is quite complementary to other
proposed methods. For example, potentially Earth-like exoplanets that
demonstrate a likely glint component in their phase-resolved spectra could
then become important targets for glint mapping (Lustig-Yaeger et al., 2018).
Alternatively, an exoplanet that shows a reflected-light polarization peak at
the Brewster angle geometry (which occurs at phase angles near 100–110°;
McCullough, 2006; Williams & Gaidos, 2008; Stam, 2008; Zugger et al., 2010)
would then be a strong candidate for follow-up crescent phase observations for
phase-dependent PCA studies. Overall, we envision that a large suite of
characterization approaches would be extensively applied to any future exo-
Earth candidate in an attempt to sort fact from fiction.

##  5 Conclusions

Detection of ocean glint signatures offers an important path to directly
constraining the surface habitability of any potentially Earth-like exoplanet.
Absent detailed, phase-resolved observations of the distant Earth, high-
fidelity, phase-dependent spectral models of Earth provide the best means for
understanding glint detectability. We explored how phase-dependent PCA could
reveal a glint component — with its characteristic reddening — over a range of
orbital inclinations using three-dimensional spectral models of Earth. Pairing
these spectral models with a direct imaging instrument noise model, we
explored the detectability of raw brightness enhancements due to glint. Our
findings are summarized as follows:

  * •

Application of spectral PCA to phase-dependent simulated observations of Earth
yields detections of a strong glint component across a wide range of viewing
geometries. For observations spanning gibbous through crescent phases, glint
is the dominant component explaining phase-dependent variations for orbital
inclinations larger than 60°and this component could likely be detected even
with low-SNR spectroscopy or photometry. At an inclination of 45°, only less-
extreme crescent phases can be accessed. Here, glint does not provide a
characteristic reddening and, as a result, PCA does not detect a strong glint
contribution.

  * •

A glint component can be detected in as few as two phase-dependent
observations for exo-Earth twins with orbital inclinations larger than 60°.

  * •

Exposure time required to detect a glint brightness enhancement is likely to
be driven by challenging crescent-phase observations. Such a detection is best
achieved in the I, Y, and J photometric bands, where limited water opacity
provides strong surface sensitivity for Earth. In Y band, the crescent-phase
glint enhancement could be detected for a typical HabEx target in less than
200 hr for orbital inclinations above 60°. The LUVOIR-A concept achieves
detections at this same distance in 3–40 hr, depending on orbital inclination,
and could detect the glint signature for a typical (more distant) target
within 200 hr for inclinations above 45°.

We thank N. Cowan for comments on an early draft of this work and for
insightful discussions on various approaches to principal component analysis.
J. Lustig-Yaeger kindly provided feedback on comparisons between work
presented here and previous results in the literature. Also, we thank M.
Marley for encouraging our writing on the utility of apparent albedo. DJR
acknowledges support from the Hooper Undergraduate Research Award program at
Northern Arizona University. TDR gratefully acknowledges support from NASA’s
Exoplanets Research Program (No. 80NSSC18K0349), Exobiology Program (No.
80NSSC19K0473), and Habitable Worlds Program (No. 80NSSC20K0226), as well as
the Nexus for Exoplanet System Science and NASA Astrobiology Institute Virtual
Planetary Laboratory (No. 80NSSC18K0829). Some results in this paper have been
derived using the HEALPix package (Górski et al., 2005).

## References

  * Abe & Matsui (1988) Abe, Y., & Matsui, T. 1988, Journal of Atmospheric Sciences, 45, 3081, doi: 10.1175/1520-0469(1988)045<3081:EOAIGH>2.0.CO;2
  * Aumann et al. (2003) Aumann, H. H., Chahine, M. T., Gautier, C., et al. 2003, IEEE Transactions on Geoscience and Remote Sensing, 41, 253, doi: 10.1109/TGRS.2002.808356
  * Battersby et al. (2018) Battersby, C., Armus, L., Bergin, E., et al. 2018, Nature Astronomy, 2, 596, doi: 10.1038/s41550-018-0540-y
  * Beer et al. (2001) Beer, R., Glavich, T. A., & Rider, D. M. 2001, Appl. Opt., 40, 2356, doi: 10.1364/AO.40.002356
  * Brack (1993) Brack, A. 1993, Origins of Life and Evolution of the Biosphere, 23, 3, doi: 10.1007/BF01581985
  * Cowan et al. (2012) Cowan, N. B., Abbot, D. S., & Voigt, A. 2012, ApJ, 752, L3, doi: 10.1088/2041-8205/752/1/L3
  * Cowan et al. (2009) Cowan, N. B., Agol, E., Meadows, V. S., et al. 2009, ApJ, 700, 915, doi: 10.1088/0004-637X/700/2/915
  * Fauchez et al. (2018) Fauchez, T., Arney, G., Kumar Kopparapu, R., & Domagal Goldman, S. 2018, AIMS Geosciences, 4, 180, doi: 10.3934/geosci.2018.4.180
  * Fresnel (1821) Fresnel, A. 1821, in Annales de Chimie et Physique, Vol. 17, 101–112
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485, doi: 10.1007/s11214-006-8315-7
  * Gaudi et al. (2018) Gaudi, B. S., Seager, S., Mennesson, B., et al. 2018, Nature Astronomy, 2, 600, doi: 10.1038/s41550-018-0549-2
  * Gaudi et al. (2020) Gaudi, B. S., Seager, S., Mennesson, B., et al. 2020, The Habitable Exoplanet Observatory (HabEx) Mission Concept Study Final Report.  https://arxiv.org/abs/2001.06683
  * Goldblatt et al. (2013) Goldblatt, C., Robinson, T. D., Zahnle, K. J., & Crisp, D. 2013, Nature Geoscience, 6, 661, doi: 10.1038/ngeo1892
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Hall et al. (1995) Hall, D. K., Riggs, G. A., & Salomonson, V. V. 1995, Remote Sensing of Environment, 54, 127, doi: 10.1016/0034-4257(95)00137-P
  * Hart (1978) Hart, M. H. 1978, Icarus, 33, 23, doi: 10.1016/0019-1035(78)90021-0
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108, doi: 10.1006/icar.1993.1010
  * Kopparapu et al. (2020) Kopparapu, R. K., Wolf, E. T., & Meadows, V. S. 2020, Characterizing Exoplanet Habitability (University of Arizona Press), 449–476.  http://www.jstor.org/stable/j.ctv105bb62.23
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131, doi: 10.1088/0004-637X/765/2/131
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., et al. 2013, A&A, 554, A69, doi: 10.1051/0004-6361/201321042
  * Livengood et al. (2011) Livengood, T. A., Deming, L. D., A’Hearn, M. F., et al. 2011, Astrobiology, 11, 907, doi: 10.1089/ast.2011.0614
  * Lustig-Yaeger et al. (2018) Lustig-Yaeger, J., Meadows, V. S., Tovar Mendoza, G., et al. 2018, AJ, 156, 301, doi: 10.3847/1538-3881/aaed3a
  * Madhusudhan & Burrows (2012) Madhusudhan, N., & Burrows, A. 2012, ApJ, 747, 25, doi: 10.1088/0004-637X/747/1/25
  * Matsui & Abe (1986) Matsui, T., & Abe, Y. 1986, Nature, 322, 526, doi: 10.1038/322526a0
  * McCullough (2006) McCullough, P. R. 2006, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0610518
  * Meadows & Barnes (2018) Meadows, V. S., & Barnes, R. K. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 57, doi: 10.1007/978-3-319-55333-7_57
  * Meadows & Crisp (1996) Meadows, V. S., & Crisp, D. 1996, J. Geophys. Res., 101, 4595, doi: 10.1029/95JE03567
  * Meadows et al. (2018) Meadows, V. S., Reinhard, C. T., Arney, G. N., et al. 2018, Astrobiology, 18, 630, doi: 10.1089/ast.2017.1727
  * Mennesson et al. (2019) Mennesson, B., Kennedy, G., Ertel, S., et al. 2019, Bulletin of the American Astronomical Society, 51, 324
  * Pallé et al. (2003) Pallé, E., Goode, P. R., Yurchyshyn, V., et al. 2003, Journal of Geophysical Research (Atmospheres), 108, 4710, doi: 10.1029/2003JD003611
  * Pierrehumbert & Gaidos (2011) Pierrehumbert, R., & Gaidos, E. 2011, ApJ, 734, L13, doi: 10.1088/2041-8205/734/1/L13
  * Qiu et al. (2003) Qiu, J., Goode, P. R., Pallé, E., et al. 2003, Journal of Geophysical Research (Atmospheres), 108, 4709, doi: 10.1029/2003JD003610
  * Quanz et al. (2018) Quanz, S. P., Kammerer, J., Defrère, D., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10701, Optical and Infrared Interferometry and Imaging VI, ed. M. J. Creech-Eakman, P. G. Tuthill, & A. Mérand, 107011I, doi: 10.1117/12.2312051
  * Ramirez & Kaltenegger (2014) Ramirez, R. M., & Kaltenegger, L. 2014, ApJ, 797, L25, doi: 10.1088/2041-8205/797/2/L25
  * Riggs et al. (1999) Riggs, G. A., Hall, D. K., & Ackerman, S. A. 1999, Remote Sensing of Environment, 68, 152, doi: 10.1016/S0034-4257(98)00107-2
  * Roberge & Moustakas (2018) Roberge, A., & Moustakas, L. A. 2018, Nature Astronomy, 2, 605, doi: 10.1038/s41550-018-0543-8
  * Robinson (2018) Robinson, T. D. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Springer), 67, doi: 10.1007/978-3-319-55333-7_67
  * Robinson et al. (2014) Robinson, T. D., Ennico, K., Meadows, V. S., et al. 2014, ApJ, 787, 171, doi: 10.1088/0004-637X/787/2/171
  * Robinson et al. (2010) Robinson, T. D., Meadows, V. S., & Crisp, D. 2010, ApJ, 721, L67, doi: 10.1088/2041-8205/721/1/L67
  * Robinson & Reinhard (2020) Robinson, T. D., & Reinhard, C. T. 2020, Earth as an Exoplanet (University of Arizona Press), 379–416.  http://www.jstor.org/stable/j.ctv105bb62.21
  * Robinson et al. (2016) Robinson, T. D., Stapelfeldt, K. R., & Marley, M. S. 2016, PASP, 128, 025003, doi: 10.1088/1538-3873/128/960/025003
  * Robinson et al. (2011) Robinson, T. D., Meadows, V. S., Crisp, D., et al. 2011, Astrobiology, 11, 393, doi: 10.1089/ast.2011.0642
  * Rugheimer et al. (2013) Rugheimer, S., Kaltenegger, L., Zsom, A., Segura, A., & Sasselov, D. 2013, Astrobiology, 13, 251, doi: 10.1089/ast.2012.0888
  * Schwieterman et al. (2015) Schwieterman, E. W., Robinson, T. D., Meadows, V. S., Misra, A., & Domagal-Goldman, S. 2015, ApJ, 810, 57, doi: 10.1088/0004-637X/810/1/57
  * Seager (2013) Seager, S. 2013, Science, 340, 577, doi: 10.1126/science.1232226
  * Shields et al. (2013) Shields, A. L., Meadows, V. S., Bitz, C. M., et al. 2013, Astrobiology, 13, 715, doi: 10.1089/ast.2012.0961
  * Stam (2008) Stam, D. M. 2008, A&A, 482, 989, doi: 10.1051/0004-6361:20078358
  * Thomas & Stamnes (1999) Thomas, G., & Stamnes, K. 1999, Radiative Transfer in the Atmosphere and Ocean, Cambridge Atmospheric and Space Science Series (Cambridge University Press).  https://books.google.com/books?id=ZVJfQgAACAAJ
  * Tinetti et al. (2006a) Tinetti, G., Meadows, V. S., Crisp, D., et al. 2006a, Astrobiology, 6, 34, doi: 10.1089/ast.2006.6.34
  * Tinetti et al. (2006b) —. 2006b, Astrobiology, 6, 881, doi: 10.1089/ast.2006.6.881
  * Turbet et al. (2018) Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86, doi: 10.1051/0004-6361/201731620
  * von Paris et al. (2010) von Paris, P., Gebauer, S., Godolt, M., et al. 2010, A&A, 522, A23, doi: 10.1051/0004-6361/201015329
  * Waters et al. (2006) Waters, J. W., Froidevaux, L., Harwood, R. S., et al. 2006, IEEE Transactions on Geoscience and Remote Sensing, 44, 1075, doi: 10.1109/TGRS.2006.873771
  * Way et al. (2016) Way, M. J., Del Genio, A. D., Kiang, N. Y., et al. 2016, Geophys. Res. Lett., 43, 8376, doi: 10.1002/2016GL069790
  * Whitmire et al. (1991) Whitmire, D. P., Reynolds, R. T., & Kasting, J. F. 1991, in Bioastronomy: The Search for Extraterrestial Life — The Exploration Broadens, ed. J. Heidmann & M. J. Klein, Vol. 390, 173–178, doi: 10.1007/3-540-54752-5_210
  * Williams & Gaidos (2008) Williams, D. M., & Gaidos, E. 2008, Icarus, 195, 927, doi: 10.1016/j.icarus.2008.01.002
  * Wolf & Toon (2010) Wolf, E. T., & Toon, O. B. 2010, Science, 328, 1266, doi: 10.1126/science.1183260
  * Wolfe (2006) Wolfe, R. E. 2006, in Science and Instruments, Vol. 1, Earth Science Satellite Remote Sensing, ed. J. J. Qu, W. Gao, M. Kafatos, R. E. Murphy, & V. V. Salomonson (Springer Berlin Heidelberg), 50–73
  * Wordsworth & Pierrehumbert (2013) Wordsworth, R., & Pierrehumbert, R. 2013, Science, 339, 64, doi: 10.1126/science.1225759
  * Zugger et al. (2010) Zugger, M. E., Kasting, J. F., Williams, D. M., Kane, T. J., & Philbrick, C. R. 2010, ApJ, 723, 1168, doi: 10.1088/0004-637X/723/2/1168
  * Zugger et al. (2011) —. 2011, ApJ, 739, 12, doi: 10.1088/0004-637X/739/1/12


11institutetext: University of St Andrews, St Andrews, Scotland
11email<EMAIL_ADDRESS>
# A Ptolemaic Partitioning Mechanism

Richard Connor

###### Abstract

For many years, exact metric search relied upon the property of triangle
inequality to give a lower bound on uncalculated distances. Two exclusion
mechanisms derive from this property, generally known as pivot exclusion and
hyperplane exclusion. These mechanisms work in any proper metric space and are
the basis of many metric indexing mechanisms. More recently, the Ptolemaic and
four-point lower bound properties have been shown to give tighter bounds in
some subclasses of metric space.

Both triangle inequality and the four-point lower bound directly imply
straightforward _partitioning_ mechanisms: that is, a method of dividing a
finite space according to a fixed partition, in order that one or more classes
of the partition can be eliminated from a search at query time. However, up to
now, no partitioning principle has been identified for the Ptolemaic
inequality, which has been used only as a filtering mechanism.

Here, a novel partitioning mechanism for the Ptolemaic lower bound is
presented. It is always better than either pivot or hyperplane partitioning.
While the exclusion condition itself is weaker than Hilbert (four-point)
exclusion, its calculation is cheaper. Furthermore, it can be combined with
Hilbert exclusion to give a new maximum for exclusion power with respect to
the number of distances measured per query.

###### Keywords:

Metric Search Partitioning Ptolemaic inequality Supermetric space

##  1 Background and Related Work

This article concerns querying a large finite space $(S,d)$ which is a subset
of an infinite metric space $(U,d)$.

In most general terms, querying the space $(S,d)$ with query $q\in U$ is the
task of finding a subset $\\{s\leftarrow S\,|\,d(q,s)\leq t\\}$, for some
value $t$ which gives a suitable size of solution set.

Table 1: Notation used throughout Symbols | Meaning
---|---
$(U,d)$ | An infinite metric space with domain $U$ and distance $d$
$(S,d)$ | A large finite space $S\subset U$ over which search is performed
$P$ | A small reference set $P\subset U$, usually $P\subset S$
$u,u_{0},u_{1},\dots$ | Elements of the infinite domain $U$
$s,s_{0},s_{1},\dots$ | Elements of the finite domain $S$
$p,p_{0},p_{1},\dots$ | Elements of the reference set $P$
$M$ | A fixed radius used to define a partition with a given $p\in P$
$q,t$ | A query $q\in U$ associated with a numeric query threshold $t$
$\mathcal{P}$ | A partition of $S$ defined according to distances to $P$
$\mathcal{S}$ | A class of $\mathcal{P}$ which may be excluded given a particular $q,t$
$\mathcal{Q}$ | A subset of $U$ defined according to a particular $q,t$
$\mathbb{R}^{n}$ | An $n$-dimensional real domain
$\ell_{2}$ | The Euclidean distance metric
$\tau$ | A numeric parameter of the Ptolemaic partitioning mechanism
| ($\tau\geq 0.5$, typically $\tau\approx 1$)

It is generally assumed that $|S|$ is large and the cost of applying the
function $d$ is high, and so the simple solution of applying $d(q,s)$ to all
$s\in S$ is intractable [9].

Table 1 gives a summary of these and other notations used thoughout the
article.

###  1.1 Filtering and Partitioning

All metric search solutions rely upon algebraic properties of $(U,d)$. A
relatively small set of distinguished reference points
$P=\\{p_{0},\dots,p_{m}\\}$ (typically, $P\subset S$) is used to avoid direct
calculation of $d(q,s)$, after the distances $d(s,P)$ and $d(q,P)$ have been
calculated. $d(s,P)$ is calculated ahead of query time, during a pre-
processing phase. Two types of usage are distinguished as follows:

filtering:

given a query $q\in U$, a specific datum $s\in S$, and the distances $d(q,P)$
and $d(s,P)$, it may be possible to determine that $d(q,s)>t$ for some $t$
without having to calculate $d(q,s)$.

partitioning:

given a partition $\mathcal{P}$ of $S$ determined at pre-processing time with
respect to $d(S,P)$, and the distances $d(q,P)$, it may be possible to
determine that some classes of $\mathcal{P}$ do not contain any elements $s$
such that $d(q,s)\leq t$.

Both types of mechanism have their place in metric search, see [1, 9] for many
examples. Filtering approaches however imply linear-time solutions, whereas
partitioning can be used to construct an indexing mechanism, typically where a
very large data set is recursively partitioned, in order to achieve a sub-
linear search time.

For filtering, the algebraic properties are required to give a lower-bound on
the distance $d(q,s)$ with reference to the sets of distances $d(q,P)$ and
$d(s,P)$. For partitioning, a further requirement is to identify a partition
that can be determined at pre-processing time, of which one or more classes
may be excluded at query time according to $d(q,P)$.

Table 2 shows partitioning mechanisms which derive from various known lower-
bound properties. The contribution of this paper is a novel partitioning
mechanism for Ptolemaic inequality, shown in bold type in the table. Until
now, such a partitioning mechanism has been missing from the literature.

Table 2: Five different partition functions and their corresponding exclusion conditions. In all cases the Partition Criterion is used to form a distinguished subset of $S$ at pre-processing time, for all $s\in S$. The Exclusion Condition is evaluated with respect to the query $q$ and a query radius $t$. Row 5 summarises the novel contribution of this paper. | Underlying property | Partition Criterion | Exclusion Condition
---|---|---|---
1 | triangle inequality | $d(s,p)\leq M$ | $d(q,p)>M+t$
2 | triangle inequality | $d(s,p)\geq M$ | $d(q,p)<M-t$
3 | triangle inequality | $d(s,p_{0})\leq d(s,p_{1})$ | $d(q,p_{0})-d(q,p_{1})>2t$
4 | four-point lower bound | $d(s,p_{0})\leq d(s,p_{1})$ | $\frac{d(q,p_{0})^{2}-d(q,p_{1})^{2}}{d(p_{0},p_{1})}>2t$
5 | Ptolemaic inequality |  $d(s,p_{0})\leq d(s,p_{1})$ | $d(q,p_{0})-d(q,p_{1})>t/\tau$
|  | $\land\,\,d(s,p_{1})\geq\tau d(p_{0},p_{1})$ |

The remainder of this section introduces some necessary preliminaries. In
Section 2 the underlying geometry of the partition mechanism is given, and
Section 3 gives its formal definition. Section 4 gives a quantitative analysis
of its value.

###  1.2 Subclasses of Metric Space

Properties (1-3) listed in Table 2 are possessed by all proper metric spaces.
Property (4) is found only in _supermetric_ spaces [2], which include all
spaces which are isometrically embeddable in Hilbert space111see eg
https://en.wikipedia.org/wiki/Hilbert_space, while property (5) is found in
any Hadamard space222see eg https://en.wikipedia.org/wiki/Hadamard_space.

Any Hilbert-embeddable space is also a Hadamard space; although Hadamard
spaces are a little more general, it is not clear that any practical non-
Hilbert spaces fall in this category. Details of Hilbert spaces are elaborated
in [3]; in this context it is sufficient to know that the following classes of
metric space are members of both classes: Euclidean, Cosine, Jensen-Shannon,
Quadratic Form, Triangular, and Mahalanobis spaces. Furthermore, the square
root of any proper metric gives a space in both classes. The partition
mechanism described here is thus applicable to any of these spaces.

###  1.3 Ptolemaic and Four-point Lower Bounds

The Ptolemaic inequality was identified for use as a distance lower-bound for
certain metric spaces in [5], and used further in a number of studies for
example [6, 7]. For any four objects $u_{0},u_{1},u_{2},u_{3}\in U$, the
Ptolemaic inequality states:

$d(u_{0},u_{2})\cdot d(u_{1},u_{3})\leq d(u_{0},u_{1})\cdot
d(u_{2},u_{3})+d(u_{1},u_{2})\cdot d(u_{3},u_{0})$

In $(\mathbb{R}^{n},\ell_{2})$ this is more simply stated as the product of
the lengths of the diagonals of any quadrilateral being no greater than the
sum of the products of the pairs of opposing sides. Given this property, a
lower bound on the distance $d(q,s)$ can be determined whenever, for two
reference values $p_{0},p_{1}$, all the distances
$d(s,p_{0}),d(s,p_{1}),d(q,p_{0}),d(q,p_{1})$ and $d(p_{0},p_{1})$ are known.
This lower bound is much tighter than those available via simple triangle
inequality, and has been used to great effect for filtering objects during
search, particularly in the context of a very expensive distance function [7].
The mechanisms used to incorporate this lower bound into metric search
techniques include the Ptolemaic pivot table, the Ptolemaic PM-Tree, and the
Ptolemaic M-Index [6]. In all cases, the inequality is used as an extra
filtering mechanism superimposed onto an existing filtering or partitioning
structure.

The four-point lower bound property, and the Hilbert exclusion mechanism, were
first identified in [2], and investigated further in [3, 4]. Any supermetric
space $(U,d)$ has the four-point property: for any four objects
$u_{0},u_{1},u_{2},u_{3}\in U$, there exists a tetrahedron with vertices
$u_{0}^{\prime},u_{1}^{\prime},u_{2}^{\prime},u_{3}^{\prime}\in\mathbb{R}^{3}$
where the distances between pairs of points are preserved, i.e.
$d(u_{i},u_{j})=\ell_{2}(u_{i}^{\prime},u_{j}^{\prime})$.

The four-point property thus implies the Ptolemaic property, but not vice-
versa.

The four-point lower-bound property applies to the case where five of the six
edge lengths of a tetrahedron are known. In this case, two adjacent faces of
the tetrahedron can be constructed. A lower bound of the final side length is
obtained by notionally rotating these faces around their common edge to
minimise the final edge length, which occurs when a planar tetrahedron is
formed.

###  1.4 Projection into 2-dimensional Space

Figure 1: Four objects $p_{0},p_{1},s_{0},s_{1}$ selected from a supermetric
space are projected onto a 2D plane according to the known distances
$K,A,B,C,D$. Although X is not known, it is known that a tetrahedron with
these four vertices exists in 3 dimensions. By the four-point lower bound,
$d(s_{0},s_{1})\geq X$. By the Ptolemaic lower bound in 2 dimensions,
$X\geq\frac{AD-BC}{K}$

Together these properties imply that if the Ptolemaic inequality is applied to
a quadrilateral in two dimensions, when that quadrilateral has been formed
according to five distances measured among four objects in a supermetric
space, then the inequality applies also to the original space. Figure 1 shows
an example of this.

The figure shows a projection in $(\mathbb{R}^{2},\ell_{2})$ of four objects
$p_{0},p_{1},q$ and $s$ selected from a supermetric space $(U,d)$. All
distances other than $d(q,s)$ have been calculated in $(U,d)$. The projections
of objects $p_{0}$ and $p_{1}$ are plotted at the points $(0,0)$ and
$(0,d(p_{0},p_{1}))$ respectively333This choice is arbitrary, any two points
which preserve $d(p_{0},p_{1})$ could be used. The projections of objects $q$
and $s$ are plotted at the unique points above the X-axis which preserve their
distances from $p_{0}$ and $p_{1}$. The supermetric properties imply that the
tetrahedron $p_{0},p_{1},q,s$ must exist in $(\mathbb{R}^{3},\ell_{2})$,
therefore the unknown distance $d(q,s)$ is lower-bounded by the sixth edge of
the planar tetrahedron plotted in $(\mathbb{R}^{2},\ell_{2})$.

By the four-point lower bound property, $d(q,s)\geq X$. By the Ptolemaic
lower-bound property, $X\geq\frac{AD-BC}{K}$. Therefore, in the original
supermetric space, $d(q,s)\geq\frac{AD-BC}{K}$. For the rest of this article,
only 2D projections like these are considered, where two distinguished
reference objects $p_{0},p_{1}$ are used to form a planar projection of the
rest of the data set, and rely on the Ptolemaic property with the context of
planar quadrilaterals. This restricts the outcome to Hilbert-embeddable
spaces, although as noted this is not a significant practical restriction.

It is worth noting that while the derivation and correctness of the mechanism
rely upon the existence of the 2D projection, the projection itself does not
require to be calculated. As shown in Table 2, the calculations required are
restricted to simple calculations over distances measured in the original
space.

##  2 The Underlying Geometry

Partitioning mechanisms differ from filtering in that, for each possibility of
exclusion, it is necessary to identify two subsets of the universal space:

  1. 1.

a _static_ subset $\mathcal{S}$, which can be identified and indexed during
the pre-processing of the finite data set, and

  2. 2.

a _dynamic_ subset $\mathcal{Q}$, which is identified only after the query and
(typically) its associated search radius become apparent.

Exclusion of $\mathcal{S}$ can be performed when every element of
$\mathcal{Q}$ is separated by at least the search radius from every element of
$\mathcal{S}$. In the following section objects denoted by $s$ and $q$ are
referred to, representing elements of $\mathcal{S}$ and $\mathcal{Q}$
respectively.

###  2.1 2D Geometry

(a) Four points $p_{0},p_{1},s,q$ on a 2D plane.

(b) For any $q^{\prime}$ in the shaded area, $d(q^{\prime},s)\geq\tau|A-B|$.

Figure 2:  In 2(a), line segments among the points are annotated by their
lengths. $s$ is chosen according to the parameter $\tau=C/K$, and $q$ is a
point on the circle defined by $p_{0},p_{1}$ and $s$. $q$ is the unique point
with $X=\tau(B-A)$. The regions shown in the figure represent regions of the
original, potentially non-Euclidean, space; the inequalities established are
generally applicable to the original space.

Figure 2(a) shows four points $p_{0},p_{1},q$ and $s$ drawn on a plane. These
points represent the 2D projections of two reference (or pivot) values
$p_{0},p_{1}$, a query value $q$ and a potential solution value $s$. The
figure is annotated with line segments labelled $A-D$, $K$ and $X$, where the
labels represent the lengths of the respective lines. $K$ is the inter-pivot
distance, and $X$ is a lower bound of the unknown distance $d(q,s)$.

For the moment, values have been chosen such that

  * •

$C=D$

  * •

the parameter $\tau$ defines the ratio $C/K$

  * •

point $q$ lies on the same circle as $p_{0},p_{1}$ and $s$

The Ptolemaic inequality states

$BC\leq AD+KX$

so in this case:

$X\geq\tau(B-A)$

The boundary of this region defines a hyperbola with foci $p_{0},p_{1}$ and
semi-major axis $X/2\tau$, as shown in Figure 2(b). It follows that any point
within the shaded region is at least distance $X$ from the point $s$.

As $q,p_{0},p_{1}$ and $s$ are co-circular, $q$ is the unique closest point on
the (left-hand) hyperbola to $s$ and $X=\tau(B-A)$. The line segment $sq$ is
therefore perpendicular to the tangent of the hyperbola at $q$. As the
gradient of the tangent is negative, it follows that any point above and to
the right of $s$ is further than $X$ from any point to the left of the
hyperbola, as illustrated in Figure 3(a).

(a) The shaded area to the right of the central axis contains points which are
at least $d(q,s)$ from that on the left.

(b) The shaded area here also contains contains points which are at least
$d(q,s)$ from that on the left.

Figure 3: The shaded area to the left of each central axis denotes the locus
defined by $d(q^{\prime},p_{0})-d(q^{\prime},p_{1})>\tau t$ for any
$q^{\prime}\in U$. In both cases, for any $s^{\prime}\in S$ in the shaded area
to the right of the central axis, $d(q^{\prime},s^{\prime})>t$.

However, the static partition can be extended to include more of the finite
search space, by including any value $s^{\prime}\in S$ to the right of the
central axis where also $d(s^{\prime},p_{0})\geq d(s,p_{0})$, as illustrated
in Figure 3(b). This not only increases the cardinality of the potentially
excluded subset, but also avoids the requirement to calculate the 2D
projection.

The static and dynamic classes represented in Figure 3(b) are now formally
defined as

$\displaystyle\mathcal{S}\quad$ $\displaystyle=\quad\\{s^{\prime}\leftarrow
S\,|\,d(s^{\prime},p_{0})\geq d(s^{\prime},p_{1})\quad\land\quad
d(s^{\prime},p_{0})\geq\tau d(p_{0},p_{1})\\}$ $\displaystyle\mathcal{Q}\quad$
$\displaystyle=\quad\\{q^{\prime}\leftarrow
U\,|\,d(q^{\prime},p_{1})-d(q^{\prime},p_{0})>t/\tau\\}$

with the property that it is impossible for any element
$q^{\prime}\in\mathcal{Q}$ to be within distance $t$ of any element
$s^{\prime}\in\mathcal{S}$.

The validity of this extension seems evident from the illustration in Figure
3(b), but of course needs to be demonstrated for the general case. A full
justification of the correctness is included in Appendix 0.A.

##  3 The Partition Mechanism

The addition of the criterion $d(s^{\prime},p_{0})\geq\tau d(p_{0},p_{1})$ to
the static partition allows further exclusion potential relying on the normal
triangle inequality method, i.e. if $d(q,p_{0})<\tau d(p_{0},p_{1})-t$. This
extension to the exclusion criterion is illustrated in Figure 4(a).

(a) The locus $d(p_{0},q)<d(p_{0},s)-t$ can be included in $\mathcal{Q}$.

(b) The locus $d(s,p_{\\{0,1\\}})<\tau d(p_{0},p_{1})$ forms the final
partition class.

Figure 4:

Furthermore, when also including the symmetric opposite criteria, the static
partition now defines three subclasses as shown in Figure 4(b). It may further
be noted that the third of these subclasses may also independently excluded if
$d(q,p_{0})$ or $d(q,p_{1})\geq\tau K+t$, again relying only on triangle
inequality.

So finally, according to the geometry established in Section 2.1, a static
partition $\\{\mathcal{S}_{1},\mathcal{S}_{2},\mathcal{S}_{3}\\}$ of $S$ can
be established for a pair of reference points $p_{0}$ and $p_{1}$ with
$K=d(p_{0},p_{1})$ and a given value of $\tau$ as follows:

$\displaystyle\mathcal{S}_{1}\quad$ $\displaystyle=\quad\\{s\leftarrow
S\,|\,d(s,p_{0})\geq d(s,p_{1})\quad\land\quad d(s,p_{0})\geq\tau K\\}$
$\displaystyle\mathcal{S}_{2}\quad$ $\displaystyle=\quad\\{s\leftarrow
S\,|\,d(s,p_{0})<d(s,p_{1})\quad\land\quad d(s,p_{1})\geq\tau K\\}$
$\displaystyle\mathcal{S}_{3}\quad$ $\displaystyle=\quad\\{s\leftarrow
S\,|\,d(s,p_{0})<\tau K\quad\land\quad d(s,p_{1})<\tau K\\}$

These static regions are illustrated on the 2D plane in Figure 4(b).

For a given query object $q$ with threshold $t$, where $A=d(q,p_{0})$ and
$B=d(q,p_{1})$, these classes can be excluded from a search as follows:

$\displaystyle\mathcal{S}_{1}\quad$ $\displaystyle:\quad
B-A>t/\tau\quad\lor\quad A<\tau K-t$ $\displaystyle\mathcal{S}_{2}\quad$
$\displaystyle:\quad A-B\geq t/\tau\quad\lor\quad B<\tau K-t$
$\displaystyle\mathcal{S}_{3}\quad$ $\displaystyle:\quad A\geq\tau
K+t\quad\lor\quad B\geq\tau K+t$

Note that it is possible for the exclusion of region $\mathcal{S}_{3}$ to
occur in conjunction with that of $\mathcal{S}_{3}$ or $\mathcal{S}_{3}$. The
mechanism resulting from these definitions is now evaluated in Section 4.

##  4 Evaluation

Before proceeding with a full quantitative evaluation, it is interesting to
view graphical representations based on a sample from a particular data set,
in order to give a more pragmatic view of the Ptolemaic partition mechanism in
comparison with hyperplane (hyperbolic) and Hilbert (four-point) partition
mechanisms.

(a) A scatter plot based on randomly selected reference points. $p_{0}$ is
plotted at $(0,0)$ and $p_{1}$ at $(0,d(p_{0},p_{1}))$. The data set is
plotted according to the distance of each value from $p_{0}$ and $p_{1}$.

(b) The boundary of the class $\mathcal{S}_{1}$ is plotted in red. Values
lying to the right of this boundary can be excluded for queries lying to the
left of either black boundary.

Figure 5: Graphical view of the Ptolemaic partition mechanism based on a 2D
projection. A value of $\tau=1.3$ has been arbitrarily selected, along with a
query threshold of 0.3, which is the mean nearest-neighbour distance within a
set of 100k objects.

(a)

(b)

Figure 6: Equivalent graphical views of hyperplane and Hilbert partition
mechanisms. In each case, queries falling to the left of the black line can be
used to exclude the subset of data falling to the right of the red line. Again
a query threshold of 0.3 has been used.

Figure 5(a) shows a scatter plot of 1,000 generated values from a
10-dimensional Euclidean space, each projected onto a 2D plane according to
their distances from two randomly generated pivot values.

Figure 5(b) shows the same projection superimposed with one of the partition
boundaries of the Ptolemaic partition mechanism with a $\tau$ value of 1.3,
and a query threshold of 0.3. The boundary of the static region
$\mathcal{S}_{1}$ is shown in red; those points lying to the right of the red
boundary are thus subset to exclusion when either
$d(q,p_{1})-d(q,p_{0})>1.3\cdot 0.3$, or if $d(q,p_{1})<1.3\cdot
d(p_{0},p_{1})-0.3$. The boundaries of these regions are shown by solid and
dotted black lines respectively; every point to the left of either boundary
represents a value for $q^{\prime}\in\mathcal{Q}$ which allows exclusion of
the static class.

Figure 6 shows the same plot with boundaries for standard hyperplane and
Hilbert exclusion, in 6(a) and 6(b) respectively.

It is clear, at least in this case, that the Ptolemaic mechanism always
excludes a smaller subset, while the probability of the exclusion being
possible is higher. From Figure 6(a) it is evident that the dimensionality of
the data set is starting to challenge hyperplane exclusion, while both
Ptolemaic and Hilbert mechanisms continue to remain effective. Finally, while
it is not possible to judge the relative efficacy of Ptolemaic vs. Hilbert
from these diagrams, it can be observed that neither is a proper subset of the
other, and it is therefore possible to use both Ptolemaic and Hilbert with
respect to the same pair of reference points. This would allow a hybrid
mechanism, more effective that either in isolation, based on the same dynamic
measurements of $d(q,p_{0})$ and $d(q,p_{1})$.

###  4.1 Quantitative Evaluation

Quantitative evaluation is performed over sets of uniformly generated
Euclidean data, from between 8 and 20 dimensions. 50k data objects are used
and 1k non-intersecting queries are evaluated. The threshold used for each
query corresponds to the 5nn distance as pre-calculated over the data.

Experiments were performed over Ptolemaic, Hilbert, and hyperplane mechanisms.
For each experiment, a fixed number of reference points was used, and each of
the $n\choose 2$ pairs of reference points was used to construct a partition
over the space. The single outcome is the mean proportion per query of values
that were successfully excluded, this value being between 0 and 1. For the
majority of the experiments 10 reference points used, this giving 45 different
partitions. Thus all results given correspond to the proportion of the data
that can be successfully excluded at cost of only 10 distance calculations per
query.

All experiments were performed using MatLab, and the code is available from
the author444note to reviewers, will be published!.

###  4.2 Choosing $\tau$

Figure 7: Probability of successful exclusion for differing values of $\tau$
in different dimensions. Note that when $\tau=0.5$ the mechanism is identical
to traditional hyperplane (Hyperbolic) exclusion.

First, different values for $\tau$ are examined. As mentioned, when $\tau=0.5$
the mechanism reverts to simple hyperplane exclusion; while a value of less
than $0.5$ is technically possible, there is no value in such a choice. As
$\tau$ gets large, then ever fewer data will be present in the partition which
may be excluded, and again the mechanism will become useless. Early tests
showed that a value somewhere around $1$ is usually close to optimal, although
for specific reference point pairs an optimum values of between around 0.8 and
1.2 were observed.

It would in fact be possible to optimise $\tau$ based on each particular pair
of reference points, which we have not yet investigated thoroughly. In this
experiment a fixed value of $\tau$ is applied to all partitions, which is
possibly more realistic for many scenarios.

Figure 7 shows the results of various values of $\tau$ when applied to data of
between 8 and 20 dimensions. As can be seen there is a general trend of larger
values being better as dimensions increase, but only within quite a small
margin; while there is clearly an element of noise in this experiment, the
best value in each case is either 1.0 or 1.1, although further investigation
is warranted. For further experiments described over the data of different
dimensions, the best value of $\tau$ found in this experiment was used.

###  4.3 Evaluation over High Dimensional Data

Figure 8:  Performance of different exclusion mechanisms as dimensionality
increases, measured as the probability of a non-solution being excluded based
on the distances among query and reference values.

Having picked a value for $\tau$, outcomes for some different mechanisms over
data ranging from 8 to 20 dimensions are given. Four mechanisms are used:
Hyperplane (Hyperbolic) exclusion; Hilbert exclusion; Ptolemaic exclusion, and
finally a combination of Hilbert and Ptolemaic exclusion.

For the combination, for each pair of reference points five different subsets
were identified at pre-processing time: two for Hilbert, and three for
Ptolemaic, and all possible exclusions were attempted for each query. The
observation here is that it is possible for different data to be excluded by
each mechanism. As the essential query-time cost of performing the exclusions
is the cost of the two distance operations and some relatively cheap
arithmetic, taking the union of all possible exclusions makes practical sense
as the distance calculations are amortised. As Hilbert exclusion always allows
exclusion from a superset of queries identified by Hyperbolic exclusion, there
is no point in combining that mechanism also.

Figure 8 shows the outcome. As can be seen, while the performance of
Hyperbolic exclusion falls rapidly away after around 8-10 dimensions, both
Ptolemaic and Hilbert perform much better into the higher dimensional range.
Hilbert always performs better than Ptolemaic, which is not very surprising as
the four-point lower bound property is stronger then the Ptolemaic inequality,
and technically applies to a smaller subset of metric spaces. What is more
interesting, however, is that the combination of Ptolemaic and Hilbert gives a
strictly better result than Hilbert alone; that is, the data sets identified
for exclusion by the two mechanisms are not in a strict subset relation. Again
it is noted that the inherent query-time cost of the joint mechanism is very
similar to the cost of just one, as in all cases the query to pivot distances
calculated are reused in both mechanisms.

The final plots in the graph show the use of 20 and 50 pivot values for the
combination mechanism. Although only doubling the number of query-to-pivot
distances required, 20 pivots gives $20\choose 2$ i.e. 190 partitions to
apply, and as can be seen the increase makes for a much higher exclusion
ratio. Similarly, 50 pivots gives 1,225 partitions. The important observation
however is that there is a clear degree of orthogonality in the randomly
selected partitions, allowing almost perfect exclusion in 12 and 16 dimensions
respectively.

##  5 Conclusions and Future Work

This paper fills a significant gap in the literature, that is a set partition
that can be used as an exclusion mechanism for the Ptolemaic inequality; for
some years, other distance lower-bounds have had known mechanisms and in this
sense the Ptolemaic inequality has been an outlier.

In its simplest form, the mechanism is quantitatively much better than
traditional hyperplane partitioning, and not quite as good as Hilbert
partitioning. This is almost inevitable, as the class of spaces to which the
inequalities can be applied are in a strict subset relation. Should this
mechanism have been identified before Hilbert exclusion it would have been
deserving of significant excitement, but this is nowadays tempered by the
existence of the more effective Hilbert exclusion over essentially the same
subclass of metric spaces.

However, it is the case that the individual data objects which the new
mechanism excludes are not a proper subset of those identified by Hilbert
exclusion, and as shown the two mechanisms may operate in conjunction to give
a unified mechanism which, for the same cost of distance calculations against
reference points, gives a better exclusion outcome than either in isolation.
Particularly in high-dimensional spaces, this therefore gives a further
increment in the limit of dimensionality for which exact search can be
effective. While the “rule of thumb” used to be that 8-10 dimensions was the
effective limit for exact search [8], with the combined mechanism 16
dimensions can be effectively searched while avoiding almost all explicit
distance calculations.

Some further avenues are worth exploring. First, it is feasible to calculate
individual $\tau$ values customised to each particular pair of pivot points,
rather than to choose a single value for the whole set. This would be expected
to give significant, if incremental, improvement in performance.

Finally, there are many other contexts beyond a simple recursive decomposition
of a large data set where such mechanisms can be used. It is therefore of
potential value in its own right for this previous gap in knowledge to be
filled.

##  6 Acknowledgements

The author would like to sincerely thank the anonymous reviewers for their
thorough and helpful comments on the submitted version of this article.

## References

  * [1] Chávez, E., Navarro, G., Baeza-Yates, R., Marroquín, J.L.: Searching in metric spaces 33(3), 273–321. https://doi.org/10.1145/502807.502808, https://dl.acm.org/doi/10.1145/502807.502808
  * [2] Connor, R., Vadicamo, L., Cardillo, F.A., Rabitti, F.: Supermetric search with the four-point property. In: Proceedings of 9th International Conference on Similarity Search and Applications (SISAP 2016). pp. 51–64. LNCS, Springer (2016)
  * [3] Connor, R., Cardillo, F.A., Vadicamo, L., Rabitti, F.: Hilbert exclusion: improved metric search through finite isometric embeddings. ACM Transactions on Information Systems (TOIS) 35(3), 1–27 (2016)
  * [4] Connor, R., Vadicamo, L., Cardillo, F.A., Rabitti, F.: Supermetric search. Information Systems 80, 108–123 (2019)
  * [5] Hetland, M.L.: Ptolemaic indexing. Journal of Computational Geometry 6, 165–184–165–184 (2015)
  * [6] Hetland, M.L., Skopal, T., Lokoč, J., Beecks, C.: Ptolemaic access methods: Challenging the reign of the metric space model. Information Systems 38(7), 989–1006 (2013). https://doi.org/https://doi.org/10.1016/j.is.2012.05.011, https://www.sciencedirect.com/science/article/pii/S0306437912000786
  * [7] Lokoč, J., Hetland, M.L., Skopal, T., Beecks, C.: Ptolemaic indexing of the signature quadratic form distance. In: Proceedings of the Fourth International Conference on SImilarity Search and APplications. p. 9–16. SISAP ’11, Association for Computing Machinery, New York, NY, USA (2011). https://doi.org/10.1145/1995412.1995417, https://doi.org/10.1145/1995412.1995417
  * [8] Weber, R., Schek, H.J., Blott, S.: A quantitative analysis and performance study for similarity-search methods in high-dimensional spaces. In: Proc. of 24th VLDB. vol. 98, pp. 194–205. Morgan Kaufmann (1998)
  * [9] Zezula, P., Amato, G., Dohnal, V., Batko, M.: Similarity search: the metric space approach, vol. 32. Springer Science & Business Media (2006)

## Appendix

##  Appendix 0.A Justification of Correctness

Figure 9: The shaded area to the left of each central axis denotes the locus
defined by $d(q^{\prime},p_{0})-d(q^{\prime},p_{1})>\tau t$ for any
$q^{\prime}\in U$. In both cases, for any $s^{\prime}\in S$ in the shaded area
to the right of the central axis, $d(q^{\prime},s^{\prime})>t$.

Figure 9 is reproduced from the text of Section 3. In the following we define
the shaded regions as:

$\displaystyle\mathcal{S}\quad$ $\displaystyle=\quad\\{s^{\prime}\leftarrow
S\,|\,d(s^{\prime},p_{0})\geq d(s^{\prime},p_{1})\quad\land\quad
d(s^{\prime},p_{0})\geq\tau d(p_{0},p_{1})\\}$ $\displaystyle\mathcal{Q}\quad$
$\displaystyle=\quad\\{q^{\prime}\leftarrow
U\,|\,d(q^{\prime},p_{1})-d(q^{\prime},p_{0})>t/\tau\\}$

We now show that for any $q^{\prime}\in\mathcal{Q}$ and
$s^{\prime}\in\mathcal{S},d(q^{\prime},s^{\prime})>d(q,s)$555Note that
$s\in\mathcal{S}$ but $q\notin\mathcal{Q}$..

There are two circles to consider, illustrated in Figure 10(a). $\mathcal{C}$
is centred around a point on the hyperbola with radius $d(q,s)$, and
$\mathcal{C}^{\prime}$ is centred around $p_{0}$ with radius $d(p_{0},s)$. So
far we have shown that $q$ is the nearest point on the hyperbola to $s$, and
therefore $s$ lies on $\mathcal{C}$ when it is centred at $q$. It is
sufficient to show that $\mathcal{C}$ does not intersect with the boundary of
$\mathcal{S}$ when its centre is at any point on the hyperbola containing $q$.

(a) Circle $\mathcal{C}$ and arc $\mathcal{C^{\prime}}$ must not intersect
when $\mathcal{C}$ is centred at or below $q$.

(b) The gradient of the hyperbola at $q^{\prime}$ is equal to the gradient of
$\mathcal{C^{\prime}}$ at $s^{\prime\prime}$.

Figure 10:

If $\mathcal{C}$ is centred around $q$, as illustrated in Figure 10(a), it
crosses the central axis. However as $p_{0},q,s,$ and $p_{1}$ are co-circular,
then (a) its radius is smaller than that of $\mathcal{C}^{\prime}$, and (b)
line $p_{0}s$ has a steeper gradient that line $qs$. Together these ensure
that no element of $\mathcal{S}$ is within distance $d(q,s)$ of $q$.

Considering the case where the centre of $\mathcal{C}$ is above $q$, it is
sufficient to observe that, as $q$ is the closest point on the hyperbola to
$s$, then the segment $sq$ is perpendicular to the tangent of the hyperbola at
$q$. The gradient of the tangent is negative, and so as the centre of
$\mathcal{C}$ moves upwards from $q$, the distance from the centre to the
nearest point on the boundary becomes increasingly greater than $d(q,s)$.

Now consider $\mathcal{C}$ as its centre adopts some position $q^{\prime}$ on
the hyperbola between $q$ and the line segment $p_{0}p_{1}$, as depicted in
Figure 10(b). Define $s^{\prime}$ as the point on $\mathcal{C}^{\prime}$
intersecting the line $p_{0}q^{\prime}$. As long as the gradient of the
hyperbola at any point between $q$ and $q^{\prime}$ is less than the gradient
of the arc at any point between $s$ and $s^{\prime}$, then it is impossible
for any intersection to occur between $\mathcal{C}$ and
$\mathcal{C}^{\prime}$.

The point $x$ is marked as the point on the central axis closest to the point
$q^{\prime}$, i.e. where the points $p_{0},q^{\prime},x$ and $p_{1}$ are co-
circular. The point $s^{\prime\prime}$ is marked as the point on
$\mathcal{C}^{\prime}$ such that the lines segments $q^{\prime}x$ and
$p_{0}s^{\prime\prime}$ are parallel. The gradient of the tangent of the
hyperbola at $q^{\prime}$ is thus equal to the gradient of the tangent of
$\mathcal{C}^{\prime}$ at $s^{\prime\prime}$.

It is now sufficient to note that the gradient of $q^{\prime}x$ is less than
that of $p_{0}x$, due to the co-circularity of $p_{0},q^{\prime},x$ and
$p_{1}$. Therefore $s^{\prime\prime}$ is below $s^{\prime}$ on
$\mathcal{C}^{\prime}$, in all cases where $q$ and $x$ are above the line
segment $p_{0}p_{1}$. Therefore $\mathcal{C}$ and $\mathcal{C}^{\prime}$ can
never intersect, for any value of $\tau$.