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

Abstract

1 INTRODUCTION

The diffuse cloud of stars observed around the Milky Way (MW) 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, White & Guiderdoni 1993; Sommer-Larsen, Götz & Portinari 2003; Chan et al. 2015). The two haloes emerge alongside each other, sharing the formation mechanism, i.e. a combination of the accretion on to 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, Belokurov & Evans 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 (MSTO) stars (e.g. Sesar, Jurić & Ivezić 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 data sets, 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 −2.5, to a much steeper one, consistent with a power-law index of ≈−4.

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). None the less, 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, Mao & Wechsler 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 sightlines 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 (SDSS) 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 can 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 blue horizontal branch (BHB) stars with spectra, Das, Williams & Binney (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, Evans & Williams 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, Wirtanen & Janes (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 MW and have been used to describe the Galactic halo with unwavering success over the last 50 yr (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 the 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 ≈21 600 RRLs out to ≈20 kpc from the Sun, with constant completeness of ≈ 20 per cent and purity ≈ 90 per cent.

Armed with this unprecedented data set, 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 ≈30 kpc. Finally, we also allow the halo to be (i) arbitrary oriented, (ii) triaxial and (iii) offset 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 N-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-fitting 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 G 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 values 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, Belokurov & Torrealba 2017). The relevant GaiaSource quantities used here (other than the sky coordinates RA, Dec.), are:

• N_obs, the number of times a source has crossed a CCD in the Gaia's focal plane.

• F_G, the flux (electron per second) measured in the G band averaging over N_obs single flux measurements.

• |$\sigma _{{\rm F}_{G}}$|⁠, the standard deviation of the N_obs flux measurements.

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

• 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 |$=\sqrt{\mu ^2_\alpha \cos ^2{\delta } + \mu ^2_\delta }$|⁠, where μ_α and μ_δ are the proper motions measured along RA and Dec., respectively.

2.2 RR Lyrae in Gaia DR1

The selection of variables through the AMP parameter suffers from the limitations that AMP is a time-averaged information and does not allow us 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 Section 2.2.6.

2.2.1 The Gaia + 2MASS sample

As mentioned, GDR1 reports photometric information only in the G band. We derived a colour (J − G) 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, Finkbeiner & Davis (1998) and the transformation A_G = 2.55E(B − V) for the G band (see Belokurov et al. 2017) and A_J = 0.86E(B − V) for the J band (Fitzpatrick 1999).

2.2.2 Auxiliary RR Lyrae data sets

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, 2014) and the Stripe82 (S82, Sesar et al. 2010) catalogues. These samples allowed us to identify the optimal selection criteria, analyse the completeness and the contamination of the catalogue¹ and estimate the RRL absolute magnitude in the G band. The CSS contains about 22 700 type-ab RRL stars distributed over a large area of the sky (about 33 000 deg² between 0° < α < 360° and −75° < δ < 65°) and extended up to a distance of 70 kpc. The completeness of this sample is constant (at ∼ 65 per cent) only for 13 < V < 15, while it quickly decreases outside this range. Most importantly, as shown in fig. 13 of Drake et al. (2013), for objects fainter than V ∼ 15, the completeness is a strong function of the number of observations and thus varies appreciably across the sky. SDSS's Stripe82 covers a 2.5°-wide and 100°-long patch of sky aligned with the celestial equator and contains ‘only’ 483 RRLs. However, the sample is very pure (with less than < 1 per cent 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 Section 2.2.3) and to estimate the absolute magnitude in the G band (see Section 2.2.4), while the high quality of S82 sample is ideal to analyse the completeness and contamination of our final sample (see Sections 2.2.5 and 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 Section 2.2.2) in order to maximize completeness of the sample and its spatial uniformity, while keeping the level of contamination low (see Section 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 b = −10° and 10°. Our limit in latitude (|b| = 10°) 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 section 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.

Fig. 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 N_obs–AMP (left-hand panel), colour–magnitude (middle panel) and AEN–G (right-hand panel) planes. The bona-fide RRLs occupy a well-defined strip in colour, thus we excluded all the stars with the J − G colour index greater than −0.4 and lower than −0.95 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 < −0.7 (see Section 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, G = 17.1 is the faintest magnitude that we can reach to obtain a sample with spatially uniform completeness (see Section 2.2.5 for further details). The number of bright stars with G < 10, 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 Sections 2.2.5 and 2.2.6), we decided to put the bright magnitude limit at G = 10.

Figure 1.

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

Open in new tabDownload slide

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 AEN > 0.65 as shown by the horizontal black line. This cut likely removes contaminant extragalactic objects since they typically have AEN ≈ 2 (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 PM greater than 50 mas yr⁻¹. Given differences in the light curves and its sampling, the significance of AMP (equation 1) might depend on the number of photometric measurements N_obs. With this in mind, we impose N_obs > 30: the focal plane of Gaia has an array of 9 × 7 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 data set (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 (l, b) = (280.47°, −32.89°) with an angular radius of 9° for the Large Magellanic Cloud (LMC) and the other centred on (l, b) = (302.80°, −44.30°) with a radius of 7° 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 Section 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: l = [167°, 189°] b = [16°, 22°] for S1 and l = [160°, 190°] b = [63°, 73°] 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 30 arcmin, 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 Section 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 Sections 3 and 5.

The final sample contains about 21 600 RRLs that can be used to have a direct look at the distribution of stars in the Galactic halo (Section 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 per cent 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 per cent in Deason et al. 2011). The summary of the applied selection cuts can be found in Table 1.

2.2.4 Absolute magnitude and distance estimate

Figure 2.

Distribution of G-band absolute magnitudes for bona-fide RRLs in the GCSS catalogue (Gaia+CSS, see Section 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).

Open in new tabDownload slide

2.2.5 Completeness

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 SDSS (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 16 < r < 18, where r is the r-band magnitude. In this magnitude range, the SDR7 can be considered 100 per cent complete.² We chose the r band because the peak of the filter response is almost coincident in wavelength with the one of the Gaia G band (see Section 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 b > 10°.

Figure 3.

Completeness analysis of the sources in GaiaSource (Section 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 (N_Gaia) and stellar sources in the SDSS DR7 (N_SDSS, 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 16 < G < 18 cut, and the SDSS sources using the 16 < r < 18 cut (further details on the text).

Open in new tabDownload slide

We estimated the completeness (and the contamination) using the RRLs in our auxiliary catalogues: CSS and S82 (Section 2.2.2). In particular, S82 represents a complete (∼100  per cent) and pure (∼99  per cent) 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 Section 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 (−0.65 red triangles, −0.70 blue circles and −0.75 green diamonds) in the range of magnitudes 15 < G < 17. 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 per cent for AMP > −0.65 to about 30 per cent for AMP > −0.75, but it is reasonably constant up to about G = 17.5 then it abruptly decreases to 0 at G ≈ 18, so we decided to conservatively cut our sample at G = 17.1 (vertical black line in Fig. 4, see Section 2.2.3).

Figure 4.

Completeness of our samples of RRLs as a function of distance from the Sun (D_⊙) or G magnitude. The conversion from D_⊙ to G has been obtained from equation (2) assuming a constant absolute magnitude M_RRL = 0.525. Different symbols indicate completeness for samples obtained using different AMP cuts: red triangles for AMP > −0.65, blue circles for AMP > −0.70 and green diamonds for AMP > −0.75. 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 G magnitude range of 15 < G < 17, while the vertical black line marks the G = 17.1 faint magnitude limit of our final sample (see Section 2.2.5).

Open in new tabDownload slide

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 N_obs. 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 (equation 1) improves its effectiveness in selecting RRLs. However, the increase is not dramatic as all the differences are within 10 per cent. 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 Table 1).

2.2.6 Contamination

Figure 5.

Contamination of the RRL sample as a function of the AMP cut. The contamination is estimated using the S82 catalogue (see Section 2.2.2) in the range of G magnitudes 15 < G < 17. The error bars have been calculated by propagating Poisson uncertainties from the number of stars in AMP bins.

Open in new tabDownload slide

3 DISTRIBUTION OF THE RRLs IN THE INNER STELLAR HALO

3.1 A first Gaia look at the stellar halo

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 G magnitude intervals: 10 < G < 15.5, corresponding to heliocentric distances between ∼1 and 10 kpc (equation 2), and 15.5 < G < 17.1, i.e. with heliocentric distances between 10 and 20 kpc (assuming that all RRLs have absolute magnitude M_G = 0.525). We stress that Fig. 6 represents the first all-sky view of the distribution of the RRLs in the inner halo.

Figure 6.

Star-count all-sky maps in the Galactic coordinates (l, b) for the RRLs in our sample (Section 2.2.3) in the magnitude intervals 10 < G < 15.5 (upper panel) and 15.5 < G < 17.1 (lower panel). The ‘holes’ between l = −50° and −100° are due to the mask used to eliminate the contribution of the LMC and SMC (see Section 2.2.3).

Open in new tabDownload slide

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 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 (l ≈ 0) and there are not evident asymmetries with respect to the Galactic plane. The second magnitude range covers the Galactic distances between 18 and 28 kpc in the side of the Galaxy containing the Sun and between 3 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 l = −50°. In both magnitude intervals, an overdense band of stars at low Galactic latitude (|b| < 15°–20°) can be seen running all around the Galaxy. A discussion of the nature of these structures can be found in Section 5.3.

3.2 Setting the frame of reference

3.3 Density distribution of the halo RRLs

3.3.1 Meridional plane

Fig. 7 shows the number density of the RRLs of our sample in the Galactic meridional plane R–Z_g. 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 Z_g. 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 Section 5.3. The density at higher Z_g 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 q = 0.6 elliptical contours (overplotted in Fig. 7) out to R ≈ 15-20 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 overdense region at Z_g > 10 kpc, which does not seem to have a counterpart below the Galactic plane (see Section 3.3.4 for further details).

Figure 7.

Number density in the Galactic R-Z_g 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 (equation 10). The black contours are plotted with spacing Δlog ρ = 0.4 from log ρ = −1.2 to 1.6, where ρ is in units of kpc⁻³; the white dashed lines represent elliptical contours with an axial ratio q = 0.6.

Open in new tabDownload slide

3.3.2 Density 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 (p = 1 and q ≠ 1) shells and dividing this number by the shell volume corrected by the coverage of our sample (equation 10). The profiles are shown in Fig. 8 for q = 1.0, 0.6 and 0.3. Independently of the assumed value of q, the density of RRLs follows a single power law (SPL) 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 Section 5.4 for a detailed comparison with results obtained in the literature).

Figure 8.

1D number density profiles of the halo RRLs as functions of the elliptical radius r_e (equation 9) assuming that the halo is stratified on spheroids with p = 1 and different values of q : q = 1.0 (blue diamonds), q = 0.6 (green dots) and q = 0.3 (red squares). The black dashed line shows, for comparison, an SPL with index α = −2.7. The Poissonian errors are smaller than the size of the symbols. The density normalization is arbitrary and different for each profile.

Open in new tabDownload slide

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

3.3.3 Halo flattening

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 q ≈ 0.6. However, in order to more rigorously study the halo flattening, we estimate the density in the r_e-θ plane (see equations 4 and 9). In practice, for a given value of q, we scan the density as a function of r_e at fixed θ. If the RRLs are truly stratified on similar spheroids with axial ratio q, the density is independent of θ, so the iso-density contours in the r_e–θ plane are vertical stripes. If the assumed value of q is smaller than the true value q_true, r_e is underestimated, the estimated density is a monotonic decreasing function of θ and the iso-density contours in the r_e–θ plane are bent in the direction of |θ| increasing with r_e (provided that the density is a decreasing function of r_e). The isodenisty contours are bent in the opposite direction, if one assumes q > q_true. The shape of the iso-density contours in the r_e–θ 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 r_e–θ plane of the RRLs in our sample are shown in Fig. 9 for q = 0.75, 0.55 and 0.35, assuming p = 1. Below |θ| = 20° (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 overestimated in all the panels. At higher latitudes (|θ| > 20°), the contours give a direct indication on the flattening of the halo: in the right-hand panel of Fig. 9 (q = 0.35) the iso-density contours are significantly inclined in a way that implies q_true > 0.35. For r_e < 20 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 (q = 0.55), but beyond 20 kpc, the contours start to bend so that q_true > 0.55. The last two iso-density contours in the left-hand panel (q ≈ 0.75) 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; Das et al. 2016). In Section 4.4, we perform a comprehensive model fitting analysis of the RRL data set and compare our results to those found in the literature.

Figure 9.

Number density of RRLs in the elliptical radius (r_e)–Galactocentric latitude (θ) space for q = 0.75 (left-hand panel), q = 0.55 (middle panel) and q = 0.35 (right-hand panel), assuming p = 1. The density is normalized to the maximum value of each panel. The contours show the normalized density levels (0.0002, 0.00055, 0.001, 0.002, 0.004, 0.02, 0.04, 0.2), while the dashed lines indicate |θ| = 20°.

Open in new tabDownload slide

3.3.4 Vertical asymmetries

Contours of the RRL density shown in Figs 7 and 9 display an overall symmetry between the Northern and Southern Galactic hemispheres, however at Z_g > 10 kpc and around θ ≈ 60° an overdensity 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 Z_g slabs to understand whether the overdensity 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 X_g-Y_g plane integrated along three different Z_g 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 per cent) 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 Z_g slabs (middle panels), the difference between the number of stars in the region above and below the Galactic plane is small (less than 5 per cent) 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 Z_g (right-hand panels) show a significant difference in the total number of stars (about 20 per cent) that cannot be explained by Poisson fluctuations only. Moreover, the excess of stars above the Galactic plane is strongly clustered in the regions between X_g ≈ (5, 15) kpc and Y_g ≈ (0, 10) kpc.

Figure 10.

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

Open in new tabDownload slide

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 Section 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 Section 4.3.5) of different halo models and find that the Magellanic Clouds mask introduces a difference of about the 7 per cent in the number of objects above and below the Galactic plane for |Z_g| > 10 kpc, significantly less than the value of 20 per cent obtained here. Therefore, the structure seen at high positive Z_g appears to be genuine, and is most likely related to the ‘Virgo overdensity’ (Jurić et al. 2008; Vivas et al. 2016). Across all Z_g > 0 slabs, portions of the Virgo Cloud are visible as an excess of stars at positive X_g. 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 overdensity at negative X_g (and Z_g).

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

Figs 7 and 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 per cent 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 (Section 2.2.3). Fig. 9 suggests that the substructure can be effectively eliminated with a selection in angle θ (see equation 7).

In particular, at |θ| = 20°, 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 |θ| > 20°: this subsample contains 13 713 objects and covers approximately 44 per cent of the halo volume within a sphere of radius 28 kpc (see Table 1). We stress that this cut is based on values of θ obtained from equations (2) and (7) assuming the same value of the absolute magnitude, M = M_RRL = 0.525, for all the stars (see Section 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 Z_g. We found that the results obtained for the sample with |θ| > 20° (see Section 4.4) are qualitatively similar to the results obtained using a sample with |b| > 30° (the lowest Galactic latitude used in most of the previous works, e.g. Deason et al. 2011) or a sample with |Z_g| > 6 kpc (Fermani & Schönrich 2013 uses |Z_g| > 4 kpc to cut disc stars). The cuts on b and Z_g 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 |θ| = 20°.

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 Section 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

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 α_inn < −3 or with α_out > −3 imply haloes 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 an SPL with α_inn < −3 (see Section 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 p = 1 and q = 1 in equation (9), so that r_e is just the spherical radius D_g (see definition 5) in the Galactic frame of reference;

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

• disc–plane axisymmetric (DP): we set q = 1, 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

  \begin{equation} {\rm {{\it q}}}({\rm {{\it r}}}_{\rm e})={\rm {{\it q}}}_\infty - ({\rm {{\it q}}}_\infty -{\rm {{\it q}}}_0) {\rm exp}\left[ 1- \frac{\sqrt{{\rm {{\it r}}}_{\rm e}^2+{\rm {{\it r}}}_{\rm eq}^2}}{{\rm {{\it r}}}_{\rm eq}} \right], \end{equation}

  (15)

  so q varies from q₀ at the centre to the asymptotic value q_∞ at large radii and the variation is tuned by the exponential scalelength r_eq;

• (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 (equation 9), we transform the Galactic coordinates (see equation 4) of each star by applying a rotation matrix R(γ, β, η) 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;

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

4.2.3 Complete halo model

Each complete halo model is defined by a model for the density law (SPL, DPL, CPL, BPL and EIN) plus a model for the shape of the iso-density surfaces (SH, DN, DP and TR) and any combination of geometrical variants (qv, tl and off). For instance, an SPL-SH model has a spherical distribution of stars with an SPL density profile, while an EIN-TR^{qv, tl} model is a triaxial tilted model with varying flattening along the Z-axis of symmetry and an EIN density profile (see Table 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 and EIN), geometrical models (SH, DN, DP and TR) and model variants (qv, tl and off). See Table 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

4.3.2 Likelihood of a single star

4.3.3 The total likelihood

4.3.4 Sampling of the parameter space

Concerning the parameter r_eb, we assume, for both the BPL and DPL density models, that the prior distribution is uniform between |${\rm {{\it r}}}^{{\rm min}}_{\rm e}$| and |${\rm {{\it r}}}^{{\rm max}}_{\rm e}$| which represent the minimum and maximum elliptical radii (equation 9) of stars in our sample. For most of our models, |${\rm {{\it r}}}^{{\rm min}}_{\rm e}$| remains constant at about 2.4 kpc, while |${\rm {{\it r}}}^{{\rm max}}_{\rm e}$| ranges between roughly 31 and 38 kpc.

We analysed the posterior probability of the parameters, |$\vec{\mu }$|⁠, for a given halo model using the equation (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 emcee³ (Foreman-Mackey et al. 2013). The technical details of our approach can be found in Appendix B.

The final best-fitting 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 and magnitude limits) as our clean sample of RRLs (Section 4.1). This algorithm distributes the stars using a combination of the density laws (Section 4.2.1) and according to the assumptions on the properties of the iso-density ellipsoidal surfaces (Section 4.2.2). The absolute magnitude of the stars in the mock catalogues are distributed using the best-fitting 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 Section 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 p = 1 and q = 1, we can detect a break in the density profile within about r_eb = 28 kpc (equation 21). For more flattened models (q around 0.5/0.6), this range extends up to about r_eb = 35 kpc. However, in this latter case most of the information for regions where r_e > 28 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 Section 3.3.4 and Fig. 9. In conclusion, we are confident that we are able to robustly detect significant deviations from an 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 data set.

4.4 Results

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

4.4.1 Model comparison

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 |$\dim (\vec{\mu})$|⁠, such that for two models with the same likelihood, the one with more parameters is penalized. The best-fitting parameters together with a comparison of the maximum likelihood and the BIC for a representative sample of halo models can be found in Table 3. These quantities have been calculated taking the maximum of the 14 400 likelihood estimates obtained with the final MCMC sampling (see Section 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

We tested different density profiles assuming a DN geometrical model: the SPL (equation (11) with α_inn = α_out), the DPL (equation 11), the BPL (equation 12), the CPL (equation 11 with α_inn = 0) and the EIN profile (equation 13).

The logarithmic maximum likelihoods and the BICs obtained for these different models are very similar as shown in Table 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 an SPL is just a constant (as shown by the black dashed line in the four panels). The best-fitting broken radius for the BPL is r_eb ≈ 22 kpc: in the inner part the best-fitting slope is practically the same as the SPL model then it decreases to a slightly larger slope of about 2.9.

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 Table 3).

Open in new tabDownload slide

Figure 12.

Radial profile of the density slope (∂(ln ρ)/∂(ln r_e)) for different density law models: BPL (first row), DPL (second row), CPL (third row) and EIN (fourth 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 (b > 10°, dotted line) and below (b < −10°, solid line) the Galactic plane. The dashed black line shows the best-fitting slope obtained for the SPL model. In all the cases, we assumed a DN geometrical halo model (Table 3).

Open in new tabDownload slide

Concerning the DPL, the posterior distributions of the slopes are compatible with α_out = α_inn, such that the final density profile is again an SPL with a slope compatible with the best-fitting SPL model (second row panel in Fig. 12), moreover the posterior distribution of r_e is almost uniform in the prior range interval (Table 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 (third row panel in Fig. 12).

Finally, for the EIN the fit favours very large values both for n (≈38) and for r_eb (≈390 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 Table 4). In principle, the overdensity of stars described in Section 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 (b < −10°) and dashed (b > 10°) 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 r_e ≈ 20 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 subsamples 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 an SPL density law: the RRLs follows an SPL with an exponent that ranges between 2.6, assuming an SH model, to 3 assuming a TR^qv halo. These results agree with the estimate of the density profiles obtained in Section 3.3.2 counting stars in ellipsoidal shells (Fig. 8).

4.4.3 Iso-density surfaces

As shown by the BIC comparison in Table 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-fitting 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 an SPL, see Table 3). The SH and the SP models predict too low a fraction of stars at lower G and a significant excess around G = 15, 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, q = 1 and p = 1), DN (orange, q = 0.58 and p = 1), DP (green, q = 1 and p = 1.59) and TR (red, q = 0.65 and p = 1.27). In all cases, we assumed an SPL for the density profile of the RRLs (with parameters given in Table 3).

Open in new tabDownload slide

As a further refinement, we tested the DN^qv and TR^qv models. Allowing q to vary, we obtained a better description of the distribution of RRLs as shown by the decrease of the BICs in Table 3. The variation of the flattening is quite significant as it is shown in Fig. 14. The halo is largely flattened (q ≈ 0.57) 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 (r_e ≈ 30 kpc). These results are in agreement with what we have seen directly in the distribution of stars in the r_e–θ plane in Section 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.

Figure 14.

Halo flattening, q, as a function of the elliptical radius r_e (equation 9). The black solid line shows the median of the functional-form distribution of q obtained for our best halo model (SPL-TR^qv; last row in Table 3), while the dark and light green bands indicate the 68 per cent and 95 per cent confidence intervals. The dashed line indicates the best-fitting q obtained for the SPL-TR halo model (eighth row in Table 3). The orange curve shows the best-fitting 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-fitting 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 innermost radii). In the work of Xue et al. (2015), q is a function of the spherical radius D_g rather than the elliptical radius r_e. In order to compare their and our results, the orange line has been calculated along the Galactic vertical direction where X_g = Y_g = 0, D_g = Z_g (equation 5) and r_eq(r_e) = Z_g (equation 9). In this case, our estimate of the flattening q at r_e should be compared with their estimates calculated at D_g = qr_e.

Open in new tabDownload slide

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 Section 4.4.4 for further details.

Open in new tabDownload slide

In the (non-tilted) TR and the DP models, the elongated axis makes an (anticlockwise) angle γ ≈ −20° 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 Φ_bar⁴ is uncertain (Antoja et al. 2014), depending on the works it ranges from about −10° (e.g. Robin et al. 2012) to about −45° (e.g. Benjamin et al. 2005). Given these large uncertainties, the angle γ can be considered compatible with Φ_bar. 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.

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 Section 4.2.2). We added these refinements to the SPL-TR model.

For the offset model, the best-fitting parameters are

• X_off = 0.39 ± 0.05  kpc,

• Y_off = −0.17 ± 0.06  kpc,

• Z_off = −0.01 ± 0.01  kpc,

and the size of the displacement vector with respect to the Galactic Centre is D_off = 0.43 ± 0.07. 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 overfitting of some local substructure at small elliptical radii. Moreover, both the posterior distribution of Y_off and Z_off are compatible (within 3σ and 1σ, 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 X_⊙ = 7.61 ± 0.05 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-fitting results for the (anticlockwise) angles describing the orientation of the halo (see Section 4.2.2) are

• γ = −20.1° ± 2.7° (rotation around the Z-axis),

• β = 5.7° ± 0.8° (rotation around the new Y-axis),

• η = 3.13° ± 0.5° (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 (Table 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 (ΔBIC ≈ −40), 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 TR^qv, TR^tl, TR^off and TR^{qv, tl, off}. 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 Section 5.5).

5 DISCUSSION

5.1 Best halo model

We have explored the properties of the Galactic stellar halo using two independent methods. First, we analysed the density of the RR Lyrae candidates binned along the Z and R dimensions. This direct method was followed by an approach where the probability of observing each individual star in the space spanned by celestial coordinates and the apparent magnitude was described by a 3D model, whose parameters were optimized using MCMC likelihood sampling. Additionally, we augmented each method with a model comparison exercise. The results of the direct analysis (Section 3) and the maximum-likelihood fit (Section 4) are highly compatible. Namely, (i) the density of the RRL stars in the inner 25 kpc of the halo appears to be stratified on ellipsoids with a pronounced flattening along the Galactic Z-axis (Fig. 7), and (ii) the radial density profile can be described by an SPL (Fig. 8). While – as evidenced by Table 3 – this simple model gives an adequate description of the RRLs density distribution, it can be further refined by adding a mild elongation in the Galactic plane in combination with a radially varying flattening (along the Z-axis). Based on the BIC analysis (Section 4.4.1) as well as the direct comparison between the properties of the data and the best-fitting models, we chose to highlight this SPL–TR^qv (last row in Table 3) as the best-performing model amongst the family examined.

The distribution of the RRLs in the inner halo as described by the above best-fitting model has the following properties: the radial density is an SPL power law with an exponent α = 2.96 ± 0.05 and the iso-density surfaces of triaxial shape. The minor axis of the density ellipsoids is aligned with the Galactic Z-axis. We find strong evidence for an evolution of the vertical density flattening: the inner parts are squashed with (q ≈ 0.57), however the flattening becomes less pronounced (q ≈ 0.75) at the border of the Galactic volume analysed in this work (r_e ≈ 30 kpc). Note that the signs of the evolution of the halo flattening can be observed in the data directly as demonstrated in Fig. 9. In the plane of Galactocentric latitude and elliptical radius, the iso-density contours should appear vertical for the correct (and constant) value of q. Clearly, the contours are slanted in the left and the right-hand panels of the figure, where the halo flattening is over- and underestimated correspondingly. However, in the middle panel (corresponding to our best-fitting fixed q model), for a large range of radii, the iso-density contours do look vertical. Beyond 20 kpc, the contours start to bend away from the Galactic Centre, signalling a change in the halo shape. Fig. 16 illustrates the effect of the radial halo shape evolution with the θ-r_e views of the mock RRL realizations for two of our best-fitting (in each class) models. Note that the middle panel in the bottom row of the figure displays the same trend as the middle panel of Fig. 9, albeit with a lot less noise. With regards to the halo's major axis, in the Galactic plane, we measure a mild elongation (p = 1.27 ± 0.03) in the direction which is rotated with respect to the Galactic Y-axis by γ = −21.3° ± 2.6° (anticlockwise).

Figure 16.

Same as in Fig. 9, but for mock stellar haloes. Top: density maps for the best-fitting SPL-DN model. Bottom: density maps for the best-fitting SPL-DN^qv (see Table 3). The maps have been obtained by averaging over 1000 halo mock realizations. 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 MW. As the middle panel of Fig. 9 demonstrates, the contours indeed start vertical for r_e < 20 kpc, but at larger radii, there is a noticeable bending away from the centre. Our best-fitting model with varying q displays exactly the same behaviour (see middle panel of bottom row).

Open in new tabDownload slide

The measured RRL (as represented by the clean sample, see Section 4.1) distribution can be compared to the best-fitting model in Figs 17–19. Fig. 17 presents all-sky maps in both Galactic (top panels) and equatorial (bottom) coordinates. Fig. 18 gives the number density in the meridional plane, and, finally, Fig. 19 shows the number density in the r_e-θ space. In each figure, the rightmost column displays the map of residuals in the corresponding projection.

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^qv, see Section 5.1 and Table 3) and the right-hand panels show the relative data-model residuals. The model maps have been obtained integrating the star density distribution (equation 21) through the magnitude interval 10 < G < 17.1 assuming a constant absolute magnitude M = M_RRL = 0.525.

Open in new tabDownload slide

Figure 18.

Number density in the Galactic R-Z_g plane. Left-hand panel: density map for the RRLs in our clean sample (Section 4.1); middle panel: density map for our best model (SPL-TR^qv, see Section 5.1) and right-hand panel: relative data-model residuals. The black iso-density contours are plotted with interval of Δlog ρ = 0.4 from log ρ = −3 to 4, where ρ has the dimension kpc⁻³. The model map has been obtained averaging the maps of 1000 different mock catalogues made assuming the best-fitting parameters obtained for the SPL-TR^qv halo model (see Table 3).

Open in new tabDownload slide

Figure 19.

Same as in Fig. 18, but for the number density in the Galactic r_e-θ plane. The r_e in the plot have been obtained assuming q = 0.6 and p = 1 (see Section 3.3.3).

Open in new tabDownload slide

Reassuringly, the residual maps indicate that the model gives a reasonable description of the RRL distribution in the G2M sample. Overall, the MW's inner halo appears smooth, but some systematic model-data mismatches are visible. The most prominent (with residuals in some pixels corresponding to a density excess of order of 200 per cent) is the overdensity of stars visible in Fig. 17 at Galactic longitudes around −70° or RA of about 190°. This structure is also discernible in Fig. 18 for Z_g > 15 kpc and at about r_e = 35 kpc and above θ = 50° in Fig. 19. Most likely, this large stellar cloud is related to the Virgo overdensity, as discussed in Section 3.3.4.

Amongst other differences, a mild deficiency of stars is visible at the edge of the Galactic volume (available to us) in the anticentre direction, i.e. around l = 180°. This mild depletion is especially noticeable below the Galactic plane, however the amplitude of these outer Galaxy residuals is much smaller as compared to the mismatch caused by the Virgo overdensity. We argue that this deficiency of stars may be a hint of the steepening of the radial density profile, but, as discussed in Section 4, it does not obey the smooth parametric model considered here. A more robust study of the outer halo behaviour ought to be possible in the near future, when the coverage of the Galactic volume is extended well beyond the current radial limit (see Section 5.5).

5.2 Possible sources of bias in the maximum-likelihood analysis

5.2.1 Assumption on the absolute magnitude

Our halo fitting method (Section 4.4) relies on the assumption that all RRLs in the sample have the same absolute magnitude M_RRL = 0.525. However, as Fig. 2 demonstrates, this is not strictly true as indicated by the Gaussian-like wings in the M_G distribution. To understand the possible bias introduced by this simplification, we assumed that the pdf of the absolute magnitude P(M) in equa-tion (22) is represented by the best-fitting double Gaussian (blue line in Fig. 2). With this addition to the model, we repeated the fit of our best halo model (Section 4.3.4), this time marginalizing the likelihood over M. Perhaps unsurprisingly, the final posterior distributions of the model parameters are totally compatible with the ones obtained assuming a constant absolute magnitude.

5.2.2 Dust

Regions with high dust extinction add severe uncertainties in the study of distribution of stars in the Galaxy. In these regions, it is difficult to recover the faintest stars and the completeness could be much lower compared to the rest of the sky. Note however, that applying the cut on θ (Section 4.1), we eliminated most of the regions with high Galactic dust extinction: the 85 per cent of the stars in our ‘clean’ sample has E(B − V) < 0.25 and only the 1 per cent are located in regions with E(B − V) > 0.8. To be sure that our results are not influenced by the regions with high reddening, we repeated the fit of our best halo model (Section 5.1) masking the regions with E(B − V) > 0.4. We did not find any appreciable differences with respect to the unmasked analysis.

5.2.3 Contamination from the Galactic disc

The RRL sample was modelled (Section 4.3.4) using equa-tion (31) in the logarithmic likelihood in equation (25): we considered f as a free parameter, but fixed the disc's stellar density profile assuming values for R_d and Z_d. We tested a wide range of thin and thick disc models using R_d values between 2 and 4 kpc and Z_d values between 0.1 and 1 kpc. In general, we did not find significant differences with respect to the results shown in Table 3, except for the SH and DP models. In these two cases, the models try to reproduce the vertical flattening observed in the data with a high percentage of disc stars, but the final fits are poor as evidenced by the presence of large residuals. For all the other models considered, we found that the final posterior distribution for f  is lower than 3 per cent and is consistent with 0. Therefore, we conclude that our ‘clean’ sample (Section 4.1) does not harbour stars coming from the Galactic disc and the observed flattening along the vertical direction is a genuine property of the distribution of RRLs in the halo.

5.3 Halo substructures

The main focus of this work is the study of the overall distribution of RRLs in the inner halo. Reinforcing the principal assumption behind the modelling, the distribution of RRLs in our sample appears smooth. None the less, several substructures are present, of which the following two are the most obvious: (i) a highly flattened distribution of stars at low Galactic latitudes and (ii) an excess of stars at high Z_g. The first substructure is clearly evident in Fig. 7 and in Fig. 9 below θ = 20°. In our investigation of this substructure, we have attempted two approaches: first, we tried to exclude the regions most affected, i.e. those close to the disc from influencing the halo fit, and, second, we aimed to model the excess as a contamination coming from the Galactic disc. In particular, we repeated the fit using our original catalogue of about 22 600 RRLs (see Table 1) and the halo+disc stellar pdf defined in equation (31). The disc in the model was a realistic double exponential distribution: we fixed the radial (R_d) and vertical (Z_d) thin and thick disc scalelengths to the classical values found in literature (e.g. Das et al. 2016). Fig. 20 gives the density distribution in the R–Z_g plane for a disc+halo model in which the disc has R_d = 2.7 kpc and Z_d = 0.2 kpc. As illustrated by the residual map (right-hand panel), the halo+disc model is able to explain the flattened distribution of stars but only within about 10 kpc. As a consequence, the posterior distribution of the halo parameters favours a more flattened density (q ≈ 0.4), thus concealing the evidence of the radial variation of q, as measured and discussed above (Table 3) as part of the analysis of the ‘clean’ sample of RRLs (Section 4.1). The residuals suggest that the two-component model gives a poor fit with a significant mismatch both at low Z_g for R > 10 kpc and at high Z_g for R < 10 kpc. Note that these results do not change if different values for the disc scalelengths are assumed (R_d between 2 and 4 kpc and Z_d between 0.1 and 1 kpc). As a further test, we repeated the fit leaving the scalelengths of the disc model as free parameters. This test yielded unrealistically large values of R_d ≈ 25 kpc and Z_d ≈ 3 kpc. This can be compared to the results of the study of the thick disc's RR Lyrae by Mateu et al. (2011) who find the disc to be ‘antitruncated’, i.e. having a flatter density profile in the outer parts. In conclusion, it seems that this flat and elongated structure cannot be modelled simply by adding a double exponential disc with the properties expected for the MW. Instead, we conjecture that the observed excess could possibly be related to the Monoceros Ring (Newberg et al. 2002). Notice, however that it is still debated (see de Boer, Belokurov & Koposov 2017, and references therein) whether Monoceros represents a portion of the accreted halo (e.g. Sollima et al. 2011) or a structure related to the Galactic disc such as a warp and/or a flare in the outer disc (e.g. López-Corredoira & Molgó 2014).

Figure 20.

Left-hand panel: number density of stars in the Galactic R-Z_g plane for a multicomponent model with a halo and a disc; the black iso-density contours are plotted with interval of Δlog ρ = 0.4 from log ρ = −3 to 4, where ρ are in unit of kpc⁻³. Right-hand panel: relative residuals (data model/model) between the density in the R-Z_g plane for our complete sample of RRLs as shown in Fig. 7 and the model shown in the left-hand panel. The model map has been obtained averaging the maps of 1000 different mock catalogues made assuming the best-fitting parameters obtained for the disc + SPL-DN halo model (see Section 5.3). The stars in the disc follow a double exponential distribution with a radial scalelength R_d = 2.7 kpc and a vertical scaleheight Z_d = 0.2 kpc.

Open in new tabDownload slide

Another, a more prosaic, explanation could simply mean that the excess of stars at low latitude is due to Gaia artefacts such as the cross-match failures (discussed in Belokurov et al. 2017). Such conclusion is supported by the fact that no survey so far has reported a clear detection of RRL in the Monoceros Ring. Additionally, this structure is not observed in a sample of metal-poor K-Giants [Fe/H] < −1 from the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST) survey (Liu et al. 2017). Of course, Monoceros Ring is indeed a prominent feature of the stellar density maps in the direction of the Galactic anticentre. However, so far it has only been traced by typical disc stellar populations such as metal-rich MS stars (Yanny et al. 2003; Jurić et al. 2008; Morganson et al. 2016) or M giants (see e.g. Crane et al. 2003).

The second substructure is a significant overdensity of stars above the Galactic plane: it is visible in the density map in Fig. 7 for Z_g > 15 kpc and in Fig. 9 for θ > 50° and large elliptical radii; it is the cause of the mismatch between the star counts above and below the Galactic plane for |Z_g| > 10 kpc (right-hand panels in Fig. 10). Moreover, it is also traced by the significant increases of the residuals in the data-model comparison in Figs 17 and 18. The spatial location and shape of this structure indicate that it is most likely related to the Virgo overdensity (Jurić et al. 2008), a large and diffuse halo substructure detected with a variety of tracers, but most notably with horizontal branch stars, such as BHBs and RRL (see Duffau et al. 2006; Vivas & Zinn 2006; Deason et al. 2011; Vivas et al. 2016). As far as the quality of the halo models is concerned, by repeating the fit using the stars only below the Galactic plane, we find that the Virgo Cloud does not strongly influence our final results (see Section 4.4).

5.4 Comparison with other works

Thanks to the all-sky view, stable completeness and substantial purity, we have been able to test Galactic stellar halo models with a degree of complexity previously unexplored. Indeed, our best halo model (SPL-TR^qv, see Section 5.1) cannot be directly compared with results from other works. In fact, to date, only Deason et al. (2011) attempted to fit a triaxial model to the stellar halo. They found that the triaxial model does increase the quality of the model of the distribution of A-coloured stars, however they considered this to be an overfit, mainly due to the presence of substructures such as Monoceros and the Sagittarius stream. We have taken care to remove stars that may be related to Monoceros in our modelling (Sections 3 and 4.1). With regards to the Sagittarius stream, we have not found any worrying levels of contamination by its members in the volume covered in this analysis. Indeed, a recent work by Liu et al. (2017) shows that the residuals due to the Sagittarius stream are in regions not sampled in our work. Comparing the two triaxial models, the best-fitting p in Deason et al. (2011) is smaller than 1, while we infer p ≈ 1.3. Pérez-Villegas, Portail & Gerhard (2017) found that, due to the gravitational effect of the bar/bulge, a triaxial distribution of RRLs arises also from an initially oblate distribution of stars. This effect is restricted to a very small Galactic region (D_g < 5 kpc), however gravitational effects at large scale (e.g. the influence of the LMC and/or of the spiral arms) could also cause halo triaxiality outside the innermost region of the Galaxy. In this context, in the future, it will be interesting to study the trend of p as a function of the elliptical radius as already done in this and previous works for the vertical flattening q (see below).

Perhaps not totally unexpectedly, we have found strong evidence for the evolution of the shape of the stellar halo as a function of the position in the Galaxy. Both parametric (Section 4.4.3) and non-parametric (Section 3.3.3) analyses require the flattening of the halo along the Galactic Z-axis to decrease from the inner to the outer regions, i.e. q to increase at larger distances. Most recently, similar claims have been put forward in the works of Xue et al. (2015), Das et al. (2016), Das & Binney (2016) and Liu et al. (2017) (see Section 3.3.3 and Fig. 14). The results of Xue et al. (2015) have been obtained by employing the same functional form q(r_e) used in this work (equation 15), while Das & Binney (2016) used a non-parametric measure of q(r_e) estimated by the iso-density contours of their best-fitting dynamical model. Finally, Das et al. (2016) estimated the flattening using both the methods described above. Fig. 14 shows the results of these works compared to our best-fitting flattening profile: the non-parametric radial profile of Das et al. (2016) is compatible with our results, while the results found in Das & Binney (2016) point towards a distribution of stars that is systematically more spherical. While the best-fitting parametric functional forms of both Xue et al. (2015) and Das et al. (2016) are compatible with our measurement of the flattening trend for r_e > 20 kpc, there are several important differences. First of all, these studies postulated the dependence of q on the spherical radius. Here, to make our models self-consistent, we chose to express q as a function of the elliptical radius r_e (see equations 9 and 15). Second, the largest discrepancies are at small radii (r_e < 15 kpc): in our best-fitting model q is almost constant at ∼0.6, while in the other works the halo flattening keeps increasing (q decreasing) towards the centre. For example, around the Galactic Centre (q ≈ 0.4 in Das et al. 2016 and q ≈ 0.2 in Xue et al. 2015). Note however that the samples used in both these studies are actually depleted in stars within 10 kpc from the Galactic Centre, and thus the reported behaviour of q there is an extrapolation rather than a measurement. It is interesting to report that, in a recent work, Xu et al. (2018) analysed the distribution of a sample of metal-poor K Giants with a novel non-parametric technique and found that the halo becomes almost spherical (q > 0.8) already at small elliptical radii (r_e ≈ 20 kpc). Curiously, Deason et al. (2011) found no evidence for the variation of the halo flattening within 40 kpc from the Galactic Centre. The reason of these discrepancies is not clear, but we note that both the sky coverage and the halo tracer population is different compared to our study. In conclusion, we are still a long way from a general convergence about the details of the flattening of the stellar halo, even for the same halo tracers. In this respect, we stress again that the superior sky coverage of our study (and of Gaia in general) should assure very robust results, especially in the study of the halo flattening.

Our best SPL-DN model (see Table 3) is highly compatible with the previous results relying on RRLs (see Table 4) such as Miceli et al. (2008) and Vivas & Zinn (2006), but partially in contrast with the result of Sesar et al. (2013) who found a similar slope outside of 16 kpc but a flatter radial density profile (α ≈ −1) in the very inner part of the halo. However, the distribution of their RRLs in this region is very clumpy and cannot be properly fitted by a simple smooth power-law model. Similarly, Watkins et al. (2009) found that the overall distribution of RRLs cannot be described by a smooth model, but we note that their results are based on a drastically lower sky coverage with respect to this and other RRL-based studies.

Comparing to inference with tracers other than RRL, our best SPL-DN model is in good agreement with the work of Liu et al. (2017) based on red giant branch (RGB) stars and to Jurić et al. (2008) who used MS and MSTO stars. The studies employing other tracers such as BHB and K-Giants (Das et al. 2016; Xue et al. 2015; Deason et al. 2011) derive similar flattening of the halo, but usually report a steeper profile for the density law (α > −4) beyond 20–30 kpc. Note however that the work of De Propris, Harrison & Mares (2010) found a more shallow profile (α ≈ −2.5) up to 100 kpc. Some of the differences can simply be due to different stars used to trace the Galactic halo. In particular the RRLs tend to be more metal-rich as compared to the BHBs – an important difference given that there is evidence of a metallicity gradient in the halo (Carollo et al. 2007). While we are not sensitive to the changes in the radial density profile outside of 25 kpc, there have been claims that such ‘break’ is an artefact of not taking the shape evolution into account (see e.g. Liu et al. 2017; Xue et al. 2015). Further discussion can be found in Section 4.4.3). Table 4 summarizes the properties of the Galactic stellar halo as reported in the literature.

As a final note, it is important to mention that the distribution of our RRLs is nicely in agreement (in terms of both density profile and halo shape) with the one estimated by Pietrukowicz et al. (2015) for a sample of genuine RRL stars in the Galactic bulge. This result seems to strengthen the hypothesis that the RRL population in the bulge is consistent with being the inward extension of the Galactic stellar halo, as also found in N-body simulations by Pérez-Villegas et al. (2017).

5.5 Waiting for Gaia DR2

The Gaia collaboration announced the new data release (GDR2) for 2018 April.⁵ The GDR2 will have a series of important improvements with respect to the data we used in this work:

• the sky-coverage will be more uniform and the number of flux measurements (N_obs) for each star will increase, thus making our simple variability estimator AMP (equation 1) more precise and more robust;

• the photometry for a larger sample of variable stars will be used to fine-tune the selection cuts of RRLs without relying only on external auxiliary data sets (Section 2.2.2);

• most importantly, photometric colours G_BP and G_RP will be used for RRL colour selection without cross-matching Gaia with other surveys (as 2MASS in this work).

We envisage applying the method presented here to the GDR2, shedding the necessity to rely on shallow 2MASS data and thus extending our analysis out to Galactocentric distances of about 100 kpc. This unprecedented volume coverage, combined with stellar PMs will allow us to shed light on the substructures hidden in the Galactic halo and thus decipher the Galaxy's formation history (see e.g. Helmi et al. 2017).

6 SUMMARY

In this work, we presented the very first use of the Gaia DR1 photometric catalogue to study the properties of the Galactic stellar halo as traced by RR Lyrae. The key points of our study are listed below.

• We used the method proposed by Belokurov et al. (2017) to mine the Gaia DR1 for variable stars in advance of the official variable object release by the Collaboration. It is based on the estimate of the object's variability, AMP (equation 1) derived from the Gaia's mean flux and its associated error.

• Our principal selection cuts are those based on the AMP statistic and the colour index, J − G obtained from the cross-match between Gaia and 2MASS. It is the depth of the latter survey which governs the reach of our most distant RR Lyrae. Our sample extends over the whole sky and contains 21 643 stars within a sphere of 20 kpc centred on the Sun, covering an unprecedented fraction of the volume ( ∼ 58 per cent) of the inner halo (R < 28 kpc). While the overall completeness of our sample cannot compete with the levels attained by dedicated RRL surveys, it is stable across the sky and the magnitude range explored – the property most important for robust measurements of the stellar density behaviour in the halo. Additionally, we demonstrate that the sample's contamination is close to zero, and is not expected to exceed ∼10 per cent.

• Assuming a constant absolute magnitude for all RRLs in our sample (M_RRL = 0.525), we analysed their density distribution in the inner Galaxy. To that end, Fig. 6 presents the first all-sky stellar halo maps. These heliocentric slices reveal a relatively smooth nearby halo, with prominent interquadrant asymmetries starting to become noticeable beyond 10 kpc from the Sun.

• We point out two main substructures superimposed on the otherwise smooth distribution of the Galactic RRLs. One is visible as a highly flattened distribution of stars close to the disc plane, at odds with the more ellipsoidal distribution observed at high Z_g. It contains a little less than half of the stars in our sample, but given its large radial extent it cannot be explained as a contamination of the Galactic disc. It may be related to the Monoceros Ring (Newberg et al. 2002), but most likely it is caused by artefacts in the Gaia data at low Galactic latitude (see discussion in Belokurov et al. 2017). At high Galactic latitudes, a large overdensity of stars is evident, clustered in a region with Z_g ≈ 15 kpc and at Galactic radius between 5 kpc and 10 kpc, as illustrated in Fig. 10. We argue that this is likely a portion of the well-known Virgo overdensity (Jurić et al. 2008). Its counterpart at negative Z is the Hercules–Aquila Cloud, also discernible in the figure, on the other side of the Galaxy, beyond the bulge.

• The RRL density distribution is stratified on ellipsoidal surfaces flattened along the direction perpendicular to the Galactic disc. The density decreases regularly following a power law, without evidence of an inner core or a change of slope within the radial range accessible, as demonstrated in Fig. 8. For the first time, we are able to show the evolution of the halo flattening with radius directly in the data (see Fig. 9). The MW's stellar halo is rather flat in the inner parts but becomes more spherical at R ∼ 20 kpc. A similar property has been glimpsed in the sample of RGB halo stars considered by Liu et al. (2017).

• To mitigate possible biases associated with the low-latitude excess of the RRL candidates, we extracted a subset of 13 713 stars from the original sample. We applied to this ‘clean’ sample a maximum-likelihood fitting method, testing a large variety of halo models. The parameters of our final best model are in good agreement with the results obtained by the non-parametric analysis of the RRLs distribution. The halo is mildly triaxial with its major and intermediate axes in the Galactic plane (p = 1.27) and a significant flattening along the Galactic Z-axis which varies from q ≈ 0.57 in the centre to q ≈ 0.75 at the edge of the radial range analysed. The halo's major axis is rotated in the anticlockwise direction by a angle γ ≈ −21° with respect to the direction of the Galactic Y-axis. The density slope is well approximated by an SPL with exponent α = −2.96.

   We have demonstrated the power of the Gaia data for the Galactic stellar halo exploration. The wealth of the information available in GDR1 is unexpected and signals that a paradigm shift in the halo studies is imminent. One only has to wait for the GDR2 to be unleashed.

Acknowledgements

The authors wish to thank the members of the Cambridge Streams club for stimulating discussions and the anonymous referee for useful comments that improved this manuscript.

This work has made use of data from the European Space Agency mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement number 308024. VB, DE and SK acknowledge financial support from the ERC. SK also acknowledges the support from the STFC (grant ST/N004493/1).

Footnotes

REFERENCES