Mapping the Milky Way with LAMOST I: method and overview

The probability density distribution (PDF) ${p}_{\mathrm{ph}}$ $(D| c,m,l,b)$ represents the probability of finding a star at distance D at given c, m, l and b. The probability of finding a star in a small volume (D, $D+{\rm{\Delta }}D$) is then written as

$$\begin{eqnarray}&&\begin{array}{l}\Pr ([D,D+{\rm{\Delta }}D]| c,m,l,b)=\\ \,\,\,\,\quad {\displaystyle \int }_{D}^{D+{\rm{\Delta }}D}{p}_{\mathrm{ph}}(D| c,m,l,b){dD}\\ \doteq {p}_{\mathrm{ph}}(D| c,m,l,b){\rm{\Delta }}D.\end{array}\end{eqnarray} \tag{ 1 }$$

On the other hand, $\Pr ([D,D+{\rm{\Delta }}D]| c,m,l,b)$ is also equivalent to the fraction of stars, which can be calculated from the underlying stellar density profile, located within the small volume. Thus, it can also be written as

$$\begin{eqnarray}&&\Pr ([D,D+{\rm{\Delta }}D]| c,m,l,b)={\nu }_{\mathrm{ph}}(D| c,m,l,b)\\ &&\displaystyle \frac{{\rm{\Omega }}{D}^{2}{\rm{\Delta }}D}{{\displaystyle \int }_{0}^{\infty }{\nu }_{\mathrm{ph}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{dD}},\end{eqnarray} \tag{ 2 }$$

where ${\nu }_{\mathrm{ph}}$ is the volume stellar density measured from the photometric stars (i.e., the complete dataset), ${\rm{\Omega }}{D}^{2}{\rm{\Delta }}D$ is the volume element between D and $D+{\rm{\Delta }}D$, and Ω is the solid angle associated with the line-of-sight. Combining Equations (1) with (2), we have

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{ph}}(D| c,m,l,b){\rm{\Delta }}D=\\ \,\displaystyle \frac{{\nu }_{\mathrm{ph}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{\rm{\Delta }}D}{{\displaystyle \int }_{0}^{\infty }{\nu }_{\mathrm{ph}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{dD}}.\end{array}\end{eqnarray} \tag{ 3 }$$

Similarly, for the spectroscopic data, Equation (3) becomes

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{sp}}(D| c,m,l,b){\rm{\Delta }}D=\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\quad \displaystyle \frac{{\nu }_{\mathrm{sp}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{\rm{\Delta }}D}{{\displaystyle \int }_{0}^{\infty }{\nu }_{\mathrm{sp}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{dD}},\end{array}\end{eqnarray} \tag{ 4 }$$

where ${\nu }_{\mathrm{sp}}$ represents the stellar density measured from counting the spectroscopic survey stars only. Within the small region around c and m, the selection function is assumed to be flat with either the color index or magnitude. Thus, it does not change the probability of finding a star in both the photometric and spectroscopic samples, i.e.,

$$\begin{eqnarray}&&{p}_{\mathrm{ph}}(D| c,m,l,b)={p}_{\mathrm{sp}}(D| c,m,l,b).\end{eqnarray} \tag{ 5 }$$

Then, Equations (3) and (4) can be combined together via Equation (5)

$$\begin{eqnarray}&&{\nu }_{\mathrm{ph}}(D| c,m,l,b)={\nu }_{\mathrm{sp}}(D| c,m,l,b){S}^{-1}(c,m,l,b),\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&S(c,m,l,b)=\displaystyle \frac{{\int }_{0}^{\infty }{\nu }_{\mathrm{sp}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{dD}}{{\int }_{0}^{\infty }{\nu }_{\mathrm{ph}}(D| c,m,l,b){\rm{\Omega }}{D}^{2}{dD}}\end{eqnarray} \tag{ 7 }$$

is the selection function at c and m along (l, b).

Then the stellar density profile for all photometric stars along (l, b) can be derived by integrating over c and m

$$\begin{eqnarray}&&\begin{array}{l}{\nu }_{\mathrm{ph}}(D| l,b)=\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\iint {\nu }_{\mathrm{sp}}(D| c,m,l,b){S}^{-1}(c,m,l,b){dcdm}.\end{array}\end{eqnarray} \tag{ 8 }$$

Now consider a stellar sub-population C selected from spectroscopic data based on, for example, the metallicity, luminosity or age. Following Bayes' theorem, the probability of finding a star at D given C, c, m, l and b is written as

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{ph}}(D| C,c,m,l,b)=\\ \,\,\,\,\,\,\,\,\,\quad \displaystyle \frac{{p}_{\mathrm{ph}}(C| D,c,m,l,b){p}_{\mathrm{ph}}(D| c,m,l,b)}{{p}_{\mathrm{ph}}(C| c,m,l,b)}.\end{array}\end{eqnarray} \tag{ 9 }$$

If there is no special selection function for C members at given c and m along (l, b), then we would expect that ${p}_{\mathrm{ph}}(C| D,c,m,l,b)$, i.e. the probability about C at given D, c, m, l and b, should be the same as that for the spectroscopic data, ${p}_{\mathrm{sp}}(C| D,c,m,l,b)$. Subsequently, by integrating over D, we obtain that ${p}_{\mathrm{ph}}(C| c,m,l,b)={p}_{\mathrm{sp}}(C| c,m,l,b)$. Finally, applying Equation (5) to this expression, all the terms of ${p}_{\mathrm{ph}}$ on the right hand side of Equation (9) can be replaced with ${p}_{\mathrm{sp}}$. Applying Bayes' theorem again to the right hand side, we obtain

$$\begin{eqnarray}&&{p}_{\mathrm{ph}}(D| C,c,m,l,b)={p}_{\mathrm{sp}}(D| C,c,m,l,b).\end{eqnarray} \tag{ 10 }$$

Therefore, similar to Equation (6), we infer that

$$\begin{eqnarray}&&\begin{array}{l}{\nu }_{\mathrm{ph}}(D| C,c,m,l,b)=\\ \,\,\,\,\,\,\,\,\,\,\quad {\nu }_{\mathrm{sp}}(D| C,c,m,l,b){S}^{-1}(C,c,m,l,b),\end{array}\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&S(C,c,m,l,b)=\displaystyle \frac{{\int }_{0}^{\infty }{\nu }_{\mathrm{sp}}(D| C,c,m,l,b){\rm{\Omega }}{D}^{2}{dD}}{{\int }_{0}^{\infty }{\nu }_{\mathrm{ph}}(D| C,c,m,l,b){\rm{\Omega }}{D}^{2}{dD}}.\end{eqnarray} \tag{ 12 }$$

Consequently,

$$\begin{array}{ccc}{\nu }_{\mathrm{ph}}(D| C,l,b) & = & \,\,\,\,\,\,\,\,\,\,\,\,\,\iint {\nu }_{\mathrm{sp}}(D| C,c,m,l,b){S}^{-1}\\ & & (C,c,m,l,b){dcdm}.\end{array} \tag{ 13 }$$

In Equations (8) and (13), ${\nu }_{\mathrm{ph}}$, as the density profile for the complete dataset, should be determined. We can estimate it by measuring ${\nu }_{\mathrm{sp}}$ from the spectroscopic data. Usually, along a given line-of-sight, spectroscopic observations can only target a few hundred to a few thousand stars. This implies that ${\nu }_{\mathrm{sp}}$ has to be established from a relatively small dataset. In order to deal with such a situation, we derive ${\nu }_{\mathrm{sp}}$ by using kernel density estimation (KDE). The KDE method can account for both uncertainties in the distance and the relatively small number of stars in a given small region around (c, m). At (l, b), the PDF of D is contributed by all the stars around c and m in terms of

$$\begin{eqnarray}&&{p}_{\mathrm{sp}}(D| c,m,l,b)=\displaystyle \frac{1}{{n}_{\mathrm{sp}}(c,m,l,b)}\underset{i}{\overset{{n}_{\mathrm{sp}}(c,m,l,b)}{\sum }}{p}_{i}(D),\end{eqnarray} \tag{ 14 }$$

where ${p}_{i}(D)$ is the PDF of the distance estimate for the ith star and ${n}_{\mathrm{sp}}(c,m,l,b)$ is the number of spectroscopic stars at a given c, m, l and b. Bringing Equation (14) into Equation (4), we obtain

$$\begin{eqnarray}&&{\nu }_{\mathrm{sp}}(D| c,m,l,b)=\displaystyle \frac{1}{{\rm{\Omega }}{D}^{2}}\underset{i}{\overset{{n}_{\mathrm{sp}}(c,m,l,b)}{\sum }}{p}_{i}(D).\end{eqnarray} \tag{ 15 }$$

Note that $\int {\nu }_{\mathrm{sp}}(c,m,l,b){\rm{\Omega }}{D}^{2}{dD}={n}_{\mathrm{sp}}(c,m,l,b)$, so they are then canceled in Equation (15).

In principle, a KDE determined ν is a continuous function and thus can return a value at any D. However, if there are very few stars nearby D or the given D is beyond the farthest star, the derived ν is not well constrained by observations and thus may no longer be reliable. Therefore, in practice, the stellar density value is sensible only at the positions with sufficient star samples. For convenience, we only adopt ν for the distance at which the spectroscopic stars are located.

It is worth noting that the stellar density derived for a given (l, b) is averaged over the solid angle of the field-of-view, Ω. Therefore, the derived density at the distance where a star is located does not equal the stellar density for the exact 3D spatial position of the star, but is only an approximation of the latter.

Theoretically, S is defined by Equations (7) and (12). Practically, it is not possible to calculate S directly from ${\nu }_{\mathrm{ph}}$ and ${\nu }_{\mathrm{sp}}$. Notice that the integration of ν is actually equivalent to the total number of stars at given c, m, l and b. Therefore, S can be evaluated from

$$\begin{eqnarray}&&S(c,m,l,b)=\displaystyle \frac{{n}_{\mathrm{sp}}(c,m,l,b)}{{n}_{\mathrm{ph}}(c,m,l,b)}.\end{eqnarray} \tag{ 16 }$$

For the case of sub-population C, ${p}_{\mathrm{ph}}(C| c,m,l,b)={p}_{\mathrm{sp}}(C| c,m,l,b)$, and we then have

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{\mathrm{ph}}(C,c,m,l,b)}{{n}_{\mathrm{ph}}(c,m,l,b)}=\displaystyle \frac{{n}_{\mathrm{sp}}(C,c,m,l,b)}{{n}_{\mathrm{sp}}(c,m,l,b)}.\end{eqnarray} \tag{ 17 }$$

Hence, we infer that $S(C,c,m,l,b)=S(c,m,l,b)$.

In order to improve the precision of the computation in Equations (8) and (13), we calculate S⁻¹ instead of S in the rest of the paper.

Before the tests, we construct the selection function to mimic the real spectroscopic survey. The selection of mock spectroscopic stars is based on the infrared (IR) bands J and K. First, we do not induce any selection in the color index $J-K$. Second, the selection function for K is divided into two parts at K = 13 mag. At $K\lt 13\,\mathrm{mag}$, we arbitrarily select stars, while at $K\gt 13$ mag, we apply a flat selection function to keep the mock spectroscopic stars evenly distributed along K. This mimics the strategy that the targeting is more biased to bright rather than faint sources. We denote such a selection function as T1.

We randomly draw 50 groups of mock spectroscopic samples according to T1 and derive the density profile for each group of samples. As a sample, the mock spectroscopic stars from one of the 50 groups are highlighted with red circles in Figure 1, in which the black dots indicate the $J-K$ vs. K distribution for the complete dataset. The right panel shows the normalized distributions of K for the mock spectroscopic stars (red line) and the complete dataset (black line). It can be seen that the red line shows a trend that is similar to the black one at $K\lt 13$ and then becomes flat at $K\gt 13$. This reflects the definition of selection function T1. The normalized distributions of $J-K$ in the bottom panel demonstrate that even though we do not apply any selection function to the color index, the normalized distribution of $J-K$ for the mock spectroscopic stars (red line) is not exactly same as that for the full mock dataset (black line). This is a natural result of the selection function in K. There are more red stars in the regime $K\gt 13$ than in $K\lt 13$. The selection strategy for $K\gt 13$ reduces the sampling rate of fainter and redder stars in the mock spectroscopic stars. This leads to a slightly different distribution in $J-K$ with the complete dataset.

Zoom In Zoom Out Reset image size

Figure 2 shows how the selection function S is calculated for the same group of mock spectroscopic stars as in Figure 1. In the left panel, color represents the number density of stars for the complete sample, i.e., ${n}_{\mathrm{ph}}(J-K,K,l,b)$. The size of each bin is ${\rm{\Delta }}(J-K)=0.1$ and ${\rm{\Delta }}K=0.25$. The plot exhibits a high density at $J-K\sim 0.8$ and $K\sim 14$, which is most likely contributed by K/M dwarf stars. Applying the selection function T1, the map of the number density of mock spectroscopic stars, ${n}_{\mathrm{sp}}(J-K,K,l,b)$, is displayed in the middle panel. The right panel shows the map of ${S}^{-1}(J-K,K,l,b)$, which is obtained by dividing values in the left by those in the middle panel of Figure 2. It is apparent that ${S}^{-1}(J-K,K,l,b)$ is roughly separated into two platforms; one is at $K\lt 13$ with a value of about 5 and the other is at $K\gt 13$ with a value of about 10. This pattern implies that the selection function T1 does not induce selection effects in $J-K$ at each K bin, although the overall distributions in $J-K$ (as shown in the bottom panel of Figure 1) are slightly different than the complete dataset.

Zoom In Zoom Out Reset image size

In order to mimic a more realistic situation, i.e., study the stellar density via a selected stellar population, we select four different populations defined in Table 1 for tests. Table 1 also lists the median numbers and scatter of the 50 groups of arbitrarily selected mock spectroscopic samples based on T1 in the 3rd column. It can be noted that for the sub-population T1S2, which is selected with ${M}_{K}\lt -2$, only about 40 stars are involved in the determination of the stellar density profile. This is a good sample to test the performance of the density estimates for a small dataset.

Table 1.  The Test Populations Selected from the Galaxia Mock Data

Population	Criteria	Number for T1	Number for T2
S1	all spectral types	2230±29	1943±22
S2	${M}_{K}\lt -2$	40±5	23±3
S3	[Fe/H]$\lt -0.3$	979±23	998±22
S4	[Fe/H]$\gt -0.3$	1244±29	946±17

Figure 3 shows the distributions of star counts without a selection function correction for the complete set (black lines) and for one of the 50 groups of the mock spectroscopic samples (red lines). The four panels show the cases for the four different sub-populations, which demonstrate that the selection function T1 does distort the spatial distributions of the mock spectroscopic stars compared to those of the complete samples.

Zoom In Zoom Out Reset image size

The panels in Figure 4 show the derived median stellar densities and their 1σ dispersions (blue circles with error bars) over the 50 groups of mock spectroscopic stars at the spatial positions in which the mock spectroscopic stars are located. In contrast, the corresponding stellar densities for the complete dataset, i.e., the "true" profiles, are superposed as black solid lines in the panels. In the bottom of each panel, the relative ${\rm{\Delta }}\nu /\nu$ is also displayed. It is seen that the derived stellar density profiles from the mock spectroscopic data are perfectly consistent with the "true" profiles. This confirms that when the selection function is very simple, the stellar density can be reconstructed well. In particular for the case T1S2, the accuracy of derived stellar density is still very high even though the sample contains only about 40 stars (see the top-right panel of Figure 4).

Zoom In Zoom Out Reset image size

We now change to a more complicated selection function. Because the LAMOST survey selects its targets using optical-band photometry, the bright/faint limiting magnitude cuts in an optical band lead to a sloped cut in K magnitude. Thus a rectangle region selected in the optical color-magnitude diagram turns out to be a wedge region in IR bands.

Figure 5 provides the selection function for a sample plate² of the LAMOST survey. The middle panel shows that the spectroscopic stars are distributed in a wedge region in the $J-K$ vs. K diagram. Thus the selection function S⁻¹ in the right panel also displays a similar shape.

Zoom In Zoom Out Reset image size

To mimic such a selection function, we randomly select stars satisfying the following criterion

$$\begin{eqnarray}&&-3(J-K)+12\lt K\lt -3(J-K)+15.\end{eqnarray} \tag{ 18 }$$

We denote this selection function as T2. Figure 6 shows a randomly drawn sample based on T2 in the $J-K$ vs. K diagram, in which the selected spectroscopic stars (red circles) are distributed in a wedge area. The median numbers of stars and their scatters over the 50 groups of arbitrarily drawn mock spectroscopic stars based on T2 are separately listed in the last column of Table 1 for the four sub-populations.

Zoom In Zoom Out Reset image size

Figure 7 features the density of stars in the color-magnitude diagram for the same sample group of mock spectroscopic data as in Figure 6. It can be seen that the selected stars exhibit a tilted distribution in the $J-K$ vs. K plane, as does the map of S⁻¹ according to the definition of T2.

Zoom In Zoom Out Reset image size

Figure 8 displays the spatial distributions of star counts without selection function correction for the four sub-populations from one of the 50 groups of mock spectroscopic stars (red lines). In contrast, the black lines indicate the distributions of stars for the corresponding complete datasets. It can be seen that for T2S1, T2S2 and T2S3, distances further than 30 kpc are significantly under-sampled.

Zoom In Zoom Out Reset image size

Figure 9 shows the derived median stellar density profiles and their 1σ dispersions over the 50 randomly drawn groups (blue circles with error bars). In general, the derived profiles are quite consistent with their corresponding "true" values within an uncertainty of ${\rm{\Delta }}\nu /\nu \lt 1$. However, a few systematic biases of about ${\rm{\Delta }}\nu /\nu \sim 0.5-1$ appear in the derived density profiles: one occurs at the distance of $\sim 20-30$ kpc for T2S1, T2S2 and T2S3; the other is at $D\sim 2.5$ kpc for T2S4.

Zoom In Zoom Out Reset image size

Considering that properties of Galactic spatial structures, such as the scale height/length of the disks and the power-law index of the stellar halo, are mostly measured in logarithmic density, such systematic biases may not lead to severe effects. For instance, a simple exponential fit within 0.3–2.5 kpc gives a scale height of 219 pc with the true density profile of T2S4 (black line in the bottom-right panel of Fig. 9), while it decreases to 191 pc obtained from the corresponding derived stellar density profile (blue circles in the same panel). This means that the derived density based on T2 may induce a systematically smaller exponential scale height of only ∼13

For a realistic survey such as LAMOST, many plates are overlapped with each other, and the selection function for different plates varies. Then the possible systematic bias of the stellar density measurement in one plate, as shown in Figure 9, may be compensated by other overlapped plates. Consequently, cumulation of the systematic bias from each individual plate may not induce an overall systematic bias to the final density profiles, but may increase random uncertainty in the resulting density profiles.

In general, the targeting strategy of the LAMOST survey is rather simple. Carlin et al. (2012) designed an elegant targeting algorithm, which tries to make the distribution of target stars flat in color-magnitude planes. However, in practice, this technique has not been fully used for a few reasons. First, the dynamical range of magnitudes observed in a single LAMOST plate is limited to only about 3 magnitudes so that bright spectra will not saturate or cross-talk between fibers. Second, the actual limiting magnitude of $r\sim 18$ mag is brighter than the designed goal by at least 2 mag. Both situations lead to insufficient source stars for targeting. For instance, the bright plates used in the LAMOST survey cover $14\lt r\lt 16.8$ mag, which may only contain less than 200 stars per square degree at high Galactic latitudes. Such a sparse stellar sampling is not sufficient to apply the selection function designed by Carlin et al. (2012) for targeting, because the targeting algorithm needs far more source stars than available fibers to achieve a flat distribution of targets in the color-magnitude plane.

Finally, the LAMOST survey separates the targets into different plates with different ranges of magnitudes in each line-of-sight. The VB plates cover $9\lt r\lt 14$, the B plates cover $14\lt r\lt 16.8$, the M plates cover $r\lt 17.8$ and the F plates cover $r\lt 18.5$. For the VB and B plates, no specific selection function is applied, i.e., the stars are randomly selected. For the M and F plates, we only apply the selection function in r magnitude, and leave the selection in color index to be arbitrary.

This simplified selection function is applied to all of the sky area observed by LAMOST except for the Galactic anti-center region, which covers $150^\circ \lt l\lt 210^\circ$ and $| b| \lt 30^\circ$. For the main survey regions, the UCAC4 catalog (Zacharias et al. 2013) with $r\lt 14$ mag and PanSTARRS-1 catalog (Schlafly et al. 2012) with $r\gt 14$ are adopted as the source catalogs for targeting. For the Galactic anti-center region, the selection strategy follows Yuan et al. (2015), and uses the 2MASS, UCAC4 and Xuyi survey catalogs (Liu et al. 2014b) as the sources for targeting.

Because the different source catalogs are based on different photometric systems, they cannot be unified unless they are cross-calibrated with each other. To avoid the complicated calibrations between at least four different systems, we finally adopt the 2MASS catalog, which covers most of the LAMOST observed stars, as the photometric dataset. Indeed, for LAMOST K giant stars, about 96% of them are within $K\lt 14.3$ mag, which is the limiting magnitude of the 2MASS catalog. Subsequently, we use the $J-K$ vs. K map to derive the selection function S for each observed plate in the LAMOST survey. The selection function S⁻¹ for all LAMOST DR3 plates and the python codes used for measurement of the stellar density profiles can be found at https://github.com/liuchaonaoc/LAMOST_density.

In order to derive the density profiles for the stellar disk and halo of the Galaxy, we select RGB stars. Although RGB stars are not the standard candles that red clump stars are, the accuracy of distance using the former ones is still around 20%–30% (Liu et al. 2014a, Carlin et al. 2015), which is sufficient for mapping Galactic structure in the outer regions. Moreover, RGB stars are intrinsically brighter and occur in stellar populations with a broader range of metallicity than red clump stars. Therefore, they are suitable tracers for the outer disk and halo, which are dominated by metal-poor stars.

For LAMOST survey data, we adopt stellar parameters, i.e., the effective temperature, surface gravity and metallicity, provided by the LAMOST pipeline (Wu et al. 2011, 2014; Luo et al. 2015). K giant stars, including the RGB and red clump stars, are selected based on criteria proposed by Liu et al. (2014a). Red clump stars are identified by Wan et al. (2015) and Tian et al. (2016). The RGB stars are then selected from K giants by excluding the red clump stars.

The absolute magnitudes at K_s band, M_Ks, for the RGB stars are estimated following the method suggested by Carlin et al. 2015. Interstellar extinction for individual RGB stars is estimated using the Rayleigh-Jeans color excess (RJCE) approach introduced by Majewski et al. (2011) and later revised to incorporate WISE bandpasses by Zasowski et al. (2013). Subsequently, the distances to RGB stars are estimated by combining the apparent K_s magnitude, the absolute magnitude M_Ks and the RJCE extinction A_Ks.

The RJCE extinction is compared with the 3D extinction map provided by Green et al. (2014, 2015), who derived the 3D dust map from PanSTARRS-1 data.

Figure 10 shows the RJCE extinction vs. distance for the RGB stars located at ($150^\circ \lt l\lt 160^\circ$, $| b| \lt 2^\circ$), ($170^\circ \lt l\lt 190^\circ$, $| b| \lt 2^\circ$) and ($200^\circ \lt l\lt 210^\circ$, $| b| \lt 2^\circ$) using circles, crosses and pluses, respectively. The Green et al. mean extinction-distance relations at similar lines-of-sight are superposed in the same figure. We see that the RJCE extinctions for individual RGB stars are roughly consistent with Green's 3D extinction map within about 0.2 mag, which may lead to only about 10% uncertainty in distance. Therefore, the distances for the RGB stars should not be substantially affected by extinction, even in the Galactic mid-plane.

Zoom In Zoom Out Reset image size

In the following subsections we demonstrate two samples: (1) RGB stars with all ranges of metallicity as tracers for the Milky Way structures, especially for the Galactic disk, and (2) metal-poor RGB stars as tracers for the stellar halo.

We select 21 954 RGB stars with ${M}_{{Ks}}\lt -3.5\,\mathrm{mag}$ and distances larger than 0.5 kpc as probes to map the Galaxy. Because these samples are dominated by metal-rich stars, they mostly trace the Galactic disk. The volume completeness for these RGB stars is within 40 kpc. The volume completeness in this paper is defined such that the range of luminosity for the stars does not significantly change within the volume. Plate by plate, we apply the method to derive stellar density and finally obtain density values for the distances at which sample stars are located.

Table 3 lists the first three rows of the samples. The columns in the table are described in Table 2. The full RGB data can be found at https://github.com/liuchaonaoc/LAMOST_density.

Table 2.  Explanation of Data Columns in the Catalogs used for Disk and Halo Density Mapping

Column name	Type	Unit	Comments
obsid	integer		The identifier of the LAMOST spectra
RA	double	degree	Right ascension in epoch J2000.0
DEC	double	degree	Declination in epoch J2000.0
${T}_{\mathrm{eff}}$	float	K	Effective temperature from LAMOST DR3 catalog
log g	float	dex	Surface gravity from LAMOST DR3 catalog
[Fe/H]	float	dex	Metallicity from LAMOST DR3 catalog
vlos	float	km s−1	Line-of-sight velocity from LAMOST DR3 catalog
MKs	float	mag	Ks-band absolute magnitude from Carlin et al. (2015)
MKs lower error	float	mag	−1σerror of MKs
MKs upper error	float	mag	+1σerror of MKs
AKs	float	mag	Ks-band interstellar extinction derived from RJCE (Majewski et al. 2011,
			Zasowski et al. (2013))
Ks	float	mag	Ks-band apparent magnitude from 2MASS catalog
Distance	float	kpc	Distance estimated from MKs, AKs and Ks
Distance lower error	float	kpc	−1σerror of distance
Distance upper error	float	kpc	+1σerror of distance
Z	float	kpc	Vertical distance from the star to the Galactic mid-plane
R	float	kpc	Galactocentric radius in cylindrical coordinates
r	float	kpc	Galactocentric radius in spherical coordinates
$\mathrm{ln}\nu$	float	ln kpc−3	Logarithmic stellar density at the position of the star

Table 3.  The Disk-Like RGB Stars Selected from the LAMOST DR3 Catalog

obsid	RA	DEC	${T}_{\mathrm{eff}}$	log g	[Fe/H]	${v}_{\mathrm{los}}$	MK	AK	K	Distance	Z	R	r	$\mathrm{ln}\nu$
	(°)	(°)	(K)	(dex)	(dex)	(km s−1)	(mag)	(mag)	(mag)	(kpc)	(kpc)	(kpc)	(kpc)	$\mathrm{ln}$(pc ${}^{-3}$)
75606117	337.46744	+6.01526	4829	1.11	−1.66	−3	$-{3.96}_{-0.44}^{+0.42}$	0.148	10.646	${7.79}_{-1.37}^{+1.75}$	−5.22	8.26	9.77	−11.40
269507178	3.66002	+55.72871	4159	1.47	−0.47	−93	$-{3.82}_{-0.36}^{+0.24}$	0.252	9.721	${4.55}_{-0.48}^{+0.82}$	−0.51	10.86	10.88	−6.25
51211049	257.70419	+18.38338	4207	2.64	−0.51	−33	$-{3.32}_{-0.42}^{+0.38}$	0.072	9.687	${3.86}_{-0.62}^{+0.82}$	1.97	5.82	6.14	−7.51

Notes: Here we only list the first three rows; the complete table is in the online file (https://github.com/liuchaonaoc/LAMOST_density).

Note that the samples contain about 20% of stars that have been observed more than twice. Because the duplicated stars were observed in different plates and selection effects for each observation can be considered separately, the duplication would not affect stellar density estimation. In fact, these data can be used to test the internal uncertainty, or the precision, of the distance and stellar density estimations. The left panel of Figure 11 shows the distribution of relative scatter in the distance estimates (solid line) of stars observed multiple times. The scatter is defined by the standard deviation of the distance. The relative scatter is the scatter divided by the mean distance estimate. The distribution of relative scatter in distance is fitted with a Gaussian (dashed line) with the best fit parameter of σ = 0.038, while the median value of the distribution is 0.028. This means that the internal uncertainty of distance estimates for the RGB stars is around 3%–4%.

Zoom In Zoom Out Reset image size

The right panel of Figure 11 shows the distribution of relative scatter in stellar density (solid line) for the duplicated stars. The definitions of scatter and relative scatter are similar to those for the distance estimates. The dashed line stands for the best fit Gaussian with σ = 0.37. The median value for the distribution is 0.25. This is better than the performance of the mock data test described in sections 3.1 and 3.2, which may reach ${\rm{\Delta }}\nu /\nu \gt 0.5$ in the worst cases. The difference in the derived stellar densities for multiply-observed stars originates from two possible channels: (1) the different selection corrections from the different plates observing the same star and (2) the different averaged sky areas since the plates observing the same star may not be completely overlapped with each other. The small dispersion of the relative scatter implies that these differences are very small for most of the stars.

Figure 12 demonstrates the averaged stellar density map in the R–Z plane³. The bin size is 0.5 × 0.5 kpc. Because this population is dominated by metal-rich stars, it displays a remarkable tomographic map of the Galactic disk from $R\sim 4$ to ∼20 kpc. Two prominent features are seen in this figure. First, the Monoceros ring, unveiled by Newberg et al. (2002) and later highlighted in SDSS star counting (Jurić et al. 2008) and metallicity distribution mapping (Ivezić et al. 2008), does not show up in our map of the whole Galactic outer disk with LAMOST RGB stars. Second, the iso-density contours (black lines) show moderate north-south asymmetry in the outer part of the disk starting from $R\sim 12$ kpc. The stellar density above the mid-plane is larger than that below at given radii. This is essentially consistent with Xu et al. (2015). The asymmetry about the mid-plane may also be the result of the warp. If that is the case, the stellar disk may bend up to some extent in the Galactic anti-center region between $150^\circ \lt l\lt 210^\circ$, in which the LAMOST disk data are mostly concentrated. Both features are very interesting and will be further discussed in an upcoming paper (Wang et al. in preparation).

Zoom In Zoom Out Reset image size

In Figure 13, black circles show the stellar surface density derived by integrating the volume density over Z in each R slice. In order to avoid issues due to a lack of data in the southern Galactic hemisphere, we use $| Z|$ instead of Z. That is, the mean stellar density at given $| Z|$ is averaged over the values at both +Z and −Z. Only $| Z| \lt 40\,\mathrm{kpc}$ is included in the integration. The width of each R slice is 1 kpc. For some bins without stellar density estimates, interpolated values are used. The red and blue dashed lines are best-fit exponential disk and power-law halo profiles, respectively, and the green dot-dashed line is the sum of the two profiles. We find that the surface density profile can be fitted well with an exponential disk with scale length of $1.6\pm 0.1\,\mathrm{kpc}$ and a power law halo with index of $2.4\pm 0.3$. Note that this power-law index is not equivalent to that derived from spherical coordinates, since it is derived from the surface density as a function of cylindrical radius R.

Zoom In Zoom Out Reset image size

The fraction of stars contributed by the halo, as derived from our models of the surface density, is indicated as the black line which is scaled by the right y-axis in the figure. It shows that, at R = 14 kpc, about half of the surface density is contributed by the disk. Therefore, no evidence is found to support that there is a truncation at around 14 kpc. Moreover, at R = 19 kpc, the disk still contributes ∼10% to the total stellar surface density, which cannot be negligible because the data points (black circles) are clearly larger than the halo model (blue lines) around the radius. In other words, the disk extends to as far as 19 kpc. Beyond this radius, the observed surface density smoothly transitions to a pure halo profile. The perfect fitting with an exponential disk model implies that the disk displays neither truncation, nor a break, nor an up-bending feature.

The extended larger size of the disk is not likely to be an artificial effect due to incorrect interstellar extinction, according to the discussion in section 4.2. On the other hand, Carlin et al. (2015) pointed out that the distances may be systematically overestimated for α-enhanced giant stars by about 15%. However, Hayden et al. (2015) show that it is the α-low giant stars that dominate the outer disk. Consequently, the systematic bias in distance due to the α abundance may not significantly affect the distance for the outer disk RGB stars. Therefore, we suggest that the largely extended disk without truncation should be real.

The Galactic disk measured in this work is significantly larger than that in some other works, which claimed that the Milky Way stellar disk truncates at around 14–15 kpc (e.g. Robin et al. 1992, Minniti et al. 2011). Carraro (2015) pointed out that the lines-of-sight investigated by these works are very close to $b\sim 0$ and may not be aligned with the disk plane in the outskirts due to the existence of the warp. In this case, at some radius the in-plane lines-of-sight no longer probe the disk, leading to a false truncation at that radius. It is worth noting that Carraro et al. (2010) and Feast et al. (2014) found some young stars at around 20 kpc from the Galactic center, which is in agreement with our suggestion that the disk extends to such radii. We also point out that Xu et al. (2015) suggested the disk may extend all the way to large radii. More quantitative analysis will appear in an upcoming work by Wang et al. (in preparation).

In this section we select 5171 halo-like RGB stars with [Fe/H]$\lt -1$ and ${M}_{{Ks}}\lt -4$ from the LAMOST DR3 catalog to map the stellar density profile for the stellar halo. The cut in K_s-band absolute magnitude ensures that the data are roughly complete within 50 kpc from the Galactic center. The first three rows of the sample are listed in Table 4. The columns in the table are explained in Table 2. The full catalog can be found at https://github.com/liuchaonaoc/LAMOST_density. We then derive stellar densities at the spatial positions at which the halo-like RGB stars are located.

Table 4.  The Halo-Like RGB Stars Selected from the LAMOST DR3 Catalog

obsid	RA	Dec	${T}_{\mathrm{eff}}$	$\mathrm{log}g$	[Fe/H]	${v}_{\mathrm{los}}$	MK	AK	K	Distance	Z	R	r	$\mathrm{ln}\nu$
	(deg)	(deg)	(K)	(dex)	(dex)	(km s−1)	(mag)	(mag)	(mag)	(kpc)	(kpc)	(kpc)	(kpc)	$\mathrm{ln}$(pc ${}^{-3}$)
75606117	337.46744	+6.01526	4829	1.11	−1.66	−3	$-{3.96}_{-0.44}^{+0.42}$	0.148	10.646	${7.79}_{-1.37}^{+1.75}$	−5.22	8.26	9.77	−11.81
342809095	263.45129	+18.10553	3870	−0.24	−2.17	−22	$-{6.12}_{-0.02}^{+0.02}$	0.053	11.280	${29.47}_{-0.27}^{+0.27}$	12.49	21.36	24.74	−13.47
318815093	233.31128	+12.40625	4386	0.53	−2.32	172	$-{5.32}_{-0.16}^{+0.28}$	0.287	10.854	${15.04}_{-1.82}^{+1.15}$	11.45	3.53	11.98	−12.43

Notes: Here we only list the first three rows; the complete table is in the online file (https://github.com/liuchaonaoc/LAMOST_density).

Figure 14 shows the averaged stellar density map in the R–Z plane for the halo population. The metal-poor RGB stars probe the shape of the halo well within a Galactocentric radius of 50 kpc. The contours of $\mathrm{ln}\nu =-12.5$ and −13.5 in the figure show clear oblate shapes within 20 kpc. The iso-density contour at −14.5 displays a slightly larger axis ratio, although it is still oblate. Beyond 30 kpc, the contours of −15.5 and −16.5 are roughly spherical. Although Xue et al. (2015) claimed that a single power-law halo with variable axis ratio can fit their data better than the broken power-law, this work is the first time, to the best of our knowledge, that the variable axis ratio from inner to outer halo is directly illustrated in the tomographic map.

Zoom In Zoom Out Reset image size

We further demonstrate how the broken power-law appears with the assumption of a constant axis ratio. Figure 15 shows $\mathrm{ln}\nu$ as a function of r for the points shown in Figure 14 with different fixed values of the axis ratio. The top panel presumes the axis ratio q = 1.0. Thus r is the same as the Galactocentric radius in spherical coordinates. The $\mathrm{ln}\nu$ profile is fitted with a broken power-law. Within $\sim 40\pm 2.8\,\mathrm{kpc}$, the best fit power-law index is $-2.9\pm 0.1$, but beyond this it is down to $-5.0\pm 1.2$. In the bottom panel, q is set to be 0.75. In this case, the broken power-law profile shown in the top panel almost disappears. Instead, it shows that the two power-law indices in the inner and outer halos are similar, i.e., a single power-law may be better to fit the data. This means that the presumption of the axis ratio can significantly change the radial profile of $\mathrm{ln}\nu$. More quantitative study about the structure of the stellar halo with LAMOST RGB samples can be found in Xu et al. (2017).

Zoom In Zoom Out Reset image size

In this work, we derive the density profile along a line-of-sight for each LAMOST plate. This implies that the stellar density estimate is smoothed over the solid angle of 20 square degrees, since the field-of-view of the LAMOST telescope is 5 degrees in diameter. Note that the stellar density profile may dramatically change from one end to another in such a large sky area, especially when the distance is large. As a consequence, the plate-wide smoothing density profile may smear out or blur the structure of the Milky Way to some extent. In this section we investigate the influence of plate-wide smoothing on the density profile with toy models.

We adopt a toy model of the vertical stellar density profile with two exponential components: the thin disk component has a scale height of 0.3 kpc and the thick disk has a scale height of 0.9 kpc. At Z = 0, the thick disk component is 15% of the thin one. The normalizer, i.e. the stellar density at Z = 0, is set to 0.05 pc⁻³. Such a configuration is similar to the Milky Way disk in the solar neighborhood.

Considering that LAMOST plates usually overlap each other in most of the surveyed sky areas, we draw mock stars from the toy model in the partly superposed LAMOST-like plates, i.e. circles with 5-degree-field-of-view, centered at b = 0°, 1°, 2°, $\ldots$, 12°, 22°, $\ldots$, 82° with l = 180°. Obviously, these plates correspond to different vertical heights above the Galactic mid-plane; the larger the distance, the larger the vertical height.

For each plate, at a given Galactocentric radius, we randomly draw 50 mock stars evenly distributed within the circular area of the plate. Note that, even though the mock stars have the same radius, their Z values are different because they are distributed at various Galactic latitudes within the plate. According to section 2, the stellar densities at the positions of the mock stars are derived by averaging over the model density covering the area of the plate at the given radius. Thus, although the Z values are different for the mock stars drawn from the same plate at a given radius, the stellar densities assigned to them are the same. This is a similar but simplified way to determine the stellar density as compared to the real observed data applied in section 4. Although uncertainties in the distance estimates and the selection function corrections are ignored, the simplified method for stellar density determination is sufficient for investigating the effect induced by plate-wide smoothing.

5.1.1. The smoothing effect in the vertical disk profile

We first investigate the effect in the vertical density profile. Figure 16(a) shows the result of the vertical density profile test at R = 12 kpc. The black solid line indicates the "true" profile from the toy model. The red filled circles represent the averaged density values over the mock stars at different Z bins. The red dashed line displays the best-fit double-exponential profile with red symbols. It is seen that the derived stellar densities perfectly follow the "true" stellar density profile and the best-fit parameters are quite consistent with the "true" values with uncertainties less than 20%.

Zoom In Zoom Out Reset image size

We then move to a larger distance at R = 20 kpc, where the effect of plate-wide smoothing should be more substantial. The result in Figure 16(b) shows that the shape of the best-fit density profile (red dashed line) derived from the "observed" data (red filled circles) is in agreement with the "true" profile (black solid line). Similar to the test at R = 12 kpc, the uncertainty of the derived model parameters is less than 20% without substantial systematic bias. In both cases at R = 12 and 20 kpc, plate-wide smoothing due to the large field-of-view does not significantly change the result.

5.1.2. The smoothing effect in the Monoceros ring-like substructure

5.1.3. Summary

The selection function used in this work is determined by comparison between the color-magnitude diagram from the photometric data and that from the final spectroscopic catalog. The selection function in the latter one is not only composed of targeting selection, but also considers selection effects induced during observations and data reductions. The original selection function is always altered during observations and data processing. For instance, one line-of-sight is designed to be observed multiple times with different ranges of magnitudes so that a larger range of magnitude can be covered. But due to weather, some designed plates may never be observed. Then the post-observation selection function for the line-of-sight would never be identical to the designed one. The site conditions, which allows more observation time in winter but less in summer, can also induce selection function in the sky area (Yao et al. 2012). Another instance is that the data reduction may lose spectra with very low signal-to-noise ratio, which may more frequently occur for fainter sources. Therefore, only using the originally designed targeting selection for the correction cannot take into account effects from observations and data reductions.