Chemical Mapping of the Milky Way with The Canada–France Imaging Survey: A Non-parametric Metallicity–Distance Decomposition of the Galaxy

Over the course of the coming decade, our view of the cosmos will be revolutionized by a series of unprecedented new surveys. Perhaps the most exciting of these in the immediate future is the Gaia satellite, a cornerstone of the European Space Agency's science strategy, which will survey the astrometric sky, taking measurements of the minute motions of about a billion stars in our Milky Way and the Local Group to understand the detailed formation history of our Galaxy. Although Gaia's precision measurements will have profound implications for many areas of astrophysics, their most obvious application will be for studies of the formation, and subsequent dynamical, chemical, and star formation evolution of the Milky Way. Indeed, in this new era, many of the questions of the formation history of galaxies, which we normally associate with high-redshift studies, will be addressed with unprecedented spatial detail by looking directly at the remnants of the structures that made our Galaxy.

The Gaia mission will provide the kinematic dimensions (particularly proper motions) that are largely absent from existing surveys and will bring about a phenomenal increase in the data quality and quantity for the nearby Galaxy. The position of every object in the sky brighter than $G\sim 20.5$ mag (over 1 billion objects) will be mapped with a positional accuracy reaching micro-arcseconds for the brightest stars. A spectrometer will provide radial velocity information and abundances for stars brighter than $G\sim 16$ mag. In addition, a photometer will measure the spectral energy distribution with sufficient resolution to estimate stellar metallicities at G = 15 mag to ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]=0.1$–0.2 dex (for FGKM stars; Liu et al. 2012). For the brightest stars ($G\lt 12$), atmospheric information and interstellar extinction will also be derived. Thus Gaia will undoubtedly provide the foundation for much of the next generation of research in Galactic and stellar astronomy, themselves the foundation for much of the rest of astrophysics.

One of the most exciting problems that Gaia will be able to contribute to is the unveiling of the dark matter distribution in the Milky Way, both on global ($\sim 100\,\mathrm{kpc}$) and small scales ($\lt 1\,\mathrm{kpc}$). The key observables that Gaia will bring to this analysis are excellent proper motion measurements—in $\mathrm{mas}\,{\mathrm{yr}}^{-1}$—of individual stars in situ in the halo. This angular displacement information must be coupled with reasonable distance measurements, to have access to the physical transverse velocities (in, e.g., $\,\mathrm{km}\,{{\rm{s}}}^{-1}$), and of course to know where the stars are in three-dimensional space. Although Gaia will also measure stellar parallaxes, such measurements will not be available for the vast majority of the surveyed halo stars, which have faint magnitudes (see Figure 1).

Zoom In Zoom Out Reset image size

In situ halo stars (say with distance $D\gt 10\,\mathrm{kpc}$) with $G\gt 18$ will not have useful Gaia parallax measurements. A further problem is that these faint stars are predominantly A-, F-, and G-type main-sequence (MS) dwarfs. The Gaia spectrophotometer will not give useful astrophysical parameters for such stars: as explained in detail in Bailer-Jones et al. (2013), at $G=19$ the metallicity uncertainty is expected to be ${\rm{\Delta }}[\mathrm{Fe}/{\rm{H}}]=0.6\mbox{--}0.74$, while the surface gravity uncertainty is expected to be ${\rm{\Delta }}\mathrm{log}g=0.37\mbox{--}0.51$. Adopting the Ivezić et al. (2008) metallicity-dependent photometric parallax calibration, even with perfect g- and r-band photometry, such metallicity uncertainties would typically incur $\gt 25$% distance errors. It is hence imperative to obtain alternative distance measurements to enable halo science with Gaia. This is especially critical given the low density of bright halo tracers (∼7 per deg² up to $G=18$ in the Besançon model simulation shown in Figure 1).

Thus, to enable much of the next generation of exciting halo science, we need to be able to measure the distances for most stars in the Gaia catalog, and, if possible, to measure distances for even fainter stars, since additional information can of course be extracted from the faint stellar populations associated with the stars detected by Gaia.

Fortunately, MS stars, which are the most numerous halo sources in the Gaia catalog, have a relatively well defined color–luminosity relation that can be exploited to derive their distances based only on multiband photometry. This is a consequence of the fact that the MS locus is not extremely sensitive to metallicity (see Equation (3) below), and the effect of age is to depopulate the bluer stars while maintaining the shape of the redder MS. Using Sloan Digital Sky Survey (SDSS) data, Jurić et al. (2008) exploited this property to derive distances for 48 million stars out to $\sim 20\,\mathrm{kpc}$ using effectively just r-band magnitudes and $(r-i)$ colors. In a subsequent landmark study, Ivezić et al. (2008, hereafter I08) demonstrated a tight correlation between the spectroscopically measured metallicity of MS stars (Lee et al. 2008) and their $(u-g)$, $(g-r)$ colors (see Figure 2 for the CFIS-u version of this color–color diagram). Systematic calibration errors are small compared to typical random errors from photometric uncertainties, allowing measurements of metallicity from ugr photometry that, for sufficiently large samples, are comparably precise to spectroscopic measurements, but much cheaper. However, to keep the random error from exceeding 0.3 dex (after which it becomes difficult to cleanly discriminate different Galactic populations), a maximum uncertainty of ∼0.03 mag in u is required. In the SDSS, this occurs at $u\sim 19.3$, which translates to an effective distance threshold for turn-off stars of $\sim 5\mbox{--}10\,\mathrm{kpc}$. This is despite the fact that the SDSS g-band photometry is substantially deeper, reaching $g\sim 20.7$ with similar uncertainty. Thus for the purpose of measuring photometric metallicities of MS stars, the SDSS u-band is really ∼2.7 mag too shallow for its g-band. Indeed, the u-band depth was the limiting factor in the I08 analysis, affecting their sample size and the discriminating power of the data.

Zoom In Zoom Out Reset image size

I08 also showed that the photometric metallicity measurement allowed in turn for the photometric distance to be refined. They set their absolute magnitude calibration M_r to depend on color $(g-i)\equiv x$ and metallicity $[\mathrm{Fe}/{\rm{H}}]$, so that

$$\begin{eqnarray}&&{M}_{r}(x,[\mathrm{Fe}/{\rm{H}}])={M}_{r}^{0}(x)+{\rm{\Delta }}{M}_{r}([\mathrm{Fe}/{\rm{H}}]),\end{eqnarray} \tag{ 1 }$$

where the color term was fitted to be

$$\begin{eqnarray}{M}_{r}^{0}(x) & = & -5.06+14.32x+14.32x-12.97{x}^{2}\\ & & +6.127{x}^{3}-1.267{x}^{4}+0.0967{x}^{5},\end{eqnarray} \tag{ 2 }$$

and the metallicity correction term is

$$\begin{eqnarray}&&{\rm{\Delta }}{M}_{r}({[\mathrm{Fe}/{\rm{H}}])=4.50-1.11[\mathrm{Fe}/{\rm{H}}]-0.18[\mathrm{Fe}/{\rm{H}}]}^{2}.\end{eqnarray} \tag{ 3 }$$

From Equation (3) it can be appreciated that the absolute magnitude is highly sensitive to metallicity. This is especially important when dealing with Galactic halo populations: if we were to assume a disk-like metallicity of $[\mathrm{Fe}/{\rm{H}}]=-0.5$ for a halo star with true metallicity $[\mathrm{Fe}/{\rm{H}}]=-1.5$, we would incur a 0.75 mag error in M_r (i.e., a 41% distance error), rendering any tomographic analysis invalid.

Here we use u-band data from the new Canada–France Imaging Survey (CFIS, which we present in detail in Ibata et al. 2017; hereafter Paper I) to greatly improve on the SDSS u-band photometry and thereby probe the MS populations of the Milky Way out to much larger distances than was possible with the SDSS. CFIS is a large community program at the Canada–France–Hawaii Telescope that was organized to obtain u- and r-band photometry needed to measure photometric redshifts for the Euclid mission (Laureijs et al. 2011; Racca et al. 2016). As we show in Paper I, at a limiting uncertainty of 0.03 mag for point sources, CFIS-u reaches approximately three magnitudes deeper than the SDSS u-band data. Although the final CFIS-u survey area will be 10,000$\,{\deg }^{2}$, covering most of the northern hemisphere at $| b| \gt 25^\circ$, here we present our first metallicity analysis, based on ∼2900$\,{\deg }^{2}$, mostly contained in the declination range of $18^\circ \lt \delta \lt 45^\circ$ (see Figure 1 of Paper I).

The outline of this paper is as follows. In Section 2, we describe the methods that we use to measure photometric metallicity, while the effect of contamination from giants and subgiants is examined in Section 3, and the survey completeness in Section 4. In Section 5, we explore the distributions in metallicity and distance throughout the Milky Way, and present a new algorithm to deconvolve the survey into sub-populations in Section 6. We conclude with a discussion of our findings in Section 7.

In this section, we will lay out the procedure by which photometric metallicities are measured. In our first attempts, this was achieved by combining CFIS-u with SDSS $g,r,$ and i, but we found that the accuracy of the SDSS measurements in these bands limited the depth we could attain. Since the release of the Pan-STARRS1 catalog (PS1, Chambers et al. 2016) in 2016 December, we now have access to more accurate $g,r,$ and i measurements, which significantly improve the depth and the size of the sample of stars for which we are able to measure good metallicities. Nevertheless, after extensive experimentation, we realized that by combining the SDSS and PS1 measures, we could obtain still deeper and more reliable photometry. We found, in particular, that metallicity outliers were often stars that had inconsistent SDSS and PS1 values (this will be detailed further below). To combine the SDSS and PS1 magnitudes, we first shifted the PS1 magnitudes onto the SDSS system using the simple linear transformations derived by Tonry et al. (2012), and calculated their uncertainty-weighted mean fluxes f and flux uncertainties $\delta f$ as

$$\begin{eqnarray}f & = & ({f}_{\mathrm{SDSS}}/\delta {f}_{\mathrm{SDSS}}^{2}+{f}_{\mathrm{PS}1}/\delta {f}_{\mathrm{PS}1}^{2})/S,\\ \delta f & = & \sqrt{1/S},\end{eqnarray} \tag{ 4 }$$

where $S\equiv \delta {f}_{\mathrm{SDSS}}^{-2}+\delta {f}_{\mathrm{PS}1}^{-2}$. Unless stated otherwise, the $g,r,$ and i magnitudes we refer to below are these combined PS1+SDSS values, on the SDSS system, while the $u$ magnitudes are from CFIS, calibrated as explained in Paper I.

CFIS-u includes a large number of stars with high-quality spectroscopy obtained as part of the SDSS-Segue project (Yanny et al. 2009). Selecting those objects in the SDSS-DR10 spectroscopic catalog that also have a less than 10% chance of being photometric variables in Pan-STARRS1 (according to the variability analysis of Hernitschek et al. 2016), and cross-matching against the current CFIS-u catalog gives a total of 65403 fiducial stars that we can use to define the photometric metallicity procedure. The motivation for removing the PS1 photometric variables is that this should help in reducing the noise in the metallicity relation, since the photometric measurements are not coeval.

We use the SDSS spectroscopic data set to make the $u-g$, $g-r$ color–color plot shown in Figure 2(a), where each circle is color-coded by spectroscopic metallicity. We use this sample as a training set to construct the photometric metallicity relation. The training sample is selected to retain only those objects classified as dwarfs by the Segue pipeline, with $3\lt \mathrm{log}g\lt 5$. Following I08, we further limit the stars to ${(g-r)}_{0}\gt 0.2$ to avoid contamination by blue horizontal branch stars, and also to ${(g-r)}_{0}\lt 0.6$ to keep preferentially MS dwarfs over giants (see Figure 1 in I08). (Extinction corrections were derived using the Schlegel et al. 1998 dust maps.) As we show in Figure 2, by applying a further cut to keep ${(u-g)}_{0}\lt 1.35$ (also used by I08), the metallicity bias of the ${(g-r)}_{0}\,=\,0.6$ limit is minimized. Finally, we also select stars to lie within the purple line boundary in Figure 2, to avoid extrapolation away from the color–color parameter region inhabited by the Segue stars. We fit the training sample with a two-dimensional Legendre polynomial, using only those stars with good CFIS-u measurements ($\delta u\lt 0.03$). The model is allowed terms in up to x³ and y³, which with cross-terms gives a total of 10 parameters. It can be seen that there is a close correspondence of the color representation of the SDSS spectroscopic [Fe/H] values (Figure 2(a)) and the corresponding model value interpolated from the photometry (Figure 2(b)). Using this model, we can now effectively place a CFIS-u survey star onto this color–color plane and read off the photometric metallicity, as would be measured by Segue. We note that the spectroscopic metallicity measure that we use is the "adopted" or FEH_ADOP value from the Segue Stellar Parameter Pipeline (SSPP; Lee et al. 2008). The other metallicity measures provided by the SSPP performed marginally less well in the tests below.

For comparison, in Figure 3, we display the theoretically expected behavior of a young (panel a) and an old (panel b) stellar population in the $(u-g)$ versus $(g-r)$ plane as a function of metallicity according to the PARSEC models (Bressan et al. 2012). This resembles closely the observed distribution (Figure 2), and one can appreciate that the derived metallicity should be independent of age. However, there are areas of this plane that are not occupied by old stars and it may be possible in future work to use this fact to identify distant younger stars in the halo.

Zoom In Zoom Out Reset image size

We notice that a small amount of contamination is caused by stars whose positions in the ${(g-i)}_{0}$ versus ${(g-r)}_{0}$ color–color plane is unusual. These objects are possibly stars whose photometric measurements are poor, variables that were not filtered out with the PS1 variability criterion, or possibly just unusual stars. We removed these objects by constructing the color index

$$\begin{eqnarray}&&{(g-i)}_{\mathrm{diff}}\equiv {(g-i)}_{0}-1.36{(g-r)}_{0}+0.023,\end{eqnarray} \tag{ 5 }$$

which is narrowly distributed around zero (with a standard deviation of 0.024 mag between $g=17$ and 18), and selected only those stars with $| {(g-i)}_{\mathrm{diff}}| \leqslant 0.1+3\delta {(g-i)}_{\mathrm{diff}}$ (after the quality cuts detailed below are applied, this selection removes only 0.3% of the sources in the final catalog).

In Figure 4(a), we show the result of applying the two-dimensional interpolation function of Figure 2 (which we will call "Method 1") to the CFIS-u stars with SDSS/Segue counterparts, this time using no spectroscopic information (i.e., we do not use the Segue surface gravity estimate) to cull the sample (see Section 3). The correspondence between the predicted photometric $[\mathrm{Fe}/{\rm{H}}]$ (plotted on the ordinate) and the spectroscopically measured value is good over most of the range for these stars, although clearly the method performs less well below $[\mathrm{Fe}/{\rm{H}}]\sim -2$. The color coding of the points represents their ${(u-{u}_{\mathrm{SDSS}})}_{0}$ colors; it can be appreciated that most of the strong outliers have high values of $| {(u-{u}_{\mathrm{SDSS}})}_{0}|$. Note, however, that we do not cull on ${(u-{u}_{\mathrm{SDSS}})}_{0}$ at this stage, since most CFIS-u stars do not have good SDSS u-band measurements. In Figure 4(b), we show the corresponding distribution of metallicity differences; fitting this distribution with a Gaussian (using a $2.5\sigma$-clipping algorithm) gives the function shown in pink, which has $\sigma =0.23$ dex. Note that the average uncertainty of the spectroscopic metallicity measurements for this sample is $\delta [\mathrm{Fe}/{\rm{H}}]=0.04$, so almost all of the scatter is due to the intrinsic color variation of MS stars of identical metallicity, plus the photometric uncertainties of the SDSS, PS1, and CFIS-u surveys.

Zoom In Zoom Out Reset image size

During our data exploration, we constructed versions of Figure 4 that were color-coded with the quantity ${(g-i)}_{\mathrm{diff}}$ instead of ${(u-{u}_{\mathrm{SDSS}})}_{0}$. These showed clearly that $[\mathrm{Fe}/{\rm{H}}]$ also depends on the ${(g-i)}_{\mathrm{diff}}$ parameter, meaning that stellar metallicity (unsurprisingly) is not only a function of ${(u-g)}_{0}$ and ${(g-r)}_{0}$, but also of ${(g-i)}_{0}$. We therefore refitted the training sample with a three-dimensional set of colors for each star, namely ${(u-g)}_{0}$, ${(g-r)}_{0}$, and ${(g-i)}_{0}$. The fitting function we used was a Legendre polynomial with up to cubic terms in each variable, and in the cross-terms (20 parameters). The result is shown in Figure 4(c). The residuals displayed in Figure 4(d) are now significantly better ($\sigma =0.20$ dex), meaning that this procedure (Method 2) should be preferred over Method 1, especially since good i-band measurements exist for almost all CFIS-u stars.

However, Figures 4(c) (and 4(a)) also shows another interesting property: the stars on the upper side of the sequence typically have higher values of ${(u-{u}_{\mathrm{SDSS}})}_{0}$ than those on the bottom side of the sequence (i.e., the upper envelope shows more green points, while blue points are more common on the lower side). Thus there is also residual metallicity information in the ${(u-{u}_{\mathrm{SDSS}})}_{0}$ quantity. This is also not surprising: the transmission curves of the SDSS u-band filter and the new CFHT u-band filter are very different (as we show in Figure 2 of Paper I), and their magnitude difference encodes information about the ∼3800–4000 Å interval, which notably contains the metallicity-sensitive Ca ii H+K lines (3968.5 and 3933.7 Å). To harness the metallicity information in this additional color, we implemented a four-dimensional metallicity fit using the four variables: ${(u-{u}_{\mathrm{SDSS}})}_{0}$, ${(u-g)}_{0}$, ${(g-r)}_{0}$, and ${(g-i)}_{0}$. The fit is allowed up to cubic terms in all variables, and includes all cross-terms up to third order (for a total of 34 parameters).

The result of this new fit (Method 3) is shown in Figure 4(e), and the corresponding distribution of ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{Photometric}}\,\mbox{--}{[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{Spectroscopic}}$ is displayed in Figure 4(f). Since we are now assuming we have a reasonable SDSS u-band measurement, we cull the sample to keep those stars within $-0.15\lt {(u-{u}_{\mathrm{SDSS}})}_{0}\lt 0.1$ (an approximately $3\sigma$ interval around the mean of ${(u-{u}_{\mathrm{SDSS}})}_{0}$). The residuals are better for this 4D fit ($\sigma =0.19$ dex) than for Method 1, and it can be seen that the color distribution in Figure 4(e) does not show an obvious bias with respect to metallicity, contrary to what was seen in Figures 4(a) and (c). This procedure can be used for the brighter stars with $\delta {u}_{\mathrm{SDSS}}\lesssim \ 0.05$ mag (in our final metallicity catalog, 59% of stars that pass all the quality criteria have have $\delta {u}_{\mathrm{SDSS}}\lt 0.05$ mag).

Some of the sky areas we observed with CFIS-u have not been covered by the SDSS, but they do contain PS1 photometry. We therefore recomputed the "Method 2" photometric metallicity calibration using PS1 as the source for the ${g}_{{\rm{P}}1},{r}_{{\rm{P}}1},{i}_{{\rm{P}}1}$ magnitudes. The result is shown in Figure 5, and the scatter in the photometric metallicity measurements turns out to be only marginally worse in this case ($\sigma =0.21$ dex). We noticed that the significantly different g-band transmission curves between the SDSS and PS1 (see, e.g., Tonry et al. 2012) contain some metallicity information, as can be seen from the color distribution in Figure 5(a). However, we found that fitting the ${g}_{{\rm{P}}1}\mbox{--}{g}_{\mathrm{SDSS}}$ information does not improve the scatter more than what was obtained through "Method 3" above, so we will ignore it henceforth.

Zoom In Zoom Out Reset image size

With these fits, it is worth re-examining the photometric accuracy required to derive a good metallicity measurement. Selecting stars in a narrow interval around ${(g-r)}_{0}\,=\,0.3$, we find an approximately linear relation in ${(u-g)}_{0}$ versus $[\mathrm{Fe}/{\rm{H}}]$ with slope $d{(u-g)}_{0}/d[\mathrm{Fe}/{\rm{H}}]=0.15$. This means that even with perfect $g,r$ photometry, a u-band uncertainty of $\delta u=0.03$ will cause a $\delta [\mathrm{Fe}/{\rm{H}}]=0.2$ dex random metallicity error. We set this value as the maximum u-band uncertainty that can be tolerated: i.e., where the random error becomes equal to the intrinsic scatter in the photometric metallicity determination procedure. At present, the number of CFIS-u stars with this photometric accuracy and that are also in the SDSS DR13 point-source catalog is $5.6\times {10}^{6}$. As we show in Figure 5 of Paper I, CFIS-u is approximately 3 mag deeper than the SDSS at this uncertainty limit, i.e., we can measure stars that are a factor of roughly four times further away at the necessary ${\rm{S}}/{\rm{N}}$ than was previously possible with the SDSS. For reference, we list the interpolation functions employed in Methods 1 and 2 in Table 1.

Table 1.  Photometric Metallicity Interpolation Functions

term	Method 1	Method 1	Method 2	Method 2
	Polynomial pi	Coefficient ai	Polynomial pi	Coefficient ai
1	1	−43.331	1	−15.0144
2	x	94.081	x	50.8150
3	y	−37.547	y	−596.8324
4	$\tfrac{1}{2}(3{x}^{2}-1)$	−9.450	z	511.4958
5	$\tfrac{1}{2}(3{y}^{2}-1)$	−39.941	$\tfrac{1}{2}(3{x}^{2}-1)$	−11.1208
6	${xy}$	−40.615	$\tfrac{1}{2}(3{y}^{2}-1)$	7.6707
7	$\tfrac{1}{2}(5{x}^{3}-3x)$	3.686	$\tfrac{1}{2}(3{z}^{2}-1)$	9.1469
8	$\tfrac{1}{2}(5{y}^{3}-3y)$	−41.192	${xy}$	−1.7758
9	$\tfrac{1}{2}(3{x}^{2}-1)y$	−31.657	${xz}$	−19.3549
10	$\tfrac{1}{2}(3{y}^{2}-1)x$	116.244	${yz}$	−73.5756
11	⋯	⋯	$\tfrac{1}{2}(5{x}^{3}-3x)$	3.8139
12	⋯	⋯	$\tfrac{1}{2}(5{y}^{3}-3y)$	−208.2240
13	⋯	⋯	$\tfrac{1}{2}(5{z}^{3}-3z)$	84.2878
14	⋯	⋯	$\tfrac{1}{2}(3{x}^{2}-1)y$	−19.8719
15	⋯	⋯	$\tfrac{1}{2}(3{x}^{2}-1)z$	−6.4548
16	⋯	⋯	$\tfrac{1}{2}(3{y}^{2}-1)x$	22.5696
17	⋯	⋯	$\tfrac{1}{2}(3{y}^{2}-1)z$	749.5947
18	⋯	⋯	$\tfrac{1}{2}(3{z}^{2}-1)x$	7.3065
19	⋯	⋯	$\tfrac{1}{2}(3{z}^{2}-1)y$	−582.5109
20	⋯	⋯	${xyz}$	68.5820

Note. We list here the multi-dimensional Legendre polynomials p_i and the fitted coefficients a_i that were used to interpolate metallicity from photometry, for both Methods 1 and 2. For clarity, the polynomials are defined in terms of $x\equiv {u}_{0}-{g}_{0}$, $y\equiv {g}_{0}-{r}_{0}$, and $z\equiv {g}_{0}-{i}_{0}$. The g-, r-, and i-band values should be on the SDSS system, while the u-band should be on the CFIS system. The interpolated result is ${\sum }_{i}{a}_{i}{p}_{i}$. The $(u-g,g-r)$ vertices of the polygon (purple line in Figure 2) within which these interpolation functions have been fitted are $(1.350,0.600)$, $(1.170,0.600)$, $(1.105,0.576)$, $(0.962,0.522)$, $(0.888,0.489)$, $(0.779,0.431)$, $(0.715,0.390)$, $(0.617,0.287)$, $(0.600,0.200)$, $(0.990,0.200)$, $(1.072,0.275)$, $(1.116,0.304)$, $(1.190,0.339)$, $(1.350,0.405)$.

Download table as:  ASCIITypeset image

We note here in passing that we refrained from using the z-band photometry from SDSS or PS1, because the MS stars of interest here are typically blue and faint, and thus likely to have large z-band uncertainties.

The metallicity error distributions shown previously in Figure 4 to prove the efficacy of the photometric metallicity method developed in the previous section included the giant stars in common between CFIS-u and SDSS. However, the metallicity calibration procedures used SDSS dwarf stars as a training sample, so it is useful to check how well the measurements perform on stars of other luminosity classes. This is reported in Figure 6, which displays the spectroscopic versus photometric $[\mathrm{Fe}/{\rm{H}}]$ measurements (using Method 2) for stars with $\mathrm{log}g\lt 3.5$. The color of the points encodes surface gravity, as measured from the spectra. It transpires that the photometric metallicity is biased (by $+0.16$ dex), and the errors are larger ($\sigma =0.25$ dex) than for dwarfs, but nevertheless, it is clear that our measurements will retain useful chemical information on the giants and subgiants that will contaminate our sample.

Zoom In Zoom Out Reset image size

Identifying the stars that are not dwarfs is clearly very challenging from photometry alone. Efforts reported in the literature include the use of narrowband filters that measure the strength of gravity-sensitive lines (see, e.g., Majewski et al. 2000). We will return to the issue of dwarf-giant discrimination in a later contribution in this series. However, for the present purposes, it is useful to estimate the extent of the contamination problem. To this end, we present in Figure 7 the distribution of luminosity classes as a function of photometric distance in a $100\,{\deg }^{2}$ simulation toward the Galactic pole using the Besançon model (Robin et al. 2003).

Zoom In Zoom Out Reset image size

The top row (panels a–c of Figure 7) shows the luminosity class distribution as a function of distance and metallicity, adopting photometric selection criteria that are essentially identical to those of I08 ($15\lt {u}_{\mathrm{CFHT}}\lt 22;$ $14\lt g\lt 22;$ $14\lt r\lt 22;$ $0.2\lt g-r\lt 0.6;$ $0.6\lt u-g\lt 1.35$). We will call this color and magnitude selection the "wide cut." To match the observations, the abscissae show photometric distances derived from Equation (1) (i.e., assuming that the stars belong to the MS), using the Besançon model photometry (the Besançon model was calculated for CFHT filters, and $g,r$ were converted into SDSS magnitudes using the transformations given in Regnault et al. 2009). While the fraction of dwarfs is $\gtrsim 80$% over most of the distance intervals, it is striking that for metal-poor stars (panel a) the sample will be highly contaminated by giants. The reason for this is that the model predicts that the density of genuinely nearby metal-poor MS stars is very low. Distant very metal-poor giants that pass the color cuts unfortunately end up heavily contaminating the counts at derived photometric distances $\lesssim 7\,\mathrm{kpc}$.

To overcome this problem, we also consider a stricter selection to $0.2\lt g-r\lt 0.35$ in panels (d–f) of Figure 7, which we will call the "narrow cut." This color selection removes the potential giant contamination in metal-poor stars that falsely appear to be nearby according to their photometric distances. However, the stricter color interval removes sensitivity to old metal-rich stars, as shown in Figure 8. We are therefore forced to work with both the "wide" and "narrow" cuts, but we will keep their respective limitations in mind when interpreting the findings.

Zoom In Zoom Out Reset image size

It is relevant to note here that the Besançon model uses a shallow density law for the stellar halo (with a power-law exponent of −2.44). From fitting the density of blue horizontal branch and blue straggler stars in the SDSS, Deason et al. (2011) find instead that the stellar halo is best reproduced by a double power-law component, with an outer power-law of exponent approximately −4.6 beyond a break at a (Galactocentric) distance of $\sim 27\,\mathrm{kpc}$. Fitting a single power law to counts of halo red giant branch stars in the SDSS, Xue et al. (2015) find an exponent of −4.2 ± 0.1. If these analyses are correct, then the Besançon model significantly overestimates the number of giants and subgiants in our sample, and the contamination is much lower than Figure 7 would suggest.

The distance of the MS stars of interest increases with magnitude, as described by Equation (1), so the more distant objects suffer from larger uncertainties, and may be less likely to be detected. The resulting incompleteness could give rise to a distance-dependent bias in any study of (for example) the density of stars in the survey. It is important therefore to quantify this effect.

Ideally, the completeness of the sample would be measured by comparison to a much deeper survey, but this is not possible at present. Instead, we examine the counts as a function of magnitude using the SDSS "Stars" catalog; these point sources are known to be complete to $\gt 95$% to $g=22.2$ (Adelman-McCarthy et al. 2006).

We undertake the completeness analysis in the large region $120^\circ \lt {\rm{R}}.{\rm{A}}.\,\lt \,220^\circ$, $25^\circ \lt \mathrm{decl}.\,\lt \,33^\circ$, which encompasses the north Galactic pole, and which is fully covered by the SDSS, PS1, and CFIS. Figure 9(a) shows the $g-r$, $g$ color–magnitude diagram (CMD) of this region, with dashed and dotted lines marking the limits of the "wide" and "narrow" cuts, respectively. The light-colored histograms in panels (b) and (c) show the magnitude distribution of the SDSS point sources for $0.2\leqslant {(g-r)}_{0}\lt 0.6$ and $0.2\leqslant {(g-r)}_{0}\lt 0.35$, respectively. The dark-colored histograms are the corresponding CFIS-u counts, matched to the SDSS point sources. To ${g}_{0}=22.2$ (and within the corresponding color intervals), CFIS-u detects $\gt 99$% of all SDSS point sources.

Zoom In Zoom Out Reset image size

Imposing the additional color cut $0.6\leqslant {(u-g)}_{0}\leqslant 1.35$ alters (b) and (c) into the distributions displayed in (e) and (f), respectively (the color cut is shown on the $u-g$, $g$ CMD in panels (d) and (g).

To measure metallicities to ∼0.2 dex requires reliable, accurate photometry, and by extensive experimentation, we have found that this can be achieved by the following.

• 1.  
  Requiring that the SDSS and PS1 $g$ and $r$ measures of individual stars agree to within 3 standard deviations (the PS1 values are first transformed onto the SDSS system using the linear equations of Tonry et al. 2012).
• 2.  
  Requiring that the ${(g-i)}_{\mathrm{diff}}$ quantity remains within the bounds $| {(g-i)}_{\mathrm{diff}}| \leqslant 0.1+3\delta {(g-i)}_{\mathrm{diff}}$.
• 3.  
  Requiring an upper limit of 0.03 mag uncertainty in all of $u,g,r$.

With these additional criteria, we obtain the magnitude distributions shown in dark histograms in panels (h) and (i) (for the "wide" and "narrow" cuts, respectively). The ratio between these distributions and the corresponding light-colored histograms define the completeness functions for these populations. We find that the completeness remains above 50% until $g=21.1$, and will use this value henceforth as the limiting cutoff of the metallicity sample. Given the $T=12\,\mathrm{Gyr}$ models in Figure 8, this corresponds to a distance limit of approximately $30\,\mathrm{kpc}$ for a metal-poor star with $[\mathrm{Fe}/{\rm{H}}]\,=-2.3$ and $\sim 22\,\mathrm{kpc}$ for a star of $[\mathrm{Fe}/{\rm{H}}]=-0.5$.

It should be possible to explore further in distance by relaxing the $u,g,r$ photometric accuracy, at the cost of a lower metallicity accuracy. We leave this issue to a subsequent contribution, however.

We will now employ CFIS to examine how the stellar populations vary as a function of height above the Galactic plane. To this end, we consider the $1.9\times {10}^{5}$ stars that lie toward the North Galactic Cap (NGC, which we define to be the sky above $b\gt 70^\circ$) and that pass the quality criteria listed in Section 4, and we also consider an intermediate latitude sample with $1.0\times {10}^{6}$ stars in the range of $30^\circ \lt | b| \lt 70^\circ$.

In Figure 10(a), we show how the metallicity varies in the NGC sample as a function of vertical distance z above the plane from 0.63 to $31.6\,\mathrm{kpc}$ (although for clarity we have transformed z into an equivalent distance modulus from 9 to 17.5). As we step out in distance, the increasing volume leads to larger numbers of stars per bin, but beyond $\sim 10\,\mathrm{kpc}$ the density begins to fall faster than the volume increases, leading to a diminution of the counts. In Figure 10(b), we have normalized the distribution in each distance interval to the peak metallicity value. This nicely shows the progression of metallicity with distance, from the inner thin disk that peaks in the fifth most metal-rich bin ($[\mathrm{Fe}/{\rm{H}}]=-0.45$, maroon) to the thick disk that peaks in the seventh bin ($[\mathrm{Fe}/{\rm{H}}]=-0.65$, red) to the halo that peaks in the fourteenth bin ($[\mathrm{Fe}/{\rm{H}}]=-1.35$, green).

Zoom In Zoom Out Reset image size

Panels (d) and (e) of Figure 10 show the same information as (a) and (b), respectively, but for the intermediate latitude sample. Since the heliocentric distances remain equal, we display a closer vertical distance range from $8\lt {\mathrm{DM}}_{z}\lt 16.5$. While the number of metal-rich (disk) stars is substantially higher than in the NGC sample, it is interesting to note the striking similarity of the normalized metallicity–distance distributions (b and e). We argue in Paper I, however, that the metal-rich component observed at intermediate latitude is probably dominated by the outer disk population (Haywood et al. 2013), with a negligible contribution from the thick disk beyond $R=11\,\mathrm{kpc}$.

Beyond $\sim 6\,\mathrm{kpc}$ ($m-M=13.9$), the populations are predominantly metal-poor, yet interestingly, there remains a significant metal-rich tail at these high extraplanar distances, which we shall attempt to quantify shortly. However, it is first useful to examine whether or not there is a significant variation of the stellar populations with distance. This is explored in Figure 11, where we show in panels (a) and (b) the density profiles at $14.5\lt {\mathrm{DM}}_{z}\lt 17.5$ for the NGC and intermediate latitude samples, respectively. The stellar populations are displayed in five different metallicity slices, chosen to have approximately equal counts, and hence similar noise properties. The profiles are also normalized, so that their peak values equal unity. The distance profiles of the metallicity samples in each panel are strikingly similar, with the exception of the most metal-rich selection ($-1.0\lt [\mathrm{Fe}/{\rm{H}}]\lt 0.0$) in the NGC region. This similarity indicates that the mix of stellar populations stays approximately constant with distance for the halo population, which dominates the counts at ${\mathrm{DM}}_{z}\gtrsim 13$ (we shall return to this point below), and implies a lack of a vertical metallicity gradient in the halo. The underlying reason for the discordant metal-rich profile in Figure 11(a) is currently unclear.

Zoom In Zoom Out Reset image size

Given the opportunity afforded by this powerful new data set, we decided to investigate whether the metallicity–distance distributions discussed above could be decomposed simply into sub-populations with a minimum of assumptions. The test hypothesis that we investigate is that the NGC and the "intermediate latitude" areas of sky possess three distinct stellar populations that each have a different density profile. Allowing for only two distinct components gives a very poor fit to the present data set: the resulting residuals map in the metallicity–distance plane possess large coherent clumps, indicative of an insufficiently flexible model. Of course, we could also have chosen to employ four or more components, and could examine the Bayesian evidence for the additional parameters, but we feel that it is beyond the scope of the present work and such a study will be deferred to a later contribution.

We introduce a mild prior that the density falls off monotonically with distance z. Based on the discussion above pertaining to Figure 11, we assume that the metallicity distribution function (MDF) of each population does not vary with z, and we further assume that each MDF falls off monotonically from a single peak.

We developed an algorithm to fit the 3D distributions shown in Figure 10(a) that uses a Markov chain Monte Carlo (MCMC) procedure to search for the binned MDFs and density functions. The input data are counts at each of 34 × 25 = 850 independent bins, and the algorithm attempts to reproduce these by adapting the MDFs of the three populations (24 parameters for each MDF; not 25 since we impose the requirement that the MDFs are normalized), and the density at each of the 34 distance bins. Thus the total number of adjustable parameters is 3 × (24 + 34) = 174. We employ the same MCMC driver software presented in Ibata et al. (2011), which uses the affine sampler method of Goodman & Weare (2010). A total of 10⁷ MCMC iterations were run. (Some further technical details on the likelihood function and plausible convergence to the global optimum solution are provided in the Appendix.)

As mentioned above in Sections 3 and 4, the I08 color cut (our "wide cut") may suffer from contamination by giants in the low-distance, metal-poor bins. And, as we showed in Figure 7, this problem can be alleviated by imposing a more stringent $g-r$ selection ("narrow cut"), but at the cost of a much smaller sample size and a loss of sensitivity to metal-rich stars. In order to avoid the bias of the "wide cut," we imposed a window function on the fits, effectively ignoring those stars with $[\mathrm{Fe}/{\rm{H}}]\lt -1.5$ and distance modulus $\lt 13.5$. This window function is not necessary for the "narrow cut" and was not used when fitting to that sample.

Figure 12 shows the resulting best-fit solution for the "wide cut" sample toward the NGC. The three populations identified by the algorithm are displayed in Figure 12(a); this has picked out a very peaked metal-rich population (red line), but which contains a non-negligible metal-poor tail. The intermediate population (green line) shows a similar behavior, but displaced toward metal-poor values. Finally, a halo-like population with peak metallicity $[\mathrm{Fe}/{\rm{H}}]=-1.35$ (blue line) is also identified. The lighter lines in this panel (and also panels c, e, and f) show the 99% confidence intervals found by the MCMC parameter exploration. Marginalizing over distance gives the counts shown in Figure 12(b), where one can appreciate that the fitted model gives an excellent representation of the data. The distance distributions of the three automatically identified populations are shown in Figure 12(c); here one sees that the metal-rich population (red line) is dominant until a distance modulus of ${\mathrm{DM}}_{z}\sim 11$, and that the metal-poor halo-like distribution (blue line) becomes dominant at ${\mathrm{DM}}_{z}\sim 13$ and beyond. Marginalizing the data and model over metallicity (Figure 12(d)) demonstrates that the model also works very well in distance.

Zoom In Zoom Out Reset image size

To convert these number counts as a function of distance into a measure of population density, we need to correct for the observational selection function imposed on the sample. To this end, we use the PARSEC models (Bressan et al. 2012) to derive artificial catalogs in SDSS filters of the stellar populations in each of our 25 metallicity bins. The catalogs are shifted in magnitude to each distance interval, and we apply the same photometric selection criteria to the artificial catalogs as to the data. At this stage, we use completeness functions measured in Section 4 to randomly filter out entries in the artificial catalog, thus allowing us to account for the effect of the photometric quality selection criteria. Also, since we are analyzing the population density as a function of extraplanar distance z (and not heliocentric distance), we need to take the Galactic latitude of the survey stars into account; this is implemented by calculating the correction functions for the appropriate latitude of the stars.

By comparing the number of artificial stars that are finally detected in a metallicity–distance interval compared to the number initially generated, we derive the factors necessary to correct for the members of the stellar population that lie outside of the photometric selection window. We assume ages of 5, 10, and $12\,\mathrm{Gyr}$ for Populations A, B, and C, respectively. Note that for this calculation the age of the population is not very important: for instance, at $z=5\,\mathrm{kpc}$, the difference in correction between a $5\,\mathrm{Gyr}$ model and a $10\,\mathrm{Gyr}$ model is only 20%. In Figure 12(e), we display the relative density of the three populations, where we have corrected for the missing stars in this way. The density of the three models decreases with distance, with the metal-rich population showing the fastest fall-off.

The profile of the metal-poor Population C (blue line) does not appear to follow a power-law behavior. This can be seen when we convert the abscissa to Galactocentric values (dark cyan line), which is clearly not straight in this log–log representation. Indeed, it appears to possess a break at approximately $20\,\mathrm{kpc}$, with the inner region following an approximate power law of exponent $\alpha =-4.2$ (cyan dashed line). Interestingly, the alternative log-linear view of this population in Figure 12(f) shows that it is close to exponential over a very large range in z, but possessing a huge scale height of $3.5\,\mathrm{kpc}$. As mentioned above, for the "wide cut," we expect the metal-poor, low-distance bins to be more contaminated by giants and subgiants, and this together with low number statistics probably accounts for the upturn in Population C in the first few distance bins.

An exponential function with scale height of $0.388\,\mathrm{kpc}$ (dashed red line) can be seen to follow the Population A profile closely out to $\sim 2.5\,\mathrm{kpc}$, equivalent to ∼6.5 scale heights. Beyond that distance, there is effectively no information about this population (as can be appreciated from Figure 12(c)). Population B has a larger scale height of $0.857\,\mathrm{kpc}$ (green dashed line) and extends out smoothly to beyond $\sim 7.5\,\mathrm{kpc}$.

Finally, Figure 10(c) shows the residuals between the data and the fitted model. Comparison to Figure 10(a) shows that these residuals are typically at the ∼2% level. Thus one can obtain a remarkably good fit to the stellar populations in the NGC out to $\sim 30\,\mathrm{kpc}$ using a very simple three-component model where the stellar populations do not change with z.

In Figure 13, we show an identical analysis for the intermediate latitude sample ($30^\circ \lt | b| \lt 70^\circ$), while the residuals of this model from the data are shown in Figure 10(f). This decomposition of a completely spatially independent sample identifies almost exactly the same stellar populations (compare panel a of both Figures 12 and 13) and the resulting density profiles are similar (panels f).

Zoom In Zoom Out Reset image size

We now check the decomposition using the "narrow cut," which should not suffer as much contamination in the low metallicity and low-distance bins. In Figure 14(a), we display the distance–metallicity distribution, this time for the CFIS-u observations at $| b| \gt 30^\circ$. Qualitatively, the distribution is similar to the "wide cut" distribution for $30^\circ \lt | b| \lt 70^\circ$ shown previously in Figure 10(d), but possessing a smaller fraction of metal-rich stars, just as expected from inspecting the PARSEC models (Figure 8). The corresponding decomposition is shown in Figure 15 (note that, in this case, we do not need to employ a window function to fit the data). While we now expect the profiles of the metal-rich components to be less secure than for the "wide cut," the metal-poor components should be more reliable. Interestingly, however, we again discover a clear halo component with an exponential profile of scale height $3.1\,\mathrm{kpc}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some caution is needed not to over-interpret these data. Our main concern is the poorly known giant and sub-giant contamination in our samples. If this contamination is a constant factor, as suggested by the Besançon model (Figure 7), then the results discussed above would hold without any further correction. But any complex contamination profile in distance would introduce errors into the results quoted above. We expect that this issue will be resolved by checking our results against Gaia measurements, and bootstrapping outwards in distance.

Also, while it is straightforward to attribute the various components at the NGC to the thin disk, thick disk, and halo, the interpretation is likely more complex at lower latitudes, especially toward the Galactic anticenter. This is because the thick disk is essentially absent in the outer disk (Hayden et al. 2015), while the properties of the thin disk itself are changing rapidly at $7\lt R\lt 9\,\mathrm{kpc}$. These issues are discussed further in Paper I.

In Paper I, we introduced the u-band component of the CFIS, a community project on the Canada–France–Hawaii Telescope that aims to secure part of the ground-based photometry necessary to measure photometric redshifts for the Euclid space mission. CFIS was designed to contribute significant stand-alone science in addition to being essential for the success of Euclid. It is composed of an excellent image quality r-band survey over 5000$\,{\deg }^{2}$ whose main scientific driver is gravitational lensing, while the u-band survey of 10,000$\,{\deg }^{2}$ aims primarily to study Galactic archeology. The contribution to Galactic science will be achieved by greatly improving the metallicities and distances of faint stars in the Milky Way, thus providing an important complement to the SDSS, PS1, and Gaia surveys. The present analysis is based on approximately one-third of the final u-band area (∼2900 ${\deg }^{2}$).

The main aim of the present contribution has been to lay out in detail the procedure we use to measure accurate metallicities using CFIS u-band photometry together with $g,r,i$-band photometry derived from the SDSS and PS1 surveys. Our method is a variant of the technique developed by I08, which was greatly limited by the photometric quality of the SDSS in the u-band. By training our fitting functions (in multi-dimensional color space) on a sample of spectroscopic stars from the Segue survey, we find a scatter of 0.2 dex between the photometric and spectroscopic $[\mathrm{Fe}/{\rm{H}}]$ measurements, covering a metallicity range between solar and $[\mathrm{Fe}/{\rm{H}}]\sim -2.5$. This opens up the possibility of mapping out the chemical properties of distant stellar populations in the Milky Way (especially in its halo) with an unprecedentedly large sample of stars. The metallicity also allows for improved photometric distance measurements that will be substantially better than Gaia parallaxes for faint distant stars (which of course are the most numerous halo tracers), and will even allow Gaia's proper motion measurements for faint stars to be converted into physically more useful transverse velocities. As we discussed in Section 1, this is essential to enable a wide range of halo science questions to be addressed with Gaia.

A significant concern with this photometric metallicity method, and with the resulting photometric parallaxes, is that unresolved binaries can introduce biases into the analysis. Such pairs will of course appear brighter than isolated stars at the same distance, and one mistakenly attributes a closer distance to them. The simulations performed by I08 suggest that the worst-case binary configuration as far as metallicity determination is concerned would lead to a low-metallicity primary having its metallicity overestimated by 0.2 dex. In Figure 4, one sees such a scatter toward higher photometric metallicity, so binaries may be one of the contributors to that (slight) bias. The effect of stellar multiplicity on photometric parallax determinations is discussed in detail in Jurić et al. (2008), who modeled the consequence of various binary fractions on the derived scale heights of the Galaxy, and found that with a 100% binary fraction the scale height is underestimated by 25%. This bias must also be present in our analysis, but it remains very difficult to quantify due to the unknown binary fraction and how this property varies spatially through the components of the Galaxy.

Already, with approximately one-third of the final survey area, the CFIS-u survey provides substantially better statistics on the metallicity and distance properties of distant Galactic halo stars than has been available from previous surveys. For instance, CFIS-u has allowed us to measure good metallicities (with approximately 0.2 dex uncertainty) for $\gt {10}^{6}$ stars beyond a heliocentric distance of $4.7\,\mathrm{kpc}$.

Examining the spatial distribution of the survey stars, we find that beyond $z=8\,\mathrm{kpc}$ above the disk, and out to the limit of the survey at about $30\,\mathrm{kpc}$, the stellar populations retain an approximately constant metallicity distribution with z, implying that the population is dynamically well-mixed. This stands in contrast to what is observed in NGC 891 (Ibata et al. 2009), the only external galaxy where it has been possible to measure metallicity variations in the halo component at comparable distances.

The greatly enhanced statistics of metallicity and distance measurements at large extraplanar heights allow us to consider undertaking a decomposition of the Milky Way without employing traditional fitting methods that rely on analytic density models. To this end, we developed a non-parametric decomposition algorithm that has almost complete liberty to alter the density profile of the populations, but is subject to the reasonable constraint that the corresponding metallicity distributions have a single peak. We stress that this method was developed to allow these excellent quality data to speak for themselves, with virtually no a priori assumptions applied to the modeling, and, in particular, with no analytical profiles assumed or imposed on the Galactic sub-components (we fitted power-law and exponential profiles to the solutions, not to the data). The decomposition into three populations with MDFs that are invariant with extraplanar distance z clearly identifies a thin disk, a thick disk and a halo component toward the NGC, and in our intermediate latitude sample $30^\circ \lt | b| \lt 70^\circ$. These are recovered when using similar selection criteria to those adopted by I08 (our "wide cut"), as well as when using a stricter selection (the "narrow cut") that should be less affected by contamination from giants and subgiants. We refrain from extending the decomposition to lower latitude due to the complex behavior of the outer Galactic disk, which is discussed in Paper I.

Curiously, the halo population possesses a close to exponential profile (with a scale height of $\sim 3\,\mathrm{kpc}$) over the distances currently probed. This stands in contrast with earlier work that had found a more gentle fall-off of the halo population with distance: for instance Robin et al. (2000) fitted a power-law exponent of −2.44 to star counts over a sample of (by today's standards) small fields, while in their analysis of blue horizontal branch stars and blue stragglers in the SDSS, Deason et al. (2011) found a shallow power-law slope of approximately −2.3 inside a break radius of $\sim 27\,\mathrm{kpc}$. But perhaps the most surprising discrepancy comes from the comparison with the analysis of Jurić et al. (2008), who studied the density profile of halo MS stars in the SDSS and derived a power-law index between −2.5 and −3, which they say is "in excellent agreement with Besançon program values" (which adopts the Robin et al. 2000 profile). The disagreement with our results is all the more striking since we also analyze MS stars in the northern sky, which almost all lie within the SDSS footprint (in our case, over an area of 2900$\,{\deg }^{2}$ versus 6500$\,{\deg }^{2}$ analyzed by Jurić et al. 2008).

The improvements in the present work include the use of more accurate $g,r,$ and i photometry (being the combination of SDSS with the more precise PS1), but most importantly we now have access to much better u-band photometry from CFIS, which opens up the dimension of metallicity for a large number of halo stars out to $\sim 30\,\mathrm{kpc}$. The discrimination afforded by metallicity is important, as we show in Figure 16, which compares our observations and the Besançon model at $b\gt 80^\circ$. There we select stars with $[\mathrm{Fe}/{\rm{H}}]\lt -1.5$, which we have shown should be completely dominated by halo stars at distances $\gtrsim 5\,\mathrm{kpc}$, and we further impose the "narrow cut" color selection, so as to eliminate worries about contamination by giants. In Figure 16, we show the Besançon model version (with Robin et al. 2003 Galactic parameters) that Jurić et al. (2008) find good agreement with. The model underpredicts the counts at ${g}_{0}\lesssim 19.5$, yet overpredicts the fainter stars. As one can see in Figure 9(a), the halo component becomes particularly important at ${g}_{0}\gtrsim 18.5$ (notice the vertical feature between $0.2\leqslant {(g-r)}_{0}\leqslant 0.35$, ${g}_{0}\gt 18.5$). This is clear evidence that the density profile of the dominant halo MS population is much steeper over the heliocentric distance range of 5–$30\,\mathrm{kpc}$ than deduced by those earlier studies.

Zoom In Zoom Out Reset image size

It is possible that this approximately exponential inner halo structure was formed from the heating of the Galactic disk by minor mergers (as seen in the models of Purcell et al. 2010), which may also explain the presence of (a small number of) metal-rich stars at high extraplanar distances. At the end of the Purcell et al. (2010) simulations, the thicker component that formed in the merger had a scale height of 4–$7\,\mathrm{kpc}$, similar to our findings. Firm conclusions on this possibility will require a combined kinematic analysis with Gaia proper motions. Indeed, we expect the full power of the CFIS data for Galactic archeology science to be realized when they are coupled with these proper motion measurements.

In a future contribution, we expect to be able to recalibrate the photometric distance–metallicity relation for MS stars that was presented here, using bright Gaia stars with well-measured trigonometric parallaxes. It will be fascinating to test whether photometric distance accuracies of ∼5% can be achieved, as claimed by I08, since that will greatly enhance the power of the dynamical analyses that can be undertaken with these numerous halo tracers.

We thank the staff of the Canada–France–Hawaii Telescope for taking the CFIS data and their extraordinary support throughout the project. We are especially indebted to Todd Burdulis for the care and dedication given to planning and observing this survey. R.A.I. and N.F.M. gratefully acknowledge support from a "Programme National Cosmologie et Galaxies" grant.

This work is based on data obtained as part of the CFIS, a CFHT large program of the National Research Council of Canada and the French Centre National de la Recherche Scientifique. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA Saclay, at the Canada–France–Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS website is http://www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE).

This section aims to explain in more detail the procedure adopted in Section 6 to fit the metallicity–distance information. The MCMC algorithm we employed explores the parameter space, trying to find the optimal fitting parameters θ that maximize the log-likelihood function

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }(\theta )=\displaystyle \sum _{i=1}^{n}-\displaystyle \frac{{({D}_{i}-{M}_{i}(\theta ))}^{2}}{2\,\delta {{D}_{i}}^{2}}+\mathrm{ln}{{ \mathcal L }}_{\mathrm{priors}}(\theta ),\end{eqnarray} \tag{ 6 }$$

where D_i is the ith data point (out of n) with corresponding uncertainty $\delta {D}_{i}$, and ${M}_{i}(\theta )$ is the model calculated at the position of datum i.

As explained in Section 6, each of the three population models possesses 34 density parameters and 24 metallicity parameters for a total of 174 parameters. In general, finding the optimal configuration of a nonlinear model with 174 parameters is of course very challenging. However, the task we are confronted with here is substantially easier than the general case due to several simplifying properties of the problem. First, we assumed that the metallicity distribution of any given population does not vary with distance. This, combined with the fact that at large extraplanar distances there is only a single population present (the halo), means that the algorithm can rapidly converge on the halo MDF. Furthermore, close to the Sun, another population is dominant (the thin disk), which again greatly simplifies the task of finding the best population decomposition.

Uniform priors were adopted for all parameters. We penalized very heavily negative densities, negative MDFs and multiple peaks in the MDF. To favor plausible solutions with a monotonically decreasing density profile, we also imposed a prior such that if the model density ${\rho }_{j+1}$ in the distance bin $j+1$ is greater than the density ${\rho }_{j}$ in bin j, we add $-10({\rho }_{j+1}-{\rho }_{j})$ to $\mathrm{ln}{{ \mathcal L }}_{\mathrm{priors}}$ (using ρ in units of $\mathrm{stars}\,{\mathrm{kpc}}^{-3}$).

The initial values of the population density parameters in the MCMC runs were assigned (arbitrarily) a uniform value of 100 counts in each bin. Given the constraint that the MDFs have to be unimodal, we found it convenient to start off the MCMC runs with broad Gaussian MDFs. However, we found that the solutions converged to the same results, within the uncertainties, for all the initial MDF guesses and initial density values we tried.

Nevertheless, it is very difficult to demonstrate conclusively that the solutions presented in Section 6 are indeed the optimal most-likely solutions over the entire huge parameter space. So instead, we will examine a much simpler model, and show that this exhibits the same behavior as the full (essentially non-parametric) method.

To this end, we assume that the MDFs of the stellar populations can each be described with a skewed Gaussian function of the form

$$\begin{eqnarray}\mathrm{MDF}([\mathrm{Fe}/{\rm{H}}]) & = & \left(1-\mathrm{erf}\left(\displaystyle \frac{\alpha ([\mathrm{Fe}/{\rm{H}}]-\mu )}{\sqrt{2}\sigma }\right)\right)\\ & & \times \exp \left\{-\displaystyle \frac{{([\mathrm{Fe}/{\rm{H}}]-\mu )}^{2}}{2{\sigma }^{2}}\right\},\end{eqnarray} \tag{ 7 }$$

where α is a skewness shape parameter, μ is a location parameter, and σ is a width parameter. Given the results presented in Section 6, the density profiles of the populations are modeled as vertical exponentials:

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}\exp \{-| z| /s\},\end{eqnarray} \tag{ 8 }$$

where ${\rho }_{0}$ is the density normalization and s is the scale height.

Thus this simpler model is the sum of three such populations, each with five parameters (α, μ, σ, ${\rho }_{0}$, s), for a total of 15 parameters. The software and likelihood function we used to explore this parameter space was essentially identical to that used for the full model. As before, we adopt uniform priors on all parameters, but require σ, ${\rho }_{0}$, and s to be positive. We ran 100 simulations, each starting with random initial values chosen uniformly in the ranges of $\alpha \in [-1,1]\,\mathrm{dex}$, $\mu \in [-2.5,0]\,\mathrm{dex}$, $\sigma \in [0.2,0.5]\,\mathrm{dex}$, ${\rho }_{0}\in [{10}^{3},{10}^{5}]\,\mathrm{stars}\,{\mathrm{kpc}}^{-3}$ and $s\in [0.2,1]\,\mathrm{kpc}$. The simulations were run for ${1.1}^{6}$ iterations, and we discarded the first 10⁵ burn-in iterations. At each iteration, the high, medium, and low-metallicity solutions were assigned to Populations "A," "B," and "C," respectively. The resulting MDFs (together with their uncertainties) are shown in the top panel in Figure 17, and they can be seen to be very similar (but not identical) to the non-parametric solutions in Figure 12(a). The correlations between the density parameters are displayed in the "corner plot" below. The recovered scale-heights of the three populations can be seen to have similar values to the exponential functions overlaid on the profiles in Figure 12(f). These similarities demonstrate that our results with the non-parametric method are close to global optimal solutions for simple models.

Zoom In Zoom Out Reset image size