Validations and Corrections of the SFD and Planck Reddening Maps Based on LAMOST and Gaia Data

Precise correction of dust reddening is fundamental to obtain the intrinsic parameters of celestial objects. The Schlegel et al. (SFD) and the Planck 2D extinction maps are widely used for the reddening correction. In this work, using accurate reddening determinations of about two million stars from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) data release 5 (DR5) and Gaia DR2, we check and calibrate the SFD and Planck maps in the middle and high Galactic latitudes. The maps show similar precision in reddening correction. We find small yet significant spatially dependent biases for the four maps, which are similar between the SFD and Planck2014-R maps, and between the Planck2014-Tau and Planck2019-Tau maps. The biases show a clear dependence on the dust temperature and extinction for the SFD and Planck2014-R maps. While those of the Planck2014-Tau and Planck2019-Tau maps have a weak dependence on the dust temperature, they both strongly depend on the dust spectral index. Finally, we present corrections of the SFD and Planck extinction maps within the LAMOST footprint, along with empirical relations for corrections outside the LAMOST footprint. Our results provide important clues for the further improvement of the Galactic all-sky extinction maps and lay an significant foundation for the accurate extinction correction in the era of precision astronomy.


INTRODUCTION
Interstellar dust is formed from the condensation of heavy elements which are expelled into the space by stellar winds or explosions. In the ultraviolet (UV), optical, and near-infrared (NIR) bands, dust absorbs and scatters the light passing through it. Therefore, correction for the effect of dust extinction is fundamental to reveal the intrinsic properties of observed objects. The existing extinction maps serve as a straight and convenient tool to correct for extinction.
The two-dimensional (2D) extinction map of Schlegel et al. (1998, hereafter SFD) is most widely used for dust corrections. The SFD map is based on the Infrared Astronomical Satellite (IRAS) 100 µm and the Dust Infrared Background Experiment (DIRBE) 100 yuanhb@bnu.edu.cn and 240 µm data. The far-infrared (FIR) emission is firstly modeled by a modified blackbody to obtain the dust temperature, which is then used to determine the dust column density and finally transferred to reddening by calibrating against 389 elliptical galaxies. Unfortunately, the map suffers from moderate systematic uncertainties as the dust temperatures have a low spatial resolution as well as the fact that the cosmic infrared background (CIB) and the zodiacal light have not been fully removed.
Many works have been carried out to validate and calibrate the SFD map. By comparing the surface number density and the average colors of galaxies from the Sloan Digital Sky Survey (SDSS; York et al. 2000) in the regions with different extinctions, Yahata et al. (2007) indicated that, toward low-reddening lines of sight, the results show the opposite of the effect of Galactic dust that the number density of galaxies increases with E(B − V ) regardless of extinction correction with the SFD map.
These unexpected effects were explained by the pollution of the sky brightness by extragalactic infrared emission of order 0.01 magnitudes in the SFD emission map. Using 151,637 passively evolving galaxies from the SDSS, Peek & Graves (2010) have presented corrections to the SFD map in the high Galactic latitude region. They found that the SFD map under-predicts extinction in regions of low dust temperature. The result is further confirmed by a large-scale reddening map from neutral hydrogen emission (Lenz et al. 2017). By measuring the reddening of a sample of over 50, 000 QSOs from the SDSS DR7, Wolf (2014) concluded that there is a non-linearity in the SFD map (see also Figure  10 of Li et al. 2017), which may be caused by pollution from unresolved extragalactic infrared emission. In addition, Schlafly et al. (2010), Schlafly & Finkbeiner (2011) and Yuan et al. (2013) measured the stellar reddening and found that the SFD extinction map overestimates E(B − V ) by about 14 percent.
Although many works have checked the errors of the SFD map, further investigations are still needed due to the limitations of previous works. The limitations include large statistical uncertainties caused by small sample sizes, systematic calibration errors in the photometric data used, and invalid assumptions adopted. Thanks to the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al. 2012) and Gaia (Gaia Collaboration et al. 2016), reddening values of an enormous number of stars can be estimated with high accuracy. The LAMOST surveys Deng et al. 2012;Liu et al. 2014) have obtained high-quality spectra and precise atmospheric parameters including effective temperature, surface gravity, and metallicity for millions of stars. Their reddening values E(G BP − G RP ) can be accurately measured from the star-pair method (Yuan et al. 2013) by combining with the Gaia photometric data.
In this work, we first fit empirical temperature-and . Then by comparing with the SFD dust map, we produce a correction map for it and investigate the potential error sources. Since Planck Collaboration et al. (2014) (hereafter Planck2014) applied a similar method as SFD to build a new extinction map (hereafter Planck2014-R) and provided another one based on 353 GHz absorption (hereafter Planck2014-Tau), we inspect and correct them at the same time. The map from Irfan et al. (2019) (hereafter Planck2019-Tau), which is an update of Planck2014-Tau, is also considered.
The paper is organized as follows. We describe our data and the method of calculating E(B − V ) in Section 2. The details of the four extinction maps and the way to build the correction factor k maps are introduced in Section 3. We present the final correction factor k maps and discuss potential error sources of the extinction maps in Section 4. We summarize in Section 5.

ACCURATE REDDENING ESTIMATION FROM
LAMOST AND GAIA DATA We use a star-pair method (Yuan et al. 2013) to measure reddening values of individual stars. The method assumes that stars of the same stellar atmospheric parameters should have the same intrinsic colors. The method requires a control sample of stars that are unreddened or of very low reddening. We can then obtain the intrinsic colors of a reddened star from its pairs out of the control sample having the same atmospheric parameters. In this work, we calculate values of E(G BP −G RP ) using the same star-pair algorithm of Yuan et al. (2015).

LAMOST and Gaia data
The stellar parameters used in this analysis are from the LAMOST Data Release 5 (DR5; Luo et al. 2015). The LAMOST spectra have a resolution of about R 1800 and cover a wavelength range of 3700 -9000Å. The LAMOST DR5 has provided stellar atmospheric parameters (effective temperatures T eff , surface gravities log g, and metallicities [Fe/H]) for over 5 million stars, with the LAMOST Stellar Parameter Pipeline (LASP; Wu et al. 2011). The typical uncertainties of the stellar parameters are about 110 K, 0.2 dex and 0.15 dex for T eff , log g and [Fe/H], respectively (Luo et al. 2015). Note the internal uncertainties are significantly lower (e.g., Niu et al. 2021).
To obtain accurate reddening values of LAMOST stars, we adopt the photometric data from the Gaia DR2 (Collaboration et al. 2018). The Gaia DR2 provides high precision G magnitudes of ∼ 1.7 billion sources and G BP and G RP photometric measurements of ∼ 1.4 billion targets. The Gaia G, G BP , and G RP filters cover wavelength ranges of 330 -1050, 330 -680, and 630 -1050 nm, respectively. The calibration uncertainties are 2, 5, and 3 mmag for G, G BP , and G RP , respectively (Evans et al. 2018).
To validate and correct 2D reddening maps, we only use stars at medium and high Galactic latitudes, where dust layers can be treated as thin screens and the reddening values of the observed stars are the same as those at infinite distances. Stars from the LAMOST DR5 are first cross-matched with the Gaia DR2. We then select stars with the following criteria: parallax > 0, Galactic latitude |b| > 20 • , distance to the Galactic plane |Z| > 200 pc, effective temperature T eff ∈ [4000, 8000] K, LAMOST spectral signal-to-noise ratio (SNR) at g-band SNR g > 20. This yields 2,083,741 sources as our sample.

Reddening determination
As we mentioned above, the star-pair method is used to estimate the reddening values for our sample stars. For a given sample star, its control stars are selected via the following criteria: ∆T eff < T eff * (0.00003 * T eff ) 2 K, ∆log g < 0.5 dex, ∆[Fe/H] < 0.3 dex, distance to the Galactic plane |Z| > 1 kpc and extinction value given by the SFD map E(B − V ) SF D < 0.02 mag. Then by using the same star-pair algorithm described in Yuan et al. (2015), we obtain E(G BP −G RP ) and E(G BP −G) for our sample. Stars of resulted reddening E(G BP − G RP ) < −0.05 mag are excluded. This leads to 2,067,805 stars in our sample, which we denote as 'LAMOST sample' in this work. Reddening coefficient R is usually used to convert reddening values of a given color to E(B − V ). The stellar extinction depends on the convolution of stellar spectral energy distribution (SED) and the extinction curve. Therefore, the reddening coefficient R is physically related to the T eff and the E(B − V ). This effect is particularly strong for the very broad Gaia passbands. Therefore, we use binary quadratic functions to fit R GBP−GRP and R GBP−G as functions of T eff and E(B − V ). A sample of stars is collected via the following criteria: • |b| > 15 • ; • |Z| > 300 pc; • If T eff > 4500 K, SNR g > 15. Otherwise, SNR g > 20; • E(B − V ) > 0.05 mag.
The selection criteria are different from those for selecting stars to correct reddening maps mentioned in Section 2.1. To make the fitting more reliable, here we exclude stars in very low reddening regions, and restrict SNR g more tightly. This yields about 700,000 sources for fitting. The exact number varies slightly for different extinction maps. In order to build grids, E(B − V ) values are equally divided into 8 bins from 0 to 0.8 mag, and T eff are also equally divided into 8 bins from 4000 to 8000 K. Then for each grid, the medians of E(B − V ), T eff , and R are estimated. The error of R is also calculated by the following formula: where N is the number of sources in the grid. The grids with N < 10 or R err > 0.03 are discarded in the fitting. We fit the relations for the four maps separately, considering their different systematics in E(B − V ). As to be mentioned in Section 3.1, the Planck2014-R map is unreliable in high extinction regions. Therefore, for this map we only choose grids whose E(B − V ) values are in the range of [0.05, 0.5] mag, and in the range of [0.05, 0.7] mag for other maps. The fitting coefficients are listed in Table 1. The results are plotted in Figure 1 It is because that the Planck2014-R map tends to underestimate reddening in high extinction regions.
Via the above relations, the E(B − V ) LAMOST values can be computed: The typical errors of E(B − V ) LAMOST for individual stars are around 0.01-0.02 mag.

The SFD and Planck maps
The SFD map traces dust reddening via thermal emission in the far-infrared based on the IRAS data. The IRAS 100 µm map was calibrated to match the COBE/DIRBE data. The zodiacal light and CIB contaminations were removed. The DIRBE 100 and 240 µm data were then used to estimate the dust temperature at a spatial resolution of 1.3 • and to transform the IRAS 100 µm flux to dust optical depth. By using a sample of 389 elliptical galaxies, the dust map was finally normalized to E(B − V ) reddening map. The SFD reddening map has a spatial resolution of 6.1 arcmin.
The Planck2014-R and Planck2014-Tau maps are based on the Planck 350, 550, and 850 µm data from the HFI 2013 delivery maps (Planck Collaboration et al. 2014) and the IRAS 100 µm data. Dust parameters, including the dust optical depth τ 353 , dust temperature T obs , and spectral index of the dust emission β obs , were obtained by the χ 2 fit of the observed dust SED with a modified blackbody (MBB) model. The dust radiance R was also calculated by integrating the MBB fit. Finally, Planck2014-R and Planck2014-Tau reddening maps were respectively obtained by transforming R and τ 353 to E(B − V ) using extinction measurements of SDSS quasars. Both maps have a resolution of 5 arcmin. Planck2014 recommended that their Planck2014-R reddening map should only be used to estimate reddening in lines of sight where E(B − V ) < 0.3 mag. However, to make a more robust comparison, we slightly enlarged this limit by adopting sightlines of reddening E(B − V ) <0.5 mag. Finally, Planck2014 indicated an offset of −0.003 mag for the SFD map. We thus corrected the SFD maps by adding the offset.
Based on the Planck Release 2 353, 545, and 857 GHz maps as well as the IRAS 100 µm data, Irfan et al. (2019) determined the MBB model parameters of thermal dust emission by a new sparsity-based, parametric method and presented the Planck2019-Tau map. By adopting the new method, they are able to produce full-resolution MBB parameter maps without smoothing and account for the CIB without removing thermal dust emission through over-smoothing. We use the same coefficient for the Planck2014-Tau map to convert τ 353 to E(B − V ). The Planck2019-Tau map has a resolution of 5 arcmin.

Correction Factor
To validate and correct the SFD and Planck maps, we first obtain the E(B − V ) map values for the individual stars in the LAMOST sample from the SFD and Planck maps, and then compare these reddening values to those derived from the star-pair method E(B − V ) LAMOST . Sample stars are then divided into different subfields (pixels) by the HEALPix grid (Gorski et al. 2005) at Nside = 64 (with a resolution of about 1 • ). There are typically 100 stars in one pixel. For each pixel, we as- where k is the correction factor and is estimated by averaging the ratios of E(B − V ) LAMOST and E(B − V ) map after 3σ clipping. Its uncertainty is derived by, where E(B − V ) fit = k ×E(B − V ) map and N is the total number of sources for a given pixel. Figure 3 shows the comparisons of E(B − V ) LAMOST and E(B − V ) map in six selected sight-lines. For a given sightline, the standard deviations between different maps are similar. The typical standard deviations are between 0.01 -0.03 mag, increasing towards high extinction sight-lines.

RESULT AND DISCUSSION
4.1. Removing bad pixels Figure 4 shows the spatial variations of k err for the SFD and Planck reddening maps. Excluding pixels in the very low extinction regions of the maps, most k err values are less than 0.1, only a few pixels have k err > 0.1 due to their small numbers of stars or large deviations between E(B − V ) LAMOST and E(B − V ) f it . Figure 5 plots the spatial and histogram distributions of star numbers in the individual pixels for the SFD map as an example. There are about 100 stars for a typical pixel. Only 887 pixels have less than 10 stars. The results for the three Planck maps are similar. We thus exclude these pixels in the following analyses. Figure 6 shows Std fit against E(B − V ) map for the SFD and Planck reddening maps, where Std fit is the standard deviation of the reddening fit residuals Values of k err are also color-coded. Only pixels having over 10 stars are plotted. Std fit increases with E(B − V ) map . Secondorder polynomials are applied to fit the binned median values. The pixels that have Std fit values larger than the fitted curves by 2σ deviation are excluded. There are about 400 pixels in each panel of Figure 6 that have very large values of Std fit and consequently k err > 0.1. Note not all pixels eliminated from the four extinction maps are the same. A total of 739 pixels are discarded from the four extinction maps.
We have checked the reddening-distance profiles of these pixels. About half (374) of them have dust clouds at large distances to the Galactic plane Z (Figure 7). For these pixels, there are lots of stars locating in front of the dust cloud. Thus the reddening values of these stars E(B − V ) LAMOST are much smaller than those from the reddening maps E(B − V ) map , causing large values of Std fit . Figure 8 shows the spatial distribution of all the discarded pixels for the four dust maps. The number of rejected times are also marked. It must be admitted that not all pixels with associated dust cloud are well eliminated in the four extinction maps. However, due to the small number of such pixels, they will not affect the statistical results obtained in the following analyses. The

Precisions of the SFD and Planck maps
Before we discuss the possible factors that affect the spatial variations of k, we first check and validate reddening correction precisions of the SFD and Planck maps. As there are systematics of E(B − V ) values between the four reddening maps, we compare E(G BP − G RP ) values here. For each pixel of each map, we estimate the dispersion of the differences between E(G BP − G RP ) LAMOST and E(G BP − G RP ) map , where E(G BP − G RP ) map is calculated using the corresponding reddening coefficient. For each map, the median dispersion values are plotted against E(G BP −G RP ) LAMOST in Figure 12. Dispersion values increase with the reddening values of the pixels. The curves of the four maps are very close, suggesting that the four reddening maps can achieve similar precisions in reddening correction. In the low-extinction region (E(G BP − G RP ) < 0.1 mag), the Planck2014-Tau map is slightly worse. In the relatively high-extinction region, the Planck2019-Tau map is slightly better. The typical dispersion is 0.016, 0.022, and 0.045 mag at E(G BP −G RP ) LAMOST = 0.01, 0.1, and 0.7 mag, respectively. Note the dispersion values include measurement uncertainties of E(G BP − G RP ) LAMOST , therefore are upper limits of the reddening correction precisions.

Properties correlated with the correction factor
In Figure 10, the spatial distributions of the correction factor k show clear patterns. We explore here the possible origins of the patterns.

Reddening values
The panels in the first row of Figure 13 Table  2.
The overall relations between k and E(B − V ) map are similar for the four reddening maps (particularly between the SFD and Planck2014-R maps, and between the Planck2014-Tau and Planck2019-Tau maps). All show dips at E(B − V ) map ∼ 0.03 mag when E(B −V ) map < 0.1 mag. The trends become flat around 1 when E(B −V ) map > 0.1 mag, as expected considering that reddening-dependent coefficients are used (see Section 2.2) to compute the E(B − V ) LAMOST values and then the k values. The dips in low reddening regions are likely caused by the stronger over-estimation of reddening compared to that in the high reddening regions.
We have explored the relations between k and E(B − V ) map for five different sky areas in the remaining rows of Figure 13. For a given sky area, the relations are similar between the SFD and Planck2014-R maps, and between the Planck2014-Tau and Planck2019-Tau maps. However, for a given map, the relations are different between different sky areas, suggesting that reddening is not the main factor for the spatial variations of k.

Dust Temperature and spectral index
Since the SFD and Planck maps are all based on the MBB model of dust thermal emissions, the correction factor k could be correlated to the model parameters, i.e. dust temperature T d and spectral index β. Note that T d and β parameters are usually fitted simultaneously and suffer strong anti-correlation in the MBB model (e.g., see Figure 7 of Planck2014). The SFD map adopted a constant value of β = 2. While the Planck2014-R map was based on dust radiance (R), which is the integrated intensity (See Equation 9 in Planck2014) and thus does not suffer from any degeneracy in the fit parameters. The possible effect of spectral index β on the Planck2014-R map should be very weak. Therefore, we first investigate the effect of T d on k for the SFD and Planck2014-R maps.
The upper-left and upper-right panels of Figure 14 plot k as functions of E(B − V ) and T d for the SFD and Planck2014-R maps, respectively. A strong dependence on T d is found for the SFD map, particularly in the low extinction regions (E(B − V ) < 0.3 mag), where k decreases when T d increases. It suggests that the SFD map underestimates/overestimates reddening in low/high dust temperature regions. One possible explanation is that the SFD T d map "underestimates/overestimates" T d in high/low dust temperature regions. The result is consistent with Peek & Graves (2010) and Lenz et al. (2017).
For the Planck2014-R map, the dependence on E(B − V ) and T d is more clearly visible, in both the low and high extinction regions. It suggests that the Planck2014-R map underestimates/overestimates reddening in low/high dust temperature regions, and indicates that the Planck dust radiance map underestimates/overestimates the integrated intensity in low/high dust temperature regions.
We also investigate the dependence of k on T d for the Planck2014-Tau and Planck2019-Tau maps. Very weak relations are found. However, a strong dependence on β is found for the above two maps. To show the dependence more clearly, we first perform a normalization of the k values with respect to the red lines in the top panels of Figure 13 to get rid of the effect of reddening. The normalized k values are named k hereafter. The upper-left panel of Figure 15 shows k as functions of T d and β for the Planck2014-Tau map. We also divide the pixels into four different reddening ranges, and their results are plotted in other upper panels. It can be seen that for a given T d , k increases as β increases. Similar results are found for the Planck2019-Tau map, as shown in Figure 16.

Fitting of the correction factor
To provide corrections of the SFD and Planck maps for regions outside the footprints in the current work, we fit the correction factor k as functions of E(B − V ) and T d for the SFD and Planck2014-R maps, and the normalized correction factor k as functions of T d and β for the Planck2014-Tau and Planck2019-Tau maps.
Two-dimensional fourth-order polynomials of 15 free parameters are adopted in the fitting with the Python procedure SciPy. The resultant coefficients are listed in Table 3. The fitting residuals are shown in Figure 14 for the SFD and Planck2014-R maps, Figure 15 for the Planck2014-Tau map, and Figure 16 for the Planck2019-Tau map, respectively. One can see that global trends with E(B − V ), T d or β are successfully removed. Figure 17 shows spatial distributions of k/k , k fit /k fit , and ∆k/∆k for the SFD and Planck maps. The fluctuations are smaller in the residual maps. However, the fitting residuals still show similar spatial patterns, despite the fact that they show no dependence on E(B − V ), T d or β. Figure 18 compares fitting residuals between different maps. The correlation coefficients are between 0.63 -0.90. The result indicates that there are likely other factors (e.g., dust sizes/reddening laws) that contribute the spatial variations of the correction factors. Such possibilities will be explored in the future work.

SUMMARY
Combining high-precision LAMOST DR5 spectroscopic data and Gaia DR2 photometric data, we have calculated the color excess E(G BP − G RP ) for about two million well selected middle and high Galactic latitude stars by using the star-pair technique. With empirical temperature-and reddening-dependent coefficients, values of E(B − V ) LAMOST are further obtained to check and correct the SFD, Planck2014-R, Planck2014-Tau, and Planck2019-Tau reddening maps. By comparing E(B − V ) LAMOST with E(B − V ) map , the following results are found: 1. On one hand, the overall agreements between the E(B − V ) LAMOST and E(B − V ) map maps are excellent (Figure 9). On the other hand, the small yet significant discrepancies show clear spatiallydependent patterns (Figure 10). The patterns for the SFD and Planck2014-R maps are similar, while those for the Planck2014-Tau and Planck2019-Tau maps are similar.
2. The four reddening maps can achieve similar precisions in reddening correction (Figure 12). In the low-extinction region (E(G BP − G RP ) < 0.1 mag), the Planck2014-Tau map is slightly worse. In the relatively high-extinction region, the Planck2019-Tau map is slightly better.
3. For a given sky area, the k and E(B − V ) map relations are similar between the SFD and Planck2014-R maps, and between the Planck2014-Tau and Planck2019-Tau maps. However, for a given map, the relations are different among different sky areas (Figure 13), suggesting that reddening is not the main factor for the spatial variations of k.
4. For the SFD and Planck2014-R maps, the dependence of k on E(B − V ) and T d is clearly visible (Figure 14). k decreases when T d increases. It suggests that the two maps underestimates/overestimates reddening in low/high dust temperature regions. One possible explanation is that the SFD T d map "underestimates/overestimates" T d in high/low dust temperature regions, consistent with Peek & Graves (2010) and Lenz et al. (2017).
While the Planck dust radiance map possibly underestimates/overestimates the integrated intensity in low/high dust temperature regions.
5. For the Planck2014-Tau and Planck2019-Tau maps, the dependence of k on T d is very weak, but very strong on β. For a given T d , the normalized correction factor k increases as β increases ( Figure 15 and Figure 16).
The k maps and their errors are publicly available 1 and can be used to perform corrections of the SFD and Planck maps. For regions outside the footprints in the current work, relations of k as functions of E(B − V ) and T d for the SFD and Planck2014-R maps, and k as functions of T d and β for the Planck2014-Tau and Planck2019-Tau maps, can be used to correct global trends with E(B − V ), T d or β to some extent (Table 3 and Figure 17). It should be noticed that, applications of the empirical correction relations are limited by the range of E(B − V ) for fitting. For the convenience of use, a python routine is provided 2 for such purposes. However, the fitting residuals between different maps still show similar spatial patterns and good correlations (Figure 18), indicating that there are likely other factors (e.g., dust sizes/reddening laws) that contribute the spatial variations of the correction factors. Such possibilities will be explored in the future work.
Our results provide important clues for the further improvement of the Galactic all-sky extinction maps and lay an important foundation for the accurate extinction correction in the era of precision astronomy.

ACKNOWLEDGMENTS
We acknowledge the referee for his/her valuable comments to improve the clarity and quality of the manuscript. This work is supported by the National Key Basic R&D Program of China via 2019YFA0405503, the National Natural Science Foundation of China through the projects NSFC 12173007, 12173034 and 11603002, and Beijing Normal University grant No. 310232102. We acknowledge the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A08 and CMS-CSST-2021-A09. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/ gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.    Table 3. Fitting coefficients of k and k . The function form is: z = ax 4 + by 4 +cyx 3 +dxy 3 +ex 2 y 2 +f x 3 +gy 3 +hyx 2 +ixy 2 +jx 2 +ky 2 +lxy+mx+ny+o. For the SFD and Planck2014-R maps, z is k, x is T d , and y is E(B − V ). For the Planck2014-Tau and Planck2019-Tau maps, z is k , x is T d , and y is β.       Figure 6). Dust clouds far above the Galactic disk (Z ∼500 pc) are found as jumps of k in the panels.