# Gaussian Process Model for the Local Stellar Velocity Field from Gaia Data
Release 2

Patrick Nelson1 and Lawrence M. Widrow1
1Department of Physics, Engineering Physics, and Astronomy, Queen’s
University, Kingston, K7L 3N6, Canada
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We model the local stellar velocity field using position and velocity
measurements for 4M stars from the second data release of Gaia. We determine
the components of the mean or bulk velocity in $\sim 27k$ spatially-defined
bins. Our assumption is that these quantities constitute a Gaussian process
where the correlation between the bulk velocity at different locations is
described by a simple covariance function or kernel. We use a sparse Gaussian
process algorithm based on inducing points to construct a non-parametric,
smooth, and differentiable model for the underlying velocity field. We
estimate the Oort constants $A$, $B$, $C$, and $K$ and find values in
excellent agreement with previous results. Maps of the velocity field within
$2\,{\rm kpc}$ of the Sun reveal complicated substructures, which provide
clear evidence that the local disk is in a state of disequilibrium. We present
the first 3D map of the divergence of the stellar velocity field and identify
regions of the disk that may be undergoing compression and rarefaction.

###### keywords:

Galaxy: kinematics and dynamics – Galaxy: structure

††pubyear: 2022††pagerange: Gaussian Process Model for the Local Stellar
Velocity Field from Gaia Data Release 2–Gaussian Process Model for the Local
Stellar Velocity Field from Gaia Data Release 2

##  1 Introduction

A common strategy for understanding the dynamics of the Milky way is to assume
that phenomena such as the bar, spiral arms, and warp can be understood as
perturbative departures from an equilibrium state (See, for example, Binney
(2013); Sellwood (2013)). A further assumption is that stars orbit in the mean
gravitational field of gas, dark matter, and the other stars. The stellar
components of the Galaxy are then described by a phase space distribution
function (DF), $f({\bf x},\,{\bf v},t)$, which obeys the collisionless
Boltzmann equation (CBE) coupled to gravity via Poisson’s equation. The
equilibrium/perturbation split then carries over to the DF and gravitational
potential $\Psi$. For a disk galaxy such as the Milky Way, the equilibrium
state exhibits symmetry about both the its midplane and rotation axis. Any
departure from these symmetries therefore signals a departure from
equilibrium.

The second Gaia data release (Gaia DR2, Gaia Collaboration et al. (2018))
includes measurements of positions and velocities for over seven million stars
thereby vastly increasing the number of stars for which all phase space
coordinates are known. This data provides an estimate of the present-day
stellar DF near the Sun in the sense that each star can be represented as a
phase space probability distribution function $P_{i}$ that reflects
observational uncertainties in its position and velocity:

$f({\bf x},\,{\bf v},\,t)=\sum_{i}P_{i}\left({\bf x}-{\bf x}_{i},{\bf v}-{\bf
v}_{i}\right)~{}.$ (1)

In the limit of perfect observations, $P_{i}$ reduces to a six-dimensional
delta function. In general, interpretation of the DF involves the construction
of a set of observables

$O=\int d^{3}x^{\prime}d^{3}v^{\prime}f({\bf x}^{\prime},\,{\bf
v}^{\prime},\,t){\cal F}({\bf x^{\prime}},\,{\bf v^{\prime}})$ (2)

where ${\cal F}$ is designed to pick out particular features of the DF and
also account for the selection function of the survey. For example, an
estimate of the number density $n({\bf x})$ is obtained by setting ${\cal
F}={\cal W}({\bf x}-{\bf x^{\prime}})$ where ${\cal W}$ is a localized window
function:

$n({\bf x})=\int d^{3}x^{\prime}d^{3}v^{\prime}f({\bf x}^{\prime},\,{\bf
v}^{\prime},\,t){\cal W}({\bf x}-{\bf x^{\prime}})~{}.$ (3)

In this paper, we focus on the first velocity moment of the DF:

${\bf V}({\bf x})=n^{-1}\int d^{3}x^{\prime}d^{3}v^{\prime}f({\bf
x}^{\prime},\,{\bf v}^{\prime},\,t){\cal W}({\bf x}-{\bf x^{\prime}})\,{\bf
v}~{}.$ (4)

For a system in equilibrium, ${\bf V}$ is a function of Galactocentric radius
$R$ and distance from the midplane $|z|$.

Studies of the bulk velocity in pre-Gaia surveys such as the Sloan Extension
for Galactic Understanding and Exploration (SEGUE) (Yanny et al., 2009), the
Radial Velocity Experiment (RAVE) (Steinmetz et al., 2006), and Large Sky Area
Multi-Object Fiber Spectroscopic Telescope (LAMOST) (Cui et al., 2012)
revealed bulk motions than signaled a Galaxy in disequilibrium. For example,
Widrow et al. (2012) determined the bulk vertical velocity $V_{Z}$ and its
dispersion $\sigma_{z}$ as functions of $z$ for SEGUE stars and found that
they showed features asymmetric about $z=0$. Williams et al. (2013) determined
the full 3D velocity field from some 70k red clump stars in the RAVE survey
and found vertical motions that varied with Galactocentric radius $R$ as well
as evidence for radial flows. Similarly, Carlin et al. (2013) and Pearl et al.
(2017) used LAMOST spectroscopic velocities with proper motions from the PPMXL
catalog (Roeser et al., 2010) for 400k stars to map out both vertical and
radial bulk motions. Finally, Schönrich et al. (2019) and Friske & Schönrich
(2019) found a wavelike pattern in the mean vertical velocity as a function of
angular momentum about the spin axis of the Galaxy, which served as a proxy
for Galactocentric radius.

Unfortunately, a clear picture of bulk motions near the Sun from these and
other studies never emerged. The surveys each had their own complicated
selection functions, as illustrated in Figure 1 of Carlin et al. (2013).
Furthermore, systematic distance errors could masquerade as velocity-space
substructure (Schönrich, 2012; Carrillo et al., 2018). These problems were
alleviated by Katz, D. et al. (2019) who presented a series of velocity field
maps using data from Gaia DR2 and straight-foward binning procedures. For
example, they constructed face-on maps by first dividing the sample into
vertical slabs of width $\Delta z=400\,{\rm pc}$ and then computing the mean
velocity in $xy$-cells of $200\,{\rm pc}$ by $200\,{\rm pc}$. Similar binning
procedures were used to construct maps of the velocity field in the $R-z$ and
$R-\phi$ planes.

The Gaia kinematic maps are simple and intuitive though the connection with
theory, say through the CBE, requires further analysis. In particular, the
mean velocities in the Gaia maps differ from the moments that appear in the
continuity and Jeans equations since the former are averaged over cells of a
finite size that must strike a balance between improved spatial resolution and
reduced statistical fluctuations. The situation is only exacerbated when one
goes to compute gradients of the velocity field.

In this paper, we take a different approach based on Gaussian process (GP)
regression, a Bayesian, (almost) nonparametric approach to fitting data from
the field of machine learning (See, for example, Rasmussen & Williams (2006)).
GP regression has been used extensively in data science and engineering though
only recently has it been applied to problems in astrophysics. An example that
has some similarities with the problem at hand is the construction of dust and
extinction maps in the Galaxy (see, for example, Sale & Magorrian (2018)).
(For a different approach for constructing smooth maps of streaming motions,
see Khanna et al. (2022).)

Our provisional starting point is the assumption that the velocities of stars
in the Galaxy constitute, to a good approximation, a Gaussian process.
Formally, a Gaussian process is a collection of random variables such that any
subset of these variables is Gaussian. Stellar velocities would seem to be a
good candidate for a GP since the velocity distributions for the
differentstellar components (thin and thick disks, stellar halo) are well-
described by anisotropic Maxwellians (Binney & Tremaine, 2008).

The fundamental object in GP regression is the covariance matrix, which
describes correlations in the output variables for different values of the
inputs. In our case, the outputs are the components of the velocity field and
the inputs are different positions within the Galaxy. The covariance matrix is
constructed from a kernel function, which depends of a set of hyperparameters
that control the strength and scale length of correlations in the velocity
field. Thus, although the model for the velocity field itself is
nonparametric, the covariance function is parametric. The driver in GP
regression is optimization of the likelihood function, defined below, over the
space of hyperparameters. Note that GP regression is different from kernel
smoothing or kernel density estimation, where one produces a smooth model by
convolving the data with a window function. In GP regression, we use the data
to select a model from the space of nonparametric functions in the GP prior.

To evaluate the GP likelihood function, one must invert an $N$ by $N$
covariance matrix where $N$ is the number of observations. This calculation is
an $O(N^{3})$ operation that requires $O(N^{2})$ rapid-access memory (RAM)
capabilities. Clearly, it is unfeasible to apply GP regression directly to the
Gaia data, where $N$ in DR2 already exceeds $10^{6}$. In this work, we apply a
two-fold strategy for handling the large-$N$ bottleneck. First, we bin stars
into cells of size $(\Delta x,\,\Delta y,\,\Delta z)=(125,\,125,\,\,50)\,{\rm
pc}$. The factor of $2.5$ difference in $\Delta z$ is meant to account for the
difference between the radial scale length and the vertical scale height of
the disk. Note that the width of the bins in $z$ are a factor of $8$ smaller
than those used in constructing the Katz, D. et al. (2019) face-on maps. Our
binning procedure leads to roughly $27k$ cells above our threshold of $20$
stars. We then use a sparse GP algorithm that deploys $M$ inducing points to
approximate the likelihood function. This algorithm reduces the computational
complexity to $O(NM^{2})$ and the RAM requirements to $O(NM)$.

Our GP analysis leads to a single model from which other properties of the
velocity field can be derived. By contrast Katz, D. et al. (2019) use
different binning strategies to explore different facets of the velocity
field. Moreover, the Gaia maps are somewhere between a projection, where one
integrates out one spatial dimension, and a two-dimensional slice, where one
takes a narrow range in one of the dimensions to get the velocity field on a
two-dimensional surface. By contrast, our GP model can be queried to give the
inferred velocity field with uncertainties at any point in the sample region.
Furthermore, since the model is differentiable, We can use it to infer
quantities such as the Oort constants at the position of the Sun. In addition,
derivatives of ${\bf V}$ can be used to make statements about the dynamical
state of the Galaxy. This feature allows us to map out variations in Oort
functions in the vicinity of the Sun. It also allows us to construct a 3D map
of the divergence, which is proportional to the total time derivative of the
stellar number density. In so doing, we are able to identify regions near the
Sun where the number density of stars may be undergoing compression and
rarefaction.

An outline of our paper is as follows. In Section 2 we describe the sample
used in our analysis as well as our binning strategy. In Section 3 we provide
a brief introduction to Gaussian processes and outline our sparse GP
algorithm. We also summarize results of tests performed with mock data. The
results of our analysis for Gaia DR2 are discussed in Sections 4 and 5. We
present estimates for the generalized Oort constants and maps of the velocity
field as well as a three-dimensional map of the divergence of the velocity
field. We discuss possible extensions of this work in Section 6 and provide a
brief summary of our results and some concluding remarks in Section 7.

##  2 Preliminaries

###  2.1 6D phase space catalog

The Gaia DR2 Radial Velocity survey contains 6D phase-space measurements for
over 7 million stars (Gaia Collaboration et al., 2018; Katz, D. et al., 2019).
In this work, we use the gaiaRVdelpepsdelsp43 catalog, which was constructed
to correct for systematic errors in the Gaia parallaxes and uncertainties
(Schönrich et al., 2019)111See https://doi.org/10.5281/zenodo.2557803.

Following Katz, D. et al. (2019) we use both Cartesian coordinates
$(X,Y,Z,U,V,W)$ and Galactic cylindrical coordinates
$(R,\phi,Z,V_{R},V_{\phi},V_{Z})$. The Cartesian coordinate system has the
Galactic centre at the origin and the Sun on the $-X$-axis; $U,\,V$ and $W$
are the usual velocity components where positive values at the position of the
Sun indicate motion toward the Galactic centre, the direction of Galactic
rotation, and the North Galactic Pole, respectively. We use $8.27\,{\rm kpc}$
for the Sun’s distance to the Galactic centre and $20\,{\rm pc}$ for its
distance to the Galactic midplane (see Schönrich et al. (2019) and references
therein) so that $(X_{\odot},\,Y_{\odot},\,Z_{\odot})=(-8.27,\,0,\,0.02)\,{\rm
kpc}$. We further set
$\left(U_{\odot},\,V_{\odot},\,W_{\odot}\right)=\left(11.1,\,250,\,7.24\right)\,{\rm
km}\,{\rm s}^{-1}$ (Schönrich, 2012). Our cylindrical coordinate system also
has the Galactic centre at the origin, the Sun at $\phi=0^{\circ}$, and $\phi$
increasing in the direction of Galactic rotation. This system is left-handed
in the sense that increasing $\phi$ corresponds to the direction toward the
South rather than North Galactic pole. We note that at the position of the Sun
$V_{R}=-U$ and $V_{\phi}=V$; $V_{Z}=W$.

We calculate the $U,\,V,W$ velocity components and their errors from
astrometric observations using the method described in Johnson & Soderblom
(1987). We then convert to $V_{R}$, $V_{\phi}$ and $V_{Z}$. The velocity
components, but not their errors, were provided with the gaiaRVdelpepsdelsp43
catalog. Note that in that catalog, the Galactocentric cylindrical velocity
components were given as $(U_{g},\,V_{g},\,W_{g})=(-V_{R},\,V_{\phi},\,V_{z})$
and the Cartesian coordinates were given as $(x,\,y,\,z)=(-X,\,Y,\,Z)$.

We apply the following restrictions to the sample:

  * •

color: $G_{bp}-G_{rp}<1.5$

  * •

magnitude: $G<14.5$ and $G_{bp},\,G,\,G_{rp}>0$

  * •

parallax signal to noise: $p/\sigma_{p}>4$

  * •

parallax uncertainty: $\sigma_{p}<0.1\,{\rm mas}$ with $\sigma_{p}$ given by
the Gaia pipeline

  * •

excess $B-R$ flux: $1.172<E_{bprp}<1.3$

These quality cuts are recommended in Schönrich et al. (2019) to ensure
minimal systematic biases in derived kinematic quantities. They leave 4584106
stars from the original catalog of 6606247. Finally, we apply the following
kinematic cuts to remove high-velocity outliers:

  * •

proper motion: $\mu_{\alpha},\,\mu_{\delta}<400\,{\rm mas\,yr^{-1}}$

  * •

proper motion error:
$\epsilon_{\mu_{\alpha}},\,\epsilon_{\mu_{\delta}}<20\,{\rm mas}\,yr^{-1}$

  * •

Galactocentric speed: $|{\bf V}|<600\,{\rm km\,s}^{-1}$

These cuts are similar to those implemented by Williams et al. (2013) with the
modification that our velocity cut is in terms of the speed in the local frame
of rest whereas they apply a cut on the radial velocity. These cuts eliminate
only about 400 stars.

###  2.2 Binning

As discussed above, computation time and RAM requirements make it unfeasible
to apply GP regression directly to data sets with much more than $10^{4}$
entries. For this reason, we bin stars so that input data for our GP analysis
are the mean velocity components in cells. In a sense, velocity components
averaged over stars in a cell are closer to Gaussian process than the stellar
velocities themselves. Stars in the region near the Sun can come from the thin
disk, the thick disk, or the stellar halo. The velocity distributions for
these components are roughly Maxwellian. Thus, the stellar velocities are
drawn from what might be better described as a mixture of Gaussians. On the
other hand, the average velocity of some large number of stars in a cell will
be approximately Gaussian due to the central limit theorem.

We set up a Cartesian grid of cells with size $(\Delta x,\,\Delta y,\,\Delta
z)=(125,\,125,\,\,50)\,{\rm pc}$ for the region $4<X({\rm kpc})<12$,
$-4<Y({\rm kpc})<4$, and $-2<Z({\rm kpc})<2$ and keep only those cells with
more than 20 stars. Mean values and uncertainties for the $V_{R},\,V_{\phi}$,
and $V_{Z}$ components are then calculated by a standard least squares
algorithm. In the end, we are left with mean velocities for 27305 cells
representing the observations of 3972825 stars.

##  3 Gaussian Process Regression

###  3.1 Overview of Gaussian processes

We begin with a brief review of GP regression. More thorough discussions can
be found in numerous resources such as the excellent book by (Rasmussen &
Williams, 2006). Suppose we have observations of a real scalar process $f$
such that

$y_{i}=f({\bf x}_{i})+\epsilon_{i}~{}~{}~{}~{}~{}i=\\{1,\,2,\,\dots N\\}$ (5)

where ${\bf x}_{i}$ is the input vector for the $i$’th observation, $y_{i}$ is
the scalar output and $\epsilon_{i}$ is additive noise for that observation.
For the case at hand, ${\bf x}_{i}=(X_{i},\,Y_{i},\,Z_{i})$ is the position
vector of the $i$’th cell and $f$ is $V_{R}$, $V_{\phi}$, or $V_{Z}$. If $f$
is a Gaussian Process then it is completely specified by a mean function
$m({\bf x})$ and covariance function $k({\bf x},{\bf x}^{\prime})$ such that

$\displaystyle m({\bf x})$ $\displaystyle=\mathbb{E}\left[f({\bf x})\right]$
(6) $\displaystyle k({\bf x},\,{\bf x}^{\prime})$
$\displaystyle=\mathbb{E}\left[\left(f({\bf x})-m({\bf x})\right)\left(f({\bf
x}^{\prime})-m({\bf x}^{\prime})\right)\right]$ (7)

where $\mathbb{E}$ denotes expectation value. Furthermore, the probability
distribution function of ${\mathbf{f}}\equiv\left[f_{1},\,\dots,f({\bf
x}_{n}\right]$ is given by

$P(\mathbf{f}\,\lvert\,\mathbf{X})=\mathcal{N}(\mathbf{f}\,\lvert\,\boldsymbol{\mu},\mathbf{K})\,$
(8)

where $\mathbf{X}$, $\boldsymbol{\mu}$, and $\mathbf{K}$ are aggregate vectors
of the input vectors ${\bf x}_{i}$, the mean functions $m_{i}\equiv m({\bf
x}_{i})$, and the kernel functions $K_{ij}\equiv
k(\mathbf{x}_{i},\mathbf{x}_{j})$, respectively. As usual, ${\cal N}$ denotes
a multivariate Gaussian. Equation 8 constitutes a Gaussian process prior on
$f$. For simplicity, we assume $m=0$.

The goal of GP regression is to infer $\mathbf{f}_{*}$ at new inputs ${\bf
X}_{*}$ given the data $\left\\{\mathbf{X},\,\mathbf{y}\right\\}$. If the
noise $\epsilon$ is identical, independent, and Gaussian with dispersion
$\sigma_{n}$, then the joint probability for $\mathbf{y}$ and $\mathbf{f}_{*}$
is itself Gaussian

$P(\mathbf{y},\,\mathbf{f}_{*}|\mathbf{X},\,\mathbf{X}_{*})=\mathcal{N}\left(0,\,\begin{bmatrix}\mathbf{K}+\sigma_{n}^{2}{I}&\mathbf{K}_{*}\\\
\mathbf{K}_{*}^{\mathsf{T}}&\mathbf{K}_{**}\end{bmatrix}\right)\ ,$ (9)

where $\mathbf{K}\equiv K(\mathbf{X},\mathbf{X})$, $\mathbf{K}_{*}\equiv
K(\mathbf{X},\mathbf{X}_{*})$ and $\mathbf{K}_{**}\equiv
K(\mathbf{X}_{*},\mathbf{X}_{*})$. The quantity of interest is the conditional
probability

$\displaystyle P(\mathbf{f}_{*}|\mathbf{y},\,\mathbf{X},\,\mathbf{X}_{*})$
$\displaystyle=\frac{P(\mathbf{y},\,\mathbf{f}_{*}|\mathbf{X},\,\mathbf{X}_{*})}{P(\mathbf{y})}$
(10) $\displaystyle={\cal N}\left(\bar{\mathbf{f}}_{*},\,C\right)$ (11)

where

$\bar{\mathbf{f}}_{*}=\mathbf{K}_{*}^{T}\left(\mathbf{K}+\sigma^{2}\mathbf{I}\right)^{-1}\mathbf{y}$
(12)

and

$\mathbf{C}=\mathbf{K}_{**}-\mathbf{K}_{*}^{T}\left(\mathbf{K}+\sigma^{2}\mathbf{I}\right)^{-1}\mathbf{K}_{*}$
(13)

(see Rasmussen & Williams (2006), section 2.2 and appendix A).

###  3.2 Kernel function

A key ingredient of GP regression is the kernel or covariance function. Though
there is considerable flexibility in choosing a kernel, for most problems it
is constructed by taking sums and/or products of standard kernels, which
themselves are constructed from elementary functions of the input variables.

In the present analysis, we allow for different kernels for $V_{R}$,
$V_{\phi}$, and $V_{Z}$ since they are modeled as independent scalar
functions. Each of the kernels uses a three-dimensional radial basis function
(RBF) plus a term proportional to the identity matrix that accounts for
unknown noise

$k({\bf x}_{i},\,{\bf
x}_{j})=\sigma_{f}^{2}e^{-\tilde{r}_{ij}^{2}/2}+\sigma_{n}^{2}\delta_{ij}$
(14)

where $\sigma_{f}^{2}$ is the signal variance and

$\tilde{r}_{ij}\equiv\frac{(X_{i}-X_{j})^{2}}{l_{x}^{2}}+\frac{(Y_{i}-Y_{j})^{2}}{l_{y}^{2}}+\frac{(Z_{i}-Z_{j})^{2}}{l_{z}^{2}}$
(15)

is a dimensionless pseudo-distance between ${\bf x}_{i}$ and ${\bf x}_{j}$.
The RBF kernel is positive definite, differentiable, and maximumal for ${\bf
x}_{i}={\bf x}_{j}$. These are all features that lead to realistic models for
the bulk velocity field of the Galactic disk. Through trial and error, we find
that the fit for $V_{\phi}$ is significantly improved by including in $k$ the
product of two linear kernels. The linear kernel is given by

$k_{\rm lin}({\bf x},\,{\bf
x}^{\prime})=\sum_{i=x,y,z}\mu_{i}x_{i}x_{i}^{\prime}~{}.$ (16)

It expands the space of priors on $f$ to include linear functions of the
inputs ${\bf x}$. Unlike the RBF kernel, it is non-stationary in the sense
that it depends on the absolute positions of the data points rather than the
distance between pairs of data points. For $V_{\phi}$, we add

$k_{\rm
lin2}=\left(\sum_{i=x,y,z}\mu_{i}x_{i}x_{i}^{\prime}\right)\left(\sum_{j=x,y,z}\nu_{j}x_{j}x_{i}^{\prime}\right)$
(17)

to the RBF and noise kernels in Equation 14.

The kernels defined above depend on a number of free parameters commonly
referred to as hyperparameters, since they parameterize the covariance
function rather $f$ itself. The hyperparameters determine characteristics of
functions in the prior of $f$. For example, $l_{x}$, $l_{y}$, and $l_{z}$
control the correlation lengths along the three coordinate axes. Proper choice
of hyperparameters is crucial for obtaining a good model of the data. If the
length scales are too small, then the model will tend to overfit the data,
thereby attributing small-scale bumps and wiggles to $f$ rather than noise.
Conversely, if the length scales are too large, then the model will tend to
miss important features represented in the data.

In very simple problems, one can find suitable hyperparameters by trial and
error as illustrated in Chapter 2 of Rasmussen & Williams (2006). A principled
approach, and one that is suitable for more complex problems, is to maximize
the marginal likelihood function $p(\mathbf{y}|{\mathbf{X}},\theta)$ where

$\displaystyle\log(p(\mathbf{y}|{\mathbf{X}},\theta)$
$\displaystyle=-\frac{n}{2}\log(2\pi)-\frac{1}{2}\log|K+\sigma_{n}^{2}I|$
$\displaystyle-\frac{1}{2}\mathbf{y}^{T}(K+\sigma_{n}^{2})^{-1}\mathbf{y}~{}.$
(18)

Here, $\theta$ represents the hyperparameters.

###  3.3 Sparse GP regression via inducing points

In any optimization scheme, one must evaluate the marginal likelihood function
in Equation 3.2 a large number of times for different choices of the
hyperparameters $\theta$. Each likelihood call involves inversion of the
$N\times N$ matrix $K$, which is an $O(N^{3})$ operation that requires
$O(N^{2})$ RAM. This process becomes unfeasible for $N$ much greater than
$10^{4}$. Fortunately, there are a number of algorithms that allow one to
estimate the marginal likelihood using $M<N$ inputs. These algorithms,
generally referred to as sparse GP regression, reduce the computational
complexity to $O(NM^{2})$ and the RAM requirement to $O(NM)$. In this work, we
follow the method described in (Bauer et al., 2016). We denote the full
covariance matrix ($\mathbf{K}$ in Equation 8) as $\mathbf{K}_{ff}$, the
covariance matrix for the inducing points as $\mathbf{K}_{uu}$ and the cross
covariance matrix as $\mathbf{K}_{fu}$ In sparse GP regression, one replaces
$\mathbf{K}_{ff}$ in Equation 3.2 with $\mathbf{Q}_{\rm
ff}\equiv\mathbf{K}_{\rm fu}\mathbf{K}_{\rm uu}^{-1}\mathbf{K}_{\rm uf}$ and
includes an additional term given by ${\rm Tr}(\mathbf{K}_{\rm
ff}-\mathbf{Q}_{\rm ff})/2\sigma_{n}^{2}$ (Titsias, 2009; Bauer et al., 2016).
Note that optimization requires that we compute the gradient of the marginal
likelihood with respect to the hyperparameters.

We implement our GP algorithm using GPy, a Gaussian process package written in
Python by the University of Sheffield machine learning group (GPy, 2012). We
apply the sparse GP module with a heteroscedastic Gaussian likelihood, which
allows us to incorporate uncertainties in the mean velocity components for
individual cells. Optimization is carried out by applying the VarDTC inference
method (Titsias, 2009) using the GPy optimization module and the L-BFGS-B
algorithm from the software package SciPy (Jorge Nocedal, 2006; Virtanen et
al., 2020).

###  3.4 Mock data tests

Figure 1: Results of a mock data test. Top left: the GP fit of the mock data.
Top right: the confidence region of the model. Bottom left: the generating
field for the mock data. Bottom right: the difference between the mock field
and the model, scaled by the dispersion as calculated in the GP regression.
Here and throughout, the black dot indicates the position of the Sun and
velocities are given in $\rm km\,s^{-1}$. Figure 2: Histogram of the
normalized residual $\delta$ as defined in the text. The histogram includes
residuals from all cells used in fitting the mock data and is itself
normalized so that the area under the curve is unity. The black dashed line is
the normal distribution ${\cal{N}}(0,1)$. The red dashed line is the normal
distribution fit ${\cal{N}}(0.11,0.79)$.

Before turning to Gaia DR2, we test our algorithm on a mock data sample. This
sample is constructed by replacing measured velocities in the
gaiaRVdelpepsdelsp43 catalog with velocities drawn from an analytic function
that is chosen so that the resulting maps are qualitatively similar to ones
found with the real data. For brevity, we present results for a single generic
component with velocity field. The analytic function is given by

$V(X,Y,Z)=2\left(11\cos{g_{1}}-3\cos{g_{2}}-2\cos{g_{3}}-Z^{2}-10\right)$ (19)

where $g_{1}=Z+(X-8)/6$, $g_{2}=0.52(Y-2)$, and $g_{3}=X-1.3Y-11.3$. The
velocity field for $Z=0$ is shown in the upper left panel of Figure 1.

The mock data is analysed using the same procedure that will be applied to the
real data. We first determine mean values for $V$ and errors about the mean
for the $27k$ cells. We then optimize the likelihood function over the
hyperparameters using the inducing point algorithm described above with
$M=3000$ and the RBF plus noise kernel function. Armed with optimal
hyperparameters, we determine the posterior of the model velocity field
$\bar{V}$ using Equation 12. The result for the $XY$-plane is shown in the
lower left panel of Figure 1. We also determine the covariance matrix of the
posterior using Equation 13. The square root of its diagonal elements provides
an estimate of the dispersion, $\epsilon_{V}$, which is shown in the upper
right panel of Figure 1. As expected, the dispersion rises sharply toward the
edges of the region where we have data. Finally, we show the normalized
residual $\delta\equiv(V-\bar{V})/\epsilon_{V}$ in the $XY$ plane. In Figure
2, we plot a histogram of $\delta$ over the sample volume. If the agreement
between model and mock data were perfect, we would expect the histogram to be
well-approximated by the normal distribution, ${\cal N}(0,1)$. Instead, we
find that the histogram is fit by the Gaussian ${\cal N}(0.11,0.79)$. This
indicates that systematic errors are roughly a factor of 10 smaller than
statistical ones and that uncertainties are over-estimated by $\sim 20\%$.
Given that mock data was drawn from an analytic function rather than one
generated from the GP prior, these small discrepancies are not unexpected.

##  4 Results

Table 1: Optimal hyperparameters from our sparse GP analysis of Gaia data. Units are ${\rm kpc}$ for $l_{i}$, ${\rm km\,s}^{-1}$ for $\sigma_{f}$ and ${\rm km\,s}^{-1}{\rm kpc}^{-2}$ for $\mu_{i}$ and $\nu_{i}$. | $V_{R}$ | $V_{\phi}$ | $V_{Z}$
---|---|---|---
Kernel | RBF | RBF $+$ lin2 | RBF
$l_{x}$ | 1.06 | 1.13 | 1.89
$l_{y}$ | 2.72 | 1.75 | 4.18
$l_{z}$ | 0.536 | 0.242 | 0.484
$\sigma_{f}$ | 12.4 | 10.4 | 3.31
$\mu_{x},\,\nu_{x}$ |  | 0.236, 0.228 |
$\mu_{y},\,\nu_{y}$ |  | 0.207, 0.209 |
$\mu_{z},\,\nu_{z}$ |  | 3.86, 3.71 |

As with the mock data, we analyse the Gaia DR2 sample by first determining the
mean values and uncertainties for the three velocity components in our $27k$
cells. Hyperparameters from our optimization of the likelihood function are
given in Table 1. We see that $l_{y}/l_{x}\sim 1.5-2.5$ and $l_{x}/l_{z}\sim
2-5$. The hierarchy of scales is expected given that the radial length scale
of the disk is a factor of $5-10$ times larger than the scale height and that
variations in the disk tend to be stronger in the radial direction than the
azimuthal direction.

###  4.1 Velocity field in components

Figure 3: Left: Components of the velocity field $V_{R}$, $V_{\phi}$, and
$V_{z}$ in the $z=0$ plane. Right: One-sigma confidence regions of the GP
models. Figure 4: Left: Components of the velocity field in the $\phi=0$
plane. Right: One-sigma confidence regions of the GP models.

Figures 7 and 6 show model predictions for the velocity field in the $z=0$ and
$\phi=0^{\circ}$ planes. These figures are produced by querying the model at
the appropriate points using Equation 12. By and large, the maps are in good
agreement with those found in Katz, D. et al. (2019) (See their Figures 10 and
11) though they are clearly much smoother. The GP analysis tends to pick out
features on the length scales characterized by $l_{x}$, $l_{y}$, and $l_{z}$.
We also show the statistical uncertainties in the mean for the three
components, $\epsilon_{i}$, as calculated from Equation 13. In general the
uncertainties are less than $1\,{\rm km\,s}^{-1}$ within the sample volume,
but rise rapidly as one approaches the edge of the volume. Our maps can also
be compared to those in Khanna et al. (2022) who combine data from Gaia EDR3
with radial velocity measurements from a number of spectroscopic surveys. They
construct parametric models of the velocity field in heliocentric coordinates
so a direct comparison is difficult.

The upper left panels of Figures 3 and 4 point to several regions of radial
bulk flows. In particular, there is inward bulk motion just beyond the Solar
circle that is approximately independent of $\phi$ across the sample volume.
This motion is primarily in the region $|z|<500\,{\rm pc}$ and is thus likely
associated with the thin disk. On the other hand, the outward radial flow just
inside the Solar circle is mainly found at positive $\phi$ but across the full
range in $z$ of the survey. There also appears to be another outward flow at
$R>10\,{\rm kpc}$, $\phi\simeq-10^{\circ}$ and positive $z$.

The middle-left panels of Figures 3 and 4 show $V_{\phi}$. The dominant flow
here, and indeed for ${\bf V}$ in general, is the motion about the Galactic
centre. We see that contours of constant $V_{\phi}$ roughly follow contours of
constant $R$ in Figure 7. As discussed below, the average of $V_{\phi}$ over
$\phi$ at $z=0$ yields the rotation curve near the Solar circle. The main
feature in Figure 4 is a decrease of $V_{\phi}$ as one moves away from the
midplane of the Galaxy. This trend is easily explained by an increase in
asymmetric drift due to a larger contribution from dynamically warmer stars in
the thick disk. As with $V_{R}$, we see that $V_{\phi}$ is not perfectly
symmetric about the midplane. Consider, for example, the ridge in $V_{\phi}$
between $X=-6$ and $X=-8$ in Figure 7, which coincides with the peak seen in
Figure 6 at $R=7\,{\rm kpc}$. The peak and its outer slopes are shifted by a
small, but non-negligible, amount above the midplane.

The bottom left panels of Figures 3 and 4 show our results for the vertical
bulk motion, $V_{Z}$. In the midplane we see a clear trend of increasing
$V_{Z}$ as one moves across the local patch of the Galaxy in the direction of
increasing $R$ and $\phi$. The view in the $Rz$ plane shows a mix of bending
and breathing motions. The disk is moving downward inside the Solar circle and
upward outside the Solar circle and in addition it appears to be expanding
away from the midplane. The generation of bending and breathing waves through
internal disk dynamics and interactions with satellites such as the
Sagittarius dwarf have been studied in N-body simulations by numerous authors
including (Gómez et al., 2013; Chequers et al., 2018; Laporte et al., 2019;
Poggio et al., 2021; Bennett et al., 2022) and Thulasidharan et al. (2022).
Observational studies using data from various surveys in Widrow et al. (2012);
Williams et al. (2013); Carlin et al. (2013); Carrillo et al. (2018); Katz, D.
et al. (2019); Wang et al. (2019); López-Corredoira et al. (2020).

###  4.2 Velocity vectors

Figure 5: Map of the velocity field in $\phi-Z$ plane along a cylindrical
surface at $R=8.27$. Arrows represent the components $V_{\phi}$ and $V_{Z}$
where, for clarity, we subtract $220\,{\rm km\,s}^{-1}$ from $V_{\phi}$. The
background color map shows $V_{R}$ while black contour indicates the curve
where $V_{R}=0$.

To aid our understanding of the velocity field we present velocity vector maps
in three projections. Similar vector field maps were presented in Pearl et al.
(2017) using data from LAMOST and, most recently, by Fedorov et al. (2022)
using data from Gaia DR2. We begin with Figure 5, which presents the vector
field in the $\phi Z$ plane at $R=R_{0}$, that is, along a curved cylindrical
surface that includes a part of the Solar circle. To help visualize the
velocity field, we have subtracted off the vector $220$ ${\rm km}\,{\rm
s}^{-1}$$\hat{\mathbf{\phi}}$. Radial velocities are shown as a color map.
Thus, all three components of the velocity field are represented in the
figure.

As already noted in the $V_{\phi}$ panel of Figure 6, the dominant feature of
the map is the decrease in $V_{\phi}$ due to asymmetric drift as one moves
away from the midplane. Velocities in the midplane are about $10\,{\rm
km\,s}^{-1}$ higher than they are at $|Z|=500\,{\rm pc}$. Though the dominant
flow is in the azimuthal direction, we do see a clear downward motion for
$\phi<0^{\circ}$, in agreement with Figure 7.

Figure 6: Velocity vector map in the $RZ$-plane for $\phi=0^{\circ}$ with
arrows representing $V_{R}$ and $V_{Z}$. The background color map shows
$V_{\phi}$ with the black contour indicating the curve where $V_{\phi}=220$
${\rm km}\,{\rm s}^{-1}$.

In Figure 6 we show the velocity field in the $RZ$-plane at $\phi=0^{\circ}$,
which passes through the position of the Sun. We see that there are three
distinct regions defined primarily by the sign of $V_{R}$. Specifically, we
find an inward flow centered on $(R,\,z)=(8.5,-0.2){\rm kpc}$, an outward and
downward flow inside the Solar circle, and an outward and upward flow for
$R>10\,{\rm kpc}$ and $z\sim 500\,{\rm pc}$. The results are in good agreement
with Figure 10 of Pearl et al. (2017) for regions where the samples overlap.

####  4.2.1 XY vectorfield

Figure 7: Velocity vector map in midplane. Arrows represent the $V_{X}$ and
$V_{Y}$ components with the vector $230\,{\rm km\,s}^{-1}\hat{\mathbf{\phi}}$
subtracted for clarity. The background colormap shows $V_{Z}$ and the black
contour indicates the curve where $V_{Z}=0$. Arcs show the model for nearby
spiral arms from Reid et al. (2019). In order from the inner to the outer disc
they are: Scutum (cyan), Sagittarius (magenta), Local Arm (blue), and Perseus
(black).

In Figure 7 we show the velocity field in the midplane with vertical bulk
motions are shown as the background color map. The map paints an image of
several flows in the vicinity of the Sun, which can be describe by a
combination of expansion or compression and shear. For example, roughly along
the Solar circle, we find radial compression and azimuthal shear — the
dominant motion inside the Solar circle is outward and toward positive $\phi$
while the motion just beyond the Solar circle is inward and toward negative
$\phi$.

It is tempting to attribute flows in the midplane to the presence of spiral
arms. With this in mind, we overlay the spiral arm model of Reid et al.
(2019). In general, one expects to see motion toward the Galactic centre on
the leading side of a spiral arm and motion away from the Galactic centre on
the trailing side (Kawata et al., 2014; Grand et al., 2015; Hunt et al.,
2015). This pattern is indeed seen along the Local Arm in Figure 7 though the
connection between flows in our map and the other arms is more tenuous. For
example, there is little evidence for changes in the velocity field along the
Sagittarius arm. This may be because Sagittarius is a relatively minor arm of
the Galaxy and mainly a region rich in star formation but having relatively
little mass (Churchwell et al., 2009). Analysis of density enhancement
(Benjamin et al., 2005) and Galactocentric rotation (Kawata et al., 2018) show
few of the expected features for the other arms.

##  5 Velocity Gradient

In the previous section, we showed that our GP model could be queried to to
give the bulk velocity field ${\bf V}$ at any position in the sample region.
Since the GP solution is continuous and differentiable, we can also use it to
infer the gradient of the velocity field, a $3\times 3$-tensor whose
components are $\partial V_{i}/\partial X_{j}$. These components can be
computed either by replacing the kernel functions in Equations 12 and 13 with
their gradients or by a finite difference scheme. We have confirmed that the
two methods agree for the RBF kernel. In what follows we use finite
differencing. Similar results are presented in Fedorov et al. (2021).

###  5.1 Oort constants

Gradients of the velocity field played a central role in Oort’s seminal work
on the rotation of the Galaxy Oort (1927). The basic idea was to expand the
velocity field near the Sun as a Taylor series in position. Though Oort
considered only axisymmetric flows in the $X-Y$ plane, the method was extended
to general three dimensional flows by Ogrodnikoff (1932); Milne (1935);
Chandrasekhar (1942) and Ogorodnikov (1965).

To linear order in position, the velocity field near the Sun can be written as
a Taylor series

${\bf V}={\bf V}_{\odot}+\mathbf{H}\cdot{\boldsymbol{x}}+{\cal
O}({\boldsymbol{x}}^{2})~{}.$ (20)

If we restrict ourselves to the projection of the velocity field onto the
midplane of the Galaxy, then ${\mathbf{H}}$ reduces to the $2\times 2$ matrix

$\mathbf{H}=\left(\begin{array}[]{c@{\;\;}c}\partial v_{x}/\partial x&\partial
v_{y}/\partial x\\\ \partial v_{x}/\partial y&\partial v_{y}/\partial
y\end{array}\right)_{{\bf
x}=0}\equiv\left(\begin{array}[]{c@{\;\;}c}K+C&A-B\\\
A+B&K-C\end{array}\right).$ (21)

where the second equality defines the Oort constants $A$, $B$, $C$, and $K$.
These constants measure, respectively, the azimuthal shear, vorticity, radial
shear, and divergence of the velocity field in the midplane of the disk. In
Galactocentric polar coordinates $(R,\varphi)$ the Oort constants are given by

$\displaystyle 2A$ $\displaystyle=$
$\displaystyle\frac{v_{\varphi}}{R}-\frac{\partial v_{\varphi}}{\partial
R}-\frac{1}{R}\frac{\partial v_{R}}{\partial\varphi}$ (22) $\displaystyle 2B$
$\displaystyle=$ $\displaystyle-\frac{v_{\varphi}}{R}-\frac{\partial
v_{\varphi}}{\partial R}+\frac{1}{R}\frac{\partial v_{R}}{\partial\varphi}$
(23) $\displaystyle 2C$ $\displaystyle=$
$\displaystyle-\frac{v_{R}}{R}+\frac{\partial v_{R}}{\partial
R}-\frac{1}{R}\frac{\partial v_{\varphi}}{\partial\varphi}$ (24)
$\displaystyle 2K$ $\displaystyle=$
$\displaystyle\frac{v_{R}}{R}+\frac{\partial v_{R}}{\partial
R}+\frac{1}{R}\frac{\partial v_{\varphi}}{\partial\varphi}$ (25)

Chandrasekhar (1942).

The usual method for determining the Oort constants from astrometric data is
based on the observation that the proper motion of a star in the direction of
Galactic longitude, $\mu_{l}$, and the line-of-sight velocity, $v_{\rm los}$,
can be expressed in terms of sine and cosine functions of $l$. The Oort
constants appear as coefficients in this truncated Fourier series and can
therefore be determined by standard statistical methods. (See, for example,
Feast & Whitelock (1997); Olling & Dehnen (2003); Bovy (2017); Vityazev et al.
(2018); Wang et al. (2021)

Table 2: Comparison of measured values for the Oort constants from the literature and from this work in units of ${\rm km}\,{\rm s}^{-1}\,{\rm kpc}^{-1}$. For the GP model we provide the inferred values at the exact position of the Sun and average values over a spherical volume of radius $500\,{\rm pc}$. Origin | $A$ | $B$ | $C$ | $K$
---|---|---|---|---
Oort (1927) |  ~19 |  ~-24 |  |
Feast (1997) | 14.8 $\pm$ 0.8 | -12.4 $\pm$ 0.6 |  |
Olling (2003) | 15.9 $\pm$ 2 | -16.9 $\pm$ 2 | -9.8 $\pm$ 2 |
Bovy (2017) | 15.3 $\pm$ 0.4 | -11.9 $\pm$ 0.4 | -3.2 $\pm$ 0.4 | -3.3 $\pm$ 0.6
Vityazev (2018) | 16.3 $\pm$ 0.1 | -11.9 $\pm$ 0.1 | -3.0 $\pm$ 0.1 | -4.0 $\pm$ 0.2
Li (2019) | 15.1 $\pm$ 0.1 | -13.4 $\pm$ 0.1 | -2.7 $\pm$ 0.1 | -1.7 $\pm$ 0.2
Wang (2021) | 16.3 $\pm$ 0.9 | -12.0 $\pm$ 0.8 | -3.1 $\pm$ 0.5 | -1.3 $\pm$ 1.0
This work (p) | 16.2 $\pm$ 0.2 | -11.7 $\pm$ 0.2 | -3.1 $\pm$ 0.2 | -3.0 $\pm$ 0.2
This work (v) | 15.2 $\pm$ 0.8 | -12.4 $\pm$ 0.9 | -2.9 $\pm$ 0.8 | -2.3 $\pm$ 0.5

In this paper, we compute the Oort constants directly from our GP model for
the velocity field using the above equations. Our results along with a
selection of values from the literature are given in Table 2. We find
excellent agreement with recent measurements from Bovy (2017); Vityazev et al.
(2018); Li et al. (2019) and Wang et al. (2021). As noted above, the values
from the literature all determine the Oort constants by fitting $\mu_{l}$ and
$v_{\rm los}$ to a low-order Fourier series in $l$. They do, however, use data
from different surveys and with different geometric selection functions and
sample sizes. For example, both Bovy (2017) and Li et al. (2019) consider
large samples of stars within $500\,{\rm pc}$ of the Sun. Bovy (2017)
considers main sequence stars from the Tycho–Gaia Astrometric (TGAS) catalog
(Michalik, Daniel et al., 2015) while Li et al. (2019) considers all stars
from Gaia DR2. Wang et al. (2021) uses a relatively small sample of $A$-stars
from LAMOST with a similar range in distance from the Sun. Finally, Vityazev
et al. (2018) uses stars from TGAS but with the larger reach of $1.5\,{\rm
kpc}$.

Our model allows us to predict the Oort constants at a single point, namely
the position of the Sun. These values are given in the second to last line in
Table 2. To allow for a closer comparison to literature values, we also
include values for the Oort constants averaged over a spherical volume of
radius $500\,{\rm pc}$.

###  5.2 Oort functions

Figure 8: Galactic parameters as a function of radius. Top: Oort parameters
$A$ (left) and $B$ (right). Bottom left: $V_{\phi}$. Bottom right: Angular
velocity $\Omega=A-B$. The model is queried along the line passing from the
Galactic centre through the Sun. The red dashed lines show the $1\,\sigma$
uncertainties. Figure 9: Oort constants $A$, $B$, $C$, and $K$ calculated from
Equations 22 \- 25. Colormap units are in ${\rm km}\,{\rm s}^{-1}\,{\rm
kpc}^{-1}$. The spiral arms model of Reid et al. (2019) is also over-plotted
for the local divergence ($K$). In order of inner to the outer disc the arms
are: Scutum (cyan), Sagittarius (magenta), Local Arm (blue), and Perseus
(black).

In Oort’s 1927 original work, stars are assumed to follow circular orbits in
the midplane of the Galaxy. Under this assumption, $C$ and $K$ are identically
zero and $A$ and $B$ can be written in terms of circular speed and its
gradient at the position of the Sun. They therefore provide a direct probe of
the gravitational potential and hence matter distribution in the Galaxy.
Olling & Merrifield (1998) extended the idea of Oort constants to Oort
functions with

$\displaystyle A(R)$
$\displaystyle=\frac{1}{2}\left(\frac{v_{c}(R)}{R}-\frac{dv_{c}}{dR}\right)$
(26) $\displaystyle B(R)$
$\displaystyle=\frac{-1}{2}\left(\frac{v_{c}(R)}{R}+\frac{dv_{c}}{dR}\right)$
(27)

where $v_{c}$ is the circular speed. They then derived $A(R)$ and $B(R)$ from
mass models of the Milky Way and Oort constant measurements.

In Figure 8 we show $A(R)$ and $B(R)$ over the sample region. Note that our
results for $A(R)$ and $B(R)$ use mean $V_{\phi}$ rather than $v_{c}$ and
therefore include the effects of asymmetric drift. Nevertheless, they show the
same broad trends as seen in Olling & Merrifield (1998) and more recently,
Fedorov et al. (2021). $A(R)$ is a decreasing function of $R$ with a "bump"
just inside the Solar circle. $B(R)$ is an increasing function of $R$. In the
bottom panels of Figure 8 we show the rotation curve $V_{\phi}$ and angular
velocity curve $\Omega=|A-B|$, which show departures from a flat rotation
curve at the few ${\rm km\,s}^{-1}$ level.

The expressions for the Oort constants in Equations 22-25 can similarly be
extended to function of $R$ and $\varphi$ to yield maps of Oort functions in
the Galactic midplane. These are shown in Figure 9. Were the Galaxy in an
axisymmetric steady state, we would find that $C=K=0$ and $A$ and $B$ depend
only on $R$. Though $C$ and $K$ are about a factor of five smaller than $A$
and $B$ they are clearly nonzero. Moreover, all four functions show variations
in $\varphi$ though the dominant gradients are in the $R$-direction. The
prominent regions of shear (lower left panel showing $C$) and compression
(lower right panel showing $K$) provide another way of visualizing the in-
plane flows already seen in Figure 7.

###  5.3 Divergence of the velocity field

Figure 10: Divergence of the velocity field, calculated using cylindrical
coordinates. Regions with $|\nabla\cdot\mathbf{V}|>5$ ${\rm km}\,{\rm
s}^{-1}\,{\rm kpc}^{-1}$are plotted, positive in red and negative in blue.
Coloring and opacity is scaled with amplitude of the divergence.

As mentioned above, the Oort constants are constructed from $X-Y$ derivatives
of the in-plane components of the velocity field. More generally, ${\bf H}$ in
Equation 20 is a $3\times 3$ tensor. This tensor can be written as the sum of
a six-component symmetric tensor $M^{+}$ and a three-component anti-symmetric
tensor $M^{-}$

$M^{\pm}=\frac{1}{2}\left(\frac{\partial v_{i}}{\partial
x_{k}}\pm\frac{\partial v_{k}}{\partial x_{i}}\right)$ (28)

where as before, the components are evaluated at the position of the Sun
(Ogrodnikoff, 1932; Milne, 1935; Ogorodnikov, 1965; Tsvetkov & Amosov, 2019;
Fedorov et al., 2021). Note that the coefficients of $M^{\pm}$ include the
Oort constants. For example, $A=M^{+}_{12}$ and $B=M^{-}_{21}$. Recently, the
nine parameters of the velocity gradient method have been estimated using data
from Gaia DR2 (Vityazev et al., 2018; Bobylev & Bajkova, 2021).

Here, we focus on the divergence of the velocity field
$\mathbf{\nabla}\cdot{\bf V}$, which corresponds to the trace of $M^{+}$.
Recall that the zeroth moment of the collisionless Boltzmann equation yields
the continuity equation

$\frac{\partial\rho}{\partial t}+\mathbf{\nabla}\cdot\left(\rho{\bf
V}\right)=0$ (29)

where $\rho$ is the stellar density. This equation may be written as

$\frac{1}{\rho}\frac{d\rho}{dt}=-\mathbf{\nabla}\cdot{\bf V}$ (30)

where $d/dt=\partial/\partial t+{\bf V}\cdot\mathbf{V}$ is the total
derivative. Thus, a region of positive divergence indicates that either the
stellar distribution is expanding, or stars are being transported into a
region of lower density. Either way, $\mathbf{\nabla}\cdot{\bf V}$ is a clear
signal of disequilibrium in the disk.

In Figure 10 we show the divergence in three dimensions from six different
viewing angles. We find a prominent region of negative divergence
corresponding to compression in the velocity field running through the Solar
Neighborhood and roughly aligned with Galactic azimuth. This is in good
agreement with the results for $K$ in the Galactic plane as seen in the lower-
right panel of Figure 9. This region is bracketed by regions of positive
divergence, one below the midplane and inside the Solar circle and the other
above the midplane and outside the Solar circle. The amplitude of the
divergence is $5-10\mbox{${\rm km}\,{\rm s}^{-1}\,{\rm kpc}^{-1}$}$, which
implies that the time-scale for the logarithm of the density to change is of
order $100-200\,{\rm Myr}$.

##  6 Discussion

Our discussion of the divergence of the stellar velocity field illustrates how
a GP model can provide a link between theory and data. In particular, the
divergence gives us a direct indication of regions in the disk where the total
time derivative of $\ln{\rho}$ is nonzero. Clearly, if we have $\rho$ itself,
then we can estimate the convective derivative term ${\bf
V}\cdot\mathbf{\nabla}\rho$ and hence infer the partial time derivative of
$\rho$. In this way, we gain information about time-dependent phenomena in the
disk from a kinematic snapshot. Models of $\rho$ for particular tracer
populations require detailed knowledge of the selection function but are
certainly accessible given the wealth of astrometric data available.

In a similar fashion, GP models for components of the velocity dispersion
tensor would allow one to study the first moments of the CBE, namely the Jeans
equations. In contrast with the continuity equation, these equations involve
gradients of the gravitational potential. Thus, we will typically have two
unknown quantities, one involving time-derivatives of the moments and the
other involving the potential. Nevertheless, one should be able to infer
something about the dynamics of the stellar disk, and in particular,
departures from equilibrium by modelling moments of the DF via GP regression.

These considerations suggest an extension of the present work where we treat
the output as a single three-components vector rather than three independent
scalars. Doing so would allow for correlations between the different
components of the velocity field. The framework for extending GP regression to
vector outputs is already well-established (Alvarez et al., 2011).

##  7 Conclusions

In this paper, we present a GP model for the mean or bulk stellar velocity
field in the vicinity of the Sun using astrometric measurements from Gaia DR2.
The model is nonparametric in the sense that there are no prior constraints on
the functional form of the velocity field. Instead, one specifies the
functional form of the kernel function, which through a set of
hyperparameters, controls properties of the prior on the velocity field such
as its coherence length.

The main challenge in applying GP regression to Gaia data comes from the
large-$N$ requirements in both computing time and RAM. Fortunately, there is a
large effort within the machine learning community on addressing these
problems (See, for example, Hensman et al. (2013)). In this work, we have just
started to exploit methods developed in that field.

Our model provides smooth, differentiable versions of the velocity field maps
found in Katz, D. et al. (2019) and elsewhere. The Gaia maps used different
binning schemes to compute different projections of the velocity field. In our
case, we compute a single GP model (or more precisely, three independent
models for each of the cylindrical velocity field components) from which
properties of the velocity field could be derived. As a check of the model, we
confirm that the values for the Oort constants derived by differentiating the
velocity field agree extremely well with recent determinations in the
literature.

A GP model for the velocity field can thus provide a starting point for making
contact between astrometric data and dynamics via the continuity and Jeans
equations. Already, our results for the divergence of the velocity field have
provided a unique perspective on departures in the disk from an equilibrium
state.

## Acknowledgements

It is a pleasure to thank Dan Foreman-Mackey, Hamish Silverwood, and Lauren
Anderson for helpful conversations on Gaussian processes. We acknowledge the
financial support of the Natural Sciences and Engineering Research Council of
Canada. LMW is grateful to the Kavli Institute for Theoretical Physics at the
University of California, Santa Barbara for providing a stimulating
environment during a 2019 program on galactic dynamics. His research at the
KITP was supported by the National Science Foundation under Grant No. NSF
PHY-1748958.

## Data Availability

Data used in this paper is available through Zenodo
(https://doi.org/10.5281/zenodo.2557803)

## References

  * Alvarez et al. (2011) Alvarez M. A., Rosasco L., Lawrence N. D., 2011, arXiv e-prints, p. arXiv:1106.6251
  * Bauer et al. (2016) Bauer M., van der Wilk M., Rasmussen C. E., 2016, arXiv e-prints, p. arXiv:1606.04820
  * Benjamin et al. (2005) Benjamin R. A., et al., 2005, ApJ, 630, L149
  * Bennett et al. (2022) Bennett M., Bovy J., Hunt J. A. S., 2022, ApJ, 927, 131
  * Binney (2013) Binney J., 2013, New Astron. Rev., 57, 29
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Bobylev & Bajkova (2021) Bobylev V. V., Bajkova A. T., 2021, Astronomy Letters, 47, 607
  * Bovy (2017) Bovy J., 2017, MNRAS:Letters, 468, L63
  * Carlin et al. (2013) Carlin J. L., et al., 2013, ApJ, 777, L5
  * Carrillo et al. (2018) Carrillo I., et al., 2018, MNRAS, 475, 2679
  * Chandrasekhar (1942) Chandrasekhar S., 1942, Principles of stellar dynamics. Dover Pheonix Editions
  * Chequers et al. (2018) Chequers M. H., Widrow L. M., Darling K., 2018, MNRAS, 480, 4244
  * Churchwell et al. (2009) Churchwell E., et al., 2009, PASP, 121, 213
  * Cui et al. (2012) Cui X.-Q., et al., 2012, Research in Astronomy and Astrophysics, 12, 1197
  * Feast & Whitelock (1997) Feast M., Whitelock P., 1997, MNRAS, 291, 683
  * Fedorov et al. (2021) Fedorov P. N., Akhmetov V. S., Velichko A. B., Dmytrenko A. M., Denischenko S. I., 2021, MNRAS, 508, 3055–3067
  * Fedorov et al. (2022) Fedorov P. N., Akhmetov V. S., Velichko A. B., Dmytrenko A. M., Denyshchenko S. I., 2022, arXiv e-prints, p. arXiv:2205.09381
  * Friske & Schönrich (2019) Friske J. K. S., Schönrich R., 2019, MNRAS, 490, 5414
  * GPy (2012) GPy since 2012, GPy: A Gaussian process framework in python, http://github.com/SheffieldML/GPy
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gómez et al. (2013) Gómez F. A., Minchev I., O’Shea B. W., Beers T. C., Bullock J. S., Purcell C. W., 2013, MNRAS, 429, 159
  * Grand et al. (2015) Grand R. J. J., Bovy J., Kawata D., Hunt J. A. S., Famaey B., Siebert A., Monari G., Cropper M., 2015, MNRAS, 453, 1867
  * Hensman et al. (2013) Hensman J., Fusi N., Lawrence N. D., 2013, Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, pp 282 – 290
  * Hunt et al. (2015) Hunt J. A. S., Kawata D., Grand R. J. J., Minchev I., Pasetto S., Cropper M., 2015, MNRAS, 450, 2132
  * Johnson & Soderblom (1987) Johnson D. R. H., Soderblom D. R., 1987, AJ, 93, 864
  * Jorge Nocedal (2006) Jorge Nocedal S. J. W., 2006, Numerical Optimization. Springer, https://doi.org/10.1007/978-0-387-40065-5
  * Katz, D. et al. (2019) Katz, D. et al., 2019, A&A, 622, A205
  * Kawata et al. (2014) Kawata D., Hunt J. A. S., Grand R. J. J., Pasetto S., Cropper M., 2014, MNRAS, 443, 2757
  * Kawata et al. (2018) Kawata D., Baba J., Ciucă I., Cropper M., Grand R. J. J., Hunt J. A. S., Seabroke G., 2018, MNRAS, 479, L108
  * Khanna et al. (2022) Khanna S., Sharma S., Bland-Hawthorn J., Hayden M., 2022, arXiv e-prints, p. arXiv:2204.13672
  * Laporte et al. (2019) Laporte C. F. P., Minchev I., Johnston K. V., Gómez F. A., 2019, MNRAS, 485, 3134
  * Li et al. (2019) Li C., Zhao G., Yang C., 2019, ApJ, 872, 205
  * López-Corredoira et al. (2020) López-Corredoira M., Garzón F., Wang H. F., Sylos Labini F., Nagy R., Chrobáková Ž., Chang J., Villarroel B., 2020, A&A, 634, A66
  * Michalik, Daniel et al. (2015) Michalik, Daniel Lindegren, Lennart Hobbs, David 2015, A&A, 574, A115
  * Milne (1935) Milne E. A., 1935, MNRAS, 95, 560
  * Ogorodnikov (1965) Ogorodnikov K. F., 1965, Dynamics of stellar systems. Pergamon Press
  * Ogrodnikoff (1932) Ogrodnikoff K., 1932, Z. Astrophys., 4, 190
  * Olling & Dehnen (2003) Olling R. P., Dehnen W., 2003, ApJ, 599, 275–296
  * Olling & Merrifield (1998) Olling R. P., Merrifield M. R., 1998, MNRAS, 297, 943
  * Oort (1927) Oort J. H., 1927, Bull. Astron. Inst. Netherlands, 3, 275
  * Pearl et al. (2017) Pearl A. N., Newberg H. J., Carlin J. L., Smith R. F., 2017, ApJ, 847, 123
  * Poggio et al. (2021) Poggio E., Laporte C. F. P., Johnston K. V., D’Onghia E., Drimmel R., Grion Filho D., 2021, MNRAS, 508, 541
  * Rasmussen & Williams (2006) Rasmussen C. E., Williams C. K., 2006, Gaussian Processes for Machine Learning. MIT Press, http://www.gaussianprocess.org/gpml/chapters/RW.pdf
  * Reid et al. (2019) Reid M. J., et al., 2019, ApJ, 885, 131
  * Roeser et al. (2010) Roeser S., Demleitner M., Schilbach E., 2010, AJ, 139, 2440
  * Sale & Magorrian (2018) Sale S. E., Magorrian J., 2018, MNRAS, 481, 494
  * Schönrich (2012) Schönrich R., 2012, MNRAS, 427, 274
  * Schönrich et al. (2019) Schönrich R., McMillan P., Eyer L., 2019, MNRAS, 487, 3568
  * Sellwood (2013) Sellwood J. A., 2013, in Oswalt T. D., Gilmore G., eds, , Vol. 5, Galactic Structure and Stellar Populations. Springer Science, p. 923, doi:10.1007/978-94-007-5612-0_18
  * Steinmetz et al. (2006) Steinmetz M., et al., 2006, AJ, 132, 1645
  * Thulasidharan et al. (2022) Thulasidharan L., D’Onghia E., Poggio E., Drimmel R., Gallagher J. S. I., Swiggum C., Benjamin R. A., Alves J., 2022, A&A, 660, L12
  * Titsias (2009) Titsias M., 2009, Journal of Machine Learning Research - Proceedings Track, 5, 567
  * Tsvetkov & Amosov (2019) Tsvetkov A. S., Amosov F. A., 2019, Astronomy Letters, 45, 462
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vityazev et al. (2018) Vityazev V. V., Popov A. V., Tsvetkov A. S., Petrov S. D., Trofimov D. A., Kiyaev V. I., 2018, Astronomy Letters, 44, 236
  * Wang et al. (2019) Wang H.-F., et al., 2019, MNRAS, 491, 2104
  * Wang et al. (2021) Wang F., Zhang H.-W., Huang Y., Chen B.-Q., Wang H.-F., Wang C., 2021, MNRAS, 504, 199
  * Widrow et al. (2012) Widrow L. M., Gardner S., Yanny B., Dodelson S., Chen H.-Y., 2012, ApJ, 750, L41
  * Williams et al. (2013) Williams M. E. K., et al., 2013, MNRAS, 436, 101
  * Yanny et al. (2009) Yanny B., et al., 2009, AJ, 137, 4377