Extinction of Taurus, Orion, Perseus and California Molecular Clouds Based on the LAMOST, 2MASS and Gaia surveys I: Three-dimensional Extinction and Structure

The optical and near-infrared (NIR) bands are selected for measuring the clouds’ extinction. Because the visual extinction is much larger than the near-infrared (the V band extinction is about ten times that in the K band (Wang & Chen, 2019)), the two wavelength ranges are expected to trace both the diffuse and dense clouds. Besides, the ratio of the optical-to-NIR extinction is an indicator of the extinction law to be studied in our next work.

The optical photometric data are taken from the space telescope Gaia for its high precision and full sky coverage. The adopted Gaia EDR3 data contains 1.5 billion sources whose photometry is accurate to mmag-level in the $G_{\rm BP}$, $G$, $G_{\rm RP}$ bands centering at 532, 673, and 797 nm respectively (Jordi et al., 2010). The very wide band, $G$, is not used in this work because the extinction coefficient of this band varies strongly with stellar effective temperature and the extinction itself (Danielski et al., 2018). The NIR data are taken from the PSC of 2MASS that surveyed the whole sky in the $J$, $H$, $K_{\rm s}$ band, which bring about the photometry over more than 500 million sources with the average uncertainty better than 0.03mag (Skrutskie et al., 2006).

The stellar parameters are taken from the LAMOST DR8 (Luo et al., 2015). The LAMOST DR8 provides the effective temperature ($T_{\rm{eff}}$), surface gravity (log $g$) and metallicity ($\rm\left[Fe/H\right]$) for more than six-million stars. For the duplicated sources, only the parameters derived from the spectrum with the highest signal-to-noise ratio are kept. Then about 6.4 million stars are kept and cross-matched with the Gaia EDR3 and 2MASS PSC datasets respectively by a radius of $1\arcsec$. The distance is a key parameter to investigate the 3D extinction of the clouds. Here, we take the geometric distance provided by Bailer-Jones et al. (2021) which contains the distances and their uncertainties for 1.47 billion stars.

The data quality is further restricted as following to exclude the poor measurements:

1. 1.

   $\rm\left[Fe/H\right]$ is within [-1.0, 0.5].

2. 2.

   The error of $T_{\rm{eff}}$, $\rm\left[Fe/H\right]$ and ${\rm log}\ g$ is smaller than 150K, 0.15 dex and 0.3 dex respectively.

3. 3.

   The signal-to-noise ratio in the $g$ band (the parameter “snrg” in the LAMOST DR8 catalog) is larger than 10.

4. 4.

   The photometric error is $<0.05$mag in the Gaia EDR3 bands and $<0.1$mag in the 2MASS bands.

With these limitations, $\sim$4.6 million are left for calculating the intrinsic color index $C_{\rm G_{BP},G_{RP}}^{0}$ while $\sim$4.3 million stars are kept for $C_{\rm J,K_{s}}^{0}$. Furthermore, only dwarf stars are retained for three reasons: (1) metallicity has less effect on the intrinsic colors of dwarfs than of giants (see Figure 1 in Zhao et al. 2020); (2) giants in the LAMOST catalog are mostly at a distance $>$ 1kpc, much further than the studied MCs; and (3) stellar parameters of dwarfs are more accurately determined from the LAMOST observations. In practice, the dwarf stars are selected from the Kiel diagram in Figure 1 by ${\rm log}\ g>3.8$ for $T_{\rm{eff}}>6600K$ or $-0.0644\times T_{\rm{eff}}^{2}+0.457\times T_{\rm{eff}}+3.09<{\rm log}\ g$ for $T_{\rm{eff}}<6600K$. This reduces the sample to $\sim$ 3.2 million and 3.1 million stars for $G_{\rm BP}-G_{\rm RP}$ and $J-K_{\rm s}$ respectively.

Lombardi et al. (2010) defined the boundary of TMC in $(l,b)=([165\arcdeg,180\arcdeg],[-10\arcdeg,-20\arcdeg]$). Later, Bialy et al. (2021) recognized an ellipse substructure in the 3D space called the Tau Ring centering at $(l,b,d)=(179\arcdeg.5,\ -14\arcdeg.2,\ 179{\rm pc})$ with a semimajor and a semiminor axis of 39pc and 26pc respectively, and a projection radius of about $10\arcdeg$ at $d=218$pc. Including all these structures, we extend the boundaries a little to cover a larger area where the extinction is comparatively low, and select the boundary as $(l,b)=([160\arcdeg,195\arcdeg],[10\arcdeg,-40\arcdeg]$.

The area selection for the other three MCs follows the same rule as for the TMC, i.e. we select an area that is extended a little from the previously defined. Lombardi et al. (2010) delimited the PMC to $155^{\circ}<l<165^{\circ}$ and $-25^{\circ}<b<-15^{\circ}$, and Rezaei Kh. & Kainulainen (2022) defined $155^{\circ}<l<170^{\circ}$ and $-14^{\circ}<b<-6^{\circ}$ for CMC which agrees with Lombardi et al. (2010). For OMC that is composed of three major structures, Lombardi et al. (2011) suggested $203^{\circ}<l<217^{\circ}$ and $-21^{\circ}<b<-17^{\circ}$ for Orion A, $201^{\circ}<l<210^{\circ}$ and $-17^{\circ}<b<-5^{\circ}$ for Orion B and $188^{\circ}<l<201^{\circ}$ and $-18^{\circ}<b<-7^{\circ}$ for $\lambda$ Orionis. Taking all these into consideration, our selected regions of the four MCs are listed in the first and second column of Table 1. There are some overlaps in sightlines between some clouds such as TMC and CMC. Section 4 will introduce the method to further separate them in the 3D space.

By summing up the four clouds region in Table 1, the whole region to be studied is within $(l,b)=([130\arcdeg,220\arcdeg],[-50\arcdeg,+20\arcdeg])$, a total area of about 5000 deg². The area of each cloud and the density of tracing stars are displayed in Figure 2.

The principle in deriving the 3D extinction in the area of the TMC, OMC, PMC and CMC is that the extinction along a sightline is a non-decreasing function of distance. With the distance by Bailer-Jones et al. (2021) and the color excess calculated above, a compromise between the number of tracers and the spatial resolution yields a step of 0.2° in both longitude and latitude. Due to the non-uniform extinction within a selected area of 0.2° squared and the error in distance and color excess, the variation of color excess with distance is scattering and sometimes not monotonically increasing as shown by the dots in Figure 4.

The isotonic regression in the SCIKIT-LEARN package for PYTHON (Pedregosa et al., 2011) is applied to trace the general tendency. For a given set of observations ($x_{1}$, $y_{1}$), ($x_{2}$, $y_{2}$),…,($x_{n}$, $y_{n}$), isotonic regression is a non-parametric regression that seeks a weighted least-square fit $\hat{y_{i}}$ for all $i$ in a monotonic model. It solves the following problem:

where $y_{i}$ is the calculated color excess of star $i$, $\hat{y_{i}}$ is the fitted color excess along the given line of sight at the distance of star $i$, and $w_{i}$ is the weight. Typically, the weights are equal to 1 for all $i$. However, because the stars closer to the center of the given sightline are more representative, an exponential kernel function is adopted to describe $w_{i}$ as following:

where $\theta$ is the angular distance of star $i$ to center of the selected sightline, $\gamma$ is a scale parameter, and $\theta_{0}$ defines the radius of the selected bin size. Meanwhile, the number of stars that have non-zero weights should be more than 10 in one sightline to ensure the credibility of the model fitting. Considering the existence of some sightlines where there are less than 10 objects within a 0.2° circle, a larger $\theta_{0}$ is selected. The scale parameter $\gamma$ is set to be equal to the resolution. Figure 4 shows the fitting results with three test values of $\gamma$ in three sightlines as examples. After testing the model with several values of resolution (i.e. $\gamma$), $0.2^{\circ}$ and $1.0^{\circ}$ (5 times of $\gamma$) are selected respectively for $\gamma$ and $\theta_{0}$.

The extinction is sliced every 25pc from 100pc to 600pc by the above non-parameter fitting, and the results are presented in Figure 5. In combination with the continuity in the sky area, it can be seen that the four MCs appear in order of distance, i.e. TMC, PMC, OMC and CMC. The distance extents of each molecular cloud listed in Table 1 are judged by observing the main extinction structure (high extinction areas) within its boundary in Figure 5. It should be noted that the extension is partly caused by the substructures that appear at different distance instead of the thickness of the cloud. The details for each substructure will be discussed later.

The results are compared with Green et al. (2019) and Leike et al. (2020) in the same area and distance range. For comparing with Green et al. (2019), the integrated extinction map up to 600pc from this work and Green et al. (2019) is displayed in Figure 6. The value of $E_{\rm{B,V}}^{\rm{Green}}$ in Green et al. (2019) is converted to $E_{\rm G_{BP},G_{RP}}$ with the factor $E_{\rm{B,V}}^{\rm{Green}}/E_{\rm{G_{BP},G_{RP}}}=0.71\pm 0.01$ suggested by Sun et al. (2021a). It can be seen that the two results are roughly identical with a mean difference and its standard deviation of $0.007\pm 0.093$mag respectively. Nevertheless, it is obvious that the extinction value is comparatively smaller in the high-extinction regions (the blue part in Figure 6), such as the dense MC regions marked by the black lines in the right panel in Figure 6. This can be explained by the relatively shallow depth of the LAMOST survey that has a limiting magnitude of about $15$ to $17$mag in the $g$ band, which is unable to detect the stars with a large extinction. The largest extinction in this work is about $2.0$mag in $E_{\rm G_{BP},G_{RP}}$ and $0.8$mag in $E_{\rm J,K_{s}}$, equivalent to $A_{V}\sim 5$mag, occurring in some Galactic plane areas. This can be taken as the limiting depth of extinction in this work. For comparing with Leike et al. (2020), the integrated extinction map up to 350pc from this work and Leike et al. (2020) is displayed in Figure 7, where the factor that converts the $G$ band extinction to $E_{\rm G_{BP},G_{RP}}$ is 1.89 (Wang & Chen, 2019). It can be seen that the results from this paper is slightly larger than that from Leike et al. (2020) with a mean difference of $0.04\pm 0.05$mag.

Though the distance-sliced extinction gives the rough range of each cloud, the 25pc step is quite large. The distance to various parts of a molecular cloud can be more accurately determined. The basic method to determine the distance to the cloud is the same as that in Zhao et al. (2020) for supernova remnants (SNRs). It assumes that the extinction will present a sharp increase (i.e. a ‘jump’) at the distance of the dense cloud due to its high dust grain density, and consequently the distance of the cloud is recognized in the variation of the interstellar extinction along the distance. Such model is used to analyze the distance to the extended sources (see e.g. Chen et al. 2017 for MCs and Zhao et al. 2020 for SNRs). Zucker et al. (2019, 2020) used the change of reddening along with parallax to derive the distance to molecular clouds as well.

Consistent with the above analysis, a circle area with a radius of 1.0° is taken for the distance determination. For a given area, the extinction in terms of color excess is the function of distance:

where $E(d)$ is the total color excess along the line of sight until $d$, $E^{\rm fgd}(d)$ is the extinction from the foreground diffuse medium, and $E^{\rm MC}(d)$ is the contribution by the molecular cloud. The $E^{\rm MC}(d)$ is described by a Gaussian error function:

where $\delta E$, $d^{\rm{MC}}$, and $\sigma$ represent the extinction ‘jump’, the distance to and the half thickness superposed onto the distance error of the specific cloud region in the given sightline respectively. However, instead of a two-order polynomial function used in Chen et al. (2017); Zhao et al. (2020) or exponential function in Sun et al. (2021b), a constant is assumed to describe the foreground extinction, i.e.:

This modification is to match the close distance of the clouds, which implies very small foreground extinction. Technically, there are inadequate stars to determine the variation of the foreground extinction within the small distance range, in particular for TMC at $\sim 150$pc.

In addition, only the stars with $E_{\rm G_{BP},G_{RP}}>$0.15mag, i.e. $>5\sigma$, are selected to ensure the ‘jump’ is significant in the given sightline. Because the ‘jump’ is likely to appear at the edge of the cloud, we include the stars at an extended distance range of each molecular cloud (as shown in Table 1) both at the close and the far side by 100pc to guarantee that the jump at the edge of the cloud can be detected.

An MCMC analysis is performed to find the best set of parameters in Equation 4 under the priors of uniform distribution for all parameters, which maximize the likelihood defined as :

where $\theta$ is the parameter to be determined, i.e. $E_{0}$, $\delta E$, $d^{MC}$ and $\sigma$; $E_{i}$ and $E(d_{i}\mid\theta)$ are the color excess calculated from the blue-edge method and the equation parameter to be fitted; $\sigma_{E_{i}}$ is the error of the color excess; $x$ is the data used to fit the equation, including $E_{i}$, $d_{i}$ and $\sigma_{E_{i}}$; and $n$ is the total number of stars in each pixel. We now consider the uncertainty of the distance and combine the distance uncertainty together with the uncertainty of the derived color excess as Chen et al. (2019) (c.f. their Equation 2). With the assistance of the distance from Bailer-Jones et al. (2021), the $\sigma_{E_{i}}$ in the likelihood function is given by:

where $\sigma_{E^{obs}_{i}}$ is the uncertainty of the observed color index, $\sigma_{E^{BlueEdge}_{i}}$ is the uncertainty imported by the blue-edge method as discussed in Section 3, $(E_{i}\frac{\sigma_{d_{i}}}{d_{i}})$ results from the distance uncertainty, which is only an approximation under the assumption that the dust opacity is constant along the line of sight. The distance uncertainties are simply adopted as $\sigma_{d_{i}}=\frac{d_{hi}-d_{lo}}{2}$ where $d_{hi}$ and $d_{lo}$ are the upper and lower bounds in the Bailer-Jones et al. (2021) distance catalog.

Stars that lie within the ranges of both the angular radius and the extended distance of the specific MC are adopted to fit Equation 4, 5 and 6. The MCMC procedure (Foreman-Mackey et al., 2013) is performed to fit the parameters in the model. The ‘burn-in’ chain has 50 walkers and 500 steps to stabilize the chain for final Monte Carlo simulation. Then 3000 steps with 50 walkers are run to estimate the final parameters and their errors. The median values (50th percentile) of the final chain are taken as the best estimation, and the uncertainties equal to the 16th and 84th percentile. Although the input distance range is extended, the result is required to satisfy the following condition so that the ‘jump’ is still in the distance range of the molecular cloud:

where $d^{\rm MC}_{\rm lower}$ and $d^{\rm MC}_{\rm upper}$ mark the lower and upper limit of the MC’s distance, $d^{\rm MC}$ and $\sigma$ mark the distance and the half thickness with error of the specific cloud region.

We adopt the Gelman$-$Rubin statistic to determine whether the fit reaches to convergence. The MCMC fitting is regarded to be converged if the square root of the R hat $(\sqrt{\hat{R}})$ for all the parameters are smaller than 1.01 (Vehtari et al., 2019). We first run the chain 2 times for each sightline and the parameter set with the smaller $\sqrt{\hat{R}}$ for the $d^{\rm{MC}}$ is adopted. If the value of $\sqrt{\hat{R}}$ is smaller than 1.01 for every parameter, this set of parameters is believed to converge and taken as the final results for the fitting in the given sightline, which leads to that about 55$\%$ of the fittings are converged. Moreover, the integrated autocorrelation time $\tau$ is used to check the effective number of independent samples for the converged fittings (typically the chains longer than about 50$\tau$ are sufficient (Foreman-Mackey et al., 2013)). Figure 8 shows the convergence results of the distance ($d^{\rm{MC}}$). The right panel in Figure 8 indicates that almost all (about 99$\%$) of the converged fittings have sufficient independent samples, which strengthens the validity of results selected by R hat.

Figure 9 displays four examples, each for one MC, of the fitting result, i.e. the $\hat{R}$ for the parameter distance ($d^{\rm{MC}}$) and the distribution of the posterior samples of $d^{\rm{MC}}$. It can be seen that the model well follows the trend of the observational points in the expected distance range and the fitting is well converged.

In some sightline, one-cloud model is not perfectly optimized to fit the MCs, there are two or more clouds in some sightlines. We test the two-cloud system with four models as an example, the case sightline is the fourth example in Figure 9 which transverses both the TMC and CMC clouds. In Figure 10, four models are run: (1) The distance range is limited from 100pc to 350pc and one-cloud model is used, then the distance is derived and converged, which detects the TMC cloud; (2) the distance is limited from 300pc to 700pc and one-cloud model is used, then we obtain the distance of convergence for the CMC cloud; (3) the distance is limited from 100 to 700pc and one-cloud model is used, then the result will be similar to the result from 300pc to 700pc, which detects the distance of CMC with the higher extinction and more star tracers; and (4) the distance is again limited from 100pc to 700pc but two-cloud model is used, then two distances, i.e. the distance to the TMC and CMC cloud, are detected, which completely coincide with the one-cloud model with a preselected distance range. This proves that the one-cloud model with a preselected distance range yields the same results as a two-cloud model with no preselection of distance. Also, from the third model which is a two-cloud system but fits with one-cloud model, we find that the one-cloud model will detect the highest extinction cloud or the cloud containing more star tracers if directly used to fit a two-cloud system. The clouds are complexes and may contain multiple cloudlets. Such complexes are presented within a large space range, and the sub-structures mostly disperse in distance. In our analysis, the bin-size in space is only $0.2^{\circ}$ by side, in which we consider no complexes. Indeed, a small distance separation (say less than 5pc) will be difficult to distinguish with the present accuracy.

With the distance and the sliced extinction determined, the 2D extinction map can be decomposed into the individual map for each MC. In Figure 11, 12, 13 and 14 for the four clouds, the left and right panels are for $E_{\rm G_{BP},G_{RP}}$ and $E_{\rm J,K_{s}}$ respectively, while the middle panels show the distance structure with the background contour map of $E_{\rm G_{BP},G_{RP}}$. The blanks denote the position with no reliable result.