# The Photo-Astrometric Vertical Tracer Density of the Milky Way II: Results
from Gaia.

Andrew Everall,1 Vasily Belokurov1, N. Wyn Evans1, Douglas Boubert2,3, Robert
J. J. Grand4,5,6
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
2Magdalen College, University of Oxford, High Street, Oxford OX1 4AU, UK
3Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, Parks
Road, Oxford OX1 3PU, UK
4Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, Postfach
1317, D-85741 Garching, Germany
5Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, E-38205 La
Laguna, Tenerife, Spain
6Departamento de Astrofísica, Universidad de La Laguna, Av. del Astrofísico
Francisco Sánchez s/n, E-38206, La Laguna, Tenerife, Spain
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We use Gaia photometry and astrometry to estimate the vertical spatial
structure of the Milky Way at the Solar radius, formally accounting for sample
incompleteness (the selection function) and parallax measurement uncertainty.
Our results show impressive precision demonstrating the power of the Gaia
data. However, systematic errors dominate the parameter value uncertainties.
We thoroughly test and quantify the impacts of all systematic uncertainties.
The vertical tracer density is modelled as a sum of two exponential profiles
for the thin and thick discs, together with a spherically symmetric power-law
for the stellar halo. We constrain the thin disc scale height as
${h_{\mathrm{Tn}}=260\pm 3\,(\mathrm{stat})\pm
26\,\mathrm{pc}\,(\mathrm{sys})}$ and thick disc ${h_{\mathrm{Tk}}=693\pm
7\,(\mathrm{stat})\pm 121\,\mathrm{pc}\,(\mathrm{sys})}$. For the halo, we
obtain a power law profile with $n_{\mathrm{H}}=3.543\pm
0.023\,(\mathrm{stat})\pm 0.259\,(\mathrm{sys})$. We infer a local stellar
mass density for non-compact object stars of
${\rho_{\mathrm{local}}^{*}=3.66\pm 0.03\,(\mathrm{stat})\pm 0.52\times
10^{-2}\,\mathrm{M}_{\odot}/\mathrm{pc}^{3}\,(\mathrm{sys})}$ and surface
density of ${\Sigma_{\mathrm{local}}^{*}=23.17\pm 0.08\,(\mathrm{stat})\pm
2.43\,\mathrm{M}_{\odot}/\mathrm{pc}^{2}\,(\mathrm{sys})}$. We find
asymmetries above and below the disc with longer disc scale heights in the
north but a flatter halo in the south at the $\lesssim 10$ per cent level.

###### keywords:

Galaxy: stellar content, stars: statistics, Galaxy: kinematics and dynamics,
methods: data analysis, methods: statistical

††pubyear: 2021††pagerange: The Photo-Astrometric Vertical Tracer Density of
the Milky Way II: Results from Gaia.–A.2

##  1 Introduction

The 3D distribution of stars throughout the Milky Way is vital for
understanding the formation history of our Galaxy. This ‘tracer density’ is
also a key ingredient in methods attempting to estimate the distribution of
dark matter in the Milky Way, with important implications for both
cosmological models and direct detection experiments (e.g., Read, 2014).

There is a rich history of research into the structure of the Milky Way from
Herschel 1785’s star-gage method to Kapteyn 1922’s lens-shaped model to recent
discoveries of asymmetries (e.g., Widrow et al., 2012) and halo substructures
(e.g., Belokurov et al., 2006). However, across the vast majority of
historical studies of Milky Way structure, the missing ingredient has been
directly measured distances. As a result, we have typically been dominated by
the statistical uncertainties regarding the distances to observed sources.
Kapteyn 1922 even says “I know of no more depressing thing in the whole domain
of astronomy than to pass from the consideration of the accidental errors of
our star places to that of their systematic errors.”

The Hipparcos mission (Perryman et al., 1997) dramatically improved the
situation, providing milli-arcsecond precision parallax measurements and
propelling the field forward (e.g. Creze et al., 1998). However, the Hipparcos
catalogue contains only $118\,000$ stars so the sample size is limited.

Countering the limited size of Hipparcos, large scale photometric surveys such
as SDSS (Gunn et al., 1998) and 2MASS (Skrutskie et al., 2006) measured
precise photometry for tens to hundreds of millions of sources. Using
photometric colours and stellar evolution models, luminosities of stars can be
estimated which are used to infer distances to stars. This has been used to
infer the structure of the Milky Way (e.g. Robin et al., 2003; Bilir et al.,
2006a, b; Jurić et al., 2008).

This picture has entirely changed with the advent of Gaia. Its primary aim is
to measure the spatial and velocity distribution of over a billion stars in
the Galaxy (Gaia Collaboration et al., 2016). To do this, Gaia and DPAC111DPAC
is the Gaia Data Processing and Analysis Consortium who we have to thank for
producing the exquisite quality of data. have now provided parallax
measurements for 1,467,744,818 sources (Gaia Collaboration et al., 2021a) with
precisions down to $10^{-2}$ mas (Lindegren et al., 2021a). Given this quality
of data, one would be forgiven for thinking a detailed 3D map of the stellar
components of the Milky Way would be a straightforward task.

For reasons related to the completeness of Gaia astrometry and awkwardness of
distance uncertainties, inferring the true spatial distribution of stars from
Gaia is a difficult statistical problem. Rix et al. 2021 provides a detailed
introduction and discussion of the importance and challenges of estimating and
applying selection functions in source density models. To avoid these issues,
some studies using Gaia data to infer the structure of the Milky Way only work
with subsets and infer distances from photometry (Deason et al., 2019; Iorio
et al., 2018). Some elect to not use Gaia data at all (Mateu & Vivas, 2018;
Fukushima et al., 2019; Dobbie & Warren, 2020).

In Everall et al. 2021a, henceforth Paper I, we introduce and validate a
method to attack these problems head-on, accounting for the selection function
of the Gaia astrometry sample and parallax uncertainties. In this work, we
leverage Paper I’s method to estimate the scale height of the Milky Way thin
and thick discs, the radial profile of the halo and the local number density
of stars for each component. We model the impacts of multiple systematic
uncertainties to estimate the effect these may have on our final results.

In Section 2, we introduce our high latitude sample, which is extracted from
Gaia early data release 3 (EDR3). We describe the cuts used to remove a small
number of contaminants. Our model for the vertical tracer populations is then
fit to the Gaia data and we describe the results in Section 4. There are
various simplifications and approximations used in the method and model which
could, in principle, bias the parameter fits. These are discussed and tested
in Section 5. We explain how our tests are used to quantify statistical and
systematic uncertainties in Section 6. Finally, we discuss the results in
comparison with the literature values in Section 7.

##  2 Data

Figure 1: Cuts on Gaia and unWISE colour photometry are used to remove
extragalactic sources from our sample. Left: Galaxies are removed using a cut
on phot_bp_rp_excess_factor after correcting for colour dependence,
$C^{*}<1.8$. Red points show the SDSS spectroscopically classified galaxies
which clearly extend to high excess flux levels. Middle: Quasars are removed
with colour-colour cuts on Gaia $G_{\mathrm{BP}}-G_{\mathrm{RP}}$ and unWISE
$W_{1}-W_{2}$ shown by the blue dashed lines. The SDSS quasars (blue points)
are clearly clustered in the region of colour-colour space beyond these cuts.
Right: The parallax SNR distribution of the galaxy and quasar samples are
nearly Gaussian unit-variance distributed with a small enhancement at high
$\varpi$ due to a small number of stars which are incorrectly removed from the
sample.

Our initial sample of Gaia sources consists of all objects in EDR3 with
$|b|>80^{\circ}$, published parallax with $\mathrm{RUWE}<1.4$ and published
$G$-band apparent magnitude with $G>5$. Brighter sources saturate the Gaia
CCDs which significantly affects the reliability of astrometric solutions. Our
sample is extracted with the following query which returns 673 926 sources in
the Galactic north and 702 599 in the south.

⬇

select ra, dec, parallax, parallax_error, phot_g_mean_mag from
gaiaedr3.gaia_source

where (b<-80 or b>80)

and parallax is not NULL

and phot_g_mean_mag>5

and ruwe<1.4

A recurring challenge with Gaia astrometry is the zero-point parallax offset,
which leads to a small bias for any individual source but can significantly
bias models fit to an entire population (e.g. see Everall et al., 2019). We
apply the zero-point correction recommended in Lindegren et al. 2021b for
sources with 5 and 6 parameter astrometric solutions. Many other groups have
attempted to measure the zero point parallax offset from Cepheid variables
(Riess et al., 2021), Red Clump stars (Huang et al., 2021), eclipsing binaries
(Stassun & Torres, 2021; Ren et al., 2021) and quasars (Groenewegen, 2021)
(although the Lindegren et al., 2021b model was constructed using quasars so
it is unsurprising that these results match well). The conclusions are that
for the majority of sources, the parallax offset is reduced to under
$10\mu$as. Zinn 2021 and Riess et al. 2021 find the parallaxes of sources
brighter than $G=10.8$ are overestimated by $\sim 15\mu$as after the
correction, so we adjust the offset for the small portion of our sample with
$G<10.8$. We test and discuss the effect of any residual offset in Section 5.

Parallax errors in Gaia are found to be typically underestimated when
considering globular clusters (Vasiliev & Baumgardt, 2021) and wide binaries
(El-Badry et al., 2021). We use the model from Equation (16) of El-Badry et
al. 2021 to revise the parallax errors of our Gaia sample, as this is
appropriate for uncrowded fields which broadly applies to our sample. Close
binary systems can bias the measured parallax for individual sources
(Belokurov et al., 2020b; Penoyre et al., 2020) however this will
predominantly affect nearby sources with specific orbital parameters and
therefore not have a significant effect on our results.

The Gaia G-band apparent magnitude also has some small systematic bias for
sources with 6-parameter astrometric solutions. We apply the apparent
magnitude correction recommended in Riello et al. 2021 to the sources where
$G_{\mathrm{BP}}-G_{\mathrm{RP}}$ colour is available. One issue this raises
is that the $G$-band apparent magnitude used for the data is subtly different
from the measurements used to derive the Gaia selection function. However, the
magnitude correction is at most $-0.025$ mag which is much smaller than our
0.2 mag resolution of the selection function. Therefore this inconsistency
will have a negligible effect on the results.

As we are only using objects at high Galactic latitude, there is likely to be
a sizable contamination from extragalactic sources (both quasars and distant
galaxies). If left in the sample, these would bias the inferred distribution
of stars towards larger distances.

Classifiers have been constructed to determine the probability of a source
being extragalactic based on Gaia astrometry and photometry complemented with
other surveys (Bailer-Jones et al., 2019; Shu et al., 2019). The issue is that
these classifications are not $100\%$ pure and will likely remove dim stars
with low parallaxes which are misclassified as extragalactic. This is
particularly clear in Fig. 10 of Bailer-Jones et al. 2019 where the ‘quasar’
population is dominated by the LMC, SMC and particular scans. The most
prominent scans are the same as those found in Appendix B of Boubert et al.
2020 which were caused by missing calibration data in the Gaia photometric
processing pipeline. To avoid introducing a bias to our data when removing
extragalactic sources, we avoid selecting on apparent magnitude and
astrometry.

Galaxies have an extended flux distribution on the sky. Due to the larger
window size used to measure BP and RP on-board Gaia, galaxies will typically
produce an excess flux in these bands over the $G$-band (see Fig. 21 Riello et
al., 2021). The flux ratio between the combined BP and RP measurements and the
$G$-band is published as phot_bp_rp_excess_factor in the Gaia archive (Evans
et al., 2018). The published excess flux has some residual colour-dependence
which needs correcting. We use the formula provided in Section 6 of Riello et
al. 2021 to estimate the corrected flux excess $C^{*}$. Galaxies are selected
as sources with $C^{*}>1.8$. The distribution of sources in excess flux vs
$G_{\mathrm{BP}}-G_{\mathrm{RP}}$ is shown in the left panel Fig. 1 with the
red dashed line showing the Galaxy cut.

Quasars are well distinguished using the WISE photometry’s $W_{1}-W_{2}$
colour (e.g., Shu et al., 2019). We crossmatch our sample with the unWISE
sample which has improved resolution over the original WISE catalogue (Lang,
2014). Taking the nearest object within 2 arcseconds correcting for proper
motions with the Gaia epoch set to 2016 and unWISE to 2010 produces a
successful match for 88% of sources in our sample. Quasars are removed from
our sample using the colour-colour cut

$W_{1}-W_{2}>0.5\quad\&\quad G_{\mathrm{BP}}-G_{\mathrm{RP}}<0.7(W_{1}-W_{2})$
(1)

which is shown by the blue dashed line in the middle panel of Fig. 1.

These cuts select 2,933 galaxies and 50,726 quasars with 553 sources
classified as both a galaxy and quasar. However, this does not tell us how
successful our selection has been. For this, we crossmatch with
spectroscopically classified sources in SDSS-IV (Blanton et al., 2017). We
again use a proper motion corrected crossmatch for sources within 2 arcseconds
with the SDSS epoch set at 2000. In this case, only 1.8% of our sample receive
SDSS spectra, the vast majority of which are in the northern field. The
objects classified as galaxies and quasars by SDSS are shown as the red and
blue points respectively in the left and middle panels of Fig. 1.

Of those with successful crossmatches, 8 900 are classified as galaxies or
quasars by SDSS whilst our cuts select 8 275 sources, of which 8 114 are
classified as extragalactic by both. This implies that our selection criteria
correctly classifies 91.2% of extragalactic sources with only 1.7% of Milky
Way sources incorrectly classified as extragalactic. The remaining 8.8% of
missing sources account for $\sim 0.3\%$ of our final sample, so we consider
this completeness to be sufficient.

Extragalactic sources are far too distant for Gaia parallax measurements
therefore the measured parallax signal to noise will be distributed as
$\varpi/\sigma_{\varpi}\sim\mathcal{N}(0,1)$. We show this distribution in the
right hand panel of Fig. 1 for galaxies (red), quasars (blue) and the
remainder of the sample (purple). The extragalactic sources are close to
normally distributed. The galaxy sample has a small amount of stellar
contamination which marginally enhances the $+\varpi$ wing, but overall this
shows that our classification has performed well.

Figure 2: Top: The number density of sources in HEALPix pixels across the
north (left) and south (right) regions of the sky with $|b|>80^{\circ}$ is
mostly uniform. The two clear exceptions are NGC 288 at the south Galactic
pole and the Sculptor dwarf spheroidal at $(l,b)=(288^{\circ},-83^{\circ}$)
both appearing in the upper right panel. Bottom: After masking these
contributions we are left with the bottom panels which are almost uniform with
a slight number density gradient from towards the Galactic center at the top
to the anticenter at the bottom. Figure 3: The effect of removing sources with
$\varpi-\sigma_{\varpi}>10^{\frac{22-G}{5}}$ is shown as a function of
absolute magnitude ($M_{G}=G+5\log_{10}(\varpi)-10$, left), inverse parallax
(middle) and parallax signal-to-noise (right) for samples in the south (top)
and north (bottom) regions. The cut conservatively removes sources which are
likely to be intrinsically dimmer than the maximum absolute magnitude of the
model ($M_{G}=12$). The removed sources (blue histograms) don’t extend beyond
$1/\varpi=400$ pc and all have parallax SNR greater than unity so we can be
confident in their high absolute magnitudes. Some sources with $M_{G}>12$ will
pass this very cautious cut but we expect these will be dominated by the
number of bright sources with well measured parallax in the Solar
neighbourhood. Our remaining samples after applying this cut are shown by the
red histograms.

The number density of sources in pixels around the north and south Galactic
poles is shown in the top panels of Fig. 2. For the most part, the
distribution is reasonably smooth and noise dominated which is good when
fitting a smooth model. However, the south field has two significant
overdensities. The overdensity close to the south Galactic pole is the
globular cluster NGC 288 which sits at a distance of approximately $9$ kpc
from the Sun with a scale radius of $\sim 3$ arcminutes (Vasiliev & Baumgardt,
2021). The other overdensity at slightly higher latitudes east of the Galactic
Centre direction is the Sculptor dwarf spheroidal at
$l=288^{\circ},b=-83^{\circ}$ with a half-light radius of $\sim 11.3$
arcminutes (McConnachie, 2012).

To prevent these objects from contaminating our smooth models, we mask the
regions of the sky occupied by the structure out to four scale radii. We then
renormalise the pixels by the fraction of the area which remains unmasked.
This is the same treatment that we apply to pixels sitting on the edge of the
$10^{\circ}$ radius fields. The resulting source density after masking NGC 288
and Sculptor is given by the bottom panels in Fig. 2, showing no further
significant residual substructure. The gradient of the source density from the
Galactic Centre (top of the figure) to the outer galaxy can now be seen. This
shows the cylindrical radius dependence of the Milky Way distribution of stars
which is not factored into our model, but we discuss its impact in Section 5.

The absolute magnitude model defined in Section 3.2 of Paper I is limited by
$M_{G}<12$ in order to avoid use of uncertain stellar evolutionary models. The
Gaia sample may still contain sources dimmer than this limit, which are
nonetheless near enough that Gaia is able to detect them. The issue is that we
cannot directly measure absolute magnitude and parallax error is large enough
for many sources that they will be scattered to that region of absolute
magnitude space independent of their true brightness. Our compromise is to cut
out sources which are likely to be fainter than $M_{G}=12$ by $1\sigma$
uncertainty in parallax. In other words, removing all sources with greater
than $84\%$ likelihood of $M_{G}>12$. This means only keeping sources with

$\varpi-\sigma_{\varpi}<10^{\frac{22-G}{5}}.$ (2)

The effect of this cut is shown in Fig. 3. The left panel shows the naive
absolute magnitude distribution calculated with $s=1/\varpi$. The cut removes
a large fraction of objects which fall outside the boundary. Importantly, from
the middle and right panels, all of the sources removed from the sample are
measured with $1/\varpi$ within $400$pc of the Sun with a parallax SNR greater
than $1.4$. Any error in this cut will introduce a dependence of the selection
function on measured parallax and parallax error. However, given the high
parallax SNR of the removed sources, we expect that this dependence should be
negligibly small. An added benefit of the cut we have placed here is that it
will likely remove sources with poor astrometric solutions as classified by
Rybizki et al. 2021 and Gaia Collaboration et al. 2021b which are typically
fainter than $M_{G}=12$. This cut removes a further 13 792 and 13 731 sources
from the north and south fields respectively.

After all of the cleaning, we are left with 633 289 north and 640 072 south
sources in our sample. We emphasise that, through all of these cuts, we remove
less than $11\%$ of the sample with published parallax, $G$ apparent magnitude
and $\mathrm{RUWE}<1.4$. By comparison, a cut on $\varpi>0$ alone (which is a
serious crime, according to Luri et al. 2018) removes over $15\%$ and a
signal-to-noise cut of $\varpi/\sigma_{\varpi}>4$ removes over $61\%$. We are
modelling the vast majority of Gaia sources using the reliable astrometric and
photometric data that is available.

(a) No Truncation

(b) $s<160$ kpc

Figure 4: Fitted models and sample number densities per unit $z$ (kpc, left), $M_{G}$ (middle) and $G$ (right) for the thin disc (orange) and thick disc (green) and halo (purple) and their sum total (black). Lines show the median model fits with shaded regions providing the $1^{\mathrm{st}}-99^{\mathrm{th}}$ percentile range of the posterior fits to the Gaia data. In most cases the posterior is so tightly constrained that the uncertainties cannot be picked out in these plots. a: For the infinite halo model there is qualitative agreement between south (top) and north (bottom) disc samples with a steeper northern halo profile. Due to the large total normalisation of the infinite halo within $b>80^{\circ}$, the halo dominates the absolute magnitude profile and it sits directly under the total profile. b: The model with halo truncated such that $s<160$ kpc also has similar north and south profiles with a marginally steeper south halo. For both models the red histograms in the left panels show the distribution $\sin(|b|)/\varpi$ which is significantly different to the fit model due to a combination of the selection function and parallax error which we have demonstrated need to be treated properly ($\sim 14$ per cent of the sample has negative parallax and cannot even be plotted). Red histograms in the right column show the $G$ distribution of the data which agrees very well with the product of our model with the selection function (black dotted line). At the bright end the model slightly overestimates the data which is likely because our model does not truncate at the tip of the red giant branch (see Fig. 3 Paper I).  Figure 5: The posterior distributions for the north (red solid contours) and south (blue solid contours) sample fits show a small but significant disagreement across most parameters suggesting a weak asymmetry. Constraining the model to $s<160$ kpc (dashed contours) has a small impact on disc parameters however the halo model is much more significantly affected. Notably, the halo power-law index, which pushes close to the lower bound for an un-truncated model, is fit with a significantly steeper profile when the truncation is applied. The truncated model is better suited to the Milky Way for which the halo will not extend indefinitely. Component | Parameter | North | South | North ($s<160$kpc) | South ($s<160$kpc)
---|---|---|---|---|---
Thin disc | $w$ | ${2.24}_{-0.05}^{+0.05}\times 10^{5}$ | ${2.11}_{-0.03}^{+0.04}\times 10^{5}$ | ${1.97}_{-0.05}^{+0.05}\times 10^{5}$ | ${1.98}_{-0.04}^{+0.04}\times 10^{5}$
| $h_{\mathrm{Tn}}$ | ${0.279}_{-0.002}^{+0.002}$ | ${0.255}_{-0.002}^{+0.002}$ | ${0.269}_{-0.002}^{+0.003}$ | ${0.250}_{-0.002}^{+0.002}$
| $f_{G}$ | ${1.16}_{-0.04}^{+0.04}\times 10^{-2}$ | ${1.32}_{-0.04}^{+0.04}\times 10^{-2}$ | ${1.25}_{-0.04}^{+0.04}\times 10^{-2}$ | ${1.36}_{-0.04}^{+0.04}\times 10^{-2}$
| $M_{\mathrm{TO}}$ | 3.1 |  |  |
| $\alpha_{3}$ | -0.6 |  |  |
Thick disc | $w$ | ${6.32}_{-0.05}^{+0.05}\times 10^{5}$ | ${6.63}_{-0.04}^{+0.04}\times 10^{5}$ | ${6.16}_{-0.05}^{+0.05}\times 10^{5}$ | ${6.53}_{-0.04}^{+0.04}\times 10^{5}$
| $h_{\mathrm{Tk}}$ | ${0.766}_{-0.009}^{+0.010}$ | ${0.711}_{-0.005}^{+0.005}$ | ${0.706}_{-0.007}^{+0.007}$ | ${0.683}_{-0.005}^{+0.005}$
| $f_{G}$ | ${5.26}_{-0.16}^{+0.17}\times 10^{-3}$ | ${5.83}_{-0.14}^{+0.15}\times 10^{-3}$ | ${5.21}_{-0.16}^{+0.16}\times 10^{-3}$ | ${5.87}_{-0.15}^{+0.15}\times 10^{-3}$
| $M_{\mathrm{TO}}$ | 3.1 |  |  |
| $\alpha_{3}$ | -0.73 |  |  |
Halo | $w$ | ${8.64}_{-1.07}^{+1.69}\times 10^{6}$ | ${3.26}_{-1.38}^{+3.11}\times 10^{8}$ | ${3.19}_{-0.05}^{+0.05}\times 10^{6}$ | ${3.39}_{-0.05}^{+0.06}\times 10^{6}$
| $n_{\mathrm{H}}$ | ${3.254}_{-0.046}^{+0.041}$ | ${3.005}_{-0.003}^{+0.004}$ | ${3.705}_{-0.023}^{+0.023}$ | ${3.382}_{-0.020}^{+0.020}$
| $f_{G}$ | ${6.04}_{-0.10}^{+0.10}\times 10^{-3}$ | ${4.86}_{-0.09}^{+0.09}\times 10^{-3}$ | ${5.91}_{-0.09}^{+0.09}\times 10^{-3}$ | ${4.66}_{-0.09}^{+0.09}\times 10^{-3}$
| $M_{\mathrm{TO}}$ | 3.1 |  |  |
| $\alpha_{3}$ | -0.64 |  |  |
Shared | $\alpha_{1}$ | ${-0.1111}_{-0.0005}^{+0.0005}$ | ${-0.1065}_{-0.0004}^{+0.0004}$ | ${-0.1123}_{-0.0005}^{+0.0005}$ | ${-0.1075}_{-0.0004}^{+0.0004}$
| $\alpha_{2}$ | ${-0.2685}_{-0.0023}^{+0.0023}$ | ${-0.2763}_{-0.0022}^{+0.0022}$ | ${-0.2705}_{-0.0023}^{+0.0022}$ | ${-0.2785}_{-0.0023}^{+0.0022}$
Table 1: The median and $16^{\mathrm{th}}-84^{\mathrm{th}}$ percentile ranges
of the posterior parameter distributions from fits to the Gaia samples are
shown for the north and south samples for an infinitely extending model and
for a halo truncated at $s=160$ kpc. Across all parameters there is a
significant asymmetry between the results of fitting to the north and south
samples. The south disc profiles are steeper than the north with a smaller
scale height, however, the southern halo is significantly shallower than the
northern halo even pushing up against the prior boundary for the un-truncated
model.

##  3 Method overview

The method and model used are described in detail in Sections 2,3 and 4 of
Paper I. Here we will briefly recall the important details.

A three component model is used with a thin disc, thick disc and halo where
each component is the product of a spatial and absolute magnitude distribution

$f(l,b,\varpi,M_{G})=\sum_{c=\\{\mathrm{Tn},\mathrm{Tk},\mathrm{H}\\}}w_{c}\,\nu_{c}(l,b,\varpi,\bm{\psi}_{\nu})\,\phi_{c}(M_{G},\bm{\psi}_{\phi}).$
(3)

$w_{c}$ is the total number of stars in the given component within
$b>80^{\circ}$ or $b<-80^{\circ}$ for the north and south samples. As
discussed in Paper I, this is a significant assumption as the disc is known to
have formed over an extended period of time (e.g. Snaith et al., 2015) and the
halo is made up of multiple stellar populations (e.g. Helmi et al., 2018;
Belokurov et al., 2018; Belokurov et al., 2020a). Bovy et al. (2012) and
Mackereth et al. (2017) have demonstrated that the thin and thick discs are
not distinct in spatial structure and metallicity and both can be modelled by
a single continuous distribution of profiles. We discuss such a model and the
challenges it presents in Section 7.2, however, for this work we apply the
simpler model composed of distinct thin and thick disc profiles.

The thin and thick discs are exponential profiles in $|z|$ with scale heights
$h_{\mathrm{Tn}}$ and $h_{\mathrm{Tk}}$. The halo is a spherically symmetric
power law profile with exponent $-n_{\mathrm{H}}$.

The absolute magnitude distributions are four-part exponential profiles, one
for the giants with $M_{G}<M_{\mathrm{TO}}$ (the turn-off magnitude) and three
for main sequence dwarfs. Exponents, $\alpha_{1},\alpha_{2}$ of the lower and
upper main sequence respectively are fit but constrained to take the same
values for all components. The middle component of the main sequence is
constrained by the continuity condition on $M_{G}$ and the IMF. The giant
exponent, $\alpha_{3}=-0.60,-0.77,-0.64$ is fixed for the thin disc, thick
disc and halo respectively. $M_{\mathrm{TO}}=3.1$ is also fixed for all
components. The free and fixed parameters are all listed in Table 1. Priors on
all parameters are listed in Table 1 of Paper I.

The model is fit to the data with the likelihood function introduced in
Section 2 of Paper I. This is a Poisson likelihood function (Appendix B,
Everall & Das, 2020) accounting for the selection function of the observatory
and marginalising over parallax uncertainty.

The selection function for the Gaia DR2 source catalogue was developed in
Boubert & Everall 2020 making use of the calibrated Gaia scanning law (Boubert
et al., 2020; Boubert et al., 2021). In Everall & Boubert 2021 this was
updated with the EDR3 nominal scanning law and selection functions were
evaluated for subsets of Gaia EDR3 using the method developed in Boubert &
Everall 2021. We use the product of the selection functions for the Gaia EDR3
source catalogue and subset with measured parallax and $\mathrm{RUWE}<1.4$.
This accurately describes the incompleteness of the samples we are modelling.

To optimize the results, we first run MCMC on the priors using emcee (Foreman-
Mackey et al., 2013) with 44 walkers, 100 step burn-in and 100 steps of
sampling. Ten samples are drawn from the prior for gradient descent with
L-BFGS-B (Zhu et al., 1997) implemented in scipy. The maximum likelihood
estimate with the highest likelihood is used to start a second MCMC process
with a small Gaussian ball around the parameter values. 44 walkers and 5000
steps are used with the final 2500 steps at 5 step intervals taken for the
posteriors.

This method is used for all fits in Sections 4 and 5.

##  4 Results

| $\mu^{\mathrm{North}}$ | $\mu^{\mathrm{South}}$ | $\mu$ | $\sigma_{\mathrm{sys}}^{\mathrm{N}/\mathrm{S}}$ | $\sigma_{\mathrm{sys}}^{\mathrm{Tests}}$ |
---|---|---|---|---|---|---
$\rho^{*}_{\mathrm{local}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{3}\right)$ | ${3.28}\pm{0.03}\times 10^{-2}$ | ${4.05}\pm{0.04}\times 10^{-2}$ | ${3.66}\pm{0.03}\times 10^{-2}$ | $\pm{0.39}\times 10^{-2}$ | $\pm{0.34}\times 10^{-2}$ | ($z_{\odot}$)
$\Sigma^{*}_{\mathrm{local}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{2}\right)$ | ${21.58}\pm{0.08}$ | ${24.77}\pm{0.10}$ | ${23.17}\pm{0.08}$ | $\pm{1.59}$ | $\pm{1.84}$ | ($z_{\odot}$)
$\log_{10}\left(M^{*}_{\mathrm{Halo}}\,\left(\mathrm{M}_{\odot}\right)\right)$ | ${8.97}\pm{0.01}$ | ${8.74}\pm{0.01}$ | ${8.86}\pm{0.01}$ | $\pm{0.12}$ | $\pm{0.15}$ | ($R$)
$\rho_{\mathrm{Tn}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{3}\right)$ | ${2.86}\pm{0.03}\times 10^{-2}$ | ${3.56}\pm{0.03}\times 10^{-2}$ | ${3.21}\pm{0.03}\times 10^{-2}$ | $\pm{0.35}\times 10^{-2}$ | $\pm{0.38}\times 10^{-2}$ | ($z_{\odot}$)
$\rho_{\mathrm{Tk}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{3}\right)$ | ${4.19}\pm{0.13}\times 10^{-3}$ | ${4.92}\pm{0.11}\times 10^{-3}$ | ${4.57}\pm{0.13}\times 10^{-3}$ | $\pm{0.34}\times 10^{-3}$ | $\pm{2.17}\times 10^{-3}$ | ($A_{V}$,$\sigma_{G}$)
$\rho_{\mathrm{H}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{3}\right)$ | ${2.13}\pm{0.03}\times 10^{-5}$ | ${1.51}\pm{0.02}\times 10^{-5}$ | ${1.82}\pm{0.03}\times 10^{-5}$ | $\pm{0.31}\times 10^{-5}$ | $\pm{0.22}\times 10^{-5}$ | ($R$)
$\rho_{\mathrm{Tk}}/\rho_{\mathrm{Tn}}$ | ${0.147}\pm{0.005}$ | ${0.138}\pm{0.003}$ | ${0.141}\pm{0.005}$ | $\pm{0.000}$ | $\pm{0.075}$ | ($A_{V}$,$\sigma_{G}$)
$\rho_{\mathrm{H}}/\rho_{\mathrm{Tn}}$ | ${7.46}\pm{0.14}\times 10^{-4}$ | ${4.25}\pm{0.07}\times 10^{-4}$ | ${5.85}\pm{0.14}\times 10^{-4}$ | $\pm{1.60}\times 10^{-4}$ | $\pm{1.10}\times 10^{-4}$ | ($R$)
$\Sigma_{\mathrm{Tn}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{2}\right)$ | ${15.36}\pm{0.15}$ | ${17.82}\pm{0.13}$ | ${16.59}\pm{0.15}$ | $\pm{1.22}$ | $\pm{2.22}$ | ($z_{\odot}$)
$\Sigma_{\mathrm{Tk}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{2}\right)$ | ${5.92}\pm{0.13}$ | ${6.72}\pm{0.11}$ | ${6.33}\pm{0.13}$ | $\pm{0.38}$ | $\pm{1.67}$ | ($\sigma_{G}$,$A_{V}$)
$\Sigma_{\mathrm{H}}\,\left(\mathrm{M}_{\odot}/\mathrm{pc}^{2}\right)$ | ${0.294}\pm{0.003}$ | ${0.225}\pm{0.002}$ | ${0.259}\pm{0.003}$ | $\pm{0.034}$ | $\pm{0.019}$ | ($R$)
$h_{\mathrm{Tn}}\,(\mathrm{kpc})$ | ${0.269}\pm{0.003}$ | ${0.250}\pm{0.002}$ | ${0.260}\pm{0.003}$ | $\pm{0.009}$ | $\pm{0.024}$ | ($\sigma_{G}$,$A_{V}$)
$h_{\mathrm{Tk}}\,(\mathrm{kpc})$ | ${0.706}\pm{0.007}$ | ${0.682}\pm{0.005}$ | ${0.693}\pm{0.007}$ | $\pm{0.010}$ | $\pm{0.121}$ | ($R$)
$n_{\mathrm{H}}$ | ${3.705}\pm{0.023}$ | ${3.382}\pm{0.020}$ | ${3.543}\pm{0.023}$ | $\pm{0.160}$ | $\pm{0.204}$ | ($R$)
$f^{\mathrm{G}}_{\mathrm{Tn}}$ | ${1.25}\pm{0.04}\times 10^{-2}$ | ${1.36}\pm{0.04}\times 10^{-2}$ | ${1.31}\pm{0.04}\times 10^{-2}$ | $\pm{0.04}\times 10^{-2}$ | $\pm{0.84}\times 10^{-2}$ | ($M_{\mathrm{TO}}^{\mathrm{Tn}}$,$\sigma_{G}$)
$f^{\mathrm{G}}_{\mathrm{Tk}}$ | ${5.21}\pm{0.16}\times 10^{-3}$ | ${5.87}\pm{0.15}\times 10^{-3}$ | ${5.55}\pm{0.16}\times 10^{-3}$ | $\pm{0.29}\times 10^{-3}$ | $\pm{1.05}\times 10^{-3}$ | ($\Delta\varpi$,$A_{V}$)
$f^{\mathrm{G}}_{\mathrm{H}}$ | ${5.91}\pm{0.10}\times 10^{-3}$ | ${4.67}\pm{0.09}\times 10^{-3}$ | ${5.28}\pm{0.10}\times 10^{-3}$ | $\pm{0.61}\times 10^{-3}$ | $\pm{0.64}\times 10^{-3}$ | ($\Delta\varpi$)
$\alpha_{1}$ | ${-0.1123}\pm{0.0005}$ | ${-0.1075}\pm{0.0004}$ | ${-0.1099}\pm{0.0005}$ | $\pm{0.0024}$ | $\pm{0.0050}$ | ($\Delta\varpi$)
$\alpha 2$ | ${-0.2705}\pm{0.0023}$ | ${-0.2785}\pm{0.0022}$ | ${-0.2745}\pm{0.0023}$ | $\pm{0.0034}$ | $\pm{0.0244}$ | ($\Delta\varpi$)
Table 2: Transformed results from the model fits to Gaia data are given for
the north, south and combined samples along with one standard deviation
uncertainties. While the statistical uncertainties for each sample are
incredibly tight for most parameters, the systematic uncertainties due to
north-south asymmetry and model oversimplifications are much larger. It is
important to consider these additional systematics when using our results. The
systematic uncertainties should be added in quadrature to the statistical
uncertainty.

The model is independently fit to the northern and southern Gaia samples. This
halves the sample size in either fit but means we can draw a comparison
between the Milky Way structure above and below the disc. The method is
described in detail in Section 2 of Paper I.

The resultant model is shown in Fig. 4(a). Solid lines and shaded regions show
the median and $1^{\mathrm{st}}-99^{\mathrm{th}}$ percentile ranges for the
fits to the individual components and sum total. We evaluate this by drawing
1000 samples from the MCMC posterior, evaluating the model and taking the
percentiles as a function of $z$, $M_{G}$ and $G$. Each of the three Galaxy
components are well constrained with the thin disc dominating the model for
$z<0.5$ kpc, the thick disc being the main contribution for $0.5<z<5$ kpc and
the halo taking over at large distances. The thin and thick disc profiles are
qualitatively very similar between the north and south samples however the
halo profile in the south fit declines much more slowly with distance.

In the left panels of Fig. 4(a) red histograms show the number density of
stars as a function of $z=\sin(b)/\varpi$, which provides a biased estimate of
height above the Milky Way disc. The distribution is significantly lower than
our model at large scale heights both due to the selection function and
because parallax uncertainty scatters measurements to either larger positive
or negative observed parallax. Faint sources in Gaia have typical parallax
uncertainties $\sigma_{\varpi}>0.1$ mas and so measuring $0<\varpi<0.1$ mas
(which corresponds to $z>10$ kpc) is unlikely due purely to measurement noise.
Some of these sources are scattered up in parallax and down in distance
generating the excess of sources with measured $\varpi\sim 0.4$ kpc. This can
also be seen in Figure 5 of Paper I where the imposed Gaia-like selection
function and parallax uncertainties have the same effect on the naive
distribution of $z=\sin(b)/\varpi$. The point we are making here is that one
must account for both parallax uncertainty and the Gaia selection function to
obtain an unbiased model of the Milky Way distribution of stars.

Unlike several previous works such as Jurić et al. 2008 and Mateu & Vivas
2018, our model extends to infinity so we require $n_{\mathrm{H}}>3$ to keep
the model normalisable. However, this is unphysical and other studies have
shown that the halo drops off much steeper beyond $r\gtrsim 50$ kpc (Deason et
al., 2014) or $r\sim 160$ kpc (Fukushima et al., 2019). Our model is dominated
by information from the inner, shallower component of their profiles. This
leads to an overestimate of the overall halo normalisation, which consequently
are untrustworthy.

To obtain a more realistic halo normalisation, we rerun the fits truncating
the parallax integral and halo normalisation with $s<160$ kpc (i.e.
$1/s>6.25\,\mu$as) and changing the halo exponent prior to
$n_{\mathrm{H}}\sim\mathrm{U}[2,7.3]$. The spatial and absolute magnitude
profiles are shown in Fig. 4(b). In this case, the north and south halo
profiles are both significantly steeper.

The right hand panels of Figs. 4(a) and 4(b) show the apparent magnitude
distribution marginalised over position on the sky and distance. We weight the
total distribution by the selection function which produces the black dotted
line. This sits directly on the red histograms which give the apparent
magnitude distribution of the Gaia data. The model slightly overestimates the
apparent magnitude distribution at the bright end ($G\lesssim 7$) which we
expect is due to the truncation of the absolute magnitude distribution at the
tip of the red giant branch which can be seen at $M_{G}\sim-3$ in Fig. 5 of
Paper I but which we do not account for in our model.

The posteriors on each parameter are shown in Fig. 5 for the north and south
samples (blue and red respectively) with distance truncated fits shown with
dashed contours. Across all parameters there are systematic differences
between the results from the north and south samples. For the thin and thick
disc parameters these differences are small. However, in the case of the halo,
the effect is far more substantial. Transitioning from an infinite to a
truncated halo also significantly modifies the halo parameters with small
knock-on effects to the disc. Given previous work (Deason et al., 2014;
Fukushima et al., 2019), we consider the truncated model to be the more
appropriate and will use those fits for our final results.

The posterior median, $16^{\mathrm{th}}$ and $84^{\mathrm{th}}$ percentiles
for all components and parameters in each of the runs are given in Table 1.

###  4.1 Stellar mass density

Our parameterisation, in particular the component normalisation ($w_{c}$), is
specific to this sample as it is the total number of source with $M_{G}<12$
within the cone $|b|>80^{\circ}$. The local stellar mass density
($\rho^{*}_{\mathrm{local}}$), local surface density
($\Sigma^{*}_{\mathrm{local}}$) and halo total stellar mass
($M_{\mathrm{Halo}}^{*}$) are more generally interesting to the Galactic
dynamics community and can be estimated from our results as we will explain
here.

The number density of sources in the Solar neighbourhood with $M_{G}<12$ is
given by $w_{c}\cdot\nu_{c}(s=0)$ where subscript $c$ refers to each of the
three Milky Way components. We can inflate this to include main sequence
sources with $M_{G}>12$ using the isochrones from Section 3 of Paper I and
IMF. The isochrones translate $M_{G}=12$ to a minimum initial mass of sources
in our sample for each component, giving
$\mathcal{M}_{\mathrm{ini,min}}=0.177,0.147,0.115\,\mathrm{M}_{\odot}$ for the
thin, thick disc and halo respectively. The maximum initial mass of stars
before they reach the post-AGB evolution phase – eventually leading to a
compact object remnant and thus disappearing from our sample – is
$\mathcal{M}_{\mathrm{ini,max}}=1.083,0.980,0.801\,\mathrm{M}_{\odot}$. To get
the total pre-compact object local number density of sources, we inflate our
local number density by a factor

$X_{c}=\frac{\int_{0}^{\mathcal{M}_{\mathrm{ini,max}}}\xi(\mathcal{M}_{\mathrm{ini}})\mathrm{d}\mathcal{M}_{\mathrm{ini}}}{\int_{\mathcal{M}_{\mathrm{ini,min}}}^{\mathcal{M}_{\mathrm{ini,max}}}\xi(\mathcal{M}_{\mathrm{ini}})\mathrm{d}\mathcal{M}_{\mathrm{ini}}}$
(4)

where $\xi(\mathcal{M}_{\mathrm{ini}})$ is the IMF (we use Kroupa 2001). This
gives ${X_{c}=3.167,2.785,2.398}$ for the three components.

To estimate the local stellar mass density, we need the mean mass of sources
in the population. We can use the IMF again for this however we need to
account for stellar mass loss. We use the three component isochrones to
transform from initial mass to current stellar mass,
$\mathcal{M}(\mathcal{M}_{\mathrm{ini}})$. The stellar evolution models do not
extend all the way to zero mass so we assume any stars with
$\mathcal{M}_{\mathrm{ini}}<0.1$ experience negligible mass loss in their
lifetimes such that
${\mathcal{M}(\mathcal{M}_{\mathrm{ini}})=\mathcal{M}_{\mathrm{ini}}}$.The
mean mass of all non-compact object stars is

$\displaystyle\langle\mathcal{M}\rangle_{c}=\frac{\int_{0}^{\mathcal{M}_{\mathrm{ini,max}}}\,\mathcal{M}(\mathcal{M}_{\mathrm{ini}})\,\xi(\mathcal{M}_{\mathrm{ini}})\,\mathrm{d}\mathcal{M}_{\mathrm{ini}}}{\int_{0}^{\mathcal{M}_{\mathrm{ini,max}}}\xi(\mathcal{M}_{\mathrm{ini}})\,\mathrm{d}\mathcal{M}_{\mathrm{ini}}}$

which gives $\langle\mathcal{M}_{c}\rangle=0.174,0.168,0.155$ for the three
components. Finally, the local mass density of non-compact object stars is

$\displaystyle\rho^{*}_{\mathrm{local},c}=w_{c}\,\nu_{c}(s=0)\cdot
X_{c}\cdot\langle\mathcal{M}\rangle_{c}.$ (5)

A critical assumption we have made is that any stars born with an initial mass
larger than $\mathcal{M}_{\mathrm{ini,max}}$ will not appear in our sample. In
reality, the White Dwarf sequence extends up to $M_{G}\sim 8$ (see Rix et al.,
2021) and so there may be many White Dwarfs in our sample. However, these will
be dominated by the main sequence dwarfs of the same absolute magnitude and
will only provide a severely sub-dominant contribution to the number density
(see Fig. 2 Gaia Collaboration et al., 2021b).

We estimate the surface densities by integrating our components with respect
to $z$ with the mathematical working given in Appendix A. Since our power-law
halo has $n>3$, the total halo stellar mass is not well normalised at $r=0$.
We estimate the total halo mass, $M^{*}_{\mathrm{Halo}}$, by integrating our
halo profile for $r>1$ kpc and taking $n_{H}=2$ (uniform density) inside. As a
result, the halo mass is largely dominated by stars inside the Solar radius
and is an extrapolation of the local halo stellar mass density so this should
be taken with caution.

The means and standard deviations of these parameters using the $s<160$ kpc
fits are $\mu^{\mathrm{North}}$ and $\mu^{\mathrm{South}}$ in Table 2.

##  5 Systematic Errors

(a) Thin disc

(b) Thick disc

(c) Halo

Figure 6: Posterior parameter fits to the mock sample from Paper I for the thin disc parameters (a, top left), thick dick (b, top right) and halo (c, bottom) under alterations to the data which could introduce systematic errors. The purple solid “Good” contours in all panels show the posteriors from Paper I, fit to the sample without any imposed systematics and black dot and lines show the input parameters used to generate the sample. Adding a $z_{\odot}$ offset (blue dotted) has only a marginal impact on most parameter estimates. Gaia-like magnitude error ($\sigma_{G}$, red dot-dashed) leads us to underestimate the thin disc giant fraction. Extinction from Bayestar ($A_{V}$, orange solid) biases the model towards an overly steep halo whilst a $-10\,\mu$as parallax offset ($\Delta\varpi$, cyan dotted) has the opposite effect. Component | Parameter | Input | Good fit |  $z_{0}=21$pc | $A_{V}$ | $\sigma_{G}$ | $\Delta\varpi$
---|---|---|---|---|---|---|---
Thin disc | $w$ | ${1.20}\times 10^{4}$ | ${9.59}_{-1.35}^{+1.38}\times 10^{3}$ | ${1.01}_{-0.17}^{+0.17}\times 10^{4}$ | ${8.93}_{-1.55}^{+1.65}\times 10^{3}$ | ${8.76}_{-1.27}^{+1.52}\times 10^{3}$ | ${9.56}_{-1.37}^{+1.52}\times 10^{3}$
| $h_{\mathrm{Tn}}$ | ${0.300}$ | ${0.281}_{-0.015}^{+0.015}$ | ${0.296}_{-0.018}^{+0.017}$ | ${0.275}_{-0.018}^{+0.017}$ | ${0.272}_{-0.015}^{+0.017}$ | ${0.283}_{-0.015}^{+0.016}$
| $f_{G}$ | ${4.50}\times 10^{-3}$ | ${3.76}_{-1.30}^{+1.43}\times 10^{-3}$ | ${3.20}_{-1.27}^{+1.37}\times 10^{-3}$ | ${3.11}_{-1.49}^{+1.60}\times 10^{-3}$ | ${1.67}_{-0.98}^{+1.28}\times 10^{-3}$ | ${3.28}_{-1.36}^{+1.43}\times 10^{-3}$
| $M_{\mathrm{TO}}$ | 3.1 |  |  |  |  |
| $\alpha_{3}$ | -0.6 |  |  |  |  |
Thick disc | $w$ | ${4.30}\times 10^{4}$ | ${4.19}_{-0.19}^{+0.20}\times 10^{4}$ | ${3.97}_{-0.21}^{+0.24}\times 10^{4}$ | ${3.95}_{-0.18}^{+0.19}\times 10^{4}$ | ${4.09}_{-0.17}^{+0.19}\times 10^{4}$ | ${3.98}_{-0.17}^{+0.17}\times 10^{4}$
| $h_{\mathrm{Tk}}$ | ${0.900}$ | ${0.812}_{-0.045}^{+0.052}$ | ${0.846}_{-0.064}^{+0.072}$ | ${0.772}_{-0.050}^{+0.058}$ | ${0.780}_{-0.042}^{+0.050}$ | ${0.798}_{-0.045}^{+0.052}$
| $f_{G}$ | ${5.40}\times 10^{-3}$ | ${5.83}_{-0.66}^{+0.69}\times 10^{-3}$ | ${5.99}_{-0.72}^{+0.77}\times 10^{-3}$ | ${6.35}_{-0.76}^{+0.74}\times 10^{-3}$ | ${6.28}_{-0.65}^{+0.69}\times 10^{-3}$ | ${6.40}_{-0.70}^{+0.76}\times 10^{-3}$
| $M_{\mathrm{TO}}$ | 3.1 |  |  |  |  |
| $\alpha_{3}$ | -0.77 |  |  |  |  |
Halo | $w$ | ${9.45}\times 10^{5}$ | ${8.81}_{-0.45}^{+0.52}\times 10^{5}$ | ${8.97}_{-0.51}^{+0.69}\times 10^{5}$ | ${8.34}_{-0.43}^{+0.49}\times 10^{5}$ | ${9.05}_{-0.51}^{+0.59}\times 10^{5}$ | ${1.03}_{-0.07}^{+0.07}\times 10^{6}$
| $n_{\mathrm{H}}$ | ${3.740}$ | ${3.812}_{-0.066}^{+0.068}$ | ${3.795}_{-0.081}^{+0.076}$ | ${3.889}_{-0.073}^{+0.075}$ | ${3.778}_{-0.068}^{+0.068}$ | ${3.661}_{-0.060}^{+0.063}$
| $f_{G}$ | ${3.50}\times 10^{-3}$ | ${3.48}_{-0.15}^{+0.15}\times 10^{-3}$ | ${3.49}_{-0.17}^{+0.16}\times 10^{-3}$ | ${3.53}_{-0.16}^{+0.16}\times 10^{-3}$ | ${3.59}_{-0.15}^{+0.16}\times 10^{-3}$ | ${3.92}_{-0.16}^{+0.16}\times 10^{-3}$
| $M_{\mathrm{TO}}$ | 3.1 |  |  |  |  |
| $\alpha_{3}$ | -0.64 |  |  |  |  |
Shared | $\alpha_{1}$ | ${-0.110}$ | ${-0.110}_{-0.002}^{+0.002}$ | ${-0.110}_{-0.002}^{+0.002}$ | ${-0.111}_{-0.002}^{+0.002}$ | ${-0.108}_{-0.002}^{+0.002}$ | ${-0.105}_{-0.002}^{+0.002}$
| $\alpha_{2}$ | ${-0.250}$ | ${-0.252}_{-0.008}^{+0.009}$ | ${-0.252}_{-0.009}^{+0.009}$ | ${-0.257}_{-0.009}^{+0.009}$ | ${-0.246}_{-0.009}^{+0.009}$ | ${-0.228}_{-0.009}^{+0.009}$
Table 3: The results of systematics due to shifting the Sun from the Galactic
plane ($z_{\odot}=21$pc), introducing dust extinction ($A_{V}$), adding
apparent magnitude error ($\sigma_{G}$) and introducing a $-10\,\mu$as zero
point parallax offset ($\Delta\varpi$) are tested for the mock sample. “Good
fit” provides the results from Paper I where the parameters are fit to data
drawn from the same model with no systematics imposed. Most tests have only a
marginal impact on parameter estimates with $\lesssim 2$ standard deviation
offsets from the input parameters. Magnitude error significantly affects thin
disc giant fraction whilst $A_{V}$ and $\Delta\varpi$ have a stronger impact
on the halo parameters.

There are various aspects of the model which may lead to systematic errors in
the posterior parameter fits. In most cases, these originate from
simplifications to make the optimization computationally tractable. Here, we
address some of the aspects which are capable of biasing the results and test
the significance of their impact on the inferred parameters.

The tests are all performed using mock catalogues. The first four systematic
tests (Sections 5.1-5.4) use the exact same sample as Section 4 of Paper I,
but resampling from the selection function where apparent magnitudes are
altered. In Fig. 6 and Table 3 we have provided the posteriors of the “SF &
$\sigma_{\varpi}$” fit from Paper I for comparison and labeled it “Good” as
this was fit under ideal circumstances where the data correctly represents the
model. The tests in Sections 5.5 and 5.6 use re-sampled catalogues applying
the same method as Section 4 of Paper I including parallax uncertainties from
the Astrometric Spread Function (Everall et al., 2021b). There is a level of
statistical error in population sampling which affects the posteriors for
tests when a new catalogue is generated.

###  5.1 Solar vertical offset

In our model, we assume that the Sun sits directly on the mid-plane of the
Milky Way and as such we have a symmetric view of the Galaxy towards the north
and south. In fact, the Sun is slightly vertically offset from the Galactic
plane to the north by $\sim 14-21$pc (Binney et al., 1997; Joshi, 2007;
Widmark & Monari, 2019; Bennett & Bovy, 2019). As a result, our model assumes
the distribution of stars in the south is closer than it actually is and in
the north, too far away. This may impact the inferred scale height of the
discs.

To test the significance of this assumption, we use our mock sample and
introduce a vertical shift to the effective Solar position. This is done for
stars assuming the sample is entirely in the north. The vertical position is
changed for all stars such that the new coordinate, $z^{\prime}$ is given by

$z^{\prime}=z-z_{\odot},$ (6)

with $z_{\odot}=21$ pc – towards the upper end of estimates of the Solar
position offset from the Galactic plane. This reduces both the latitude and
distance of sources and therefore also reduces the apparent magnitudes

$\displaystyle\tan b^{\prime}$ $\displaystyle=\frac{s\tan b-z_{\odot}\sec
b}{s}$ (7) $\displaystyle s^{\prime}$
$\displaystyle=s\left[1+\frac{z_{\odot}}{s}\left(\frac{z_{\odot}}{s}-2\sin
b\right)\right]^{\frac{1}{2}}.$ (8)

Our latitude cut is applied on the updated latitudes, $|b^{\prime}|>80$. We do
not use the southern population, as this requires re-sampling the mock outside
the original selection bounds which would be more complicated to interpret.
For the north sample, this cut simply removes some sources from the original
data set. The source apparent magnitudes are then recomputed from their
original absolute magnitudes and the new distances, after which the selection
function is applied to the sample and finally observed parallaxes are re-
sampled from the expected uncertainties. The sample size is reduced by $\sim
1.8$ per cent over the original.

The results of the parameter fits to the new sample are given in Table 3 and
shown by the blue dashed contours in Fig. 6. The shifts of parameters from the
true values are marginally significant in some cases. Specifically, the scale
height of the thin disc is slightly increased which may be considered counter
intuitive given that we are effectively pushing the Sun closer to sources.
However, pairing this with the increased weight of the thin disc and reduced
weight of the thick disc it suggests the thin disc is taking on some thick
disc sources. Overall, these results suggest that the simplification to the
model of setting $z_{\odot}=0$ is only likely to have a marginal effect on
parameter estimates.

###  5.2 Dust extinction

Extinction due to inter-stellar dust causes stars to appear dimmer than they
would otherwise be at a given distance and absolute magnitude. This is one of
the motivations behind narrowing the sample to high-latitude regions. In these
areas, the effects of dust extinction are small for any individual star.
However, since this systematically affects all sources in the same direction,
there can still be a sizable affect on the model parameter estimates.

Why don’t we use the published dust maps to de-redden the sources in the Gaia
samples in the first place? At first glance, this is an appealing suggestion,
but there is a subtle issue here which would lead to an underestimate in
uncertainties. Because the Green et al. 2019 extinction map is evaluated using
Gaia parallax information, as is our model, there would be a double counting
of information. The formally correct way to handle this problem is to
simultaneously fit the structure of the Milky Way and the extinction map. One
immediate challenge is a strong degeneracy between extinction, distance and
absolute magnitude of sources. This is a significantly more complex problem
and well beyond the scope of this paper.

Nonetheless, we can gauge the impact of extinction by applying a Milky Way
extinction map to the mock catalogue, re-sample the selection function from
the new observed apparent magnitudes and fit the model parameters to this.

The most detailed 3D extinction map to date for the Milky Way is that of Green
et al. 2019. This uses apparent magnitudes from a wide range of pass-bands
throughout the optical and infra-red to estimate stellar reddening, whilst
Gaia parallaxes are used to provide distance information to the model. Using
the dustmaps Python module (Green, 2018), we take a single sample of the
extinction parameter for each source. As a proxy for the reddening vector
component in the Gaia $G$ band, we use the Pan-STARRS’ $g$-filter value of
$3.518$ (Table 1 of Green et al., 2019). The mean extinction for sources in
the selected sample is $\langle\delta G\rangle\sim 0.01$.

The addition of stellar extinction from Green et al. 2019 has a marginally
significant effect on parameter estimates shown by the orange contours in the
posteriors in Fig 6 and more quantitatively in the $A_{V}$ column of Table 3.
The scale height and normalisation of the discs and halo are pushed down and
the halo is too steep. Stars further away will be more obscured by dust and
less likely to be included in our sample. Our method doesn’t account for this
so we instead fit a marginally steeper model than is actually the case. At
high latitudes, this is a small effect but if we widened our on-sky sample,
this could dramatically impact the results.

###  5.3 Magnitude uncertainty

Figure 7: The distribution of $G$-band flux error amplitude (uncertainty per
observation) of all sources in Gaia EDR3 is shown by the log-normalised grey-
scale histograms. The yellow line and shaded regions provide the median and
$16^{\mathrm{th}}-84^{\mathrm{th}}$ percentile range in 0.1mag bins. The
median is used with the scanning law to estimate the expected apparent
magnitude error for sources in Gaia as a function of position on the sky and
apparent magnitude.

The method used to fit the model assumes no apparent magnitude measurement
uncertainty. Of course this is not the case and Gaia has uncertainties on all
apparent magnitude measurements. To estimate the systematic effect of this
uncertainty, we re-sample the mock apparent magnitudes from Gaia-like
uncertainties and apply the method to the new sample.

For any source in the Gaia catalogue, the apparent magnitude in the $G$-band
is estimated up to nine times whenever it is scanned by the nine columns of
CCDs in the field of view. The set of all apparent magnitude measurements,
which can number in the hundreds, is used to estimate the magnitude
uncertainty. We reverse engineer this process to estimate the apparent
magnitude uncertainty per observation. The amplitude of apparent magnitude
measurement uncertainty is given by

$\xi_{F_{G}}=\sqrt{N}\frac{\sigma_{F_{G}}}{F_{G}}$ (9)

where $F_{G}$ is the measured source flux and $\xi_{F_{G}}$ is the flux error
per observation (Eq. 2 Belokurov et al., 2017). We estimate this for all
sources in Gaia EDR3 and take the median as a function of apparent magnitude
shown in Fig. 7.

To estimate apparent magnitude errors for sources in the mock sample, we then
replicate Gaia’s observations for those sources again. The error per
observation is taken from the median in Fig. 7 and the number of scans of the
source is given by the number of scans of that position on the sky in the EDR3
nominal scanning law222Gaia EDR3 nominal scanning law:
http://cdn.gea.esac.esa.int/Gaia/gedr3/auxiliary/commanded_scan_law/. Ideally,
on average, $62/7$ observations are taken with each scan as there are nine CCD
columns but in one of the seven rows, a CCD is replaced by a wave-guide sensor
(Gaia Collaboration et al., 2016). However, Gaia is not 100% efficient and not
all observations are successfully recorded or make it through the data
processing pipeline. To account for this, we also multiply by the efficiency
at a given magnitude taken from Boubert & Everall 2020. Whilst this efficiency
is estimated from Data Release 2 (DR2), it should give a rough approximation
of the behaviour in EDR3. Eq. 9 is now reversed to estimate the flux error

$\sigma_{F_{G}}=\sqrt{N(l,b,G)}\,\xi_{F_{G}}(G)\,F_{G}(G).$ (10)

where $N(l,b,G)$ is the product of the number of scans at the given position
on the sky and $62/7$ times the observation efficiency. This is used to sample
an observed flux $F_{G}^{\prime}\sim\mathcal{N}(F_{G},\sigma_{F_{G}})$. The
mean magnitude change for all sources in our sample is $\langle|\delta
G|\rangle\sim 0.003$. Finally, we apply the selection function to the newly
estimated magnitudes which replicates the fact that the selection function was
estimated as a function of measured apparent magnitude.

The introduction of apparent magnitude error has marginal effects on some
parameter posteriors, shown by the red dot-dashed contours in Fig. 6 and
included in Table 3 as $\sigma_{G}$, with the most significant being the thin
disc dwarf fraction. This may be explained by a blurring of the sharp dwarf-
giant absolute magnitude boundary which, given the low giant fraction in the
thin disc, will lead to more dwarfs being estimated as giants than vice versa
and reduce the dwarf fraction.

###  5.4 Parallax offset

As discussed in Section 2, a significant amount of work has been devoted to
constraining the zero-point parallax offset of the Gaia astrometry sample.
However, most tests are applied on sources at the bright end of the Gaia
magnitude range. At the faint end, the correction from Lindegren et al. 2021b
reduces the parallax bias to $\sim$ a few micro-arcseconds as a function of
apparent magnitude. However, as can be seen in the third panel of Fig. 2 in
Lindegren et al. 2021b, there are variations over the sky of $\sim 10\,\mu$as
towards the north and south Galactic poles.

To test the impact of a residual parallax offset, we subtract $10\,\mu$as from
the parallax measurements in our mock sample and rerun the fits without
correcting for this. The posteriors are shown in the final column of Table 3
and with cyan dotted contours in Fig. 6. The parallax bias has no significant
impact on the thin disc and thick disc parameters however the impact on the
halo is considerable. Because this is tested using the same sample as the
“Good” fits, we are comparing with those posteriors rather than the input
parameters. In the Fig. 6(c) we can see how much the parallax offset shifts
the halo parameter posteriors away from the Good results towards a shallower,
more extended halo profile. This is not too surprising. A constant parallax
offset is small compared to the true source parallax for nearby sources, but
becomes much more significant with increasing distance.

A similar effect was found in Everall et al. 2019 when measuring the local
velocity ellipsoid with Gaia DR2. As we stray further from the Solar
neighbourhood, the impact of a negative parallax offset becomes more
significant, systematically overestimating the distances to sources. In our
case, this causes an overestimate of the radial extent of the halo.

(a) Thin disc

(b) Thick disc

(c) Halo

Figure 8: Posterior distributions for parameters fit to mock samples with
systematic differences to the type of model assumed in the fitting procedure.
Ground truth input parameters for the mock sample are shown by the black lines
and points for the thin disc (a, top left sub-figure), thick disc (b, top
right sub-figure) and halo parameters (c, bottom sub-figure). The fits to the
mock sample with a shifted thin disc turn-off
($M_{\mathrm{TO}}^{\mathrm{Tn}}$, blue contours) overestimate the population
of thin disc giants which is unsurprising as the fits assume a giant is any
source with $M_{G}<3.1$. Orange dashed contours are the fits to a mock sample
with cylindrical radius dependence and an oblate halo ($q=0.9$, orange
dashed). This significantly impacts all parameters producing extended thin and
thick discs with overestimated scale heights and an overly-steep halo.

###  5.5 Turn-off magnitude

When setting up the model, we fixed some parameter values. The most notable is
the main-sequence turn-off which determines the absolute magnitude at which
the population transitions from dwarf-dominated to giant-dominated. For all
populations, we set this to $M_{G}=3.1$ motivated by the theoretical
isochrones, however, the turn-off magnitude is a function of stellar age. We
used a thin disc model with $\tau=6.9$ Gyr, but this is only approximately the
mean age of the thin disc. The disc has formed over a long period of time and
so is made up of sources with a wide range of ages (Snaith et al., 2015;
Fantin et al., 2021; Katz et al., 2021). This leads to a main sequence turn-
off which is extended over a range of magnitudes.

Ruiz-Lara et al. 2020 demonstrated that there was a large burst of star
formation in the thin disc $\sim 5.9$ Gyr ago which, by inspection of
isochrones, corresponds to a turn-off magnitude of $M_{\mathrm{TO}}=2.9$. To
test this, we generate a new mock Gaia-like catalogue with the same input
parameter values as described in Section 4.1 of Paper I except that
$M_{\mathrm{TO}}=2.9$ for the thin disc population. We then refit this
incorrectly assuming a fixed $M_{\mathrm{TO}}=3.1$ for all components.

The results are provided in Table 4 and shown by the blue solid contours in
Fig. 8. This systematic error increases the giant fraction of the thin disc,
whilst not significantly affecting any other parameter. By moving
$M_{\mathrm{TO}}$ lower and fitting with a higher value, we are classifying
many dwarfs as giants in the model. Therefore, the effect on $f_{\mathrm{G}}$
is unsurprising but it is reassuring to see that the remaining parameters are
not sensitive to small changes in the absolute magnitude distribution.

Component | Parameter | Input | $M_{\mathrm{TO}}^{\mathrm{Tn}}=2.9$ |  $R$-free, $q=0.9$
---|---|---|---|---
Thin disc | $w$ | ${1.20}$ | ${1.15}_{-0.13}^{+0.13}\times 10^{4}$ | ${1.41}_{-0.17}^{+0.15}\times 10^{4}$
| $h_{\mathrm{Tn}}$ | ${0.300}$ | ${0.294}_{-0.013}^{+0.013}$ | ${0.323}_{-0.015}^{+0.014}$
| $f_{G}$ | ${4.50}$ | ${7.34}_{-1.39}^{+1.56}\times 10^{-3}$ | ${4.93}_{-1.01}^{+1.12}\times 10^{-3}$
Thick disc | $w$ | ${4.30}$ | ${4.34}_{-0.23}^{+0.25}\times 10^{4}$ | ${5.05}_{-0.41}^{+0.49}\times 10^{4}$
| $h_{\mathrm{Tk}}$ | ${0.900}$ | ${0.901}_{-0.055}^{+0.061}$ | ${1.049}_{-0.085}^{+0.099}$
| $f_{G}$ | ${5.40}$ | ${5.99}_{-0.73}^{+0.75}\times 10^{-3}$ | ${5.58}_{-0.69}^{+0.71}\times 10^{-3}$
Halo | $w$ | ${9.45}$ | ${8.78}_{-0.48}^{+0.59}\times 10^{5}$ | ${8.39}_{-0.40}^{+0.48}\times 10^{5}$
| $n_{\mathrm{H}}$ | ${3.740}$ | ${3.791}_{-0.070}^{+0.070}$ | ${3.957}_{-0.089}^{+0.086}$
| $f_{G}$ | ${3.50}$ | ${3.53}_{-0.15}^{+0.15}\times 10^{-3}$ | ${3.58}_{-0.16}^{+0.16}\times 10^{-3}$
Shared | $\alpha_{1}$ | ${-0.110}$ | ${-0.108}_{-0.002}^{+0.002}$ | ${-0.111}_{-0.002}^{+0.002}$
| $\alpha_{2}$ | ${-0.250}$ | ${-0.249}_{-0.008}^{+0.009}$ | ${-0.238}_{-0.008}^{+0.008}$
Table 4: We provide the median parameter estimate with
$16^{\mathrm{th}}-84^{\mathrm{th}}$ percentiles for the mock sample fits with
thin disc $M_{\mathrm{TO}}=2.9$ and the cylindrical radius dependent sample
with halo oblateness $q=0.9$.

###  5.6 Galactocentric radius and an oblate halo

The model used in this work has no dependence on Galactocentric cylindrical
radius or azimuth. Due to the complexity of the model integration, we made the
simplifying assumption that all sources have the same cylindrical radius as
the Sun, $R_{\odot}$. For a detailed discussion of this, see Appendix A of
Paper I.

For the thin and thick disc profiles, this means approximating the Milky Way
disc as a uniform sheet. In the most extreme case, where sources are against
the edge of the cone with $l=0^{\circ}$ or $180^{\circ}$, the cylindrical
radius is incorrect by

$\delta R=z/\tan 80^{\circ}=0.176\,z.$ (11)

Whilst this maximum offset is significant, it does not provide much
information on how the disc profile will affect the results. For that, we
examine the mean cylindrical radius offset integrating over the disc profile.
Using radial scalelengths of $L_{\mathrm{Tn}}=2.6$ kpc and
$L_{\mathrm{Tk}}=3.6$ kpc for the thin and thick disc respectively from Jurić
et al. 2008, we draw a sample within our $|b|>80^{\circ}$ cone at fixed $z$
for a thin and thick disc weighted by $\exp(-R/L)$ and estimate the mean $R$.
For the thin disc, we get $\delta R\sim 0.2$ pc at $z=0.3$ kpc, whilst the
thick disc produces $\delta R\sim 1.7$ pc at $z=0.9$ kpc. Therefore, the
average offset of sources from their true position is small.

The halo spatial distribution is defined as a power-law profile of
Galactocentric spherical radius. An incorrect cylindrical radius leads to an
incorrect spherical radius. The spherical radius will be incorrect by

$\delta r=\sqrt{z^{2}+R_{\odot}^{2}}-\sqrt{z^{2}+y^{2}+(R_{\odot}-x)^{2}},$
(12)

where $x,y,z$ are the standard Galactic Cartesian coordinates with $x$
positive towards the Galactic Centre. At the edges of our cone with
$l=0^{\circ}$ and $180^{\circ}$ with $z=1$ kpc, this corresponds to $\delta
r\sim+0.175$ and $-0.175$ kpc respectively. This increases to $\delta
r\sim+1.05$ and $-1.19$ kpc for $z\sim 10$ kpc. To test the impact of this on
the model fits, we sample a halo profile with $n_{\mathrm{H}}=3.724$ within
the $b>80^{\circ}$ vertical cone and use the sample to estimate the mean
spherical radius error. We find that the spherical radius of sources is
underestimated by $\sim 0.3$ per cent on average.

We also assumed our halo was spherically symmetric and ignored any flattening.
This assumption may impact the inferred steepness of the halo. Many works have
also demonstrated that the stellar halo exhibits a tri-axial distribution
(Iorio et al., 2018; Iorio & Belokurov, 2019; Naidu et al., 2021). Given the
narrow vertical cone we have used, our results will primarily be sensitive to
the halo oblateness parameter so we only test that in this work.

We test the impact of placing all sources at the Solar radius and assuming a
spherical halo by regenerating our mock sample with the correct cylindrical
radius with an oblate halo $q=0.9$ (Mateu & Vivas, 2018). Using the same
parameters as discussed in Section 4.1 of Paper I, we generate a mock
catalogue with thin and thick disc scale lengths of $L_{\mathrm{Tn}}=2.6$ and
$L_{\mathrm{Tk}}=3.6$ kpc (Jurić et al., 2008). We also include the
cylindrical radius dependence of the spherical radius for the halo model.

The results are shown in Table 4 and orange dashed contours of Fig. 8. This
has had a more significant impact on the fits with the posteriors $\sim 2-3$
sigma off the input parameter values in several cases. The disc scale heights
and halo power law profile all have significant offsets from the true input
parameters. It appears that the dominant effect is a level of source confusion
between the components. The sample is no longer exactly representative of our
assumed vertical exponential profile and power law halo but instead a
marginalisation over this with radius. This leads to overestimated disc
normalisation and scale heights and an overly steep halo.

##  6 Statistical and Systematic uncertainties

We have produced fits to the observed data around the northern and southern
Galactic poles and to mock samples to test the effects of limitations in the
model. To provide results which are informative and usable, we will quantify
what our results mean for the true model parameters and their statistical and
systematic uncertainties.

We do this by assuming the posterior distributions for all parameters in Table
2 are drawn from independent normal distributions. This enables us to
parameterise all posteriors with a mean ($\mu$) and uncertainties ($\sigma$).

The first two columns of Table 2 provide the means and standard deviation
uncertainties of the MCMC posteriors for the north and south fits to the Gaia
data. We model the combined north/south posteriors as being drawn from a
normal distribution with mean $\mu$, uncertainty
$\sigma_{\mathrm{sys,N/S}}^{2}$ convolved with an additional normal
distribution, $\mathcal{N}(0,\sigma^{2})$ providing the standard deviation
uncertainty for each sample. Therefore the likelihood of the posteriors is

$\displaystyle\log\mathcal{L}$
$\displaystyle=\sum_{i,\mathrm{North}}\log\mathcal{N}\left(x_{i}\,|\,\mu,\sigma_{\mathrm{North}}^{2}+\sigma_{\mathrm{sys,N/S}}^{2}\right)$
$\displaystyle+\sum_{i,\mathrm{South}}\log\mathcal{N}\left(x_{i}\,|\,\mu,\sigma_{\mathrm{South}}^{2}+\sigma_{\mathrm{sys,N/S}}^{2}\right)$
(13)

where $x_{i}$ are the posterior samples provided by the MCMC chains. We then
maximise the log-likelihood with respect to $\mu$,
$\sigma_{\mathrm{sys,N/S}}$. The results are given in the third and forth
columns of Table 2. The method we are using here is similar to inflating
systematics until the chi-squared reaches a ‘reasonable’ value. However, we
avoid defining an arbitrary chi-squared target by instead maximising the
Gaussian log-likelihood.

Given a sample drawn from an equally weighted sum of Gaussian distributions
with the same mean but different variance, the sample variance will be the
mean of the individual component variances. Therefore, the statistical error
for the Gaia data fits is the root-mean-square of the north and south fits

$\sigma_{\mathrm{stat}}=\sqrt{\frac{\sigma_{\mathrm{North}}^{2}+\sigma_{\mathrm{South}}^{2}}{2}}.$
(14)

This is given as the statistical error on $\mu$ in Table 2.

We play a similar game with the results of our systematic test runs from
Section 5. However, in this case, we know the true parameters because we
provided the input parameters. The likelihood is given by

$\displaystyle\log\mathcal{L}$
$\displaystyle=\sum_{i,\mathrm{test}}\log\mathcal{N}\left(x_{i}\,|\,\mu_{\mathrm{true}},\sigma_{\mathrm{test}}^{2}+\sigma_{\mathrm{sys,test}}^{2}\right)$
(15)

where $\sigma_{\mathrm{test}}$ is the statistical uncertainty of the fit given
by the standard deviation of the posterior and $\mu_{\mathrm{true}}$ is the
input parameter. We maximise this with respect to $\sigma_{\mathrm{sys,test}}$
to estimate the systematic uncertainty contribution from the given test. We
then rescale the systematic errors by the measured Milky Way parameters,
$\mu/\mu_{\mathrm{true}}$, to estimate the systematic error on our fits to the
Gaia data. In Table 2, we provide $\sigma_{\mathrm{sys,tests}}$ which is the
maximum systematic error for the given parameter from our tests. We also state
the test(s) which dominate the systematic uncertainty contribution. Where more
than one test is listed, it is because they provided a similar systematic
uncertainty to within $10$ per cent.

We recommend that anyone using our results should take the root mean square
sum of all quoted uncertainties to obtain the total uncertainty on each
parameter.

##  7 Discussion

Here, we interpret the results on the structural parameters of the Milky Way
thin disc, thick disc and halo given in Table 2, comparing them with previous
work as well as considering future developments.

###  7.1 Results

The most striking thing to notice when examining our results is the comparison
between statistical and systematic uncertainties. In general, the total
systematic uncertainty is more than an order of magnitude greater than
statistical uncertainty. In some cases, it is over two orders of magnitude
larger. This demonstrates two things. First, Gaia has ushered in an era where
it is necessary to model systematic errors once considered insignificant.
Rigorous systematic analysis of the kind we have performed is essential to
provide accurate and reliable results. Secondly, the precision which can be
achieved with Gaia data is impressive. We used a deliberately constrained
sample of objects on the sky consisting of less than $0.1\%$ of the entire
Gaia catalogue and yet the precision on most parameters is more than an order
of magnitude better than anything in the literature.

We infer a local stellar mass density for pre-compact object stars with of
${\rho_{\mathrm{local}}^{*}=3.66\pm 0.03\,(\mathrm{stat})\pm 0.52\times
10^{-2}\,\mathrm{M}_{\odot}/\mathrm{pc}^{3}\,(\mathrm{sys})}$. This is smaller
than the value of $\rho_{\mathrm{local}}^{*}\approx 4.2\times
10^{-2}\,\mathrm{M}_{\odot}/\mathrm{pc}^{3}\,$ as derived (without errors
bars) in Flynn et al. 2006, using the Hipparcos and Tycho surveys, together
with the Catalogue of Nearby Stars).

We compute a surface density of ${\Sigma_{\mathrm{local}}^{*}=23.17\pm
0.08\,(\mathrm{stat})\pm
2.43\,\mathrm{M}_{\odot}/\mathrm{pc}^{2}\,(\mathrm{sys})}$. This is
significantly smaller than Bovy et al. 2012 who estimate $30\pm
1\,M_{\odot}/\mathrm{pc}^{2}$, as well as Flynn et al. 2006 who estimate
$35.5M_{\odot}/\mathrm{pc}^{2}$.

We expect that the most significant difference is that their works include
compact objects in the stellar mass estimates. In particular, Bovy et al. 2012
use the initial mass function to infer the contribution from all sources
similar to our work in Section 4.1. However, we account for stellar mass loss
and only include stars with mass low enough that they would not have evolved
into a compact object or gone supernova. Bovy et al. 2012 extrapolate to
higher mass stars which will have evolved to a compact object without
accounting for mass loss. This means their results will significantly
overestimate the total stellar mass density for evolved stars which have
undergone significant mass loss. We do not extrapolate our results to include
compact objects as there is significant uncertainty over how much of the
initial mass is kept in the final compact object remnant.

Our relative thick-to-thin disc local density ratio sits between the values of
Mackereth et al. 2017 and Jurić et al. 2008, although the systematic
uncertainties on this due to extinction and magnitude error are quite large.

Ample past research has been dedicated to estimating the scale heights of the
thin and thick disc. There is some discrepancy between studies, with thin disc
estimates in the range $h_{\mathrm{Tn}}\sim 120-300$pc and thick disc in the
range $h_{\mathrm{Tk}}\sim 500-1900$pc (Kuijken & Gilmore, 1989; Bilir et al.,
2006a; Jurić et al., 2008; Ak et al., 2008; de Jong et al., 2010; Mateu &
Vivas, 2018; Dobbie & Warren, 2020). We constrain the thin disc scale height
as ${h_{\mathrm{Tn}}=260\pm 3\,(\mathrm{stat})\pm
26\,\mathrm{pc}\,(\mathrm{sys})}$ and thick disc ${h_{\mathrm{Tk}}=693\pm
7\,(\mathrm{stat})\pm 121\,\mathrm{pc}\,(\mathrm{sys})}$. Our estimates are
broadly in agreement with Jurić et al. 2008, de Jong et al. 2010 and Mateu &
Vivas 2018 with reasonably strong constraints on the thin disc scale height
($\pm 26$ pc), but the thick disc scale height is dominated by systematic
uncertainty ($\pm 121$ pc) due to the cylindrical radius dependence.

The power law profile of the halo has received substantial attention with
typical estimates in the range $n_{\mathrm{H}}\sim 2.5-4.4$ (Saha, 1985; Gould
et al., 1996; Yanny et al., 2000; Newberg & Yanny, 2006; Jurić et al., 2008;
Smith et al., 2009; de Jong et al., 2010; Cohen et al., 2017; Iorio et al.,
2018; Mateu & Vivas, 2018; Hernitschek et al., 2018). Our model sits in the
middle of these estimates with $n_{\mathrm{H}}=3.542\pm
0.023\,(\mathrm{stat})\pm 0.259\,(\mathrm{sys})$. Recent works have suggested
the halo is better represented by a broken power law distribution (Deason et
al., 2011; Thomas et al., 2018; Fukushima et al., 2019). Whilst we have
focused in the inner halo by the definition of Fukushima et al. 2019 and
truncated at $s=160$ kpc, there is a wide range of distances inferred for the
truncation, e.g. $25$ kpc (Watkins et al., 2009) and $42$ kpc for (Cohen et
al., 2017). In reality we expect we are covering both sides of the break
especially considering many stars from other author’s samples will likely have
made it into the Gaia astrometry sample. Another issue, extensively discussed
in the literature, which can impact halo fits is accreted substructure (Bell
et al., 2008). Our current understanding of the halo is that it is composed of
stars from from GES (a major merger event $\sim 8$ Gyr ago Helmi et al., 2018;
Belokurov et al., 2018), the “Splash” (in-situ stars kicked up by the merger
event, Belokurov et al., 2020a; Grand et al., 2020) and other accreted
substrucuture such as Sagittarius and smaller streams. We masked problematic
regions of the southern field in Section 2, however, there are likely to be
more diffuse substructures which are not so easy to mask.

Another notable feature of our results is the north-south asymmetry across
several parameters. We find the northern thin disc scale height is larger than
the south at $\lesssim 10$ per cent. Dobbie & Warren 2020 found a similar
asymmetry, although they claimed a much larger $25$ per cent difference. We
also find that the southern halo is significantly shallower with a smaller
power law exponent than the north which was also seen by Hernitschek et al.
2018. These effects may be caused by dynamical instabilities which
asymmetrically excite the disc (Widrow et al., 2012; Antoja et al., 2018) and
diffuse halo substructure such as Sagittarius which may contribute many more
stars to the southern high latitude field (Vasiliev et al., 2021).

We used our halo local mass density and profile to estimate the total halo
stellar mass, using a flat uniform density for $r<1$ kpc in order to prevent
the integral from diverging, obtaining $M_{\mathrm{Halo}}^{*}\sim 7.2\times
10^{8}\,\mathrm{M}_{\odot}$. This is quite a rough estimate of the total halo
stellar mass however our results do agree reasonably well with the broad range
of literature results $M_{\mathrm{Halo}}^{*}\sim 2-14\times 10^{8}\,M_{\odot}$
(Bell et al., 2008; Deason et al., 2011, 2019).

The dwarf fractions and absolute magnitude profiles were defined specifically
for this work and were mainly fit as nuisance parameters in order to get at
the spatial distribution of stars so we do not discuss these in detail here.

###  7.2 More general models

Our results have small statistical uncertainty compared with the dominant
systematic uncertainty. This tells us that the model we have chosen to fit is
over-constrained by the data. The solution to this is to significantly
increase the amount of freedom in the model until our systematic and
statistical uncertainties are comparable.

Some generalizations of our model are obvious: inclusion of radial dependence;
provision of free parameters on the radial disc profile; allowing a free Solar
vertical position $z_{\odot}$; introduction of halo oblateness as a free
parameter. We have not provided these freedoms due to numerical complications
in the parallax error integral discussed in Section 2.1 and Appendix A of
Paper I.

We could also move away from the simplistic two-component disc model towards a
continuous distribution of disc profiles

$\rho(z,h_{z})=\Sigma(h_{z})\,\frac{\exp\left(-|z|/h_{z}\right)}{2h_{z}}$ (16)

where $\rho$ is the number density of stars as a function of position and
scale height and $\Sigma(h_{z})$ is the surface density as a function of scale
height. The distinct thin and thick disc model we have applied is the simple
case with

$\Sigma(h_{z})=\Sigma_{\mathrm{Tn}}\,\delta(h_{z}-h_{\mathrm{Tn}})+\Sigma_{\mathrm{Tk}}\,\delta(h_{z}-h_{\mathrm{Tk}})$
(17)

where $\Sigma_{\mathrm{Tn}}$ and $\Sigma_{\mathrm{Tk}}$ are the surface
densities of the thin and thick discs respectively. Motivated by spectroscopic
data Bovy et al. (2012) proposed a model where

$\Sigma(h_{z})=\mathcal{N}\exp\left(-h_{z}/c\right)$ (18)

such that the contribution of each exponential profile to the disc surface
density declines exponentially with a scale length $c$. The normalisation is

$\mathcal{N}=\frac{\Sigma_{\mathrm{disc}}}{c\left(\exp(-h_{\mathrm{min}})-\exp(-h_{\mathrm{max}})\right)}$
(19)

where $\Sigma_{\mathrm{disc}}$ is the total disc surface density from
integrating over scale heights from $h_{\mathrm{min}}$ to $h_{\mathrm{max}}$.
The total density profile as a function of $z$ generated by this scale height
distribution is given by

$\nu(z)=\int_{h_{\mathrm{min}}}^{h_{\mathrm{max}}}\frac{\Sigma_{\mathrm{disc}}}{4h_{z}c}\,\exp\left(-\frac{h_{z}}{c}-\frac{z}{h_{z}}\right)\,\mathrm{d}h_{z}.$
(20)

There is no analytic solution to this integral such that we cannot write down
a simple formula for the vertical density profile. This is a significant issue
for the parallax integration where our numerical integration requires a
unimodal integrand as discussed extensively in Paper I.

One further complication for this approach is that the continuous density
profile is correlated with the age and metallicity of the stellar population
which determine the luminosity function.

Although beyond the scope of this work, continuous disc profile models would
be a worthwhile avenue to pursue, whether that entails finding forms of
$\Sigma(h_{z})$ which produce analytic density profiles or developing
numerical methods which can fit more general models.

Alternatively, we could go for a much more data-driven approach and fit the
source density at nodes with a smooth model such as a Gaussian Process to
enable correlations between neighbouring points.

More ambitiously still, we could leverage the BP and RP photometry provided in
Gaia EDR3 for 1.5 billion sources. Rather than using our simple magnitude
model, we could directly infer the population of the HR diagram as a function
of position in the Milky Way, from which the star formation histories and
metallicity distributions could also be inferred. This will require selection
functions for the BP and RP photometry samples which have not yet been
produced.

###  7.3 Extragalactic Component

In this work, we filtered extragalactic sources from our sample using cuts on
colour and excess flux. However, another option is to add an additional
component to the model for sources at infinite distance.

The spatial distribution of the extragalactic sources would simply be

$\nu_{\mathrm{EG}}(l,b,s)\mathrm{d}V=\frac{1}{2\pi(1-\sin(b_{\mathrm{min}}))}\,\delta(1/s)\,s^{2}\,\mathrm{d}l\,\mathrm{d}\sin(b)\,\mathrm{d}s$
(21)

where extragalactic sources have zero parallax and are uniformly distributed
across the sky for $b>b_{\mathrm{min}}$ or $b<-b_{\mathrm{min}}$ with
${b_{\mathrm{min}}=80^{\circ}}$ in this work.

However, this spatial model needs an apparent magnitude distribution for all
extragalactic sources to which the selection function can be applied. This
adds significantly more complexity, as the apparent magnitude distribution is
dependent on the distance and luminosity distribution which are different for
quasars and galaxies. For this reason, we have not chosen to model the
extragalactic population in this work. However, adding this additional
component would be an interesting and worthwhile route forward, measuring the
population of galaxies and quasars as a function of apparent magnitude with
Gaia data.

##  8 Conclusions

We used the Gaia Early Data Release 3 (EDR3) photometry and astrometry to
model the vertical distribution of stars in the Milky Way at the Solar radius.
Our sample includes the majority of stars with measured parallax in Gaia
within $10^{\circ}$ of the Galactic north and south Poles. Our method formally
accounts for parallax measurement uncertainty and the Gaia selection function
(see the companion Paper I for the algorithmic details).

We represent the vertical density of the thin and thick discs by exponentials
with scale heights ${h_{\mathrm{Tn}}}$ and ${h_{\mathrm{Tk}}}$ respectively.
The stellar halo density is a power-law of spherical radius, i.e.,
$\rho\propto r^{-n_{\mathrm{H}}}$. We thoroughly test possible sources of
systematic uncertainty in our approach, in particular from oversimplifications
of the model. This enables us to quantify the systematic uncertainty
associated with all parameter estimates.

We find the scale height of the thin disc is ${h_{\mathrm{Tn}}=260\pm
3\,(\mathrm{stat})\pm 9\pm 24\,\mathrm{pc}\,(\mathrm{sys})}$, Here, the two
levels of systematic error correspond to north-south asymmetry about the
Galactic plane and simplifying model assumption (particularly the treatment of
extinction and the assumption of halo spherical symmetry). The scale height of
the thick disc is ${h_{\mathrm{Tk}}=693\pm 7\,(\mathrm{stat})\pm 10\pm
121\,\mathrm{pc}\,(\mathrm{sys})}$. Here, the larger systematic error is
controlled by our assumption that all sources have the same cylindrical polar
radius as the Sun. For the stellar halo, we are able to constrain a power law
profile of $n_{\mathrm{H}}=3.542\pm 0.023\,(\mathrm{stat})\pm 0.160\pm
0.204\,(\mathrm{sys})$.

We infer a local stellar mass density for non-compact object stars
${\rho_{\mathrm{local}}^{*}=3.66\pm 0.03\,(\mathrm{stat})\pm 0.39\pm
0.34\times 10^{-2}\,\mathrm{M}_{\odot}/\mathrm{pc}^{3}\,(\mathrm{sys})}$ and
surface density ${\Sigma_{\mathrm{local}}^{*}=23.17\pm
0.08\,(\mathrm{stat})\pm 1.59\pm
1.84\,\mathrm{M}_{\odot}/\mathrm{pc}^{2}\,(\mathrm{sys})}$. Whilst these
values are lower than previous estimates (Flynn et al., 2006; Bovy et al.,
2012), this discrepancy may be explained by the absence of any contribution
from compact object remnants to the total stellar mass. We have not included
this due to the uncertain correction for stellar mass loss, itself not well
accounted for in previous works.

We also find a north-south asymmetry with respect to the Galactic plane. The
thin and thick disc scale heights are larger in the north, and the halo
profile is shallower in the south. However, this asymmetry is only at the
$\lesssim 10$ percent level, much less than the 25 percent claimed by Dobbie &
Warren 2020.

The impressive information content of the Gaia data produces parameter
estimates with significantly improved precision over previous studies, even
for our sample using only a small region of the sky. However, systematics now
completely dominate the error budgets meaning that we need better models to
fully realise the potential of the Gaia data.

As discussed in Paper I and Section 3, the model we have applied does not
represent our current understanding of the Milky Way from dynamical and
chemical information. Therefore this work does not produce new insight into
the physics governing the formation and evolution of the Galaxy but it does
introduce a novel approach through which new insights can be found with far
greater accuracy and precision than ever before.

The approach taken here demonstrates the power of information available from
Gaia which has yet to be unlocked. There is a substantial prize if we can
control the systematic uncertainties involved with modelling the Gaia data. .

## Acknowledgements

AE thanks the Science and Technology Facilities Council of the United Kingdom
for financial support. DB thanks Magdalen College for his fellowship and the
Rudolf Peierls Centre for Theoretical Physics for providing office space and
travel funds. RG acknowledges financial support from the Spanish Ministry of
Science and Innovation (MICINN) through the Spanish State Research Agency,
under the Severo Ochoa Program 2020-2023 (CEX2019-000920-S).

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

AE is very grateful to Eugene Vasiliev who provided valuable assistance on
numerical integration methods and who’s comments on the paper content lead to
many significant improvements.

Whilst not detailed in the paper, the authors made use of the AuriGaia mock
catalogues (Grand et al., 2018) which helped to motivate the model and
approach taken.

## Data Availability

The data underlying this article are publicly available from the European
Space Agency’s Gaia archive (https://gea.esac.esa.int/archive/). The selection
function implementations used in this work are described in Boubert & Everall
2020 and Everall & Boubert 2021 and made publicly accessible through the
Python package selectionfunctions
(https://github.com/gaiaverse/selectionfunctions).

The code used to fit the model and produce all figures is made publicly
available as a GitHub repository (https://github.com/aeverall/mwtrace.git).

## References

  * Ak et al. (2008) Ak T., Bilir S., Ak S., Eker Z., 2008, New Astron., 13, 133
  * Antoja et al. (2018) Antoja T., et al., 2018, Nature, 561, 360
  * Bailer-Jones et al. (2019) Bailer-Jones C. A. L., Fouesneau M., Andrae R., 2019, MNRAS, 490, 5615
  * Bell et al. (2008) Bell E. F., et al., 2008, ApJ, 680, 295
  * Belokurov et al. (2006) Belokurov V., et al., 2006, ApJ, 642, L137
  * Belokurov et al. (2017) Belokurov V., Erkal D., Deason A. J., Koposov S. E., De Angeli F., Evans D. W., Fraternali F., Mackey D., 2017, MNRAS, 466, 4711
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Belokurov et al. (2020a) Belokurov V., Sanders J. L., Fattahi A., Smith M. C., Deason A. J., Evans N. W., Grand R. J. J., 2020a, MNRAS, 494, 3880
  * Belokurov et al. (2020b) Belokurov V., et al., 2020b, MNRAS, 496, 1922
  * Bennett & Bovy (2019) Bennett M., Bovy J., 2019, MNRAS, 482, 1417
  * Bilir et al. (2006a) Bilir S., Karaali S., Ak S., Yaz E., Hamzaoğlu E., 2006a, New Astron., 12, 234
  * Bilir et al. (2006b) Bilir S., Karaali S., Güver T., Karataş Y., Ak S. G., 2006b, Astronomische Nachrichten, 327, 72
  * Binney et al. (1997) Binney J., Gerhard O., Spergel D., 1997, MNRAS, 288, 365
  * Blanton et al. (2017) Blanton M. R., et al., 2017, AJ, 154, 28
  * Boubert & Everall (2020) Boubert D., Everall A., 2020, MNRAS, 497, 4246
  * Boubert & Everall (2021) Boubert D., Everall A., 2021, arXiv e-prints, p. arXiv:2111.04126
  * Boubert et al. (2020) Boubert D., Everall A., Holl B., 2020, MNRAS, 497, 1826
  * Boubert et al. (2021) Boubert D., Everall A., Fraser J., Gration A., Holl B., 2021, MNRAS, 501, 2954
  * Bovy et al. (2012) Bovy J., Rix H.-W., Hogg D. W., 2012, ApJ, 751, 131
  * Cohen et al. (2017) Cohen J. G., Sesar B., Bahnolzer S., He K., Kulkarni S. R., Prince T. A., Bellm E., Laher R. R., 2017, ApJ, 849, 150
  * Creze et al. (1998) Creze M., Chereul E., Bienayme O., Pichon C., 1998, A&A, 329, 920
  * Deason et al. (2011) Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903
  * Deason et al. (2014) Deason A. J., Belokurov V., Koposov S. E., Rockosi C. M., 2014, ApJ, 787, 30
  * Deason et al. (2019) Deason A. J., Belokurov V., Sanders J. L., 2019, MNRAS, 490, 3426
  * Dobbie & Warren (2020) Dobbie P. S., Warren S. J., 2020, The Open Journal of Astrophysics, 3, 5
  * El-Badry et al. (2021) El-Badry K., Rix H.-W., Heintz T. M., 2021, MNRAS,
  * Evans et al. (2018) Evans D. W., et al., 2018, A&A, 616, A4
  * Everall & Boubert (2021) Everall A., Boubert D., 2021, arXiv e-prints, p. arXiv:2111.04127
  * Everall & Das (2020) Everall A., Das P., 2020, MNRAS, 493, 2042
  * Everall et al. (2019) Everall A., Evans N. W., Belokurov V., Schönrich R., 2019, MNRAS, 489, 910
  * Everall et al. (2021a) Everall A., Evans N. W., Belokurov V., Vasiliev E., Boubert D., Grand R., 2021a, in prep
  * Everall et al. (2021b) Everall A., Boubert D., Koposov S. E., Smith L., Holl B., 2021b, MNRAS, 502, 1908
  * Fantin et al. (2021) Fantin N. J., et al., 2021, arXiv e-prints, p. arXiv:2103.14721
  * Flynn et al. (2006) Flynn C., Holmberg J., Portinari L., Fuchs B., Jahreiß H., 2006, MNRAS, 372, 1149
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fukushima et al. (2019) Fukushima T., et al., 2019, PASJ, 71, 72
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2021a) Gaia Collaboration et al., 2021a, A&A, 649, A1
  * Gaia Collaboration et al. (2021b) Gaia Collaboration et al., 2021b, A&A, 649, A6
  * Gould et al. (1996) Gould A., Bahcall J. N., Flynn C., 1996, ApJ, 465, 759
  * Grand et al. (2018) Grand R. J. J., et al., 2018, MNRAS, 481, 1726
  * Grand et al. (2020) Grand R. J. J., et al., 2020, MNRAS, 497, 1603
  * Green (2018) Green G., 2018, The Journal of Open Source Software, 3, 695
  * Green et al. (2019) Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D., 2019, ApJ, 887, 93
  * Groenewegen (2021) Groenewegen M., 2021, arXiv e-prints, p. arXiv:2106.08128
  * Gunn et al. (1998) Gunn J. E., et al., 1998, AJ, 116, 3040
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hernitschek et al. (2018) Hernitschek N., et al., 2018, ApJ, 859, 31
  * Herschel (1785) Herschel W., 1785, Philosophical Transactions of the Royal Society of London Series I, 75, 213
  * Huang et al. (2021) Huang Y., Yuan H., Beers T. C., Zhang H., 2021, ApJ, 910, L5
  * Iorio & Belokurov (2019) Iorio G., Belokurov V., 2019, MNRAS, 482, 3868
  * Iorio et al. (2018) Iorio G., Belokurov V., Erkal D., Koposov S. E., Nipoti C., Fraternali F., 2018, MNRAS, 474, 2142
  * Joshi (2007) Joshi Y. C., 2007, MNRAS, 378, 768
  * Jurić et al. (2008) Jurić M., et al., 2008, ApJ, 673, 864
  * Kapteyn (1922) Kapteyn J. C., 1922, ApJ, 55, 302
  * Katz et al. (2021) Katz D., Gomez A., Haywood M., Snaith O., Di Matteo P., 2021, arXiv e-prints, p. arXiv:2102.02082
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kuijken & Gilmore (1989) Kuijken K., Gilmore G., 1989, MNRAS, 239, 605
  * Lang (2014) Lang D., 2014, AJ, 147, 108
  * Lindegren et al. (2021a) Lindegren L., et al., 2021a, A&A, 649, A2
  * Lindegren et al. (2021b) Lindegren L., et al., 2021b, A&A, 649, A4
  * Luri et al. (2018) Luri X., et al., 2018, A&A, 616, A9
  * Mackereth et al. (2017) Mackereth J. T., et al., 2017, MNRAS, 471, 3057
  * Mateu & Vivas (2018) Mateu C., Vivas A. K., 2018, MNRAS, 479, 211
  * McConnachie (2012) McConnachie A. W., 2012, AJ, 144, 4
  * Naidu et al. (2021) Naidu R. P., et al., 2021, arXiv e-prints, p. arXiv:2103.03251
  * Newberg & Yanny (2006) Newberg H. J., Yanny B., 2006, in Journal of Physics Conference Series. pp 195–204 (arXiv:astro-ph/0507671), doi:10.1088/1742-6596/47/1/024
  * Penoyre et al. (2020) Penoyre Z., Belokurov V., Wyn Evans N., Everall A., Koposov S. E., 2020, MNRAS, 495, 321
  * Perryman et al. (1997) Perryman M. A. C., et al., 1997, A&A, 500, 501
  * Read (2014) Read J. I., 2014, Journal of Physics G Nuclear Physics, 41, 063101
  * Ren et al. (2021) Ren F., Chen X., Zhang H., de Grijs R., Deng L., Huang Y., 2021, ApJ, 911, L20
  * Riello et al. (2021) Riello M., et al., 2021, A&A, 649, A3
  * Riess et al. (2021) Riess A. G., Casertano S., Yuan W., Bowers J. B., Macri L., Zinn J. C., Scolnic D., 2021, ApJ, 908, L6
  * Rix et al. (2021) Rix H.-W., et al., 2021, arXiv e-prints, p. arXiv:2106.07653
  * Robin et al. (2003) Robin A. C., Reylé C., Derrière S., Picaud S., 2003, A&A, 409, 523
  * Ruiz-Lara et al. (2020) Ruiz-Lara T., Gallart C., Bernard E. J., Cassisi S., 2020, Nature Astronomy, 4, 965
  * Rybizki et al. (2021) Rybizki J., Green G., Rix H.-W., Demleitner M., Zari E., Udalski A., Smart R. L., Gould A., 2021, arXiv e-prints, p. arXiv:2101.11641
  * Saha (1985) Saha A., 1985, ApJ, 289, 310
  * Shu et al. (2019) Shu Y., Koposov S. E., Evans N. W., Belokurov V., McMahon R. G., Auger M. W., Lemon C. A., 2019, MNRAS, 489, 4741
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Smith et al. (2009) Smith M. C., et al., 2009, MNRAS, 399, 1223
  * Snaith et al. (2015) Snaith O., Haywood M., Di Matteo P., Lehnert M. D., Combes F., Katz D., Gómez A., 2015, A&A, 578, A87
  * Stassun & Torres (2021) Stassun K. G., Torres G., 2021, ApJ, 907, L33
  * Thomas et al. (2018) Thomas G. F., et al., 2018, MNRAS, 481, 5223
  * Vasiliev & Baumgardt (2021) Vasiliev E., Baumgardt H., 2021, arXiv e-prints, p. arXiv:2102.09568
  * Vasiliev et al. (2021) Vasiliev E., Belokurov V., Erkal D., 2021, MNRAS, 501, 2279
  * Watkins et al. (2009) Watkins L. L., et al., 2009, MNRAS, 398, 1757
  * Widmark & Monari (2019) Widmark A., Monari G., 2019, MNRAS, 482, 262
  * Widrow et al. (2012) Widrow L. M., Gardner S., Yanny B., Dodelson S., Chen H.-Y., 2012, ApJ, 750, L41
  * Yanny et al. (2000) Yanny B., et al., 2000, ApJ, 540, 825
  * Zhu et al. (1997) Zhu C., Byrd R. H., Lu P., Nocedal J., 1997, ACM Trans. Math. Softw., 23, 550–560
  * Zinn (2021) Zinn J. C., 2021, AJ, 161, 214
  * de Jong et al. (2010) de Jong J. T. A., Yanny B., Rix H.-W., Dolphin A. E., Martin N. F., Beers T. C., 2010, ApJ, 714, 663

##  Appendix A Surface Density Integrals

Here we show the maths for evaluating the surface density of each component of
our model from the number density of sources.

###  A.1 Discs

The number density of sources in the each disc component is

$\nu_{c}=\frac{w_{c}\,\tan^{2}(b_{\mathrm{min}})}{2\pi\,h_{c}^{3}}\exp\left(-\frac{|z|}{h_{c}}\right)$
(22)

as given in Eq. 23 of Paper I where subscript $c$ refers to the disc
component, such that

$\int_{0}^{2\pi}\int_{b_{\mathrm{min}}}^{\pi/2}\int_{0}^{\infty}\nu_{c}\,s^{2}\,\cos(b)\,\mathrm{d}s\,\mathrm{d}b\,\mathrm{d}l=w_{c},$
(23)

is the total number of sources in the cone. The surface number density is

$\displaystyle\Sigma_{c}$
$\displaystyle=\int_{-\infty}^{\infty}\,\nu_{c}\,\mathrm{d}z$ (24)
$\displaystyle=2\int_{0}^{\infty}\,\frac{w_{c}\,\tan^{2}(b_{\mathrm{min}})}{2\pi\,h_{c}^{3}}\exp\left(-\frac{z}{h_{c}}\right)\,\mathrm{d}z$
(25) $\displaystyle=\frac{w_{c}\,\tan^{2}(b_{\mathrm{min}})}{\pi\,h_{c}^{2}}.$
(26)

We can also write this relative to the local number density of stars

$\nu_{c}(s=0)=\frac{w_{c}\,\tan^{2}(b_{\mathrm{min}})}{2\pi\,h_{c}^{3}}$ (27)

such that

$\Sigma_{c}=2\,h_{c}\,\nu_{c}(s=0).$ (28)

###  A.2 Halo

The number density of sources in the halo is

$\nu_{\mathrm{H}}=\mathcal{N}_{\mathrm{H}}\,r^{-n_{\mathrm{H}}}=\mathcal{N}_{\mathrm{H}}\,\left(R_{0}^{2}+z^{2}\right)^{-n_{\mathrm{H}}/2}$
(29)

where

$\frac{2\pi\,\mathcal{N}_{\mathrm{H}}}{R_{0}^{n_{\mathrm{H}}-3}}\mathcal{I}=w_{\mathrm{H}}$
(30)

and $\mathcal{I}$ is the dimensionless integral

$\mathcal{I}=\int_{\sin\left(b_{\mathrm{min}}\right)}^{1}\int_{0}^{\frac{s_{\mathrm{max}}}{R_{0}}}\,x^{2}\,\left(1+\sin^{2}(b)\,x^{2}\right)^{-n_{\mathrm{H}}/2}\,\mathrm{d}x\,\mathrm{d}\sin(b).$
(31)

$s_{\mathrm{max}}$ is the upper limit of the distances used in the model as
discussed in Section 4.

The surface density integral is

$\Sigma_{\mathrm{H}}=2\mathcal{N}_{\mathrm{H}}\,\int_{0}^{s_{\mathrm{max}}}\left(R_{0}^{2}+z^{2}\right)^{-n_{\mathrm{H}}/2}\,\mathrm{d}z.$
(32)

where we have made the approximation that the distance upper bound,
$s_{\mathrm{max}}$ approximately corresponds to an upper bound on $z$ due to
the narrow cone we have used. Since the number density of sources at $160$kpc
is far smaller than in the solar neighbourhood, the impact of this assumption
is negligible.

The surface density is then

$\displaystyle\Sigma_{\mathrm{H}}$
$\displaystyle=\frac{2\mathcal{N}_{\mathrm{H}}}{R_{0}^{n-1}}\int_{0}^{\frac{s_{\mathrm{max}}}{R_{0}}}\left(1+x^{2}\right)^{-n_{\mathrm{H}}/2}\,\mathrm{d}x$
(33)
$\displaystyle=\frac{w_{\mathrm{H}}}{\pi\,\mathcal{I}\,R_{0}^{2}}\mathcal{I}_{z}$
(34)

where

$\mathcal{I}_{z}=\int_{0}^{\frac{s_{\mathrm{max}}}{R_{0}}}\,\left(1+x^{2}\right)^{-n_{\mathrm{H}}/2}\mathrm{d}x.$
(35)

Again, we can write this in terms of the local source density

$\nu_{\mathrm{H}}(s=0)=\frac{w_{\mathrm{H}}\,R_{0}^{-3}}{2\pi\,\mathcal{I}}$
(36)

such that

$\Sigma_{\mathrm{H}}=2R_{0}\,\mathcal{I}_{z}\,\nu_{\mathrm{H}}(s=0)$ (37)