Mapping the Milky Way with LAMOST II: the stellar halo

TL;DR

This study measures the density and shape of the Milky Way's stellar halo using LAMOST K giants, revealing a single power-law profile and radius-dependent flattening, along with identifying tidal streams and an unknown overdensity.

Contribution

It provides a model-independent measurement of the stellar halo's density profile and flattening variation with radius using a large spectroscopic dataset.

Findings

Halo density follows a power law with index ~5.

Halo flattening varies from 0.64 to 0.96 with radius.

Identifies tidal streams and an unknown overdensity.

Abstract

The radial number density and flattening of the Milky Way's stellar halo is measured with $\mathrm{5351}$ metal-poor ([Fe/H]$<-1$) K giants from LAMOST DR3, using a nonparametric method which is model independent and largely avoids the influence of halo substucture. The number density profile is well described by a single power law with index $5.03^{+0.64}_{-0.64}$, and flattening that varies with radius. The stellar halo traced by LAMOST K giants is more flattened at smaller radii, and becomes nearly spherical at larger radii. The flattening, $q$, is about 0.64, 0.8, 0.96 at $r=15$, 20 and 30 kpc (where $r=\sqrt{R^2+\left[Z/q\left(r\right)\right]^2}$), respectively. Moreover, the leading arm of the Sagittarius dwarf galaxy tidal stream in the north, and the trailing arm in the south, are significant in the residual map of density distribution. In addition, an unknown overdensity is identified in the residual map at (R,Z)=(30,15) kpc.

Tables

Table 1. Table 1: Incomplete list of recent stellar halo profile fits

reference	origin	tracer	sample size	distance(kpc)	model	parameters
Iorio (2017)	GAIA+2MASS	RR Lyrae	21600		Triaxial	n=2.96, p=1.27, $q=f\left(r\right)$, $q_0=0.57$, $q_{inf}=0.84$,
						$r_0=12.2kpc$
Das et al. (2016)	SEGUE2	BHB			BPL	$n_{in}=3.61$, $n_{out}=4.75$, $r_{break}=29.87$, q=0.72
					SPL	n=4.65, $q=f\left(r_{GC}\right)$, $q_0=0.39$, $q_{inf}$=0.81, $r_0=7.32Kpc$
X15	SEGUE2	K giants	1757		BPL	$n_{in}=2.8$, $n_{out}=4.3$, $r_{break}=29$, q=0.77
					Einasto	n=2.3, $r_{eff}=18$, q=0.77
					SPL	n=4.4, $q=f\left(r_{GC}\right)$, $q_0=0.3$, $q_{inf}$=0.9, $r_0=9Kpc$
Pila-Diez et al. (2015)	CFHTS & INT	near MSTO			SPL	n=4.3, q=0.79
					Triaxial	n=4.28, q=0.77, $\omega =0.87$
					BPL	$n_{in}=2.4$, $n_{out}=4.8$, $r_{break}=19$, q=0.77
					BPLq	$n_{in}=3.3$, $n_{out}=4.9$, $q_{in}=0.7$, $q_{out}=0.88$
Deason et al. (2011)	SDSS DR8	BS,BHB	$\sim 20000$		BPL	$n_{in}=2.3$, $n_{out}=4.6$, $r_{break}=27$, q=0.6
Deason et al. (2014)	SDSS DR9	BS,BHB			BPL	$n_{outer}=6-10$, $r_{break}=50$
Watkins et al. (2009)	Stripe82	RRly	417		BPL	$n_{in}=2.4$, $n_{out}=4.5$, $r_{break}=25$
Sesar et al. (2011)	CFHTLS	near MSTO	27544		BPL	$n_{in}=2.62$, $n_{out}=3.8$, $r_{break}=28$, q=0.7
Juric et al. (2008)	SDSS	MS			SPL	$n=-2.8$, q=0.64
Bell et al. (2008)	SDSS	MS	4 million		SPL	,
Siegel et al. (2002)	Kapteyn		70000		SPL	$n=2.75,q=0.6$
Robin et al. (2000)	PB				SPL	$n=2.44,q=0.76$

Table 2. Table 2: Best-fit parameters r 𝑟 r , q 𝑞 q and their uncertainties for each iso-density contour

$\ln \left(\nu \right)$	number	$r$	${\sigma }_r$	$q$	${\sigma }_q$
-11.50	165	12.095	0.969	0.468	0.068
-11.75	200	11.735	0.893	0.546	0.078
-12.00	203	12.595	1.034	0.511	0.074
-12.25	221	12.988	1.083	0.543	0.081
-12.50	246	13.904	0.842	0.565	0.064
-12.75	255	13.203	0.761	0.681	0.077
-13.00	244	15.467	0.750	0.609	0.057
-13.25	261	15.656	0.676	0.683	0.058
-13.50	268	17.375	0.708	0.654	0.051
-13.75	298	17.572	0.641	0.728	0.052
-14.00	271	17.361	0.629	0.803	0.060
-14.25	261	19.748	0.636	0.783	0.051
-14.50	267	19.857	0.610	0.796	0.049
-14.75	224	21.496	0.613	0.814	0.047
-15.00	221	23.059	0.614	0.837	0.045
-15.25	216	23.185	0.578	0.903	0.047
-15.50	192	24.076	0.554	0.974	0.047
-15.75	140	25.141	0.628	0.999	0.050
-16.00	106	29.935	0.602	0.823	0.034
-16.38	176	29.989	0.612	0.910	0.046
-16.88	81	31.349	0.641	1.043	0.052
-17.38	52	36.076	0.677	1.010	0.045

Equations

$$p_{ph}(D|c,m,l,b)=p_{sp}(D|c,m,l,b),$$

$${\nu }_{ph}(D|c,m,l,b)={\nu }_{sp}(D|c,m,l,b)S^{-1}(c,m,l,b),$$

$${\nu }_{ph}(D|l,b)=\iint {\nu }_{sp}(D|c,m,l,b)S^{-1}(c,m,l,b)dcdm.$$

$$S(c,m,l,b)=\frac{n_{sp}(c,m,l,b)}{n_{ph}(c,m,l,b)},$$

$${\nu }_{sp}(D|c,m,l,b)=\frac{1}{\Omega D^2}\sum _i^{n_{sp}(c,m,l,b)}{p}_i(D),$$

$$\ln (\nu )\propto n/2ln[R^2+(Z/q{)}^2]$$

$$\ln (\nu )\propto nln(Z)-nln(q)$$

$$\ln (\nu )\propto nln(R)$$

$$\begin{array}[]{ll}R=r\cos(\eta),\\ Z=rq\sin(\eta),\end{array}$$

$$r_{GC}=\sqrt{R^2+Z^2}.$$

$$r_{GC}=r\sqrt{1+(q^2-1)si{n}^2(\eta )}.$$

$$\sin (\eta )=\frac{r_{GC}\sin (\theta )}{aq}.$$

$$r_{GC}=rq\sqrt{\frac{1}{q^2-(q^2-1)si{n}^2(\theta )}}.$$

Full text

††pubyear: 2017††pagerange: Mapping the Milky Way with LAMOST II: the stellar halo–A

1 Introduction

The spatial distribution of the halo tracers, specifically the radial number density profile and the flattening of the tracer population, is crucial for understanding the formation history and dark matter distribution in the Milky Way. However, there is little consensus on the radial density and flattening profile of the stellar halo.

In the early studies, based on a few hundred objects (i.e. globular clusters, RR Lyrae variables, blue horizontal branch stars, and K dwarfs, etc.), radial number density profile of the Milky Way halo was usually described as a single power law (SPL) with index of $2.5-3.5$ and the flattening either kept constant $q=\mathrm{0.5-1}$ (Harris, 1976; Hawkins, 1984; Sommer-Larsen, 1987; Preston et al., 1991; Soubiran, 1993; Wetterer et al., 1996; Reid et al., 1998; Gould et al., 1998; Vivas, 2006) or varied with radius Preston et al. (1991). However, other studies explored a deviation from the single power law. For instance, Saha (1985) found the density profile traced by RR Lyrae stars followed a broken power law (BPL) with an index of 3 within 25 kpc and an index of 5 beyond 25 kpc.

With the advent of large sky suveys, the halo density profile can be studied using an expanded sample with larger sky coverage and out to larger distances. Measurements from recent literature are summarized in Table 1. These studies show that more complicated, multi-parameter models of the global number density profile are often required to fit the expanded datasets. Some researchers found that the density profile follows a single power law with a constant flattening (Robin et al., 2000; Siegel et al., 2002; Bell et al., 2008; Juric et al., 2008), while other studies claimed that a broken power law or equivalent Einasto profile (Einasto, 1965) with a constant flattening is more suitable to describe the observed data (Sesar et al., 2011; Deason et al., 2011). Deason et al. (2014) found that there are two breaks located at 27 kpc and 50 kpc, respectively, in the density profile. Xue et al. (2015, herafter X15) clarified that a broken power law or Einasto profile with constant flattening can fit the density profile of SEGUE K giants or, alternatively, a single power law with variable flattening can also fit the data well . They attributed the apparent radial change in power law index to the variation of flattening. Pila-Diez et al. (2015) also found the halo is more flattened within 19 kpc, and becomes less flattened beyond 19 kpc.

The results of previous works listed in Table 1 come with some caveats. First, some tracers (such as RR Lyrae stars) are quite rare objects in the Galactic halo, with total sample size of only a few hundreds. Second, different tracers may trace different stellar populations. For example, since on the horizontal branch is correlated with metallicity, BHBs and RR Lyraes represent different stellar populations. Finally, for all populations, the density profile is not well measured at large distances, preventing tight constraints of the shape of the outer halo.

Furthermore, the power law index and flattening parameters in the density models are degenerate to some degree. X15 found that the SDSS K giants can be empirically characterized using several models with constant or variable flattening, but it is impossible to distinguish which one is correct due to degeneracy between the models. Therefore, it is crucial to determine power law index and flattening independently.

An important consideration in determining the density profile of the stellar halo is that a large fraction of halo stars are associated with substructure. Watkins et al. (2009) indicated that 60% of the stars in their sample belong to substructure. Sesar et al. (2010) found that 20% of the stars within 30 kpc belongs to substructure, and beyond that distance the fraction is larger. X15 found that the substructure does affect the halo density profile. However, it is difficult to identify each individual member star of a given substructure. Fitting the density profile of the smooth component alone is often accomplished by cutting out regions of the sky that contain significant density substructure to avoid their influence. The problem of this approach is that most significant substructures (such as the Sgr stream and Virgo overdensity) are very extended. As a consequence, if we exclude the Sgr stream and Virgo overdensity according to the Bell et al. (2008) criteria, most stars higher than $b=60^{\circ}$ are masked. This makes it difficult to constrain $q$, since $q$ is sensitive to density comparisons between the low and high Galactic latitudes.

Liu et al. (2017 hereafter Paper I) developed a statistical method to derive the stellar density profile from the spectroscopic survey data. They qualitatively found that the metal-poor ([Fe/H]$<$-1) halo-like RGB stars have a quite oblate shape within a Galatocentric radius of about 20 kpc. When looking at radii larger than 30 kpc, the shape of the halo density profile becomes nearly spherical. This is evidence that the flattening of the stellar halo varies with radius, rather than following a broken power law.

In this work, we quantitatively investigate this with more carefully selected K giant samples and a well defined approach. A nonparametric technique is adopted to measure the density profile of the Galactic halo with The Large sky Area Multi-Object fiber Spectroscopic Telescope (LAMOST) Survey metal poor ([Fe/H]$<-1$) K giants as far as 35 kpc. Using this method, we fit the iso-density surface of LAMOST K giants directly without presuming functional forms that describe the density and variation of flattening. In addition, the influence of the outliers caused by substructure will be eliminated by fitting the median value of the iso-density surface if the substructure do not dominate in a particular bin. So the advantage of the nonparameteric method is that it can estimate the power law index and the flattening independently, while avoiding the influence of substructures.

In next section, we describe the K giant sample and the stellar density estimation. The nonparametric technique and the halo profile are presented in section 3. Finally, in section 4 we discuss the substructures that were not included in our fit, and draw brief conclusions.

2 LAMOST K giants and the stellar density

2.1 LAMOST K Giants

LAMOST is a 4-meter, quasi-meridian, reflecting Schmidt telescope which was designed for a large spectroscopic survey(Cui et al., 2012; Zhao et al., 2012; Deng et al., 2012). The stellar parameters are estimated in the LAMOST pipeline (Wu et al. 2011, 2014).

The LAMOST Data Release 3 (DR3) catalogue contains over 4 million stellar spectra, of which 350,000 stars are identified as K giants (Liu et al., 2014). Contamination from red clump stars is removed using the technique for identifying red clump stars provided by Wan et al. (2015), Tian et al. (2016), and Wan et al. (2017). The distances of these K giants are derived based on the fiducial isochrones in SDSS extinction-corrected $g-r$ $\it vs.$ absolute r band magnitude, $M_{r}$, Pan-STARRS1 $g_{p1}$ and $r_{p1}$, and LAMOST [Fe/H], using Bayesian method described in Xue et al. (2014). The photometry for the K giants is corrected for extinction based on $E(B-V)$ of Schlegel et al. (1998), and the coeficients in Table 6 of Schlafly et al. (2011). K giants with $E(B-V)>0.25$ are eliminated to avoid the influence of poor extinction correction in areas of high extinction; the typical $g-r$ error caused by the extinction correction is about 0.017 using this criterion (Xue et al. 2014). The errors in colour, magnitude, and metallicity are propagated into the distance uncertainty. The systematic distance bias which is caused by the difference between SDSS $g$ and $r$ and Pan-STARRS1 $g_{p1}$ and $r_{p1}$ is calibrated using $\sim$ 400 common K giants, as shown in Figure 1. Including all of these factors, the distance estimates for LAMOST K giants have a typical uncertainty of about 15% .

To eliminate contamination from disk stars, we cull K giants with [Fe/H]$<-1$. Because the spectroscopic selection effects are corrected based on the photometry from The Two Micron All-Sky Survey (2MASS, Cutri et al. 2003), we impose a faint magnitude limit of $K<14^{m}.3$. In addition, the absolute $r_{p1}$ band magnitude, $M_{r,p1}<-0.5$, is restricted in order to maintain the completeness of $M_{r,p1}$ up to distances of 35 kpc. We also constrain $E(B-V)<0.25$ in our K giant samples, to avoid unreliable reddening at low latitudes. After all of these cuts, there are 5351 K giants with [Fe/H]$<-1$, $M_{r,p1}<-0.5$, $K<14^{m}.3$, and $E(B-V)<0.25$ in our sample. Figure 2 shows the sky coverage of the selected K giant sample. The red points are eliminated from our sample because their $E(B-V)$ is larger than 0.25, which excludes the stars in the range of $-10^{\circ}<b<10^{\circ}$. But our sample still has a sizeable sky coverage in the range $b=-60^{\circ}\sim 90^{\circ}$, which is sufficient for constraining the global structure of the Galactic stellar halo. Figure 3 illustrates the distribution of the K giants in the $M_{r}$ vs. distance plane. The total sample, with $M_{r,p1}$ in the range $[-4^{m},-0.5^{m}]$ is barely complete up to a distance of 35 kpc (The completeness here is defined as a volume within which the coverage of the absolute magnitude of the samples should be similar at any place.). Beyond 35 kpc, we begin to lose intrinsically fainter stars. In addition, the number of available stars declines dramatically. Therefore, in order to keep a relatively large samples in a relatively large distance, we select a lower limit of $M_{r,p1}=-0.5^{m}$. A lower limit fainter than this would mean that we may not have complete coverage out to 35 kpc, which is critical for investigating the breaks of power law index.

2.2 Estimation of the stellar density

Paper I described a statistical method to derive the stellar density profile along each line-of-sight using spectroscopic survey data. Selection effects due to target selection and observational biases are taken into account in the density profile estimation. Here we briefly describe the method and suggest that the readers to refer to Paper I for more details about the approach.

We first assume that the photometric survey data is complete for the stars brighter than its limiting magnitude. Then we assume that, in a given line-of-sight which is defined by the Galactic coordinates ($l$, $b$) with solid angle $\Omega$, the probability of finding a star at distance $D$ in the spectroscopic survey data is approximately same as it is in the photometric data in a small region in the –magnitude plane, i.e.,

[TABLE]

where $c$ and $m$ stand for the index and apparent magnitude, respectively.

The two probability density functions can be expressed in terms of the stellar densities derived from the photometric data and the spectroscopic data, respectively. Note that the photometric data is presumed to be complete, while the spectroscopic data is not. Then, based on Eq. (1) the two densities, denoted as $\nu_{ph}$ and $\nu_{sp}$ respectively, are associated with a factor $S$, which represents the selection function of the spectroscopic data:

[TABLE]

Integrating over $c$ and $m$, we obtain the stellar density profile along the given line-of-sight:

[TABLE]

Paper I proved that Eq. (3) is held not only for the whole dataset, but also for a selected stellar population.

The selection function $S(c,m,l,b)$ can be derived from

[TABLE]

where $n_{sp}(c,m,l,b)$ and $n_{ph}(c,m,l,b)$ are the star counts of the spectroscopic and photometric data in the -magnitude plane, respectively.

To obtain $\nu_{ph}$ we need only estimate $\nu_{sp}$ from the observed spectroscopic data according to Eq. (3). We apply a kernel density estimation (KDE) to estimate $\nu_{sp}$ along a line-of-sight considering that 1) the number of spectroscopic stars is usually limited along a line-of-sight and 2) there is uncertainty of the distance. Then we can write

[TABLE]

where $p_{i}(D)$ is the probability density function of $R$ for the $i$th star. Note that in this form the derived $\nu_{ph}$ is a continuous function of $D$.

For a large survey like LAMOST, thousands of “plates" (used here to mean an observation with in this case 4000 fibers along a line-of-sight, even though the fiber positioners do not require production of a physical metal plate as in other surveys) have been observed. In principle, we can derive $\nu_{ph}$ for all plates and we can obtain $\nu_{ph}(D)$ at any arbitrarily selected distance $D$. However, at the position where no spectroscopic stars are sampled, the derived stellar density may suffer from large uncertainty. Hence, in practice, we only use the stellar density at the position in which a spectroscopic star is located, so that the precision of stellar density is preserved.

The advantage of this method is that no specific analytical form of the stellar density profile of the Milky Way is required. In other words, the stellar density estimation does not depend on the selection of a stellar density model of the Galaxy. For this reason, in this work, we are able to map the stellar halo density profile in a non-parametric manner.

Paper I did performance assessments and found that the typical error in the derived stellar density is less than $\sigma_{\nu}/\nu=0.4$. When the samples are very sparse along a line-of-sight, the derived stellar density profile may have some systematic bias. However, this systematic bias can be blurred when we solve for the Galactic structure by combining multiple lines-of-sight. But even in the worst case, when the systematic bias reaches $\sigma_{\nu}/\nu\sim 0.5$–$1$, it can only produce deviations of about 10% in the structural parameter estimates, according to the simulations presented in Paper I. Therefore, we infer that the accuracy of the stellar density estimates based on Paper I should be sufficient for our the study of the stellar halo.

We then estimate the stellar density at the positions of the selected K giant samples with [Fe/H]$<-1$, $M_{r,p1}<-0.5$, $K<14^{m}.3$, and $E(B-V)<0.25$. Because these samples are complete out to 35 kpc, the stellar density estimates should be reliable within this volume.

2.3 Number Density Map

We explore the stellar number density map as a function of ($R,Z$), where $R$ is radial distance on the Galactic plane ($R=\sqrt{X^{2}+Y^{2}}$) and $Z$ is the vertical distance to the Galactic mid-plane. The number density of K giants ($\nu$) at each spatial position is calculated using the method described in Paper I. In Figure 4, a tomographic map is constructed by calculating median for all values of the density of K giants ($median(\nu)$) at each R-Z pixel. In this map, there is no obvious exponential disc component near the Galactic plane due to the low metallicity selection. The isodensity contours look more like oblate ellipses rather than spherical within 30 kpc. Although the sample selection, distance estimates and interstellar extinction correction are different, our density map is quite consistent with the results of Paper I (see their figure 14). This indicates that our estimate of the halo density profile is quite robust.

In light of paper I, the stellar halo may be oblate in inner region and spherical in outer region. This can be directly tested by comparing the stellar density profile along $ln(R)$ with that along $ln(Z)$. Now consider that the profile follows a power-law:

[TABLE]

When the density profile is measured along $ln(Z)$, it becomes

[TABLE]

And when the density profile is along $ln(R)$, it reads

[TABLE]

If the stellar halo has a constant q, then it is easy to infer that the stellar density profiles along $ln(R)$ and $ln(Z)$ should be parallel with each other. If q increases with increasing radius, then the slope in equation 7 would not be n any more and hence the two density profiles should not be parallel any more.

Figure 5 shows the result of the test. We select the data located within $\theta>70^{\circ}(\theta=arctan(|Z|/R))$ in Figure 4 to represent the stellar density profile roughly along ln(Z) (displayed as red symbols in Figure 5) and select the data located at $\theta<20^{\circ}$ to stand for the profile along $ln(R)$ (black symbols in the figure). We fit them with single power-law and obtain that the best fit power index for the subset with $\theta<20^{\circ}$ is -4.78, while for the subset with $\theta>70^{\circ}$ is -2.89. It is obvious that the two density profiles are not parallel with each other, meaning that the flattening q must be variable with radius. More quantification of the variation is discusssed in next section.

Figure 6 shows the relative dispersion of the number density in each (R,Z) pixel. The relative dispersion is defined by (MAD of $\nu$)/median($\nu$). MAD of $\nu$ is median of absolute deviation defined as $median|\nu-median(\nu)|$. The relative dispersion is quite uniform in the range $10<r_{GC}<20$ kpc, with a mean relative dispersion of $0.42$. In the region where $r_{GC}>20$ kpc, the relative dispersions have larger scatter, but the mean value is still $0.44$. The dispersion in number density is caused by the intrinsic scatter of the number density due to the non-axisymetric substructures as well as the error induced by correction for selection effects and distance estimation. From the simulation in paper I, the uncertainty of the stellar density caused by selection effect correction can be as large as $\sigma_{\nu}/\nu=0.5-1$ in the worst case for a single plate, which will induce similar level random error on the iso-density surface. The relative dispersion in this work is not larger than the estimation in the paper I. This means that the density error caused by selection effect correction is well controlled.

Figure 7 shows that the distribution of the median [Fe/H] of our K giant samples is quite uniform within 25 kpc, while beyond 25 kpc the distribution of median [Fe/H] is more complicated because of overdensities.

3 Probing the number density profile using a nonparametric method

3.1 Method description

Previous studies explored the radial number density and flattening profile of the Milky Way’s stellar halo by fitting to different parameterized models. X15 shows the number density profile traced by SEGUE K giants can be described using the Einasto profile, a broken power law, or a single power law with constant or variable flattening comparably well, which indicates the degeneracy of the models. Thanks to the large sample and wide sky coverage of LAMOST K giants, it is possible to apply a direct approach to study the structure of the Galaxy with minimal assumptions. Unlike the popular model-dependent method, we devise a nonparametric method to probe the radial number density and flattening profiles of the Galactic stellar halo. We study the structure of the halo by fitting the outline of each iso-density contour to the density surface, directly.

Since many studies show the Galactic halo is oblate, we trace the halo density profile along the semi-major axis ($r$) of an ellipse. Figure 8 is a schematic cartoon to illustrate the geometry. $r$ is the semi-major axis of the ellipse defined as $r=\sqrt{R^{2}+(Z/q)^{2}}$, where flattening $q$ is the ratio of semi-minor to semi-major axis. $r_{GC}$ is the Galactocentric distance ($r_{GC}=\sqrt{R^{2}+Z^{2}}$). If the density profile is characterized by power law, then for each iso-density surface the stellar number density is given by $\nu=\nu_{0}(1/r)^{n}$. Therefore, $\ln(\nu)$ has linear relationship with $ln(r)$ i.e. $\ln(\nu)=\ln(\nu_{0})-n\;ln(r)$. For a given $\ln(\nu)$, the locus of R-Z pixels also satisfies an elliptical function.

The semi-major axis, $r$, and flattening, $q$, are determined by fitting the ellipse to each iso-density contour. We find it difficult to constrain $r$ and $q$ in Galacticcentric cylindrical coordinate, because $q$ is very sensitive to the slope of the iso-density surface, and the slope is too steep at low latitude. Therefore, we perform our analysis in Galactocentric polar coordinates:

[TABLE]

and

[TABLE]

Substituting Equation (9) (elliptic parameter equation) into Equation (10), we find

[TABLE]

We would like to solve for $r_{GC}$ in terms of the $\theta$ given in Figure 8 instead of $\eta$ used in Equation  (9). Since $\sin(\theta)=Z/r_{GC}$,

[TABLE]

Therefore, we can fit the isodensity contour with Equation (13):

[TABLE]

3.2 Fitting to the data

Following Paper I, we obtain $\nu$ at the position of each individual K giant. We divide the K-giant sample into 22 $\ln(\nu)$ bins. The binsize is $0.25$ in the range of $-16.375<\ln(\nu)\leq-11.5$, and increases to $0.5$ in the range of $-17.375\leq\ln(\nu)\leq-16.375$. The iso-density surfaces and the numbers of stars belonging to them are given in Table 2. The distribution of K giants on each iso-density surface and the fitting result is exhibited in Figure 9. The value of each $\ln(\nu)$ bin is given in the X-axis labels.

In general, we should fit the distribution of K giants for each iso-density surface in $(r_{GC},\sin\theta)$ space. In fact, we fit the median value of the distribution of the iso-density surface instead of fitting all data points to reduce the infuence of the outliers, many of which could result from halo substructure. The fitting is based on equation ( 13), with free parameters $r$ and $q$, using a Markov chain Monte Carlo (MCMC, Foreman-Mackey et al. 2012) method. The median values have larger scatter in the last few panels, because there are fewer K giants in each $\sin\theta$ bin. Also in the last three panels, the presence of the leading arm Sagittarius in the north strongly influences the last $\sin\theta$ bin; the clump of ourliers aroung $(r_{GC},\sin\theta)=(40,0.95)$, which is primarily composed of Sgr leading tail stars, represents more than half of stars in this bin. So the last bin in the three panels is not used the fitting process. The best-fit parameters $r$ and $q$ are listed in Table 2.

Note that the bins with $\ln(\nu)>-11.5$ are not used in the analysis. In these bins, the distribution of K giants can be fit by an ellipse or a straight line, because minor axis of the ellipse is small that it is comparable to the scatter of the distribution.

In order to demonstrate the performance of the model fitting, we take a case of $\ln(\nu)=-13.75$ as an example. The left panel of Figure 10 shows the deviations from the fit for the data in the $\ln(\nu)=-13.75$ bin ($r_{GC}-r_{GC,fit}$). The distribution of $r_{GC}-r_{GC,fit}$ is well described by Gaussian with $\sigma=2.8$, with the exception that the right wing is slightly over-populated. This may be due to the presence of density substructures in the halo. Furthermore, we transfer the fitting result back to $R-Z$ space in the right panel of Figure 10. As expected, most K giants in this $\ln(\nu)$ bin are roughly distributed in the shape of an ellipse. The outline of the fitting relation passing through the median values of each $\sin\theta$ bin of this $\ln(\nu)$ bin is shown in the $R-Z$ plane. In most cases, the best fit ellipes describe the iso-density surface well; some outliers appear at radii where the halo K giant sample is sparser.

Figure 11 shows the median values and outline of the elliptical relationship $r=\sqrt{R^{2}+(Z/q)^{2}}$ for each $\ln(\nu)$ bin, which shows that the ellipses are more oblate at small distances to the Galactic centre and become more spherical at larger distances.

The relationship between $r$ and $q$, as a function of $\ln(\nu)$, is shown in Figure 12. The upper panel shows the relationship between $r$ and $q$, and the lower panel shows the relationship between $ln(r)$ and $\ln(\nu)$. The error bars indicate the fitting errors derived from MCMC. We empirically fit the $r,q$ relation with the rational polynomial (black line): $q=\frac{p_{1}r^{2}+p_{2}r}{r^{2}+p_{3}r+p_{4}}$, with $p_{1}$=0.887, $p_{2}$=5.15, $p_{3}$=-7.418, and $p_{4}$=315.79. The red points show results from X15, for comparison. X15 defines $q$ as function of $r_{GC}$ (see Equation 5 of X15). We calculate $r$ and $q$ for each individual K giant of our sample with X15’s function and parameters and plot it in red on our figure. The value of $r$ is different for stars with same $R_{GC}$ but different latitude, so the relationship between $r$ and $q$ is a band in this diagram. Our fit is steeper than X15’s result within 25 kpc. Beyond 25 kpc, we find a more spherical halo with $q$ larger than 0.9. We didn’t fit isodensity contours beyond $r>35$ kpc, because our sample is incomplete at those distances.

The lower panel of Figure 12 shows the relationship between $\ln(\nu)$ and $\ln(r)$. We fit the relation with a power law. It can be well fitted with a single power law, with best fit power index of $n=-5.03^{+0.64}_{-0.64}$. We also tried to fit the $\ln(\nu)$ and $\ln(r)$ relation with broken power law, but the break radius was larger than 55 kpc, which is beyond the sampling range of our data.

The error bars in the upper and lower panels of Figure 12 indicate the fitting errors. There are other errors which are not included in the fitting errors: 1) The selection correction is a key point of this work. In the worst case, the error of $\ln(\nu)$ can be 0.5-1 for a single plate which introduces a similar fluctuation in the isodensity surface. 2) The random and systematic error in determining the distance is a source of error in the selection correction. 3) Unresolved substructure will introduce bias into the isodensity surface; substructure which does not belong to “the smooth component” will make the photometric catalogue denser, and that will in turn cause the number density function to be larger as well. 4) The remaining metal poor disc stars may weakly affect the density at low latitude. We tried to fit the median of the iso-density surface along $\sin\theta$ to reduce the effect of fluctuations in the isodensity surface and of substructure. In the future we will take steps to update of the selection correction to include errors in the estimation of distance.

3.3 Residuals and substructure

The residual map in Figure 13 is constructed by calculating residual, $(\nu-\nu_{fitting})/\nu_{fitting}$, for each star and then finding the median residual in each $R-Z$ bin. The bin size is still $2\times 2$ kpc as in Figures 4. $\nu_{fitting}$ is obtained from the power law density profile $\nu_{fitting}=(1/r)^{n}$. The value of $r$ and $q$ for each star is obtained by solving the set of equations: $r=\sqrt{R^{2}+(Z/q)^{2}}$ and $q=\frac{p_{1}r^{2}+p_{2}r}{r^{2}+p_{3}r+p_{4}}$. The residual map within 10 kpc can be ignored, because there is no estimate of $q$ and $r$ within 10 kpc, and the extrapolation of $q$ in this range is not reliable. Although we show the residual map up to 60 kpc, we should keep in mind that the map is only complete to 35 kpc for stars with $-4<M_{r}<-0.5$. Also, $q$ is extrapolated as a constant beyond $r=35$ kpc. We assume that the power law index is unchanged for $r>35$ kpc.

Figure 13 shows that the number density distribution of K giants is well fit within 25 kpc. The median of distribution of residuals is less than 20% in the region where $10<r_{GC}<25$ kpc. Beyond $r_{GC}=25$ kpc , the number density residual shows some prominent features. The significant feature in the location of $Z>30$ kpc and $R<25$ kpc may be associated with the leading arm of the Sagittarius tidal stream in the northern Galactic hemisphere. The scattered overdensity in the southern sky may associated with the trailing arm of Sgr. We show this in detail in section 3.4 and Figure 15. In addition to the known features, we detect one prominent overdensity around $(R,Z)=(30,15)$ kpc. We label the overdensity with black rectangle in Figure 13.

3.4 The overdensities related to the Sagittarius stream

The prominent overdensities in north and south Galactic hemisphere shown in the residual map may be associated with the leading and trailing arms of Sgr. In this section, we validate the correlation between the identified overdensities in the residual map and the Sgr tidal stream by selecting Sgr candidates from our K giant samples using criteria from Belokurov et al. (2014), and then overplotting the Sgr candidates on the residual map.

Figure 14 shows the LAMOST K giants within $\pm 12^{\circ}$ of the equator of the Sgr coordiante system (within $\pm 12^{\circ}$ of the Sgr dwarf orbital plane). The figure shows the distribution of $V_{gsr}$ 111$V_{gsr}=rv+9.58*\cos(b)*\cos(l)+(220+10.52)*\cos(b)*\sin(l)+7.01*\sin(b)$ (Xu et al. 2015) and distance modulus as a function of angle along the tidal stream respectively. The coordinate of ($\Lambda_{\odot}$, $B_{\odot}$) is a heliocentric coordinate system whose equator is aligned with the Sgr trailing tail. The longitude $\Lambda_{\odot}$ increases in the direction of Sgr motion, and the latitude $B_{\odot}$ points to the North Galactic pole, which is the same as defined in Belokurov et al. (2014).

The paths of the Sgr leading and trailing arms along the distance modulus, as described in Belokurov et al. (2014), are shown in the upper panel of Figure 14 as grey curves. The north-leading-arm path starts from ($\Lambda_{\odot}$, (m-M))=$(37.1^{\circ},17.83)$ and ends at $(117.2^{\circ},17.3)$. Compared with the leading arm path of Belokurov et al. (2014), the leading arm traced by LAMOST K giants with $V_{gsr}\sim-100$ kms*-1* (blue points) extends to ($\Lambda_{\odot}$, (m-M))=$(160^{\circ},16)$, which is consistent with that traced by the stream’s sub-giant branch (SGB) stars reported by Belokurov et al. (2006). The path of the Sagittarius trailing arm in the north, defined by Belokurov et al. (2014), starts from ($\Lambda_{\odot}$, (m-M))=$(138.9^{\circ},19.11)$ and ends at $(191.7^{\circ},19.49)$, the distance of trailing arm in the north is further than 60 kpc. The path of Sgr trailing arm in the south, as traced by Belokurov et al. (2014) at $220^{\circ}<\Lambda_{\odot}<300^{\circ}$, is observed in LAMOST K giants in the region $220^{\circ}<\Lambda_{\odot}<250^{\circ}$ and $16.5<(m-M)<18.5$. This identification is confirmed by the obvious clump at ($\Lambda_{\odot},V_{gsr}$)=$(245^{\circ},-150$ kms${}^{-1})$ in lower panel of Figure 14.

Note that the Sgr trailing arm in the north is farther than 60 kpc, which is beyond our detection range, so it won’t be considered seriously in this work.

We select Sgr northern leading arm and southern trailing arm candidates along the path of $V_{gsr}$ and distance modulus as a function of $\Lambda_{\odot}$ as defined by Belokurov et al. (2014). The northern leading arm candidates are selected using $V_{gsr}$ within 50 kms*-1* and $(m-M)$ within 1 mag of the Sgr leading arm paths in Figure 14. The southern trailing arm candidates are also selected along the corresponding path depicted in Figure 14. Because the clump of LAMOST K giants at ($\Lambda_{\odot},V_{gsr}$)=$(245^{\circ},-150$ kms${}^{-1})$ is slightly shifted from the grey curve of the Sgr southern trailing arm in the lower panel of Figure 14, the distance between $V_{gsr}$ of LAMOST K giants and the path centre is required to be greater than -70 kms*-1* and smaller than 30 kms*-1* within the longitude range $230^{\circ}<\Lambda_{\odot}<255^{\circ}$.

The selected Sgr northern leading arm candidates are overplotted in the left panel of Figure 15, while the southern trailing arm candidates are labelled in the right panel. The background is part of residual map of Figure 13. The Sgr candidates lie exactly in the regions where the observational number density is higher than the best fit model.

We also tried to select candidates associated with the overdensity $(R,Z)=(30,15)$ kpc, within the distance range $20<R<35$ kpc, $5<Z<25$ kpc, as depicted by the rectangle in Figure 13. The candidates are selected with $(\nu-\nu_{fitting})/\nu_{fitting}>0.5$. To determine whether the overdensity is associated with Sgr, we plot the candidates in the Sgr coordinate system in Figure 16. It is seen that there are two separated sub-groups of candidates. The first one, with $\Lambda_{\odot}<80^{\circ}$ and $B_{\odot}<-20^{\circ}$, is far from the equator of the Sgr stream coordinate system. The second one, with $110^{\circ}<\Lambda_{\odot}<180^{\circ}$ and $-20^{\circ}<B_{\odot}<40^{\circ}$, passes through the Sgr equator. Figure 17 shows the distribution of distance modulus and $V_{gsr}$ for candidates associated with the overdensity at $(R,Z)=(30,15)$ kpc as a function of distance along Sgr longitude. In the upper panel of Figure 17, the overdensity candidates are distributed along the extension of the line tracing the Sgr northern leading arm. In the lower panel of Figure 17, the $V_{gsr}$ of part of the candidates overlap the grey curve of Sgr northern leading arm. However, many other the candidates are surprisingly distributed closer to the grey curve of north-trialing-arm. The $V_{gsr}$ distribution of the overdensity at $(R,Z)=(30,15)$ kpc is clearer in Figure 24, showing the metal poor subsample, in APPENDIX A. The distance to the overdensity at $(R,Z)=(30,15)$ kpc is about 25 kpc, while the northern trialing arm is further than 60 kpc; the distances are not consistent with each other. So the overdensity cannot be a piece of northern trialing arm. It might be part of Sgr in other wrap, or, it might be a new subsucture. Because the origin of the overdensity is beyond the scope of this paper, we will pursue this subject in the further work.

4 Discussion and Conclusion

Carollo et al.(2007) found that the stellar halo can be separated into the inner and outer halos based on density profile, kinematic, and metallicity. The boundary between the inner and outer halos is about 15 kpc. The inner halo is more oblate than the outer halo. In addition, the inner and outer halos may have different formation histories. The accretion of the inner halo is thought to have happened 10 Gyr ago, and part of inner halo may have formed in situ. The accretion of the outer halo happened much more recently than that of inner halo; most of the outer halo is thought to be composed of accretion debris (Abadi et al. 2006, Bland-Hawthorn et al. 2016).

In later studies, the density profile of the halo has been mapped with different tracers (Deason et al.2011, Sesar et al. 2010). They found the inner halo has a shallower density profile than outer halo. The simulation of Deason et al. (2013) indicates that if there is a break in the power law index, then the strength of the break depends on the accretion time and the size of the accretion debris; the location of the break depends on the average apocentre of the accreted halo component.

We presnt a direct method to detect the structure of the Galactic halo. The method relies on relatively complete sky coverage (including areas of both high Galactic latitude and low Galactic latitude) and a large sample of spectroscopic tracers. In our work, we started from the defination of density distribution and fit the contours of the iso-density surface. We did not predetermine whether the flattening would change with Galactic radius, whether there is break on the power law index, or where is the possible break might occur. We determined the density profile and flattening by fitting the iso-density surface in the $(r_{GC},\sin\theta)$ plane. We did not eliminate whole areas of the sky occupied by stellar streams; we fit medians of the iso-density surface to avoid the influence of outliers which could be related to substructures. Because we used the median value of ($r_{GC}$,$\sin\theta$) when we fit the iso-density surface in the $r_{GC}-\sin\theta$ plane, the affect of the tidal streams was reduced, though not entirely removed.

We obtained the flattening and power law index independently, breaking the degeneracy between flattening and power law index. We find that the flattening $q$ changes with Galactic radius, which is consistent with result of the SPL model of X15; $q$ changes rapidly from $r=15$ kpc to 25 kpc. We find no evidence of a break of power law between 15 kpc and 35 kpc. Our data does not extend as far as X15’s SDSS K giants, but our sample is much larger. Because of the large sample and the fact that we adopted a non-parametric method, the constraint on the flattening is more independent and reliable in the distance range ($r_{GC}<$ 35 kpc).

To show the difference between our result and previous results more clearly, we compare them further. From table 1, some of previous works support SPL (Robin, et al. 2008, Juric et al. 2008, Bell et al. 2008). Meanwhile, some of others support BPL (Watkins et al. 2009, Sesar et al. 2011, Deason et al. 2011, Pila-Diez et al. 2015). The broken position changes from 19 kpc to 28 kpc, but the power law index of inner and outer halo are quite similar with each other. The density can be fitted either by SPL with varying q or BPL in the work of X15 and Das et al. (2016). Other recent works support triaxial halo model (Pila-Diez et al. 2015, Iorio et al. 2017).

The typical results of previous works are compared with our result in Figure 18. Because our analysis is under the assumption that the halo is axisymmetric. It’s hard to compare with results of fitting by triaxial halo models. From the comparison in figure 18, our result is much steeper than the other results fitted with either BPL (Deason et al. 2011, blue lines) and SPL (Juric et al. 2008, brown lines) within the range of 25 kpc. The density profile modelled by BPL with segmental flattening (Pila-Diez et al. 2015, green lines) is steeper than Deason et al. (2011) and Juric et al. (2008) but flatter than our result. The result from X15 (red points and lines) and Das et al. (2016, light blue points and lines) fitting by SPL with varying q are similar to this work.

Compared to previous works, we adopted different tracers and apply different approach, which may lead to systematic difference in the slope of the power-law model. However, the assumption of a constant q with value smaller than 1 and the selection effect of excluding the low Galactic halo stars may be the main reasons that the slope is much flatter in many of the previous works.

Figure 5 and the discussion in section 2.3 provide two exterme cases of the stellar density profile, one is along $ln(R)$ and the other along $ln(Z)$. It shows that when q is variable the two profiles have different slopes. If an inproper assumption that q is a constant is made, then the observed stellar density profile should be somewhere in between the two extreme profiles shown in Figure 5. Consequently, the best fit slope of the data should be substantially flatter than the slope derived from the profile along $ln(R)$ (the black symbols). Note that the parameter r in our model represent for the length of the semi-major axis, which is equivalent to R when q<1. Therefore, the derived slope with the assumption of a constant q smaller than 1 would lead to a flatter slope of the power-law than our result. Moreover, almost all previous works excluded the low Galactic latitude halo stars to aviod disc contamination. This leads to a systematics bias toward the density profile along $ln(Z)$, which is flatter than that along $ln(R)$ (see Figure 5). Thus, again, the selection effects prefers to a flatter power-law. Combining the two situations together, it is not surprise that many previous works using constant q derived much flatter slope, whileas X15 and Das et al.(2016), who adopted a variable q, obtain steeper slope similar to this work.

Our work measure the iso-density of the halo directly and model independent. The most crucial ingredients relies on the selection correction. We plan to revisit this technique in the future, with an updated selection correction and combining Pan-STARRS and UCAC4 photometry, so that the completeness will be extended to a fainter magnitude limit. Also, in LAMOST DR4, the $\alpha$ abundance will be provided so that we can calculate better distance estimates for the giant stars. The proposed improvements will allow us to better study the structure of the halo, particularly for $r_{GC}>35$ kpc.

Acknowledgements

This work is supported by the Strategic Priority Research Program “The Emergence of Cosmological Structures” of the Chinese Academy of Sciences, Grant No. XDB09000000 and the National Key Basic Research Program of China 2014CB845700. CL acknowledges the National Natural Science Foundation of China grant No. 11373032 and 11333003. HJN acknowledges support from the US NSF grant AST 16-15688. 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. We appreciate the anonymous referee for his/her the very helpful comments and suggestions.

Appendix A Checking our result with metal-poor subsample

In the main body of this paper, we use LAMOST K giants with [Fe/H]$<-1$ to trace the density profile of the halo. This criterion, and the exclusion of the $\sin\theta=0.1$ bin when we fit the density profile of halo in Figure  9, removes most of the disc stars. But our sample might still include contamination from a metal-poor thick disc component (X15). Although we don’t expect it to inflence our result much, we check our result with more metal poor sample which is more likely a pure halo sample in this Appendix.

We limit the subsample of K giants to [Fe/H]$<-1.5$ in this analysis; the lower metallicity cut results in a sample of 2484 K giants, which is about half of the previously analyzed sample. We fit each iso-density surface for the free parameters $r$ and $q$ with the same method used in the larger sample. Figure  19 shows the flattening variation and density profile obtained from the metal-poor subsample. The upper panel shows how the flattening changes with the radius. The axis ratio increases from 15 kpc to 30 kpc, but the uptrend is a little more gentle than that of the larger sample. The flattening q is about 0.62, 0.76, and 0.98 for $r=15,20,$ and $30$ kpc, respectively. The power law index $n=-5.19^{0.85}_{0.85}$ is consistent with that of larger sample.

Because the subsample of [Fe/H]$<-1.5$ is only half of total sample, it can only trace the isodensity surface within $r<30$ kpc. Futher than that, there are not enough stars in each $\ln(\nu)$ bin. Besides, the scatter in the distribution of axial ratio and power law index observed in Figure 19 is larger than that of Figure 12.

We follow the same method used in the larger sample to calculate the residual map of the stellar density. Figure  20 is the residual map for the metal-poor subsample. The results are consistent with that of the larger sample. The median of residual is less than 27% in the range $r_{GC}<25$ kpc. Beyond 25 kpc, the overdensities of the Sgr northern leading arm, Sgr southern trialing arm, and the overdensity at $(R,Z)=(30,15)$ kpc are still significant in the map.

To study the relationship between the overdensities and the largest substructure, Sagiattarius, we replot Figure 14 with the metal poor subsample in Figure  21. The patterns are quite similar to those of Figure  14. In Figure 22, the candidates in the metal-poor subsample also overlap the overdensities in the residual map, which are consistent with those in the larger sample.

We also select the candidates associated with the overdensity at (R,Z)=(30,15) kpc. The candidate selection criteria are the same as that of the larger sample and the selected candidates are shown in Figures 23 and  24. The metal-poor subsample also shows two groups in Figure 23. The $V_{gsr}$ distribution of metal-poor subsample candidates, in the lower panel of Figure 24, shows quite clearly that some of the candidates follow the path of Sgr northern leading arm, and some of the candidates follow the path of Sgr northern trailing arm, in agreement with the findings from the larger sample.