# A Blueprint for the Milky Way’s Stellar Populations. V. 3D Local Dust
Extinction

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  Anirudh Chiti  Department of Astronomy & Astrophysics,
University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637, USA  Kavli
Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637,
USA

(Received November 23, 2023; Accepted March 15, 2024)

###### Abstract

Using a grid of empirically calibrated synthetic spectra developed in our
previous study, we construct an all-sky 3D extinction map from the large
collection of low-resolution XP spectra in Gaia DR3. Along each line of sight,
with an area ranging from $0.2$ to $13.4$ deg2, we determine both the
reddening and metallicity of main-sequence stars and model the foreground
extinction up to approximately $3$ kpc from the Sun. Furthermore, we explore
variations in the total-to-selective extinction ratio in our parameter search
and identify its mean systematic change across diverse cloud environments in
both hemispheres. In regions outside the densest parts of the clouds, our
reddening estimates are validated through comparisons with previous reddening
maps. However, a notable discrepancy arises in comparison to other independent
work based on XP spectra, which can be attributed to systematic offsets in
their metallicity estimates. On the other hand, our metallicity scale exhibits
reasonable agreement with the high-resolution spectroscopic abundance scale.
We also assess the accuracy of the XP spectra by applying our calibrated
models, and we confirm an increasing trend of flux overestimation at shorter
wavelengths below $400$ nm.

Unified Astronomy Thesaurus concepts: Interstellar dust extinction (837);
Milky Way Galaxy (1054); Stellar abundances (1577); Gaia (2360)

††journal: The Astrophysical Journal Supplement Series

##  1 Introduction

The availability of a vast amount of observational data for Galactic stars has
significantly advanced our understanding of the Milky Way’s structure and
evolution (e.g., Ivezić et al., 2012; Bland-Hawthorn & Gerhard, 2016; Helmi,
2020). Photometric surveys are particularly noteworthy in this regard, as they
can offer the largest and most representative sample of stars through
observing a much greater number of stars than spectroscopic surveys. As
demonstrated in this series of papers (An & Beers, 2020, 2021a, 2021b; An et
al., 2023, hereafter Paper I through Paper IV, respectively), photometric data
can serve as a valuable resource for constraining the fundamental properties
of stars, including effective temperatures ($T_{\rm eff}$) and metallicities
([Fe/H]), thereby facilitating the construction of comprehensive phase-space
distributions of Galactic stars on a global scale.

In this regard, the set of low-resolution spectra for $220$ million objects
obtained using the Blue and Red Photometer (BP/RP; hereafter XP) on board the
Gaia spacecraft is a new exciting addition to its data release 3 (DR3; Gaia
Collaboration et al., 2023b). As described in De Angeli et al. (2023), the XP
spectra were taken using the BP/RP photometers, which cover the wavelength
ranges [330, 680] nm and [640, 1050] nm, respectively, resulting in a
composite spectrum that spans from the near-ultraviolet (UV) to near-infrared
wavelengths. The XP data product contains information on almost all sources
brighter than $G=17.65$, making it an invaluable resource across the entire
sky, including low-latitude regions near the Galactic plane. Unlike the Gaia
mission’s medium-resolution Radial Velocity Spectrometer (RVS; Katz et al.,
2023), however, the spectral resolution of XP spectra varies between
$30<\lambda/\Delta\lambda<100$ in BP and $70<\lambda/\Delta\lambda<100$ in RP,
depending on the brightness of the source, the position on the detector, and
the wavelength (Carrasco et al., 2021), making it challenging to derive
accurate stellar parameters using traditional absorption-line analysis
methods. Nonetheless, XP spectra essentially provide approximately ten times
higher resolution than multiwavelength broadband photometry, allowing for the
precise estimation of fundamental stellar parameters.

In this study, we utilize the XP data to investigate the foreground reddening
of stars, by taking advantage of its finely sampled flux measurements across a
broad wavelength range. To accomplish this, it is essential to establish a
well-calibrated relationship that links stellar spectral energy distributions
with stellar parameters. To address this need, we employ a suite of
empirically calibrated isochrones in this study, as developed in our previous
paper of this series. Specifically, in Paper IV, we utilized a comprehensive
dataset consisting of precise observations of main-sequence (MS) stars in
Galactic clusters and those with spectroscopic measurements, spanning a broad
range of $T_{\rm eff}$ and [Fe/H], which facilitated the development of an
empirical calibration technique for theoretical synthetic spectra. Notably,
our calibration incorporated data from numerous photometric passbands obtained
from various photometric surveys, which was a crucial step toward creating
accurate spectral templates over a wide range of wavelengths.

Over the past year, a number of studies have also attempted to determine
stellar parameters from XP spectra. Andrae et al. (2023a) estimated the
metallicity and foreground extinction for around 500 million sources with XP
spectra using the General Stellar Parameterizer from Photometry (GSP-Phot).
However, the authors acknowledged biases in their metallicity estimates due to
the systematic differences between the purely theoretical models they have
adopted and the observed data. This situation has been improved in a
subsequent work (Andrae et al., 2023b), by using a machine-learning algorithm
to train empirical models based on spectroscopic data from the Apache Point
Observatory Galactic Evolution Experiment (APOGEE; Majewski et al., 2017) in
the Sloan Digital Sky Survey (SDSS; Abdurro’uf et al., 2022), augmented by
spectra for a list of metal-poor stars (Li et al., 2022).

In parallel to Andrae et al. (2023b), Zhang et al. (2023b) developed an
empirical forward model connecting the observed XP spectra to stellar
parameters, determined through medium-resolution spectroscopy from the Large
Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al.,
2012). By incorporating near-infrared photometry from the Two Micron All Sky
Survey (2MASS; Skrutskie et al., 2006) and the Wide-field Infrared Survey
Explorer (WISE; Wright et al., 2010) and employing Bayesian statistics, they
constrained astrophysical parameters ($T_{\rm eff}$, [Fe/H], distance, and
extinction) for the majority of stars with XP spectra, achieving uncertainties
of $\sim 0.15$ dex in [Fe/H] and $\sim 0.03$ mag in $E(B\,-\,V)$. Unlike the
approach adopted in Andrae et al. (2023a), who derived astrophysical
parameters based on purely theoretical models, both Andrae et al. (2023b) and
Zhang et al. (2023b) relied solely on empirical modeling to relate these
parameters to observables (XP spectra). This methodology is similar to earlier
studies that have utilized deterministic models (e.g., Ivezić et al., 2008) or
machine-learning techniques (e.g., Whitten et al., 2021).

In this work, we aim to improve the accuracy of extinction estimates by using
our empirically calibrated set of theoretical models. Our hybrid approach
distinguishes itself from the aforementioned studies in that it does not
exclusively rely on either empirical or theoretical relationships to establish
connections between observables and physical parameters. In Paper IV, we made
use of precisely known astrophysical parameters obtained from star clusters
and individual stars with spectroscopic measurements, similar to the approach
employed by Andrae et al. (2023b) or Zhang et al. (2023b). However, in
contrast to these approaches, we employed them to empirically calibrate
theoretical models, thereby allowing us to account for the inherent systematic
errors in these models while still benefiting from their predictive powers.
Given that the relationship linking observables to physical parameters serves
as the primary source of systematic errors in parameter estimation, our
calibrated models should help to mitigate such errors by combining the
strengths of both empirical and theoretical modeling.

The structure of this paper is organized as follows. In Section 2, we provide
a concise overview of XP spectra and present our selection of the XP MS
sample. Section 3 describes our methodology for estimating stellar parameters
using a set of four alternative models, allowing us to assess systematic
errors in our approach. Additionally, we offer a quality assessment of the XP
spectra through the application of best-fit models to individual samples. In
Section 4, we construct 3D models of foreground extinction while assuming a
constant total-to-selective extinction ratio ($R_{V}\equiv A_{V}/E(B\,-\,V)$),
and compare these models to previously published work in the literature, as
discussed in Section 5. In Section 6, we explore variations in the line-of-
sight $R_{V}$ values across the sky by allowing $R_{V}$ to vary during our
parameter search. Finally, we summarize our findings and discuss the
limitations of the current work in Section 7. Additionally, the appendices
include a validation of our technique and instructions for accessing our data
products.

##  2 Gaia XP Spectra

We obtained time-averaged mean XP spectra from the Gaia DR3 archive, which are
provided in the form of a continuous representation in basis
functions111http://cdn.gea.esac.esa.int/Gaia/gdr3/Spectroscopy/. To calibrate
the raw (internally calibrated, continuously represented) mean spectra, we
employed a Python library suite, GaiaXPy (De Angeli et al.,
2023)222https://gaia-dpci.github.io/GaiaXPy-website/. See Montegriffo et al.
(2023) for more information on the external calibration of the XP spectra. We
binned the spectra on a grid of [336, 1011] nm with a sampling interval of $9$
nm, which is wider than the default case ($2$ nm) to ensure efficient
processing while retaining sufficient information for our parameter estimates.

We matched the XP sources to the main Gaia DR3 catalog using the Gaia source
identifier (source_id). Our analysis focuses specifically on MS stars, as our
empirically calibrated isochrones are currently only applicable to the
parameter space covered by these stars. To ensure that our parameter estimates
are not sensitive to our age assumptions in the underlying models, we
restricted our analysis to objects with $3.5<M_{G}<7.5$ and $0.4<G_{\rm
BP}\,-\,G_{\rm RP}<2.0$, where $M_{G}$ is computed using Gaia parallaxes. In
this initial selection of the sample, we adopted foreground extinction from
our 3D extinction map in an iterative manner, with extinction coefficients in
Gaia passbands provided by Casagrande & VandenBerg (2018). Only those having
parallax uncertainties less than $20\%$ were kept, but the majority of the
objects in the XP catalog are sufficiently bright to have reliable parallax
measurements. The above selection reduced the total number of XP sources in
our sample from $220$ million to $\sim 80$ million.

##  3 Stellar Parameter Estimation

###  3.1 Strategy

Table 1: A Family of Solutions on Stellar Parameters Designation | Distance Prior | Empirical Model Calibration in Paper IV
---|---|---
Case A | Gaia DR3 | Stellar sequences + individual spectroscopic samples
Case B | Gaia DR3 | Stellar sequences
Case C | $\cdots$ | Stellar sequences + individual spectroscopic samples
Case D | $\cdots$ | Stellar sequences

As outlined in Table 1 and further explained below, we have generated four
distinct sets of parameter estimates in this study. They were derived from two
sets of empirical calibrations, taking into account whether they were computed
with or without the constraint from Gaia parallaxes. We refer to each of these
sets as Cases A–D, respectively. While Case A and Case B are preferred, due to
the additional constraints from Gaia parallaxes, we have included Case C and
Case D to explore systematic uncertainties in our parameter estimates by
analyzing the differences among these solutions. We also note that both Case A
and Case B are considered equally valid; therefore, we take the average of
these two solutions to derive our reddening map, while considering their
differences as a realistic measure of systematic errors in our models.

###  3.2 Recap: Empirically Calibrated Isochrones

In this series of papers, we derived stellar parameters, including mass
($M_{*}$) and [Fe/H], of a star by utilizing stellar isochrones with empirical
corrections (see below). Given our primary focus on the lower MS, $M_{*}$ is
directly correlated with $T_{\rm eff}$. Stellar isochrones play a crucial role
in our methodology, as they not only establish a fundamental
$M_{*}$-luminosity relation for a given age and chemical abundances but also
provide essential information on stellar radii necessary for placing synthetic
spectra onto the absolute flux scale. In instances where individual
spectroscopic measurements of $T_{\rm eff}$ and surface gravity ($\log{g}$)
for cluster stars in our calibration sample were unavailable, we relied on
absolute magnitudes (luminosities) of the isochrones to infer these
parameters. These isochrones were also employed to construct calibration
samples from diverse spectroscopic surveys, all tied to the same $T_{\rm eff}$
scale.

We continued to employ YREC evolutionary models (Sills et al., 2000), each
intricately linked to a synthetic spectrum generated by the MARCS atmosphere
models (Gustafsson et al., 2008). Our isochrones cover the range
$-3\leq{\rm[Fe/H]}\leq+0.4$, and for this range we assumed a monotonic
relationship between [Fe/H] and $\alpha$-element abundance ([$\alpha$/Fe]):
([Fe/H], [$\alpha$/Fe])$=(-3.0,0.4)$, $(-2.0,0.3)$, $(-1.0,0.3)$,
$(-0.75,0.25)$, $(-0.5,0.2)$, $(-0.3,0.0)$, $(0.4,0.0)$. Furthermore, our
models incorporate a linear relationship between age and metallicity, as
described by ([Fe/H], age)$=(-3.0,13\ {\rm Gyr})$, $(-1.2,13\ {\rm Gyr})$,
$(-0.3,4\ {\rm Gyr})$, $(0.4,4\ {\rm Gyr})$.

We adopted the linear relationship between [$\alpha$/Fe] and [Fe/H] owing to
the predominant influence of $T_{\rm eff}$ and overall metallicity
(represented by [Fe/H]) on stellar continuum flux, with the impact of surface
gravity being relatively modest. In contrast, the contribution of
[$\alpha$/Fe] is minimal. From an observational standpoint, although modern
spectroscopic surveys offer valuable [$\alpha$/Fe] information for a large
number of field stars, the challenge arises from the limited number of
benchmark clusters in our calibration sample with known [$\alpha$/Fe], making
it challenging to precisely determine the dependence of stellar flux on
[$\alpha$/Fe]. As described above, we have also opted for a simple linear
relationship of stellar age with [Fe/H] because of the relative absence of age
information for field stars, both in the calibration and when applying the
models to observational databases. To mitigate the influence of our assumed
age, especially near the MS turnoff, we focused our analysis on the lower MS
by restricting the sample to $M_{G}>3.5$ mag. Some stars near the MS turnoff
could be excluded owing to poor model fits resulting from inaccurate ages,
leading to a reduced number of stars in our final sample.

Our empirical corrections were formulated based on our working hypothesis that
YREC interior models predict accurate relations among luminosity, $T_{\rm
eff}$, and $M_{*}$, while deviations observed from the models are entirely
attributed to how fluxes are computed in the stellar atmosphere. In our
earlier calibration exercise (e.g., An et al., 2013), the empirical correction
procedure involved aligning our theoretical isochrones with the observed
sequences of both globular and open clusters in broadband photometry. This was
achieved by defining a table of magnitude corrections in each filter as a
function of $T_{\rm eff}$ and [Fe/H]. The resulting set of empirically
corrected isochrones has been extensively utilized in Papers I–III. In Paper
IV, we extended the calibration by directly adjusting synthetic spectra,
rather than making corrections on an individual filter basis. Additionally, we
expanded our calibration sample by incorporating extensive spectroscopic
observations of individual stars from the Sloan Extension for Galactic
Understanding and Exploration (SEGUE; Yanny et al., 2009; Rockosi et al.,
2022), APOGEE, and the Galactic Archaeology with HERMES (GALAH; Buder et al.,
2021).

In this study, we employed two sets of isochrones, each featuring distinct
empirical correction schemes. The first set, introduced in Paper IV, is based
on calibrations using both cluster sequences and individual stars with
spectroscopic measurements. The second set relies exclusively on the cluster
sequences, newly constructed as part of the current work, following the
correction procedure outlined in Paper IV. Corrections for both sets are
defined within the ranges $-3\leq{\rm[Fe/H]}\leq+0.4$ and $4000$ K $\leq
T_{\rm eff}\leq 7000$ K. However, the incorporation of spectroscopic samples
into the correction procedure depends on our assumed age (see above), leading
to systematic differences in the derived [Fe/H] for hot stars (see Figure 4 in
Paper IV). Inconsistent metallicity scales among diverse spectroscopic surveys
and the cluster data pose additional challenges. Nevertheless, we regard the
spectroscopic samples as valuable resources for extending our calibration, and
we incorporated both model sets in the subsequent analysis to address
potential systematic uncertainties in our methodology.

###  3.3 Parameter Estimation

We obtained the optimal parameter set, including $\\{$[Fe/H], $M_{*}$,
$E(B\,-\,V)$, and/or distance modulus $(m\,-\,M)_{0}$, and/or $R_{V}\\}$, for
each XP spectrum utilizing empirically corrected isochrones. To match the
resolution of the input XP spectra (see § 2), we adjusted synthetic spectra
accordingly. Given [Fe/H] and stellar mass, we calculated and minimized the
total $\chi^{2}$ of a fit employing the MPFIT package (Markwardt, 2009), for a
robust nonlinear least-squares curve fitting. This process yielded the best-
fitting foreground reddening, and/or $R_{V}$, and/or $(m\,-\,M)_{0}$. In our
base case (§ 4), we employed the average Galactic extinction curve from
Fitzpatrick (1999), corresponding to the standard value of $R_{V}\equiv
A_{V}/E(B\,-\,V)=3.1$, where $A_{V}$ represents the foreground extinction in
the $V$ band. Conversely, in scenarios allowing for the variation of $R_{V}$
in our parameter estimation (§ 6), we adopted the $R_{V}$-dependent extinction
curve model from Fitzpatrick et al. (2019).

Besides the two different methods for calibrating isochrones, there are two
distinct choices regarding whether Gaia parallaxes should be incorporated into
the fitting process as a prior. When our sample includes nearby stars with
reliable parallaxes, incorporating the Gaia prior results in a more stringent
constraint on parameter estimates compared to the other case. On the other
hand, the alternative approach, independently determining a stellar distance
through the fitting process, can produce a more internally consistent set of
parameters, notwithstanding potential systematic errors in the models. As
outlined in Table 1, we explore systematic uncertainties in our derived
parameters by examining all four families of solutions.

For each XP spectrum, we calculated the $\chi^{2}$ distribution within a model
grid, with grid intervals set at $\Delta{\rm[Fe/H]}=0.05$ dex and $\Delta
M_{*}=0.02\ M_{\odot}$. The process involved fitting a paraboloid to the
$\chi^{2}$ distribution to locate the minimum $\chi^{2}$, thereby determining
the best-fitting set of parameters for a star. To gauge parameter
uncertainties, we considered $\Delta\chi^{2}=3.53$ and $4.72$ from the minimum
$\chi^{2}$ for $3$ and $4$ degrees of freedom, respectively. When
incorporating Gaia parallaxes, we propagated parallax uncertainties and
combined them with the estimated uncertainties.

Figure 1: Model fits illustrating comparisons for various effective
temperatures ($T_{\rm eff}$, left) and metallicities ([Fe/H]; right). Best-
fitting estimates for $T_{\rm eff}$ and [Fe/H] are shown, with high-resolution
spectroscopic estimates from GALAH indicated in parentheses. The left panel
depicts all cases at solar metallicity.

Figure 1 illustrates our model fits to XP spectra for a subset of stars
characterized by minimal foreground reddening, $E(B\,-\,V)<0.05$ mag. The left
panel illustrates model comparisons (based on Case A) as a function of $T_{\rm
eff}$, while the right panel shows the case as a function of [Fe/H]. Although
individual absorption lines are not distinctly visible in the XP spectra,
molecular absorption bands from CH ($430$ nm) and MgH ($515$ nm) become
evident in low-$T_{\rm eff}$ or high-[Fe/H] stars. Notably, strong absorption
is clearly observed at wavelengths shorter than $400$ nm, primarily driven by
the Ca II H and K lines ($395$ nm). These features, along with the overall
spectral shape, exhibit variations with $T_{\rm eff}$ and [Fe/H], providing
valuable constraints in our model fitting. Nonetheless, although our models
exhibit a satisfactory overall fit above $360$ nm, systematic deviations
emerge at shorter wavelengths. In the following section, we delve into a
comprehensive discussion of these deviations, attributing them to inherent
characteristics of the XP spectra. It is important to note, however, that
these discrepancies do not markedly impact the overall fit, given the large
uncertainties associated with the affected pixels.

As our study primarily aims to derive the average reddening structure along
each line of sight, we fully utilized our estimated uncertainties to
effectively assign weights to each $E(B\,-\,V)$ measurement. Hence, we applied
a set of moderate-quality criteria to exclude only instances of stochastic
failures, requiring that $\chi^{2}_{\rm min}/\nu<5$, $\sigma(E(B\,-\,V))<0.3$
mag and $\sigma(\rm[Fe/H])<1.5$ dex, where $\nu$ represents the number of
degrees of freedom. We retained stars with metallicities falling within our
model grid ($-3<{\rm[Fe/H]}<0.4$). Furthermore, we required that the distance
moduli obtained from our pure spectrophotometric approach (Case C or Case D)
are within $1$ mag compared to those derived from Gaia parallaxes. This
criterion was applied to exclude objects that are less likely to belong to the
MS. Alongside the initial selection of the XP MS sample, as detailed in § 2,
these requirements served as the basis for our sample selection, resulting in
$2.6\times 10^{7}$, $2.3\times 10^{7}$, $2.8\times 10^{7}$, and $2.5\times
10^{7}$ stars in our final sample from Case A through Case D, respectively.

Our XP sample may also include a significant number of unresolved binaries
and/or blends, which exhibit systematically higher luminosities and/or redder
colors than single stars. For details on the impact of such binaries on
photometric metallicity estimates, we refer interested readers to An et al.
(2013). In the current study, we examined the potential bias by selecting
stars from the XP sample associated with two open clusters, M67 and NGC 2516,
both known to contain a noticeable number of unresolved binaries in the
$G_{\rm BP}\,-\,G_{\rm RP}$ color-magnitude diagram. Our analysis revealed
that their $E(B\,-\,V)$ values are consistently higher than those of resolved
stars. When Gaia parallaxes are employed in the parameter estimation (Case A
and Case B), these binaries exhibit systematically larger $E(B\,-\,V)$ values,
with differences of up to $0.15$ mag, while the purely spectrophotometric
solutions (Case C and Case D) show differences of less than $\Delta
E(B\,-\,V)\sim 0.05$ mag. The greater bias found in the former cases arises
from the inherent discrepancy between the actual (observed) and expected
(modeled) luminosities of the binaries. Nonetheless, the collective influence
of binary populations in the sample is anticipated to be modest: When
considering the ensemble of objects (both single and binary populations) in
these clusters, the overall bias ranges from $0.005$ to $0.04$ mag in Case A
and Case B and remains below $0.01$ mag in Case C and Case D.

###  3.4 Systematic Deviations and Uncertainties in XP Spectra

Figure 2: Model deviations and uncertainties in XP MS Sample. Top panel: mean
model deviations of the XP MS sample as a function of wavelength. Weighted
mean differences from all stars with $E(B\,-\,V)<0.1$ are shown by a
histogram, along with weighted standard deviations displayed by error bars.
Bottom panel: weighted standard deviations (error bars in the top panel),
overlaid on top of median fractional uncertainties in the XP spectra.

In the top panel of Figure 2, the mean model deviations of the XP MS sample as
a function of wavelength are shown by a histogram. It illustrates the weighted
mean differences of the best-match models from Case A for all stars with
$E(B\,-\,V)<0.1$ mag. Error bars indicate the weighted standard deviations of
the differences, but uncertainties in the mean differences are too small to
display, due to the large number of stars in the sample. Our best-match models
are within $3\%$ of the XP spectra at $400$ nm $<\lambda<1000$ nm, with an rms
difference of only $1.4\%$. The mild bump at $\sim 640$ nm is associated with
the boundary between the BP and RP spectra. However, the mean model deviation
gradually increases toward shorter wavelengths, reaching $40\%$ at the blue
edge of the spectra ($336$ nm). The large systematic deviation is accompanied
by a large scatter of the differences, as shown by the error bars.

The systematic errors of the XP spectra evident in the top panel of Figure 2
align with earlier validation work conducted by Montegriffo et al. (2023).
Their validation process incorporated a set of high-precision spectra from
$111$ stars in the Gaia spectrophotometric standard stars survey (Pancino et
al., 2012), with an additional $60$ stars from the ‘passband validation
library’ (Pancino et al., 2021). Furthermore, they made use of a dataset
consisting of $348$ stars from the Hubble Space Telescope’s Next Generation
Spectral Library (NGSL; Heap & Lindler, 2016). For wavelengths above $370$ nm,
the first two datasets displayed excellent agreement with the XP data, with
their smoothed median differences falling within $3\%$ (see their Figure 23).
Importantly, their work unveiled an intriguing comparison using the NGSL,
which was not part of the external calibration for the XP spectra, unlike the
other two datasets. In the range [400, 1000] nm, the smoothed median
difference from the NGSL was similar to the other two. However, at wavelengths
below $400$ nm, the deviation from the NGSL data exhibited a sharp increase,
reaching $\sim 10\%$ at $350$ nm—a value that escalates with higher $G_{\rm
BP}\,-\,G_{\rm RP}$ colors. The average deviation observed in the NGSL sample
closely mirrors our findings obtained through the calibrated models.
Therefore, it becomes evident that the substantial systematic deviation at
shorter wavelengths primarily arises from errors inherent in the XP spectra
themselves, rather than being a consequence of inaccurate isochrone models.

In the bottom panel of Figure 2, we compare the weighted standard deviations
derived in our study (red solid histogram) with the median flux uncertainties
extracted from the XP spectra (blue shaded histogram). It is noteworthy that,
within the range [400, 1000] nm, the above two uncertainty estimates closely
align with each other, suggesting that our parameter search is reasonably
accurate, given the flux uncertainties in the XP spectra. Below $360$ nm,
however, the observed scatter surpasses the estimated uncertainties, which may
indicate an underestimation of flux uncertainties in the XP spectra.

Table 2: Mean Model Deviations

Wavelength | Fractional DifferencesaaFractional differences of the flux, $(f_{\rm obs}-f_{\rm model})/f_{\rm model}$, as shown in the top panel of Figure 2, where $f_{\rm obs}$ and $f_{\rm model}$ represent the monochromatic flux in the XP spectra and the best-fitting model, respectively. | Median Frac.
---|---|---
(nm) | Weighted Mean | Weighted Std. | UncertaintybbMedian fractional uncertainties calculated using the cataloged values of flux uncertainties in the XP spectra.
336 | 0.412 | 0.476 | 0.266
345 | 0.087 | 0.189 | 0.162
354 | 0.128 | 0.240 | 0.185
363 | -0.075 | 0.103 | 0.165
372 | 0.058 | 0.121 | 0.175
381 | 0.065 | 0.169 | 0.251
390 | 0.018 | 0.112 | 0.190
399 | -0.034 | 0.056 | 0.087
408 | -0.011 | 0.043 | 0.058

Note. — Only a portion of this table is shown here to demonstrate its form and
content. A machine-readable version of the full table is available.

A machine-readable version of the data presented in Figure 2 is available in
Table 2. In the following analysis, we employed XP data spanning the
wavelength range of [336, 1011] nm. Despite the presence of significant
systematic errors in the blue and red edges of the XP spectra, alongside
underestimated flux uncertainties at these wavelengths, the substantial flux
uncertainties still underweight data below $400$ nm and above $1000$ nm,
respectively. We provide further evidence of this insensitivity in Appendix A,
where we demonstrate that including UV data ($\lambda<400$ nm) has minimal
impact on our parameter search, based on comparisons to spectroscopic
parameter estimates from GALAH and APOGEE.

##  4 3D Reddening Cube

###  4.1 Multiresolution HEALPix Scheme

Figure 3: Multiresolution HEALPix scheme. (a) The number density of the XP MS
sample in the Galactic coordinate system (Case A). The number of stars per
$0.2$ deg2 area is shown. (b) A multiresolution HEALPix map created using the
data shown in panel (a), in the same coordinate system. The map is composed of
four distinct sets of HEALPix cells, each with varying surface areas: $0.2$
deg2 ($28,152$ cells), $0.8$ deg2 ($24,802$), $3.4$ deg2 ($3872$), and $13.4$
deg2 ($114$). (c) Distribution of the number of stars per cell, derived from
the multiresolution grid scheme shown in panel (b).

Panel (a) of Figure 3 presents the number density distribution of the XP MS
sample in the Galactic coordinate system. For this visualization, we employed
the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) scheme (Górski
et al., 2005), which utilizes diamond-shaped cells, ensuring equal surface
area coverage on the celestial sphere. The number of cells and their
respective areas are determined by a single parameter, denoted as $N_{\rm
side}$. In this illustration, we selected $N_{\rm side}=128$, resulting in a
pixel area equivalent to $0.2$ deg2. Given that our analysis is confined to MS
dwarfs, our parameter estimation is limited to approximately $3$ kpc from the
Sun (see below), encompassing a substantial number of stars near the Galactic
plane.

Figure 4: Example of the multiresolution HEALPix scheme. Left panel: number
density distribution of the XP MS sample in the Galactic coordinate system in
a patch of the sky. The number of stars per unit HEALPix cell with an area of
$0.2$ deg2 ($N_{\rm side}=128$) is depicted. Right panel: mean reddening for
the same region of the sky as in the left panel. We note that the region with
high foreground extinction in the right panel (e.g.,
$(l,b)\sim(20\arcdeg,5\arcdeg$)) is underrepresented by the XP MS sample.
Larger multiresolution HEALPix cells are adopted in this area to maintain a
sufficient number of stars for modeling the foreground extinction distribution
as a function of distance.

Panel (a) of Figure 3 reveals a nonuniform number density distribution along
the Galactic plane, primarily due to substantial, patchy foreground extinction
in the magnitude-limited star catalog. A more detailed illustration of this
issue is provided in Figure 4. It displays the number density (left panel) and
foreground reddening (right panel) in a $30\arcdeg\times 30\arcdeg$ patch
centered around $(l,b)\sim(25\arcdeg,15\arcdeg$), which includes a high-
contrast strip of dense clouds near the Aquila Rift. As evidenced in the lower
right area, regions with high extinction (right panel) correspond to a lower
number density of stars from the XP MS sample (left panel), resulting in a
limited number of stars within each HEALPix cell for precise modeling of the
distance vs. reddening distribution.

To address this challenge, we employed multiresolution HEALPix maps developed
by Martinez-Castellanos et al. (2022), to effectively redistribute and adjust
the sizes of HEALPix cells, thereby ensuring a comparable number of stars
within each merged HEALPix cell. We initially divided the sky into HEALPix
cells with $N_{\rm side}=128$, but we adjusted the level of pixelization to
ultimately arrive at HEALPix cells as large as having $N_{\rm side}=16$, by
ensuring that each merged cell can accommodate a maximum of $1000$ stars. As
depicted in the right panel of Figure 4, the adoption of the multiresolution
HEALPix approach resulted in larger HEALPix cells in regions characterized by
high foreground extinction, while still maintaining a sufficiently large
number of stars.

The multiresolution mapping grid constructed in this manner is shown in panel
(b) of Figure 3. It consists of $56,940$ merged cells, of which $N_{\rm side}$
and the number of cells correspond to $N_{\rm side}$=$16$ ($114$ cells), $32$
($3872$), $64$ ($24,802$), and $128$ ($28,152$). The surface area ranges from
$0.2$ to $13.4$ deg2. As illustrated in panel (c), each multiresolution cell
contains between $\sim 200$ and $\sim 1000$ stars, with a median count of
$\sim 400$ stars, providing a sufficient number for modeling the reddening
structure, as demonstrated below. This multiresolution grid is employed
throughout the paper unless stated otherwise.

###  4.2 Modeling the Cumulative Reddening

Figure 5: Reddening of stars against distance modulus in a single,
multiresolution HEALPix cell located along $l=149.8\arcdeg$ and
$b=4.8\arcdeg$. The four panels display different sets of $E(B\,-\,V)$
measurements, Case A through Case D (see Table 1). Individual fits to the
simulated data using a double logistic function in equation (1) are
illustrated by green lines, and the best-fitting dust model is depicted by a
red solid line.

Figure 5 illustrates the line-of-sight distribution of reddening toward
$l=149.8\arcdeg$ and $b=4.8\arcdeg$, displaying the distribution of
$E(B\,-\,V)$ of individual stars as a function of distance modulus. Four
alternative solutions for $E(B\,-\,V)$ measurements are shown in separate
panels (Cases A–D). In this low-latitude example, the observed distribution in
all four cases exhibits characteristic jumps in the reddening at distances of
approximately $250$ and $800$ pc from the Sun ($(m\,-\,M)_{0}\approx 7$ and
$9.5$, respectively), representing nearby dust clouds and those possibly
associated with large-scale spiral arm(s), respectively.

Figure 5 also showcases systematic differences among the four alternative
solutions. First, it is apparent that the uncertainties in $E(B\,-\,V)$ are
significantly smaller in Case A and Case B compared to Case C and Case D,
mainly because of the additional constraint applied in the parameter
estimation. The smaller uncertainties in both axes lead to a tighter
correlation in Case A and Case B, while in Case C and Case D, the distribution
appears more dispersed. Despite this limitation, we note that, since
metallicity and reddening are two primary factors that influence the shape of
stellar spectra, the purely spectrophotometric solutions (Case C and Case D)
enable us to obtain an essentially more internally consistent result, provided
that our calibrated synthetic spectra accurately reflect the shape of stellar
spectra. There are also a larger number of nearby stars in Case C and Case D,
as opposed to the other cases. This can be attributed to the tight parallax
constraints imposed on the nearest stars in Case A and Case B, which, in turn,
could lead to the exclusion of such stars with substantial $\chi^{2}$ values
owing to a small systematic error in the models and/or the presence of
unresolved star pairs.

Secondly, Case A and Case C generally contain more stars than the cases using
isochrones with sequence-based calibration (Case B and Case D) by
approximately $10$–$20\%$. Furthermore, it is noteworthy in Figure 5 that the
$E(B\,-\,V)$ estimates derived from the sequence-based solutions (Case B and
Case D) exhibit values approximately $\Delta E(B\,-\,V)=0.02$ mag higher than
those obtained in the other cases, highlighting the distinction between the
two calibration methods. This discrepancy indicates that the difference
between the two calibration approaches extends beyond merely having additional
stars in the calibration, and reflects the unique information on stellar
properties provided by each calibration sample. Specifically, our assumptions
regarding stellar age, which are made for cluster sequences (with direct
information from an MS turnoff) and individual field stars (without reliable
age estimates), could have a profound impact on our calibration. This would
leave bluer stars near the MS turnoff unconstrained in Case B and Case D,
since isochrones with sequence-based calibration have a redder MS turnoff than
the other cases. The two approaches also exhibit a significant discrepancy in
the derived metallicities near the MS turnoff, as demonstrated in Paper IV.
Although a systematic difference between the two calibrations can be reduced
by selecting stars fainter than $M_{G}\approx 4.5$ mag, we kept stars with
$M_{G}>3.5$ mag as this severely limits the available stars in the XP MS
sample.

In order to construct a 3D reddening cube, we employed a logistic function
that effectively captures the observed distribution of $E(B\,-\,V)$ on a
purely empirical basis:

$E(B\,-\,V)=\epsilon_{0}+\sum_{i=1}^{2}A_{i}(1+e^{-k_{i}(\mu-\mu_{i})})^{-1},$
(1)

where $\mu$ indicates the distance modulus. Here we employed a double logistic
function with $i=\\{1,2\\}$ to represent the observed two-step features, such
as shown in Figure 5. Both $k_{i}$ and $\mu_{i}$ jointly determine the shape
of the cumulative $E(B\,-\,V)$ distribution, whereas $A_{i}$ serves as a
normalization factor.

In equation (1), $\epsilon_{0}$ signifies the zero-point offset in our
$E(B\,-\,V)$ measurements. To determine the systematic offsets for each
dataset, we employed $E(B\,-\,V)$ of stars at high Galactic latitudes
($|b|>60\arcdeg$) with minimal foreground extinction. The $E(B\,-\,V)$
distribution in these high-latitude regions shows systematic variations across
different datasets: we computed the mode of the distribution using a bin size
of $0.002$ mag and found $E(B\,-\,V)=\\{0.011$, $0.035$, $-0.021$, and
$0.001\\}$ mag for Case A through Case D, respectively. We took these values
as $\epsilon_{0}$ for each case and kept them constant throughout the
subsequent modeling.

We fit equation (1) to the observed distribution of $E(B\,-\,V)$ against
$(m\,-\,M)_{0}$ along each line of sight utilizing the MPFIT package. To
enhance the constraints on models involving a limited number of stars, we used
$E(B\,-\,V)$ measurements from $10^{4}$ stars in the surrounding area to
obtain shape parameters $\\{k_{1},\mu_{1},k_{2},\mu_{2}\\}$ in equation (1),
assuming that the structural properties of dust clouds remain constant over a
moderate angular scale (typically with a radius of $1\arcdeg$–$3\arcdeg$). For
these stars, we binned the $(m\,-\,M)_{0}$ vs. $E(B\,-\,V)$ data using bin
sizes of $0.10$ and $0.05$ mag for $(m\,-\,M)_{0}$ and $E(B\,-\,V)$,
respectively, and then used equation (1) to model the mode of the $E(B\,-\,V)$
distribution for a given $(m\,-\,M)_{0}$ with equal weights. This approach was
particularly useful for tracing the ridgeline at small $(m\,-\,M)_{0}$, given
the relatively small number of available stars in that range.

In the subsequent modeling in the central, multiresolution HEALPix cell, we
sought the best-fitting normalization factors $\\{A_{1},A_{2}\\}$ in equation
(1), after a couple of iterations with a $3\sigma$ outlier rejection. To take
into account heteroscedastic uncertainties in $E(B\,-\,V)$ and $(m\,-\,M)_{0}$
in both axes, we repeated all of the above fitting processes using a simple
Monte Carlo simulation, in which we assumed a normal distribution of
measurement uncertainties in these parameters based on our estimated
uncertainties and those from Gaia. The green lines in Figure 5 show individual
models obtained from $100$ realizations in such simulations, while the red
solid line represents their mean trend. To characterize this mean trend, we
consistently applied equation (1) to fit the median $E(B\,-\,V)$ from
individual models as a function of $(m\,-\,M)_{0}$.

A couple of caveats in our modeling should be noted, including the neglect of
correlations between $E(B\,-\,V)$ and $(m\,-\,M)_{0}$ and the internal scatter
originating from the finite size of our adopted HEALPix map scheme. However,
as evident in Figure 5, the double logistic function effectively captures the
two-layered dust structure in all families of solutions. Notably, the mean
models obtained from Case C and Case D exhibit a more gradual increase in
$E(B\,-\,V)$ at $(m\,-\,M)_{0}\sim 9.5$, due to the degraded precision in both
distance and reddening estimation in the pure spectrophotometric solutions.
While the range of $E(B\,-\,V)$ and $(m\,-\,M)_{0}$ remains comparable across
all cases, this comparison leads us to the conclusion that both Case A and
Case B, relying on Gaia parallaxes, outperform the other cases.

Figure 6: Same as Figure 5, but displaying Case A measurements for selected
lines of sight. Sight lines directed toward the Galactic plane are displayed
on the left, while additional sight lines situated at various Galactic
latitudes within the same longitude range are displayed on the right. Note
that distinct coverages are employed, with each one representing a single
multiresolution HEALPix cell.

Figure 6 displays the same $E(B\,-\,V)$ vs. $(m\,-\,M)_{0}$ distribution as
Figure 5 but for some lines of sight, providing a glimpse of the systematic
change in Galactic dust structure and the quality of our model fits. Only
$E(B\,-\,V)$ derived from Case A is shown. The left panels display data taken
near the Galactic plane, while the right panels show observations at
$|b|\approx 30\arcdeg$, where reddening is significantly less pronounced. Our
models generally exhibit good agreement with the observed structure in areas
with moderate cumulative extinction ($E(B\,-\,V)\lesssim 0.5$). Some regions
near the Galactic plane exhibit a two-step functional form similar to that
observed in Figure 5. In contrast, for most high-latitude samples, a single
logistic function suffices, which is unsurprising given that these lines of
sight do not traverse dense clouds in the local spiral arms.

Figure 7: Range of minimum (top panel) and maximum (bottom panel) distance
modulus determined by the number of stars accessible in Case A.

Figure 7 shows the range of the valid distance modulus of the sample employed
in Case A. Within the most densely obscured regions along the Galactic plane,
distances are confined to $d\lesssim 1$ kpc, reflecting the limited number of
stars within the magnitude-limited XP MS sample. Beyond $|b|>10\arcdeg$, our
sample encompasses volumes extending, on average, up to $\sim 3$ kpc, with
certain areas exhibiting even greater depth.

###  4.3 3D Reddening Map

Figure 8: Average reddening of the XP MS sample shown in a Mollweide
projection of the Galactic coordinate system, represented by the
multiresolution HEALPix. Each panel corresponds to one of the four different
solutions obtained in this study (Table 1). Figure 9: Differential reddening
of the XP MS sample based on Case A in selected distance bins. Reddening maps
are displayed in a Mollweide projection of the Galactic coordinate system
using the multiresolution HEALPix.

Figure 8 displays the distribution of average $E(B\,-\,V)$ in the Galactic
coordinate system for the four sets of parameter estimation utilized in our
analysis (Table 1). In all panels, a prominent large-scale structure of clouds
in the local volume is observed along the Galactic plane. However, finer
structures become more evident when these maps are divided into narrower
distance bins, as shown in Figure 9, on the same multiresolution grid as
depicted in panel (b) of Figure 3. It presents cross sections of our 3D
extinction map from Case A within selected distance intervals. The $\Delta
E(B\,-\,V)$ indicates an increment of reddening induced by dust clouds in a
given distance bin, calculated using our best-fitting models in equation (1)
along each line of sight.

The majority of dark clouds and cloud complexes identified in previous CO
surveys (e.g., Dame et al., 2001) are discernible in Figure 9; see also
Schlafly et al. (2014) for distances to individual clouds. Focusing solely on
the prominent structures among these, the extended molecular cloud complex in
the constellation of Ophiuchus ($\langle d\rangle\sim 100$ pc) becomes visible
in the nearest distance bin (panel (a)), situated toward the Galactic center
at $(l,b)\sim(0^{\circ},+20^{\circ})$. This complex overlaps with the Aquila
Rift ($\langle d\rangle\sim 100$ pc), centered at
$(l,b)\sim(30^{\circ},+5^{\circ})$, which gains prominence in the larger
distance bin (panel (b)). At $(l,b)\sim(300^{\circ},-15^{\circ})$, the
Chamaeleon cloud complex is visible (Voirin et al., 2018, $\langle
d\rangle\sim 200$ pc). In the direction of the Galactic anticenter, the
Taurus–Perseus–Auriga cloud complex ($\langle d\rangle\sim 100$–$300$ pc) is
also discernible within these volumes, located at
$(l,b)\sim(170^{\circ},-10^{\circ})$. In panels (b) and (c), the Orion Nebula
($\langle d\rangle\sim 400$ pc) centered approximately at
$(210^{\circ},-20^{\circ})$ and the Cepheus and Polaris Flares ($\langle
d\rangle\sim 400$ pc) at $(l,b)\sim(110^{\circ},+20^{\circ})$ come into play.
Beyond this point, a predominantly flattened structure of dust prevails along
the Galactic plane at greater distances.

Figure 10: Difference in modeled $E(B\,-\,V)$ between different cases of
parameter estimates (Table 1) in selected distance bins.

To gain insight into the systematic variations among different $E(B\,-\,V)$
maps, we compared different solutions with each other in Figure 10 at three
representative distance bins ($251$, $630$, and $1584$ pc). As a base case, we
chose Case A, and subtracted cumulative (rather than differential)
$E(B\,-\,V)$ of the other maps from it. Overall, the differences are minimal
when comparing Case A with Case B, suggesting that incorporating Gaia
parallaxes as a constraint in our parameter estimates leads to internal
consistency, irrespective of the choice of calibration. On the contrary, the
solutions derived from the pure spectrophotometric approach (Case C and Case
D) display substantial deviations from Case A. Specifically, Figure 10 reveals
noticeable discrepancies in the direction of the Aquila Rift and the Cepheus
Flare at a distance of $251$ pc. These discrepancies can be attributed to the
large uncertainties in $E(B\,-\,V)$ and $(m\,-\,M)_{0}$ for Case C and Case D,
leading to a smearing of the $E(B\,-\,V)$ vs. $(m\,-\,M)_{0}$ relation (see
Fig. 5).

In addition, large deviations are evident at $630$ and $1584$ pc along the
Galactic plane, where there is substantial foreground extinction. The
discrepancy in $E(B\,-\,V)$ is most evident with Case D, indicating that this
difference arises from a redder MS turnoff and the stronger correlation in our
parameter estimates when the Gaia’s prior on distance is absent. This may also
indicate that our modeling encounters challenges when attempting to accurately
reproduce the observed relationship between $E(B\,-\,V)$ and $(m\,-\,M)_{0}$
at low Galactic latitudes. Because our modeling approach involves deriving
structural parameters using stars in the adjacent HEALPix cells (§ 4.2), the
problem may be attributed to the presence of compact, dense molecular clouds,
causing $E(B\,-\,V)$ values to change rapidly over relatively short angular
and heliocentric distance intervals.

Despite the substantial differences illustrated in Figure 10, we note that the
models from Case A or B are favored over those of Case C and Case D, primarily
due to more precise measurements of $E(B\,-\,V)$ and $(m\,-\,M)_{0}$. Hence,
while the observed differences among different solutions serve as upper limits
for the uncertainties in our estimates, the systematic difference between Case
A and Case B may represent a more realistic measure of the true internal
accuracy of our $E(B\,-\,V)$ estimates. Within the range of $|b|\lesssim
10\arcdeg$, the $E(B\,-\,V)$ difference between Case A and Case B tends to
increase with distance. At a distance of $1584$ pc, approximately $1\%$ of the
area exhibits a deviation greater than $\Delta E(B\,-\,V)=0.1$ mag. Above this
latitude threshold, however, the median discrepancy between the two cases
remains minimal at all distances: after pixelating the maps into HEALPix cells
with a base resolution of $N_{\rm side}=128$, effective $2\sigma$
uncertainties, calculated from the $15.9$th and $84.1$st quartiles, amount to
$0.022$, $0.010$, and $0.012$ mag at distances of $251$, $630$, and $1584$ pc,
respectively. Taking these factors into consideration, it is evident that the
internal accuracy of our $E(B\,-\,V)$ models is approximately on the order of
$\sim 0.005$ mag at $|b|\gtrsim 10\arcdeg$.

Figure 11: Average $E(B\,-\,V)$ from Case A and Case B (Table 1) in the
Galactic coordinate system as a function of heliocentric distance.

Both Case A and Case B rely on parameter estimates with constraints from Gaia
parallaxes but diverge in terms of their model calibration. Nonetheless, their
extinction models are equally valid, as we do not favor any specific
calibration method for our isochrones (also see discussions in Paper IV).
Hence, we merged the extinction models from both cases by computing their
averages in bins of $\Delta(m\,-\,M)_{0}=0.02$ mag, spanning from
$(m\,-\,M)_{0}=5$ mag to $(m\,-\,M)_{0}=12.4$ mag. The resulting composite 3D
reddening map is displayed in Figure 11 in selected distance bins. We
incorporated the case at a distance of $100$ pc in panel (a) to ensure
completeness, although our models exhibit reduced accuracy for distances less
than $\sim 250$ pc owing to the limited number of nearby stars in the XP MS
sample.

##  5 Comparisons with Previous Extinction Maps

###  5.1 Schlegel et al. (1998)

Figure 12: Comparison to the $E(B\,-\,V)$ map of Schlegel et al. (1998). Top
panel: differences from the averaged $E(B\,-\,V)$ map at $d=1584$ pc (panel
(e) in Fig. 11), in the sense of their $E(B\,-\,V)$ minus our values. Bottom
panel: the solid line represents the median difference in $E(B\,-\,V)$, with
error bars indicating the standard deviation.

In Figure 12, we present a comparison between our averaged reddening map,
developed in the preceding section (Fig. 11), and the widely used extinction
map from Schlegel et al. (1998, hereafter SFD98). The latter map offers a
cumulative measure of the dust column along each line of sight, utilizing an
all-sky far-infrared dust emission map. In this comparison, we selected our
map at a distance of $1584$ pc, or $(m\,-\,M)_{0}=11$ mag, to ensure an
adequate assessment, leveraging reddening measurements from stars located at
sufficiently large distances. We applied a scaling factor of $0.86$ to the
SFD98 map to correct for the systematic difference from stellar locus
measurements (Schlafly & Finkbeiner, 2011). Given that the SFD98 dataset has a
finer spatial resolution ($\sim 6\arcmin$) in comparison to our baseline map
with $N_{\rm side}=128$ (having a mean spacing of $30\arcmin$), we computed
the median $E(B\,-\,V)$ from SFD98 for our XP MS sample within each HEALPix
cell.

The result, presented in Figure 12, reveals an excess in the SFD98 dust map
near the Galactic plane. This surplus is attributed to a significant volume of
dust located beyond the distance covered by our reddening map. This variation
is also clearly evident in the bottom panel, where the median difference in
$E(B\,-\,V)$ is illustrated against Galactic latitude, organized into bins of
$\Delta b=2.5\arcdeg$; it exceeds $\Delta E(B\,-\,V)=0.1$ mag for
$|b|<10\arcdeg$. A similar surplus in $E(B\,-\,V)$ is also evident in the
direction of the Magellanic Clouds, at $(l,b)=(280\arcdeg,-33\arcdeg)$ and
$(303\arcdeg,-44\arcdeg)$. While SFD98 appropriately accounts for the dust
present within these galaxies, our localized dust map does not encompass them.
Despite the inherent limitations in this comparison, however, Figure 12
demonstrates a rapid convergence of both maps for $|b|>10\arcdeg$.
Specifically, at $|b|=30\arcdeg$, the difference in $E(B\,-\,V)$ is merely
$0.003$ mag, and it remains within $0.01$ mag at larger $|b|$.

While this level of agreement already lends substantial support to the
accuracy of the zero-point in our $E(B\,-\,V)$ measurements, there is a
tantalizing possibility that the small residual discrepancy observed at high
Galactic latitudes ($\sim 0.005$ mag) can be attributed to a zero-point offset
within SFD98 itself. Based on more accurate measurements of dust temperature
and column density using Planck far-infrared observations, Planck
Collaboration et al. (2014) revealed that the original $E(B\,-\,V)$ estimates
in SFD98 are smaller by approximately $0.003$–$0.006$ mag in the diffuse
interstellar medium ($E(B\,-\,V)<0.04$ mag). This discrepancy is of the same
order as the difference observed at high Galactic latitudes in Figure 12.

###  5.2 Green et al. (2019)

Figure 13: Comparison to the $E(B\,-\,V)$ map of Green et al. (2019). Top
panels: comparisons are shown with respect to our averaged $E(B\,-\,V)$ maps
at selected distance bins ($d=251$, $630$, and $1584$ pc), in the sense of
their $E(B\,-\,V)$ minus our values. Bottom panels: median deviations as a
function of Galactic latitude in intervals of $2.5\arcdeg$. Error bars denote
their standard deviations.

A more rigorous comparison of our extinction map can be facilitated by
utilizing the 3D extinction map presented in Green et al. (2019). They
employed data from Pan-STARRS 1 (PS1) and 2MASS photometry in conjunction with
Gaia DR2 parallaxes and constructed a high-resolution extinction map by
subdividing HEALPix cells down to a level of $N_{\rm side}=1024$. The top
panels of Figure 13 illustrate a comparison between our averaged extinction
map and their ‘SFD-like’ $E(B\,-\,V)$ values in the Galactic coordinate system
at selected distance bins, $(m\,-\,M)_{0}=7$, $9$, and $11$ mag. For
consistency, we scaled their reddening values by a factor of $0.86$, similar
to our treatment of the SFD98 map. To ensure a consistent spatial resolution,
we degraded their high-resolution map to match our base HEALPix resolution of
$N_{\rm side}=128$.

As illustrated in Figure 13, the extinction map provided by Green et al.
(2019) is not available for regions with decl. less than $-30\arcdeg$ owing to
limitations in PS1 photometry. Nevertheless, it is evident that regions with
noticeable discrepancies are more prominent in areas characterized by
relatively high values of $E(B\,-\,V)$, although the discrepancy does not
follow a consistent pattern. In the bottom panels, the median deviations with
respect to Galactic latitude are shown as solid lines in bins of $2.5\arcdeg$.
To provide a sense of the dispersion, the standard deviations of these
differences are represented by error bars. At a distance of $251$ pc, the
median deviation reaches up to $\Delta E(B\,-\,V)=0.03$ mag, while it stays
below $0.015$ mag at $630$ pc. At a distance of $1584$ pc, a more pronounced
deviation is evident near the Galactic plane, but the median difference
remains modest at $|b|>10\arcdeg$, rarely exceeding $\Delta E(B\,-\,V)=0.02$
mag. The observed discrepancies may arise from differences in the adopted
distances (Gaia DR2 vs. DR3), extinction laws, fundamental stellar color
relationships, and the input data. Moreover, in regions characterized by high
extinction, the lower angular resolution of our map could also influence the
observed differences.

###  5.3 Zhang et al. (2023b)

Figure 14: Same as in Figure 13, but displaying a comparison with Zhang et al.
(2023b).

Figure 14 provides a comparison to more recent, all-sky extinction
measurements presented in Zhang et al. (2023b). Similar to ours, they utilized
XP spectra to derive line-of-sight extinction. Following their
recommendations, we applied a set of quality flags, including quality_flags
$<8$, teff_confidence $>0.5$, feh_confidence $>0.5$, $\sigma(T_{\rm eff})<500$
K, $\sigma({\rm[Fe/H]})<0.5$ dex, and a maximum $20\%$ uncertainty in
distance. In each of the HEALPix cells with $N_{\rm side}=128$, we computed a
weighted mean of their reddening estimates (denoted as ‘$E$’), which closely
aligns with the $E(B\,-\,V)$ scale used in our work, differing by $1\%$–$3\%$.
The comparisons depicted in Figure 14 focus on three specific distance bins,
where stars from the Zhang et al. catalog are situated within
$|\Delta(m\,-\,M)_{0}|<0.1$ mag of the given distance modulus. There are
numerous missing pixels in these maps (gray colors), which are primarily due
to the lack of measurements that pass the aforementioned quality flags in
Zhang et al. within this narrow distance interval.333The situation can
potentially be improved using extinction models, but Zhang et al. (2023b) only
made their estimates of individual stars’ reddening public. There is a recent
study (Edenhofer et al., 2023) that utilized Zhang et al. estimates to
construct a 3D extinction model, but it is limited to $1.25$ kpc.

Unlike the earlier comparison with Green et al. (2019), their reddening
estimates reveal substantial discrepancies from ours throughout all distance
bins. At high Galactic latitudes ($|b|>30\arcdeg$), the difference in
reddening amounts to up to $\Delta E(B\,-\,V)\approx 0.05$ mag. Given the
agreement of our reddening measurements at these latitudes with both SFD98 and
Green et al. (2019), the observed deviation is plausibly attributed to an
inherent offset in Zhang et al.. Further concern arises from the presence of
positive residuals originating from the Large and Small Magellanic Clouds,
even within these local volumes (panel (e)): $\Delta E(B\,-\,V)=0.019\pm
0.006$ mag and $0.044\pm 0.009$ mag, respectively. This suggests a potential
misclassification of bright stars belonging to the Magellanic Clouds as local
Milky Way stars in their study.

Figure 15: Comparison of $E(B\,-\,V)$ and [Fe/H] estimates between Zhang et
al. (2023b) and Case A for the identical XP sample from a region with a
$3\arcdeg$ radius, centered at $(l,b)=(20\arcdeg,20\arcdeg)$. An additional
[Fe/H] comparison is presented in the bottom right panel along a line of sight
with low extinction at $(l,b)=(300\arcdeg,45\arcdeg)$. Unity lines are
represented by dashed lines, while zero differences are displayed in the
bottom left panel.

A substantial body of evidence regarding the systematic offsets emerges from a
star-by-star comparison of $E(B\,-\,V)$ and [Fe/H] as illustrated in Figure
15. To facilitate this comparison, we restricted our analysis to stars within
a $3\arcdeg$ radius centered at $(l,b)=(20\arcdeg,20\arcdeg)$, utilizing the
same quality flags for data selection. As demonstrated in the top left panel,
the $E(B\,-\,V)$ values in Zhang et al. are approximately $20\%$ higher than
our own estimates. As illustrated in the top right panel, this discrepancy is
substantiated by the systematic deviation in their metallicity scale, which is
caused by a significant correlation between $E(B\,-\,V)$ and [Fe/H] in the
parameter estimation (bottom left panel). An additional comparison of [Fe/H]
is shown in the bottom right panel toward a line of sight at
$(300^{\circ},45^{\circ})$, indicating that the discrepancy in the metallicity
scale persists even at small foreground extinction ($E(B-V)<0.1$). In Appendix
A, we validate our parameter estimation by comparing it with metallicity
estimates from GALAH and APOGEE. The reasonably good agreement implies that
there is a substantial discrepancy in the $E(B\,-\,V)$ estimates in Zhang et
al.

##  6 Systematic Changes in the Extinction Laws

Figure 16: Mean metallicity estimates on the sky in the Galactic coordinate
system assuming $R_{V}=3.1$, for (a) Case A and (b) Case B.

Figure 16 illustrates the mean metallicity distribution of the XP MS sample
from Case A (panel (a)) and Case B (panel (b)). Apart from the scale
difference in metallicity between these maps ($\Delta{\rm[Fe/H]}\approx 0.05$
dex), there are small-scale systematic variations in [Fe/H]. For instance, the
region near $l\sim 130\arcdeg$ and $b\sim 5\arcdeg$ exhibits relatively higher
[Fe/H] compared to the surrounding areas, while an opposite trend is observed
near $l\approx 355\arcdeg$ and $b\approx 15\arcdeg$. These variations appear
to be more pronounced in regions with higher $E(B\,-\,V)$, but the spatial
trend does not appear to be entirely correlated with it (see Fig. 8). In
addition to these small-scale angular variations, stars in the first and
second Galactic quadrants generally exhibit higher mean metallicity compared
to the opposite side of the Galaxy. On the contrary, all-sky metallicity
estimates obtained using the General Stellar Parametriser from Spectroscopy
(GSP-Spec; Recio-Blanco et al., 2023) reveal a more symmetric distribution of
stars observed with the Gaia RVS (Gaia Collaboration et al., 2023a, see their
Fig. 4). Within a distance of $500$ pc, their metallicity distribution appears
nearly uniform, while exhibiting a flattened structure in the $1$ kpc $<d<2$
kpc bin.

The observed systematic trend can be interpreted as indicative of systematic
changes in the extinction curve along various lines of sight within the Milky
Way. In our previous analysis, we held $R_{V}$ fixed at $3.1$ (Fitzpatrick,
1999) during the parameter estimation. However, as initially studied by
Cardelli et al. (1989) and as corroborated by numerous subsequent studies,
including Fitzpatrick (1999), it is well established that extinction curves
can vary substantially, ranging from $R_{V}\sim 2$ to $6$. This variability
has been attributed to differences in the distributions of dust grain sizes
within the clouds along different lines of sight.

Indeed, the lines of sight displaying systematically lower [Fe/H] in Figure 16
are associated with dense molecular clouds, such as Ophiuchus, located near
$(l,b)\approx(355\arcdeg,15\arcdeg)$, Taurus $(170\arcdeg,-15\arcdeg)$, Orion
$(210\arcdeg,-20\arcdeg)$, and Chamaeleon $(300\arcdeg,-15\arcdeg)$. Previous
studies (e.g., Vrba & Rydgren, 1984; Cardelli et al., 1989; Whittet et al.,
2001) found that the denser parts of these regions exhibit higher $R_{V}$
($3.5$–$6$) than in the diffuse interstellar medium, which have been
attributed to the increase in grain size. On the contrary, stars within the
Cassiopeia and Cepheus clouds, encompassing the region
$110\arcdeg<l<140\arcdeg$ and $0\arcdeg<b<20\arcdeg$, consistently exhibit
higher [Fe/H] values. These stars, located on the periphery of the Cepheus
flare shell, an ancient supernova remnant, are associated with a cloud complex
that interacts with supernova shocks (Grenier et al., 1989). Earlier estimates
of $R_{V}$ in the adjacent region suggest a value of $R_{V}=2.9$, which is
slightly lower than the canonical value (Zdanavičius et al., 2002).

To account for the influence of varying extinction laws along different lines
of sight within the local volume, we employed the $R_{V}$ dependence formula
outlined in Fitzpatrick et al. (2019). However, the use of their relation
resulted in a mild reduction in the sample size with valid parameter estimates
that passed our quality criteria. For sight lines exhibiting moderate to heavy
foreground extinction ($\langle E(B\,-\,V)\rangle\sim 0.1$–$1$), the number of
stars remains within $5\%$–$10\%$ for both Case A and Case B. However, the
magnitude of the reduction reaches $\sim 5\%$–$30\%$ for Case C and Case D
owing to some stars falling outside the model grid. This is related to the
fact that the extinction curve at the standard $R_{V}=3.1$ in Fitzpatrick et
al. (2019) differs by up to $10\%$ on the short-wavelength side compared to
the original curve in Fitzpatrick (1999), which was utilized in the above
analysis. Hence, the derived extinction also increases by this amount when
adopting the standard $R_{V}=3.1$ curve from Fitzpatrick et al. (2019).
Despite this, the decrease in the number of stars indicates that the
extinction curve may not properly leverage the overall spectral energy
distribution of stars. Consequently, we concluded that the original extinction
curve presented by Fitzpatrick (1999) provides a more internally consistent
result for our set of calibrated models.

Therefore, we utilized the extinction law from Fitzpatrick (1999) at
$R_{V}=3.1$, while incorporating the $R_{V}$ dependence described in
Fitzpatrick et al. (2019). This was accomplished by comparing
$A_{\lambda}/E(B\,-\,V)$ at $R_{V}=3.1$ in both cases and then scaling the
$A_{\lambda}/E(B\,-\,V)$ relations at arbitrary $R_{V}$. We conducted a repeat
analysis of Case A and Case B by restricting our parameter search to $1\leq
R_{V}\leq 6$.444While Fitzpatrick et al. (2019) originally formulated the
$R_{V}$ dependence over the range $2.5\leq R_{V}\leq 6$, we extended the lower
limit in our modeling to account for uncertainties in our parameter estimates.
Nonetheless, in our mapping, most lines of sight exhibit $R_{V}$ values
greater than 2.5. We found $\epsilon_{0}=\\{0.015$, $0.029\\}$ mag for Case A
and Case B, respectively, for the zero-point in equation (1).

Figure 17: Mean and statistical properties of $R_{V}$ estimates from stars
within $(m\,-\,M)_{0}\leq 12$. Only regions with significant reddening
($E(B\,-\,V)\geq 0.1$ at $(m\,-\,M)_{0}=12$) are included. (a) Mean $R_{V}$
estimates from Case A and Case B. (b) Differences in $R_{V}$ between Case A
and Case B. (c) Half of the difference between the $15.9$th and $84.1$st
percentiles of $R_{V}$ from Case A ($\sigma(R_{V})$). (d) Significance of the
deviation of $R_{V}$ from the fiducial case ($R_{V}=3.1$) in relation to the
error in the mean $R_{V}$ ($\epsilon\equiv\sigma(R_{V})/\sqrt{N}$).

Panel (a) of Figure 17 shows the distribution of mean $R_{V}$ values estimated
using stars within $(m\,-\,M)_{0}\leq 12$ from both Case A and Case B. Here we
included only regions with significant reddening ($E(B\,-\,V)\geq 0.1$ at
$(m\,-\,M)_{0}=12$), as $R_{V}$ can be reliably obtained in such areas. As
shown in panel (b), the median $R_{V}$ estimates between the two cases differ
by $\Delta R_{V}\lesssim 0.1$ within $|b|\lesssim 30\arcdeg$, beyond which our
estimates become less accurate owing to the small amount of foreground
extinction. On the contrary, the variation in $R_{V}$ observed along each line
of sight is significant, as demonstrated in panel (c), where we present half
of the difference between the $15.9$th and $84.1$st percentiles of $R_{V}$
from Case A. In the vicinity of the Galactic plane, this effective $1\sigma$
deviation, $\sigma(R_{V})$, typically falls within a range of $0.5$–$1.0$ in
each multiresolution HEALPix cell. Nonetheless, the error in the mean $R_{V}$,
computed as $\epsilon\equiv\sigma(R_{V})/\sqrt{N}$, where $N$ indicates the
number of stars in each cell, is still about $3$ times larger than the
difference between Case A and Case B (panel (b)). This suggests that a
significant dispersion in $R_{V}$ can be attributed, at least in part, to
actual star-by-star variations.

The mean $R_{V}$ values shown in panel (a) in Figure 17 range from $\approx 3$
to $4.5$, with a bulk average of $3.7$, which is larger than the canonical
$R_{V}=3.1$ (e.g., Cardelli et al., 1989; Fitzpatrick et al., 2019). If one
takes $\epsilon$ as an effective $1\sigma$ uncertainty, our estimates reveal a
significant change in $R_{V}$ across the sky (panel (d)). The systematic trend
between metallicity in Figure 16 (evaluated at $R_{V}=3.1$) and $R_{V}$ in
Figure 17 is not entirely correlated with each other. Nonetheless, as noted
above, areas where metallicities appear to be underestimated tend to have
relatively high $R_{V}$ values (such as the Ophiuchi cloud), while the
opposite trend is observed in the low-$R_{V}$ region (most strikingly in the
Cassiopeia and Cepheus cloud complexes).

In earlier studies, Schlafly et al. (2016) and Zhang et al. (2023a)
investigated large-scale variations in $R_{V}$ using spectroscopic data from
APOGEE and LAMOST, respectively. Their 2D mapping is inherently limited by
target availability, with the latter study covering a significantly larger
area, yet restricted to the northern celestial hemisphere. Nevertheless, it is
possible to make a qualitative comparison for the specific regions that are
mentioned above. First, the regions with the highest $R_{V}$ in Zhang et al.
(2023a) exhibit $R_{V}\sim 4$ and are observed near the Ophiuchi cloud, in
agreement with our findings. Their map also reveals elevated $R_{V}$ values in
patchy regions near the Orion and Taurus cloud complexes, consistent with our
findings. Moreover, the Cassiopeia and Cepheus clouds in Schlafly et al.
(2016) exhibit low $R_{V}$ ($\sim 3.0$), and there is also an indication of
low $R_{V}$ in Zhang et al. (2023a) near the border of the cloud complexes, a
region not fully covered in their study. However, we note that their samples
include giants and extend beyond our distance limit, reaching out to $4$–$5$
kpc. Consequently, their projected $R_{V}$ map may be influenced by distant
stars that are not included in our sample.

Figure 18: Mean metallicity estimates from varying $R_{V}$ in the parameter
search, for (a) Case A and (b) Case B.

The revised metallicity maps, which are based on varying $R_{V}$, are
presented in Figure 18. In contrast to Figure 16, they reveal enhanced mean
metallicities toward the Galactic center and are more symmetric with respect
to both the Galactic plane and the Galactic prime meridian in both Case A
(panel (a)) and Case B (panel (b)). In Figure 18, the highest-metallicity
regions are located near the Aquila Rift ($20\arcdeg$,$15\arcdeg$) and the
Lupus dark clouds ($340\arcdeg$,$10\arcdeg$). Along with the high-metallicity
regions near Taurus and Orion, they appear to follow the Gould Belt, a
projected ring of a large number of bright stars, OB associations, and a
significant amount of interstellar material, inclined by $\sim$20∘ with
respect to the Galactic plane (Zari et al., 2018, and references therein).
Although this ring-like structure is not immediately evident in the RVS sample
(Gaia Collaboration et al., 2023a), the overall mean metallicity distribution
in Figure 18, characterized by higher metallicity along the Galactic plane,
aligns with their findings, suggesting that the approach presented in this
section produces more physically plausible results.

Figure 19: Differences in $A_{V}$ derived from Case A when assuming a fixed
value of $R_{V}=3.1$, compared to using a median $R_{V}$ (see text). The
comparisons are illustrated at three representative distances.

Despite this, the $E(B\,-\,V)$ data cube constructed using this approach does
not necessarily surpass our base case with $R_{V}=3.1$ due to the limited
precision of our $R_{V}$ measurements for individual stars, which would
inevitably increase the uncertainty in the foreground extinction.
Alternatively, one can compute a median value for $R_{V}$ along each line of
sight and utilize it in estimating foreground extinction ($A_{V}$). The
comparisons, presented in Figure 19 at three select distances, highlight
differences in $A_{V}$ for two scenarios. The first is derived from Case A,
assuming a fixed value of $R_{V}=3.1$, while the second assumes a median
$R_{V}$ computed from each cell. As demonstrated here, the most substantial
discrepancies become evident in regions where $R_{V}$ notably deviates from
the canonical value (as depicted in Figure 17), effectively tracing regions
with high dust columns. In regions with the highest $E(B\,-\,V)$ values, the
difference in $A_{V}$ can reach a few tenths of a magnitude.

##  7 Summary and Caveats

In this study, we utilized the extensive dataset of low-resolution spectra
from Gaia DR3 to construct a 3D map of foreground dust reddening across the
entire sky. This endeavor was made possible through the application of our
newly developed, empirically calibrated stellar isochrones, as detailed in
Paper IV. Our approach involved modeling foreground extinction using a
straightforward, two-step function along each line of sight, yielding results
that exhibit a good agreement with widely used extinction maps by Schlegel et
al. (1998) and Green et al. (2019). However, our map reveals a significant
discrepancy, particularly at low Galactic latitudes, compared to the study by
Zhang et al. (2023b), who provided reddening estimates based on the same XP
spectra but through an empirical forward modeling technique. We attribute
these differences to systematic errors in their parameter estimation, a
conclusion further supported by the alignment of our metallicity scale with
high-resolution spectroscopic results from the GALAH survey and, to some
extent, SDSS/APOGEE (see Appendix A).

While our mapping of foreground extinction provides comprehensive coverage in
both the northern and southern hemispheres, there are certain limitations in
our approach. First, our current results rely on calibrated isochrones for MS
stars, which restricts this study to a relatively nearby volume compared to
using giant stars. On average, our map extends up to $3$ kpc from the Sun, but
the physical depth of the reddening cube varies across the sky owing to our
sample being primarily constrained by the brightness limit of the XP spectra
and the availability of accurate Gaia parallaxes. Consequently, some areas
with nearby dust clouds have a limited distance range in our data cube.

Secondly, our reddening cube is constructed by modeling cumulative reddening
as a function of distance using a simple two-step logistic function along each
line of sight. This choice of model is primarily empirical, driven by the
observed distribution of stars, rather than being motivated by the physical
and structural properties of interstellar dust. Due to the trade-off between
the number of stars within each HEALPix cell and the spatial resolution of the
reddening map, our chosen mapping scheme did not offer HEALPix cell sizes
large enough to generate a purely empirical representation of cumulative
extinction.

Lastly, we presented two different cases of a reddening cube, one assuming a
fixed value of $R_{V}$ and the other varying $R_{V}$ in our parameter search.
Although the latter is more physically motivated, the limited precision of our
measurements prevents us from accurately tracing fine-structure variations of
$R_{V}$ along each line of sight. We have provided only mean $R_{V}$ estimates
along each sight line, but adopting these mean $R_{V}$ values for subsequent
analyses may potentially introduce bias into the results if the true $R_{V}$
of stars significantly deviates from its mean value for a given sight line.

The 3D reddening map presented in this paper holds significant value for our
ongoing research within this series of papers. In our previous work, we
utilized photometric data to infer the metallicities of stars, but the low-
latitude regions near the Galactic plane posed challenges, as they were
largely obscured by substantial extinction in optical imaging surveys. High-
resolution spectroscopic surveys are particularly noteworthy in this regard,
as they can offer the most detailed information on individual stars, such as
their 3D kinematics and complete elemental abundances, even in the presence of
large foreground extinction. Nevertheless, due to inherent limitations in
spectroscopic surveys, such as the relatively stronger target selection bias
and a smaller number of stars compared to photometry, it is crucial to
establish a robust foundation for understanding the various Galactic stellar
populations through the development of photometric metallicity mapping, as
demonstrated in our series of works. Eventually, this will allow us to delve
deeper into spatial distributions, dynamics, and chemical properties of the
various (sub)structures that make up the Milky Way and gain a better
understanding of its formation and evolution history.

The authors gratefully acknowledge Young Sun Lee for his careful review of an
early draft of the manuscript and for providing valuable feedback. D.A. also
thanks Seo-Won Chang for meticulously examining the extinction cubes in the
GitHub library. D.A. acknowledges support provided by the National Research
Foundation of Korea grant funded by the Ministry of Science and ICT (No.
2021R1A2C1004117). T.C.B. acknowledges partial support from grant PHY
14-30152, Physics Frontier Center/JINA Center for the Evolution of the
Elements (JINA-CEE), and from OISE-1927130: The International Research Network
for Nuclear Astrophysics (IReNA), awarded by the US National Science
Foundation. A.C. is supported by a Brinson Prize Fellowship at the University
of Chicago/KICP. 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 work has made use of the Python
package GaiaXPy, developed and maintained by members of the Gaia Data
Processing and Analysis Consortium (DPAC), and in particular, Coordination
Unit 5 (CU5), and the Data Processing Centre located at the Institute of
Astronomy, Cambridge, UK (DPCI). Software: GNU parallel (Tange 2018).

##  Appendix A Comparisons to High-resolution Spectroscopic Parameter
Estimates

Figure 20: Comparisons of metallicity estimates in this study with high-
resolution spectroscopic measurements from GALAH (top panels) and APOGEE
(bottom panels). In the right panels, the same suite of comparisons is shown,
but when our parameter estimates are restricted to $390$–$993$ nm, excluding
UV portion of the spectra. The solid line represents a unity line.

In this appendix, we present a comparison between our metallicity estimates
derived from XP spectra and spectroscopic measurements obtained through
traditional line analysis methods. Figure 20 illustrates these comparisons for
two distinct datasets: GALAH DR3 (Buder et al., 2021) and APOGEE/SDSS DR17
(Majewski et al., 2017; Abdurro’uf et al., 2022), both based on high-
resolution spectroscopy ($R\approx 28,000$ and $R\approx 22,500$,
respectively). To conduct this comparison, we performed a cross-match between
our catalogs and those of GALAH and APOGEE, using a search radius of
$1\arcsec$. We specifically included our metallicity measurements with the
following criteria: $\sigma(E(B\,-\,V))<0.3$, $3.5<M_{r}<7.5$, and
$\sigma(\rm[Fe/H])<1.5$ as a minimum quality requirement, in alignment with
the primary sample criteria for our work. In both datasets, we limited our
analysis to stars with $\log{g}>3.5$ and $E(B\,-\,V)<0.1$ in Schlegel et al.
(1998). Additionally, we applied specific selection cuts for each
spectroscopic dataset, ensuring that aspcapflag is not set for APOGEE. For the
GALAH sample, we required that the stellar parameter quality flag (flag_sp)
and the overall iron abundance quality flag (flag_fe_h) are not set.

In Figure 20, the first and second columns present data corresponding to Case
A and Case B, respectively, in our study. These columns illustrate that the
systematic differences between the two calibrations lead to an approximately
$0.1$ dex difference in [Fe/H]. Despite these systematic distinctions between
the two cases, these comparisons reveal that our metallicity estimates align
well with those obtained from the GALAH survey, with differences within
$\Delta{\rm[Fe/H]}\lesssim 0.1$ at ${\rm[Fe/H]}>-1$. While our estimates also
exhibit reasonable agreement with the APOGEE dataset, it is worth noting that
in Case A there is a tendency for the underestimation of [Fe/H] by
approximately $0.15$ dex. Such a scale difference in metallicity is not
uncommon when comparing different spectroscopic abundance analyses.
Considering the inherent systematics between these two spectroscopic samples,
we interpret these results as satisfactory evidence that our metallicity scale
remains reasonable and in line with existing spectroscopic measurements.

The final column in Figure 20 illustrates the effect of incorporating UV flux
in the XP spectra on our metallicity measurements, specifically for the
comparisons with Case A. In essence, there are negligible differences between
the first and third columns, indicating that the exclusion of UV data has
minimal impact on the parameter estimates. This outcome aligns with
expectations from the fact that the uncertainty in flux measurements increases
at wavelengths shorter than $400$ nm. While it is widely acknowledged that UV
flux plays a crucial role in deriving stellar parameters from broadband
photometry, our results indicate that even low-resolution spectra at
$\lambda\gtrsim 400$ nm can offer a wealth of information about photospheric
metal abundance and foreground extinction.

##  Appendix B Comparisons to Cluster $E(B\,-\,V)$ Estimates in the Literature

Figure 21: Comparisons with $E(B\,-\,V)$ estimates from the literature. (a)
Open cluster samples with fundamental parameters from An et al. (2007, 2019).
(b) A homogeneous set of fundamental parameters for open cluster samples from
Dias et al. (2021). (c) Compilation of cluster $E(B\,-\,V)$ for Galactic
globular clusters as reported in Harris (1996). In panel (c), only the lower
limit in our $E(B\,-\,V)$ estimates is displayed, as the cluster distances
exceed the maximum threshold of our $E(B\,-\,V)$ cube. In all panels, the red
solid line corresponds to the best-fitting regression to the data, with a
zero-point and a slope displayed at the top, whereas dotted lines represent
the unity and zero lines.

Figure 21 provides a comparison of our $E(B\,-\,V)$ estimates, obtained from
the averaged reddening cube at $R_{V}=3.1$, with those from the literature. In
panel (a), we display open cluster samples with precise fundamental parameters
in An et al. (2007, 2019). These clusters, namely Praesepe, the Pleiades, M67,
NGC 6791, and NGC 2516 from left to right, have been essential resources in
our prior tests for empirically calibrated isochrones in the
Johnson–Cousins–2MASS photometric systems. Additionally, M67 and NGC 6791 have
played a pivotal role in the empirical calibration of isochrones utilized in
this work (Paper III). Panel (b) displays a comparison to a homogeneous set of
fundamental parameter estimates for open cluster samples in Dias et al.
(2021). Their estimates relied on fitting isochrones to the observed $G_{\rm
BP}$ and $G_{\rm RP}$ photometric data in Gaia DR2. We assumed $R_{V}=3.1$ to
convert their $A_{V}$ to $E(B\,-\,V)$, and we restricted their sample to those
at $|b|>10\arcdeg$.

In the right panel, we present a comparison to a compilation of cluster
parameters for Galactic globular clusters from Harris (1996). In this
comparison, we display our $E(B\,-\,V)$ estimates as lower limits because the
cluster distances extend beyond the maximum threshold of our $E(B\,-\,V)$
model. Nevertheless, our $E(B\,-\,V)$ lower limits exhibit a strong
correlation with the values from Harris (1996). This is because it is unlikely
to encounter additional dust clouds beyond our distance limit at
$|b|>10\arcdeg$.

The error bars in both panels (a) and (b) of Figure 21 reflect the combined
uncertainties, calculated as the quadrature sum of two components: half of the
difference between Case A and Case B solutions and the propagated uncertainty
from the adopted distance. The solid line in each comparison represents a
linear fit to the data after applying a $3\sigma$ clipping, accounting for
errors in both axes. In the right panel, we also made an attempt to fit the
data using the maximum $E(B\,-\,V)$ value within our models along each line of
sight (i.e., lower limits in $E(B\,-\,V)$ for distant clusters), as
demonstrated by the solid line. The resulting best-fitting lines reveal a
zero-point offset of $0.01$–$0.02$ mag. Our $E(B\,-\,V)$ estimates for open
clusters are systematically smaller by $5\%$ compared to the values in Dias et
al. (2021).

While Dias et al. (2021) also adopted the extinction law by Fitzpatrick et al.
(2019), as in this study, the systematic difference is further increased to
approximately $16\%$ when comparing our estimates to the reddening cube with
varying $R_{V}$. Therefore, other factors, such as the correlation between
metallicity and reddening estimates, may be responsible for the observed
systematic difference. Alternatively, a finite size of the angular resolution
in our reddening map could be the cause if the open clusters are
preferentially situated within extra dust clouds. In any case, a star-by-star
comparison, such as shown in Figure 15, would be invaluable in determining the
origin of the systematic difference in $E(B\,-\,V)$.

##  Appendix C Accessing Extinction Data Cubes

The primary data products resulting from this study can be accessed through
the following link.555https://github.com/deokkeunan/Galactic-extinction-map
Within this repository, we have made available two sets of data cubes. The
first set, described in Section 4, corresponds to $E(B\,-\,V)$ cubes derived
for $R_{V}=3.1$. The second set, discussed in Section 6, includes data cubes
resulting from a parameter search involving varying values of $R_{V}$. Each of
these data cubes contains the mean $E(B\,-\,V)$ values from Case A and Case B,
using the multiresolution HEALPix scheme in the Galactic coordinate system.
Additionally, we have included an example script to facilitate the utilization
of these data sets.

## References

  * Abdurro’uf et al. (2022) Abdurro’uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35. doi:10.3847/1538-4365/ac4414
  * 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. (2023) An, D., Beers, T. C., Lee, Y. S., et al. 2023, ApJ, 952, 66. doi:10.3847/1538-4357/acd5cb (Paper IV)
  * 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. (2019) An, D., Pinsonneault, M. H., Terndrup, D. M., et al. 2019, ApJ, 879, 81. doi:10.3847/1538-4357/ab23ed
  * An et al. (2007) An, D., Terndrup, D. M., Pinsonneault, M. H., et al. 2007, ApJ, 655, 233. doi:10.1086/509653
  * Andrae et al. (2023a) Andrae, R., Fouesneau, M., Sordo, R., et al. 2023a, A&A, 674, A27. doi:10.1051/0004-6361/202243462
  * Andrae et al. (2023b) Andrae, R., Rix, H.-W., & Chandra, V. 2023b, ApJS, 267, 8. doi:10.3847/1538-4365/acd53e
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529. doi:10.1146/annurev-astro-081915-023441
  * Buder et al. (2021) Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150. doi:10.1093/mnras/stab1242
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245. doi:10.1086/167900
  * Carrasco et al. (2021) Carrasco, J. M., Weiler, M., Jordi, C., et al. 2021, A&A, 652, A86. doi:10.1051/0004-6361/202141249
  * Casagrande & VandenBerg (2018) Casagrande, L. & VandenBerg, D. A. 2018, MNRAS, 479, L102. doi:10.1093/mnrasl/sly104
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197. doi:10.1088/1674-4527/12/9/003
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792. doi:10.1086/318388
  * De Angeli et al. (2023) De Angeli, F., Weiler, M., Montegriffo, P., et al. 2023, A&A, 674, A2. doi:10.1051/0004-6361/202243680
  * Dias et al. (2021) Dias, W. S., Monteiro, H., Moitinho, A., et al. 2021, MNRAS, 504, 356. doi:10.1093/mnras/stab770
  * Edenhofer et al. (2023) Edenhofer, G., Zucker, C., Frank, P., et al. 2023, arXiv:2308.01295. doi:10.48550/arXiv.2308.01295
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63. doi:10.1086/316293
  * Fitzpatrick et al. (2019) Fitzpatrick, E. L., Massa, D., Gordon, K. D., et al. 2019, ApJ, 886, 108. doi:10.3847/1538-4357/ab4c3a
  * Gaia Collaboration et al. (2023a) Gaia Collaboration, Recio-Blanco, A., Kordopatis, G., et al. 2023a, A&A, 674, A38. doi:10.1051/0004-6361/202243511
  * Gaia Collaboration et al. (2023b) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023b, A&A, 674, A1. doi:10.1051/0004-6361/202243940
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759. doi:10.1086/427976
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., et al. 2019, ApJ, 887, 93. doi:10.3847/1538-4357/ab5362
  * Grenier et al. (1989) Grenier, I. A., Lebrun, F., Arnaud, M., et al. 1989, ApJ, 347, 231. doi:10.1086/168112
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951. doi:10.1051/0004-6361:200809724
  * Harris (1996) Harris, W. E. 1996, AJ, 112, 1487. doi:10.1086/118116
  * Heap & Lindler (2016) Heap, S. R. & Lindler, D. 2016, The Science of Calibration, 503, 211
  * Helmi (2020) Helmi, A. 2020, ARA&A, 58, 205. doi:10.1146/annurev-astro-032620-021917
  * Ivezić et al. (2012) Ivezić, Ž., Beers, T. C., & Jurić, M. 2012, ARA&A, 50, 251. doi:10.1146/annurev-astro-081811-125504
  * Ivezić et al. (2008) Ivezić, Ž., Sesar, B., Jurić, M., et al. 2008, ApJ, 684, 287. doi:10.1086/589678
  * Katz et al. (2023) Katz, D., Sartoretti, P., Guerrier, A., et al. 2023, A&A, 674, A5. doi:10.1051/0004-6361/202244220
  * Li et al. (2022) Li, H., Aoki, W., Matsuno, T., et al. 2022, ApJ, 931, 147. doi:10.3847/1538-4357/ac6514
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94. doi:10.3847/1538-3881/aa784d
  * Markwardt (2009) Markwardt C. B. 2009 Astronomical Data Analysis Software and Systems XVIII ed D. A. Bohlender, D. Durand and P. Dowler (San Francisco, CA: ASP) 251
  * Martinez-Castellanos et al. (2022) Martinez-Castellanos, I., Singer, L. P., Burns, E., et al. 2022, AJ, 163, 259. doi:10.3847/1538-3881/ac6260
  * Montegriffo et al. (2023) Montegriffo, P., De Angeli, F., Andrae, R., et al. 2023, A&A, 674, A3. doi:10.1051/0004-6361/202243880
  * Pancino et al. (2012) Pancino, E., Altavilla, G., Marinoni, S., et al. 2012, MNRAS, 426, 1767. doi:10.1111/j.1365-2966.2012.21766.x
  * Pancino et al. (2021) Pancino, E., Sanna, N., Altavilla, G., et al. 2021, MNRAS, 503, 3660. doi:10.1093/mnras/stab766
  * Planck Collaboration et al. (2014) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11. doi:10.1051/0004-6361/201323195
  * Recio-Blanco et al. (2023) Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2023, A&A, 674, A29. doi:10.1051/0004-6361/202243750
  * 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
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525. doi:10.1086/305772
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103. doi:10.1088/0004-637X/737/2/103
  * Schlafly et al. (2014) Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29. doi:10.1088/0004-637X/786/1/29
  * Schlafly et al. (2016) Schlafly, E. F., Meisner, A. M., Stutz, A. M., et al. 2016, ApJ, 821, 78. doi:10.3847/0004-637X/821/2/78
  * Sills et al. (2000) Sills, A., Pinsonneault, M. H., & Terndrup, D. M. 2000, ApJ, 534, 335. doi:10.1086/308739
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163. doi:10.1086/498708
  * Tange (2018) Tange, O. 2018, GNU Parallel 2018, v1, Zenodo, doi: 10.5281/zenodo.1146014
  * Voirin et al. (2018) Voirin, J., Manara, C. F., & Prusti, T. 2018, A&A, 610, A64. doi:10.1051/0004-6361/201731153
  * Vrba & Rydgren (1984) Vrba, F. J. & Rydgren, A. E. 1984, ApJ, 283, 123. doi:10.1086/162281
  * Whitten et al. (2021) Whitten, D. D., Placco, V. M., Beers, T. C., et al. 2021, ApJ, 912, 147. doi:10.3847/1538-4357/abee7e
  * Whittet et al. (2001) Whittet, D. C. B., Gerakines, P. A., Hough, J. H., et al. 2001, ApJ, 547, 872. doi:10.1086/318421
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868. doi:10.1088/0004-6256/140/6/1868
  * 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
  * Zari et al. (2018) Zari, E., Hashemi, H., Brown, A. G. A., et al. 2018, A&A, 620, A172. doi:10.1051/0004-6361/201834150
  * Zdanavičius et al. (2002) Zdanavičius, J., Straižys, V., & Corbally, C. J. 2002, A&A, 392, 295. doi:10.1051/0004-6361:20020944
  * Zhang et al. (2023a) Zhang, R., Yuan, H., & Chen, B. 2023a, ApJS, 269, 6. doi:10.3847/1538-4365/acf764
  * Zhang et al. (2023b) Zhang, X., Green, G. M., & Rix, H.-W. 2023b, MNRAS, 524, 1855. doi:10.1093/mnras/stad1941