The first allsky view of the Milky Way stellar halo with Gaia+2MASS RR Lyrae
Abstract
We exploit the first Gaia data release to study the properties of the Galactic stellar halo as traced by RR Lyrae. We demonstrate that it is possible to select a pure sample of RR Lyrae using only photometric information available in the Gaia+2MASS catalogue. The final sample contains about 21600 RR Lyrae covering an unprecedented fraction () of the volume of the Galactic inner halo ( kpc). We study the morphology of the stellar halo by analysing the RR Lyrae distribution with parametric and nonparametric techniques. Taking advantage of the uniform allsky coverage, we test halo models more sophisticated than usually considered in the literature, such as those with varying flattening, tilt and/or offset of the halo with respect to the Galactic disc. A consistent picture emerges: the inner halo is well reproduced by a smooth distribution of stars settled on triaxial density ellipsoids. The shortest axis is perpendicular to the Milky Way’s disc, while the longest axis forms an angle of with the axis connecting the Sun and the Galactic centre. The elongation along the major axis is mild (), and the vertical flattening is shown to evolve from a squashed state with in the centre to a more spherical at the outer edge of our dataset. Within the radial range probed, the density profile of the stellar halo is well approximated by a single powerlaw with exponent . We do not find evidence of tilt or offset of the halo with respect to the Galaxy’s disc.
keywords:
galaxies: individual (Milky Way) – Galaxy: structure – Galaxy: stellar content – Galaxy: halo – stars: varaibles (RR Lyrae)1 Introduction
The diffuse cloud of stars observed around the Milky Way and known as the stellar halo is the alter ego of the much more massive structure, whose presence is inferred indirectly: the dark matter (DM) halo. The DM halo dominates the Galactic mass budget and, according to the currently favoured theories of structure formation, holds clues to a number of fundamental questions in astrophysics. These include, amongst others, the properties of the DM particles (see e.g Davé et al., 2001; Governato et al., 2004; Lovell et al., 2014), the nature of gravity itself (see e.g. Milgrom, 1983; Cabré et al., 2012), as well as the coupling between the DM and the baryons (e.g. Kauffmann et al., 1993; SommerLarsen et al., 2003; Chan et al., 2015). The two halos emerge alongside each other, sharing the formation mechanism, i.e. a combination of the accretion onto the Galaxy and the subsequent relaxation and phase mixing. Thus, there ought to exist a bond between them, which can be exploited to reveal the properties of the dark halo through the study of the luminous one. For example, by positing the continuity of the phasespace flow, the DM halo can be mapped out if the stellar halo spatial shape is known and complemented by stellar kinematics (see e.g. Jeans, 1915; Helmi, 2008; Posti et al., 2015; Williams & Evans, 2015a).
Leaving the ideas of James Jeans aside, could the structural parameters of the stellar halo alone inform our understanding of the mass assembly of the Galaxy? At our disposal are the numbers pertaining to the slope of the stellar halo’s radial density profile and its vertical flattening (see e.g. Jurić et al., 2008; Bell et al., 2008; Deason et al., 2011; Xue et al., 2015). The radial profile has so far been measured with a variety of stellar tracers. Studies based on Main Sequence TurnOff stars (e.g. Sesar et al., 2011; PilaDíez et al., 2015), Blue Straggler and Horizontal Branch stars (e.g. Deason et al., 2011) and RR Lyrae (e.g. Sesar et al., 2007; Watkins et al., 2009) seem to favour a “broken” profile. According to these datasets, somewhere between 20 and 30 kpc from the MW centre, the density slope changes from a relatively shallow one, as described by powerlaw index of approximately , to a much steeper one, consistent with a powerlaw index of .
In an attempt to interpret the observed radial density profile, Deason et al. (2013) conjecture that the presence or absence of a break is linked to the details of the stellar halo accretion history. In their exposition, a prominent break can arise if the stellar halo is dominated by the debris from a i) single, ii) early and iii) massive accretion event. This hypothesis appears to be supported by semianalytic MW stellar halo models (see e.g. Bullock & Johnston, 2005; Amorisco, 2017b) but is yet to be fully tested with Cosmological zoomin simulations (see, however Pillepich et al., 2014). Nonetheless, a consistent picture is now emerging: in line with other pieces of evidence, the “broken” MW stellar halo appears to be the telltale sign of an earlypeaked and subsequently quiescent accretion history. Note also that such destitute state of the stellar halo is not permanent but rather transient (see e.g. Deason et al., 2016; Amorisco, 2017a). In agreement with simulations, the MW stellar halo is destined to transform dramatically with the dissolution of the debris from the Sgr dwarf and the Magellanic Clouds.
While the radial density profile can be gauged based on the data from a limited number of 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 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 Acoloured 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 spectroscopicallyconfirmed K giants to detect a noticeable change of flattening with radius. Furthermore, they argue that if the halo shape is allowed to vary with radius, then a break in the radial density profile is not required. Finally, to add to the puzzle, based on a set of BHB stars with spectra, Das et al. (2016) report both evolving halo shape and the break in the radial density. Pinning down the shape of the stellar halo is important both for the dynamical inference of the shape of the DM halo (see e.g. Williams & Evans, 2015b; Bowden et al., 2016) and for our understanding of the response of the DM distribution to the presence of baryons (see Kazantzidis et al., 2004; Gnedin et al., 2004; Duffy et al., 2010; Abadi et al., 2010).
Looking at some of the earliest halo studies, which inevitably had to rely on much more limited samples of tracers, it is worth pointing out that, strikingly, glimpses of the variation of the halo shape were already caught by Kinman et al. (1966). This pioneering work took advantage of perhaps the most reliable halo tracer, the RR Lyrae stars (RRLs, hereafter). These old and metalpoor pulsating stars suffer virtually no contamination from other populations of the Milky Way and have been used to describe the Galactic halo with unwavering success over the last 50 years (see e.g. Hawkins, 1984; Saha, 1984; Wetterer & McGraw, 1996; Ivezić et al., 2000; Vivas & Zinn, 2006; Catelan, 2009; Watkins et al., 2009; Sesar et al., 2010; Akhter et al., 2012; Soszyński et al., 2014; Torrealba et al., 2015; Soszyński et al., 2016).
While deep, widearea 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 PanSTARRS1 (PS1, Hernitschek et al., 2016; Sesar et al., 2017), they have yet to be used to model the Galactic halo globally. In case of CSS, this might be due to the varying completeness of the sample. For PTF and PS1, this is sadly due to the public unavailability of the data. To remedy this, here we attempt to extract an allsky sample of RRLs from the Gaia Data Release 1 (GDR1, Gaia Collaboration et al., 2016b) data. Our primary goal is to use the thus procured RRL candidates to model the global properties of the MW stellar halo. Therefore, we are not concerned with maximizing the completeness but instead strive to achieve homogeneous selection efficiency and reasonably high purity. While GDR1 does not contain any explicit variability information for stars across the sky, Belokurov et al. (2017) and Deason et al. (2017) show that likely variable objects can be extracted from the GaiaSource table available as part of GDR1. We build on these ideas and combine Gaia and Two Micron All Sky Survey (2MASS) photometry (and astrometry) to produce a sample of RRLs out to kpc from the Sun, with constant completeness of and purity .
Armed with this unprecedented dataset, we simultaneously extract the radial density profile as well as the shape of the Galactic stellar halo. Furthermore, taking advantage of the stable completeness and the allsky view provided by Gaia+2MASS, we explore whether the density slope and the shape evolve with radius out to kpc. Finally, we also allow the halo to be i) arbitrary oriented ii) triaxial and iii) 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 selfconsistent 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 “phasemixed”, as suggested by cosmological body simulations (Helmi et al., 2011).
The paper is organized as follows. In Section 2 we describe the Gaia data as well as the method used to select an allsky sample of RRL candidates from a crossmatch 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 bestfit 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 allsky scanning space observatory, currently collecting multiepoch photometric and astrometric measurements of about a billion stars in the Galaxy. More details on the Gaia mission and on GDR1 can be found in Gaia Collaboration et al. (2016a). In the first data release, the information available for most faint sources is limited to basic properties, such as positions on the sky and fluxes in the broad Gaia band, which covers most of the visible spectra from approximately 400 nm to 10000 nm (Jordi et al., 2010).
In this work, we used the table GaiaSource released as part of the GDR1 (Gaia Collaboration et al., 2016b). GaiaSource contains a number of auxiliary pieces of information, which provide plenty of added value to the GDR1. For example, the errors on the mean flux measurements can be used to separate constant and variable sources, and even gauge the amplitude of the variability (see, e.g. Belokurov et al., 2017; Deason et al., 2017). Moreover, the quality of the astrometric fit, encapsulated by the socalled astrometric excess noise, contains information regarding the morphology of the source, and can be used to separate stars from galaxies (see Koposov et al., 2017). The relevant GaiaSource quantities used here (other than the sky coordinates RA, Dec), are:

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

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

, the standard deviation of the flux measurements;

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

AEN, the astrometric excess noise, which measures strong deviations from the best astrometric solution. The AEN should be large for objects whose behaviour deviates from that of pointlike sources, as, for example, unresolved stellar binaries or galaxies (see Lindegren et al. (2012) for details).
Additionally, relying on the crossmatch between Gaia and 2MASS, we calculated:

PM, the total proper motion of each object. PM, where and are the proper motions measured along RA and Dec, respectively.
2.2 RR Lyrae in Gaia Dr1
We use the parameters provided in GaiaSource to select RRLs from the GDR1. Following the method outlined in Belokurov et al. (2017) and Deason et al. (2017), we defined the quantity
(1) 
which can be used as a proxy for the amplitude of the stellar variability. Indeed, for variable stars with well sampled light curves, is proportional to the amplitude of the flux oscillation, while for nonvariable stars it is just a measure of photoncount Poisson errors. Thus, it is possible to set a threshold value for AMP that will select only variable candidates. For instance, Belokurov et al. (2017) showed that most variable stars like Cepheids and RRLs have .
The selection of variables through the AMP parameter suffers from the limitations that AMP is a timeaveraged information and does not allow to distinguish between various types of variable objects: RRLs, Cepheids or Mira variables. However, these different classes of pulsating stars populate a welldefined strip in the colourmagnitude 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 crossmatch failures). We discuss their importance in Sect. 2.2.6.
2.2.1 The Gaia + 2MASS sample
As mentioned, GDR1 reports photometric information only in the band. We derived a colour () for each source by crossmatching GaiaSource with the 2MASS survey data (Skrutskie et al., 2006) using the nearestneighbour 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 allsky coverage. The observed magnitudes have been corrected for extinction due to interstellar dust using the maps of Schlegel et al. (1998) and the transformation for the band (see Belokurov et al., 2017) and for the band (Fitzpatrick, 1999).
2.2.2 Auxiliary RR Lyrae datasets
In order to extract a reliable sample of RRL stars from the G2M catalogue, we must apply adhoc selection criteria. To this aim, we used two samples of bonafide RRLs: the CSS (Drake et al., 2013; Drake et al., 2014) and the Stripe 82 (S82, Sesar et al. 2010) catalogues. These samples allowed us to identify the optimal selection criteria, analyse the completeness and the contamination of the catalogue^{1}^{1}1The completeness indicates the fraction of recovered true RRLs as a function of the apparent magnitude, while the contamination is an estimate of the fraction of spurious objects (non RRLs) that “pollute” our sample. and estimate the RRL absolute magnitude in the band. The CSS contains about 22700 typeab RRL stars distributed over a large area of the sky (about 33,000 deg between and ) and extended up to a distance of 70 kpc. The completeness of this sample is constant (at ) only for , while it quickly decreases outside this range. Most importantly, as shown in Fig. 13 of Drake et al. (2013), for objects fainter than , the completeness is a strong function of the number of observations and thus varies appreciably across the sky. SDSS’s Stripe 82 covers a 2.5wide and long patch of sky aligned with the celestial equator and contains “only” 483 RRLs. However, the sample is very pure (with less than of contaminants) and complete up to a distance of 100 kpc.
The large number of stars in CSS is useful to define the selection criteria (see Sec. 2.2.3) and to estimate the absolute magnitude in the band (see Sec. 2.2.4), while the high quality of S82 sample is ideal to analyse the completeness and contamination of our final sample (see Sec. 2.2.5 and Sec. 2.2.6). A crossmatching of the CSS and S82 catalogues with G2M using an aperture of 1 arcsec led to the two samples CSS+G2M (GCSS hereafter) and S82+G2M (GS82 hereafter).
2.2.3 RR Lyrae selection cuts
In this section we describe how the final sample of RRLs was obtained from the G2M catalogue. The selection was driven by the properties of the bona fide RRLs in the GCSS and GS82 catalogues (see Sec. 2.2.2) in order to maximise completeness of the sample and its spatial uniformity, while keeping the level of contamination low (see Sec. 2.2.5).
First of all, in order to exclude a region likely dominated by the Galactic disc, we removed all the stars in the G2M catalogue located between the Galactic latitudes and . Our limit in latitude () is lower with respect to other works in literature (e.g. Deason et al. 2011; Das et al. 2016), in which, however, the choice was mainly motivated by the limited skycoverage of the used survey. Given the unprecedented skycoverage of our data sample, we decided to push forward the study of the halo structure exploring also the region at low Galactic latitude that is usually not well sampled by other surveys. We are aware that our final sample could be polluted by Galactic disc contaminants, but in the following sections we carefully analyse the level of contamination and all our results are obtained taking into account the possible biases due to stars of the Galactic disc.
Figure 1 shows the distribution of the G2M stellar density (yellowbluepurple colourmaps), a randomly selected subsample of RRLs in GCSS (red points) and in GS82 (orange squares) in the AMP (lefthand panel), colormagnitude (middle panel) and AEN (righthand panel) planes. The bona fide RRLs occupy a well defined strip in colour, thus we excluded all the stars with the colour index greater than and lower than as shown by the vertical black lines in the middle panel of Fig. 1. It is worth noting that most of the “normal” stars occupy this colour interval, therefore this cut mostly eliminates artefacts. The lefthand panel shows that the genuine RRLs are almost uniformly distributed in the AMP range of the G2M sample, however the contamination by spurious objects increases rapidly for AMP (see Sec. 2.2.5 and Fig. 5), thus we only retained stars with variability amplitudes above this value. With regards to completeness, the faint magnitude limit plays an important role. According to our analysis, is the faintest magnitude that we can reach to obtain a sample with spatially uniform completeness (see Sec. 2.2.5 for further details). The number of bright stars with , corresponding to RR Lyrae with distances less than 1 kpc from the Sun, is very small compared to the number of objects in our final catalogue. Therefore, instead of extending our completeness/contamination analysis at the very bright magnitudes (see Sec. 2.2.5 and 2.2.6), we decided to put the bright magnitude limit at .
The selection criteria described above involving colour, AMP and magnitude have the largest impact on the definition of our sample of RR Lyrae. However, we also applied a few minor refinements. The righthand panel of Fig. 1 shows that most of the bona fide RRLs have a very small value of AEN, so we excluded all sources with as shown by the horizontalblack line. This cut likely removes contaminant extragalactic objects since they typically have (Belokurov et al., 2017) and some of the eclipsing binaries that survive the colour selection. Additionally, to further clean the sample from possible nearby Galactic disc contaminants, we cull all the stars with a total proper motion, PM, greater than mas/yr. Given differences in the lightcurves and its sampling, the significance of AMP (Eq. 1) might depend on the number of photometric measurements . With this in mind, we impose : the focal plane of Gaia has an array of 97 CCDs, so all objects with less than 3 complete Gaia transits are excluded.
It would be useful to have an estimate of the photometric metallicity to retain only genuine metal poor stars from the halo and effectively exclude metal rich contaminants from the Galactic disc. However, the photometric metallicity estimate requires a basic knowledge of the shapes of the RRL light curves which is not available in our dataset (see e.g. Jurcsik & Kovacs 1996).
Finally, we masked a few regions of the sky. First we removed the area near the Magellanic Clouds using two circular apertures: one centred on with an angular radius of for the Large Magellanic Cloud (LMC) and the other centred on with a radius of for the Small Magellanic Cloud (SMC). By inspecting the sky distribution of the stars in our RRL sample, we noticed the presence of two extended structures (S1 and S2, hereafter), that were not connected to any known halo substructures, but are likely objects instead produced by Gaia crossmatch failures (see Sec. 2.2.5) that “survived” our selection cuts. We decided to mask these sky regions as well by removing all the stars in the following boxes for S1 and for S2. The final sample has been further cleaned to exclude “hot pixels” using a simple medianfilter method. We first built a skymap using pixels of , then we replaced the number of stars in each pixel by the median of the starcounts calculated in a squared window of four pixels. Finally, we calculated the ratio between the original skymap 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 medianfilter has been gauged to reveal smallscale features, since the most evident largescale structures have been already removed (LMC, SMC, S1 and S2). The “spotted” hotpixels correspond to known globular clusters (e.g. M3 and M5) or are connected to “remnants” of crossmatch failure structures (see Sec. 2.2.5). In order to fully exploit the allsky capacity of our sample, we decided not to exclude a priori any portion of the sky containing known substructures (e.g. unlike Deason et al. 2011 which masked out the Sagittarius stream). The analysis of the most significant substructures found in this work can be found in Sec. 3 and 5.
The final sample contains about 21600 RRLs that can be used to have a direct look at the distribution of stars in the Galactic halo (Sec. 3). The final number of object is similar to the one in the CSS catalogue, however we cover a larger area of the sky (almost allsky). Our sample populates 58% of the halo spherical volume within the Galactocentric distance of about 28 kpc which represents a significant improvement in volume fraction as compared to previous works (e.g. 20% in Deason et al. 2011). A summary of the applied selection cuts can be found in Tab. 1.
Selection cuts  

b [deg]  
[mag]  (10, 17.1) 
AMP  (0.7, 0.4) 
(0.95, 0.4)  
AEN  
PM [mas/yr]  
 [deg]  
Structure cuts  
LMC  
SMC  
S1  
S2  
21643 (13713)  
[%]  58 (44) 
2.2.4 Absolute magnitude and distance estimate
Despite photometric variability, RRLs have an almost constant absolute magnitude and, having the apparent magnitudes , we can directly estimate the heliocentric distances through
(2) 
where is the absolute magnitude in the Gaia band. To estimate we used the RRLs in GCSS, as they have heliocentric distances estimated from the periodluminosity relation. The resulting distribution is shown in Fig. 2: it has a well defined peak at and a small dispersion. We fit this distribution with a Gaussian function obtaining a mean of and a dispersion of , comparable to the uncertainties of the Vband absolute magnitude of RRLs (see e.g. Vivas & Zinn 2006). The above Gaussian function perfectly describes the data in the central parts, but the distribution shows broader wings for and . The objects that populate the wings could be Oosteroff type II RRLs and objects influenced by the Blazhko effect (Drake et al. 2013 and references therein). A better fit can be obtained using a double Gaussian model, where two Gaussian components peak at about the same absolute magnitude of the single Gaussian fit (see Fig. 2). Given the small dispersion around the mean, we decided to set the absolute magnitude for all the RRLs in our sample to a single value . An analysis of the effect of this approximation on the results can be found in Sec. 5.2.1.
2.2.5 Completeness
Before studying the properties of the stellar halo (as traced by RRLs), it is fundamental to consider the completeness of our sample of RRLs.
First of all, we checked that the scanning law of Gaia does not cause an intrinsic decrease of the completeness at low Galactic latitude. We compared the number of objects in GaiaSource with the number of stellar sources in the Data Release 7 of the Sloan Digital Sky Survey (SDR7, Abazajian et al. 2009) in a series of stripes at fixed Galactic longitudes, selected using the footprint of SDR7. The top panel of Fig. 3 shows the position of the stripes in Galactic coordinates. We selected all stellar sources in the SDR7 with , where is the band magnitude. In this magnitude range the SDR7 can be considered 100% complete^{2}^{2}2http://classic.sdss.org/dr7/products/general/completeness.html. We chose the band because the peak of the filter response is almost coincident in wavelength with the one of the Gaia band (see Sec. 2.1), therefore the two magnitudes are directly comparable. The bottom panel of Fig. 3 shows the ratio between the number of stars in Gaia and SDR7 in bins of Galactic latitude for individual SDSS stripes and considering all the stripes together. The ratio does not show significant variations as a function of b with values between 0.9 and 1.0. We conclude that there is no evidence for strong intrinsic completeness variations in the Gaia catalogue for .
We estimated the completeness (and the contamination) using the RRLs in our auxiliary catalogues: CSS and S82 (Sec. 2.2.2). In particular, S82 represents a complete (100 %) and pure (99 %) catalogue of RRLs located up to 100 kpc from the Sun, so it is perfect to test both the contamination and the completeness of our sample. We compared the number of stars in the original S82 sample with the ones contained in the GS82 after the selection cuts described in Sec. 2.2.3, in bins of heliocentric distance. Fig. 4 shows the level of completeness as a function of magnitude/distance assuming different lower limits in AMP ( red triangles, blue circles, green diamonds) in the range of magnitudes . The results are in agreement with the distancebased estimate of completeness as shown by the dashed lines. The level of completeness is relatively low, ranging from about 15% for to about 30% for but it is reasonably constant up to about then it abruptly decreases to 0 at so we decided to conservatively cut our sample at (vertical black line in Fig. 4, see Sec. 2.2.3).
We also used the GCSS catalogue and found no significant variation of the completeness as a function of the Galactic sky coordinates l and b, although we found a mild trend for increasing . This is expected since the larger the number of flux measurements the better the sampling of the light curves and, as a consequence, the AMP cut (Eq. 1) improves its effectiveness in selecting RRLs. However, the increase is not dramatic as all the differences are within 10%. In conclusion, for the purpose of this work, we considered the completeness of our catalogue of RRLs uniform across our sky coverage and the considered magnitude range (see Tab. 1).
2.2.6 Contamination
We estimated the contamination of spurious sources as
(3) 
where and indicate the number of stars in our samples and in GS82 catalogue after the selection cuts in Sec. 2.2.3. The S82 sample is pure (Sesar et al., 2010), so all the stars in our sample that are not present in the GS82 catalogue are likely contaminants. As before, we considered the magnitude range for different lower limits of AMP. For AMP threshold lager than the level of contamination is lower than 10% (Fig. 5) with a mild increase toward low Galactic latitudes, where the contribution of the Galactic disc is larger. For AMP cut lower than the contamination fraction rapidly increases (about 25% for ). We expect that the contaminants of our sample are eclipsingbinaries in the Galactic disc and possible instrumental artefacts. As shown and discussed in detail in Belokurov et al. (2017) some regions of the Gaia allsky map are crossed by sharp strips with a large excess of objects. These strips are similar in width to the field of view of Gaia and are regularly spaced on the sky. Belokurov et al. (2017) propose that these are spurious features due to failures in the crossmatch procedure of Gaia, so that at some epochs the flux of a star comes from a different object. The spurious measurement of the flux increases and moves the stars in the region of high AMP “polluting” our samples. We found that for the contaminants are mostly related to crossmatch failures. Unlike the disc contaminants, the crossmatch failures have a complicated and poorly understood spatial distribution, so these structures are difficult to be taken in account in the study of the properties of the Galactic halo. For this reason, we decided to use as the lower limit in our selection cut (see Sec. 2.2.3).
3 Distribution of the RRLs in the inner stellar Halo
3.1 A first Gaia look at the stellar halo
Fig. 6 shows the distribution of the RRLs of our final sample as a function of the Galactic coordinates (l, b). We split the stars in the sample into magnitude intervals: , corresponding to heliocentric distances between 1 kpc and 10 kpc (Eq. 2), and , i.e. with heliocentric distances between 10 kpc and 20 kpc (assuming that all RRLs have absolute magnitude ). We stress that Fig. 6 represents the first allsky view of the distribution of the RRLs in the inner halo.
The first magnitude range covers a portion of the Galaxy mostly located in the side of the Galaxy containing the Sun between the Galactic radii 1.4 kpc and 18 kpc. The distribution of stars in this region is quite regular: as expected most of the stars are in the direction of the Galactic centre () and there are not evident asymmetries with respect to the Galactic plane. The second magnitude range covers the Galactic distances between 18 kpc and 28 kpc in the side of the Galaxy containing the Sun and between 3 kpc and 12 kpc in the other side. In this distance range the distribution of the RRLs is less regular with clear asymmetries: in particular, an excess of stars is evident at high Galactic latitude around . In both magnitude intervals an overdense band of stars at low Galactic latitude () can be seen running all around the Galaxy. A discussion of the nature of these structures can be found in Sec. 5.3.
3.2 Setting the frame of reference
We set a lefthanded Cartesian reference frame centred in the Galactic centre and such that the Galactic disc lies in the plane , the Sun lies on the positive Xaxis and the Sun rotation velocity is . The actual vertical position of the Sun with respect to the galactic disc is uncertain, but it is estimated to be be smaller than 50 pc (Karim & Mamajek 2017 and reference therein), and thus negligible for the purpose of this work. In this Cartesian reference frame, the coordinates of an object with Galactic longitude l, Galactic latitude b and at a distance from the Sun are
(4) 
where is the distance of the Sun from the Galactic centre. In this work we assume (Bovy et al., 2012), but the main results of our work are unchanged for other values of , within the observational uncertainties (between 7.5 kpc, e.g. Francis & Anderson 2014, and 8.5 kpc, e.g. Schönrich 2012).
For a star with Galactocentric Cartesian coordinates we define the distance from the Galactic centre
(5) 
the Galactocentric cylindrical radius
(6) 
the Galactocentric latitude
(7) 
and the Galactocentric longitude
(8) 
Assuming that the stellar halo is stratified on concentric ellipsoids, it is useful to introduce another Cartesian frame , aligned with the ellipsoid principal axes, and to define the elliptical radius
(9) 
where p and q are, respectively, the YtoX and ZtoX ellipsoid axial ratios. In general we will allow the to differ from in both orientation and position of the origin. When the origin of is the Galactic centre and , as we will assume in this section, the system can be identified with without loss of generality.
3.3 Density distribution of the halo RRLs
Given the allsky RRL sample illustrated in Fig. 6, we can compute their volume number density distribution . In particular we define the number density of halo RRLs in a cell centred at , where is a coordinate vector, such as for instance , as
(10) 
where is the number of stars observed in the cell, V is the volume of the cell and is the fraction of this volume accessible to our analysis, which depends on the sky coverage, on the selection cuts and the mask that we applied (see Sec. 2.2.3), and on the offset between the Sun and the Galactic centre. We estimate numerically as follows. First, we define a Cartesian Galactic grid sampling each of the axes with 300 points linearly spaced between and ( kpc is the maximum value of , given our magnitude limit at ), so the reference Galactic volume (a cube of 56 kpc side centred on the Galactic centre) is sampled with 27 million points with a density of about 125 points per kpc. Each of the points of the grid is a “probe” of the Galactic volume. We assign to each of them Galactic coordinates, heliocentric coordinates and observational coordinates (l, b, , ). A secondary (nonuniform) grid is built by applying the selection cuts in Tab. 1 to the points on the primary grid, so we end up with a grid sampling the complete Galactic volume and a grid sampling the portion of the volume accessible to our analysis. Given a cell in a certain coordinates space (e.g. R, ), we define as the ratio between the number of points of the secondary grid and the number of points of the primary grid contained in the cell. In summary, when we estimate the number density in a 1D or 2D space we build three binned maps: the first contains the number of stars in each bin/cell (), the second the total volume of each bin/cell (V) and the third the fraction of this volume that we are actually sampling (). These three quantities are inserted in Eq. 10 to estimate the stellar number density. We tested this method with both simple analytic distributions and mock catalogues (see Sec. 4.3.5).
3.3.1 Meridional plane
Fig. 7 shows the number density of the RRLs of our sample in the Galactic meridional plane . The shape of the isodensity contours clearly shows the presence of two components: a spheroidal component and a discy component. The discy component causes the flattening of the contours at low . The nature of this component is uncertain: it could represent the disc RRL stars, it could be a lowlatitude substructure of the MW halo or, finally, it could be due to nonRRL contaminants from the Galactic disc (both genuine variables such as eclipsing binaries as well as artefacts, e.g. due to crossmatch failures). A more detailed analysis of this lowlatitude substructure can be found in Sec. 5.3. The density at higher is less contaminated by the discy component and it represents more directly the density behaviour of the RRLs in the halo. The isodensity contours nicely follow the the elliptical contours (overplotted in Fig. 7) out to kpc. At larger radii the contours tend to be rounder. The overall density distribution looks reasonably symmetric with respect to the Galactic plane, although there is an overdense region at kpc, which does not seem to have a counterpart below the Galactic plane (see Sec. 3.3.4 for further details).
3.3.2 Density profile
Fig. 7 demonstrates that the RRLs in the inner halo are consistent with being stratified on spheroids (except in the region close to the Galactic plane), thus we estimate the 1D RRL density profile by counting the stars in spheroidal (, ) shells and dividing this number by the shell volume corrected by the coverage of our sample (Eq. 10). The profiles are shown in Fig. 8 for , and . Independently of the assumed value of q, the density of RRLs follows a single power law with no significant evidence of change of slope out to an elliptical radius of 35 kpc. It must be noted that the most distant region is only sampled by a low number of stars located in a small portion of the total volume, so it is possible that a change of slope starting near the edge of our elliptical radial range could be missed with this analysis. As the elliptical radius depends on the assumed axial ratios, the results of this analysis are not straightforwardly comparable with previous works (see Sec. 5.4 for a detailed comparison with results obtained in the literature).
The slope of this power law is similar in the three cases and very close to (see Fig. 8), but, based on this simple analysis, all values of in the range are consistent with the data.
3.3.3 Halo flattening
The distribution of RRLs in the meridional Galactic plane (Fig. 7) suggests that the inner halo should be reasonably well represented by a spheroidal stratification with . However, in order to more rigorously study the halo flattening, we estimate the density in the plane (see equations 4 and 9). In practice, for a given value of , we scan the density as a function of at fixed . If the RRLs are truly stratified on similar spheroids with axial ratio , the density is independent of , so the isodensity contours in the plane are vertical stripes. If the assumed value of q is smaller than the true value , is underestimated, the estimated density is a monotonic decreasing function of and the isodensity contours in the plane are bent in the direction of increasing with (provided that the density is a decreasing function of ). The isodenisty contours are bent in the opposite direction, if one assumes . The shape of the isodensity contours in the plane is a very efficient and direct diagnostic of the evolution of q as function of the elliptical radius.
The number density maps in the plane of the RRLs in our sample are shown in Fig. 9 for , and , assuming . Below (indicated by the whitedashed 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 () the contours give a direct indication on the flattening of the halo: in the righthand panel of Fig. 9 () the isodensity contours are significantly inclined in a way that implies . For kpc a flattening of about 0.55 gives a good description of the data as shown by the vertical isodensity contours in the middle panel of Fig. 9 (), but beyond 20 kpc the contours start to bend so that . The last two isodensity contours in the lefthand panel () look vertical enough to assert that at the outer radii the halo becomes more spherical. This analysis, together with the recent works of Liu et al. (2017), shows the first direct evidence of a change of shape of the stellar halo going from the inner to the outer halo. Moreover the unique allsky view of our sample allows us to confirm that this trend is symmetric with respect to the Galactic plane. A variation of the halo flattening was also proposed in previous works (e.g. Xue et al. 2015 and Das et al. 2016). In Sec. 4.4 we perform a comprehensive model fitting analysis of the RRL dataset and compare our results to those found in the literature.
3.3.4 Vertical asymmetries
Contours of the RRL density shown in Fig. 7 and Fig. 9 display an overall symmetry between the Northern and Southern Galactic hemispheres, however at kpc and around an 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 slabs to understand whether the overdensity is compatible with the Poisson starcount fluctuations or it is due to a genuine halo substructure.
Fig. 10 shows starcount maps in the plane integrated along three different slabs. Close to the plane (lefthand panels) the two maps look very similar and the difference between the total number of stars is small (less than 2%) and compatible within the Poisson errors. In these maps, an elongated component (also visible in the other density maps, Figs. 7 and 9) can be seen stretching from the Galactic centre to the Sun and beyond: this component appears symmetric with respect to the Galactic plane but is present only in the side of the Galaxy containing the Sun. In the intermediate slabs (middle panels) the difference between the number of stars in the region above and below the Galactic plane is small (less than 5%) as in the previous case, but the maps look less symmetric. In particular the peaks of the star counts are now apart by about 10 kpc. Finally, the slabs at the highest (righthand panels) show a significant difference in the total number of stars (about 20%) that can not be explained by Poisson fluctuations only. Moreover the excess of stars above the Galactic plane is strongly clustered in the regions between kpc and kpc.
Note that some of the differences between the star counts in the regions above and below the Galactic plane could be due to the mask used to eliminate the contribution of the Magellanic Clouds (see Sec. 2.2.3). Indeed, some mismatch is expected given that we have excluded a relatively large region of the halo volume below the Galactic plane. In order to quantify the differences introduced by the mask, we produce mock catalogues (see Sec. 4.3.5) of different halo models and find that the Magellanic Clouds mask introduces a difference of about the 7% in the number of objects above and below the Galactic plane for kpc, significantly less than the value of 20% obtained here. Therefore, the structure seen at high positive appears to be genuine, and is most likely related to the “Virgo overdensity” (Jurić et al., 2008; Vivas et al., 2016). Across all slabs, portions of the Virgo Cloud are visible as an excess of stars at positive . Virgo’s counterpart underneath the disc is the HerculesAquilla Cloud (see Belokurov et al., 2007; Simion et al., 2014), discernible in middle bottom panel as strong overdensity at negative (and ).
4 Model fitting
In this section we present the stellar halo models and we compare them with the observed sample of RRLs.
4.1 Clean sample
Fig. 7 and Fig. 9 show the presence of a highly “flattened” structure close to the Galactic plane. The properties of this structure are clearly at odds with the distribution of stars at high Galactic latitude that are more likely a “genuine” tracer of the halo population. The “flattened” component contains about 35% of the RRLs in our sample, so it must be taken into account to infer the properties of the stellar halo. Therefore, we built a “clean” sample of RRLs eliminating all the stars belonging to the substructure from our original catalogue (Sec. 2.2.3). Fig. 9 suggests that the substructure can be effectively eliminated with a selection in angle (see Eq. 7).
In particular, at there is a transition between a very flattened component and a more spheroidal structure. Therefore, we define our clean sample of halo RRLs as the stars with : this subsample contains 13713 objects and covers approximately 44% of the halo volume within a sphere of radius 28 kpc (see Tab. 1). We stress that this cut is based on values of obtained from Eq. 2 and Eq. 7 assuming the same value of the absolute magnitude, , for all the stars (see Sec. 2.2.4).
We also tried to exclude the flattened structure with alternative cuts, e.g. using a higher cut on the Galactic latitude b or a direct cut on . We found that the results obtained for the sample with (see Sec. 4.4) are qualitatively similar to the results obtained using a sample with (the lowest Galactic latitude used in most of the previous works, e.g. Deason et al. 2011) or a sample with kpc (Fermani & Schönrich 2013 uses kpc to cut disc stars). The cuts on b and significantly reduce the number of tracers in the inner part of the halo and the final number of stars is smaller, in both cases, with respect to the one obtained cutting the original sample at .
In the following subsections we present the method used to fit halo models to this subsample of RRLs and the final results.
4.2 Halo models
As in Sec. 3.1, we assume here that the number density of the halo RRLs is stratified on ellipsoids with axial ratios p and q. Therefore, a halo model is defined by a functional form for the density profile (number density as a function of the elliptical radius) and a “geometrical” model for the ellipsoidal isodensity surfaces (in practice, characterized by values of p and q, and the orientation of the principal axes with respect to the Galactic disc).
4.2.1 Density profiles
We consider five families of number density profiles: double powerlaw (DPL), single powerlaw (SPL), cored powerlaw (CPL), broken powerlaw (BPL) and Einasto profiles (EIN). The DPL profile has a number density
(11) 
where and indicate the inner and outer powerlaw slopes and is the scale length. The number density of the SPL is given by Eq. 11 when , while the number density of the CPL has , in which case represents the length of the inner core. The BPL number density profile is given by a piecewise function:
(12) 
The EIN profile (Einasto, 1965) is given by
(13) 
where for (Graham et al., 2006). The steepness of the EIN profile, , changes continuously as a function of tuned by the parameter n,
(14) 
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 powerlaw profiles with or with imply halos with infinite total mass, however our study focuses only on a limited radial range and we do not exclude a priori any solution. We anticipate that our best density model is a SPL with (see Sec. 4.4), therefore there should be a physical radius, outside our radial range, beyond which the profile becomes, either abruptly (e.g. BPL or an exponential truncation) or gently (e.g. DPL), steeper.
4.2.2 Isodensity ellipsoidal surfaces
Concerning the isodensity ellipsoidal surfaces we define four different models:

discnormal axisymmetric (DN): we set , the axis of symmetry is normal to the Galactic disc, and q is a freeparameter;

discplane axisymmetric (DP): we set , the axis of symmetry is within the Galactic plane making an (anticlockwise) angle with respect to the Galactic Yaxis, and p is a freeparameter,.

triaxial (TR): both p and q are considered freeparameters, the Zaxis 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 isodensity ellipsoidal surfaces:

qvarying (qv): q depends on the elliptical radius as
(15) so q varies from at the centre to the asymptotic value at large radii and the variation is tuned by the exponential scale length ;

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

offset (off), in this model the elliptical radius of each star is estimated with respect to a point offset by (, , ) with respect to the Galactic centre.
The specific functional form of Eq. 15 is empirical: the same expression was adopted by in Xue et al. (2015) and Das et al. (2016), where, however, q is a function of the spherical radius (). We decided to maintain the dependence on given that this approach is selfconsistent with our assumption that the RRLs are stratified on ellipsoidal surfaces. If follows that for our models with a varying q, the elliptical radius of a star with coordinates is the root of
(16) 
where is defined in Eq. 15. In our fitting code we solve this equation numerically with a NewtonRaphson root finder.
4.2.3 Complete halo model
Each complete halo model is defined by a model for the density law (SPL, DPL, CPL, BPL, EIN) plus a model for the shape of the isodensity surfaces (SH, DN, DP, TR) and any combination of geometrical variants (qv, tl, off). For instance, a SPLSH model has a spherical distribution of stars with a single powerlaw density profile, while a model is a triaxial tilted model with varying flattening along the Zaxis of symmetry and an Einasto density profile (see Tab. 2 for a reference on the various model labels). In the next sections when we discuss the properties or the results of density models (e.g. SPL models), unless otherwise stated, we are implicitly referring to all the models that share the same density model whatever the geometrical and geometrical variants properties are. The same applies when we focus on geometrical models only. We fit our data with all possible combinations of density laws (SPL, DPL, CPL, BPL, EIN), geometrical models (SH, DN, DP, TR) and model variants (qv, tl, off). See Tab. 3 for a summary of results obtained with a sample of such models.
4.3 Comparing models with observations
4.3.1 Density of stars in the observed volume
Given a certain halo model (Sec. 4.2), the normalised number density of stars in an infinitesimal volume is given by
(17) 
where is the vector containing all the model parameters (e.g. for an SPL+discnormal axisymmetric model, see Sec. 4.2) and is defined such that
(18) 
It is useful to define the normalised star number density
(19) 
within the infinitesimal projected (onto the sky) volume element centred at (,l,b), where is the observed magnitude, l and b are the Galactic longitude and latitude, and
(20) 
is the determinant of the Jacobian matrix (Appendix A). The number density depends also on the additional parameter , which is the absolute magnitude needed to pass from the observable variables to the Cartesian variables (see Eqs. 2 and 4). Substituting Eq. 20 in Eq. 19 we obtain
(21) 
4.3.2 Likelihood of a single star
From the normalised density function (Eq. 21) we can define the expected rate function for finding a star with (, l, b) given the absolute magnitude and a halo model with parameters
(22) 
In Eq. 22, is a constant, represents the probability density function (pdf) of the absolute magnitude of the stars, while the function is the selection function that takes into account the incomplete coverage of the Galactic volume. It is a function of l, b and and returns a result in the Boolean domain . In particular, is always equal to 1 except for the points (, l, b) that are outside the volume covered by our clean sample (see Tab. 1).


Density parameters  Isodensity surfaces parameters  

/n  (kpc)  p  q()  (kpc)  (kpc)  (degree)  
Single powerlaw  SPL  U[0,20]  
Cored powerlaw  CPL  U[0,20]  U[0.01,10]  
Double powerlaw  DPL  U[0,20]  U[0,20]  U[,]  
Broken powerlaw  BPL  U[0,20]  U[0,20]  U[,]  
Einasto  EIN  U[0,100]  U[0.1,500]  
Spherical  SH  
Discnormal axsym  DN  U[0.1,10]  
Discplane axsym  DP  U[0.1,10]  U[80,80]  
Triaxial  TR  U[0.1,10]  U[0.1,10]  U[80,80]  
qvar  qv  U[0.1,100]  
Offset  off  U[0,10]  
Tilted  tl  U[80,80] 
In this work, we decided to set the absolute magnitude of the RRLs to a single value, (Sec. 2.2.4), so we are assuming that in Eq. 22 is a Dirac delta. As a consequence, we can marginalize Eq. 22 over to obtain
(23) 
and we define the normalised rate function as
(24) 
The normalised rate function in Eq. 24 is the pdf of stars, i.e. the likelihood per star, at a certain position (, l, b) for halo model parameters .
4.3.3 The total likelihood
Consider a sample of stars with coordinates (, l, b) and a halo model with parameters . Given the pdf of the stellar distribution (Eq. 24) the logarithmic likelihood is
(25) 
where is the number of stars in our sample. Plugging Eq. 21 and Eq. 23 into Eq. 24 we can write the logarithmic likelihood as
(26) 
where
(27) 
is the normalisation integral. Notice that the numerator of Eq. 26 is evaluated only in regions of the sky where (i.e. where we observe stars).
4.3.4 Sampling of the parameter space
Considering the Bayes’s law, the logarithmic posterior probability of the parameters of an halo model given the data is
(28) 
where is the probability of the data, , given the parameters (See sec. 4.3.3) and represents the prior probability of (Tab. 2). In Eq. 28 we omit the Bayesian evidence term that is defined as the integral of the likelihood over the whole parameter space. This term is negligible in the determination of the best set of parameters for a given model (see below).
Using Eq. 26 we can write the posterior probability as
(29) 
In order to explore the parameter space without biases, we decided to use very large priors for most of our parameters (see Tab. 2). In particular, the very large ranges for and approximate uniform infinite priors.
Concerning the parameter , we assume, for both the BPL and DPL density models, that the prior distribution is uniform between and which represent the minimum and maximum elliptical radii (Eq. 9) of stars in our sample. For most of our models, remains constant at about 2.4 kpc, while ranges between roughly 31 kpc and 38 kpc.
We analysed the posterior probability of the parameters, , for a given halo model using the Eq. 29 and sampled the parameters space with Goodman & Weare’s Affine Invariant Markov chain Monte Carlo (mcmc, Goodman & Weare 2010), making use of the Python module emcee^{3}^{3}3https://github.com/dfm/emcee (ForemanMackey et al., 2013). The technical details of our approach can be found in Appendix B.
The final bestfit values of the model parameters have been estimated using the 50th percentile of the posterior distributions and the 16th and 84th percentiles have been used to estimate the 1 uncertainties.
4.3.5 Tests on mock catalogues
In order to test our fitting method, we developed a simple Python script to generate mock catalogues with the same properties (e.g. number of stars, magnitude limits) as our clean sample of RRLs (Sec. 4.1). This algorithm distributes the stars using a combination of the density laws (Sec. 4.2.1) and according to the assumptions on the properties of the isodensity ellipsoidal surfaces (Sec. 4.2.2). The absolute magnitude of the stars in the mock catalogues are distributed using the bestfit double Gaussian functional form shown in Fig. 2. The final mock catalogues do not include Gaia and 2MASS photometric and skycoordinates errors because the uncertainties they cause on the estimate of the distance are negligible.
We applied our fitting method to mock catalogues finding that we are able to recover the input parameters for all possible model combinations, including in the sample a fraction of Galactic disc contaminants (see Sec. 5.2.3).
Analysing the mock catalogues, we can estimate what properties of the halo we are able to measure and constrain using the stars in our catalogue. We found that we are able to detect a core in the density profile (CPL density model) all the way down to 100 pc. However for a very small core size most of the information comes from a relative small region where the transition between the inner core and the outer density profile takes place. The detection of a change of slope depends on the halo flattening: for and we can detect a break in the density profile within about kpc (Eq. 21). For more flattened models (q around ) this range extends up to about kpc. However, in this latter case most of the information for regions where kpc comes from a small number of stars at very high Galactic latitude, therefore the fit can be easily biased by the presence of substructures such as the one highlighted in Sec. 3.3.4 and Fig. 9. In conclusion, we are confident that we are able to robustly detect significant deviations from a SPL within a range of elliptical radii from less than 1 kpc to about 30 kpc. Finally, we verified that the method described above can recover a variation of the halo flattening given a realistic mock dataset.
4.4 Results
In this section we present the main results obtained applying the method described in Sec. 4.3 to the RRLs in our clean sample (Sec. 4.1).
Model  
Density law  Surface  Parameters  
SPL  SH  868  1688  
SPL  DP 

405  779  
SPL  DN 

81  124  
BPL  DN 

79  141  
DPL  DN 

81  143  
CPL  DN 

81  134  
EIN  DN 

87  145  
SPL 

66  113  
SPL  TR 

22  23  
SPL 

0  0 
4.4.1 Model comparison
In addition to estimating the parameters for a single model, it is important to compare the results of different models to determine which of them gives the best description of the data. The most robust way to perform a model comparison is through the ratio of the Bayesian evidences. Under the assumption that the posterior distributions are almost Gaussian, the Bayesian evidence can be approximated by the Bayesian information criterion (BIC, Schwarz 1978) defined as
(30) 
where is the number of objects in the sample and is the maximum value of the likelihood (Eq. 26).
The BIC is often used to compare different models with different dimensions in parameter space and the model with the lowest BIC is preferred. The BIC is similar to the maximum likelihood criterion, but it includes a penalty depending on the number of free parameters, such that for two models with the same likelihood, the one with more parameters is penalised. The bestfit parameters together with a comparison of the maximum likelihood and the BIC for a representative sample of halo models can be found in Tab. 3. These quantities have been calculated taking the maximum of the likelihood estimates obtained with the final mcmc sampling (see Sec. 4.3.4).
Other than the BIC, we also compared the ability of the different models to describe the observed properties of the RRLs (e.g. distribution of the stars in the sky).
4.4.2 RRLs density law
We tested different density profiles assuming a DN geometrical model: the SPL (Eq. 11 with ), the DPL (Eq. 11), the BPL (Eq. 12), the CPL (Eq. 11 with ) and the EIN profile (Eq. 13).
The logarithmic maximum likelihoods and the BICs obtained for these different models are very similar as shown in Tab. 3, moreover Fig. 11 shows that the predicted fraction of RRLs in bins of magnitude, for the different density models, are practically coincident. Therefore, the DPL, BPL, CPL and EIN models do not offer any significant improvement in the description of the RRLs distribution with respect to the simpler SPL. Fig. 12 shows the comparison between the (elliptical) radial profiles of the density slope for different halo models. The slopes have been calculated as the logarithmic (elliptical) radial derivative of the logarithmic density, therefore the result for a SPL is just a constant (as shown by the blackdashed line in the four panels). The bestfit broken radius for the BPL is kpc: in the inner part the bestfit slope is practically the same as the SPL model then it decreases to a slightly larger slope of about 2.9.
Concerning the DPL, the posterior distributions of the slopes are compatible with , such that the final density profile is again an SPL with a slope compatible with the bestfit SPL model (2throw panel in Fig. 12), moreover the posterior distribution of is almost uniform in the prior range interval (Tab. 2). Similarly, the core radius of the CPL model is very small and the final CPL model is compatible with nocore and the outer slope is the same of the SPL model (3throw panel in Fig. 12).
Finally, for the EIN the fit favours very large values both for n () and for ( kpc), so that in the analysed radial range the variation of the slope is minimal and the density profiles effectively mimic an SPL. However, since the BIC of the EIN model is higher than the BIC of the SPL models, there is no evidence in favour of an EIN profile for the inner part of the stellar halo.
The break radius of the BPL model is compatible with the one obtained in Xue et al. (2015) after the subtraction of halo substructures, although the change of slope is much more significant in their work with respect to our result (see Tab. 4). In principle, the overdensity of stars described in Sec. 3.3.4 might make it difficult to detect the break, so we repeated the maximum likelihood analysis using the RRLs only above and below the Galactic plane. The results are shown in Fig. 12 by the solid () and dashed () coloured curve: in general the stars above the Galactic plane prefer a mean slope of about 2.65 while the ones below have a mean slope of about 2.78. Concerning the change of slope above the Galactic plane, the posterior distribution of the parameters of the BPL, DPL and CPL models are such that these models are similar to the SPL profile. Below the Galactic plane both the BPL and the DPL favour a change of slope, however the one of the BPL model is much more significant going from about 2.73 in the inner part to about 3.3 beyond kpc. It is evident that the stars below the Galactic plane have a steeper decrease of the density with respect to the stars above the Galactic plane. This result could be due to the excess of stars at high latitude, however in both the 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 a SPL density law: the RRLs follows an SPL with an exponent that ranges between 2.6, assuming a SH model, to 3 assuming a halo. These results agree with the estimate of the density profiles obtained in Sec. 3.3.2 counting stars in ellipsoidal shells (Fig. 8).
4.4.3 Isodensity surfaces
As shown by the BIC comparison in Tab. 3, the SH and DP models give significantly worse fits with respect to the DN and TR models. This result is confirmed by a comparison of the properties of the RRLs in our sample with the ones expected for the bestfit models. For example, Fig. 13 shows the fraction of RRLs in our sample in different observed magnitude bins (black points) compared to the predicted distribution for different halo models (assuming a SPL, see Tab. 3). The SH and the SP models predict too low a fraction of stars at lower and a significant excess around , while the DN and the TR models give a good match to the data.
As a further refinement we tested the and models. Allowing q to vary, we obtained a better description of the distribution of RRLs as shown by the decreases of the BICs in Tab. 3. The variation of the flattening is quite significant as it is shown in Fig. 14. The halo is largely flattened () in the very inner part and then it becomes more spherical reaching a flattening of about 0.75 at the border of the Galactic volume analysed in this work ( kpc). These results are in agreement with what we have seen directly in the distribution of stars in the  plane in Sec. 3.3.3 and Fig. 9. Essentially the same result is obtained by considering the stars only above or only below the Galactic plane. Fig. 15 shows the comparison between the magnitude distribution of the stars in our sample and the one expected for the varying q model.
In the (nontilted) TR and the DP models the elongated axis makes an (anticlockwise) angle with respect to the positive part of the Galactic Yaxis. 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 ^{4}^{4}4 is defined as the anticlockwise angle between the direction of the bar and the Galactic Xaxis. is uncertain (Antoja et al., 2014), depending on the works it ranges from about (e.g. Robin et al. 2012) to about (e.g. Benjamin et al. 2005). Given these large uncertainties, the angle can be considered compatible with . As a consequence, the elongated axis (Yaxis) 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 Sec. 4.2.2). We added these refinements to the SPLTR model.
For the offset model, the bestfit parameters are

kpc,

kpc,

kpc,
and the size of the displacement vector with respect to the Galactic centre is . The offset we found is small and quite insignificant compared to the radial extent of the halo. Indeed these small values could also be due to the overfitting of some local substructure at small elliptical radii. Moreover, both the posterior distribution of and are compatible (within and , respectively) with 0 and most of the displacement is along the Galactic Xaxis. This axis connects the Sun to the Galactic centre, therefore the offset in this direction can be also interpreted as changing the distance of the Sun from the Galactic centre. Considering our results, the estimated Sun position is kpc, which is totally compatible with the Sundistance 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 bestfit results for the (anticlockwise) angles describing the orientation of the halo (see Sec. 4.2.2) are

(rotation around the Zaxis),

(rotation around the new Yaxis),

(rotation around the final Xaxis).
The angle represents the tilt of the elongated axis of the halo with respect to the Galactic Yaxis and it is compatible with the orientation found for the nontilted DN and TR discplane models (Tab. 3). The angles and measure the tilt of the halo meridional plane with respect to the plane of the Galactic disc. Even though the BIC indicates a mild improvement with respect to the nontilted triaxial model (), the tilt is quite small.
Finally, we compared the ability of tilted and nontilted models to reproduce the observed distribution of stars. Fig. 15 shows the comparison between data and different models for the magnitude distribution of the stars. The models are , , and . Note that all of them are capable to give a good match to the data.
In conclusion, we conservatively assume that there is no need for a significant tilt to describe the RRL distribution in the inner halo. A more detailed study about this interesting properties is postponed to future works in which we will extend the coverage of the Galactic volume (see Sec. 5.5).