Relating the Diverse Merger Histories and Satellite Populations of Nearby Galaxies

A hallmark of the current Λ cold dark matter (ΛCDM) paradigm is that the halos of dark matter in which galaxies reside assemble hierarchically (e.g., White & Rees 1978). As a consequence of this hierarchical growth, dark matter and galaxy formation models alike predict that galaxies should accrue extensive populations of lower-mass dwarf "satellite" galaxies, tidally interacting with many of them (e.g., Somerville & Kolatt 1999; Bullock et al. 2001; van den Bosch 2002; Wechsler et al. 2002). The largest of these accreted galaxies experience strong dynamical friction and can merge with the central galaxy before they tidally disrupt, strongly impacting the properties and evolution of the central galaxy (e.g., Barnes & Hernquist 1991; Hopkins et al. 2006). Modern cosmological simulations, built on a ΛCDM framework and incorporating these processes, are now able to reproduce many of the most well-constrained large-scale observables, such as the bright end of the galaxy luminosity function (e.g., Bell et al. 2003) and the cosmic star formation rate (SFR; Madau et al. 1998). Perhaps the most important benchmarks of these simulations are the canonical galaxy scaling relations—for example, the Kennicutt-Schmidt star formation relation (e.g., Kennicutt 1998; Orr et al. 2018), the Tully Fisher relation (e.g., Tully & Fisher 1977; Dalcanton et al. 1997; Bell & de Jong 2001; Schaye et al. 2015), and the galaxy star-forming main sequence (Noeske et al. 2007; Schaye et al. 2015; Jarrett et al. 2017)—which operate across the large range of scales that govern galaxies' properties. Yet, despite this impressive progress in explaining the formation and evolution of relatively massive galaxies, the regime of small-scale galaxy formation has remained problematic for simulations.

This difficulty in modeling small-scale galaxy formation is due, in large part, to the sensitivity of low-mass galaxy properties to stellar feedback and other baryonic processes. Consequently, this results in a magnification of small differences in physical prescriptions on predicted galaxy properties in the low-mass regime (e.g., Bullock & Boylan-Kolchin 2017). As processes such as stellar feedback can operate on both extremely small spatial and short timescales, simultaneously achieving the required spatial and temporal resolution for the complex associated physics is extremely difficult (e.g., Hopkins et al. 2014). Thus, modeling galaxy formation at the smallest scales, and therefore investigating the model-based links between small- and large-scale galaxy formation, has proved challenging to date.

The persistent difficulty of modeling small-scale galaxy formation has motivated two decades of careful theoretical and observational study of dwarf galaxies. Difficult to detect observationally owing to their intrinsic faintness, dwarf galaxies are both the lowest-mass and most dark-matter-dominated galaxies in the universe (e.g., see the recent review by Wechsler & Tinker 2018). Consequently, their properties are critical benchmarks for our models of galaxy formation. Dwarf galaxies are the observational bedrock of our understanding of small-scale cosmology (e.g., see the recent review by Bullock & Boylan-Kolchin 2017).

While they are predicted to also exist in isolation, until recently nearly all of the low-mass dwarf galaxies readily accessible for observational study were "satellites" of Milky Way (MW)-mass galaxies (e.g., see McConnachie 2012, for an overview of the Local Group (LG) satellite populations)—existing within an environment that is dominated by a much larger central galaxy. This is an important consideration, as the large-scale environments of galaxies like the MW are not static, "closed" boxes—they are complex ecosystems that frequently experience substantial mergers with other systems throughout their lives.

Are galaxies like the MW expected to have experienced similar ensemble merger histories over the course of their lives? Quite the opposite! Theoretical modeling has long predicted that MW-mass galaxies should have a wide range in their merger histories, and particularly in the largest mergers that they have experienced (e.g., Purcell et al. 2010; Deason et al. 2015a), with many galaxies experiencing mergers with progenitors more massive than the Large Magellanic Cloud (LMC; Stewart et al. 2008). Study of their stellar halos has largely confirmed this predicted diversity by observing a wide range of stellar halo masses and metallicities around galaxies with masses comparable to the MW, with correlations between them showing that stars accreted from the largest past merger dominate the properties of the halo stellar populations (Bell et al. 2017; Harmsen et al. 2017; D'Souza & Bell 2018a; Smercina et al. 2020). Meanwhile, recent evidence suggests that the "classical" satellite populations of nearby MW-analogs are equally diverse (Geha et al. 2017; Smercina et al. 2018; Bennet et al. 2019; Carlsten et al. 2021).

Recent evidence from the LG suggests that the diversity in both of these fundamental properties of galactic systems may be related. Studies of the MW's satellites using Gaia data and detailed orbital modeling postulate that a number of its satellites may have been brought in during the infall of the LMC system (e.g., Gaia Collaboration et al. 2018; Kallivayalil et al. 2018; Jahn et al. 2019; Erkal & Belokurov 2020; Pardy et al. 2020; Patel et al. 2020)—possibly including some of the MW's "classical" satellites, such as Carina and Fornax. Additionally, deep resolved star formation histories (SFHs) of the M31's dwarf satellites with the Hubble Space Telescope (HST) have recently revealed that nearly 50% share a remarkably similar timescale for the "shutdown" of their star formation, ∼3–6 Gyr ago (Weisz et al. 2019). This is similar to the predicted first infall time of the massive progenitor galaxy whose merger with M31 likely formed M31's massive and structured stellar halo (M32p; D'Souza & Bell 2018b; Hammer et al. 2018). This suggests that the merger histories of the MW and M31 may need to be accounted for in order to explain the properties of their present-day satellite populations.

This represents an important new lens through which to view and test predictions for MW-like systems from galaxy formation models. Is the infall of satellites during mergers with galaxies such as the LMC or M32p surprising? In a dark matter only (DMO) view, the overall number of surviving subhalos within a central halo will primarily be a function of the mass of the central halo (e.g., Kravtsov et al. 2004), similar to what has long been seen in the idealized case of rich galaxy clusters (e.g., Bahcall & Cen 1993). However, due to the self-similarity of CDM, dwarf galaxies are predicted to host their own populations of satellites prior to group infall (e.g., Deason et al. 2015b; Dooley et al. 2017a, 2017b), and there has been direct confirmation of these "satellites-of-satellites" in both the Local and nearby groups (e.g., Deason et al. 2014, 2015a; Smercina et al. 2017). Therefore, these large satellites are entirely expected to bring in their own populations of low-mass satellites. Indeed, it is expected that the distribution of satellite infall times in MW-like halos should cluster around massive accretion events (e.g., Deason et al. 2015b; D'Souza & Bell 2021). Yet, because of the strong halo mass dependence on the number of preexisting low-mass satellites, and the fact that all galaxies experience accretions, the default DMO prediction is that the satellite populations of MW-mass galaxies should not, overall, correlate strongly with their late-time accretion history (Bose et al. 2020).

However, the baryonic physics involving infalling satellites' interactions with circumgalactic gas, tidal disruption by the central galaxy, and internal stellar feedback complicates matters significantly (e.g., Wetzel et al. 2015). Understanding the importance of merger history in predicting the properties of present-day satellite populations requires both extremely high resolution and a sample of MW-mass systems with diverse merger histories—a set of constraints that are only now being achieved in modern simulations (e.g., Simpson et al. 2018; Garrison-Kimmel et al. 2019; Monachesi et al. 2019; Samuel et al. 2020; Engler et al. 2021). Given this theoretical progress, it is therefore important to observationally explore how merger history relates to satellite populations, and compare it to our theoretical predictions.

In this paper, we leverage the recent wealth of observational insight into the merger histories and satellite populations of MW-mass galaxies to empirically investigate a connection between them in nearby systems, with the goal of providing an observational benchmark for galaxy formation models. The paper is structured as follows: in Section 2 we first outline the existing observational data for both satellite populations (Section 2.1) and stellar halo measurements (Section 2.2) around MW-like galaxies, and the intuition that has been gained from each. In Section 3 we present the results of our initial comparison between the most dominant merger and the number of classical satellites in eight nearby systems, which exhibit a clear relationship. In Section 4 we quantify the form and significance of this relationship (Section 4.1), accounting for an intrinsic scaling between satellite populations and galaxy mass (Section 4.2). In Section 5, we demonstrate the evolution of the MW's and M81's satellite populations along this relation. In Section 6 we compare the observed scaling relation to theoretical predictions from the Feedback In Realistic Environments (FIRE) simulations, followed by an in-depth discussion of the implications for galaxy formation and possible origins of the relation in Section 7. We summarize our conclusions in Section 8.

Addressing a possible connection between the satellite populations and merger histories of MW-mass galaxies requires measurement of two notoriously difficult-to-measure properties: (1) global satellite populations, and (2) global stellar halo properties. Efforts to study both of these properties for MW-mass systems have remained a substantial focus of the field over the last two decades. However, they have been historically kept in separate "intellectual boxes." In this paper, we combine these insights to ask the question: are galaxies' merger histories and satellite populations connected? As the vast majority of existing data has been obtained in MW-mass systems, this is where we concentrate our efforts. In this section we summarize the current state of the field regarding the satellites (Section 2.1) and stellar halos (Section 2.2) of nearby MW-mass galaxies.

2.1. The Diverse Satellite Population of MW-mass Galaxies

Tensions in the number and properties of observed MW satellites, relative to model predictions—e.g., the "Missing Satellites" and "Too Big to Fail" problems (Klypin et al. 1999; Moore et al. 1999; Boylan-Kolchin et al. 2011)—constitute some of the largest hurdles for the ΛCDM paradigm (e.g., Bullock & Boylan-Kolchin 2017). Solutions to these problems—often focusing on the impact of baryonic processes such as reionization (e.g., Bullock et al. 2001) and stellar feedback (e.g., Brooks et al. 2013)—have used the satellite population of the MW as a benchmark. Just this one system has been used to direct the scope of some of the most important problems in galaxy formation.

Motivated by this lack of context, the field's focus has shifted to surveying the satellite populations of nearby MW-analogs (i.e., central galaxies in the Local Volume with stellar mass M_⋆ ∼ 3–20 × 10¹⁰ M_⊙). Through a combination of wide-field integrated light searches, targeted resolved star studies, and spectroscopic surveys, these efforts have greatly enhanced our understanding of satellite galaxy populations. In all, eight MW-mass systems have now been surveyed to the depth of the MW's "classical" satellite population (M_V ≲ − 9): the MW and M31 (compiled by McConnachie 2012), M81 (compiled by Karachentsev & Kudrya 2014), M101 (Danieli et al. 2017; Bennet et al. 2019; Carlsten et al. 2019), Centaurus A (Crnojević et al. 2019), M94 (Smercina et al. 2018), M83 (e.g., Müller et al. 2015, 2017; Carrillo et al. 2017), and M104 (Carlsten et al. 2021).

With this newfound access to a true sample of satellite populations in MW-like galactic systems has come the realization that these populations are incredibly diverse (Geha et al. 2017; Smercina et al. 2018; Bennet et al. 2019; Carlsten et al. 2021). This diversity, particularly the discovery of the sparse satellite population around the "lonely giant" M94, is a particularly powerful constraint on small-scale galaxy formation, including the low-mass end of the stellar mass–halo mass (SMHM) relation (Smercina et al. 2018; Carlsten et al. 2021). Figure 1 (left) shows the satellite V-band luminosity functions (LFs) within 150 kpc projected galactic radius for the eight MW-mass galaxies that have been studied down to M_V ≲ − 9. We note that most of the satellites in these eight systems now have existing distance estimates (most from HST color–magnitude diagrams) and therefore their line-of-sight (LOS) distances can largely be constrained to within ±1 Mpc of the central galaxies, except in the cases of M101 and M104 (described below). Given this, we apply a 1 Mpc constraint, in addition to the 150 kpc projected radius, following Smercina et al. (2018). This is important for the MW and M31 where we take the projected satellite LFs, calculated from random LOSs within each. Another notable consequence is that NGC 5253 is excluded as an M83 satellite, due to its likely foreground status. ⁵

Zoom In Zoom Out Reset image size

There has been significant back-and-forth discussion in the literature about M101's satellite population, regarding which satellite candidates are truly associated with M101 versus in the background, as well as their integrated brightnesses (Müller et al. 2017; Danieli et al. 2017; Bennet et al. 2019; Carlsten et al. 2021). Bennet et al. (2019) determined that a number of the satellites considered by Danieli et al. (2017) were actually background objects, while Carlsten et al. (2021) added three satellites (eight versus five) compared to Bennet et al. (2019), due to a claimed discrepancy in their V-band luminosities. For this paper, we adopt Bennet et al. (2019) as M101's fiducial LF, given its similarity to the other surveys around nearby systems, but incorporate the Carlsten et al. (2021) LF into our uncertainties. The Carlsten et al. (2021) LF for M104 is also relatively uncertain, as surface brightness fluctuation (SBF) distance estimates are inconclusive for a number of its satellites. Given the richness of its satellite population, we consider it to be "measured," despite a range of faint-end LFs. The M104 LF shown in Figure 1 is an average of the lower- and upper-bound LFs in each magnitude bin.

2.2. Inferring a Galaxy's Most Dominant Merger from Its Stellar Halo Properties

Through these efforts, a sample of eight systems —the MW, M31, M81, M83, M94, M101, M104, and Cen A—now exists for which both complete satellite LFs and accreted masses have been measured. These data sets present a brand new and unprecedented opportunity to explore a connection between these two fundamental properties of nearby well-studied galactic systems.

First, we will define several quantities that will be used in this analysis. As discussed in Section 2.2, the most dominant merger a galaxy has experienced is well estimated by its total accreted stellar mass. However, in a number of systems, the most dominant merger is still in progress, in the form of a nearby massive satellite. In this case, the galaxy "ecosystems" (i.e., dark matter halos) have already begun to merge, but the galaxies have yet to coalesce. The most familiar example of this is the MW itself: as described in Section 1, a number of MW satellites may have originally been satellites of the (not-yet-disrupted) LMC and Small Magellanic Cloud (SMC). Additionally, recent evidence suggests that the LMC's massive halo has already had a substantial global impact on the MW's halo (Conroy et al. 2021), suggesting that they have indeed begun to merge on a large scale.

We contend that the most massive merged ecosystem is likely the important property to consider for satellite populations, and while the accreted mass is useful as a proxy for identifying past such mergers, the most massive existing satellite is an equally powerful proxy of current ongoing mergers, even if the galaxies have not yet interacted significantly. For example, Smercina et al. (2020) showed that in systems such as M81, which is in the early stages of a close merger with its massive satellite M82, the merger is not advanced enough to have redistributed the accreted material into the stellar halo (the estimate of the current accreted mass in M81's stellar halo is ∼1/20 of its estimated total future accreted mass). Yet given M81's modest accreted mass, M82 undoubtedly represents the most massive ecosystem merger that M81 has experienced. In light of this, we instead adopt a revised metric for the mass of the most dominant ecosystem merger, M_⋆,Dom, where the proxy is the larger of the total accreted mass or the mass of the most massive satellite. We adopt uniform conservative 0.3 dex uncertainties on M_⋆,Dom (following Bell et al. 2017; Harmsen et al. 2017; D'Souza & Bell 2018a), reflecting the large uncertainties on both the inferred stellar mass of nearby galaxies and model-dependent accreted mass estimates.

For satellites, we only have access to the "full" classical satellite populations (i.e., M_V < − 9 satellites out to the virial radius) of the MW and M31. We therefore follow recent studies (e.g., Smercina et al. 2018; Bennet et al. 2019; Carlsten et al. 2021) and restrict our attention to satellites within a projected galactic radius of 150 kpc, down to an absolute V-band magnitude of M_V < − 9 (see Figure 1). Though there is diversity in both survey area and depth, all eight systems are considered complete within these radius and brightness criteria.

Figure 2 shows M_⋆,Dom plotted against N_Sat for the eight available MW-mass systems, color-coded by whether M_⋆,Dom represents a past merger (i.e., accreted mass; blue) or an existing satellite (yellow). ⁷ A clear and surprisingly tight relationship is visible, spanning the entire three orders of magnitude in dominant merger mass: systems that have experienced larger mergers host more satellites.

Zoom In Zoom Out Reset image size

In this section, we will quantify the significance of the observed N_Sat–M_⋆,Dom relation. In Section 4.1 we define our adopted model to describe the relation, followed by the results of an MCMC analysis. In Section 4.2 we consider and account for a possible relationship between galaxy stellar mass and number of satellites, showing that the relationship remains in the residuals of an imposed M_⋆,Gal–N_Sat relation.

4.1. Defining the Model

For the purposes of simplicity, we will parameterize the relationship between M_⋆,Dom and N_Sat as a simple linear model, given by

$$\begin{eqnarray}&&N(\alpha ,{N}_{9.5};{M}_{\star ,\mathrm{Dom}})=\alpha ({\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}-9.5)+{N}_{9.5},\end{eqnarray} \tag{ 1 }$$

where we have normalized the model to the MW's value of ${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}=9.5$ (see Table 1). To estimate the probability of the data given our assumed model, we adopt the appropriate likelihood function for a linear model with intrinsic scatter (following Hogg et al. 2010):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & = & -\displaystyle \frac{1}{2}\sum _{i=1}^{n}\displaystyle \frac{{\left[{N}_{i,\mathrm{Sat}}-{N}_{i}(\alpha ,{N}_{9.5};{M}_{\star ,\mathrm{Dom}})\right]}^{2}}{{\sigma }_{i}^{2}}\\ & & +\mathrm{ln}({\sigma }_{i}^{2}),\end{array}\end{eqnarray} \tag{ 2 }$$

where the variance ${\sigma }_{i}^{2}$ is expressed as

$$\begin{eqnarray}&&{\sigma }_{i}^{2}={(\alpha \,{\sigma }_{{\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom},i}})}^{2}+{\sigma }_{{N}_{\mathrm{Sat},i}}^{2}+{\sigma }_{0}^{2}.\end{eqnarray} \tag{ 3 }$$

${\sigma }_{i}^{2}$, defined for each system i, takes into account uncertainties on both M_⋆,Dom and N_Sat, ${\sigma }_{{\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom},i}}$ and ${\sigma }_{{N}_{\mathrm{Sat},i}}$, respectively, as well as an intrinsic scatter term, which is assumed to be the same for all systems, ${\sigma }_{0}^{2}$.

Table 1. Nearby Galaxy Group Properties

Galaxy	${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Gal}}$	Dom. Satellite	${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}.\mathrm{Sat}}$	${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Acc}}$	NSat	References
	(M⊙)		(M⊙)	(M⊙)
Milky Way	10.79	LMC	9.50	8.95	7 ± 1	1, 2, 9, 15, 16
Messier 31	11.01	M33	9.51	10.3	12 ± 2	2, 10, 16, 25
Messier 81	10.75	M82	10.45	9.06	11 ± 2	2, 11, 14, 18, 20, 22
Messier 83	10.72	KK 208	∼7.81	8.62	4 ± 1	2, 4, 17, 26
Messier 94	10.72	M94-dw1	5.99	8.70	2 ± 0	21, 24, 27
Messier 101	10.77	NGC 5474	9.13	7.91	5${}_{-1}^{+3}$	2, 3, 6, 8, 12, 18, 20
Messier 104	11.30	dw1240-1118	∼7.66	11.0	15${}_{-4}^{+3}$	5, 6, 13, 23
Centaurus A	11.05	ESO324-024	∼8.13	10.6	14 ± 3	2, 7, 10, 19

Notes. The central galaxy stellar mass, name and mass of the most massive existing satellite, estimated total accreted stellar mass, and number of satellites within 150 kpc of the central with M_V < −9, for each of the eight nearby MW-mass galaxies composing our observational sample. The larger of the current dominant satellite mass and the accreted mass is bolded, which represents M_⋆,Dom for each system. Systems with no published stellar mass for the dominant satellite are denoted by a "∼" symbol; in these cases, we convert the published M_V by assuming a 1 M_⋆/L_V mass-to-light ratio. References: (1) Bell et al. (2008); (2) Bell et al. (2017); (3) Bennet et al. (2019); (4) Carrillo et al. (2017); (5) Cohen et al. (2020); (6) Carlsten et al. (2021); (7) Crnojević et al. (2019); (8) Danieli et al. (2017); (9) Deason et al. (2019); (10) D'Souza & Bell (2018b); (11) Harmsen et al. (2017); (12) Jang et al. (2020); (13) Jardel et al. (2011); (14) Karachentsev & Kudrya (2014); (15) Mackereth & Bovy (2020); (16) McConnachie (2012); (17) Müller et al. (2017); (18) Querejeta et al. (2015); (19) Rejkuba et al. (2011); (20) Sheth et al. (2010); (21) Smercina et al. (2018); (22) Smercina et al. (2020); (23) Tempel & Tenjes (2006); (24) Trujillo et al. (2009); (25) van der Marel et al. (2012); (26) M. Cosby et al. (2022, in preparation); (27) K. Gozman et al. (2022, in preparation).

Download table as:  ASCIITypeset image

We compute the posterior probability distributions for each model parameter given the data using the Python-based Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013). Figure 3 (left panel) shows data plotted against the likeliest model in blue, as well as the 16%–84% confidence regions. We also show the marginalized single and joint posterior distributions for the model parameters given the data (right panel, again in blue). For comparison, Figure 3 also shows results from the FIRE simulation suite; we defer description and discussion of this comparison to Section 6.

Zoom In Zoom Out Reset image size

The eight galaxies are best described by a tight relation, with intrinsic scatter well within the uncertainties (in fact, consistent with zero scatter in our best-fit model). The model with the highest likelihood of describing the observations is:

$$\begin{eqnarray}&&{N}_{{\rm{Sat}}}={4.73}_{-0.65}^{+0.76}({{\rm{log}}}_{10}{M}_{\star ,{\rm{Dom}}}-9.5)+{7.53}_{-0.60}^{+0.62}.\end{eqnarray} \tag{ 4 }$$

4.2. Stellar Mass Considerations

An important consideration in explaining the emergence of this relation is the range of stellar masses of the observed galaxy sample. New evidence suggests that the stellar mass of central galaxies may correlate with how many satellites they host (Carlsten et al. 2021), reflecting the nonzero slope of the SMHM relation at the 'MW mass'—i.e., galaxies with higher stellar masses should live, on average, in more massive DM halos and thus are predicted in hierarchical assembly to host more subhalos. Of the eight observed galaxies, only the MW and M31 have reliable estimates of their halo masses, inferred from the motions of individual stars (e.g., Cautun et al. 2020; Deason et al. 2021, and many others) or their full satellite populations (e.g., Watkins et al. 2010; Patel et al. 2018); thus, we cannot directly investigate the role of halo mass in setting N_Sat. We can, however, follow precedents set in the literature (e.g., Carlsten et al. 2021) and use galaxy stellar mass as a (scattered) proxy for halo mass to assess the possible intrinsic impact of galaxy mass on satellite number.

In Figure 4, we first show the stellar mass of each galaxy (see Table 1) against the estimated dominant merger mass (left). M_⋆,Dom is not expected to be particularly well correlated with galaxy mass (e.g., Bell et al. 2017; D'Souza & Bell 2018a), and we find this to be true for our observed sample. Using dark matter halo catalogs from Illustris-TNG and the abundance matching method, Carlsten et al. (2021) found an average relation between the stellar mass of the central galaxy, ${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Gal}}$, and the number of M_V < − 9 satellites within a 150 kpc projected radius, N_Sat—the same definition used in this paper. In Figure 4 (center), we show this average relation, and its large estimated scatter, along with the M_⋆,Gal and N_Sat for the eight observed galaxies (as well as galaxies from the FIRE simulation, which will be discussed in Section 6). The stellar masses assumed for this galaxy sample vary relatively substantially throughout the literature. Comparing available methods for stellar mass estimation (e.g., Bell et al. 2003; Querejeta et al. 2015), we estimate a typical uncertainty of 0.2 dex on each stellar mass. It is worth noting that the subsequent full possible 0.4 dex stellar mass range considered for each galaxy is more than half of the full stellar mass range of 0.58 dex for the sample. This is an unfortunate limitation of stellar mass comparison that cannot be avoided.

Zoom In Zoom Out Reset image size

To assess the potential impact of galaxy mass on satellite number, we evaluate the residual N_Sat for each galaxy from its predicted average, 〈N_Sat〉, estimated from its stellar mass (see Table 1) and the Carlsten et al. (2021) relation. Given the small dynamic range in stellar mass, yet steep M_⋆,Gal–N_Sat relation, small changes in M_⋆,Gal have a large change in expected N_Sat. It is important to incorporate these uncertainties from the stellar mass of the main galaxy. Accordingly, we run 10,000 Monte Carlo tests to convert the individual points to probability density space. For each test, we assess M_⋆,Dom and N_Sat − 〈N_Sat〉 for each of the eight systems as random draws from a Gaussian random variable with width equal to the estimated uncertainty in each direction.

We show the results of this Monte Carlo experiment in Figure 4 (right). Due to the large overlapping uncertainties, the eight individual measurements naturally cluster into three broad groups in M_⋆,Dom – $8.1\,\lesssim \,{\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}\,\lesssim \,8.9$, $9\,\lesssim \,{\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}\,\lesssim \,9.8$, and $10.2\,\lesssim \,{\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}\,\lesssim \,11$. For each of these three M_⋆,Dom bins, we calculate the average N_Sat − 〈N_Sat〉. A clear positive correlation is visible between M_⋆,Dom and N_Sat − 〈N_Sat〉 across the three bins.

We do note that the distribution of N_Sat − 〈N_Sat〉 for the eight observed galaxies appears shifted to negative N_Sat − 〈N_Sat〉, relative to the abundance-matched Illustris-TNG distribution. It is possible that this is a statistical effect, though our randomized treatment of the uncertainties should yield a relatively robust distribution. The other possibility is a difference in assumptions between the observations and the abundance matching applied to Illustris-TNG. The most likely of these seems to be a modestly lower stellar mass normalization (likely reflecting a lower stellar mass-to-light ratio) for the halo occupation model used on TNG halos relative to the nearby galaxy sample.

Regardless of the large uncertainties or potential mass-to-light differences, it seems clear that a strong relationship between number of satellites and mass of the dominant merger persists, even after taking into careful account a possible intrinsic relationship with stellar mass of the central. ⁸

Throughout the previous sections, we have established that a strong relationship exists between the number of satellites around nearby MW-mass systems and the mass of their dominant mergers. Of the eight observed systems considered in this work, two galaxies have estimates of both their first and second most dominant mergers, as well as usable evidence for which of their satellites may have existed prior to the most dominant event—the MW and M81. Both galaxies are currently experiencing their most dominant merger—the MW with the LMC, and M81 with M82—but their accreted masses (i.e., limit on their previous dominant mergers) have also been measured from their stellar halos (see Table 1). Recent work has shown promise in differentiating existing and newcomer satellites in both systems. ⁹ In this section, we will use this recent work to sketch the evolution of the MW and M81 satellite systems from their past to current dominant mergers, showing that their satellite systems evolve along the relation.

5.1. Recent MW Satellites

5.2. Recent M81 Satellites

Though no comparable-quality dynamical evidence exists for the satellites of M81 (and likely will not for the foreseeable future), rich existing resolved star data sets do exist for many M81 satellites, along with corresponding SFH analyses. M82 and NGC 3077 are both commonly assumed to be recent additions to the M81 system, as evidenced by M81's ongoing interaction with both (e.g., Yun et al. 1994; Okamoto et al. 2015; Smercina et al. 2020). Additionally, its third most massive satellite, NGC 2976, is thought to have interacted with M81 as recently as ∼1 Gyr ago (Williams et al. 2010) and therefore is also very likely to be a recent addition on a similar timescale. Given their very comparable infall times, we consider NGC 3077 and NGC 2976 to be past M82 satellites.

It has been long known that M81's global satellite population exhibits more active and recent star formation than is found in the MW's or M31's satellites (Weisz et al. 2008). Using the intuition built from the MW's and M31's satellites, both of which show clustering in their τ₉₀ timescales, we analyze the broad SFH properties of all eight available low-mass satellites within 150 kpc of M81 to assess possible clustering around the infall time of M82. SFHs for each of these satellites were calculated as part of the ANGST survey (Weisz et al. 2011). We use these published SFHs to evaluate τ₅₀ and τ₉₀ for each satellite and estimate uncertainties from the 16th–84th percentile uncertainties on the cumulative SFHs calculated by Weisz et al. (2011). Figure 5 (right) shows the estimated τ₅₀ and τ₉₀ for each of the eight satellites, color-coded by their V-band absolute magnitude, also presented in Table 2.

Zoom In Zoom Out Reset image size

Table 2. M81 Satellites

Galaxy	Rproj	MV	τ50	τ90
	(kpc)		(Gyr)	(Gyr)
BK3N	11	−9.59	3.9${}_{-0.5}^{+6.2}$	3.2${}_{-2.7}^{+0.5}$
KDG61	31	−13.87	12.9${}_{-10.9}^{+0.3}$	1.6${}_{-1.0}^{+2.0}$
FM1	61	−11.46	13.0${}_{-2.4}^{+0.3}$	4.1${}_{-1.3}^{+8.5}$
BK5N	72	−11.33	13.0${}_{-5.5}^{+0.0}$	10.8${}_{-9.3}^{+1.7}$
IKN	82	−11.51	11.7${}_{-3.5}^{+1.1}$	10.3${}_{-4.0}^{+0.8}$
KDG64	101	−13.43	9.3${}_{-6.3}^{+3.7}$	1.8${}_{-0.5}^{+2.7}$
KK77	102	−12.84	9.8${}_{-2.0}^{+3.2}$	3.4${}_{-2.0}^{+5.32}$
F8D1	119	−13.14	10.6${}_{-4.0}^{+4.0}$	1.8${}_{-0.1}^{+10.7}$

Notes. 50% and 90% star formation timescales for faint M81 satellites, calculated from their resolved SFHs (Weisz et al. 2011). Satellites are ordered by projected radii, calculated assuming a distance to M81 of 3.6 Mpc (Radburn-Smith et al. 2011), and absolute V-band magnitudes are given (from Karachentsev et al. 2000, 2001; Sharina et al. 2005; Georgiev et al. 2009). Satellites with τ₉₀ < 4 Gyr are bolded.

Download table as:  ASCIITypeset image

Of these eight satellites, two satellites have ancient τ₉₀ > 9 Gyr. The remaining six have τ₉₀ ≤ 4 Gyr, which is comparable to the timescale of M82's infall (e.g., Yun et al. 1994). While the uncertainties on these timescales are large (due to model uncertainty in stellar populations shallower than the main-sequence turnoff), three of these satellites—F8D1, KDG61, and KDG64—had very recent star formation (τ₉₀ ≲ 2 Gyr). We assert that these three satellites of M81, along with NGC 3077 and NGC 2976, very likely originated in the M82 system prior to infall into M81's halo. The three satellites with τ₉₀ between 3–4 Gyr—FM1, KK77, and BK3N—also have quite recent implied infall times, and could be more recent, given their uncertainties.

Given the existing evidence, we take NGC 3077, NGC 2976, F8D1, KDG61, and KDG64 to be likely former M82 satellites, and FM1, KK77, and BK3N to be possible former M82 satellites.

5.3. Evolution of the MW and M81 Satellite Systems

Having established in the previous sections that the satellite populations of MW-mass galaxies tightly correlate with the masses of the dominant mergers they have experienced, it is necessary to contextualize this result. How does this empirical relation compare to theoretical predictions?

To first order, CDM predicts that the number of subhalos a central halo hosts should scale with its virial mass. This can be intuited from the much richer satellite populations of massive galaxy clusters, relative to MW-mass groups. However, we are considering MW-mass galaxies, which should occupy halos with a much smaller range of virial masses than the difference between clusters and groups. Particularly in this narrower range of halo mass, where low-mass satellites hold essentially all of the statistical weight, there is considerable precedent for the deviation of satellite galaxy properties from CDM predictions when accounting for the complex impact of baryons (e.g., the "Missing Satellites" and "Too Big to Fail" problems; e.g., Macciò et al. 2010; Font et al. 2011; Brooks et al. 2013; Wetzel et al. 2016). Given this, it is most fruitful to compare the observed M_⋆,Dom–N_Sat relation against hydrodynamic simulations.

Recent results from the Illustris TNG50 simulation reproduce the predictions from CDM; mainly, a scattered relationship between N_Sat and virial mass (Engler et al. 2021). However, Engler et al. (2021) also found hints of a correlation between later assembly times and higher satellite number, albeit with high scatter. This may indicate that systems with later-time, larger mergers have more satellites—qualitatively similar to the trend seen in the observations. However, TNG50 does not resolve the more numerous low-mass classical satellites probed by the observations, making a fair comparison difficult. Ideally, the observations would be compared to a simulation with high enough resolution to capture the entire satellite population.

Currently, the highest-resolution baryonic simulation that meets the necessary criteria to compare complete satellite populations and merger histories in a sample of MW-mass galaxies is the FIRE (Hopkins et al. 2014, 2018) simulation suite. FIRE has resolved both the satellite galaxy populations (down to M_⋆ > 10⁵ M_⊙; Garrison-Kimmel et al. 2019) and stellar halo properties (Sanderson et al. 2018) of a sample of MW-mass galaxies, in both isolated and LG-like environments, with virial masses ranging from 0.8 − 2 × 10¹² M_⊙. ¹⁰

To directly compare FIRE to the eight observed systems, we take the satellite catalogs for each simulated system in FIRE (Samuel et al. 2020), assess the number of satellites with M_⋆ > 4 × 10⁵ M_⊙ (comparable to M_V < − 9; Smercina et al. 2018) and projected radius < 150 kpc for 1000 randomly drawn 3D orientations, and measure the average N_Sat across these 1,000 different sightlines—identical to the technique we used to construct the MW's and M31's satellite LFs (Section 2.1). We take M_⋆,Dom to be the larger of the peak stellar mass of each galaxies' largest z = 0 satellite or the total accreted mass—identical to the definition of M_⋆,Dom for the observed sample. We provide accreted masses, dominant satellite masses, and sightline-averaged satellite counts for the FIRE galaxies in Table 3 in Appendix C.

Table 3. FIRE Data

Galaxy	M⋆,Gal	${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}.\mathrm{Sat}}$	${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Acc}}$	NSat
	(1010 M⊙)	(M⊙)	(M⊙)
m12b	7.3	8.72	10.05	6${}_{-1}^{+2}$
m12c	5.1	9.91	9.59	14${}_{-2}^{+1}$
m12f	6.9	8.22	9.82	6 ± 1
m12i	5.5	8.10	9.47	6 ± 1
m12m	10.0	8.94	9.87	6${}_{-1}^{+2}$
m12z	1.8	8.81	9.14	6 ± 1
Romeo	5.9	9.84	9.20	10${}_{-2}^{+1}$
Juliet	3.4	9.30	9.48	11${}_{-3}^{+2}$
Romulus	8.0	8.69	10.21	4 ± 2
Remus	4.0	8.78	9.64	8${}_{-2}^{+1}$
Thelma	6.3	8.56	9.67	9 ± 1
Louise	2.3	8.56	9.13	8${}_{-2}^{+1}$

Notes. Relevant data for the FIRE sample. All data presented here are from the most current highest-resolution runs, and all include the subgrid physics detailed in Garrison-Kimmel et al. (2019). References: galaxy stellar mass (M_⋆,Gal)—Samuel et al. (2020) and Santistevan et al. (2020); dominant satellite mass (M_⋆,Dom.Sat), total accreted stellar mass (M_⋆,Acc), and average number of satellites within projected 150 kpc (N_Sat)—Samuel et al. (2020), this work. We ask that any works that wish to use the updated accreted masses also cite Sanderson et al. (2018).

Download table as:  ASCIITypeset image

We combine the LG-analog "ELVIS on FIRE" (treated as individual systems) and isolated "m12" suites, as our observational comparison set combines a number of different environments. ¹¹ Following Section 4, we use emcee to calculate the posterior probability distributions for the parameters of the model defined in Equation (1) given this mock-observed FIRE sample. We show the results of this analysis in Figure 3, where the results for the FIRE simulation are shown in orange. As with the observations (Section 4), we show the model with highest likelihood of describing the FIRE systems and the 16%–84% confidence region, as well as the marginalized single and joint posterior distributions for the model parameters (right panel). The difference between the FIRE and observed samples is striking. Unlike the observed systems, FIRE shows shows no correlation and large intrinsic scatter at a given M_⋆,Dom. In Appendix D we compare FIRE to the Auriga simulation to show that this result is, in general, independent of the differences between them.

We also show the FIRE systems in Figure 4 for comparison, with stellar masses taken from Samuel et al. (2020) and Santistevan et al. (2020; see Table 3). Surprisingly, the FIRE systems exhibit a positive correlation between M_⋆,Dom and M_⋆,Gal, in contrast with the observations and expectations (e.g., Bell et al. 2017). ¹² Conversely, FIRE does does not exhibit any discernible correlation between N_Sat and M_⋆,Gal. This lack of a trend between M_⋆,Gal and N_Sat in FIRE has been previously noted and is thought to be related to the smaller range of halo masses in FIRE (0.9–1.7×10¹² M_⊙) than, e.g., the Carlsten et al. (2021) analysis (0.8–8×10¹² M_⊙), which may allow for the more effective tidal disruption from more massive central galaxies to dominate any scaling of satellite number with virial mass (Samuel et al. 2020). We will discuss this further in Section 7.3. The FIRE N_Sat residuals exhibit essentially pure scatter in N_Sat − 〈N_Sat〉 for a fixed M_⋆,Dom, which in turn is entirely consistent with the TNG-abundance matching results. It is, however, completely unlike the correlation visible in the observed sample, despite a very comparable stellar mass range as seen in the left panel. We observe no statistical difference between the LG-like ELVIS hosts and the isolated m12 hosts. In addition, as M104 is also more massive (by a factor of ∼2 in M_⋆,Gal) than the range of galaxies produced in FIRE, we also conduct the comparison analysis between the observed sample excluding M104 and FIRE, and find identical results.

It is important to note that FIRE produces only half the logarithmic dynamic range in dominant merger mass experienced by MW-mass galaxies compared to the observations—${\mathrm{log}}_{10}{M}_{\odot }\,\sim \,9\mbox{--}10.4$ versus ${\mathrm{log}}_{10}{M}_{\odot }\,\sim \,8.3\mbox{--}11$. For reference, this lack of diversity is borne out in the comparable-resolution Auriga simulation as well (see Appendix D), though neither simulation imposes a selection constraint on the merger history. It is unclear whether this lack of diversity is related somehow to the by-design choice of specific halos in these cosmological zoom-in simulations thought to be capable of hosting an MW-mass disk-like galaxy (as cosmological simulations typically show a broader range more consistent with the observations; Deason et al. 2015b; D'Souza & Bell 2018a), or possibly indicates a more fundamental problem with halo occupation. An important first step toward better understanding this tension will likely be to work toward producing comparable-resolution simulations that encompass a larger range of halo masses.

This being said, it is entirely clear that even if FIRE were to encompass a broader range of halo masses, and were it shown that such a broader halo mass range is more consistent with the observed galaxy sample considered here, the considerable scatter in N_Sat at a fixed M_⋆,Dom (see Figure 3) would persist. This large scatter is present in other baryonic simulations (e.g., TNG50 and Auriga; Engler et al. 2021, Appendix D), regardless of the halo mass range used. This suggests that, while there may well be an inevitable scaling of satellite number with virial mass, the tight empirical M_⋆,Dom–N_Sat relation is not forthcoming in current galaxy formation models.

This brings us to a last, important consideration: despite the lack of a tight M_⋆,Dom–N_Sat relation, the satellite count normalization (given by N_9.5 in our adopted model) is almost identical in both FIRE and the observations. Rather than indicative of specific calibration on the part of the models to match the MW, we contend that this reflects the status of the MW's satellites as a long-standing benchmark for the community. This suggests that benchmarking against the MW, irrespective of merger history, has allowed high-resolution zoom-in simulations to relatively faithfully reproduce the MW's satellites while failing to fully capture the drivers of the scaling relation between MW-mass galaxies' satellite number and the mass of their dominant merger.

Our direct comparison of the satellite populations and dominant mergers of eight nearby systems benefited from impressive observational progress on both the satellite and stellar halo fronts. The next decade will undoubtedly see this sample of galaxies grow, encompassing systems in different environments and with different merger histories. These future observations will allow us to better understand and interpret the observational findings of this paper; mainly, the empirical relation between satellite number and dominant merger mass. Acknowledging the likely future refinement of our results, we aim to evaluate the scientific insight gained from the existence of this relation, and discuss what we find to be a number of important implications for galaxy formation.

7.1. Sketching the Problem

7.2. Virial Radius

7.3. Enhanced Tidal Disruption

7.4. Group-scale Feedback

7.5. Time of Accretion

7.6. Moving Forward

Motivated by the infall of MW satellites associated with the LMC and the broad diversity in the satellite populations of nearby MW-mass galaxies, throughout this paper we have, for the first time, investigated a possible relationship between these diverse satellite populations and their equally diverse merger histories as traced by their most dominant mergers. For this effort, we define a new metric for the mass of a galaxies' most dominant merger, M_⋆,Dom, as the larger of either its total accreted mass M_⋆,Acc, or the mass of its most massive current satellite, M_⋆,Dom.Sat. Using this new definition of M_⋆,Dom, and the satellite populations of eight nearby galaxies, we find the following:

• 1.  
  We find a strong observed correlation between M_⋆,Dom and the number of M_V < −9 satellites within a projected 150 kpc (N_Sat) of the eight nearby MW-mass galaxies in which these two properties have been best measured.
• 2.  
  Assuming a simple linear model between ${\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}$ and N_Sat, in concert with MCMC sampling methods, we find an unexpectedly tight relation, which remains even after accounting for possible intrinsic relationships with galaxy stellar mass. We conclude that the number of satellites around an MW-mass galaxy is a strong function of its largest merger partner.
• 3.  
  Using recent evidence from orbital and SFH modeling, we empirically demonstrate that both the MW's and M81's satellite populations have evolved along the observed M_⋆,Dom–N_Sat relation over the course of their current dominant mergers with the LMC and M82, respectively. Consequently, we estimate that more than 50% of M81's current satellites within a projected 150 kpc were accreted during its recent merger with M82.
• 4.  
  We compare the observed relation to theoretical MW-like systems from the FIRE simulation, assessing average satellite populations within projected radii <150 kpc across many random orientations. We find the simulated galaxies' behavior to be completely different from the observations. Instead of the tight observed relationship, FIRE exhibits no discernible relationship between M_⋆,Dom and N_Sat, with a large intrinsic scatter in N_Sat at a fixed M_⋆,Dom. However, it does produce the correct number of satellites for an MW-like M_⋆,Dom. Given the similarity between the results from DMO simulations, semi-analytic models, and these high-resolution hydrodynamic simulations, this tension suggests that, while current simulations can reproduce the small-scale structure observed around the MW, they are not fully capturing the diversity of merger histories essential to understanding the formation of other MW-mass galaxies and their satellites.
• 5.  
  The M_⋆,Dom–N_Sat relation directly implies that galaxies like M94 or M83, with difficult-to-explain sparse satellite populations, may represent the "baseline" present-day satellite population of an MW-mass galaxy, in the absence of any larger merger throughout its life.
• 6.  
  We speculate that some combination of heavily time-dependent satellite radial distributions, enhanced tidal disruption, and early-time group-scale feedback in the vicinity of a massive central may help to explain the emergence of the M_⋆,Dom–N_Sat relation. The "time of accretion" may be an important factor to consider when exploring the validity of these different physical mechanisms.

M_⋆,Dom and N_Sat are proxies for two of the most fundamental aspects of hierarchical galaxy formation theory—galaxy merger histories and the buildup of satellite galaxy populations. Though the infall of satellites during merger events may feel intuitive, the tight scaling relation described here cannot be explained as a natural consequence of traditional hierarchical structure formation. The solution must be sought in the context of galaxy formation physics, in which any trend is currently either highly scattered or nonexistent. While future observational facilities will deepen our understanding of the relationship between these proxies, as it currently stands, these results, which demonstrate a surprisingly acute difference between their predicted and observed connection, constitute an important tension that must be understood and addressed if we are to achieve a holistic understanding of galaxy formation. The ability of models to reproduce this relationship, including the general diversity in the observed merger histories of MW-like galaxies, may serve as a guide for successfully coupling small- and large-scale galaxy physics, including stellar feedback and environmental processes.

We thank the referee for a careful and thoughtful review that improved this paper. We also thank Andrew Wetzel, Monica Valluri, Jeremy Bailin, Mario Mateo, JD Smith, Julianne Dalcanton, Oleg Gnedin, and Ian Roederer for valuable discussions and feedback.

A.S. was supported by NASA through grant No. GO-14610 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. E.F.B. was partly supported by the National Science Foundation through grant 2007065 and by the WFIRST Infrared Nearby Galaxies Survey (WINGS) collaboration through NASA grant NNG16PJ28C through subcontract from the University of Washington.

This research uses cosmological zoom-in baryonic simulations of MW-mass galaxies from the Feedback In Realistic Environments (FIRE) simulation project, run using the Gizmo gravity plus hydrodynamics code in meshless finite-mass mode (Hopkins 2015) and the FIRE-2 physics model (Hopkins et al. 2018). We thank the FIRE team for helping provide us with electronic versions of the FIRE satellite catalogs, and for valuable advice on working with the data.

This research has made use of the following astronomical databases and catalogs: the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). The VizieR catalog access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in Ochsenbein et al. 2000. The Extragalactic Distance Database (Tully et al. 2009).

Software: Matplotlib (Hunter 2007), NumPy (Harris et al. 2020; van der Walt et al. 2011), Astropy (Astropy Collaboration et al. 2018), SciPy (Virtanen et al. 2020), emcee (Foreman-Mackey et al. 2013), corner (Foreman-Mackey 2016).

In addition to the eight galaxies considered in this paper, estimates of the satellite populations exist for two other galaxies with well-measured M_⋆,Dom: NGC 4565 and M51. NGC 4565 has an estimated accreted mass from the GHOSTS survey (Harmsen et al. 2017) that is substantially more massive than any existing satellite. M51 is currently interacting with its well-studied massive satellite, M51b, which has a similar mass to M81's dominant satellite, M82. It is unlikely that M51's current accreted mass exceeds M51b's stellar mass, so M51b's mass can be assumed to represent M_⋆,Dom. Both galaxies' satellite populations were studied by Carlsten et al. (2021), but they have large uncertainties stemming from inconclusive SBF distance measurements, similar to M104. However, both have significantly fewer confirmed satellites than M104, and much larger uncertainties on the faint-end of their LFs. Therefore, we treat them as lower limits, compared to the rest of the sample. Figure 6 shows M_⋆,Dom plotted against N_Sat, identical to Figure 2, but with the lower limits on N_Sat for NGC 4565 and M51 included.

Zoom In Zoom Out Reset image size

N_Sat limits for both NGC 4565 and M51 are consistent with the relation, though M51 lies particularly low relative to the other galaxies with its dominant merger mass (e.g., M31 and M81). If M51's census were complete, it would be a strong outlier from the N_Sat–M_⋆,Dom relation, and it will be important to understand in detail in future efforts. In light of this, we consider two possible factors to be particularly important to reflect upon and investigate in the future.

The first, as mentioned above, is that the census of M51's satellite population is incomplete, and possibly highly incomplete due to its merger stage. Watkins et al. (2015) showed that M51's interaction with M51b has distributed tidal debris to radii of at least a projected 50 kpc from its center (<100 kpc in 3D radius). This debris has a surface brightness of μ_V < 26 mag arcsec ⁻², which is brighter than the majority of the faint classical satellites in the LG (McConnachie 2012). Carlsten et al. (2020) showed that at least half of M51's possible satellites lie in this region dominated by bright tidal debris. Other recent studies of the radial distribution of classical satellites around MW-like hosts also suggest that this projected radial range (i.e., ≲50 kpc) should encompass a substantial fraction of the central's satellite population (e.g., Samuel et al. 2020). It is entirely possible that this bright in situ contamination uniquely impacts the completeness of current satellite searches around M51, relative to other galaxies with much lower surface brightness accreted halos.

A second possibility is that M51 is truly deficient in satellites within its central 150 kpc, due to the current stage of its merger with M51b. D'Souza & Bell (2021) showed that during a merger between an MW-mass galaxy and a massive satellite, the orbit of the massive satellite influences the bulk distribution of the satellite population. When the massive satellite is close to its first pericenter, its satellites are also close to pericenter and are unusually centrally concentrated (e.g., the MW's; see Samuel et al. 2020). At later merger stages, the newly arrived satellites tend to be toward their apocenters (for a timescale of a few gigayears), which often lie well outside the virial radius. As M51's merger with M51b is in a relatively late stage, it is very possible that it is temporarily deficient in satellites within its inner 150 kpc. If this is the case, one would expect a significant number of satellites at substantially larger distance from M51—a hypothesis that is straightforward to test with future observations.

Independent of the importance of either of these two factors, it is clear that M51's satellite population will be a critical case study to better understand the emergence of the M_⋆,Dom–N_Sat relation.

Here we assess the impact of normalizing M_⋆,Dom to M_⋆,Gal—a "dominant merger ratio." We adopt a slightly modified version of Equation (1):

$$\begin{eqnarray}&&\begin{array}{l}N(\alpha ,{N}_{9.5};{M}_{\star ,\mathrm{Dom}})=\alpha ({\mathrm{log}}_{10}{M}_{\star ,\mathrm{Dom}}-9.5+{\mathrm{log}}_{10}{M}_{\star ,\mathrm{Gal}})\\ \quad +{N}_{9.5}.\end{array}\end{eqnarray} \tag{ B1 }$$

We repeat our MCMC analysis from Section 4.1, using this updated model. Figure 7 shows the results of this analysis. While uncertainties on M_⋆,Dom/M_⋆,Gal are higher than for M_⋆,Dom alone, we find results that are very comparable with those presented in Figure 3: a tight positive relation, with intrinsic scatter consistent with zero.

Zoom In Zoom Out Reset image size

We also repeat this analysis of the dominant merger ratio with the FIRE simulations. We find a more significant correlation between M_⋆,Dom/M_⋆,Gal and N_Sat than M_⋆,Dom alone: while the distribution of slopes is very broad, the average slope is positive. However, the intrinsic scatter is still far higher than the observed relation. The increased correlation observed in Figure 7 could be driven at least in part by the strong correlation between M_⋆,Dom and M_⋆,Gal in FIRE (see Figure 4). It is also important to note that the dynamic range of M_⋆,Dom/M_⋆,Gal for the FIRE galaxies is even smaller than M_⋆,Dom: nearly all of the FIRE galaxies exhibit M_⋆,Dom/M_⋆,Gal in between the values for the MW (0.05) and M31 (0.2), compared to the range M_⋆,Dom/M_⋆,Gal = 0.008–0.5 for the observed sample.

While Sanderson et al. (2018) measured accreted masses for the FIRE sample, there were some differences in the mass resolution of the galaxies and the subgrid physics used between different simulated galaxies. The runs from which the satellite catalogs adopted in this paper were created (Samuel et al. 2020) to include some of the same systems from Sanderson et al. (2018) that have since been re-simulated at higher resolution, and all now include a subgrid prescription for turbulent metal diffusion not present in the original simulations. This affects both the isolated "m12" suite, which achieves the now-standard FIRE initial mass resolution of 7070 M_⊙ for star and gas particles, and the LG-analogous "ELVIS on FIRE" suite, which was re-simulated at twice the resolution. See Garrison-Kimmel et al. (2019) for more nuanced details of the updated simulations.

With these updates to the models, there is concern that the orbital properties of mergers may be slightly different than those in previous runs, and therefore differences could be present in the distribution of accreted material relative to the measurements of Sanderson et al. (2018). To address this, we use the star particle catalogs for each FIRE galaxy to estimate updated accreted masses. We use the same scheme as Sanderson et al. (2018) to effectively "tag" accreted particles; mainly, that they initially formed at a distance greater than 30 kpc from the host (d_form > 30 kpc) and do not belong to an identified satellite at z = 0. This selection scheme effectively separates accreted and in situ material, accurately accounts for the steep observed accreted density profiles and consequent dominant reservoir of accreted material at small radii, and produces accreted masses that are consistent with estimates using the observational techniques described in Section 2.2.

In Table 3 we present the updated accreted masses, dominant merger masses, and sightline-averaged satellite counts for galaxies in the FIRE simulation, which were used in our comparison to the observations. A description of the calculation of N_Sat can be found in Section 6. We stress that while future works are fully free to use and distribute the updated accreted stellar masses presented here, Sanderson et al. (2018) should also be cited, as we adopt the particle-tagging method devised in that work for identifying accreted star particles.

Currently, two fully hydrodynamic, high-resolution cosmological zoom-in simulations meet the necessary criteria to compare complete satellite populations and merger histories in a sample of MW-mass galaxies: the FIRE (Hopkins et al. 2018) and Auriga simulations (Grand et al. 2017). FIRE achieves higher resolution for its satellite galaxies, which is why it was chosen as the default simulation for model–observation comparison throughout the main body of the paper. However, Auriga also achieves impressively high resolution, and has simulated an even larger sample of MW-like systems. In this section, we compare these two simulations directly to show their very similar behavior in the context of comparing satellite populations and merger histories.

Each simulation has resolved both the satellite galaxy populations (Simpson et al. 2018; Garrison-Kimmel et al. 2019) and stellar halo properties (Sanderson et al. 2018; Monachesi et al. 2019) of a sample of MW-mass galaxies, with virial masses ranging from approximately 0.8–2 × 10¹² M_⊙. While FIRE resolves dwarf galaxies down to M_⋆ > 10⁵ M_⊙, Auriga only resolves satellites to M_⋆ > 4 × 10⁵ M_⊙. When comparing simulations to each other, we consider dwarf galaxies in FIRE to a limit of M_⋆ > 4 × 10⁵ M_⊙ to put the simulations on equal footing, noting that there may be other differences between them owing to their modest differences in resolution.

In Figure 8 we show the estimated mass of the most dominant merger (either accreted mass or largest existing satellite) plotted against the number of satellites within 300 kpc with stellar mass above 4 × 10⁵ M_⊙, for both the Auriga and FIRE simulations. As in Section 6, for the FIRE runs, we combine the LG analog and isolated "m12" suites, as our observational comparison set combines a number of different environments.

Zoom In Zoom Out Reset image size

For each simulation suite, we run emcee to evaluate the posterior probability distributions for each of the model parameters defined in Equation (1), identical to our approach in Sections 4 and 6. The "best" model, with parameters corresponding to the likeliest value for each marginalized posterior distribution, is shown in Figure 8 for both FIRE and Auriga, along with their respective 16%–84% confidence regions. Both models are consistent with no intrinsic correlation, with nearly flat slopes and large intrinsic scatter. We show the positions of the MW and M31 for comparison, adopting the M_⋆,Dom values given in Table 1 and R < 300 kpc satellite populations compiled by Garrison-Kimmel et al. (2019). While M31 is a relative outlier in its satellite populations, compared to both simulations, both FIRE and Auriga closely bracket the MW and M31 in M_⋆,Dom–N_Sat parameter space and possess a very similar range of dominant merger masses. This comparison highlights that different high-resolution models, with different subgrid prescriptions for baryonic physics, consistently meet the benchmark set by satellite populations in the LG, and yet fail to capture the emergence of a tight scaling relation between M_⋆,Dom and N_Sat.