# Fast Active Monocular Distance Estimation from Time-to-Contact

Levi Burner , Nitin J. Sanket, Cornelia Fermüller, Yiannis Aloimonos
Perception and Robotics Group, University of Maryland Institute for Advanced
Computer Studies
University of Maryland, College Park
Corresponding author. The support of the NSF under awards DGE-1632976 and OISE
2020624 is gratefully acknowledged.

###### Abstract

Distance estimation is fundamental for a variety of robotic applications
including navigation, manipulation and planning. Inspired by the mammal’s
visual system, which gazes at specific objects (active fixation), and
estimates when the object will reach it (time-to-contact), we develop a novel
constraint between time-to-contact, acceleration, and distance that we call
the $\tau$-constraint. It allows an active monocular camera to estimate depth
using time-to-contact and inertial measurements (linear accelerations and
angular velocities) within a window of time.

Our work differs from other approaches by focusing on patches instead of
feature points. This is, because the change in the patch area determines the
time-to-contact directly. The result enables efficient estimation of distance
while using only a small portion of the image, leading to a large speedup.

We successfully validate the proposed $\tau$-constraint in the application of
estimating camera position with a monocular grayscale camera and an Inertial
Measurement Unit (IMU). Specifically, we test our method on different real-
world planar objects over trajectories 8-40 seconds in duration and 7-35
meters long. Our method achieves 8.5 cm Average Trajectory Error (ATE) while
the popular Visual-Inertial Odometry methods VINS-Mono and ROVIO achieve 12.2
and 16.9 cm ATE respectively. Additionally, our implementation runs 27$\times$
faster than VINS-Mono’s and 6.8$\times$ faster than ROVIO’s. We believe these
results indicate the $\tau$-constraints potential to be the basis of robust,
sophisticated algorithms for a multitude of applications involving an active
camera and an IMU.

Figure 1: Top: Five seconds of sequence 9’s camera trajectory along with the
fixated scene patch used to estimate 3D distance using a monocular camera and
an IMU. Bottom: $l_{2}$ error for Seq. 9 of our method compared to VINS-Mono,
ROVIO, and AprilTag 3. Averaged over all sequences, our method is 50% and 30%
more accurate and 27$\times$ and 6.8$\times$ faster than VINS-Mono and ROVIO
while also using only a small portion of the image (initially 2.5%). Time is
color coded with the hot colormap.

##  1 Introduction

Early researchers in computer vision were fascinated by living being’s ability
to control their movement in order to gather information about their
environment. The process was named “Active Perception” [1, 3, 6, 11], and
numerous mechanisms for exploiting activeness were developed over the past
thirty years. In [18], four forms of activeness are described, three of which
require physical movement. It is in this sense, that we use the word “active”,
_e.g_. to mean “movement at will” or “moving” rather than active illumination
like RGB-D cameras. The goal of all forms is to reduce the computational load
on the system by distributing the task at hand between action and perception
[29].

Because, active vision systems are usually moving in some way, they must
accelerate to produce changes to this movement. Despite this, utilizing
observer acceleration (measured using an Inertial Measurement Unit or IMU) to
facilitate visual computations has not received much attention in the Active
Vision literature. To fill this gap, we introduce a new mathematical
constraint relating (a) time-to-contact ($\tau$) (the ratio of camera velocity
to scene distance) with (b) scene depth and (c) observer acceleration. This
constraint makes it possible for an active monocular camera to estimate scene
depth by accelerating in any direction. We call this the $\tau$-constraint.

We demonstrate the utility of the constraint in a series of experiments
involving fixation on a single object (i.e, tracking) while the active
(moving) monocular camera is accelerating. Fixation is key to the method’s
efficiency because it can ignore most of the incoming pixels. Such object-
centered representations have been argued to be useful for a wide variety of
tasks ranging from manipulation to navigation [6, 4, 12].

The key concept of fusing inertial measurements with camera observations has
been extensively studied in the fields of Structure from Motion (SfM) and
Simultaneous Localization and Mapping (SLAM). Traditionally, this fusion has
been achieved in Visual Inertial Odometry (VIO) frameworks using a series of
Bayesian filters or Factor graphs [17, 30].

It is important to stress that our implementation is not a fully-functional
VIO implementation, and is only intended to empirically demonstrate the
efficacy of the $\tau$-constraint. Nevertheless, it naturally results in a
trajectory estimate. For this reason, we compare our method with the popular
state-of-the-art VIO methods VINS-Mono and ROVIO as well as the fiducial
marker based pose estimation method, AprilTag 3 [27, 20]. To provide a firm
grounding to our results, we also compare our results with millimeter accurate
ground truth from a Vicon motion capture system.

We formally define our problem statement next along with a summary of our
contributions.

###  1.1 Problem Definition and Contributions

We solve the mathematical problem of estimating distance to a scene patch,
using time-to-contact, and acceleration measurements in closed form.

To test our result, we implement a method of estimating time-to-contact from
the change of a patch’s area in the image plane and couple it with
acceleration measurements through the novel $\tau$-constraint.

A list of our contributions are given below:

  * •

A closed form solution for estimating the 3D position of the camera given a
short history of time-to-contact and acceleration.

  * •

A computationally efficient method of 3D position estimation utilizing the
$\tau$-constraint, which demonstrates viability with inexpensive sensors.

  * •

Comparisons against the popular VIO methods, VINS-Mono and ROVIO, as well as
the fiducial marker based posed estimation method, AprilTag 3, to show the
efficacy of the novel constraint in a real-world setting.

We believe that our constraint will open doors for a multitude of applications
using an active camera and an IMU.

##  2 Related Work

A multitude of prior works from both the computer vision and robotics
literature deal with time-to-contact and the fusion of camera and IMU
measurements to obtain relative camera pose (odometry). However, to the best
of our knowledge, none of these works provided a closed form solution for the
estimation of distance from time-to-contact and inertial measurements.

#### Time-to-contact

One of the earliest works to discuss how time-to-contact $\tau$ can be used in
control tasks is [22]. It provides a mathematical analysis on how time-to-
contact, $\tau$, could be used by humans driving a vehicle to judge when to
start braking and how to control ongoing braking. The idea was further
developed to create a “General Tau Theory” to understand how human’s might
perform control for a wide variety of tasks [23].

Because of its intuitive formulation, $\tau$ has also been the subject of many
studies in computer vision and robotics. [31] showed how to land a spacecraft
using event cameras by computing $\tau$ from the divergence of optical flow.
[32] fuses information from a depth camera and $\tau$ to compute “time-to-
impact” which can in-turn be used to dodge dynamic obstacles without prior
knowledge of scene geometry or obstacles. In robotics, most methods that
perform optical flow based control use initial height estimates either
implicitly or explicitly, as remarked in [15]. To this end, [15] proposed to
fuse control effort with flow divergence or $\tau$ with an extended Kalman
filter in order to recover depth. Another strategy exploits the instability
that manifests at certain heights when performing $\tau$ control with fixed
gain feedback to estimate depth [10]. In the context of self driving cars,
BinaryTTC [2] proposed a network to predict per-pixel $\tau$ with arbitrarily
fine quantization to achieve real-time performance. Finally, it is important
to note that $\tau$ can be efficiently computed when the scene under
consideration is planar in nature [16].

#### Visual Inertial Odometry (VIO)

One of the earliest works [26] for real-time VIO fused sparse feature tracks
with IMU measurements using a Multi-State Constraint Kalman Filter. The
pipeline was then made more robust by ROVIO [8] which used patch tracking and
enforced photometric consistency with IMU measurements using an Iterative
Extended Kalman Filter. ROVIO required robust initialization, accurate
intrinsic and extrinsics, and IMU calibration for accurate odometry
estimation. To this end, VINS-Mono [27] introduced a nonlinear optimization
based sliding window estimator with pre-integrated IMU factors and an
initialization-free formulation. Note that in the context of classical VIO,
direct methods operate by directly incorporating constraints on pixel
intensities into the optimization formulation whereas indirect methods utilize
features and enforce geometric constraints defined on these features from
multiple views. Fundamental to these methods is the underlying proof that
tracked points and inertial measurements allow pose to be recovered in closed
form [24, 25]. VIO is called loosely-coupled if both the sensors are used to
predict pose separately. On the other hand, tightly-coupled formulations
obtain pose by combining measurements from both sensors directly in an
optimization problem.

Recent deep learning based VIO methods have achieved better accuracy than
classical approaches. VINet [9] presented the first supervised method to
estimate odometry using a CNN + LSTM architecture by combining visual and
inertial data. This was later improved by DeepVIO [13] using supervision from
a stereo camera for better robustness. The limiting factor of these approaches
is a lack of speed and generalization across various compute platforms [28].
As will be shown, the $\tau$-constraint can be formulated as a loss, which may
be interesting for future work in deep VIO.

Our classical implementation is based on a tightly coupled optimization of the
$\tau$-constraint with respect to acceleration and time-to-contact. Time-to-
contact is determined from a directly formulated patch tracker. The results
are combined with a linear observer to limit the effects of noise.

##  3 Derivation of the $\tau$-constraint

Figure 2: Definition of camera coordinates for a scene point $\mathbf{X}$
being imaged at image point $\mathbf{x}$. The moving camera “perceives” the
movement of the scene point as optical flow $\mathbf{u}$. Figure 3: The patch
being tracked over time changes in shape and location and is modelled by an
affine transformation. Time is color coded as yellow to red. Figure 4: System
overview of our method to estimate camera position using the
$\tau$-constraint.

This work incorporates two concepts: time-to-contact and visual-inertial
fusion. We mathematically derive a formulation for the relation between $\tau$
at a scene point, the depth at that scene point, and acceleration of the
camera, which we call the $\tau$-constraint. This allows us to estimate metric
distance and relative camera pose an order of magnitude faster than standard
VIO methods by ignoring pixels outside of a local neighborhood of the scene
point.

###  3.1 Organization

Below, we derive the $\tau$-constraint starting from an affine homography of a
scene patch. We show how the affine homography determines the planar flow and
then rearrange the planar flow’s parameters to find time-to-contact, and its
generalization as frequency-of-contact. Once these preliminaries are in place,
the $\tau$-constraint is stated and proven.

It is important to note that time-to-contact can be estimated with methods
other than planar flow and thus the $\tau$-constraint stands on its own.
However, building up to the statement of the $\tau$-constraint through planar
flow is more illuminating.

###  3.2 Estimating Time-to-Contact with Affine Flow

As illustrated in Fig. 2, the optical flow $\mathbf{u}(t)$ at an image
location $\mathbf{x}=\begin{bmatrix}x&y&1\end{bmatrix}^{T}$ is given as

$\mathbf{u}_{\mathbf{x}}=K\frac{\dot{\mathbf{X}}_{\mathbf{x}}}{Z_{\mathbf{x}}}-K\frac{\dot{Z}_{\mathbf{x}}}{Z_{\mathbf{x}}}\mathbf{x},$
(1)

where $\dot{\mathbf{X}}_{\mathbf{x}}$ is the velocity of the scene point
corresponding to pixel $\mathbf{x}$ in the camera coordinate frame.

Without loss of generality, let us assume the camera is calibrated so that the
focal length is of unit length and the projection matrix is
$K=\mathbb{I}_{3\times 3}$.

Under rigid motion, $\dot{\mathbf{X}}_{\mathbf{x}}$ is related to the camera
translational velocity $\mathbf{V}$ and rotational velocity $\mathbf{\Omega}$
as
$\dot{\mathbf{X}}_{\mathbf{x}}=-\mathbf{\Omega}\times\mathbf{X}_{\mathbf{x}}-\mathbf{V}$.
Due to de-rotation from IMU measurements (rotation compensated image frames)
we proceed with the assumption that $\mathbf{\Omega}=0$ and without loss of
generality refer to $\dot{\mathbf{X}}_{\mathbf{x}}$ as $\dot{\mathbf{X}}$.

Let us assume that the scene patch we are looking at is planar. A simple
formulation to model the change of patch area over time is the affine
homography transformation. Because the affine homography is linear, its
derivative defines an affine flow field. Thus, we can model the optical flow
of the rotation compensated image patch corresponding to the planar scene
patch as

$\mathbf{u}_{\mathbf{x}}=A\mathbf{x}=\begin{bmatrix}a_{1}&a_{2}&a_{3}\\\
a_{4}&a_{5}&a_{6}\\\ 0&0&0\end{bmatrix}\mathbf{x},$ (2)

where $A(t)$ is a time-varying affine transformation of the homogeneous image
location $\mathbf{x}$.

Let us derive A in terms of physical quantities. Let the scene points
$\mathbf{X}_{\mathbf{x}}$ lie on the plane
$1/Z_{\mathbf{x}}=\begin{bmatrix}{n_{x},n_{y},n_{z}}\end{bmatrix}\begin{bmatrix}X_{\mathbf{x}}/Z_{\mathbf{x}},Y_{\mathbf{x}}/Z_{\mathbf{x}},1\end{bmatrix}^{T}=\mathbf{n}^{T}\mathbf{x}$.
By substituting this expression for $1/{Z_{\mathbf{x}}}$ into equation (2), we
obtain:

$\begin{split}\mathbf{u}_{\mathbf{x}}&=\left[\dot{\mathbf{X}}-\dot{Z}\mathbf{x}\right]\mathbf{n}^{T}\mathbf{x}\\\
&=\dot{\mathbf{X}}\mathbf{n}^{T}\mathbf{x}-\dot{Z}n_{z}\mathbf{x}-\underbrace{\dot{Z}\mathbf{x}\left(n_{x}x+n_{y}y\right)}_{\gamma}\\\
&=\underbrace{\begin{bmatrix}\dot{X}n_{x}-\dot{Z}n_{z}&\dot{X}n_{y}&\dot{X}n_{z}\\\
\dot{Y}n_{x}&\dot{Y}n_{y}-\dot{Z}n_{z}&\dot{Y}n_{z}\\\ 0&0&0\\\
\end{bmatrix}}_{A(t)\coloneqq}\mathbf{x}-\begin{bmatrix}\gamma_{x}\\\
\gamma_{y}\\\ 0\end{bmatrix}\end{split}$ (3)

The assumption of an affine homography, as used by the tracker, implies an
affine flow-field. Thus, $\gamma$ will always be measured as zero, i.e., the
quadratic terms in the flow field are not known.

By rearranging the components of $A$ to recover
$\dot{\mathbf{X}}\mathbf{n}^{T}$ we can then write
$\dot{\mathbf{X}}/Z_{\mathbf{x}}$ in terms of the affine flow parameters $A$
and an image location $\mathbf{x}$. Let
$\eta\coloneqq((a_{4}a_{3}/a_{6})-a_{1})$, then

$\dfrac{\dot{\mathbf{X}}}{Z_{\mathbf{x}}}=\begin{bmatrix}a_{4}a_{3}/a_{6}&a_{2}&a_{3}\\\
a_{4}&a_{2}a_{6}/a_{3}&a_{6}\\\
\eta(a_{4}/a_{6})&\eta(a_{2}/a_{3})&\eta\end{bmatrix}\mathbf{x}.$ (4)

The third element of $\dot{\mathbf{X}}/Z_{\mathbf{x}}$ is the inverse of time-
to-contact. Thus, by tracking a patch by an affine homography, we can estimate
the affine flow parameters, which determine time-to-contact. Because an affine
homography defines the location of every point on the planar scene patch in
image space at all times, we assume $\dot{\mathbf{X}}/Z$ is defined in
reference to a particular scene point in the tracked patch for the rest of the
paper.

###  3.3 Frequency-of-Contact as a Generalization of Time-to-Contact

Now, we define the frequency-of-contact as

$\mathbf{F}(t)\coloneqq\dfrac{\dot{\mathbf{X}}(t)}{Z(t)}$ (5)

This is the generalization of time-to-contact to all three dimensions because,
time-to-contact $\tau$, is classically defined only for camera $Z$ direction
as $\tau(t)\coloneqq-1/F_{z}(t)$.

Frequency-of-contact is more often well-defined than time-to-contact because
it avoids division by zero when the camera is stationary. $\mathbf{F}(t)$ is
ill-defined when $Z(t)=0$. However, this means that the object is on the image
plane and cannot be “seen”. For this reason, and without loss of generality,
for the rest of the paper we assume $Z$ is bounded away from zero, i.e.,
$Z(t)\geq\epsilon>0\,\,|\,\,\forall t$.

###  3.4 The $\tau$-constraint

To allow concise expression of the $\tau$-constraint we define the following
quantities and an operator that is the double integral of vectors of
functions.

$\begin{split}&\mathbf{E}(t)\coloneqq\begin{bmatrix}\displaystyle\int_{0}^{t}\begin{bmatrix}F_{X}(\lambda)\\\
F_{Y}(\lambda)\end{bmatrix}\Phi_{F_{Z}}(\lambda)d\lambda\\\
\Phi_{F_{Z}}(t)-1\end{bmatrix}-t\mathbf{F}(0)\\\
&\Phi_{F_{Z}}(t)\coloneqq\exp\left(\int_{0}^{t}F_{Z}(\lambda)d\lambda\right)\\\
&\mathbf{\Delta}\\{\mathbf{f}\\}(t)\coloneqq\int_{0}^{t}\left(\int_{0}^{\lambda}\mathbf{f}(\lambda_{2})d\lambda_{2}\right)d\lambda\end{split}$
(6)

$\mathbf{E}(t)$ equals the position change of the observer due to acceleration
divided by the scene depth $Z(0)$. Thus, the magnitude of $\mathbf{E}(t)$ is
proportional to the applied acceleration (action) and so we call it the
“action’s effect”.

######  Theorem 3.1 ($\tau$-constraint).

If $\mathbf{F}(t)$ and $\ddot{\mathbf{X}}(t)$ are known for all $t$ in a
closed interval and $Z(t)\geq\epsilon>0$, then the following constraint
between depth $Z(0)$, frequency-of-contact $\mathbf{F}$, and acceleration
$\ddot{\mathbf{X}}$ holds for each point on the scene patch.

$\begin{split}Z(0)\mathbf{E}(t)=\mathbf{\Delta}\\{\ddot{\mathbf{X}}\\}(t)\end{split}$
(7)

###### Proof.

The proof proceeds by solving a differential equation whose solution is
determined by $\mathbf{F}$ and initial conditions $\mathbf{X}(0)$. Next,
$\mathbf{X}(t)$ is expanded using the fundamental theorem of calculus so that
$\mathbf{X}$ is determined by $\mathbf{F}(0)$, $\mathbf{X}(0)$, and
$\ddot{\mathbf{X}}(0)$. Setting the two expressions for $\mathbf{X}(t)$ equal,
leads directly to the $\tau$-constraint.

The definition of frequency-of-contact Eq. 5 can be interpreted as a time-
varying linear differential equation.

$\dot{\mathbf{X}}(t)=\mathbf{F}(t)Z(t)$ (8)

The solution to the above is given by

$\begin{split}&\mathbf{X}(t)=\begin{bmatrix}X(0)\\\ Y(0)\\\
0\end{bmatrix}+\begin{bmatrix}\displaystyle\int_{0}^{t}\begin{bmatrix}F_{X}(\lambda)\\\
F_{Y}(\lambda)\end{bmatrix}\Phi_{F_{Z}}(\lambda)d\lambda\\\
\Phi_{F_{Z}}(t)\end{bmatrix}Z(0)\end{split}$ (9)

Now, from the fundamental theorem of calculus, $\mathbf{X}(t)$ is given by

$\begin{split}\mathbf{X}(t)&=\mathbf{X}(0)+t\underbrace{\mathbf{F}(0)Z(0)}_{\dot{\mathbf{X}}(0)}+\underbrace{\int_{0}^{t}\left(\int_{0}^{\lambda}\ddot{\mathbf{X}}(\lambda_{2})d\lambda_{2}\right)d\lambda}_{\mathbf{\Delta}\\{\ddot{\mathbf{X}}\\}(t)}\\\
\end{split}$ (10)

Eq. 10 satisfies the differential equation defined by $\mathbf{F}$ in Eq. 8,
given equal initial conditions. Thus, if the solutions to both systems are
unique, they are equal. The complete proof is given in the supplementary
material. An intuitive explanation is that $\mathbf{X}$ and its derivatives
define all the terms of both of the above expressions for $\mathbf{X}$. Thus,
the solution of Eq. 8 equals Eq. 10.

Setting Eq. 9 equal to Eq. 10 and canceling terms, leads directly to the
$\tau$-constraint.

∎

##  4 Fusion of Interial Measurements and Frequency-of-Contact

To tightly couple IMU measurements and frequency-of-contact from a camera in a
sliding window manner, we setup a least squares problem over the
$\tau$-constraint. The solution to this problem gives us the current depth of
a point in the scene and gravitational acceleration in a fixed frame (camera
frame at time $t=0$, commonly called start-of-service frame).

To re-iterate, as illustrated in Fig. 4, we make our system rotation invariant
by integrating the gyroscope measurements $\mathbf{\Omega}$ to get an estimate
of $R$ that rotates the measured acceleration $\mathbf{a}^{m}_{c}(t)$ back to
a fixed frame $\mathbf{a}^{m}(t)=R(t)\mathbf{a}^{m}_{c}(t)$. The IMU measures
the resultant of gravitational acceleration and linear acceleration. Thus,
$\mathbf{a}^{m}(t)=-\ddot{X}+\bm{g}$.

Note that instead of using the IMU measurements, $\ddot{\mathbf{X}}$ could
also be determined from the control effort exerted by an active observer
(sometimes called “efference copies”).

Let us suppose a history of $\mathbf{a}^{m}(t)$ and $\mathbf{F}(t)$ are
available over a time interval $[0,T]$ and define them to be zero if
$t\not\in[0,T]$. Similarly, let $\mathbf{1}_{[0,T]}$ be the indicator function
over the time window $[0,T]$. Now, without loss of generality, we consider the
problem along only the $Z$ axis. The problem can also be solved using motion
along the $X$ or $Y$ axis by replacing the subscripts of $Z$ with $X$ or $Y$.
The fusion from two sources is done through the following optimization problem
over the aforementioned time window.

$\operatorname*{argmin}_{Z(0),g_{Z}}\|{E_{Z}Z(0)+\Delta\\{a^{m}_{Z}\\}+\Delta\\{\mathbf{1}_{[0,T]}\\}g_{Z}}\|_{2}^{2}$
(11)

Equation (11) is represented by the $\tau$-constraint element in Fig. 4.

###  4.1 Efficient Computation of $Z$ Distance

Now, we can re-write Eq. 11 as

$\operatorname*{argmin}_{y}y^{T}Qy+c^{T}y$ (12)
$\begin{split}Q&\coloneqq\begin{bmatrix}\left\|E_{Z}\right\|_{2}^{2}&\langle
E_{Z},\Delta\\{\mathbf{1}_{[0,T]}\\}\rangle_{2}\\\ \langle
E_{Z},\Delta\\{\mathbf{1}_{[0,T]}\\}\rangle_{2}&\left\|\Delta\\{\mathbf{1}_{[0,T]}\\}\right\|_{2}^{2}\end{bmatrix}\\\
c&\coloneqq 2\begin{bmatrix}\langle
E_{Z},\Delta\\{a^{m}_{Z}\\}\rangle_{2}&\langle\Delta\\{\mathbf{1}_{[0,T]}\\},\Delta\\{a^{m}_{Z}\\}\rangle_{2}\end{bmatrix}^{T}\\\
y&\coloneqq\begin{bmatrix}Z(0)&g_{Z}\end{bmatrix}^{T}.\\\ \end{split}$ (13)

Which, if $Q\succ 0$, has the following unique solution:

$\begin{split}\begin{bmatrix}Z^{*}(0)\\\
g^{*}_{Z}\end{bmatrix}&=-\frac{1}{2}Q^{-1}c\\\ \end{split}$ (14)

By the positive semi-definite form of the problem given by Eq. 11, $Q\succeq
0$. Also, the problem is well-posed ($Q\succ 0$) if and only if the
acceleration changes during the considered time window. The proof is given in
the supplemental material.

Thus, given motion along an axis, scale and gravitational acceleration can be
recovered efficiently, by integrating time-to-contact and acceleration
measurements and solving a 2-dimensional linear system.

###  4.2 Filtering of Depth

Now, we extend the above derivation for trajectory estimation. Suppose
measurements of $\mathbf{F(t)}$ and $\mathbf{a}^{m}(t)$ are available over a
time window $[t-T,t],T>0$. Then, as long as the acceleration is not constant
in the time window, Eq. 14 can provide estimates of $g^{*}_{z}(t)$ and $Z(t)$.
Further, an estimate of $\dot{Z}(t)$ can be found by multiplying depth with
the frequency of contact in the Z direction, i.e., $\dot{Z}(t)=F_{Z}(t)Z(t)$.

Then, all three quantities can be fused in a Luenberger observer [14] to
obtain smooth trajectory estimates. Let $\hat{Z}$ and $\hat{\dot{Z}}$ be the
estimated quantities, then the Luenberger observer for estimating depth and
relative speed is

$\begin{bmatrix}\dot{\hat{Z}}(t)\\\ \dot{\hat{\dot{Z}}}(t)\\\
\end{bmatrix}=\begin{bmatrix}\hat{\dot{Z}}(t)\\\ a^{m}_{z}(t)-g^{*}_{z}(t)\\\
\end{bmatrix}+L\begin{bmatrix}Z(t)-\hat{Z}(t)\\\
\dot{Z}(t)-\hat{\dot{Z}}(t)\\\ \end{bmatrix}$ (15)

This classic method from the control literature has the property that the
estimation error converges to zero globally and exponentially fast [14].

If multiple estimates for $Z(t)$ are available from Eq. 14 due to motion along
multiple axis, then results are averaged. In the case that an estimate of
$Z(t)$ and $g^{*}_{z}(t)$ is not available, due to lack of acceleration change
(jerk) that is required by Eq. 14, then dead reckoning based on integration of
$\dot{Z}$ is used since time-to-contact is always available during fixation.

Finally, we recover the three dimensional trajectory by multiplying $\hat{Z}$
with the current image location of the point for which time-to-contact is
being estimated, _e.g_.

$\begin{bmatrix}\hat{X}(t)&\hat{Y}(t)\end{bmatrix}^{T}\coloneqq\begin{bmatrix}x(t)&y(t)\end{bmatrix}^{T}\hat{Z}(t)$
(16)

(a) Sequence 1

(b) Sequence 3

(c) Sequence 4

(d) Sequence 7

Figure 5: Five second sample camera trajectories for various sequences that
are used to evaluate our novel $\tau$-constraint that uses a monocular camera
and an IMU to evaluate the depth of a scene patch. Fixation is illustrated by
the highlighted scene region. Time is color coded with the hot colormap.

##  5 Experiments

We designed an experiment to evaluate the $\tau$-constraint and determine if
it is a promising method for future research.

We created ten sequences with five distinct scenes in an indoor setting. Each
scene contains a planar object to fixate on as shown in Fig. 1 and Fig. 5. For
each scene, two recordings were made, one with an AprilTag and one without.

The sequences were recorded with an Intel${}^{\text{\textregistered}}$
RealSense${}^{\text{TM}}$ D435i camera using the built-in IMU and the left
grayscale camera [19]. We used the grayscale camera because the D435i’s IMU is
hardware time-stamped to the grayscale imager. The D435i was configured to
record at 90 frames per second at 848 $\times$ 480 px. resolution. The IMU was
configured to record gyroscope measurements at 400 Hz and the acceleration at
250 Hz.

We tracked the fixated patch using gyroscope measurements for rotation
stabilization and by continually fitting an affine homography using the
ubiquitous patch tracking method from [5]. The frequency-of-contact
$\mathbf{F}$ was then recovered using Eq. 4. The affine tracker was
initialized by tracking a 100$\times$100 pixel patch sub-sampled to 4000
pixels. However, note that the patch size in pixels changes dramatically
during fixation as the tracked patch grows and shrinks in area as illustrated
by Fig. 3.

Then, the $\tau$-constraint was fused with the IMU measurements using a 2
second signal history and a 100 Hz sampling rate using linear interpolation.
If the the average power of the bias corrected acceleration along an axis was
below 2 $\text{m}/\text{s}^{2}$ the resulting depth from the $\tau$-constraint
was not used. If no observations of depth were available the depth estimate
was forward propagated in time using the frequency of contact measurements
(which are always available). The Luenberger’s gain was set to
$L=\text{diag}(2,20)$ for all the sequences. Also, it is important to note,
that these parameters were not tuned using the sequences considered in this
paper.

Our implementation is prototyped in Python 3.8 using standard scientific
Python libraries with all processing done on a CPU. Numba is used to
accelerate a few critical sections such as the affine tracker and the
computation of integrals over time [21]. Two threads on a
Intel${}^{\text{\textregistered}}$ Core${}^{\text{TM}}$ i7-6820HQ laptop
processor were used to perform all computations.

The open-source VINS-Mono and ROVIO implementations provided by their
respective authors were used for comparison [27, 8, 7]. VINS-Mono was
configured to output poses at 90 Hz (the camera frame-rate), however it
produced poses only at 80 Hz. ROVIO was configured to produce poses at 90 Hz.

As another source of comparison we used the standard AprilTag 3 [20] library
with the 36h11 tag family to detect AprilTags and their corner locations. The
corners were used to solve the Perspective-n-Point problem to recover the
tag’s location in the current frame. The integrated gyroscope measurements
were used to rotate each AprilTag pose measurement back to the fixed
orientation used by the $\tau$-constraint.

We also measured the ground truth trajectory using a Vicon motion capture
system with 8 Vantage V8 cameras at 200 Hz.

We use the Trajectory Evaluation toolbox from [33] to align all trajectories
to the Vicon ground truth. Average Trajectory Error (ATE) is calculated as
described in [33] as

$\text{ATE}(\mathbf{X},\mathbf{X}^{v})=\left(\frac{1}{N}\sum_{n=0}^{N-1}\|\mathbf{X}_{n}-\mathbf{X}^{v}_{n}\|_{2}^{2}\right)^{1/2}.$
(17)

Seq. | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
Sequence Duration (s) | 14.99 | 26.14 | 32.22 | 36.23 | 16.35 | 16.27 | 8.02 | 32.15 | 26.75 | 40.07
Distance Traveled (m) | 15.71 | 29.65 | 22.14 | 34.75 | 15.56 | 15.78 | 7.30 | 26.75 | 21.40 | 35.36
Method |  |  |  |  | ATE (cm) $\downarrow$ |  |  |  |  |
AprilTag 3 [20] | 2.88 | - | 2.66 | - | 3.77 | - | 0.66 | - | 2.55 | -
VINS-Mono [27] | 5.39 | 8.79 | 14.21 | 15.37 | - | 6.11 | 1.15 | 18.45 | 13.07 | 4.34
ROVIO [8, 7] | 7.55 | 9.89 | 11.86 | 33.23 | 29.96 | 2.84 | 0.69 | 3.93 | 16.62 | 3.79
$\tau$-constraint (ours) | 8.06 | 6.91 | 11.96 | 10.27 | 16.80 | 7.21 | 10.70 | 4.34 | 2.38 | 3.24
$\tau$-constraint (ours) (30 Hz LK output) | 7.75 | 6.59 | 12.43 | 15.04 | 17.40 | 8.18 | 5.91 | 4.75 | 3.63 | 3.49
$\tau$-constraint (ours) (15 Hz LK output) | 9.81 | 7.59 | 16.17 | 19.29 | 14.96 | 8.70 | 5.86 | 5.37 | 3.47 | 3.48
Table 1: Sequence duration (second), path length (meters), and each method’s
accuracy in centimeters of Average Trajectory Error (ATE).

Here $\mathbf{X}^{v}$ is the Vicon’s trajectory estimate and $\mathbf{X}$ is
the trajectory under evaluation after alignment. $N$ is the number of Vicon
measurements in the considered time window. The trajectory under evaluation
was linearly interpolated in order to evaluate the ATE.

Since the different trajectory estimation methods require different amounts of
time to initialize, after which they produce position estimates, we perform
the alignment and ATE calculation over the sequences’ time period where
measurements were available from all sources. The duration and path length for
each sequence were also calculated over this time period.

Because the accuracy of the numerical methods used to evaluate integrals is
affected by the sampling rate, in additional experiments we reduced the output
of the LK tracker from 90 Hz to 30 Hz and 15 Hz. The results illustrate the
effect of low-framerate input on the $\tau$-constraint’s accuracy.

##  6 Results

The ATE results across all the 10 sequences are given in Tab. 1 along with the
path length and duration of the sequence. The proposed $\tau$-constraint
achieves a lower ATE than VINS-Mono in six out of nine comparable sequences
and a lower ATE than ROVIO in 5 out of 10 comparable sequences. ATE averaged
over all sequences was 8.5 cm for our method, 16.9 cm for ROVIO, and 2.8 cm
for AprilTag 3. The average ATE for VINS-Mono, averaged over all sequences
except sequence 5 was 12.2 cm. The VINS-Mono result is omitted for sequence 5
because its trajectory estimate diverged.

It is no surprise that the AprilTag method for pose estimation routinely
achieved the best ATE. This is because the AprilTag system uses the known size
of the visual fiducial to estimate depth. Regardless, the $\tau$-constraint
performs comparably to the AprilTag method in some sequences. This is
particularly noticeable in Sequence 9. In Fig. 1 the instantaneous $l_{2}$
error of each method in Sequence 9 is plotted for comparison.

Lowering the sample rate going from the LK-tracker to the $\tau$-constraint
from 90 to 30 and 15 Hz had a relatively small effect on accuracy. In a few
sequences, the ATE error actually decreased. This is probably because the
bandwidth of the cameras trajectory is limited to a few hertz since the camera
was moved by hand and because decimation can have a similar effect as low-pass
filtering.

While the ATE errors are promising, they do not indicate that our method is
better than existing VIO methods. Such a claim would require the development
of a full VIO stack around the $\tau$-constraint and comparisons on existing
VIO datasets. Developing a full VIO stack around the $\tau$-constraint is far
beyond the scope of this paper. Regardless, the ATE errors are very promising
and indicate that development of the $\tau$-constraint is a promising
direction for future work to enhance the capabilities of current VIO and SLAM.
We believe this would help achieve a new state-of-the-art.

The speed of our method was much higher than VINS-Mono and ROVIO. Our Python
implementation achieved an average performance of 7.14 $\times$ realtime or
642 frames per second (fps). VINS-Mono’s C++ implementation ran at 0.26
$\times$ realtime or 23.6 fps. ROVIO’s C++ implementation ran at 1.05 $\times$
realtime or 94.5 fps.

##  7 Limitations and Future Work

The main limitation of the $\tau$-constraint is that there must be
acceleration for the constraint to be useful. However, this is a fundamental
constraint of all monocular geometric methods that are not coupled with a
position or velocity sensor. A related limitation is that our method of using
least-squares with the $\tau$-constraint introduces the additional condition
that jerk (change of acceleration) must be non-zero at some point during the
measurement window. In practice we found sinusoidal accelerations with
approximately 2 $\text{m}/\text{s}^{2}$ of power were necessary to get good
measurements. Future work should investigate this carefully from a thorough
mathematical and empirical perspective. Additionally, it may be of interest to
the robotics community to investigate actively controlling the agent in order
to use the $\tau$-constraint.

For our method to be deployable in-the-wild, one would need to fixate (track)
on multiple patches and actively switch between them as necessary to obtain
robust pose for fast movement when one or more patches enter and leave the
field of view. Such a method would result in a fully functional VIO system
built upon the $\tau$-constraint. Based on our preliminary results, we believe
this to be a promising avenue for future research.

##  8 Conclusion

A new mathematical constraint between time-to-contact ($\tau$), acceleration,
and depth, called the $\tau$-constraint, was presented. To validate the
practical efficacy of the constraint, a series of experiments in a trajectory
estimation task were conducted using inexpensive and commercially available
sensors. The results were particularly promising. Our method ran 6.8$\times$
to 27$\times$ faster than two state-of-the-art VIO methods that would
typically be used for the same task, while achieving 30% to 50% less absolute
trajectory error. This speedup is largely due to using only a small portion
(initially 2.5%) of the image. These results indicate that further development
of the $\tau$-constraint, in theory and practice, is a promising direction for
VIO, VI-SLAM, active perception, and robotics.

## References

  * [1] John Aloimonos, Isaac Weiss, and Amit Bandyopadhyay.  Active vision.  International Journal of Computer Vision, 1(4):333–356, 1988.
  * [2] Abhishek Badki, Orazio Gallo, Jan Kautz, and Pradeep Sen.  Binary TTC: A temporal geofence for autonomous navigation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12946–12955, 2021.
  * [3] Ruzena Bajcsy.  Active perception.  Proceedings of the IEEE, 76(8):966–1005, 1988.
  * [4] Ruzena Bajcsy, Yiannis Aloimonos, and John K. Tsotsos.  Revisiting active perception.  Autonomous Robots, 42:177–196, 2018.
  * [5] Simon Baker and Iain Matthews.  Lucas-Kanade 20 years on: A unifying framework.  International Journal of Computer Vision, 56(3):221–255, 2004.
  * [6] Dana H. Ballard.  Animate vision.  Artificial Intelligence, 48(1):57–86, 1991.
  * [7] Michael Bloesch, Michael Burri, Sammy Omari, Marco Hutter, and Roland Siegwart.  Iterated extended Kalman filter based visual-inertial odometry using direct photometric feedback.  The International Journal of Robotics Research, 36(10):1053–1072, 2017.
  * [8] Michael Bloesch, Sammy Omari, Marco Hutter, and Roland Siegwart.  Robust visual inertial odometry using a direct EKF-based approach.  In 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 298–304, 2015.
  * [9] Ronald Clark, Sen Wang, Hongkai Wen, Andrew Markham, and Niki Trigoni.  VINet: Visual-inertial odometry as a sequence-to-sequence learning problem.  In Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, volume 31, pages 3995–4001, 2017.
  * [10] Guido C H E de Croon.  Monocular distance estimation with optical flow maneuvers and efference copies: a stability-based strategy.  Bioinspiration & Biomimetics, 11(1):016004, 2016.
  * [11] Cornelia Fermüller.  Navigational preliminaries.  In Yiannis Aloimonos, editor, Active Perception, chapter 3, pages 103–150. Lawrence Erlbaum Assoc., Hillsdale, NJ, 1993.
  * [12] Cornelia Fermüller and Yiannis Aloimonos.  Tracking facilitates 3-D motion estimation.  Biological Cybernetics, 67(3):259–268, 1992.
  * [13] Liming Han, Yimin Lin, Guoguang Du, and Shiguo Lian.  DeepVIO: Self-supervised deep learning of monocular visual inertial odometry using 3D geometric constraints.  In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 6906–6913. IEEE, 2019.
  * [14] João P. Hespanha.  Linear Systems Theory.  Princeton University Press, 2nd edition, 2018.
  * [15] Hann Woei Ho, Guido CHE de Croon, and Qiping Chu.  Distance and velocity estimation using optical flow from a monocular camera.  International Journal of Micro Air Vehicles, 9(3):198–208, 2017\.
  * [16] Berthold K.P. Horn, Yajun Fang, and Ichiro Masaki.  Time to contact relative to a planar surface.  In 2007 IEEE Intelligent Vehicles Symposium, pages 68–74, 2007\.
  * [17] Vadim Indelman, Stephen Williams, Michael Kaess, and Frank Dellaert.  Information fusion in navigation systems via factor graph based incremental smoothing.  Robotics and Autonomous Systems, 61(8):721–738, 2013.
  * [18] Nitin Jagannatha Sanket.  Active Vision Based Embodied-AI Design for Nano-UAV Autonomy.  PhD thesis, University of Maryland, College Park, 2021.
  * [19] Leonid Keselman, John Iselin Woodfill, Anders Grunnet-Jepsen, and Achintya Bhowmik.  Intel(R) RealSense(TM) stereoscopic depth cameras.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pages 1267–1276, 2017.
  * [20] Maximilian Krogius, Acshi Haggenmiller, and Edwin Olson.  Flexible layouts for fiducial tags.  In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 1898–1903, 2019.
  * [21] Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert.  Numba: A LLVM-based Python JIT compiler.  In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, pages 7:1–7:6, 2015.
  * [22] David N. Lee.  A theory of visual control of braking based on information about time-to-collision.  Perception, 5(4):437–459, 1976.
  * [23] David N Lee, Reinoud J Bootsma, Mike Land, David Regan, and Rob Gray.  Lee’s 1976 paper.  Perception, 38(6):837–858, 2009.
  * [24] Agostino Martinelli.  Vision and IMU data fusion: Closed-form solutions for attitude, speed, absolute scale, and bias determination.  IEEE Transactions on Robotics, 28(1):44–60, 2012.
  * [25] Agostino Martinelli.  Closed-form solution of visual-inertial structure from motion.  International Journal of Computer Vision, 106(2):138–152, 2014\.
  * [26] Anastasios I. Mourikis and Stergios I. Roumeliotis.  A multi-state constraint Kalman filter for vision-aided inertial navigation.  In Proceedings of the 2007 IEEE International Conference on Robotics and Automation, pages 3565–3572, 2007.
  * [27] Tong Qin, Peiliang Li, and Shaojie Shen.  VINS-Mono: A robust and versatile monocular visual-inertial state estimator.  IEEE Transactions on Robotics, 34(4):1004–1020, 2018.
  * [28] Nitin J Sanket, Chahat Deep Singh, Cornelia Fermüller, and Yiannis Aloimonos.  PRGFlow: Unified SWAP-aware deep global optical flow for aerial robot navigation.  Electronics Letters, 57(16):614–617, 2021.
  * [29] Nitin J. Sanket, Chahat Deep Singh, Kanishka Ganguly, Cornelia Fermüller, and Yiannis Aloimonos.  GapFlyt: Active vision based minimalist structure-less gap detection for quadrotor flight.  IEEE Robotics and Automation Letters, 3(4):2799–2806, 2018.
  * [30] Fendy Santoso, Matthew A. Garratt, and Sreenatha G. Anavatti.  Visual–inertial navigation systems for aerial robotics: Sensor fusion and technology.  IEEE Transactions on Automation Science and Engineering, 14(1):260–275, 2017.
  * [31] Olaf Sikorski, Dario Izzo, and Gabriele Meoni.  Event-based spacecraft landing using time-to-contact.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, pages 1941–1950, 2021.
  * [32] Celyn Walters and Simon Hadfield.  EVReflex: Dense time-to-impact prediction for event-based obstacle avoidance.  In 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 1304–1309, 2021.
  * [33] Zichao Zhang and Davide Scaramuzza.  A tutorial on quantitative trajectory evaluation for visual(-inertial) odometry.  In 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 7244–7251, 2018.


#  The precession and inclinational variation across the plane of the Milky
Way revealed by open clusters

Zhihong He (何治宏) School of Physics and Astronomy, China West Normal
University, No. 1 Shida Road, Nanchong 637002, China

###### Abstract

This article presents a study of the geometry and motion of the Galactic disk
using open clusters in the Gaia era. The findings suggest that the inclination
$\theta_{i}$ of the Galactic disk increases gradually from the inner to the
outer disk, with a shift in orientation at the Galactocentric radius of
approximately 6 $\pm$ 1 kpc. Furthermore, the study reveals that the inclined
orbits may be elliptical rather than circular, however, more observations are
needed to confirm this. An analysis of the vertical motion along the
Galactocentric radius reveals that the disk has warped with precession, and
that the line-of-nodes shifts at different radii, aligning with the results
from the classical Cepheids. Although there is uncertainty for
precession/peculiar motion in Solar orbit, after considering the uncertainty,
the study derives a median value of $\dot{\phi}_{LON}$ = 6.8 km s-1 kpc-1in
the Galaxy. This value for the derived precession in the outer disk is lower
than those in the literature due to the systematic motion in Solar orbit
($\theta_{i}$ = 0.6∘). The study also finds that the inclinational variation
of the disk is significant and can cause systematic motion, with the variation
rate $\dot{\theta}_{i}$ decreasing along the Galactic radius with a slope of
-8.9 $\mu$as yr-1 kpc-1. Moreover, the derived $\dot{\theta}_{i}$ in Solar
orbit is 59.1 $\pm$ 11.2sample $\pm$ 7.7${}_{V_{Z\odot}}$ $\mu$as yr-1, which
makes it observable for high precision astrometry.

Galaxy: stellar content - star clusters: general -Galaxy: warp

##  1 Introduction

When observing galaxies, astronomers often find that the disk is bent or
tilted (Jarrett et al. 2003) \- this is known as the warp structure. It is
believed that this phenomenon is caused by a variety of sources of
gravitational influence, such as interactions with neighboring galaxies or
dark matter halos (e.g. Hunter & Toomre 1969; Sparke & Casertano 1988; Quinn
et al. 1993; Shen & Sellwood 2006). The motion of stars and gas in a galaxy’s
disk can cause distortion in its shape by moving in different directions or
orbits. To study the warp structure of galaxies, researchers use observations
of the positions and movements of stars and gas clouds (e.g. Sancisi 1976;
Bosma 1981; Thilker et al. 2005) along with simulations and models of galactic
dynamics (e.g. Pringle 1992; Velazquez & White 1999; Robin et al. 2003; Roškar
et al. 2010). A deeper understanding of the warp structure is crucial to
comprehend the evolution of galaxies and the underlying processes that drive
their formation and transformation in the universe.

The warped structure in the outer disk of our Milky Way has been known for
some time, first detected through the observation of HI gas (Kerr 1957; Burke
1957). Studies have shown that the vertical angle increases to 3 degrees at a
Galactocentric radius of around 16 kpc (Burton 1988). Similarly, molecular CO
clouds observed in the first and second Galactic quadrant have also
demonstrated similar warp features in the outer arm (Dame & Thaddeus 2011; Sun
et al. 2015). However, the uncertainty of gas kinematic distances has made it
difficult to study the geometric structure of the Galactic warp. It was only
through the observation of classical Cepheids (Chen et al. 2019; Skowron et
al. 2019a, hereafter CCs) in recent years that a clear warp in the outer disk
has become measurable. This allowed for the direct measurement of the
Galactocentric radius (RGC) and scale height of the warp. Findings indicate
that both the gas and young star disk warp upwards in the first and/or second
quadrants, and downwards in the third and/or fourth quadrants. Moreover, the
absolute vertical amplitude reaches 0.3 to 5 kpc above the Galactic plane
beyond RGC = 10 to 30 kpc (e.g. Gum et al. 1960; López-Corredoira et al. 2002;
Nakanishi & Sofue 2003; Levine et al. 2006; Voskes & Butler Burton 2006; Chen
et al. 2019; Skowron et al. 2019a, b; Lemasle et al. 2022), highlighting the
significant dimensions of the warp structure.

Furthermore, the orbital inclination of the warp can lead to a systematic
velocity in the vertical direction, particularly close to the line-of-node
(LON) where the warp surface intersects with the Galactic plane. During the
ascending semicycle, this systematic velocity often displays an upward trend
while the opposite is observed in the descending semicycle. Multiple
statistical studies have identified this trend using various stellar tracers
(López-Corredoira et al. 2014; Liu et al. 2017; Poggio et al. 2017),
especially with the improved accuracy of proper motion data from later Gaia
releases (Gaia Collaboration et al. 2018, 2021, 2022b). Additionally, the
vertical velocity of stars in low-latitude regions near the Galactic plane can
be determined by its proper motion $\mu_{b}$ and is less influenced by the
radial velocity (e.g. Poggio et al. 2018, 2020; Romero-Gómez et al. 2019; Wang
et al. 2020; Li et al. 2020; Cheng et al. 2020; Li et al. 2023; Dehnen et al.
2023).

The question of whether the Milky Way’s warped disk is precessing has sparked
debate. Poggio et al. (2020) and Cheng et al. (2020) proposed that there is
precession of around 10.9 and 13.6 km s-1 kpc-1in the outer Galactic disk,
respectively, while others argue that there is no significant precession (e.g.
Wang et al. 2020; Chrobáková & López-Corredoira 2021). The debate stems from
the fact that the precession velocity calculation depends on different tracers
and warp models. Recently, Dehnen et al. (2023) introduced a new perspective
on precession in the outer disk of the Milky Way based on classical Cepheids,
which takes into account the unsteady inclination (known as the inclinational
variation rate $\dot{\theta}_{i}$ in this context) for the first time. This
study showed a precession rate that gradually decreases from 12 (RGC = 12 kpc)
to 6 (RGC = 14 kpc) km s-1 kpc-1without any inclination variations. However,
the previously mentioned studies only indicate precession (or lack thereof) in
the Galactic outer disk (RGC > 9 kpc), without any evidence of warping in the
inner Galaxy.

Star clusters are essential in Galactic studies because they provide more
precise distance and motion estimates compared to individual stellar tracers.
However, only a few open clusters (OCs) located in the southern warp have been
identified in previous study (Cantat-Gaudin et al. 2020). This is primarily
due to the scarcity of distant OCs, making it difficult to explore the
geometry of the Milky Way warp on a larger scale. Recently, we discovered
nearly 1500 reliable open clusters through machine learning and visual
inspection based on Gaia DR3 (He et al. 2023b). This new result represents a
substantial improvement over previous studies, which identified only one-forth
of the present sample’s distant OCs located more than 4 kpc away. This immense
dataset presents an opportunity to explore the wider disk via OCs’ study that
were previously unattainable. The primary objective of this study is to
utilize this extensive sample of Gaia DR3-based OCs to gather vital geometric
and kinematic information throughout the Galactic disk. By doing so, it may
provide us with further insight into the Milky Way’s warp.

Chapter 2 introduces the selection process and the warp structure inferred
from open clusters and classical Cepheids. Chapters 3 and 4 present the
geometry and kinematic results, respectively. We analyze how orbital
inclination, precession, and inclinational variation changes at different
radii. The study presents an in-depth analysis of these findings, providing
necessary insights into the warp structure traced by the sample. Chapter 5
summarizes the paper’s main findings and presents potential future research
opportunities. Finally, Appendix A presents error analysis of the OC sample,
and Appendix B show a sketch map of elliptical inclined plane.

##  2 OC sample and warp structure

Our study utilized a large sample of 2017 open clusters from Cantat-Gaudin et
al. (2020), previously used to study the structure of the Milky Way disk
within approximately 4 kpc of the Solar system in the Gaia DR2 era. We then
cross-matched the member stars of these clusters in Gaia DR3, selecting only
clusters with 20 or more member stars (a total of 1837 clusters).
Additionally, we obtained new star clusters from Gaia DR2 and EDR3 with a
sample size of over 500, including 615 in He et al. (2021, 2022), 628 in
Castro-Ginard et al. (2022), and 1656 in He et al. (2023a). We cross-matched
the member stars of these clusters in Gaia DR3 to obtain additional line-of-
sight velocity information. We removed 138 duplicate clusters (from He et al.
2021, 2022; Castro-Ginard et al. 2022) and 746 clusters with less than 20
member stars (mostly in He et al. 2023a). In total, these clusters contain
3852 open clusters, among which 92$\%$ are within 4 kpc. Our study’s sample
provides a more extensive dataset for analyzing the structure of the Milky
Way’s disk.

Recently, a machine learning algorithm was employed by He et al. (2023b) to
identify and verify 2085 star clusters/candidates, out of which 1488 were
reliable OCs (Type 1) verified through visual inspection. A total of 944 Type
1 OCs were located beyond 4 kpc, some with large extinction values (AV > 5
mag). Each star cluster underwent isochrone fitting and visual inspection,
resulting in a more in-depth analysis of distant OCs’ geometric structure. The
coordinates of the above nearby OC samples, combined with Type 1 distant OCs,
were transformed into the Galactocentric coordinate system (more details in
Appendix A). A significant number of distant star clusters demonstrate the
twisted/spiral structures of the Milky Way, as shown in Figure 1. This is the
second celestial body type, after classical Cepheids (as shown in Figure 1,
grey symbols), to directly trace the twisted geometric structure. The highest
scale height can reach over 1.5 kpc, both in the southern ($\theta$ < 0) and
northern ($\theta$ > 0) Galaxy. However, the edge-on image from the anti-
Galactocentric line shown in Figure 1 only shows the warp in the outer disk.
Are there any other tilted or warped structures in the inner regions?

Figure 1:  The distribution of open clusters in the XY plane is depicted in
both the face-on (a.) and edge-on (b.) views, with Galactocentric coordinates
($\phi,\theta$) also displayed. In the left panel, the black arrow indicates
the rotation direction of the Milky Way, and the position of the Solar system
is (-8.15, 0) kpc, according to Reid et al. (2019). Red dots represent distant
open clusters from He et al. (2023b), while blue crosses denote other Gaia-
based open clusters (Cantat-Gaudin et al. 2020; Castro-Ginard et al. 2022; He
et al. 2021, 2022, 2023a). For comparison, the gray pluses signify the CC
sample (Skowron et al. 2019a), showing a similar intuitive plot of the warp
structure (based on CCs) from Chen et al. (2019); Lemasle et al. (2022);
except for the southern warp at $\phi\sim~{}270^{\circ}$, which still lack of
OC samples there.

##  3 Disk inclination

To examine the spatial distribution of open clusters (OCs) inside and outside
the Galactic disk, we present an edge-on view of the plane at different radii
using a 2 kpc bin (Figure 2), from the Galactic anticenter direction. The
traced disc’s degree of tilt by the cluster increases as the radius gets
larger, demonstrating that the Milky Way’s disc is tilted at various radii.
Although the tilt is not very apparent at 6 kpc, a trend of tilt angle change
can be observed. Remarkably, the distribution of Cepheid variables exhibits
the same tilt pattern, and their distances obtained through the period-
luminosity relation indicate consistency with the geometric distribution
derived from Gaia’s parallax measurements. Our analysis suggests that the
tilted disc structure is not unique to outer disk, as previously thought, but
rather extends to other regions. Such findings could help to refine our
understanding of the Milky Way’s structure.

Figure 2:  The edge-on views display Gaia-based open clusters (colored dots)
and classical Cepheids (black dots) from Skowron et al. (2019a) located at
varying Galactocentric radii around RGC = 4, 6, 8, 10, 12, and 14 kpc, all
within a width range of $\pm$ 1 kpc. The gray dashed line shows the linear
regression of the inclination angle on the line of sight anti-Galactic center.
Young (< 108 yr), middle-aged (108 \- 109 yr), and old (> 108 yr) OCs are
identified by the blue, magenta, and red symbols, respectively. As the radius
increases, the thickness of the OC/CC disks gradually increase, displaying
flaring characteristics. Additionally, compared to younger open clusters, the
older ones tend to have higher scale heights, as stated in Cantat-Gaudin et
al. (2020). It is essential to note that the LON are not situated at Y = 0
kpc, and the fitted $\theta_{i}$ is visible in Figure 3-a.

To shed more light on our findings, we categorized the OC samples into three
sub-groups based on their age: young (YOCs < 100 Myr), middle-aged (MOCs:
0.1-1 Gyr), and old (OOCs > 1 Gyr) OCs. Interestingly, the inclination angles
observed by different age tracers show only slight variations, indicating that
YOCs and MOCs are highly consistent. On the other hand, OOCs usually show
varying degrees of tilt; however, the trend of the tilt angle change is
consistent with the other samples. We attribute this inconsistency among older
clusters to radial migration or incomplete sub-samples. As presented in Figure
2, there is a consistent inclination angle of about 0.6 degrees on the orbit
where the Sun is located. This value remains consistent for YOCs, MOCs, and
CCs, with only a slightly larger inclination in OOCs. However, in the case of
older clusters, the member stars in the main sequence exhibit a tendency to be
fainter compared to their younger counterparts. This difference in brightness
has the potential to result in notable variations in the observed distribution
between the older OCs and the younger ones, particularly under heavy
extinctions. Consequently, the apparent strong pattern observed in older
samples, particularly within the 3 to 5 kpc range of RGC (Figure 2), may not
accurately reflect the true inclination.

To calculate the Milky Way disc’s inclination angle ($\theta_{i}$) and
minimize uncertainties caused by bin selection and LON positions, we fitted OC
sample values at different radii using a 0.5 kpc step and obtained the LON
position from the OC kinematics (Section 4). During this process, we
discovered that the disk did not share a cohesive plane at different radii and
that the scale height of the Solar system varied in different literature.
Additionally, The lopsided warps were presented in various radius,
particularly in the outer arm regions. We believe this is due to the objects’
orbits not being circular but instead elliptical (discussed in Section 5).
Nevertheless, to simplify our calculation, we used a circular orbit
approximation and added $d_{Z_{0}}$ to eliminate orbital ellipticity/Z⊙ impact
on the scale height:

$Z(R_{GC},\phi_{GC},\phi_{LON})_{obs}=R_{GC}\times{\sin(\phi_{LON}-\phi_{GC})}\times\sin\theta_{i}+d_{Z_{0}}$
(1)

Here, $\phi_{LON}$ and $\phi_{GC}$ represent the position of LON and OC
samples in Galactocentric coordinates, and $Z_{obs}$ denotes the OCs’ scale
height. Figure 3-a displays the inclination angle of the disk at different
radii, and the error bars reflect the uncertainty in using OC samples of
varying age limits. From the inner disk to the outer disk, the tilt angle
gradually increases, with the absolute value of the tilt angle being
approximately 0∘ near RGC = 6 kpc, and the tilt direction changes. The tilt
angle experiences a slight decrease near RGC = 10 kpc, followed by a rapid
increase to 3 degrees at RGC = 14 kpc. Unfortunately, the insufficient OC
samples within RGC = 4 kpc leaves the orbital tilt in the Milky Way bulge
unknown. Nevertheless, the limited information we do have suggests a negative
trend in the inclination angle at that location.

The consistent geometric features of tilt disks displayed by different sub-
samples, as illustrated in Figure 2 and Figure 3, indicate that the samples
utilized in this study are adequately representative for investigating the
warp structures, especially for YOCs and MOCs. The degree of tilt at most
radii appears to be less reliant on the sample selection. Additionally, as
previously stated, orbital tilt is not solely reflected in geometric features;
it also displays a vertical systemic velocity in kinematics, which is
unaffected by extinction in the Galactic inner disk. As shown in Figure 4, a
notable negative value is presented in the vertical velocity in the inner
Galactic disc. This value is still maintained despite uncertainties regarding
the Sun’s peculiar motion, orbital precession, and inclinational variation in
the Galaxy. When combined with their geometric features, these compelling
pieces of evidence demonstrate that the orbit of the most inner regions of the
Milky Way is also tilted, and the direction of this tilt is opposite to that
of the outer disk.

Figure 3:  The derived parameters from Equation 1 and 2 are presented with
error bars accounting for uncertainties in OC sample age limitations. The
parameters include: (a) inclination angle $\theta_{i}$, (b) LON positions
$\phi_{LON}$, (c) precession rate $\dot{\phi}_{LON}$, (d) inclinational
variation rate $\dot{\theta}_{i}$, (e) differential scale height
d${}_{z_{0}}$, and (f) uncertainty (from V⊙) of $\dot{\phi}_{LON}$, plotted
against Galactocentric radius RGC. The typical uncertainty resulting from V⊙
for $\dot{\theta}_{i}$ is 4.3 $\mu$as yr-1.  Figure 4:  The edge-on views of
OCs in the inner (left upper panel) and outer (right upper panel) Galactic
disk, along with their corresponding vertical velocity distributions
(displayed in the lower panels) are presented. A linear regression is applied
(gray dashed line) to demonstrate the smaller inclination angles of the inner
disk for clear visualization. Even though the inner disk does not exhibit the
same significant tilt as the outer disk, a noticeable negative VZ is still
observed, consistent with the direction of the tilt-induced velocity
deviation. In the right lower panel, the black dashed line and cyan dashed
line represent the vertical velocity when precession rate $\dot{\phi}_{LON}$ =
0 and 6.8 km s-1 kpc-1, respectively. It can be observed that the
$\dot{\phi}_{LON}$ = 0 does not correspond to the observed VZ distribution, as
supported by the study conducted by Poggio et al. (2020) utilizing proper
motions of giant stars.

##  4 Galactic Motion

As mentioned above, a tilted disc can cause a vertical velocity component that
reaches its maximum near the LON and gradually decrease towards zero at the
highest point of the warp 90 degrees away from the northern/southern
solstices. Moreover, a positive LON rotation (precession, $\dot{\phi}_{LON}$)
can weaken the VZ⊙ systematical velocity caused by this tilted orbit. Figure 4
displays the measured results at RGC = 13 kpc, revealing that if only tilt and
no precession existed, the observed vertical velocity would be higher.
However, the observed velocity is consistent with a precession rate of
$\dot{\phi}_{LON}$ = 6.8 km s-1 kpc-1, which is lower than what some previous
studies (Poggio et al. 2020; Cheng et al. 2020) have reported. Nonetheless,
this outcome is independent of any warping or kinematic model, confirming the
presence of the precession of the Milky Way.

The Galactic coordinate system exhibits hints of rotation, ranging from 0.05
to 0.8 mas yr-1, based on various objects and observations (e.g. Miyamoto &
Zhu 1998; Zhu 2000; Bobylev & Bajkova 2019). Lindegren et al. (2018)
identified in Gaia data that the rotation is not more than 0.15 mas yr-1,
showing 0.1 mas yr-1 when several VLBA sources were compared (Lindegren 2020).
Such a magnitude of frame rotation can also affect the proper motion $\mu_{b}$
or vertical velocity in Solar orbit systematically. For instance, a 50 $\mu$as
yr-1 rate may cause 0 ($\phi_{LON}$ = 0∘) to $\sim\pm$ 2 km s-1($\phi_{LON}$ =
$\mp$ 90∘) vertical velocity around the Solar orbit. Additionally, it is
reasonable to consider the variation of inclination since the inclination
itself needs to increase/decrease to the current positions. Nevertheless,
previous investigations on warping disks did not consider its effect, except
for recent work Dehnen et al. (2023). In this study, the precession
$\dot{\phi}_{LON}$ and mean (since the orbit is not a rigid body)
inclinational variation $\dot{\theta}_{i}$ are described by the simple
equation:

$V_{Z}(R_{GC},\phi_{GC},\theta_{i})_{obs}=\frac{|\theta_{i}|}{\theta_{i}}\times(\frac{R_{GC}\times(\dot{\phi}_{GC}-\dot{\phi}_{LON})}{\sqrt{\frac{1+\tan^{-2}(\theta_{i})}{\cos(\phi_{LON}-\phi_{GC})}-1}})+\dot{\theta_{i}}\times
R_{GC}\times\sin(\phi_{LON}-\phi_{GC})\times\cos(\theta_{i})-V_{Z\odot}$ (2)

The most significant uncertainty when calculating the precession rate is VZ⊙
due to the 0.6-degree inclination near the Solar orbit, along with potential
uncertainties from Solar peculiar motion, precession, and inclinational
variation near the Sun. To compensate, an uncertainty of VZ⊙ = 1 $\pm$ 1 km
s-1was added. Using the OC kinematics, LON positions of all radii were fitted,
leading to a wide distribution around $\phi$ = 170∘ (Figure 3-b). The LON
positions in the outer Galactic disk ranges from $\phi$ = 145∘ (RGC = 8 kpc)
to 180∘ (RGC = 14 kpc), consistent with the LON positions from CC samples
(Chen et al. 2019), whereas the inner regions show an opposite trend.
Additionally, our study revealed that the potential systematic movement
associated with the rotation of the LON around the line-of-solstice is
negligible, with an error bar of approximately $\sim$10 $\mu$as yr-1.

In the inner disk, despite a large sub-sample size, the small absolute value
of the tilt angle (especially at RGC = 5 to 7 kpc, where the $\theta_{i}$ is
close to zero) results in a significant influence of VZ⊙ uncertainty on the
outcome. A slight change in VZ⊙ can cause a large differential
$\dot{\phi}_{LON}$ (Figure 3-f). Nevertheless, considering the uncertainty, a
precession rate greater than zero ( > 5 km s-1 kpc-1or greater) can be
observed at RGC = 4 kpc, indicating the presence of precession characteristics
in the inner disk, similar to the outer disk. Figure 3-d shows inclinational
variation across the Milky Way disk, where $\dot{\theta}_{i}$ is higher in the
inner disk and decreases with a slope of -8.9 $\mu$as yr-1 kpc-1 towards the
outer disk. The $\dot{\theta}_{i}$ in the region where the Solar system is
located is 59.1 $\pm$ 11.2sample $\pm$ 7.7${}_{V_{Z\odot}}$ $\mu$as yr-1,
consistent with hints of previous findings (Lindegren et al. 2018; Bobylev &
Bajkova 2019; Lindegren 2020). This inclinational variation across the disk
suggests that the rotation of the LON reflects not only the motion of the
Solar system but also the interior motion of the entire Galaxy.

##  5 Discussion and Conclusion

The findings presented in Figure 3-e reveal that the $d_{Z_{0}}$ values are
not zero, suggesting that the median height of the tilted Milky Way disk is
not on the same plane. This is likely an indication that the orbital motion of
the disk in the Milky Way is a tilted ellipse plane instead of a circular
orbit. Upon closer inspection, it was observed that $d_{Z_{0}}$ steadily
increases from RGC = 4 to 9 kpc, which could be attributed to the varying
position between the two foci at different orbital radii (as depicted in the
Sketch map of Appendix B). However, relying solely on ellipse orbit, the
observed trend would suggest that eccentricity of the orbits in the inner disk
are greater than 0.2 and could even be as high as 0.3, which is not
significant enough to serve as evidence. High eccentricity may cause the
lopsided warp, but it also give rise to systematic radial motion near the
solar orbit. Therefore, it is necessary to conduct further research to gain a
more profound understanding of this phenomenon.

This study presents observational evidence for precession and inclinational
variation spanning both the inner and outer disks for the first time. The
analysis of different-aged OC samples produced significant results, revealing
that the precession rates of the Galactic disk are lower than previously
reported in the literature. This deviation is attributable to the fact that
systematic vertical velocity, which affects the local standard of rest in
Solar orbit, has been taken into account. Additionally, a global inclinational
variation spread across the Galactic disk may suggest the disk galaxy is
presently undergoing a shift in inclination. However, further observations and
simulations are necessary to examine its origin better. It should be noted
that the considerable rotation of the Solar orbit may impact high precision
astrometry, necessitating additional analysis to identify any potential
effects of precession and inclinational variation on the coordinate system.

##  6 Acknowledgements

This work has made use of data from the European Space Agency (ESA) mission
GAIA (https://www.cosmos.esa.int/gaia), processed by the GAIA Data Processing
and Analysis Consortium
(DPAC,https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the
DPAC has been provided by national institutions, in particular the
institutions participating in the GAIA Multilateral Agreement. This work is
supported by "Young Data Scientists" project of the National Astronomical Data
Center, CAS (NADC2023YDS-07); and Fundamental Research Funds of China West
Normal University (CWNU, No.21E030).

##  Appendix A Data processing and error analysis

We employed the Gaia DR3 data to derive the cluster position and velocity
components from the member stars of 5340 open clusters. To do so, we first
calculated the weighted means of the member stars’ parameters and then
adjusted the parallax by adding 0.017 mas (parallax zero point, Lindegren et
al. 2021; Fabricius et al. 2021). By considering the Solar system’s position:
[RGC⊙, Z⊙] are [8.15 kpc, 5.5 pc], and Solar peculiar motion: [U⊙, V⊙, W⊙] are
[10.6, 10.7, 7.6] km s-1, and Vϕ0 = 236 km s-1 (Reid et al. 2019), we
transformed these values into a coordinate system centered on the Galactic
center (RGC, $\phi$, Z) and a three-dimensional velocity (V${}_{R_{GC}}$, Vϕ,
VZ). It is worth noting that 1269 of the 5340 clusters did not have radial
velocities, and we estimated the median value by selecting OCs in ($l\pm$ 10∘,
$b\pm$ 5∘, $\varpi\pm$ 0.05 mas) and set the error to maximum of 50 km s-1and
velocity dispersion of the sample.

It is worth mentioning that the uncertainty of the vertical velocity component
V⊙ had the most substantial impact on our results. Given the variation in W⊙
values across different studies (e.g., 7.6, 8.6, and 6.9 km s-1from VLBA
observations, Gaia data, and others, respectively. Reid et al. 2019; Gaia
Collaboration et al. 2022a; Semczuk et al. 2023), we included the effect of
this uncertainty in our calculation (Section 4). Furthermore, we added
$d_{Z_{0}}$ to reduce the impact of the inclinational ecliptic orbit and scale
height uncertainty in the Solar system (Section 3). To fit for parameters in
Figure 3, we selected an open cluster sample with err${}_{V_{Z}}$< 1 km s-1and
err${}_{V_{\phi}}$ < 2.5 km s-1, which is equivalent to the typical intrinsic
velocity dispersion of OCs (around 1 km s-1or less, Mermilliod et al. 2009),
without placing any distance-based constraints. However, we accounted for
distance errors while estimating velocity errors, resulting in the exclusion
of some star clusters with large distance errors in the velocity selections.
Furthermore, we employed a bin size of 2 kpc while fitting the data,
effectively minimizing the effects of distance errors (distance error < 500 pc
for 93$\%$ of the clusters). Figure 5 depicts the histograms of the derived
velocity errors.

Figure 5: Statistical distribution of velocity component errors presented in
the Galactocentric coordinate system, using 1 km s-1bin size for VZ and 2.5 km
s-1bin size for Vϕ and V${}_{R_{GC}}$.

##  Appendix B Sketch map of elliptical inclined galactic plane

Figure 6 depicts the sketch map of the scale height in the elliptical inclined
galactic plane. The corresponding median scale height follows the trend of
$d_{Z_{0}}$ in RGC = 4 to 11 kpc, as seen in Figure 3-e.

Figure 6:  Sketch map of scale height in the elliptical inclined galactic
plane. The left panel displays the elliptical plane at different semi-major
axis, represented by the function: $x^{2}/a^{2}+y^{2}/b^{2}=1$, where $f_{1}$
(galactic center) and $f_{2}$ are the foci, and the eccentricity is set to
0.2. Right panel demonstrates the differential scale height of the $f_{2}$ and
semi-minor axis, with the inclination angles ($\theta_{i}$) set at 0.5 (blue),
0 (magenta), and -0.5 (black) degrees.

## References

  * Bobylev & Bajkova (2019) Bobylev, V. V., & Bajkova, A. T. 2019, Astronomy Letters, 45, 208
  * Bosma (1981) Bosma, A. 1981, AJ, 86, 1791
  * Burke (1957) Burke, B. F. 1957, AJ, 62, 90
  * Burton (1988) Burton, W. B. 1988, in Galactic and Extragalactic Radio Astronomy, ed. K. I. Kellermann & G. L. Verschuur, 295–358
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1
  * Castro-Ginard et al. (2022) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2022, A&A, 661, A118
  * Chen et al. (2019) Chen, X., Wang, S., Deng, L., et al. 2019, Nature Astronomy, 3, 320
  * Cheng et al. (2020) Cheng, X., Anguiano, B., Majewski, S. R., et al. 2020, ApJ, 905, 49
  * Chrobáková & López-Corredoira (2021) Chrobáková, Ž., & López-Corredoira, M. 2021, ApJ, 912, 130
  * Dame & Thaddeus (2011) Dame, T. M., & Thaddeus, P. 2011, ApJ, 734, L24
  * Dehnen et al. (2023) Dehnen, W., Semczuk, M., & Schönrich, R. 2023, MNRAS, arXiv:2305.09343
  * Fabricius et al. (2021) Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2022a) Gaia Collaboration, Drimmel, R., Romero-Gomez, M., et al. 2022a, arXiv e-prints, arXiv:2206.06207
  * Gaia Collaboration et al. (2022b) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022b, arXiv e-prints, arXiv:2208.00211
  * Gum et al. (1960) Gum, C. S., Kerr, F. J., & Westerhout, G. 1960, MNRAS, 121, 132
  * He et al. (2023a) He, Z., Liu, X., Luo, Y., Wang, K., & Jiang, Q. 2023a, ApJS, 264, 8
  * He et al. (2023b) He, Z., Luo, Y., Wang, K., et al. 2023b, arXiv e-prints, arXiv:2305.10269
  * He et al. (2022) He, Z., Li, C., Zhong, J., et al. 2022, ApJS, 260, 8
  * He et al. (2021) He, Z.-H., Xu, Y., Hao, C.-J., Wu, Z.-Y., & Li, J.-J. 2021, Research in Astronomy and Astrophysics, 21, 093
  * Hunter & Toomre (1969) Hunter, C., & Toomre, A. 1969, ApJ, 155, 747
  * Jarrett et al. (2003) Jarrett, T. H., Chester, T., Cutri, R., Schneider, S. E., & Huchra, J. P. 2003, AJ, 125, 525
  * Kerr (1957) Kerr, F. J. 1957, AJ, 62, 93
  * Lemasle et al. (2022) Lemasle, B., Lala, H. N., Kovtyukh, V., et al. 2022, A&A, 668, A40
  * Levine et al. (2006) Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881
  * Li et al. (2023) Li, X., Wang, H.-F., Luo, Y.-P., et al. 2023, ApJ, 943, 88
  * Li et al. (2020) Li, X. Y., Huang, Y., Chen, B. Q., et al. 2020, ApJ, 901, 56
  * Lindegren (2020) Lindegren, L. 2020, A&A, 633, A1
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4
  * Liu et al. (2017) Liu, C., Tian, H.-J., & Wan, J.-C. 2017, in Formation and Evolution of Galaxy Outskirts, ed. A. Gil de Paz, J. H. Knapen, & J. C. Lee, Vol. 321, 6–9
  * López-Corredoira et al. (2014) López-Corredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A&A, 572, A101
  * López-Corredoira et al. (2002) López-Corredoira, M., Cabrera-Lavers, A., Garzón, F., & Hammersley, P. L. 2002, A&A, 394, 883
  * Mermilliod et al. (2009) Mermilliod, J. C., Mayor, M., & Udry, S. 2009, A&A, 498, 949
  * Miyamoto & Zhu (1998) Miyamoto, M., & Zhu, Z. 1998, AJ, 115, 1483
  * Nakanishi & Sofue (2003) Nakanishi, H., & Sofue, Y. 2003, PASJ, 55, 191
  * Poggio et al. (2020) Poggio, E., Drimmel, R., Andrae, R., et al. 2020, Nature Astronomy, 4, 590
  * Poggio et al. (2017) Poggio, E., Drimmel, R., Smart, R. L., Spagna, A., & Lattanzi, M. G. 2017, A&A, 601, A115
  * Poggio et al. (2018) Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21
  * Pringle (1992) Pringle, J. E. 1992, MNRAS, 258, 811
  * Quinn et al. (1993) Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523
  * Romero-Gómez et al. (2019) Romero-Gómez, M., Mateu, C., Aguilar, L., Figueras, F., & Castro-Ginard, A. 2019, A&A, 627, A150
  * Roškar et al. (2010) Roškar, R., Debattista, V. P., Brooks, A. M., et al. 2010, MNRAS, 408, 783
  * Sancisi (1976) Sancisi, R. 1976, A&A, 53, 159
  * Semczuk et al. (2023) Semczuk, M., Dehnen, W., Schönrich, R., & Athanassoula, E. 2023, MNRAS, 519, 902
  * Shen & Sellwood (2006) Shen, J., & Sellwood, J. A. 2006, MNRAS, 370, 2
  * Skowron et al. (2019a) Skowron, D. M., Skowron, J., Mróz, P., et al. 2019a, Science, 365, 478
  * Skowron et al. (2019b) —. 2019b, Acta Astron., 69, 305
  * Sparke & Casertano (1988) Sparke, L. S., & Casertano, S. 1988, MNRAS, 234, 873
  * Sun et al. (2015) Sun, Y., Xu, Y., Yang, J., et al. 2015, ApJ, 798, L27
  * Thilker et al. (2005) Thilker, D. A., Bianchi, L., Boissier, S., et al. 2005, ApJ, 619, L79
  * Velazquez & White (1999) Velazquez, H., & White, S. D. M. 1999, MNRAS, 304, 254
  * Voskes & Butler Burton (2006) Voskes, T., & Butler Burton, W. 2006, arXiv e-prints, astro
  * Wang et al. (2020) Wang, H. F., López-Corredoira, M., Huang, Y., et al. 2020, ApJ, 897, 119
  * Zhu (2000) Zhu, Z. 2000, Ap&SS, 271, 353