# A Blueprint for the Milky Way’s Stellar Populations. IV.
Spectrum-based Empirical Corrections and the Galactic Starburst Sequence

Deokkeun An  Department of Science Education, Ewha Womans University, 52
Ewhayeodae-gil, Seodaemun-gu, Seoul 03760, Republic of Korea;
<EMAIL_ADDRESS>Timothy C. Beers  Department of Physics and Astronomy and
JINA Center for the Evolution of the Elements, University of Notre Dame, Notre
Dame, IN 46556, USA  Young Sun Lee  Department of Physics and Astronomy and
JINA Center for the Evolution of the Elements, University of Notre Dame, Notre
Dame, IN 46556, USA  Department of Astronomy and Space Science, Chungnam
National University, Daejeon 34134, Republic of Korea  Thomas Masseron
Instituto de Astrofísica de Canarias (IAC), E-38200 La Laguna, Tenerife, Spain
Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206 La
Laguna, Tenerife, Spain

###### Abstract

We continue our series of papers on phase-space distributions of stars in the
Milky Way based on photometrically derived metallicities and Gaia astrometry,
with a focus on the halo-disk interface in the local volume. To exploit
various photometric databases, we develop a method of empirically calibrating
synthetic stellar spectra based on a comparison with observations of stellar
sequences and individual stars in SDSS, SMSS, and PS1, overcoming band-
specific corrections employed in our previous work. In addition, photometric
zero-point corrections are derived to provide an internally consistent
photometric system with a spatially uniform metallicity zero point. Based on
our phase-space diagrams, we find a strikingly narrow sequence in the
rotational velocity ($v_{\phi}$) versus metallicity ([Fe/H]) space for a
sample of high proper-motion stars ($>25$ mas yr-1), which runs along the Gaia
Sausage/Enceladus (GSE) and the Splash sub-structures, and is connected to the
disk. Notably, a rapid increase of $v_{\phi}$ from a nearly zero net rotation
to $200\ \rm km~{}s^{-1}$ within a narrow metallicity interval
($-0.8<{\rm[Fe/H]}<-0.2$) suggests that these stars were formed on a short
gas-depletion time scale. Based on measurements of a scale height and length,
we argue that they are distinct from stars dynamically heated by mergers, and
are the relics formed during the starburst when the young Milky Way
encountered the gas-rich GSE merger. The chain of high proper-motion stars,
which we dub the Galactic Starburst Sequence (GSS), provides evidence that the
post-merger metal-enriched gas settled onto the disk.

Unified Astronomy Thesaurus concepts: Milky Way Galaxy (1054); Milky Way
stellar halo (1060); Milky Way dynamics (1051); Milky Way formation (1053);
Milky Way evolution (1052); Stellar abundances (1577); Stellar populations
(1622)

††journal: the AAS Journals††facilities: AAVSO, PS1, SkyMapper,
Sloan††software: GNU parallel (Tange 2021)

##  1 Introduction

Photometric survey databases are useful resources for studying stellar
populations and structures of the Milky Way Galaxy. The size of spectroscopic
samples has grown rapidly in recent years, but photometric surveys still cover
a significantly larger volume of space and thereby can provide the least-
biased sample of Galactic stars. Multi-band observations are particularly
useful, since they can be used to constrain fundamental stellar parameters,
such as effective temperature ($T_{\rm eff}$) and metallicity ([Fe/H]), with
sufficient accuracy. The overall shape of a spectral energy distribution as
traced by multi-band photometry depends on $T_{\rm eff}$, and ultra-violet
excess provides information about a star’s metallicity. Importantly,
photometric data can be combined with all-sky, high-precision astrometry from
Gaia, providing rich information on chemical and kinematical properties of
stars, as we demonstrated in a series of papers (An & Beers, 2020, 2021a,
2021b, hereafter Paper I, II, and III, respectively).

To obtain a clear view of Galactic stellar populations, one needs to establish
accurate relationships between photometry and fundamental stellar parameters
on an empirical basis or by using theoretical predictions. To take advantage
of each method, we adopted a hybrid approach in previous papers in this series
to calibrate theoretical isochrones of the main sequence using observations of
well-studied Galactic globular and open clusters. The models were taken from
YREC (Sills et al., 2000), and were combined with MARCS (Gustafsson et al.,
2008) synthetic spectra, in order to convert $T_{\rm eff}$ and luminosities
into photometric colors and magnitudes. Differences of the models from
observations typically amount to a few hundredths of a magnitude for warm
($T_{\rm eff}>5000$ K) stars, but they become as large as a few tenths of a
magnitude for cooler stars. The model offsets are also a function of
metallicity, in that more metal-rich cluster sequences tend to exhibit larger
color deviations.

The observed offsets from the models are systematic in nature and cannot be
simply reconciled by adjusting input cluster parameters, which implies that
they may originate from errors in the input physics and physical parameters.
To overcome these difficulties, we took the observed offsets as empirical
correction factors that one needs to apply to our specific choice of
theoretical stellar models. When models are used with empirical corrections,
we obtain distances from SDSS photometry that are consistent with Gaia
parallaxes, and our photometrically derived metallicities ([Fe/H]) are also in
overall agreement with spectroscopic measurements in SDSS (Paper I).

While we developed a method of deriving photometric metallicities from SDSS,
this set of color-$T_{\rm eff}$ corrections is only valid for observations
taken in the SDSS $ugriz$ filter set. Other photometric surveys also adopt
filter sets similar to that of SDSS, but their transmission curves are not
exactly the same as each other, leading to non-negligible differences in
magnitudes. In this sense, direct calibration of synthetic spectra can serve
as an alternative way of establishing such relations in various filter
passbands. Importantly, it enables us to combine various photometric survey
databases and produce chemo-kinematical phase-space maps over the entire
celestial hemisphere in an internally consistent manner.

One of the goals in this study is to generalize our empirical-correction
procedure, and construct a set of corrected synthetic spectra in order to
generate magnitudes in any given filter set with high confidence. An essential
requirement to achieve this goal is to finely sample flux for a set of
calibration stars over a wide range of wavelength and stellar parameters using
multi-band photometry, which has become practical in the era of massive
photometric surveys. The basic idea of calibrating model fluxes in theoretical
stellar spectra, as opposed to making corrections on individual color indices
in the models, was introduced by Lejeune et al. (1997, 1998). However, the
current work is based on a significantly larger set of photometric,
spectroscopic, and astrometric data, which were unavailable then.

The other goal of this work is to probe multi-dimensional data cubes
constructed based on the revised metallicity estimates. In addition to
chemical information from photometry, we exploit kinematic data from Gaia, as
in the previous papers of this series. The majority of main-sequence stars in
our sample are too faint to have radial velocity measurements. However, along
the great circle perpendicular to the direction of disk rotation ($l=0\arcdeg$
and $180\arcdeg$), rotational velocities ($v_{\phi}$) in the cylindrical
coordinate in the rest frame of the Galaxy do not depend on radial velocity.
This enables to derive $v_{\phi}$ by utilizing the proper motions of stars
near the Galactic prime meridian and construct phase-space diagrams in
$v_{\phi}$ and [Fe/H], which can subsequently be used to characterize
kinematical and chemical properties of individual populations.

Because our method relies on calibration of theoretical models for main-
sequence stars, our current approach is effectively limited to a local volume
($d<6$ kpc). Nonetheless, in our previous papers, we demonstrated the
usefulness of such data by providing an unbiased, global perspective on local
stellar populations, including Gaia Sausage/Enceladus (GSE; Belokurov et al.,
2018; Helmi et al., 2018) and the Splash (Belokurov et al., 2020). In
particular, we applied Gaussian mixture models in Paper III to isolate
individual stellar populations, and evaluated their fractional contributions
as a function of distance from the Galactic plane ($Z$) and Galactocentric
distance ($R_{\rm GC}$). In this work, we present evidence for yet another
stellar population, which appears to have been formed during a starburst
episode in the early evolution of the Galaxy, possibly induced by the GSE
merger.

This paper is organized into two parts, which describe each of the above two
subjects: the revised calibration of theoretical models (§ 2) and its
application to photometric databases (§ 3). We summarize our results in § 4.

##  2 Spectrum-based Empirical Corrections on Theoretical Models

Stellar metallicities presented in this study are computed using a set of
theoretical stellar isochrones with revised empirical corrections. As
described below, our new color-$T_{\rm eff}$-[Fe/H] relations, which convert
theoretically predicted quantities ($T_{\rm eff}$ and luminosity) into
observables (colors and magnitudes), hinge on both stellar sequences and field
stars with spectroscopic metallicity estimates; see An et al. (2009, 2013),
Papers I and II, and references therein, for more information on our previous
model corrections. The same set of base theoretical models (Sills et al.,
2000; Gustafsson et al., 2008) with identical model parameters, including an
age-metallicity relation and $\alpha$-element abundance mixtures, are adopted
in this work. In contrast to the models adopted in the previous papers of this
series, this newer version of calibration utilizes photometry in various
filter passbands over a wide range of wavelength, and thereby enables fine-
tuning of synthetic stellar spectra.

###  2.1 Calibration Samples

Table 1: A Summary of Calibration Samples Photometry | Bandpasses | Stellar sequences | Individual stars (spectroscopic sample)
---|---|---|---
SDSS | $ugriz$ |  An et al. (2008), Gaia double | SEGUE, GALAH
SMSS | $uvgriz$ | Gaia double | SEGUE, GALAH
PS1 | $grizy$ |  Bernard et al. (2014), Gaia double | SEGUE, GALAH
APASS | $BV$ | Gaia double | SEGUE, GALAH
Stetson | $BVI_{C}$ | Stetson’s standard photometry |  $\cdots$
Table 2: Fundamental Parameters and Uncertainties of Stellar Sequences Cluster/Sequence | [Fe/H]aaUniform uncertainties except for M67 and NGC 6791. |  $(m\,-\,M)_{0}$aaUniform uncertainties except for M67 and NGC 6791. |  $E(B\,-\,V)$aaUniform uncertainties except for M67 and NGC 6791. | ageaaUniform uncertainties except for M67 and NGC 6791. | Min $M_{r}$bbA minimum $M_{r}$ to select MS stars. |
---|---|---|---|---|---|---
Name | (dex) | (mag) | (mag) | (Gyr) | (mag) | ReferencesccReferences for [Fe/H], $(m\,-\,M)_{0}$, and $E(B\,-\,V)$, respectively: (1) Kraft & Ivans (2003); (2) Carretta et al. (2000); (3) Paper II; (4) An et al. (2019) and references therein.
M15 | $-2.42\pm 0.10$ | $15.25\pm 0.15$ | $0.100\pm 0.020$ | $13.0\pm 2.6$ | $19.4$ | 1/1/1
M92 | $-2.38\pm 0.10$ | $14.64\pm 0.15$ | $0.020\pm 0.004$ | $13.0\pm 2.6$ | $18.5$ | 1/2/1
M13 | $-1.60\pm 0.10$ | $14.38\pm 0.15$ | $0.020\pm 0.004$ | $13.0\pm 2.6$ | $18.5$ | 1/2/1
M3 | $-1.50\pm 0.10$ | $15.02\pm 0.15$ | $0.010\pm 0.002$ | $13.0\pm 2.6$ | $19.0$ | 1/1/1
Gaia double (blue MS) | $-1.30\pm 0.10$ |  $\cdots$ ddParallaxes for individual field stars in Gaia EDR3, along with their foreground extinctions in Schlegel et al. (1998). |  $\cdots$ ddParallaxes for individual field stars in Gaia EDR3, along with their foreground extinctions in Schlegel et al. (1998). | $13.0\pm 2.6$ | $4.0$ | 3/./.
M5 | $-1.26\pm 0.10$ | $14.46\pm 0.15$ | $0.030\pm 0.006$ | $12.0\pm 2.4$ | $18.5$ | 1/2/1
Gaia double (red MS) | $-0.40\pm 0.10$ |  $\cdots$ ddParallaxes for individual field stars in Gaia EDR3, along with their foreground extinctions in Schlegel et al. (1998). |  $\cdots$ ddParallaxes for individual field stars in Gaia EDR3, along with their foreground extinctions in Schlegel et al. (1998). | $12.0\pm 2.4$ | $4.0$ | 3/./.
M67 | $+0.00\pm 0.01$ | $9.61\pm 0.03$ | $0.041\pm 0.004$ | $4.0\pm 0.5$ | $13.0$ | 4/4/4
NGC 6791 | $+0.37\pm 0.07$ | $13.06\pm 0.06$ | $0.120\pm 0.020$ | $9.0\pm 1.0$ | $17.6$ | 4/4/4

As summarized in Table 1, we adopt both cluster sequences and a set of
individual stars with spectroscopic metallicity estimates for the calibration
of models. As shown in the first two columns, a total $19$ passbands from
various photometric surveys are utilized in this study: $ugriz$ photometry
from the Sloan Digital Sky Survey (SDSS) DR14, $uvgriz$ from the SkyMapper Sky
Survey (SMSS) DR2 (Onken et al., 2019), $grizy$ from the Pan-STARRS1 surveys
(PS1, Chambers et al., 2016), and $BV$ from the AAVSO Photometric All-Sky
Survey (APASS) DR10 (Henden et al., 2018). The $gri$ photometry in APASS is
not used owing to large photometric zero-point offsets (see also Tonry et al.,
2018). We add to the list standard-star photometry in $BVI_{C}$ constructed by
P. Stetson (see Stetson,
2000).111https://www.canfar.net/storage/list/STETSON/Standards. Some of these
databases use the same notation for their filter passbands ($ugriz$), but
their response functions are not identical. Below we make a distinction
between these filters by specifying the survey names.

Table 2 lists the stellar sequences adopted in this study. As in our previous
exercise, we employ a set of well-studied Galactic globular and open clusters
(M15, M92, M13, M3, M5, M67, and NGC 6791) over a wide range of metallicity
($-2.4\leq{\rm[Fe/H]}\leq 0.4$). We use fiducial sequences from An et al.
(2008), which were derived from the SDSS imaging data. Zero-point corrections
(An et al., 2013) are applied to tie An et al. cluster photometry to SDSS DR14
(Abolfathi et al., 2018). Fiducial sequences in PS1 for the same set of
clusters are taken from Bernard et al. (2014). For Stetson’s photometry, we
only take data in $BVI_{C}$, since the $U$-band and $R$-band are less well-
defined than the others. At the time of this writing, cluster photometry in
the SMSS and APASS filter systems is not yet available. In future work, we
plan to improve the calibration of synthetic spectra by incorporating such
data whenever they become available.

In addition to clusters, we employ the Gaia double sequence, which appears on
a color-magnitude diagram from stars with large proper motions (Gaia
Collaboration et al., 2018a). As demonstrated in Paper II, each of the
sequences represents two dominant populations in the local halo – GSE and
Splash – and has [Fe/H]$\approx-1.3$ and $-0.4$, respectively (see also
Sahlholdt et al., 2019). This dynamically defined group of stars also provides
a powerful constraint on the shape of a sequence, bridging the gap between
globular and open clusters at intermediate metallicities. We follow the
procedure developed in Paper II to extract individual sequences from the Gaia
double sequence (see Appendix A). In short, this technique relies on
metallicity-sensitive $u$-band photometry to separate the two chemically
distinct populations. For SDSS and SMSS, we use their $u$-band data. For PS1
and APASS, we use SDSS $u$. Additional cuts on kinematics further help to
isolate each of the populations. As in our previous work, we compute
$v_{\phi}$ using Gaia’s proper motions and parallaxes, but without radial
velocity measurements, within $\pm 30\arcdeg$ along the Galactic prime
meridian. To construct a clean sequence, we use objects with good astrometry,
having less than $20\%$ uncertainty in parallax and $30\%$ uncertainty in
proper motion. We also apply cuts on $E(B\,-\,V)<0.1$ and $|b|>20\arcdeg$. As
in Paper II, we impose $-150<v_{\phi}\leq-50\ \rm km~{}s^{-1}$ and
$120<v_{\phi}\leq 150\ \rm km~{}s^{-1}$ on the sample to derive the blue
(metal-poor) and red (metal-rich) main sequences, respectively, for which we
compute weighted median colors in bins of $0.5$ mag in $M_{r}$ ($\Delta
M_{r}=0.2$ mag in $u\,-\,r$ CMDs).

For individual calibration stars with spectroscopic metallicities, we utilize
the Sloan Extension for Galactic Understanding and Exploration (SEGUE, Yanny
et al., 2009; Rockosi et al., 2022) and the Galactic Archaeology with HERMES
(GALAH) survey (Buder et al., 2021), as they are among the largest and most
uniform spectroscopic data sets in the Northern and Southern Hemisphere,
respectively. Specifically, we adopt metallicity estimates from a rerun of the
updated SEGUE Stellar Parameter Pipeline (SSPP; Lee et al., 2008a, b),
performed by one of the coauthors (Y. S. Lee). We apply cuts based on a
signal-to-noise ratio (SNR) of the spectra ($>30$), uncertainty in $T_{\rm
eff}<200$ K, and $\log{g}>3.5$ to select main-sequence stars with high-quality
parameter estimates. For the GALAH sample, we also require that a stellar
parameter quality flag (flag_sp) and an overall iron abundance quality flag
(flag_fe_h) are not set.

The advantage of the SEGUE and GALAH samples is corroborated by the
availability of photometric data in various passbands (such as SMSS $uv$; see
Table 1), while a stable baseline provided by the cluster sample makes the
sequence-based approach reliable. All SEGUE stars are covered by SDSS and PS1
imaging surveys, but can only be matched to objects in SMSS near the
equatorial region. Likewise, the majority of GALAH stars have good matches to
SMSS photometry, but only a small fraction of its survey area overlaps with
SDSS. Primary stellar sources (type$=6$) in the SDSS are kept, with a set of
minimal quality flags in the $r$-band measurements to ensure that sources do
not have issues such as de-blending, interpolation, and saturation. Similarly,
primary detections in PS1 from its stacked imaging catalog are taken. We
select point-like sources by imposing a maximum $0.05$ mag difference in
$i$-band photometry between a point-spread function and Kron magnitudes. For
SMSS, we apply cuts on a number of photometric quality flags to only retain
good photometric measurements: class_star $>0.9$, flags $<3$, nch_max $=1$,
prox $>7.5$, ngood_min $>1$, and nimaflags $=0$ in each passband.

Photometric catalog objects are matched with Gaia Early Data Release 3 (EDR3,
Gaia Collaboration et al., 2021) and Data Release 3 (DR3, Gaia Collaboration
et al., 2022) using a $1\arcsec$ search radius. Zero-point corrections on
parallax by the Gaia team (Lindegren et al., 2021) are adopted. In the
following calibration, objects with good parallaxes ($\sigma_{\pi}/\pi<0.2$)
are used; about $68\ \%$ of them have parallaxes within $2\%$ from those
inferred based on Bayes’ theorem (Bailer-Jones et al., 2021). The foreground
reddening values in Schlegel et al. (1998) are adopted, except for clusters,
along with extinction coefficients at $R_{V}=3.1$ in Schlafly & Finkbeiner
(2011) for SDSS, PS1, and Johnson-Cousins (those listed as ‘Landolt’) bands.
For SMSS, we adopt values in Wolf et al. (2018). These coefficients assume
$14\%$ reduction in the original $E(B\,-\,V)$ in Schlegel et al. (1998).

###  2.2 Spectrum-based Corrections

####  2.2.1 Scope

Because the overall shape of a spectrum is influenced by $T_{\rm eff}$ and
metallicity, and by $\log{g}$ to a lesser degree, we aim to provide an
empirically calibrated set of synthetic stellar spectra over a wide range of
$T_{\rm eff}$ and [Fe/H], which in turn can be employed to derive such
quantities in other stars. For this purpose, theoretical flux ratios between
different filter passbands (color indices) are computed as a function of
wavelength, and any deviation from the calibration data are attributed to
systematic errors in the models. All flux ratios are referenced to the SDSS or
SMSS $r$-band among different filter passbands, because of the large amount of
flux collected in this passband, and its relatively weak metallicity
sensitivity in the $T_{\rm eff}$ regime considered in this study ($4000<T_{\rm
eff}<7000$ K). Its bolometric corrections are also less prone to systematic
errors.

We derive synthetic stellar colors from MARCS model atmospheres (Gustafsson et
al., 2008); see An et al. (2009) for more information on the construction of
the MARCS model library. Filter-response curves in SDSS $ugriz$ are found on
the project webpage.222https://www.sdss.org/instruments/camera/#Filters
References for other filter transmissions are Tonry et al. (2012, PS1
$grizy$), Bessell et al. (2011, SMSS $uvgriz$), and Bessell & Murphy (2012,
Johnson-Cousins $BVI_{C}$). APASS $BV$ are taken from the SVO filter profile
service.333http://svo2.cab.inta-csic.es/svo/theory/fps3 When deriving a Vega
magnitude, we use a Vega model from the HST CALSPEC library (Bohlin et al.,
2014), with the suggested flux re-scaling from Riello et al. (2021). Effective
wavelengths ($\lambda_{\rm eff}$) of the filter passbands are computed for
each model, as they have a mild dependence on the underlying stellar spectrum
(primarily on $T_{\rm eff}$).

Isochrones generated using YREC (Sills et al., 2000) are combined with MARCS
synthetic colors. For a given [Fe/H], differences in colors between
observational data and models are computed as a function of absolute magnitude
in the $r$-band ($M_{r}$). We use $M_{r}$ to infer $T_{\rm eff}$ from a
$M_{r}$-$T_{\rm eff}$ relation in the isochrones, because spectroscopic
$T_{\rm eff}$ measurements are not typically available for individual stars in
a stellar sequence. The infrared flux method (IRFM) is another useful way of
deriving $T_{\rm eff}$ from photometry, but our approach makes a specific
prediction on both $T_{\rm eff}$ and $\log{g}$ of a star, which are necessary
ingredients for generating precise synthetic model colors. There also exist
$T_{\rm eff}$ scale differences of a few hundred kelvin between photometric
and spectroscopic approaches (e.g., Pinsonneault et al., 2004). For this
reason, instead of taking spectroscopic $T_{\rm eff}$, we compute $T_{\rm
eff}$ photometrically using stellar isochrones in the same manner as for the
stellar sequences. This significantly reduces systematic differences in
color-$T_{\rm eff}$ relations from our heterogeneous data sets, and makes
color-$T_{\rm eff}$-[Fe/H] relations internally more consistent.

Figure 1: The $T_{\rm eff}$-[Fe/H] space covered by the calibration sample.
The orange solid rectangles exhibit main sequences of clusters and the Gaia
double sequence with $u$-band measurements in SDSS, while more extended
sequences with deeper $gr$ photometry are shown by orange-dashed lines. The
number density of SEGUE and GALAH stars is shown by the gray histogram and
blue contours, respectively. Only those having SNR $>30$, $\log{g}>3.5$, and
photometric $T_{\rm eff}$ errors less than $200$ K are retained for these
individual spectroscopic metallicity measurements.

Figure 1 displays a $T_{\rm eff}$-[Fe/H] space covered by the stellar
sequences and spectroscopic targets in our sample, where the $T_{\rm eff}$ are
determined from our isochrones. Stellar sequences in SDSS are shown by orange
boxes with a width of $\pm 0.1$ dex in [Fe/H]. The dotted lines represent the
same sequences, but without valid $u$-band measurements, demonstrating the
necessity of deeper $u$-band photometry. Stellar sequences have discrete
metallicities, while the spectroscopic sample (the gray-scale histogram for
SEGUE and contours for GALAH targets) fills up the remaining space. Our
spectroscopic sample is heavily weighted by metal-rich stars, with a
significantly lower number in the metal-poor regime ([Fe/H]$<-1$). The upper
right side is not covered by both samples, owing to the increased metal
content and relatively old ages of the stars. The lower left part is
depopulated because of the lack of faint metal-poor main-sequence stars in the
spectroscopic surveys.

####  2.2.2 Model Comparisons

For each star in the spectroscopic sample, an isochrone is generated by
interpolating the model grid at the star’s metallicity ([Fe/H]), and
differences in flux are computed in various filter passbands, as done for the
stellar sequences. However, to evenly sample cluster sequences and the
spectroscopic targets, we bin each of the SEGUE and GALAH samples in [Fe/H]
and compute mean flux offsets as a function of $T_{\rm eff}$ (see Appendix A).
In accordance with the metallicity of the observed stellar sequences (Table
2), the central metallicities are set to [Fe/H]$=\\{-2.4$, $-1.6$, $-1.3$,
$-0.4$, $0.0$, $+0.3\\}$. To match this binning and keep the average
metallicity of a subset of stars as close as possible to the central
metallicity values, each spectroscopic sample is divided into
[Fe/H]=$\\{(-2.9,-2.1)$, $(-1.8,-1.4)$, $(-1.5,-1.1)$, $(-0.5,-0.3)$,
$(-0.1,0.1)$, $(+0.2,+0.4)\\}$. A large width ($0.8$ dex) is set in the lowest
metallicity bin to compensate for the small number of stars in the sample.
Because the flux difference from our models changes mildly with metallicity,
our adopted bin sizes have little impact on the calibration. Nonetheless, mean
flux offsets from GALAH are taken only at [Fe/H]$=\\{-0.4$, $0.0$, $+0.3\\}$
due to the lack of metal-poor stars in the sample.

Figure 2: Magnitude differences between observed data and isochrones as a
function of $\lambda_{\rm eff}$ of filter passbands. A $T_{\rm eff}=5800$ K
case is shown above to display the overall quality of a comparison based on
stellar sequences (globular, open clusters, and Gaia’s double sequences; open
symbols) and spectroscopic measurements of individual stars from SEGUE and
GALAH (filled symbols). Comparisons are made in six metallicity bins at
[Fe/H]$=-2.4$, $-1.6$, $-1.3$, $-0.4$, $0.0$, and $0.3$, respectively (from
top left to bottom right panels). The red solid line represents moving
averages of the difference as a function of $\lambda_{\rm eff}$ (see text).

Figure 3: Same as in Figure 2, but at $T_{\rm eff}=4500$ K. Comparisons in the
two lowest metallicity bins are not shown due to the lack of observational
data in this $T_{\rm eff}$ bin.

Figures 2 and 3 show color differences in magnitude between models and
observations, as a function of $\lambda_{\rm eff}$, at two selected $T_{\rm
eff}$ ($5800$ K and $4500$ K, respectively). In each panel, differences are
shown for the sequences and an ensemble of spectroscopic targets in open and
filled symbols, respectively, and different symbols are used to indicate
references for photometry. All flux differences are registered to the SDSS
$r$-band, or the SMSS $r$-band for the GALAH sample, owing to a relatively
small number of cross-matches with SDSS objects in the Southern sky. The
differences from the model at these two $r$ passbands are nearly the same, and
are defined to be zero (i.e., their bolometric corrections are assumed to be
correct).444Similarly, because cluster fiducial sequences in PS1 (Bernard et
al., 2014) are not directly tied to SDSS photometry, color indices are
registered using the PS1 $r$-band, instead of the fiducial SDSS $r$-band.
Nevertheless, differences from model colors in the PS1 $r$-band are negligible
(see Appendix A). Thus, the effect of switching between SDSS $r$-band and PS1
$r$-band is not significant.

Error bars for the stellar sequences in Figures 2 and 3 represent propagated
uncertainties from photometry and input parameters ([Fe/H], $(m\,-\,M)_{0}$,
$E(B\,-\,V)$, and age; Table 2). For the fiducial sequences in SDSS and PS1,
constant uncertainties of $0.02$ mag are assumed in the color indices.
Similarly, a $2\%$ error is adopted for the mean colors of Stetson’s $BVI_{C}$
cluster sequences as a conservative limit. Photometry of Gaia’s double
sequence is collected across a large area on the sky, and therefore an
observed scatter is taken as uncertainties in the mean colors, since it
represents a sum of random and systematic zero-point errors, unless propagated
photometric uncertainties are larger. A comparison of APASS photometry to
Stetson’s standard photometry for the sample clusters (except NGC 6791)
reveals an rms dispersion of $0.04$ mag in $BV$. Thus, it is added in
quadrature to a total error budget for the APASS-based double sequence.

For the spectroscopic sample, the error bars in Figures 2 and 3 also indicate
the quadratic sum of random and systematic uncertainties. The random component
includes uncertainties in photometry, spectroscopic [Fe/H], and Gaia parallax.
For the systematic uncertainty, a $0.1$ dex in [Fe/H] is assumed to take into
account a scale difference between our models and the spectroscopic
determinations. A $20\%$ uncertainty in age is adopted for all stars. Flux
differences from these systematic uncertainties are typically less than $0.01$
mag in $gizyBV$, but are as high as $0.05$ mag in $uv$ at high metallicities.
In APASS, $0.04$ mag uncertainty in photometry is further incorporated into
the final uncertainty (see above).

In Figures 2 and 3, there are systematic differences between the two classes
of samples. For instance, the stellar sequences tend to show smaller flux
deviations from the models than the spectroscopic sample; at [Fe/H]$\geq-0.4$,
our result indicates that the spectroscopic sample is fainter at
$\lambda<4000$ Å, but brighter at $\lambda>7000$ Å than the sequences. Indeed,
the cluster sequences in PS1 $grizy$ (Bernard et al., 2014) exhibit the
smallest differences overall. Such a difference can be caused by inconsistent
metallicity scales, but other sources of errors, such as adopted ages, may
also contribute to the observed offset, although the absolute model deviation
changes monotonically with age, without modifying the observed wavelength-
dependent offsets.

In Figures 2 and 3, the red line indicates an average model deviation as a
function of $\lambda_{\rm eff}$. Mean magnitude differences between models and
observational data are computed by linearly interpolating values at three
adjacent filter $\lambda_{\rm eff}$. They are smoothed by applying a boxcar
average with a width of $1000$ Å, which is comparable to the FWHM of a broad
filter passband. Average differences indicate that our models greatly over-
estimate flux below $5000$ Å, by up to $20\%$, while they under-estimate flux
at longer wavelengths, by $10\%$ at the most. The $BV$ photometry from APASS
(open and filled downward triangles) exhibits consistently larger fluxes than
the other calibration sample by $0.05$–$0.1$ mag, but the differences from the
mean line are within our estimated $1\sigma$ uncertainties. We suspect all-sky
photometric zero-point errors (at roughly $4\%$ levels) as a likely source of
the systematics (see above).

####  2.2.3 Construction of a Correction Cube

Figure 4: Mean magnitude differences at $\lambda=4000$ Å, as a function of
[Fe/H], at selected $T_{\rm eff}$. Each point represents one of the
metallicity groups in this work. The blue line represents the best-fitting
$3$rd-order polynomial. Constant offsets are assumed beyond the metallicity
range covered by the data.

Figure 5: Mean magnitude offsets of the original models from observational
data in [Fe/H] versus $T_{\rm eff}$ at selected wavelengths.

Figure 4 displays model differences as a function of [Fe/H]. Each data point
represents one of the metallicity groups in this work. As shown by the blue
lines, a $3^{\rm rd}$-order polynomial function is used to depict the observed
trend. Beyond the metallicity range covered by the sample, a constant offset
is assumed in the model deviation. Figure 5 displays slices of these mean
model deviations at some selected wavelengths. Red colors indicate that models
over-predict the flux, while the blue colors show regions with under-predicted
flux. In this way, we construct a three-dimensional data cube of model
deviations as a function of $T_{\rm eff}$, [Fe/H], and wavelength.

Figure 6: Same as in Figure 5, but assuming $400$ K cooler $T_{\rm eff}$ in
the YREC isochrones.

To first order, the model deviations change monotonically with wavelength:
models are brighter than observations at shorter wavelengths, while the sign
reverses at longer wavelengths. This suggests an offset in the $T_{\rm eff}$
scale of the isochrones as a major source of the systematic mismatch. However,
the required amount of offset must be very large, by about $400$ K, even for
warm stars (Figure 6); it is even larger for cooler stars. Even if there are
systematic differences between different approaches of determining $T_{\rm
eff}$, such as the IRFM and spectroscopic determination from
excitation/ionization balance, this is beyond the accepted range of errors in
the models. Therefore, it seems that the observed offsets originate from a
combination of various sources of systematic errors, such as incorrect input
physics or under-estimated line absorption in the models. Boundary conditions
in stellar-interior models may also be incompatible with the atmosphere models
in this study. On the observational side, an inconsistent metallicity scale,
incorrect assumptions on elemental-abundance ratios, or errors in the assumed
age could be responsible for the systematic $T_{\rm eff}$ offsets.

The model differences are highly non-linear in the [Fe/H] versus $T_{\rm eff}$
space, and no simple correction of the theoretical isochrones is possible to
remedy the problem. For this reason, we take the observed flux offsets in
Figure 5 as a correction matrix for our choice of stellar isochrones and
synthetic spectra. More specifically, we employ a semi-empirical approach to
correct synthetic spectra based on observations, while keeping stellar-
interior models intact. Unlike in our previous work, the data cube in Figure 5
provides a continuous function of magnitude correction in wavelength. Thus,
the corrected synthetic spectra can be applied to any filter sets in the
wavelength range covered by our calibration sample.

####  2.2.4 Comparison with Previous Calibration

In comparison to purely theoretical models, the net result of our empirical
corrections is redder colors or a higher photometric $T_{\rm eff}$, due to
over-estimation of the model flux at shorter wavelengths and under-estimation
at longer wavelengths. Apart from this fundamental change in the models, the
revised calibration also differs from our earlier versions of the empirical
corrections. The biggest change is the inclusion of individual spectroscopic
targets in the sample, which inevitably modifies the metallicity scale of the
models. As shown earlier in Paper I, isochrones calibrated using fiducial
clusters produce photometric metallicities that are in agreement with
spectroscopic estimates in SEGUE within $\Delta{\rm[Fe/H]}\sim 0.1$ dex at
[Fe/H]$>-1.5$, but the difference amounts to $\sim 0.4$ dex at [Fe/H]$=-2$, in
the sense of a lower metallicity from our cluster-based approach. Such a
difference is a direct consequence of a systematic offset in the metallicity
scale between the cluster sequences and the SEGUE stars.

Figure 7: Comparison between photometric metallicity estimates using
calibrated models in Paper III (‘seq’: sequence-based) and this work (‘comb’:
combined from both cluster sequences and spectroscopic samples). A logarithmic
number density is shown with arbitrary contour levels. The open diamonds show
a weighted median difference in bins of $\Delta M_{r}=0.2$ and
$\Delta{\rm[Fe/H]}=0.2$ dex, respectively, with error bars indicating its
uncertainty. Only those having $\sigma_{\pi}/\pi<0.1$ and $E(B\,-\,V)<0.1$ are
included in the above comparisons. The comparison in the bottom panel is
further restricted to $4.5<M_{r}<7.5$.

Figure 7 compares the photometric metallicities from Paper III with those in
this work, based on the revised calibration. The comparisons are shown for the
solutions based on Gaia parallaxes using SDSS $ugriz$ photometry. Only stars
having distance uncertainties less than $10\%$ and $E(B\,-\,V)<0.1$ are
included. Other constraints are the same as for the main sample, as described
in the next section, except in the top panel, which displays a metallicity
difference over the full range of $M_{r}$. The large deviations for bright
stars are evident, which are caused by the lack of hot stars in the
calibration sample (Figure 1). In our subsequent analysis, including the
bottom panel of Figure 7, we adopt $4.5<M_{r}<7.5$ to avoid regions with
potentially large calibration errors.

The weighted median difference in Figure 7 indicates that the two calibration
versions agree at high metallicity ([Fe/H]$>-1$), but that photometric
metallicities from this work become larger for metal-poor stars, amounting to
$0.35$ dex at [Fe/H]$=-2$. Because the GALAH sample is confined to metal-rich
stars in our calibration, the systematic trend highlights an inconsistent
match of the metallicity scales between the cluster sequences and SEGUE stars.
The SSPP estimates have been checked thoroughly using clusters and high-
resolution spectroscopic abundance determinations, and the overall agreement
is impressive (Rockosi et al., 2022). Nonetheless, such comparisons are biased
to giants and main-sequence turn-off stars in the metal-poor regime, due to
the lack of metal-poor main-sequence dwarfs in the SEGUE sample with a
sufficiently high SNR. Therefore, the difference may originate from an
internally inconsistent metallicity scale between dwarfs and giants in the
SEGUE sample.

In summary, because neither metallicity scale is preferred over the other, the
revised isochrones obtained in the current experiment should be taken as an
alternative to our earlier cluster-based calibration. More precisely, a
metallicity distribution of metal-poor stars obtained in this study based on
the hybrid sample is hinged on an intermediate metallicity scale between SEGUE
and the cluster-based work, set by their relative weights to the final
calibration sample. The sense is that it is skewed toward more metal-rich
stars below [Fe/H]$=-1$ compared to previous work of this series.

##  3 Chemo-kinematical Properties of the Local Halo

In this section, we apply our newly calibrated set of models to large
photometric catalogs, and provide new insights for Galactic stellar
populations in the local volume. In addition to chemical information from
photometry, we exploit kinematical data from Gaia to generate phase-space maps
at various distances from the Galactic plane and Galactocentric distances. As
in our previous papers, we restrict our analysis to a strip within $\pm
30\arcdeg$ from the Galactic prime meridian ($l=0\arcdeg$ and $180\arcdeg$),
where a conversion from transverse motions into $v_{\phi}$ is reliable (see
Paper III for more details). Below, we first inspect metallicity distributions
and phase-space diagrams to validate our new calibration (§ 3.1), and present
distributions of scale heights and lengths for each group of stars in bins of
$v_{\phi}$ and [Fe/H] (§ 3.2). Based on phase-space diagrams of high proper
motion stars, we demonstrate that our data reveal yet another stellar
population formed during a period of Galactic starburst activity (§ 3.3).

###  3.1 Validation of Photometric Metallicity Estimates

In the following applications, we use Gaia EDR3 as a main source catalog,
since $v_{\phi}$ in the rest frame of the Galaxy is a major ingredient of our
phase-space diagrams. We combine Gaia EDR3 astrometric data with photometry in
SDSS, SMSS, and PS1 using a $1\arcsec$ match radius. As there is a little
overlap between SDSS and SMSS, two combined catalogs — SDSS $\cap$ Gaia and
SMSS $\cap$ Gaia, respectively — are created, in addition to a master
photometric catalog from all survey data. PS1 $grizy$ photometry is added to
each data set, when there exist SDSS or SMSS photometry, to better constrain
the stellar parameters.

For each set of catalogs, we estimate metallicities for individual stars
through a grid search of calibrated models, imposing Gaia parallaxes in this
computation. For the uncertainty in metallicity, we propagate parallax
uncertainties by assuming $1\sigma$ lower/higher parallaxes and taking
differences in the metallicity estimation. The final uncertainties in
metallicity are a quadratic sum of this propagated uncertainty and the one
from a $\Delta\chi^{2}=2.3$ boundary for two degrees of freedom. Only
solutions showing at minimum fair model matches (judged from a reduced
$\chi^{2}$ of the fit less than $3$ and detected in at least five photometric
passbands) are kept. See Paper I for detailed information on model comparisons
and quality cuts on the sample.

Figure 8: A global metallicity map using the Mollweide projection in the
Galactic coordinate system. Weighted mean metallicities are shown for stars at
$1.2<d<3$ kpc from (a) SDSS, (b) SMSS, and (c) SDSS $\cup$ SMSS, where PS1
photometry is used if available. All metallicity estimates are derived using
Gaia parallaxes. Only stars with $\sigma_{\pi}/\pi<0.3$ and $|b|>20\arcdeg$
are shown, while high-extinction regions with $E(B\,-\,V)>0.1$ are excluded.
Each HEALPix area is $3.3\ {\rm deg}^{2}$ at $N_{\rm side}=32$.

Figure 8 shows metallicity maps of stars in the local volume ($1.2<d<3$ kpc)
in the Galactic coordinate system (Mollweide projection), based on our
metallicity estimates for individual stars in each of the three combined
catalogs. Panel (a) displays a map from SDSS $\cap$ Gaia, which mostly covers
the Northern Galactic Hemisphere, while panel (b) shows the Southern
Hemisphere from SMSS $\cap$ Gaia. Panel (c) displays a full coverage map from
SDSS, SMSS, and PS1. All three maps are smoothed using a median filter with a
$2\arcdeg$ radius. There are some overlaps between SDSS $\cap$ Gaia and SMSS
$\cap$ Gaia, mainly along the celestial equator, but they are not significant
overall. For the same reason, a weak discontinuity appears in the overlapping
areas. Low latitude regions ($|b|<20\arcdeg$) and some areas are excluded to
avoid large cumulative extinction, $E(B\,-\,V)>0.1$, which would otherwise
affect our photometric parameter estimates.

Reassuringly, Figure 8 reveals that more metal-rich stars are found near the
Galactic plane, from which a mean metallicity gradient is evident from the low
to high Galactic latitudes, as expected from a simple population gradient.
This exercise proves not only that our technique can be used to determine
metallicities of stars precisely from multi-band photometry, but also that we
can use corrected synthetic spectra to combine data in various filter sets to
generate an internally consistent all-sky metallicity maps. Our calibration
procedure is currently valid for main-sequence stars, and has lower precision
for giant stars. Consequently, the above mapping in this study based on main-
sequence stars probes a local volume out to $\sim 6$ kpc from the Sun. Giant
stars are excluded in our sample using color-magnitude relations based on Gaia
parallaxes, but a purely photometric approach can also be employed to tag such
stars, as demonstrated in Paper III.

Our photometric technique is a sensitive probe of photometric zero-point
errors. It is particularly useful for large photometric surveys because it is
non-trivial to have an internally consistent photometric zero point across
large areas on the sky. In Appendix B, we demonstrate the existence of
spatially correlated photometric zero-point errors in SDSS $u$ and SMSS $uv$
based on our corrected models. The size of photometric zero-point offsets is a
few hundredth of magnitude level, but is as high as $0.1$ mag in some areas.
By inverting the problem, zero-point offsets in photometry can be derived to
make a uniform mean metallicity of nearby stars on the sky. This backward
design on photometric zero-point corrections improves the quality of the
metallicity mapping and somewhat narrows the gap in model deviations between
the cluster sequences and the SEGUE sample. For this reason, we iterate the
calibration procedure (§ 2) using zero-point corrected photometry in SDSS $u$
and SMSS $uv$. In the following analysis, including Figure 8, all input
photometry is corrected for the spatially correlated zero-point offsets. These
steps closely parallel similar exercises in zero-point corrections for the
SMSS $uv$ in Huang et al. (2021, 2022).

Figure 9: Phase-space diagrams based on photometric metallicity and distance
estimates from SDSS $\cap$ PS1 photometry. (a) Stars at $3.5<|Z|\leq 6$ kpc
and $8<R_{\rm GC}\leq 11$ kpc. (b) A full sample of stars at the same $|Z|$ as
in panel (a), but at $11<R_{\rm GC}\leq 14$ kpc. (c) Same as in panel (b), but
from stars having proper motions less than $10$ mas yr-1. (d) Same as in panel
(b), but for stars having $>10$ mas yr-1.

More quantitative comparisons with previous studies can be made using phase-
space diagrams such as shown in Figure 9, in which the number density of stars
is displayed as a function of $v_{\phi}$ and [Fe/H]. Positive values of
$v_{\phi}$ indicate that stars are moving in the same direction as the
Galactic disk. The displayed data are taken at $3.5<|Z|\leq 6$ kpc, where the
contribution of disk stars can be minimized (see Paper III). Since Gaia
parallaxes have large uncertainties for our main-sequence star sample at this
distance, we use metallicity and distance estimates from fully photometric
solutions based on SDSS $\cap$ PS1 photometry, without relying on parallaxes.

The phase-space diagrams in panels (a) and (b) of Figure 9 are taken from
$8<R_{\rm GC}\leq 11$ kpc and $11<R_{\rm GC}\leq 14$ kpc, respectively, in
Galactocentric spherical coordinates, which reveal complex sub-structures of
local stars. The high metallicity, fast-rotating clump at $\langle
v_{\phi}\rangle\approx 150\ \rm km~{}s^{-1}$ and
$\langle{\rm[Fe/H]}\rangle\approx-0.5$ in panel (b) represents thick-disk
stars. Its mean $v_{\phi}$ and [Fe/H] are similar to earlier results in the
literature (see Yan et al., 2019, and references therein). On the other hand,
metal-rich ([Fe/H]$>-1$) stars in panel (a) have a skewed $v_{\phi}$
distribution, owing to an increased contribution of the Splash in the inner
Galactic region, as demonstrated in Paper III.

In panels (a) and (b) of Figure 9, a group of metal-poor stars is seen at
[Fe/H] $\approx-2.2$, which exhibits a slow net prograde rotation ($\langle
v_{\phi}\rangle\approx 70\ \rm km~{}s^{-1}$). Along with a more metal-rich
([Fe/H] $\approx-1.5$) counterpart in the inner Galactic region, it was seen
as one of the main constituents of the Galactic halo in our series of papers.
In Paper I, we considered both of these components as “inner and outer halos”
in the dual halo paradigm (Carollo et al., 2007, 2010; Beers et al., 2012),
but we called them “metal-poor and metal-rich halos” in Paper II. On the other
hand, we made a presumption in Paper III that they constitute a main body of
GSE, owing to the lack of other analogous structures known in the same phase
space.

However, recent evidence suggests that the slowly rotating metal-poor stars
are likely a separate entity from GSE. Belokurov & Kravtsov (2022) used
aluminium abundances to separate in situ stars in the halo from accreted
stars, and found a large number of metal-poor in situ stars ([Fe/H] $<-1$).
Although they span an extreme range in $v_{\phi}$ (from nearly $\sim-150\ \rm
km~{}s^{-1}$ to $\sim 300\ \rm km~{}s^{-1}$) and their spectroscopic sample is
limited to [Fe/H] $>-1.5$, their approximate mean $v_{\phi}\sim 100\ \rm
km~{}s^{-1}$ and metallicity ([Fe/H] $<-1$) suggest that these stars (dubbed
“Aurora”) are a part of the structure seen in our previous work. Given its
lower metallicity than the Splash, it is likely an old (primordial) in situ
halo that formed before the GSE merger at $z=1$–$2$.

In support of this view, we employ a simple proper-motion cut in the sample to
separate stars in GSE from the metal-poor halo distribution, as demonstrated
in panels (c) and (d) of Figure 9. They show the same phase-space diagrams of
stars as in panel (b), but have different ranges of proper motion, $<10$ mas
yr-1 and $>10$ mas yr-1, respectively. In panel (c), both the metal-poor in
situ halo and thick disk are seen, connected by a narrow band of stars, which
we assigned to the Metal-Weak Thick Disk (MWTD) in a Gaussian mixture model in
Paper III. On the other hand, an elongated structure along the $v_{\phi}=0\
\rm km~{}s^{-1}$ line stands out from the high proper-motion sample in panel
(d), which encompasses a wide range of metallicity ($-3<{\rm[Fe/H]}<-0.5$).
The chemical and kinematical properties of these high proper-motion stars are
similar to GSE in the original works of Belokurov et al. (2018) and Helmi et
al. (2018).

Although our $v_{\phi}$ measurement is strongly correlated with proper
motions, the above separation can be understood by the highly radial orbits of
GSE stars, which contrast with a nearly isotropic velocity distribution of
metal-poor in situ halo stars. We also note that the separation does not
change appreciably even if the heliocentric distances of the sample are
further narrowed down, indicating that smaller (negative) $v_{\phi}$ is not
merely caused by systematically shorter distances. In summary, all of these
chemical and kinematical properties of individual populations in our phase-
space diagrams are consistent with those found from previous (mostly
spectroscopic) studies, which essentially validates our photometric [Fe/H] and
$v_{\phi}$ estimates.

###  3.2 Scale-Height and Scale-Length Distributions

To study structural properties of each group of stars, we compute a scale
length and height in bins of $v_{\phi}$ and [Fe/H], with bin sizes set to
$25~{}\rm km~{}s^{-1}$ and $0.2$ dex, respectively. In each bin, a number
density of stars is computed in bins of $500$ pc in $Z$ and $R$, respectively,
and is fit using the following exponential function:

$\rho=\rho_{0}e^{-R/L}e^{-Z/H},$ (1)

where $L$ and $H$ represent a scale length and a scale height, respectively,
and $\rho_{0}$ is a normalization. Although it has been known that bulk halo
stars on average follow a power-law density profile, we use equation (1) as a
proxy in the local volume near the Galactic plane in order to carry a relative
comparison between different populations. Given the cone-shaped survey data
volume, a scale height is computed using stars at $1<|Z|<4$ kpc within $\pm 1$
kpc from the solar radius in each Galactic hemisphere. A scale length is
derived from stars across $7<R<12$ kpc, but tailored at $1<|Z|<3$ kpc to
include stars at nearly equal $Z$ over the $R$ range. Uncertainties in stellar
density are computed by taking a standard deviation from the best-fitting
model, and the regression is repeated to estimate final uncertainties in the
scale measurement.

Figure 10: Distribution of scale height based on the SDSS $\cap$ Gaia sample.
Results from the Northern and Southern Hemispheres are shown in the top and
bottom panels, respectively. Scale heights are computed within a $2$ kpc-wide
zone centered at the Sun ($7.34\leq R\leq 9.34$ kpc). Notable features seen in
Panel (a) are marked by circled numbers in both panels: \raisebox{-.9pt} {1}⃝
GSE, \raisebox{-.9pt} {2}⃝ Splash, and \raisebox{-.9pt} {3}⃝ metal-poor in
situ halo. Figure 11: Distribution of scale length based on the SDSS $\cap$
Gaia sample in the Northern (top) and Southern Hemisphere (bottom). Scale
lengths are computed using stars at $7<R<12$ kpc and $1<|Z|<3$ kpc. The
circled numbers are the same as in Figure 10.

Figure 10 displays a distribution of scale height in each Galactic hemisphere
from the SDSS $\cap$ Gaia sample. Only those measurements with a fractional
distance uncertainty less than $30\%$ are included. A $0.2$ dex $\times\ 25\
\rm km~{}s^{-1}$ pixel is subdivided into $4\times 4$ sub-pixels, which are
smoothed using a $5$-point boxcar average, in order to have a smoothed, global
look at changes in the structural properties. Likewise, Figure 11 shows a
distribution of scale length in both hemispheres, binned and smoothed in the
same manner as in Figure 10. A requirement that stars lie in the stripe along
the Galactic prime meridian severely limits the number of stars available in
each data set, resulting in $\approx 4\times 10^{6}$ and $2\times 10^{6}$
stars in the Northern and Southern Galactic Hemispheres, respectively.

In the top panels of Figure 10 and 11, three notable features are seen in the
Northern Galactic Hemisphere, which are marked by the circled numbers. The
intermediate scale-height and scale-length valley is particularly evident,
extending from a region mainly occupied by GSE stars (‘\raisebox{-.9pt} {1}⃝’)
to a region populated by the Splash (‘\raisebox{-.9pt} {2}⃝’). This valley
appears even more dramatic, as it contrasts with a scale-height “highland” at
[Fe/H] $\approx-1.8$ and $v_{\phi}\approx+80\ \rm km~{}s^{-1}$
(‘\raisebox{-.9pt} {3}⃝’). For reference, the circled numbers are also marked
in the following figures, including the bottom panels of Figures 10 and 11.

The intermediate scale-height valley directly shows that GSE stars are
distributed farther from the Galactic plane than disk stars, but not as far as
metal-poor in situ halo stars. This observation can be understood by a low
inclination, highly radial orbit of the GSE progenitor, which penetrated deep
into the primordial Galaxy. In addition, a comparison between the Northern and
Southern Hemispheres in Figures 10 and 11 clearly demonstrates a larger amount
of debris in the Northern Hemisphere in the local volume, owing to a pileup of
stars at the last apocenter of a highly eccentric orbit of the GSE progenitor
(Naidu et al., 2020).

Figure 12: Distribution of scale height (top) and scale length (bottom) based
on the SMSS $\cap$ Gaia sample. The circled numbers are the same as in Figure
10.

The asymmetric stellar distribution in the halo between the Northern and
Southern Galactic Hemispheres can be checked using SMSS $\cap$ Gaia. SDSS
mainly covers the Northern Galactic Hemisphere, but has a limited coverage in
the South. On the other hand, SMSS covers almost the opposite side of the
celestial hemisphere (see Figure 8) with a total of $\approx 6\times 10^{5}$
stars in the Northern Galactic Hemisphere and $\approx 10^{6}$ stars in the
South in our sample. Figure 12 displays scale-height and scale-length
distributions from SMSS $\cap$ Gaia in the Southern Galactic Hemisphere,
generated following the same steps used for Figures 10 and 11. The coherent
structure in the parameter space covered by GSE is not seen as clearly as in
the Northern Galactic Hemisphere from SDSS $\cap$ Gaia, but resembles the
distributions in the South, highlighting the highly eccentric orbit of the GSE
merger. On the other hand, owing to the limited overlap along the Galactic
prime meridian, the SDSS $\cap$ Gaia sample is biased toward the Galactic
anticenter direction in the Southern Hemisphere, while the sample from SMSS
$\cap$ Gaia is more populated toward the Galactic center. Therefore, a trace
of the intermediate scale-height valley in the top panel could be real, and
may capture the GSE debris in the direction toward the Galactic center. Deeper
SMSS photometry in future data releases would be useful to explore this volume
in more detail.

In the top panel of Figure 10, a vertical trough with intermediate scale
heights lies at [Fe/H] $\sim-0.8$, extending from $v_{\phi}\sim 0\ \rm
km~{}s^{-1}$ to $v_{\phi}\sim 100\ \rm km~{}s^{-1}$. Originally, Belokurov et
al. (2020) defined a region occupied by the Splash as $-0.7<{\rm[Fe/H]}<-0.2$
and $-150<v_{\phi}<+100\ \rm km~{}s^{-1}$, which overlaps with this trough.
The fact that the intermediate scale-height valley from GSE stretches out to a
region populated by Splash stars supports a previous claim that the GSE merger
has dynamically heated stars in the primordial disk of the Milky Way (Bonaca
et al., 2017; Belokurov et al., 2020). In other words, because a nearly in-
plane collision with a dwarf galaxy would leave behind heated stars confined
to a disk plane, this apparent coincidence supports a view on a causal
connection between the Splash and GSE. In addition, since such stars
originated from the primordial disk of the Galaxy, which was smaller in size
than the current stellar disk, this also naturally explains the short scale
length of these stars (Figure 11).

If the small scale length is an indication that stars were born in the inner
region, and then were displaced by mergers, the vast stretch of the lowland
with small scale lengths in Figure 11 implies that dynamical heating took
place over a wider range of metallicity than previously considered by
Belokurov et al. (2020). At [Fe/H] $>-1$, stars with the lowest scale length
are located at $v_{\phi}\sim+100\ \rm km~{}s^{-1}$. At lower metallicities,
the nearly horizontal trough ($v_{\phi}\sim+150\ \rm km~{}s^{-1}$ at [Fe/H]
$\approx-1.5$) coincides with the MWTD (see Paper III). If the lowland in the
top panel of Figure 11 is fenced by a contour line at $L=1.1$ kpc and is taken
as an approximate extent of the dynamical heating from these mean lines, it
covers $-2<{\rm[Fe/H]}<0$. Therefore, even if the GSE progenitor was a major
source of the dynamical heating of a primordial disk, it seems unlikely that
the orbital properties of stars in this region are altered by a single massive
merger event. Instead, our result indicates that such a dynamical heating
process on the primordial disk was in operation, even before the GSE merger,
driven by more numerous minor mergers, as predicted by numerical simulations
of the early universe (e.g., Grand et al., 2020).

The other interesting feature is the “highland” in the scale-height
distribution (‘\raisebox{-.9pt} {3}⃝’ in Figure 10), which coincides with the
metal-poor in situ halo (see Figure 9). Its scale height reaches $H\sim 1.1$
kpc, which is significantly higher than those of heated stars ($\sim 0.8$
kpc); therefore, it is unlikely that it formed through the same dynamical
heating process. Moreover, the highland exhibits a gradient in the scale
length, as shown in the top panel of Figure 11, in the sense that more metal-
rich stars are more strongly concentrated in the inner region of the Galaxy.
This implies that more active star formation took place in the deeper
potential well, while the Milky Way has grown by chaotic coalescence of
numerous small gas-rich dwarf galaxies in the early universe. It is unclear,
however, how the Milky Way attained the net angular momentum in the same
direction as that of the Galactic disk at this stage. Nonetheless, its small
prograde net rotation ($\langle v_{\phi}\rangle\sim+60\ \rm km~{}s^{-1}$)
suggests that it had maintained a puffy disk-like structure owing to turbulent
nature of the gas-rich minor mergers.

###  3.3 The Galactic Starburst Sequence (GSS)

Figure 13: A phase-space diagram of stars with large proper motions (PM) from
SDSS $\cap$ Gaia.

Our phase-space diagrams also reveal a long and narrow sequence of stars when
a simple cut on proper motion is made. All four panels in Figure 13 are drawn
from the same SDSS $\cap$ Gaia sample, but with different cuts on a minimum
value of proper motions ($20$, $40$, $60$, and $70$ mas yr-1, respectively).
The top left panel is dominated by disk populations (mostly thick-disk stars),
but a low $v_{\phi}$, low [Fe/H] tail begins to show up when a mild cut on
proper motion ($20$ mas yr-1) is imposed. At $>40$ mas yr-1, a striking elbow-
like feature emerges from these diagrams, where halo and disk stars form a
narrow, continuous sequence over a wide range of [Fe/H] and $v_{\phi}$. This
feature is characterized by two joint, orthogonal branches; the horizontal arm
is nearly parallel to the $v_{\phi}=0\ \rm km~{}s^{-1}$ line over a wide range
of metallicity ($-2\lesssim{\rm[Fe/H]}\lesssim-0.4$), while the vertical arm
has a narrow metallicity range ($-0.8\lesssim{\rm[Fe/H]}\lesssim-0.2$). The
sequence passes through GSE and the Splash, and is eventually connected to the
disk. Although our $v_{\phi}$ estimates depend on distance, the elbow-like
feature becomes stronger with higher proper-motion cuts, but changes little
with distance from the Sun.

Figure 14: Same as in Figure 13, but from SMSS $\cap$ Gaia. Figure 15: (a)
Same as in Figure 13, but a phase-space diagram of stars for which Gaia’s
metallicity estimates are available. A vertical distance range ($1<|Z|<2$ kpc)
and a proper-motion cut ($>25$ mas yr-1) are adjusted to emphasize the
observed sequence more clearly. (b) Same as in panel (a), but based on three
dimensional motions in Gaia.

As shown in Figure 14, the sequence persists even when SMSS $\cap$ Gaia is
used, indicating that it is present in both hemispheres. Aside from our
photometric metallicities, the same sequence can also be seen in Figure 15,
which is based on spectroscopic metallicities in Gaia DR3 from the Radial
Velocity Spectrometer (RVS) spectra ($\lambda/\Delta\lambda=11,500$) using the
General Stellar Parameteriser-Spectroscopy (GSP-Spec) module (Recio-Blanco et
al., 2022). Because these spectroscopic observations are available for bright
stars ($G<14$ mag), the sample is limited to relatively nearby stars
($1<|Z|<2$ kpc); the lower $|Z|$ cut is made to exclude numerous disk stars.
The proper-motion cut is lowered to $25$ mas yr-1 to retain as many stars as
possible along the sequence. On the other hand, the constraint on Galactic
latitudes is lifted, as the spectroscopic metallicities are only weakly
dependent on foreground reddening. While the left panel in Figure 15 is based
on the projected $v_{\phi}$ as for our photometric samples along the prime
meridian, the right panel is based on full three-dimensional space motions
from both radial velocities and proper motions in Gaia. The sequence appears
almost identical in both panels, which supports the existence of this coherent
structure. Besides, this also validates our approach for obtaining $v_{\phi}$
from proper motions only.

The continuous chain of stars in the phase-space diagram implies that they
formed through successive enrichment of metals along this pathway, and somehow
it is connected to the GSE merger. Strikingly, the vertical arm of the
sequence has a very narrow metallicity range, indicating that these stars
formed in a short gas-depletion time scale without sufficient time for
successive enrichment of metals (but with $\alpha$-element enhancement from
core-collapse supernovae), while star-forming clouds were falling or changing
the orientation of an angular momentum vector rapidly from nearly zero to
$\sim 200\ \rm km~{}s^{-1}$ in $v_{\phi}$. In other words, our result
indicates that the young Milky Way went through a phase of starburst activity.
For this reason, we refer to this structure as the Galactic Starburst Sequence
(GSS). Notably, the full GSS is not seen in panel (d) of Figure 9, although it
is also made using high proper-motion stars. This difference is manifested by
a high central concentration of Splash stars in the Galaxy, as shown in Paper
III, and more clearly demonstrated by the short scale lengths in Figure 11.
Thus, the GSS is not a dominant structure at large Galactocentric distances.

Figure 16: Contours of equal scale height (red solid) and length (blue
dotted), plotted on top of a phase-space diagram of stars with large proper
motions ($>60$ mas yr-1; bottom left panel in Figure 13). The scale-height
contour delineates $H=1.0$ kpc (top panel of Figure 10), and the scale-length
contour has $L=1.1$ kpc (top panel of Figure 11).

To better understand the properties of the GSS, we combine the above results
in Figure 16. The red solid contour shows a region with a scale height $H=1.0$
kpc from the top panel of Figure 10, while the blue dashed contour delineates
a scale length $L=1.1$ kpc from the top panel of Figure 11. They are overlaid
on top of contours of stars with high proper motions ($>60$ mas yr-1) in
Figure 13. First, the width of the vertical arm of the GSS is narrower than
the metallicity range of the heated population, the extent of which can be
delineated by the blue dotted enclosure with a small scale length. This
indicates that the Splash, which can be defined as a group of metal-rich stars
with halo-like kinematics, is in fact composed of two distinct groups of
stars. The first is a group of dynamically heated stars, and the other is the
starburst population. Heated stars formed in the inner region of the
primordial Galaxy and then were displaced to the current location in the halo,
while starburst stars formed in a top-down fashion, as indicated by a weak
positive correlation between [Fe/H] and $v_{\phi}$.

Secondly, both GSE and the GSS lie along the intermediate scale-height valley,
as shown by the red solid line in Figure 16. The simplest explanation is that
the GSE merger has triggered a starburst in the mixture of gas from the
primordial Milky Way and gas donated by the GSE merger, because mergers were
likely gas rich at high redshifts. According to this scenario, the radially
biased orbit of GSE should be responsible for the relatively high transverse
motions of stars in the GSS, which contrast with a nearly isotropic velocity
distribution of metal-poor in situ halo stars. In other words, these stars
formed out of metal-enriched gas having similar orbital properties with the
GSE progenitor. Therefore, they are expected to show up in the high proper-
motion sample, even though their orbits have evolved significantly over time
from halo-like to disk-like orbits.

##  4 Summary and Discussions

Thanks to Gaia and large-scale photometric and spectroscopic databases, it is
now possible to perform accurate comparisons of theoretical models to
extensive, high-quality data. In this work, we use models generated using YREC
and MARCS, and quantify model deviations as a function of $T_{\rm eff}$,
metallicity, and wavelength. The current approach relies on a comparison with
multi-wavelength ultra-violet, optical, and near-infrared color-$T_{\rm eff}$
relations derived from Galactic cluster sequences, Gaia’s double sequence, and
a sample of spectroscopic data in SEGUE and GALAH, using photometric data in a
number of broad-band filters from SDSS, SMSS, and APASS and standard-star
photometry in the literature. Mean flux deviations are derived as a function
of wavelength, which amount to up to $\sim 20\%$ at $<4000$ Å, but
significantly less in longer wavelengths. We find that no single factor can
remove the observed offsets, but it is more likely a problem arising from a
combination of various sources of errors in the models and/or observational
data. Subsequently, we define the model offset as an empirical correction
function for our specific choice of models.

When combined with proper-motion measurements from Gaia, our technique can be
used to construct phase-space ($v_{\phi}$ versus [Fe/H]) diagrams of stars to
provide a global perspective on the stellar populations in the Milky Way. In
this way, we identify a long and narrow sequence of stars in a phase-space
diagram, which we call the GSS. The GSS overlaps with GSE, and runs along a
valley with intermediate scale heights, suggesting that GSE has likely
triggered star formation. It also passes through the Splash, showing rapid
evolution of $v_{\phi}$ of star-forming clouds within a narrow metallicity
range, which testifies to a starburst event in the young Milky Way. The wide
metallicity range of dynamically heated stars as traced by small scale-length
regions in a phase-space diagram indicates that the Splash is likely composed
of two stellar populations with distinct origins – dynamically heated and
starburst populations.

The double sequence of stars with high transverse motions ($>200\ \rm
km~{}s^{-1}$) in Gaia (Gaia Collaboration et al., 2018a) is another
manifestation of the GSS. As shown in Paper II, the two sequences are
separated by a few tenths of a magnitude in $griz$ colors, and more strongly
when the $u$-band is included, indicating a metallicity offset by $\sim 1$
dex. According to the corresponding metallicity ranges, stars on the blue
sequence belong to GSE, while stars in the red sequence mostly constitute the
Splash or, more specifically, the vertical arm of the GSS. Our new perspective
into the phase-space diagram of these stars reveals two consecutive modes of
chemo-kinematical evolution: a rapid chemical enrichment along the orbit of
GSE, followed by starburst in rapidly evolving orbits of gas clouds.

Previously, Necib et al. (2020) identified a group of $\sim 100$ stars in the
solar vicinity, dubbed “Nyx”, which show coherent radial ($134\ \rm
km~{}s^{-1}$) and azimuthal ($130\ \rm km~{}s^{-1}$) motions. Their
metallicities peak at [Fe/H] $=-0.55$, with a dispersion of $0.13$ dex. In
addition, they are mostly confined to a plane with a maximum vertical distance
of $Z=1.7$ kpc, which is comparable to our scale-height measurement ($H\approx
0.7$ kpc) for stars having similar [Fe/H] and $v_{\phi}$ as Nyx stars (Figure
10). The mean metallicity and a small scale height suggest its close
connection to stars in the vertical arm of the GSS. Nonetheless, most of Nyx
stars have $100\lesssim v_{\phi}<250\ \rm km~{}s^{-1}$, with a few extreme
cases, and therefore trace only the upper half of the vertical arm. Necib et
al. (2020) postulated that the Nyx stream is a remnant of a disrupted dwarf
galaxy. However, a more recent analysis in Zucker et al. (2021) argues against
their extragalactic origin, based on the fact that a subset of these stars are
indistinguishable from thick-disk stars in the elemental-abundance space,
having larger $\alpha$-element abundances than those of accreted stars at a
given metallicity.

In other studies, age-metallicity relations also reveal a chain of metal-rich
stars, establishing a link between the halo and the disk (Haywood et al.,
2013; Nissen et al., 2020; Xiang & Rix, 2022). Recently, Xiang & Rix (2022)
used $\alpha$-element-rich stars with low orbital angular momenta, and
demonstrated the existence of a narrow and continuous age-metallicity relation
in the Milky Way. According to their analysis, the Milky Way achieved a high
metallicity floor ([Fe/H] $\sim-1$) about $13$ Gyr ago, with successive metal
enrichment over the following $\sim 5$ Gyr. The majority of stars are found in
a narrow metallicity range ($-0.7<{\rm[Fe/H]}<-0.2$), which leads to its
possible connection to the vertical arm of the GSS. Importantly, its narrow
and continuous channel of stars reinforces the physical entity of the GSS.

In addition, recent numerical simulations in the literature also predict the
existence of starburst populations in the Milky Way-size galaxies (Cooper et
al., 2015; Bignone et al., 2019; Grand et al., 2020; Renaud et al., 2021). For
instance, in Grand et al. (2020), a merger with a gas-rich GSE-like progenitor
induces starburst activity owing to the increased compression of gas clouds in
both galaxies. The star-formation rate suddenly increases by a factor of two,
which lasts less than $1$ Gyr. Interestingly, these simulations reproduce some
of the key features observed in our data, in that both starburst and heated
populations produced by the merger event are present in a narrow metallicity
bin, but cover a wide range of $v_{\phi}$. Furthermore, the heated population
shows a correlation between [Fe/H] and $v_{\phi}$ in their simulation, while
there is essentially no such dependence for the starburst population
(considering an excessively large scatter in $v_{\phi}$ at a given
metallicity). This may well be explained by a narrow range in [Fe/H] of the
vertical arm of the GSS. Another intriguing aspect is that the starburst
population in the simulation exhibits a centrally concentrated, rotationally
supported disk. This prediction is also in agreement with the fact that the
high proper-motion stars are found along the intermediate scale-height valley.

The vertical and horizontal arms of the GSS can be attributed to stars formed
during a merger-induced star formation and those accreted from the GSE
progenitor, respectively. In our phase-space diagram in Figures 13 and 14, we
count each group of stars using a simple box criterion: $-50<v_{\phi}\leq 200\
\rm km~{}s^{-1}$ and $-0.9<{\rm[Fe/H]}\leq 0$ for the starburst, and
$-60<v_{\phi}\leq+60\ \rm km~{}s^{-1}$ and $-2<{\rm[Fe/H]}\leq-0.9$ for the
accreted population, respectively. We find that the ratios between these two
populations are $1.2$–$4.4$, depending on the proper motion cut made in each
subset. These estimates qualitatively agree with Grand et al. (2020), in that
the total mass of stars formed during a merger-induced star formation in the
host galaxy exceeds the stellar mass of accreted stars. Strikingly, they are
quantitatively in agreement with a ratio of $\sim 2$–$5$ predicted by a more
recent numerical simulation in Orkney et al. (2022).

At high $v_{\phi}$, the GSS is connected to the disk, suggesting that disk
stars formed in part from gas clouds left behind after the starburst episode.
Interestingly, the mean metallicity of the tip of the GSS ([Fe/H]
$\approx-0.4$) is coincident with the low-metallicity end of the thin disk, if
the disk stars are divided into the low and high $\alpha$-element sequences on
the [$\alpha$/Fe]-[Fe/H] plane (e.g., Hayden et al., 2015, and references
therein). This implies that old stars in the low-$\alpha$ sequence were formed
out of a mixture of gas clouds in the protogalactic disk and those accreted by
the GSE merger, possibly through the classical “two-infall” model (Chiappini
et al., 1997; Brook et al., 2007; Buck, 2020). The lower metallicity of old
thin-disk stars than those of the most metal-rich thick-disk stars and their
intermediate [$\alpha$/Fe] can be understood by the dilution of the metal-rich
clouds in the protogalactic disk and a boost in [$\alpha$/Fe] by core-collapse
supernovae during the starburst.

Our finding naturally explains the G-dwarf problem as well — the apparent
excess of metal-rich stars in the local disk compared to the prediction from a
simple closed-box chemical model (e.g., Greener et al., 2021). The above
interpretation is further supported by recent findings on in situ stars that
lie between the GSE and the low-$\alpha$ sequence (Myeong et al., 2022) and
the fraction of Splash stars that exhibit intermediate [$\alpha$/Fe] (Lee, A.
et al. 2022, submitted).

The extensibility of calibrated synthetic spectra presented in this work
enables accurate predictions of stellar magnitudes for filter passbands in
various photometric surveys, such as the Javalambre/Southern Photometric Local
Universe Survey (J/S-PLUS; Cenarro et al., 2019; Mendes de Oliveira et al.,
2019) and the Legacy Survey of Space and Time (LSST; Ivezić et al., 2019). Our
metallicity-mapping technique can serve as a useful resource for studying the
demographics of stellar populations that are yet to be discovered in the local
universe.

We thank Donald M. Terndrup and Christopher A. Onken for carefully checking
our manuscript and providing a number of helpful comments. D.A. thanks Anirudh
Chiti for a useful discussion. T.C.B. and Y.S.L. thank Andreia Carrillo for an
enlightening discussion on the Galactic starburst. D.A. acknowledges support
provided by the National Research Foundation (NRF) of Korea grant funded by
the Ministry of Science and ICT (No. 2021R1A2C1004117). T.C.B. and Y.S.L.
acknowledge partial support for this work from grant PHY 14-30152; Physics
Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), and
OISE-1927130: The International Research Network for Nuclear Astrophysics
(IReNA), awarded by the US National Science Foundation. Y.S.L. acknowledges
support from the NRF of Korea grant funded by the Ministry of Science and ICT
(NRF-2021R1A2C1008679). T.M. acknowledges financial support from the Spanish
Ministry of Science and Innovation (MICINN) through the Spanish State Research
Agency, under the Severo Ochoa Program 2020-2023 (CEX2019-000920-S) as well as
support from the ACIISI, Consejería de Economía, Conocimiento y Empleo del
Gobierno de Canarias and the European Regional Development Fund (ERDF) under
grant with reference PROID2021010128. Funding for the Sloan Digital Sky Survey
IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of
Energy Office of Science, and the Participating Institutions. SDSS-IV
acknowledges support and resources from the Center for High Performance
Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV
is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration including the Brazilian Participation
Group, the Carnegie Institution for Science, Carnegie Mellon University,
Center for Astrophysics — Harvard & Smithsonian, the Chilean Participation
Group, the French Participation Group, Instituto de Astrofísica de Canarias,
The Johns Hopkins University, Kavli Institute for the Physics and Mathematics
of the Universe (IPMU) / University of Tokyo, the Korean Participation Group,
Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik
Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-
Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für
Extraterrestrische Physik (MPE), National Astronomical Observatories of China,
New Mexico State University, New York University, University of Notre Dame,
Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State
University, Shanghai Astronomical Observatory, United Kingdom Participation
Group, Universidad Nacional Autónoma de México, University of Arizona,
University of Colorado Boulder, University of Oxford, University of
Portsmouth, University of Utah, University of Virginia, University of
Washington, University of Wisconsin, Vanderbilt University, and Yale
University. This work presents results from the European Space Agency (ESA)
space mission Gaia. Gaia data are being processed by the Gaia Data Processing
and Analysis Consortium (DPAC). Funding for the DPAC is provided by national
institutions, in particular the institutions participating in the Gaia
MultiLateral Agreement (MLA). The Gaia mission website is
https://www.cosmos.esa.int/gaia. The Gaia archive website is
https://archives.esac.esa.int/gaia. This paper makes use of data from the
AAVSO Photometric All Sky Survey, whose funding has been provided by the
Robert Martin Ayers Sciences Fund and from the NSF (AST-1412587).

##  Appendix A Construction of the Gaia Double Sequence and Model Comparisons

Figure 17: Gaia’s double sequence in selected passbands. The colored histogram
shows a number density of stars with $-150<v_{\phi}\leq+150\ \rm km~{}s^{-1}$.
The blue and red lines are our derived fiducial sequences. Figure 18:
Comparisons of Stetson’s photometry with theoretical models for selected
clusters. Figure 19: Comparison between the observed and original model colors
for the SEGUE sample in the SDSS passbands. The above example shows
differences of extinction-corrected colors using stars at
$-0.5<{\rm[Fe/H]}<-0.3$ as a function of photometric $T_{\rm eff}$. The solid
line shows average differences in moving boxes, and the error bars indicate
uncertainties of the mean values. Figure 20: Same as in Figure 19, but in the
SMSS filters. To avoid confusion, SDSS $r$ is shown as $r^{\prime}$. Figure
21: Same as in Figure 19, but in the PS1 filters. To avoid confusion, SDSS $r$
is shown as $r^{\prime}$. Figure 22: Same as in Figure 19, but in the APASS
filters. To avoid confusion, SDSS $r$ is shown as $r^{\prime}$.

Here, a number of complementary plots are provided for some selected data
sets, which are used in our model calibration (§ 2). Figure 17 displays Gaia’s
double sequences in SDSS, SMSS, PS1, and APASS using color indices containing
$u$ or $B$ passbands, which best separate the blue and red sequences of stars
having $-150<v_{\phi}\leq 150\ \rm km~{}s^{-1}$. Figure 18 displays
comparisons of the original models (without empirical corrections) with
Stetson’s cluster photometry. Model comparisons with the SEGUE sample
($-0.5<{\rm[Fe/H]}<-0.3$) are shown in Figures 19–22 in SDSS, SMSS, PS1, and
APASS photometry, respectively. See Paper II for more information on the
construction of Gaia’s double sequence and model comparisons of the cluster
and the Gaia double sequence in the SDSS filter set.

##  Appendix B Photometric Properties of SDSS, SMSS, and PS1

###  B.1 Random and Zero-Point Uncertainties

Figure 23: Statistical properties of SMSS (left) and PS1 photometry (right)
from comparisons with SDSS. Top: A best-fitting Gaussian distribution of
magnitude offsets in each passband, after subtracting a color-dependent
difference between different passbands. Standard deviations are shown in the
inset. Bottom: A distribution of the ratio between a standard deviation of a
magnitude difference [$\sigma(\Delta{\rm mag})$] and a propagated uncertainty
[$\sigma({\rm phot})$] in each filter passband. The mean ratios are shown in
the inset. Both SMSS $u$\- and $v$ bands are compared to SDSS $u$, and PS1 $y$
to SDSS $z$.

Since observations in various filter passbands are combined to derive a set of
parameters in this study, it is necessary to adopt realistically sized
photometric uncertainties in the $\chi^{2}$ statistics. Here, we compare SMSS
and PS1 photometry with SDSS, and compute a mean magnitude difference and a
dispersion to infer the size of the true uncertainties. The comparison with
SMSS is limited to narrow regions, since the imaging stripes of SDSS overlaps
only a little with SMSS footprints along the celestial equator. The overlap
with PS1 is more extensive, as both surveys covered the Northern Hemisphere.

Figure 23 shows statistical properties of the comparison of SDSS photometry
with SMSS (left) and PS1 (right), respectively. The top panels display a
distribution of a magnitude difference of bright ($14.5<r<18$ mag) stars in
each passband. To take into account non-negligible color terms between
different filter passbands (i.e., color transformations), the magnitude
differences in each $1\arcdeg$-wide strip in R.A. are fit using a third-order
polynomial as a function of $g\,-\,r$ in $0.3<g\,-\,r<1.1$, and its mean trend
is removed. The SMSS $v$ is compared to SDSS $u$, and PS1 $y$ to SDSS $z$.

A weighted median difference is computed in each
$(\Delta\alpha,\Delta\delta)=(1\arcdeg,1\arcdeg$) region, and an ensemble of
these differences are fit using a Gaussian function, as shown by solid lines
in the top panels of Figure 23. Standard deviations of the best-fitting
Gaussian functions provide a measure of the spatial variation of the
photometric zero points across the sky. They are $\sigma_{\rm zp}^{\rm
SMSS}(u,v,g,r,i,z)=\\{13,14,4,4,4,6\\}$ mmag and $\sigma_{\rm zp}^{\rm
PS1}(g,r,i,z,y)=\\{5,5,5,6,8\\}$ mmag, respectively, which are comparable to
the quoted systematic uncertainties by each survey.

The bottom panels of Figure 23 show ratios between the standard deviation of a
magnitude difference and a propagated uncertainty for each filter passband in
each $(\Delta\alpha,\Delta\delta)=(1\arcdeg,1\arcdeg$) patch. Again, SMSS $v$
is compared to SDSS $u$, and PS1 $y$ to SDSS $z$. Some of the brightest
objects are rejected in this comparison, owing to unrealistically small
uncertainties in the SDSS PSF $i$ magnitude ($<0.005$ mag). As shown in the
bottom left panel, the ratios for the SMSS and SDSS $grz$ are near unity,
indicating that photometric uncertainties are comparable to the observed
scatter, while photometric uncertainties in $uv$ are likely under-estimated.
As shown in the bottom right panel, the differences between propagated
uncertainties and observed dispersions are even larger for PS1 passbands.

Based on the above comparisons, an uncertainty ‘floor’ ($\sigma_{f}$) is
computed in each passband in order to make a median of a standard deviation
equal to a median of propagated uncertainties from all $1\arcdeg\times
1\arcdeg$ patches. Assuming that quoted uncertainties in SDSS are correct, we
find $\sigma_{f}^{\rm SMSS}(u,v,g,r,i,z)=\\{42,42,10,9,9,10\\}$ mmag and
$\sigma_{f}^{\rm PS1}(g,r,i,z,y)=\\{43,44,42,44,43\\}$ mmag for SMSS and PS1,
respectively. In all cases, zero-point uncertainties ($\sigma_{\rm zp}$) are
overwhelmed by the uncertainty floors ($\sigma_{f}$). We add both
uncertainties in quadrature to the original photometric uncertainties in SMSS
and PS1, and use them throughout this work.

###  B.2 Re-calibration of SMSS $uv$-Band Photometry

By design, the zero point of SMSS DR2 photometry was set based on $griz$
photometry in the PS1 system calculated from all-sky Gaia photometry (Onken et
al., 2019). The wavelengths of SMSS $griz$ passbands overlap with Gaia and PS1
passbands, so in principle one can tie them together without losing
information on the properties of stellar spectra. On the contrary, photometric
zero points in $u$-band and $v$-band remain largely unconstrained due to the
lack of short-wavelength passbands in Gaia or PS1.

Figure 24: Panel (a): A projection map of mean metallicities in the Galactic
coordinate system based on the original photometry. The color scheme shows a
mean metallicity of stars at $1.2<d<3$ kpc from a generalized histogram in
each pixel. Panel (b): Zero-point corrections on $u$-band and $v$-band
photometry, in order to match the peak metallicities of nearby stars
($0.5<d<1$ kpc) to a reference value (see text). Panel (c): Same as in panel
(a), but based on the re-calibrated $uv$ photometry. Panel (d): Zero-point
corrections in the $u$-band in Huang et al. (2022). In all panels, the
Galactic center is at the center, and the North Galactic Pole is to the top.
Areas at low Galactic latitudes ($|b|<20\arcdeg$) and large cumulative
extinction [$E(B\,-\,V)>0.1$] are excluded. Note that the above metallicity
maps are not exactly the same as in Figure 8, owing to different ways of
estimating metallicities with and without Gaia parallaxes, respectively.

Panel (a) of Figure 24 shows the mean metallicity distribution of stars
($1.2<d<3$ kpc) from the SMSS in the Galactic coordinate system, which
demonstrates the necessity for a second-order photometric zero-point
correction. Here, we use metallicities from a fully photometric solution, as
it is more prone to photometric errors than the case based on the Gaia
parallax. The mean metallicity in each pixel is derived from a generalized
histogram of photometric metallicities, which accounts for an uncertainty in
metallicity by taking it as a standard deviation of a normal probability
distribution. We implement HEALPix (Górski et al., 2005) in the Galactic
coordinate system, for which we set a resolution parameter $N_{\rm side}$ to
$32$, corresponding to an equal area of $3.3\ {\rm deg}^{2}$ for each pixel.
Although the bright survey limit restricts the sample to relatively nearby
stars, nearly uniform metallicities of stars are contrary to what is expected
in the local volume.

To obtain zero-point corrections on SMSS $uv$ photometry, we assume that
nearby stars ($0.5<d<1$ kpc) have the same metallicity in every direction, as
they are mostly thin-disk stars. Any deviation in the mean metallicity is
attributed entirely to a zero-point error in the $u$\- and $v$-band, because
of their larger zero-point uncertainties (Figure 23) and stronger
sensitivities on metallicity than other passbands. We also assume the same
amount of offset in both passbands. The sensitivity of our metallicity
estimate on zero point is estimated using high-latitude stars
($|b|>60\arcdeg$) by comparing with a case assuming an arbitrary $0.06$ mag
offset in the $u$\- and $v$-band, from which we find
$\Delta{\rm[Fe/H]}/\Delta{\rm mag}=-2.3$. Using this factor, the metallicity
map is forced to match a reference metallicity ([Fe/H]$=-0.28$), which is
taken from the average metallicity in the above volume. As shown in panel (b)
of Figure 24, we adopt $N_{\rm side}=16$ (a pixel area of $13.4\ {\rm
deg}^{2}$) for the zero-point correction map, because a higher spatial
resolution results in the loss of some pixels with small numbers of stars,
while information on spatial dependence is lost on a lower resolution map. In
a new metallicity map based on zero-point corrections [panel (c)], a global
change of the mean metallicity from low- to high-latitude regions is evident,
as expected from a simple population gradient from the disk to the halo.

For comparison, panel (d) in Figure 24 shows the zero-point corrections in
Huang et al. (2021, 2022), who used spectroscopic estimates from the GALAH
survey and Gaia parallaxes to derive photometric offsets in each of the SMSS
filters. It shows $u$-band corrections, but similar patterns and amplitudes
are seen in the $v$-band. The median difference from our map is negligible
($0.004$ mag and $0.001$ mag in $u$\- and $v$-band, respectively, in the sense
of our study minus their values), and a standard deviation amounts to $0.025$
mag in both bands. However, high-latitude regions ($|b|>-50\arcdeg$) are only
sparsely populated by GALAH targets, and therefore their correction functions
are weakly constrained in the polar regions. In this respect, our correction
map provides more complete information for our chemo-kinematic sample along
the Galactic prime meridian, which passes through the high-latitude regions.

###  B.3 Re-calibration of SDSS $u$-Band Photometry

Figure 25: Same as in panels (a)–(c) in Figure 24, but from SDSS photometry.
Note that the above metallicity maps [panels (a) and (c)] are not exactly same
as in Figure 8, owing to different ways of estimating metallicities with and
without Gaia parallaxes, respectively.

In the same manner as for the SMSS $u$ and $v$ photometry, SDSS $u$-band
photometry is re-calibrated for a small, but significant offset in our
metallicity map. Panel (a) in Figure 25 shows a mean metallicity distribution
from stars at $1.2<d<3$ kpc. There are strips with distinctly lower or higher
metallicities than surrounding areas (e.g., the stripe along $l=30\arcdeg$;
see also Figure 1 in Paper I). These strips are generally parallel to the
scanning direction of the SDSS imaging footprints (each $2.5\arcdeg$ wide and
$\sim 120\arcdeg$ long), suggesting that the spatially correlated offsets in
metallicity are induced by photometric zero-point errors of the metallicity-
sensitive $u$-band.

The zero-point correction map in panel (b) is derived in the same way as for
the SMSS $u$ and $v$ passbands. As there are more stars available for the
construction of a metallicity distribution function in SDSS, we adopt a finer
pixel size for HEALPix, $N_{\rm side}=32$. The metallicity sensitivity
$\Delta{\rm[Fe/H]}/\Delta{\rm mag}=-2.6$ is obtained from a case with a $0.06$
mag offset in $u$, and is used to make a uniform metallicity distribution of
nearby stars ($0.5<d<1$ kpc) at the ensemble average ([Fe/H]$=-0.25$). In this
way, only the relative zero-point offsets are rectified, while the global mean
remains intact. The revised metallicity map shown in panel (c) is
significantly more smooth than in panel (a), and shows no artificial
metallicity structures.

## References

  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42. doi:10.3847/1538-4365/aa9e8a
  * An & Beers (2020) An, D. & Beers, T. C. 2020, ApJ, 897, 39. doi:10.3847/1538-4357/ab8d39 (Paper I)
  * An & Beers (2021a) An, D. & Beers, T. C. 2021a, ApJ, 907, 101. doi:10.3847/1538-4357/abccd2 (Paper II)
  * An & Beers (2021b) An, D. & Beers, T. C. 2021b, ApJ, 918, 74. doi:10.3847/1538-4357/ac07a4 (Paper III)
  * An et al. (2013) An, D., Beers, T. C., Johnson, J. A., et al. 2013, ApJ, 763, 65. doi:10.1088/0004-637X/763/1/65
  * An et al. (2015) An, D., Beers, T. C., Santucci, R. M., et al. 2015, ApJ, 813, L28. doi:10.1088/2041-8205/813/2/L28
  * An et al. (2008) An, D., Johnson, J. A., Clem, J. L., et al. 2008, ApJS, 179, 326. doi:10.1086/592090
  * An et al. (2009) An, D., Pinsonneault, M. H., Masseron, T., et al. 2009, ApJ, 700, 523. doi:10.1088/0004-637X/700/1/523
  * An et al. (2019) An, D., Pinsonneault, M. H., Terndrup, D. M., et al. 2019, ApJ, 879, 81. doi:10.3847/1538-4357/ab23ed
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., et al. 2021, AJ, 161, 147. doi:10.3847/1538-3881/abd806
  * Beers et al. (2012) Beers, T. C., Carollo, D., Ivezić, Ž., et al. 2012, ApJ, 746, 34. doi:10.1088/0004-637X/746/1/34
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., et al. 2018, MNRAS, 478, 611. doi:10.1093/mnras/sty982
  * Belokurov et al. (2020) Belokurov, V., Sanders, J. L., Fattahi, A., et al. 2020, MNRAS, 494, 3880. doi:10.1093/mnras/staa876
  * Belokurov & Kravtsov (2022) Belokurov, V. & Kravtsov, A. 2022, MNRAS, 514, 689. doi:10.1093/mnras/stac1267
  * Bernard et al. (2014) Bernard, E. J., Ferguson, A. M. N., Schlafly, E. F., et al. 2014, MNRAS, 442, 2999. doi:10.1093/mnras/stu1081
  * Bessell et al. (2011) Bessell, M., Bloxham, G., Schmidt, B., et al. 2011, PASP, 123, 789. doi:10.1086/660849
  * Bessell & Murphy (2012) Bessell, M. & Murphy, S. 2012, PASP, 124, 140. doi:10.1086/664083
  * Bignone et al. (2019) Bignone, L. A., Helmi, A., & Tissera, P. B. 2019, ApJ, 883, L5. doi:10.3847/2041-8213/ab3e0e
  * Bohlin et al. (2014) Bohlin, R. C., Gordon, K. D., & Tremblay, P.-E. 2014, PASP, 126, 711. doi:10.1086/677655
  * Bonaca et al. (2020) Bonaca, A., Conroy, C., Cargile, P. A., et al. 2020, ApJ, 897, L18. doi:10.3847/2041-8213/ab9caa
  * Bonaca et al. (2017) Bonaca, A., Conroy, C., Wetzel, A., et al. 2017, ApJ, 845, 101. doi:10.3847/1538-4357/aa7d0c
  * Brook et al. (2007) Brook, C., Richard, S., Kawata, D., et al. 2007, ApJ, 658, 60. doi:10.1086/511056
  * Buck (2020) Buck, T. 2020, MNRAS, 491, 5435. doi:10.1093/mnras/stz3289
  * Buder et al. (2021) Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150. doi:10.1093/mnras/stab1242
  * Carollo et al. (2010) Carollo, D., Beers, T. C., Chiba, M., et al. 2010, ApJ, 712, 692. doi:10.1088/0004-637X/712/1/692
  * Carollo et al. (2007) Carollo, D., Beers, T. C., Lee, Y. S., et al. 2007, Nature, 450, 1020. doi:10.1038/nature06460
  * Carretta et al. (2000) Carretta, E., Gratton, R. G., Clementini, G., et al. 2000, ApJ, 533, 215. doi:10.1086/308629
  * Casagrande et al. (2019) Casagrande, L., Wolf, C., Mackey, A. D., et al. 2019, MNRAS, 482, 2770. doi:10.1093/mnras/sty2878
  * Cenarro et al. (2019) Cenarro, A. J., Moles, M., Cristóbal-Hornillos, D., et al. 2019, A&A, 622, A176. doi:10.1051/0004-6361/201833036
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560
  * Chiappini et al. (1997) Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765. doi:10.1086/303726
  * Cooper et al. (2015) Cooper, A. P., Parry, O. H., Lowing, B., et al. 2015, MNRAS, 454, 3185. doi:10.1093/mnras/stv2057
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, A&A, 616, A10. doi:10.1051/0004-6361/201832843
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, A&A, 616, A1. doi:10.1051/0004-6361/201833051
  * 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. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv:2208.00211
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759. doi:10.1086/427976
  * Grand et al. (2020) Grand, R. J. J., Kawata, D., Belokurov, V., et al. 2020, MNRAS, 497, 1603. doi:10.1093/mnras/staa2057
  * Greener et al. (2021) Greener, M. J., Merrifield, M., Aragón-Salamanca, A., et al. 2021, MNRAS, 502, L95. doi:10.1093/mnrasl/slab012
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951. doi:10.1051/0004-6361:200809724
  * 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
  * Haywood et al. (2013) Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2013, A&A, 560, A109. doi:10.1051/0004-6361/201321397
  * Helmi (2020) Helmi, A. 2020, ARA&A, 58, 205. doi:10.1146/annurev-astro-032620-021917
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85. doi:10.1038/s41586-018-0625-x
  * Helmi et al. (1999) Helmi, A., White, S. D. M., de Zeeuw, P. T., et al. 1999, Nature, 402, 53. doi:10.1038/46980
  * Henden et al. (2018) Henden, A. A., Levine, S., Terrell, D., et al. 2018,
  * Huang et al. (2021) Huang, Y., Yuan, H., Li, C., et al. 2021, ApJ, 907, 68. doi:10.3847/1538-4357/abca37
  * Huang et al. (2022) Huang, Y., Yuan, H., Li, C., et al. 2022, ApJ, 924, 141. doi:10.3847/1538-4357/ac425d
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111. doi:10.3847/1538-4357/ab042c
  * Ivezić et al. (2008) Ivezić, Ž., Sesar, B., Jurić, M., et al. 2008, ApJ, 684, 287. doi:10.1086/589678
  * Kraft & Ivans (2003) Kraft, R. P. & Ivans, I. I. 2003, PASP, 115, 143. doi:10.1086/345914
  * Lee et al. (2008a) Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008a, AJ, 136, 2022. doi:10.1088/0004-6256/136/5/2022
  * Lee et al. (2008b) Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008b, AJ, 136, 2050. doi:10.1088/0004-6256/136/5/2050
  * Lejeune et al. (1997) Lejeune, T., Cuisinier, F., & Buser, R. 1997, A&AS, 125, 229. doi:10.1051/aas:1997373
  * Lejeune et al. (1998) Lejeune, T., Cuisinier, F., & Buser, R. 1998, A&AS, 130, 65. doi:10.1051/aas:1998405
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4. doi:10.1051/0004-6361/202039653
  * Mendes de Oliveira et al. (2019) Mendes de Oliveira, C., Ribeiro, T., Schoenell, W., et al. 2019, MNRAS, 489, 241. doi:10.1093/mnras/stz1985
  * Myeong et al. (2022) Myeong, G. C., Belokurov, V., Aguado, D. S., et al. 2022, arXiv:2206.07744
  * Naidu et al. (2020) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48. doi:10.3847/1538-4357/abaef4
  * Naidu et al. (2021) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2021, ApJ, 923, 92. doi:10.3847/1538-4357/ac2d2d
  * Necib et al. (2020) Necib, L., Ostdiek, B., Lisanti, M., et al. 2020, Nature Astronomy, 4, 1078. doi:10.1038/s41550-020-1131-2
  * Nissen et al. (2020) Nissen, P. E., Christensen-Dalsgaard, J., Mosumgaard, J. R., et al. 2020, A&A, 640, A81. doi:10.1051/0004-6361/202038300
  * Onken et al. (2019) Onken, C. A., Wolf, C., Bessell, M. S., et al. 2019, PASA, 36, e033. doi:10.1017/pasa.2019.27
  * Orkney et al. (2022) Orkney, M. D. A., Laporte, C. F. P., Grand, R. J. J., et al. 2022, arXiv:2206.09246
  * Pinsonneault et al. (2004) Pinsonneault, M. H., Terndrup, D. M., Hanson, R. B., et al. 2004, ApJ, 600, 946. doi:10.1086/379925
  * Recio-Blanco et al. (2022) Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2022, arXiv:2206.05541
  * Renaud et al. (2021) Renaud, F., Agertz, O., Read, J. I., et al. 2021, MNRAS, 503, 5846. doi:10.1093/mnras/stab250
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3. doi:10.1051/0004-6361/202039587
  * Rockosi et al. (2022) Rockosi, C. M., Lee, Y. S., Morrison, H. L., et al. 2022, ApJS, 259, 60. doi:10.3847/1538-4365/ac5323
  * Sahlholdt et al. (2019) Sahlholdt, C. L., Casagrande, L., & Feltzing, S. 2019, ApJ, 881, L10. doi:10.3847/2041-8213/ab321e
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103. doi:10.1088/0004-637X/737/2/103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525. doi:10.1086/305772
  * Sills et al. (2000) Sills, A., Pinsonneault, M. H., & Terndrup, D. M. 2000, ApJ, 534, 335. doi:10.1086/308739
  * Stetson (2000) Stetson, P. B. 2000, PASP, 112, 925. doi:10.1086/316595
  * Tonry et al. (2012) Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99. doi:10.1088/0004-637X/750/2/99
  * Tonry et al. (2018) Tonry, J. L., Denneau, L., Flewelling, H., et al. 2018, ApJ, 867, 105. doi:10.3847/1538-4357/aae386
  * Wolf et al. (2018) Wolf, C., Onken, C. A., Luvaul, L. C., et al. 2018, PASA, 35, e010. doi:10.1017/pasa.2018.5
  * Xiang & Rix (2022) Xiang, M. & Rix, H.-W. 2022, Nature, 603, 599. doi:10.1038/s41586-022-04496-5
  * Yan et al. (2019) Yan, Y., Du, C., Liu, S., et al. 2019, ApJ, 880, 36. doi:10.3847/1538-4357/ab287d
  * Yanny et al. (2009) Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377. doi:10.1088/0004-6256/137/5/4377
  * Zucker et al. (2021) Zucker, D. B., Simpson, J. D., Martell, S. L., et al. 2021, ApJ, 912, L30. doi:10.3847/2041-8213/abf7cd