The Profile of the Galactic Halo from Pan-STARRS1 3π RR Lyrae

Stellar density models can take various functional forms. We first describe what the stellar density models we use for evaluating the RRab sample have in common.

In what follows we will assume that our models are characterized by a set of parameters denoted as ${\boldsymbol{\theta }}$ and that the density ρ_RRL is ellipsoidal, allowing for a halo flattening q along the Z direction. Oblate density distributions have q < 1, spherical have q = 1, and prolate have q > 1.

The density is a function of right-handed Cartesian coordinates (X, Y, Z), which we evaluate through the Galactic longitude, Galactic latitude, and heliocentric distance (l, b, D), so its dimension is kpc⁻³:

$$\begin{eqnarray}X & = & {R}_{\odot }-D\cos l\cos b\\ Y & = & -D\sin l\cos b\\ Z & = & D\sin b.\end{eqnarray} \tag{ 1 }$$

This reference frame is centered at the Galactic center. The Galactic disk is in the (X, Y) plane, with the X-axis pointing to the Sun and the Z-axis to the north Galactic pole. R_⊙ denotes the distance of the Galactic center from the Sun, in this work assumed to be 8 kpc, and the main results of our work should not change for other values of R_⊙ within the assumed observational uncertainties.

The vertical position of the Sun with respect to the Galactic disk is uncertain, but it is estimated to be smaller than 50 pc (Iorio et al. 2017; Karim & Mamajek 2017) and thus negligible for the purpose of this work.

Caution must be taken when comparing our work to others: some papers use a left-handed system instead (e.g., Iorio et al. 2017), where the Y-axis is flipped with respect to our definition.

With Equation (1), the galactocentric distance R_gc is then defined as ${R}_{\mathrm{gc}}=\sqrt{{X}^{2}+{Y}^{2}+{Z}^{2}}$, and the flattening-corrected radius is defined as ${r}_{q}=\sqrt{{X}^{2}+{Y}^{2}+{(Z/q)}^{2}}$, where q gives the halo flattening along the Z direction as a minor-to-major-axis ratio. This describes an oblate stellar halo that is stratified on concentric ellipsoids, where X, Y, Z are the ellipsoid principal axes.

Following a number of previous studies, we presume that the overall radial density profile of the halo can be described by a power law or an Einasto profile, with the density stratified on concentric ellipsoidal surfaces of constant r_q in all cases.

3.1.1. Power-law Profile

A simple power-law halo model ρ_halo is widely used (e.g., Sesar et al. 2013b) to describe the distribution of the halo stars:

$$\begin{eqnarray}&&{\rho }_{\mathrm{halo}}(X,Y,Z)={\rho }_{\odot \mathrm{RRL}}{({R}_{\odot }/{r}_{q})}^{n}.\end{eqnarray} \tag{ 2 }$$

For a power-law profile, the shape of the density profile is described by the parameter n. Larger values of n indicate a steeper profile.

The free parameters are ${\boldsymbol{\theta }}=(n,q)$. Here ${\rho }_{\odot \mathrm{RRL}}$ is the number density of RR Lyrae at the position of the Sun, R_⊙ is the distance of the Sun from the Galactic center, and r_q is the flattening-corrected radius. As we are not interested in absolute numbers, we are not fitting for ρ_⊙RRL.

Others presume a broken power law (BPL; e.g., Xue et al. 2015), where inner and outer power-law indices are defined. The change in the power-law index then occurs by a step function at the break radius. As our sample starts at a galactocentric radius of 20 kpc, and the break radius is found to be around or below 20 kpc (e.g., Xue et al. 2015), we cannot compare to the results by Xue et al. (2015). However, in order to compare to the findings by Deason et al. (2014), who find a BPL with three ranges of subsequently steepening slope, where one of the breaks is occurring within the distance range present in our sample, we fit a BPL:

$$\begin{eqnarray}&&{\rho }_{\mathrm{halo}}(X,Y,Z)=\left\{\begin{array}{l}{\rho }_{\odot \mathrm{RRL}}{({R}_{\odot }/{r}_{q})}^{{n}_{\mathrm{inner}}},{\rm{if}}\,{\rm{}}{r}_{q}\leqslant {r}_{\mathrm{break}}\\ {\rho }_{\odot \mathrm{RRL}}{r}_{\mathrm{break}}^{{n}_{\mathrm{outer}}-{n}_{\mathrm{inner}}}{({R}_{\odot }/{r}_{q})}^{{n}_{\mathrm{outer}}},{\rm{else.}}\,{\rm{}}\end{array}\right.\end{eqnarray} \tag{ 3 }$$

3.1.2. Einasto Profile

The Einasto profile (Einasto 1965; Einasto & Haud 1989) is the 3D analog to the Sérsic profile (Sérsic 1963) for surface brightnesses and has been used to describe the halo density distribution (Merritt et al. 2006; Deason et al. 2011; Sesar et al. 2011; Xue et al. 2015; Iorio et al. 2017) and dark matter halos (Merritt et al. 2006; Navarro et al. 2010). It is given by

$$\begin{eqnarray}&&\gamma (r)\equiv -\displaystyle \frac{d\mathrm{ln}\rho (r)}{d\mathrm{ln}r}\propto {r}^{\alpha },\end{eqnarray} \tag{ 4 }$$

where the steepness of the Einasto profile, α, changes continuously as a function of the effective radius r_eff,

$$\begin{eqnarray}&&\alpha =-\displaystyle \frac{{d}_{n}}{n}{\left(\displaystyle \frac{{r}_{q}}{{r}_{\mathrm{eff}}}\right)}^{1/n}.\end{eqnarray} \tag{ 5 }$$

This can be rearranged to

$$\begin{eqnarray}&&{\rho }_{\mathrm{halo}}({r}_{q})\equiv {\rho }_{0}\exp \{-{d}_{n}[{({r}_{q}/{r}_{\mathrm{eff}})}^{1/n}-1]\},\end{eqnarray} \tag{ 6 }$$

where ρ₀ is the (here irrelevant) normalization, r_eff is the effective radius, and n is the concentration index. The parameter d_n is a function of n, where for n ≥ 0.5 a good approximation is given by ${d}_{n}\approx 3n-1/3+0.0079/n$ (Graham et al. 2006).

For an Einasto profile, the shape of the density profile is described by the parameter n. This profile allows for a nonconstant fall-off without the need for imposing a discontinuous break radius: density distributions with steeper inner profiles and shallower outer profiles are generated by large values of n, whereas small values of n account for a shallower inner and steeper outer profile. The parameter r_eff describes the radius of the inner core of the profile.

The free parameters of an Einasto profile with a constant flattening q are ${\boldsymbol{\theta }}=({r}_{\mathrm{eff}},n,q)$.

3.1.3. Profiles with Varying Flattening

The models described so far assume a constant flattening q. However, Preston et al. (1991) found evidence for a decrease in the flattening with increasing radius. Carollo et al. (2007, 2010) find evidence that at least the innermost part of the halo is quite flattened.

We thus increase the complexity of the model by allowing for a nonconstant flattening of the halo, parameterized by the galactocentric radius. To describe such radial variations of the stellar halo's flattening, we consider the functional form for $q({R}_{\mathrm{gc}})$ as

$$\begin{eqnarray}&&q({R}_{\mathrm{gc}})={q}_{\infty }-({q}_{\infty }-{q}_{0})\exp \left(1-\displaystyle \frac{\sqrt{{R}_{\mathrm{gc}}^{2}+{r}_{0}^{2}}}{{r}_{0}}\right),\end{eqnarray} \tag{ 7 }$$

with q₀ being the flattening at the center, q_∞ being the flattening at large galactocentric radii, and r₀ being the exponential scale radius over which the change of flattening occurs.

Thus, the flattening q now varies from q₀ at the center to the asymptotic value q_∞ at large radii, and the variation is tuned by the exponential scale length r₀.

All other equations to describe the radial profile given above apply from the previously described Einasto and power-law profile, replacing only q with q(R_gc), and thus replacing the fitting parameter q with three fitting parameters q₀, q_∞, and r₀.

In general, a selection function describes the fraction of stars that are targeted, as a function of, e.g., position, distance, or magnitude.

We introduce the selection function for two reasons: to correct for the noncomplete volume sampling naturally occurring during a survey, and to remove known overdensities to build a "clean" sample of RRab, eliminating all the stars belonging to the substructures from our original catalog. Both cosmological models and observations imply that a good portion of halo stars, at least beyond 20 kpc, are in substructures. Especially the prominent ones, such as the Sagittarius stream and the Virgo overdensity, can and will affect the fits of smooth models, as pointed out also by Deason et al. (2011).

After we remove those known substructures, it is, of course, still possible that there are previously unknown substructures, as well as that the smooth component of the halo is also structured but at a level that is below our resolution.

Our selection function ${ \mathcal S }(l,b,D)$ is binary $[0,1]$ so that ${ \mathcal S }$ is always equal to 1 except for the points (l, b, D) that are excluded. The predicted density of stars is then simply the product of the underlying density distribution with the selection function, suggesting that one constrains this underlying density by forward modeling of the observations.

The RRab candidates from Sesar et al. (2017) were selected uniformly from the set of objects in the PS1 3π survey in the area and apparent magnitude range available for this survey. The selection completeness and purity are uniform over a wide range of apparent magnitude up to a flux-averaged r-band magnitude of 20 mag (Sesar et al. 2017), which is described later on in Equation (12).

Starting from the 44,403 RRab stars in the sample of Sesar et al. (2017), we exclude known overdensities in (l, b, D). Among the largest overdensities are the Sagittarius stream, dwarf galaxies such as Draco dSph, and globular clusters. A complete list can be found in Table 1. Also, we cut out sources too close to the Galactic plane ($| b| \lt 10^\circ$), or too close to the Galactic center (R_gc ≤ 20 kpc), as we want to avoid regions with many overdensities such as streams as mostly found within 20 kpc, we want to excise the Galactic bulge, and additionally the RRab sample is relatively sparse toward the Galactic disk.

Table 1.  Removed Overdensities within R_gc > 20 kpc

Name	${R}_{\mathrm{gc}}$ a	Remove	Remove	Remove	Remove	Remove	Remove	Removed
	(center)	l min (deg)	l max (deg)	b min (deg)	b max (deg)	D min (kpc)	D max (kpc)	Sources
Bootes III dSph	37.87	32	34.1	74.5	75.4	45.9	46.5	3
Sextans dSph	45.14	242	245	41	44	60	120	99
NGC 292 Bootes I dSph	45.88	357	359	69	70	55	70	4
UMa 1 dSph	60.17	150	160	54	54.6	90	120	4
Draco dSph	80.70	84	87	33.5	35.5	65	100	191
UMi dSph	48.23	100	110	40	50	60	80	53
NGC 7089 M2	25.08	53	53.5	−36	−35.5	11	12.5	4
NGC 6626 M28	71.90	7.8	8	−6	−5.3	5.3	5.8	3b
Pal 3	45.25	240	240.2	41.8	42	80	100	3
Laevens 3	45.56	63.58	63.602	−21.2	−21.13	55	62	2
NGC 2419	56.06	178	183	24	26	76	84	8b
NGC 6293	82.37	357	359	7	9	9	10	16b
NGC 6402 M14	93.26	21	21.8	14.5	15.2	8	10	6
NGC 6171 M107	109.44	2.8	3.8	22.1	23.7	5.5	8	7b
Pisces Overdensity	51.67	87.3	87.4	−58.2	−57.9	79	82	1
RR10c	25.49	186.37	186.38	51.5	51.6	41.2	41.3	1b

Notes. Overdensities are grouped by dwarf galaxies, globular clusters, and others (Pisces overdensity, and the single RRab star RR10). In each group, they are ordered by R_gc.

^aThe center of the removed overdensity. ^bThese sources are also removed by other cuts. ^cThis RRab is a member of the Orphan stream (Sesar et al. 2013a).

Download table as:  ASCIITypeset image

From the 33,378 sources we exclude in total, 6575 are within ±10° of the Galactic plane, 26,951 are within 20 kpc of the Galactic center, 5960 are in the Sgr stream, and 578 are in other overdensities as listed in Table 1; as those regions partially overlap, the numbers stated here would add up to 35,484.

The selection function ${ \mathcal S }(l,b,D)$ is thus composed of

$$\begin{eqnarray}&&{ \mathcal S }(l,b,D)={{ \mathcal S }}_{\mathrm{RRL}}(l,b,D)\times {{ \mathcal S }}_{\mathrm{area}}(l,b,D),\end{eqnarray} \tag{ 8 }$$

where ${{ \mathcal S }}_{\mathrm{RRL}}(l,b,D)$ describes the selection cuts of the sample introduced by the survey and Sesar et al. (2017), itself leading to the 44,403 RRab stars, and ${{ \mathcal S }}_{\mathrm{area}}(l,b,D)$ describes area cuts to exclude overdensities.

The area and depth of the PS1 sample of RRab lead to

$$\begin{eqnarray}{{ \mathcal S }}_{\mathrm{RRL}}(l,b,D)=\left\{\begin{array}{cl}1, & {\rm{if}}\,{\rm{}}\delta \gt -30^\circ \,\mathrm{and}\ {D}_{\min }\lt D\lt {D}_{\max }\\ 0, & {\rm{else.}}\,{\rm{}}\end{array}\right.\end{eqnarray} \tag{ 9 }$$

The spatial cuts to geometrically excise bulge and thick-disk stars beyond a galactocentric distance of 20 kpc are

$$\begin{eqnarray}{{ \mathcal S }}_{\mathrm{bulge},\mathrm{disc}}(l,b,D)=\left\{\begin{array}{cl}1, & {\rm{if}}\,{\rm{}}| b| \geqslant 10^\circ \,\mathrm{and}\ {R}_{\mathrm{gc}}\geqslant 20\,\mathrm{kpc}\\ 0, & {\rm{else.}}\end{array}\right.\end{eqnarray} \tag{ 10 }$$

The spatial cuts to geometrically excise the Sagittarius (Sgr) stream are based on our previous work describing the Sgr stream's 3D geometry as traced by PS1 RRab stars (Hernitschek et al. 2017). To each star in the sample, we can assign a probability that it is associated with the Sgr stream, p_sgr (Hernitschek et al. 2017, see Equation (11) therein). We excise sources with p_sgr > 0.2 as members of the Sgr stream, leading to a selection function of

$$\begin{eqnarray}{{ \mathcal S }}_{\mathrm{sgr}}(l,b,D)=\left\{\begin{array}{cl}1, & {\rm{if}}\,{\rm{}}{p}_{\mathrm{sgr}}(l,b,D)\lt 0.2\\ 0, & {\rm{else.}}\,{\rm{}}\end{array}\right.\end{eqnarray} \tag{ 11 }$$

Additional spatial cuts are used to remove all stars in the boxes listed in Table 1 in the Appendix, in order to excise known overdensities. This results in ${{ \mathcal S }}_{\mathrm{other}\mathrm{overdensities}}(l,b,D)$.

Taking into account that the RRab sample is not complete, with the completeness varying with magnitude, another term for the selection function needs to be introduced.

Sesar et al. (2017) find that the RRab selection function is approximately constant at ∼90% for a flux-averaged r-band magnitude r_F ≲ 20 mag, after which it steeply drops to zero at r_F ∼ 21.5 mag. Writing r_F as r_F(D), the selection function characterizing the distance-dependent completeness is

$$\begin{eqnarray}&&{{ \mathcal S }}_{c}({r}_{{\rm{F}}})=L-\displaystyle \frac{L}{1+\exp (-k({r}_{{\rm{F}}}-{x}_{0}))},\end{eqnarray} \tag{ 12 }$$

with (Sesar et al. 2017)

$$\begin{eqnarray}&&L=0.91\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&k=4.0\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{x}_{0}=20.6\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{r}_{{\rm{F}}}=2.05\mathrm{log}(D)+11.\end{eqnarray} \tag{ 16 }$$

In addition, we estimated the distance-dependent purity to supplement the overall sample purity that was given as 90% by Sesar et al. (2017). Using the RRab sample within SDSS S82, as done by Sesar et al. (2017) to estimate the distance-dependent completeness of our RRab sample, we find the purity staying stable at a level of 98% to 95% over a range from 15 to more than 20 mag in the r band. In contrast, over the same magnitude range, the completeness drops from 91% to 80%. The faintest RRab in S82 (which we use as the validation set; see Sesar et al. 2017) is found at r_F = 20.58 mag, and there are in total only two sources in this faintest 0.5 mag bin. The 10 faintest RRab stars within S82 span a distance range from 85 to 102 kpc. This means that for sources fainter than 20.5 mag, the purity cannot be estimated in this way. For sources beyond D = 90 kpc, we adopted a purity of 94%. There is no SDSS source within S82 that was not picked up by PS1. The different distance dependency of purity and completeness reflects that it is easier to lose objects (i.e., not to classify them as RRab stars) than to get spurious sources into the catalog of PS1 RRab stars, given the rigorous definition adopted to consider a star as RRab (Sesar et al. 2017). Although the effect of the purity is negligible, as the effect of a dropping completeness at large distances dominates, and we cannot determine the purity beyond D = 90 kpc, we included it as part of the selection function, ${{ \mathcal S }}_{p}(D)$.

We end up with a selection function

$$\begin{eqnarray}{ \mathcal S }(l,b,D) & = & {{ \mathcal S }}_{\mathrm{RRL}}(l,b,D)\times {{ \mathcal S }}_{c}(D)\times {{ \mathcal S }}_{p}(D)\times {{ \mathcal S }}_{\mathrm{area}}(l,b,D)\\ & = & {{ \mathcal S }}_{\mathrm{RRL}}(l,b,D)\times {{ \mathcal S }}_{c}(D)\\ & & \times {{ \mathcal S }}_{\mathrm{bulge},\mathrm{disc}}(l,b,D)\\ & & \times {{ \mathcal S }}_{\mathrm{Sgr}}(l,b,D)\times {{ \mathcal S }}_{\mathrm{other}\mathrm{overdensities}}(l,b,D).\end{eqnarray} \tag{ 17 }$$

The overdensities listed in Table 1 are chosen in the following way: based on a list of dwarf galaxies within 3 Mpc by McConnachie (2012), its update from 2014,⁷ and a list of currently known halo streams by Grillmair & Carlin (2016), we select overdensities that could show up in a survey that covers the position and distance cuts of PS1 3π. We check each overdensity to see whether it appears in the RRab sample, and if so, we cut it out by defining a selection box in (l, b, D). We end up with the cuts described in Table 1.

After excising stars using ${{ \mathcal S }}_{\mathrm{area}}(l,b,D)$, the sample reduces to 11,025 RRab stars, which we call the "cleaned sample." The original sample and the cleaned sample are shown in Figure 2.

Out of these sources, 679 lie beyond a galactocentric distance of 80 kpc, and 101 beyond 100 kpc, in contrast to 1093 sources beyond 80 kpc, and 238 beyond 100 kpc in the original sample.

We now incorporate this selection function in fitting a parameterized model for the stellar density of the halo.

With the models ${\rho }_{\mathrm{RRL}}({ \mathcal D }| {\boldsymbol{\theta }})$ and the selection function ${ \mathcal S }$ at hand, we can directly calculate the likelihood of the data ${ \mathcal D }$ given the model ρ_RRL, the fitting parameters ${\boldsymbol{\theta }}$, and the selection function ${ \mathcal S }$ following Bovy et al. (2012).

The normalized unmarginalized log likelihood for the ith star with the observables ${{ \mathcal D }}_{i}$ is then

$$\begin{eqnarray}&&\mathrm{ln}p({{ \mathcal D }}_{i}| {\boldsymbol{\theta }})=\displaystyle \frac{{\rho }_{\mathrm{RRL}}({{ \mathcal D }}_{i}| {\boldsymbol{\theta }})| {\boldsymbol{J}}| { \mathcal S }({l}_{i},{b}_{i},{D}_{i})}{\int \int \int {\rho }_{\mathrm{RRL}}(l,b,D| {\boldsymbol{\theta }})| {\boldsymbol{J}}| { \mathcal S }(l,b,D){dldbdD}},\end{eqnarray} \tag{ 18 }$$

where the normalization integral is over the observed volume. The Jacobian term $| {\boldsymbol{J}}| ={D}^{2}\cos b$ reflects the transformation from (X, Y, Z) to (l, b, D) coordinates.

We evaluate the logarithmic posterior probability of the parameters ${\boldsymbol{\theta }}$ of the halo model, given the full data ${ \mathcal D }$ and a prior $p({\boldsymbol{\theta }})$, $\mathrm{ln}p({\boldsymbol{\theta }}| { \mathcal D })=\mathrm{ln}p({ \mathcal D }| {\boldsymbol{\theta }})+\mathrm{ln}p({\boldsymbol{\theta }})$, with

$$\begin{eqnarray}&&\mathrm{ln}p({ \mathcal D }| {\boldsymbol{\theta }})=\displaystyle \sum _{i}\mathrm{ln}p({{ \mathcal D }}_{i}| {\boldsymbol{\theta }})\end{eqnarray} \tag{ 19 }$$

being the marginal log likelihood for the full data set.

To determine the best-fit parameters and their uncertainties, we sample the posterior probability over the parameter space with Goodman & Weare's affine-invariant Markov chain Monte Carlo (Goodman & Weare 2010), making use of the Python module emcee (Foreman-Mackey et al. 2013).

The final best-fit values of the model parameters have been estimated using the median of the posterior distributions; the uncertainties have been estimated using the 15.87th and 84.13th percentiles. For a parameter whose probability distribution function (pdf) can be well described by a Gaussian distribution, the difference between the 15.87th and 84.13th percentile is equal to 1σ.

The calculation of the normalization integral in Equation (18) is complicated by the presence of the selection function ${ \mathcal S }$, leading to the fact that in some regions of the integrated space the integrand function is not continuous and shows an abrupt decrease to 0. For this reason, the classical multidimensional quadrature methods in Python are not able to give robust results. We decided to calculate the integral instead on a fine regular grid that is (Δl = 1°) ×(Δb = 1°) × (δD = 1 kpc) wide.

3.3.1. Model Priors

We now lay out the "pertinent range," across which the model priors are given. We set a different prior distribution $p({\boldsymbol{\theta }})$ for each of the five following cases:

Power-law model:

$$\begin{eqnarray}\mathrm{ln}p({\boldsymbol{\theta }}) & = & \mathrm{Uniform}(1.0\lt n\lt 6.0)\\ & & +\mathrm{Uniform}(0.1\lt q\lt 1.0).\end{eqnarray} \tag{ 20 }$$

BPL model:

$$\begin{eqnarray}\mathrm{ln}p({\boldsymbol{\theta }}) & = & \mathrm{Uniform}(1.0\lt n\lt 6.0)\\ & & +\mathrm{Uniform}(0.1\lt q\lt 1.0)\\ & & +\mathrm{Uniform}(\mathrm{log}({R}_{\min })\lt \mathrm{log}({r}_{\mathrm{break}})\lt \mathrm{log}({R}_{\max })),\end{eqnarray} \tag{ 21 }$$

where R_min, R_max give the galactocentric distance range available in the sample.

Einasto profile:

$$\begin{eqnarray}\mathrm{ln}p({\boldsymbol{\theta }}) & = & \mathrm{Uniform}(\mathrm{log}(0.01)\lt \mathrm{log}({r}_{\mathrm{eff}})\lt \mathrm{log}(50))\\ & & +\mathrm{Uniform}(\mathrm{log}(0.01)\lt \mathrm{log}({r}_{0})\lt \mathrm{log}(50))\\ & & +\mathrm{Uniform}(0.5\lt n\lt 20.0)\\ & & +\mathrm{Uniform}(0.1\lt q\lt 1.0).\end{eqnarray} \tag{ 22 }$$

Power-law model with q(R_gc):

$$\begin{eqnarray}\mathrm{ln}p({\boldsymbol{\theta }}) & = & \mathrm{Uniform}(\mathrm{log}(0.01)\lt \mathrm{log}({r}_{0})\lt \mathrm{log}(50))\\ & & +\mathrm{Uniform}(1.0\lt n\lt 5.0)\\ & & +\mathrm{Uniform}(0.1\lt {q}_{0}\lt 1.0)\\ & & +\mathrm{Uniform}(0.1\lt {q}_{\infty }\lt 1.0).\end{eqnarray} \tag{ 23 }$$

Einasto profile with q(R_gc):

$$\begin{eqnarray}\mathrm{ln}p({\boldsymbol{\theta }}) & = & \mathrm{Uniform}(\mathrm{log}(0.01)\lt \mathrm{log}({r}_{\mathrm{eff}})\lt \mathrm{log}(50))\\ & & +\mathrm{Uniform}(\mathrm{log}(0.01)\lt \mathrm{log}({r}_{0})\lt \mathrm{log}(50))\\ & & +\mathrm{Uniform}(\mathrm{log}(0.01)\lt \mathrm{log}({r}_{0})\lt \mathrm{log}(50))\\ & & +\mathrm{Uniform}(0.5\lt n\lt 20.0)\\ & & +\mathrm{Uniform}(0.1\lt {q}_{0}\lt 1.0)\\ & & +\mathrm{Uniform}(0.1\lt {q}_{\infty }\lt 1.0).\end{eqnarray} \tag{ 24 }$$

In order to test the methodology for fitting the density as discussed in Section 3, we created mock data samples of RR Lyrae stars in the Galactic halo, which should have the same properties as the observed sample of RRab stars, using a combination of a density law and assumptions on the selection function imposed by both PS1 3π and our selection cuts (Section 3.2). In detail, we first sampled ∼50,000 stars from mock halos generated with an underlying density given by a power law, Einasto profile, power law with q(R_gc), or Einasto profile with q(R_gc). We then applied a 3% distance uncertainty, superimposed the sample with faint and far Gaussian blobs away from the regions excluded by the selection function to simulate unknown overdensities, added the RRab known as members of the Sgr stream, and then applied the selection function. After that, the sample has ∼12,000 sources, and we randomly sample 11,025 sources to match the cleaned observed sample.

An example of a simulated distribution of halo RR Lyrae is shown in the top panel of Figure 17.

We then run the same analysis code on this sample as for the PS1 3π RRab sample. This enables us to estimate which halo properties we are able to identify and constrain with our approach.

We find results that are consistent with the input model within reasonable uncertainties, which means that we are able to recover the input parameters for all models in their assumed parameter range, and compare well with results we got from the PS1 3π data.

The one- and two-dimensional projections of the pdfs for fitting one of these mock halos, along with the parameters used to generate the mock halo, are given in Figure 3.

Zoom In Zoom Out Reset image size

Based on the five models described above in Section 3.1 and the selection function as described in Section 3.2, we apply our likelihood approach (Section 3.3) in order to constrain the best-fit model parameters.

We estimate those best-fit model parameters for the complete cleaned RRab sample, which spans 3/4 of the sky and contains 8917 sources. Figure 4 compares the observed number density of RR Lyrae stars with the density predicted by best-fit models.

Zoom In Zoom Out Reset image size

Table 2 summarizes the best-fit parameters of our five halo density models. For each model, we give the type of the density model, its best-fitting parameters along with their 1σ uncertainties estimated as the 15.87th and 84.13th percentiles, and the maximum log likelihood $\mathrm{ln}({{ \mathcal L }}_{\max })$. We also give the BIC, a measure for model comparison described in Section 4.2.

The one- and two-dimensional projections of the pdf for each model are given in Figures 5–9.

Table 2.  Best-fitting Halo Models

Density Model	Best-fit Parameters	$\mathrm{ln}({{ \mathcal L }}_{\max })$	BIC	ΔBIC
Power-law model	$q={0.918}_{-0.014}^{+0.016}$, $n={4.40}_{-0.04}^{+0.05}$	−157625	315269	203
BPL model	${r}_{\mathrm{break}}={38.7}_{-0.58}^{+0.69}$, $q={0.908}_{-0.006}^{+0.008}$,
	${n}_{\mathrm{inner}}={4.97}_{-0.05}^{+0.02}$, ${n}_{\mathrm{outer}}={3.93}_{-0.04}^{+0.05}$	−214222	428464	113398
Einasto profile	reff = 1.07 ± 0.10 kpc,
	q = 0.923 ± 0.007, n = 9.53+0.27−0.28	−157685	315388	322
Power-law model with q(Rgc)	${r}_{0}={25.0}_{-1.7}^{+1.8}\,\mathrm{kpc}$, ${q}_{0}={0.773}_{-0.016}^{+0.017}$,
	${q}_{\infty }={0.998}_{-0.001}^{+0.002}$, n = 4.61 ± 0.03	−157524	315066	0
Einasto profile with q(Rgc)	${r}_{0}={26.7}_{-2.0}^{+2.2}\,\mathrm{kpc}$, q0 = 0.779 ± 0.018,
	${q}_{\infty }={0.998}_{-0.002}^{+0.001}$, ${r}_{\mathrm{eff}}={1.04}_{-0.13}^{+0.25}\,\mathrm{kpc}$,
	$n={8.78}_{-0.30}^{+0.33}$	−157582	315182	116

Note. Summary of our best-fitting halo density models. For each model, we give the type of the density model, its best-fitting parameters along with their 1σ uncertainties estimated as the 15.87th and 84.13th percentiles, the maximum log likelihood $\mathrm{ln}({{ \mathcal L }}_{\max })$, and the BIC. ΔBIC gives the difference between the BIC of the best-fit model, the power-law model with q(R_gc), and the model used.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the power-law and BPL models, the pdf shows an almost Gaussian-like distribution with no covariance between the model parameters q and n. For the Einasto profile, as for the power law, the concentration index n and the flattening parameter q show an almost Gaussian distribution with no covariance. The concentration index n is covariant with the effective radius parameter, r_eff. The pdf of the power-law model with q(R_gc) shows covariance, and the pdf is strongly distorted from a Gaussian distribution. For the Einasto profile with q(R_gc), the pdf is more complex and skewed. The fitting parameters r₀, q₀, and q_∞ show a covariance, but their marginalizations have a Gaussian-like appearance.

Among models with constant flattening, the distribution of the sources is reasonably well fit by a power-law model with $n={4.40}_{-0.04}^{+0.05}$ and a halo flattening of $q={0.918}_{-0.014}^{+0.016}$. Allowing for a break in the power-law profile, we find a break radius of ${r}_{\mathrm{break}}={38.7}_{-0.58}^{+0.69}$, a halo flattening of $q\,={0.908}_{-0.006}^{+0.008}$, and the inner and outer slopes ${n}_{\mathrm{inner}}={4.97}_{-0.05}^{+0.02}$ and ${n}_{\mathrm{outer}}={3.93}_{-0.04}^{+0.05}$, respectively. The distance distribution is fit comparably well by a model with an Einasto profile with $n={9.53}_{-0.28}^{+0.27}$, an effective radius r_eff = 1. 07 ± 0.10 kpc, and a halo flattening of q₀ = 0.923 ± 0.007. If we allow for a radius-dependent flattening q(R_gc), we find the best-fit parameters for a power-law model with q(R_gc) as ${r}_{0}={25.0}_{-1.8}^{+1.7}\,\mathrm{kpc}$, n = 4.61 ± 0.03, ${q}_{0}={0.773}_{-0.016}^{+0.017}$, and ${q}_{\infty }={0.998}_{-0.001}^{+0.002}$. The best-fit parameters for an Einasto profile with q(R_gc) are ${r}_{0}={26.7}_{-2.0}^{+2.2}\,\mathrm{kpc}$, q₀ = 0.779 ± 0.018, ${q}_{\infty }={0.998}_{-0.002}^{+0.001}$, ${r}_{\mathrm{eff}}={1.04}_{-0.13}^{+0.25}\,\mathrm{kpc}$, and $n={8.78}_{-0.30}^{+0.33}$.

We find here ${q}_{0}\lt {q}_{\infty }$ for both models with variable flattening, indicating that the inner halo is more flattened than the outer halo. Assuming a constant flattening q instead, its best-fit value is also consistent among the power-law and Einasto profile models.

For all five models, the best-fit values along with their 1σ uncertainties are summarized in Table 2.

Our results confirm that if a varying flattening is assumed, the halo profile has an r₀ close to 20 kpc and the inner halo is more flattened than the outer. This is also consistent with results by Carollo et al. (2007, 2010), as well as Xue et al. (2015), Das & Binney (2016), and Iorio et al. (2017). For a BPL, we cannot confirm the Deason et al. (2014) result of a steepening found beyond 65 kpc. We discuss our results in comparison with previous attempts in more detail in Section 5.1.

We have estimated the best-fitting parameters for each model. In addition to that, it is important to compare the results of different models to determine which of them gives the best description of the data.

The most reliable way would be to compute the ratio of the Bayesian evidence, which is defined as the integral of the likelihood over all of the parameter space, for each model in order to compare them. Especially in higher-dimensional parameter spaces, like the ones we deal with here, this turns out to be too computationally expensive. However, under the assumption that the posterior distributions are almost Gaussian, an approximation can be used, called the BIC (Schwarz 1978).

The BIC takes into account both the statistical goodness of fit and the number of parameters that have to be estimated to achieve this particular degree of fit, by imposing a penalty for increasing the number of parameters in order to avoid overfitting. The BIC is defined as

$$\begin{eqnarray}&&\mathrm{BIC}=\dim ({\boldsymbol{\theta }})\mathrm{ln}(N)-2\mathrm{ln}({{ \mathcal L }}_{\max }),\end{eqnarray} \tag{ 25 }$$

where ${\boldsymbol{\theta }}$ are the model parameters, N is the number of objects in the sample, and ${{ \mathcal L }}_{\max }$ is the maximum likelihood, where we defined the likelihood function in Equation (19) as $\mathrm{ln}p({ \mathcal D }| {\boldsymbol{\theta }})$.

Using the BIC for selecting a best-fit model, the model with lowest BIC is preferred.

We have computed the BIC for all of our models and show them in Table 2 along with the best-fit parameters.

According to the BIC, we find the best-fit model to be the power law with q(R_gc), followed by the Einasto profile with q(R_gc), the constant-flattening power law, the constant-flattening Einasto profile, and finally the BPL. As the values of BIC in Table 2 indicate, allowing for flattening variations makes for distinctly better fits to the distribution of the RRab stars.

However, attention should be paid to the shape of the posterior distribution. When calculating the BIC, it is assumed that the posterior distributions are reasonably comparable to a Gaussian. As we see from Figures 5–9, the power-law model and the Einasto profile have posterior distributions that compare well to a Gaussian distribution, whereas for the cases with q(R_gc) the posterior distributions are somewhat distorted and show also a covariance between parameters.

Another issue is whether a difference in BIC is significant. A rating of the strength of the evidence against the model with the higher BIC value is given in Kass & Raftery (1995): a ΔBIC > 10 indicates very strong evidence against the model with the higher BIC.

In Section 4.1, we estimated best-fit parameters for the complete cleaned RRab sample, which spans 3/4 of the sky. Here we estimate them on smaller parts of the sky. This will help us to resolve and identify possible local variations in the best-fit model, especially in the halo flattening q and steepness n. We also look for previously unknown overdensities that we might find owing to the spatial extent and depth of the RRab sample.

4.3.1. Fitting Hemispheres and Pencil Beams

We now fit the halo profile for both the north and south Galactic hemisphere independently, in order to explore what the effects on our models—of rather restrictive functional form—are. The north hemisphere contains 6880 sources, whereas the south hemisphere contains only 4145 sources because of the PS1 3π survey footprint.

The results of this fitting attempt are summarized in Table 3. What we find is that the steepness parameters n of all best-fit hemisphere models compare well for both the north and south Galactic hemisphere and also compare well with the fit for the complete halo. When taking a look at the flattening-related parameters, q, q₀, q_∞, r_eff, we find that for models with constant flattening (both the power-law and Einasto profile models) q_south < q < q_north. In the case of models with q(R_gc), we find that the value of parameter q₀ is smaller for the south than for the north hemisphere, q_0,south < q₀ < q_0,north, whereas the value of the parameter q_∞ is similar for both hemispheres. Furthermore, we find that r_0,north > r_0,south > r₀ for both the power law with q(R_gc) and the Einasto profile with q(R_gc).

Table 3.  Best-fitting Halo Models for Each Hemisphere

Density Model	Best-fit Parameters North Galactic Hemisphere	Best-fit Parameters South Galactic Hemisphere
Power-law model	$q={0.925}_{-0.009}^{+0.010}$, n = 4.36 ± 0.03	$q={0.852}_{-0.011}^{+0.010}$, n = 4.40 ± 0.04
Einasto profile	${r}_{\mathrm{eff}}={1.11}_{-0.10}^{+0.09}\,\mathrm{kpc}$,	${r}_{\mathrm{eff}}={1.18}_{-0.11}^{+0.09}\,\mathrm{kpc}$
	$q={0.934}_{-0.010}^{+0.009}$, $n={9.59}_{-0.26}^{+0.30}$	$q={0.851}_{-0.011}^{+0.013}$, $n={9.10}_{-0.28}^{+0.31}$
Power-law model with q(Rgc)	r0 = 29.2 ± 4.4 kpc, ${q}_{0}={0.831}_{-0.017}^{+0.031}$,	${r}_{0}={18.8}_{-1.6}^{+1.4}\,\mathrm{kpc}$, ${q}_{0}={0.515}_{-0.058}^{+0.027}$
	${q}_{\infty }={0.997}_{-0.001}^{+0.004}$, $n={4.53}_{-0.06}^{+0.04}$	${q}_{\infty }={0.998}_{-0.002}^{+0.004}$, $n={4.88}_{-0.03}^{+0.06}$
Einasto profile with $q({R}_{\mathrm{gc}})$	${r}_{0}={31.9}_{-3.1}^{+3.9}\,\mathrm{kpc}$, ${q}_{0}={0.837}_{-0.017}^{+0.036}$,	${r}_{0}={20.9}_{-2.0}^{+2.3}\,\mathrm{kpc}$, ${q}_{0}={0.545}_{-0.066}^{+0.040}$
	${q}_{\infty }={0.998}_{-0.003}^{+0.001}$, ${r}_{\mathrm{eff}}={1.00}_{-0.11}^{+0.v}\,\mathrm{kpc}$,	${q}_{\infty }={0.998}_{-0.005}^{+0.001}$, ${r}_{\mathrm{eff}}={1.14}_{-0.12}^{+0.16}\,\mathrm{kpc}$
	$n={9.07}_{-0.28}^{+0.35}$	$n={7.57}_{-0.23}^{+0.40}$

Note. Summary of our best-fitting halo density models. For each model, we give the type of the density model and its best-fitting parameters along with their 1σ uncertainties estimated as the 15.87th and 84.13th percentiles.

Download table as:  ASCIITypeset image

The results of finding q_south < q < q_north for models with constant flattening and q_0,south < q₀ < q_0,north, ${q}_{\infty ,\mathrm{south}}\,\sim {q}_{\infty ,\mathrm{north}}\sim {q}_{\infty ,\mathrm{south}}$, r_0,north > r_0,south > r₀ in the case of a radius-dependent flattening are consistent: by definition of q(R_gc) (Equation (7)), q₀ is the flattening at center, q_∞ is the flattening at large galactocentric radii, and r₀ is the exponential scale radius over which the change of flattening occurs. A larger r₀ means that the flattening of the inner halo, where we find ${q}_{0}\lt {q}_{\infty }$, is in force out to a larger radius than for a smaller r_eff, thus leading to a larger part of the halo being more flattened.

The generalized result is thus that the south Galactic hemisphere is somewhat more flattened than the north Galactic hemisphere.

We also tried fitting models to the data in disjoint pencil beams (Δl = 30°) × (Δb = 60°), to further understand possible local variations in the best-fit model, especially in the halo flattening q and steepness n.

The angular source number density for the cleaned RRab sample, given per (Δl = 30°) × (Δb = 60°) bin, is shown in Figure 10.

Zoom In Zoom Out Reset image size

The resulting best-fit parameters for the power-law model, power-law model with q(R_gc), Einasto profile, and Einasto profile with q(R_gc) are shown in Figures 13–16 and are given in Tables 4–7 along with their 1σ uncertainties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Parameters for the Power-law Model on Δl = 30°, Δb = 60° Bins

l	b	Sources	q	n
0	−90	158	$({0.307}_{-0.077}^{+0.126})$	${3.53}_{-0.19}^{+0.18}$
0	−30	532	${0.912}_{-0.068}^{+0.057}$	${4.23}_{-0.11}^{+0.12}$
0	30	327	${0.973}_{-0.039}^{+0.020}$	${3.914}_{-0.127}^{+0.131}$
30	−90	582	${0.557}_{-0.055}^{+0.052}$	${3.79}_{-0.09}^{+0.10}$
30	−30	1099	${0.896}_{-0.051}^{+0.050}$	${4.89}_{-0.08}^{+0.06}$
30	30	289	${0.943}_{-0.053}^{+0.04}$	${4.25}_{-0.14}^{+0.16}$
60	−90	529	${0.714}_{-0.063}^{+0.061}$	$({3.39}_{-0.09}^{+0.08})$
60	−30	899	${0.967}_{-0.041}^{+0.024}$	${4.71}_{-0.08}^{+0.08}$
60	30	325	${0.919}_{-0.054}^{+0.048}$	${4.60}_{-0.14}^{+0.13}$
90	−90	260	${0.813}_{-0.076}^{+0.076}$	${3.86}_{-0.14}^{+0.14}$
90	−30	428	${0.949}_{-0.058}^{+0.037}$	${4.72}_{-0.12}^{+0.12}$
90	30	247	${0.927}_{-0.056}^{+0.0476}$	${4.862}_{-0.129}^{+0.094}$
120	−90	172	${0.568}_{-0.085}^{+0.081}$	${4.02}_{-0.17}^{+0.19}$
120	−30	312	${0.935}_{-0.084}^{+0.046}$	${4.92}_{-0.09}^{+0.06}$
120	30	232	${0.960}_{-0.042}^{+0.027}$	${4.80}_{-0.17}^{+0.12}$
150	−90	74	${0.447}_{-0.159}^{+0.168}$	${4.04}_{-0.28}^{+0.32}$
150	−30	215	${0.952}_{-0.065}^{+0.036}$	${4.71}_{-0.17}^{+0.14}$
150	30	277	${0.919}_{-0.054}^{+0.048}$	${4.82}_{-0.13}^{+0.11}$
180	−90	161	${0.940}_{-0.061}^{+0.041}$	${4.88}_{-0.13}^{+0.09}$
180	−30	318	${0.967}_{-0.051}^{+0.026}$	$({4.30}_{-0.12}^{+0.12})$
180	30	377	$({0.730}_{-0.067}^{+0.068})$	$({2.83}_{-0.11}^{+0.10})$
210	−90	98	${0.640}_{-0.099}^{+0.100}$	${4.61}_{-0.25}^{+0.22}$
210	−30	387	${0.942}_{-0.077}^{+0.042}$	${4.95}_{-0.06}^{+0.04}$
210	30	402	${0.934}_{-0.052}^{+0.042}$	${3.84}_{-0.11}^{+0.11}$
240	−30	292	${0.959}_{-0.054}^{+0.030}$	${4.90}_{-0.11}^{+0.07}$
240	30	476	${0.933}_{-0.043}^{+0.039}$	${4.79}_{-0.11}^{+0.11}$
270	−30	20	${0.818}_{-0.193}^{+0.128}$	${4.74}_{-0.38}^{+0.20}$
270	30	521	${0.967}_{-0.033}^{+0.023}$	${4.70}_{-0.11}^{+0.11}$
300	−30	2	$({0.458}_{-0.191}^{+0.318})$	$({2.20}_{-0.93}^{+1.66})$
300	30	520	${0.992}_{-0.012}^{+0.006}$	${3.44}_{-0.10}^{+0.10}$
330	−30	102	${0.953}_{-0.066}^{+0.036}$	${4.18}_{-0.28}^{+0.29}$
330	30	390	${0.976}_{-0.031}^{+0.018}$	${3.55}_{-0.11}^{+0.11}$

Note. Best-fit values for the power-law model, when carried out on Δl = 30°, Δb = 60° bins. The table gives the bin limits from (l, b) to (l + Δl, b + Δb), the number of sources contained, as well as the best-fit model parameters. Only bins containing sources are listed. Values in brackets are unreliable for reasons mentioned in Section 4.3. The (l, b) distribution of these table values is depicted in Figure 13.

Download table as:  ASCIITypeset image

Table 5.  Best-fit Parameters for the Einasto Profile on Δl = 30°, Δb = 60° Bins

l	b	Sources	reff	q	n
0	−90	158	${1.26}_{-0.10}^{+0.09}$	$({0.539}_{-0.095}^{+0.094})$	${10.4}_{-0.7}^{+0.7}$
0	−30	532	${1.22}_{-0.09}^{+0.09}$	${0.901}_{-0.071}^{+0.063}$	${8.90}_{-0.52}^{+0.60}$
0	30	327	${1.25}_{-0.10}^{+0.10}$	${0.972}_{-0.034}^{+0.020}$	${10.0}_{-0.6}^{+0.7}$
30	−90	582	${1.27}_{-0.10}^{+0.09}$	${0.629}_{-0.048}^{+0.046}$	${11.6}_{-0.6}^{+0.6}$
30	−30	1099	${1.18}_{-0.11}^{+0.10}$	${0.865}_{-0.051}^{+0.052}$	${6.83}_{-0.26}^{+0.30}$
30	30	289	${1.23}_{-0.10}^{+0.10}$	${0.949}_{-0.050}^{+0.034}$	${9.23}_{-0.58}^{+0.63}$
60	−90	529	${1.34}_{-0.10}^{+0.10}$	${0.830}_{-0.052}^{+0.050}$	$({12.8}_{-0.6}^{+0.6})$
60	−30	899	${1.18}_{-0.10}^{+0.10}$	${0.966}_{-0.047}^{+0.026}$	${7.62}_{-0.33}^{+0.37}$
60	30	325	${1.20}_{-0.10}^{+0.09}$	${0.929}_{-0.049}^{+0.043}$	${8.07}_{-0.48}^{+0.54}$
90	−90	260	${1.25}_{-0.09}^{+0.09}$	${0.896}_{-0.063}^{+0.060}$	${10.3}_{-0.6}^{+0.7}$
90	−30	428	${1.19}_{-0.09}^{+0.10}$	${0.943}_{-0.071}^{+0.041}$	${7.51}_{-0.39}^{+0.46}$
90	30	247	${1.17}_{-0.10}^{+0.10}$	${0.923}_{-0.057}^{+0.045}$	${6.97}_{-0.45}^{+0.51}$
120	−90	172	${1.24}_{-0.10}^{+0.09}$	${0.665}_{-0.070}^{+0.071}$	${9.52}_{-0.65}^{+0.69}$
120	−30	312	${1.184}_{-0.101}^{+0.099}$	${0.926}_{-0.076}^{+0.053}$	${6.70}_{-0.35}^{+0.41}$
120	30	232	${1.18}_{-0.10}^{+0.10}$	${0.964}_{-0.040}^{+0.025}$	${7.356}_{-0.457}^{+0.488}$
150	−90	74	${1.21}_{-0.10}^{+0.10}$	${0.694}_{-0.111}^{+0.117}$	${8.68}_{-0.73}^{+0.76}$
150	−30	215	${1.20}_{-0.10}^{+0.10}$	${0.956}_{-0.064}^{+0.032}$	${7.56}_{-0.50}^{+0.56}$
150	30	277	${1.19}_{-0.11}^{+0.09}$	${0.914}_{-0.049}^{+0.049}$	${7.36}_{-0.45}^{+0.53}$
180	−90	161	${1.17}_{-0.09}^{+0.11}$	${0.938}_{-0.059}^{+0.043}$	${6.69}_{-0.51}^{+0.58}$
180	−30	318	${1.23}_{-0.10}^{+0.09}$	${0.972}_{-0.044}^{+0.021}$	$({9.17}_{-0.52}^{+0.54})$
180	30	377	${1.37}_{-0.09}^{+0.09}$	${0.931}_{-0.048}^{+0.043}$	$({13.8}_{-0.6}^{+0.7})$
210	−90	98	${1.19}_{-0.10}^{+0.10}$	${0.666}_{-0.088}^{+0.093}$	${7.84}_{-0.71}^{+0.71}$
210	−30	387	${1.17}_{-0.10}^{+0.11}$	${0.944}_{-0.073}^{+0.041}$	${6.48}_{-0.36}^{+0.37}$
210	30	402	${1.26}_{-0.10}^{+0.09}$	${0.961}_{-0.040}^{+0.026}$	${10.8}_{-0.56}^{+0.62}$
240	−30	292	${1.16}_{-0.10}^{+0.09}$	${0.961}_{-0.048}^{+0.028}$	${6.93}_{-0.39}^{+0.47}$
240	30	476	${1.19}_{-0.10}^{+0.10}$	${0.931}_{-0.040}^{+0.041}$	${7.43}_{-0.36}^{+0.41}$
270	−30	20	${1.19}_{-0.10}^{+0.10}$	${0.814}_{-0.194}^{+0.129}$	${6.68}_{-1.01}^{+0.99}$
270	30	521	${1.18}_{-0.10}^{+0.10}$	${0.966}_{-0.034}^{+0.023}$	${7.80}_{-0.42}^{+0.45}$
300	−30	2	${1.20}_{-0.10}^{+0.10}$	$({0.453}_{-0.186}^{+0.330})$	${7.30}_{-0.97}^{+1.03}$
300	30	520	${1.30}_{-0.09}^{+0.09}$	${0.993}_{-0.010}^{+0.005}$	$({12.2}_{-0.7}^{+0.6}$
330	−30	102	${1.21}_{-0.10}^{+0.10}$	${0.947}_{-0.073}^{+0.037}$	${8.36}_{-0.76}^{+0.77})$
330	30	390	${1.28}_{-0.09}^{+0.10}$	${0.981}_{-0.027}^{+0.014}$	$({11.4}_{-0.6}^{+0.6})$

Note. Best-fit values for the Einasto profile, when carried out on ${\rm{\Delta }}l=30^\circ$, ${\rm{\Delta }}b=60^\circ$ bins. The table gives the bin limits from (l, b) to (l + Δl, b + Δb), the number of sources contained, as well as the best-fit model parameters. Only bins containing sources are listed. Values in brackets are unreliable for reasons mentioned in Section 4.3. The (l, b) distribution of these table values is depicted in Figure 14.

Download table as:  ASCIITypeset image

Table 6.  Best-fit Parameters for the Power-law Model with q(R_gc) on Δl = 30°, Δb = 60° Bins

l	b	Sources	r0	q0	q∞	n
0	−90	159	${18.2}_{-2.4}^{+2.7}$	$({0.204}_{-0.020}^{+0.051})$	${0.789}_{-0.111}^{+0.122}$	${4.99}_{-0.32}^{+0.14}$
0	−30	517	${39.5}_{-8.3}^{+21.9}$	${0.996}_{-0.418}^{+0.210}$	$({0.211}_{-0.261}^{+0.099})$	${3.77}_{-0.19}^{+0.188}$
0	30	325	${32.2}_{-3.5}^{+20.7}$	${0.999}_{-0.410}^{+0.200}$	${0.991}_{-0.042}^{+0.022}$	${3.91}_{-0.15}^{+0.15}$
30	−90	577	${34.6}_{-3.8}^{+4.5}$	${0.472}_{-0.024}^{+0.027}$	${0.998}_{-0.001}^{+0.002}$	${4.99}_{-0.03}^{+0.14}$
30	−30	1086	${39.4}_{-9.1}^{+19.7}$	${0.878}_{-0.487}^{+0.110}$	${0.841}_{-0.089}^{+0.068}$	${4.91}_{-0.096}^{+0.063}$
30	30	287	${39.1}_{-16.0}^{+16.2}$	${0.855}_{-0.418}^{+0.110}$	${0.999}_{-0.058}^{+0.042}$	${4.40}_{-0.17}^{+0.18}$
60	−90	529	${31.1}_{-3.5}^{+4.0}$	${0.432}_{-0.047}^{+0.050}$	${0.999}_{-0.025}^{+0.012}$	${4.57}_{-0.209}^{+0.209}$
60	−30	903	${34.8}_{-3.3}^{+19.8}$	${0.996}_{-0.393}^{+0.255}$	${0.960}_{-0.054}^{+0.028}$	${4.71}_{-0.09}^{+0.09}$
60	30	326	${34.7}_{-10.2}^{+19.8}$	${0.821}_{-0.382}^{+0.110}$	${0.999}_{-0.060}^{+0.045}$	${4.83}_{-0.180}^{+0.169}$
90	−90	264	${37.8}_{-5.9}^{+4.8}$	${0.590}_{-0.065}^{+0.074}$	${0.999}_{-0.045}^{+0.019}$	${4.53}_{-0.23}^{+0.22}$
90	−30	431	${39.8}_{-5.8}^{+17.9}$	${0.976}_{-0.423}^{+0.194}$	${0.977}_{-0.089}^{+0.041}$	${4.75}_{-0.13}^{+0.12}$
90	30	249	${39.4}_{-4.4}^{+18.0}$	${0.919}_{-0.364}^{+0.223}$	${0.997}_{-0.059}^{+0.047}$	${4.94}_{-0.14}^{+0.09}$
120	−90	179	${37.0}_{-7.9}^{+5.9}$	${0.485}_{-0.069}^{+0.062}$	${0.984}_{-0.108}^{+0.084}$	${4.99}_{-0.24}^{+0.10}$
120	−30	314	${39.7}_{-12.2}^{+20.1}$	${0.992}_{-0.478}^{+0.135}$	$({0.407}_{-0.271}^{+0.085})$	${4.72}_{-0.13}^{+0.10}$
120	30	234	${34.8}_{-4.3}^{+16.8}$	${0.974}_{-0.391}^{+0.210}$	${0.999}_{-0.046}^{+0.032}$	${4.91}_{-0.17}^{+0.13}$
150	−90	77	${39.7}_{-6.0}^{+22.3}$	${0.497}_{-0.151}^{+0.270}$	${0.992}_{-0.171}^{+0.416}$	${4.99}_{-0.65}^{+0.82}$
150	−30	218	${38.0}_{-7.0}^{+15.6}$	${0.978}_{-0.318}^{+0.302}$	${0.998}_{-0.068}^{+0.036}$	${4.77}_{-0.15}^{+0.16}$
150	30	280	${39.4}_{-7.4}^{+21.4}$	${0.940}_{-0.497}^{+0.120}$	${0.977}_{-0.075}^{+0.049}$	${4.87}_{-0.17}^{+0.14}$
180	−90	166	${39.7}_{-25.2}^{+7.3}$	${0.997}_{-0.451}^{+0.044}$	${0.585}_{-0.309}^{+0.211}$	${4.60}_{-0.55}^{+0.26}$
180	−30	327	${39.6}_{-5.5}^{+17.0}$	${0.998}_{-0.357}^{+0.236}$	${0.996}_{-0.042}^{+0.022}$	${4.46}_{-0.12}^{+0.12}$
180	30	379	${39.7}_{-15.5}^{+8.7}$	${0.999}_{-0.157}^{+0.099}$	$({0.511}_{-0.139}^{+0.116})$	$({2.43}_{-0.23}^{+0.23})$
210	−90	102	${14.1}_{-10.2}^{+10.2}$	${0.934}_{-0.219}^{+0.192}$	${0.203}_{-0.182}^{+0.228}$	${2.92}_{-0.71}^{+0.65}$
210	−30	401	${37.2}_{-6.4}^{+21.1}$	${0.992}_{-0.397}^{+0.191}$	${0.899}_{-0.110}^{+0.062}$	${4.99}_{-0.06}^{+0.03}$
210	30	412	${39.2}_{-23.0}^{+9.3}$	${0.808}_{-0.148}^{+0.107}$	${0.995}_{-0.053}^{+0.028}$	${4.07}_{-0.16}^{+0.15}$
240	−30	299	${26.3}_{-4.26}^{+18.2}$	${0.998}_{-0.450}^{+0.219}$	${0.993}_{-0.085}^{+0.034}$	${4.99}_{-0.099}^{+0.060}$
240	30	485	${32.2}_{-11.3}^{+14.4}$	${0.830}_{-0.324}^{+0.093}$	${0.999}_{-0.043}^{+0.032}$	${4.97}_{-0.14}^{+0.10}$
270	−30	21	${39.1}_{-14.2}^{+11.8}$	${0.961}_{-0.230}^{+0.135}$	${0.723}_{-0.211}^{+0.373}$	${4.99}_{-0.544}^{+0.270}$
270	30	524	${39.2}_{-2.9}^{+19.3}$	${0.920}_{-0.456}^{+0.167}$	${0.998}_{-0.045}^{+0.027}$	${4.87}_{-0.12}^{+0.13}$
300	−30	2	${39.6}_{-13.2}^{+14.0}$	${0.969}_{-0.187}^{+0.342}$	${0.506}_{-0.231}^{+0.322}$	${3.53}_{-1.13}^{+1.53}$
300	30	520	${39.6}_{-1.8}^{+15.8}$	${0.997}_{-0.376}^{+0.314}$	${0.999}_{-0.017}^{+0.006}$	${3.50}_{-0.10}^{+0.10}$
330	−30	97	${38.6}_{-8.1}^{+22.7}$	${0.990}_{-0.425}^{+0.139}$	${0.971}_{-0.143}^{+0.061}$	${3.97}_{-0.270}^{+0.277}$
330	30	388	${34.92072}_{-2.63}^{+16.9}$	${0.998}_{-0.444}^{+0.243}$	${0.993}_{-0.036}^{+0.019}$	$({3.47}_{-0.12}^{+0.13})$

Note. Best-fit values for the power-law model with $q({R}_{\mathrm{gc}})$, when carried out on ${\rm{\Delta }}l=30^\circ$, Δb = 60° bins. The table gives the bin limits from (l, b) to (l + Δl, b + Δb), the number of sources contained, as well as the best-fit model parameters. Only bins containing sources are listed. Values in brackets are unreliable for reasons mentioned in Section 4.3. The (l, b) distribution of these table values is depicted in Figure 15.

Download table as:  ASCIITypeset image

Table 7.  Best-fit Parameters for the Einasto Profile with q(R_gc) on Δl = 30°, Δb = 60° Bins

l	b	Sources	r0	q0	q∞	reff	n
0	−90	158	${23.4}_{-4.8}^{+1.0}$	$({0.200}_{-0.01}^{+0.06})$	${0.987}_{-0.230}^{+0.031}$	${1.18}_{-0.10}^{+0.10}$	${5.15}_{-0.38}^{+2.2}$
0	−30	532	${13.2}_{-0.04}^{+6.5}$	${0.621}_{-0.032}^{+0.303}$	${0.995}_{-0.163}^{+0.017}$	${1.21}_{-0.08}^{+0.11}$	${8.80}_{-0.49}^{+0.73}$
0	30	327	${25.3}_{-10.6}^{+4.4}$	${0.840}_{-0.107}^{+0.070}$	${0.999}_{-0.045}^{+0.004}$	${1.25}_{-0.11}^{+0.08}$	${9.55}_{-0.61}^{+0.82}$
30	−90	582	${35.6}_{-3.9}^{+2.4}$	${0.431}_{-0.033}^{+0.034}$	${1}_{-0.049}^{+0.005}$	${1.19}_{-0.12}^{+0.08}$	${6.53}_{-0.22}^{+0.82}$
30	−30	1099	${15.9}_{-1.6}^{+6.7}$	${0.997}_{-0.189}^{+0.031}$	$({0.762}_{-0.048}^{+0.156})$	${1.19}_{-0.11}^{+0.09}$	${7.42}_{-0.60}^{+0.15}$
30	30	289	${38.3}_{-19.7}^{+1.7}$	${0.800}_{-0.04}^{+0.12}$	${0.998}_{-0.075}^{+0.006}$	${1.21}_{-0.10}^{+0.10}$	${8.56}_{-0.41}^{+0.97}$
60	−90	529	${31.7}_{-4.3}^{+2.6}$	${0.418}_{-0.049}^{+0.026}$	${0.999}_{-0.034}^{+0.002}$	${1.198}_{-0.095}^{+0.102}$	${8.50}_{-0.46}^{+0.80}$
60	−30	899	${27.3}_{-13.0}^{+6.0}$	${0.999}_{-0.113}^{+0.011}$	${0.989}_{-0.131}^{+0.006}$	${1.21}_{-0.12}^{+0.08}$	${7.85}_{-0.32}^{+0.42}$
60	30	325	${37.3}_{-20.3}^{-1.5}$	${0.803}_{-0.060}^{+0.099}$	${0.999}_{-0.088}^{+0.008}$	${1.18}_{-0.091}^{+0.105}$	${7.60}_{-0.403}^{+0.766}$
90	−90	260	${36.1}_{-8.59}^{+0.246}$	${0.538}_{-0.066}^{+0.054}$	${0.999}_{-0.047}^{+0.004}$	${1.20}_{-0.10}^{+0.09}$	${8.37}_{-0.55}^{+0.77}$
90	−30	428	${37.7}_{-23.7}^{+4.9}$	${0.998}_{-0.229}^{+0.023}$	${0.999}_{-0.189}^{+0.021}$	${1.23}_{-0.14}^{+0.06}$	${7.64}_{-0.29}^{+0.64}$
90	30	247	${36.6}_{-22.5}^{+4.6}$	${0.860}_{-0.122}^{+0.077}$	${0.999}_{-0.112}^{+0.014}$	${1.193}_{-0.113}^{+0.086}$	${6.794}_{-0.354}^{+0.790}$
120	−90	172	${35.5}_{-10.4}^{+0.3}$	${0.399}_{-0.069}^{+0.064}$	${0.999}_{-0.141}^{+0.016}$	${1.20}_{-0.12}^{+0.08}$	${6.56}_{-0.16}^{+0.20}$
120	−30	312	${14.2}_{-0.2}^{+20.5}$	${0.998}_{-0.316}^{+0.026}$	${0.999}_{-0.305}^{+0.029}$	${1.19}_{-0.10}^{+0.10}$	${6.79}_{-0.28}^{+0.65}$
120	30	232	${37.9}_{-23.8}^{+6.1}$	${0.957}_{-0.147}^{+0.015}$	${0.999}_{-0.087}^{+0.009}$	${1.18}_{-0.09}^{+0.11}$	${7.34}_{-0.42}^{+0.69}$
150	−90	74	${39.7}_{-11.5}^{+1.1}$	${0.429}_{-0.088}^{+0.080}$	${0.993}_{-0.233}^{+0.026}$	${1.18}_{-0.10}^{+0.10}$	${6.85}_{-0.38}^{+0.23}$
150	−30	215	${38.4}_{-25.4}^{+12.9}$	${0.999}_{-0.652}^{+0.077}$	${0.996}_{-0.109}^{+0.08}$	${1.19}_{-0.10}^{+0.10}$	${7.67}_{-0.56}^{+0.51}$
150	30	277	${14.0}_{-0.2}^{+8.7}$	${0.997}_{-0.199}^{+0.032}$	${0.916}_{-0.086}^{+0.051}$	${1.18}_{-0.10}^{+0.10}$	${7.62}_{-0.61}^{+0.57}$
180	−90	161	${32.5}_{-16.0}^{+3.7}$	${0.997}_{-0.118}^{+0.010}$	${0.749}_{-0.080}^{+0.183}$	${1.22}_{-0.13}^{+0.08}$	${7.45}_{-0.83}^{+0.75}$
180	−30	318	${12.2}_{-0.9}^{+8.5}$	${0.994}_{-0.434}^{+0.031}$	${0.997}_{-0.087}^{+0.004}$	${1.22}_{-0.11}^{+0.107}$	$({9.30}_{-0.56}^{+0.65})$
180	30	377	${34.4}_{-6.1}^{+0.8}$	$({0.328}_{-0.058}^{+0.034})$	${0.998}_{-0.060}^{+0.005}$	${1.30}_{-0.11}^{+0.09}$	$({11.4}_{-0.6}^{+0.8})$
210	−90	98	${38.5}_{-24.9}^{+6.2}$	${0.6}_{-0.2}^{+0.1}$	${0.748}_{-0.164}^{+0.070}$	${1.22}_{-0.12}^{+0.08}$	${7.48}_{-0.67}^{+0.10}$
210	−30	387	${32.2}_{-18.4}^{+3.4}$	${0.997}_{-0.297}^{+0.032}$	${0.999}_{-0.152}^{+0.016}$	${1.16}_{-0.10}^{+0.11}$	${6.62}_{-0.32}^{+0.41}$
210	30	402	${31.8}_{-9.2}^{+3.1}$	${0.688}_{-0.079}^{+0.067}$	${0.999}_{-0.040}^{+0.004}$	${1.25}_{-0.11}^{+0.08}$	${9.81}_{-0.46}^{+0.82}$
240	−30	292	${36.7}_{-22.4}^{+2.7}$	${0.997}_{-0.174}^{+0.014}$	${0.996}_{-0.162}^{+0.014}$	${1.19}_{-0.12}^{+0.08}$	${6.98}_{-0.32}^{+0.58}$
240	30	476	${29.0}_{-13.0}^{+3.57}$	${0.786}_{-0.061}^{+0.095}$	${0.999}_{-0.062}^{+0.006}$	${1.18}_{-0.10}^{+0.10}$	${6.92}_{-0.25}^{+0.68}$
270	−30	20	${18.8}_{-4.10}^{+4.52}$	${0.992}_{-0.390}^{+0.045}$	$({0.544}_{-0.261}^{+0.250})$	${1.20}_{-0.11}^{+0.08}$	${7.10}_{-1.02}^{+0.98}$
270	30	521	${38.0}_{-23.0}^{+3.49}$	${0.939}_{-0.062}^{+0.043}$	${0.999}_{-0.083}^{+0.009}$	${1.17}_{-0.08}^{+0.12}$	${7.83}_{-0.40}^{+0.60}$
300	−30	2	${13.0}_{-1.3}^{+9.8}$	${0.984}_{-0.722}^{+0.193}$	${0.999}_{-0.733}^{+0.197}$	${1.20}_{-0.10}^{+0.09}$	${7.65}_{-1.27}^{+0.74}$
300	30	520	${26.0}_{-11.6}^{+6.7}$	${0.919}_{-0.087}^{+0.058}$	${0.999}_{-0.019}^{+0.001}$	${1.32}_{-0.14}^{+0.09}$	$({11.9}_{-0.7}^{+0.9})$
330	−30	102	${23.2}_{-9.2}^{+9.9}$	${0.988}_{-0.248}^{+0.018}$	${0.999}_{-0.190}^{+0.021}$	${1.19}_{-0.07}^{+0.12}$	${8.51}_{-0.86}^{+0.65}$
330	30	390	${27.6}_{-7.71}^{+6.78}$	${0.748}_{-0.066}^{+0.090}$	${0.999}_{-0.029}^{+0.001}$	${1.26}_{-0.10}^{+0.12}$	${10.3}_{-0.5}^{+0.1}$

Note. Best-fit values for the Einasto profile with $q({R}_{\mathrm{gc}})$, when carried out on ${\rm{\Delta }}l=30^\circ$, ${\rm{\Delta }}b=60^\circ$ bins. The table gives the bin limits from (l, b) to (l + Δl, b + Δb), the number of sources contained, as well as the best-fit model parameters. Only bins containing sources are listed. Values in brackets are unreliable for reasons mentioned in Section 4.3. The (l, b) distribution of these table values is depicted in Figure 16.

Download table as:  ASCIITypeset image

The fitting procedure also works well with small pieces of the sky. As an example, we show the fitted models for two small patches on the sky, 240° < l < 270°, −30° < b < 30° and 30° < l < 60°, −90° < b < −30° (see Figure 11). To illustrate the fitting performance further, in Figure 12 we give the posterior probability distribution in the case of fitting a power law with varying flattening q(R_gc) (Equation (7)) to a 30° × 60° patch of mock data. The posterior distribution is comparable to those when fitting the same halo profile to the full cleaned sample (see Figure 3 for comparison). However, the width of the posterior probability distribution increases compared to the full cleaned sample. This is also reflected in the 1σ intervals given along with the best-fit parameters in the top right part of the figure.

Starting with the results of fitting the power-law model to the (Δl = 30°) × (Δb = 60°) patches (Figure 13), we find that the flattening parameter q is homogeneous over almost the complete sky. There are a few exceptions, i.e., for 0° < l < 30°, −90° < b < −30°, the resulting flattening parameter q is suspiciously small. However, within that region, there are only 22 sources, which makes a reliable fit difficult.

Again for 300 < l < 330, −30 < b < 30, the resulting q, as well as here the power-law index n, is small. Since there are only two sources within that region, we obviously have to exclude that fit. Outside of these regions, the resulting q and n are relatively homogeneous, with a trend to smaller n near the edges of the survey (see white empty region at l > 240° in the figures) and at very high latitudes.

For the Einasto profile (Figure 14) we also find regions on the sky where the fitting parameters are considerably deviating. For the parameter n, this is especially the case for 180° < l < 210°, −30° < b < 90°, as well as at some regions at high latitudes. In those cases, the best-fitting n is sometimes much higher and sometimes much smaller than for the power-law model; however, this is a result of the different definition of n in both models (see Equation (2) versus Equation (6), and the steepness of the Einasto profile not being constant but changing continuously as given by Equation (5)), and the fitted profiles look comparable.

In the case of a power law with q(R_gc), as shown in Figure 15, the best-fit values for q₀, q_∞ are more similar than for the fit of the complete halo. In general, as is clearly visible in Figure 15, the distribution of the halo-flattening parameters q₀, q_∞, and the power-law index n shows more scatter than for a power law with constant halo flattening. This might be caused by the model tending to overfit the data, a problem common to higher-dimensional models, overreacting to fluctuations in the underlying data set that should be fitted.

In the case of an Einasto profile with q(R_gc), as shown in Figure 16, the best-fit values for q₀, q_∞ are again more similar than for the fit of the complete halo. We find about the same deviations as reported for the other models, such as unreliable fits at 0° < l < 30°, −90° < b < −30°, and 300 < l < 330, −30 < b < 30, due to the small number of sources within those regions.

Within 180° < l < 210°, the best-fit value of n is influenced by the presence of outskirts of the Sagittarius stream, which were not fully removed by our cuts. Within this region, a small number of stars from the stream appear to be present, and in general, the number of sources in this region of the sky is small after applying our cuts on overdensities. This is also the case for 300° < l < 330°, −30° < b < 90°. A higher-dimensional model is more affected by this than a lower-dimensional one; compare the extreme cases of the two-dimensional power-law model and the five-dimensional Einasto profile with q(R_gc).

Again, in those cases, the best-fitting n is much higher than for the case of a power-law model; however, this is a result of the different definition of n in each model (see Equation (2) versus Equation (6), and the steepness of the Einasto profile being not constant but changing continuously as given by Equation (5)), and the fitted profiles look comparable.

We give the mean and variance for the best-fit parameters on Δl = 30°, Δb = 60° bins for all four models in Table 8.

Table 8.  Mean and Variance for Best-fit Parameters on Δl = 30°, Δb = 60° Bins

Density Model	Parameter Mean and Variance
Power-law model	mean(n) = 4.25, Var(n) = 0.443,
	mean(q) = 0.840, Var(q) = 0.0333
Einasto profile	mean(reff) = 1.22, Var(reff) = 0.00252,
	mean(q) = 0.873, Var(q) = 0.0191
	mean(n) = 3.85, Var(n) = 8.76
Power-law model with q(Rgc)	mean(r0) = 36.0, Var(r0) = 36.3,
	mean(q0) = 0.856, Var(q0) = 0.0443,
	mean(q∞) = 0.859, Var(q∞) = 0.0559,
	mean(n) = 4.45, Var(n) = 0.441
Einasto profile with q(Rgc)	mean(r0) = 29.1, Var(r0) = 81.5,
	mean(q0) = 0.790, Var(q0) = 0.058,
	mean(q∞) = 0.958, Var(q∞) = 0.011,
	mean(reff) = 1.21, Var(reff) = 0.001,
	mean(n) = 7.961, Var(n) = 2.01

Note. Mean and variance for the best-fit parameters of all four halo models on Δl = 30°, Δb = 60° bins.

Download table as:  ASCIITypeset image

4.3.2. Density Residuals and Their Significance

Additionally, we compared the best-fit model to the PS1 3π RR Lyrae sample by calculating the residuals of that model. In Figure 17, we give density plots in the Cartesian reference frame (X, Y, Z) (see Equation (1)) for the best-fit model, as well as residuals for the observed cleaned sample of PS1 3π RRab stars. Densities are each color-coded according to the legend.

Zoom In Zoom Out Reset image size

The first row of Figure 17 shows a realization of a mock "cleaned sample" of 11,025 sources (the same number of sources as in the observed cleaned sample), sampled from the best-fit model, a power law with q(R_gc) with r₀ = 25.0 kpc, q₀ = 0.773, ${q}_{\infty }=0.998$, and n = 4.61, with the applied selection function.

This mock sample looks very comparable to the observed cleaned sample, Figure 2(b). The positions of the Galactic plane and Sgr stream, as removed by the selection function, are indicated. The dashed circle represents the R_gc > 20 kpc cut. Sources within the circle but farther away than 20 kpc are seen as a result of projection effects; the distinctly higher density just after 20 kpc shows the stars that are no longer affected by this distance cut.

We calculated the number density of our observed cleaned sample (given in Figure 2) at each (X, Y, Z), using a nearest-neighbor-based adaptive Bayesian density estimator (Ivezić et al. 2005; Sesar et al. 2013b), yielding ln (ρ_obs). The result is shown in the second row of Figure 17.

We then applied the same estimation of the 3D number density to 10 realizations of mock samples from the best-fit model; the resulting mean density is given in the third row of the figure as $\mathrm{ln}\langle {\rho }_{\mathrm{model}}\rangle$.

The logarithmic residuals of the best-fit model were calculated by subtracting the ln model mean number density (third row) from the observed number density (second row), yielding ${\rm{\Delta }}\mathrm{ln}\langle \rho \rangle =\mathrm{ln}({\rho }_{\mathrm{obs}})-\mathrm{ln}\langle {\rho }_{\mathrm{model}}\rangle$, as given in the last row of this figure. A ${\rm{\Delta }}\mathrm{ln}\langle \rho \rangle \lt 0$ indicates that the best-fit model overestimates the number densities, whereas a ${\rm{\Delta }}\mathrm{ln}\langle \rho \rangle \gt 0$ means that it underestimates the number density.

We find that the best-fit model leads, as expected, to a ${\rm{\Delta }}\mathrm{ln}\langle \rho \rangle \sim 0$ over wide ranges, but it also shows regions where the model underestimates the number density (yellow to red). This can be due to selection effects from the PS1 3π observing strategy, but it can also be an indicator for unknown structure and overdensities, as well as a more distorted halo shape. Also, our finding that the flattening is different for both hemispheres points toward a halo structure that is more complex than just an ellipsoid.

There are also regions showing slightly negative, near-zero residuals. As we draw the same number of stars from the mock sample as were observed, the overall density is naturally slightly overpredicted if there are underpredicted regions (regions in the PS1 RRab sample containing previously unknown overdensities) in order to match the total number of sources. This leads to slightly negative residuals when comparing the observed and the mock sample. A similar behavior is shown in Sesar et al. (2013b), Figure 10. They illustrate that this behavior is also found when fitting a mock data sample consisting of an underlying halo profile with added diffuse overdensities: slightly negative residuals are found over a wide area owing to a clumpy halo, i.e., a halo with diffuse overdensities. As indicated when discussing the selection function in Section 3.2, our estimations of purity and completeness are rather uncertain beyond distances of 90 kpc.

The dark-blue regions in this plot occur as a result of edge effects when the samples become sparse at the survey's outskirts. The diffuse overdensities thus revealed are of further interest; we will discuss them in more detail in Section 4.3.3 and thus label them in Figure 17 according to Sesar et al. (2007).

Also, as the relative sparseness of our cleaned RRab sample (11,025 sources within almost 3/4 of the sky and an extent of 20 kpc < R_gc < 131 kpc) introduces local number density fluctuations even for a smooth underlying density distribution, we have to estimate the significance of these overdensities. To do so, we carry out the following approach:

• 1.  
  We bootstrap the observed RRab sample N = 50 times. We estimate the density of each of these bootstrapped samples, using the density estimator by Ivezić et al. (2005), resulting in ρ_obs,i for i = 1 ... N. We fit each of these bootstrapped samples, sample each of them 10 times, and get the mean model density using the density estimator. This yields $\langle {\rho }_{\mathrm{model}}{\rangle }_{i}$ for i = 1 ... N, and further ${\rm{\Delta }}{(\mathrm{ln}\langle \rho \rangle )}_{i}\equiv \mathrm{ln}({\rho }_{\mathrm{obs},i})-\mathrm{ln}\langle {\rho }_{\mathrm{model},i}\rangle$ for i = 1...N.
• 2.  
  From the above, we can construct the variance $\sigma ({\rm{\Delta }}(\mathrm{ln}\langle \rho \rangle ))\equiv \mathrm{Var}\langle {\rm{\Delta }}{(\mathrm{ln}\langle \rho \rangle )}_{i}\rangle$.
• 3.  
  The 3D significance is then ${\rm{\Delta }}(\mathrm{ln}\langle \rho \rangle )/\sigma ({\rm{\Delta }}(\mathrm{ln}\langle \rho \rangle ))$.

The resulting variance and significance are shown in Figure 18, each projected using the mean. Per definition, the significance is 0 where ${\rm{\Delta }}(\mathrm{ln}\langle \rho \rangle )=0$.

Zoom In Zoom Out Reset image size

We find a significance of ∼20 to >50 at regions that coincide with the lower row of panels in Figure 17, and the variance is small and does not exceed 0.04–0.08 within these overdense regions. We count this as a strong indicator of these overdensities being real and not caused by Poisson number density fluctuations.

4.3.3. Overdensities

We compare the overdensities found by us with those discovered previously by Sesar et al. (2007, 2010).

In their studies they analyzed the spatial distribution of candidate RR Lyrae stars discovered by SDSS Stripe 82 along the Celestial Equator. They had used 634 RR Lyrae candidates from SDSS Stripe 82 and 296 RR Lyrae candidates from Ivezić et al. (2000) in their 2007 analysis (Sesar et al. 2007), and later on they cleaned the SDSS Stripe 82 sample of RR Lyrae (Sesar et al. 2010), using then 366 highly probable RRab stars.

In Figure 19, we plot the overdensities in an $({\rm{R}}.{\rm{A}}.,D)$ projection similar to Sesar et al. (2007) (see their Figure 13; see also Figure 11 in Sesar et al. 2010), using our full range in decl. and highlighting the region covered by their analysis. Uppercase letters denote overdensities found in the SDSS sample of Sesar et al. (2007, 2010), numbers denote overdensities found in their analysis of the Ivezić et al. (2000) sample (not numbered in Sesar et al. 2007), and lowercase letters denote overdensities we found in regions not covered by the analysis of Sesar et al. (2010).

Zoom In Zoom Out Reset image size

We can recover most of the overdensities found by Sesar et al. (2007, 2010), i.e., we recover their overdensities A, B, C, E, F, G, I, J, L. Among them, Sesar et al. (2010) claim that they do not find overdensities I and L they had found in their previous analysis and attribute this to their then better, cleaned sample of RR Lyrae stars. However, we find the overdensities I and L, where especially L stands out. We could verify that some overdensities found in Sesar et al. (2007) will disappear in a more cleaned sample, as shown in Sesar et al. (2010): consistent with Sesar et al. (2010), we do not find the overdensities D, H, K, and M. However, the overdensity D appears very small in Sesar et al. (2007), so we cannot say for sure if we are able to identify anything at this position. Our sample does not cover exactly the same extent in R.A. and decl., so the overdensities D, H, K, M as given in Sesar et al. (2007) lie in regions we do not cover. We have checked adjacent slices in decl. to see whether we might detect those overdensities, but we cannot find them. So our conclusion regarding those four overdensities is that either they are not real, as assumed by Sesar et al. (2007), or they have a small extent in decl.

The left wedge of Figure 19 compares to the left wedge of Figure 13 in Sesar et al. (2007), where they used RR Lyrae from Ivezić et al. (2005). We label these overdensities (1), (2), (3).

In regions not covered by Sesar et al. (2007), we detect many new overdensities out to D ≳ 100 kpc, continuing the overall distribution of overdensities found before. We label them by lowercase letters.

The strongest clump in the left wedge, (2), stems from stars belonging to the Sgr stream not being fully excised by our cuts. The same holds for the small and sparse overdensity C being part of the stream's trailing arm (Ivezić et al. 2003; Sesar et al. 2007).

Sesar et al. (2010) claim that the overdensity J is most probably a stellar stream, the Pisces oversensitivity (see also Watkins et al. 2009). In contrast to Sesar et al. (2007, 2010), we find that the overdensities J and L might be connected.

We now compare our results—especially halo steepness and flattening—to earlier findings from the many other groups that have attacked this important and interesting problem. In previous efforts to determine the halo shape, RR Lyrae and BHB stars have often been used as tracers because they are found in old populations, have precise distances, and are bright enough to be observed at radii out to ∼100 kpc (see, e.g., Xue et al. 2008, 2011; Sesar et al. 2010; Deason et al. 2011, 2013, 2014).

The major issue with BHB stars is potential confusion with blue stragglers and with QSOs. Samples of BHB stars must be carefully vetted to ensure that contamination is minimalized. This is not easy, especially as one moves farther out to fainter objects (see, e.g., Deason et al. 2012), where, in an effort to build a sample of BHB stars beyond 80 kpc, of 48 candidates selected photometrically, after the acquisition of low spectral resolution VLT-FORS2 spectra, only 7 turned out to be bona fide BHB stars. RR Lyrae, on the other hand, can be identified and verified with just photometric light curves, available from the application of modern machine-learning techniques to the databases of the large multiepoch photometric surveys that have been carried out over the past decade.

It is important to remember that our sample selection procedure that we apply to the Pan-STARRS1 3π database takes advantage of the experience gained by attempting to use the SDSS (Sesar et al. 2007, 2010) and then to analyze the more difficult PTF (Cohen et al. 2017) with few observations taken with a random cadence. Our PS1 sample of RRab stars (Sesar et al. 2017) is unique in that it contains a large number (44,403) of RR Lyrae in total, of which 17,452 lie within the radial range of R_gc from 20 up to 130 kpc, all with highly precise photometry. This yields a sample of high purity and completeness exceeding 80% out to R_gc = 80 kpc. Out of these sources, 11,025 are found outside of dense regions such as stellar streams, the Galactic disk and bulge, or globular clusters and stellar streams.

Furthermore, Sesar et al. (2017) have quantified this completeness with extensive testing throughout the entire radial range. Confusion with QSOs and with blue stragglers is eliminated by requiring a light curve characteristic of an RR Lyrae star.

First glimpses of the variation of the halo shape were already caught by Kinman et al. (1966), based on RR Lyrae stars as halo tracers. As a first attempt, Preston et al. (1991) argued that the flattening changes from q = 0.5 at 1 kpc to q ∼ 1 at 20 kpc. However, later work by Sluis & Arnold (1998) shows a constant flattening of q ∼ 0.5 without any evidence for a radius-dependent flattening. Recent work by De Propris et al. (2010) utilizing 666 BHB stars from the 2dF QSO Redshift Survey states that the halo is approximately spherical with a power-law index of ∼2.5 out to 100 kpc. Sesar et al. (2010), who studied main-sequence turnoff stars from the Canada–France–Hawaii Telescope Legacy Survey, find that the flattening is approximately constant at q ∼ 0.7 out to 35 kpc.

Carollo et al. (2007, 2010) found that the inner halo is highly flattened with axis ratios of q ∼ 0.6, whereas the outer halo is more spherical with axis ratios of q ∼ 0.9. In contrast to us, they include stars as close as 2 kpc in their fitting and thus get a more pronounced flattening within about 5–10 kpc.

Others also find evidence that at least the innermost part of the halo is quite flattened: Sesar et al. (2010) fit the Galactic halo profile based on ∼5000 RR Lyrae stars from the recalibrated LINEAR data set, spanning 5 kpc < D < 30 kpc over ∼8000 deg² of the sky. They find for their best-fit model an oblate ellipsoid with an axis ratio of q = 0.63 and a double power-law model with q = 0.65, n_inner = 1, n_outer = 2.7, and r_break = 16 kpc.

Hence, based on the different distances span by the aforementioned work, there is much evidence that the innermost part of the halo is quite flattened, that the outer part of the halo is more spherical, and that our results confirm that.

As a reason for the smaller minor-to-major-axis ratio q found for the inner halo, Deason et al. (2011) state that inner-halo stars possess generally high orbital eccentricities and exhibit a modest prograde rotation around the Galactic center. In contrast, stars in the outer halo exhibit a much more spherical spatial distribution as they cover a wide range of orbital eccentricities and show a retrograde rotation about the Galactic center.

For the density slope of the halo profile, many studies (e.g., Sesar et al. 2007, 2011; Watkins et al. 2009; Deason et al. 2011) find that it shows a shift from a relatively shallow one, as described by n ∼ 2.5, to a much steeper one outside of about 20–30 kpc that is consistent with n ∼ 4. The earliest evidence for that is from Saha (1985), who found that RR Lyrae are well described by a BPL with n ∼ 3 out to 25 kpc and n ∼ 5 beyond.

Subsequently, Sesar et al. (2011) used 27,544 near-turnoff MS stars out to 35 kpc selected from the Canada–France–Hawaii Telescope Legacy Survey to find a flattening of the stellar halo of 0.7 and a density distribution consistent with a BPL with an inner slope of 2.62 and an outer slope of 3.8 at the break radius of 28 kpc, or an equally good Einasto profile with a concentration index of 2.2 and an effective radius of 22.2 kpc.

Xue et al. (2015) probe the Galactic halo at $10\,\mathrm{kpc}\,\lt {R}_{\mathrm{gc}}\lt 80\,\mathrm{kpc}$ using 1757 stars from the SEGUE K-giant Survey. The majority of their sources are found at R_gc < 30 kpc, whereas in our sample 1093 RRab stars exist beyond a galactocentric distance of 80 kpc, and 238 beyond 100 kpc. They find that they can fit their sample by an Einasto profile with n = 3.1, r_eff = 15 kpc, and q = 0.7 and by an equally flattened BPL with n_in = 2.1, n_out = 3.8, and r_break = 18 kpc (this is something we had not applied), and when fitting by an Einasto profile with q(R_gc), they find that the halo is considerably more flattened as q changes from 0.55 ± 0.02 at 10 kpc to 0.8 ± 0.03 at large radii.

Bell et al. (2008) used ∼4 million color-selected MS turnoff stars from DR5 of the SDSS out to 40 kpc and find a best-fit flattening of the stellar halo of 0.5–0.8, and the density profile of the stellar halo is approximately described by a power law with index of 2–4.

Other estimates of the power-law index, or slope, of the halo give break radii or effective radii of ∼20–30 kpc and power-law slopes of n ∼ 3 (e.g., Deason et al. 2011, 2014; Sesar et al. 2011; Xue et al. 2015). For example, Sesar et al. (2011) fit the Galactic halo within heliocentric distances of <35 kpc, steeper at R_gc > 28 kpc, n_inner = 2.62, n_outer = 3.8, or a best-fit Einasto profile with n = 2.2, r_e = 22.2 kpc, and q = 0.7, where they found no evidence that it changes across the range of probed distances. Subsequently, Deason et al. (2014) found a very steep outer halo profile with a power law of 6 beyond 50 kpc, and yet steeper slopes of 6–10 at larger radii.

Our findings of n = 4.40–4.61 for a power-law model or n = 8.78–9.53 for an Einasto profile (keep in mind the different definitions of n) for our sample starting at R_gc = 20 kpc are thus in good agreement with most previous results, assuming no break radius and thus a constant-density slope or the slope variation as introduced by the Einasto profile (see Equation (5)). We claim that the estimate of the power-law parameter from De Propris et al. (2010), n ∼ 2.5, is too shallow, as well as that the estimate of the power-law parameter beyond 65 kpc from Deason et al. (2014), n = 6–10, is too steep, as we do not see such a drop from our fit or from our data.

We cannot verify results showing a break radius near 20 kpc, as our sample starts at 20 kpc and thus only a small number of sources would be found, if at all, within the break radius. However, the extent of our sample enables us to check the finding by Deason et al. (2014). They find a BPL with three ranges of subsequently steepening slope: 2.5 for 10 kpc < R_gc < 25 kpc, 4.5 for 25 $\mathrm{kpc}\lt {R}_{\mathrm{gc}}\lt 65\,\mathrm{kpc}$, and 10 for $65\,\mathrm{kpc}\lt {R}_{\mathrm{gc}}\lt 100\,\mathrm{kpc}$. Whereas in the distance range 25 kpc < R_gc < 65 kpc their profile agrees well with our results, we cannot probe the inner part, as our sample starts at 20 kpc, and our data argue against a significantly steeper profile at R_gc > 65 kpc as long as the PS1 selection function is applied. We thus fit a BPL to our sample, allowing for a break radius beyond 20 kpc, and find a slope of 3.93 beyond 39 kpc all the way out to the limit of our sample. There is no evidence for a steepening of the halo profile beyond 65 kpc as found by Deason et al. (2014), or any other indication that there is a truncation or break in the halo profile within the range we probe by the RRab sample. We roughly confirm their power-law slope of 4.8 for the region 25 kpc < R_gc < 45 kpc, as we find a power-law slope of 4.97 for 20 kpc < R_gc < 39 kpc. The small change in slope near 39 kpc may be related to how, in detail, the Sgr stream is removed. Beyond 40 kpc we find a much less steep power-law slope than do Deason et al. (2014).

Deason et al. (2013) interpret the presence or absence of a break as linked to the details of the stellar accretion history. They state that a prominent break can arise if the stellar halo is dominated by the debris from an accretion event that is massive, single, and early.

Xue et al. (2015) and Slater et al. (2016) find a halo profile that is shallower than our best-fit models; it is difficult to know whether this difference results from methodological differences or some intrinsic difference in the distribution between RR Lyrae and giants (Slater et al. 2016) or K giants (Xue et al. 2015).

Iorio et al. (2017) very recently carried out an attempt to map the Galactic inner halo with R_gc < 28 kpc based on a sample of 21,600 RR Lyrae from the Gaia and 2MASS surveys. They found that the best-fit model to describe the halo distribution is a power law with n = 2.96 ± 0.05, and flattening is present, resulting in an triaxial ellipsoid.

In Table 9, we give the best-fit parameters for the Galactic halo as found in other work, along with the distance range over which they were estimated. These models are visualized, together with our best-fit models, in Figure 20.

Table 9.  Model Parameters for Selected Other Surveys

Paper	Density Model	Parameter Mean and Variance	Distance Range
Watkins et al. (2009)	BPL	ninner = 2.4, ninner = 4.5,	5 kpc < Rgc < 115 kpc
		rbreak = 23 kpc, q = 1 (assumes no flattening)
Deason et al. (2014)	BPL with three segments	α1 = 2.5, α2 = 4.5, αout = 10,	10 kpc < Rgc < 100 kpc
		rc = 25 kpc, rbreak = 65 kpc, q = 1 (assumes no flattening)
Cohen et al. (2015)	power law	n = 3.8 ± 0.3, q = 1 (assumes no flattening)	50 kpc < Rgc < 115 kpc
Xue et al. (2015)	power law with q(Rgc)	α = 4.2 ± 0.1, q0 = 0.2 ± 0.1, ${q}_{\infty }=0.8\pm 0.3$, r0 = 6 ± 1 kpc	10 kpc < Rgc < 80 kpc
Slater et al. (2016)	power law	n = 3.5 ± 0.2, q = 1 (assumes no flattening)	20 kpc < Rgc < 80 kpc

Note. Mean and, as far as available, variance for the best-fit parameters found in other work. The models are shown, together with our best-fit models, in Figure 20. Some authors call the exponent of their power-law model α instead of n. We have kept their notation.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Figure 20 shows a remarkably different slope for the models based on giants (the yellow, orange, and green lines in this figure) in contrast to those from using horizontal branch and RRab stars (all other lines in this figure, including our halo fits). We interpret this as a sign that older stellar populations (RR Lyrae and BHB stars) are distributed in a more concentrated way than giants that should span a wider age spread, thus giving information on the assembly process of the halo: by definition, very early accretion contains only old stellar populations, and these could be more concentrated because only more bound orbits accreted onto the young and lighter proto–Milky Way, whereas more giants originate from later accretions of dwarf galaxies that had prolonged star formation and therefore formed more stars but not more RR Lyrae.

We find that models of stellar accretions support these ideas, especially Font et al. (2008) stating that the larger spread in ages found in the outer halo results from the late assembly of those stars compared to those in the inner halo (Bullock & Johnston 2005; Font et al. 2006). The outer halos in the models studied by Font et al. (2008) tend to show a larger spread in the ages and metallicities of their stellar populations than the inner halos, and this suggests that the outer halo should have a significantly larger fraction of intermediate-age versus old stars than the inner halo.

In our subsequent analysis, we found that even after removing all known prominent substructures, the halo might be more irregular than only being influenced by a flattened halo or flattened inner and outer halo ("dichotomy" of the halo). We find that the south Galactic hemisphere is somewhat more flattened than the north Galactic hemisphere. From calculating the residuals of our best-fit model, we find that there is some deviation of the real halo structure from our best-fit model.