Saying Hallo to M94's Stellar Halo: Investigating the Accretion History of the Largest Pseudobulge Host in the Local Universe

Resolved stellar populations are powerful tools for measuring halo properties, but at the distance of M94, old low- and intermediate-mass MS stars, containing the bulk of the stellar mass, are too faint to detect at the depths of HSC. Therefore, we rely on measuring RGB stars whose metallicities are a strong function of their color (Streich et al. 2014). They are also a more luminous population, hence we are able to resolve RGB stars within two magnitudes of the TRGB, which is 2–3 orders of magnitude brighter than their MS counterparts and still trace the underlying stellar population.

Unfortunately, at the depth achieved by HSC, resolved stellar sources are contaminated by galaxies because HSC can detect, but not resolve, distant galaxies (see Figure 2, top left), so it is important to have a systematic method to achieve star–galaxy separation. We first cut sources by FWHM with the criterion FWHM < 075 to select only point and compact sources. This is shown at the top right of Figure 2, which depicts a color–magnitude diagram (CMD) of all unresolved "point" sources in our sample. In shallower data sets, such cuts are sufficient owing both to the relatively lower number of galaxies and their typically larger FWHMs at brighter magnitudes. In our deep data set, compact and unresolved galaxies still outnumber stars, being particularly apparent at colors g − i ∼ 0.

Zoom In Zoom Out Reset image size

Stars occupy an empirical "stellar locus" in color–color diagrams that that can be distinct from galaxies. We use the g − r/r − i stellar locus (e.g., Covey et al. 2007; Ivezic et al. 2007; High et al. 2009; Davenport et al. 2014) to select stars in our survey and remove contaminants like galaxies and quasars. This is shown as the red line in Figure 2 at the bottom left. We chose sources that are within their photometric $\sqrt{{\sigma }_{g-r}^{2}+{\sigma }_{r-i}^{2}}$ errors of this line. This leaves us with our final list of stellar candidates, whose CMD is shown at the bottom right of Figure 2. The RGB sequence is clearer in this CMD, and unresolved galaxies with g − i ∼ 1 are much reduced in number (compare the top-right and bottom-right panels of Figure 2), but compact galaxies with g − i ∼ 0 are still numerous, primarily because it is at this point where distant background galaxies and stars have similar gri colors.

Using this final selection, we then narrow our focus to just the region around the RGB.

We wish to obtain a more robust color distribution that will better reflect the metallicity of the population, so we calculate the "Q-color" of candidate RGB stars following Monachesi et al. (2016a). The intention of this metric is to give the color of an RGB star, referenced to a magnitude 0.5 mag below the TRGB. We achieve this by rotating the (color, TRGB magnitude) point of (2, 25.036) around a −22° angle, making the metal-poor RGB vertical.

From the distribution of stellar sources between 12 and 30 kpc from the center of the galaxy, we see that the RGB branch on the CMD looks to have two different regions within it, one centered around a Q-color of ∼1.5 and the other centered around ∼2. We also fit a four-component Gaussian mixture model to the peaks of the 1D Q-color distribution to obtain the peak colors of the two regions. A Q-color CMD and a 1D histogram of the Q-colors are shown in Figure 3. Each fitted Gaussian peak is overlaid in yellow, with the sum of all the peaks shown in green.

Zoom In Zoom Out Reset image size

A four-component Gaussian mixture model describes the features of the Q-color distribution reasonably well. The RGB region of M94 has two clear and distinct populations corresponding to RGB stars at two different characteristic colors of Q = (g − i)_TRGB+0.5 = 1.64 and 2.22. In addition, a blue peak accounts for the population of background galaxies; this component has little impact on this study as we select against it in color–magnitude space. The fourth broad peak models (likely imperfectly) the foreground star and compact background population, which has a wide range in color, extending well beyond the colors expected of RGB stars. It is clear from this consideration that fore- and background objects do contaminate our RGB selections, requiring our later use of background subtraction for quantitative investigation of the density and color profiles of RGB stars.

Bearing this in mind, we focus on the two RGB stellar populations for the remainder of this paper. These regions are displayed on top of the full stellar source CMD in Figure 4, with best-fit isochrones overlaid. We will refer to the blue-most region as bRGB and the redder one as rRGB throughout. The best-fit isochrones for each region have an [M/H] = −1.4 for bRGB and [M/H] = −0.4 for rRGB assuming an age of 10 Gyr. The metallicity for bRGB is derived from generating a set of PARSEC (Bressan et al. 2012) isochrones with metallicities between 0.0001 and 0.02 in 0.00005 increments, finding their TRGB + 0.5 color, and comparing to the peak Q-color of the region that has been transformed to a g − i color. The metallicity for rRGB stars is derived star-by-star because the metal-richer isochrones are more tilted and curved than the metal-poorer ones, and so a conversion to Q-color becomes mildly magnitude dependent. For each star, we find the g − i color difference at each i-band magnitude for each star in the RGB region. The best-fit metallicity is then taken from the model that minimizes the difference between the g − i color of the RGB stars and the isochrone model. Thereby, the median metallicity of the rRGB stars is taken from the median of the individual metallicity measurements of all the rRGB stars. ¹⁰

Zoom In Zoom Out Reset image size

We note that it would be possible to fit the RGB colors with different combinations of age and metallicity. In particular, younger, slightly more metal-rich isochrones can also lead to a good fit (we estimate an [M/H] ∼ −1.07/−0.17 and −1.25/−0.3 for the bRGB/rRGB of a 4 and 7 Gyr population, respectively). This degeneracy does not affect the main conclusions of this paper, with the degree of systematic error introduced by the choice of stellar population being within the final quoted error bars. Nonetheless, we argue that our choice of a fiducial age for the halo stellar population of 10 Gyr is appropriate. Indeed, magnetohydrodynamics simulations such as those used in Monachesi et al. (2019) find that the median age for accreted stellar populations in the stellar halo of Milky Way–mass galaxies is 9.4 Gyr at a galactocentric distance of 30 kpc. Other studies such as by Font et al. (2006) and McCarthy et al. (2012) also find a lower limit on the age of accreted halo stars to be ≳10 Gyr.

In Figure 5 we show the spatial distribution of stars in bRGB and rRGB color coded by their best-fit metallicities. Superimposed on a roughly uniform background of fore- and background contaminants, we see clear concentrations of stars at M94's distance. This visualization shows the striking difference between the two regions in both spatial distribution and metallicity: bRGB is clearly composed of a more metal-poor population, and its stars are more widely distributed than those in the rRGB population, whose stars are clustered in a dense circular region at small radii until ∼30 kpc. We later calculate the brightness profile of each population, showing that the rRGB population follows an exponential profile that directly follows the outer disk profile from Watkins et al. (2016); we therefore attribute the rRGB population to M94's disk in what follows. We attribute the bRGB population, which is metal poorer, more diffuse, appears mildly asymmetric, and shows hints of possible streams, to M94's stellar halo.

Zoom In Zoom Out Reset image size

In order to calculate the color and surface brightness profiles accurately, it is important to determine the radius at which we are dominated by the roughly uniform background of contaminants apparent in Figures 3 and 5. We examine Hess CMDs covering a series of radial bins from the center of the galaxy, as shown in Figure 6.

Zoom In Zoom Out Reset image size

From the visual inspection, we conclude that the bRGB region starts to be dominated by background galaxies after ∼50 kpc, and the rRGB region after ∼30 kpc. After these thresholds, the density counts drop drastically in each respective region and we cannot see a clear RGB feature. These background thresholds are used for the rest of the analysis.

Before we can measure these profiles however, we need to account for completeness and crowding in our sample. As illustrated in Figure 5, there is a gap in the very central part of the galaxy due to crowding effects such that the HSC reduction pipeline was unable to extract sources reliably. There are also sources that overlap in high-density regions, especially around the first few radial bins away from the galaxy's center. The HSC source detection pipeline is not able to locate and categorize every source with full accuracy, especially those that are blended together, making it necessary for us to account for this source of error.

In order to quantify the completeness of our sample rigorously, we ran a suite of artificial star tests (ASTs). In order to avoid the injected artificial stars from causing crowding, we inject a total of 674,980 artificial stars split into five runs of 134,996 stars each. Artificial stars were drawn from uniform grids in position, in i-band magnitude between i = 22.5 and i = 28, and color between g − i = −2 and g − i = 5. The r-band magnitudes were then calculated assuming each artificial star lies on the g − r/r − i stellar locus. At the input artificial source density of ∼10 arcmin⁻², these artificial sources do not measurably effect the measurement of the source properties by the pipeline (e.g., Smercina et al. 2020). The pipeline was then run on each of these five artificial star runs in the exact setup used to process the original observations, resulting in a list of recovered artificial stars. For our purposes, an input artificial star was declared recovered if its recovered i-band magnitude was within 0.5 mag and its recovered position within 03 of its input magnitude and position, respectively.

Our particular interest is in calculating the completeness in each of our two RGB regions (bRGB and rRGB), in the same regions of the image used to create the star count profiles (e.g., Figure 7). We choose a set of radial bins, 5 kpc wide from 10 to 15 kpc and 2 kpc wide thereafter. Artificial stars in the color–magnitude regions of bRGB and rRGB were analyzed for each spatial region, and the fraction of those recovered was weighted by the ratio of the number of input stars to the expected luminosity function of RGB stars with their appropriate best-fit isochrones (from Figure 4), as a function of i-band magnitude. In this way, for each spatial region of interest, a single completeness for each bRGB and rRGB star was determined with which to scale the observed number of bRGB and rRGB stars, respectively, in each region.

Zoom In Zoom Out Reset image size

For studying resolved stellar populations in the context of mergers, typically one would like to measure only the stars along the minor axis of a galaxy since these are relatively free of contamination from in situ stars (see Monachesi et al. 2016b; Smercina et al. 2020). Since M94 is not as highly inclined as galaxies such as M81, its minor axis traces stars in the plane of the disk, so we once again used the same circular radial regions as referenced in Section 3.3 for our analysis. We choose to analyze radial profiles out to a 70 kpc radius. For the bRGB (halo) population we adopt a background region from a 50 to 70 kpc projected radius, yielding a background density of 0.21 arcmin⁻². The rRGB (disk) is detectable only at a much smaller radius, and a background density of 0.13 arcmin⁻², determined between 30 and 70 kpc, was adopted.

For each radial bin, we calculate the RGB-selected and background-subtracted source density, dividing by the completeness fraction found from our AST analysis. The conversion between the number of RGB stars in the observed magnitude range to stellar mass assumes an [M/H] = −1.4 (−0.4) isochrone for the bRGB (rRGB) population and a 44% mass loss between the initial and present day masses, which is compatible with a Chabrier initial mass function. When calculating the surface brightnesses, we adopt a V-band mass to light ratio of 1.25, appropriate for ∼10 Gyr-old populations with the range of relevant metallicities from Bruzual & Charlot (2003)'s single-burst models.

Formal uncertainties from counting statistics are negligible. Instead, we estimate uncertainties in our density profiles that incorporate larger-scale variations (owing to halo substructure or variations in the background) by assessing the variation of our profiles as a function of azimuth. We divide the data set into radial slices of 45° width and repeat the same background estimate and mass density calculations as before for both RGB regions. For each radial bin, we perform a bootstrap resampling analysis and chose one of the quadrant mass densities at that radius, with replacement. We repeat this 1000 times and find the 16th and 84th percentile values of the resulting distribution, which we use, respectively, as the upper and lower error bars in Figure 7.

We also overplot the mass density from Watkins et al. (2016), who measured M94's inner disk. This paper assumes elliptical annuli around M94's center, while we used circular annuli, so we adjust their given radii for this and convert it to an equivalent circular radii, assuming an axis ratio of b/a = 0.79 (Watkins et al. 2016). This is shown in Figure 7. The rRGB profile agrees closely with the profile of the outer disk from Watkins et al. (2016) within 20 kpc, the range where the measurements overlap. We can also see that bRGB's profile is much flatter than rRGB's. These stark differences leads us to affirm that bRGB and rRGB represent two different stellar populations in the galaxy. rRGB shows an exponential drop in its profile and agrees with Watkins et al. (2016)'s inner disk profile, strongly suggesting that its stars are actually in the disk of M94, while bRGB represents the halo population.

To measure the total accreted stellar mass, we fit a three-parameter power-law model with intrinsic scatter following Hogg et al. (2010) from 16 to 42 kpc (Figure 8). We use a power law of the form ${\rm{\Sigma }}=a{\mathrm{log}}_{10}(\tfrac{x}{{R}_{0}})+{{\rm{\Sigma }}}_{0}$ with a likelihood of

$$\begin{eqnarray*}&&\mathrm{log}{ \mathcal L }=-0.5\displaystyle \sum _{i=0}^{n}\displaystyle \frac{{({{\rm{\Sigma }}}_{i}-{\rm{\Sigma }}(x,{{\rm{\Sigma }}}_{0}))}^{2}}{{\sigma }^{2}}+{log}({\sigma }^{2}),\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{\sigma }^{2}={\sigma }_{{\rm{\Sigma }}}^{2}+V,\end{eqnarray*}$$

with σ_Σ being the error in the mass density in each radial bin and V being the variance of the intrinsic scatter. R₀ = 37.9 kpc is a normalization factor that we take to be the mean of the radial bins used for fitting. We find best-fit parameters of $a=-{2.61}_{-0.17}^{+0.16}$ and $\mathrm{log}{{\rm{\Sigma }}}_{0}={3.68}_{-0.03}^{+0.03}{M}_{\odot }\,{\mathrm{kpc}}^{-2}$. Following Harmsen et al. (2017), we integrate this profile from 10 to 40 kpc and calculate an accreted stellar mass of M_*,10–40 = 9.25 × 10⁷ M_⊙. Harmsen et al. (2017) uses simulated stellar halo data from Bullock & Johnston (2005) to estimate that the total accreted halo mass, M_acc, is about three times that of the mass between 10 and 40 kpc, M_*,10–40, with variations of ∼40% from halo to halo. Accounting for this uncertainty in our error calculation, we extrapolate to the total accreted mass following this method and get a total accreted stellar mass of ${M}_{* }={2.77}_{-1.00}^{+1.54}\times {10}^{8}{M}_{\odot }$.

Zoom In Zoom Out Reset image size

We calculate a median g − i color and metallicity profile for the bRGB population using the same bins as the mass density profile. We first calculate a 2D histogram of all the sources and then estimate the area and sources in the background region (everything between 50 and 70 kpc) to get the background source color distribution. For each radial bin, we find the cumulative distribution of the background Q-color weighted by the ratio of the bin area to background area, which we use to get the total number of background sources that correspond to a certain color. We then find the cumulative distribution of bRGB stars in each bin and subtract the cumulative background distribution normalized to the area of the radial bin and find the median of this background-subtracted distribution, converting back to g − i color. The resulting color profile is shown in Figure 9. We estimate a metallicity gradient of slope ∼ −0.005 dex kpc⁻¹ between 16 and 40 kpc.

Zoom In Zoom Out Reset image size

The color gradient is d(g − i)/dr = −0.0028 mag/kpc. Assuming 10 Gyr-old stellar populations, and recalling that RGB color is a very weak function of age, the g − i color should correspond to the metallicity of the population, shown schematically on the right-hand side of Figure 9. While the mapping to metallicity is nonlinear, in the color range relevant to M94's bRGB population a metallicity gradient of d[M/H]/dr ∼ −0.005 dex kpc⁻¹ is implied, if the g − i color gradient is interpreted as being due to metallicity alone. To compare this with other halos requires translating to the HST passbands used in Monachesi et al. (2016a) and Harmsen et al. (2017). A linear fit to g − i versus F606W – F814W for stars with HST measurements in the rough color range around g − i ∼ 1 gives a fit of g − i ∼ 1.91 (F606W − F814W) − 0.64; correspondingly, the implied color gradient in F606W – F814W is d(F606W − F814W)/dr = −0.0015 mag/kpc. This color gradient is rather mild, similar to galaxies like the Milky Way, NGC 253, or NGC 4945.

We follow numerous past works and use M94's stellar halo properties to make inferences about its merger history (e.g., Bell et al. 2017). We show M94's halo in context of the observed stellar halo metallicity–mass relation for 14 nearby Milky Way–stellar mass galaxies in Figure 10, which was first discussed in Harmsen et al. (2017). Since RGB color is affected primarily by [M/H] and not [Fe/H] (Streich et al. 2014), we use the α-agnostic [M/H] instead of converting our measurement to [Fe/H]. M94 follows the halo metallicity–mass correlation well and contains one of the least-massive and most metal-poor halos compared to these other disk galaxies, with only M101 harboring a smaller, more-anemic halo. This generally agrees with the general expectations from theoretical accretion-only models wherein the halo is accreted from the merger of its last satellite progenitor, and that a smaller halo is formed through the aggregation of a few less-massive satellites instead of being dominated by a single very-massive merger (Monachesi et al. 2019). This also points to M94 likely having a more quiescent and quieter merger history compared to the majority of local Milky Way–stellar mass galaxies (Deason et al. 2015a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

D'Souza & Bell (2018) found a correlation between the mass of a galaxy's dominant merger, M_dom, and the halo's total accreted stellar mass, M_acc, (see also Zhu et al. 2022). This assumes that only a fraction, frac_Dom, of the accreted mass contributes to the stellar halo by the dominant merging partner. For the majority of galaxies that fall a slight ways below the 1-to-1 line in the M_dom−M_acc plane, frac_Dom is less than 1. In most cases, the accreted halo mass M_acc is greater than the mass of the dominant progenitor M_dom. This is especially true at the lower-mass end of M_acc and indicates that the halo contains contributions from a number of other merger events, consistent with our idea of the hierarchical buildup of galaxies. In order to estimate the mass of M94's dominant merger partner, we follow D'Souza & Bell (2018) and correlate the dominant merger mass with the total accreted stellar mass for galaxies taken from the Illustris simulation suite (Vogelsberger et al. 2014a, 2014b; Genel et al. 2014), chosen to have DM halo masses of $11.7\leqslant \mathrm{log}{M}_{\mathrm{DM}}\leqslant 12.15$. We fit this distribution with a line, and infer the mass of M94's dominant merger from our accreted mass calculated in Section 4.1, as shown in Figure 11. The uncertainty in the dominant mass estimate is taken as the sum in quadrature of the intrinsic scatter in the distribution and the uncertainty of the measured accreted mass. The resulting dominant stellar mass is estimated to be $\mathrm{log}{M}_{\mathrm{dom}}={8.15}_{-0.22}^{+0.29}{M}_{\odot }$. The Small Magellanic Cloud's (SMC) stellar mass is 4.6 × 10⁸ M_⊙ (McConnachie 2012); M94's most-massive merger partner is likely to be less massive than even the SMC.

Thus we see that M94 is characterized by an inactive merger history which gives rise to its small, metal-poor halo. Two other nearby galaxies, M83 (M_halo = 4.2 × 10⁸ M_⊙; M. Cosby et al. 2023, in preparation) and M101 (M_halo = 8.2 × 10⁷ M_⊙; Jang et al. 2020), exhibit this same type of quiet merger history. Yet all three of these galaxies have very different bulge masses. While M94's pseudobulge is massive (almost 10^10.5 M_⊙; Trujillo et al. 2009), there is no evidence for a classical bulge in either M101 (Kormendy et al. 2010) or M83 (de Jong et al. 2008; Fisher & Drory 2010); both galaxies host small (≤10^9.5 M_⊙) pseudobulges. These three galaxies—with outwardly very similar merger histories—have such wildly dissimilar bulge properties, suggesting that the mass of a galaxy's dominant merger is not simply correlated with the mass or properties of its bulge.

We connect M94 to the broader population of galaxies with measured halos and bulges in Figure 12, with galaxies marked as having either a classical bulge or a pseudobulge. If we take the halo stellar mass to act as a proxy for the mass of the galaxy's dominant merger, then we see that almost every other pseudobulge-hosting galaxy in our local universe underwent a larger merger than M94, and yet hosts a less-massive bulge. This adds evidence to the picture that mergers are not the primary drivers of pseudobulge creation, and any role that they play in determining its properties would be complicated and intertwined with other, secular processes that are ongoing within the galaxy (Bell et al. 2017; Gargiulo et al. 2019, 2022). Indeed, gas flows in so-called "oval distortions" like the kind M94 has can be a major impetus for the formation of pseudobulges (Speights et al. 2019) and can also fuel star-forming rings and enhance star formation in general (Athanassoula & Palous 1984; Mulder & Combes 1996; Li et al. 2016). M94 currently has about equal fractions of intermediate-age and old-age stars in its pseudobulge, found by taking spectra within its central 1'' × 11 region (Fernandes et al. 2004; Mason et al. 2015), indicating semirecent star formation, potentially from gas driven inward from a Lindblad resonance (Wong & Blitz 2001; Trujillo et al. 2009).

Zoom In Zoom Out Reset image size

What makes M94 even more interesting are features such as its warped disk, with its H i disk misaligned with its optical counterpart and differential alignment of its inner disk (oval distortion) and outer disk, and its highly noncircular motions (de Blok et al. 2008; Walter et al. 2008; Trujillo et al. 2009). Trujillo et al. (2009) rules out merger as an explanation for M94's differential inclination, citing merger simulations that give too low of a resulting inclination to match the observations. Instead, Trujillo et al. (2009) note that the inclination discrepancy in this case goes away if the "inner disk" is an oval distortion, so it has a completely different ellipticity than the outer disk. And while merger is a common explanation for irregularities in galactic morphology, given our inferences regarding M94's quiet merger history, the presence of these features indicate that the gas morphology of a galaxy can be strongly distorted by processes that involve no merger with an established galaxy. M94 would not be the only curious example—galaxies such as M83 have been found to also harbor a warped and disturbed disk (Barnes et al. 2014; Heald et al. 2016), yet have no sign of a recent large merger event (M. Cosby et al. 2023, in preparation).

The measure of M94's accreted and dominant merger masses has also been used to gain insight into the satellite populations of Milky Way–stellar mass galaxies. Satellite luminosity functions have now been measured for a substantial number of Local Volume galaxies, finding a large range in the number of satellites in each group (Smercina et al. 2018; Müller et al. 2019; Carlsten et al. 2020; see Carlsten et al. 2022 for an overview), where M94 is a particularly extreme example, having only two classical satellites (Smercina et al. 2018). The number of satellites appears to vary with the host galaxy's stellar mass (Carlsten et al. 2020) as would be expected with a scaling with overall DM halo mass. Yet, considerable scatter remains, and has been claimed to correlate with several different measures that relate to a group's merger history: "tidal index," a measure of environmental density (D. Sand, private communication), magnitude gap (Kim et al. 2022), and the mass of the galaxy's most-massive past merger (Smercina et al. 2022). While simulations do not yet show such behavior (Smercina et al. 2022), such correlations may broadly reflect the delivery of satellites in large merger events (Deason et al. 2015b; D'Souza & Bell 2021). Owing to M94's low number of satellites, it has particular constraining power in such correlations, and M94's combination of a low halo mass (and thus low-mass dominant merger) and low satellite number are consistent with the view that satellite number correlates with merger mass (Smercina et al. 2022).