Modelling the Milky Way. I -- Method and first results fitting the thick disk and halo with DES-Y3 data

TL;DR

This paper introduces MWFitting, a new method for modeling the Milky Way's stellar components using DES-Y3 data, successfully recovering known structures and revealing a break in the halo density profile.

Contribution

The paper presents MWFitting, a novel technique for fitting Galactic components with Hess Diagrams, validated on synthetic data, and applied to DES-Y3 data to analyze the Milky Way's structure.

Findings

Recovered known stellar over-densities in DES data.

Detected a break in the Milky Way halo density profile.

Validated method achieves better than 3% precision on synthetic data.

Abstract

We present MWFitting, a method to fit the stellar components of the Galaxy by comparing Hess Diagrams (HDs) from TRILEGAL models to real data. We apply MWFitting to photometric data from the first three years of the Dark Energy Survey (DES). After removing regions containing known resolved stellar systems such as globular clusters, dwarf galaxies, nearby galaxies, the Large Magellanic Cloud and the Sagittarius Stream, our main sample spans a total area of $\sim$2,300 deg$^2$ distributed across the DES footprint. We further explore a smaller subset ($\sim$ 1,300 deg$^2$) that excludes all regions with known stellar streams and stellar overdensities. Validation tests on synthetic data possessing similar properties to the DES data show that the method is able to recover input parameters with a precision better than 3\%. Based on the best-fit models, we create simulated stellar catalogues covering the whole DES footprint down to $g = 24$ magnitude. Comparisons of data and simulations provide evidence for a break in the power law index describing the stellar density of the Milky Way (MW) halo. Several previously discovered stellar over-densities are recovered in the residual stellar density map, showing the reliability of MWFitting in determining the Galactic components. Simulations made with the best-fitting parameters are a promising way to predict MW star counts for surveys such as LSST and Euclid.

Tables

Table 1. Table 1: The MW model adopted in this work includes the bulge (as a triaxial truncated spheroid component), the thin disk (as an exponential model in the radial direction and a squared hyperbolic secant model in the vertical direction), the thick disk (as an exponential model in both R 𝑅 R and z 𝑧 z directions), and the halo (modeled with a double power-law profile). The columns list: the formula describing each MW component (first), free parameters (second), a description of each component (third), units (forth), initial value (fifth), and the best-fit value with errors (last column) for both samples ( † for raw sample and ‡ ‡ \ddagger for refined sample).

Formula	Symbol	Meaning	Unit	Initial value	Fixed/Best-fit value⋆
Bulge1
${\rho }^{\mathrm{bulge}}={\rho }_{GC}^{\mathrm{bulge}}\frac{\exp \left(-a^2/a_m^2\right)}{{\left(1+a/a_0\right)}^{1.8}}$ $\text{with }{\rho }^{\mathrm{bulge}}(0,0,0)={\rho }_{\mathrm{GC}}^{\mathrm{bulge}}$ $\text{with }a={\left(x^{\prime 2}+y^{\prime 2}/{\eta }^2+z^2/{\zeta }^2\right)}^{1/2}$ $\text{and }{x}^{\prime },y^{\prime }\text{ rotated by }{\varphi }_0\text{. w.r.t. }x,y$
	${\rho }_{\mathrm{GC}}^{\mathrm{bulge}}$	space density at GC	M⊙ pc-3	400	fixed
	$a_m$	scale length	pc	2500	fixed
	$a_0$	truncation scale length	pc	95	fixed
	$\eta$, $\zeta$	1:$\eta$:$\zeta$ scale ratios	-	0.68, 0.31	fixed
	${\varphi }_0$	angle w.r.t. Sun-GC line	deg (∘)	15	fixed
Thin disk
${\rho }^{\mathrm{thin}}=A^{\mathrm{thin}}{\text{sech}}^2(h/h_z^{\mathrm{thin}})\times$ $\exp \left(R/h_R^{\mathrm{thin}}\right)$ $\text{with }{h}_z^{\mathrm{thin}}=h_{z,0}^{\mathrm{thin}}+{\left(1+t/t_{\mathrm{incr}}^{\mathrm{thin}}\right)}^{\alpha }$ ${\text{and }{\int }_{h=-\infty }^{+\infty }{\rho }^{\mathrm{thin}}dz|}_{\odot }={\Sigma }_{\odot }^{\mathrm{thin}}$	${\Sigma }_{\odot }^{\mathrm{thin}}$	local mass surface density	M⊙ pc-2	55.412	fixed
	$h_R^{\mathrm{thin}}$	thin disk scale length	pc	29132	fixed
	$R_{\max }^{\mathrm{thin}}$	maximum radius	kpc	15	fixed
	$h_{z,0}^{\mathrm{thin}}$	scale height for	pc	94.72	fixed
		youngest stars
	$t_{\mathrm{incr}}^{\mathrm{thin}}$	timescale for increase in $h_z$	Gyr	5.552	fixed
	$\alpha$	exponent for increase in $h_z$	-	1.673	fixed
Thick disk
${\rho }^{\mathrm{thick}}=A^{\mathrm{thick}}\exp (h/h_z^{thick})\times$ $\exp \left(R/h_R^{\mathrm{thick}}\right)$ ${\text{with }{\int }_{h=-\infty }^{+\infty }{\rho }^{\mathrm{thick}}dz|}_{\odot }={\Sigma }_{\odot }^{\mathrm{thick}}$	$h_z^{\mathrm{thick}}$	scale height	pc	925	$925\pm {40}^{†}$
					$910\pm {46}^{‡}$
	$h_R^{\mathrm{thick}}$	thick disk scale length	pc	2667	$2667\pm {95}^{†}$
					$2631\pm {121}^{‡}$
	${\Sigma }_{\odot }^{\mathrm{thick}}$	local mass surface density	${10}^{-3}$M⊙ pc-2	3.89	$3.89\pm {0.65}^{†}$
					$3.97\pm {0.74}^{‡}$
	$R_{\max }^{\mathrm{thick}}$	maximum radius (fixed)	kpc	15	fixed
Halo
${\rho }^{\mathrm{halo}}=f{\rho }_{\odot }^{halo}{\left(\frac{r_{\odot }}{r_{obl}}\right)}^n$ $\text{with }{\rho }^{halo}(R_{\odot },0,z_{\odot })={\rho }_{\odot }^{halo},$ $r_{obl}=\sqrt{R^2+\left(z/q^2\right)}$ $\text{if }{r}_{obl}\le r_{br},n=n_1,f=1$ $\text{else }\left(r_{obl}>r_{br}\right),n=n_2,f={\left(r_{\odot }/r_{br}\right)}^{n_1-n_2}$	$n_1$	inner exponent	-	1.82	$1.82\pm {0.08}^{†}$
					$1.86\pm {0.11}^{‡}$
	$n_2$	outer exponent	-	4.14	$4.14\pm {0.05}^{†}$
					$4.24\pm {0.06}^{‡}$
	$r_{BR}$	break radius	kpc	18.52	$18.52\pm {0.27}^{†}$
					$18.59\pm {0.49}^{‡}$
	$q$	axial ratio $z/x$	-	0.75	$0.75\pm {0.01}^{†}$
		(oblateness)			$0.74\pm {0.02}^{‡}$
	${\rho }_{\odot }^{halo}$	local mass space density	${10}^{-5}\text{M}\text{<sub><span>⊙</span></sub>}$pc-3	3.31	$3.31\pm {0.20}^{†}$
					$3.51\pm {0.26}^{‡}$
Dust layer
${\rho }^{\mathrm{dust}}=A^{\mathrm{dust}}\exp \left(h/h_z^{\mathrm{dust}}\right)$ $\text{with }{\int }_{\ell =0}^{+\infty }{\rho }^{\mathrm{dust}}d\ell =A_V^{\infty }$
	$A_V^{\infty }$	total extinction at infinity	-	4	fixed
	$h_z^{\mathrm{dust}}$	dust scale height	pc	1105	fixed
Others
	$R_{\odot }$	Sun’s distance to the GC	kpc	8.1226	fixed
	$z_{\odot }$	Sun’s height above plane	pc	20.87	fixed
1 Parameters from Vanhollebeke et al. (2009)
2 Best-fit parameter from Girardi et al. (2005)
3 Adopted in Girardi et al. (2005)
4 Schlegel et al. (1998)
5 Lynga (1982)
6 Gravity Collaboration et al. (2018)
7 Bennett & Bovy (2019)
⋆ See Table 3 for more details about those parameters.

Table 2. Table 2: Results of two tests (A and B) using mwfitting  . Even though the initial guesses start far from the input values, the final parameter values are within 3% of the input values. The simulations in this table compare 100 fields and oversample the models in the same way as the comparison to real data.

Parameter	Unit	True	Initial Guess		Best-fitting		$\frac{\mid Best-True\mid }{True}$ (%)		$\frac{Best-True}{\sigma }$
		Value	A	B	A	B	A	B	A	B
ThickDisk $h_z$	$\mathrm{pc}$	925	1037.6	903.6	${923.1}_{-1.9}^{+2.3}$	${922.5}_{-1.9}^{+1.9}$	0.2	0.2	-0.8	-1.2
ThickDisk $R_e$	$\mathrm{pc}$	2666	2849	2397	${2657}_{-6}^{+6}$	${2675}_{-6}^{+6}$	0.3	0.3	-1.5	+1.5
ThickDisk $\rho$ (R=R0)	$\times {10}^{-3}$M${}_{\odot }{}^{}\mathrm{pc}_{}^{-2}$	3.90	4.23	3.57	${3.91}_{-0.03}^{+0.02}$	${3.94}_{-0.02}^{+0.02}$	0.3	1.0	+0.4	+1.7
Halo $n_1$	-	1.821	1.501	1.991	${1.861}_{-0.011}^{+0.010}$	${1.867}_{-0.010}^{+0.009}$	2.2	2.5	+3.8	+4.8
Halo $n_2$	-	4.139	4.407	4.520	${4.133}_{-0.005}^{+0.005}$	${4.124}_{-0.005}^{+0.005}$	0.1	0.4	-1.2	-3.0
Halo $r_{br}$	kpc	18.52	17.61	20.06	${18.65}_{-0.04}^{+0.04}$	${18.63}_{-0.05}^{+0.04}$	0.7	0.6	+3.0	+2.5
Halo $q$	-	0.748	0.785	0.827	${0.745}_{-0.001}^{+0.001}$	${0.748}_{-0.001}^{+0.001}$	0.4	0.0	-3.0	0.0
Halo $\rho$ (R=R0)	$\times {10}^{-5}$ M${}_{\odot }{}^{}\mathrm{pc}_{}^{-3}$	3.31	3.36	3.52	${3.40}_{-0.02}^{+0.02}$	${3.39}_{-0.02}^{+0.02}$	2.7	2.4	+4.5	+4.0

Table 3. Table 3: Best-fit parameters for the raw and refined samples. The last two columns are results from the literature. In our results, the first errors listed are the 1 ​ σ 1 𝜎 1\sigma statistical error or the standard deviation of the mean estimated by the jackknife block method (see more details in the text). The second errors are the systematic errors as discussed in Section  5 and   4.3 . They represent the ability of the pipeline to recover the true model, and the degeneracy of the parameters regarding the uncertainty of the local density of the thin disk.

Parameter	Unit	mwfitting		Juri$\acute{c}$	de Jong	Deason
		Raw sample	Refined sample	et al. 2008	et al. 2010	et al. 2011
ThickDisk $h_e$	pc	$925\pm 6\pm 40$	$910\pm 8\pm 45$	$743\pm 150$	$750\pm 70$	-
ThickDisk $R_e$	pc	$2667\pm 89\pm 34$	$2631\pm 111\pm 49$	$3261\pm 650$	$4100\pm 400$	-
ThickDisk $\rho$ (R=R0)	${10}^{-3}$ M⊙pc-2	$3.89\pm 0.09\pm 0.64$	$3.97\pm 0.12\pm 0.73$	$7.53\pm 0.75$	$5.01\pm 1.30$	-
Halo $n_1$	-	$1.82\pm 0.05\pm 0.06$	$1.86\pm 0.07\pm 0.08$	-	-	${2.3}_{-0.1}^{+0.1}$
Halo $n_2$	-	$4.14\pm 0.03\pm 0.04$	$4.24\pm 0.04\pm 0.05$	-	-	${4.6}_{-0.1}^{+0.2}$
Halo $r_{br}$	kpc	$18.52\pm 0.15\pm 0.23$	$18.59\pm 0.39\pm 0.29$	-	-	${27}_{-1}^{+1}$
Halo $q$	-	$0.75\pm 0.01\pm 0.01$	$0.74\pm 0.02\pm 0.01$	$0.64\pm 0.01$	$0.88\pm 0.03$	${0.59}_{-0.03}^{+0.02}$
Halo $\rho$ (R=R0)	${10}^{-5}\text{M}\text{<sub><span>⊙</span></sub>}p{c}^{-3}$	$3.31\pm 0.10\pm 0.17$	$3.51\pm 0.13\pm 0.23$	$2.95\pm 0.74$	$6.31\pm 0.77$	-

Table 4. Table 4: Stellar masses estimates for the MW components fit in this work, for the raw and refined samples.

Component	Estimated mass (M⊙)
	Raw sample	Refined sample
Thick disk	$3.57\pm 0.43\times {10}^6$	$3.70\pm 0.44\times {10}^6$
Halo (r100kpc)	$6.80\pm 1.04\times {10}^8$	$6.98\pm 1.56\times {10}^8$

Table 5. Table 5: The main parameters to run mwfitting   pipeline.

Parameter name	Description	Std. configuration	Unit
Input data
NSIDE	Footprint map granularity	1024	HEALPix Nside
Input as simulation	Build and fit mock catalogues	False	-
Random factor	Range to multiply initial values	2	-
	(for ‘input as simulation’ method)
Seed for input	Seed to generate random numbers	0	-
	and multiply random factor
Mock catalogues
Build mock catalogue?	Build a mock catalogue	False	-
	following the best-fitting parameters
Hess Diagrams
Nside	Size of cells for HD	16	HEALPix Nside
Minimum area	minimum coverage of the cell	0.9	Cell area
Global seed	Global seed for fields choice	1	-
Magnitude range	min, max, step	17.0;21.0;0.1	mag
Colour range	min, max, step	0.0;0.8;0.1	mag
Cell counts	Number of fields to be fitted	10	Field
Filters
Streams	Filter streams	True	-
Minimum Galactic latitude	Lowest absolute value for $b$	30	degree
Optimise
Overfactor	Oversampling the models	20	area in cell

Equations

$$-2\ln \lambda (\theta )=2\sum _{i=1}^N\left({\nu }_i(\theta )-n_i+n_i\ln \frac{n_i}{{\nu }_i(\theta )}\right)$$

$$p{s}_i(X)=n_b{\varphi }_n(X_1,...,X_n)-(n_b-1){\varphi }_{n-k}((X_1,...,X_n{)}_{[i]})$$

Full text

††pagerange: Modelling the Milky Way. I - Method and first results fitting the thick disk and halo with DES-Y3 data–Affiliations††pubyear: 2019

1 Introduction

Over the last 40 years, we have learned the utility of describing a complex system like the Milky Way (MW) through simple building blocks (e.g., Bahcall & Soneira, 1981), composed of nearly homogeneous stellar populations, smoothly distributed in space in a few components like the thin and thick disks, bulge and halo. The derivation of simple parameters for these components – such as scale lengths and heights, limiting radii, central densities, etc. – allows us to put our Galaxy in perspective by comparing it to other spiral galaxies (Courteau et al., 2011) and to galaxies produced in cosmological simulations (see, e.g. Hopkins et al., 2014; Bland-Hawthorn & Gerhard, 2016). Examining the residuals of the best-fit models enables the identification of stellar substructure such as dwarf galaxies and stellar streams (e.g., Shipp et al., 2018). Fitted models can also be used to estimate the distribution of stars in future surveys.

Our understanding of the MW has steadily advanced over the past several decades. For example, the thick disk (Gilmore & Reid, 1983) has long been proposed to explain the MW stellar population within 1-5 kpc on either side of the Galactic plane. Thick disk stars differ from those closer to the Galactic plane in kinematics, age and metalicity, being older, more metal-poor, less rotationally supported, and having typically higher [$\alpha$/Fe] at a fixed metalicity (for instance, see Reddy et al., 2006; Fuhrmann, 2008). More recently, the spatial structure of different stellar populations has been studied by Anders et al. (2014) and Bovy et al. (2016), among others, using survey data from APOGEE (Majewski et al., 2016). In brief, high [$\alpha$/Fe] stars tend to follow a double exponential density profile parallel and perpendicular to the Galactic plane, with scales of $h_{R}\simeq 2.2$ kpc and $h_{z}\simeq 1.0$ kpc, respectively (Bovy et al., 2016). The lower [$\alpha$/Fe] stars display a more complex distribution, including a metalicity gradient and disk flaring (Anders et al., 2014). Even so, the traditional description of the thin and thick disk components with double exponential profiles (or a sech${}^{2}z$ perpendicular to the disk plane) is adequate (Cabrera-Lavers et al., 2005; Jurić et al., 2008; de Jong et al., 2010).

At the outer limit of the MW, the stellar Galactic halo is roughly spherical in shape. Early studies indicated that the radial density of this component is better described by a power-law profile with index $n\sim-2.75$ than an exponential profile (Jurić et al., 2008; de Jong et al., 2010). However, more recent work has found that the stellar density drops off faster at typical distances $\simeq$ 20 kpc, suggesting that the density of the stellar halo follows a broken power-law profile (Watkins et al., 2009; Deason et al., 2011; Sesar et al., 2011; Deason et al., 2018) or another model that decreases more rapidly at large radii (Einasto, 1965; Merritt et al., 2006; Deason et al., 2011; Hernitschek, et al., 2018). These observations are not unexpected, since a power-law index of $n<-3$ is necessary at large radii in order for the integrated mass of the stellar halo to converge.

In addition to the aforementioned developments in describing the stellar content of the Galaxy, an impressive amount of work has been dedicated to determine the star formation rate (SFR, Ryan & Norris, 1991; Fuhrmann, 1998), initial mass function (IMF, Kroupa, 2001; Chabrier, 2003; Kroupa & Weidner, 2003; Wood & Mao, 2005), and Age-Metalicity Relation (AMR, Rocha-Pinto et al., 2000; Zoccali et al., 2003; Fuhrmann, 2008; Casagrande et al., 2011) for the stars in the MW, along with the modelling of stellar evolution (Bertelli et al., 1994; Girardi et al., 2000, 2002; VandenBerg et al., 2006; Marigo & Girardi, 2007; Girardi et al., 2010; Paxton et al., 2011; Spada et al., 2013) and the stellar contents of the Galaxy itself (Sharma et al. 2011; Czekaj et al. 2014111See https://model.obs-besancon.fr/modele_ref.php for a complete list of publications of the Besançon group.; Pasetto et al. 2018). Thanks to all these developments, we are now able to build a detailed structural model for the Galaxy.

To take advantage of this knowledge and the increasing number of deep wide-field astronomical surveys, we have developed mwfitting. This work aims to present the method and to show its first application to data in the Dark Energy Survey (DES; DES Collaboration 2005).

In this work we aim to:

• •

Present an efficient method to describe the structure of the Galaxy by comparing star counts to predictions of stellar population synthesis models. The comparison between data and models is made through binned colour-magnitude diagrams (i.e., Hess Diagram, HD) in specific regions in the sky. Many different models are used to predict star counts, such as the spatial distribution of stars in the MW components, the stellar IMF, SFR, and AMR. Also crucial in determining star counts are the input stellar evolutionary models that prescribe magnitudes and colours as a function of fundamental stellar parameters, such as mass, age, and metalicity.

• •

Validate the code using mock data. These tests are done to test the accuracy of mwfitting  to evaluate systematic uncertainties, and to measure the effect initial values has on recovering the input parameters.

• •

Apply mwfitting to model the Galactic thick disk and halo in DES three year (Y3) data.

• •

Show and discuss the results of the method and the implications on the Galactic model adopted.

This paper is structured as follows: in Section 2 we discuss the mwfitting. In Section 3 we briefly describe the DES Y3 data. In Section 4 we present the results of mwfitting. In Section 5 we describe a simulation based on the best fitting parameters and discussion of the results. Finally, we conclude in Section 6.

2 MWFitting package

In this paper, we adopt Trilegal222http://stev.oapd.inaf.it/cgi-bin/trilegal models to describe the stellar content of the Galaxy. Trilegal is a stellar population synthesis code, based on the Girardi et al. (2002) database of stellar isochrones, and augmented with models for brown and white dwarfs. For more details about the stellar models, we refer to Girardi et al. (2005). Note that even though several upgrades in the database of evolutionary tracks and stellar atmospheres have become available recently (see, e.g. Marigo et al., 2017), they severely reduce computational speed, and only include short-lived evolutionary phases and cool stars, which are not the subject of the present work.

The following subsections present the sequence of steps that leads to a final product of the mwfitting . Section 2.1 describes Trilegal input parameters to model a sky region with a specific Galactic model. The previous attempts to calibrate the Galactic model using Trilegal are briefly discussed in Section 2.2; the adopted Galactic model is presented in Section 2.3; in Section 2.4 we discuss the implementation of the mwfitting; in Section 2.5 we validate the mwfitting pipeline using synthetic data with known input parameters.

2.1 Trilegal parameters

The Trilegal population synthesis simulation requires input parameters such as: covered area, photometric system, magnitudes and colour ranges, 3D position of the Sun, dust distribution, IMF for single stars, binary fraction, and mass ratios of unresolved binaries. It also requires structural models, SFR, and AMR for each Galactic component (see Table 1).

Regarding the color and magnitude ranges, Trilegal models are very successful in describing the stellar evolutionary phases as the main sequence (MS), including the turn-off (MSTO), and stars in the sub and red giant branches (respectively, SGB and RGB), for stars in a wide range of masses.

Stellar evolutionary models present a poor colour-fit for low-mass stars with [Fe/H] $\geq$ -2, such as M-type stars, which is the most abundant spectral type in thin disk. See for instance (Sarajedini et al., 2007), for a discussion about the comparisons of simple stellar populations of globular clusters to theoretical models.

Based on that, we choose to exclude the red thin-disk stars (see figure 2 and discussion in de Jong et al. 2010) and keep the parameters of this component fixed. The magnitude depth of DES also favours stars farther away than those in the thin disk, which supports our choice.

2.2 Previous attempts to calibrate Trilegal

Early descriptions of the MW components and their calibrations using Trilegal are found in Groenewegen et al. (2002) and Girardi et al. (2005). Those first attempts were based on a simple trial-and-error approach, where each model parameter was set to literature values, changed by hand until a “good description” for the star counts was met for a given survey. Surveys used in these analyses compromise both deep (e.g., DMS and EIS-deep, Osmer et al., 1998; Arnouts et al., 2001), shallow (e.g., , Skrutskie et al., 2006) photometric data, and local (e.g., Hipparcos catalogue, Perryman et al., 1997).

Vanhollebeke et al. (2009) explored a different approach to calibrate the bulge’s parameters using Trilegal. They defined a likelihood function to quantitatively evaluate the goodness-of-fit between data and model (see also Eidelman et al., 2004; Dolphin, 2002) as:

[TABLE]

where $n_{i}$ is the number of observed objects in a given magnitude/colour bin i, and $\nu_{i}(\theta)$ is the number of objects predicted by the set of parameters $\theta$ that describes the model. The summation is performed over all lines-of-sight, and magnitude/colour bins included in the comparison. The authors used the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Fletcher, 1987) to maximize their likelihood and derived uncertainties from the likelihood profile, as detailed in that work.

In this context, the fitting of disk and halo parameters using the latter method requires an extra set of variables. This presents several issues:

• •

Fitting the disk (thin and thick) and halo implies the simultaneous fitting of $\sim 30$ structural parameters, with many samples across the sky. The resulting analysis is very time-consuming.

• •

Local maxima in likelihood space may be very common, and due to the high dimensionality of the problem, finding the absolute maximum may be challenging.

This is not the case when fitting the bulge, as there are fewer parameters, and there are a large set of lines-of-sight which leaves little chance for solutions to be trapped in local maxima (Vanhollebeke et al., 2009). In the present case, it is advisable to implement tests for local maxima in log-likelihood space, and check whether different starting conditions lead to the same solution. These tests imply even longer computing times.

In the following sections, we describe the implementation of mwfitting aimed at tackling the challenges discussed above (see also Girardi et al., 2012).

2.3 Galactic model adopted

Table 1 summarizes the functional form utilized for each Galactic component, the parameters that describe the component, and whether the parameter is fixed or free in the fit. We adopt an exponential model along the disk plane and a square hyperbolic secant perpendicular to it for the thin disk. The parameters of the thin disk and bulge modelled by Trilegal in this work are kept fixed at the values described in Girardi et al. (2005), with some minor tweaks as in Girardi et al. (2012). The only parameters allowed to vary are related to the thick disk and to the halo of our Galaxy. An exponential model in both radial and vertical directions describes the distribution of stars in the thick disk. The stellar halo is described by a double power-law profile, with an inner exponent, $n_{1}$, describing the stellar density of the halo out to a certain distance, $r_{BR}$ (radius of the break) and an outer exponent, $n_{2}$, for farther distances. We require that the density of halo stars is continuous at $r_{BR}$ for both exponents. Since DES covers the south Galactic cap ($b<-30^{\circ}$), it largely excludes the MW bulge. Therefore, in this analysis, we fix the parameters of the bulge component following the triaxial model presented in Binney et al. (1997).

The IMF assumed for Galactic stars is the Chabrier lognormal IMF (Chabrier, 2003) and the fraction of binaries adopted is 30%, with the mass ratio of the secondary over the primary limited between 0.7 and 1.0 (Barmina et al., 2002). The SFR and AMR are specific to each MW component. Stars in the bulge and in the thick disk follow a SFR and AMR described by Zoccali et al. (2003) and Boeche et al. (2013), respectively. Thin disk and halo stars are modelled following previous comparisons from Groenewegen et al. (2002) and Girardi et al. (2005).

2.4 The MWFitting pipeline: fitting the galaxy with Hess diagrams

The mwfitting pipeline code fits a global, multi-component model of the MW to the observed stellar density in bins of Galactic longitude and latitude, magnitude, and color. The inclusion of spatial and color-magnitude information allows us to break degeneracies between the various MW model components.

We begin by pixelizing the sky using the HEALPix333https://healpix.jpl.nasa.gov/ scheme to define individual lines-of-sight (which we call “cells”). We select cells that reside within the survey, and remove cells that are contaminated by resolved stellar populations such as globular clusters and dwarf galaxies. For each cell, we calculate the coordinates of the centre, the average reddening and reddening dispersion, the limiting magnitude (as specified by the user), the colour range, and the bin size in the color-magnitude diagram (CMD) space.

Within each cell, we calculate model HDs for each component (i.e., bulge, halo, thin, and thick disk) over a range of distance moduli, typically separated by 0.1 mag. These so-called “partial HDs” for each component and distance are stored in separate Header Data Units (HDUs) of a multi-extension FITS444https:////fits.gsfc.nasa.gov/fits_standard.html file. This data format allows the normalizations of different model components to be quickly adjusted. For example, the normalization of the stellar halo can be adjusted by a factor $f$, by multiplying all partial HDs associated with the halo by the same factor $f$. The total model-predicted MW HD can be quickly calculated from a linear combination of the individual partial HDs, incorporating typical photometric errors of the survey in each band. This method allows us to rapidly construct stellar density predictions for a wide range of MW model parameters as listed in Table 1. Variation in each parameter corresponds to varying the weight of each partial HD, which are then combined to produce a total HD in each HEALPix cell.

The Poisson log-likelihood (Eq. 1) is calculated by first comparing the total model-predicted HDs to the data in each cell and then summing the log-likelihoods over all cells. To fit the MW model to an observed data set, we apply an Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble sampler (i.e.,emcee, Foreman-Mackey et al., 2013). The free and fixed parameters of our model, along with their initial values, are listed in Table 1. We assume flat priors ranging from 0.5-2.0 times the initial value of each free model parameter. We also checked visually whether the walkers converged or not at the end of the burn-in phase, in order to inform realistic best-fitting parameters.

Since Trilegal computes a discrete distribution of points as a realization of the expected population of stars in each cell, we are left with statistical noise due to sampling a finite number of points. To mitigate this noise, we increase the number of simulated stars by an over-factor which is then taken into account while normalizing the final Hess diagram for each cell. A typical over-factor value is 30, for the magnitude, colour range, and MW components explored in this work.

The mwfitting code was developed and is currently implemented in the DES-Brazil Portal powered by Laboratório Interinstitucional de e-Astronomia (LIneA555http://www.linea.gov.br/). More details on the DES-Brazil Portal can be found in Gschwend et al. (2018) and Fausti Neto et al. (2018). The application of mwfitting  to the DES data took 23h in a SGI ICE-X FC3Y cluster with 4 compute nodes. Each node contained 48 cores and 125 GB of RAM. For more detailed or technical information, the reader is directed to the Appendix B, where the input parameters of the pipeline and details about them are described.

2.5 Validating the code with mock data

In this section we test mwfitting using mock data. We verify that we can recover the input parameters of our simulated data set when applied to an area with the same footprint as DES-Y3.

Each test utilizes 100 cells, and each cell has the same area as the unit cell designed for the real data (HEALPix pixels with nside=16), following identical footprint and coverage maps (see Section 3). The range in magnitude and colour is the same as the DES data ($17<g<23$ and $0.0<g-r<0.6$, respectively), with the same bin size in magnitude and color space (0.1 mag). Uncertainties in the magnitude of the stars were also incorporated in the synthetic data.

Table 2 lists the parameters, units, input values, best-fit values and their errors, as indicated by emcee, and the significance of the differences between the best-fit and the true value, for two trials. We run two tests with the same input parameters but different initial values for the MCMC, which we refer to as test A and B. Analyzing Table 2, we find that mwfitting is able to recover the input values of the mock data accurately, even when the initial starting points are far from the true ones. Differences between true and best-fit values are within $3\%$ of the true parameters at the maximum, and the deviations are within 3$\sigma$, with a few exceptions. The maximum differences occur for the density of the halo and its inner exponent, while the differences for the remaining parameters are all below 1%.

Inspecting the HDs, there is excellent concordance between the mock data and the best-fit model data. The overall range of differences in test A between input data and best-fit models is [-2.28%, +1.40%], in terms of star counts. Fig. 1 shows the HDs of the cell with the largest absolute difference (-2.28%), with the center located at [$l=226.41\degr$, $b=-69.02\degr$]. The panels of Fig. 1 shows the HD of the best-fit model, simulated input data (mock), absolute difference, and the Poissonian significance over the HD cells, limited by the maximum significance (given in the title of the panel). The distribution of differences and their significance values show no systematic trend in the colour-magnitude plane. Note that the best-fit HD is smoother than the mock HD distribution due to the oversampling of the model.

Test B produced similar results to test A, with star counts differences in the range [-2.25%, +2.19%]. The cell with the largest absolute difference (-2.25%) exhibits one bin in the HD diagram with maximum significance of 4.6$\sigma$. There is a general concordance in the remaining cells, with typical maximum significance $\leq 4\sigma$ in the cells of the HDs.

The differences between the recovered and true values (the last two columns of Table 2) are expected to follow a standard normal distribution, with $\mu$=0 and $\sigma$=1. However, those values appear to be somewhat higher than expected, reflecting a systematic error in the recovery of the true model greater than the uncertainty reported by emcee. In order to encompass half of the recovery errors within $\pm 0.67\sigma$ (or 50% of the area of the standard normal distribution), the uncertainties provided by emcee are increased by a factor of 6.0. In this way, we aim to incorporate realistic systematic errors, and we are assuming they are downsampled by emcee method.

3 DES Data

DES is a wide-area photometric survey covering about $5\,000$ deg2 in the southern Galactic cap (DES Collaboration, 2005). DES images were taken with the Dark Energy Camera (DECam, Flaugher et al., 2015), with a typical single-exposure (90s in $griz$ bands and 45s in $Y$ band) 10$\sigma$ limiting magnitudes of g = 23.57, r = 23.34, i = 22.78, $z$ = 22.10 and $Y$ = 20.69 for point sources (Morganson et al., 2018). The final coadded images at the end of the first 3 years of observations achieve g = 24.33, r = 24.08, i = 23.44, $z$ = 22.69 and $Y$ = 21.44 at 10$\sigma$ (DES Collaboration, 2018). DES was designed for cosmological analyses, avoiding the Galactic plane (DES Collaboration, 2018). Therefore, also considering the depth of the survey, the DES stellar sample will mostly contain stars from the Galactic thick disk and halo. In this section, we characterize the main aspects of the photometry and star/galaxy (S/G) separation in the DES.

DES-Y3 data were processed by the DES Data Management system (DESDM, Morganson et al., 2018) and include observations from the first three years of the survey. The DES catalogue studied here is the Year 3 Gold release version 2.2 (Sevilla-Noarbe, in preparation), hereafter referred to as the DES-Y3 catalogue. This catalogue is composed of the same objects as the first public data release (DES-DR1; see DES Collaboration, 2018), but contains enhanced photometric and morphological measurements and other ancillary information.

In order to identify the area covered by the DECam observations, the sky is partitioned in HEALPix pixels (nside=4096) with size equal to 52 arcsec $\times$ 52 arcsec (footprint map). Regions around stars brighter than $J=12$ in 2MASS (Skrutskie et al., 2006), globular clusters (Harris, 1996, updated 2010) and a small area close to Large Magellanic Cloud (LMC) were masked. The area covered by DECam in each band and pixel (coverage map) is also estimated by a coverage map produced from mangle (Swanson et al., 2008). The DES-Y3 catalogue lists objects located in pixels (with nside=4096) with sampled area $>50$% in g, r, i and z bands and imaged at least once in all those four filters.

The DES-Y3 Gold data is photometrically calibrated by the Forward Global Calibration Method (FGCM666https://github.com/lsst/fgcmcal, see Burke et al. 2018). A comparison between DES-Y3 and $Gaia$ DR1 (Lindegren et al., 2016) shows a mean difference of 0.0014 magnitudes with $\sigma=0.0067$ magnitudes (DES Collaboration, 2018). The PSF photometry for DES-Y3 catalogue was performed by simultaneously fitting each object in multiple exposures (single object fitting or ‘SOF’). This procedure is very similar to the multi-object PSF-fitting (‘MOF’) described in Drlica-Wagner et al. (2018).

Initially, we apply a S/G separation procedure that is similar to Shipp et al. (2018). We use the parameter EXTEND_CLASS_MASH_SOF, which is a variable designed to classify point source (star or quasi-stellar objects - QSO) or extended sources (galaxies) based on ngmix (Sheldon, 2015). We nominally adopt values from the SOF photometry and when SOF photometry is unavailable we adopt photometry from the coadded images. This criterion increases the stellar sample by including stars with good PSF-fitting in coadded images but with failures in SOF. This S/G separation is applied for objects in the full range of magnitudes. Similar to Shipp et al. (2018), the same weight-averaged SPREAD_MODEL in i band is applied as S/G classification for the sample of bright stars ($g<18$) where PSF photometry fails.

Extensive completeness assessments were carried out in the DES year 1 (DES-Y1) catalogue, assuring that the catalogue is virtually complete in the range $17<g<22$, with estimated completeness $\geq$ 95% at the faint limit (Sevilla-Noarbe et al., 2018).

The quality of the DES photometry and S/G classification is illustrated in Fig. 2, where we show a colour-colour diagram ($g-r$ vs $r-i$) for sources classified as stars and corrected for reddening following Schlegel et al. (1998). There are 13,990,013 sources within the magnitude range $17<g<23$ and the limits shown in Fig. 2, namely $0.0<g_{0}-r_{0}<1.6$ and $-0.3<r_{0}-i_{0}<1.6$. A blue plume close to $g_{0}-r_{0}\cong 0$ and $r_{0}-i_{0}\cong 0.25$ amounts to a few thousands stars, probably due to binary systems with a white dwarf and a main sequence star (Kleinman et al., 2004). A lower level of contamination by QSO’s is expected in that region of the color-color diagram ($0.0\leq g-r\leq 0.5$ and $-0.25\leq r-i\leq 0.50$), along with contamination on the redder end, which is not taken into account in the process of fitting.

In order to further decrease contamination from mis-classified galaxies, we tested alternative methods for star-galaxy separation. The best method that we found was to use the photometric redshift as a way to identify galaxies that were morphologically classified as stars. Photometric redshifts were estimated using the Directional Neighbourhood Fitting or DNF (De Vicente, Sánchez & Sevilla-Noarbe, 2016), and we refer to this work for details about the fitting of the redshift. We removed objects with DNF photo-z $z>0.55$.

To assess the stellar completeness of DES at faint magnitudes, we matched the DES stars to the SPLASH-SXDF catalogue (Mehta, et al.2019, ), using as reference for the S/G classification the tag STAR_FLAG, based on the $BzK$ color-color-diagram. The comparison between DES data and SPLASH-SXDF show DES data is $>$90% complete down to $g=23$. This confirms the estimates in Shipp et al. (2018), and we expect that this sample will have minimal contamination from galaxies and QSOs.

4 MWFitting applied to DES-Y3 stars

We partition the DES data into cells corresponding to HEALPix pixels with nside=16, covering a solid angle of $13.43$ deg2. The cells included in the analysis are those with a fill factor $\geq 80$% ($>10.74$ deg2) of its footprint. Such criterion (and others mentioned below) are identical to those adopted for the validation tests.

We choose a constant range of magnitude ($17<g<23$) and color ($0.0<g-r<0.6$) when applying mwfitting  to DES data. This constant color-magnitude selection is motivated by the uniformity of the DES footprint in this magnitude depth, and by the high confidence of the modelled stars in this color range. We bin the data in color-magnitude space with a bin size of 0.1 mag in both color and magnitude. This choice of bin size is somewhat arbitrary, and we have found that the results of our analysis are insensitive to the choice of bin size.

The stars in our sample are not reddening corrected, instead the reddening is incorporated in the models following a Gaussian distribution based on the average and dispersion of the reddening on each cell.

We exclude cells with known stellar clusters and dwarf galaxies. The list of objects includes globular clusters and dwarf galaxies discovered up-to-date (Harris 1996, 2010 edition; McConnachie 2012; Drlica-Wagner et al. 2015; Koposov et al. 2015; Kim & Jerjen 2015; Luque et al. 2018), along with nearby galaxies partially resolved into stars in the DES images and catalogues (IC5152, ESO294-G010, NGC55, NGC300, NGC1399, NGC247, IC1613, ESO410-G005). The stars from those objects represent a potential contamination to Galactic fields and these fields contained positive residuals in initial iterations of mwfitting.

Cells with any region closer than 22*∘* from the LMC centre were also masked. Nidever et al. (2019) clearly shows (see their figure 5) a significant population of LMC main-sequence stars located out to 21*∘* from the center of the LMC. Furthermore, we masked the Sagittarius Stream, removing a stripe of width equal to 20*∘* along the centre of the stream (Majewski et al., 2003).

After removing the aforementioned regions and selecting only cells with a fill factor of more than 80%, the remaining 194 cells constitute our so-called raw sample. This sample includes the stellar population of streams discovered in the DES footprint (Shipp et al., 2018) and the Eridanus-Phenix overdensity (Eri-Phe, Li et al., 2016). Since these objects cover a large area with a much lower stellar density than that of the Galaxy, we retain them in the raw sample.

4.1 With or without streams?

To explore the influence of including regions with known stellar streams and the Eri-Phe overdensity, we define a second sample removing the regions where those objects are located. The list of masked stellar streams is that described by Mateu (2017), and we refer to this work for further details. In the case of Eri-Phe over-density, the masked area has a triangular shape as shown in figure 3 of the discovery paper (Li et al., 2016). The second sample of DES data comprise 105 cells, and we refer to this sample as the refined sample.

Fig. 3 puts into perspective the footprint of the raw and refined samples using an orthonormal projection of the southern Galactic Hemisphere. The DES footprint is outlined in black. The cells included in mwfitting are displayed in green and masked cells are shown in orange. The raw and refined samples are top and bottom respectively. A significant portion of the DES footprint is masked in the refined sample.

The Sagittarius Stream stands out in both panels of Fig. 3 as a wide stripe crossing the South Galactic Pole and cells masked due to proximity to the LMC are in the lower left corner. The area sampled by DES-Y3 and compared to models amounts to 2,315 deg2 (194 cells) in the raw sample, and to 1,256 deg2 (105 cells) in the refined sample.

4.2 MWFitting configuration and errors

Before discussing the outcomes from applying mwfitting to DES data, we first discuss the emcee configuration used. We use 200 walkers along 250 steps with step length as 1% of each parameter to sample the posterior distribution. We perform initial iteration, starting with input values from the literature. In a second iteration, we redo the fit starting from previous fitting. The first 200 steps are discarded as a burn-in phase, and we examine the remaining distribution to check that the walkers have converged. We apply a Gelman-Rubin convergence diagnostic ($R_{c}\leq 1.004$) to check for convergence of the Markov chains.

The results from applying mwfitting to the raw and refined samples are listed in Table 3. We find that the errors reported from the posterior distribution are smaller than the difference of best-fit parameters when we tested the pipeline with subsets of the raw or refined sample.

Hence, we have decide to estimate the statistical errors from a jackknife resampling method (Feigelson and Babu, 2012), in addition to the systematic errors based on the emcee method. The jackknife method creates $n$ samples (where $n$ is the number of observations), replicating the initial sample in each iteration, but omitting the i-$th$ observation. The jackknife block method is similar, but instead we group the observations into $n_{b}$ data blocks with size $k$ (in our case, the blocks are a set of cells). In each $i$ subsample with $k$ size, a pseudo-value $ps_{i}$ is calculated:

[TABLE]

where $\phi_{n}$ is the statistical estimator (e. g. mean or dispersion) defined for $n$ blocks and $\phi_{n-k}((X_{1},...,X_{n})_{[i]}$ is the same estimator but for the deleted-one sample. The pseudo-values, $ps_{i}$, follow a standard normal distribution for the $\phi$ parameter with mean and standard deviation.

We adopted $k=10$ for both samples, with $n_{b}=20$ blocks in the raw and $n_{b}=10$ blocks in refined sample. Following this method, the statistical errors indicated in Table 3 bound 1$\sigma$ or 68% of the likelihood distribution of each parameter. One potential concern is that imperfect modelling of the thin disk could affect fitted parameters of the thick disk and halo. In order to assess this possible degeneracy, we run multiple fits of the halo and thick disk with the thin disk density set to 60%–110% (with bin size equal to 10%) of the benchmark value listed in Table 1 (55.41 M*⊙* pc*-2*). Assuming an uncertainty of 10% in the local surface density of the thin disk (similar to the uncertainty of Holmberg & Flynn 2004), those trials indicate a strong dependence between the densities of thin/thick disk. A decrease of 10% in the modelled density of the thin disk means an increase of the same amount in the density of the thick disk, while for the remaining parameters the dependence is much weaker. In this way, we assume an uncertainty of 10% in the local density of the thin disk, and we added the systematic dependence on the thin disk local density as an systematic error in Table 3 for all the parameters. We assume that the correlation between the parameters in Table 3 and the parameters of the thin disk (with the exception of the density of the thick disk) is much weaker than the correlation between the same parameters and the normalization of the thin disk.

Following this reasoning, the systematic errors included in Table 3 account for the ability of the pipeline to recover input values, and the dependence of the parameters on the local density of the thin disk. The best-fit parameters for the raw and refined samples agree within 1$\sigma$ and have quite similar errors.

4.3 MWFitting results

There is a general agreement between our results and previous works (see Table 3), even though our uncertainties are smaller in most of the cases.

The vertical and radial scale of the thick disk are consistent within $\sim 1\sigma$ given the estimate and uncertainty from Jurić et al. (2008), and the density normalization of the thick disk is within 1$\sigma$ of the estimate and uncertainty from de Jong et al. (2010).

The large differences in the density of the thick disk reported by previous works may be related to the different methods used to estimate the total stellar mass. Different IMFs heavily influence the number of low-mass stars, most of which are not sampled by the HDs in this work. Different approaches in selecting stars also impact the estimation of the total stellar mass. Likewise, we point out there is a discrepancy by a factor of $\sim$2 in the local halo stellar density between the estimations of Jurić et al. (2008) and de Jong et al. (2010).

Comparing our measurements of the Galactic halo to the literature, the best-fit values of oblateness ($q$) are between the results of Jurić et al. (2008) and Deason et al. (2011) and that of de Jong et al. (2010). Regarding the inner and outer exponents of the double power law describing the halo density, we find that estimates from Deason et al. (2011) are steeper than ours, but that the two results are consistent to within 20%. This relative discrepancy could be due to many factors: minor tweaks in the stellar evolutionary models, the different regions sampled (SDSS imaged most the northern hemisphere while DES samples the southern hemisphere), or the other model parameters adopted. Similar explanations could account for the difference between our results and the single power-law fit by Hernitschek, et al. 2018 (n=$4.40^{+0.05}_{-0.04}$), in addition to the fact that they use RR Lyrae stars from Pan-STARRS1 between 20 kpc $\leq R_{GC}\leq$ 131 kpc, which extend to much larger distances than our sample ($r_{GC}\leq 60$ kpc).

Our model indicates a closer power-law break radius than that indicated by Deason et al. (2011); however, our best-fit break radius is consistent with the larger range of fits in the literature. In order to illustrate the range of distances for the radius of the break in previous works, we cite a few examples using diverse methods. For example, Watkins et al. (2009) use a sample of RR Lyrae stars in Stripe-82 region sampled by SDSS, finding a break radius of 23 kpc. Pila-Díez, et al. (2015) fit F stars from fields of MENeaCS and CCCP projects determining a power-law break at 20 kpc from the Galactic center, and in a more recent work Xue et al. (2015) modelled giant stars from SDSS/SEGUE-2 found a closer break radius than our value ($18\pm 1$ kpc).

In a more recent work, Deason, et al. (2018) determined the orbital properties of a sample of MS and BHB halo stars using position, kinematic properties and metalicites from Gaia DR2 and SDSS. Adding the Galactic gravitational potential, they derive the apocenter of the star’s orbits, addressing the break of the halo to a “pile-up” effect where the stars with eccentricity $e>0.9$ slow-down near the most distant part of the orbit. After excluding stars from the disk, the average apocenter derived for MS and BHB stars are 16$\pm$6 kpc and 20$\pm$7 kpc, respectively, in excellent agreement with our fit (see also Bullock & Johnston 2005; Deason, et al. 2013 and references therein).

We also fit two alternative models for the Galactic halo: an Einasto profile and a single power law (in both cases the thick disk was modeled with the same exponential profile). The best-fit model with halo modeled by an Einasto profile yielded a lower likelihood ($-2\ln\lambda=440{,}340$) than the model with the double power law ($-2\ln\lambda=218{,}355$), both following Eq. 1. The best-fitting model for thick disk and halo with a single power law resulted in an even lower likelihood ($-2\ln\lambda\simeq 1{,}000{,}000$). These conclusions are quite similar to those of Deason et al. (2011).

N-body simulations (e.g., Bullock & Johnston, 2005) show that an excess of stars in the central region of the halo (similar to the observed double power law of the halo) can be related to accretion events. Thus, the features of the halo’s stellar profile provide clues about the epoch, number an characteristics of past accretion events. Following this reasoning, the observed features of the best-fit halo model strongly favour a massive accretion event where the stars from the accreted satellite dominates the Galactic halo out to the break radius. We posit that these stars may be associated with the merger of the Gaia-Enceladus-Sausage galaxy (Belokurov, et al., 2018; Helmi et al., 2018). The stellar mass derived from integrating the best-fit model for the $raw\ sample$ from the Galactic center out to the break radius results in a mass of $\simeq 3.6\times 10^{8}\,M_{\odot}$. This estimate can be used as an upper limit for the current stellar mass of the Gaia-Enceladus-Sausage galaxy, excluding possible globular clusters associated to the former satellite. In a comparison to a recent work, Mackereth & Bovy (2020) selected a mono-abundance population of halo ($-3<$[Fe/H]$<-1$ and $0.0<$[Mg/Fe]$<0.4$) red giant stars from APOGEE-DR14, and they found a estimation of the current mass of stars with high eccentricity ($e>0.7$) associated to Gaia-Enceladus-Sausage galaxy ($3\pm 1\times 10^{8}\,M_{\odot}$) very close to our estimate above.

5 Simulating the stellar contents of DES-Y3

With the best-fit parameters, we produce a simulated stellar catalogue matched to DES-Y3 with limiting magnitude of $g=24$ and in the colour range $0<g-r<0.6$. We compare these simulations to the real data to study the stellar distribution in DES-Y3, to highlight asymmetries in the Galactic components (such as flares and warps in the disk), and to reveal stellar substructures.

Fig. 4 compares the star counts as a function of g magnitude in DES-Y3 to simulations using the best-fit models. The regions where the stars are sampled exclude areas containing dwarf galaxies, globular clusters, stellar streams, the Sagittarius Stream and Eridanus-Phoenix over-density, and regions with high reddening ($b\leq-30^{\circ}$). The magnitude bins in this figure are twice the size of the magnitude bins in the fitting, in order to sample a smooth histogram.

The distribution of DES-Y3 stars in Fig. 4 is shown as a blue line, while the distribution of stars in the simulations using the best-fit parameters from the raw and refined samples are shown as thick and thin green lines, respectively. Grey lines sample the distribution of modelled stars belonging to the bulge and disk (dotted line) or to the halo (dashed line), both following the best-fitting model for raw sample.

An initial look at Fig. 4 reveals a high level of consistency between the two best-fit models. The differences between both models are $<\ 2$% in general. These models are reasonably similar to the data, agreeing within 5% in the magnitude range $17<g<23$. The histogram shows an excess in the DES-Y3 data with respect to both best-fit models between $20<g<22$, with an excess in the modelled stars of a few percent between $18<g<19.5$. The discrepancies between data and model in Fig. 4 may be improved in several ways: better models for the evolution of metal-poor stars (population of the halo), minor tweaks in the halo’s SFH, additional components in the Galactic halo model (e.g., the Gaia-Enceladus galaxy), a potential metalicity gradient in the halo, slight changes to the Sun-Galactic center distance, or in any other parameter taken into account in the Trilegalmodels, such as the interstellar extinction. We are investigating the possible causes for the observed excess of stars, in order to develop an improved model for the Galaxy. Interestingly, the difference between data and models is dependent on the region of the sky examined, with better agreement found including only Galactic fields at higher latitude.

Fig. 5 shows the same distribution of stars, but in the $g\times g-r$ CMD space, with bins in magnitude and color equal to 0.2 and 0.02, respectively. We note that this bin size is different from that used in the fitting, but is equal to that in Fig. 4 for the $g$ magnitude. In order to highlight subtle differences in color, we oversample the $g-r$ color range with bin sizes equal to 0.02.

The left panel of Fig. 5 shows the distribution of DES-Y3 stars, similar to the blue line in Fig. 4. The best-fit model for the raw and refined samples are shown in the central and right panels representing the green lines in Fig. 4. Analogous to the observed magnitude distribution at the faint end in Fig. 4, the fainter end of the first panel in Fig. 5 shows a decreasing number of sources below $g=23$, where the dashed line delimits the bound of the stars compared to models in the fitting. The panels share the same colorbar, indicated on the right of the figure.

The three panels of Fig. 5 exhibit strong similarities down to $g\lesssim 23$. The thick disk leaves its main imprint by the plume of MSTO metal-rich stars at $g<19$ and $g-r\simeq 0.4$. There is a smooth transition between the crowding of MSTO stars of the thick disk and the MSTO stars of the halo, which starts at $g\simeq 19$ but in a bluer region. This transition is seen in the Fig. 4 as a distribution of stars slightly more flat ($18<g<19$) than the preceding or subsequent range. The MSTO halo’s stars are concentrated in a large range of magnitudes centered at $g\simeq 21$, whose density smoothly decreases towards the fainter end.

An excess of DES-Y3 stars in the range $21<g<22$ is seen on the left panel of Fig. 5, similar to Fig. 4, but with the additional information that the excess stars are concentrated near the MSTO of halo stars. The most populated bin in the left panel of Fig. 5 ($g-r$=0.33, $g=21.25$) contains 25% more stars when compared to the same bin in the central panel.

Even though they are not included in this comparison, the estimation of star counts fainter than $g=23$ is important for future surveys such as the Rubin Observatory LSST (LSST Science Collaboration et al., 2009) and Euclid (Sartoris et al., 2016), where S/G classification will be important. For example, at $g=24$, Fig. 4 shows that the expected number of halo stars in the models is roughly double the number of stars in the data. Realistic simulations for future large and deep surveys must consider and account for this incompleteness.

5.1 Poissonian significance maps

Fig. 6 shows the Poissonan significance maps generated for both samples of the DES-Y3 data using the best-fit model parameters. Given the steep decrease of stars at faint g-magnitudes in Fig. 4, we restrict the sample to stars with $g<23.5$. The significance of each $7\times 7$ arcmin2 pixel is taken as the residual star counts (difference between the DES-Y3 catalogue and the model catalogue) divided by the square root of modelled star counts. Both maps are smoothed by a Gaussian kernel with $\sigma$ = 7 arcmin, resulting in a minimum significance of -1.67 for refined sample and -1.69 for raw sample. In order to highlight under/overdensities as blue/reddish colours, and white colour as a perfect agreement between models and data, the significance range is set to [-1.67,1.67]. Pixels with significance higher than 1.67 (mainly known globular clusters and dwarf galaxies) are saturated with that maximum value.

Many known Galactic substructures are enhanced in this residual map, attesting to the accuracy of the mwfitting model. We label the most significant stellar over-densities on both panels of Fig. 6. For instance, the stripe roughly parallel to $l=180^{\circ}$ is the Sagittarius Stream. The over-density associated with SMC (SMCNOD) in the anti-LMC side (Pieres et al., 2017; Mackey et al., 2018) is also evident. Although we are not using a $matched\ filter$, a technique commonly applied to highlight fainter substructures as streams (e.g., Odenkirchen et al., 2003), a few streams are noticeable in Fig. 6. The ATLAS stream (Koposov et al., 2014; Shipp et al., 2018), a track of stars close to Galactic Pole (indicated in Fig. 6), is a good example of such a structure, as well as the Phoenix stream (Balbinot, et al., 2016), a long track of stars seemingly pointing toward the Phoenix dwarf galaxy. Other visible features are the Indus stellar stream (just below Tuc II), and the Tuc III stream, centered on the dwarf galaxy Tuc III.

The regions at the lowest Galactic latitudes between $240^{\circ}<l<270^{\circ}$ presents smooth and flat over-densities (with the exception of the region close to LMC) in both panels of Fig. 6, which may indicate that there is room for improvement in the thin disk model. The region at $b<-30^{\circ}$, $220^{\circ}<l<240^{\circ}$ in DES-Y3 footprint exhibits a strong excess of stars close to NGC1904, which may be the result of disk flaring or the Southern extension of the Monoceros Ring (Newberg et al., 2002).

The Eridanus-Phoenix over-density (Li et al., 2016) is a very large over-density of stars between $270^{\circ}<l<330^{\circ}$ and $-40^{\circ}<b<-70^{\circ}$, populating a triangle with vertices close of LMC, SMC and Fornax dwarf galaxy, seen on both panels of Fig. 6. Subtracting the stars in the modelled catalog, the Eridanus-Phoenix cloud contains an over-density of 4756 (4755) stars within the range ($17<g<22$ and $0.0<g-r<0.6$) when compared to the best-fit of the raw (refined) sample. Accounting for stars more massive than 0.1 M*⊙* in a Chabrier mass function (Chabrier et al., 2000) for a disk-like IMF stars, those values correspond to an object with $\simeq 1.6\times 10^{4}$ ($\simeq 1.5\times 10^{4}$) M*⊙* for the raw (refined) sample. These mass estimations represent a decrease in mass of at least by factor of five compared to the estimates in Li et al. (2016).

Even though the best-fit parameters for both samples agree within 1$\sigma$, there are slight differences regarding the two panels of Fig 6. For instance, ATLAS and Phoenix streams seem to be more continuous in the left panel with best-fitting parameters from raw sample than with the refined sample.

5.2 Milky Way stellar mass

We calculate the stellar masses of the halo and thick disk MW components with the best-fit parameters (Table 3) and list them in Table 4. These mass estimations only include field stars following from a smooth model for the Galactic components, and therefore exclude the mass from globular clusters, dwarf galaxies, and streams.

The bulge parameters are kept fixed, and the model described in Table 1 amounts to a stellar mass of $1.28\times 10^{10}$ M*⊙* or 21.4% of the total stellar mass of the Galaxy ($5.97\pm 0.99\times 10^{10}$ M*⊙, following the adopted model here). This agrees with mass estimates from the literature, where estimates of the stellar bulge mass range from 10-20% of the MW stellar mass (Licquia & Newman, 2015; Portail et al., 2017). Our model includes a thin disk (with fixed parameters) and has a stellar mass of $4.62\times 10^{10}$ M⊙, which is within 1$\sigma$ of the estimation by Licquia & Newman 2015 ($5.17\pm 1.11\times 10^{10}$ M⊙) and within 2$\sigma$ of McMillan 2011 ($5.54\pm 0.63\times 10^{10}$ M⊙*). The thick disk has a small contribution to total disk mass, with a ratio of stellar masses in the thin and thick disks as $\cong 13\,000:1$.

The halo mass is estimated by integrating the double power-law profile from the Galactic Centre out to 100 kpc. Based on our fits, we find that the stellar halo contributes 1.1% of the Galactic stellar mass, while the disks contribute with $\simeq$80% of the total. Our estimate of the total stellar halo mass is within the range estimated by Deason et al. 2011 ($2-10\times 10^{8}$ M*⊙), while being more massive than estimated by Bell, et al. 2008, where the latter authors found a halo with an integrated stellar mass out to 40 kpc of $3.7\pm 1.2\times 10^{8}$ M⊙*.

6 Concluding Remarks

We have developed a new code to fit the stellar components of the MW. In this first paper, we concentrate on fitting the thick disk and the halo due to the location of the DES footprint in the south Galactic cap. We list below our main conclusions from this work.

• •

This work presents mwfitting, a pipeline constructed to fit structural parameters for the Galactic components with Trilegal stellar population synthesis models.

• •

The mwfitting  pipeline is validated with synthetic catalogues. We successfully recovered the input parameters (with a maximum deviation $\leq$ 3%) using the same oversampling factor and a footprint smaller than the real data (see Table 2).

• •

Our main goal in this work is to model the halo and the thick disk components by applying the mwfitting pipeline to data from DES-Y3 Gold catalogue. We defined two different samples based on known stellar over-densities. Both samples excluded cells populated by dwarf galaxies, globular clusters and cells close to the LMC. In the refined sample, we further excluded cells where stellar streams and Eridanus-Phoenix over-density are located.

• •

Table 3 lists the results for both samples, with statistical uncertainties determined by jackknife resampling and the emcee method. The systematic uncertainties are sampled by the ability of the pipeline in recovering the true parameters based on simulations and the uncertainties in the local density of the thin disk. Results from both samples agree within a confidence level of 68% (1$\sigma$).

• •

The distribution of DES-Y3 stars presents a reasonable agreement (within $\leq 5\%$ in number of stars in each bin) with our models down to $g=23$. The distribution of stars in the DES-Y3 catalogue and in the models both peak close to $g=22.0$. Fainter than $g=23$, there is a decrease in the number of stars, that we interpret as a result of the S/G classification schema applied here, coupled with the relative scarcity of stars in the outer MW halo.

• •

CMDs comparing DES-Y3 stars and both simulations reasonably agree down to $g=23$, suggesting that the double power law is a good description of the Galactic halo, at least at this depth.

• •

The star counts in the stellar halo is crucial for predicting the density of faint stars with $g-r\lesssim 1$, which will be sampled in future surveys such as the Rubin Observatory LSST and Euclid.

• •

Simulations over the entire DES-Y3 footprint based on our best-fit models were produced. Both simulations agree well with the data. Residual maps highlight many over-densities associated with globular clusters, dwarf galaxies, clouds, and streams in the DES footprint.

• •

We found a mass ratio between the thin and thick disks equal to $\simeq$ 13000:1, while the halo amounts to 1.1% of the total MW stellar mass.

Future work with mwfitting will include data from other wide-field surveys to extend the analysis to both the north and south Galactic hemispheres and will include improvements to the modeling for the Galactic halo (e.g., tri-axial models).

Acknowledgements

The authors are grateful to James Binney for many useful suggestions and comments. We thank Cecilia Mateu for the discussion about the location of the streams on the eastern part of the DES footprint.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020, and the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) do e-Universe (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Affiliations

1* Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil *

2* Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil *

3* Osservatorio Astronomico di Padova, INAF, Vicolo dell’Osservatorio 5, I-35122 Padova, Italy *

4* Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands *

5* Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS - 91501-970, Brazil *

6* Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain *

7* George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA *

8* LSST 933 North Cherry Avenue, Tucson, Arizona 85721, USA *

9* Koninklijke Sterrenwacht van België, Ringlaan 3, B–1180 Brussels, Belgium *

10* Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA *

11* Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA *

12* Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 16, CH-8093 Zurich, Switzerland *

13* Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049 Madrid, Spain *

14* CNRS, UMR 7095, Institut d’Astrophysique de Paris, F-75014, Paris, France *

15* Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d’Astrophysique de Paris, F-75014, Paris, France *

16* Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK *

17* Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, CA 94305, USA *

18* SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA *

19* Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA *

20* National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA *

21* Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain *

22* Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India *

23* Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA *

24* Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA *

25* Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain *

26* Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain *

27* Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA *

28* Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA *

29* Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA *

30* Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA *

31* Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA *

32* Department of Physics, The Ohio State University, Columbus, OH 43210, USA *

33* Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA *

34* Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia *

35* Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain *

36* Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA *

37* Brookhaven National Laboratory, Bldg 510, Upton, NY 11973, USA *

38* School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK *

39* Brandeis University, Physics Department, 415 South Street, Waltham MA 02453 *

40* Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859, Campinas, SP, Brazil *

41* Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831 *

42* Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth, PO1 3FX, UK *

43* Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA *

44* Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile *

Appendix A MWFitting pipeline inputs

This appendix describes input parameters of the pipeline when submitting MWFitting through the science portal, intended to guide LIneA users. Table 5 lists the name, description, standard configuration and units of the parameters that the user should use to reproduce our results. These parameters are not related to the components of the MW models, but needed to run the code.