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

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.


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,b;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 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 < λ/∆λ < 100 in BP and 70 < λ/∆λ < 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 tra-An et al. ditional 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 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 eff , [Fe/H], dis-tance, and extinction) for the majority of stars with XP spectra, achieving uncertainties of ∼ 0.15 dex in [Fe/H] and ∼ 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 totalto-selective extinction ratio (R V ≡ 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.

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 functions1 .To calibrate the raw (internally calibrated, continuously represented) mean spectra, we employed a Python library suite, GaiaXPy (De Angeli et al. 2023) 2 .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 BP − G 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 ∼ 80 million.

Strategy
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.

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 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 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 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).
We adopted the linear relationship between [α/Fe] and [Fe/H] owing to the predominant influence of T 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 [α/Fe] is minimal.From an observational standpoint, although modern spectroscopic surveys offer valuable [α/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 [α/Fe], making it challenging to precisely determine the dependence of stellar flux on [α/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 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 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 ≤ [Fe/H] ≤ +0.4 and 4000 K ≤ T eff ≤ 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.

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 χ 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 ≡ 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 χ 2 distribution within a model grid, with grid intervals set at ∆[Fe/H] = 0.05 dex and ∆M * = 0.02 M ⊙ .The process involved fitting a paraboloid to the χ 2 distribution to locate the minimum χ 2 , thereby determining the bestfitting set of parameters for a star.To gauge parameter uncertainties, we considered ∆χ 2 = 3.53 and 4.72 from the minimum χ 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 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 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 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 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 χ 2 min /ν < 5, σ(E(B − V )) < 0.3 mag and σ([Fe/H]) < 1.5 dex, where ν represents the number of degrees of freedom.We retained stars with metallicities falling within our model grid (−3 < [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 × 10 7 , 2.3 × 10 7 , 2.8 × 10 7 , and 2.5 × 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 BP − G 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 ∆E(B − V ) ∼ 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.

Systematic Deviations and Uncertainties in 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 < λ < 1000 nm, with an rms difference of only 1.4%.The mild bump at ∼ 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 li-An et al. brary' (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 ∼ 10% at 350 nm-a value that escalates with higher G BP − G 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 ex-tracted 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.
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 (λ < 400 nm) has minimal impact on our parameter search, based on comparisons to spectroscopic parameter estimates from GALAH and APOGEE.

Multiresolution HEALPix Scheme
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 cov-

An et al.
erage on the celestial sphere.The number of cells and their respective areas are determined by a single parameter, denoted as N side .In this illustration, we selected N side = 128, resulting in a pixel area equivalent to 0.2 deg 2 .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.
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 • × 30 • patch centered around (l, b) ∼ (25 • , 15 • ), 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 side = 128, but we adjusted the level of pixelization to ultimately arrive at HEALPix cells as large as having N 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 side and the number of cells correspond to N side =16 (114 cells), 32 (3872), 64 (24, 802), and 128 (28, 152).The surface area ranges from 0.2 to 13.4 deg 2 .As illustrated in panel (c), each multiresolution cell contains between ∼ 200 and ∼ 1000 stars, with a median count of ∼ 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.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 ≈ 7 and 9.5, respectively), representing nearby dust clouds and those possibly associated with large-scale spiral arm(s), respectively.

Modeling the Cumulative Reddening
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 χ 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 sequencebased 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 sequencebased solutions (Case B and Case D) exhibit values approximately ∆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 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: where µ 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 µ i jointly determine the shape of the cumulative E(B − V ) distribution, whereas A i serves as a normalization factor.
In equation ( 1), ϵ 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 • ) 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 ϵ 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 , µ 1 , k 2 , µ 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 • -3 • ).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.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.
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σ 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 els display data taken near the Galactic plane, while the right panels show observations at |b| ≈ 30 • , 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 ) ≲ 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 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 ≲ 1 kpc, reflecting the limited number of stars within the magnitude-limited XP MS sample.Beyond |b| > 10 • , our sample encompasses volumes extending, on average, up to ∼ 3 kpc, with certain areas exhibiting even greater depth.

3D Reddening Map
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.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 1).
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| ≲ 10 • , the E(B − V ) difference between Case A and Case B tends to increase with distance.At a distance of 1584 pc, approxi-mately 1% of the area exhibits a deviation greater than ∆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 side = 128, effective 2σ uncertainties, calculated from the 15.9th and 84.1st 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 ∼ 0.005 mag at |b| ≳ 10 • .
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 ∆(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 ∼ 250 pc owing to the limited number of nearby stars in the XP MS sample.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 reso-lution (∼ 6 ′ ) in comparison to our baseline map with N side = 128 (having a mean spacing of 30 ′ ), 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 ∆b = 2.5 • ; it exceeds ∆E(B − V ) = 0.1 mag for |b| < 10 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 (∼ 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.006mag 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.

Green et al. (2019)
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 con-structed a high-resolution extinction map by subdividing HEALPix cells down to a level of N 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 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 • 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 • .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 ∆E(B − V ) = 0.03 mag, while it stays below 0.015 mag at 630 pc.At a distance of 1584 pc,  1) in the Galactic coordinate system as a function of heliocentric distance.
a more pronounced deviation is evident near the Galactic plane, but the median difference remains modest at |b| > 10 • , rarely exceeding ∆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.

Zhang et al. (2023b)
Figure 14 provides a comparison to more recent, allsky extinction measurements presented in Zhang et al. (2023b).Similar to ours, they utilized XP spectra to derive line-of-sight extinction.Following their rec-ommendations, we applied a set of quality flags, including quality flags < 8, teff confidence > 0.5, feh confidence > 0.5, σ(T eff ) < 500 K, σ([Fe/H]) < 0.5 dex, and a maximum 20% uncertainty in distance.In each of the HEALPix cells with N 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  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 • ), the difference in reddening amounts to up to ∆E(B − V ) ≈ 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)): ∆E(B − V ) = 0.019 ± 0.006 mag and 0.044 ± 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.
Star-by-star comparisons of reddening and [Fe/H] are illustrated in Figure 15.To facilitate this comparison, we restricted our analysis to stars within a 3 • radius centered at (l, b) = (20 • , 20 • ), utilizing the same quality flags for data selection.The linear fit to the data indicates that their E estimates are 13% larger than our values (orange line in the top left panel), which is still significantly larger than expected from the different definitions of reddening parameters.While there is a significant correlation between reddening and [Fe/H] in the parameter estimation (bottom left panel), the discrepancy cannot be attributed to a systematic deviation in their metallicity scale, as the two independent methods yield essentially consistent metallicity estimates (top right panel).This is further supported by the comparison in the bottom right panel, showing [Fe/H] along a line of sight at (300 • , 45 • ) with small foreground extinction (E(B − V ) < 0.1).In Appendix A, we further validate our parameter estimation by comparing it with metallicity estimates from GALAH and APOGEE.

SYSTEMATIC CHANGES IN THE EXTINCTION LAWS
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 (∆[Fe/H] ≈ 0.05 dex), there are small-scale systematic variations in [Fe/H].For instance, the region near l ∼ 130 • and b ∼ 5 • exhibits relatively higher [Fe/H] compared to the surrounding areas, while an opposite trend is observed near l ≈ 355 • and b ≈ 15 • .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 (Fitz-  patrick 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 ∼ 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) ≈ (355 • , 15 • ), Taurus (170 • , −15 • ), Orion (210 • , −20 • ), and Chamaeleon (300 • , −15 • ).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 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 (⟨E(B − V )⟩ ∼ 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 ∼ 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 λ /E(B − V ) at R V = 3.1 in both cases and then scaling the A λ /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 ≤ R V ≤ 6. 4 We found ϵ 0 = {0.015,0.029} mag for Case A and Case B, respectively, for the zero-point in equation (1).
Panel (a) of Figure 17 shows the distribution of mean R V values estimated using stars within (m − M ) 0 ≤ 12 from both Case A and Case B.Here we included only regions with significant reddening (E(B − V ) ≥ 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 ∆R V ≲ 0.1 within |b| ≲ 30 • , 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.9th and 84.1st percentiles of R V from Case A. In the vicinity of the Galactic plane, this effective 1σ deviation, σ(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 ϵ ≡ σ(R V )/ √ 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 ≈ 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 ϵ as an effective 1σ 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 ∼ 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 (∼ 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.
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 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 ∼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.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.

SUMMARY AND CAVEATS
In this study, we utilized the extensive dataset of lowresolution 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).In addition, our metallicity estimates show good agreement with high-resolution spectroscopic results from the GALAH survey and, to some extent, SDSS/APOGEE (see Appendix A).
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.Both our method and that of Zhang et al. yield consistent metallicity estimates, suggesting that the reddening differences arise from systematic factors other than the metallicity scale.
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.Con-sequently, 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 finestructure 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.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 Ga-iaXPy, 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).

A. COMPARISONS TO HIGH-RESOLUTION SPECTROSCOPIC PARAMETER ESTIMATES
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 ≈ 28, 000 and R ≈ 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 ′′ .We specifically included our metallicity measurements with the following criteria: σ(E(B − V )) < 0.3, 3.5 < M r < 7.5, and σ([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 ∆[Fe/H] ≲ 0.1 at [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 λ ≳ 400 nm can offer a wealth of information about photospheric metal abundance and foreground extinction.

B. COMPARISONS TO CLUSTER E(B − V ) ESTIMATES IN THE LITERATURE
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. (2007An et al. ( , 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 BP and G 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 • .
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 • .(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.
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σ 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.02mag.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 ).

C. ACCESSING EXTINCTION DATA CUBES
The primary data products resulting from this study can be accessed through the following link. 5Within 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.

Figure 1 .
Figure 1.Model fits illustrating comparisons for various effective temperatures (T eff , left) and metallicities ([Fe/H]; right).Best-fitting estimates for T 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 2 .
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.

ab
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 deg 2 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 deg 2 (28, 152 cells), 0.8 deg 2 (24, 802), 3.4 deg 2 (3872), and 13.4 deg 2 (114).(c) Distribution of the number of stars per cell, derived from the multiresolution grid scheme shown in panel (b).

Figure 5
Figure 5 illustrates the line-of-sight distribution of reddening toward l = 149.8• and b = 4.8 • , displaying

Figure 4 .
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 deg 2 (N 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) ∼ (20 • , 5 • )) 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.andCase 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 ≈ 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:

Figure 5 .
Figure 5. Reddening of stars against distance modulus in a single, multiresolution HEALPix cell located along l = 149.8• and b = 4.8 • .The four panels display different sets of E(B − V ) measurements, Case A through Case D (see Table1).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 7 .
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.
However, finer structures become more evident when these maps are divided into narrower distance bins, as shown in Figure9, on the same multiresolution grid as depicted in panel (b) of Figure3.It presents cross sections of our 3D extinction map from Case A within selected distance intervals.The ∆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 Figure9; see alsoSchlafly 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 (⟨d⟩ ∼ 100 pc) becomes visible in the nearest distance bin (panel (a)), situated toward the Galactic center at (l, b) ∼ (0 • , +20 • ).This complex overlaps with the Aquila Rift (⟨d⟩ ∼ 100 pc), centered at (l, b) ∼ (30 • , +5 • ), which gains prominence in the larger distance bin (panel (b)).At (l, b) ∼ (300 • , −15 • ), the Chamaeleon cloud complex is visible (Voirin et al. 2018, ⟨d⟩ ∼ 200 pc).In the direction of the Galactic anticenter, the Taurus-Perseus-Auriga cloud complex (⟨d⟩ ∼ 100-300 pc) is also discernible within these volumes, located at (l, b) ∼ (170 • , −10 • ).In panels (b) and (c), the Orion Nebula (⟨d⟩ ∼ 400 pc) centered approximately at (210 • , −20 • ) and the Cepheus and Polaris Flares (⟨d⟩ ∼ 400 pc) at (l, b) ∼ (110 • , +20 • ) come into play.Beyond this point, a predominantly flattened structure of dust prevails along the Galactic plane at greater distances.
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).

Figure 8 .
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 (Table1).

Figure 9 .
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. 5. COMPARISONS WITH PREVIOUS EXTINCTION MAPS 5.1.Schlegel et al. (1998)

Figure 10 .
Figure 10.Difference in modeled E(B − V ) between different cases of parameter estimates (Table 1) in selected distance bins.ure 12 demonstrates a rapid convergence of both maps for |b| > 10 • .Specifically, at |b| = 30 • , 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 (∼ 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.006mag 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 Figure12.

Figure 11 .
Figure11.Average E(B − V ) from Case A and Case B (Table1) in the Galactic coordinate system as a function of heliocentric distance.

Figure 12 .
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.

Figure 13 .Figure 14 .
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 • .Error bars denote their standard deviations.

Figure 15 .
Figure 15.Comparison of reddening and metallicity estimates between Zhang et al. (2023b) and Case A for the identical XP sample from a region with a 3 • radius, centered at (l, b) = (20 • , 20 • ).An additional [Fe/H] comparison is presented in the bottom right panel along a line-of-sight with low extinction at (l, b) = (300 • , 45 • ).Unity lines are represented by dashed lines, while null differences are displayed in the bottom left panel.
Figure 16.Mean metallicity estimates on the sky in the Galactic coordinate system assuming RV = 3.1, for (a) Case A and (b) Case B.

Figure 18 .
Figure 17.Mean and statistical properties of RV estimates from stars within (m − M )0 ≤ 12.Only regions with significant reddening (E(B − V ) ≥ 0.1 at (m − M )0 = 12) are included.(a) Mean RV estimates from Case A and Case B. (b) Differences in RV between Case A and Case B. (c) Half of the difference between the 15.9th and 84.1st percentiles of RV from Case A (σ(RV )).(d) Significance of the deviation of RV from the fiducial case (RV = 3.1) in relation to the error in the mean RV (ϵ ≡ σ(RV )/ √ N ).

Figure 19 .
Figure 19.Differences in AV derived from Case A when assuming a fixed value of RV = 3.1, compared to using a median RV (see text).The comparisons are illustrated at three representative distances.
Figure20.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.

Figure 21 .
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 inHarris (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.

Table 1 .
A Family of Solutions on Stellar Parameters

Table 2 .
Mean Model Deviations 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.