<figure><img src="image_1.png" /></figure>

# The Profile of the Galactic Halo from Pan-STARRS1 $3\pi$ RR Lyrae

Nina Hernitschek¹ ID, Judith G. Cohen¹ ID, Hans-Walter Rix² ID, Branimir Sesar² ID, Nicolas F. Martin²,³ ID, Eugene Magnier⁴ ID, Richard Wainscoat⁴ ID, Nick Kaiser⁴ ID, John L. Tonry⁴ ID, Rolf-Peter Kudritzki⁴, Klaus Hodapp⁴ ID, Ken Chambers⁴ ID, Heather Flewelling⁴ ID, and William Burgett⁵ ID

¹Division of Physics, Mathematics and Astronomy, Caltech, Pasadena, CA 91125, USA; ninah@astro.caltech.edu

²Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

³Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 11 rue de l'Université, F-6700 Strasbourg, France

⁴Institute for Astronomy, University of Hawaii at Manoa, Honolulu, HI 96822, USA

⁵GMTO Corporation (Pasadena), 465 N. Halstead Street, Suite 250, Pasadena, CA 91107, USA

Received 2017 November 1; revised 2018 March 28; accepted 2018 April 18; published 2018 May 21

## Abstract

We characterize the spatial density of the Pan-STARRS1 (PS1) sample of Rrab stars to study the properties of the old Galactic stellar halo. This sample, containing 44,403 sources, spans galactocentric radii of $0.55$ kpc $\le R_{gc} \le 141$ kpc with a distance precision of 3% and thus is able to trace the halo out to larger distances than most previous studies. After excising stars that are attributed to dense regions such as stellar streams, the Galactic disk and bulge, and halo globular clusters, the sample contains ~11,000 sources within $20$ kpc $\le R_{gc} \le 131$ kpc. We then apply forward modeling using Galactic halo profile models with a sample selection function. Specifically, we use ellipsoidal stellar density models $\rho(l, b, R_{gc})$ with a constant and a radius-dependent halo flattening $q(R_{gc})$. Assuming constant flattening $q$, the distribution of the sources is reasonably well fit by a single power law with $n = 4.40_{-0.04}^{+0.05}$ and $q = 0.918_{-0.014}^{+0.016}$ and comparably well fit by an Einasto profile with $n = 9.53_{-0.28}^{+0.27}$, an effective radius $r_{eff} = 1.07 \pm 0.10$ kpc, and a halo flattening of $q = 0.923 \pm 0.007$. If we allow for a radius-dependent flattening $q(R_{gc})$, we find evidence for a distinct flattening of $q \sim 0.8$ of the inner halo at ~25 kpc. Additionally, we find that the south Galactic hemisphere is more flattened than the north Galactic hemisphere. The results of our work are largely consistent with many earlier results (e.g., Watkins et al.; Iorio et al.). We find that the stellar halo, as traced in RR Lyrae stars, exhibits a substantial number of further significant over- and underdensities, even after masking all known overdensities.

**Key words:** Galaxy: halo – stars: variables: RR Lyrae

## 1. Introduction

The Milky Way's extended stellar halo contains only a small fraction (~1%) of the Galaxy's stars but is an important diagnostic of the Milky Way's formation, dark matter distribution, and mass.

The stellar halo shows great complexity in its spatial structure, with abundant globular clusters, dwarf galaxies, and stellar streams. This makes it difficult to dissect with local spectroscopic or photometric data. While the radial density profile can be derived from data of a limited number of sightlines through the Galaxy, a sensible description of the overall stellar halo shape requires nearly complete coverage of the sky.

As stellar halos formed from disrupted satellites and still show signs of their accretion history in the form of overdensities such as streams, they are central to studies on galaxy formation such as the hierarchical galaxy formation in the ΛCMD model. The spatial distribution, as well as kinematics, metallicities, and thus ages of halo stars, enables us to get information on those merger processes, as well as to compare them to simulations from theoretical models.

Many studies were carried out within the past 50 yr to map the Galactic halo, and those studies often took advantage of RR Lyrae stars as reliable halo tracers. These old and metal-poor pulsators are ideal for this task, as they can be selected with a high purity, thus showing only very little contamination from other populations of the Milky Way. Furthermore, RRab are luminous variable stars pulsating in the fundamental mode that obey a well-defined period-luminosity relation, albeit with a

small dependence on metallicity. Thus, the mean luminosity of an RRab variable, and hence its distance, can be determined with knowledge of the light curve only. RRab stars were used by many previous studies, including those of Hawkins (1984), Saha (1984), Wetterer et al. (1996), Ivezić et al. (2000), Vivas & Zinn (2006), Jurić et al. (2008), Catelan (2009), Watkins et al. (2009), de Jong et al. (2010), Sesar et al. (2010, 2011), Deason et al. (2011), Akhter et al. (2012), Drake et al. (2014), Torrealba et al. (2015), Cohen et al. (2015), Xue et al. (2015), Soszyński et al. (2016), Bland-Hawthorn & Gerhard (2016), Iorio et al. (2017), and Cohen et al. (2017).

The key to using RRab to explore the Galactic halo is having a reliable list of RRab variables selected from a suitable multiepoch imaging survey covering a wide distance range and as much of the sky as possible. Recently the inner halo out to ~30 kpc was explored by a sample of ~5000 RRab generated from a recalibration of the LINEAR catalog by Sesar et al. (2013b) and most recently by Iorio et al. (2017) using a sample selected from the combination of the Gaia Data Release 1 (Gaia Collaboration et al. 2016) and the Two Mass All Sky Survey (2MASS; Skrutskie et al. 2006).

Drake et al. (2014) used the Catalina Real-Time Transient Surveys DR1 to select a sample of 47,000 periodic variables, of which 16,797 are RR Lyrae, and the bulk of them are at $R_{gc} < 40$ kpc. In total, the Catalina Surveys RR Lyrae Data Release 1⁶ (Drake et al. 2013a, 2013b, 2014;

<footnote n="6">⁶ http://nesssi.cacr.caltech.edu/DataRelease/RRL.html</footnote>
Torrealba et al. 2015; Drake et al. 2017) contains 43,599 RR Lyrae, of which 32,980 are RRab stars.

To reach larger distances with larger samples was very difficult in the past. One approach was to use brighter tracer stars, usually K giants and usually selected from the Sloan Digital Sky Survey (SDSS), but with larger distance uncertainties and only modest sample sizes (see, e.g., Xue et al. 2015, who probe the Galactic halo out to 80 kpc using 1757 stars from the SEGUE K-giant Survey). There have also been efforts to reach the outer halo using blue horizontal branch (BHB) stars, which to first order have a fixed luminosity similar to that of RRab (see, e.g., Deason et al. 2014), but these run into problems of confusion with much more numerous blue stragglers at the same apparent magnitude and with quasars. Prior to the present work, perhaps the most successful attempt to probe the density distribution in the outer Galactic halo was by Cohen et al. (2017), reaching out to above 100 kpc, with a small (~450) sample of RRab stars selected from the Palomar Transient Facility (PTF) database.

Here, we overcome these difficulties by using a selection of RRab from the PS1 survey, which covers the entire northern sky to a limiting magnitude such that detection of RRab out to more than 100 kpc is not difficult. Hernitschek et al. (2016) and Sesar et al. (2017) exploited the PS1 survey to create a sample of RRab that reaches far into the outer halo, which is very large (44,403 RRab), with known high purity and completeness. The details of the machine-learning techniques that were used to select this sample and the assessment of its purity and completeness as a function of distance are described in Hernitschek et al. (2016) and Sesar et al. (2017).

In this paper we exploit the PS1 RRab sample to study the Milky Way halo out to distances in excess of 100 kpc.

We develop and apply a rigorous density modeling approach for Galactic photometric surveys that enables investigation of the structure of the Galactic halo as traced by RR Lyrae stars from 20 kpc to more than 100 kpc. We fit models that characterize the radial density and flattening of the Milky Way's stellar halo, while accounting for the complex selection function resulting from both the survey itself and the selection of sources within the survey data.

In Section 2, we lay out the properties of the PS1 RRab stars. In Section 3 we present the method of fitting a series of parameterized models to the RRab stars while considering a selection function. This step is key to obtaining accurate radial profiles. In the following, first two types of parameterized models for the radial stellar density are shown in Section 3.1, followed by the description of the selection function in Section 3.2 and our approach to constrain the model parameters in Section 3.3. A test method, relying on mock data, is shown in Section 3.4. In Section 4, we present the results for the profile and flattening of the Milky Way's stellar halo, as well as findings of previously unknown halo overdensities. In Section 5, we discuss results and methodology and compare them to work by others. Finally, Section 6 summarizes the paper.

## 2. RR Lyrae Stars from the PS1 Survey

Our analysis is based on a sample of highly likely RRab stars, as selected by Sesar et al. (2017) from the Pan-STARRS1 $3\pi$ survey. In this section, we describe the pertinent properties of the PS1 $3\pi$ survey and the RR Lyrae light curves obtained, and we recapitulate briefly the process of selecting the likely

RRab, as laid out in Sesar et al. (2017). We also briefly characterize the obtained candidate sample.

The Pan-STARRS1 (PS1) survey (Kaiser et al. 2010) collected multiepoch, multicolor observations undertaking a number of surveys, among which the PS1 $3\pi$ survey (Chambers et al. 2016) is currently the largest. It has observed the entire sky north of decl. $-30^\circ$ in five filter bands ($g_{P1}$, $r_{P1}$, $i_{P1}$, $z_{P1}$, $y_{P1}$) with a $5\sigma$ single epoch depth of about 22.0, 22.0, 21.9, 21.0, and 19.8 mag in $g_{P1}$, $r_{P1}$, $i_{P1}$, $z_{P1}$, and $y_{P1}$, respectively (Stubbs et al. 2010; Tonry et al. 2012).

Starting with a sample of more than $1.1 \times 10^9$ PS1 $3\pi$ sources, Hernitschek et al. (2016) and Sesar et al. (2017) subsequently selected a sample of 44,403 likely RRab stars, of which ~17,500 are at $R_{gc} \ge 20$ kpc, by applying machine-learning techniques based on light-curve characteristics. RRab stars are the most common type of RR Lyrae, making up ~91% of all observed RR Lyrae (Smith 2004) and displaying the steep rises in brightness typical of RR Lyrae.

The identification of the RRab stars is highly effective, and the sample of RRab stars is pure (90%) and complete ($\ge 80\%$ at 80 kpc) at high galactic latitudes. The distance estimates are precise to 3%, based on newly derived period-luminosity relations for the optical/near-infrared PS1 bands (Sesar et al. 2017). Overall, this results in the widest (3/4 of the sky) and deepest (reaching $>120$ kpc) sample of RR Lyrae stars to date, allowing us to observe them globally across the Milky Way. Out of these sources, 1093 exist beyond a galactocentric distance of 80 kpc, and 238 exist beyond 100 kpc.

In the subsequent analysis, we refer to this sample (Sesar et al. 2017) as “RRab stars.”

The left panels of Figure 1 show the source density of the PS1 sample of RRab stars for different distance bins $0\text{ kpc} < D \le 20\text{ kpc}$, $20\text{ kpc} < D \le 50\text{ kpc}$, and $50\text{ kpc} < D \le 120\text{ kpc}$. The right panels of the same figure show the sample after a cleaning to remove overdensities was applied; the details of this cleaning are described later.

Figure 2 is based on the same data but shown in the Cartesian reference frame $(X, Y, Z)$ for an easier comparison with subsequent plots of halo models, as well as to highlight the individual effects of removing certain overdensities.

While the sample covers the entire sky above a decl. $\delta > -30^\circ$, which enables a view of halo substructure like the Sagittarius stream (Hernitschek et al. 2017), in this paper we focus on stars away from the Galactic plane and center, as well as away from known large overdensities like the Sagittarius stream. Details of the process of removing these overdensities are given in Section 3.2.

## 3. Density Fitting

In this section we lay out a forward-modeling approach to describe the spatial distribution of the stellar halo using a set of flexible but ultimately smooth and symmetric functions.

We presume that the stellar halo distribution can be sensibly approximated by a spheroidal distribution with a parameterized radial profile. Similar approaches were carried out by, e.g., Sesar et al. (2013b), Xue et al. (2015), Cohen et al. (2015), and Iorio et al. (2017), but all with either a smaller sample size than in our analysis or probing a smaller distance range.

The number of halo parameters depends on the complexity of the model assumed for the stellar halo distribution. The mathematics of this approach essentially follows Bovy et al. (2012) and Rix & Bovy (2013).
<figure><img src="image_2.png" /><figcaption>cleaned RRab sample</figcaption></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Figure 1.** Density of the uncleaned and cleaned PS1 $3\pi$ sample of RRab stars, shown in Galactic coordinates $(l, b)$ for different heliocentric distance bins. The logarithmic source number density is given within $5\deg^2$ wide bins, in units of deg$^{-2}$. This bin size was chosen to reduce Poisson noise. White cells are empty, and dark-blue cells have 1 source per deg$^2$. Starting from a sample of 44,403 sources (Sesar et al. 2017), containing overdensities like globular clusters and streams and affected by sample incompleteness near the Galactic plane and apocenter (here shown in the left column as “uncleaned sample”), we construct a sample of 11,025 sources outside of such known overdensities. To do so, we apply the selection cuts described in Section 3.2, to geometrically excise such overdensities. The largest overdensities removed are the Sagittarius stream (we remove sources associated with the Sgr stream according to Hernitschek et al. 2017), as well as the thick disk (we remove sources within $|b| < 10^\circ$) and close to the Galactic center and the bulge (we remove sources within $R_{gc} \le 20$ kpc). Showing the source density in three different distance bins shows major overdensities, as well as how excising such overdense regions affects the cleaned sample.

A number of very different models have been proposed for the density profile of the stellar halo. We denote the spatial number density here as $\rho_{\text{RRL}}(l, b, D)$ and the general form of the models as $\rho_{\text{RRL}}(D|\theta)$, where $\theta$ are the model parameters

(see Section 3.1) and $D = (l, b, D)$ are the observables with Galactic coordinates and the heliocentric distance $D$.

An approach for fitting the spatial density profiles of the RRab sample must account for the fact that the observed star
<figure><img src="image_9.png" /><figcaption>cleaned RRab sample</figcaption></figure>

<figure><img src="image_10.png" /></figure>

**Figure 2.** PS1 3π sample of RRab stars, shown in the Cartesian reference frame (X, Y, Z) as given in Equation (1). This reference frame is centered at the Galactic center, the Galactic disk is placed in the (X, Y) plane, the X-axis is pointing to the Sun, and the Z-axis is pointing to the north Galactic pole. The logarithmic source number density in each projection is given for 1 kpc² wide bins. Starting from a sample of 44,403 sources (Sesar et al. 2017), containing overdensities like globular clusters and streams and affected by sample incompleteness near the Galactic plane and apocenter (here shown in the top panel as “uncleaned sample”), we construct a sample of 11,025 sources outside of such known overdensities. To do so, we apply the selection cuts described in Section 3.2, to geometrically excise such overdensities. The largest overdensities removed are the Sagittarius stream (we remove sources associated with the Sgr stream according to Hernitschek et al. 2017), as well as the thick disk (we remove sources within |b| < 10°) and close to the Galactic center and the bulge (we remove sources within R<sub>gc</sub> ≤ 20 kpc). The effects of removing those sources are clearly visible in the bottom panels and are each labeled. The dashed circle here represents the 20 kpc cut. Sources within the circle but farther away are seen owing to projection effects; the distinctly higher density just after 20 kpc shows the stars that are no longer affected by this distance cut.
counts do not reflect the underlying stellar distribution but are strongly shaped by selection effects both from the survey itself and from selection cuts we chose while preparing the sample. We denote the spatial selection function as S(l, b, D) (see Section 3.2).

To properly take all of these effects into account, we need to use forward modeling. In what follows we fit stellar density models to the data by generating the expected observed distribution of stars in the RRab sample, based on our model for the selection function and the halo density models. This predicted distribution is then automatically compared to the observed star counts to calculate the likelihood of the observed RRab star counts.

### 3.1. Stellar Density Models

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 $\theta$ and that the density $\rho_{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^{-3}$:

$$ 
\begin{aligned}
X &= R_{\odot} - D \cos l \cos b \\
Y &= -D \sin l \cos b \\
Z &= D \sin b.
\end{aligned}
\qquad (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_{\odot}$ 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_{\odot}$ 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_{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 $\rho_{halo}$ is widely used (e.g., Sesar et al. 2013b) to describe the distribution of the halo stars:

$$ \rho_{\text{halo}}(X, Y, Z) = \rho_{\odot\text{RRL}}(R_{\odot}/r_q)^n. \qquad (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 $\theta = (n, q)$. Here $\rho_{\odot\text{RRL}}$ is the number density of RR Lyrae at the position of the Sun, $R_{\odot}$ 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 $\rho_{\odot\text{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:

$$ \rho_{\text{halo}}(X, Y, Z) = \begin{cases} \rho_{\odot\text{RRL}}(R_{\odot}/r_q)^{n_{\text{inner}}}, & \text{if } r_q \le r_{\text{break}} \\ \rho_{\odot\text{RRL}}r_{\text{break}}^{n_{\text{outer}}-n_{\text{inner}}}(R_{\odot}/r_q)^{n_{\text{outer}}}, & \text{else.} \end{cases} \qquad (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

$$ \gamma(r) = - \frac{d \ln \rho(r)}{d \ln r} \propto r^{\alpha}, \qquad (4) $$

where the steepness of the Einasto profile, $\alpha$, changes continuously as a function of the effective radius $r_{eff}$,

$$ \alpha = - \frac{d_n}{n} \left( \frac{r_q}{r_{\text{eff}}} \right)^{1/n}. \qquad (5) $$

This can be rearranged to

$$ \rho_{\text{halo}}(r_q) \equiv \rho_0 \exp \left\{ -d_n \left[ (r_q / r_{\text{eff}})^{1/n} - 1 \right] \right\}, \qquad (6) $$

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

**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}$.

<sup>a</sup> The center of the removed overdensity.

<sup>b</sup> These sources are also removed by other cuts.

<sup>c</sup> This RRab is a member of the Orphan stream (Sesar et al. 2013a).

The free parameters of an Einasto profile with a constant flattening $q$ are $\theta = (r_{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_{gc})$ as

$$q(R_{gc}) = q_{\infty} - (q_{\infty} - q_{0})\exp\left(1 - \frac{\sqrt{R_{gc}^{2} + r_{0}^{2}}}{r_{0}}\right), \quad (7)$$

with $q_0$ being the flattening at the center, $q_\infty$ being the flattening at large galactocentric radii, and $r_0$ being the exponential scale radius over which the change of flattening occurs.

Thus, the flattening $q$ now varies from $q_0$ at the center to the asymptotic value $q_\infty$ at large radii, and the variation is tuned by the exponential scale length $r_0$.

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_0$, $q_\infty$, and $r_0$.

### 3.2. Selection Function

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 $S(l, b, D)$ is binary [0, 1] so that $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\pi$ 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| < 10^\circ$), or too close to the Galactic center ($R_{gc} \le 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.

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 $S(l, b, D)$ is thus composed of

$$S(l, b, D) = S_{\text{RRL}}(l, b, D) \times S_{\text{area}}(l, b, D), \quad (8)$$

where $S_{\text{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 $S_{\text{area}}(l, b, D)$ describes area cuts to exclude overdensities.

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

$$S_{\text{RRL}}(l, b, D) = \begin{cases} 1, & \text{if } \delta > -30^{\circ} \text{ and } D_{\text{min}} < D < D_{\text{max}} \\ 0, & \text{else.} \end{cases} \quad (9)$$

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

$$S_{\text{bulge,disc}}(l, b, D) = \begin{cases} 1, & \text{if } |b| \ge 10^{\circ} \text{ and } R_{\text{gc}} \ge 20 \text{ kpc} \\ 0, & \text{else.} \end{cases} \quad (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_{\text{sgr}}$ (Hernitschek et al. 2017, see Equation (11) therein). We excise sources with $p_{\text{sgr}} > 0.2$ as members of the Sgr stream, leading to a selection function of

$$S_{\text{sgr}}(l, b, D) = \begin{cases} 1, & \text{if } p_{\text{sgr}}(l, b, D) < 0.2 \\ 0, & \text{else.} \end{cases} \quad (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 $S_{\text{other 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 \lesssim 20$ mag, after which it steeply drops to zero at $r_F \sim 21.5$ mag. Writing $r_F$ as $r_F(D)$, the selection function characterizing the distance-dependent completeness is

$$S_c(r_F) = L - \frac{L}{1 + \exp(-k(r_F - x_0))}, \quad (12)$$

with (Sesar et al. 2017)

$$L = 0.91 \quad (13)$$

$$k = 4.0 \quad (14)$$

$$x_0 = 20.6 \quad (15)$$

$$r_F = 2.05 \log(D) + 11. \quad (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, $S_p(D)$.

We end up with a selection function

$$\begin{aligned} S(l, b, D) &= S_{\text{RRL}}(l, b, D) \times S_c(D) \times S_p(D) \times S_{\text{area}}(l, b, D) \\ &= S_{\text{RRL}}(l, b, D) \times S_c(D) \\ & \quad \times S_{\text{bulge,disc}}(l, b, D) \\ & \quad \times S_{\text{sgr}}(l, b, D) \times S_{\text{other overdensities}}(l, b, D). \end{aligned} \quad (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,<fnref n="7" /> 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\pi$. 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 $S_{\text{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.

### 3.3. Constraining Model Parameters

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

<footnote n="7">⁷ http://edd.ifa.hawaii.edu/describe_columns.php?table=kmcc</footnote>
mock halo power law profile with varying flattening $q(R_{gc})$

<figure><img src="image_11.png" /><figcaption>**Figure 3.** One- and two-dimensional projections of the posterior probability distributions (pdf) of parameters ($r_0, q_0, q_\infty, n$) of the power law with varying flattening $q(R_{gc})$ (Equation (7)) fitted to a mock sample, used to test the methodology for fitting the halo density profile. The blue lines and squares mark the maximum likely value of each parameter. The best-fit parameters are given along with their $1\sigma$ intervals in the top right part of the figure. The parameters used for generating the mock sample are indicated by dark red lines and squares and also given in the right part of the figure.</figcaption></figure>

The normalized unmarginalized log likelihood for the ith star with the observables $D_i$ is then

$$ \ln p(D_i|\theta) = \frac{\rho_{RRL}(D_i|\theta)|J|S(l_i, b_i, D_i)}{\iint \rho_{RRL}(l, b, D|\theta)|J|S(l, b, D)dldbdD}, \quad (18) $$

where the normalization integral is over the observed volume. The Jacobian term $|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 $\theta$ of the halo model, given the full data $D$ and a prior $p(\theta)$, $\ln p(\theta|D) = \ln p(D|\theta) + \ln p(\theta)$, with

$$ \ln p(D|\theta) = \sum_i \ln p(D_i|\theta) \quad (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
<figure><img src="image_12.png" /><figcaption>**Figure 4.** Comparison between the observed distance distribution of the cleaned samples and the predicted distributions by the best-fit models, with the number density shown in a log plot. The black histogram shows the galactocentric distance distribution of our cleaned sample of 11,025 RRab stars, whereas the gray histogram gives the distance distribution of the full data set of 44,403 RRab stars from Sesar et al. (2017). Removed overdensities are highlighted with dashed lines and are listed in Table 1. The overplotted solid lines represent the best-fit model for each of the five halo profiles. As a result of the selection function, these models do not follow a straight line in the log plot, but drop much more rapidly, especially beyond a galactocentric distance of 80 kpc. For comparison, dashed lines, in the same color as the solid lines, represent each $\rho_{halo} \times S(l, b, D)$, where S is the selection function as given in Equation (17). We see that each our five models can fit the distance distribution properly, and our assumption about the selection function S represents the true selection effects and overdensity cuts. The best-fit parameters for each of the models are given in Table 2.</figcaption></figure>

<table><caption>**Table 2**
Best-fitting Halo Models</caption><thead><tr><th>Density Model</th><th>Best-fit Parameters</th><th>ln(L<sub>max</sub>)</th><th>BIC</th><th>ΔBIC</th></tr></thead><tbody><tr><td>Power-law model</td><td>q = 0.918<sup>+0.016</sup><sub>-0.014</sub>, n = 4.40<sup>+0.05</sup><sub>-0.04</sub></td><td>-157625</td><td>315269</td><td>203</td></tr><tr><td>BPL model</td><td>r<sub>break</sub> = 38.7<sup>+0.69</sup><sub>-0.58</sub>, q = 0.908<sup>+0.008</sup><sub>-0.006</sub>,<br>n<sub>inner</sub> = 4.97<sup>+0.02</sup><sub>-0.05</sub>, n<sub>outer</sub> = 3.93<sup>+0.05</sup><sub>-0.04</sub></td><td>-214222</td><td>428464</td><td>113398</td></tr><tr><td>Einasto profile</td><td>r<sub>eff</sub> = 1.07 ± 0.10 kpc,<br>q = 0.923 ± 0.007, n = 9.53<sup>+0.27</sup><sub>-0.28</sub></td><td>-157685</td><td>315388</td><td>322</td></tr><tr><td>Power-law model with q(R<sub>gc</sub>)</td><td>r<sub>0</sub> = 25.0<sup>+1.8</sup><sub>-1.7</sub> kpc, q<sub>0</sub> = 0.773<sup>+0.017</sup><sub>-0.016</sub>,<br>q<sub>∞</sub> = 0.998<sup>+0.002</sup><sub>-0.001</sub>, n = 4.61 ± 0.03</td><td>-157524</td><td>315066</td><td>0</td></tr><tr><td>Einasto profile with q(R<sub>gc</sub>)</td><td>r<sub>0</sub> = 26.7<sup>+2.2</sup><sub>-2.0</sub> kpc, q<sub>0</sub> = 0.779 ± 0.018,<br>q<sub>∞</sub> = 0.998<sup>+0.001</sup><sub>-0.002</sub>, r<sub>eff</sub> = 1.04<sup>+0.25</sup><sub>-0.13</sub> kpc,<br>n = 8.78<sup>+0.33</sup><sub>-0.30</sub></td><td>-157582</td><td>315182</td><td>116</td></tr></tbody></table>

**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\sigma$ uncertainties estimated as the 15.87th and 84.13th percentiles, the maximum log likelihood $\ln(L_{max})$, and the BIC. $\Delta$BIC gives the difference between the BIC of the best-fit model, the power-law model with $q(R_{gc})$, and the model used.

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\sigma$.

The calculation of the normalization integral in Equation (18) is complicated by the presence of the selection function 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 multi-dimensional quadrature methods in Python are not able to
power law profile

<figure><img src="image_13.png" /><figcaption>**Figure 5.** One- and two-dimensional projections of the posterior probability distributions of parameters (q, n) of the power law (Equation (2)) fitted to the cleaned sample. The blue lines and squares mark the median value of each parameter. The best-fit parameters are given along with their 1σ intervals in the top right part of the figure. Both the power-law index n and the flattening parameter q show an almost Gaussian distribution with no covariance.</figcaption></figure>

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(\theta)$ for each of the five following cases:

Power-law model:

$$ \ln p(\theta) = \begin{aligned} & \text{Uniform}(1.0 < n < 6.0) \\ & + \text{Uniform}(0.1 < q < 1.0). \end{aligned} \qquad (20) $$

BPL model:

$$ \ln p(\theta) = \begin{aligned} & \text{Uniform}(1.0 < n < 6.0) \\ & + \text{Uniform}(0.1 < q < 1.0) \\ & + \text{Uniform}(\log(R_{\text{min}}) < \log(r_{\text{break}}) < \log(R_{\text{max}})), \end{aligned} \qquad (21) $$

where $R_{\text{min}}$, $R_{\text{max}}$ give the galactocentric distance range available in the sample.

Einasto profile:

$$ \ln p(\theta) = \begin{aligned} & \text{Uniform}(\log(0.01) < \log(r_{\text{eff}}) < \log(50)) \\ & + \text{Uniform}(\log(0.01) < \log(r_0) < \log(50)) \\ & + \text{Uniform}(0.5 < n < 20.0) \\ & + \text{Uniform}(0.1 < q < 1.0). \end{aligned} \qquad (22) $$

Power-law model with $q(R_{\text{gc}})$:

$$ \ln p(\theta) = \begin{aligned} & \text{Uniform}(\log(0.01) < \log(r_0) < \log(50)) \\ & + \text{Uniform}(1.0 < n < 5.0) \\ & + \text{Uniform}(0.1 < q_0 < 1.0) \\ & + \text{Uniform}(0.1 < q_{\infty} < 1.0). \end{aligned} \qquad (23) $$

Einasto profile with $q(R_{\text{gc}})$:

$$ \ln p(\theta) = \begin{aligned} & \text{Uniform}(\log(0.01) < \log(r_{\text{eff}}) < \log(50)) \\ & + \text{Uniform}(\log(0.01) < \log(r_0) < \log(50)) \\ & + \text{Uniform}(\log(0.01) < \log(r_0) < \log(50)) \\ & + \text{Uniform}(0.5 < n < 20.0) \\ & + \text{Uniform}(0.1 < q_0 < 1.0) \\ & + \text{Uniform}(0.1 < q_{\infty} < 1.0). \end{aligned} \qquad (24) $$

### 3.4. Fitting Tests on Mock Data

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_{\text{gc}})$, or Einasto profile with $q(R_{\text{gc}})$. We then applied a 3% distance uncertainty, superimposed the sample with faint and far Gaussian blobs
<figure><img src="image_14.png" /><figcaption>**Figure 6.** One- and two-dimensional projections of the posterior probability distributions of parameters ($q, r_{break}, n_{inner}, n_{outer}$) of the power law (Equation (3)) fitted to the cleaned sample. The blue lines and squares mark the median value of each parameter. The best-fit parameters are given along with their $1\sigma$ intervals in the top right part of the figure. The parameters show an almost Gaussian distribution with no covariance.</figcaption></figure>

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.
<figure><img src="image_15.png" /><figcaption>**Figure 7.** One- and two-dimensional projections of the posterior probability distributions of parameters ($r_{eff}$, $q$, $n$) of the Einasto profile (Equation (6)) fitted to the cleaned sample. The blue lines and squares mark the median value of each parameter. The best-fit parameters are given along with their $1\sigma$ intervals in the top right part of the figure. As is the case 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}$.</figcaption></figure>

## 4. Results

We now present the results of applying the modeling from Section 3.1 to the cleaned sample of RRab stars as described in Section 3.2. We fitted the simplest model, a power law, to our data, as well as the Einasto model, to allow easy comparison with density profiles obtained from N-body simulations (Diemand et al. 2004; Navarro et al. 2004; Graham et al. 2006; Merritt et al. 2006). Both models are fitted with a constant halo-flattening parameter $q$, as well as with a distance-dependent flattening $q(R_{gc})$. We illustrate these results in three ways: (1) by showing the predicted distribution by the best-fit models, (2) by showing the joint posterior distribution functions of the halo model parameters of each model, and (3) by comparing the models using the Bayesian information criterion (BIC).

First, we discuss the result of fitting the complete cleaned $3\pi$ sample, in order to explore the broad trends in spatial structure. Subsequently, we split the sample into two hemispheres, as well as into relatively broad $\Delta l = 30^\circ$, $\Delta b = 60^\circ$ bins, and map the local halo structure. Finally, we calculate and analyze the residuals of the best-fit model.

### 4.1. Best-fit Model Parameters

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
power law profile with varying flattening $q(R_{gc})$

<figure><img src="image_16.png" /><figcaption>**Figure 8.** One- and two-dimensional projections of the posterior probability distributions (pdf) of parameters ($r_0, q_0, q_\infty, n$) of the power law with varying flattening $q(R_{gc})$ (Equation (7)) fitted to the cleaned sample. The blue lines and squares mark the maximum likely value of each parameter. The best-fit parameters are given along with their $1\sigma$ intervals in the top right part of the figure. The fitting parameters show a covariance, and the pdf is strongly distorted from a Gaussian distribution, including local maxima in the distribution of $r_0$.</figcaption></figure>

8917 sources. Figure 4 compares the observed number density of RR Lyrae stars with the density predicted by best-fit models.

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\sigma$ uncertainties estimated as the 15.87th and 84.13th percentiles, and the maximum log likelihood $\ln(\mathcal{L}_{\text{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.

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_{\text{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
<figure><img src="image_17.png" /><figcaption>Einasto profile with varying flattening $q(R_{gc})$</figcaption></figure>

**Figure 9.** One- and two-dimensional projections of the posterior probability distributions (pdf) of parameters ($r_{eff}$, $r_0$, $q_0$, $q_\infty$, $n$) of the Einasto profile with varying flattening $q(R_{gc})$ (Equation (7)) fitted to the cleaned sample. The blue lines and squares mark the maximum likely value of each parameter. The best-fit parameters are given along with their $1\sigma$ intervals in the top right part of the figure. The fitting parameters $r_0$, $q_0$, and $q_\infty$ show a (partially strong) covariance, and the pdf is strongly distorted from a Gaussian distribution, including local maxima in the distribution of $r_0$, $q_0$, and $q_\infty$. The best-fit model and the pdf lead to $q_0 \sim q_\infty$, thus being quite similar to the Einasto profile with a constant flattening.

parameters $r_0$, $q_0$, and $q_\infty$ 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_{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_{inner} = 4.97_{-0.05}^{+0.02}$.

and $n_{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 \pm 0.10$ kpc, and a halo flattening of $q_0 = 0.923 \pm 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}$ kpc, $n = 4.61 \pm 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
<table><caption>**Table 3**
Best-fitting Halo Models for Each Hemisphere</caption><thead><tr><th>Density Model</th><th>Best-fit Parameters North Galactic Hemisphere</th><th>Best-fit Parameters South Galactic Hemisphere</th></tr></thead><tbody><tr><td>Power-law model</td><td>q = 0.925<sup>+0.010</sup><sub>-0.009</sub>, n = 4.36 ± 0.03</td><td>q = 0.852<sup>+0.010</sup><sub>-0.011</sub>, n = 4.40 ± 0.04</td></tr><tr><td>Einasto profile</td><td>r<sub>eff</sub> = 1.11<sup>+0.09</sup><sub>-0.10</sub> kpc,<br>q = 0.934<sup>+0.009</sup><sub>-0.010</sub>, n = 9.59<sup>+0.30</sup><sub>-0.26</sub></td><td>r<sub>eff</sub> = 1.18<sup>+0.09</sup><sub>-0.11</sub> kpc<br>q = 0.851<sup>+0.013</sup><sub>-0.011</sub>, n = 9.10<sup>+0.31</sup><sub>-0.28</sub></td></tr><tr><td>Power-law model with q(R<sub>gc</sub>)</td><td>r<sub>0</sub> = 29.2 ± 4.4 kpc, q<sub>0</sub> = 0.831<sup>+0.031</sup><sub>-0.017</sub>,<br>q<sub>∞</sub> = 0.997<sup>+0.004</sup><sub>-0.001</sub>, n = 4.53<sup>+0.04</sup><sub>-0.06</sub></td><td>r<sub>0</sub> = 18.8<sup>+1.4</sup><sub>-1.6</sub> kpc, q<sub>0</sub> = 0.515<sup>+0.027</sup><sub>-0.058</sub><br>q<sub>∞</sub> = 0.998<sup>+0.004</sup><sub>-0.002</sub>, n = 4.88<sup>+0.06</sup><sub>-0.03</sub></td></tr><tr><td>Einasto profile with q(R<sub>gc</sub>)</td><td>r<sub>0</sub> = 31.9<sup>+3.9</sup><sub>-3.1</sub> kpc, q<sub>0</sub> = 0.837<sup>+0.036</sup><sub>-0.017</sub>,<br>q<sub>∞</sub> = 0.998<sup>+0.001</sup><sub>-0.003</sub>, r<sub>eff</sub> = 1.00<sup>+0.3</sup><sub>-0.11</sub> kpc,<br>n = 9.07<sup>+0.35</sup><sub>-0.28</sub></td><td>r<sub>0</sub> = 20.9<sup>+2.3</sup><sub>-2.0</sub> kpc, q<sub>0</sub> = 0.545<sup>+0.040</sup><sub>-0.066</sub><br>q<sub>∞</sub> = 0.998<sup>+0.001</sup><sub>-0.005</sub>, r<sub>eff</sub> = 1.14<sup>+0.16</sup><sub>-0.12</sub> kpc<br>n = 7.57<sup>+0.40</sup><sub>-0.23</sub></td></tr></tbody></table>

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\sigma$ uncertainties estimated as the 15.87th and 84.13th percentiles.

Einasto profile with $q(R_{gc})$ are $r_0 = 26.7_{-2.0}^{+2.2}$ kpc, $q_0 = 0.779 \pm 0.018$, $q_\infty = 0.998_{-0.002}^{+0.001}$, $r_{eff} = 1.04_{-0.13}^{+0.25}$ kpc, and $n = 8.78_{-0.30}^{+0.33}$.

We find here $q_0 < 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\sigma$ uncertainties are summarized in Table 2.

Our results confirm that if a varying flattening is assumed, the halo profile has an $r_0$ 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.

## 4.2. Comparing Models

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

$$ \text{BIC} = \text{dim}(\theta)\ln(N) - 2\ln(\mathcal{L}_{\text{max}}), \quad (25) $$

where $\theta$ are the model parameters, $N$ is the number of objects in the sample, and $\mathcal{L}_{\text{max}}$ is the maximum likelihood, where we defined the likelihood function in Equation (19) as $\ln p(\mathcal{D}|\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 $\Delta$BIC > 10 indicates very strong evidence against the model with the higher BIC.

## 4.3. Local Halo Properties

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\pi$ 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
<figure><img src="image_18.png" /><figcaption>**Figure 10.** Angular source number density for the cleaned RRab sample, given per ($\Delta l = 30^{\circ}$) $\times$ ($\Delta b = 60^{\circ}$) bin. This binning is used to fit for the local halo properties. The number density is color-coded, as well as given in numbers. Empty cells are outside the survey footprint. Away from the Galactic equator, the angular number density drops as the spanned area decreases. We find that the number of sources is increased within $30^{\circ} < l < 90^{\circ}$, $|b| < 30^{\circ}$; stars not fully excised from the Galactic bulge and the crossing Sagittarius stream account for that. Also around $240^{\circ} < l < 330^{\circ}$, $30^{\circ} < b < 90^{\circ}$ we find an increase of sources, as we can remove most but not all stars from the Sgr stream by setting geometric cuts on their angle above the plane of the stream. A significant increase near the Galactic center, where sources would fall into the bin $0^{\circ} < l < 30^{\circ}$, $-30^{\circ} < b < 30^{\circ}$, is not found, as we remove everything within $R_{gc} \le 20$ kpc and $|b| < 10^{\circ}$ as well.</figcaption></figure>

Galactic hemisphere and also compare well with the fit for the complete halo. When taking a look at the flattening-related parameters, $q$, $q_0$, $q_{\infty}$, $r_{\text{eff}}$, we find that for models with constant flattening (both the power-law and Einasto profile models) $q_{\text{south}} < q < q_{\text{north}}$. In the case of models with $q(R_{gc})$, we find that the value of parameter $q_0$ is smaller for the south than for the north hemisphere, $q_{0,\text{south}} < q_0 < q_{0,\text{north}}$, whereas the value of the parameter $q_{\infty}$ is similar for both hemispheres. Furthermore, we find that $r_{0,\text{north}} > r_{0,\text{south}} > r_0$ for both the power law with $q(R_{gc})$ and the Einasto profile with $q(R_{gc})$.

The results of finding $q_{\text{south}} < q < q_{\text{north}}$ for models with constant flattening and $q_{0,\text{south}} < q_0 < q_{0,\text{north}}$, $q_{\infty,\text{south}} \sim q_{\infty,\text{north}} \sim q_{\infty,\text{south}}$, $r_{0,\text{north}} > r_{0,\text{south}} > r_0$ in the case of a radius-dependent flattening are consistent: by definition of $q(R_{gc})$ (Equation (7)), $q_0$ is the flattening at center, $q_{\infty}$ is the flattening at large galactocentric radii, and $r_0$ is the exponential scale radius over which the change of flattening occurs. A larger $r_0$ means that the flattening of the inner halo, where we find $q_0 < q_{\infty}$, is in force out to a larger radius than for a smaller $r_{\text{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 ($\Delta l = 30^{\circ}$) $\times$ ($\Delta b = 60^{\circ}$), 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 ($\Delta l = 30^{\circ}$) $\times$ ($\Delta b = 60^{\circ}$) bin, is shown in Figure 10.

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\sigma$ uncertainties.

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^{\circ} < l < 270^{\circ}$, $-30^{\circ} < b < 30^{\circ}$ and $30^{\circ} < l < 60^{\circ}$, $-90^{\circ} < b < -30^{\circ}$ (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^{\circ} \times 60^{\circ}$ 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\sigma$ intervals given along with the best-fit parameters in the top right part of the figure.
<figure><img src="image_19.png" /><figcaption>PS1 3π RRab stars within $R_{gc} > 20$ kpc, $30^{\circ} < l < 60^{\circ}$, $-90^{\circ} < b < -30^{\circ}$ with best-fit halo profiles</figcaption></figure>

<figure><img src="image_20.png" /></figure>

**Figure 11.** Comparison between the observed distance distribution of the cleaned samples and the predicted distributions by the best-fit models for $30^{\circ} \times 60^{\circ}$ patches on the sky, with the number density shown in a log plot. We find that the fitting procedure also works reliably with small pieces on the sky. The black histogram shows the galactocentric distance distribution of our cleaned sample of RRab stars within the given patch on the sky, whereas the gray histogram gives the distance distribution of all RRab stars within the given patch on the sky. Removed overdensities are highlighted with dashed lines and are listed in Table 1. The overplotted solid lines represent the best-fit model for each of the four halo profiles. As a result of the selection function, these models do not follow a straight line in the log plot, but drop much more rapidly, especially beyond a galactocentric distance of 80 kpc. For comparison, dashed lines, in the same color as the solid lines, represent each $\rho_{halo} \times S(l, b, D)$, where $S$ is the selection function as given in Equation (17). All color-coding and lines are comparable to those in Figure 4.

Starting with the results of fitting the power-law model to the ($\Delta l = 30^{\circ}) \times (\Delta b = 60^{\circ}$) 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^{\circ} < l < 30^{\circ}$, $-90^{\circ} < b < -30^{\circ}$, 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^{\circ}$ 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^{\circ} < l < 210^{\circ}$, $-30^{\circ} < b < 90^{\circ}$, 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.
mock halo power law profile with varying flattening $q(R_{gc})$
on ($\Delta l=30^{\circ}$) $\times$ ($\Delta b=60^{\circ}$)

<figure><img src="image_21.png" /><figcaption>**Figure 12.** One- and two-dimensional projections of the posterior probability distributions (pdf) of parameters ($r_0, q_0, q_\infty, n$) of the power law with varying flattening $q(R_{gc})$ (Equation (7)) fitted to a 30° $\times$ 60° patch of the cleaned sample. The blue lines and squares mark the maximum likely value of each parameter. The same patch on the sky as in the bottom panel of Figure 11 was chosen. 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 is increased in the case shown here, compared to the full cleaned sample. This is also reflected in the $1\sigma$ intervals given along with the best-fit parameters in the top right part of the figure. The best-fit parameters are given along with their $1\sigma$ intervals in the top right part of the figure. The covariance of the parameters is comparable to those found when fitting the complete cleaned sample.</figcaption></figure>

In the case of a power law with $q(R_{gc})$, as shown in Figure 15, the best-fit values for $q_0, q_\infty$ 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_0, q_\infty$, and the power-law index $n$ shows more scatter than for a power law with constant halo flattening. This might be caused
<figure><img src="image_22.png" /><figcaption>**Figure 13.** Angular distribution for the best-fit power-law parameters *q* and *n*, respectively, given per ($\Delta l = 30^{\circ}$) $\times$ ($\Delta b = 60^{\circ}$) bin. The number density is color-coded, as well as given in numbers. Empty cells are outside of the survey area. Some cells show a large deviation of *q* or *n* from the mean or from the expected value. These cells are highlighted with a thick frame. Reasons for those deviations are discussed in Section 4.3. The values for the cells in these plots are given in Table 4.</figcaption></figure>

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_0$, $q_\infty$ 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^{\circ} < l < 30^{\circ}$, $-90^{\circ} < b < -30^{\circ}$, and $300^{\circ} < l < 330^{\circ}$, $-30^{\circ} < b < 30^{\circ}$, due to the small number of sources within those regions.

Within $180^{\circ} < l < 210^{\circ}$, 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^{\circ} < l < 330^{\circ}$, $-30^{\circ} < b < 90^{\circ}$. 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 $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ bins for all four models in Table 8.

### 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.

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_0 = 25.0$ kpc, $q_0 = 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
<figure><img src="image_23.png" /><figcaption>**Figure 14.** Angular distribution for the best-fit Einasto profile parameters *q* and *n*, given per ($\Delta l = 30^{\circ}$) $\times$ ($\Delta b = 60^{\circ}$) bin. The number density is color-coded, as well as given in numbers. Empty cells are outside of the survey area. The Einasto profile parameter $r_{\text{eff}}$ was neglected here for the sake of clarity, and as we compare mainly the results on oblateness and steepness of the halo profile. Some cells show a large deviation of *q* or *n* from the mean or from the expected value. These cells are highlighted with a thick frame. Reasons for those deviations are discussed in Section 4.3. The values for the cells in these plots are given in Table 5.</figcaption></figure>

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(\rho_{\text{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 $\ln\langle\rho_{\text{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 $\Delta\ln\langle\rho\rangle = \ln(\rho_{\text{obs}}) - \ln\langle\rho_{\text{model}}\rangle$, as given in the last row of this figure. A $\Delta\ln\langle\rho\rangle < 0$ indicates that the best-fit model overestimates the number densities, whereas a $\Delta\ln\langle\rho\rangle > 0$ means that it underestimates the number density.

We find that the best-fit model leads, as expected, to a $\Delta\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).
<figure><img src="image_24.png" /><figcaption>**Figure 15.** Angular distribution for the best-fit parameters $q_0$, $q_\infty$, and $n$, respectively, of a power-law model with $q(R_{gc})$. The distribution is given on a ($\Delta l = 30^\circ$) $\times$ ($\Delta b = 60^\circ$) grid. The number density is color-coded, as well as given in numbers. Empty cells are outside of the survey area. The power-law parameter $r_0$ was neglected here for the sake of clarity, and as we compare mainly the results on oblateness and steepness of the halo profile. Some cells show a large deviation of $q_0$, $q_\infty$, or $n$ from the mean or from the expected value for the parameter in case. These cells are highlighted with a thick frame. Reasons for those deviations are discussed in Section 4.3. The values for the cells in this plot are given in Table 6.</figcaption></figure>

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 $\rho_{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
<figure><img src="image_25.png" /><figcaption>**Figure 16.** Angular distribution for the best-fit parameters $q_0, q_\infty$, and $n$ of an Einasto profile with $q(R_{gc})$. The distribution is given on a ($\Delta l = 30^\circ$) $\times$ ($\Delta b = 60^\circ$) grid. The number density is color-coded, as well as given in numbers. Empty cells are outside of the survey area. The power-law parameter $r_0$ was neglected here for the sake of clarity, and as we compare mainly the results on oblateness and steepness of the halo profile. Some cells show a large deviation of $q_0, q_\infty$, or $n$ from the mean or from the expected value for the value in case. These cells are highlighted with a thick frame. Reasons for those deviations are discussed in Section 4.3. The values for the cells in this plot are given in Table 7.</figcaption></figure>

estimator. This yields $\langle \rho_{\text{model}} \rangle_i$ for $i = 1...N$, and further
$\Delta(\ln\langle\rho\rangle)_i \equiv \ln(\rho_{\text{obs},i}) - \ln\langle\rho_{\text{model},i}\rangle$ for $i = 1...N$.

2. From the above, we can construct the variance
   $\sigma(\Delta(\ln\langle\rho\rangle)) \equiv \text{Var}\langle\Delta(\ln\langle\rho\rangle)_i\rangle$.

3. The 3D significance is then $\Delta(\ln\langle\rho\rangle)/\sigma(\Delta(\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 $\Delta(\ln\langle\rho\rangle) = 0$.

We find a significance of ~20 to >50 at regions that
coincide with the lower row of panels in Figure 17, and the
<table><caption>**Table 4**
Best-fit Parameters for the Power-law Model on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ Bins</caption><thead><tr><th>l</th><th>b</th><th>Sources</th><th>q</th><th>n</th></tr></thead><tbody><tr><td>0</td><td>-90</td><td>158</td><td>(0.307<sup>+0.126</sup><sub>-0.077</sub>)</td><td>3.53<sup>+0.18</sup><sub>-0.19</sub></td></tr><tr><td>0</td><td>-30</td><td>532</td><td>0.912<sup>+0.057</sup><sub>-0.068</sub></td><td>4.23<sup>+0.12</sup><sub>-0.11</sub></td></tr><tr><td>0</td><td>30</td><td>327</td><td>0.973<sup>+0.020</sup><sub>-0.039</sub></td><td>3.914<sup>+0.131</sup><sub>-0.127</sub></td></tr><tr><td>30</td><td>-90</td><td>582</td><td>0.557<sup>+0.052</sup><sub>-0.055</sub></td><td>3.79<sup>+0.10</sup><sub>-0.09</sub></td></tr><tr><td>30</td><td>-30</td><td>1099</td><td>0.896<sup>+0.050</sup><sub>-0.051</sub></td><td>4.89<sup>+0.06</sup><sub>-0.08</sub></td></tr><tr><td>30</td><td>30</td><td>289</td><td>0.943<sup>+0.04</sup><sub>-0.053</sub></td><td>4.25<sup>+0.16</sup><sub>-0.14</sub></td></tr><tr><td>60</td><td>-90</td><td>529</td><td>0.714<sup>+0.061</sup><sub>-0.063</sub></td><td>(3.39<sup>+0.08</sup><sub>-0.09</sub>)</td></tr><tr><td>60</td><td>-30</td><td>899</td><td>0.967<sup>+0.024</sup><sub>-0.041</sub></td><td>4.71<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>60</td><td>30</td><td>325</td><td>0.919<sup>+0.048</sup><sub>-0.054</sub></td><td>4.60<sup>+0.13</sup><sub>-0.14</sub></td></tr><tr><td>90</td><td>-90</td><td>260</td><td>0.813<sup>+0.076</sup><sub>-0.076</sub></td><td>3.86<sup>+0.14</sup><sub>-0.14</sub></td></tr><tr><td>90</td><td>-30</td><td>428</td><td>0.949<sup>+0.037</sup><sub>-0.058</sub></td><td>4.72<sup>+0.12</sup><sub>-0.12</sub></td></tr><tr><td>90</td><td>30</td><td>247</td><td>0.927<sup>+0.0476</sup><sub>-0.056</sub></td><td>4.862<sup>+0.094</sup><sub>-0.129</sub></td></tr><tr><td>120</td><td>-90</td><td>172</td><td>0.568<sup>+0.081</sup><sub>-0.085</sub></td><td>4.02<sup>+0.19</sup><sub>-0.17</sub></td></tr><tr><td>120</td><td>-30</td><td>312</td><td>0.935<sup>+0.046</sup><sub>-0.084</sub></td><td>4.92<sup>+0.06</sup><sub>-0.09</sub></td></tr><tr><td>120</td><td>30</td><td>232</td><td>0.960<sup>+0.027</sup><sub>-0.042</sub></td><td>4.80<sup>+0.12</sup><sub>-0.17</sub></td></tr><tr><td>150</td><td>-90</td><td>74</td><td>0.447<sup>+0.168</sup><sub>-0.159</sub></td><td>4.04<sup>+0.32</sup><sub>-0.28</sub></td></tr><tr><td>150</td><td>-30</td><td>215</td><td>0.952<sup>+0.036</sup><sub>-0.065</sub></td><td>4.71<sup>+0.14</sup><sub>-0.17</sub></td></tr><tr><td>150</td><td>30</td><td>277</td><td>0.919<sup>+0.048</sup><sub>-0.054</sub></td><td>4.82<sup>+0.11</sup><sub>-0.13</sub></td></tr><tr><td>180</td><td>-90</td><td>161</td><td>0.940<sup>+0.041</sup><sub>-0.061</sub></td><td>4.88<sup>+0.09</sup><sub>-0.13</sub></td></tr><tr><td>180</td><td>-30</td><td>318</td><td>0.967<sup>+0.026</sup><sub>-0.051</sub></td><td>(4.30<sup>+0.12</sup><sub>-0.12</sub>)</td></tr><tr><td>180</td><td>30</td><td>377</td><td>(0.730<sup>+0.068</sup><sub>-0.067</sub>)</td><td>(2.83<sup>+0.10</sup><sub>-0.11</sub>)</td></tr><tr><td>210</td><td>-90</td><td>98</td><td>0.640<sup>+0.00</sup><sub>-0.099</sub></td><td>4.61<sup>+0.22</sup><sub>-0.25</sub></td></tr><tr><td>210</td><td>-30</td><td>387</td><td>0.942<sup>+0.042</sup><sub>-0.077</sub></td><td>4.95<sup>+0.04</sup><sub>-0.06</sub></td></tr><tr><td>210</td><td>30</td><td>402</td><td>0.934<sup>+0.042</sup><sub>-0.052</sub></td><td>3.84<sup>+0.11</sup><sub>-0.11</sub></td></tr><tr><td>240</td><td>-30</td><td>292</td><td>0.959<sup>+0.040</sup><sub>-0.054</sub></td><td>4.90<sup>+0.07</sup><sub>-0.11</sub></td></tr><tr><td>240</td><td>30</td><td>476</td><td>0.933<sup>+0.039</sup><sub>-0.043</sub></td><td>4.79<sup>+0.11</sup><sub>-0.11</sub></td></tr><tr><td>270</td><td>-30</td><td>20</td><td>0.818<sup>+0.128</sup><sub>-0.193</sub></td><td>4.74<sup>+0.20</sup><sub>-0.38</sub></td></tr><tr><td>270</td><td>30</td><td>521</td><td>0.967<sup>+0.023</sup><sub>-0.033</sub></td><td>4.70<sup>+0.11</sup><sub>-0.11</sub></td></tr><tr><td>300</td><td>-30</td><td>2</td><td>(0.458<sup>+0.318</sup><sub>-0.191</sub>)</td><td>(2.20<sup>+1.66</sup><sub>-0.93</sub>)</td></tr><tr><td>300</td><td>30</td><td>520</td><td>0.992<sup>+0.006</sup><sub>-0.012</sub></td><td>3.44<sup>+0.10</sup><sub>-0.10</sub></td></tr><tr><td>330</td><td>-30</td><td>102</td><td>0.953<sup>+0.036</sup><sub>-0.066</sub></td><td>4.18<sup>+0.29</sup><sub>-0.28</sub></td></tr><tr><td>330</td><td>30</td><td>390</td><td>0.976<sup>+0.018</sup><sub>-0.031</sub></td><td>3.55<sup>+0.11</sup><sub>-0.11</sub></td></tr></tbody></table>

Note. Best-fit values for the power-law model, when carried out on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ bins. The table gives the bin limits from $(l, b)$ to $(l + \Delta l, b + \Delta 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.

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 $(R.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.

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

Note. Best-fit values for the Einasto profile, when carried out on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ bins. The table gives the bin limits from $(l, b)$ to $(l + \Delta l, b + \Delta 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.

and highlighting the region covered by their analysis. Upper-case 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).

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

Note. Best-fit values for the power-law model with $q(R_{gc})$, when carried out on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ bins. The table gives the bin limits from $(l, b)$ to $(l + \Delta l, b + \Delta 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.

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 \gtrsim 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.

## 5. Discussion

Before discussing possible implications of the results, it is worth discussing some potential sources of bias.

Carrying out the described modeling of the halo profile for the complete cleaned sample, we find that among models with constant flattening the distribution of the sources is reasonably well fit by a power-law model with $n = 4.40^{+0.05}_{-0.04}$ and a halo flattening of $q = 0.918^{+0.016}_{-0.014}$. The distance distribution is fit comparably well by a model with an Einasto profile with $n = 9.53^{+0.27}_{-0.28}$, an effective radius $r_{eff} = 1.07 \pm 0.10$ kpc, and a halo flattening of $q = 0.923 \pm 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}$ kpc, $n = 4.61 \pm 0.03$, $q_0 = 0.773^{+0.017}_{-0.016}$, and $q_{\infty} = 0.998^{+0.002}_{-0.001}$. The best-fit parameters for an Einasto profile with $q(R_{gc})$ are $r_0 = 26.7^{+2.2}_{-2.0}$ kpc, $q_0 = 0.779 \pm 0.018$, $q_{\infty} = 0.998^{+0.001}_{-0.002}$,
<table><caption>**Table 7**
Best-fit Parameters for the Einasto Profile with $q(R_{gc})$ on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ Bins</caption><thead><tr><th>l</th><th>b</th><th>Sources</th><th>r<sub>0</sub></th><th>q<sub>0</sub></th><th>q<sub>∞</sub></th><th>r<sub>eff</sub></th><th>n</th></tr></thead><tbody><tr><td>0</td><td>-90</td><td>158</td><td>23.4<sup>+1.0</sup><sub>-4.8</sub></td><td>(0.200<sup>+0.06</sup><sub>-0.01</sub>)</td><td>0.987<sup>+0.031</sup><sub>-0.230</sub></td><td>1.18<sup>+0.10</sup><sub>-0.10</sub></td><td>5.15<sup>+2.2</sup><sub>-0.38</sub></td></tr><tr><td>0</td><td>-30</td><td>532</td><td>13.2<sup>+6.5</sup><sub>-0.04</sub></td><td>0.621<sup>+0.303</sup><sub>-0.032</sub></td><td>0.995<sup>+0.17</sup><sub>-0.163</sub></td><td>1.21<sup>+0.11</sup><sub>-0.08</sub></td><td>8.80<sup>+0.73</sup><sub>-0.49</sub></td></tr><tr><td>0</td><td>30</td><td>327</td><td>25.3<sup>+4.4</sup><sub>-10.6</sub></td><td>0.840<sup>+0.070</sup><sub>-0.107</sub></td><td>0.999<sup>+0.004</sup><sub>-0.045</sub></td><td>1.25<sup>+0.08</sup><sub>-0.11</sub></td><td>9.55<sup>+0.82</sup><sub>-0.61</sub></td></tr><tr><td>30</td><td>-90</td><td>582</td><td>35.6<sup>+2.4</sup><sub>-3.9</sub></td><td>0.431<sup>+0.034</sup><sub>-0.033</sub></td><td>1<sup>+0.005</sup><sub>-0.049</sub></td><td>1.19<sup>+0.08</sup><sub>-0.12</sub></td><td>6.53<sup>+0.82</sup><sub>-0.22</sub></td></tr><tr><td>30</td><td>-30</td><td>1099</td><td>15.9<sup>+6.7</sup><sub>-1.6</sub></td><td>0.997<sup>+0.031</sup><sub>-0.189</sub></td><td>(0.762<sup>+0.156</sup><sub>-0.048</sub>)</td><td>1.19<sup>+0.09</sup><sub>-0.11</sub></td><td>7.42<sup>+0.15</sup><sub>-0.60</sub></td></tr><tr><td>30</td><td>30</td><td>289</td><td>38.3<sup>+1.7</sup><sub>-19.7</sub></td><td>0.800<sup>+0.12</sup><sub>-0.04</sub></td><td>0.998<sup>+0.006</sup><sub>-0.075</sub></td><td>1.21<sup>+0.10</sup><sub>-0.10</sub></td><td>8.56<sup>+0.97</sup><sub>-0.41</sub></td></tr><tr><td>60</td><td>-90</td><td>529</td><td>31.7<sup>+2.6</sup><sub>-4.3</sub></td><td>0.418<sup>+0.026</sup><sub>-0.049</sub></td><td>0.999<sup>+0.002</sup><sub>-0.034</sub></td><td>1.198<sup>+0.102</sup><sub>-0.095</sub></td><td>8.50<sup>+0.80</sup><sub>-0.46</sub></td></tr><tr><td>60</td><td>-30</td><td>899</td><td>27.3<sup>+6.0</sup><sub>-13.0</sub></td><td>0.999<sup>+0.011</sup><sub>-0.113</sub></td><td>0.989<sup>+0.006</sup><sub>-0.131</sub></td><td>1.21<sup>+0.08</sup><sub>-0.12</sub></td><td>7.85<sup>+0.42</sup><sub>-0.32</sub></td></tr><tr><td>60</td><td>30</td><td>325</td><td>37.3<sup>-1.5</sup><sub>-20.3</sub></td><td>0.803<sup>+0.099</sup><sub>-0.060</sub></td><td>0.999<sup>+0.008</sup><sub>-0.088</sub></td><td>1.18<sup>+0.105</sup><sub>-0.091</sub></td><td>7.60<sup>+0.766</sup><sub>-0.403</sub></td></tr><tr><td>90</td><td>-90</td><td>260</td><td>36.1<sup>+0.246</sup><sub>-0.589</sub></td><td>0.538<sup>+0.054</sup><sub>-0.066</sub></td><td>0.999<sup>+0.004</sup><sub>-0.047</sub></td><td>1.20<sup>+0.09</sup><sub>-0.10</sub></td><td>8.37<sup>+0.77</sup><sub>-0.55</sub></td></tr><tr><td>90</td><td>-30</td><td>428</td><td>37.7<sup>+4.9</sup><sub>-23.7</sub></td><td>0.998<sup>+0.023</sup><sub>-0.229</sub></td><td>0.999<sup>+0.021</sup><sub>-0.189</sub></td><td>1.23<sup>+0.06</sup><sub>-0.14</sub></td><td>7.64<sup>+0.64</sup><sub>-0.29</sub></td></tr><tr><td>90</td><td>30</td><td>247</td><td>36.6<sup>+4.6</sup><sub>-22.5</sub></td><td>0.860<sup>+0.077</sup><sub>-0.122</sub></td><td>0.999<sup>+0.014</sup><sub>-0.112</sub></td><td>1.193<sup>+0.086</sup><sub>-0.113</sub></td><td>6.794<sup>+0.790</sup><sub>-0.354</sub></td></tr><tr><td>120</td><td>-90</td><td>172</td><td>35.5<sup>+0.3</sup><sub>-10.4</sub></td><td>0.399<sup>+0.064</sup><sub>-0.069</sub></td><td>0.999<sup>+0.016</sup><sub>-0.041</sub></td><td>1.20<sup>+0.08</sup><sub>-0.12</sub></td><td>6.56<sup>+0.20</sup><sub>-0.16</sub></td></tr><tr><td>120</td><td>-30</td><td>312</td><td>14.2<sup>+20.5</sup><sub>-0.2</sub></td><td>0.998<sup>+0.026</sup><sub>-0.316</sub></td><td>0.999<sup>+0.029</sup><sub>-0.305</sub></td><td>1.19<sup>+0.10</sup><sub>-0.10</sub></td><td>6.79<sup>+0.65</sup><sub>-0.28</sub></td></tr><tr><td>120</td><td>30</td><td>232</td><td>37.9<sup>+6.1</sup><sub>-23.8</sub></td><td>0.957<sup>+0.015</sup><sub>-0.147</sub></td><td>0.999<sup>+0.009</sup><sub>-0.087</sub></td><td>1.18<sup>+0.11</sup><sub>-0.09</sub></td><td>7.34<sup>+0.69</sup><sub>-0.42</sub></td></tr><tr><td>150</td><td>-90</td><td>74</td><td>39.7<sup>+1.1</sup><sub>-11.5</sub></td><td>0.429<sup>+0.080</sup><sub>-0.088</sub></td><td>0.993<sup>+0.026</sup><sub>-0.033</sub></td><td>1.18<sup>+0.10</sup><sub>-0.10</sub></td><td>6.85<sup>+0.33</sup><sub>-0.38</sub></td></tr><tr><td>150</td><td>-30</td><td>215</td><td>38.4<sup>+12.9</sup><sub>-25.4</sub></td><td>0.999<sup>+0.077</sup><sub>-0.652</sub></td><td>0.996<sup>+0.08</sup><sub>-0.109</sub></td><td>1.19<sup>+0.10</sup><sub>-0.10</sub></td><td>7.67<sup>+0.51</sup><sub>-0.56</sub></td></tr><tr><td>150</td><td>30</td><td>277</td><td>14.0<sup>+8.7</sup><sub>-0.2</sub></td><td>0.997<sup>+0.032</sup><sub>-0.199</sub></td><td>0.916<sup>+0.051</sup><sub>-0.086</sub></td><td>1.18<sup>+0.10</sup><sub>-0.10</sub></td><td>7.62<sup>+0.57</sup><sub>-0.61</sub></td></tr><tr><td>180</td><td>-90</td><td>161</td><td>32.5<sup>+3.7</sup><sub>-16.0</sub></td><td>0.997<sup>+0.010</sup><sub>-0.118</sub></td><td>0.749<sup>+0.183</sup><sub>-0.080</sub></td><td>1.22<sup>+0.08</sup><sub>-0.13</sub></td><td>7.45<sup>+0.75</sup><sub>-0.83</sub></td></tr><tr><td>180</td><td>-30</td><td>318</td><td>12.2<sup>+8.5</sup><sub>-0.9</sub></td><td>0.994<sup>+0.031</sup><sub>-0.434</sub></td><td>0.997<sup>+0.004</sup><sub>-0.087</sub></td><td>1.22<sup>+0.107</sup><sub>-0.11</sub></td><td>(9.30<sup>+0.65</sup><sub>-0.56</sub>)</td></tr><tr><td>180</td><td>30</td><td>377</td><td>34.4<sup>+0.8</sup><sub>-6.1</sub></td><td>(0.328<sup>+0.034</sup><sub>-0.058</sub>)</td><td>0.998<sup>+0.005</sup><sub>-0.060</sub></td><td>1.30<sup>+0.09</sup><sub>-0.11</sub></td><td>(11.4<sup>+0.8</sup><sub>-0.6</sub>)</td></tr><tr><td>210</td><td>-90</td><td>98</td><td>38.5<sup>+6.2</sup><sub>-24.9</sub></td><td>0.6<sup>+0.1</sup><sub>-0.2</sub></td><td>0.748<sup>+0.070</sup><sub>-0.164</sub></td><td>1.22<sup>+0.08</sup><sub>-0.12</sub></td><td>7.48<sup>+0.10</sup><sub>-0.67</sub></td></tr><tr><td>210</td><td>-30</td><td>387</td><td>32.2<sup>+3.4</sup><sub>-18.4</sub></td><td>0.997<sup>+0.032</sup><sub>-0.297</sub></td><td>0.999<sup>+0.016</sup><sub>-0.152</sub></td><td>1.16<sup>+0.11</sup><sub>-0.10</sub></td><td>6.62<sup>+0.41</sup><sub>-0.32</sub></td></tr><tr><td>210</td><td>30</td><td>402</td><td>31.8<sup>+3.1</sup><sub>-9.2</sub></td><td>0.688<sup>+0.067</sup><sub>-0.079</sub></td><td>0.999<sup>+0.004</sup><sub>-0.040</sub></td><td>1.25<sup>+0.08</sup><sub>-0.11</sub></td><td>9.81<sup>+0.82</sup><sub>-0.46</sub></td></tr><tr><td>240</td><td>-30</td><td>292</td><td>36.7<sup>+2.7</sup><sub>-27.4</sub></td><td>0.997<sup>+0.014</sup><sub>-0.174</sub></td><td>0.996<sup>+0.014</sup><sub>-0.162</sub></td><td>1.19<sup>+0.08</sup><sub>-0.12</sub></td><td>6.98<sup>+0.58</sup><sub>-0.32</sub></td></tr><tr><td>240</td><td>30</td><td>476</td><td>29.0<sup>+3.57</sup><sub>-13.0</sub></td><td>0.786<sup>+0.095</sup><sub>-0.061</sub></td><td>0.999<sup>+0.006</sup><sub>-0.062</sub></td><td>1.18<sup>+0.10</sup><sub>-0.10</sub></td><td>6.92<sup>+0.68</sup><sub>-0.25</sub></td></tr><tr><td>270</td><td>-30</td><td>20</td><td>18.8<sup>+4.52</sup><sub>-4.10</sub></td><td>0.992<sup>+0.045</sup><sub>-0.390</sub></td><td>(0.544<sup>+0.250</sup><sub>-0.261</sub>)</td><td>1.20<sup>+0.08</sup><sub>-0.11</sub></td><td>7.10<sup>+0.98</sup><sub>-1.02</sub></td></tr><tr><td>270</td><td>30</td><td>521</td><td>38.0<sup>+3.49</sup><sub>-23.0</sub></td><td>0.939<sup>+0.043</sup><sub>-0.062</sub></td><td>0.999<sup>+0.009</sup><sub>-0.083</sub></td><td>1.17<sup>+0.12</sup><sub>-0.08</sub></td><td>7.83<sup>+0.60</sup><sub>-0.40</sub></td></tr><tr><td>300</td><td>-30</td><td>2</td><td>13.0<sup>+9.8</sup><sub>-1.3</sub></td><td>0.984<sup>+0.193</sup><sub>-0.722</sub></td><td>0.999<sup>+0.197</sup><sub>-0.733</sub></td><td>1.20<sup>+0.09</sup><sub>-0.10</sub></td><td>7.65<sup>+0.74</sup><sub>-1.27</sub></td></tr><tr><td>300</td><td>30</td><td>520</td><td>26.0<sup>+6.7</sup><sub>-11.6</sub></td><td>0.919<sup>+0.058</sup><sub>-0.087</sub></td><td>0.999<sup>+0.001</sup><sub>-0.019</sub></td><td>1.32<sup>+0.09</sup><sub>-0.14</sub></td><td>(11.9<sup>+0.9</sup><sub>-0.7</sub>)</td></tr><tr><td>330</td><td>-30</td><td>102</td><td>23.2<sup>+9.9</sup><sub>-9.2</sub></td><td>0.988<sup>+0.018</sup><sub>-0.248</sub></td><td>0.999<sup>+0.021</sup><sub>-0.190</sub></td><td>1.19<sup>+0.12</sup><sub>-0.07</sub></td><td>8.51<sup>+0.65</sup><sub>-0.86</sub></td></tr><tr><td>330</td><td>30</td><td>390</td><td>27.6<sup>+6.78</sup><sub>-7.71</sub></td><td>0.748<sup>+0.090</sup><sub>-0.066</sub></td><td>0.999<sup>+0.001</sup><sub>-0.029</sub></td><td>1.26<sup>+0.12</sup><sub>-0.10</sub></td><td>10.3<sup>+0.1</sup><sub>-0.5</sub></td></tr></tbody></table>

Note. Best-fit values for the Einasto profile with $q(R_{gc})$, when carried out on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ bins. The table gives the bin limits from $(l, b)$ to $(l + \Delta l, b + \Delta 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.

$r_{eff} = 1.04^{+0.25}_{-0.13}$ kpc, and $n = 8.78^{+0.33}_{-0.30}$. Allowing for a break in the power-law profile, we find a break radius of $r_{break} = 38.7^{+0.69}_{-0.58}$, a halo flattening of $q = 0.908^{+0.008}_{-0.006}$, and the inner and outer slopes $n_{inner} = 4.97^{+0.02}_{-0.05}$ and $n_{outer} = 3.93^{+0.05}_{-0.04}$, respectively.

From these fits, we find that two robust effects emerge: there is evidence for the stellar halo being distinctly more flattened at small radii ($q \sim 0.8$) and more spherical at large radii ($q \sim 1$). The flattening is consistent among all halo profiles we explored. There is no evidence for a steepening of the halo profile beyond 65 kpc as found by Deason et al. (2014), neither from our fits nor from our data.

We also fitted the halo profile for both the north and south Galactic hemisphere independently, in order to look for possible local variations in the best-fit model, especially in the halo flattening $q$ and steepness $n$. The north hemisphere contains 6880 sources, while the south hemisphere contains only 4145 sources because of the PS1 $3\pi$ survey footprint. We find that the steepness parameters $n$ of all best-fit hemisphere models compare well for both the north and south Galactic hemispheres and also compare well with the fit for the

complete halo. We further find that for models with constant flattening (thus the power-law and Einasto profile models) $q_{south} < q < q_{north}$, and the same applies for $q_0$ in the case of models with $q(R_{gc})$. For those models, we find that the value of the parameter $q_{\infty}$ is similar for the north and south hemispheres and the complete halo. However, we find that $r_{0,north} > r_{0,south} > r_0$ for both the power law with $q(R_{gc})$ and the Einasto profile with $q(R_{gc})$.

The results $q_{south} < q < q_{north}$ for models with constant flattening and $q_{0,south} < q_0 < q_{0,north}$, $q_{\infty,south} \sim q_{\infty,north} \sim q_{\infty,south}$, $r_{0,north} > r_{0,south} > r_0$ in the case of a radius-dependent flattening are consistent: by definition of $q(R_{gc})$ (see Equation (7)), $q_0$ is the flattening at the center, whereas $q_{\infty}$ is the flattening at large galactocentric radii, and $r_0$ is the exponential scale radius for the flattening. A larger $r_0$ means that the flattening of the inner halo, where we find $q_0 < q_{\infty}$, in effect extends out to a larger galactocentric radius than for a smaller $r_{eff}$, thus leading to a larger fraction of the halo being more flattened.

Thus, the generalized result suggests that the south Galactic hemisphere is somewhat more flattened than the north Galactic hemisphere.
<table><caption>**Table 8**
Mean and Variance for Best-fit Parameters on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ Bins</caption><thead><tr><th>Density Model</th><th>Parameter Mean and Variance</th></tr></thead><tbody><tr><td>Power-law model</td><td>mean(n) = 4.25, Var(n) = 0.443,<br>mean(q) = 0.840, Var(q) = 0.0333</td></tr><tr><td>Einasto profile</td><td>mean(r<sub>eff</sub>) = 1.22, Var(r<sub>eff</sub>) = 0.00252,<br>mean(q) = 0.873, Var(q) = 0.0191<br>mean(n) = 3.85, Var(n) = 8.76</td></tr><tr><td>Power-law model with q(R<sub>gc</sub>)</td><td>mean(r<sub>0</sub>) = 36.0, Var(r<sub>0</sub>) = 36.3,<br>mean(q<sub>0</sub>) = 0.856, Var(q<sub>0</sub>) = 0.0443,<br>mean(q<sub>&infin;</sub>) = 0.859, Var(q<sub>&infin;</sub>) = 0.0559,<br>mean(n) = 4.45, Var(n) = 0.441</td></tr><tr><td>Einasto profile with q(R<sub>gc</sub>)</td><td>mean(r<sub>0</sub>) = 29.1, Var(r<sub>0</sub>) = 81.5,<br>mean(q<sub>0</sub>) = 0.790, Var(q<sub>0</sub>) = 0.058,<br>mean(q<sub>&infin;</sub>) = 0.958, Var(q<sub>&infin;</sub>) = 0.011,<br>mean(r<sub>eff</sub>) = 1.21, Var(r<sub>eff</sub>) = 0.001,<br>mean(n) = 7.961, Var(n) = 2.01</td></tr></tbody></table>

**Note.** Mean and variance for the best-fit parameters of all four halo models on $\Delta l = 30^{\circ}$, $\Delta b = 60^{\circ}$ bins.

Our finding that a best-fit model requires a flattened, and especially a varying-flattened, halo with a smaller *q* (minor-to-major-axis ratio) in the inner parts and a larger value, *q* ~ 1, in the outskirts is supported by many other studies.

In our subsequent analysis, we found that the halo might be more irregular than only being influenced by a flattened halo or flattened inner and outer halo (“dichotomy” of the halo). 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.

Using our best-fit halo model, we then continued by computing the residuals in order to find local deviations from the smooth halo described by the best-fit model. We found striking overdensities and compare them to the ones discovered by Sesar et al. (2007, 2010). Additionally, we find new overdensities in regions of the sky not accessible to Sesar et al. (2007, 2010).

After describing the outcome of our study and its scientific relevance, we now briefly discuss possible sources of bias in the maximum likelihood analysis. Removing known overdensities such as streams, globular clusters, and dwarf galaxies, as well as the Milky Way disk and bulge, from our sample, thus producing the “cleaned sample,” was crucial. Our cuts on the disk and bulge are fairly broad. For the Sgr stream, we tried to find a good trade-off between removing most of the stream and not removing too many background stars, and thus we decided to remove sources based on their probability for being associated with the Sgr stream as shown in our previous work (Hernitschek et al. 2017, see Equation (11) therein).

Another crucial point is our assumption on the cleanness of the RRab sample, as well as on its distance precision. For both we refer to Sesar et al. (2017), who claim, based on extensive testing, a high-latitude sample purity of 90%, a sample completeness of 90% within 60 kpc and ≥80% at 80 kpc, and a distance precision of 3%. To account for the distance-dependent completeness, we introduced a term in our selection function.

Regions with high dust extinction can add severe uncertainties in the study of the distribution of stars in the Galaxy. We account for that with our cut on the region around the Galactic disk, $|b| \ge 10^{\circ}$.

The halo fit can also be influenced by up to now unknown overdensities. On the other hand, we can use the halo fit to identify such overdensities, as well as to infer the distortion of the halo from an ellipsoid that is flattened in the Z direction.

## 5.1. Comparison with Previous Results

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\pi$ 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
<figure><img src="image_26.png" /><figcaption>**Figure 17.** Density plots in the Cartesian reference frame (X, Y, Z) for the best-fit model, as well as its residuals. Densities are each color-coded according to the legend. Our results are described in detail in Section 4.3. First row: realization of a mock “cleaned sample” from the best-fit model with the selection function applied. This sample consists of 11,025 sources, the same number of sources as in the observed cleaned sample. Second row: number density of the observed cleaned sample at each (X, Y, Z), given in Figure 2, using a nearest-neighbor approach. Third row: mean model density from applying the same estimation of the number density to 10 realizations of mock samples from the best-fit model, where a single mock sample looks like those given in the first row. Last row: logarithmic residuals of the best-fit model. A $\Delta \ln\langle\rho\rangle < 0$ indicates that the best-fit model overestimates the number densities, whereas a $\Delta \ln\langle\rho\rangle > 0$ means that it underestimates the number density. The green color ($\Delta \ln\langle\rho\rangle = 0$) corresponds to the density predicted by our best-fit halo model. We find $\Delta \ln\langle\rho\rangle \sim 0$ over wide ranges, but also regions where the model underestimates the number density (yellow to red) and thus these regions are overdensities. The dark-blue regions are edge effects when the samples become sparse at the survey's outskirts. The overdensities are of further interest; we label them according to Sesar et al. (2007).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure 18.** In order to estimate the significance of the overdensities shown in Figure 17, we calculate their variance and significance as described in Section 4.3.2.</figcaption></figure>

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 \sim 1$ at 20 kpc. However, later work by Sluis & Arnold (1998) shows a constant flattening of $q \sim 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 $\sim 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 \sim 0.7$ out to 35 kpc.

Carollo et al. (2007, 2010) found that the inner halo is highly flattened with axis ratios of $q \sim 0.6$, whereas the outer halo is more spherical with axis ratios of $q \sim 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 $\sim 5000$ RR Lyrae stars from the recalibrated LINEAR data set, spanning $5 \text{ kpc} < D < 30 \text{ kpc}$ over $\sim 8000 \text{ deg}^2$ 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_{\text{inner}} = 1$, $n_{\text{outer}} = 2.7$, and $r_{\text{break}} = 16 \text{ 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 \sim 2.5$, to a much steeper one outside of about 20–30 kpc that is consistent with $n \sim 4$. The earliest evidence for that is from Saha (1985), who found that RR Lyrae are well
<figure><img src="image_28.png" /><figcaption>**Figure 19.** Plot of the overdensities in an (R.A., *D*) projection similar to that in Sesar et al. (2007). The green color ($\Delta \ln\langle\rho\rangle = 0$) corresponds to the model density, yellow and red regions are overdensities, and blue regions are underdensities, analogous to Figure 17. 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. (2005) 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). A detailed discussion is given in Section 4.3.3.</figcaption></figure>

described by a BPL with $n \sim 3$ out to 25 kpc and $n \sim 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 kpc $< R_{gc} < 80$ 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 \pm 0.02$ at 10 kpc to $0.8 \pm 0.03$ at large radii.

Bell et al. (2008) used $\sim 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 $\sim 20–30$ kpc and power-law slopes of $n \sim 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 \sim 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$ kpc $< R_{gc} < 65$ kpc, and $10$ for $65$ kpc $< R_{gc} < 100$ 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
<figure><img src="image_29.png" /><figcaption>**Figure 20.** Comparison of two of our best-fit models, the power law, and the Einasto profile, each with radius-dependent $q(R_{gc})$, to best-fit models of other works. Their model parameters are given in Table 9, whereas the parameters of our best-fit models are given in Table 2. Most best-fit models compare well to ours within their distance range. However, the models by Xue et al. (2015) and Slater et al. (2016) are slightly shallower, and the model by Deason et al. (2014) shows a BPL shape that we find neither from our fits nor from our data.</figcaption></figure>

<table><caption>**Table 9**
Model Parameters for Selected Other Surveys</caption><thead><tr><th>Paper</th><th>Density Model</th><th>Parameter Mean and Variance</th><th>Distance Range</th></tr></thead><tbody><tr><td>Watkins et al. (2009)</td><td>BPL</td><td>n<sub>inner</sub> = 2.4, n<sub>inner</sub> = 4.5,<br>r<sub>break</sub> = 23 kpc, q = 1 (assumes no flattening)</td><td>5 kpc &lt; R<sub>gc</sub> &lt; 115 kpc</td></tr><tr><td>Deason et al. (2014)</td><td>BPL with three segments</td><td>α<sub>1</sub> = 2.5, α<sub>2</sub> = 4.5, α<sub>out</sub> = 10,<br>r<sub>c</sub> = 25 kpc, r<sub>break</sub> = 65 kpc, q = 1 (assumes no flattening)</td><td>10 kpc &lt; R<sub>gc</sub> &lt; 100 kpc</td></tr><tr><td>Cohen et al. (2015)</td><td>power law</td><td>n = 3.8 ± 0.3, q = 1 (assumes no flattening)</td><td>50 kpc &lt; R<sub>gc</sub> &lt; 115 kpc</td></tr><tr><td>Xue et al. (2015)</td><td>power law with q(R<sub>gc</sub>)</td><td>α = 4.2 ± 0.1, q<sub>0</sub> = 0.2 ± 0.1, q<sub>∞</sub> = 0.8 ± 0.3, r<sub>0</sub> = 6 ± 1 kpc</td><td>10 kpc &lt; R<sub>gc</sub> &lt; 80 kpc</td></tr><tr><td>Slater et al. (2016)</td><td>power law</td><td>n = 3.5 ± 0.2, q = 1 (assumes no flattening)</td><td>20 kpc &lt; R<sub>gc</sub> &lt; 80 kpc</td></tr></tbody></table>

**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.

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 \pm 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.

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.

## 6. Summary

We used a sample of 44,403 PS1 RRab stars from Sesar et al. (2017) in order to determine the spatial structure of the Galactic halo using Pan-STARRS1 $3\pi$ RR Lyrae. We excluded known overdensities, among them the Sagittarius stream, dwarf galaxies such as the Draco dSph, and globular clusters. Also, we cut out sources too close to the Galactic plane ($|b| < 10^\circ$) or too close to the Galactic center ($R_{gc} \le 20$ kpc). We end up with a sample of 11,025 RRab in the Galactic stellar halo, called the “cleaned sample.” Each RRab star has a highly precise distance (3% uncertainty). The sample is very pure and with high completeness. Each star has a photometric light curve that resembles that of an RRab, guaranteeing a very low level of interlopers.

A forward-modeling approach using different density models for the Galactic halo profile, a selection function of the sample describing the aforementioned cuts to exclude overdensities, and the distance-depending completeness was then applied to this sample.

Our basic result is that the stellar Galactic halo, when described purely by RRab stars outside of known overdensities, can be characterized by a power law with an exponent of $n = 4.61 \pm 0.03$ and a varying flattening ($q(R_{gc})$) with a more oblate inner halo and an almost spherical outer halo, as described by the parameters $r_0 = 25.0_{-1.7}^{+1.8}$ kpc, $q_0 = 0.773_{-0.016}^{+0.017}$, and $q_\infty = 0.998_{-0.001}^{+0.002}$. From our halo fits, we find that three robust effects emerge: there is evidence for the stellar halo being distinctly more flattened at small radii ($q \sim 0.8$) and more spherical at large radii ($q \sim 1$), and the flattening is consistent among all halo profiles we explored. We have no indication that there is a truncation or break in the halo profile within the range we probe by the RRab sample. As discussed in Section 5.1, broadly speaking the results of our work are largely consistent with most earlier work. However, we do not find a BPL halo claimed by Deason et al. (2014), in particular, we cannot reproduce their extreme power law of $n = 6-10$ beyond 65 kpc, but we do confirm their power-law slope of 4.8 for the region $25 \text{ kpc} < R_{gc} < 45 \text{ kpc}$.

Further, we claim that the estimate of the power-law parameter given by De Propris et al. (2010), $n \sim 2.5$, is too small to agree with our results.

To explore further, we fitted the halo profile for both the north and south Galactic hemisphere independently, in order to look for possible local variations in the best-fit model, especially in the halo flattening $q$ and steepness $n$. Our generalized result suggests that the south Galactic hemisphere is somewhat more flattened than the north Galactic hemisphere.

The final step in our analysis of the structure of the outer halo of the Milky Way was to compute the residuals from our results as compared to the smooth halo described by the best-fit model. This difference map was used to find local deviations from the smooth halo described by the best-fit model. We found striking overdensities and compare them to the ones discovered earlier by Sesar et al. (2007, 2010). Additionally, we find new overdensities that are in regions of the sky not accessible to Sesar et al. (2007, 2010).

The overarching goals of studies of the Milky Way's outer halo are to determine the total mass of the Galaxy and, to the extent possible, the origin and importance of substructure from which one can infer clues regarding the accretion history of the formation of the Galaxy. Having established the spatial profile of the Galactic halo and evaluated at least partially the local deviations from smooth structure as present but not overwhelming over the regime from $R_{gc} = 20$ to 130 kpc, the next step toward these goals is a study of the kinematics of the outer halo. We need 6D information, i.e., positions, distances, proper motions, and radial velocities. We already have the first two of this list. However, unfortunately the accuracy of Gaia proper motions of RRab in the outer halo at the large $R_{gc}$ we have probed is poor in Gaia DR1, and even after the end of the Gaia mission, it is still not as accurate as might be desired. Determination of the radial velocity distribution as a function of $R_{GC}$ out to these large distances will require a massive dedicated spectroscopic program at a large telescope. Such an effort has been initiated.

H.-W.R. acknowledges funding from the European Research Council under the European Unions Seventh Framework Programme (FP 7) ERC grant agreement no. 321035.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant no. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under grant no. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), and the Los Alamos National Laboratory.

## ORCID iDs

Nina Hernitschek https://orcid.org/0000-0003-1681-0430
Judith G. Cohen https://orcid.org/0000-0002-8039-4673
Hans-Walter Rix https://orcid.org/0000-0003-4996-9069
Branimir Sesar https://orcid.org/0000-0002-0834-3978
Nicolas F. Martin https://orcid.org/0000-0002-1349-202X
Eugene Magnier https://orcid.org/0000-0002-7965-2815
Richard Wainscoat https://orcid.org/0000-0002-1341-0952
Nick Kaiser https://orcid.org/0000-0001-6511-4306
John L. Tonry https://orcid.org/0000-0003-2858-9657
Klaus Hodapp https://orcid.org/0000-0003-0786-2140
Ken Chambers https://orcid.org/0000-0001-6965-7789
Heather Flewelling https://orcid.org/0000-0002-1050-4056
William Burgett https://orcid.org/0000-0003-4401-9582

## References

Akhter, S., Da Costa, G. S., Keller, S. C., et al. 2012, ApJ, 756, 23
Bell, E. F., Zucker, D. B., Belokurov, V., et al. 2008, ApJ, 690, 295
Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529
Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148
Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931
Carollo, D., Beers, T. C., Chiba, M., et al. 2010, ApJ, 712, 692

Carollo, D., Beers, T. C., Lee, Y. S., et al. 2007, Nature, 450, 1020

Catelan, M. 2009, in IAU Symp. Ser. 258, The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse (Cambridge: Cambridge Univ. Press), 209

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Cohen, J. G., Sesar, B., Banholzer, S., et al. 2015, in IAU Symp. 317, The General Assembly of Galaxy Halos: Structure, Origin and Evolution, ed. A. Bragaglia et al. (Cambridge: Cambridge Univ. Press), 91

Cohen, J. G., Sesar, B., Bahnolzer, S., et al. 2017, ApJ, 849, 150

Das, P., & Binney, J. 2016, MNRAS, 460, 1725

Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903

Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012, MNRAS, 425, 2840

Deason, A. J., Belokurov, V., Koposov, S. E., et al. 2014, ApJ, 787, 30

Deason, A. J., Van der Marel, R. P., Guhathakurta, P., et al. 2013, ApJ, 766, 24

de Jong, J. T. A., Yanny, B., Rix, H.-W., et al. 2010, ApJ, 714, 663

De Propris, R., Harrison, C. D., & Mares, P. J. 2010, ApJ, 719, 1582

Diemand, J., Moore, B., & Stadel, J. 2004, MNRAS, 352, 535

Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013a, ApJ, 763, 3

Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013b, ApJ, 765, 154

Drake, A. J., Djorgovski, S. G., Catelan, M., et al. 2017, MNRAS, 469, 368

Drake, A. J., Graham, M. J., Djorgovski, S. G., et al. 2014, ApJS, 213, 9

Einasto, J. 1965, TrAlm, 5, 87

Einasto, J., & Haud, U. 1989, A&A, 23, 89

Font, A. S., Johnston, K. V., Bullock, J. S., et al. 2006, ApJ, 638, 585

Font, A. S., Johnston, K. V., Ferguson, A. M. N., et al. 2008, ApJ, 673, 215

Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2

Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, Vol. 5 (Berkeley, CA: Mathematical Sciences Publishers)

Graham, A. W., Merritt, D., Moore, B., et al. 2006, AJ, 132, 2701

Grillmair, C. J., & Carlin, J. L. 2016, in Tidal Streams in the Local Group and Beyond, Vol. 420, ed. H. J. Newberg & J. L. Carlin (Basel: Springer International), 87

Hawkins, M. R. S. 1984, MNRAS, 206, 443

Hernitschek, N., Rix, H.-W., Schlafly, E. F., et al. 2016, ApJ, 817, 73

Hernitschek, N., Rix, H.-W., Sesar, B., et al. 2017, ApJ, 850, 96

Iorio, G., Belokurov, V., Erkal, D., et al. 2017, MNRAS, 474, 2142

Ivezić, Ž, Goldston, J., Finlator, K., et al. 2000, AJ, 120, 963

Ivezić, Ž, Lupton, R. H., Anderson, S., et al. 2003, arXiv:astro-ph/0301400

Ivezić, Ž, Vivas, A. K., Lupton, R. H., & Zinn, R. 2005, AJ, 129, 1096

Jurić, M., Ivezić, Ž, Brooks, A., et al. 2008, ApJ, 673, 864

Kaiser, N., Burgett, W., Chambers, K., et al. 2010, Proc. SPIE, 7733, 77330E

Karim, M. T., & Mamajek, E. E. 2017, MNRAS, 465, 472

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Kinman, T. D., Wirtanen, C. A., & Janes, K. A. 1966, ApJS, 13, 379

McConnachie, A. W. 2012, AJ, 144, 4

Merritt, D., Graham, A. W., Moore, B., et al. 2006, AJ, 132, 2685

Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039

Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21

Preston, G. W., Shectman, S. A., & Beers, T. C. 1991, ApJ, 375, 121

Rix, H.-W., & Bovy, J. 2013, A&ARv, 21, 61

Saha, A. 1984, ApJ, 283, 580

Saha, A. 1985, ApJ, 289, 310

Schwarz, G. 1978, AnSta, 6, 461

Sérsic, J. L. 1963, BAAA, 6, 41

Sesar, B., Grillmair, J. C., Cohen, J. G., et al. 2013a, ApJ, 776, 26

Sesar, B., Hernitschek, N., Mitrović, S., et al. 2017, AJ, 153, 204

Sesar, B., Ivezić, Ž, Grammer, S. H., et al. 2010, ApJ, 798, 717

Sesar, B., Ivezić, Ž, Lupton, R. H., et al. 2007, ApJ, 134, 2236

Sesar, B., Ivezić, Ž, Scott, J. S., et al. 2013b, AJ, 146, 21

Sesar, B., Jurić, M., & Ivezić, Ž 2011, ApJ, 731, 4

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 116

Slater, C. T., Nidever, D. L., Munn, J. M., et al. 2016, ApJ, 832, A206

Sluis, A. P. N., & Arnold, R. A. 1998, MNRAS, 297, 732

Smith, H. A. 2004, RR Lyrae Stars (Cambridge: Cambridge Univ. Press)

Soszyński, I., Udalski, A., & Szymański, M. K. 2016, AcA, 66, 131

Stubbs, C. W., Doherty, P., Cramer, C., et al. 2010, ApJS, 191, 376

Tonry, C. W., Stubbs, K. R., Lykke, K. R., et al. 2012, ApJ, 750, 99

Torrealba, G., Catelan, M., Drake, A. J., et al. 2015, MNRAS, 446, 2251

Vivas, A. K., & Zinn, R. 2006, AJ, 132, 714

Watkins, L. L., Evans, N. W., Belokurov, V., et al. 2009, MNRAS, 398, 1757

Wetterer, C. J., McGraw, J. T., Hess, T., et al. 1996, AJ, 112, 742

Xue, X.-X., Rix, H.-W., Yanny, B., et al. 2011, ApJ, 738, 79

Xue, X.-X., Rix, H.-W., Zhao, G., et al. 2008, ApJ, 684, 1143

Xue, X.-X., Rix, H.-W., Zhibo, M., et al. 2015, ApJ, 809, 144