First all-sky view of the Galactic halo

The first all-sky view of the Milky Way stellar halo with Gaia+2MASS RR Lyrae

G. Iorio, V. Belokurov, D. Erkal, S. E. Koposov, C. Nipoti, F. Fraternali
Dipartimento di Fisica e Astronomia, Università di Bologna, via Gobetti 93/2, I-40129, Bologna, Italy
INAF - Osservatorio Astronomico di Bologna, via Gobetti 93/3, I-40129, Bologna, Italy
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
Department of Physics, University of Surrey, Guildford, GU2 7XH, UK
McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213, USA
Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands
Accepted XXX. Received YYY; in original form ZZZ

We exploit the first Gaia data release to study the properties of the Galactic stellar halo as traced by RR Lyrae. We demonstrate that it is possible to select a pure sample of RR Lyrae using only photometric information available in the Gaia+2MASS catalogue. The final sample contains about 21600 RR Lyrae covering an unprecedented fraction () of the volume of the Galactic inner halo ( kpc). We study the morphology of the stellar halo by analysing the RR Lyrae distribution with parametric and non-parametric techniques. Taking advantage of the uniform all-sky coverage, we test halo models more sophisticated than usually considered in the literature, such as those with varying flattening, tilt and/or offset of the halo with respect to the Galactic disc. A consistent picture emerges: the inner halo is well reproduced by a smooth distribution of stars settled on triaxial density ellipsoids. The shortest axis is perpendicular to the Milky Way’s disc, while the longest axis forms an angle of with the axis connecting the Sun and the Galactic centre. The elongation along the major axis is mild (), and the vertical flattening is shown to evolve from a squashed state with in the centre to a more spherical at the outer edge of our dataset. Within the radial range probed, the density profile of the stellar halo is well approximated by a single power-law with exponent . We do not find evidence of tilt or offset of the halo with respect to the Galaxy’s disc.

galaxies: individual (Milky Way) – Galaxy: structure – Galaxy: stellar content – Galaxy: halo – stars: varaibles (RR Lyrae)
pubyear: 2017pagerange: The first all-sky view of the Milky Way stellar halo with Gaia+2MASS RR LyraeB

1 Introduction

The diffuse cloud of stars observed around the Milky Way and known as the stellar halo is the alter ego of the much more massive structure, whose presence is inferred indirectly: the dark matter (DM) halo. The DM halo dominates the Galactic mass budget and, according to the currently favoured theories of structure formation, holds clues to a number of fundamental questions in astrophysics. These include, amongst others, the properties of the DM particles (see e.g Davé et al., 2001; Governato et al., 2004; Lovell et al., 2014), the nature of gravity itself (see e.g. Milgrom, 1983; Cabré et al., 2012), as well as the coupling between the DM and the baryons (e.g. Kauffmann et al., 1993; Sommer-Larsen et al., 2003; Chan et al., 2015). The two halos emerge alongside each other, sharing the formation mechanism, i.e. a combination of the accretion onto the Galaxy and the subsequent relaxation and phase mixing. Thus, there ought to exist a bond between them, which can be exploited to reveal the properties of the dark halo through the study of the luminous one. For example, by positing the continuity of the phase-space flow, the DM halo can be mapped out if the stellar halo spatial shape is known and complemented by stellar kinematics (see e.g. Jeans, 1915; Helmi, 2008; Posti et al., 2015; Williams & Evans, 2015a).

Leaving the ideas of James Jeans aside, could the structural parameters of the stellar halo alone inform our understanding of the mass assembly of the Galaxy? At our disposal are the numbers pertaining to the slope of the stellar halo’s radial density profile and its vertical flattening (see e.g. Jurić et al., 2008; Bell et al., 2008; Deason et al., 2011; Xue et al., 2015). The radial profile has so far been measured with a variety of stellar tracers. Studies based on Main Sequence Turn-Off stars (e.g. Sesar et al., 2011; Pila-Díez et al., 2015), Blue Straggler and Horizontal Branch stars (e.g. Deason et al., 2011) and RR Lyrae (e.g. Sesar et al., 2007; Watkins et al., 2009) seem to favour a “broken” profile. According to these datasets, somewhere between 20 and 30 kpc from the MW centre, the density slope changes from a relatively shallow one, as described by power-law index of approximately , to a much steeper one, consistent with a power-law index of .

In an attempt to interpret the observed radial density profile, Deason et al. (2013) conjecture that the presence or absence of a break is linked to the details of the stellar halo accretion history. In their exposition, a prominent break can arise if the stellar halo is dominated by the debris from a i) single, ii) early and iii) massive accretion event. This hypothesis appears to be supported by semi-analytic MW stellar halo models (see e.g. Bullock & Johnston, 2005; Amorisco, 2017b) but is yet to be fully tested with Cosmological zoom-in simulations (see, however Pillepich et al., 2014). Nonetheless, a consistent picture is now emerging: in line with other pieces of evidence, the “broken” MW stellar halo appears to be the tell-tale sign of an early-peaked and subsequently quiescent accretion history. Note also that such destitute state of the stellar halo is not permanent but rather transient (see e.g. Deason et al., 2016; Amorisco, 2017a). In agreement with simulations, the MW stellar halo is destined to transform dramatically with the dissolution of the debris from the Sgr dwarf and the Magellanic Clouds.

While the radial density profile can be gauged based on the data from a limited number of sight-lines through the Galaxy, the shape of the stellar halo requires a much more complete coverage of the sky. So far, much of the halo modelling has relied on the Sloan Digital Sky Survey data, which is biased towards the Northern celestial hemisphere. It is therefore possible that the incomplete view has troubled the efforts to simultaneously infer the details of the radial density evolution and the shape of the halo. For example, using A-coloured stars, Deason et al. (2011) measure substantial flattening of the stellar halo in the direction perpendicular to the Galactic disc plane, but no evidence for the change of the shape with radius. On the other hand, Xue et al. (2015) use a sample of spectroscopically-confirmed K giants to detect a noticeable change of flattening with radius. Furthermore, they argue that if the halo shape is allowed to vary with radius, then a break in the radial density profile is not required. Finally, to add to the puzzle, based on a set of BHB stars with spectra, Das et al. (2016) report both evolving halo shape and the break in the radial density. Pinning down the shape of the stellar halo is important both for the dynamical inference of the shape of the DM halo (see e.g. Williams & Evans, 2015b; Bowden et al., 2016) and for our understanding of the response of the DM distribution to the presence of baryons (see Kazantzidis et al., 2004; Gnedin et al., 2004; Duffy et al., 2010; Abadi et al., 2010).

Looking at some of the earliest halo studies, which inevitably had to rely on much more limited samples of tracers, it is worth pointing out that, strikingly, glimpses of the variation of the halo shape were already caught by Kinman et al. (1966). This pioneering work took advantage of perhaps the most reliable halo tracer, the RR Lyrae stars (RRLs, hereafter). These old and metal-poor pulsating stars suffer virtually no contamination from other populations of the Milky Way and have been used to describe the Galactic halo with unwavering success over the last 50 years (see e.g. Hawkins, 1984; Saha, 1984; Wetterer & McGraw, 1996; Ivezić et al., 2000; Vivas & Zinn, 2006; Catelan, 2009; Watkins et al., 2009; Sesar et al., 2010; Akhter et al., 2012; Soszyński et al., 2014; Torrealba et al., 2015; Soszyński et al., 2016).

While deep, wide-area samples of RRLs now exist, for example provided by the Catalina Sky Survey (CSS, Drake et al., 2013), Palomar Transient Factory (PTF, Sesar et al., 2014) and Pan-STARRS1 (PS1, Hernitschek et al., 2016; Sesar et al., 2017), they have yet to be used to model the Galactic halo globally. In case of CSS, this might be due to the varying completeness of the sample. For PTF and PS1, this is sadly due to the public unavailability of the data. To remedy this, here we attempt to extract an all-sky sample of RRLs from the Gaia Data Release 1 (GDR1, Gaia Collaboration et al., 2016b) data. Our primary goal is to use the thus procured RRL candidates to model the global properties of the MW stellar halo. Therefore, we are not concerned with maximizing the completeness but instead strive to achieve homogeneous selection efficiency and reasonably high purity. While GDR1 does not contain any explicit variability information for stars across the sky, Belokurov et al. (2017) and Deason et al. (2017) show that likely variable objects can be extracted from the GaiaSource table available as part of GDR1. We build on these ideas and combine Gaia and Two Micron All Sky Survey (2MASS) photometry (and astrometry) to produce a sample of RRLs out to kpc from the Sun, with constant completeness of and purity .

Armed with this unprecedented dataset, we simultaneously extract the radial density profile as well as the shape of the Galactic stellar halo. Furthermore, taking advantage of the stable completeness and the all-sky view provided by Gaia+2MASS, we explore whether the density slope and the shape evolve with radius out to kpc. Finally, we also allow the halo to be i) arbitrary oriented ii) triaxial and iii) off-set from the nominal MW centre.

The analysis of the density distribution of the stars in our sample is based on the fit of density models, rather than on the fit of full dynamical models (see e.g. Das & Binney 2016; Das et al. 2016). The main reason behind this choice is that we do not have any kinematic information, so the use of self-consistent dynamical models does not add any significant improvement to our study. Moreover, the knowledge of the spatial density distribution of the stellar halo is a useful piece of information not only if the halo is stationary, but also if it is not “phase-mixed”, as suggested by cosmological -body simulations (Helmi et al., 2011).

The paper is organized as follows. In Section 2 we describe the Gaia data as well as the method used to select an all-sky sample of RRL candidates from a cross-match between Gaia and the 2MASS. Here, we also give the estimates of the purity and completeness of the resulting sample. In Section 3 we show and discuss the spatial distribution of the selected RRLs. Section 4 presents the details of the maximum likelihood approach employed to fit the data with different halo density models and the final results of this analysis. In Section 5 the best-fit halo model is discussed together with the possible biases that can affect our results. The summary of the results can be found in Section 6.

2 The RR Lyrae Sample

In this Section we describe the method used to select a sample of RRLs from GDR1.

2.1 Gaia Data release 1

Gaia is an all-sky scanning space observatory, currently collecting multi-epoch photometric and astrometric measurements of about a billion stars in the Galaxy. More details on the Gaia mission and on GDR1 can be found in Gaia Collaboration et al. (2016a). In the first data release, the information available for most faint sources is limited to basic properties, such as positions on the sky and fluxes in the broad Gaia band, which covers most of the visible spectra from approximately 400 nm to 10000 nm (Jordi et al., 2010).

In this work, we used the table GaiaSource released as part of the GDR1 (Gaia Collaboration et al., 2016b). GaiaSource contains a number of auxiliary pieces of information, which provide plenty of added value to the GDR1. For example, the errors on the mean flux measurements can be used to separate constant and variable sources, and even gauge the amplitude of the variability (see, e.g. Belokurov et al., 2017; Deason et al., 2017). Moreover, the quality of the astrometric fit, encapsulated by the so-called astrometric excess noise, contains information regarding the morphology of the source, and can be used to separate stars from galaxies (see Koposov et al., 2017). The relevant GaiaSource quantities used here (other than the sky coordinates RA, Dec), are:

  • , the number of times a source has crossed a CCD in the Gaia’s focal plane;

  • , the flux (electron per second) measured in the band averaging over single flux measurements;

  • , the standard deviation of the flux measurements;

  • , the mean magnitude in the Gaia band (van Leeuwen et al., 2016) calculated from ;

  • AEN, the astrometric excess noise, which measures strong deviations from the best astrometric solution. The AEN should be large for objects whose behaviour deviates from that of point-like sources, as, for example, unresolved stellar binaries or galaxies (see Lindegren et al. (2012) for details).

Additionally, relying on the cross-match between Gaia and 2MASS, we calculated:

  • PM, the total proper motion of each object. PM, where and are the proper motions measured along RA and Dec, respectively.

2.2 RR Lyrae in Gaia Dr1

We use the parameters provided in GaiaSource to select RRLs from the GDR1. Following the method outlined in Belokurov et al. (2017) and Deason et al. (2017), we defined the quantity


which can be used as a proxy for the amplitude of the stellar variability. Indeed, for variable stars with well sampled light curves, is proportional to the amplitude of the flux oscillation, while for non-variable stars it is just a measure of photon-count Poisson errors. Thus, it is possible to set a threshold value for AMP that will select only variable candidates. For instance, Belokurov et al. (2017) showed that most variable stars like Cepheids and RRLs have .

The selection of variables through the AMP parameter suffers from the limitations that AMP is a time-averaged information and does not allow to distinguish between various types of variable objects: RRLs, Cepheids or Mira variables. However, these different classes of pulsating stars populate a well-defined strip in the colour-magnitude diagram. Therefore, as we show below, one can overcome this problem by applying selection cuts both in AMP and in colour obtaining a fairly clean sample of RRLs. Note that, high values of AMP are also expected for contaminants (e.g. eclipsing binaries) and artefacts (e.g. spurious variations related to Gaia cross-match failures). We discuss their importance in Sect. 2.2.6.

2.2.1 The Gaia + 2MASS sample

Figure 1: Diagnostic diagrams for the selection criteria of our RRL sample. Left-hand panel: distribution of stars in the AMP- space; middle panel: distribution of stars in the colour-magnitude diagram ( magnitude from Gaia, magnitude from 2MASS); right-hand panel: distribution of stars in the AEN- magnitude space. The colour-maps show the distribution of objects in the G2M catalogue, while points show a a randomly selected subsample of bona fide RRLs from GCSS (red circles, 5% of the original sample) and GS82 (orange squares, 35% of the original sample) catalogues (see Sec. 2.2.2). The horizontal-black and vertical-black lines show the selection cuts used to obtain the final sample of RRLs (see Tab. 1), while the green arrows indicate the regions used to obtain the final sample.

As mentioned, GDR1 reports photometric information only in the band. We derived a colour () for each source by cross-matching GaiaSource with the 2MASS survey data (Skrutskie et al., 2006) using the nearest-neighbour method with an aperture of 10 arcsec obtaining the final GaiaSource + 2MASS sample of stars (G2M, hereafter). We chose 2MASS mainly due to its uninterrupted all-sky coverage. The observed magnitudes have been corrected for extinction due to interstellar dust using the maps of Schlegel et al. (1998) and the transformation for the band (see Belokurov et al., 2017) and for the band (Fitzpatrick, 1999).

2.2.2 Auxiliary RR Lyrae datasets

In order to extract a reliable sample of RRL stars from the G2M catalogue, we must apply ad-hoc selection criteria. To this aim, we used two samples of bona-fide RRLs: the CSS (Drake et al., 2013; Drake et al., 2014) and the Stripe 82 (S82, Sesar et al. 2010) catalogues. These samples allowed us to identify the optimal selection criteria, analyse the completeness and the contamination of the catalogue111The completeness indicates the fraction of recovered true RRLs as a function of the apparent magnitude, while the contamination is an estimate of the fraction of spurious objects (non RRLs) that “pollute” our sample. and estimate the RRL absolute magnitude in the band. The CSS contains about 22700 type-ab RRL stars distributed over a large area of the sky (about 33,000 deg between and ) and extended up to a distance of 70 kpc. The completeness of this sample is constant (at ) only for , while it quickly decreases outside this range. Most importantly, as shown in Fig. 13 of Drake et al. (2013), for objects fainter than , the completeness is a strong function of the number of observations and thus varies appreciably across the sky. SDSS’s Stripe 82 covers a 2.5-wide and -long patch of sky aligned with the celestial equator and contains “only” 483 RRLs. However, the sample is very pure (with less than of contaminants) and complete up to a distance of 100 kpc.

The large number of stars in CSS is useful to define the selection criteria (see Sec. 2.2.3) and to estimate the absolute magnitude in the band (see Sec. 2.2.4), while the high quality of S82 sample is ideal to analyse the completeness and contamination of our final sample (see Sec. 2.2.5 and Sec. 2.2.6). A cross-matching of the CSS and S82 catalogues with G2M using an aperture of 1 arcsec led to the two samples CSS+G2M (GCSS hereafter) and S82+G2M (GS82 hereafter).

2.2.3 RR Lyrae selection cuts

In this section we describe how the final sample of RRLs was obtained from the G2M catalogue. The selection was driven by the properties of the bona fide RRLs in the GCSS and GS82 catalogues (see Sec. 2.2.2) in order to maximise completeness of the sample and its spatial uniformity, while keeping the level of contamination low (see Sec. 2.2.5).

First of all, in order to exclude a region likely dominated by the Galactic disc, we removed all the stars in the G2M catalogue located between the Galactic latitudes and . Our limit in latitude () is lower with respect to other works in literature (e.g. Deason et al. 2011; Das et al. 2016), in which, however, the choice was mainly motivated by the limited sky-coverage of the used survey. Given the unprecedented sky-coverage of our data sample, we decided to push forward the study of the halo structure exploring also the region at low Galactic latitude that is usually not well sampled by other surveys. We are aware that our final sample could be polluted by Galactic disc contaminants, but in the following sections we carefully analyse the level of contamination and all our results are obtained taking into account the possible biases due to stars of the Galactic disc.

Figure 1 shows the distribution of the G2M stellar density (yellow-blue-purple colour-maps), a randomly selected subsample of RRLs in GCSS (red points) and in GS82 (orange squares) in the -AMP (left-hand panel), color-magnitude (middle panel) and AEN- (right-hand panel) planes. The bona fide RRLs occupy a well defined strip in colour, thus we excluded all the stars with the colour index greater than and lower than as shown by the vertical black lines in the middle panel of Fig. 1. It is worth noting that most of the “normal” stars occupy this colour interval, therefore this cut mostly eliminates artefacts. The left-hand panel shows that the genuine RRLs are almost uniformly distributed in the AMP range of the G2M sample, however the contamination by spurious objects increases rapidly for AMP (see Sec. 2.2.5 and Fig. 5), thus we only retained stars with variability amplitudes above this value. With regards to completeness, the faint magnitude limit plays an important role. According to our analysis, is the faintest magnitude that we can reach to obtain a sample with spatially uniform completeness (see Sec. 2.2.5 for further details). The number of bright stars with , corresponding to RR Lyrae with distances less than 1 kpc from the Sun, is very small compared to the number of objects in our final catalogue. Therefore, instead of extending our completeness/contamination analysis at the very bright magnitudes (see Sec. 2.2.5 and 2.2.6), we decided to put the bright magnitude limit at .

The selection criteria described above involving colour, AMP and magnitude have the largest impact on the definition of our sample of RR Lyrae. However, we also applied a few minor refinements. The right-hand panel of Fig. 1 shows that most of the bona fide RRLs have a very small value of AEN, so we excluded all sources with as shown by the horizontal-black line. This cut likely removes contaminant extragalactic objects since they typically have (Belokurov et al., 2017) and some of the eclipsing binaries that survive the colour selection. Additionally, to further clean the sample from possible nearby Galactic disc contaminants, we cull all the stars with a total proper motion, PM, greater than  mas/yr. Given differences in the light-curves and its sampling, the significance of AMP (Eq. 1) might depend on the number of photometric measurements . With this in mind, we impose : the focal plane of Gaia has an array of 97 CCDs, so all objects with less than 3 complete Gaia transits are excluded.

It would be useful to have an estimate of the photometric metallicity to retain only genuine metal poor stars from the halo and effectively exclude metal rich contaminants from the Galactic disc. However, the photometric metallicity estimate requires a basic knowledge of the shapes of the RRL light curves which is not available in our dataset (see e.g. Jurcsik & Kovacs 1996).

Finally, we masked a few regions of the sky. First we removed the area near the Magellanic Clouds using two circular apertures: one centred on with an angular radius of for the Large Magellanic Cloud (LMC) and the other centred on with a radius of for the Small Magellanic Cloud (SMC). By inspecting the sky distribution of the stars in our RRL sample, we noticed the presence of two extended structures (S1 and S2, hereafter), that were not connected to any known halo substructures, but are likely objects instead produced by Gaia cross-match failures (see Sec. 2.2.5) that “survived” our selection cuts. We decided to mask these sky regions as well by removing all the stars in the following boxes for S1 and for S2. The final sample has been further cleaned to exclude “hot pixels” using a simple median-filter method. We first built a sky-map using pixels of , then we replaced the number of stars in each pixel by the median of the star-counts calculated in a squared window of four pixels. Finally, we calculated the ratio between the original sky-map and the one processed with the median filter and all the objects in pixels with a ratio larger than 10 were removed. The properties of the median-filter has been gauged to reveal small-scale features, since the most evident large-scale structures have been already removed (LMC, SMC, S1 and S2). The “spotted” hot-pixels correspond to known globular clusters (e.g. M3 and M5) or are connected to “remnants” of cross-match failure structures (see Sec. 2.2.5). In order to fully exploit the all-sky capacity of our sample, we decided not to exclude a priori any portion of the sky containing known substructures (e.g. unlike Deason et al. 2011 which masked out the Sagittarius stream). The analysis of the most significant substructures found in this work can be found in Sec. 3 and 5.

The final sample contains about 21600 RRLs that can be used to have a direct look at the distribution of stars in the Galactic halo (Sec. 3). The final number of object is similar to the one in the CSS catalogue, however we cover a larger area of the sky (almost all-sky). Our sample populates 58% of the halo spherical volume within the Galactocentric distance of about 28 kpc which represents a significant improvement in volume fraction as compared to previous works (e.g. 20% in Deason et al. 2011). A summary of the applied selection cuts can be found in Tab. 1.

Selection cuts
|b| [deg]
[mag] (10, 17.1)
AMP (-0.7, -0.4)
(-0.95, -0.4)
PM [mas/yr]
|| [deg]
Structure cuts
21643 (13713)
[%] 58 (44)
Table 1: Summary of the selection cuts used to obtain the final sample of RRLs from the G2M catalogue. The description of parameters used in the sample selection and the details on the cut substructures can be found in Sec. 2.1. and are the sky angular distances from the LMC and SMC, respectively. is the Galactocentric latitude (defined in Eq. 7) and was estimated by assuming an RRL absolute magnitude of . The refers to the subsample used in the likelihood analysis in Sec. 4. The bottom part of the table gives a summary of the whole sample. is the number of stars in the sample and is the fraction of the spherical volume of the halo sampled by our stars between Galactocentric distance 0 kpc and 28 kpc.

2.2.4 Absolute magnitude and distance estimate

Figure 2: Distribution of band absolute magnitudes for bona fide RRLs in the GCSS catalogue (Gaia+CSS, see Sec. 2.2.2). The coloured lines show the fit with different functions: single Gaussian (red) and double Gaussian (blue, with blue dashed lines showing individual Gaussian components of the fit).

Despite photometric variability, RRLs have an almost constant absolute magnitude and, having the apparent magnitudes , we can directly estimate the heliocentric distances through


where is the absolute magnitude in the Gaia band. To estimate we used the RRLs in GCSS, as they have heliocentric distances estimated from the period-luminosity relation. The resulting distribution is shown in Fig. 2: it has a well defined peak at and a small dispersion. We fit this distribution with a Gaussian function obtaining a mean of and a dispersion of , comparable to the uncertainties of the V-band absolute magnitude of RRLs (see e.g. Vivas & Zinn 2006). The above Gaussian function perfectly describes the data in the central parts, but the distribution shows broader wings for and . The objects that populate the wings could be Oosteroff type II RRLs and objects influenced by the Blazhko effect (Drake et al. 2013 and references therein). A better fit can be obtained using a double Gaussian model, where two Gaussian components peak at about the same absolute magnitude of the single Gaussian fit (see Fig. 2). Given the small dispersion around the mean, we decided to set the absolute magnitude for all the RRLs in our sample to a single value . An analysis of the effect of this approximation on the results can be found in Sec. 5.2.1.

2.2.5 Completeness

Figure 3: Completeness analysis of the sources in GaiaSource (Sec. 2.1) as a function of Galactic latitude. Top panel: regions of the sky considered in this analysis, each stripe has a constant Galactic longitude. Bottom panel: ratio between the number of stars in the GaiaSource () and stellar sources in the SDSS DR7 (, Abazajian et al. 2009) in bins of Galactic latitude. The lines refer to the ratio obtained for stars located in regions of the same colour shown in the top panel. The dots and the black line indicate the ratio obtained considering the stars in all the “stripes”. The stars in Gaia have been selected using the cut, and the SDSS sources using the cut (further details on the text).

Before studying the properties of the stellar halo (as traced by RRLs), it is fundamental to consider the completeness of our sample of RRLs.

First of all, we checked that the scanning law of Gaia does not cause an intrinsic decrease of the completeness at low Galactic latitude. We compared the number of objects in GaiaSource with the number of stellar sources in the Data Release 7 of the Sloan Digital Sky Survey (SDR7, Abazajian et al. 2009) in a series of stripes at fixed Galactic longitudes, selected using the footprint of SDR7. The top panel of Fig. 3 shows the position of the stripes in Galactic coordinates. We selected all stellar sources in the SDR7 with , where is the -band magnitude. In this magnitude range the SDR7 can be considered 100% complete222 We chose the band because the peak of the filter response is almost coincident in wavelength with the one of the Gaia band (see Sec. 2.1), therefore the two magnitudes are directly comparable. The bottom panel of Fig. 3 shows the ratio between the number of stars in Gaia and SDR7 in bins of Galactic latitude for individual SDSS stripes and considering all the stripes together. The ratio does not show significant variations as a function of b with values between 0.9 and 1.0. We conclude that there is no evidence for strong intrinsic completeness variations in the Gaia catalogue for .

We estimated the completeness (and the contamination) using the RRLs in our auxiliary catalogues: CSS and S82 (Sec. 2.2.2). In particular, S82 represents a complete (100 %) and pure (99 %) catalogue of RRLs located up to 100 kpc from the Sun, so it is perfect to test both the contamination and the completeness of our sample. We compared the number of stars in the original S82 sample with the ones contained in the GS82 after the selection cuts described in Sec. 2.2.3, in bins of heliocentric distance. Fig. 4 shows the level of completeness as a function of magnitude/distance assuming different lower limits in AMP ( red triangles, blue circles, green diamonds) in the range of magnitudes . The results are in agreement with the distance-based estimate of completeness as shown by the dashed lines. The level of completeness is relatively low, ranging from about 15% for to about 30% for but it is reasonably constant up to about then it abruptly decreases to 0 at so we decided to conservatively cut our sample at (vertical black line in Fig. 4, see Sec. 2.2.3).

Figure 4: Completeness of our samples of RRLs as a function of distance from the Sun () or magnitude. The conversion from to has been obtained from Eq. 2 assuming a constant absolute magnitude . Different symbols indicate completeness for samples obtained using different AMP cuts: red-triangles for , blue circles for and green diamonds for The error bars were calculated using the number of stars in each magnitude range and Poisson statistics. The dashed lines indicate the completeness estimated in the magnitude range of , while the vertical black line marks the faint magnitude limit of our final sample (see Sec. 2.2.5).
Figure 5: Contamination of the RRLs sample as a function of the AMP cut. The contamination is estimated using the S82 catalogue (see Sec. 2.2.2) in the range of magnitudes . The error-bars have been calculated by propagating Poisson uncertainties from the number of stars in AMP bins.

We also used the GCSS catalogue and found no significant variation of the completeness as a function of the Galactic sky coordinates l and b, although we found a mild trend for increasing . This is expected since the larger the number of flux measurements the better the sampling of the light curves and, as a consequence, the AMP cut (Eq. 1) improves its effectiveness in selecting RRLs. However, the increase is not dramatic as all the differences are within 10%. In conclusion, for the purpose of this work, we considered the completeness of our catalogue of RRLs uniform across our sky coverage and the considered magnitude range (see Tab. 1).

2.2.6 Contamination

We estimated the contamination of spurious sources as


where and indicate the number of stars in our samples and in GS82 catalogue after the selection cuts in Sec. 2.2.3. The S82 sample is pure (Sesar et al., 2010), so all the stars in our sample that are not present in the GS82 catalogue are likely contaminants. As before, we considered the magnitude range for different lower limits of AMP. For AMP threshold lager than the level of contamination is lower than 10% (Fig. 5) with a mild increase toward low Galactic latitudes, where the contribution of the Galactic disc is larger. For AMP cut lower than the contamination fraction rapidly increases (about 25% for ). We expect that the contaminants of our sample are eclipsing-binaries in the Galactic disc and possible instrumental artefacts. As shown and discussed in detail in Belokurov et al. (2017) some regions of the Gaia all-sky map are crossed by sharp strips with a large excess of objects. These strips are similar in width to the field of view of Gaia and are regularly spaced on the sky. Belokurov et al. (2017) propose that these are spurious features due to failures in the cross-match procedure of Gaia, so that at some epochs the flux of a star comes from a different object. The spurious measurement of the flux increases and moves the stars in the region of high AMP “polluting” our samples. We found that for the contaminants are mostly related to cross-match failures. Unlike the disc contaminants, the cross-match failures have a complicated and poorly understood spatial distribution, so these structures are difficult to be taken in account in the study of the properties of the Galactic halo. For this reason, we decided to use as the lower limit in our selection cut (see Sec. 2.2.3).

3 Distribution of the RRLs in the inner stellar Halo

3.1 A first Gaia look at the stellar halo

Figure 6: Star-count all-sky maps in the Galactic coordinates (l, b) for the RRLs in our sample (Sec. 2.2.3) in the magnitude intervals (upper panel) and (lower panel). The “holes” between and are due to the mask used to eliminate the contribution of the Large and Small Magellanic Clouds (see Sec. 2.2.3).

Fig. 6 shows the distribution of the RRLs of our final sample as a function of the Galactic coordinates (l, b). We split the stars in the sample into magnitude intervals: , corresponding to heliocentric distances between 1 kpc and 10 kpc (Eq. 2), and , i.e. with heliocentric distances between 10 kpc and 20 kpc (assuming that all RRLs have absolute magnitude ). We stress that Fig. 6 represents the first all-sky view of the distribution of the RRLs in the inner halo.

The first magnitude range covers a portion of the Galaxy mostly located in the side of the Galaxy containing the Sun between the Galactic radii 1.4 kpc and 18 kpc. The distribution of stars in this region is quite regular: as expected most of the stars are in the direction of the Galactic centre () and there are not evident asymmetries with respect to the Galactic plane. The second magnitude range covers the Galactic distances between 18 kpc and 28 kpc in the side of the Galaxy containing the Sun and between 3 kpc and 12 kpc in the other side. In this distance range the distribution of the RRLs is less regular with clear asymmetries: in particular, an excess of stars is evident at high Galactic latitude around . In both magnitude intervals an over-dense band of stars at low Galactic latitude () can be seen running all around the Galaxy. A discussion of the nature of these structures can be found in Sec. 5.3.

3.2 Setting the frame of reference

We set a left-handed Cartesian reference frame centred in the Galactic centre and such that the Galactic disc lies in the plane , the Sun lies on the positive X-axis and the Sun rotation velocity is . The actual vertical position of the Sun with respect to the galactic disc is uncertain, but it is estimated to be be smaller than 50 pc (Karim & Mamajek 2017 and reference therein), and thus negligible for the purpose of this work. In this Cartesian reference frame, the coordinates of an object with Galactic longitude l, Galactic latitude b and at a distance from the Sun are


where is the distance of the Sun from the Galactic centre. In this work we assume (Bovy et al., 2012), but the main results of our work are unchanged for other values of , within the observational uncertainties (between 7.5 kpc, e.g. Francis & Anderson 2014, and 8.5 kpc, e.g. Schönrich 2012).

For a star with Galactocentric Cartesian coordinates we define the distance from the Galactic centre


the Galactocentric cylindrical radius


the Galactocentric latitude


and the Galactocentric longitude


Assuming that the stellar halo is stratified on concentric ellipsoids, it is useful to introduce another Cartesian frame , aligned with the ellipsoid principal axes, and to define the elliptical radius


where p and q are, respectively, the Y-to-X and Z-to-X ellipsoid axial ratios. In general we will allow the to differ from in both orientation and position of the origin. When the origin of is the Galactic centre and , as we will assume in this section, the system can be identified with without loss of generality.

3.3 Density distribution of the halo RRLs

Given the all-sky RRL sample illustrated in Fig. 6, we can compute their volume number density distribution . In particular we define the number density of halo RRLs in a cell centred at , where is a coordinate vector, such as for instance , as


where is the number of stars observed in the cell, V is the volume of the cell and is the fraction of this volume accessible to our analysis, which depends on the sky coverage, on the selection cuts and the mask that we applied (see Sec. 2.2.3), and on the off-set between the Sun and the Galactic centre. We estimate numerically as follows. First, we define a Cartesian Galactic grid sampling each of the axes with 300 points linearly spaced between and ( kpc is the maximum value of , given our magnitude limit at ), so the reference Galactic volume (a cube of 56 kpc side centred on the Galactic centre) is sampled with 27 million points with a density of about 125 points per kpc. Each of the points of the grid is a “probe” of the Galactic volume. We assign to each of them Galactic coordinates, heliocentric coordinates and observational coordinates (l, b, , ). A secondary (non-uniform) grid is built by applying the selection cuts in Tab. 1 to the points on the primary grid, so we end up with a grid sampling the complete Galactic volume and a grid sampling the portion of the volume accessible to our analysis. Given a cell in a certain coordinates space (e.g. R, ), we define as the ratio between the number of points of the secondary grid and the number of points of the primary grid contained in the cell. In summary, when we estimate the number density in a 1D or 2D space we build three binned maps: the first contains the number of stars in each bin/cell (), the second the total volume of each bin/cell (V) and the third the fraction of this volume that we are actually sampling (). These three quantities are inserted in Eq. 10 to estimate the stellar number density. We tested this method with both simple analytic distributions and mock catalogues (see Sec. 4.3.5).

3.3.1 Meridional plane

Figure 7: Number density in the Galactic plane for the RRLs in our sample. The number density is calculated dividing the number of stars found in a cylindrical ring by the volume of the ring corrected for the non-complete volume sampling of the data (Eq. 10). The black contours are plotted with spacing from to , where is in units of ; the white dashed lines represent elliptical contours with an axial ratio .

Fig. 7 shows the number density of the RRLs of our sample in the Galactic meridional plane . The shape of the iso-density contours clearly shows the presence of two components: a spheroidal component and a discy component. The discy component causes the flattening of the contours at low . The nature of this component is uncertain: it could represent the disc RRL stars, it could be a low-latitude substructure of the MW halo or, finally, it could be due to non-RRL contaminants from the Galactic disc (both genuine variables such as eclipsing binaries as well as artefacts, e.g. due to cross-match failures). A more detailed analysis of this low-latitude substructure can be found in Sec. 5.3. The density at higher is less contaminated by the discy component and it represents more directly the density behaviour of the RRLs in the halo. The iso-density contours nicely follow the the elliptical contours (overplotted in Fig. 7) out to kpc. At larger radii the contours tend to be rounder. The overall density distribution looks reasonably symmetric with respect to the Galactic plane, although there is an over-dense region at kpc, which does not seem to have a counterpart below the Galactic plane (see Sec. 3.3.4 for further details).

3.3.2 Density profile

Figure 8: 1D number density profiles of the halo RRLs as functions of the elliptical radius (Eq. 9) assuming that the halo is stratified on spheroids with and different values of q: (blue diamonds), (green dots) and (red squares). The black-dashed line shows, for comparison, a single power law with index . The Poissonian errors are smaller than the size of the symbols. The density normalisation is arbitrary and different for each profile.

Fig. 7 demonstrates that the RRLs in the inner halo are consistent with being stratified on spheroids (except in the region close to the Galactic plane), thus we estimate the 1D RRL density profile by counting the stars in spheroidal (, ) shells and dividing this number by the shell volume corrected by the coverage of our sample (Eq. 10). The profiles are shown in Fig. 8 for , and . Independently of the assumed value of q, the density of RRLs follows a single power law with no significant evidence of change of slope out to an elliptical radius of 35 kpc. It must be noted that the most distant region is only sampled by a low number of stars located in a small portion of the total volume, so it is possible that a change of slope starting near the edge of our elliptical radial range could be missed with this analysis. As the elliptical radius depends on the assumed axial ratios, the results of this analysis are not straightforwardly comparable with previous works (see Sec. 5.4 for a detailed comparison with results obtained in the literature).

The slope of this power law is similar in the three cases and very close to (see Fig. 8), but, based on this simple analysis, all values of in the range are consistent with the data.

3.3.3 Halo flattening

Figure 9: Number density of RRLs in the elliptical radius () - Galactocentric latitude () space for (left-hand panel), (middle panel) and (right-hand panel), assuming . The density is normalized to the maximum value of each panel. The contours show the normalised density levels , while the dashed lines indicate .

The distribution of RRLs in the meridional Galactic plane (Fig. 7) suggests that the inner halo should be reasonably well represented by a spheroidal stratification with . However, in order to more rigorously study the halo flattening, we estimate the density in the plane (see equations 4 and 9). In practice, for a given value of , we scan the density as a function of at fixed . If the RRLs are truly stratified on similar spheroids with axial ratio , the density is independent of , so the isodensity contours in the plane are vertical stripes. If the assumed value of q is smaller than the true value , is underestimated, the estimated density is a monotonic decreasing function of and the isodensity contours in the plane are bent in the direction of increasing with (provided that the density is a decreasing function of ). The isodenisty contours are bent in the opposite direction, if one assumes . The shape of the iso-density contours in the plane is a very efficient and direct diagnostic of the evolution of q as function of the elliptical radius.

The number density maps in the plane of the RRLs in our sample are shown in Fig. 9 for , and , assuming . Below (indicated by the white-dashed lines) the contours are nearly horizontal, because the density is dominated by the highly flattened discy component (see Fig. 7) for which q is over-estimated in all the panels. At higher latitudes () the contours give a direct indication on the flattening of the halo: in the right-hand panel of Fig. 9 () the iso-density contours are significantly inclined in a way that implies . For kpc a flattening of about 0.55 gives a good description of the data as shown by the vertical iso-density contours in the middle panel of Fig. 9 (), but beyond 20 kpc the contours start to bend so that . The last two iso-density contours in the left-hand panel () look vertical enough to assert that at the outer radii the halo becomes more spherical. This analysis, together with the recent works of Liu et al. (2017), shows the first direct evidence of a change of shape of the stellar halo going from the inner to the outer halo. Moreover the unique all-sky view of our sample allows us to confirm that this trend is symmetric with respect to the Galactic plane. A variation of the halo flattening was also proposed in previous works (e.g. Xue et al. 2015 and Das et al. 2016). In Sec. 4.4 we perform a comprehensive model fitting analysis of the RRL dataset and compare our results to those found in the literature.

3.3.4 Vertical asymmetries

Figure 10: Number counts of stars in the plane for different -slabs: left-hand panels with , middle panels with and right-hand panels with . The top and bottom panels show the layers above () and below () the Galactic plane, respectively. The colour maps are normalised to the maximum number in each column of plots. is the total number of stars found on each -slab. The black dots show the position of the Sun (at ), while the crosses indicate the maximum of the star counts. The maps have been smoothed with a spline kernel.

Contours of the RRL density shown in Fig. 7 and Fig. 9 display an overall symmetry between the Northern and Southern Galactic hemispheres, however at kpc and around an over-density of RRLs above the disc plane is evident. Thanks to the unprecedented sky coverage of our sample we can directly analyse the distribution of stars in different -slabs to understand whether the over-density is compatible with the Poisson star-count fluctuations or it is due to a genuine halo substructure.

Fig. 10 shows star-count maps in the plane integrated along three different -slabs. Close to the plane (left-hand panels) the two maps look very similar and the difference between the total number of stars is small (less than 2%) and compatible within the Poisson errors. In these maps, an elongated component (also visible in the other density maps, Figs. 7 and 9) can be seen stretching from the Galactic centre to the Sun and beyond: this component appears symmetric with respect to the Galactic plane but is present only in the side of the Galaxy containing the Sun. In the intermediate -slabs (middle panels) the difference between the number of stars in the region above and below the Galactic plane is small (less than 5%) as in the previous case, but the maps look less symmetric. In particular the peaks of the star counts are now apart by about 10 kpc. Finally, the slabs at the highest (right-hand panels) show a significant difference in the total number of stars (about 20%) that can not be explained by Poisson fluctuations only. Moreover the excess of stars above the Galactic plane is strongly clustered in the regions between kpc and kpc.

Note that some of the differences between the star counts in the regions above and below the Galactic plane could be due to the mask used to eliminate the contribution of the Magellanic Clouds (see Sec. 2.2.3). Indeed, some mismatch is expected given that we have excluded a relatively large region of the halo volume below the Galactic plane. In order to quantify the differences introduced by the mask, we produce mock catalogues (see Sec. 4.3.5) of different halo models and find that the Magellanic Clouds mask introduces a difference of about the 7% in the number of objects above and below the Galactic plane for kpc, significantly less than the value of 20% obtained here. Therefore, the structure seen at high positive appears to be genuine, and is most likely related to the “Virgo over-density” (Jurić et al., 2008; Vivas et al., 2016). Across all slabs, portions of the Virgo Cloud are visible as an excess of stars at positive . Virgo’s counterpart underneath the disc is the Hercules-Aquilla Cloud (see Belokurov et al., 2007; Simion et al., 2014), discernible in middle bottom panel as strong over-density at negative (and ).

4 Model fitting

In this section we present the stellar halo models and we compare them with the observed sample of RRLs.

4.1 Clean sample

Fig. 7 and Fig. 9 show the presence of a highly “flattened” structure close to the Galactic plane. The properties of this structure are clearly at odds with the distribution of stars at high Galactic latitude that are more likely a “genuine” tracer of the halo population. The “flattened” component contains about 35% of the RRLs in our sample, so it must be taken into account to infer the properties of the stellar halo. Therefore, we built a “clean” sample of RRLs eliminating all the stars belonging to the substructure from our original catalogue (Sec. 2.2.3). Fig. 9 suggests that the substructure can be effectively eliminated with a selection in angle (see Eq. 7).

In particular, at there is a transition between a very flattened component and a more spheroidal structure. Therefore, we define our clean sample of halo RRLs as the stars with : this subsample contains 13713 objects and covers approximately 44% of the halo volume within a sphere of radius 28 kpc (see Tab. 1). We stress that this cut is based on values of obtained from Eq. 2 and Eq. 7 assuming the same value of the absolute magnitude, , for all the stars (see Sec. 2.2.4).

We also tried to exclude the flattened structure with alternative cuts, e.g. using a higher cut on the Galactic latitude b or a direct cut on . We found that the results obtained for the sample with (see Sec. 4.4) are qualitatively similar to the results obtained using a sample with (the lowest Galactic latitude used in most of the previous works, e.g. Deason et al. 2011) or a sample with kpc (Fermani & Schönrich 2013 uses kpc to cut disc stars). The cuts on b and significantly reduce the number of tracers in the inner part of the halo and the final number of stars is smaller, in both cases, with respect to the one obtained cutting the original sample at .

In the following subsections we present the method used to fit halo models to this subsample of RRLs and the final results.

4.2 Halo models

As in Sec. 3.1, we assume here that the number density of the halo RRLs is stratified on ellipsoids with axial ratios p and q. Therefore, a halo model is defined by a functional form for the density profile (number density as a function of the elliptical radius) and a “geometrical” model for the ellipsoidal iso-density surfaces (in practice, characterized by values of p and q, and the orientation of the principal axes with respect to the Galactic disc).

4.2.1 Density profiles

We consider five families of number density profiles: double power-law (DPL), single power-law (SPL), cored power-law (CPL), broken power-law (BPL) and Einasto profiles (EIN). The DPL profile has a number density


where and indicate the inner and outer power-law slopes and is the scale length. The number density of the SPL is given by Eq. 11 when , while the number density of the CPL has , in which case represents the length of the inner core. The BPL number density profile is given by a piecewise function:


The EIN profile (Einasto, 1965) is given by


where for (Graham et al., 2006). The steepness of the EIN profile, , changes continuously as a function of tuned by the parameter n,


The EIN profile is the only density law among those considered here that assures a halo model with a finite total mass for any choice of the parameters. The power-law profiles with or with imply halos with infinite total mass, however our study focuses only on a limited radial range and we do not exclude a priori any solution. We anticipate that our best density model is a SPL with (see Sec. 4.4), therefore there should be a physical radius, outside our radial range, beyond which the profile becomes, either abruptly (e.g. BPL or an exponential truncation) or gently (e.g. DPL), steeper.

4.2.2 Iso-density ellipsoidal surfaces

Concerning the iso-density ellipsoidal surfaces we define four different models:

  • spherical (SH): we set and in Eq. 9, so that is just the spherical radius (see definition 5) in the Galactic frame of reference;

  • disc-normal axisymmetric (DN): we set , the axis of symmetry is normal to the Galactic disc, and q is a free-parameter;

  • disc-plane axisymmetric (DP): we set , the axis of symmetry is within the Galactic plane making an (anticlockwise) angle with respect to the Galactic Y-axis, and p is a free-parameter,.

  • triaxial (TR): both p and q are considered free-parameters, the Z-axis of symmetry is coincident with the normal to the Galactic plane and the X and Y axes are within the plane making an (anticlockwise) angle with respect to the Galactic X and Y axes, respectively.

Given the unprecedented sky coverage of our sample, we also tested more complex models for the iso-density ellipsoidal surfaces:

  • q-varying (qv): q depends on the elliptical radius as


    so q varies from at the centre to the asymptotic value at large radii and the variation is tuned by the exponential scale length ;

  • tilt (tl), in this model we assume that the principal axes of the ellipsoids are tilted with respect to the Galactic plane, in practice, before calculating the elliptical radii (Eq. 9), we transform the Galactic coordinates (see Eq. 4) of each star by applying a rotation matrix following a ZYX formalism so that is the rotation angle around the original Z-axis, the one around the new Y-axis and finally is the rotation angle around the final X-axis, all the rotation are defined in the anticlockwise direction;

  • offset (off), in this model the elliptical radius of each star is estimated with respect to a point offset by (, , ) with respect to the Galactic centre.

The specific functional form of Eq. 15 is empirical: the same expression was adopted by in Xue et al. (2015) and Das et al. (2016), where, however, q is a function of the spherical radius (). We decided to maintain the dependence on given that this approach is self-consistent with our assumption that the RRLs are stratified on ellipsoidal surfaces. If follows that for our models with a varying q, the elliptical radius of a star with coordinates is the root of


where is defined in Eq. 15. In our fitting code we solve this equation numerically with a Newton-Raphson root finder.

4.2.3 Complete halo model

Each complete halo model is defined by a model for the density law (SPL, DPL, CPL, BPL, EIN) plus a model for the shape of the iso-density surfaces (SH, DN, DP, TR) and any combination of geometrical variants (qv, tl, off). For instance, a SPL-SH model has a spherical distribution of stars with a single power-law density profile, while a model is a triaxial tilted model with varying flattening along the Z-axis of symmetry and an Einasto density profile (see Tab. 2 for a reference on the various model labels). In the next sections when we discuss the properties or the results of density models (e.g. SPL models), unless otherwise stated, we are implicitly referring to all the models that share the same density model whatever the geometrical and geometrical variants properties are. The same applies when we focus on geometrical models only. We fit our data with all possible combinations of density laws (SPL, DPL, CPL, BPL, EIN), geometrical models (SH, DN, DP, TR) and model variants (qv, tl, off). See Tab. 3 for a summary of results obtained with a sample of such models.

4.3 Comparing models with observations

4.3.1 Density of stars in the observed volume

Given a certain halo model (Sec. 4.2), the normalised number density of stars in an infinitesimal volume is given by


where is the vector containing all the model parameters (e.g. for an SPL+disc-normal axisymmetric model, see Sec. 4.2) and is defined such that


It is useful to define the normalised star number density


within the infinitesimal projected (onto the sky) volume element centred at (,l,b), where is the observed magnitude, l and b are the Galactic longitude and latitude, and


is the determinant of the Jacobian matrix (Appendix A). The number density depends also on the additional parameter , which is the absolute magnitude needed to pass from the observable variables to the Cartesian variables (see Eqs. 2 and 4). Substituting Eq. 20 in Eq. 19 we obtain


4.3.2 Likelihood of a single star

From the normalised density function (Eq. 21) we can define the expected rate function for finding a star with (, l, b) given the absolute magnitude and a halo model with parameters


In Eq. 22, is a constant, represents the probability density function (pdf) of the absolute magnitude of the stars, while the function is the selection function that takes into account the incomplete coverage of the Galactic volume. It is a function of l, b and and returns a result in the Boolean domain . In particular, is always equal to 1 except for the points (, l, b) that are outside the volume covered by our clean sample (see Tab. 1).

Density parameters Iso-density surfaces parameters
/n (kpc) p q() (kpc) (kpc) (degree)
Single power-law SPL U[0,20]
Cored power-law CPL U[0,20] U[0.01,10]
Double power-law DPL U[0,20] U[0,20] U[,]
Broken power-law BPL U[0,20] U[0,20] U[,]
Einasto EIN U[0,100] U[0.1,500]
Spherical SH
Disc-normal axsym DN U[0.1,10]
Disc-plane axsym DP U[0.1,10] U[-80,80]
Triaxial TR U[0.1,10] U[0.1,10] U[-80,80]
q-var qv U[0.1,100]
Offset off U[0,10]
Tilted tl U[-80,80]
Table 2: Prior distribution (Eq. 28) of the halo model parameters. Each row refers to a given component of the halo models, each column indicates a single parameter or a group of parameters if they share the same prior distribution. The second column indicates the labels where we indicate the halo models through the text (Sec. 4.2.3). The third column refers to the parameter n for the Einasto profile and to the parameter for the other density models. The U indicates a uniform distribution within the values shown inside the square brackets; a indicates a Dirac delta distribution. The prior range of the parameter for the BPL and DPL density models is not the same in all models, depending on the minimum () and maximum () elliptical radius of the stars in the sample.

In this work, we decided to set the absolute magnitude of the RRLs to a single value, (Sec. 2.2.4), so we are assuming that in Eq. 22 is a Dirac delta. As a consequence, we can marginalize Eq. 22 over to obtain


and we define the normalised rate function as


The normalised rate function in Eq. 24 is the pdf of stars, i.e. the likelihood per star, at a certain position (, l, b) for halo model parameters .

4.3.3 The total likelihood

Consider a sample of stars with coordinates (, l, b) and a halo model with parameters . Given the pdf of the stellar distribution (Eq. 24) the logarithmic likelihood is


where is the number of stars in our sample. Plugging Eq. 21 and Eq. 23 into Eq. 24 we can write the logarithmic likelihood as




is the normalisation integral. Notice that the numerator of Eq. 26 is evaluated only in regions of the sky where (i.e. where we observe stars).

4.3.4 Sampling of the parameter space

Considering the Bayes’s law, the logarithmic posterior probability of the parameters of an halo model given the data is


where is the probability of the data, , given the parameters (See sec. 4.3.3) and represents the prior probability of (Tab. 2). In Eq. 28 we omit the Bayesian evidence term that is defined as the integral of the likelihood over the whole parameter space. This term is negligible in the determination of the best set of parameters for a given model (see below).

Using Eq. 26 we can write the posterior probability as


In order to explore the parameter space without biases, we decided to use very large priors for most of our parameters (see Tab. 2). In particular, the very large ranges for and approximate uniform infinite priors.

Concerning the parameter , we assume, for both the BPL and DPL density models, that the prior distribution is uniform between and which represent the minimum and maximum elliptical radii (Eq. 9) of stars in our sample. For most of our models, remains constant at about 2.4 kpc, while ranges between roughly 31 kpc and 38 kpc.

We analysed the posterior probability of the parameters, , for a given halo model using the Eq. 29 and sampled the parameters space with Goodman & Weare’s Affine Invariant Markov chain Monte Carlo (mcmc, Goodman & Weare 2010), making use of the Python module emcee333 (Foreman-Mackey et al., 2013). The technical details of our approach can be found in Appendix B.

The final best-fit values of the model parameters have been estimated using the 50th percentile of the posterior distributions and the 16th and 84th percentiles have been used to estimate the 1- uncertainties.

4.3.5 Tests on mock catalogues

In order to test our fitting method, we developed a simple Python script to generate mock catalogues with the same properties (e.g. number of stars, magnitude limits) as our clean sample of RRLs (Sec. 4.1). This algorithm distributes the stars using a combination of the density laws (Sec. 4.2.1) and according to the assumptions on the properties of the iso-density ellipsoidal surfaces (Sec. 4.2.2). The absolute magnitude of the stars in the mock catalogues are distributed using the best-fit double Gaussian functional form shown in Fig. 2. The final mock catalogues do not include Gaia and 2MASS photometric and sky-coordinates errors because the uncertainties they cause on the estimate of the distance are negligible.

We applied our fitting method to mock catalogues finding that we are able to recover the input parameters for all possible model combinations, including in the sample a fraction of Galactic disc contaminants (see Sec. 5.2.3).

Analysing the mock catalogues, we can estimate what properties of the halo we are able to measure and constrain using the stars in our catalogue. We found that we are able to detect a core in the density profile (CPL density model) all the way down to 100 pc. However for a very small core size most of the information comes from a relative small region where the transition between the inner core and the outer density profile takes place. The detection of a change of slope depends on the halo flattening: for and we can detect a break in the density profile within about kpc (Eq. 21). For more flattened models (q around ) this range extends up to about kpc. However, in this latter case most of the information for regions where kpc comes from a small number of stars at very high Galactic latitude, therefore the fit can be easily biased by the presence of substructures such as the one highlighted in Sec. 3.3.4 and Fig. 9. In conclusion, we are confident that we are able to robustly detect significant deviations from a SPL within a range of elliptical radii from less than 1 kpc to about 30 kpc. Finally, we verified that the method described above can recover a variation of the halo flattening given a realistic mock dataset.

4.4 Results

In this section we present the main results obtained applying the method described in Sec. 4.3 to the RRLs in our clean sample (Sec. 4.1).

Density law Surface Parameters
SPL SH -868 1688
-405 779
-81 124
, ,
-79 141
, ,
-81 143
, kpc,
-81 134
, kpc,
-87 145
, ,
-66 113
-22 23
, ,
0 0
Table 3: Summary of properties and results for a sample of families of halo models. For each family we report the assumed density law (Sec. 4.2.1) and geometry of the iso-density surfaces (Sec. 4.2.2). See Tab. 2 for a reference on the model labels. For each fitted parameter we report the median and the uncertainties estimated as the 16th and 84th percentiles of the posterior distribution, the values in parentheses indicate the value for which we obtain the maximum value of the likelihood (Eq. 26). The last two columns indicate the logarithmic likelihood and BIC differences between the best model of the family and the best of of all the presented models. In the SPL models , see Sec. 4.2.1. The (anticlockwise) angle indicates the tilt of the X and Y axes of symmetry of the halo with respect to the Galactic X and Y axes, see Sec. 4.2.2.

4.4.1 Model comparison

In addition to estimating the parameters for a single model, it is important to compare the results of different models to determine which of them gives the best description of the data. The most robust way to perform a model comparison is through the ratio of the Bayesian evidences. Under the assumption that the posterior distributions are almost Gaussian, the Bayesian evidence can be approximated by the Bayesian information criterion (BIC, Schwarz 1978) defined as


where is the number of objects in the sample and is the maximum value of the likelihood (Eq. 26).

The BIC is often used to compare different models with different dimensions in parameter space and the model with the lowest BIC is preferred. The BIC is similar to the maximum likelihood criterion, but it includes a penalty depending on the number of free parameters, such that for two models with the same likelihood, the one with more parameters is penalised. The best-fit parameters together with a comparison of the maximum likelihood and the BIC for a representative sample of halo models can be found in Tab. 3. These quantities have been calculated taking the maximum of the likelihood estimates obtained with the final mcmc sampling (see Sec. 4.3.4).

Other than the BIC, we also compared the ability of the different models to describe the observed properties of the RRLs (e.g. distribution of the stars in the sky).

4.4.2 RRLs density law

Figure 11: Top panel: the black points show the fraction of stars in magnitude bins, while the curves show the same fraction expected for halo models with different density laws: SPL (blue), BPL (orange), DPL (green), CPL (red) and EIN (magenta). Error bars on the data points and on the model distributions indicate Poisson uncertainties. Bottom panel: relative residual (Data-Model/Model). In all the cases we assumed a DN geometrical halo model (see Tab. 3).
Figure 12: Radial profile of the density slope ()) for different density law models: BPL (1th-row), DPL (2th-row), CPL (3th-row) and EIN (4th-row). The bands represent the posterior distribution between 16th and 84th percentiles. The coloured lines represent the median of the posterior obtained fitting only the stars above (, dotted line) and below (, solid line) the Galactic plane. The dashed-black line shows the best-fit slope obtained for the SPL model. In all the cases we assumed a DN geometrical halo model (Tab. 3).

We tested different density profiles assuming a DN geometrical model: the SPL (Eq. 11 with ), the DPL (Eq. 11), the BPL (Eq. 12), the CPL (Eq. 11 with ) and the EIN profile (Eq. 13).

The logarithmic maximum likelihoods and the BICs obtained for these different models are very similar as shown in Tab. 3, moreover Fig. 11 shows that the predicted fraction of RRLs in bins of magnitude, for the different density models, are practically coincident. Therefore, the DPL, BPL, CPL and EIN models do not offer any significant improvement in the description of the RRLs distribution with respect to the simpler SPL. Fig. 12 shows the comparison between the (elliptical) radial profiles of the density slope for different halo models. The slopes have been calculated as the logarithmic (elliptical) radial derivative of the logarithmic density, therefore the result for a SPL is just a constant (as shown by the black-dashed line in the four panels). The best-fit broken radius for the BPL is kpc: in the inner part the best-fit slope is practically the same as the SPL model then it decreases to a slightly larger slope of about 2.9.

Concerning the DPL, the posterior distributions of the slopes are compatible with , such that the final density profile is again an SPL with a slope compatible with the best-fit SPL model (2th-row panel in Fig. 12), moreover the posterior distribution of is almost uniform in the prior range interval (Tab. 2). Similarly, the core radius of the CPL model is very small and the final CPL model is compatible with no-core and the outer slope is the same of the SPL model (3th-row panel in Fig. 12).

Finally, for the EIN the fit favours very large values both for n () and for ( kpc), so that in the analysed radial range the variation of the slope is minimal and the density profiles effectively mimic an SPL. However, since the BIC of the EIN model is higher than the BIC of the SPL models, there is no evidence in favour of an EIN profile for the inner part of the stellar halo.

The break radius of the BPL model is compatible with the one obtained in Xue et al. (2015) after the subtraction of halo substructures, although the change of slope is much more significant in their work with respect to our result (see Tab. 4). In principle, the over-density of stars described in Sec. 3.3.4 might make it difficult to detect the break, so we repeated the maximum likelihood analysis using the RRLs only above and below the Galactic plane. The results are shown in Fig. 12 by the solid () and dashed () coloured curve: in general the stars above the Galactic plane prefer a mean slope of about 2.65 while the ones below have a mean slope of about 2.78. Concerning the change of slope above the Galactic plane, the posterior distribution of the parameters of the BPL, DPL and CPL models are such that these models are similar to the SPL profile. Below the Galactic plane both the BPL and the DPL favour a change of slope, however the one of the BPL model is much more significant going from about 2.73 in the inner part to about 3.3 beyond kpc. It is evident that the stars below the Galactic plane have a steeper decrease of the density with respect to the stars above the Galactic plane. This result could be due to the excess of stars at high latitude, however in both the sub-samples the minimum BIC is still obtained for the SPL model. Therefore, the change of slope is not highly significant and the fact that the DPL and the BPL show a different behaviour could indicate that the change of slope could be due to some local artefact (such as a local decrease of the completeness).

The comparisons discussed above refer only to DN halo models, however we obtained similar results also for DP and TR halo models. In particular, we found that the SPL density model always provides the lowest BIC value, independent of the assumed halo geometry.

In conclusion, there is no significant evidence of deviation from a SPL density law: the RRLs follows an SPL with an exponent that ranges between 2.6, assuming a SH model, to 3 assuming a halo. These results agree with the estimate of the density profiles obtained in Sec. 3.3.2 counting stars in ellipsoidal shells (Fig. 8).

4.4.3 Iso-density surfaces

As shown by the BIC comparison in Tab. 3, the SH and DP models give significantly worse fits with respect to the DN and TR models. This result is confirmed by a comparison of the properties of the RRLs in our sample with the ones expected for the best-fit models. For example, Fig. 13 shows the fraction of RRLs in our sample in different observed magnitude bins (black points) compared to the predicted distribution for different halo models (assuming a SPL, see Tab. 3). The SH and the SP models predict too low a fraction of stars at lower and a significant excess around , while the DN and the TR models give a good match to the data.

Figure 13: As in Fig. 12 but for halo models SH (blue, and ), DN (orange, and ), DP (green, and ) and TR (red, and ). In all cases we assumed a SPL for the density profile of the RRLs (with parameters given in Tab. 3).

As a further refinement we tested the and models. Allowing q to vary, we obtained a better description of the distribution of RRLs as shown by the decreases of the BICs in Tab. 3. The variation of the flattening is quite significant as it is shown in Fig. 14. The halo is largely flattened () in the very inner part and then it becomes more spherical reaching a flattening of about 0.75 at the border of the Galactic volume analysed in this work ( kpc). These results are in agreement with what we have seen directly in the distribution of stars in the - plane in Sec. 3.3.3 and Fig. 9. Essentially the same result is obtained by considering the stars only above or only below the Galactic plane. Fig. 15 shows the comparison between the magnitude distribution of the stars in our sample and the one expected for the varying q model.

In the (non-tilted) TR and the DP models the elongated axis makes an (anticlockwise) angle with respect to the positive part of the Galactic Y-axis. It is interesting to compare the orientation of the halo principal axes to the one of the Galactic bar. The orientation of the Galactic bar 444 is defined as the anticlockwise angle between the direction of the bar and the Galactic X-axis. is uncertain (Antoja et al., 2014), depending on the works it ranges from about (e.g. Robin et al. 2012) to about (e.g. Benjamin et al. 2005). Given these large uncertainties, the angle can be considered compatible with . As a consequence, the elongated axis (Y-axis) of the halo is almost perpendicular to the orientation of the Galactic bar. It is unclear if this correlation is a coincidence, due to the large errors, or if it reflects a real link between the two structures. For the purpose of this work, we note that the fact that the elongated axis is almost perpendicular to the Galactic bar rules out the hypothesis that the fit has been influenced by some contaminants coming from the inner Galactic bar.

Figure 14: Halo flattening, q, as a function of the elliptical radius (Eq. 9). The black solid-line shows the median of the functional-form distribution of q obtained for our best halo model (SPL-TR; last row in Tab. 3) while the dark and light green bands indicate the 68% and 95% confidence intervals. The dashed line indicates the best-fit q obtained for the SPL-TR halo model (8th row in Tab. 3). The orange curve shows the best-fit functional form of q found in Xue et al. (2015) using a sample of K giants. The other curves show the non-parametric estimate of q from Das & Binney (2016) (magenta line) using a sample of K giants stars and from Das et al. (2016) (red line) using a sample of BHB stars. Concerning the best-fit parametric functional forms of Das et al. (2016), we note that they found a profile that is very similar to the one of Xue et al. (2015)(practically coincident, except at the very inner radii). In the work of Xue et al. (2015), q is a function of the spherical radius rather than the elliptical radius . In order to compare their and our results the orange line has been calculated along the Galactic vertical direction where , (Eq. 5) and (Eq. 9). In this case our estimate of the flattening q at should be compared with their estimates calculated at .

4.4.4 Tilt and offset

As a final analysis, we considered halo models that can be tilted with respect to the Galactic plane and offset from the Galactic centre (see Sec. 4.2.2). We added these refinements to the SPL-TR model.

For the offset model, the best-fit parameters are

  • kpc,

  • kpc,

  • kpc,

and the size of the displacement vector with respect to the Galactic centre is . The offset we found is small and quite insignificant compared to the radial extent of the halo. Indeed these small values could also be due to the over-fitting of some local substructure at small elliptical radii. Moreover, both the posterior distribution of and are compatible (within and , respectively) with 0 and most of the displacement is along the Galactic X-axis. This axis connects the Sun to the Galactic centre, therefore the offset in this direction can be also interpreted as changing the distance of the Sun from the Galactic centre. Considering our results, the estimated Sun position is kpc, which is totally compatible with the Sun-distance estimates in literature (e.g. McMillan & Binney 2010).

For the tilt model, we take the triaxial model and allow the axes of symmetry to be arbitrarily oriented. The best-fit results for the (anticlockwise) angles describing the orientation of the halo (see Sec. 4.2.2) are

  • (rotation around the Z-axis),

  • (rotation around the new Y-axis),

  • (rotation around the final X-axis).

The angle represents the tilt of the elongated axis of the halo with respect to the Galactic Y-axis and it is compatible with the orientation found for the non-tilted DN and TR disc-plane models (Tab. 3). The angles and measure the tilt of the halo meridional plane with respect to the plane of the Galactic disc. Even though the BIC indicates a mild improvement with respect to the non-tilted triaxial model (), the tilt is quite small.

Finally, we compared the ability of tilted and non-tilted models to reproduce the observed distribution of stars. Fig. 15 shows the comparison between data and different models for the magnitude distribution of the stars. The models are , , and . Note that all of them are capable to give a good match to the data.

In conclusion, we conservatively assume that there is no need for a significant tilt to describe the RRL distribution in the inner halo. A more detailed study about this interesting properties is postponed to future works in which we will extend the coverage of the Galactic volume (see Sec. 5.5).

Figure 15: As in Fig. 11 but for triaxial halo models with varying q (blue), halo tilt (green), centre offset (yellow) and a combination of the previous models (red), see Sec. 4.4.4 for further details.

5 Discussion

5.1 Best halo model

Figure 16: Same as in Fig. 9 but for mock stellar halos. Top: density maps for the best-fit SPL-DN model. Bottom: density maps for the best-fit SPL-DN (see Tab. 3). The maps have been obtained by averaging over 1000 halo mock realisations. Note that for a halo with a fixed flattening, for the (almost) correct value of q (middle panel) the iso-density contours remain vertical across the whole range of elliptical radii. This is at odds with the RRL distribution in the Milky Way. As the middle panel of Fig. 9 demonstrates, the contours indeed start vertical for kpc, but beyond that, there is a noticeable bending away from the centre. Our best-fit model with varying q displays exactly the same behaviour (see middle panel of bottom row).
Figure 17: Top panels: all-sky star number density maps in Galactic coordinates. Bottom panels: sky-maps density plots in equatorial coordinates. The left-hand panels show the data, while the middle panels show our best model (SPL-TR, see Sec. 5.1 and Tab. 3), the right-hand panels show the relative data-model residuals. The model maps have been obtained integrating the star density distribution (Eq. 21) through the magnitude interval assuming a constant absolute magnitude