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

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

###### 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 ($\sim 60\%$) of the volume of the Galactic
inner halo ($\text{R}<28$ kpc). We study the morphology of the stellar halo by
analysing the RR Lyrae distribution with parametric and non-parametric
techniques. Taking advantage of the uniform all-sky coverage, we test halo
models more sophisticated than usually considered in the literature, such as
those with varying flattening, tilt and/or offset of the halo with respect to
the Galactic disc. A consistent picture emerges: the inner halo is well
reproduced by a smooth distribution of stars settled on triaxial density
ellipsoids. The shortest axis is perpendicular to the Milky Way’s disc, while
the longest axis forms an angle of $\sim 70^{\circ}$ with the axis connecting
the Sun and the Galactic centre. The elongation along the major axis is mild
($\text{p}=1.27$), and the vertical flattening is shown to evolve from a
squashed state with $\text{q}\approx 0.57$ in the centre to a more spherical
$\text{q}\approx 0.75$ at the outer edge of our dataset. Within the radial
range probed, the density profile of the stellar halo is well approximated by
a single power-law with exponent $\alpha=-2.96$. 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)

††pubyear: 2017††pagerange: The first all-sky view of the Milky Way stellar
halo with Gaia+2MASS RR Lyrae–B

##  1 Introduction

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

Leaving the ideas of James Jeans aside, could the structural parameters of the
stellar halo alone inform our understanding of the mass assembly of the
Galaxy? At our disposal are the numbers pertaining to the slope of the stellar
halo’s radial density profile and its vertical flattening (see e.g. Jurić et
al., 2008; Bell et al., 2008; Deason et al., 2011; Xue et al., 2015). The
radial profile has so far been measured with a variety of stellar tracers.
Studies based on Main Sequence Turn-Off stars (e.g. Sesar et al., 2011; Pila-
Díez et al., 2015), Blue Straggler and Horizontal Branch stars (e.g. Deason et
al., 2011) and RR Lyrae (e.g. Sesar et al., 2007; Watkins et al., 2009) seem
to favour a “broken” profile. According to these datasets, somewhere between
20 and 30 kpc from the MW centre, the density slope changes from a relatively
shallow one, as described by power-law index of approximately $-2.5$, to a
much steeper one, consistent with a power-law index of $\approx-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).
Nonetheless, a consistent picture is now emerging: in line with other pieces
of evidence, the “broken” MW stellar halo appears to be the tell-tale sign of
an early-peaked and subsequently quiescent accretion history. Note also that
such destitute state of the stellar halo is not permanent but rather transient
(see e.g. Deason et al., 2016; Amorisco, 2017a). In agreement with
simulations, the MW stellar halo is destined to transform dramatically with
the dissolution of the debris from the Sgr dwarf and the Magellanic Clouds.

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

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

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

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

The analysis of the density distribution of the stars in our sample is based
on the fit of density models, rather than on the fit of full dynamical models
(see e.g. Das & Binney 2016; Das et al. 2016). The main reason behind this
choice is that we do not have any kinematic information, so the use of self-
consistent dynamical models does not add any significant improvement to our
study. Moreover, the knowledge of the spatial density distribution of the
stellar halo is a useful piece of information not only if the halo is
stationary, but also if it is not “phase-mixed”, as suggested by cosmological
$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-fit halo model is discussed together with the possible
biases that can affect our results. The summary of the results can be found in
Section 6.

##  2 The RR Lyrae Sample

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

###  2.1 Gaia Data release 1

Gaia is an all-sky scanning space observatory, currently collecting multi-
epoch photometric and astrometric measurements of about a billion stars in the
Galaxy. More details on the Gaia mission and on GDR1 can be found in Gaia
Collaboration et al. (2016a). In the first data release, the information
available for most faint sources is limited to basic properties, such as
positions on the sky and fluxes in the broad Gaia $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 value to the GDR1. For example,
the errors on the mean flux measurements can be used to separate constant and
variable sources, and even gauge the amplitude of the variability (see, e.g.
Belokurov et al., 2017; Deason et al., 2017). Moreover, the quality of the
astrometric fit, encapsulated by the so-called astrometric excess noise,
contains information regarding the morphology of the source, and can be used
to separate stars from galaxies (see Koposov et al., 2017). The relevant
GaiaSource quantities used here (other than the sky coordinates RA, Dec), are:

  * •

$N_{\text{obs}}$, the number of times a source has crossed a CCD in the Gaia’s
focal plane;

  * •

$\text{F}_{G}$, the flux (electron per second) measured in the $G$ band
averaging over $N_{\text{obs}}$ single flux measurements;

  * •

$\sigma_{\text{F}_{G}}$, the standard deviation of the $N_{\text{obs}}$ flux
measurements;

  * •

$G$, the mean magnitude in the Gaia $G$ band (van Leeuwen et al., 2016)
calculated from $\text{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
$\mu_{\alpha}$ and $\mu_{\delta}$ 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

$\text{AMP}\equiv\log_{10}\left(\sqrt{N_{\text{obs}}}\frac{\sigma_{\text{F}_{G}}}{\text{F}_{G}}\right),$
(1)

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

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

####  2.2.1 The Gaia \+ 2MASS sample

Figure 1: Diagnostic diagrams for the selection criteria of our RRL sample.
Left-hand panel: distribution of stars in the AMP-$N_{\text{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 a randomly selected
subsample of bona fide RRLs from GCSS (red circles, 5% of the original sample)
and GS82 (orange squares, 35% of the original sample) catalogues (see Sec.
2.2.2). The horizontal-black and vertical-black lines show the selection cuts
used to obtain the final sample of RRLs (see Tab. 1), while the green arrows
indicate the regions used to obtain the final sample.

As mentioned, GDR1 reports photometric information only in the $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 et al. (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 datasets

In order to extract a reliable sample of RRL stars from the G2M catalogue, we
must apply ad-hoc selection criteria. To this aim, we used two samples of
bona-fide RRLs: the CSS (Drake et al., 2013; Drake et al., 2014) and the
Stripe 82 (S82, Sesar et al. 2010) catalogues. These samples allowed us to
identify the optimal selection criteria, analyse the completeness and the
contamination of the catalogue111The completeness indicates the fraction of
recovered true RRLs as a function of the apparent magnitude, while the
contamination is an estimate of the fraction of spurious objects (non RRLs)
that “pollute” our sample.  and estimate the RRL absolute magnitude in the $G$
band. The CSS contains about 22700 type-ab RRL stars distributed over a large
area of the sky (about 33,000 deg2 between $0^{\circ}<\alpha<360^{\circ}$ and
$-75^{\circ}<\delta<65^{\circ}$) and extended up to a distance of 70 kpc. The
completeness of this sample is constant (at $\sim 65\%$) 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\sim 15$, the
completeness is a strong function of the number of observations and thus
varies appreciably across the sky. SDSS’s Stripe 82 covers a 2.5∘-wide and
$100^{\circ}$-long patch of sky aligned with the celestial equator and
contains “only” 483 RRLs. However, the sample is very pure (with less than
$<1\%$ 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 $G$ band (see
Sec. 2.2.4), while the high quality of S82 sample is ideal to analyse the
completeness and contamination of our final sample (see Sec. 2.2.5 and Sec.
2.2.6). A cross-matching of the CSS and S82 catalogues with G2M using an
aperture of 1 arcsec led to the two samples CSS+G2M (GCSS hereafter) and
S82+G2M (GS82 hereafter).

####  2.2.3 RR Lyrae selection cuts

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

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

Figure 1 shows the distribution of the G2M stellar density (yellow-blue-purple
colour-maps), a randomly selected subsample of RRLs in GCSS (red points) and
in GS82 (orange squares) in the $N_{\text{obs}}$-AMP (left-hand panel), color-
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 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,
$G=17.1$ 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 $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 Sec.
2.2.5 and 2.2.6), we decided to put the bright magnitude limit at $G=10$.

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 $\text{AEN}>0.65$ as shown by the horizontal-black line. This
cut likely removes contaminant extragalactic objects since they typically have
$\text{AEN}\approx 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 proper motion, PM, greater than $50$ mas/yr. Given differences in
the light-curves and its sampling, the significance of AMP (Eq. 1) might
depend on the number of photometric measurements $N_{\text{obs}}$. With this
in mind, we impose $N_{\text{obs}}>30$: the focal plane of Gaia has an array
of 9$\times$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 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
$(\text{l},\text{b})=(280.47^{\circ},-32.89^{\circ})$ with an angular radius
of $9^{\circ}$ for the Large Magellanic Cloud (LMC) and the other centred on
$(\text{l},\text{b})=(302.80^{\circ},-44.30^{\circ})$ with a radius of
$7^{\circ}$ for the Small Magellanic Cloud (SMC). By inspecting the sky
distribution of the stars in our RRL sample, we noticed the presence of two
extended structures (S1 and S2, hereafter), that were not connected to any
known halo substructures, but are likely objects instead produced by Gaia
cross-match failures (see Sec. 2.2.5) that “survived” our selection cuts. We
decided to mask these sky regions as well by removing all the stars in the
following boxes $\text{l}=[167^{\circ},189^{\circ}]$
$\text{b}=[16^{\circ},22^{\circ}]$ for S1 and
$\text{l}=[160^{\circ},190^{\circ}]$ $\text{b}=[63^{\circ},73^{\circ}]$ 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 Sec. 2.2.5). In order to fully exploit the
all-sky capacity of our sample, we decided not to exclude a priori any portion
of the sky containing known substructures (e.g. unlike Deason et al. 2011
which masked out the Sagittarius stream). The analysis of the most significant
substructures found in this work can be found in Sec. 3 and 5.

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

Selection cuts
---
|b| [deg] | $>10$
$G$ [mag] | (10, 17.1)
AMP | (-0.7, -0.4)
$J-G$ | (-0.95, -0.4)
$N_{\text{obs}}$ | $>30$
AEN | $<0.65$
PM [mas/yr] | $<50$
|$\theta$| [deg] † | $>20$
Structure cuts
LMC | $\text{D}_{\text{LMC}}>9^{\circ}$
SMC | $\text{D}_{\text{SMC}}>7^{\circ}$
S1 | $\text{l}\notin[167^{\circ},189^{\circ}]\vee\text{b}\notin[16^{\circ},22^{\circ}]$
S2 | $\text{l}\notin[160^{\circ},190^{\circ}]\vee\text{b}\notin[63^{\circ},73^{\circ}]$
$N_{\text{stars}}$ | 21643 (13713†)
$f_{\text{V}}$ [%] | 58 (44†)
Table 1: Summary of the selection cuts used to obtain the final sample of RRLs
from the G2M catalogue. The description of parameters used in the sample
selection and the details on the cut substructures can be found in Sec. 2.1.
$\text{D}_{\text{LMC}}$ and $\text{D}_{\text{SMC}}$ are the sky angular
distances from the LMC and SMC, respectively. $\theta$ is the Galactocentric
latitude (defined in Eq. 7) and was estimated by assuming an RRL absolute
magnitude of $M_{\text{RLL}}=0.525$. The $\dagger$ refers to the subsample
used in the likelihood analysis in Sec. 4. The bottom part of the table gives
a summary of the whole sample. $N_{\text{stars}}$ is the number of stars in
the sample and $f_{\text{V}}$ is the fraction of the spherical volume of the
halo sampled by our stars between Galactocentric distance 0 kpc and 28 kpc.

####  2.2.4 Absolute magnitude and distance estimate

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

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

$\log\left(\frac{\text{D}_{\odot}}{\text{kpc}}\right)=\frac{G-M_{G}}{5}-2,$
(2)

where $M_{G}$ is the absolute magnitude in the Gaia $G$ band. To estimate
$M_{G}$ we used the RRLs in GCSS, as they have heliocentric distances
estimated from the period-luminosity relation. The resulting $M_{G}$
distribution is shown in Fig. 2: it has a well defined peak at $M_{G}\approx
0.5$ and a small dispersion. We fit this distribution with a Gaussian function
obtaining a mean of $0.525$ and a dispersion of $0.090$, comparable to the
uncertainties of the V-band absolute magnitude of RRLs (see e.g. Vivas & Zinn
2006). The above Gaussian function perfectly describes the data in the central
parts, but the distribution shows broader wings for $M_{G}<0.3$ and
$M_{G}>0.7$. 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 $M_{\text{RRL}}=0.525$. An analysis of the effect of this
approximation on the results can be found in Sec. 5.2.1.

####  2.2.5 Completeness

Figure 3: Completeness analysis of the sources in GaiaSource (Sec. 2.1) as a
function of Galactic latitude. Top panel: regions of the sky considered in
this analysis, each stripe has a constant Galactic longitude. Bottom panel:
ratio between the number of stars in the GaiaSource ($N_{\textit{Gaia}}$) and
stellar sources in the SDSS DR7 ($N_{\text{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).

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
$16<r<18$, where $r$ is the $r$-band magnitude. In this magnitude range the
SDR7 can be considered 100%
complete222http://classic.sdss.org/dr7/products/general/completeness.html. 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 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
$\text{b}>10^{\circ}$.

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 ($\sim$100 %) and pure ($\sim$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 ($-0.65$ red triangles, $-0.70$ blue circles, $-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% for $\text{AMP}>-0.65$ to about 30% for $\text{AMP}>-0.75,$ but it is
reasonably constant up to about $G=17.5$ then it abruptly decreases to 0 at
$G\approx 18,$ so we decided to conservatively cut our sample at $G=17.1$
(vertical black line in Fig. 4, see Sec. 2.2.3).

Figure 4: Completeness of our samples of RRLs as a function of distance from
the Sun ($\text{D}_{\odot}$) or $G$ magnitude. The conversion from
$\text{D}_{\odot}$ to $G$ has been obtained from Eq. 2 assuming a constant
absolute magnitude $M_{\text{RRL}}=0.525$. Different symbols indicate
completeness for samples obtained using different AMP cuts: red-triangles for
$\text{AMP}>-0.65$, blue circles for $\text{AMP}>-0.70$ and green diamonds for
$\text{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 Sec. 2.2.5). Figure 5: Contamination of the RRLs sample as a
function of the AMP cut. The contamination is estimated using the S82
catalogue (see Sec. 2.2.2) in the range of $G$ magnitudes $15<G<17$. The
error-bars have been calculated by propagating Poisson uncertainties from the
number of stars in AMP bins.

We also used the GCSS catalogue and found no significant variation of the
completeness as a function of the Galactic sky coordinates l and b, although
we found a mild trend for increasing $N_{\text{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 (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

$\text{Contamination}=\frac{N_{\rm S}-N_{\text{GS82}}}{N_{S}},$ (3)

where $N_{\rm S}$ and $N_{\text{GS82}}$ 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 $15<G<17$ for different lower limits of AMP.
For AMP threshold lager than $-0.7$ 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 $-0.7$ the
contamination fraction rapidly increases (about 25% for $\text{AMP}=-0.8$). We
expect that the contaminants of our sample are eclipsing-binaries in the
Galactic disc and possible instrumental artefacts. As shown and discussed in
detail in Belokurov et al. (2017) some regions of the Gaia all-sky map are
crossed by sharp strips with a large excess of objects. These strips are
similar in width to the field of view of Gaia and are regularly spaced on the
sky. Belokurov et al. (2017) propose that these are spurious features due to
failures in the cross-match procedure of Gaia, so that at some epochs the flux
of a star comes from a different object. The spurious measurement of the flux
increases $\sigma_{\text{F}_{G}}$ and moves the stars in the region of high
AMP “polluting” our samples. We found that for $\text{AMP}>-0.8$ the
contaminants are mostly related to cross-match failures. Unlike the disc
contaminants, the cross-match failures have a complicated and poorly
understood spatial distribution, so these structures are difficult to be taken
in account in the study of the properties of the Galactic halo. For this
reason, we decided to use $\text{AMP}=-0.7$ as the lower limit in our
selection cut (see Sec. 2.2.3).

##  3 Distribution of the RRLs in the inner stellar Halo

###  3.1 A first Gaia look at the stellar halo

Figure 6: Star-count all-sky maps in the Galactic coordinates (l, b) for the
RRLs in our sample (Sec. 2.2.3) in the magnitude intervals $10<G<15.5$ (upper
panel) and $15.5<G<17.1$ (lower panel). The “holes” between
$\text{l}=-50^{\circ}$ and $\text{l}=-100^{\circ}$ are due to the mask used to
eliminate the contribution of the Large and Small Magellanic Clouds (see Sec.
2.2.3).

Fig. 6 shows the distribution of the RRLs of our final sample as a function of
the Galactic coordinates (l, b). We split the stars in the sample into $G$
magnitude intervals: $10<G<15.5$, corresponding to heliocentric distances
between $\sim$1 kpc and 10 kpc (Eq. 2), and $15.5<G<17.1$, i.e. with
heliocentric distances between 10 kpc 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.

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
($\text{l}\approx 0$) 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
$\text{l}=-50^{\circ}$. In both magnitude intervals an over-dense band of
stars at low Galactic latitude ($\text{|b|}<15^{\circ}-20^{\circ}$) can be
seen running all around the Galaxy. A discussion of the nature of these
structures can be found in Sec. 5.3.

###  3.2 Setting the frame of reference

We set a left-handed Cartesian reference frame
$(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}})$ centred in the
Galactic centre and such that the Galactic disc lies in the plane
$(\text{X}_{\text{g}},\text{Y}_{\text{g}})$, the Sun lies on the positive
X-axis and the Sun rotation velocity is $\dot{\text{Y}}_{\text{g}}>0$. 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 $\text{D}_{\odot}$ from the
Sun are

$\left\\{\begin{array}[]{lll}\text{X}_{\text{g}}=\text{X}_{\text{g},\odot}-\text{D}_{\odot}\cos\text{b}\cos\text{l},\\\
\text{Y}_{\text{g}}=\text{D}_{\odot}\cos\text{b}\sin\text{l},\\\
\text{Z}_{\text{g}}=\text{D}_{\odot}\sin\text{b},\end{array}\right.$ (4)

where $\text{X}_{\text{g},\odot}$ is the distance of the Sun from the Galactic
centre. In this work we assume $\text{X}_{\text{g},\odot}=8\ \text{kpc}$ (Bovy
et al., 2012), but the main results of our work are unchanged for other values
of $\text{X}_{\text{g},\odot}$, 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
$(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}})$ we define the
distance from the Galactic centre

$\text{D}_{\text{g}}=\sqrt{\text{X}^{2}_{\text{g}}+\text{Y}^{2}_{\text{g}}+\text{Z}^{2}_{\text{g}}},$
(5)

the Galactocentric cylindrical radius

$\text{R}=\sqrt{\text{X}^{2}_{\text{g}}+\text{Y}^{2}_{\text{g}}},$ (6)

the Galactocentric latitude

$\theta=\arctan\frac{\text{Z}_{\text{g}}}{\text{R}},$ (7)

and the Galactocentric longitude

$\phi=\arctan{\frac{\text{Y}_{\text{g}}}{\text{X}_{\text{g}}}}.$ (8)

Assuming that the stellar halo is stratified on concentric ellipsoids, it is
useful to introduce another Cartesian frame $(\text{X},\text{Y},\text{Z})$,
aligned with the ellipsoid principal axes, and to define the elliptical radius

$\text{r}_{\text{e}}=\sqrt{\text{X}^{2}+\text{Y}^{2}\text{p}^{-2}+\text{Z}^{2}\text{q}^{-2}},$
(9)

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

###  3.3 Density distribution of the halo RRLs

Given the all-sky RRL sample illustrated in Fig. 6, we can compute their
volume number density distribution $\rho$. In particular we define the number
density of halo RRLs in a cell $\Delta\vec{\epsilon}$ centred at
$\vec{\epsilon}$, where $\vec{\epsilon}$ is a coordinate vector, such as for
instance $(\text{R},\text{Z}_{\text{g}})$, as

$\rho(\vec{\epsilon},\Delta\vec{\epsilon})\approx
N(\vec{\epsilon},\Delta\vec{\epsilon})f^{-1}_{\text{V}}(\vec{\epsilon},\Delta\vec{\epsilon})\text{V}^{-1}(\vec{\epsilon},\Delta\vec{\epsilon}),$
(10)

where $N(\vec{\epsilon},\Delta\vec{\epsilon})$ is the number of stars observed
in the cell, V is the volume of the cell and $f_{\text{V}}$ is the fraction of
this volume accessible to our analysis, which depends on the sky coverage, on
the selection cuts and the mask that we applied (see Sec. 2.2.3), and on the
off-set between the Sun and the Galactic centre. We estimate $f_{\text{V}}$
numerically as follows. First, we define a Cartesian Galactic grid sampling
each of the axes
$(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}})$ with 300
points linearly spaced between $-\text{D}_{\text{g,max}}$ and
$\text{D}_{\text{g,max}}$ ($\text{D}_{\text{g,max}}\simeq 28$ kpc is the
maximum value of $\text{D}_{\text{g}}$, given our magnitude limit at
$G=17.1$), 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 kpc3. 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, $\alpha$, $\delta$). A
secondary (non-uniform) grid is built by applying the selection cuts in Tab. 1
to the points on the primary grid, so we end up with a grid sampling the
complete Galactic volume and a grid sampling the portion of the volume
accessible to our analysis. Given a cell in a certain coordinates space (e.g.
R, $\text{Z}_{\text{g}}$), we define $f_{\text{V}}$ 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 ($N$), the second the total volume of each bin/cell (V)
and the third the fraction of this volume that we are actually sampling
($f_{\text{V}}$). These three quantities are inserted in Eq. 10 to estimate
the stellar number density. We tested this method with both simple analytic
distributions and mock catalogues (see Sec. 4.3.5).

####  3.3.1 Meridional plane

Figure 7: Number density in the Galactic $\text{R}-\text{Z}_{\text{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 (Eq. 10). The black
contours are plotted with spacing $\Delta\log\rho=0.4$ from $\log\rho=-1.2$ to
$\log\rho=1.6$, where $\rho$ is in units of $\text{kpc}^{-3}$; the white
dashed lines represent elliptical contours with an axial ratio $\text{q}=0.6$.

Fig. 7 shows the number density of the RRLs of our sample in the Galactic
meridional plane $\text{R}-\text{Z}_{\text{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 $\text{Z}_{\text{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 Sec. 5.3.
The density at higher $\text{Z}_{\text{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 the $\text{q}=0.6$
elliptical contours (overplotted in Fig. 7) out to $\text{R}\approx 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 over-dense region at $\text{Z}_{\text{g}}>10$ kpc, which
does not seem to have a counterpart below the Galactic plane (see Sec. 3.3.4
for further details).

####  3.3.2 Density profile

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

Fig. 7 demonstrates that the RRLs in the inner halo are consistent with being
stratified on spheroids (except in the region close to the Galactic plane),
thus we estimate the 1D RRL density profile by counting the stars in
spheroidal ($\text{p}=1$, $\text{q}\neq 1$) 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 $\text{q}=1.0$, $\text{q}=0.6$ and
$\text{q}=0.3$. 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
$\alpha=-2.7$ (see Fig. 8), but, based on this simple analysis, all values of
$\alpha$ in the range $-3<\alpha<-2.5$ are consistent with the data.

####  3.3.3 Halo flattening

Figure 9: Number density of RRLs in the elliptical radius
($\text{r}_{\text{e}}$) - Galactocentric latitude ($\theta$) space for
$\text{q}=0.75$ (left-hand panel), $\text{q}=0.55$ (middle panel) and
$\text{q}=0.35$ (right-hand panel), assuming $\text{p}=1$. The density is
normalized to the maximum value of each panel. The contours show the
normalised density levels $(0.0002,0.00055,0.001,0.002,0.004,0.02,0.04,0.2)$,
while the dashed lines indicate $|\theta|=20^{\circ}$.

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 $\text{q}\approx 0.6$. However, in order to more
rigorously study the halo flattening, we estimate the density in the
$\text{r}_{\text{e}}-\theta$ plane (see equations 4 and 9). In practice, for a
given value of $q$, we scan the density as a function of $\text{r}_{\text{e}}$
at fixed $\theta$. If the RRLs are truly stratified on similar spheroids with
axial ratio $q$, the density is independent of $\theta$, so the isodensity
contours in the $\text{r}_{\text{e}}-\theta$ plane are vertical stripes. If
the assumed value of q is smaller than the true value
$\text{q}_{\text{true}}$, $\text{r}_{\text{e}}$ is underestimated, the
estimated density is a monotonic decreasing function of $\theta$ and the
isodensity contours in the $\text{r}_{\text{e}}-\theta$ plane are bent in the
direction of $|\theta|$ increasing with $\text{r}_{\text{e}}$ (provided that
the density is a decreasing function of $\text{r}_{\text{e}}$). The isodenisty
contours are bent in the opposite direction, if one assumes
$\text{q}>\text{q}_{\text{true}}$. The shape of the iso-density contours in
the $\text{r}_{\text{e}}-\theta$ 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 $\text{r}_{\text{e}}-\theta$ plane of the RRLs
in our sample are shown in Fig. 9 for $\text{q}=0.75$, $\text{q}=0.55$ and
$\text{q}=0.35$, assuming $\text{p}=1$. Below $|\theta|=20^{\circ}$ (indicated
by the white-dashed lines) the contours are nearly horizontal, because the
density is dominated by the highly flattened discy component (see Fig. 7) for
which q is over-estimated in all the panels. At higher latitudes
($|\theta|>20^{\circ}$) the contours give a direct indication on the
flattening of the halo: in the right-hand panel of Fig. 9 ($\text{q}=0.35$)
the iso-density contours are significantly inclined in a way that implies
$\text{q}_{\text{true}}>0.35$. For $\text{r}_{\text{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 ($\text{q}=0.55$), but
beyond 20 kpc the contours start to bend so that
$\text{q}_{\text{true}}>0.55$. The last two iso-density contours in the left-
hand panel ($\text{q}\approx 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 and Das
et al. 2016). In Sec. 4.4 we perform a comprehensive model fitting analysis of
the RRL dataset and compare our results to those found in the literature.

####  3.3.4 Vertical asymmetries

Figure 10: Number counts of stars in the
$\text{X}_{\text{g}}-\text{Y}_{\text{g}}$ plane for different
$\text{Z}_{\text{g}}$-slabs: left-hand panels with
$0<|\text{Z}_{\text{g}}/\text{kpc}|<5$, middle panels with
$5<|\text{Z}_{\text{g}}/\text{kpc}|<10$ and right-hand panels with
$|\text{Z}_{\text{g}}/\text{kpc}|>10$. The top and bottom panels show the
layers above ($\text{Z}_{\text{g}}>0$) and below ($\text{Z}_{\text{g}}<0$) the
Galactic plane, respectively. The colour maps are normalised to the maximum
number in each column of plots. $N_{\text{tot}}$ is the total number of stars
found on each $\text{Z}_{\text{g}}$-slab. The black dots show the position of
the Sun (at $\text{Z}_{\text{g}}=0$), while the crosses indicate the maximum
of the star counts. The maps have been smoothed with a spline kernel.

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

Fig. 10 shows star-count maps in the $\text{X}_{\text{g}}-\text{Y}_{\text{g}}$
plane integrated along three different $\text{Z}_{\text{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%) 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 $\text{Z}_{\text{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%) 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
$\text{Z}_{\text{g}}$ (right-hand panels) show a significant difference in the
total number of stars (about 20%) that can not be explained by Poisson
fluctuations only. Moreover the excess of stars above the Galactic plane is
strongly clustered in the regions between $\text{X}_{\text{g}}\approx(5,15)$
kpc and $\text{Y}_{\text{g}}\approx(0,10)$ 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 $|\text{Z}_{\text{g}}|>10$ kpc, significantly less than the
value of 20% obtained here. Therefore, the structure seen at high positive
$\text{Z}_{\text{g}}$ appears to be genuine, and is most likely related to the
“Virgo over-density” (Jurić et al., 2008; Vivas et al., 2016). Across all
$\text{Z}_{\text{g}}>0$ slabs, portions of the Virgo Cloud are visible as an
excess of stars at positive $\text{X}_{\text{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 over-
density at negative $\text{X}_{\text{g}}$ (and $\text{Z}_{\text{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

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 $\theta$ (see Eq. 7).

In particular, at $|\theta|=20^{\circ}$ 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 $|\theta|>20^{\circ}$: 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 $\theta$ obtained from Eq. 2 and Eq. 7 assuming the same
value of the absolute magnitude, $M=M_{\text{RRL}}=0.525$, 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
$\text{Z}_{\text{g}}$. We found that the results obtained for the sample with
$|\theta|>20^{\circ}$ (see Sec. 4.4) are qualitatively similar to the results
obtained using a sample with $|\text{b}|>30^{\circ}$ (the lowest Galactic
latitude used in most of the previous works, e.g. Deason et al. 2011) or a
sample with $|\text{Z}_{\text{g}}|>6$ kpc (Fermani & Schönrich 2013 uses
$|\text{Z}_{\text{g}}|>4$ kpc to cut disc stars). The cuts on b and
$\text{Z}_{\text{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
$|\theta|=20^{\circ}$.

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

###  4.2 Halo models

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

####  4.2.1 Density profiles

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

$\rho(\text{r}_{\text{e}})\propto\left(\frac{\text{r}_{\text{e}}}{\text{r}_{\text{eb}}}\right)^{-\alpha_{\text{inn}}}\left(1+\frac{\text{r}_{\text{e}}}{\text{r}_{\text{eb}}}\right)^{-(\alpha_{\text{out}}-\alpha_{\text{inn}})},$
(11)

where $\alpha_{\text{inn}}$ and $\alpha_{\text{out}}$ indicate the inner and
outer power-law slopes and $\text{r}_{\text{eb}}$ is the scale length. The
number density of the SPL is given by Eq. 11 when
$\alpha_{\text{inn}}=\alpha_{\text{out}}$, while the number density of the CPL
has $\alpha_{\text{inn}}=0$, in which case $\text{r}_{\text{eb}}$ represents
the length of the inner core. The BPL number density profile is given by a
piecewise function:

$\rho(\text{r}_{\text{e}})\propto\begin{cases}\text{r}_{\text{e}}^{-\alpha_{\text{inn}}}&\text{r}_{\text{e}}\leq\text{r}_{\text{eb}}\\\
\text{r}_{\text{eb}}^{\alpha_{\text{out}}-\alpha_{\text{inn}}}\text{r}_{\text{e}}^{-\alpha_{\text{out}}}&\text{r}_{\text{e}}\geq\text{r}_{\text{eb}}\end{cases}.$
(12)

The EIN profile (Einasto, 1965) is given by

$\rho(\text{r}_{\text{e}})\propto\text{exp}\left[-d_{\text{n}}\left(\left(\frac{\text{r}_{\text{e}}}{\text{r}_{\text{eb}}}\right)^{\frac{1}{\text{n}}}-1\right)\right],$
(13)

where $d_{\text{n}}=3\text{n}-0.3333+0.0079/\text{n}$ for $\text{n}\geq 0.5$
(Graham et al., 2006). The steepness of the EIN profile,
$\alpha_{\text{EIN}}$, changes continuously as a function of
$\text{r}_{\text{e}}$ tuned by the parameter n,

$\alpha_{\text{EIN}}=-\frac{d_{\text{n}}}{\text{n}}\left(\frac{\text{r}_{\text{e}}}{\text{r}_{\text{eb}}}\right)^{\frac{1}{\text{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 power-law profiles with $\alpha_{\text{inn}}<-3$ or with
$\alpha_{\text{out}}>-3$ 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
$\alpha_{\text{inn}}<-3$ (see Sec. 4.4), therefore there should be a physical
radius, outside our radial range, beyond which the profile becomes, either
abruptly (e.g. BPL or an exponential truncation) or gently (e.g. DPL),
steeper.

####  4.2.2 Iso-density ellipsoidal surfaces

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

  * •

spherical (SH): we set $\text{p}=1$ and $\text{q}=1$ in Eq. 9, so that
$\text{r}_{\text{e}}$ is just the spherical radius $\text{D}_{\text{g}}$ (see
definition 5) in the Galactic frame of reference;

  * •

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

  * •

disc-plane axisymmetric (DP): we set $\text{q}=1$, the axis of symmetry is
within the Galactic plane making an (anticlockwise) angle $\gamma$ 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 $\gamma$ 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

$\text{q}(\text{r}_{\text{e}})=\text{q}_{\infty}-(\text{q}_{\infty}-\text{q}_{0})\text{exp}\left[1-\frac{\sqrt{\text{r}_{\text{e}}^{2}+\text{r}_{\text{eq}}^{2}}}{\text{r}_{\text{eq}}}\right],$
(15)

so q varies from $\text{q}_{0}$ at the centre to the asymptotic value
$\text{q}_{\infty}$ at large radii and the variation is tuned by the
exponential scale length $\text{r}_{\text{eq}}$;

  * •

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 $R(\gamma,\beta,\eta)$ following
a ZYX formalism so that $\gamma$ is the rotation angle around the original
Z-axis, $\beta$ the one around the new Y-axis and finally $\eta$ is the
rotation angle around the final X-axis, all the rotation are defined in the
anticlockwise direction;

  * •

offset (off), in this model the elliptical radius of each star is estimated
with respect to a point offset by ($\text{X}_{\text{off}}$,
$\text{Y}_{\text{off}}$, $\text{Z}_{\text{off}}$) 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 ($\text{D}_{\text{g}}$). We decided to
maintain the dependence on $\text{r}_{\text{e}}$ given that this approach is
self-consistent with our assumption that the RRLs are stratified on
ellipsoidal surfaces. If follows that for our models with a varying q, the
elliptical radius $\text{r}_{\text{e}}$ of a star with coordinates
$(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}})$ is the root of

$\text{r}^{2}_{\text{e}}-\text{X}^{2}_{\text{g}}-\text{Y}^{2}_{\text{g}}\text{p}^{-2}-\text{Z}^{2}_{\text{g}}\text{q}(\text{r}_{\text{e}})^{-2}=0,$
(16)

where $\text{q}(\text{r}_{\text{e}})$ is defined in Eq. 15. In our fitting
code we solve this equation numerically with a Newton-Raphson root finder.

####  4.2.3 Complete halo model

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

###  4.3 Comparing models with observations

####  4.3.1 Density of stars in the observed volume

Given a certain halo model (Sec. 4.2), the normalised number density of stars
in an infinitesimal volume
$\text{dV}_{\text{g}}=\text{d}\text{X}_{\text{g}}\text{d}\text{Y}_{\text{g}}\text{d}\text{Z}_{\text{g}}$
is given by

$\tilde{\rho}(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}}|\vec{\mu})\propto\frac{\text{d}N}{\text{dV}_{\text{g}}},$
(17)

where $\vec{\mu}$ is the vector containing all the model parameters (e.g.
$\vec{\mu}=(\alpha_{\text{inn}},\text{q})$ for an SPL+disc-normal axisymmetric
model, see Sec. 4.2) and $\tilde{\rho}$ is defined such that

$\int\tilde{\rho}(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}}|\vec{\mu})\text{dV}=1.$
(18)

It is useful to define the normalised star number density

$\tilde{\nu}(m,\text{l},\text{b}|M,\vec{\mu})=|\text{J}|\tilde{\rho}(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}}|m,\text{l},\text{b},M,\vec{\mu})$
(19)

within the infinitesimal projected (onto the sky) volume element
$\text{d}S=\text{d}m\text{dl}\text{db}$ centred at ($m$,l,b), where $m$ is the
observed magnitude, l and b are the Galactic longitude and latitude, and

$|\text{J}|=\frac{\ln 10}{5}\text{D}^{3}_{\odot}(m,M)\cos\text{b}$ (20)

is the determinant of the Jacobian matrix
$J=\left[\frac{\text{dV}_{\text{g}}}{\text{d}S}\right]$ (Appendix A). The
number density $\tilde{\nu}$ depends also on the additional parameter $M$,
which is the absolute magnitude needed to pass from the observable variables
$(m,\text{l},\text{b})$ to the Cartesian variables (see Eqs. 2 and 4).
Substituting Eq. 20 in Eq. 19 we obtain

$\tilde{\nu}(m,\text{l},\text{b}|M,\vec{\mu})=\frac{\ln
10}{5}\tilde{\rho}(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}}|m,\text{l},\text{b},M,\vec{\mu})\text{D}^{3}_{\odot}(m,M)\cos\text{b}.$
(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 ($m$, l, b) given the absolute magnitude $M$
and a halo model with parameters $\vec{\mu}$

$\lambda(m,\text{l},\text{b},M|\vec{\mu})=A\tilde{\nu}(m,\text{l},\text{b}|M,\vec{\mu})W(m,\text{l},\text{b}|M_{\text{RRL}})P(M).$
(22)

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

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

In this work, we decided to set the absolute magnitude of the RRLs to a single
value, $M_{\text{RRL}}=0.525$ (Sec. 2.2.4), so we are assuming that $P(M)$ in
Eq. 22 is a Dirac delta. As a consequence, we can marginalize Eq. 22 over $M$
to obtain

$\lambda(m,\text{l},\text{b}|M_{\text{RRL}},\vec{\mu})=A\tilde{\nu}(m,\text{l},\text{b}|M_{\text{RRL}},\vec{\mu})W(m,\text{l},\text{b}|M_{\text{RRL}}),$
(23)

and we define the normalised rate function as

$\tilde{\lambda}(m,\text{l},\text{b}|M_{\text{RRL}},\vec{\mu})=\frac{\lambda(m,\text{l},\text{b}|M_{\text{RRL}},\vec{\mu})}{\int\lambda(m,\text{l},\text{b}|M_{\text{RRL}},\vec{\mu})\text{d}m\
\text{dl}\ \text{db}}.$ (24)

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

####  4.3.3 The total likelihood

Consider a sample of stars $D$ with coordinates ($m$, l, b) and a halo model
with parameters $\vec{\mu}$. Given the pdf of the stellar distribution
$\tilde{\lambda}(\text{m},\text{l},\text{b})$ (Eq. 24) the logarithmic
likelihood is

$\ln\mathcal{L}=\sum^{N_{\text{s}}}_{i=1}\ln\tilde{\lambda}(m_{i},\text{l}_{i},\text{b}_{i}|\vec{\mu}),$
(25)

where $N_{\text{s}}$ is the number of stars in our sample. Plugging Eq. 21 and
Eq. 23 into Eq. 24 we can write the logarithmic likelihood as

$\ln\mathcal{L}=\sum^{N_{\text{s}}}_{i=1}\ln\frac{\tilde{\rho}(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}}|m_{i},\text{l}_{i},\text{b}_{i},M_{\text{RRL}},\vec{\mu})\text{D}^{3}_{\odot}(m_{i},M_{\text{RRL}})\cos\text{b}}{\text{V}_{c}(\vec{\mu})}$
(26)

where

$\text{V}_{c}(\vec{\mu})=\int^{90^{\circ}}_{-90^{\circ}}\cos\text{b}\
\text{db}\int^{360^{\circ}}_{0^{\circ}}\text{dl}\int^{G_{\text{max}}}_{G_{\text{min}}}\tilde{\rho}\text{D}^{3}_{\odot}W\text{d}m$
(27)

is the normalisation integral. Notice that the numerator of Eq. 26 is
evaluated only in regions of the sky where $W=1$ (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 $\vec{\mu}$ of an halo model given the data is

$\ln P\left(\vec{\mu}|D\right)=\ln
P\left(D|\vec{\mu}\right)+\ln\Pi\left(\vec{\mu}\right)=\ln\mathcal{L}+\ln\Pi\left(\vec{\mu}\right),$
(28)

where $P\left(D|\vec{\mu}\right)=\mathcal{L}$ is the probability of the data,
$D$, given the parameters $\vec{\mu}$ (See sec. 4.3.3) and
$\Pi\left(\vec{\mu}\right)$ represents the prior probability of $\vec{\mu}$
(Tab. 2). In Eq. 28 we omit the Bayesian evidence term $P(D)$ 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

$\begin{split}\ln&P\left(\vec{\mu}|D\right)=\\\
&\sum^{N_{\text{s}}}_{i=1}\ln\frac{\tilde{\rho}(\text{X}_{\text{g}},\text{Y}_{\text{g}},\text{Z}_{\text{g}}|m_{i},\text{l}_{i},\text{b}_{i},M_{\text{RRL}},\vec{\mu})\text{D}^{3}_{\odot}(m_{i},M_{\text{RRL}})\cos\text{b}}{\text{V}_{c}(\vec{\mu})}\\\
&+\ln\Pi\left(\vec{\mu}\right).\end{split}$ (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 $\alpha_{\text{inn}}$ and $\alpha_{\text{out}}$ approximate
uniform infinite priors.

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

We analysed the posterior probability of the parameters, $\vec{\mu}$, 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
emcee333https://github.com/dfm/emcee (Foreman-Mackey et al., 2013). The
technical details of our approach can be found in Appendix B.

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

####  4.3.5 Tests on mock catalogues

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

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

Analysing the mock catalogues, we can estimate what properties of the halo we
are able to measure and constrain using the stars in our catalogue. We found
that we are able to detect a core in the density profile (CPL density model)
all the way down to 100 pc. However for a very small core size most of the
information comes from a relative small region where the transition between
the inner core and the outer density profile takes place. The detection of a
change of slope depends on the halo flattening: for $\text{p}=1$ and
$\text{q}=1$ we can detect a break in the density profile within about
$\text{r}_{\text{eb}}=28$ kpc (Eq. 21). For more flattened models (q around
$0.5/0.6$) this range extends up to about $\text{r}_{\text{eb}}=35$ kpc.
However, in this latter case most of the information for regions where
$\text{r}_{\text{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 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 | $\Delta\ln(\mathcal{L}_{\text{max}})$ | $\Delta\text{BIC}$
SPL | SH | $\alpha_{\text{inn}}=2.61\pm 0.01\ (2.61)\dagger$ | -868 | 1688
SPL | DP |  |  $\alpha_{\text{inn}}=2.64\pm 0.01\ (2.64)\dagger$,
---
$\text{p}=1.59\pm 0.03\ (1.59)$, $\gamma=-22.0\pm 1.7\ \text{deg}\ (-22.3\
\text{deg})\ddagger$
-405 | 779
SPL | DN |  |  $\alpha_{\text{inn}}=2.71\pm 0.01\ (2.71)\dagger$,
---
$\text{q}=0.58\pm 0.01\ (0.58)$
-81 | 124
BPL | DN |  |  $\alpha_{\text{inn}}=2.70\pm 0.01\ (2.70)$, $\alpha_{\text{out}}=2.90\pm 0.11\ (2.87)$,
---
$\text{r}_{\text{eb}}=22.3^{+3.3}_{-2.8}\ (22.2)$ kpc, $\text{q}=0.58\pm 0.01\
(0.58)$
-79 | 141
DPL | DN |  |  $\alpha_{\text{inn}}=2.70\pm 0.03\ (2.71)$, $\alpha_{\text{out}}=2.74^{+0.04}_{-0.02}\ (2.72)$,
---
$\text{r}_{\text{eb}}=21.1^{+5.7}_{-9.9}\ (23.1)$ kpc, $\text{q}=0.58\pm 0.01\
(0.58)$
-81 | 143
CPL | DN |  |  $\alpha_{\text{out}}=2.72\pm 0.02\ (2.71)$, $\text{r}_{\text{eb}}=0.03^{+0.03}_{-0.01}\ (0.03)$ kpc,
---
$\text{q}=0.58\pm 0.01\ (0.58)$
-81 | 134
EIN | DN |  |  $\text{n}=37.5^{+3.5}_{-4.6}\ (40)$, $\text{r}_{\text{eb}}=391.1^{+84.6}_{-169.8}\ (474)$ kpc,
---
$\text{q}=0.58\pm 0.01\ (0.58)$
-87 | 145
SPL | $\text{DN}^{\it qv}$ |  |  $\alpha_{\text{inn}}=2.93\pm 0.05\ (2.93)\dagger$,
---
$\text{q}_{0}=0.52\pm 0.02\ (0.52)$, $\text{q}_{\infty}=0.74\pm 0.05\ (0.75)$,
$\text{r}_{\text{eq}}=14.8\pm 1.9\ (15.0)$ kpc
-66 | 113
SPL | TR |  |  $\alpha_{\text{inn}}=2.71\pm 0.01\ (2.71)\dagger$,
---
$\text{p}=1.27\pm 0.03\ (1.27)$, $\gamma=-21.1\pm 2.6\ \text{deg}\ (-21.1\
\text{deg})\ddagger$
$\text{q}=0.65\pm 0.01\ (0.65)\ddagger$
-22 | 23
SPL | $\text{TR}^{\it qv}$ |  |  $\alpha_{\text{inn}}=2.96\pm 0.05\ (2.96)\dagger$,
---
$\text{p}=1.27\pm 0.03\ (1.27)$, $\gamma=-21.3\pm 2.6\ \text{deg}\ (-21.3\
\text{deg})\ddagger$
$\text{q}_{0}=0.57\pm 0.02\ (0.57)$, $\text{q}_{\infty}=0.84\pm 0.06\ (0.84)$,
$\text{r}_{\text{eq}}=12.2^{+2.4}_{-1.8}\ (12.2)$ kpc
0 | 0
Table 3: Summary of properties and results for a sample of families of halo
models. For each family we report the assumed density law (Sec. 4.2.1) and
geometry of the iso-density surfaces (Sec. 4.2.2). See Tab. 2 for a reference
on the model labels. For each fitted parameter we report the median and the
uncertainties estimated as the 16th and 84th percentiles of the posterior
distribution, the values in parentheses indicate the value for which we obtain
the maximum value of the likelihood $\mathcal{L}_{\text{max}}$ (Eq. 26). The
last two columns indicate the logarithmic likelihood and BIC differences
between the best model of the family and the best of of all the presented
models. $\dagger$ In the SPL models $\alpha_{\text{inn}}=\alpha_{\text{out}}$,
see Sec. 4.2.1. $\ddagger$ The (anticlockwise) angle $\gamma$ indicates the
tilt of the X and Y axes of symmetry of the halo with respect to the Galactic
X and Y axes, see Sec. 4.2.2.

####  4.4.1 Model comparison

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

$\text{BIC}=\dim(\vec{\mu})\ln\left(N_{\text{S}}\right)-2\ln\left(\mathcal{L}_{\text{max}}\right),$
(30)

where $N_{\text{S}}$ is the number of objects in the sample and
$\mathcal{L}_{\text{max}}$ is the maximum value of the likelihood (Eq. 26).

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

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

####  4.4.2 RRLs density law

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

We tested different density profiles assuming a DN geometrical model: the SPL
(Eq. 11 with $\alpha_{\text{inn}}=\alpha_{\text{out}}$), the DPL (Eq. 11), the
BPL (Eq. 12), the CPL (Eq. 11 with $\alpha_{\text{inn}}=0$) and the EIN
profile (Eq. 13).

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

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

Finally, for the EIN the fit favours very large values both for n ($\approx
38$) and for $\text{r}_{\text{eb}}$ ($\approx 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 Tab. 4). In principle, the over-density of stars described in Sec. 3.3.4
might make it difficult to detect the break, so we repeated the maximum
likelihood analysis using the RRLs only above and below the Galactic plane.
The results are shown in Fig. 12 by the solid ($\text{b}<-10^{\circ}$) and
dashed ($\text{b}>10^{\circ}$) 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 $\text{r}_{\text{e}}\approx 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 sub-
samples the minimum BIC is still obtained for the SPL model. Therefore, the
change of slope is not highly significant and the fact that the DPL and the
BPL show a different behaviour could indicate that the change of slope could
be due to some local artefact (such as a local decrease of the completeness).

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

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

####  4.4.3 Iso-density surfaces

As shown by the BIC comparison in Tab. 3, the SH and DP models give
significantly worse fits with respect to the DN and TR models. This result is
confirmed by a comparison of the properties of the RRLs in our sample with the
ones expected for the best-fit models. For example, Fig. 13 shows the fraction
of RRLs in our sample in different observed magnitude bins (black points)
compared to the predicted distribution for different halo models (assuming a
SPL, see Tab. 3). The SH and the SP models predict too low a fraction of stars
at lower $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, $\text{q}=1$ and
$\text{p}=1$), DN (orange, $\text{q}=0.58$ and $\text{p}=1$), DP (green,
$\text{q}=1$ and $\text{p}=1.59$) and TR (red, $\text{q}=0.65$ and
$\text{p}=1.27$). In all cases we assumed a SPL for the density profile of the
RRLs (with parameters given in Tab. 3).

As a further refinement we tested the $\text{DN}^{\it qv}$ and $\text{TR}^{\it
qv}$ 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 ($\text{q}\approx 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
($\text{r}_{\text{e}}\approx 30$ kpc). These results are in agreement with
what we have seen directly in the distribution of stars in the
$\text{r}_{\text{e}}$-$\theta$ plane in Sec. 3.3.3 and Fig. 9. Essentially the
same result is obtained by considering the stars only above or only below the
Galactic plane. Fig. 15 shows the comparison between the magnitude
distribution of the stars in our sample and the one expected for the varying q
model.

In the (non-tilted) TR and the DP models the elongated axis makes an
(anticlockwise) angle $\gamma\approx-20^{\circ}$ 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 $\Phi_{\text{bar}}$444$\Phi_{\text{bar}}$ is defined as the
anticlockwise angle between the direction of the bar and the Galactic X-axis.
is uncertain (Antoja et al., 2014), depending on the works it ranges from
about $-10^{\circ}$ (e.g. Robin et al. 2012) to about $-45^{\circ}$ (e.g.
Benjamin et al. 2005). Given these large uncertainties, the angle $\gamma$ can
be considered compatible with $\Phi_{\text{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.

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

####  4.4.4 Tilt and offset

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

For the offset model, the best-fit parameters are

  * •

$\text{X}_{\text{off}}=0.39\pm 0.05$ kpc,

  * •

$\text{Y}_{\text{off}}=-0.17\pm 0.06$ kpc,

  * •

$\text{Z}_{\text{off}}=-0.01\pm 0.01$ kpc,

and the size of the displacement vector with respect to the Galactic centre is
$\text{D}_{\text{off}}=0.43\pm 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 over-fitting of some local substructure at
small elliptical radii. Moreover, both the posterior distribution of
$\text{Y}_{\text{off}}$ and $\text{Z}_{\text{off}}$ are compatible (within
$3-\sigma$ and $1-\sigma$, 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 $\text{X}_{\odot}=7.61\pm 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-fit results for the (anticlockwise)
angles describing the orientation of the halo (see Sec. 4.2.2) are

  * •

$\gamma=-20.1^{\circ}\pm 2.7^{\circ}$ (rotation around the Z-axis),

  * •

$\beta=5.7^{\circ}\pm 0.8^{\circ}$ (rotation around the new Y-axis),

  * •

$\eta=3.13^{\circ}\pm 0.5^{\circ}$ (rotation around the final X-axis).

The angle $\gamma$ represents the tilt of the elongated axis of the halo with
respect to the Galactic Y-axis and it is compatible with the orientation found
for the non-tilted DN and TR disc-plane models (Tab. 3). The angles $\beta$
and $\eta$ 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 ($\Delta\text{BIC}\approx-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^{\it qv}$, $TR^{\it tl}$, $TR^{\it off}$ and $TR^{\it 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 Sec. 5.5).

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

##  5 Discussion

###  5.1 Best halo model

Figure 16: Same as in Fig. 9 but for mock stellar halos. Top: density maps for
the best-fit SPL-DN model. Bottom: density maps for the best-fit SPL-DNqv (see
Tab. 3). The maps have been obtained by averaging over 1000 halo mock
realisations. Note that for a halo with a fixed flattening, for the (almost)
correct value of q (middle panel) the iso-density contours remain vertical
across the whole range of elliptical radii. This is at odds with the RRL
distribution in the Milky Way. As the middle panel of Fig. 9 demonstrates, the
contours indeed start vertical for $\text{r}_{\text{e}}<20$ kpc, but beyond
that, there is a noticeable bending away from the centre. Our best-fit model
with varying q displays exactly the same behaviour (see middle panel of bottom
row). Figure 17: Top panels: all-sky star number density maps in Galactic
coordinates. Bottom panels: sky-maps density plots in equatorial coordinates.
The left-hand panels show the data, while the middle panels show our best
model (SPL-TRqv, see Sec. 5.1 and Tab. 3), the right-hand panels show the
relative data-model residuals. The model maps have been obtained integrating
the star density distribution (Eq. 21) through the magnitude interval
$10<G<17.1$ assuming a constant absolute magnitude $M=M_{\text{RRL}}=0.525$.
Figure 18: Number density in the Galactic $\text{R}-\text{Z}_{\text{g}}$
plane. Left-hand panel: density map for the RRLs in our clean sample (Sec.
4.1); middle panel: density map for our best model (SPL-TRqv, see Sec. 5.1);
right-hand panel: relative data-model residuals. The black iso-density
contours are plotted with interval of $\Delta\log\rho=0.4$ from $\log\rho=-3$
to $\log\rho=4$, where $\rho$ has the dimension $\text{kpc}^{-3}$. The model
map has been obtained averaging the maps of 1000 different mock catalogues
made assuming the best-fit parameters obtained for the SPL-TRqv halo model
(see Tab.3). Figure 19: Same as in Fig. 18 but for the number density in the
Galactic $\text{r}_{\text{e}}-\theta$ plane. The $\text{r}_{\text{e}}$ in the
plot have been obtained assuming $\text{q}=0.6$ and $\text{p}=1$ (see Sec.
3.3.3).

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
3-D 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 (Sec. 3) and the maximum likelihood fit (Sec.
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 a single power-law (SPL, Fig. 8). While - as
evidenced by Tab. 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 (Sec. 4.4.1) as well
as the direct comparison between the properties of the data and the best-fit
models, we chose to highlight this SPL-TRqv (last row in Tab. 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-
fit model has the following properties: the radial density is an SPL power-law
with an exponent $\alpha=2.96\pm 0.05$ and the iso-density surfaces of
triaxial shape. The minor axis of the density ellipsoids is aligned with the
Galactic Z. We find strong evidence for an evolution of the vertical density
flattening: the inner parts are squashed with ($\text{q}\approx 0.57$),
however the flattening becomes less pronounced ($\text{q}\approx 0.75$) at the
border of the Galactic volume analysed in this work
($\text{r}_{\text{e}}\approx 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 panels
of the Figure, where the halo flattening is over- and under-estimated
correspondingly. However, in the middle panel (corresponding to our best-fit
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, signaling a change in the halo shape. Fig. 16 illustrates the effect
of the radial halo shape evolution with the $\theta-\text{r}_{\text{e}}$ views
of the mock RRL realisations for two of our best-fit (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 ($\text{p}=1.27\pm 0.03$) in the direction which is rotated with
respect to the Galactic Y-axis by $\gamma=-21.3^{\circ}\pm 2.6^{\circ}$
(anticlockwise).

The measured RRL (as represented by the clean sample, see Sec. 4.1)
distribution can be compared to the best-fit model in Figs. 17, 18 and 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
$\text{r}_{\text{e}}-\theta$ space. In each Figure, the right-most column
displays the map of residuals in the corresponding projection.

Reassuringly, the residual maps indicate that the model gives a reasonable
description of the RRLs distribution in the G2M sample. Overall, the Milky
Way’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%) 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 $\text{Z}_{\text{g}}>15\ \text{kpc}$ and at
about $\text{r}_{\text{e}}=35$ kpc and above $\theta=50^{\circ}$ in Fig. 19.
Most likely, this large stellar cloud is related to the Virgo over-density, as
discussed in Sec. 3.3.4.

Amongst other differences, is a mild deficiency of stars visible at the edge
of the Galactic volume (available to us) in the anti-centre direction, i.e.
around l=$180^{\circ}$. 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 over-density. We argue
that this deficiency of stars may be a hint of the steepening of the radial
density profile, but as discussed in Sec. 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 Sec.
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 (Sec. 4.4) relies on the assumption that all RRLs in
the sample have the same absolute magnitude $M_{\text{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 Eq. 22 is represented by the best-fit double
Gaussian (blue line in Fig. 2). With this addition to the model, we repeated
the fit of our best halo model (Sec. 4.3.4), this time marginalising 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 $\theta$ (Sec.
4.1), we eliminated most of the regions with high Galactic dust extinction:
the 85% of the stars in our “clean” sample has $E(B-V)<0.25$ and only the 1%
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 (Sec. 5.1) masking the regions with $E(B-V)>0.4$. We did not find
any appreciable differences with respect to the un-masked analysis.

####  5.2.3 Contamination from the Galactic disc

The final results of the fitting analysis (Sec. 4.4) have been obtained
considering the presence of only one smooth stellar component, however disc
RRLs (if any) and other disc variables, including artefacts (see Sec. 2.2.5)
could pollute our sample. Although the cut on Galactic b (Sec. 2.2.3) and
$\theta$ (Sec. 4.1) is employed to limit such contamination, we repeated the
modelling with the addition of an extra disc component. Thus, this new model
has contribution from an ellipsoidal halo with normalised density
$\tilde{\rho}_{h}$ (Sec. 4.2) as well as a double-exponential disc with
normalised density
$\tilde{\rho}_{d}\propto\text{exp}[-\text{R}/\text{R}_{\text{d}}]\text{exp}[-|\text{Z}|/\text{Z}_{\text{d}}]$.
Here $\text{R}_{\text{d}}$ and $\text{Z}_{\text{d}}$ are the radial and the
vertical disc scale-lengths. Considering the two components, the global pdf of
the stellar distribution is

$\tilde{\lambda}(m,\text{l},\text{b}|f,\vec{\mu_{\text{h}}},\vec{\mu_{\text{d}}})=(1-f)\tilde{\lambda}_{\text{h}}(m,\text{l},\text{b}|\vec{\mu}_{\text{h}})+f\tilde{\lambda}_{d}(m,\text{l},\text{b}|\vec{\mu}_{\text{d}}),$
(31)

where $\tilde{\lambda}_{\text{h}}$ and $\tilde{\lambda}_{d}$ represent the pdf
of the halo and the disc stars, while $f$ is the the disc-to-total stellar
number ratio. The pdfs that appear in Eq. 31 above have been already
marginalised over the absolute magnitude distribution of the halo and disc
stars. For the halo we used a Dirac delta centred on $M=M_{\text{RRL}}=0.525$
(Sec. 2.2.4), while for the disc a uniform distribution between $M=-2$ and
$M=5$. The absolute magnitude distribution of the stars in the disc could be
much more complicated with respect to a uniform distribution. However, making
use of mock catalogues we found that our assumption is a good choice to
minimise the bias also in the presence of a truly complicated multi-peak
distribution.

The RRL sample was modelled (Sec. 4.3.4) using Eq. 31 in the logarithmic
likelihood in Eq. 25: we considered $f$ as a free parameter, but fixed the
disc’s stellar density profile assuming values for $\text{R}_{\text{d}}$ and
$\text{Z}_{\text{d}}$. We tested a wide range of thin and thick disc models
using $\text{R}_{\text{d}}$ values between 2 and 4 kpc and
$\text{Z}_{\text{d}}$ values between 0.1 and 1 kpc. In general, we did not
find significant differences with respect to the results shown in Tab. 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 an 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% and is consistent with 0.
Therefore, we conclude that our “clean” sample (Sec. 4.1) does not harbor
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. Nonetheless, 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 $\text{Z}_{\text{g}}$. The first substructure is
clearly evident in Fig. 7 and in Fig. 9 below $\theta=20^{\circ}$. 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 22600 RRLs (see Tab. 1) and the
halo+disc stellar pdf defined in Eq. 31. The disc in the model was a realistic
double exponential distribution: we fixed the radial ($\text{R}_{\text{d}}$)
and vertical ($\text{Z}_{\text{d}}$) thin and thick disc scale-lengths to the
classical values found in literature (e.g. Das et al. 2016). Fig. 20 gives the
density distribution in the $\text{R}-\text{Z}_{\text{g}}$ plane for a
disc+halo model in which the disc has $\text{R}_{\text{d}}=2.7$ kpc and
$\text{Z}_{\text{d}}=0.2$ kpc. As illustrated by the residual map (right
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
($\text{q}\approx 0.4$), thus concealing the evidence of the radial variation
of q, as measured and discussed above (Tab.3) as part of the analysis of the
“clean” sample of RRLs (Sec. 4.1). The residuals suggest that the two-
component model gives a poor fit with a significant mismatch both at low
$\text{Z}_{\text{g}}$ for $\text{R}>10$ kpc and at high $\text{Z}_{\text{g}}$
for $\text{R}<10$ kpc. Note that these results do not change if different
values for the disc scale lengths are assumed ($\text{R}_{\text{d}}$ between 2
and 4 kpc and $\text{Z}_{\text{d}}$ between 0.1 and 1 kpc). As a further test,
we repeated the fit leaving the scale lengths of the disc model as free
parameters. This test yielded unrealistically large values of
$\text{R}_{\text{d}}\approx 25$ kpc and $\text{Z}_{\text{d}}\approx 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 “anti-truncated”, i.e. having a
flatter density profile in the outer parts. In conclusion, it seems that this
flat and elongated structure can not be modelled simply by adding a double
exponential disc with the properties expected for the Milky Way. 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 et al., 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
$\text{R}-\text{Z}_{\text{g}}$ plane for a multi-component model with an halo
and a disc; the black iso-density contours are plotted with interval of
$\Delta\log\rho=0.4$ from $\log\rho=-3$ to $\log\rho=4$, where $\rho$ are in
unit of $\text{kpc}^{-3}$. Right-hand panel: relative residuals (Data-
Model/Model) between the density in the $\text{R}-\text{Z}_{\text{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-fit parameters obtained for
the disc + SPL-DN halo model (See 5.3). The stars in the disc follow a double
exponential distribution with a radial scale length $\text{R}_{\text{d}}=2.7$
kpc and a vertical scale height $\text{Z}_{\text{d}}=0.2$ kpc.

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 $[\text{Fe/h}]<-1$ from the 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 anti-centre. 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 over-density of stars above the
Galactic plane: it is visible in the density map in Fig. 7 for
$\text{Z}_{\text{g}}>15$ kpc and in Fig. 9 for $\theta>50^{\circ}$ and large
elliptical radii; it is the cause of the mismatch between the star counts
above and below the Galactic plane for $|\text{Z}_{\text{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 Fig. 17 and Fig.
18. The spatial location and shape of this structure indicate that it is most
likely related to the Virgo over-density (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 Sec. 4.4).

###  5.4 Comparison with other works

Authors |  | Tracer
---
(1)
| Range [kpc]
---
(2)
| Slope ($\alpha$)
---
(3)
| Axial ratios
---
(4)
| Catalogue
---
(5)
Das et al. (2016) | BHB | 10-70 | $\approx-4.7\pm 0.3$ | $\text{q}=0.39\pm 0.09\xrightarrow{}0.81\pm 0.05$ | SegueII
Xue et al. (2015) | K-Giants | 10-80 | $-4.2\pm 0.1$ | $\text{q}=0.2\pm 0.1\xrightarrow{}0.80\pm 0.03$ | SDSS
Sesar et al. (2013) | RRLs | 5-30 | $\approx-2.4$ | $\text{q}\approx 0.65$ | LINEAR
Deason et al. (2011) | BHB | 10-80 | $-2.3\pm 0.1\xrightarrow{\text{R}_{\text{b}}\approx 27\ \text{kpc}}-4.6\pm 0.2$ | $\text{q}=0.59\pm 0.02$ | SDSS
De Propris et al. (2010) | BHB | 10-100 | $\approx-2.5$ | $\text{q}\approx 1.0$ | 2dF QSO
Watkins et al. (2009) | RRLs | 5-110 | $-2.4\xrightarrow{\text{R}_{\text{b}}\approx 25\ \text{kpc}}-4.5$ |  q=1 (assumed) | Stripe82
Jurić et al. (2008) | MSTO | 0-20 | $\approx-2.8$ | $\text{q}\approx 0.6$ | SDSS
Miceli et al. (2008) | RRLs | 3-30 | $-3.15\pm 0.07$ |  $\text{q}=0.5\xrightarrow{}1.0$ (assumed) | SegueII
Vivas & Zinn (2006) | RRLs | 4-60 | $\approx-2.8$ | $\text{q}\approx 0.6$ | QUEST-1
This work | RRLs | 1.5-28 | $-2.96\pm 0.05$ |  | $\text{q}=0.57\pm 0.02\xrightarrow{}0.84\pm 0.06$
---
$\text{p}=1.27\pm 0.03$
Gaia DR1 + 2MASS
Table 4: Galactic halo properties found in a sample of literature works. 1:
used halo tracers, Blue horizontal branch stars (BHB), K-Giants stars
(K-Giants), RR Lyrae stars (RRLs), stars near the main-sequence turn-off
(MSTO). 2: range of sampled Galactocentric distances. 3: exponent of the
power-law describing the density profile of the stars in the halo; when two
values are present they represent the exponents in the inner part and in the
outer part of the halo separated by the transition cylindrical radius
$\text{R}_{\text{b}}$. 4: results on the shape of the halo iso-density
surfaces, q represents the the Z-to-X axial ratio, p the Y-to-X axial ratio
(if not shown $\text{p}=1$); when two values are present they represent the
axial ratio in the inner part and in the outer part of the halo. 5: Survey or
Catalogue used, SDSS (Abazajian et al., 2009), Stripe82 (Sesar et al., 2010),
SegueII (Yanny et al., 2009), 2dF QSO (Shanks et al., 2000), QUEST-1 (Vivas et
al., 2004), LINEAR (Stokes et al., 2000), Gaia DR1 (Gaia Collaboration et al.,
2016b), 2MASS (Skrutskie et al., 2006).

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-TRqv, see
Sec. 5.1) can not 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-colored 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 (Sec. 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-fit p in Deason et al. (2011) is smaller than 1, while we infer
$\text{p}\approx 1.3$. Pérez-Villegas et al. (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 ($\text{D}_{\text{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 (Sec. 4.4.3) and non-parametric (Sec. 3.3.3)
analysis 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 Sec. 3.3.3 and Fig. 14). The results of Xue et al. (2015) have
been obtained by employing the same functional form
$\text{q}(\text{r}_{\text{e}})$ used in this work (Eq. 15), while Das & Binney
(2016) used a non parametric measure of $\text{q}(\text{r}_{\text{e}})$
estimated by the iso-density contours of their best-fit 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-fit 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 fit parametric functional forms of both Xue et al.
(2015) and Das et al. (2016) are compatible with our measurement of the
flattening trend for $\text{r}_{\text{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 $\text{r}_{\text{e}}$ (see
Eqs. 9 and 15). Second, the largest discrepancies are at small radii
($\text{r}_{\text{e}}<15$ kpc): in our best-fitting model q is almost constant
at $\sim 0.6$, while in the other works the halo flattening keeps increasing
(q decreasing) towards the centre. For example, around the Galactic centre,
($\text{q}\approx 0.4$ in Das et al. 2016 and $\text{q}\approx 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.
(2017) 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 ($\text{q}>0.8$) already at small elliptical radii
($\text{r}_{\text{e}}\approx 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 Tab. 3) is highly compatible with the previous
results relying on RRLs (see Tab. 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 ($\alpha\approx-1$) in the very inner part of the halo. However, the
distribution of their RRLs in this region is very clumpy and can not be
properly fitted by a simple smooth power-law model. Similarly, Watkins et al.
(2009) found that the overall distribution of RRLs can not 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 RGB stars and to
Jurić et al. (2008) who used MS main-sequence and MSTO main-sequence turn-off
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
($\alpha>-4$) beyond 20-30 kpc. Note however that the work of De Propris et
al. (2010) found a more shallow profile ($\alpha\approx-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 artifact 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 Sec. 4.4.3). Table 4 summarises 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 April
2018555https://www.cosmos.esa.int/web/gaia/release. 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_{\text{obs}}$) for each star will increase, thus making our simple
variability estimator AMP (Eq. 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 datasets
(Sec. 2.2.2);

  * •

most importantly, photometric colours $G_{\text{BP}}$ and $G_{\text{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 proper motions 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.

  1. 1.

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
(Eq. 1) derived from the Gaia’s mean flux and its associated error.

  2. 2.

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 21643 stars
within a sphere of 20 kpc centred on the Sun, covering an unprecedented
fraction of the volume ($\sim 58\%$) of the inner halo ($\text{R}<28$ kpc).
While the overall completeness of our sample can not 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 $\sim$10%.

  3. 3.

Assuming a constant absolute magnitude for all RRLs in our sample
($M_{\text{RRL}}=0.525$), we analysed their density distribution in the inner
Galaxy. To that end, Figure 6 presents the first all-sky stellar halo maps.
These heliocentric slices reveal a relatively smooth nearby halo, with
prominent inter-quadrant asymmetries starting to become noticeable beyond 10
kpc from the Sun.

  4. 4.

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 $\text{Z}_{\text{g}}$. It contains a
little less than half of the stars in our sample, but given its large radial
extent it can not 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 over-density