A profile in FIRE: resolving the radial distributions of satellite galaxies in the Local Group with simulations

TL;DR

This study uses FIRE-2 simulations to analyze the radial distribution of satellite galaxies around Milky Way-like hosts, finding broad agreement with Local Group observations and predicting additional satellites for future discovery.

Contribution

First detailed analysis of satellite radial distributions in FIRE-2 simulations, including effects of baryonic physics and host variation, matching Local Group data.

Findings

Simulations resolve satellite survival and destruction for Mstar > 10^5 Msun.

Broad agreement with Local Group satellite distributions, independent of satellite mass.

Predicted 2-10 undiscovered satellites around the Milky Way.

Abstract

While many tensions between Local Group (LG) satellite galaxies and LCDM cosmology have been alleviated through recent cosmological simulations, the spatial distribution of satellites remains an important test of physical models and physical versus numerical disruption in simulations. Using the FIRE-2 cosmological zoom-in baryonic simulations, we examine the radial distributions of satellites with Mstar > 10^5 Msun around 8 isolated Milky Way- (MW) mass host galaxies and 4 hosts in LG-like pairs. We demonstrate that these simulations resolve the survival and physical destruction of satellites with Mstar >~ 10^5 Msun. The simulations broadly agree with LG observations, spanning the radial profiles around the MW and M31. This agreement does not depend strongly on satellite mass, even at distances <~ 100 kpc. Host-to-host variation dominates the scatter in satellite counts within 300 kpc of the hosts, while time variation dominates scatter within 50 kpc. More massive host galaxies within our sample have fewer satellites at small distances, likely because of enhanced tidal destruction of satellites via the baryonic disks of host galaxies. Furthermore, we quantify and provide fits to the tidal depletion of subhalos in baryonic relative to dark matter-only simulations as a function of distance. Our simulated profiles imply observational incompleteness in the LG even at Mstar >~ 10^5 Msun: we predict 2-10 such satellites to be discovered around the MW and possibly 6-9 around M31. To provide cosmological context, we compare our results with the radial profiles of satellites around MW analogs in the SAGA survey, finding that our simulations are broadly consistent with most SAGA systems.

Figures

Tables

Table 1. Table 1: Host galaxy properties and satellite galaxy counts

Name	M200m [${10}^{12}$ M⊙]	M∗ [${10}^{10}$ M⊙]	Nsat( 50 kpc)	Nsat( 100 kpc)	Nsat( 300 kpc)
MW	$\sim$1.4	$\sim$5	1 $\pm$ 0.5 (50%)	6 $\pm$ 0.5 (10%)	13 $\pm$ 0 (0%)
M31	$\sim$1.6	$\sim$10	2 $\pm$ 0.5 (25%)	5 $\pm$ 1.0 (20%)	27 $\pm$ 0.5 (<5%)
m12m	1.6	10.0	1 $\pm$ 1.0 (100%)	7 $\pm$ 3.0 (45%)	27 $\pm$ 2.0 (10%)
m12b	1.4	7.3	0 $\pm$ 0.2 (N/A)	3 $\pm$ 1.0 (35%)	11 $\pm$ 0.7 (5%)
m12f	1.7	6.9	0 $\pm$ 0.0 (N/A)	1 $\pm$ 1.0 (100%)	16 $\pm$ 1.0 (5%)
Thelma	1.4	6.3	1 $\pm$ 0.5 (50%)	6 $\pm$ 1.2 (20%)	17 $\pm$ 1.2 (10%)
Romeo	1.3	5.9	1 $\pm$ 0.7 (70%)	4 $\pm$ 0.7 (20%)	17 $\pm$ 1.0 (5%)
m12i	1.2	5.5	0 $\pm$ 0.5 (N/A)	3 $\pm$ 1.2 (40%)	12 $\pm$ 0.5 (5%)
m12c	1.4	5.1	1 $\pm$ 0.9 (90%)	8 $\pm$ 1.5 (20%)	23 $\pm$ 1.0 (5%)
m12w	1.1	4.8	0 $\pm$ 0.5 (N/A)	5 $\pm$ 0.9 (20%)	22 $\pm$ 1.5 (10%)
Juliet	1.1	3.4	1 $\pm$ 0.5 (50%)	8 $\pm$ 1.7 (20%)	20 $\pm$ 0.5 (5%)
Louise	1.2	2.3	1 $\pm$ 0.5 (50%)	8 $\pm$ 1.5 (20%)	23 $\pm$ 0.7 (5%)
m12z	0.9	1.8	1 $\pm$ 0.5 (50%)	7 $\pm$ 0.7 (10%)	17 $\pm$ 0.7 (5%)
m12r	1.1	1.5	1 $\pm$ 1.2 (120%)	5 $\pm$ 1.7 (35%)	14 $\pm$ 2.2 (15%)
Time variation			1 $\pm$ 0.6 (60%)	5 $\pm$ 1.3 (25%)	17 $\pm$ 1.1 (5%)
Host-to-host variation			1 $\pm$ 0.5 (50%)	5 $\pm$ 2.5 (50%)	17 $\pm$ 4.7 (30%)
Total variation across hosts+time			1 $\pm$ 1.0 (100%)	5 $\pm$ 3.0 (60%)	17 $\pm$ 6.0 (35%)

Table 2. Table 2: Parameters for fits to the ratio of subhalos in baryonic versus dark matter-only simulations in Equation 1. Cumulative distributions refer to the total number of subhalos enclosed as a function of 3D distance while differential distributions refer to discrete bins in 3D distance.

Subhalo selection method	$\alpha$	$d_0$ [kpc]	$d_1$ [kpc]
Cumulative distributions
M${}_{\mathrm{peak}}{}^{}>8\times {10}^8$ M⊙	0.7	11	89
M${}_{\mathrm{bound}}{}^{}>{10}^8$ M⊙	0.8	13	106
M${}_{\mathrm{bound}}{}^{}>{10}^7$ M⊙	0.8	2	98
Differential distributions
M${}_{\mathrm{peak}}{}^{}>8\times {10}^8$ M⊙	0.8	8	78
M${}_{\mathrm{bound}}{}^{}>{10}^8$ M⊙	0.9	21	95
M${}_{\mathrm{bound}}{}^{}>{10}^7$ M⊙	0.9	0	100

Equations

$$f\left(d\right)=\{\begin{array}{ll}0& 0\le d<d_0\\ \alpha \left[1-e^{-\frac{d-d_0}{d_1}}\right]& d\ge d_0\end{array}$$

$$f_{\mathrm{bound}}<0.32(N_{\mathrm{acc}}/1000{)}^{-0.8}$$

$$f_{\mathrm{bound}}<\frac{1.79}{f(c)}\left(\frac{ϵ}{r_{\mathrm{s},0}}\right)\left(\frac{r_{\mathrm{h}}}{r_{\mathrm{s},0}}\right)$$

Full text

††pubyear: 2019††pagerange: A profile in FIRE: resolving the radial distributions of satellite galaxies in the Local Group with simulations–C

1 Introduction

Dark matter dominates the matter content of dwarf galaxies by up to several orders of magnitude, making them ideal sites for small-scale tests of the standard paradigm for structure formation: cold dark matter (CDM) with a cosmological constant ($\Lambda$). CDM makes testable predictions for both the central mass profile of dwarf galaxies and their number density and spatial distribution around more massive host galaxies. However, on such small scales, tests of CDM require highly resolved observations that are only feasible within the nearby Universe. Fortunately, the Milky Way (MW) and Andromeda (M31) galaxies that make up the Local Group (LG) are host to populations of satellite dwarf galaxies which can provide quantitative tests of CDM on small scales.

LG satellite galaxies have been a source of significant tensions within the CDM model, largely stemming from comparisons of observations to dark matter-only (DMO) simulations that lack the effects of baryonic physics. Arguably the most famous of these tensions, the “missing satellites” problem, describes a discrepancy between the number of subhalos predicted by DMO simulations compared to the smaller number of luminous satellite galaxies observed around the MW (e.g. Moore et al., 1999; Klypin et al., 1999b). However, newer simulations that include the effects of baryonic physics through hydrodynamics and sub-grid models show agreement with the number of observed satellite dwarf galaxies in the LG, in part from enhanced tidal disruption of satellites by the baryonic disks of host galaxies (e.g. Brooks et al., 2013; Sawala et al., 2016; Wetzel et al., 2016; Garrison-Kimmel et al., 2018a; Kelley et al., 2018; Simpson et al., 2018; Buck et al., 2019). N-body simulations in conjunction with semi-analytic models of galaxy formation have also yielded similar results, showing similar radial distributions for observable satellites (e.g. Macciò et al., 2010; Font et al., 2011). Simultaneously, a better understanding of observational incompleteness has also been critical in alleviating the missing satellites tension (Tollerud et al., 2008; Walsh et al., 2009; Hargis et al., 2014; Kim et al., 2018).

In addition to over-predicting the number of satellites, DMO simulations predict too many dense, massive (“too-big-to-fail") satellites (Boylan-Kolchin et al., 2011, 2012), and satellites with steeper (“cuspier") central density profiles than seen in observations (Navarro et al., 1996). Again, baryonic effects in simulations are a pathway to reconciling these problems because stellar feedback acts to redistribute the central dark matter and “core-out" the density profile of dwarfs (e.g. Mashchenko et al., 2008; Chan et al., 2015; Oñorbe et al., 2015; El-Badry et al., 2016; Dutton et al., 2016; Garrison-Kimmel et al., 2017b; Fitts et al., 2017).

Baryonic effects are also crucial for understanding the predicted phase space coordinates of satellites around simulated MW/M31-like galaxies. This phase space information can be used to infer the formation history of satellites and rigorously test CDM predictions. For example, orbit modeling of the Large Magellanic Cloud (LMC) has provided evidence that it is undergoing its first pericentric passage around the MW, and this may partly be why it is still able to form stars (Besla et al., 2007; Kallivayalil et al., 2013). The phase space distribution of LG satellites have further challenged the CDM model because the MW’s satellite galaxies appear to be arranged in a thin, planar structure that is coherently rotating, and a similar structure has been found around M31 (Lynden-Bell, 1976; Metz et al., 2007; Conn et al., 2013; Ibata et al., 2013; Pawlowski, 2018).

Furthermore, the Satellites Around Galactic Analogs (SAGA) survey is broadening our understanding of LG satellites by targeting satellites of MW analogs within 20-40 Mpc of the LG (Geha et al., 2017). Their goal is to obtain a complete census of satellites around 100 MW analogs, down to the luminosity of the Leo I dwarf galaxy (M${}_{r}<−12.3$; M${}_{*}\approx 5\times 10^{6}$ M*⊙*). This will make it possible to connect LG satellite galaxies with a large sample of satellite populations, providing a statistically robust cosmological context to interpret LG galaxy formation and evolution.

Satellite dwarf galaxies can be used to study the effects of environment on galaxy formation as well. Even before satellites accrete onto their host, they are preprocessed by interactions with other dwarf galaxies that are bound to them in small groups (e.g. Zabludoff & Mulchaey, 1998; McGee et al., 2009; Wetzel et al., 2013; Hou et al., 2014). Once they fall into their host halo, satellites can be tidally disrupted into diffuse stellar structures by their host. For instance, the Sagittarius dwarf galaxy is being actively torn apart into a stellar stream within the MW’s halo (e.g. Lynden-Bell & Lynden-Bell, 1995; Belokurov et al., 2006). As satellites orbit in the halos of their host galaxies, they are thought to be ram pressure-stripped of their gas, causing their star formation to be subsequently suppressed (e.g. Gunn & Gott, 1972; Fillingham et al., 2016). The MW and M31 may exert some of the strongest observed environmental effects on their satellite populations: most of their satellites are gas-poor and no longer forming stars, making them an interesting case study for environmental effects (e.g. Einasto et al., 1974; Mateo, 1998; Grcevich & Putman, 2009; McConnachie, 2012; Slater & Bell, 2014; Spekkens et al., 2014; Wetzel et al., 2015).

Given the unique ability to measure full 3D positions and velocities of LG satellites, and thus infer their orbital histories, the LG also provides a fertile physical testing ground for numerical evolution and disruption of subhalos in simulations. Historically, it has been difficult to use simulations to interpret observations of LG satellites because baryonic simulations have only recently begun to produce dwarf galaxies that do not suffer from numerical over-merging. Simulations of satellites undergoing tidal disruption have revealed that the most critical simulation parameters in dynamically resolving satellites are spatial and mass resolution (e.g. Carlberg, 1994; van Kampen, 1995; Moore et al., 1996; Klypin et al., 1999a; van Kampen, 2000; Diemand et al., 2007; Wetzel & White, 2010; van den Bosch & Ogiya, 2018). While current large-volume simulations can offer a sizable sample of satellite galaxies, their dwarf galaxies have limited resolution both in terms of particle mass and gravitational force softening, which curbs their usefulness in tests that require accurate tidal disruption and survival of satellites.

Understanding formation and evolution is contingent on resolving the radial distribution of satellites as a function of distance from their hosts in cosmological simulations. Higher resolution, ‘zoom-in’ simulations are now providing satellite populations that are sufficiently well-resolved for studying LG satellite populations in detail. The main questions this paper aims to answer are:

• •

Do cosmological zoom-in baryonic simulations reproduce the observed radial distributions of satellites around the MW, M31, and MW analogs?

• •

Do the radial profiles reflect physical disruption from the host galaxy and/or numerical disruption inherent in the simulations?

• •

How do radial profiles in hydrodynamic simulations differ from those in DMO simulations?

• •

If the simulations are representative of the LG, how complete are observations of dwarf galaxies out to large distances around the MW and M31?

This paper is organized as follows: in Section 2 we describe the simulations used and how satellites were selected from them, in Section 3 we describe the observational data set used, in Section 4 we present our results on radial profiles with comparisons of the hydrodynamic simulations to both observations and dark matter-only simulations, and a discussion of the conclusions and implications is given in Section 5.

2 Simulations

2.1 FIRE simulation suite

We use cosmological zoom-in baryonic simulations from the Feedback In Realistic Environments (FIRE) project111https://fire.northwestern.edu/, run with the upgraded FIRE-2 (Hopkins et al., 2018) numerical implementations of fluid dynamics, star formation, and stellar feedback. The FIRE-2 simulations use a Lagrangian meshless finite-mass (MFM) hydrodynamics code, GIZMO (Hopkins, 2015). The MFM method allows for hydrodynamic gas particle smoothing to adapt based on the density of particles while still conserving mass, energy, and momentum to machine accuracy. Gravitational forces are solved using an improved version of the N-body GADGET-3 Tree-PM solver (Springel, 2005), and the gravitational force softening of gas particles automatically adapts to match their hydrodynamic smoothing length.

The FIRE-2 simulations invoke realistic gas physics through a metallicity-dependent treatment of radiative heating and cooling over $10-10^{10}$ K, including free-free, photo-ionization and recombination, Compton, photo-electric and dust collisional, cosmic ray, molecular, metal-line, and fine-structure processes, accounting for 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe). The cosmic UVB background is included using the Faucher-Giguère et al. (2009) model, in which HI reionization occurs early ($z_{\rm reion}\sim 10$). The simulations that we use also model sub-grid diffusion of metals via turbulence (Hopkins, 2016; Su et al., 2017; Escala et al., 2018).

Star formation occurs in gas that is self-gravitating, Jeans-unstable, cold (T $<10^{4}$ K), dense ($n>1000$ cm*-3*), and molecular (following Krumholz & Gnedin 2011). Each star particle represents a single stellar population under the assumption of a Kroupa stellar initial mass function (Kroupa, 2001), and we evolve star particles according to standard stellar population models from STARBURST99 v7.0 (Leitherer et al., 1999). The simulations explicitly model several stellar feedback processes including core-collapse and Type Ia supernovae, continuous stellar mass loss, photoionization, photoelectric heating, and radiation pressure.

For all simulations, we generate cosmological zoom-in initial conditions at $z=99$ using the MUSIC code (Hahn & Abel, 2011), and we save 600 snapshots from $z=$99 to 0, with typical spacing of $\lesssim$25 Myr.

We use two suites of simulations in this paper. The first is the Latte suite of individual MW/M31-mass halos introduced in Wetzel et al. (2016). Latte consists of 7 hosts with halo masses M${}_{\rm 200m}=1-2\times 10^{12}$ M*⊙* (where ‘200m’ indicates a measurement relative to 200 times the mean matter density of the Universe), selected from a periodic volume of length 85.5 Mpc. Gas and star particles have initial masses of 7070 M*⊙, though because of stellar mass loss, at $z=0$ a typical star particle has mass $\approx$ 5000 M⊙. Dark matter particles have a mass resolution of m${}_{\rm dm}=3.5\times 10^{4}$ M⊙*. Dark matter and stars have fixed gravitational softening (comoving at $z>9$ and physical at $z<9$): $\epsilon_{\rm dm}=40$ pc and $\epsilon_{\rm star}=4$ pc (Plummer equivalent). The minimum gas resolution (inter-element spacing) and softening length reached in each simulation is $\approx 1$ pc.

In this paper we introduce two new hosts into the Latte suite: m12w and m12r. We select them using the same criteria as the Latte suite: M${}_{\rm 200m}(z=0)=1-2\times 10^{12}$ M⊙ and no neighboring halos of similar mass ($>3\times 10^{11}$ M⊙) within at least 5 R200m, to limit computational cost. However, for these two halos we impose an additional criterion: each must host an LMC-mass subhalo. Specifically, within the initial dark-matter-only simulation, we select halos that host (only) one subhalo within the following limits at $z=0$: maximum circular velocity V${}_{\rm circ,max}=92\pm 12$ km/s, distance $d=51\pm 40$ kpc, radial velocity $v_{\rm rad}=64\pm 17$ km/s, tangential velocity $v_{\rm tan}=314\pm 60$ km/s. These criteria are centered on the observed values for the LMC (e.g., Kallivayalil et al., 2013; van der Marel & Kallivayalil, 2014), though we use a wider selection window than the observational uncertainties to find a sufficient sample in our cosmological volume, which for this sample is a periodic box of length $172$ Mpc with updated cosmology to match Planck Collaboration et al. (2018): $h=0.68$, $\Omega_{\Lambda}=0.69$, $\Omega_{m}=0.31$, $\Omega_{b}=0.048$, $\sigma_{8}=0.82$, $n_{\rm s}=0.97$. The zoom-in re-simulations use the same resolution as the existing Latte suite (given the slightly different cosmology, dark-matter particles have slightly higher mass of m${}_{\rm dm}=3.9\times 10^{4}$ M⊙). While we select these halos to have LMC-like subhalos in the pilot dark-matter-only simulation, when we re-run with baryonic physics the details of the satellite orbit (in particular the orbital phase) do change. m12w’s most massive satellite has M${}_{*}=8\times 10^{8}$ M⊙ and at $z=0$ is at $d=248$ kpc, having experienced pericentric passage of 78 kpc 2.4 Gyr ago ($z=0.19$). m12r’s most massive satellite has M${}_{*}=2.8\times 10^{9}$ M⊙ and at $z=0$ is at $d=390$ kpc, having experienced pericentric passage of 30 kpc 0.7 Gyr ago at $z=0.05$. We will examine the dynamics of these LMC-like passages in upcoming work (Chapman et al., in preparation).

In addition to the Latte suite, we include one additional individual host (m12z), selected to have a slightly lower halo mass at $z=0$ and simulated at a higher mass resolution of m${}_{\rm baryon,ini}=4200$ M*⊙* (Garrison-Kimmel et al., 2018a).

We also use the ELVIS on FIRE suite of two simulations, which selected halos to mimic the separation and relative velocity of the MW-M31 pair in the LG (Garrison-Kimmel et al., 2018a). These simulations have $\approx 2\times$ better mass resolution than the Latte suite: the Romeo & Juliet simulation has m${}_{\rm baryon,ini}=3500$ M*⊙* and the Thelma & Louise simulation has m${}_{\rm baryon,ini}=4000$ M*⊙*.

All simulations assume flat $\Lambda$CDM cosmologies, with slightly different parameters across the full suite: $h=0.68-0.71$, $\Omega_{\Lambda}=0.69-0.734$, $\Omega_{m}=0.266-0.31$, $\Omega_{b}=0.0455-0.048$, $\sigma_{8}=0.801-0.82$, and $n_{\rm s}=0.961-0.97$, broadly consistent with Planck Collaboration et al. (2018).

2.2 Halo finder

We identify dark-matter (sub)halos using the ROCKSTAR 6D halo finder (Behroozi et al., 2013a). We identify halos according to their radius that encloses 200 times the mean matter density, R200m, and keep those with bound mass fraction $>0.4$ and at least 30 dark matter particles. We generate a halo catalog at each of the 600 snapshots for each simulation. We then construct merger trees using CONSISTENT-TREES (Behroozi et al., 2013b). For numerical stability, we generate halo catalogs and merger trees using only dark matter particles.

We then assign star particles to each (sub)halo in post-processing as follows (adapted from the method described in Necib et al. 2018). Given each (sub)halo’s radius, Rhalo, and Vcirc,max as returned by ROCKSTAR, we first identify all star particles whose position is within 0.8 Rhalo (out to a maximum radius of 30 kpc) and whose velocity is within 2 Vcirc,max of each (sub)halo’s center-of-mass velocity. We then keep star particles (1) whose positions are within 1.5 R90 (the radius that encloses 90 per cent of the mass of member star particles) of both the center-of-mass position of member stars and the dark matter halo center (thus ensuring the galaxy center is coincident with the halo center), and (2) whose velocities are within 2 $\sigma_{\rm vel}$ (the velocity dispersion of member star particles) of the center-of-mass velocity of member stars. We then iteratively repeat (1) and (2) until M*∗, the sum of the masses of all member star particles, converges to within 1 per cent. We keep all halos with at least 6 star particles and average stellar density $>300$ M${}_{\odot}\,$kpc-3*.

We examined each galaxy in our sample at $z=0$ by eye and found that this method robustly identifies real galaxies with stable properties across time; in particular, it reliably separates true galaxies from transient alignments between subhalos and stars in the stellar halos of the MW-mass hosts. All of the subhalos (within 300 kpc of their host) that we analyze are uncontaminated by low-resolution dark matter.

2.3 Satellite selection

We refer to the MW- and M31-mass galaxies in our simulations as “hosts” and their surrounding populations of dwarf galaxies with M${}_{*}>10^{5}$ M*⊙* within 300 kpc as “satellites". Each of the eight Latte+m12z simulations contains a single isolated host while each of the two ELVIS on FIRE simulations contains two hosts in a LG-like pair, with their own distinct satellite populations. This provides a total of 12 host-satellite systems to study and compare to observations. Table 1 summarizes properties of these systems. Host galaxies have stellar masses M${}_{*}\sim 10^{10-11}$ M*⊙* and dark matter halos M${}_{\rm h}=0.9-1.7\times 10^{12}$ M*⊙*. Host stellar mass is measured by computing the stellar mass enclosed by a 2D radius in the plane of the host disk and a height above and below the plane that together define a cylinder containing 90 per cent of the total stellar mass within a sphere of radius 30 kpc around the host galaxy.

Our satellite selection of M${}_{*}>10^{5}$ M*⊙* corresponds to a minimum of $\sim$20 star particles and a peak halo mass (throughout their history) of M${}_{\rm{peak}}>8\times 10^{8}$ M*⊙* (or $\gtrsim 2.3\times 10^{4}$ dark matter particles prior to infall). We expect subhalos that contain satellite galaxies with M${}_{*}\sim 10^{5}$ M*⊙* to be both resolved in the simulations (Hopkins et al., 2018) and nearly complete in observations (e.g. Koposov et al., 2007; Tollerud et al., 2008; Walsh et al., 2009; Tollerud et al., 2014; Martin et al., 2016), so we choose this as our lower mass limit to make reasonable comparisons to the MW and M31 (but see Sections 4.8 and 4.9 for further discussion on potential incompleteness in the LG). For this analysis, we consider only the total (3D) radial distance from satellite to host galaxy, leaving a complete study of the full 3D positions and the problem of satellite planes for future work. For further details on the stellar masses, velocity dispersions, dynamical masses, and star-formation histories of dwarf galaxies in our simulations, see Garrison-Kimmel et al. (2018a, 2019).

3 Observations

We use the compilation of observed stellar masses of LG satellite galaxies in Garrison-Kimmel et al. (2018a), in which they assume stellar mass-to-light ratios from Woo et al. (2008) where available, and elsewhere use M∗/L${}_{\rm V}=1.6$ (Martin et al., 2008; Bell & de Jong, 2001). We apply the same stellar mass limit and host-satellite distance limit to the MW and M31 as in our simulation satellite criteria (M${}_{*}>10^{5}$ M*⊙* and d < 300 kpc). For the satellite galaxies around the MW we take sky coordinates and distances with uncertainties from McConnachie (2012). To model the effects of uncertainties in observed distances, we sample MW satellite distances 1000 times assuming Gaussian distributions for the uncertainties to generate a median radial profile with scatter around the MW (Figure 1).

We exclude the Sagittarius dwarf spheroidal galaxy from our MW sample, because it is undergoing significant tidal interactions and we do not believe our halo finder would correctly identify it as a subhalo in the simulations. We also include two more recently discovered ultra-diffuse satellite dwarf galaxies of the MW: Crater 2 (D${}_{\odot}\sim 118$ kpc; Torrealba et al. 2016) and Antlia 2 (D${}_{\odot}\sim 130$ kpc; Torrealba et al. 2018), bringing the total number of MW satellites considered to 13. Using the nominal mass-to-light ratio of 1.6 we estimate the stellar masses of these additional galaxies to be M${}_{*}\sim 2.6\times 10^{5}$ M*⊙* for Crater 2 and M${}_{*}\sim 3.4\times 10^{5}$ M*⊙* for Antlia 2.

For the satellite galaxies around M31, we use sky coordinates where available from McConnachie (2012) and apply the same stellar mass and distance restrictions, leaving us with a total of 28 satellite galaxies. To obtain the 3D radial profiles of M31’s satellites with uncertainties, we sample 1000 line-of-sight distances from the posterior distributions published in Conn et al. (2012), where available. However, several M31 satellites do not have published distance distributions: M32, NGC205, IC10, And VI, And VII, And XXIX, LGS 3, And XXXI, and And XXXII. In the cases of M32 and NGC205, they are too close to M31 to reliably determine their distances, so we assume they have the same line-of-sight distance distribution as M31 itself. Positions on the sky, distances, and distance uncertainties for And XXXI and And XXXII are taken from their discovery paper (Martin et al., 2013). For the remaining satellites without line-of-sight distance posteriors, we sample the distances published in McConnachie (2012), assuming Gaussian distributions on the uncertainties.

Figure 1 shows the cumulative number of satellite galaxies around the MW and M31 as a function of 3D distance from the host, and the shaded regions represent estimated scatter in these profiles when we consider uncertainties in line-of-sight distance. While the sample for M31 includes 28 total satellite galaxies, when we include uncertainties the median number of satellites within 300 kpc is 27. The resulting 68 per cent scatter averaged across distance from host in LG radial profiles is $\pm$0.3 satellites for the MW and $\pm$1.2 satellites on average for M31. We discuss comparisons to the scatter in simulation profiles in Section 4.2.

Comparisons to the LG must also be understood in terms of observational completeness. However, the observational data used for comparison to the simulations in this work have not been completeness-corrected. From the Pan-Andromeda Archaeological Survey (PAndAS, McConnachie et al., 2009), the satellite population around M31 is complete to within 150 kpc (projected) of M31 down to half-light luminosities L${}_{1/2}>10^{5}$ L*⊙* (Tollerud et al., 2012). This includes our lowest satellite galaxy stellar mass limit ($10^{5}$ M*⊙*), so we think we are making a fair comparison to M31 at least within 150 kpc (where we find evidence for tidal disruption of galaxies by the host, see Sections 4.5 and 4.6). However, if we assume that our simulations are representative of the LG we find that there may be more galaxies to discover around M31 beyond 150 kpc (see Section 4.9). Given that M31 already has a somewhat high number of satellite galaxies compared to the MW, this could potentially make M31’s satellite population larger than those of the simulations used here.

Completeness around the MW is complicated by varied survey coverage and seeing through the Galactic disk (Kim et al., 2018). However, these sources of incompleteness are likely to affect only satellite galaxies fainter than classical dwarf galaxies and therefore they are unlikely to significantly change the results of this work. Of some concern is the proper identification of diffuse or low surface brightness galaxies (especially through the disk), but this is already being addressed using Gaia data to identify dynamically coherent stellar structures (like the Antlia 2 galaxy included in this work). We cannot preclude the possibility of further observational incompleteness down to our lowest stellar mass cut out to 300 kpc around the MW and M31. For this reason, we present comparisons at multiple (higher) stellar mass limits (Sections 4.1 and 4.4) and make predictions for the numbers of satellites to potentially be discovered around the MW and M31 (Sections 4.8 and 4.9).

We also compare our simulations and observations of the LG to results from the Satellites Around Galactic Analogs (SAGA) survey (Geha et al., 2017). SAGA targets MW analogs down to the luminosity of the Leo I dwarf galaxy (M${}_{r}<−12.3$; M${}_{*}\approx 5\times 10^{6}$ M*⊙*), and the initial results include the 2D radial profiles of satellite galaxies around 8 MW analogs within 20-40 Mpc of the LG. For more details on how we made this comparison, see Section 4.3.

4 Results

We analyze satellite galaxy positions in our simulations over time using halo catalogs from 11 snapshots, taken over $z=0-0.1$ ($\sim$1.3 Gyr) in steps of $\Delta z=0.01$. We do this for each of the 12 simulated hosts, providing a total of 132 radial distributions of satellite galaxies at different times to study. In the inner halo, a typical satellite can undergo a full orbit in under 1 Gyr, while it may take $\sim$3-4 Gyr for a complete orbit in the outer halo. This time baseline allows us to time-average over satellite orbits to minimize sampling noise over at least 1/4 of an orbit, which is especially important at small distances where satellites spend the least amount of time. Our choice is motivated by a compromise between sampling sufficiently across orbital histories and avoiding systematic redshift evolution (compared with $z=0$) in the satellite populations. We find that time-averaging is critical for obtaining accurate results (see Section 4.2 for results on time variation in radial profiles).

4.1 Radial profiles

Figure 2 shows the cumulative number of satellite galaxies with M${}_{*}>10^{5}$ M*⊙* as a function of 3D distance from the host, or the radial profile, for each individual host-satellite system. The solid, colored lines are the simulated median radial profile across $z=0-0.1$, and the shaded regions show the 68 per cent and 95 per cent variation over time. The median number of satellites within 300 kpc for the simulated hosts ranges from 11-27, consistent with the observed total number of M${}_{*}>10^{5}$ M*⊙* satellite galaxies within 300 kpc of the MW (median N${}_{\rm sat}=13$) and M31 (median N${}_{\rm sat}=27$) today. Hosts are ordered by stellar mass with m12m being the most massive (M${}_{*}\approx 10^{11}$ M*⊙) and m12r the least massive (M${}_{*}\approx 1.5\times 10^{10}$ M⊙*). The number of satellite galaxies does not have an obvious correlation with host mass. The hosts show a wide range of profile shapes: m12m, m12c, m12w, and Louise closely follow M31, while Thelma, Romeo, m12i, m12z, and m12r more closely follow the MW.

Figure 3 summarizes the key result of this work: the radial profiles of satellite galaxies around the 12 hosts in our simulations span the observed radial distributions of satellites in the LG. Figure 3 aggregates all of our simulated profiles at three different satellite stellar mass thresholds: M∗ $>10^{5}$ M*⊙* (left), M∗ $>10^{6}$ M*⊙* (middle), and M∗ $>10^{7}$ M*⊙* (right). In the top panels we show the median and scatter across all 132 radial profiles simultaneously. The median radial profile for all simulated hosts (blue) lies on top of the median LG observations at distances $<150$ kpc (where observational completeness is more secure), and at larger distances it lies between the MW and M31 profiles. The median for paired hosts (orange) is slightly above the total median (blue), while the median for isolated hosts (pink) is slightly below the total median. However, the paired and isolated medians are still within the total 68 per cent scatter across all the simulations.

The 68 per cent scatter in the simulations overlaps the 68 per cent scatter in MW observations (shown in Figure 1) at nearly all distances, and the MW’s median profile is always within the 95 per cent simulation scatter. M31’s median profile lies within the 95 per cent simulation scatter at nearly all distances. However, M31 appears to have a slight excess of satellites compared to the 68 per cent simulation scatter at small distances ($<$50 kpc) and large distances ($>$250 kpc) for all satellite M∗ thresholds, though the uncertainties in M31’s profile at small distances are relatively high (see Section 3 for more details). The 95 per cent scatter in simulations overlaps with the 68 per cent scatter in LG profiles at all distances (not shown here, but see Figure 1). We also a differentially-binned radial distribution for satellites with M∗ $>10^{5}$ M*⊙* in Appendix B, where we also see general agreement the LG and our simulations. We conclude that our simulation sample broadly agrees with and spans the profiles around the MW and M31.

In the bottom panels of Figure 3, we normalize the total simulation median and scatter to the MW and M31 radial profiles. To calculate the simulated-to-observed ratios, we divide the time-averaged radial profile of each of the 12 simulated hosts by 1000 sampled observational radial profiles of the MW or M31. Thus, the scatter in each of the bottom panels is from simulated host-to-host variation as well as observational uncertainties. We find that the MW ratio is consistent with unity at the 68 per cent level at nearly all distances for satellites with M∗ $>10^{5}$ M*⊙* and M∗ $>10^{6}$ M*⊙. This consistency breaks down at distances $<$150 kpc for M∗ $>10^{7}$ M⊙* given the presence of the LMC and SMC, which are currently near their pericentric passage around the MW.

The M31 ratio is consistent with unity at the 95 per cent level at most distances for satellite M∗ $>10^{5}$ M*⊙* and M∗ $>10^{7}$ M*⊙, while for M∗ $>10^{6}$ M⊙* the upper scatter in the ratio typically reaches $\sim$0.8. The M31 ratio is consistent with unity at the 68 per cent level within 50-150 kpc of the host for satellite galaxies with M∗ $>10^{5}$ M*⊙*. Beyond 50 kpc, the median M31 ratio is typically $\sim$50 per cent across the different mass thresholds, indicating that it has a somewhat large satellite galaxy population compared to our average simulation. This excess of satellite galaxies around M31 relative to the simulations is consistent at all distances, suggesting that M31 may just have more satellites overall, which may mean that its host halo mass is higher than in our simulated sample. The M31 ratio is most consistent with unity for our lowest mass bin and within 50-100 kpc. We interpret this, along with our resolution tests in Appendix A, as evidence that we are resolving our sample well even at these lower satellite masses.

Finally, to statistically test whether our simulations’ radial profiles are consistent with the LG, we perform a two sample Kolmogorov-Smirnov (KS) test between the median profiles of the LG and each simulated host’s profiles (at all 11 snapshots) for satellite galaxies with M${}_{*}>10^{5}$ M*⊙. The KS test compares the overall shape of the radial profile, and is less sensitive to the absolute number of satellites than taking a ratio between simulations and observations. We calculate the KS statistic for each of the 11 snapshots over $z=0-0.1$, and we quote the percentage of snapshots where a simulation was inconsistent with either the MW or M31. The KS test results show that a few of the simulations are inconsistent with being drawn from the same distribution as the MW at a significance level of 95 per cent: m12f (83 per cent), m12m (27 per cent), m12i (18 per cent), and m12w (9 per cent). Only m12r (9 per cent) is inconsistent with M31’s distribution, and only at one of the 11 snapshots. We also use the Anderson-Darling (AD) test to check these results and maximize sensitivity to the tails of the radial distributions. With the AD tests, we achieve essentially the same results as the KS tests. We also repeat the KS and AD tests for satellite galaxies with M${}_{*}>10^{6}$ M⊙*, and found that none of the simulated profiles are inconsistent with the MW or M31 at the 95 per cent level, possibly indicating even better agreement at higher masses and that simulations and observations are well resolved and complete in this mass range.

4.2 Scatter across hosts versus across time

Table 1 summarizes host galaxy mass, number of satellites per host within representative distances, and the scatter over time in each host’s radial profile. We quantify the scatter in radial profiles using the 68 per cent scatter around the median number of satellites with M${}_{*}>10^{5}$ M*⊙* within a given distance from their host. To understand the importance of time versus host-to-host scatter, we compare the radial profile scatter within individual hosts over time (the pink and orange shaded regions from Figure 2), scatter among hosts after their time dependence has been averaged out (the solid, colored median lines in Figure 2), and total scatter among all hosts and snapshots simultaneously (the blue shaded region of Figure 3). We quote the 68 per cent scatter about the median in absolute number of satellites and also quote scatter as a percentage relative to the median to give an idea of the fractional variation. We consider the scatter at three different distances (50, 100, and 300 kpc) to measure time dependence over the full range of the radial profiles.

First, we consider the scatter in the total number of satellite galaxies within 300 kpc. The combined scatter across all hosts and snapshots within 300 kpc is $\pm$6 satellites, or a 35 per cent variation when normalized to the median of 17 satellites. The host-to-host scatter after averaging time dependence out is $\pm$4.7 satellites (30 per cent), whereas the average scatter over time for an individual host is much lower at $\pm$1.1 satellites (5 per cent). Thus we find that total scatter at large distances is dominated by to host-to-host variations rather than time dependence.

Within 100 kpc, the combined scatter across hosts and time is $\pm$3 satellites (60 per cent), while the host-to-host scatter is $\pm$2.5 satellites (50 per cent), and the time scatter is $\pm$1.3 satellites (25 per cent). The increased fractional significance of time scatter is likely caused by the relatively small amount of time that satellites spend near pericenter of their orbits. Within 50 kpc we see that this effect is exacerbated: the combined scatter across hosts and time is $\pm$1 satellite (100 per cent), while host-to-host scatter is $\pm$0.5 satellites (50 per cent), and time scatter is $\pm$0.6 satellites (60 per cent). We conclude that at large distances ($\gtrsim$100 kpc) the total scatter across all 132 radial profiles is dominated by host-to-host variation, and at small distances ($\lesssim$50 kpc) the total scatter is dominated by time dependence from satellite orbits.

4.3 Comparison to the SAGA survey

We also compare our simulated and LG profiles to the on-the-sky projected radial profiles of 8 MW analogs in the SAGA survey. To match the SAGA luminosity limit, we select satellite galaxies for comparison to SAGA in our simulations, the MW, and M31 by requiring them to have stellar masses above the value of Leo I, M${}_{*}\approx 5\times 10^{6}$ M*⊙*. We generate 2D projections of the simulations along 1000 lines of sight for each of the 12 host-satellite systems at 11 snapshots, from which we compute the median and scatter across the simulated sample. For M31 satellites, we use only their projected on-the-sky distances from M31, assuming a line-of-sight distance to M31 of 780 kpc. For the MW, we use the 3D positions of the satellites and their line-of-sight distance uncertainties to generate 2D realizations from 1000 lines of sight as we did for the simulations.

Figure 4 (left) shows the observed 2D profiles for SAGA hosts, the MW, and M31. Most SAGA systems have fewer satellite galaxies compared to the MW and M31, which could be an effect of the broad mass selection function used in SAGA to choose MW analogs within uncertainties on the MW’s stellar mass (Geha et al., 2017). Because our simulations show only slightly higher satellite counts in our LG pairs compared with isolated hosts (see Figure 3), this implies that the SAGA selection of isolated hosts is unlikely to be a significant cause of difference as compared with the LG. The MW profile lies in the middle of the SAGA sample, and its scatter via line-of-sight averaging spans most of the range between the average SAGA profile and M31’s profile within 200 kpc. M31 still has a relatively large number of satellites compared to the SAGA observations at all distances (especially beyond 150 kpc), but two of the SAGA hosts have numbers of satellites approaching M31’s profile.

Figure 4 (right) shows the SAGA profiles compared to the simulations. The blue line is the median and the shaded regions show the 68 per cent and 95 per cent scatter in simulations. The scatter in simulated 2D profiles is mainly due to host-to-host variation and line-of-sight averaging, while time variation contributes a negligible amount of scatter in projection. At distances $>$100 kpc, three of the eight SAGA hosts are at or below the 95 per cent simulation limits. The SAGA average lies within the 68 per cent simulation scatter at most distances, though still slightly below the simulation median for distances $>$100 kpc. The best agreement between the SAGA average and the simulations is at small distances ($<$100 kpc), where they overlap the most. Overall, the simulation scatter encompasses five of the eight SAGA profiles and we find reasonable agreement among SAGA results, the LG, and our simulations in projection.

4.4 Dependence on satellite galaxy stellar mass

In this section, we examine in more detail whether our results within small distance ($\lesssim$100 kpc) depend on the stellar mass of satellite galaxies. This is a test of how our simulations compare to observations across our satellite mass range, and because higher-mass satellites are better resolved in both stellar mass and halo mass, it is also a test of dependence on resolution. Here, we assume that satellites with larger stellar masses inhabit more massive dark matter halos, but this may not always be true given scatter in the galaxy stellar mass-halo mass relation (Garrison-Kimmel et al., 2017a; Fattahi et al., 2018). We use the number of satellite galaxies within small distances as our summary statistic because this is where we expect to see the most prominent effects of tidal disruption and perhaps numerical over-merging in simulations. However, given the small numbers of satellites within 50 kpc of the hosts in both our simulations and the observations, we choose 100 kpc as the limiting distance as a reasonable trade-off between testing at small distances and obtaining reasonable statistics.

The top panel of Figure 5 shows the median number of simulated satellite galaxies within 100 kpc of their host across all hosts and snapshots as a function of the lower limit on satellite stellar mass, compared to the MW and M31. We consider satellite stellar mass limits from M${}_{*}>10^{5}$ M*⊙* up to M${}_{*}>10^{7}$ M*⊙, the highest stellar mass bin where we still have sufficient statistics in our simulations. The observed medians for the MW and M31 are $\lesssim$2$\times$ higher than the simulation median. This difference for satellites with M${}_{*}>10^{7}$ M⊙* could be caused by the presence of the LMC and SMC near their pericenters around the MW, which is not typical in a time-averaged sense. Even so, the 95 per cent simulation scatter always encompasses the observations, and the 68 per cent scatter mostly contains the MW and M31 lines.

In the bottom panels, we normalize our simulations to the MW and M31 observations, by sampling from both the simulated hosts and observational uncertainties simultaneously. In general, satellite galaxies with smaller stellar masses reside in less massive dark matter halos, so they are resolved with fewer star and dark matter particles. Therefore, in the absence of confounding numerical artifacts, we might expect our simulations to show increasingly fewer satellites relative to observations at lower stellar masses if we are reaching our resolution limit. Interestingly, we find the best agreement with observations when we include our lowest mass satellite galaxies (M${}_{*}>10^{5}$ M*⊙). The simulated-to-observed ratios are always consistent with unity at the 95 per cent level, but are only consistent with unity at the 68 per cent level when we include satellites with M${}_{*}>10^{5-6}$ M⊙*. The trend in the ratios as a function of minimum satellite stellar mass considered is relatively flat, though our simulations may not be producing as many higher-mass satellites as the LG. This is broadly consistent with results from Garrison-Kimmel et al. (2018a), who examined all satellites out to 300 kpc and found that most hosts are consistent with the MW and M31’s satellite population is only slightly larger than the simulations. Therefore, relative to observations, our simulations do not suffer from obvious over-destruction of satellites at the stellar masses that we consider.

To more thoroughly analyze numerical resolution, we test for convergence of the radial distributions of subhalos in these simulations compared to those from lower resolution simulations in Appendix A. There, we show that our (high resolution) simulations are converged to within $\sim 20$ per cent on average, and the 68 per cent (host-to-host) scatter is consistent with 100 per cent convergence at distances $>30-40$ kpc. We note that this exercise suffers from the effects of an additional disruptive effect in the low-resolution simulations: because the low-resolution host galaxies have $\sim 2\times$ larger stellar masses, their subhalos may be more easily tidally stripped or destroyed as they orbit close to the host. This may have the effect of making the high-resolution simulations appear less converged, at least at small distances from the host. We conclude that subhalos hosting the satellite galaxies in our high-resolution simulations are sufficiently resolved for tests of the satellite populations’ radial distributions. For a more nuanced discussion of convergence and additional resolution tests, see Appendix A.

4.5 Dependence on host mass

We test whether our results for satellite galaxies with M${}_{*}>10^{5}$ M*⊙* are sensitive to the stellar and halo masses of the host galaxies within 50 and 100 kpc. Figure 6 (left) shows the median number of satellite galaxies within 50 and 100 kpc of the host as a function of host stellar mass. The simulations agree well with the MW, and while the simulation trends lie below M31, the scatter for the most massive simulated host (m12m) is still consistent with M31. Within both 50 and 100 kpc of the host there is a negative trend in the number of satellites as a function of host galaxy stellar mass, and 4 hosts have no satellites at all (median over time) within 50 kpc. These 4 hosts all have stellar masses $\gtrsim 5\times 10^{10}$ M*⊙*, which is the average host stellar mass for the simulations. We interpret this and the trend lines as evidence for enhanced tidal destruction of satellites in our simulations due to the increased gravitational potential from the more massive hosts’ baryonic disks.

Figure 6 (right) shows the median number of satellites within 50 and 100 kpc as a function of host halo mass (M200m). When controlling for host halo mass, the time variation or scatter in the simulations is consistent with the MW and M31. However, M31 lies above both the simulation trends and the MW lies slightly above the 100 kpc trend line. M31 on the other hand, lies above the simulation trends, but still within the simulation scatter. The trend in Nsat(d<100 kpc) as a function of host halo mass is slightly less steep than as a function of host stellar mass for the simulations. The correlation of Nsat(d<50 kpc) with host mass is stronger for stellar mass (Pearson correlation coefficient: $r_{*}=-0.32$) than it is for halo mass ($r_{200\rm{m}}=-0.22$). The correlations of Nsat(d<100 kpc) with each type of host mass are: $r_{*}=-0.35$ and $r_{200\rm{m}}=-0.43$. Within 300 kpc (not shown) we find little to no correlation: $r_{*}=0.14$ and $r_{200\rm{m}}=0.07$. We interpret the larger correlation with host stellar mass within 50 kpc and steeper trend with host stellar mass within 100 kpc as the destructive tidal effects of the host baryonic disk manifesting at sufficiently small distances. Since host stellar mass is more correlated with satellite count within 50 kpc, we conclude that host stellar mass is a better predictor of the total number of surviving satellite galaxies within 50 kpc of the host, where we expect disk effects to be strongest.

Naively, we might expect the number of satellites at any distance to correlate positively with halo mass, and because M∗ correlates with M200m, we might also expect a similar correlation with stellar mass. Both the negative trend with host stellar mass and the lack of satellites around the more massive galactic disks suggest instead that the baryonic disk is depleting the satellite population at small distances. However, we note that because of the correlation between host M∗ and M200m in our simulations (see Figure 13), we cannot strictly disentangle the tidal effects of the separate disk and halo components of the host independently in our analysis. Despite this uncertainty, we find it physically plausible that tidal destruction of satellites can negate our initial expectations, at least for satellites closer to the host galaxy, consistent with results presented in Garrison-Kimmel et al. (2017b) and Kelley et al. (2018) that show a lack of satellites or subhalos at small distances in the presence of a disk potential. This also explains the lack of correlation between the number of satellites within 300 kpc and host halo mass (also noted in Fig 3 of Garrison-Kimmel et al. 2018a as a function of host halo virial mass): while increasing halo mass increases the number of expected satellites, the correlation of host stellar mass with host halo mass and the tidal destruction from the host disk act to cancel out this dependence, at least within the limited host mass range that we explore with our simulations.

We also note that while our simulated hosts have a wide range of stellar masses (M${}_{*}\sim 10^{10-11}$ M*⊙), they were selected over only a narrow range in host halo mass (M${}_{200\rm{m}}\sim 1-2\times 10^{12}$ M⊙*). Therefore, our sample is missing MW/M31-like host galaxies with much larger (or smaller) halo masses, but with stellar masses that scatter into our sample’s range. Hosts with more extreme halo masses like this could potentially lead to a less negative correlation of Nsat with host M∗ within 100 kpc.

4.6 Comparison with dark matter-only simulations

The seven Latte simulations also have dark matter-only (DMO) versions run with the same number of DM particles and the same force softening222However, the DMO simulations have DM particles with slightly higher masses of m${}_{\rm dm}=4.2\times 10^{4}$ M*⊙* due to the lack of baryons. We correct for this by multiplying DMO subhalo masses by $1-f_{\rm b}$ to account for the mass that would be otherwise relegated to baryons given the cosmic baryon fraction ($f_{\rm b}\equiv\Omega_{\rm b}/\Omega_{\rm m}$) of our baryonic simulations.. We compare the DMO versions to the baryonic simulations in order to investigate the effects of baryonic physics on the radial profiles of subhalos. To compare with the baryonic simulations, we find that satellite galaxies with M${}_{*}>10^{5}$ M*⊙* have typical peak dark matter halo masses M${}_{\rm peak}\gtrsim 8\times 10^{8}$ M*⊙*.

We select subhalos in the DMO and baryonic simulations by requiring them to be within 1000 kpc of their host and to have M${}_{\rm peak}>8\times 10^{8}$ M*⊙* so their halos are approximately as well-resolved as baryonic satellites down to M${}_{*}\sim 10^{5}$ M*⊙*. We then average the radial profiles of each host-subhalo system over $z=0-0.1$ using all available snapshots (67 total) for improved subhalo statistics at small distances. We compute the ratio of a host’s baryonic-to-DMO profiles for each host individually, and then examine the median and scatter across hosts. We compute the ratio as a function of distance for both cumulative and differential subhalo counts: N($<d$) and N($d_{1}<d<=d_{2}$), respectively. Figure 7 shows the cumulative ratio of baryonic-to-DMO subhalos (top) and the differential ratio (bottom). The line and shaded regions are the median and scatter showing only host-to-host variation, as the time dependence has been averaged out prior to taking the ratio.

Within $\sim$100 kpc from the hosts, there are $\lesssim$50 per cent the number of baryonic subhalos compared to DMO subhalos, and this continues to rapidly drop as distance decreases until the (median) ratio reaches zero at 10-15 kpc. As Garrison-Kimmel et al. (2017b) studied extensively using embedded disk potentials in DMO simulations of m12i and m12f, this is almost entirely due to the presence of the additional gravitational potential from the disk in the baryonic simulations. Here, we provide a more robust sample of simulations where we also time-average for each host, which is critical given how little time satellites spend near pericenter. The scatter within 100 kpc is greater than the scatter at 200-1000 kpc, due to a few hosts that have a number of baryonic subhalos closer to their number of DMO subhalos at small distances.

At large distances ($\gtrsim$200 kpc), the median ratios of baryonic-to-DMO subhalos flatten to $\sim$0.7 for the cumulative case and $\sim$0.8 for the differential case. This indicates that baryonic effects can reduce the masses of halos even at large distances from the host by $\sim$20-30 per cent as compared with DMO simulations. The overall reduction of substructure in the baryonic simulations relative to the DMO simulations is likely due to a combination of various baryonic effects such as reionization through our UV background and environmental effects like ram-pressure stripping and interactions with large scale structure such as filaments (Benítez-Llambay et al., 2013). Any of these processes may act to blow out gas from the galaxies in our baryonic simulations, shallowing their gravitational potential significantly in lower-mass galaxies like dwarfs, which in turn allows for easier removal of dark matter mass through gravitational interactions. Sawala et al. (2017) also found that the abundance of for subhalos with masses below $10^{9.5}$ M*⊙* in the APOSTLE simulations was reduced at all distances out to 200 kpc from the hosts. For the largest distances they examine, between 50 and 200 kpc, Sawala et al. (2017) found a reduction in substructure abundance of 23 per cent. This is similar to our results for the ratio of baryonic-to-DMO differential profiles between about 200 to 1000 kpc, where where we see a reduction in substructure abundance of about 20 per cent.

We provide fits to the ratio of baryonic-to-DMO subhalo counts as a function of distance that may be used to estimate the number of subhalos containing satellite galaxies (M${}_{*}>10^{5}$ M*⊙*) in other DMO simulations. We fit to the median ratio across hosts, and use the 68 per cent variation in the ratio as uncertainty on the fitted median values333The $z=0$ snapshots of baryonic m12i, m12f, and m12m are publicly available at ananke.hub.yt for comparison to individual hosts.. In Table 2, we also explore fits using other subhalo mass cuts. We fit the median of the cumulative and differential baryonic-to-DMO ratios as a function of distance ($d$):

[TABLE]

Where $\alpha$ is the asymptotic value of the ratio for infinitely large $d$, $d_{0}$ is the inner cutoff where the ratio goes to zero, and $d_{1}$ is the distance within which the ratio sharply declines. For the cumulative profile shown we find: $\alpha=0.7$, $d_{0}=11$ kpc, and $d_{1}=89$ kpc. For the differential profile shown we find: $\alpha=0.8$, $d_{0}=8$ kpc, and $d_{1}=78$ kpc. Table 2 shows these parameters for other fits using instantaneous bound halo mass for subhalo selection (not shown in Figure 7).

We find that, as expected, the fitted baryonic-to-DMO subhalo count ratios tend towards close to unity at large distances and drop to zero near the baryonic disk boundary. The fits indicate that even at arbitrarily large distances from the host, the baryonic subhalos are subject to additional destructive baryonic effects. The decline in the fitted ratios within $\sim$100 kpc is strikingly sharp: the cumulative and differential ratios both go to zero within $\sim 10$ kpc, indicating the physical boundary of intense gravitational effects from the baryonic disk. We see the same general trends in fits, for both cumulative and differential ratios, across the three different subhalo selection methods we use.

Our results agree with studies that have found that satellite survival depends on host-satellite distance at pericentric passage (e.g. D’Onghia et al., 2010; Zhu et al., 2016; Sawala et al., 2017; Garrison-Kimmel et al., 2018b; Nadler et al., 2018; Rodriguez Wimberly et al., 2019). We note that the destruction that we see is somewhat less extreme than in Garrison-Kimmel et al. (2017b), who used two of our baryonic simulations (m12i and m12f) and found no surviving subhalos at $z=0$ within $\sim$20 kpc of the host. Our results here are more robust given the larger host sample and that we time-average the profiles.

Kelley et al. (2018) examined the destructive effects of an analytical disk+bulge potential embedded in DMO simulations, where the analytical potential was allowed to realistically grow over time to match the MW’s potential at $z=0$. They found the ratio of subhalo counts that were subject to the embedded potential relative to subhalo counts that were not subject to the additional potential to be $\sim$1/3 within 50 kpc of their hosts. We find that our baryonic simulations are more efficient at destroying subhalos within 50 kpc, with a median ratio of baryonic-to-DMO subhalo counts of $\sim$1/5 at this distance. This could mean that additional baryonic effects, such as supernovae, in our simulations lead to enhanced modulation of the baryonic-to-DMO ratio. However, the simulations used in Kelley et al. (2018) were calibrated to the mass of the MW and may not capture the full effects of our wider mass range which encapsulate more massive M31-like galaxies as well.

Sawala et al. (2017) performed a similar comparison of the radial distributions of substructure in baryonic and DMO simulations, averaging over time and four hosts from the APOSTLE simulations. However, the baryonic disks of the hosts in their simulations are $\sim 2\times 10^{10}$ $M_{\odot}$, which is lower than the average disk masks of our hosts. Thus, based on our results from Section 4.5, we expect to see more substructure destroyed around our hosts than Sawala et al. (2017) found. They found a baryonic-to-DMO ratio of subhalo counts of $\sim 1/2$ at $d<10$ kpc, and $\gtrsim 3/4$ at $d>50$ kpc. By comparison, we see a much smaller median baryonic-to-DMO ratio of zero within $10$ kpc of our hosts, but the host-to-host scatter reaches as high as $\sim 1/5$ at $d<10$ kpc. At $d>50$ kpc, the scatter in our ratio varies from $\sim 1/5\textendash 1$ and at $d>100$ kpc it is $\gtrsim 1/2-1$. Newton et al. (2018) repeated this exercise and found subhalo depletion similar to Sawala et al. (2017): their baryonic-to-DMO ratio ranged from $\sim 1/2$ at small distances and rose to $\sim 4/5$ at large distances ($R_{200}$) from the host. Considering the differences in the stellar masses of the host disks between our simulations, the larger subhalo depletion we see at small distances compared to that from Sawala et al. (2017); Newton et al. (2018) is unsurprising, and we note that far from the host disk our results are more similar to each other.

4.7 Radial concentration

We further quantify satellite radial profiles using their shape, which we refer to as radial concentration. A profile with higher concentration generally has more of its satellites at small distances than at large distances from the host. We parameterize the concentration of our simulated and observed radial profiles using two metrics: $R_{90}/R_{50}$, the ratio of the distances enclosing 90 per cent and 50 per cent of the total number of satellite galaxies around a host, and $R_{90}/R_{10}$ to be sensitive to variations in satellite counts at smaller distances.

We analyze the concentration of 3D profiles considering LG satellites, baryonic simulation satellites, and DMO simulation subhalos that are within 300 kpc of their host. We measure concentration of the baryonic profiles for satellite galaxies with M${}_{*}>10^{5}$ M*⊙* around each of the 12 baryonic hosts and in the LG. DMO subhalos were selected as in Section 4.6, by requiring M${}_{\rm peak}>8\times 10^{8}$ M*⊙* for each of the 7 available DMO hosts. We also analyze the concentration of profiles in 2D projection for LG satellites, simulated baryonic satellites, and SAGA survey satellites with M${}_{*}>5\times 10^{6}$ M*⊙*. We report the concentration of each simulated host as the median over 11 snapshots from $z=0-0.1$, and the observed LG values as the median across 1000 sampled profiles.

Figure 8 (left) summarizes $R_{90}/R_{50}$ concentration measurements for the baryonic simulations, DMO simulations, the LG, and the SAGA survey. $R_{90}/R_{50}$ does not significantly differentiate baryonic simulations from DMO simulations. M31’s $R_{90}/R_{50}$ agrees with both the baryonic and DMO simulations, but the MW’s $R_{90}/R_{50}$ is slightly higher than the baryonic simulations, and is more consistent with the DMO simulations. However, we do find that $\sim$10-30 per cent of individual snapshots for half of the baryonic hosts (m12b, m12c, m12r, m12z, Romeo, and Thelma) have $R_{90}/R_{50}$ values that are at least as concentrated as the MW. This suggests that the MW has a slightly more concentrated profile shape relative to the median values for each baryonic simulation host. In 2D projection, the MW appears more highly concentrated than the baryonic simulations and M31. Most of the 2D SAGA profiles over the baryonic simulation distribution, but two SAGA systems have much higher concentration that is closer to the MW and one SAGA system has a concentration nearly twice that of the MW.

Figure 8 (right) shows $R_{90}/R_{10}$ concentration measurements for the simulations and observations. The baryonic simulations cover a broader range of values for $R_{90}/R_{10}$ than they do for $R_{90}/R_{50}$. DMO simulations tend to have systematically higher average $R_{90}/R_{10}$ than the baryonic simulations. Thus, the primary difference between baryonic and DMO profiles lies in the fraction of satellites at small distances ($\lesssim$100 kpc), where the DMO simulations have a larger fraction of their subhalos. This is consistent with the results of Section 4.6, where we show that the largest discrepancies between baryonic and DMO profiles occur within $\lesssim$100 kpc of the hosts. The $R_{90}/R_{10}$ values for the MW and M31 are consistent with the baryonic simulations and lie outside the range of DMO values. The 2D profile span an even broader range in $R_{90}/R_{10}$ than the 3D profiles. The SAGA systems are broadly consistent with the baryonic simulations, with a few more SAGA systems lying at the high end of the baryonic distribution. The MW in projection is also near the higher end of the baryonic simulations, and M31 appears much more concentrated than anything else.

The concentrations of the MW and M31 profiles generally overlap with the concentrations of the simulated baryonic profiles. Considering incompleteness in M31’s satellite population, if there are more M31 satellites to discover beyond 150 kpc, it could potentially push M31’s $R_{90}$ higher. This could increase M31’s concentration to a point where it is discrepant with the baryonic simulations. However, when using the $R_{90}/R_{50}$ metric, the MW is slightly more radially concentrated and therefore less consistent with the baryonic simulations than the DMO simulations. The MW in 2D projection appears more concentrated than most of the baryonic simulations, and M31 in projection is more concentrated than anything else. The SAGA profiles mostly overlap the projected baryonic simulation profiles, with a few SAGA systems having higher concentration more like the MW. Our results indicate that the MW may not be as much of a high-concentration outlier as previously thought: Yniguez et al. (2014) noted that the MW had a larger concentration than all of their DMO simulations. Concentration depends strongly on observational completeness assumptions though, and this may hint that there are more satellites just above M${}_{*}=10^{5}$ M*⊙* remaining to be discovered at farther distances from the MW as we explore next.

4.8 Implications for incompleteness around the Milky Way

While the MW and M31 profiles agree quite well out to 150 kpc, the MW appears to have a larger proportion of its satellite galaxies at small distances than both our simulations and M31. This could be a peculiarity of the MW profile, or it may be hinting at more satellites remaining to be discovered beyond 150 kpc from the MW. For example, the difference in shape could be due to the current presence of the LMC and the SMC near their pericenters around the MW (Kallivayalil et al., 2013). While Yniguez et al. (2014) found that potential incompleteness in the census of MW satellites meant there could be $\sim$10 classical dwarf satellite galaxies remaining to be discovered, which could bring the MW into better agreement with M31.

We expect observations of MW satellites to be complete down to at least M${}_{*}\sim 10^{5}$ M*⊙* within 150 kpc and out of the plane of the disk. Beyond this distance and through the disk the completeness may be uncertain, as evidenced by the discovery of Antlia 2, which had been obscured by the MW disk. Here, we focus on implications for incompleteness without considering the effects of seeing through the MW’s disk. While our theoretical results are suggestive, a more in-depth account of observational completeness for ‘classical’ dwarf galaxies also depends on the surface brightness distribution of the population and their on-the-sky positions with respect to the Galactic plane (or any other foreground structure). Our simulations can provide more detailed predictions for these effects on the completeness of the satellite population, especially through the use of Gaia-like mocks (Sanderson et al., 2018), which we plan to pursue in future work.

To investigate potential incompleteness in observations of the MW’s satellites (that have not been completeness-corrected), we examine how many additional satellites we would expect to find around the MW based on our simulations that match the MW profile out to 150 kpc. We choose simulated profiles for comparison by requiring them to have the same number of satellites within 150 kpc as the median value for the MW, which is 10 for M${}_{*}>10^{5}$ M*⊙*. One host meets this criteria at four snapshots (m12z), and four hosts meet this criteria at a single snapshot each (m12w, m12r, Romeo, and Juliet), providing a total of 8 matched profiles.

Figure 9 (left) shows the range of simulated profiles that match the MW at 150 kpc compared to the observed MW profile, for satellites with M${}_{*}>10^{5}$ M*⊙. The simulations agree remarkably well with the MW below 150 kpc, which further strengthens our claim that if we match the profile at this distance, then we are accurately resolving survivability of satellites closer to the host. Notably, beyond $\sim$150 kpc the simulation profiles are systematically higher than the MW profile. In total, the simulation median profile has 5 more satellites than the MW median profile within 300 kpc. The lower 68 per cent (95 per cent) limits on the simulation profile imply that there may be at least 4 (2) more satellites at 150-300 kpc from the MW. If our simulations are representative of the real MW, then based on the median simulation profile, we predict that there should be 5 more satellites with M${}_{*}>10^{5}$ M⊙* within 150-300 kpc of the MW.

We expect observational completeness to be better at higher satellite stellar masses, so we repeat this exercise for satellites with M${}_{*}>10^{6}$ M*⊙* to check if the agreement between simulations and observations is indeed better. At this satellite stellar mass threshold, the MW has 4 satellites within 150 kpc. We find that 8 out of the 12 simulated hosts match the MW’s profile at 150 kpc for at least 1 snapshot out of 11, together providing a total of 27 matching profiles. Notably, m12i matches for 8 snapshots and Romeo matches for 5 snapshots.

Figure 9 (right) shows the range of simulated profiles that match the MW at 150 kpc compared to the observed MW profile, for satellites with M${}_{*}>10^{6}$ M*⊙. We find that the agreement between the simulations and observations spans the full range of distances in this satellite mass range. The 95 per cent simulation scatter almost completely encompasses the MW observational scatter below 150 kpc, and beyond that the 95 per cent simulation scatter overlaps with the upper half of the observational scatter. The lower 68 per cent (95 per cent) limits on the simulation profile imply that there may be at least 1 (0) more satellite with M${}_{*}>10^{6}$ M⊙* to be discovered within 150-300 kpc of the MW. The median simulation profile indicates that there are on average 2 satellites in this mass range remaining to be discovered around the MW. Compared to the larger number of undiscovered satellites that we predict for the lower mass range, the observations of satellites with M${}_{*}>10^{6}$ M*⊙* appear to be more complete. We find that this strengthens our conclusion that the census of MW satellite galaxies may not be complete down to M${}_{*}>10^{5}$ M*⊙*.

4.9 Incompleteness around M31

M31’s satellite population is complete down to our lowest stellar mass limit (M${}_{*}>10^{5}$ M*⊙) and within 150 kpc of the host given the uniform depth and coverage of PAndAS in this area (McConnachie et al. 2009 and see Section 3 for more discussion). However, outside of the PAndAS footprint, the completeness limit for M31’s satellite galaxies is not clear. We use our simulations as testing grounds to examine effects of this incompleteness on recovering M31’s true radial profile. For simplicity and to match M31’s profile (which has a median value of 27 satellites at 300 kpc), we select hosts from our simulations that have at least 20 (median over time) satellite galaxies within 300 kpc with M${}_{*}>10^{5}$ M⊙: m12m, m12c, m12w, Juliet, and Louise. We perform a mock survey by selecting satellites in 2D projection along 1000 lines of sight. To mimic the PAndAS footprint, we assume that our mock observations are complete down to M${}_{*}=10^{5}$ M⊙* within a projected radial distance of 150 kpc from the host, and within 150-300 kpc we assume two possible estimates of the completeness: M${}_{*}>10^{6}$ M*⊙* and M${}_{*}>10^{7}$ M*⊙*.

Figure 10 shows the results of our mock surveys compared to the true radial profiles for the 5 hosts with M31-like profiles. Comparing our simulated true median profiles (blue dashed) to the recovered profiles (blue solid, with shaded regions showing 68 per cent and 95 per cent scatter), we find that we typically recover 75-80 per cent (median) of our satellites, depending on the completeness mass. Thus, if our estimates of stellar completeness beyond 150 kpc are correct, M31 reasonably has 6-9 undetected satellites with M${}_{*}>10^{5}$ M*⊙* within 150-300 kpc of the host, which we obtain by applying 20-25 per cent incompleteness to M31’s observed profile. It is also worth noting that beyond $\sim$200 kpc, the M31 profile lies above the scatter in the selected simulations. This may indicate that M31 is more massive than our simulated hosts, or that there is something else fundamentally different about M31 compared to our simulations. This result motivates deeper PAndAS-like surveys out to greater distances around M31, which are likely to find several dwarf galaxies, based on our simulations.

5 Summary and Discussion

Using the FIRE-2 baryonic cosmological zoom-in simulations of MW- and M31-mass halos, we study the radial profiles of satellite galaxies with M${}_{*}>10^{5}$ M*⊙. We explore 12 host-satellite systems: 8 isolated MW/M31-like galaxies from the Latte suite + m12z and 4 galaxies in LG-like pairs from the ELVIS on FIRE suite, where the hosts span M${}_{200\rm{m}}=0.9-1.7\times 10^{12}$ M⊙*. To reduce noise in profiles at small distances from satellites momentarily near pericenter, we time-average the simulated radial profiles over $z=0-0.1$ ($\sim$1.3 Gyr). We compare against the 3D profiles measured around the MW and M31 (including observational uncertainties in line-of-sight distance), and against the 2D profiles of MW analogs in the SAGA survey. Our main conclusions are as follows:

• •

The radial distributions of satellite galaxies with M${}_{*}>10^{5}$ M*⊙* within 300 kpc of their host in the FIRE-2 simulations agree well with LG observations. The scatter in the simulations spans the radial profiles of the MW and M31, and the median ratio of simulated-to-observed profiles is typically $\sim$1 for the MW and $\sim$1/2 for M31. Though M31 has a relatively large satellite population, it is still within our simulation scatter.

• •

The radial concentration of the baryonic simulations generally agrees with LG observations, but the MW (and M31 in 2D projection) has a more concentrated shape than the simulations. If we examine simulations with the same number of satellites as the MW at d$<$150 kpc, we find excellent agreement with the MW down to $\sim$50 kpc. Beyond 150 kpc, the matched simulation profiles all lie above the MW profile. We predict 2-10 satellites (at 95%) with M${}_{*}>10^{5}$ M*⊙* to be discovered within 150-300 kpc from the MW.

• •

If we perform mock surveys with the same observational characteristics as PAndAS on our simulations, we recover on average 75-80 per cent of the true satellite population. Based on this, we predict there may be 6-9 undetected satellites around M31 and outside the PAndAS footprint depending on the (uncertain) completeness limit outside of the PAndAS footprint.

• •

2D projected radial profiles of satellite galaxies with M${}_{*}>5\times 10^{6}$ M*⊙* for the simulations also agree with the profiles for the 8 MW analogs from the SAGA survey. The scatter in the simulations spans a majority of SAGA profiles, though 3 SAGA hosts have fewer satellites at large distances ($>$100 kpc).

• •

The agreement we find in radial profiles does not depend strongly on satellite galaxy stellar mass. Thus, even at small distances (<100 kpc) where satellite galaxies are subject to stronger tidal forces from the host’s disk, our simulations resolve the survival and physical destruction of satellites down to our lower stellar mass limit (M${}_{*}>10^{5}$ M*⊙, with typical M${}_{\rm peak}>8\times 10^{8}$ M⊙* or $\sim 2\times 10^{4}$ DM particles).

• •

Simulated hosts with larger stellar masses have fewer satellite galaxies at small distances ($\lesssim 100$ kpc). We interpret this as caused by tidal destruction of satellite galaxies by the gravitational potential the host’s disk. We find a similar correlation with halo mass as well, which we interpret as a manifestation of more massive halos having bigger disks. We note, however, that we examined hosts only over a narrow host halo mass range M${}_{\rm h}=0.9-1.7\times 10^{12}$ M*⊙*.

• •

The variation from host-to-host scatter among the different simulations dominates over time variation at large distances ($\gtrsim$100 kpc), while time variation is the dominant contributor to scatter at small distances ($\lesssim$50 kpc).

• •

KS tests between the radial profiles of the simulations the profiles of the MW and M31 show that most of the simulated profiles are consistent with being drawn from the same underlying distribution as LG observations. However, 4 (1) of the simulations have radial profiles inconsistent with the MW (M31).

• •

Consistent with previous studies, our dark matter-only simulations have many more subhalos at small distances (<100 kpc), and hence larger concentrations in their radial profiles, than their baryonic counterparts. This corroborates the idea that the baryonic simulations have enhanced tidal destruction of satellites due to the additional disk potential present in baryonic hosts. We provide fits to the ratio of baryonic to dark matter-only subhalo counts as a function of distance, which one can use to renormalize existing DMO simulations to include baryonic effects.

We present a thorough comparison of satellite galaxy radial profiles around MW/M31-like galaxies in the FIRE-2 simulations to the LG and to MW analogs from the SAGA survey. Incorporating time dependence of the radial profile over the last 1.3 Gyr in the simulations is key to a robust comparison of the simulations with observations, because the profile at small distances from the host can be highly time-variable. Overall, we find that our simulations are generally representative of current observations. Specifically, combined with the recent results of Garrison-Kimmel et al. (2018a) and Garrison-Kimmel et al. (2019), who analyzed the same FIRE simulation suite, we see broad agreement with the population of ‘classical’ dwarf galaxies (M${}_{*}\gtrsim 10^{5}$ M⊙) in the LG across a wide range of properties: stellar masses, stellar velocity dispersion and dynamical mass profiles, star-formation histories, and now spatial distributions in terms of radial profiles. However, we emphasize that these simulations are not yet able to resolve ultra-faint dwarf galaxies, so it remains unclear how well simulations agree with the profiles of ultra-faints (especially in incorporating incompleteness). The simulations used here also do not yet include the most realistic treatment of cosmic ray physics implemented in the FIRE project, which has effects on the mass of the host and hence the survivability of satellites (Chan et al., 2018, Hopkins et al., in preparation).

Given the correlation with number of satellites at small distances with host stellar mass, we interpret the analogous dearth of subhalos in the baryonic simulations relative to the DMO simulations as primarily from tidal disruption of satellites by the baryonic disk. This agrees with a wealth of previous work that generally finds an excess of DMO subalos near the host relative to the number in baryonic simulations or DMO+analytical disk potential (Taylor & Babul, 2001; Hayashi et al., 2003; Read et al., 2006a, b; Berezinsky et al., 2006; D’Onghia et al., 2010; Peñarrubia et al., 2010; Brooks et al., 2013; Zhu et al., 2016; Errani et al., 2017; Garrison-Kimmel et al., 2017b; Sawala et al., 2017; Kelley et al., 2018).

The shape of the radial profile of satellite galaxies also has significant implications for how other satellite phenomena are measured. For example, the missing satellites problem (e.g. Moore et al., 1999; Klypin et al., 1999b) and the satellite plane problem (e.g. Pawlowski, 2018) are both sensitive to concentration of the radial profile, and the MW’s satellite distribution is often found to be unusually concentrated compared to simulations (e.g. Zentner et al., 2005; Li & Helmi, 2008; Metz et al., 2009; Yniguez et al., 2014). However, controlling for the shape of the profile proves difficult because typical metrics of radial concentration do not necessarily produce the comprehensive description of spatial distribution that is needed to interpret observations. We find that DMO simulations have systematically higher radial concentration than baryonic simulations. Other studies have reached the same conclusion by comparing DMO simulations to baryonic simulations or to DMO simulations with a semianalytic model of galaxy formation (e.g. Kang et al., 2005; Ahmed et al., 2017). This suggests that DMO simulations alone cannot be used to accurately predict the shapes of observed radial profiles which are likely affected by baryonic processes.

We also find that while the radial concentration of the M31 profile agrees with our baryonic simulations, the MW is more concentrated than the baryonic simulations when we compare profile shape with $R_{90}/R_{50}$. The MW is more concentrated than the baryonic simulations (and even most of the DMO simulations) under this metric because 50 per cent of MW satellites are within 110 kpc of the MW, but the simulations only attain this fraction of satellites within $\gtrsim$140 kpc of the host on average. This is similar to what Yniguez et al. (2014) found by comparing the number of satellites within 100 and 400 kpc of their host for LG profiles and DMO simulations: the MW has a more concentrated shape than all of their simulations and M31. If we instead match the number of simulated satellites within 150 kpc of their host to the observed number within 150 kpc of the MW, we find that simulations meeting this criteria unanimously show a larger number of satellites within 150-300 kpc than the MW. We interpret this as potential evidence for incompleteness in the MW’s satellite population at large distances, and our simulations predict there are on average 5 (at least 2) satellites with M${}_{*}>10^{5}$ M*⊙* to be discovered beyond 150 kpc from the MW.

Due to the peculiarity of the MW profile, we also use KS testing to accurately compare our simulations with observations. We find that all 12 of the simulated hosts have at least 10 snapshots matching M31’s profile, and 9 of the hosts have at least 10 snapshots matching the MW’s profile. We will examine the full three-dimensional spatial and dynamical distributions of satellite galaxies in detail and examine the satellite plane problem in our simulations in future work (Samuel et al., in preparation).

The spatial distribution of satellite galaxies correlates with attributes of the host galaxy, both in our simulations and in the LG. Importantly, the correlated host attributes are not limited to the dark matter halo properties of the host, and the spatial distribution of satellites may be most strongly correlated with the host’s baryonic features. DMO simulations are insensitive to the effects of a realistic host galaxy disk, and thus are not sufficient predictors of observed radial profiles at small distances which are the most influenced by the host’s baryonic structure. We have provided a correction to such DMO radial profiles by modeling the depletion of subhalos by the baryonic disk as a function of distance from the host.

Acknowledgements

We thank Marla Geha, Risa Wechsler, and Ethan Nadler for their helpful comments. This research made use of Astropy,444http://www.astropy.org a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018), the IPython package (Pérez & Granger, 2007), NumPy (Van Der Walt et al., 2011), SciPy (Jones et al., 2001), Numba (Lam et al., 2015), and matplotlib, a Python library for publication quality graphics (Hunter, 2007).

JS, AW, and SB were supported by NASA, through ATP grant 80NSSC18K1097 and HST grants GO-14734 and AR-15057 from the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. We performed this work in part at the Aspen Center for Physics, supported by NSF grant PHY-1607611, and at the KITP, supported NSF grant PHY-1748958. Support for SGK and PFH was provided by an Alfred P. Sloan Research Fellowship, NSF grant #1715847 and CAREER grant #1455342, and NASA grants NNX15AT06G, JPL 1589742, 17-ATP17-0214. Support for SRL was provided by NASA through Hubble Fellowship grant #HST-JF2-51395.001-A awarded by STScI. KE was supported by an NSF graduate research fellowship. MBK acknowledges support from NSF grant AST-1517226 and CAREER grant AST-1752913 and from NASA grants NNX17AG29G and HST-AR-13888, HST-AR-13896, HST-AR-14282, HST-AR-14554, HST-AR-15006, HST-GO-12914, and HST-GO-14191 from STScI. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522, by NASA through grant 17-ATP17-0067, and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. JSB was supported by NSF AST-1518291, HST-AR-14282, and HST-AR-13888. We ran simulations using the Extreme Science and Engineering Discovery Environment (XSEDE) supported by NSF grant ACI-1548562, Blue Waters via allocation PRAC NSF.1713353 supported by the NSF, and NASA HEC Program through the NAS Division at Ames Research Center.

Appendix A Resolution test

To examine the dependence of our satellite profiles on numerical resolution, we use the lower-resolution (LR) versions of the Latte simulation suite. We simulated each of the 7 Latte hosts at 8$\times$ lower mass resolution, with baryonic particle masses of m${}_{\rm bary}\sim 5.7\times 10^{4}$ M*⊙* and m${}_{\rm dm}=2.8\times 10^{5}$ M*⊙*. Furthermore, all gravitational force softenings are 2$\times$ larger.

In principle, we could compare satellites at fixed M∗ between LR and high-resolution (HR) simulations. However, we choose to compare the survival of subhalos in the baryonic simulations that are resolved with the same number of DM particles, for three reasons. First, as studied extensively in Section 4.1.4 of Hopkins et al. (2018), the stellar masses of dwarf galaxies resolved with small numbers of star particles are sensitive to numerical convergence; the lowest-mass galaxies resolved in our LR simulations ($\sim 20$ star particles) form systematically $\sim 2\times$ higher M∗ at fixed subhalo Mpeak than in our HR simulations. Thus, comparing satellites at fixed M∗ at our resolution limit mixes the numerical effects of star-formation efficiency and tidal disruption, but comparing satellites at fixed Mpeak isolates the effects of tidal disruption, which is our goal here. Second, because dwarf galaxies at these masses are so DM-dominated, the survivability of a satellite is governed more directly by the number of DM particles in its subhalo than its number of star particles. Finally, most previous works on numerical disruption of satellites (e.g. van den Bosch & Ogiya, 2018) focused on DM-only simulations and how well resolved subhalos are, so using Mpeak makes our tests more comparable to those previous works.

For the lowest-mass galaxies that we examine in this work, M${}_{*}\sim 10^{5}$ M*⊙, we find that, across the HR simulations, they are hosted by subhalos with average M${}_{\rm{peak}}\sim 8\times 10^{8}$ M⊙* ($\sim 2\times 10^{4}$ DM particles). Similarly, the HR satellites with M${}_{*}\sim 10^{6}$ M*⊙* have subhalos of M${}_{\rm{peak}}\sim 2\times 10^{9}$ M*⊙* ($\sim 7\times 10^{4}$ DM particles) and HR satellites with M${}_{*}\sim 10^{7}$ M*⊙* have subhalos of M${}_{\rm{peak}}\sim 10^{10}$ M*⊙* ($\sim 3\times 10^{5}$ DM particles). With 8$\times$ larger particle mass in the LR simulations, a LR subhalo that is resolved as well as the subhalos of our lowest-mass satellites in the HR simulations has M${}_{\rm peak}\sim 6.4\times 10^{9}$ M*⊙. We thus select subhalos in both the LR and HR simulations with M${}_{\rm{peak}}>6.4\times 10^{9}$ M⊙* and compute their radial profiles out to 1000 kpc around each host, averaging over all 67 snapshots at $z=0-0.1$ to improve statistics.

Figure 11 shows a comparison of the LR radial profiles to the HR versions. The top panel compares the cumulative profile while the bottom panel compares the differential (discrete distance bins) profile, which more directly indicates where converge occurs. The solid lines and shaded regions show the median and 68 per cent and 95 per cent host-to-host scatter over the 7 Latte simulations. Within $\sim$30 kpc, the LR simulations show a deficit of $\sim$25 per cent (in the differential) compared with HR. Beyond 30 kpc, the median ratio stays mostly between 0.8 and 1, and the 68 per cent host-to-host scatter is always consistent with 1 for both cumulative and differential profiles. We thus conclude that the radial profiles are well converged (to $\sim 20$ per cent) beyond $\sim$30 kpc, where almost all observed satellites (M${}_{*}>10^{5}$ M*⊙*) of the MW and M31 are.

Moreover, the difference between LR and HR within 30 kpc is exaggerated by the fact that the LR simulations have more massive host galaxies. As with the dwarf galaxies themselves, and as studied in detail in Hopkins et al. (2018), the stellar masses of the host galaxies are sensitive to resolution as well, with LR hosts having on average 1.7 times higher stellar mass than their HR counterparts. Using the results from Figure 6 (left), a host galaxy with 1.7$\times$ higher stellar mass will have $\sim$15 per cent fewer satellites at d$<$50 kpc, even at fixed resolution, which can account for most of the difference between LR and HR at small distances. Based on this, we plot the expected median values for the ratio of LR to HR profiles at 50 and 100 kpc as black points in the upper panel of Fig 11. These profile ratio values are within 10 per cent of unity, showing better agreement between the LR and HR simulations.

Furthermore, given the complicating effects of different host galaxy masses in LR versus HR simulations, we also note similar results from the extensive numerical tests in Hopkins et al. (2018). Specifically, in Section 4.14 and Figure 14 they compared the (differential) number of subhalos versus distance in a DM-only simulation of the same m12i host at the same resolutions that we use here, using a broadly similar subhalo selection (instantaneous M${}_{\rm bound}>10^{8}$ M⊙). Thus, while that convergence test did not include the additional tidal force of the central galaxy or any other baryonic effects, its does provide a cleaner numerical test in the DM-only regime. They found convergence to better than 20 per cent down to $d\approx 50$ kpc, consistent with the results of Figure 11.

We also can use the more rigorous criteria articulated in van den Bosch & Ogiya (2018) to test how well our lowest-mass subhalos are resolved. They stipulate that a subhalo on a circular orbit around a static, spherically symmetric host potential will suffer from numerical noise or disruption if the bound mass fraction falls below either of two limits that depend on mass and force resolution:

[TABLE]

[TABLE]

Where $N_{\rm acc}$ is the number of DM particles in the subhalo at accretion, $c$ is the NFW concentration parameter of the subhalo at accretion, $f(c)=\ln(1+c)-\frac{c}{1+c}$, $r_{\rm s,0}$ is the NFW scale radius at accretion, and $r_{\rm h}$ is the instantaneous half-mass radius of the subhalo. We have verified that nearly all of the lowest mass subhalos considered in our HR simulations ($N_{\rm acc}\sim 2\times 10^{4}$ DM particles) verify these criteria. Note however, that these criteria were generated from idealized DMO simulations that do not account for the disk potential present in our baryonic hosts, which more efficiently (and rapidly) disrupts subhalos that orbit close to the disk (see e.g. Garrison-Kimmel et al., 2017b).

Appendix B Differential radial distribution

We also examined the differentially-binned radial profiles of satellites around the hosts in our simulations. Figure 12 shows these profiles for the simulations and MW/M31 observations, considering all satellites with M${}_{*}>10^{5}$ M*⊙*. The simulation scatter (blue regions) encompasses both the MW and M31 profiles at the 68 per cent level out to about 150 kpc. Past this distance, the MW shows a known lack of satellites between $150-200$ kpc, but is otherwise consistent with the simulations at the 68 per cent level. Potential incompleteness in the MW’s satellite population is explored further in Section 4.8. M31 remains within the simulation scatter at the 95 per cent level from $150-300$ kpc. In general, the simulated differential distributions are a reasonable match to the Local Group.

Appendix C Correlation with host galaxy and halo mass

In Figure 13, we repeat the exercise of section 4.5, but this time plotting the number of satellites at small distances as a function of host halo mass and color coding my host disk mass. This illustrates that both host disk mass and halo mass simultaneously correlate with the number of satellites at small distances from the host. Therefore, the main driver of the negative trends with host mass remains uncertain in our analysis. However, other work that has systematically varied an analytical disk potential at fixed halo mass and found that the disk was the source of a reduction in DMO substructure close to the host when compared to a disk-less host halo (Kelley et al., 2018). Although our results are not necessarily definitive on their own, we conclude that it is not unreasonable for the trends we see in Nsat as a function of host mass to originate from enhanced tidal destruction of satellites due to the host’s baryonic disk.