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

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

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 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 disc and bulge as well as halo globular clusters, the sample contains sources within .
We then apply forward modeling using Galactic halo profile models with a sample selection function. Specifically, we use ellipsoidal stellar density models with a constant and a radius-dependent halo flattening . Assuming constant flattening , the distribution of the sources is reasonably well fit by a single power law with and , and comparably well by an Einasto profile with , an effective radius and a halo flattening of . If we allow for a radius-dependent flattening , we find evidence for a distinct flattening of of the inner halo at . 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. (2009), Iorio et al. (2017). 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.

\affiliation

Division of Physics, Mathematics and Astronomy, Caltech, Pasadena, CA 91125

\affiliation

Division of Physics, Mathematics and Astronomy, Caltech, Pasadena, CA 91125

\affiliation

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

\affiliation

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

\affiliation

Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany \affiliationObservatoire astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 11 rue de l’Université, F-6700 Strasbourg, France

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

Institute for Astronomy, University of Hawaiâi at Manoa, Honolulu, HI 96822, USA

\affiliation

\correspondingauthor

Nina Hernitschek

1 Introduction

The Milky Way’s extended stellar halo contains only a small fraction () 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 sight-lines through the Galaxy, a sensible description of the overall stellar halo shape requires nearly complete coverage of the sky.

As stellar haloes 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 model. The spatial distribution, as well as kinematics and metallicities and thus ages of halo stars enable 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 last 50 years 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 which obey a well defined period - luminosity relation, albeit with a small dependence on metallicity. Thus the mean luminosity of a 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 et al. (2006), Jurić et al. (2008), Catelan et al. (2009), Watkins et al. (2009), de Jong et al. (2010), Sesar et al. (2010), Deason et al. (2011) Sesar 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 et al. (2016), Iorio et al. (2017), 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 multi-epoch 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 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 (GDR1, Gaia Collaboration et al., 2016b) and 2MASS (Skrutskie et al., 2006).

Drake et al. (2014) used the Catalina Real-Time Transient Surveys (CRTS) DR1 to select a sample of 47,000 periodic variables of which 16,797 are RR Lyrae, the bulk of them are at  kpc. In total, the Catalina Surveys RR Lyrae Data Release 11 (Drake et al., 2013a, b, 2014; 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 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 () 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. (2017b) exploited the PS1 survey to create a sample of RRab which 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 which 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. (2017b).

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 as well as 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 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 overdensitis. 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. (2017b) from the Pan-STARRS1 3 survey. In this section, we describe the pertinent properties of the PS1 3 survey and the RR Lyrae light curves obtained, and recapitulate briefly the process of selecting the likely RRab, as laid out in Sesar et al. (2017b). We also briefly characterize the obtained candidate sample.

The Pan-STARRS 1 (PS1) survey (Kaiser et al., 2010) collected multi-epoch, multi-color observations undertaking a number of surveys, among which the PS1 3 survey (Chambers et al., 2016) is currently the largest. It has observed the entire sky north of declination in five filter bands () with a 5 single epoch depth of about 22.0, 22.0, 21.9, 21.0 and 19.8 magnitudes in , and , respectively (Stubbs et al., 2010; Tonry et al., 2012).

Starting with a sample of more than PS1 3 sources, Hernitschek et al. (2016) and Sesar et al. (2017b) subsequently selected a sample of 44,403 likely RRab stars, of which 17,500 are at  kpc, by applying machine-learning techniques based on light-curve characteristics. RRab stars are the most common type of RR Lyrae, making up 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 (), and complete ( at 80 kpc) at high galactic latitudes. The distance estimates are precise to , based on newly derived period-luminosity relations for the optical/near-infrared PS1 bands (Sesar et al., 2017b). Overall, this results in the widest (3/4 of the sky) and deepest (reaching  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 beyond 100 kpc.

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

The left panels of Fig. 1 show the source density of the PS1 sample of RRab stars for different distance bins , , . The right panels of the same figure show the sample after a cleaning to remove overdensities was applied; the details of this cleaning are descried later.

Fig. b is based on the same data but shown in the Cartesian reference frame 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 declination , 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, and also away from known large overdensities like the Sagittarius stream. Details of the process of removing these overdensities are given in Sec. 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), 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 et al. (2013).

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 and the general form of the models as , where are the model parameters (see Sec. 3.1) and are the observables with Galactic coordinates and the heliocentric distance .

An approach for fitting the spatial-density profiles of the RRab sample must account for the fact that the observed star counts do not reflect the underlying stellar distribution, but are strongly shaped by selection effects both from the survey itself as well as from selection cuts we chose while preparing the sample. We denote the spatial selection function as (see Sec. 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 , and that the density is ellipsoidal, allowing for a halo flattening along the direction. Oblate density distributions have , spherical have and prolate have .

The density is a function of right-handed Cartesian coordinates , that we evaluate through the Galactic longitude, Galactic latitude and heliocentric distance , so its dimension is :

 X =R⊙−Dcoslcosb (1) Y =−Dsinlcosb Z =Dsinb.

This reference frame is centered at the Galactic center. The Galactic disc is in the plane, with the axis pointing to the Sun and the axis to the North Galactic Pole. denotes the distance of the Galactic center from the Sun, in this work assumed to be 8 kpc, and main results of our work should not change for other values of within the assumed observational uncertainties.

The vertical position of the Sun with respect to the Galactic disc is uncertain, but it is estimated to be smaller than 50 pc (Karim et al., 2017; Iorio et al., 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 -axis is flipped with respect to our definition.

With Equ. (1), the Galactocentric distance is then defined as , and the flattening-corrected radius defined as where gives the halo flattening along the direction as a minor-to-major-axis ratio. This describes an oblate stellar halo that is stratified on concentric ellipsoids, where , , 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 in all cases.

Power-Law Profile

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

 ρhalo(X,Y,Z)=ρ⊙RRL(R⊙/rq)n. (2)

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

The free parameters are . Here, is the number density of RR Lyrae at the position of the Sun, the distance of the Sun from the Galactic center, and is the flattening-corrected radius. As we are not interested in absolute numbers, we are not fitting for .

Others presume a broken power-law (BPL) (e.g. Xue et al., 2015), where an inner and outer power law index 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:

 Unknown environment '% (3)

Einasto Profile

The Einasto profile (Einasto, 1965, 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) as well as dark matter haloes (Merritt et al., 2006; Navarro et al., 2010). It is given by

 γ(r)≡−dlnρ(r)dlnr∝rα (4)

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

 α=−dnn(rqreff)1/n. (5)

This can be rearranged to

 ρhalo(rq)≡ρ0exp{−dn[(rq/reff)1/n−1]}, (6)

where is the (here irrelevant) normalization, is the effective radius, is the concentration index. The parameter is a function of , where for , a good approximation is given by (Graham et al., 2006).

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

The free parameters of an Einasto profile with a constant flattening are .

Profiles with Varying Flattening

The models described so far assume a constant flattening . 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 non-constant flattening of the halo, parameterized by the Galactocentric radius. To describe such a radial variations of the stellar halo’s flattening, we consider the functional form for as:

 q(Rgc)=q∞−(q∞−q0)exp⎛⎜ ⎜⎝1−√R2gc+r20r0⎞⎟ ⎟⎠ (7)

with being the flattening at the center, being the flattening at large Galactocentric radii, and being the exponential scale radius over which the change of flattening occurs.

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

All other equations to describe the radial profile given above apply from the previously described Einasto and power law profile, replacing only with , thus replacing the fitting parameter with three fitting parameters , and .

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 non-complete 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 catalogue. 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 the smooth component of the halo is also structured, but at a level that is below our resolution.

Our selection function is binary so that is always equal to 1 except for the points 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. (2017b) were selected uniformly from the set of objects in the PS1 3 survey in the area and apparent magnitude range available for this survey. The selection completeness and purity are uniform over a wide range of apparent magnitude up to a flux-averaged r-band magnitude of 20 mag (Sesar et al., 2017b), which is described later on in Equ. (12).

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

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

The selection function is thus composed of:

 Extra open brace or missing close brace (8)

where describes the selection cuts of the sample introduced by the survey and Sesar et al. (2017b) itself leading to the 44,403 RRab stars, and describes area cuts to exclude overdensities.

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

 SRRL(l,b,D)={1,if δ>−30\arcdegandDmin

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

 Sbulge,disc(l,b,D)={1,if% |b|≥10\arcdegandRgc≥20 kpc0,else. (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, (Hernitschek et al., 2017, see Equ. (11) therein). We excise sources with as members of the Sgr stream, leading to a selection function of

 Ssgr(l,b,D)={1,if psgr(l,b,D)<0.20,else. (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 into .

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. (2017b) find that the RRab selection function is approximately constant at for a flux-averaged -band magnitude mag, after which it steeply drops to zero at mag. Writing as , the selection function characterizing the distance-dependent completeness is

 Sc(rF)=L−L1+exp(−k(rF−x0)) (12)

with (Sesar et al., 2017b):

 L =0.91 (13) k =4.0 (14) x0 =20.6 (15) rF =2.05log(D)+11. (16)

In addition, we estimated the distance-dependent purity to supplement the overall sample purity that was given as 90% by Sesar et al. (2017b). Using the RRab sample within SDSS S82, as done by Sesar et al. (2017b) 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 band. In contrast, over the same magnitude range, the completeness drops from 91% to 80%. The faintest RRab in S82 (which we use as validation set, see (Sesar et al., 2017b)) is found at =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 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., 2017b). 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 kpc, we included it as part of the selection function, .

We end up with a selection function

 S(l,b,D)= SRRL(l,b,D)×Sc(D)×Sp(D)×Sarea(l,b,D) (17) = SRRL(l,b,D)×Sc(D) ×Sbulge,disc(l,b,D) Extra open brace or missing close brace

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 20142 and a list of currently known halo streams by Grillmair et al. (2016), we select overdensities that could show up in a survey that covers the position and distance cuts of PS1 3. We check each overdensity to see if it appears in the RRab sample, and if so, cut it out by defining a selection box in . We end up with the cuts described in Table 1.

After excising stars using , the sample reduces to 11,025 RRab stars which we call the “cleaned sample”. The original and the cleaned sample are shown in Fig. b.

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 and the selection function at hand, we can directly calculate the likelihood of the data given the model , the fitting parameters and the selection function following Bovy et al. (2012).

The normalized un-marginalized log likelihood for the -th star with the observables is then

 lnp(Di|θ)=ρRRL(Di|θ)|J|S(li,bi,Di)∫∫∫ρRRL(l,b,D|θ)|J|S(l,b,D)dldbdD (18)

where the normalization integral is over the observed volume. The Jacobian term reflects the transformation from to coordinates.

We evaluate the logarithmic posterior probability of the parameters of the halo model, given the full data and a prior , with

 lnp(D|θ)=∑ilnp(Di|θ) (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 parameters space with Goodman & Weare’s Affine Invariant Markov Chain Monte Carlo (Goodman et al., 2010), making use of the Python module emcee (Foreman 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 pdf can be well-described by a Gaussian distribution, the difference between the 15.87th and 84.13th percentile is equal to 1.

The calculation of the normalization integral in Equ. (18) is complicated by the presence of the selection function , 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 give robust results. We decided to calculate the integral instead on a fine regular grid that is wide.

Model Priors

We now lay out the “pertinent range”, across which the model priors are given. We set a different prior distribution for each of the five following cases:

power law model:

 lnp(θ)= Uniform(1.0

BPL model:

 lnp(θ)= Uniform(1.0

where , give the Galactocentric distance range available in the sample.

Einasto profile:

 lnp(θ)= Uniform(log(0.01)

power law model with :

 lnp(θ)= Uniform(log(0.01)

Einasto profile with :

 lnp(θ)= Uniform(log(0.01)

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 (Sec. 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 , or Einasto profile with . We then applied a 3% distance uncertainty, superimpose the sample with faint and far Gaussian blobs away from the regions excluded by the selection function to simulate unknown overdensities, added the RRab known as members of the Sgr stream, and then applied the selection function. After that, the sample has 12,000 sources, and we randomly sample 11,025 sources to match the cleaned observed sample.

An example of a simulated distribution of halo RR Lyrae is shown in the upper panel of Fig. 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 posterior probability distributions (pdf) for fitting one of these mock halos, along with the parameters used to generate the mock halo, are given in Fig. 3.

4 Results

We now present the results of applying the modeling from Sec. 3.1 to the cleaned sample of RRab stars as described in Sec. 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 (Navarro et al., 2004; Diemand et al., 2004; Merritt et al., 2006; Graham et al., 2006). Both models are fitted with a constant halo flattening parameter as well as with a distance-dependent flattening . We illustrate these results in three ways: first by showing the predicted distribution by the best-fit models, then by showing the joint posterior distribution functions of the halo model parameters of each model; and third, we compare the models using the Bayesian information criterion (BIC).

First, we discuss the result of fitting the complete cleaned 3 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 , 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 Sec. 3.1 and the selection function as described in Sec. 3.2, we apply our likelihood approach (Sec. 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 8,917 sources. Fig. a 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 uncertainties estimated as the 15.87th and 84.13th percentiles, and the maximum log likelihood . We also give the BIC, a measure for model comparison described in Sec. 4.2.

The one- and two-dimensional projections of the posterior probability distributions (pdf) for each model are given in Fig. 5 to 9.

For the power-law and BPL model, the pdf shows an almost Gaussian-like distribution with no covariance between the model parameters and .
For the Einasto profile, as for the power law, the concentration index and the flattening parameter show an almost Gaussian distribution with no covariance. The concentration index is covariant with the effective radius parameter, .
The pdf of the power law model with shows covariance, and the pdf is strongly distorted from a Gaussian distribution.
For the Einasto profile with , the pdf is more complex and skewed. The fitting parameters , , 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 and a halo flattening of . Allowing for a break in the power-law profile, we find a break radius of , a halo flattening of , and the inner and outer slopes and , respectively. The distance distribution is fit comparably well by a model with an Einasto profile with , an effective radius and a halo flattening of . If we allow for a radius-dependent flattening , we find the best-fit parameters for a power law model with as , , , . The best-fit parameters for an Einasto profile with are , , , , .

We find here for both models with variable flattening, indicating that the inner halo is more flattened than the outer halo. Assuming a constant flattening 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 uncertainties are summarized in Tab. 2.

Our results confirm that if a varying flattening is assumed, the halo profile has a 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) and Carollo et al. (2010), as well as Xue et al. (2015), Das et al. (2016) and Iorio et al. (2017). For a BPL, we cannot confirm 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 Sec. 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 Bayesian information criterion (Schwarz, 1978, BIC).

The BIC takes into account both the statistical goodness of fit, as well as 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

 Extra open brace or missing close brace (25)

where are the model parameters, is the number of objects in the sample, and is the maximum likelihood, where we defined the likelihood function in Equ. (19) as .

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 Tab. 2 along with the best-fit parameters.

According to the BIC, we find the best-fit model to be the power law with , followed by the Einasto profile with , the constant-flattening power law, the constant-flattening Einasto profile and finally the BPL. As the values of BIC in Tab. 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 reasonable comparable to a Gaussian. As we see from Fig. 5 to 9, the power-law model as well as the Einasto profile have posterior distributions that compare well to a Gaussian distribution, whereas for the cases with , 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 et al. (1995): A indicates a 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 now 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 and steepness . We also look for previously unknown overdensities that we might find due to the spatial extent and depth of the RRab sample.

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 effects on our models – of rather restrictive functional form – are. The north hemisphere contains 6,880 sources, whereas the south hemisphere contains only 4,145 sources because of the PS1 3 survey footprint.

The results of this fitting attempt are summarized in Tab. 3. What we find is that the steepness parameters of all best-fit hemisphere models compare well for both the north and south Galactic hemisphere and also compare well with the fit for the complete halo. When taking a look at the flattening-related parameters, , , , , we find that for models with constant flattening (both the power law and Einasto profile models), . In the case of models with , we find the value of parameter being smaller for the south than for the north hemisphere, , whereas the value of the parameter is similar for both hemispheres. Furthermore, we find that for both the power law with and the Einasto profile with .

The results of finding for models with constant flattening, and , in the case of a radius-dependent flattening, are consistent: by definition of (Eq. (7)), is the flattening at center, is the flattening at large Galactocentric radii, and is the exponential scale radius over which the change of flattening occurs. A larger means that the flattening of the inner halo, where we find , is in force out to a larger radius than for a smaller , 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 , to further understand possible local variations in the best-fit model, especially in the halo flattening and steepness .

The angular source number density for the cleaned RRab sample, given per bin, is shown in Fig. 10.

The resulting best-fit parameters for the power law model, power-law model with , Einasto profile and Einasto profile with are shown in the Figures 13 to 16, and are given in the Tables Tab. 4 to 7 along with their 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, , and , (see Fig. 11). To illustrate the fitting performance further, in Fig. 12 we give the posterior probability distribution in the case of fitting a power law with varying flattening (Equ. (7)) to a patch of mock data. The posterior distribution is comparable to those when fitting the same halo profile to the full cleaned sample (see Fig. 3 for comparison). However, the width of the posterior probability distribution increases compared to the full cleaned sample. This is also reflected in the intervals given along with the best-fit parameters in the top right part of the figure.

Starting with the results of fitting the power law model to the patches, Fig. 13, we find that the flattening parameter is homogeneous over almost the complete sky. There are a few exceptions, i.e. for , , the resulting flattening parameter is suspiciously small. However, within that region, there are only 22 sources which makes a reliable fit difficult.

Again for , , the resulting , and here also the power-law index , is small. Since there are only 2 sources within that region, we obviously have to exclude that fit. Outside of these regions, the resulting and are relatively homogeneous, with a trend to smaller near the edges of the survey (see white empty region at in the figures) and at very high latitudes.

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

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

In the case of an Einasto profile with , as shown in Fig. 16, the best-fit values for , 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 , , and , , due to the small number of sources within that regions.

Within , the best-fit value of 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 appears 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 , . A higher-dimensional model is more affected by this than a lower-dimensional one; compare the extreme cases of the 2-dimensional power law model and the 5-dimensional Einasto profile with .

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

We give the mean and variance for the best-fit parameters on , bins for all four models in Tab. 8.

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 Fig. 17, we give density plots in the Cartesian reference frame (see Equ. (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 Fig. 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 with , , , , with applied selection function.
This mock sample looks very comparable to the observed cleaned sample, Fig. b(b). The position of the Galactic plane and Sgr stream, as removed by the selection function, are indicated. The dashed circle represents the cut. Sources within the circle but further away than 20 kpc are seen due to projection effects; the distinctly higher density just after 20 kpc shows the stars that are no longer affected by this distance cut.
We calculated the number density of our observed cleaned sample (given in Fig. b) at each , using a nearest-neighbor based adaptive Bayesian density estimator (Sesar et al., 2013b; Ivezić et al., 2005) yielding . The result is shown in the second row of Fig. 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 .
The logarithmic residuals of the best-fit model were calculated by subtracting the model mean number density (third row) from the observed number density (second row), yielding , as given in the last row of this figure. A indicates the best-fit model overestimates the number densities, whereas a means that it underestimates the number density.
We find that the best-fit model leads, as expected, to a over wide ranges, but 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 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), Fig. 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 due to a clumpy halo, i.e. a halo with diffuse overdensities. As indicated when discussing the selection function in Sec. 3.2, our estimation of purity and completeness are rather uncertain beyond distances of 90 kpc.
The dark blue regions in this plot occur due to 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 Sec. 4.3.3 and thus label them in Fig. 17 accordingly to Sesar et al. (2007).

Also, as the relatively sparseness of our cleaned RRab sample (11,025 sources within almost 3/4 of the sky and an extent of ) 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:

• We bootstrap the observed RRab sample times. We estimate the density of each of these bootstrapped samples, using the density estimator by Ivezić et al. (2005), resulting in for . We fit each of these bootstrapped samples, sample each of them 10 times and get the mean model density using the density estimator. This yields for , and further for .

• From the above, we can construct the variance .

• The 3D significance is then .

The resulting variance and significance are shown in Fig. 18, each projected using the mean. Per definition, the significance is 0 where .

We find a significance of to at regions that coincide with the lower row of panels in Fig. 17, and the variance is small and does not exceed 0.04 to 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.

Overdensities

We compare the overdensities found by us with those discovered previously by Sesar et al. (2007) and Sesar et al. (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 cleaned the SDSS Stripe 82 sample of RR Lyrae (Sesar et al., 2010), using then 366 highly probable RRab stars.

In Fig. 19, we plot the overdensities in a projection similar to Sesar et al. (2007) (see their Fig. 13, see also Fig. 11 in Sesar et al. (2010)), using our full range in declination and highlighting the region covered by their analysis. Upper-case letters denote overdensities found in the SDSS sample of Sesar et al. (2007) and Sesar et al. (2010), numbers denote overdensities found in their analysis of the Ivezić et al. (2000) sample (not numbered in Sesar et al. (2007)), and lower-case 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) and Sesar et al. (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 don’t find the overdensities D, H, K and M. However, the overdensitiy D appears very small in Sesar et al. (2007), so we cannot say for sure if we are able to identify anything at this position. Our sample does not cover exactly the same extent in RA and Dec, so the overdensities D, H, K, M as given in Sesar et al. (2007) lie in regions we don’t cover. We have checked adjacent slices in declination to see if we might detect those overdensities, but 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 declination.

The left wedge of Fig. 19 compares to the left wedge of Fig. 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 continuing the overall distribution of overdensities found before. We label them by lower-case 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’ trailing arm (Sesar et al., 2007; Ivezić et al., 2003a).

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

5 Discussion

Before discussing possible implications of the results, it is worth to discuss 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 and a halo flattening of . The distance distribution is fit comparably well by a model with an Einasto profile with , an effective radius and a halo flattening of . If we allow for a radius-dependent flattening , we find the best-fit parameters for a power law model with as , , , . The best-fit parameters for an Einasto profile with are , , , , . Allowing for a break in the power-law profile, we find a break radius of , a halo flattening of , and the inner and outer slopes and , respectively.

From these fits, we find two robust effects to emerge: There is evidence for the stellar halo being distinctly more flattened at small radii (), and more spherical at large radii (). 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 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 and steepness . The north hemisphere contains 6,880 sources, while the south hemisphere contains only 4,145 sources because of the PS1 3 survey footprint. We find that the steepness parameters of all best-fit hemisphere models compare well for each the north and south Galactic hemisphere 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), , and the same applying for in the case of models with . For those models, we find the value of the parameter being similar for the north and south hemisphere and the complete halo. However, we find that for both the power law with and the Einasto profile with .

The results for models with constant flattening, and , in the case of a radius-dependent flattening, are consistent: by definition of (see Eq. (7)), is the flattening at the center, whereas is the flattening at large Galactocentric radii, and is the exponential scale radius for the flattening. A larger means that the flattening of the inner halo, where we find , in effect extends out to a larger Galactocentric radius than for a smaller , 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.

Our finding that a best-fit model requires a flattened, and especially a varying-flattened, halo with a smaller (minor-to-major-axis ratio) in the inner parts and a larger value, , 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) and Sesar et al. (2010). Additionally, we find new overdensities in regions of the sky not accessible to Sesar et al. (2007) and Sesar et al. (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 disc and bulge, from our sample, thus producing the “cleaned sample”, was crucial. Our cuts on the disc and bulge are fairly broad. For the Sgr stream, we tried to find a good tradeoff between removing most of the stream and not removing too many background stars, and thus 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 Equ. (11) therein).

Another crucial point is our assumption on the cleanness of the RRab sample, as well on its distance precision. For both we refer to Sesar et al. (2017b), who claim, based on extensive testing, a high-latitude sample purity of 90%, a sample completeness of within 60 kpc and 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 disc, .

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 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 further 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 multi-epoch photometric surveys that have been carried out over the past decade.

It is important to remember that our sample selection procedure which we apply to the Pan-STARRS1 3 database takes advantage of the experience gained by attempting to use the SDSS (Sesar et al., 2007, 2010) and then to analyze the more difficult Palomar Transient Facility (Cohen et al., 2017) with few observations taken with a random cadence. Our PS1 sample of RRab stars (Sesar et al., 2017b) 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 from 20 up to 130 kpc, all with highly precise photometry. This yields a sample of high purity and completeness exceeding 80% out to . Out of these sources, 11,025 are found outside of dense regions such as stellar streams, the Galactic disc and bulge or globular clusters and stellar streams.

Furthermore, Sesar et al. (2017b) 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 a 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 tracer. As a first attempt, Preston et al. (1991) argued that the flattening changes from at 1 kpc to at 20 kpc. However, later work by Sluis et al. (1998) shows a constant flattening of without any evidence for a radius-depending flattening. A recent work by De Propris et al. (2010) utilizing 666 BHB stars from the 2dF QSO Redshift Survey states that the halo is approximately spherical with a power-law index of 2.5 out to 100 kpc. Sesar et al. (2010) who studied main-sequence turn-off (MSTO) stars from the Canada-France-Hawaii Telescope Legacy Survey find that the flattening is approximately constant at out to 35 kpc.

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

Others also find evidence that at least the innermost part of the halo is quite flattened: Sesar et al. (2010) fit the Galactic halo profile based on 5000 RR Lyrae stars from the recalibrated LINEAR data set, spanning over 8000 deg of the sky. They find for their best-fit model an oblate ellipsoid with an axis ratio of , and a double power-law model with , , , .

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

As a reason for the smaller minor-to-major axis ratio 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; Watkins et al., 2009; Sesar et al., 2011; Deason et al., 2011) find that it shows a shift from a relatively shallow one, as described by , to a much steeper one outside of about 20 to 30 kpc that is consistent with . The earliest evidence for that is from Saha (1985) who found that RR Lyrae are well described by a broken power law (BPL) with out to 25 kpc, and 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 the density distribution to be 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 using 1,757 stars from the SEGUE K-giant Survey. The majority of their sources are found at , whereas in our sample, 1,093 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 , = 15 kpc, , by an equally flattened broken power law with = 2.1, = 3.8, = 18 kpc (this is something we had not applied), and when fitting by an Einasto profile with , they find the halo being considerably more flattened as changes from at 10 kpc to at large radii.

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

Other estimates of the power law index, or slope, of the halo give break radii or effective radii of kpc, power-law slopes of (e.g. Sesar et al., 2011; Deason et al., 2011, 2014; Xue et al., 2015). For example, Sesar et al. (2011) fit the Galactic halo within heliocentric distances of kpc, steeper at , , , or a best-fit Einasto profile with , , , 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 for a power-law model, or for an Einasto profile(keep in mind the different definitions of ) for our sample starting at 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 Equ. (5)). We claim that the estimate of the power-law parameter from De Propris et al. (2010), , is too shallow, as well as that the estimate of the power-law parameter beyond 65 kpc from Deason et al. (2014), , is too steep, as we don’t see such a drop from our fit nor 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 , 4.5 for , 10 for . Whereas in the distance range , 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 significant steeper profile at 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 , as we find a power-law slope of 4.97 for . The small change in slope near 39 kpc may be related to how, in detail, the Sgr stream is removed. Beyond 40 kpc we find a much less steep power-law slope than do Deason et al. (2014).

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

Xue et al. (2015) and Slater et al. (2016) find a halo profile that is shallower than our best-fit models; it is difficult to know if this difference results from methological 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 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 , 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 Fig. 20.

Fig. 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, and thus gives 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 supporting 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 haloes 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 haloes, 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 of 44,403 PS1 RRab stars from Sesar et al. (2017b) in order to determine the spatial structure of the Galactic halo using Pan-STARRS 1 3 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 (), or too close to the Galactic center (). 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 which resembles that of a RRab, guaranteeing a very low level of interlopers.

A forward modeling approach using different density models for the Galactic halo profile, as well as a selection function of the sample describing the aforementioned cuts to exclude overdensities, as well as 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 , and a varying flattening () with a more oblate inner halo and an almost spherical outer halo, as described by the parameters , , . From our halo fits, we find three robust effects to emerge: There is evidence for the stellar halo being distinctly more flattened at small radii (), and more spherical at large radii (). 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 Sec. 5.1, broadly speaking the results of our work are largely consistent with most earlier work. However, we do not find a broken power-law halo claimed by Deason et al. (2014), in particular we cannot reproduce their extreme power law of beyond 65 kpc, but do confirm their power-law slope of 4.8 for the region .

Further, we claim that the estimate of the power-law parameter given by De Propris et al. (2010), , 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 and steepness . 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) and Sesar et al. (2010). Additionally, we find new overdensities which are in regions of the sky not accessable to Sesar et al. (2007) and Sesar et al. (2010).

The overaching 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 to 130 kpc, the next step towards 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 we have probed is poor in Gaia DR1, and even after the end of the Gaia mission, is still not as accurate as might be desired. Determination of the radial velocity distribution as a function of 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 n. [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, The 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, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory.