THE SHAPE AND PROFILE OF THE MILKY WAY HALO AS SEEN BY THE CANADA–FRANCE–HAWAII TELESCOPE LEGACY SURVEY

Studies of the Galactic stellar halo set constraints on the formation history of the Milky Way and galaxy formation processes in general. For example, contemporary simulations of galaxy formation predict that stellar halos of Milky-Way-type galaxies are assembled from inside out, with the majority of mass (50%–80%) coming from several massive (10⁸–10¹⁰ M_☉) satellites that have merged more than 9 Gyr ago, while the remaining mass comes from lower mass satellites accreted in the past 5–9 Gyr (Bullock & Johnston 2005; De Lucia & Helmi 2008). The actual fraction of massive versus less massive mergers will depend on the formation history of the galaxy in question.

A further prediction is that the stellar halos should be more centrally concentrated and should have steeper density profiles at moderate radii (>20–30 kpc) than host dark matter halos (e.g., Figure 9 in Bullock & Johnston 2005 or Figure 12 in De Lucia & Helmi 2008). According to Bullock & Johnston (2005, p. 942), "the difference in profile shapes—and the steep rollover in the light matter at moderate to large radii—is a natural consequence of embedding the light matter deep within the dark matter satellites: the satellites' orbits can decay significantly before any of the more tightly bound material is lost." Hence, they anticipate a correlation between the extent of the stellar halo (steepness of the density profile) and the extent (or mass) of satellites that built the halo: less extended (more massive) satellites will build more concentrated stellar halos. Therefore, the contribution of massive mergers to the formation of the Milky Way halo can be constrained by characterizing the stellar halo number density profile over large distances and over a wide sky area.

The large area covered by the Sloan Digital Sky Survey (SDSS; York et al. 2000), with accurate photometric measurements (∼0.02 mag; Ivezić et al. 2004) and faint flux limits (r < 22), allowed for a novel approach to studies of the stellar distribution in the Galaxy. Using a photometric parallax relation appropriate for main-sequence stars, Jurić et al. (2008, hereafter J08) estimated distances for a large number of stars and directly mapped the Galactic stellar number density to heliocentric distances of 20 kpc. They found that the halo stellar number density distribution within 20 kpc of the Sun can be fit with a two parameter, single power-law ellipsoid model

$$\begin{equation} \rho (R,Z) \propto [R^2 + (Z/q)^2]^{n/2}, \end{equation} \tag{ 1 }$$

where R and Z are the cylindrical galactocentric radius and height above the Galactic plane, respectively, n = −2.77 ± 0.2 is the power-law index, and q ≡ c/a = 0.64 ± 0.1 is the ratio of major axes in the Z- and R-direction, indicating that the halo is oblate (flattened in the Z-direction). However, additional data suggest that the J08 single power-law halo cannot be extrapolated beyond 20 kpc. A kinematic analysis by Carollo et al. (2007, 2010) suggests that the halo consists of two components with different spatial density profiles and median metallicities, with the "inner" to "outer" halo transition happening at 15–20 kpc, and the density profile becoming shallower beyond that point. On the other hand, the distribution of RR Lyrae stars from the SEKBO survey (Keller et al. 2008), and RR Lyrae and main-sequence stars from SDSS stripe 82 data seem to show a steeper density profile beyond 30 kpc (main-sequence stars can be detected up to 40 kpc in co-added SDSS stripe 82; Sesar et al. 2010).

While these studies indicate a change in the halo density profile, each has its shortcomings: kinematic studies do not give a direct measurement of the density profile but rather model it under an assumed dark matter halo potential, RR Lyrae stars are relatively sparse tracers of stellar number density (∼5 kpc⁻³ in the solar neighborhood; Sesar et al. 2010 and references therein), and the SDSS stripe 82 region covers only about 1% of the sky. Ideally, the halo stellar number density distribution should be mapped using main-sequence stars over a large fraction of the sky and to a distance of at least 50 kpc. Such studies will be enabled by next generation wide-area surveys such as the Dark Energy Survey (Lin et al. 2009), Pan-STARRS (Kaiser et al. 2002), and Large Synoptic Survey Telescope (LSST; Ivezić et al. 2008b; LSST Science Book 2009).

Meanwhile, the Canada–France–Hawaii Telescope Legacy Survey (CFHTLS) has observed 170 deg² of sky in four fields as part of the CFHTLS "wide" survey.⁶ The names, positions, and sky coverage of CFHTLS "wide" fields are listed in Table 1. Due to their relatively small sky coverage, these fields are effectively "pencil-beam" surveys for the purposes of this paper. Despite the small sky coverage (comparable to SDSS stripe 82), these data are very useful because of their depth (95% completeness at i' = 23.5 for point sources; Goranova et al. 2009), corresponding to a distance limit of ∼35 kpc for main-sequence stars, and because they observe lines of sight unexplored by other surveys. These properties allow one to study how the halo density profile changes as a function of distance and line of sight, both in the inner and outer halo.

This paper is organized as follows. In Section 2, we give an overview of CFHTLS data and describe the synthesis of SDSS ugri magnitudes from recalibrated CFHTLS u*g'r'i'i'₂ observations. The recalibration and transformation of CFHTLS data into the SDSS photometric system (Fukugita et al. 1996) is done to allow the usage of CFHTLS data with relations already defined on the SDSS system, such as the color–luminosity and photometric–metallicity relations (Ivezić et al. 2008a; Bond et al. 2010). The CFHTLS z'-band observations are not recalibrated as they were not publicly available for all CFHTLS fields at the time of writing. In Section 3, we analyze the distribution of stellar counts as a function of position and compare it to the J08 halo model. The best-fit broken power-law model is derived in Section 4. In Section 5, we study the metallicity distribution in the halo and analyze the dependence of best-fit model parameters on the adopted metallicity distribution. We finish by discussing our conclusions in Section 6.

We perform several data quality tests and post-processing steps before using the CFHTLS data in scientific analysis. First, we photometrically recalibrate the data by utilizing repeated CFHTLS observations and then transform them to the SDSS photometric system. We also investigate the performance of different types of CFHTLS magnitudes.

2.1. Overview of CFHTLS Data

2.2. Star–Galaxy Separation

We separate point-like and extended sources using the half-light radius (HLR) measured (in pixel units) for each detected source by SExtractor. For point-like sources, the HLR is independent of magnitude and depends only on image seeing (Schultheis et al. 2006). For each CFHTLS band and pointing, we remove the dependence of HLR on image seeing by subtracting the median HLR value of bright (17.5–18.5 mag) sources. Bright sources are used in this procedure because they are dominated by point-like sources (see Figure 9 in Gwyn 2008). We find that after subtraction, the distribution of HLR values of bright sources in the r' and i' (or i'₂) bands can be modeled as a ∼0.1 pixel wide Gaussian centered at zero.

We classify a source as a star if its HLR in the r' and i' (or i'₂) bands is less than 0.2 pixels (after removing the dependence of HLR on image seeing). Sources that do not satisfy this condition are classified as galaxies. To estimate the quality of this classification, we compare it to the SDSS star–galaxy classification (Lupton et al. 2002) based on deep, co-added SDSS stripe 82 data. The star–galaxy separation in co-added SDSS stripe 82 data is reliable to at least r ∼ 23 (J. Annis et al. 2011, in preparation), and for purposes of this comparison we consider it to be the ground truth.

Figure 1 shows the fraction of SDSS stars identified as stars in the CFHTLS data (completeness) and the fraction of CFHTLS stars identified as galaxies by the SDSS (contamination) as a function of r' magnitude. Using this plot, we estimate that the observed number counts will be underestimated by about 5% for r' < 21 and overestimated by about 15%–20% at the faint end (r' ∼ 22.5).

Zoom In Zoom Out Reset image size

2.3. Recalibration and Transformation of CFHTLS Data to the SDSS ugriz Photometric System

Following Padmanabhan et al. (2008), the photometric recalibration of CFHTLS data is separated into "relative" and "absolute" calibration. This process is schematically illustrated in Figure 2. For relative calibration, we use overlaps between pointings to recalibrate all pointings within a CFHTLS field onto a single photometric system (specific to that field). This calibration is then tied to the SDSS photometric system (absolute calibration). The final step in the process is the synthesis of SDSS ugri photometry from recalibrated CFHTLS u*g'r'i'i'₂ observations.

Zoom In Zoom Out Reset image size

2.3.1. Flux Extraction

The CFHTLS catalogs provide adaptive-aperture (m_auto) and fixed-aperture magnitudes (m_aper) for sources detected by SExtractor. To compare fluxes extracted by these two methods, we select stars from the W3 field and calculate m_auto − m_aper residuals, where m stands for the u*, g', r', and i' bands. The residuals are binned in m_auto bins and median values are calculated for each bin. The dependence of median values on m_auto for the u*g'r'i' bands is shown in Figure 3. In general, we find that the median difference between adaptive- and fixed-aperture magnitudes increases linearly toward fainter magnitudes, reaching ∼0.04 mag in the i' band.

This indicates that there is a problem with one or both extraction methods. As we will demonstrate in Section 2.4, the adaptive-aperture magnitudes (m_auto) are responsible for the dependence seen in Figure 3, based on the following reasons. First, measuring adaptive-aperture magnitudes is a more complex process than measuring at the fixed aperture (for details see the SExtractor v2.5 manual⁹), making adaptive-aperture magnitudes less robust. Second, a caveat in the SExtractor v2.5 manual warns of potential problems with adaptive-aperture magnitudes when the SExtractor R_min parameter is set too low. This caveat states that "when signal to noise is low, it may appear that an erroneously small aperture is taken by the algorithm. That is why we have to bound the smallest accessible aperture to R_min." Therefore, if R_min was set too low during flux extraction, inadequate (too small) apertures may have been used, making m_auto fainter and causing m_auto − m_aper residuals to be biased toward positive values. Even though the caveat states that this may happen at fainter magnitudes (low signal to noise), Figure 3 seems to indicate that brighter magnitudes may be affected as well. We therefore use fixed-aperture magnitudes in the rest of this work.

Zoom In Zoom Out Reset image size

2.3.2. Relative Calibration

Relative calibration places all pointings within a CFHTLS field onto a single photometric system (specific to that field) using repeatedly observed stars that are located in regions where pointings overlap (CFHTLS pointings overlap about 3'' in right ascension and 4'' in declination directions). Mathematically, the problem of relative calibration can be described as a χ² minimization problem (Padmanabhan et al. 2008) with

$$\begin{equation} \chi ^2 = \sum ^{n_{{\rm star}}}_i \chi ^2_i = \sum ^{n_{{\rm star}}}_i \sum _{j\in \mathcal {O}(i)} \left[ \frac{m_i - m_{j,\rm orig} - \Delta ^{{\rm rel}}_{j,f}}{\sigma _j} \right]^2, \end{equation} \tag{ 2 }$$

where n_star is the number of unique, repeatedly observed stars in the f = W1, W2, W3, W4 field, j runs over multiple observations (pointings), $\mathcal {O}(i)$, of the ith star with unknown true magnitude m_i, m_orig and σ are the magnitude and its error supplied by MegaPipe, respectively, and Δ^rel_{j, f} is the correction needed to place the jth pointing on the f field's relative photometric system. The system defined with Equation (2) is overdetermined, since the number of unknowns (n_star + n[parameters]) is smaller than the number of observations (n_obs is at least 2n_star). Equation (2) can be expressed in matrix form and can be solved using sparse matrix techniques (see Section 3.2 in Padmanabhan et al. 2008).

To solve Equation (2) for each field, we use repeatedly observed stars with 17 < m_orig < 23, where m_orig = u*, g', r', i', i'₂. The observations in this magnitude range are not saturated and their reported photometric errors are smaller than 0.2 mag. The observations are supplied to a sparse matrix inversion code written by Padmanabhan et al. (2008), and the best-fit Δ^rel_{j, f} values returned by this code are renormalized so that the median value of Δ^rel_{j, f} is equal to zero. On average, we find the rms scatter of Δ^rel_{j, f} values to be between 0.02 and 0.03 mag, reflecting the systematic uncertainty in CFHTLS magnitudes.

2.3.3. Absolute Calibration

Following relative calibration, we tie relative photometric systems to the SDSS photometric system using CFHTLS stars matched to SDSS stars. The goal is to find the best-fit Δ^abs_f value for each field, such that

$$\begin{equation} \chi ^2 = \sum _{f=W1}^{W4} \sum _{i} \left[ \frac{m_{i,\rm rel}-m_{i,\rm sdss}-C\, {\rm color} - \Delta ^{{\rm abs}}_f}{\sqrt{\sigma _i^2+\sigma _{i, \rm sdss}^2}} \right]^2 \end{equation} \tag{ 3 }$$

is minimized, where i runs over CFHTLS stars in the f field matched to SDSS stars, m_rel is the CFHTLS magnitude after relative calibration (m_rel = m_orig + Δ^rel_{j, f}), and m_sdss and σ_sdss are the SDSS point-spread function (PSF) magnitude and its error, respectively. For m_rel = u*, g', r', i', i'₂, the SDSS color and band are color = u-g, g − r, g − r, r − i, r − i and m_sdss = u, g, r, i, i, respectively (colors and magnitudes are not corrected for interstellar medium (ISM) extinction).

The C color term corrects the linear dependence of $m_{i,\rm rel}-m_{i,\rm sdss}$ residuals on color caused by differences between CFHTLS and SDSS spectral response curves (see Figure 1 in Gwyn 2008 for a comparison of spectral response curves). This term reduces the scatter in $m_{i,\rm rel}-m_{i,\rm sdss}$ residuals and improves our estimate of Δ^abs_f. Since CFHTLS and SDSS spectral response curves should not change significantly with time, the C color term is simultaneously fit for all four CFHTLS fields. Note that this term is only used when estimating Δ^abs_f; linear and other higher-order color terms that model the transformation of CFHTLS bands into SDSS bands are derived in Section 2.3.4.

To determine Δ^abs_f values, we use CFHTLS stars with m_rel > 17 matched to SDSS stars with m_sdss < 21 and with the SDSS color as specified in Table 2. The best-fit C and Δ^abs_f are also listed in Table 2. Finally, we define recalibrated CFHTLS magnitudes as

$$\begin{equation} m_{{\rm cfht}} = m_{j,f,\rm orig} + \Delta ^{{\rm rel}}_{j,f} + \Delta ^{{\rm abs}}_f = m_{j,f,\rm orig} + \Delta ^{{\rm tot}}_{j,f}, \end{equation} \tag{ 4 }$$

where Δ^tot_{j, f} is the total correction for the jth pointing in the f field. The Δ^tot_{j, f} values for CFHTLS fields are listed in Table 3.

Table 2. Δ^abs_f Values for CFHTLS "Wide" Fields

			Δabsf
CFHT Band	SDSS Color Range	C	W1	W2	W3	W4
u*	1.1 < u − g < 1.4	−0.155	0.087	0.099	0.098	0.097
g'	0.3 < g − r < 1.2	−0.161	−0.017	−0.044	−0.037	−0.043
r'	0.4 < g − r < 1.2	−0.018	−0.030	−0.019	−0.027	−0.027
i'	0.3 < r − i < 1.0	−0.081	−0.020	−0.015	−0.022	−0.030
i'2	0.3 < r − i < 1.0	0.020	0	n/a	n/a	0

Download table as:  ASCIITypeset image

Table 3. Δ^total_{j, f} Values for the CFHTLS W1–W4 Fields

Pointinga	u*	g'	r'	i'
W1+0-2	0.106	−0.034	−0.034	−0.048
W1+0-3	0.105	−0.030	−0.028	−0.040
W1+0-4*	0.077	−0.023	−0.026	0.001

Note. ^aPointings with names ending with "*" have i'₂ instead of i' observations.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

2.3.4. Transformation to SDSS Bandpasses

With CFHTLS m_cfht = u*, g', r', i', i'₂ observations calibrated onto the SDSS system, we now derive the equations that transform CFHTLS observations into SDSS m_sdss = u, g, r, i, i magnitudes. In general, the transformation from CFHTLS to SDSS bandpasses can be defined as

$$\begin{equation} m = m_{{\rm cfht}} + f\rm (color) + Z_0, \end{equation} \tag{ 5 }$$

where Z₀ is the constant term, and f(color) is some function of CFHT colors (colors and magnitudes are not corrected for ISM extinction).

To find f(color) and Z₀ for each SDSS ugri bandpass, we bin m_sdss − m_cfht residuals as a function of CFHT color, and fit polynomials to m_sdss − m_cfht medians. Here, we only use CFHTLS stars with m_cfht > 17 matched to SDSS stars with m_sdss < 20. The dependence of m_sdss − m_cfht residuals on CFHT color is shown in Figure 4.

Zoom In Zoom Out Reset image size

The best-fit polynomials shown in Figure 4 define the transformation of recalibrated CFHTLS magnitudes into SDSS magnitudes:

$$\begin{eqnarray} u &=& u^* + 1.07 - 2.375(u^*-g^\prime) + 1.954(u^*-g^\prime)^2\nonumber\\ && -\, 0.483(u^*-g^\prime)^3 \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} g &=& g^\prime + 0.05 - 0.062(g^\prime -r^\prime) + 0.365(g^\prime -r^\prime)^2 \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} r &=& r^\prime - 0.05 + 0.275(g^\prime -r^\prime) - 0.380(g^\prime -r^\prime)^2 \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} i &=& i^\prime - 0.002 + 0.092(r^\prime -i^\prime) - 0.015(r^\prime -i^\prime)^2 \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} i &=& i^\prime _2 - 0.005 + 0.145(r^\prime -i^\prime _2) - 0.280(r^\prime -i^\prime _2)^2\nonumber\\ && +\, 0.140(r^\prime -i^\prime _2)^3. \end{eqnarray} \tag{ 10 }$$

Using Equations (6)–(10), we synthesize SDSS photometry (m = u, g, r, i) for all CFHTLS observations. Note that these transformations can be used even if the relative recalibration step from Section 2.3.2 is skipped because the best-fit Δ^rel_{j, f} values are renormalized so that their median is zero.

2.4. Quality of Photometric Calibration

To estimate the quality of synthesized photometry, we use CFHT sources brighter than m = 17 that are matched to SDSS stars with m_sdss < 20, where m stands for u, g, r, and i. For each band m we calculate m − m_sdss residuals, and bin them as a function of color (not corrected for ISM extinction, color = u-g, g − r, g − r, r − i for m = u, g, r, i, respectively). The absolute value of m − m_sdss medians for the m = g, r, i bands is smaller than 0.01 mag and shows no dependence on color, for all CFHTLS fields. For the W1, W3, and W4 fields, the absolute value of the u − u_sdss medians is smaller than 0.02 mag, and shows no dependence on the u − g color, while for the W2 field the u − u_sdss medians seem to show linear dependence on the u − g color for u − g < 1.3, as illustrated in Figure 5 (bottom panel). Since all four fields were calibrated using the same procedure, and the transformations from recalibrated CFHTLS to SDSS magnitudes were derived using the data from all fields, the u − u_sdss dependence on the u − g color points to a problem with the original CFHTLS u*-band observations in the W2 field. We hypothesize that the u  −  u_sdss dependence on color may be due to incorrectly determined color-dependent air-mass term by Elixir or MegaPipe pipelines.

Zoom In Zoom Out Reset image size

The systematic uncertainty of synthesized magnitudes can be determined by comparing synthesized photometry of repeatedly observed CFHT sources; such sources can be found in regions where pointings overlap. We find that the systematic uncertainty is ∼0.03 mag (see Figure 6), which is consistent with the value cited by Gwyn (2008).

Zoom In Zoom Out Reset image size

To quantify the non-linearity of synthesized photometry, we bin m − m_sdss residuals in m = u, g, r, i magnitude bins. As shown in Figure 7, the medians of g − g_sdss residuals in g magnitude bins show no dependence on magnitude, indicating linear behavior of synthesized photometry. We have repeated this last test using adaptive-aperture magnitudes and have found a magnitude dependence in g − g_sdss residuals similar to the one shown in Figure 3. Similar results are also obtained for uri bands. These results point to the dependence shown in Figure 3 as being due to the problematic adaptive-aperture flux extraction, and justify our choice of fixed-aperture magnitudes.

Zoom In Zoom Out Reset image size

The data presented in previous sections allow us to measure the stellar number density of near main-sequence turnoff (MSTO) halo stars, and examine it as a function of position in the Galaxy. The sample we obtained from the CFHTLS extends to distances and Galactocentric radii (R_gal) nearly a factor of two greater than previous wide-area studies that used main-sequence stars (e.g., 5 kpc ≲ R_gal ≲ 15 kpc; J08). The sample also overlaps in range with studies based on RR Lyrae stars (5 kpc ≲ R_gal ≲ 110 kpc; Sesar et al. 2010), and it therefore presents an opportunity to examine the behavior of the halo density profile in the intermediate range (5 kpc ≲ R_gal ≲ 30 kpc).

We begin by selecting a sample of near-MSTO stars with the following criteria:

$$\begin{eqnarray} &&0.2 < g-r < 0.3 \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} g > 17\, \, \& \, \, 17 < r < 22.5\, \, \& \, \, i > 17 \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&5 < D/{\rm kpc} < 35, \end{eqnarray} \tag{ 13 }$$

where magnitudes and colors are corrected for ISM extinction using maps from Schlegel et al. (1998). The g − r color cut serves to select near-MSTO stars, and D is the heliocentric distance to a star computed using the photometric-parallax relation from Ivezić et al. (2008a, see Equations (A2) and (A7) in their Appendix A):

$$\begin{eqnarray} M_r &=& {-}0.56 + 14.32x - 12.97x^2 + 6.127x^3 - 1.267x^4\nonumber\\ && +\, 0.0967x^5 - 1.11\rm [Fe/H] - 0.18\rm [Fe/H]^2. \end{eqnarray} \tag{ 14 }$$

The x = g − i color in Equation (14) is computed from the more accurately measured g − r color using the stellar locus fit from J08 (see their Figure 8). The [Fe/H] is estimated from the u − g and g − r color using the photometric-metallicity relation from Bond et al. (2010, see their Equation (A1)):

$$\begin{eqnarray} \rm [Fe/H] &=& {-}13.13 + 14.09x + 28.04y -5.51xy -5.90x^2\nonumber\\ && -\,58.68y^2 + 9.14x^2y -20.61xy^2 + 58.20y^3,\nonumber\\ \end{eqnarray} \tag{ 15 }$$

where x = u − g and y = g − r. After these cuts, 13692, 7347, 6505, and 6676 stars are left in the W1, W2, W3, and W4 beams, respectively.

For each pencil beam, we bin the resulting subset in ΔDM = 0.2 mag wide bins in distance modulus,¹⁰ DM = 5log (D) − 5. We thus obtain the distribution of number counts ΔN(l, b, D) as a function of distance D for each pencil beam.

We transform the observed counts to density:

$$\begin{equation} \rho (l_i, b_i, D) = \frac{\Delta {}N(l_i, b_i, {\rm DM})}{0.2 \ln (10) D^3 \Delta \Omega \Delta {}{\rm DM}}, \end{equation} \tag{ 16 }$$

where l_i, b_i, and ΔΩ_i are the field centers and area covered by each beam as listed in Table 1, and ρ(l_i, b_i, D) is the number density in stars pc⁻³. The sight lines sampled by the CFHTLS data are shown in Figure 8.

Zoom In Zoom Out Reset image size

In panels of Figure 9, we plot the dependence of the measured number density on the distance from the Galactic center, R_gal, for each of the four CFHTLS beams. The measurements are marked by symbols with error bars and connected with a solid line for clarity. Overplotted with a dotted line on each panel is the prediction of the axisymmetric oblate halo model of J08. Plotted as open circles, and connected by red line segments, are samples within |Z| ⩽ 5 kpc of the Galactic plane. As these may be contaminated by disk stars, we leave them out of all further analyses. The residuals of the data for the J08 model are shown in Figure 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The density profile observed in W1 beam clearly stands out. While it roughly (within ∼20%) agrees with the predictions of J08 at R_gal ∼ 15 kpc, beyond that radius the observed density begins to exceed the J08 extrapolation, peaking with a factor of ∼2 excess at R_gal ∼ 28 kpc, and dropping toward the end of the observed range (R_gal ∼ 35 kpc). By comparing the location of this overdensity with the best-fit model of the Sagittarius dwarf spheroidal galaxy and its tidal tails (Law & Majewski 2010), we conclude that the excess is due to the leading and trailing arm of the Sagittarius stream (Ibata et al. 1994; Majewski et al. 2003) crossing the W1 beam. While the distribution and metallicity of stars in this beam may provide useful new constraints for the study of the stream (see Section 6), this (un)fortunate fact makes the majority of W1 data unusable for the study of the smooth halo profile.

Density profiles in W3 and W4 beams agree within ∼15% with the predictions of the J08 halo model in the 10 < R_gal/kpc < 25 and 8 < R_gal/kpc < 27 ranges, respectively. The observed agreement is nontrivial. First, the directions observed by these beams do not overlap with the SDSS data used by J08, and therefore test their model in an entirely different part of the halo (especially the W4 beam). Second, the CFHTLS data cover a significantly larger distance range than the data used by J08, thus validating the extrapolation to ∼10 kpc greater distances (and an order of magnitude change in stellar number density). And finally, while J08 did use a small (∼300 deg²) area in the southern Galactic hemisphere to construct their model, their best-fit parameters were largely determined by the ≳ 6000 deg² of data from the north. The fact that the W4 beam (b = −4184) is very well matched by the model puts a constraint to any asymmetries between the northern and southern hemispheres to ≲ 15%, out to distances of at least R_gal ∼ 25 kpc. This conclusion is also supported to distances of D ≲ 18 kpc by the analysis of north versus south SDSS III wide-area imaging (A. Bonaca et al. 2011, in preparation).

Beyond R_gal ∼ 25 kpc, in both the W3 and W4 beams, the J08 model overpredicts the observed counts. In particular, the J08 model overpredicts the density observed in the final two bins of the W3 beam by 20% and 40%, respectively. Since we have measured incompleteness to be on the order of ∼5% in almost the entire observed range (see Figure 1), and have made the bins volume-complete, this turnover cannot be a result of observational bias but an indication of a change in halo density profile beyond R_gal ≳ 25 kpc. The overprediction by the model may be even greater since the observed CFHTLS stellar counts are overestimated by 15%–20% at the faint end due to the contamination by galaxies, as shown in Figure 1.

The profile exhibited by field W2 is also unusual. As shown in Figure 9, it is significantly steeper than either the W3 or W4 profile or the J08 prediction. By itself it is well described by a single n = −4.5, q = 0.65 power law. This fit, however, is incompatible with observations from the other three directions where a shallower profile closer to n ∼ −2.8 is strongly preferred. The curious behavior appears to be a combination of two effects: (1) the overdensity created by the Monoceros stream at R_gal ≲ 25 kpc and (2) the steepening of the halo profile beyond R_gal ∼ 30 kpc.

The Monoceros stream (Newberg et al. 2002) is present in the general direction of the W2 field (l = 232°, b = 26°). For example, in an AAT/WFI survey of the anticenter region, Conn et al. (2007) detect the stream in the l = 220°, b = 15° direction at D = 11 ± 1.6 kpc, as well as at l = 240°, b = 10° with D = 13.8 kpc. Similarly, J08 are able to trace the stream in SDSS star count maps to l ∼ 230° and, importantly, show that it extends to at least b ∼ 25° of Galactic latitude where it forms a factor of ∼1.5–2 overdensity with respect to the extrapolation of the smooth stellar background before and after the stream (e.g., as seen in the top left panel of Figure 13 in J08). This is consistent with the factor of ∼2 overdensity observed in W2. Second, both of these studies estimate the Galactocentric distance and width of the stream of R_gal ∼ 18 kpc and ΔR_gal ∼ 3–4 kpc in direction of W2. Given that these characteristics are broadly consistent with the enhancement in the W2 direction for R_gal ≲ 25 kpc, and having no additional evidence to the contrary, we interpret the observed enhancement as a detection of the Monoceros stream.

Taking into account the enhancements due to the Sagittarius and Monoceros streams, a single power law remains an appropriate description of the smooth halo component to Galactocentric distances of R_gal ∼ 25–30 kpc. Beyond this limit, however, the observed profile appears to turn over rather quickly and the model needs to be modified to explain it.

To assess the character of the observed turnover, we fit a series of models of increasing complexity to the observed data. We use χ² per degree of freedom (χ²_dof) as the goodness of fit metric, and search for minima in χ² hypersurfaces using a Levenberg–Marquardt non-linear solver as implemented by the GNU Scientific Library.¹¹ To increase the likelihood of finding the true global minimum, we repeat the minimization procedure with 10,000 different initial conditions selected randomly from a plausible range of initial values of each parameter.

We begin by fitting a single J08-type power law to all admissible data points¹² of beams W3 and W4, the first eight data points of the W1 beam (those that show no contamination by the Sagittarius stream), and the last six data points of the W2 beam (those that we judge are past the influence of the Monoceros stream). We obtain ρ₀ = 1.7 × 10⁻⁶ pc⁻³, q = 0.72, n = −2.9 as the best, but less satisfactory (χ²_dof = 6.8), fit. The ρ₀ is the normalization (number of stars per pc³) for the 0.2 < g −  r < 0.3 color bin we use. For comparison, the fiducial model with J08 parameters shown in Figure 9 has χ²_dof = 9.1 when fitted to the same data.

We next increase the complexity of the model by allowing for triaxiality of the ellipsoid, parameterized by w ≡ b/a (the ratio of ellipsoid axes):

$$\begin{equation} \rho (x, y, z) \propto \left(x^2 + \frac{y}{w}^2 + \frac{z}{q}^2\right)^{\frac{n}{2}}. \end{equation} \tag{ 17 }$$

This addition makes practically no difference; the best-fit model changes only slightly (ρ₀ = 1.7 × 10⁻⁶ pc⁻³, q = 0.72, n = −2.96, w = 1.02) while χ²_dof actually increases (to ∼7) because of the extra degree of freedom.

We continue by permitting the triaxial halo ellipsoid to rotate in the X–Y (Galactic) plane by an angle ϕ:

$$\begin{eqnarray} x^{\prime } & = & \cos (\phi)\, x - \sin (\phi)\, y \nonumber \\ y^{\prime } & = & \sin (\phi)\, x + \cos (\phi)\, y \nonumber \\ \rho (x, y, z) & \propto & \left(x^{\prime 2} + \frac{y^{\prime }}{w}^2 + \frac{z}{q}^2\right)^{\frac{n}{2}}. \end{eqnarray} \tag{ 18 }$$

This five-parameter model marginally improves the fit (χ²_dof = 6.8), but converges to parameters ρ₀ = 1 × 10⁻⁶ pc⁻³, q = 1.11, n = −3.3, w = 1.3, ϕ = 230° that are strongly excluded by prior data (e.g., Chen et al. 2001, J08, and others).

Given the results of this series of experiments, a single power law is unlikely to describe the observed counts: the data require a functional form allowing for a change in the radial profile beyond R_gal ≳ 25 kpc. We therefore attempt a series of simple "broken power law" models, where the density follows one (the "inner") power law until radius R_br is reached, and the other (the "outer" power law) beyond.

We begin with a minimal extension of a single power-law model, allowing for a change of the power-law index beyond a certain radius R_br:

$$\begin{eqnarray} R_e & = & \left(x^2 + y^2 + \frac{z}{q}^2\right)^{\frac{1}{2}} \nonumber \\ \rho (x, y, z) & \propto & \left\lbrace \begin{array}{@{}c@{}}(R_e)^{n_{\rm inner}}, \;\;\;\;\;\;\; R_e < R_{\rm br} \\ (R_e)^{n_{\rm outer}}, \;\;\;\;\;\;\; R_e > R_{\rm br}. \\ \end{array} \right. \end{eqnarray} \tag{ 19 }$$

Note that as the ellipsoid is allowed to be oblate or prolate, the break radius R_br is only equal to the physical Galactocentric radius R_gal on the x-axis (along the line connecting the Galactic center and the Sun). In the vertical direction, the physical radius corresponding to R_br is reduced by a factor of q⁻¹.

The above model, with five free parameters, produces a significantly better fit to the data (χ²_dof = 3.9). The best-fit parameters for the inner power law, ρ₀ = (1.45 ± 0.05) × 10⁻⁶ pc⁻³, q = 0.70 ± 0.01, n_inner = −2.62 ± 0.04, R_br = 27.8 ± 0.8 kpc, n_outer = −3.8 ± 0.1 are in excellent agreement with the J08 model, while beyond R_br = 27.8 kpc the best-fit profile becomes steeper than the J08 model.

As shown in Figure 10, the observed profiles are better fit by the broken power law than by the J08 model, excluding regions with known tidal streams. A fit that entirely excludes the W1 and W2 fields (not just the regions with known tidal streams) does not strongly constrain the broken power-law model. In this fit, the normalization ρ₀ has high (∼50%) fractional uncertainty and strongly correlates with oblateness q and break radius R_br (correlation coefficient is ∼1). This strong correlation is caused by similar positions of the W3 and W4 beams in the x–z and y–z planes with respect to the Galactic plane (see Figure 8). In comparison, our best-fit profile that uses all four beams while excluding regions with known tidal streams shows much weaker correlation between ρ₀, q, and R_br (∼0.4).

The inner parts of the W3 field in our best-fit model do show a systematic underestimate of the counts by the model on the order of ∼15% (at an approximately 1σ–2σ level¹³), while the opposite occurs in inner parts of the W1 field. These lines of sight point toward high latitudes in the northern and southern Galactic hemispheres, and the observed difference may be a signature of slight north–south asymmetry of halo star counts recently seen in SDSS III data (A. Bonaca et al. 2011, in preparation).

To further assess the robustness of the detected break, we examined two more variants of the model: the first, where we allowed the outer halo to have an oblateness parameter q different from that of the inner halo, and the second, where we fixed the parameters of the inner halo to J08 values, and allowed only those of the outer halo to vary. Both cases resulted in a similar value of χ²_dof, as well as similar break radii (R_br ∼ 28 kpc) and outer power-law indices (n_outer ∼ −3.8). Importantly, the model with varying q_inner and q_outer produced best-fit values of 0.71 and 0.69 for the two, respectively, indicating that there is no evidence for a change in oblateness of the halo across the range of distances examined.

Based on the above analysis, we conclude that the detection of steepening of the halo density power law is robust. To explain the CFHTLS data, the power-law index needs to change from n_inner = −2.62 to n_outer = −3.8 around R_br ∼ 28 kpc. An R–Z plane visualization of the J08 power-law halo model and the broken power-law model presented in this paper is shown in Figure 11.

Zoom In Zoom Out Reset image size

We also fit Einasto's model (Einasto 1965) to our data to allow easy comparison with density profiles obtained from N-body simulations (Navarro et al. 2004; Diemand et al. 2004; Merritt et al. 2005; Graham et al. 2006). The best-fit parameters for Einasto's model are n = 2.2 ± 0.2, R_e = 22.2 ± 0.4 kpc, ρ₀ = 1.06 ± 0.05 stars pc⁻³, and q = 0.70 ± 0.01 with χ²_pdf = 4.25.

Due to contamination by galaxies at the faint end (see Figure 1), the n_outer = −3.8 ± 0.1 power-law index given above is likely somewhat shallower than what it should be. To estimate by how much n_outer is shallower due to contamination by galaxies, we determine f(r) = [1 − contamination(r)]/completeness(r) for each star, where r is a star's r-band magnitude, and completeness(r) and contamination(r) are the solid and dashed lines from Figure 1, respectively. We sum f(r) values in distance modulus bins and use them instead of raw number counts when fitting the model. With this approach, a star located in a distance modulus bin with higher galaxy contamination will contribute a value smaller than one toward the total count in that bin. Finally, we find that n_outer decreases by 0.1 to n_outer = −3.9, which is still within the uncertainty of n_outer = −3.8 ± 0.1 determined previously.

In addition to the stellar number density distribution, the metallicity distribution can also provide strong constraints for halo formation models. Thanks to CFHTLS u-band data, it is possible to compute photometric metallicity using a method developed by Ivezić et al. (2008a). Furthermore, the absolute magnitude of a star and hence its distance depend on the star's metallicity (Equation (14)) and thus it is important to study systematic errors in best-fit model parameters as a function of metallicity distribution.

Figure 12 shows the median halo metallicity in four CFHTLS beams as a function of distance from the Galactic center, R_gal. In the two beams where significant substructure is not detected (W3 and W4 beams), the median metallicity is independent of distance and averages to [Fe/H] ∼ -1.5 dex with a range of 0.1 dex. This average value is consistent with the median halo metallicity measured by Ivezić et al. (2008a) using SDSS F- and G-type main-sequence stars, and the range is consistent with the systematic uncertainty inherent to this photometric metallicity method (∼0.1 dex; Ivezić et al. 2008a).

Zoom In Zoom Out Reset image size

The average metallicity in the W1 beam is also independent of distance and averages [Fe/H] ∼ -1.5 dex. This trend is probably a coincidence since the Sagittarius tidal stream passes through the beam and the stream's metallicity is not required to be [Fe/H] ∼ -1.5 dex everywhere. On the contrary, models and observations (Chou et al. 2007; Sesar et al. 2010; Law & Majewski 2010, and references therein) suggest that the metallicity exhibits a gradient along the Sagittarius tidal stream—we simply happen to observe the stream where its metallicity is [Fe/H] ∼ -1.5 dex.

The metallicity in the W2 beam is a bit higher ([Fe/H] ∼ -1.3 dex) than in other beams, even at large Galactocentric distances (R_gal > 28 kpc) where, according to Figure 10 (top right panel), the contribution of the Monoceros stream stars should be small. Since the Ivezić et al. (2008a) photometric metallicity method depends on accurate u-band photometry, anything affecting u-band measurements will also affect the photometric metallicity estimate. As discussed in Section 2.4, the u-band measurements synthesized from CFHTLS data may be impacted by some calibration issues in this beam, so these issues are likely responsible for the apparently higher metallicity in the W2 beam (a +0.05 mag systematic offset in the u band will increase the metallicity by 0.2 dex).

Our observation that the halo metallicity is independent of Galactocentric distance goes contrary to the results of Carollo et al. (2007) and de Jong et al. (2010, hereafter deJ10), who conclude that the halo metallicity changes from [Fe/H] ∼ -1.6 dex to [Fe/H] ∼ -2.2 dex at R_gal ∼ 15 kpc. To see how our best-fit model varies when the metallicity changes from [Fe/H] ∼ -1.6 dex for R_gal ≲ 15 kpc to [Fe/H] ∼ -2.2 dex beyond R_gal ∼ 15 kpc, we iteratively modify the metallicity distribution in the CFHTLS beams to reflect the metallicity distributions shown in deJ10, Figure 7 (right, hereafter [Fe/H]_deJ10(D)). New metallicities are assigned to CFHT stars by interpolating [Fe/H] from [Fe/H]_deJ10(D) using initial heliocentric distances, D (calculated from Equations (14) and (15)). Heliocentric distances are then recalculated using new metallicity values, and the metallicities are again interpolated from [Fe/H]_deJ10(D)). This process is repeated (about 2–3 times) until the fractional difference between distances in consecutive steps dips below 0.15 (fractional distance uncertainty for main-sequence stars is ∼15%; Sesar et al. 2008). The broken power-law model is fit to modified data once convergence is achieved.

The best-fit model parameters obtained for a deJ10-like metallicity distribution are listed in Table 4. In addition, Table 4 also lists best-fit model parameters obtained for constant metallicity distributions (i.e., all stars are assumed to have fixed metallicity) and for the best-fit model presented in Section 4. Within a range of plausible metallicity distributions, the best-fit model parameters do not seem to vary greatly and average around q ∼ 0.7, n_inner ∼ −2.4, n_outer ∼ −3.8, and R_br ∼ 25 kpc.

We have recalibrated CFHTLS "wide" survey u*g'r'i'i'₂ observations and transformed them to the SDSS ugri photometric system. Using a series of tests, we demonstrated that synthesized ugri observations, when compared to SDSS observations, show no dependence on color or magnitude. The only exception to this are the u-band observations synthesized from CFHTLS W2 field u*-band observations, which show a slight linear dependence on the u − g color. Median photometric error in synthesized ugri photometry is ∼0.03 mag at the bright end and ∼0.1 mag at ∼22.5 mag.

By obtaining the synthesized ugri photometry from CFHTLS observations, we were able to use the photometric parallax and metallicity relations which allowed us to study the spatial and metallicity distributions of near-turnoff main-sequence stars in the Galactic halo to heliocentric distances of ∼35 kpc. We find that the halo number density profile becomes steeper at Galactocentric distances greater than R_gal ∼ 28 kpc, with the power-law index changing from n_inner = −2.62 ± 0.04 to n_outer = −3.8 ± 0.1. While the best-fit model parameters do change slightly depending on the adopted metallicity distribution (see Table 4), we find that a broken power-law model is required for a good fit to the data. We also find that the best-fit R_br value cannot be smaller than R_br = 22 kpc even for the most extreme assumptions about the halo metallicity distribution (i.e., constant metallicity at [Fe/H] = -2.0 dex from 5 to 35 kpc).

This study provides further evidence for the steepening of the halo density profile previously detected by Sesar et al. (2010) using main-sequence and RR Lyrae stars from the SDSS stripe 82, and by Keller et al. (2008) using RR Lyrae stars from the SEKBO survey. This result is consistent with predictions of galaxy formation simulations which find a steepening of the density profile beyond ∼30 kpc (Bullock & Johnston 2005; De Lucia & Helmi 2008; Zolotov et al. 2009).

We see no evidence of change in halo metallicity within the range of probed distances. The halo metallicity ranges between [Fe/H] ∼ -1.4 dex and [Fe/H] ∼ -1.6 dex, and averages at [Fe/H] ∼ -1.5 dex. This result runs contrary to the Carollo et al. (2007) and de Jong et al. (2010) studies which report a metallicity of [Fe/H] ∼ -1.6 dex within R_gal ≲ 15 kpc, and [Fe/H] ∼ -2.2 dex beyond. Only the in situ spectroscopic metallicities of distant main-sequence stars may provide a definitive answer to this discrepancy. With the sky density of near-MSTO stars at high Galactic latitudes of about 100 stars deg⁻², the multi-object capability over a wide field of view would be well matched to such a program.

While the total sky coverage of the four CFHTLS beams is slightly smaller than that of SDSS stripe 82 (220 deg² versus 300 deg²), the CFHTLS beams provide a much stronger constraint on the oblateness (Z to R semi-major axis ratio) of the stellar halo because they probe very different Galactic lines of sight. We find the oblateness to be q = 0.70 ± 0.01 and see no evidence of change across the range of probed distances (5 < R_gal/kpc < 35). This result is quite consistent with the oblateness of the dark matter halo, q_DM = 0.72, obtained by Law & Majewski (2010) using the positions and kinematics of Sagittarius stream stars. However, we find the minor axis of the stellar halo to be aligned with the spin vector of the Milky Way, while Law & Majewski (2010) find the minor axis of the dark matter halo to be perpendicular to the spin vector of the Milky Way.

We have detected the Sagittarius and Monoceros streams as excesses of stars in the CFHTLS W1 and W2 fields, respectively. These detections provide new constraints on models of these streams. For example, the Law & Majewski (2010) model of the Sagittarius dwarf spheroidal galaxy predicts the positions and velocities of Sagittarius stream stars in the CFHTLS W1 beam (Figure 13). The spatial distribution of Sagittarius stars is predicted to be bimodal, with a narrow peak at R_gal ∼ 20 kpc, and a broader peak at R_gal ∼ 35 kpc. The observed distribution of Sagittarius stars in this region, overplotted in Figure 13 (top panel), has only one peak at R_gal ∼ 35 kpc. Unfortunately, we do not have radial velocity measurements in this region, and cannot quantify contributions of particular streams (leading and trailing) to this peak. Multi-epoch surveys, such as the Palomar Transient Factory (PTF; Law et al. 2009) and LSST, will enable robust identification of bright tracers such as RR Lyrae stars, which can be used to map the velocity structure in this region. The resulting improvement in models of these streams may then help to further constrain the shape, orientation, and mass of the Galactic dark matter halo.

Zoom In Zoom Out Reset image size

B.S. thanks NSF grant AST-0908139 to J. G. Cohen and NSF grant AST-1009987 to S. R. Kulkarni for partial support. Ž.I. acknowledges support by NSF grants AST-0707901 and AST-1008784 to the University of Washington, by NSF grant AST-0551161 to LSST for design and development activity, and by the Croatian National Science Foundation grant O-1548-2009. M.J. acknowledges support by NASA through Hubble Fellowship grant HF-51255.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. Partial support for this work was provided by NASA through a contract issued by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, 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 of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work is based in part on data products produced at TERAPIX and the Canadian Astronomy Data Centre as part of the Canada–France–Hawaii Telescope Legacy Survey, a collaborative project of NRC and CNRS.