# [X/Fe] Marks the Spot: Mapping Chemical Azimuthal Variations in the Galactic
Disk with APOGEE

Zoe Hackshaw Department of Astronomy, The University of Texas at Austin, 2515
Speedway Boulevard Austin, TX 78712, USA<EMAIL_ADDRESS>Keith Hawkins
Department of Astronomy, The University of Texas at Austin, 2515 Speedway
Boulevard Austin, TX 78712, USA  Carrie Filion Department of Physics &
Astronomy, The Johns Hopkins University, Baltimore, MD 21218, USA  Danny Horta
Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave., New
York, NY 10010, USA  Chervin F. P. Laporte Institut de Ciències del Cosmos
(ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028
Barcelona, Spain  Christopher Carr Department of Astronomy, Columbia
University, 550 West 120th Street, New York, NY 10027, USA  Adrian M. Price-
Whelan Center for Computational Astrophysics, Flatiron Institute, 162 5th
Ave., New York, NY 10010, USA

###### Abstract

Chemical cartography of the Galactic disk provides insights to its structure
and assembly history over cosmic time. In this work, we use chemical
cartography to explore chemical gradients and azimuthal substructure in the
Milky Way disk with giant stars from APOGEE DR17. We confirm the existence of
a radial metallicity gradient in the disk of $\Delta$[Fe/H]/$\Delta$R
$\sim-0.066\ \pm 0.0004$ dex/kpc and a vertical metallicity gradient of
$\Delta$[Fe/H]/$\Delta$Z $\sim-0.164\pm 0.001$. We find azimuthal variations
($\pm 0.1$ dex) on top of the radial metallicity gradient that have been
previously established with other surveys. The APOGEE giants show strong
correlations with stellar age and the intensity of azimuthal variations in
[Fe/H]; older stellar populations show the largest deviations from the radial
metallicity gradient. Beyond iron, we show that other elements (e.g., Mg, O)
display azimuthal variations at the $\pm 0.05$ dex-level across the Galactic
disk. We illustrate that moving into the orbit-space could help constrain the
mechanisms producing these azimuthal metallicity variations. These results
suggest that the spiral arms of the Galaxy are not solely responsible for
azimuthal metallicity variations and other Galactic processes are at play.

Galactic archaeology (2178), Milky Way disk (1050), Stellar abundances (1577)

††software: astropy (Astropy Collaboration et al., 2022, 2013), astroNN (Leung
& Bovy, 2019a), gala (Price-Whelan et al., 2020), PyAstronomy (Czesla et al.,
2019), scipy (Virtanen et al., 2020), matplotlib.pyplot (Hunter, 2007), numpy
(van der Walt et al., 2011)

##  1 Introduction

The field of Galactic archaeology exists to answer long-standing questions
about the processes that drive Galactic formation and evolution (Eggen et al.,
1962; Searle & Zinn, 1978). We can use the Milky Way and its resolved Galactic
components as a laboratory to answer questions about Galactic evolution and
characterize the hierarchical formation (Davis et al., 1985) of the Milky Way.

The current investigation of the Galactic processes that drive the Milky Way’s
evolution has exploded within the last few decades due to wide-field missions
aimed at mapping the stellar content of the Milky Way such as APOGEE (Majewski
et al., 2017a) and Gaia (Gaia Collaboration et al., 2022b). Employing stars as
the key witnesses to Galactic evolution, we are able to take stellar
information (stellar parameters, kinematics, chemical abundances, etc.) and
apply them with spatial information to create information-dense maps of the
Milky Way in a process known as chemical cartography (Hayden et al., 2015).
Measuring the distribution of elements throughout the Milky Way disk can
inform us about global and secular processes that dominate over space and time
(e.g., Hawkins et al., 2015).

To postulate which Galactic phenomena are the most influential, observations
of the Galactic disk (both thick and thin) are needed to identify which
signatures and patterns prevail. One of the most prominent trends in the Milky
Way is the existence of the negative radial and vertical metallicity gradients
(an incomprehensive list includes: Mayor, 1976; Andrievsky et al., 2002;
Magrini et al., 2009; Luck & Lambert, 2011; Bergemann et al., 2014; Xiang et
al., 2015; Yan et al., 2019; Hawkins, 2022). These gradients could provide
supporting evidence for certain formation theories of the Milky Way, such as
the ‘inside-out’ formation theory (Larson, 1976). This theory suggests that
the negative radial ($\Delta$[Fe/H]/$\Delta$R) metallicity gradient in the
Milky Way could be caused by the inner Galaxy forming first and fast, leading
to a higher metallicity concentration towards the Galactic center. As time
proceeds, the outer Galaxy starts to form causing it to be metal-poor compared
to its centralized counterpart (Matteucci & Francois, 1989; Chiappini et al.,
1997).

While the metallicity gradients are some of the most prominent features
observed in the Galaxy, there are more subtle chemical characteristics that
appear to be washed out by these strong trends. Chemical azimuthal
substructure in the Milky Way has been previously identified using a variety
of different tracers such as Hii regions (Balser et al., 2011, 2015) and
Cepheids (Lemasle et al., 2008; Pedicelli et al., 2009). Recently, using
large-scale stellar surveys, azimuthal variations in [Fe/H] throughout the
disk are quantifiable on the level of $\sim 0.1$ dex (Hawkins, 2022; Poggio et
al., 2022; Imig et al., 2023).

One motivator for looking for angle-dependent chemical variations in our
Galaxy is the confirmed existence of azimuthal variations in galaxies beyond
the Milky Way. Ho et al. (2017) found [O/H] azimuthal substructure in NGC
1365, proposing that the spiral arms in this galaxy plays a large role in the
dispersion of metals throughout the galaxy. Hwang et al. (2019) show clear
azimuthal metallicity variations in SDSS IV MaNGA galaxies in Figure 13.
Kreckel et al. (2019) observed subtle azimuthal variations in 4 out of 8 in
their sample of nearby galaxies from the PHANGS–MUSE survey, with ranging
associations between the metallicity variations and the spiral arms. This
range in correlations with the spiral arms calls into question if spiral arms
are the most adequate explanation for the causes of azimuthal chemical
variations.

There are two distinct mechanisms that spiral arms can generate to induce
azimuthal metallicity variations. The first process is radial
migration/churning (Sellwood & Binney, 2002) which, in this context, is the
morphing of stellar orbits that alters the angular momentum of an orbit
without the addition of excess energy. The second process is kinematic
heating/blurring which would induce changes (heat) in the stellar orbital
parameters without increasing the angular momentum. An example of this is the
spiral arms dynamically evoking changes in the motions of stars along the
leading and trailing edges of the spiral arms (Grand et al., 2012, 2016).

In simulations, Khoperskov et al. (2018) found that azimuthal variations in
metallicity may arise from the dynamics of stellar disks alone, without the
need of radial migration to reshape the stellar population. This is due to
dynamically cooler populations in the disk showing a larger contribution to
spiral arms than dynamically hotter populations, leading to azimuthal
variations in metallicity. Khoperskov & Gerhard (2021) claimed that in various
phase-space coordinates, stars in an angular momentum overdensity caused by
the spiral arms also exhibit a mean metallicity that is systematically higher
than stars not in an angular momentum overdensity.

Using Gaia DR3 stars, Poggio et al. (2022) found that the azimuthal variations
present in the disk correlate with where the spiral arms are predicted to be
using different samples of bright stars (using effective temperature as a
proxy for age) within $\sim 4$ kpc of the Sun. Interestingly, they concluded
the inverse of the simulated results from Di Matteo et al. (2013), in that the
youngest population exhibits the most pronounced inhomogenities.

Debattista et al. (2024) quantified azimuthal metallicity variations
($\delta$[Fe/H]) in a Milky Way-like galaxy simulation and found variations
that are comprable with the magnitude of the azimuthal variations seen in the
Galaxy. Their azimuthal variations in metallicity were coincident with the
spiral density waves and present in both young and old populations of stars.
When looking at the pattern speeds of the $\delta$[Fe/H] variations, they
found that these pattern speeds matched those of the spirals, pointing to the
spiral arms as the root cause of these variations.

There are a handful of explanations, aside from spiral arms, for these
azimuthal variations that have been introduced in the literature. One
possibility for these angle-dependent trends could be secular processes, such
as influence and migration due to the Galactic bar (Di Matteo et al., 2013).
Specifically, it has been speculated that the presence of azimuthal variations
in an old stellar population is a probe of bar activity. Using $N$-body
simulations, Filion et al. (2023) found that the bar induces azimuth-dependent
trends in stellar radii, which then coincides with angle-dependent variations
in metallicity using younger to intermediate-age stars ($\sim 1-4$ Gyr).

A recent study done suggested that the azimuthal variations could arise from
interactions with a satellite galaxy, such as Sagittarius (Carr et al., 2022).
Torques from the gravitational interaction between an external satellite and
the disk of the host galaxy can cause radial migration of the stellar
population where inward-migrating stars would be more metal-poor on average
compared to in-situ populations and outward-migrating populations would be
more metal-rich than the in-situ outer Galactic population (with an assumed
negative radial metallicity gradient). This migration response will be induced
with respect to the location of the perturber, which then stimulate the
azimuthal variations between inward and outward migrating populations. With
the authors’ simulations, they posited that the influence of Sagittarius will
be the most prominent in the outer disk, a conclusion also drawn from Laporte
et al. (2018). Observationally, Hwang et al. (2019) found gas-phase azimuthal
metallicity variations in close or interacting galactic pairs, providing
evidence for merger-induced azimuthal metallicity variations.

Most likely, there will be a contribution from all of the aforementioned
dynamical processes that could cause azimuthal variations. Empirical
observations are necessary to try to disentangle the processes that are
responsible for the observed azimuthal metallicity variations in the Galactic
disk. While the azimuthal variations of [Fe/H] in the Galactic disk have been
characterized by different tracer populations such as APOGEE red giants
(Eilers et al., 2022) and LAMOST OBAF-type stars as well as Gaia (Hawkins,
2022; Poggio et al., 2022), little work has been done when looking beyond iron
towards other elements, investigating the effect that stellar population age
has on the intensity of azimuthal variations, and examining the dynamical
perspective.

Figure 1: The face-on (top) and edge-on (bottom) distribution of our thin disk
sample of 202,510 stars in which the the hexagonal bins are colored by the
logarithmic number of stars. The black contours in the top panel are the
spiral arms of the Milky Way determined by Reid et al. (2019). The orange star
in both panels represents the Sun.

Throughout this work, we aim to confirm the azimuthal variations in the
Galactic disk and investigate any correlation with the spiral arms of the
Galaxy. We additionally split our sample by age to quantify the effect that
stellar age has on the intensity of the azimuthal variations, examine
azimuthal variations in elements beyond Fe, and probe the dynamical origins of
these azimuthal metallicity variations. In Section 2, we delineate the dataset
adopted for this project. In Section 3, we explain the methods taken to
achieve our goals. In Section 4.1 we present our radial and vertical [Fe/H]
gradients, and highlight the azimuthal deviations from the radial gradient in
Section 4.2. We quantify azimuthal variations in multiple other elements in
Section 4.3. We separate our sample into distinct age bins and characterize
how stellar populations of different ages have disparities in their azimuthal
variations in Section 4.4. In Section 4.5, we link the stellar chemistry to
kinematics. We summarize our results in Section 5.

##  2 Data

The initial sample of stars came from the seventeenth data release (DR17,
Abdurro’uf et al., 2022) of the Apache Point Observatory Galactic Evolution
Experiment (APOGEE) (Majewski et al., 2017b). APOGEE is a large scale near-
infrared (15140Å < $\lambda$ < 16960Å) stellar spectroscopic survey. The
survey spans both hemispheres, consisting of a spectrograph on the 2.5-meter
Sloan Foundation Telescope at the Apache Point Observatory in New Mexico, USA
as well as the 2.5-meter Irénée du Pont Telescope at the Las Campanas
Observatory in Chile.

We made use of the astroNN value-added catalog (Leung & Bovy, 2019a)111The
astroNN package is available here: https://github.com/henrysky/astroNN of
abundances, distances, and ages for APOGEE DR17 sources. astroNN is an open
source python package developed for the neural network trained on the APOGEE
data and is designed to be a general package for deep learning in astronomy.
The APOGEE-astroNN catalog contains results from applying astroNN neural nets
on APOGEE DR17 spectra to infer stellar parameters, abundances trained with
ASPCAP DR17 (Holtzman et al., in prep.), distances retrained with Gaia eDR3
(Gaia Collaboration et al., 2021) from Leung & Bovy (2019b), and ages trained
with APOKASC-2 (Mackereth et al., 2019) in combination with low-metallicity
asteroseismic ages (Montalbán et al., 2021).

The neural network from Leung & Bovy (2019a) mimics a ‘standard’ spectroscopic
analysis by using the full wavelength range to deduce stellar parameters and
specific sections of the spectrum to derive individual elemental abundances.
astroNN takes into account incomplete and noisy training data while applying
dropout variational inference to find uncertainties on the measurements. This
catalog contains stellar parameters (Teff, log $g$, and [Fe/H]) as well as 18
individual element abundances with precisions at the $\sim 0.03$ dex level,
agreeing quite well with the traditional ASPCAP pipeline and producing a
smaller scatter with tighter uncertainties.

For the distances in this work, we adopted the weighted distance in APOGEE-
astroNN. This parameter is a weighted combination of the astroNN distance
(spectro-photometric calibrated distances from Leung & Bovy, 2019b) and Gaia
parallax (Gaia Collaboration et al., 2022a). The Galactocentric X, Y, and Z
coordinates were found as transformations from the Galactocentric cylindrical
radius, azimuth, and vertical height given by astroNN. The Galactocentric
positions and velocities were computed assuming the Sun is located at 8.125
kpc from the Galactic center (GRAVITY Collaboration et al., 2018), 20.8 pc
above the Galactic midplane (Bennett & Bovy, 2019) and has radial, rotational,
and vertical velocities of -11.1, 242, and 7.25 km/s, respectively (Schönrich
et al., 2010; Bovy et al., 2012).

To explore the effect of stellar age on the chemical azimuthal variations, we
sorted our sample into 3 distinct groups with the ages in APOGEE-astroNN. Due
to the tendency to under-predict old ages, Leung & Bovy (2019a) and Mackereth
et al. (2019) used a non-parametric Locally Weighted Scatterplot Smoothing
(LOWESS) to correct the ages. While this inconsistency mainly applies to stars
older than $8+$ Gyr, we followed the recommendation and applied the LOWESS
corrected ages for the analysis in this work.

To explore the intersection of chemistry and dynamics, we aim to quantify how
metallicity excess ($\delta$[Fe/H]) changes with respect to the dynamical
parameters radial action (Jr), angular momentum (Lz) which is proportional to
the azimuthal action (Jϕ), vertical action (Jz), orbital eccentricity (e),
maximum height above the Galactic plane (Z${}_{\text{max}}$), and total
energy. These parameters were found by integrating the stellar orbits with the
Gala code designed primarily for Galactic dynamics evaluations (Price-Whelan,
2017). The mass-model of the Milky Way used for the gala.dynamics.orbit
function was the MilkyWayPotential2022 that has been fit to the rotation curve
in Eilers et al. (2019) and incorporates the phase-space spiral in the solar
neighborhood set by Darragh-Ford et al. (2023).

Figure 2:  The [Fe/H] abundances with respect to Galactocentric radius of our
planar thin disk sample. The artifact at R $\sim 8$ kpc is an observational
effect of the over-representation of stars in the solar neighborhood. We fit a
linear model to the colored data points (column-normalized planar thin disk
sample) which is represented by the grey line. From this, we obtain a
metallicity gradient of $\Delta$[Fe/H]/$\Delta$R -0.066$\pm 0.0004$ dex/kpc
and a y-intercept of 0.546 dex for the stars in our planar thin disk sample.

With the APOGEE-astroNN catalog, we employed the following cuts to obtain a
set of stars that is well-sampled and a reliable representation of the
population we aim to characterize:

  1. 1.

The Galactocentric distances are an integral part of this work to characterize
the metallicity gradient in the Galactic disk. To ensure we have precise
positions we place an error cut on the distance measurements obtained from the
astroNN dataset, selecting stars with errors $<30\%$.

  2. 2.

To minimize systematic effects and dwarf contamination in our sample, we
removed stars outside of the effective temperature range $3500<$ Teff $<5000$
K as well as any star with log $g$ $>3.6$ dex. This selection criterion
establishes that we are sampling the true red giant section of the Color-
Magnitude Diagram (CMD).

  3. 3.

We removed any stars with an error larger than 0.08 dex on [Fe/H], [O/H], or
[Mg/H] due to these elements being the most relevant in our chemical
cartography. The limit of 0.08 dex was chosen based on previous studies that
have found the dispersion of line-to-line abundances in APOGEE is $\sim 0.08$
dex for [Fe/H] (Chen et al., 2015; Hawkins et al., 2016) We placed no
selection cut on any of the other elements to maintain a balance between the
robustness and size of our sample. For completeness we include all of the
elements in Section 11 but we refrain from drawing any conclusions about
Galactic evolution from the elements with no error cuts (elements other than
[Fe/H], [O/H], and [Mg/H]).

  4. 4.

Chemical composition and evolution varies differently among the thin disk,
thick disk, and halo. We select only the kinematic thin disk stars based on
the processes outlined in Ramírez et al. (2013) (full derivation shown in
Appendix A). In this work, the authors outlined a probabilistic approach to
assign membership probabilities for stars to thin disk, thick disk, or halo
populations based on their Galactic space velocities. We required all thin
disk stars to have a TD/D probability $>80\%$, providing us with a final
sample of 202,510 stars, shown in Figure 1.

With this full thin disk sample, we are able to characterize the vertical
metallicity gradient as a function of Galactocentric radius as well as the
radial metallicity gradient as a function of height above the plane. To
investigate the chemical gradients and azimuthal variations in the plane of
the disk, we further sub-sample the 202,510 stars in the following section.

###  2.1 Planar Thin Disk Sample

Figure 3:  The radial metallicity gradient as a function of absolute vertical
height Z above (and below) the plane. The gradient derived in this work is
represented by the black dots, compared to a variety of other studies that use
different tracers [Cheng et al. (2012) (purple triangles), Hawkins (2022)
(blue triangles), Hayden et al. (2014) (pink squares), Yan et al. (2019)
(green squares)]. Consistently, $\Delta$[Fe/H]/$\Delta$R starts off at its
most negative in the plane and shallows out with greater distances from the
disk. A vertical line is plotted at $\Delta$[Fe/H]/$\Delta$R=0 to illustrate
where the radial gradient is no longer negative. The points derived in this
work lie generally in the middle of the other studies conducted.

The density of stars in the stellar disk of the Milky Way is thought to
noticeably decline at about $\sim 15.5$ kpc away from the Galactic center
(Momany et al., 2006; Reylé et al., 2009; Carraro et al., 2010). The inner
radii of the Galactic disk are thought to be dominated by the bar, typically
believed to extend out to $\sim 3.5$ kpc (Hammersley et al., 1994; Wegg et
al., 2015; Lucey et al., 2023). Thus to minimize the chemical and dynamical
effects driven by the bar in the most central region, and the drop in the
density of stars in the furthest regions, we adopted a Galactocentric radius
cut of $3.5<R<15.5$ kpc to our final thin disk sample. We place a 0.3 kpc cut
on Z${}_{\text{max}}$ to ensure that we are only selecting bona fide thin disk
stars with orbits confined to the Galactic plane. After these final selection
cuts, we have a planar thin disk sample of 32,768 stars with the median
uncertainties for Teff, log $g$, and [Fe/H] being 33K, 0.07 dex, and 0.03 dex,
respectively.

##  3 Methods

Table 1: Radial Metallicity Gradient Parameters for $\left|\text{Z}\right|$ bins in the full thin disk sample $\left|\text{Z}\right|$ |  $\Delta$[Fe/H]/$\Delta$R |  $\sigma\Delta$[Fe/H]/$\Delta$R | N${}_{\text{stars}}$
---|---|---|---
(kpc) | (dex/kpc) | (dex/kpc) |
0.1 | -0.0672 | 0.0004 | 40434
0.3 | -0.0579 | 0.0005 | 42167
0.5 | -0.0492 | 0.0006 | 35415
0.7 | -0.0389 | 0.0007 | 25742
0.9 | -0.0274 | 0.0009 | 14874
1.1 | -0.0218 | 0.0011 | 9051
1.3 | -0.0153 | 0.0013 | 5780
1.5 | -0.0132 | 0.0015 | 3573
1.7 | -0.0126 | 0.0016 | 2298
1.9 | -0.0107 | 0.0018 | 1580

The behavior of the metallicity of stars throughout the Galactic disk can be
best characterized by a negative linear gradient (e.g. Janes, 1979; Rolleston
et al., 2000). Using our planar thin disk sample, we derived the radial and
vertical metallicity gradients stars as follows.

The gradients are initially represented as linear functions in which the
gradients are represented by $m_{R}$ and $m_{z}$ in Equations 1 and 2. In
this, [Fe/H]${}_{\text{R}}$ is the metallicity at a certain Galactocentric
radius (R) and $m_{R}=$$\Delta$[Fe/H]/$\Delta$R. To account for the
metallicity at the Galaxy’s center, we introduce the term $b_{R}$:

$\mathrm{[Fe/H]_{R}}=m_{R}\mathrm{R}+b_{R}$ (1)

Similarly, we model the vertical metallicity gradient as:

$\mathrm{[Fe/H]_{Z}}=m_{Z}\mathrm{|Z|}+b_{Z}$ (2)

where $m_{z}$ and $b_{z}$ are constants.

To find the gradients, we fit linear models to the equations above with
scipy.stats.linregress (Virtanen et al., 2020) which calculates a linear
least-squares regression for two sets of measurements (in our case
Galactocentric radius/height and [Fe/H]). This returns the slope of the line
(metallicity gradient), the standard error of the gradient, and the
y-intercept (metallicity at the Galactic center/midplane). When computing the
vertical gradient, we use our full thin disk sample without a radius or height
cut. The resulting gradients are in Section 4.1. In Tables 1 (and 2), we
compute the linear (vertical) metallicity gradients in radial (Galactic
height) bins.

We use our 1D models to create 2D metallicity residual maps to search for
potential signature of azimuthal variations. We subtract the model abundances
(in which the abundances are exactly equal to the linear gradient) from the
observed [Fe/H] of each star in our sample. If the linear gradient is the only
chemical feature of the disk, we would expect to see stochastic noise. If
there is structure in this noise, then there are other processes driving the
chemistry of the disk. We plot these residuals and investigate this azimuthal
substructure in Section 4.2. When comparing the metallicity excess to azimuth,
it is important to note that we define the line of sight from the Galactic
center through the Sun to be an azimuth of $\phi(\pi)=1$.

Figure 4:  The vertical metallicity gradient as a function of Galactocentric
radius. The gradient derived in this work is depicted by black dots, whereas
the other colored points represent the gradient as determined by different
tracers in other studies [Hayden et al. (2014) (pink squares), Hawkins (2022)
(blue traingles), Nandakumar et al. (2020) (purple triangles), Yan et al.
(2019) (green squares)]. The vertical metallicity gradient is at its most
negative closest to the Galactic center and shallows out
($\Delta$[Fe/H]/$\Delta$Z approaches 0) as distance increases.

The linear gradient is a useful approximation for [Fe/H], however this does
not necessarily hold true for all elements (see Appendix B). To attempt to
account for any non-linearity in our sample of elements, we employed a
‘running median’ method. We began by calculating the median [X/Y] abundances
in 0.2 kpc radial bins within the bounds of $\sim 6-15$ kpc. In these bins, we
then took the median [X/Y] value and following the linear method, subtracted
$Data-Median$ to quantify how the individual stars may deviate azimuthally
from the median [X/Y] abundance. However, when this approach was taken, the
radial bins did not sample azimuthal angles isotropically due to APOGEE only
probing one part of the disk. Thus, for the context of this work we only
select elements that behave monotonically with respect to Galactocentric
radius and move forward with the linear gradient method in Section 4.3.

##  4 Results and Discussion

Table 2: Vertical Metallicity Gradient Parameters for R bins in the full thin disk sample $\left|\text{R}\right|$ |  $\Delta$[Fe/H]/$\Delta$Z |  $\sigma\Delta$[Fe/H]/$\Delta$Z | N${}_{\text{stars}}$
---|---|---|---
(kpc) | (dex/kpc) | (dex/kpc) |
7.5 | -0.2307 | 0.0026 | 29634
8.5 | -0.1785 | 0.0022 | 49208
9.5 | -0.1383 | 0.0025 | 33837
10.5 | -0.0951 | 0.0026 | 22099
11.5 | -0.0632 | 0.0028 | 14729
12.5 | -0.0617 | 0.0030 | 8025
13.4 | -0.0518 | 0.0034 | 3090

In this section we present our radial and vertical metallicity gradients
(Section 4.1), as well as azimuthal [Fe/H] variations throughout the thin disk
(Section 4.2). We uncover azimuthal variations in elements beyond Fe in
Section 4.3. To investigate which mechanisms may be responsible for these
variations, we bin our sample by age to explore the age dependence on
azimuthal variations (Section 4.4) and close out the section by examining the
link between metallicity excess and several dynamical properties (Section
4.5).

###  4.1 Metallicity Gradients

Figure 5: Left: The full APOGEE planar thin disk sample colored by the
metallicity. Middle: Each position of the data points colored by a model
gradient of $\Delta$[Fe/H]/$\Delta$R = $-0.066\ +0.546$ dex/kpc. Right: the
residuals of the observed [Fe/H] abundances and the linear model abundances.
Figure 6: The azimuthal metallicity variations ($\delta$[Fe/H] in the last
panel of Figure 5) as compared to different determinations of the spiral arms.
The black contours in the top panel represent the spiral arms as derived by
Poggio et al. (2021) using main sequence stars in Gaia and the solid black
lines in the bottom panel are the spiral arms as determined by Reid et al.
(2019) using high-mass star-forming regions.

Our first step in this analysis was to measure the radial
($\Delta$[Fe/H]/$\Delta$R) and vertical ($\Delta$[Fe/H]/$\Delta$Z) metallicity
gradients in our planar thin disk sample that have been previously observed
with different tracer populations. Employing the methods outlined in Section
3, we obtained a radial metallicity gradient with our planar thin disk sample
of $\sim-0.066\ \pm 0.0004$ dex/kpc (Figure 2) as well as a vertical
metallicity gradient of $\sim-0.164\pm 0.001$ dex/kpc. Overall our gradients
match within reason to other recently calculated gradients. Hawkins (2022)
derived a radial metallicity gradient of $\sim-0.078\pm 0.001$ dex/kpc in the
Galactic disk and a vertical metallicity gradient of $\sim-0.15\pm 0.01$
dex/kpc with an OBAF-type stellar sample from the LAMOST survey. Using the
sixteenth data release of APOGEE, Eilers et al. (2022) determined a radial
gradient of $\sim-0.057\pm 0.001$ dex/kpc. Imig et al. (2023) calculated the
running median radial metallicity gradients with APOGEE DR17 ASPCAP parameters
and found $\Delta$[Fe/H]/$\Delta$R$=-0.064\pm 0.001$ dex/kpc in their
low-$\alpha$ sample at $|Z|\leq 0.25$, which is in good agreement with our
gradient of $\sim-0.066\ $ dex/kpc. We note one of the main distinctions in
our sample compared to other works in the literature is the cut we imposed on
Z${}_{\text{max}}$ as opposed to present day height above the plane.

Figure 7:  The $Data-Model$ panel for the elements we deem to behave
monotonically with respect to Galactocentric radius. The contours are colored
by the shared [X/Fe] colorbar to illustrate the varying intensities of the
deviations from the linear radial gradient. The second panel represents the
average [$\alpha$/Fe] abundance, with the individual $\alpha$-elements in the
following panels. The elements with the most saturated contours, such as
[Fe/H], showcase the most exaggerated deviation from the radial gradients. The
$\alpha$-elements appear to be loosely anti-correlated with [Fe/H], following
predictions from the difference in timescales between events that mainly
produce $\alpha$-elements (Type II SN) and events that mainly produce [Fe/H]
(Type Ia SN).

The radial metallicity gradient is illustrated in Figure 2. The metal
enhancement of the inner Galaxy as compared to the outer Galaxy could point to
the ‘inside-out’ formation theory (Larson, 1976; Kobayashi & Nakasato, 2011).
Starting at the inner Galaxy (i.e., R $\sim$3.5 kpc), the metallicity
decreases linearly $-$ for every kiloparsec moving towards the outer disk, the
global metallicity decreases by $\sim 0.07$ dex.

Figure 3 explores how $\Delta$[Fe/H]/$\Delta$R varies with height above (and
below) the Galactic plane. We separate the full thin disk sample into 0.2 kpc
steps of —Z— and plot these results along with a select few other studies that
use different tracer populations: LAMOST OBAF-type stars (Hawkins, 2022),
SEGUE main sequence turn-off stars (Cheng et al., 2012), APOGEE red giant
stars (Hayden et al., 2014), and LAMOST FGK-type stars (Yan et al., 2019).
Regardless of tracer population, all studies show a clear trend that as the
distance from the Galactic plane increases, the radial metallicity gradient
becomes more shallow, approaching zero. Physically, this means that the radial
metallicity gradient is a prominent feature in the disk, but with increasing
height from the Galactic plane the radial gradient is no longer the dominant
observed relation and the vertical metallicity gradient starts to take over.
We tabulate these results in Table 1. This work’s $\Delta$[Fe/H]/$\Delta$R vs.
—Z— trend lies generally in the center of the other studies’ trends. When
compared with Hayden et al. (2014) which used the same survey as this work,
merely an earlier data release (DR10 compared to our DR17), our work seems to
be shifted up by $\sim 0.01$ dex. This is not surprising due to the results of
Jönsson et al. (2018) showing that APOGEE data releases can have systematic
differences of less than 0.05 dex.

Recent studies have shown that there is a vertical metallicity gradient that
varies as a function of Galactocentric radius. Önal Taş et al. (2016) found
that the radial gradient is flat within 0.5 - 1 kpc of the plane and then
becomes positive greater than 1 kpc away from the plane, suggesting that there
is a vertical metallicity gradient in the Galaxy. To derive the vertical
gradient, we applied the same methodology used for the radial metallicity
gradient, obtaining $\Delta$[Fe/H]/$\Delta$Z $\sim-0.164\pm 0.001$ dex/kpc.

To quantify how the vertical gradient changes as a function of Galactocentric
radius, we follow a similar methodology to the evaluation of the radial
gradient varying with Galactic height, this time splitting our full thin disk
sample of stars into radial bins of 1 kpc (Table 2). In Figure 4, we show our
derived vertical metallicity gradient using the APOGEE red giants, denoted by
black circles. Other works that we compared to include vertical gradients
obtained from APOGEE red giant stars (Hayden et al., 2014), LAMOST OBAF-type
stars (Hawkins, 2022), GALAH FGK-type stars (Nandakumar et al., 2020), and
LAMOST FGK-type stars (Yan et al., 2019). We find that the vertical
metallicity gradient is heavily correlated with Galactocentric radius in that
$\Delta$[Fe/H]/$\Delta$Z approaches zero with increasing distance from the
Galactic center.

###  4.2 Azimuthal Variations in $\Delta$[Fe/H]/$\Delta$R

Figure 8:  This figure is similar to Panel 3 of Figure 5 increasing age groups
of $0-2,3-6,$ and $6+$ Gyr from left to right with 8730, 18493 and 5545 stars
respectively. The color bars are constant throughout all panels, thus the
youngest age group has the least amount of contrast and represents the
smallest deviations from the linear gradient. Conversely, the oldest age group
has the most saturated colors due to the larger variations from the modelled
linear gradient.

In this section, we aim to characterize the angle-dependent deviations from
the linear radial metallicity gradient in our planar thin disk sample and
investigate any correlations with these deviations and the spiral arms. The
discernible radial metallicity gradient has been recognized and characterized
for decades (e.g., Mayor, 1976; Andrievsky et al., 2002; Magrini et al., 2009;
Boeche et al., 2013, 2014; Cunha et al., 2016). With increased sample size and
precision provided by advances in instrumentation throughout the years, we are
now able to peel back the gradient to see if there is any under-lying
structure beneath this strong signature.

To look beyond the metallicity gradient, we follow the steps outlined in
Section 3 and compute the model linear $\Delta$[Fe/H]/$\Delta$R gradient
throughout the disk. We subtract this model off of the data to search for any
structure in the residuals. When applied, we find variations in the [Fe/H]
abundances that correlate with azimuthal angle.

We plot the top-down view of the Milky Way with our planar thin disk data in X
and Y colored by [Fe/H] in Panel 1 of Figure 5. We illustrate the model linear
metallicity gradient in the second panel of Figure 5 and in the final panel we
show the residuals that are found after subtracting the model from the data
(labelled $\delta$[Fe/H]).

Figure 9:  This figure shows $\delta$[Fe/H] as a function of azimuthal angle
($\phi$). For reference, $\phi(\pi)=1$ is the line of sight from the Sun
towards the Galactic center. The youngest age group, $0-2$ Gyr (blue squares),
displays the smallest deviations from zero in $\delta$[Fe/H]. The intermediate
age group, $2-6$ Gyr (green triangles), is the largest sample and has
noticable deviations from zero. The oldest age group, $6+$ Gyr (purple
circles), varies the most dramatically in $\delta$[Fe/H] azimuthally. This
depiction matches the theorized azimuthal variations, discussed further in
Section 4.4.

There is a blue-red-blue-red pattern in $\delta$[Fe/H] shown in the final
panel of Figure 5. The red bins represent areas of the Milky Way that are more
metal-rich than the model predicts and the blue bins are where the model over-
estimates the stellar metallicities. This oscillating pattern has been
observed previously with different stellar surveys, such as Gaia DR3 stars in
Figure 2 of Poggio et al. (2022), with similar results to this work. Hawkins
(2022) use LAMOST OBAF stars as well as Gaia DR3 stars, uncovering azimuthal
variations in both of these populations, showing similar patterns in which the
solar-neighborhood is more metal-rich (red) and the inner/outer galaxy are
more metal-poor (blue).

To probe different formation pathways for this oscillating pattern, we first
compare these results with the location of the spiral arms. Khoperskov et al.
(2018) (and Grand et al., 2012) postulate that metal-rich stars from the
center of the Galaxy may travel along the spiral arms, causing these filaments
in the outer Galaxy to appear more metal rich. Khoperskov et al. (2018) show
that azimuthal [Fe/H] variations in spiral galaxies would arise naturally if
there is an initial negative radial metallicity gradient due to the migration
of stars. To detect any correlation with the spiral arms, we plot our
$\delta$[Fe/H] residuals and different determinations of the spiral arms in
Figure 6.

The black solid lines in the bottom panel of Figure 6 represent the spiral
arms as determined by Reid et al. (2019) using the parallaxes and proper
motions associated with high-mass star-forming regions using a Very Long
Baseline Array, European VLBI Newtork, and the Japanese VLBI Exploration of
Radio Astrometry project. They locate multiple arm segments with pitch angles
ranging from 7∘ to 20∘ with the widths of the arms increasing with distance
from the Galactic center. The black contours in the top panel of Figure 6
represent the spiral arm locations determined by Poggio et al. (2021) using
the over-density of upper main sequence stars determined by Gaia DR3. These
results were consistent with some of the Reid et al. (2019) arms, such as the
Sagittarius-Carina spiral arm, while the geometry of arms with Galactic
longitudes from 180∘ to 270∘ were significantly different from other spiral
arm models. The metal-rich and metal-poor portions of the Milky Way in our
work are not fully encompassed by either determination of the spiral arm
locations, thus we look towards other processes that may have induced this
azimuthal substructure.

If the spiral arm locations do not obviously coincide with the azimuthal
substructure, this begs the question, what is causing these variations? Carr
et al. (2022) quantified the effect that a Sagittarius (Sgr)-like dwarf galaxy
would have on the kinematics and chemsitry of the Milky Way-like disk upon
first interaction. They did this by painting the particles in their N-body
simulation of the disk with an negative radial [Fe/H] gradient. The passage of
Sgr through the plane will cause radial rearrangement and disrupt stellar
orbits. If the heating of orbits is due to a non-axisymmetric feature, such as
Sgr, migration and mixing signatures will manifest as an approximate
quadrupole in chemical azimuthal variations across the disk. In their present
day snapshot of a simulated Milky Way being perturbed by a Sgr-like dwarf
galaxy, Carr et al. (2022) found that the maximum departure from the
$\Delta$[Fe/H]/$\Delta$R gradient occurs in the solar annulus on the side of
disk closest to the current position of Sgr, while the minimum is found
adjacent to Sgr in the outermost annuli. These azimuthal variations in
metallicity will be dependent on the model of Sgr chosen. It is interesting to
note that in our sample the solar annulus appears to be more metal-rich,
although we lack a complete view of the Galaxy to classify this definitively.

While external influences like Sgr can cause azimuthal metallicity variations,
studies have also shown that these features can arise through secular
evolution with structures like the Galactic bar. Filion et al. (2023)
simulated the effects of radial rearrangement in a barred galaxy and find
substantial changes in the radii of stars that, when paired with the negative
radial metallicity gradient, will induce azimuthal substructure similar to
what we find with the APOGEE stars (see also Di Matteo et al., 2013). The
radial rearrangement of stars due to the bar can be split into 3 zones to
characterize the dynamics of the stars. The effects driven by the bar would be
most prominent in the inner Galaxy where orbits evolve inward due to the
angular momentum loss of the stellar orbits. Their intermediate zone is
composed of stars moving both inwards and outwards producing no mean radial
evolution. In the outer zone of their simulations, orbits generally evolve
outwards (likely due to the net effect of the bar moving angular momentum
outwards) and the trends are less aligned with the bar angle. Any migration of
stars across the disk will produce azimuthal metallicity variations. In line
with their findings, we see that the inner-disk stars of our sample have a
lower metallicity than predicted by the gradient alone.

To discern the responsible mechanisms for this substructure, it is crucial to
confirm that this pattern is not unique to iron but extends to other elements
as well. While the azimuthal [Fe/H] variations have been quantified before
using Gaia (Hawkins, 2022; Poggio et al., 2022), APOGEE gives us access to a
suite of other elements. To confirm that this substructure is present in
elements aside from iron, we evaluate other azimuthal chemical variations in
the following section.

###  4.3 Azimuthal Variations in Other Elements

Here, we explore azimuthal variations in [Fe/H] and $\alpha$-elements (O, Mg,
Si, S, Ca) in our planar thin disk sample. We choose to specifically highlight
the $\alpha$-elements because the abundances of these elements in stars appear
to behave monotonically with respect to Galactocentric radius, similar to
[Fe/H]. For this methodology, we used our planar thin disk sample and plotted
the median abundances of these elements in 0.2 kpc radial bins in Figure 11.
As not all of these elements are monotonic with respect to Galactocentric
radius, we focus on the $\alpha$-elements. We fit linear gradients to these
elements and subtracted off the linear model (explained further in Section 3)
to look for structure in the residuals (similar to [Fe/H]), shown in Figure 7.

In Figure 7, we show azimuthal variations with ranging intensities for all
elements in our analysis, with variations detectable on at least the $\sim
0.05$ dex-level. The second panel in this figure represents the average
[$\alpha$/Fe] abundance, while the individual [$\alpha$/Fe] abundances are
found in the following panels. We find that the $\alpha$-elements are loosely
anticorrelated with [Fe/H], which aligns with expected Galactic chemical
evolution (e.g., Tinsley, 1980); in the disk where there is a $\delta$[Fe/H]
excess, there is a $\delta$[$\alpha$/Fe] deficit.

The anti-correlation between [Fe/H] and [$\alpha$/Fe] can be explained by the
astrophysical processes (and timescales) that are largely responsible for the
produciton of these elements. At early times, Type II SNe were able to
effectively disperse $\alpha$-elements (e.g., O, Mg, Si, S, and Ca) as well as
small amounts of Fe that were synthesized during the lifetimes of the first
generations of stars. As the Galaxy continues to evolve, Type Ia SNe will then
dominate, causing the [$\alpha$/Fe] abundances to lower due to more Fe-peak
(Fe, Cr, Mn, Co, Ni) elements being made available to the next generation of
stars (e.g., Tinsley, 1980; Greggio & Renzini, 1983; Woosley & Hoffman, 1992;
Matteucci & Recchi, 2001; Gonzalez et al., 2011). Following this sequence of
events, wherever there are [Fe/H]-excess regions, there will be
[$\alpha$/Fe]-defecit regions and vice versa.

The presence of azimuthal variations in other elements confirms the existence
of some process that is causing the observed chemical substructure throughout
the Galactic disk. When looking for the responsible process, Filion et al.
(2023) noted that focusing on younger populations ($\sim 1-2$ Gyr) in the
inner disk would help discern the role of secular evolution. Other simulations
have used exclusively older stars to identify the cause of azimuthal chemical
variations (Di Matteo et al., 2013). Consequently, we delve into the possible
age-dependencies observed in azimuthal metallicity variations in the following
section.

###  4.4 Azimuthal Metallicity Variations by Age

In this section, we aim to quantify the relationship between stellar age has
and the magnitude of the azimuthal substructure found in Section 4.2. We are
particularly interested in whether the amplitude of the azimuthal [Fe/H]
variations are stronger in older or younger populations. We additionally aim
to explore how the [Fe/H] gradient and azimuthal variations change as a
function of age. For reference, it is important to note the difference between
look-back time and present-day age. Look-back time is primarily used in
simulations, where it is the time elapsed from the ‘final’/present time back
to a previous point in time, used to analyze historical states. In
observations, we are limited to the present-day age of the object, or the
total age of the object from when it formed to the present time. Thus, we aim
to use the present-day ages of these stars and compare any trends in the ages
with simulated works.

Using present-day age, the Debattista et al. (2024) simulations found that
removing stars younger than 2 Gyr still produces azimuthal variations. This
implies that azimuthal variations are not solely primordial in origin because
older populations exhibit these variations as well. Bellardini et al. (2021)
demonstrated that azimuthal scatter increases with increasing look-back time.
This would indicate that stars that are born at earlier times will be born
with stronger azimuthal variations due to the inhomogeneous interstellar
medium (ISM) from which they form. While these two hypotheses appear to be
contradictory, we cannot probe look-back times observationally so we cannot
compare directly to Bellardini et al. (2021), but we can compare our results
to Debattista et al. (2024). Consequently, in the rest of this section we make
and test predictions about the observed relationship between azimuthal
metallicity variations and stellar present-day age.

If the azimuthal variations are exclusively natal in origin (i.e., they are a
result of variations in the gas-phase abundances at the time a population was
born), then we would predict that younger populations would exhibit strong
metallicity variations due to these populations forming more recently.
However, if azimuthal variations are exclusively the result of dynamical
processes, we would predict that older populations would show stronger
azimuthal variations because these stars would have had more time to interact
with non-axisymmetric features in the Galaxy that would cause the
rearrangement of stars, leading to the observed azimuthal metallicity
variations. To distinguish between these predictions, we applied the present-
day ages from APOGEE-astroNN to separate our sample by stellar age.

We separated the APOGEE giants into 3 distinct age groups to explore the
effects that age has on departures from the metallicity gradient. The age bins
we selected are 0-2 Gyr, 2-6 Gyr, and $6+$ Gyr, with 8730, 18493 and 5545
stars, respectively. The strongest radial metallicity gradient appears in the
middle age group with $\Delta$[Fe/H]/$\Delta$R $\sim-0.0774\pm 0.0005$
dex/kpc, outranking the youngest ($\Delta$[Fe/H]/$\Delta$R $\sim-0.0674\pm
0.0006$ dex/kpc) and oldest ($\Delta$[Fe/H]/$\Delta$R $\sim-0.0484\pm 0.0012$
dex/kpc) age groups.

The correlations we find with radial metallicity gradient and stellar age are
consistent with the findings of Wang et al. (2019) who used the LAMOST survey
to illustrate that stars of $4-6$ Gyr in age exhibit a steeper gradient than
either younger or older stars. Similar work has been done with a CoRoT and
APOGEE sample by Anders et al. (2017) where their middle-aged population of
stars had the steepest radial metallicity gradient, followed by the youngest
then oldest population. This age-gradient trend is consistent with a
systematic offset of $\sim-0.01$ dex/kpc between this work and Anders et al.
(2017).

After following the same methodology as Section 4.2 and subtracting the data
from the linear gradient model, Figure 8 shows the $\delta$[Fe/H] residuals
for each of the 3 age groups. Visually, the signatures seem to be the
strongest in the panel containing the oldest stars, while the deviations
lessen in intensity with decreasing age.

In Figure 9 we plot azimtuhal angle on the X axis and metallicity excess on
the Y axis for each of the different age bins. We show that the oldest age
group has the most extreme variations in terms of $\delta$[Fe/H], with the
youngest stars showing minimal azimuthal variations. Interestingly, there
seems to be a phase offset in the curves among the different age bins. This
phase offset is also seen in Figure 12 of Carr et al. (2022). For all age
groups, the variations are minimized at $\phi(\pi)=1$, which is most likely
due to the sample selection and the number of stars located along the line of
sight from the Sun towards the center of the Galaxy. The future releases of
SDSS V will provide better radial and azimuthal coverage to alleviate this
problem.

We find that the strongest azimuthal [Fe/H] variations are apparent in solar-
age stars. Additionally, the metallicity excess is largest in the oldest stars
(6+ Gyr) and similar in the youngest (0-2 Gyr) and intermediate aged (2-6 Gyr)
stars, with a phase offset that differentiates between them. These results are
in line with the simulated results of Debattista et al. (2024). The authors
found that azimuthal metallicity variations were still present when excluding
younger populations, whereas we find the strongest azimuthal metallicity
variations in our older populations. Although azimuthal vartiations have been
quantified in younger populations (e.g., the OBA sample in Hawkins, 2022), the
presence of the strongest azimuthal variations in our older populations
suggest an important contribution from dynamical processes to the creation of
azimuthal metallicity variations throughout the Galactic disk. Thus, in the
following section, we delve into possible dynamical origins.

###  4.5 Linking Chemistry to Dynamics

Figure 10:  The correlation between absolute metallicity excess
($\left|\delta\text{[Fe/H]}\right|$) on the y axes and different dynamical
properties on the x axes. The hexagonal bins are colored by density with the
pink stars representing the running medians in each panels with error bars
showing the standard deviation. In the Jr and eccentricity panels, we see a
clear increasing trend hinting that stars with high eccentricities and high
radial actions contribute to the variations with the largest magnitudes.

The driving mechanisms behind the chemical azimuthal substructure we see in
the Galactic disk can be predominantly categorized into either radial
migration/churning or kinematic heating/blurring. Radial migration (Sellwood &
Binney, 2002) (also sometimes known as cold torquing), hereafter referred to
as churning, is a dynamical process which altering stellar orbits through the
interactions with non-axisymmetric features in the disk. More specifically,
churning results in a change in the stellar orbits by changing the angular
momenta of the orbits without adding any excess energy (e.g., stars on
circular orbits experiencing churning will not see a change in eccentricity,
only radii) (Carr et al., 2022). Blurring (also known as kinematic heating)
heats an orbit without an increase in angular momenta of stellar orbits (i.e.,
there is an increase in eccentricity). Churning and blurring could be due to
interactions with the spiral arms of the Galaxy (Jenkins & Binney, 1990),
interactions with the bar (Schönrich & Binney, 2009; Filion et al., 2023), and
interactions with a satellite galaxy like Sgr (Carr et al., 2022).

In this section, we aim to provide the initial basis of linking the azimuthal
variations with Galactic dynamical properties to help elucidate the mechanisms
responsible for these observed chemical signatures. While the goal of this
paper is not to define the exact mechanisms responsible for azimuthal chemical
variations, this initial linking to dynamics will motivate future studies. The
dynamical properties of these stars could point to whether churning or
blurring has a greater impact on the magnitude of the azimuthal variations.
Thus we plot the determined radial action (Jr), angular momentum (Lz),
vertical action (Jz), orbital eccentricity (e), maximum height above the
Galactic plane (Z${}_{\text{max}}$), and total energy to plot against absolute
metallicity excess ($\left|\delta\text{[Fe/H]}\right|$).

Figure 10 is a six panel plot showing the aforementioned parameters versus
$\left|\delta\text{[Fe/H]}\right|$. The bins are colored by density with the
pink stars representing the running medians in each panels with error bars
showing the standard deviation. There is a slight positive trend with
$\left|\delta\text{[Fe/H]}\right|$ and Jr, suggesting that stars with a larger
metallicity excess tend to have higher radial actions. This correlation with
radial action follows the predictions of the Debattista et al. (2024)
simulation222Note that those authors suggest that Jr should be computed as a
time-averaged quantity. We do not investigate the impact of time-averaging,
and leave further investigation of this to future work.. There are no obvious
trends with the angular momentum, vertical action, maximum height, and total
energy of stars when comparing to $\left|\delta\text{[Fe/H]}\right|$.

With the first and third panels of Figure 10, we can make inferences about the
relative importance of churning and blurring assuming that there is an initial
negative metallicity gradient and no natal azimuthal variations. If churning
is the primary mechanism causing azimuthal variations, we would predict that
the largest $\left|\delta\text{[Fe/H]}\right|$ would be present at smaller
eccentricities. This is because churning would largely cause stars to migrate
from nearly circular orbits to nearly circular orbits of a different radius.
If blurring is the primary mechanism causing azimuthal variations, we would
expect the largest $\left|\delta\text{[Fe/H]}\right|$ values at larger
eccentricities (and radial actions) because blurring is generally
characterized by an increase in eccentricity.

In the first and third panels of Figure 10, we see larger
$\left|\delta\text{[Fe/H]}\right|$ at higher eccentricities and radial
actions. This implies that blurring (heating of orbits) contributes a non-
negligible amount to the mechanisms that are causing these observed azimuthal
metallicity variations. In reality, multiple Galactic mechanisms are likely
responsible for the azimuthal metallicity variations and more work is required
to pin down mechanism responsible for the variations observed in this work and
others (e.g., Hawkins, 2022; Poggio et al., 2022; Imig et al., 2023).

While the goal of this paper is not to define the exact mechanisms
responsible, this initial linking to dynamics will motivate future work on
mock observations of simulations to directly compare different dynamical
processes. We show that working in orbit-space, rather than present-day
positions, can provide valuable insight into the dynamical properties that are
driving the chemistry of the Galactic disk.

##  5 Summary

Stellar spectroscopists have been using chemical cartography to map the
chemistry of the Milky Way for the past decade (Hayden et al., 2015). With the
chemistry of stars, we are able to discern what observable patterns are
present throughout the Galactic disk and identify the key processes that
contribute to Galactic evolution. Some of the most prominent patterns
throughout the disk are the negative vertical and radial metallicity gradients
(Mayor, 1976; Andrievsky et al., 2002; Magrini et al., 2009; Luck & Lambert,
2011; Bergemann et al., 2014; Xiang et al., 2015; Yan et al., 2019; Hawkins,
2022). These gradients support theories such as ‘inside-out’ formation, which
would cause the stars in the Galactic center to form first and ignite rapid
star formation, with the outer Galaxy forming at later times. This will lead
to the inner Galaxy having more metal-rich stars and the outer Galaxy to be
populated by metal-poor stars. Hiding under the radial
$\Delta$[Fe/H]/$\Delta$R gradient are azimuthal variations on the scale of
$\sim~{}0.1$ dex (Poggio et al., 2022; Hawkins, 2022). To constrain the
properties of these azimuthal variations, we aim to confirm the [Fe/H]
azimuthal variations in APOGEE DR17, characterize how the variations interplay
with stellar age, identify if azimuthal variations exist in elements other
than [Fe/H], and attempt to link the azimuthal variations to dynamical
properties.

To accomplish these tasks, we model the linear radial metallicity gradient in
the planar thin disk and subtract this off from the actual stellar
metallicities to look for azimuthal substructure in [Fe/H] and
$\alpha$-elements. We then separate our sample by age in to 3 subsets: 0-2
Gyr, 2-6 Gyr, and 6+ Gyr. We then quantified how metallicity excess
($\delta$[Fe/H]) changes with respect to radial action (Jr), angular momentum
(Lz), vertical action (Jz), orbital eccentricity (e), maximum height above the
Galactic plane (Z${}_{\text{max}}$), and total energy. Based on these methods,
we present the following results:

  1. 1.

A radial metallicity gradient ($\Delta$[Fe/H]/$\Delta$R) of $\sim-0.066\ \pm
0.0004$ dex/kpc is found throughout the kinematic thin disk ( 3.5 $<$ R $<$
15.5 kpc and Z${}_{\text{max}}<$ 0.3 kpc) of the Milky Way. We additionally
find a vertical metallicity gradient ($\Delta$[Fe/H]/$\Delta$Z) of
$\sim-0.164\pm 0.001$. Both gradients are in fairly good agreement with
previous studies that use different tracer populations (Section 4.1).

  2. 2.

Azimuthal variations are found throughout the disk in [Fe/H]. These deviations
were quantified by modelling the linear gradient and subtracting the model
from the data to illustrate the [Fe/H]-poor to [Fe/H]-rich oscillating
azimuthal pattern in the disk. These results are consistent with other
azimuthal studies that have used different datasets or tracer populations
(Eilers et al., 2022; Hawkins, 2022). While previous studies have found a
correlation with azimuthal metallicity variations and the location of the
spiral arms, we find some deviations to that in this work. We argue this is
worth further investigation over a wider range of azimuthal angles (Hackshaw
et al. in prep) (Section 4.2).

  3. 3.

Azimuthal substructure is seen in elements beyond iron in [$\alpha$/Fe]. While
[Fe/H] exhibits the strongest variations from the gradient, the
$\alpha$-elements (O, Mg, Si, S) still show azimuthal variations on the $\sim
0.1$ dex-level. These $\alpha$-abundances manifest as anti-correlations with
[Fe/H], following the predictions by the production mechanisms between
$\alpha$ and Fe-peak elements (Section 4.3).

  4. 4.

Azimuthal substructure varies by stellar age. Older populations (6+ Gyr)
exhibit more extreme deviations from the radial metallicity gradient than seen
in younger (0-2 Gyr) and intermediate (2-6) populations, which suggests that
the azimuthal metallicity variations largely stem from dynamical processes
rather than natal processes. This age dependence can be coupled with mock
observed simulations to identify why the older populations are
disproportionately affected by the cause of these azimuthal variations
(Section 4.4).

  5. 5.

There appear to be discernible trends with absolute metallicity excess
($\left|\delta\text{[Fe/H]}\right|$) and different dynamical properties. There
is a positive trend between $\left|\delta\text{[Fe/H]}\right|$ and Jr, as well
as $\left|\delta\text{[Fe/H]}\right|$ and eccentricity. This implies that
stars with high eccentricities and high radial actions contribute to the
azimuthal variations with the largest magnitudes, hinting that blurring is an
important dynamical process in the production of azimuthal [Fe/H] variations.
Further investigation of these links can help definitively characterize the
cause of azimuthal metallicity variations (Section 4.5).

Azimuthal chemical variations are an apparent feature of the thin disk of our
Galaxy. However, the cause of this substructure is unknown, with hints lying
among the age dependence of these variations, azimuthal variations in elements
beyond iron, and the dynamics of these stars. To answer this open question, we
recommend deeper exploration of azimuthal variations with data like the
impending release of SDSS-V, providing better radial and azimuthal coverage of
the Galaxy. Including the orbit-space perspective rather than only the
positions is a key step to understanding what drives these processes.
Connecting the observable signatures in the Galactic disk with with mock
simulated observations will bring us one step closer conclusively identifying
what drives Galactic evolution.

## Acknowledgements

ZH thanks Catherine Manea, Maddie Lucey, and Malia Kao for invaluable insight
on this project.

KH & ZH acknowledges support from the National Science Foundation grant
AST-2108736 and from the Wootton Center for Astrophysical Plasma Properties
funded under the United States Department of Energy collaborative agreement
DE-NA0003843. This work was performed in part at Aspen Center for Physics,
which is supported by National Science Foundation grant PHY1607611. This work
was performed in part at the Simons Foundation Flatiron Institute’s Center for
Computational Astrophysics during KH’s tenure as an IDEA Fellow. This project
was developed in part at the Gaia Fete, hosted by the Flatiron Institute
Center for Computational Astrophysics in 2022 June. This project was developed
in part at the Gaia Hike, a workshop hosted by the University of British
Columbia and the Canadian Institute for Theoretical Astrophysics in 2022 June.
This project was developed in part at the 2023 Gaia XPloration, hosted by the
Institute of Astronomy, Cambridge University. CL acknowledges funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No. 852839).

Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P.
Sloan Foundation, the Heising-Simons Foundation, the National Science
Foundation, and the Participating Institutions. SDSS acknowledges support and
resources from the Center for High-Performance Computing at the University of
Utah. The SDSS web site is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration, including the Carnegie Institution for
Science, Chilean National Time Allocation Committee (CNTAC) ratified
researchers, the Gotham Participation Group, Harvard University, Heidelberg
University, The Johns Hopkins University, L’Ecole polytechnique fédérale de
Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-
Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für
Extraterrestrische Physik (MPE), Nanjing University, National Astronomical
Observatories of China (NAOC), New Mexico State University, The Ohio State
University, Pennsylvania State University, Smithsonian Astrophysical
Observatory, Space Telescope Science Institute (STScI), the Stellar
Astrophysics Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Illinois
at Urbana-Champaign, University of Toronto, University of Utah, University of
Virginia, Yale University, and Yunnan University.

##  Appendix A Thin Disk Sample Determination

The thick-to-thin disk probability ratios (TD/D) were found for each star
using the assumption that the space velocities U, V, and W have Gaussian
distributions (Bensby et al., 2014),

$f=k\cdot\text{exp}\Big{(}-\frac{(U_{\text{LSR}}-U_{\text{asym}})^{2}}{2\sigma_{\text{U}}^{2}}-\frac{(V_{\text{LSR}}-V_{\text{asym}})^{2}}{2\sigma_{\text{V}}^{2}}-\frac{W_{\text{LSR}}^{2}}{2\sigma_{\text{W}}^{2}}\Big{)}$
(A1)

where the normalization factor $k$ is given by

$k=\frac{1}{(2\pi)^{3/2}\sigma_{\text{U}}\sigma_{\text{V}}\sigma_{\text{W}}}$
(A2)

The characteristic velocity dispersions are represented by
$\sigma_{\text{U}},\sigma_{\text{V}}$ and $\sigma_{\text{W}}$ and
$U_{\text{asym}}$ and $V_{\text{asym}}$ are the asymmetric drifts. Finally the
thick (TD) to thin disk (D) ratio is found with

$\text{TD/D}=\frac{X_{\text{TD}}}{X_{\text{D}}}\cdot\frac{f_{\text{TD}}}{f_{\text{D}}}$
(A3)

Figure 11: The running median [X/Fe] abundances in 0.2 kpc radial bins for 15
elements in our sample. The high-fidelity elements ([Fe/H], [O/Fe], and
[Mg/Fe]) are featured on the top row. A number of these elements would not be
best represented by a linear gradient, so we leave a majority of them out of
our analysis. The elements included when searching for azimuthal variations in
the linear gradient are the approximately linear elements: Fe, O, Mg, Si, S.
To look at the $\alpha$-element behavior more holistically, we averaged the
$\alpha$-element abundances in Figure 7.

We characterize thin disk stars as any star with a thin disk probability $\geq
80\%$.

##  Appendix B Elemental Gradients

This appendix contains the [X/Fe] median abundance trends in 0.2 kpc bins in
Figure 11. Note the y-axis of the panels vary, but the elemental families show
general trends as expected, such as the Fe-peak (Cr, Mn, Fe, Co, and Ni)
exhibiting negative gradients. The $\alpha$-elements (O, Mg, Si, S, Ca) all
have positive trends. These gradients match well with shape and scale to other
studies that have looked at abundance trends throughout the Galactic disk (for
example, Figure 7 in Eilers et al., 2022).

Due to the fact that azimuths are not sampled isotropically in the 0.2 kpc
radial bins, the methodology of calculating $\Delta$[Fe/H] from the abundance
medians is inaccurate. Thus, we only choose elements with roughly linear and
monotonic shapes to complete our analysis (Fe, O, Mg, Si, S).

## References

  * Abdurro’uf et al. (2022) Abdurro’uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35, doi: 10.3847/1538-4365/ac4414
  * Anders et al. (2017) Anders, F., Chiappini, C., Minchev, I., et al. 2017, A&A, 600, A70, doi: 10.1051/0004-6361/201629363
  * Andrievsky et al. (2002) Andrievsky, S. M., Kovtyukh, V. V., Luck, R. E., et al. 2002, A&A, 392, 491, doi: 10.1051/0004-6361:20021035
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Balser et al. (2011) Balser, D. S., Rood, R. T., Bania, T. M., & Anderson, L. D. 2011, ApJ, 738, 27, doi: 10.1088/0004-637X/738/1/27
  * Balser et al. (2015) Balser, D. S., Wenger, T. V., Anderson, L. D., & Bania, T. M. 2015, ApJ, 806, 199, doi: 10.1088/0004-637X/806/2/199
  * Bellardini et al. (2021) Bellardini, M. A., Wetzel, A., Loebman, S. R., et al. 2021, MNRAS, 505, 4586, doi: 10.1093/mnras/stab1606
  * Bennett & Bovy (2019) Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417, doi: 10.1093/mnras/sty2813
  * Bensby et al. (2014) Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71, doi: 10.1051/0004-6361/201322631
  * Bergemann et al. (2014) Bergemann, M., Ruchti, G. R., Serenelli, A., et al. 2014, A&A, 565, A89, doi: 10.1051/0004-6361/201423456
  * Boeche et al. (2013) Boeche, C., Chiappini, C., Minchev, I., et al. 2013, A&A, 553, A19, doi: 10.1051/0004-6361/201219607
  * Boeche et al. (2014) Boeche, C., Siebert, A., Piffl, T., et al. 2014, A&A, 568, A71, doi: 10.1051/0004-6361/201423974
  * Bovy et al. (2012) Bovy, J., Rix, H.-W., & Hogg, D. W. 2012, ApJ, 751, 131, doi: 10.1088/0004-637X/751/2/131
  * Carr et al. (2022) Carr, C., Johnston, K. V., Laporte, C. F. P., & Ness, M. K. 2022, arXiv e-prints, arXiv:2201.04133.  https://arxiv.org/abs/2201.04133
  * Carraro et al. (2010) Carraro, G., Vázquez, R. A., Costa, E., Perren, G., & Moitinho, A. 2010, ApJ, 718, 683, doi: 10.1088/0004-637X/718/2/683
  * Chen et al. (2015) Chen, Y.-Q., Zhao, G., Liu, C., et al. 2015, Research in Astronomy and Astrophysics, 15, 1125, doi: 10.1088/1674-4527/15/8/003
  * Cheng et al. (2012) Cheng, J. Y., Rockosi, C. M., Morrison, H. L., et al. 2012, ApJ, 746, 149, doi: 10.1088/0004-637X/746/2/149
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765, doi: 10.1086/303726
  * Cunha et al. (2016) Cunha, K., Frinchaboy, P. M., Souto, D., et al. 2016, Astronomische Nachrichten, 337, 922, doi: 10.1002/asna.201612398
  * Czesla et al. (2019) Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages.  http://ascl.net/1906.010
  * Darragh-Ford et al. (2023) Darragh-Ford, E., Hunt, J. A. S., Price-Whelan, A. M., & Johnston, K. V. 2023, ApJ, 955, 74, doi: 10.3847/1538-4357/acf1fc
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371, doi: 10.1086/163168
  * Debattista et al. (2024) Debattista, V. P., Khachaturyants, T., Amarante, J. A. S., et al. 2024, arXiv e-prints, arXiv:2402.08356, doi: 10.48550/arXiv.2402.08356
  * Di Matteo et al. (2013) Di Matteo, P., Haywood, M., Combes, F., Semelin, B., & Snaith, O. N. 2013, A&A, 553, A102, doi: 10.1051/0004-6361/201220539
  * Eggen et al. (1962) Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748, doi: 10.1086/147433
  * Eilers et al. (2019) Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120, doi: 10.3847/1538-4357/aaf648
  * Eilers et al. (2022) Eilers, A.-C., Hogg, D. W., Rix, H.-W., et al. 2022, ApJ, 928, 23, doi: 10.3847/1538-4357/ac54ad
  * Filion et al. (2023) Filion, C., McClure, R. L., Weinberg, M. D., D’Onghia, E., & Daniel, K. J. 2023, MNRAS, 524, 276, doi: 10.1093/mnras/stad1832
  * Gaia Collaboration et al. (2022a) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022a, A&A, A2, doi: 10.1051/0004-6361/202243940
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2022b) Gaia Collaboration, Recio-Blanco, A., Kordopatis, G., et al. 2022b, arXiv e-prints, arXiv:2206.05534.  https://arxiv.org/abs/2206.05534
  * Gonzalez et al. (2011) Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2011, A&A, 530, A54, doi: 10.1051/0004-6361/201116548
  * Grand et al. (2012) Grand, R. J. J., Kawata, D., & Cropper, M. 2012, MNRAS, 421, 1529, doi: 10.1111/j.1365-2966.2012.20411.x
  * Grand et al. (2016) Grand, R. J. J., Springel, V., Kawata, D., et al. 2016, MNRAS, 460, L94, doi: 10.1093/mnrasl/slw086
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15, doi: 10.1051/0004-6361/201833718
  * Greggio & Renzini (1983) Greggio, L., & Renzini, A. 1983, A&A, 118, 217
  * Hammersley et al. (1994) Hammersley, P. L., Garzon, F., Mahoney, T., & Calbet, X. 1994, MNRAS, 269, 753, doi: 10.1093/mnras/269.3.753
  * Hawkins (2022) Hawkins, K. 2022, arXiv e-prints, arXiv:2207.04542, doi: 10.48550/arXiv.2207.04542
  * Hawkins et al. (2016) Hawkins, K., Masseron, T., Jofré, P., et al. 2016, A&A, 594, A43, doi: 10.1051/0004-6361/201628812
  * Hawkins et al. (2015) Hawkins, K., Kordopatis, G., Gilmore, G., et al. 2015, MNRAS, 447, 2046, doi: 10.1093/mnras/stu2574
  * Hayden et al. (2014) Hayden, M. R., Holtzman, J. A., Bovy, J., et al. 2014, AJ, 147, 116, doi: 10.1088/0004-6256/147/5/116
  * Hayden et al. (2015) Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132, doi: 10.1088/0004-637X/808/2/132
  * Ho et al. (2017) Ho, I. T., Seibert, M., Meidt, S. E., et al. 2017, ApJ, 846, 39, doi: 10.3847/1538-4357/aa8460
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hwang et al. (2019) Hwang, H.-C., Barrera-Ballesteros, J. K., Heckman, T. M., et al. 2019, ApJ, 872, 144, doi: 10.3847/1538-4357/aaf7a3
  * Imig et al. (2023) Imig, J., Price, C., Holtzman, J. A., et al. 2023, arXiv e-prints, arXiv:2307.13887, doi: 10.48550/arXiv.2307.13887
  * Janes (1979) Janes, K. A. 1979, ApJS, 39, 135, doi: 10.1086/190568
  * Jenkins & Binney (1990) Jenkins, A., & Binney, J. 1990, MNRAS, 245, 305
  * Jönsson et al. (2018) Jönsson, H., Allende Prieto, C., Holtzman, J. A., et al. 2018, AJ, 156, 126, doi: 10.3847/1538-3881/aad4f5
  * Khoperskov et al. (2018) Khoperskov, S., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 611, L2, doi: 10.1051/0004-6361/201732521
  * Khoperskov & Gerhard (2021) Khoperskov, S., & Gerhard, O. 2021, arXiv e-prints, arXiv:2111.15211.  https://arxiv.org/abs/2111.15211
  * Kobayashi & Nakasato (2011) Kobayashi, C., & Nakasato, N. 2011, ApJ, 729, 16, doi: 10.1088/0004-637X/729/1/16
  * Kreckel et al. (2019) Kreckel, K., Ho, I. T., Blanc, G. A., et al. 2019, ApJ, 887, 80, doi: 10.3847/1538-4357/ab5115
  * Laporte et al. (2018) Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286, doi: 10.1093/mnras/sty1574
  * Larson (1976) Larson, R. B. 1976, MNRAS, 176, 31
  * Lemasle et al. (2008) Lemasle, B., François, P., Piersimoni, A., et al. 2008, A&A, 490, 613, doi: 10.1051/0004-6361:200810192
  * Leung & Bovy (2019a) Leung, H. W., & Bovy, J. 2019a, MNRAS, 483, 3255, doi: 10.1093/mnras/sty3217
  * Leung & Bovy (2019b) —. 2019b, MNRAS, 489, 2079, doi: 10.1093/mnras/stz2245
  * Lucey et al. (2023) Lucey, M., Pearson, S., Hunt, J. A. S., et al. 2023, MNRAS, 520, 4779, doi: 10.1093/mnras/stad406
  * Luck & Lambert (2011) Luck, R. E., & Lambert, D. L. 2011, AJ, 142, 136, doi: 10.1088/0004-6256/142/4/136
  * Mackereth et al. (2019) Mackereth, J. T., Bovy, J., Leung, H. W., et al. 2019, MNRAS, 489, 176, doi: 10.1093/mnras/stz1521
  * Magrini et al. (2009) Magrini, L., Sestito, P., Randich, S., & Galli, D. 2009, A&A, 494, 95, doi: 10.1051/0004-6361:200810634
  * Majewski et al. (2017a) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017a, AJ, 154, 94, doi: 10.3847/1538-3881/aa784d
  * Majewski et al. (2017b) —. 2017b, AJ, 154, 94, doi: 10.3847/1538-3881/aa784d
  * Matteucci & Francois (1989) Matteucci, F., & Francois, P. 1989, MNRAS, 239, 885, doi: 10.1093/mnras/239.3.885
  * Matteucci & Recchi (2001) Matteucci, F., & Recchi, S. 2001, ApJ, 558, 351, doi: 10.1086/322472
  * Mayor (1976) Mayor, M. 1976, A&A, 48, 301
  * Momany et al. (2006) Momany, Y., Zaggia, S., Gilmore, G., et al. 2006, A&A, 451, 515, doi: 10.1051/0004-6361:20054081
  * Montalbán et al. (2021) Montalbán, J., Mackereth, J. T., Miglio, A., et al. 2021, Nature Astronomy, 5, 640, doi: 10.1038/s41550-021-01347-7
  * Nandakumar et al. (2020) Nandakumar, G., Hayden, M. R., Sharma, S., et al. 2020, arXiv e-prints, arXiv:2011.02783.  https://arxiv.org/abs/2011.02783
  * Önal Taş et al. (2016) Önal Taş, Ö., Bilir, S., Seabroke, G. M., et al. 2016, PASA, 33, e044, doi: 10.1017/pasa.2016.33
  * Pedicelli et al. (2009) Pedicelli, S., Bono, G., Lemasle, B., et al. 2009, A&A, 504, 81, doi: 10.1051/0004-6361/200912504
  * Poggio et al. (2021) Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, A&A, 651, A104, doi: 10.1051/0004-6361/202140687
  * Poggio et al. (2022) Poggio, E., Recio-Blanco, A., Palicio, P. A., et al. 2022, arXiv e-prints, arXiv:2206.14849.  https://arxiv.org/abs/2206.14849
  * Price-Whelan et al. (2020) Price-Whelan, A., Sipőcz, B., Lenz, D., et al. 2020, adrn/gala: v1.3, v1.3, Zenodo, doi: 10.5281/zenodo.4159870
  * Price-Whelan (2017) Price-Whelan, A. M. 2017, The Journal of Open Source Software, 2, doi: 10.21105/joss.00388
  * Ramírez et al. (2013) Ramírez, I., Allende Prieto, C., & Lambert, D. L. 2013, ApJ, 764, 78, doi: 10.1088/0004-637X/764/1/78
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131, doi: 10.3847/1538-4357/ab4a11
  * Reylé et al. (2009) Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819, doi: 10.1051/0004-6361/200811341
  * Rolleston et al. (2000) Rolleston, W. R. J., Smartt, S. J., Dufton, P. L., & Ryans, R. S. I. 2000, A&A, 363, 537
  * Schönrich & Binney (2009) Schönrich, R., & Binney, J. 2009, MNRAS, 396, 203, doi: 10.1111/j.1365-2966.2009.14750.x
  * Schönrich et al. (2010) Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829, doi: 10.1111/j.1365-2966.2010.16253.x
  * Searle & Zinn (1978) Searle, L., & Zinn, R. 1978, ApJ, 225, 357, doi: 10.1086/156499
  * Sellwood & Binney (2002) Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785, doi: 10.1046/j.1365-8711.2002.05806.x
  * Tinsley (1980) Tinsley, B. M. 1980, Fund. Cosmic Phys., 5, 287, doi: 10.48550/arXiv.2203.02041
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2019) Wang, C., Liu, X. W., Xiang, M. S., et al. 2019, MNRAS, 482, 2189, doi: 10.1093/mnras/sty2797
  * Wegg et al. (2015) Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050, doi: 10.1093/mnras/stv745
  * Woosley & Hoffman (1992) Woosley, S. E., & Hoffman, R. D. 1992, ApJ, 395, 202, doi: 10.1086/171644
  * Xiang et al. (2015) Xiang, M.-S., Liu, X.-W., Yuan, H.-B., et al. 2015, Research in Astronomy and Astrophysics, 15, 1209, doi: 10.1088/1674-4527/15/8/009
  * Yan et al. (2019) Yan, Y., Du, C., Liu, S., et al. 2019, ApJ, 880, 36, doi: 10.3847/1538-4357/ab287d