The Saga of M81: Global View of a Massive Stellar Halo in Formation

Given the crowding- and magnitude-dependent impacts on photometric completeness shown in Figure 4, any inferences made from RGB star counts in a given region of sky, selected as described in Section 3, must account for these effects. To calibrate RGB counts specifically, we apply the stellar selection criteria and RGB selection box (given in Table 2) to both the input AST catalog and the recovered catalog to assess the correspondence between input and recovered RGB density, as a function of true source density.

Figure 5 shows the results of this analysis, with average input RGB density (in arcmin⁻²) across the AST area plotted against the average recovered RGB density, for each of the 12 ASTs. As expected from the completeness functions in Figure 4, the lower density half of the ASTs (with ≲100 stars arcmin⁻² initial density) exhibits a smooth, linear correspondence between input and recovered RGB counts. The normalization factor is 0.48—i.e., approximately 48% of input RGB-like stars are recovered in the linear density regime. However, for large input densities (≫100 stars arcmin⁻²), the correspondence between input and recovered counts deviates nonlinearly, presumably due to the increasing impact of crowding. Though we take a somewhat nonstandard approach to AST analysis, the reliable behavior of the AST results relative to expectations, as well as the consistency with empirical comparisons to existing HST observations, suggests that this is an effective approach. We fit the results with a model of the form

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{out}}={{\rm{\Sigma }}}^{* }\,{{\rm{\Sigma }}}_{\mathrm{in}}\,{e}^{-\alpha {{\rm{\Sigma }}}_{\mathrm{in}}/{{\rm{\Sigma }}}_{0}},\end{eqnarray} \tag{ 1 }$$

with best-fit parameters of Σ^* = 0.48, Σ₀ = 750, and α = 2. In Section 5.1.1, we will use these results to correct the observed RGB density to the intrinsic density in order to measure SB and stellar mass.

Zoom In Zoom Out Reset image size

Perhaps the more nuanced measurement of the observed stellar populations is that of color and in turn estimates of abundance. Our survey is optimally geared to efficiently detecting RGB stars at colors of g − i = 1–1.5. For M81, this corresponds to limiting g-band magnitudes of ∼27. However, the most metal-rich RGB stars, i.e., those with [M/H] ≫ −0.5, will have g-band magnitudes of 28–29—substantially fainter than the depths achieved by this survey. Therefore, unless g-band observations are substantially deeper than those in the i band, any metal-rich populations that might exist will be too faint to observe in this survey and all similarly designed ground-based surveys.

However within our HSC footprint, there are 17 existing HST fields (ACS and WFC3) with high-quality stellar catalogs from the GHOSTS survey (Radburn-Smith et al. 2011; Monachesi et al. 2013). These fields span the majority of the stellar density scale probed by our ASTs, ranging from ∼5–400 stars arcmin⁻². A number of these observed fields have been used to robustly measure the minor axis color (Monachesi et al. 2016a) and SB (Harmsen et al. 2017) profiles of nearby galaxies, including M81, using RGB stars detected in the F606W/F814W filters. While covering far less of the M81 Group than our HSC observations, the GHOSTS fields are ∼1 mag deeper in the i band and can effectively probe the redder stellar populations unobserved by HSC. We show in Section 5.1.1 that this difference in photometric depth does not impact our SB analysis. However, in the event that there is a substantial subset of higher-metallicity stellar populations, the shallower color sensitivity of HSC could impact our ability to accurately measure the color, and therefore metallicity, of the intrinsic halo populations.

To assess the presence of such populations, and the impact on measured color they might effect, we use the existing GHOSTS fields in combination with our AST catalogs. We take the composite CMD of 14 of the GHOSTS fields and we convert from F606W−F814W versus F814W to g − i versus i. Following Monachesi et al. (2013) and Harmsen et al. (2017), we use an [M/H] = −1.2, 10 Gyr old PARSEC (Bressan et al. 2012) isochrone model for the conversion, though we note that the g − i/F606W−F814W color–color relation for RGB stars is relatively insensitive to metallicity for halo-like populations. Using this g − i/i converted catalog of GHOSTS sources, we matched each to the nearest star in the CMD of one of our AST catalogs (of intermediate density; 80 stars arcmin⁻²). This input star catalog was then matched to the corresponding pipeline-produced source catalog, and RGB star candidates were selected following Section 3. The resulting CMD of GHOSTS–AST matched sources is shown in Figure 6 (inset right).

Zoom In Zoom Out Reset image size

Monachesi et al. (2016a) measured color profiles along the major and minor axes of the GHOSTS sample, including M81. To measure more robust colors, which also are intrinsically better tied to population changes due to metallicity, they adopt a revised color metric, Q (see also Monachesi et al. 2013). As the isochrone model curve for a metal-poor stellar population is nearly a straight line for the upper portion of the RGB, the Q color corresponds to a CMD which has been rotated around a point 0.5 mag below the TRGB, such that the RGB is nearly vertical. We adopt the Q metric for this paper as well, for all of our color-based analysis. As we are operating in the g − i filters, we define a new Q_Col corresponding to a rotation angle of −22°. An example of this redefined color schema is shown overlaid on the CMD of GHOSTS–AST matched input stars in Figure 6 (inset right).

The Q_Col distributions of both the converted GHOSTS stellar catalog and the sources recovered from the GHOSTS stars matched to the input AST catalog are shown in Figure 6. The Q_Col distribution of HSC-observed RGB star candidates shows a distinct deficit of red sources, relative to GHOSTS, amounting to an offset of 0.08 mag bluewards in the median Q_Col. This offset represents a difference in the median measured and median intrinsic colors of the halo populations. We thus take this offset into account when reporting median color measurements, and related median metallicities, throughout the rest of the paper.

We caution that without similar extensive overlap with high-quality HST-derived stellar catalogs, it is impossible to estimate the contribution from higher-metallicity stellar populations and, thus, this "blue bias" is unable to be reliably corrected for. This effect will present an issue for all similarly designed ground-based stellar population surveys at distances ≫1 Mpc.

We define the minor axis according to the region shown in Figure 3 in red. Leveraging our large survey footprint, we define a much wider minor axis region than is covered by the GHOSTS survey, allowing for more robust averaging and inclusion of any potential substructure absent in the sparse GHOSTS measurements.

We divide the minor axis into projected annular radial regions, 2 kpc wide from 10 to 40 kpc, and wider 5 kpc bins outside 40 kpc, to account for the lower number of sources. We then calculate the average density for each radial region. Visually inspecting the CMDs in each bin, we find that at radii >65 kpc along the minor axis, the RGB was indistinguishable from a CMD composed of only background galaxies and MW foreground stars (see Figure 7). We thus consider the halo beyond 65 kpc along the minor axis to be undetected and parameterize the density of RGB candidates at >65 kpc as a uniform background equivalent to 0.53 arcmin⁻². This is different from the findings of Harmsen et al. (2017), where it was assumed that stars outside of 45 kpc were background contamination, due to the limited radial range and overall number counts. Using the extended radial range provided by HSC, as well as the improved star counts provided by our broad minor axis selection region, we detect RGB halo sources out to 65 kpc, indicating that the initial GHOSTS SB profile was oversubtracted.

In order to correct our measured RGB counts for this background contamination, as well as intrinsic incompleteness due to photometric detection and source crowding, we turn to the results of our ASTs, described in Section 4.1. Using our discovery of the nonlinearity between intrinsic and recovered RGB counts and the presence of a uniform contaminating background beyond 65 kpc, we model "corrected" RGB density as a function of observed density, with the form

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{corr}}={{\rm{\Sigma }}}^{* ^{\prime} }({{\rm{\Sigma }}}_{\mathrm{obs}}\,-\,{{\rm{\Sigma }}}_{\mathrm{BG}}){e}^{\beta {\left(\tfrac{{{\rm{\Sigma }}}_{\mathrm{obs}}-{{\rm{\Sigma }}}_{\mathrm{BG}}}{{{\rm{\Sigma }}}_{0}}\right)}^{\gamma }}.\end{eqnarray} \tag{ 2 }$$

The best fit to our AST results, shown in Figure 8, gives parameters of ${{\rm{\Sigma }}}^{* ^{\prime} }=2.09$, Σ_BG = 0.53, Σ₀ = 140, β = 2, and γ = 1.6. To corroborate the validity of this model, we also compare the RGB densities in each of the GHOSTS fields, measured with Subaru, to the densities published in Harmsen et al. (2017). Shown in Figure 8, we find that after correcting for the oversubtraction of background sources, the GHOSTS measurements agree excellently with our model for corrected RGB counts. This is a powerful and completely empirical support for our chosen model and AST approach.

Zoom In Zoom Out Reset image size

We calculate the mean RGB-selected source density in each radial region and then correct to the "intrinsic" density, accounting for background contamination and incompleteness, using Equation (2). We then convert this corrected density to equivalent V-band SB following Harmsen et al. (2017), as

$$\begin{eqnarray}&&{\mu }_{V}=-2.5{\mathrm{log}}_{10}[{{\rm{\Sigma }}}_{\mathrm{corr}}({\mathrm{arcsec}}^{-2})\cdot 2.09\times {10}^{-10}],\end{eqnarray} \tag{ 3 }$$

assuming a 10 Gyr, [M/H]  =   −1.2 isochrone model.

Uncertainties on the density measurements were carefully accounted for from three distinct sources. we assume errors on the average density in each radial section by taking the standard deviation in density across all independent ∼arcmin² regions, divided by the square root of the number of these regions. Owing to the wide minor axis selection area, and therefore large number of independent arcmin² regions in each radial section, these uncertainties are quite small—ranging from ∼20% at the smallest radii, to ∼1% at the largest radii. Second, we account for the model-based correction of the RGB density from our ASTs. Owing to the smoothness of the AST results, this is a small effect. As the high-density exponential portion of the model is both the most difficult to fit and has the greatest potential to impact mass estimates, we estimate the uncertainty by performing the fit while excluding either the highest-density or second highest-density point. The difference, which we take into account, is small—approximately 1%–5%. Lastly, we calculate Poisson uncertainties on the total number of stars in each radial bin. Owing to the angular broadening of our minor axis selection region, as well as the broadening of the radial binning, with radius, these uncertainties are relatively constant at 3%–5%.

Figure 9 shows the minor axis SB profile for M81. Our measurements are shown in blue, with the GHOSTS points shown in gray for comparison. Our measurements extend ∼50% farther along M81's minor axis than GHOSTS and reach remarkable depths of μ_V  ≃  33 mag arcsec ⁻² at 65 kpc. This is among the deepest SB profiles ever measured (e.g., compare to μ_V  ∼  32 mag arcsec ⁻², PISCeS Survey, Crnojević et al. 2016; μ_V  ∼  30 mag arcsec ⁻², Dragonfly Survey, Merritt et al. 2016).

Zoom In Zoom Out Reset image size

To quantify the halo mass, we fit a power law to the density profile. As Figures 7 and 9 show, the inner two radial bins are impacted by crowding too severe to be accounted for by our AST-based corrective model. To account for this in our density fit, we combine our Subaru measurements from 15–65 kpc with the three smaller-radius GHOSTS measurements from Harmsen et al. (2017, 11.2, 12.2, and 13.4 kpc; shown as filled gray circles in Figure 9) to fit the full density profile. We use a model of the form

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\rm{\Sigma }}(R)={\mathrm{log}}_{10}{{\rm{\Sigma }}}_{0}-\alpha {\mathrm{log}}_{10}(R/{R}_{0}),\end{eqnarray} \tag{ 4 }$$

where the densities are expressed in arcsec⁻². This best-fit model yields parameters of ${\mathrm{log}}_{10}{{\rm{\Sigma }}}_{0}$ = −0.53, R₀ = 4.19 kpc, and α = −2.59. Though this is much shallower than the density profile found by Harmsen et al. (2017; α = −3.53), it is in good agreement with the GHOSTS measurements when corrected for their initially oversubtracted background (above). Following Harmsen et al. (2017), we integrate the profile from 10–40 kpc, using elliptical annuli with the same assumed projected axis ratio of 0.61, obtaining an accreted stellar mass from 10–40 kpc of M_⋆,10−40 = 3.7 × 10⁸ M_⊙. Extrapolating to total accreted mass using the Harmsen et al. (2017) 10–40-to-total ratio of 0.32, we estimate a total accreted mass of M_⋆,Acc = 1.16 × 10⁹ M_⊙–within 2% of the GHOSTS estimate of 1.14 × 10⁹ M_⊙. Despite the change in power-law parameters, the large characteristic radius, combined with the relatively narrow 10–40 kpc radial range over which the profile is evaluated, results in nearly identical total mass estimates.

Finally, we compare our resolved star-based minor axis SB profile to integrated light measurements, which excel in the bright innermost parts of the galaxy, where resolved star measurements suffer from strong crowding. Figure 10 combines our measured profile with a near-infrared version of M81's minor axis SB profile, following Harmsen et al. (2017). In this case, we have chosen the WISE W1 (3.4 μm) profile measured as part of the WISE Enhanced Resolution Galaxy Atlas (Jarrett et al. 2012; Jarrett et al. 2013; T.H. Jarrett 2020). We have adjusted the elliptically averaged profile to a minor-axis-only version using the measured axis ratio for each elliptical annulus. Then, using the same 10 Gyr, [Fe/H] = −1.2 isochrone model which was used to convert our RGB counts to μ_V, we instead convert these counts to W1. The WISE profile agrees well with our resolved star-based profile, with the the different methods converging nicely at 10 kpc.

Zoom In Zoom Out Reset image size

We calculate the average minor axis g − i color profile using the same minor axis area and radial regions as used for the SB profile (Section 5.1.1). For detected sources in each radial bin, we convert the measured g − i color to Q_Col by rotating the CMD −22° around a point (1.62,24.8) 0.5 mag below the TRGB (see Section 4.1.2; Figure 6). We also compute Q_Col for all sources at >65 kpc, assuming these to be background contamination (following Figure 7, Section 5.1.1). For each radial bin, we scale these uniform background counts to the area of the bin. We then subtract the scaled cumulative Q_Col distribution of these background sources from the measured cumulative Q_Col distribution in each radial bin. The average Q_Col at each radius is then calculated as the median of the background-subtracted distribution, which is then rotated back to obtain the average g − i. Following the results of our ASTs (see Section 4.1.2), we add a constant 0.08 mag to the average color at each radius to account for the HSC observations' lack of sensitivity to red stellar populations observed in all GHOSTS fields.

We accounted for uncertainties on our average color measurements from two sources. We first estimate the upper and lower standard errors around the average Q_Col by separately calculating the 16%–50% and 50%–84% percentile spreads of the Q_Col distribution in each radial bin and dividing by the square root of the number of stars in each bin. Second, we estimate standard $\sqrt{N}/N$ Poisson errors at each radius, where N is the number of stars in each bin.

Figure 11 shows the minor axis color profile for M81. Our measurements are again shown in blue, and the GHOSTS points again in light gray for comparison (Monachesi et al. 2016a). While in relative agreement with the GHOSTS measurements, our Subaru profile appears to be ∼0.05 mag bluer than the GHOSTS profile, though with much less scatter. To ascertain the origin of this discrepancy, we measured the average color in each of the overlapping GHOSTS field regions. Though there is considerable scatter, due to the low number of RGB-like sources detected in regions as small as the ACS and WFC3 fields, we confirm that, on average, the median colors measured with Subaru in each GHOSTS field exhibit the same "blue bias" observed in the ASTs (Section 4.1). The average Subaru-measured colors in each GHOSTS field, corrected for the 0.08 mag "blue bias" between Subaru and GHOSTS, are shown in dark gray. These measurements agree much better with the original GHOSTS measurements, suggesting that the difference in the profiles is due to GHOSTS' narrow and more stochastic coverage of M81's minor axis, compared to the global view provided by Subaru.

Zoom In Zoom Out Reset image size

We recover the GHOSTS measurement of an approximately flat profile at R  ≳  20 kpc, g − i  ≃  1.68. However, we also observe a distinct negative color gradient for R  ≲  20 kpc, which is not due to the effects of crowding, which we observe to have no color dependence in our AST analysis (Section 4.1). This gradient smoothly connects the flat region of the profile to a single inner GHOSTS field (10 kpc), observed by Monachesi et al. (2013, 2016a), which is much redder than the outer fields. At first seemingly a "jump" in the profile, when combined with our Subaru observations, this inner field measurement appears to confirm that M81 possesses a steep minor axis color gradient within 20 kpc.

To estimate how this translates to metallicity, we use model HST–SDSS color–color tracks (Section 4.1.2) to convert our average g − i colors to metallicity, using the calibration of Streich et al. (2014). Though this conversion is somewhat uncertain, it is heartening that the outer portion (i.e., >20 kpc) of our halo profile matches the previous estimate of [M/H] = −1.2 (Durrell et al. 2010; Monachesi et al. 2013), which used deep HST data reaching the "Red Clump", almost exactly. With this metallicity calibration, we estimate that the ∼0.22 mag change in color from 10–20 kpc corresponds to a ∼0.51 dex change in [M/H], from  ∼ −1.2 to  ∼   −0.69. This yields a metallicity gradient of slope  ∼   −0.05 dex kpc⁻¹ inside 25 kpc—five times steeper than the global metallicity profile of M31, and comparable to M31's inner 25 kpc (Gilbert et al. 2014).

While this is the first observed case of such a distinct break in the color/metallicity profile of a MW-mass galaxy, galaxies with similar metallicity profiles to M81—i.e., displaying negative initial gradients, which flatten at large radii—have been observed in simulations (e.g., Monachesi et al. 2019). However, it is very rare to find even a simulated galaxy with such a sharp transition at <30 kpc. We discuss two possible origins of this steep color profile in Section 6.

In Figure 12, we present a global map of resolved RGB stars in M81's halo. Each star has been color-coded by its best-fit photometric metallicity, rather than g − i color, as metallicity is the more intuitive (while uncertain) quantity and is more directly comparable to other similar data sets. For this result, we estimate the metallicity for each individual star, using a grid of PARSEC isochrones (Bressan et al. 2012), age = 10 Gyr, ranging from [M/H] = −2 to 0 with steps of Δ[M/H] =0.05 dex. The distance in g − i color, at the given i magnitude, is evaluated for each star, for each isochrone. The best-fit metallicity is then defined as the model that minimizes the data−model g − i color residual. We display [M/H], rather than [Fe/H], so as to remain agnostic about [α/Fe]. Accounting for photometric uncertainties alone (not systematic uncertainties associated with age or different stellar evolution models), the typical [Fe/H] error is  ≤0.2 dex.

Zoom In Zoom Out Reset image size

The ongoing interaction between M81, M82, and NGC 3077 is immediately visible in the resolved star map. The NGC 3077 outskirts display an "S" shape, typical in tidally disrupting systems, while M82's debris is more compact. The tidal debris around both satellites is quite metal rich. The rest of the halo, however, is quite metal poor, comparable to M81's minor axis. Other than the interaction debris, five previously known satellite galaxies are visible (see Figure 3; IKN: Karachentsev et al. 2006; BK5N: Caldwell et al. 1998; KDG 61: Karachentseva & Karachentsev 1998; d0955+70: Chiboucas et al. 2009, 2013; d1005+68: Smercina et al. 2017), though there are no obvious substructures.

Figure 13 turns our map of resolved RGB stars into a map of stellar mass density in M81's halo. Using the method described in Section 4.1.1, we convert our RGB map to corrected RGB counts (Equation (2)). Again using a fiducial age = 10 Gyr, [Fe/H] −1.2 isochrone (following Harmsen et al. 2017), we convert RGB density to a corresponding stellar mass density, Σ_⋆ in M_⊙ kpc⁻², computed within ∼kpc × kpc pixels. We showed in Figure 10 that this method of SB/stellar mass estimation agrees well with integrated light measurements. The crowded centers of M81, M82, and NGC 3077 (see Figure 12) have been filled in with publicly available K_s-band images from the 2MASS Large Galaxy Atlas (Jarrett et al. 2003). We have clipped Σ_⋆ to >5.4 × 10³ M_⊙ kpc⁻²–equivalent to μ_V ≃ 33 mag arcsec ⁻², or one RGB star kpc⁻², the faintest limit we measure along the minor axis (e.g., Figure 9, Section 5.1.1). Combining our star count measurements with traditional near-infrared imaging, this map of stellar mass spans >4 orders of magnitude—from the dense stellar bulges at the centers of the primary galaxies, to the faintest stellar outskirts. This is among the most sensitive maps of stellar mass density ever constructed for a MW-mass galaxy.

Zoom In Zoom Out Reset image size

Figures 12 and 13 clearly indicate that there is a significant amount of metal-rich stellar material around M82 and NGC 3077. While M81's minor axis gives the properties of its past accretion history (i.e., M_Acc  = 1.16  × 10⁹ M_⊙;  see Section 5.1.1), any of the material around the two satellites that is unbound should be included in the current halo properties. To estimate how much of the material is unbound from M82 and NGC 3077, we estimate their respective tidal radii, using the basic approximation (von Hoerner 1957; King 1962)

$$\begin{eqnarray}&&{r}_{\mathrm{tid}}\simeq R{\left(\displaystyle \frac{{M}_{\star ,\mathrm{sat}}}{2{M}_{\mathrm{enc}}(R)}\right)}^{1/3},\end{eqnarray} \tag{ 5 }$$

where r_tid is the tidal radius, R is the separation between the central and the satellite adjusted for projection (i.e., $R\,=\sqrt{3}\,{R}_{\mathrm{proj}}$), M_⋆,sat is the stellar mass of the satellite, and M_enc(R) is the total mass of the central enclosed within R. To estimate M_enc(R), we adopt the familiar approximation for a flat rotation curve,

$$\begin{eqnarray}&&{M}_{{\rm{enc}}}(R)=\displaystyle \frac{{v}_{{\rm{c}}}^{2}\,R}{G},\end{eqnarray} \tag{ 6 }$$

where we have taken v_c = 230 km s⁻¹ from M81's H I rotation curve at 10 kpc (de Blok et al. 2008).

The projected separations from M81 of M82 and NGC 3077 are 39 kpc and 48 kpc, respectively, and their stellar masses are 2.8 × 10¹⁰ M_⊙ and 2.3 × 10⁹ M_⊙ (S4G; Sheth et al. 2010; Querejeta et al. 2015). Taking v_c = 230 km s⁻¹, this yields projected tidal radii of 10 kpc for M82 and 8.2 kpc for NGC 3077. Circles with radii equal to these tidal radii are shown in white on Figure 13. We then consider all material outside of these circles to be unbound. This amounts to 2.12 × 10⁸ M_⊙ around M82 and 3.36 × 10⁸ M_⊙ around NGC 3077—a total of ∼5.4 × 10⁸ M_⊙, which is a substantial fraction of M81's integral past accreted mass (∼10⁹ M_⊙). Taking a mass-weighted average metallicity of this material yields [Fe/H] ≃−0.9—significantly more metal rich than the rest of the halo.

Figure 14 combines Figures 12 and 13. The mass density map is divided into three average metallicity channels: [Fe/H] ∼ −0.8 (red), [Fe/H] ∼ −1.2 (green), and [Fe/H] ∼ −1.8 (blue). Each channel is then intensity weighted and combined into a three-channel color image. This figure highlights the visual impact that the massive and metal-rich debris around M82 and NGC 3077 has on the inferred mass and metallicity of M81's halo.

Zoom In Zoom Out Reset image size