Dark and luminous satellites of LMC-mass galaxies in the FIRE simulations

TL;DR

This study uses FIRE simulations to analyze the dark and luminous satellite populations of LMC-mass galaxies, comparing simulation results with recent Gaia observations to understand satellite abundance, properties, and associations.

Contribution

It provides new insights into the substructure depletion, satellite luminosity functions, and dark matter content of LMC-mass hosts, aligning simulations with recent observational data.

Findings

LMC-mass hosts deplete ~30% of substructure within 0.2 r200, less than Milky Way hosts.

Simulations produce 5-10 luminous satellites with M* > 10^4 M_sun, matching observed pre-infall companions.

Simulated luminosity functions are steeper than observed, indicating possible observational incompleteness.

Abstract

Within $\Lambda$CDM, dwarf galaxies like the Large Magellanic Cloud (LMC) are expected to host numerous dark matter subhalos, several of which should host faint dwarf companions. Recent Gaia proper motions confirm new members of the LMC-system in addition to the previously known SMC, including two classical dwarf galaxies ($M_\ast$ > $10^5$ M$_{\odot}$; Carina and Fornax) as well as several ultra-faint dwarfs (Car2, Car3, Hor1, and Hyd1). We use the Feedback In Realistic Environments (FIRE) simulations to study the dark and luminous (down to ultrafaint masses, $M_\ast$ ~ $6\times10^3$ M$_{\odot}$) substructure population of isolated LMC-mass hosts ($M_\text{200}$ = $1-3\times10^{11}$ M$_{\odot}$) and place the Gaia + DES results in a cosmological context. By comparing number counts of subhalos in simulations with and without baryons, we find that, within 0.2 $r_\text{200}$, LMC-mass hosts deplete ~30% of their substructure, significantly lower than the ~70% of substructure depleted by Milky Way (MW) mass hosts. For our highest resolution runs ($m_\text{bary}$ = 880 M$_{\odot}$), ~5-10 subhalos form galaxies with $M_\ast$ > $10^4$ M$_{\odot}$, in agreement with the 7 observationally inferred pre-infall LMC companions. However, we find steeper simulated luminosity functions than observed, hinting at observation incompleteness at the faint end. The predicted DM content for classical satellites in FIRE agrees with observed estimates for Carina and Fornax, supporting the case for an LMC association. We predict that tidal stripping within the LMC potential lowers the inner dark matter density of ultra faint companions of the LMC. Thus, in addition to their orbital consistency, the low densities of dwarfs Car2, Hyd1, and Hyd2 reinforce their likelihood of Magellanic association.

Tables

Table 1. Table 1: Properties of the host halo of all simulations analyzed herein for both the N-body dark matter only run ( -dmo ) and hydrodynamic run (no suffix) at z = 0 𝑧 0 z=0 . References for additional details on each simulation are shown in the rightmost column. Our primary sample, the LMC-mass m11 hosts, are shown in the top block, and our reference sample, the MW-mass m12 hosts, are shown below. Resolution refers to the baryonic mass resolution. All quantities are computed from halo catalogs generated by the ROCKSTAR halo finder (calculated using DM particles unless the item explicitly refers to stellar quantities, e.g. M ∗ subscript 𝑀 ∗ M_{\ast} ). The minimum V max subscript 𝑉 max V_{\text{max}} is the median value for subhalos of 195<N particles particles {}_{\text{particles}} <205, and serves as an effective resolution limit for substructures in each simulation. This particle number was chosen in accordance with the Hopkins et al. ( 2018b ) determination of DM convergence radii with the Power et al. ( 2003 ) criterion. The virial radius ( r 200m subscript 𝑟 200m r_{\text{200m}} ) is the radius where the average interior DM density is equal to 200 times the mean matter density in the Universe. The number of subhalos within 0.2 × \times r 200m subscript 𝑟 200m r_{\text{200m}} , 0.4 × \times r 200m subscript 𝑟 200m r_{\text{200m}} , and r 200m subscript 𝑟 200m r_{\text{200m}} of the host are selected as having V max subscript 𝑉 max V_{\text{max}} > 5 km s -1 , this cutoff value is chosen to include only confidently resolved subhalos in all simulations. Note that m11c and m11v do not have DMO versions. In the last column, we include the number of luminous satellites, defined as M ∗ subscript 𝑀 ∗ M_{\ast} > 0, though, practically, the stellar mass of the smallest resolved satellite (and hence total number of resolved satellites) varies with resolution. All M 200m subscript 𝑀 200m M_{\text{200m}} and M ∗ subscript 𝑀 ∗ M_{\ast} for non-DMO m12 simulations were taken from Samuel et al. ( 2019 ) , as they include all mass components in the calculation of M 200m subscript 𝑀 200m M_{\text{200m}} , as well as only include stellar mass associated with the galactic disk (i.e. excluding stellar mass contained in the stellar halo). References: 1 - Chan et al. ( 2018 ) ; 2 - El-Badry et al. ( 2018a ) ; 3 - Lazar et al. (in prep); 4 - Hopkins et al. ( 2018b ) ; 5 - Wetzel et al. ( 2016 ) ; 6 - Samuel et al. ( 2019 ) .

Simu-	Resolution	$M_{\text{200m}}$	$M_{\ast }$	$r_{\text{200m}}$	$V_{\text{max}}$	min. $V_{\text{max}}$	$N_{\text{sub}}$	$N_{\text{sub}}$	$N_{\text{sub}}$	$N_{\text{lum}}$	Ref.
lation	(M⊙)	(M⊙)	(M⊙)	(kpc)	(km s-1)	(km s-1)	< 0.2$r_{\text{200m}}$	< 0.4$r_{\text{200m}}$	< $r_{\text{200m}}$	< $r_{\text{200m}}$
m11c	2100	1.5e11	8.2e8	167.4	81.0	2.5	9	34	118	9	1
m11d	7070	2.8e11	4.1e9	203.9	88.8	4.0	16	76	257	11	2
m11d-dmo	7070	3.4e11		216.5	95.4	4.8	21	95	333		3
m11e	7070	1.5e11	1.4e9	166.0	84.3	4.2	10	31	129	5	2
m11e-dmo	7070	1.8e11		176.2	92.1	4.6	4	34	162		3
m11q	880	1.5e11	3.4e8	168.7	81.0	1.7	17	56	123	8	4
m11q-dmo	880	1.9e11		182.5	91.2	1.8	20	81	192		4
m11v	7070	2.9e11	2.4e9	210.5	84.0	4.0	11	62	257	6	4
m12b	7070	1.4e12	7.3e10	335.1	180.9	4.2	19	125	544	25	5
m12b-dmo	7070	1.4e12		358.8	178.8	4.7	91	343	1116		5
m12c	7070	1.4e12	5.1e10	328.4	156.4	4.0	62	285	895	42	5
m12c-dmo	7070	1.3e12		350.0	154.1	4.4	156	537	1492		5
m12f	7070	1.7e12	6.9e10	354.7	183.4	4.0	34	202	863	31	4,5
m12f-dmo	7070	1.8e12		383.7	176.0	4.6	132	519	1693		4,5
m12i	7070	1.2e12	5.5e10	314.2	161.2	4.0	28	194	681	24	4,5
m12i-dmo	7070	1.2e12		339.3	162.3	4.8	103	431	1225		4,5
m12m	7070	1.6e12	1.0e11	341.6	184.3	4.3	48	259	900	40	4,5
m12m-dmo	7070	1.5e12		364.5	171.1	4.9	116	515	1622		4,5
m12r	7070	1.1e12	1.5e10	304.3	136.8	4.1	53	250	822	28	6
m12r-dmo	7070	1.1e12		322.6	146.0	4.8	95	420	1192		6
m12w	7070	1.1e12	4.8e10	300.5	156.2	4.5	33	181	638	31	6
m12w-dmo	7070	1.2e12		328.3	158.6	4.9	125	448	1167		6

Table 2. Table 2: Average subhalo depletion ( N dmo subscript 𝑁 dmo N_{\text{dmo}} / N hydro subscript 𝑁 hydro N_{\text{hydro}} ) at V max subscript 𝑉 max V_{\text{max}} = 10 km/s for both MW-mass hosts ( M halo ∼ similar-to subscript 𝑀 halo absent M_{\text{halo}}\sim 10 12 M ⊙ ) and LMC-mass hosts ( M halo ∼ similar-to subscript 𝑀 halo absent M_{\text{halo}}\sim 10 11 M ⊙ ). Each column specifies the cutoff distance inside which subhalos are counted. Among LMC-mass hosts, we observe modest depletion ( N dmo subscript 𝑁 dmo N_{\text{dmo}} / N hydro ≈ 1.1 − 1.4 subscript 𝑁 hydro 1.1 1.4 N_{\text{hydro}}\approx 1.1-1.4 ) of subhalos at all radii, while subhalo depletion is a strong function of distance for MW-mass hosts.

	$N_{\text{dmo}}$/$N_{\text{hydro}}$	$N_{\text{dmo}}$/$N_{\text{hydro}}$	$N_{\text{dmo}}$/$N_{\text{hydro}}$
d${}_{\text{sub}}$ <	$r_{\text{200m}}$	0.4 $r_{\text{200m}}$	0.2 $r_{\text{200m}}$
MW-mass	1.54$\pm$0.06	1.98${}_{\text{-0.22}}{}^{}{}_{}{}^{\text{+0.28}}$	3.51${}_{\text{-1.72}}{}^{}{}_{}{}^{\text{+1.82}}$
LMC-mass	1.28${}_{\text{-0.1}}{}^{}{}_{}{}^{\text{+0.12}}$	1.14${}_{\text{-0.23}}{}^{}{}_{}{}^{\text{+0.35}}$	1.39${}_{\text{-0.58}}{}^{}{}_{}{}^{\text{+0.79}}$

Table 3. Table 3: Normalized cartesian components of satellite galaxy orbital angular momenta, as predicted by S17 for the first pericentric passage of the LMC (labelled ‘ t 1p subscript 𝑡 1p t_{\text{1p}} ’) and observed angular momenta as calculated form GAIA proper motions, originally tabulated in Gaia Collaboration et al. ( 2018 ) Table C.4 (labelled ‘obs’). These new measurements show Carina and Fornax as being consistent with a co-infall scenario with the LMC using criteria for Magellanic Cloud system membership as defined in S17: j x < 0 subscript 𝑗 𝑥 0 j_{x}<0 and | j x | > | j y | , | j z | subscript 𝑗 𝑥 subscript 𝑗 𝑦 subscript 𝑗 𝑧 |j_{x}|>|j_{y}|,|j_{z}| .

Name	Time	$j_x$	$j_y$	$j_z$
LMC	$t_{\text{1p}}$	-0.97 $\pm$ 0.03	0.14 $\pm$ 0.07	-0.19 $\pm$ 0.10
	obs	-0.93 $\pm$ 0.06	-0.1 $\pm$ 0.03	-0.36 $\pm$ 0.03
SMC	$t_{\text{1p}}$	-0.92 $\pm$ 0.05	0.04 $\pm$ 0.10	0.35 $\pm$ 0.08
	obs	-0.87 $\pm$ 0.06	-0.4 $\pm$ 0.04	-0.27 $\pm$ 0.04
Carina	$t_{\text{1p}}$	-0.93 $\pm$ 0.12	0.25 $\pm$ 0.07	-0.04 $\pm$ 0.20
	obs	-0.97${}_{-0.14}{}^{}{}_{}{}^{+0.13}$	0.17 $\pm$ 0.04	-0.19 $\pm$ 0.1
Fornax	$t_{\text{1p}}$	-0.92 $\pm$ 0.20	0.19 $\pm$ 0.11	0.20 $\pm$ 0.07
	obs	-0.96 ${}_{-0.22}{}^{}{}_{}{}^{+0.23}$	-0.17${}_{-0.14}{}^{}{}_{}{}^{+0.13}$	0.24 $\pm$ 0.07
Sculptor	$t_{\text{2p}}$	-0.94 $\pm$0.06	-0.00 $\pm$ 0.41	0.04 $\pm$ 0.05
	obs	0.99 $\pm$ 0.01	-0.03 ${}_{-0.07}{}^{}{}_{}{}^{+0.08}$	0.13 $\pm$ 0.01
Sagittarius	$t_{\text{1p}}$	-	-	-
	obs	0.05$\pm$ 0.02	-0.99${}_{-0.08}{}^{}{}_{}{}^{+0.09}$	-0.13${}_{-0.09}{}^{}{}_{}{}^{+0.08}$
Ursa Minor	$t_{\text{1p}}$	-	-	-
	obs	-1.0 $\pm$ 0.08	-0.09 $\pm$ 0.07	-0.02 $\pm$ 0.06
Leo I	$t_{\text{1p}}$	-	-	-
	obs	-0.5 ${}_{-0.37}{}^{}{}_{}{}^{+0.4}$	-0.59${}_{-0.48}{}^{}{}_{}{}^{+0.5}$	-0.64${}_{-0.34}{}^{}{}_{}{}^{+0.36}$
Sextans	$t_{\text{1p}}$	-	-	-
	obs	-0.39 $\pm$ 0.05	-0.58 $\pm$ 0.06	-0.71 $\pm$ 0.06
Leo II	$t_{\text{2p}}$	-0.92 $\pm$ 0.05	0.21 $\pm$ 0.21	-0.28 $\pm$ 0.15
	obs	-0.13${}_{-2.05}{}^{}{}_{}{}^{+2.26}$	0.97${}_{-2.12}{}^{}{}_{}{}^{+2.24}$	0.22 ${}_{-0.81}{}^{}{}_{}{}^{+0.88}$
Bootes I	$t_{\text{1p}}$	-	-	-
	obs	0.63${}_{-0.15}{}^{}{}_{}{}^{+0.16}$	0.71 $\pm$ 0.11	-0.31 $\pm$ 0.08
Draco	$t_{\text{1p}}$	-	-	-
	obs	0.89${}_{-0.09}{}^{}{}_{}{}^{+0.1}$	0.36 $\pm$ 0.06	-0.27 ${}_{-0.09}{}^{}{}_{}{}^{+0.08}$

Table 4. Table 4: Properties of dwarf galaxies near the Milky Way, selected as having either M ∗ subscript 𝑀 ∗ M_{\ast} > 3 × 10 4 absent superscript 10 4 \times 10^{4} M ⊙ , or begin determined a potential LMC satellite by S17, K18, or this work. All units and coordinates are in the galactocentric frame. The columns ‘possible’ and ‘confirmed’ refer to the LMC association criteria established by Sales et al. ( 2017 , ‘S17’) and Kallivayalil et al. ( 2018 , ‘K18’) . Galaxies determined to be consistent with the LMC system by our calculations of angular momenta (listed in Table 3 ) from Gaia DR2 data (Gaia Collaboration et al., 2018 ) are listed as ‘this work.’ Any M ∗ subscript 𝑀 ∗ M_{\ast} with a star marker ( ⋆ ) was calculated from the visual magnitude listed in K18 (except for Antlia 2, whose M V subscript 𝑀 𝑉 M_{V} was obtained from Torrealba et al. 2018 ) assuming a mass-to-light ratio of 2. Any V r subscript 𝑉 𝑟 V_{r} with a dagger ( † ) was converted from its originally tabulated heliocentric value in K18. Numbered superscripts refer to sources. Any row with no superscripts is from McConnachie 2012 . A superscript on the name of the galaxy means all properties came from that source. If a property has no superscript but there are others in its row, the nearest superscript to the left is its source. References: 1 - K18; 2 - S17; 3 - Ferrero et al. 2012 ; 4 - Koposov et al. 2018 ; 5 - Gaia Collaboration et al. 2018 ; 6 - Torrealba et al. 2018

Name	$M_{\ast }$	$l$	$b$	$D_{\text{MW}}$	$D_{\text{LMC}}$	$V_{\text{rad}}$ (MW)	possible	confirmed
	(M⊙)	(deg∘)	(deg∘)	(kpc)	(kpc)	(km s-1)
Antlia 2	4.30e5⋆,6	264.9	11.2	124.2	145.4	49.6
Bootes I	2.90e4	358.55	63.3	68.3	98.7	120.4
Canes Ven. I	2.30e5	74.3	79.8	218	254.1	78
Carina	3.80e5	264.85	-22.5	104	60.8	-13.8		this work
Car2	1.08e4⋆	270.01	-17.1	36.2	25.3	211.4†		K18
Car3	1.56e3⋆	270.01	-16.8	27.8	62.9	51.6†		K18
Cra2	3.26e5⋆	283.82	42.0	116.9	114.5	-79.1†
Draco	2.90e5	79.25	34.3	77.2	125.8	-74.6
Dra2	2.47e3⋆	98.31	42.9	22.32	125.4	-159†	K18
Eri3	1.08e33	274.32	-59.6	87.1	48.2		S17
Fornax	2.00e7	243.845	-67.1	146.1	114.5	-38.2		this work
Hercules	3.70e4	28.7	36.9	126	158.5	145
Hor1	3.92e33	270.92	-54.9	79.3	38.5	-30.4†		K18, S17
Hor2	1.88e3⋆	262.52	-54.1	79	39.6		S17
Hyd1	1.30e4⋆	304.51	-16.5	27.6	28.4	-51.4†		K18
Hyd2	1.42e43	295.62	30.5	125.2	113.5	134.2†	K18
Leo I	5.50e6	228.15	50.1	250.5	263.7	159.6
Leo II	7.40e5	223.75	68.6	231.2	255.2	18.6
LMC	1.50e9	290.25	-32.5	51.8	0.0	43.4
Phx2	2.25e33	323.32	-60.2	80.2	51.9		K18, S17
Ret32	2.00e3	273.9	-45.7	92	45.0		S17
Sag2	2.06e4⋆	18.92	-22.9	60.1	49.6
Sculptor	2.30e6	318.65	-80.2	86.7	65.6	77.0
Sextans I	4.40e5	250.3	43.7	83.8	94.1	37.7
SMC	4.60e8	309.45	-41.8	66.2	23.2	-4.2		K18, S17
Tuc42	2.20e3	313.3	-55.3	45.5	26.6		S17
Tuc52	5.00e2	316.3	-51.9	51.9	28.2		S17
Ursa Minor	2.90e5	96.55	45.6	75	125.4	-72.4

Full text

††pubyear: 2019††pagerange: Dark and luminous satellites of LMC-mass galaxies in the FIRE simulations–4

1 Introduction

The $\Lambda$CDM structure formation scenario predicts a nested hierarchy of dark matter (DM) halos, subhalos, and sub-subhalos at all mass scales from galaxy clusters to well below the molecular cooling limit of 106 M*⊙* (Press & Schechter, 1974; White & Rees, 1978; Tegmark et al., 1997; Springel et al., 2008). The abundance of substructure is thought to be nearly scale-free and self-similar. That is, the mass function of subhalos takes a universal form when normalized to the mass of the host (Gao et al., 2004; Kravtsov et al., 2004; Giocoli et al., 2008; Yang et al., 2011). This means that objects from the most massive halos at the centers of giant clusters to isolated dwarf galaxies should host a similar distribution of DM substructure when normalized properly111Though, disruption of subhalos and galaxies due to baryons complicates this relation, as described later. These halos and subhalos act as hosts of galaxy formation, providing potential wells in which gas can accumulate and condense into star forming regions. The fact that the relation between stellar mass $M_{*}$ and halo mass is a near power-law below Milky Way (MW) masses (Yang et al., 2003; Behroozi et al., 2010; Guo et al., 2010; Moster et al., 2013; Wright et al., 2017), together with the near invariance of subhalo abundance, means that the number of satellite galaxies normalized to the stellar mass of the central is expected to be independent of host halo mass for log${}_{\text{10}}$($M_{\ast}^{\text{host}}$ / M*⊙*) = 7.25 $-$ 11.75 (Sales et al., 2013).

The local group, consisting of the Milky Way (MW), M31, and their numerous satellite galaxies, offers an ideal testing ground for these ideas. For instance, the Large Magellanic Cloud (LMC), the largest satellite of the MW, has long been speculated to host satellites of its own (Lynden-Bell & Lynden-Bell, 1995; D’Onghia & Lake, 2008; Sales et al., 2011). Recent observational missions including DES, SMASH, PAN-STARRS, ATLAS, and Gaia have revealed numerous dwarf galaxies in the vicinity of the Magellanic system (Bechtol et al., 2015; Drlica-Wagner et al., 2015; Kim & Jerjen, 2015; Koposov et al., 2015; Torrealba et al., 2016; Laevens et al., 2015; Gaia Collaboration et al., 2018), greatly expanding the population of potential LMC-satellites. Based on the orbital properties of these dwarfs, recent works (Sales et al. 2017; Kallivayalil et al. 2018, hereafter S17 and K18) have discovered several possible associations of dwarf galaxies to the LMC (see also Deason et al., 2015; Jethwa et al., 2016; Shao et al., 2018a). According to S17 & K18, there are currently five firm candidates to LMC-system membership: the Small Magellanic Cloud (SMC), Car2, Car3, Hor1, Hyd1, as well as eight promising possible associations awaiting additional proper motion measurements: Dra2, Eri3, Hor2, Hyd2, Phx2, Ret3, Tuc4, and Tuc5. Pardy et al. (2019) have also suggested the possibility of LMC association for the classical dwarf galaxies Carina and Fornax. We expand on this claim in Section 4. Each new measurement marks a step closer to a complete picture of the pre-infall satellite population of the LMC, which will greatly inform our theories of galaxy formation.

While it is known that satellite galaxies trace substructure, it is unclear from a theoretical perspective exactly how the population of dark subhalos is mapped to a population of luminous companions. Cosmological N-Body (dark matter only, ‘DMO’) simulations (e.g. Navarro et al. 1996; Jenkins et al. 1998; Springel et al. 2005, 2008) and hydrodynamical baryonic simulations (e.g. Abadi et al. 2003; Governato et al. 2004; Sales et al. 2007; Brooks et al. 2013; Sawala et al. 2016; Wetzel et al. 2016) have been instrumental in making predictions for both central galaxies and smaller-scale subhalo and satellite mass functions, as well as in revealing discrepancies between theoretical predictions of dark matter structure and the observed stellar structure that follows it.

Tension in predictions made by numerical simulations with the observed population of galaxies include the ‘missing satellites’ problem (Klypin et al., 1999; Moore et al., 1999) and ‘too big to fail’ (Boylan-Kolchin et al. 2011b; Tollerud et al. 2014; Kirby et al. 2014; though see also Read et al. 2006; Okamoto & Frenk 2009), which address this mapping from dark to luminous substructure. Encouragingly, several solutions to such problems have been proposed, mostly invoking a combination of reionization heating, observational incompleteness, and the addition of realistic feedback modeling (Bullock et al., 2000; Somerville, 2002; Wetzel et al., 2016; Sawala et al., 2017; Wang et al., 2015; Garrison-Kimmel et al., 2018; Simpson et al., 2018). Furthermore, tides within the host potential affect subhalos. In particular, the addition of the gravitational potential of the central baryonic disk in MW-mass halos has been shown to suppress the presence and survival of dark matter subhalos compared to DMO runs (D’Onghia et al., 2010; Garrison-Kimmel et al., 2017b; Graus et al., 2018). With a smaller population of subhalos predicted, the number of eligible sites for satellite dwarf galaxies to form is also reduced, highlighting the need for a better understanding of the dark-to-luminous mass mapping predicted within $\Lambda$CDM. Because the central galaxy is the source of this effect, the actual efficacy of this effect depends on the $M_{\ast}$ - $M_{\text{halo}}$ ratio. Galaxies near the size of the MW have the highest $M_{\ast}$ - $M_{\text{halo}}$ ratio, while galaxies of both higher and lower stellar mass are observed to be hosted by fractionally larger halos, a trend which motivates exploring the magnitude of subhalo depletion by the central galaxy at different host mass scales.

Many of the above studies do not resolve satellite galaxies down to the scale of ultrafaints (UFs; M_{\ast}$$\lesssim 105 M*⊙*) - presumably, the most abundant class of galaxies in the universe. Since the satellite population of the LMC can be considered a scaled-down version of that of the MW, it is likely to be dominated by UFs, as the current (incomplete) population of LMC satellites seems to suggest. Due to their shallower gravitational potential, the effects of reionization at this scale of galaxy formation are believed to be stronger than for more massive galaxies like the MW or classical dwarf spheroidals (though feedback effects are weaker since they form proportionally fewer stars). Thus, exploring the ultrafaint population of the LMC (and of LMC-mass hosts in general) is critically important to push the limits of our knowledge of dark-to-luminous substructure mapping into much fainter scales than currently known.

A possible avenue to overcome numerical resolution limitations in cosmological simulations is to focus on the formation of isolated dwarf galaxies. The lower masses and smaller sizes expected for dwarf halos translate into smaller mass per particle and smaller gravitational softening at a fixed number of particles compared to a more massive halo. For example, Wheeler et al. (2015); Wheeler et al. (2018) were able to study an extremely high resolution population of ultrafaint satellites in simulations of dwarf host halos in the scale $M_{\text{200m}}$ $\sim 10^{10}$ M*⊙. On average, $1-2$ ultrafaints were expected above $M_{\ast}$ $\geq 3\times 10^{3}$ M⊙, making important predictions for the numbers and distribution of ultrafaint dwarfs around dwarf galaxies in the field. Promisingly, larger numbers of ultrafaints shall be expected for more massive dwarf hosts. The inferred virial mass of the LMC, $\sim 10^{11}$ M⊙* at the large end of the dwarf galaxy scale, promises to provide numerous substructures with which subhalo abundance can be studied. In addition, its recent infall to the MW ($1-3$ Gyr ago, Kallivayalil et al. 2013), means its satellite population will remain relatively undisturbed by the tidal field of the MW (Sales et al., 2011; Boylan-Kolchin et al., 2011a; Deason et al., 2015) and may offer an observational avenue to reconstruct its pre-infall satellite companions.

Previous theoretical works on local group satellite galaxies have predicted that $\sim$30% of M_{\ast}$$\lesssim 105 M*⊙* satellites of MW-mass hosts fell in as a satellite of a more massive galaxy (Wetzel et al., 2015), and that such groups of dwarf galaxies typically disperse in phase space about 5 Gyr after infall to a MW/M31 system (Deason et al., 2015). In addition, Dooley et al. (2017) predicts $\sim$ 8 large ultrafaint dwarf satellites of the LMC using reionization and abundance matching models with the Caterpillar simulations (Griffen et al., 2016). However, such predictions are all based on DMO simulations missing important phenomena expected once the baryons are self-consistently taken into account, for example, SNe feedback, gas hydrodynamics affecting morphology, and tidal effects from the central baryonic disk.

In this work, we present the first analysis of satellite galaxies down to the ultra-faint mass scale in simulations of LMC-mass hosts, using the Feedback In Realistic Environments (FIRE) project. We investigate the luminous and dark substructure of these hosts in hopes of gleaning insight into the real history of the LMC-system. How many satellites, dark and luminous, did the LMC bring with it as it fell into the MW system? What is the mass distribution of these satellites? How have they been shaped by co-evolution with the LMC, and does this differ from satellites/substructures of the MW? Are there ways of constraining membership beyond orbital phase space similarity with the LMC?

Previous works using the FIRE simulations have shown promising results, making them a compelling laboratory to investigate these questions. Recent relevant work includes the observational agreement in classical satellites of MW-mass hosts (Wetzel et al., 2016; Garrison-Kimmel et al., 2018), the star formation histories of satellite galaxies in various local-group inspired environemnts (Garrison-Kimmel et al., 2019), radial profiles of classical satellites (Samuel et al., 2019), the origin and evolution of the mass-metallicity relation (Ma et al., 2016), chemical abundance distribution of dwarf galaxies (Escala et al., 2018), as well as gas kinematics and morphologies in dwarf galaxies (El-Badry et al., 2018a, b).

This paper is organized as follows: In Section 2 we descibe the FIRE simulations, incuding the feedback model, the zoom-in technique, halo finding, and resolution. We also describe our sample of LMC-mass zoom-in runs. In Section 3 we describe the effect of the central galaxy population of the abundance of subhalos as compared to DMO runs. In Section 4 we expand on the analysis of S17, K18, and Pardy et al. (2019) by analyzing Gaia orbital angular momenta of satellites in Gaia Collaboration et al. (2018). We compare the updated observationally inferred pre-infall mass function of the LMC to the satellite populations in our FIRE runs. We also investigate the structural kinematics of both observed and simulated dwarf satellites, with a focus on ultrafaint dSphs, and consideration of infall time and tidal stripping. We present a summary of this work and concluding remarks in Section 5.

2 Simulations

We analyze a sample of five cosmological zoom-in simulations of LMC-mass host galaxy systems from the Feedback In Realistic Environments project222http://fire.northwestern.edu (FIRE). These runs, implemented in the updated FIRE-2 scheme (Hopkins et al., 2018b), use the fully conservative cosmological hydrodynamic code GIZMO333http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html (Hopkins, 2015), a multi-method gravity plus hydrodynamics code, in its meshless finite-mass (“MFM”) mode. This is a mesh-free Lagrangian finite-volume Godunov method which automatically provides adaptive spatial resolution while maintaining conservation of mass, energy, and momentum, and excellent shock-capturing and conservation of angular momentum, capturing advantages of both smoothed-particle hydrodynamics (SPH) and Eulerian adaptive mesh refinement (AMR) schemes (for extensive tests, see Hopkins 2015). Gravity is solved with an improved version of the Tree-PM solver from GADGET-3 (Springel, 2005), with fully-adaptive (and fully-conservative) gravitational force softenings for gas (so hydrodynamic and force softenings are always self-consistently matched), following Price & Monaghan (2007).

FIRE-2 implements a variety of methods for cooling, star formation, and stellar feedback processes. Heating and cooling rates are calculated across 10 - 1010K, including CLOUDY ionization states for free-free, photoionization & recombination, Compton scattering, photoelectric, metal-line, molecular, fine structure, dust collisional, uniform cosmic ray heating, from a spatially uniform, redshift-dependent UV background (Faucher-Giguère et al., 2009). Local self-shielding is accounted for using a Sobolev approximation. Stars are formed in accordance with Hopkins et al. (2013), requiring gas to be locally self-gravitating, self-shielding (following Krumholz & Gnedin, 2011), Jeans unstable, and with density $n_{H}>n_{\text{crit}}=1000$ cm*-3*, with all conditions met. Global star formation efficiency is naturally self-regulated by feedback processes, with good observational agreement (Orr et al., 2018). All newly formed star particles inherit mass and metallicity from their progenitor gas particles.

Stellar feedback quantities are tabulated from the stellar population model STARBUST99 (Leitherer et al., 1999), assuming a Kroupa (2001) IMF, including supernova Type Ia, II, and stellar winds, as detailed in Hopkins et al. (2018b) and Hopkins et al. (2018a). Radiative feedback is modelled with the Locally Extincted Background Radiation in Optically-thin Networks (LEBRON) algorithm (Hopkins et al., 2018b), accounting for absorbed photon momentum, photo-ionization, and photo-electric heating.

The zoom-in technique (Katz & White, 1993; Oñorbe et al., 2015) is implemented by first simulating a large, low-resolution cosmological box in which a convex Lagrangian region (at initial redshift) is then defined. The region contains all particles within $\sim$5$r_{\text{vir}}$ at $z=0$ and not containing a halo of similar mass to that of the primary, and is then reinitialized with higher resolution. This process is repeated to refine the Lagrangian region until the intended resolution is reached, at which point it is re-simulated with dark matter, gas, and star particles, buffered by a region of lower resolution dark matter particles. Initial conditions444http://www.tapir.caltech.edu/~phopkins/publicICs/ are generated with the MUSIC code (Hahn & Abel, 2011), which implements a second-order Lagrangian perturbation theory to $z=99$.

Structure and substructure are identified using an updated version of the ROCKSTAR halo finder (Behroozi et al., 2013a), which implements 6+1 dimensional phase-space analysis to determine the particles that are gravitationally bound and assign them to (sub)halos. ROCKSTAR assigns mass to halos and subhalos using spherical overdensity calculations relative to a specified threshold, such as the critical density or the average matter density of the universe. Such quantities necessarily become ambiguous and unreliable when quantifying masses of substructure, because the mass density of such embedded objects is by definition above the chosen threshold. To circumvent this ambiguity, we will mostly refer to the maximum circular velocity as an analog for subhalo mass, because it is a more robust and well-defined metric. For numerical stability, we run ROCKSTAR using only dark-matter particles, and we assign star particles in post-processing using an iterative procedure. We first select star particles within a DM halo out to 0.8 $r_{\text{halo}}$ and with velocities less than 2 $V_{\text{max}}$ of the (sub)halo’s center of mass velocity, and we iteratively compute stellar position and velocity (making sure the star particles and halo are coincident) until the total mass of star particles, $M_{\ast}$, converges to <1%. See Samuel et al. (2019) for a more detailed description.

2.1 Resolution Convergence

All simulations analyzed herein contain no low-resolution particles within $\sim$ 3 $r_{\text{200m}}$ from the host center. However, the 7070 M*⊙* resolution runs naturally use fewer particles to represent objects at a given mass than do the higher resolution runs, leading to convergence issues. For this reason, we have determined a minimum value for $V_{\text{max}}$, above which subhalos have $\gtrsim$200 DM particles, as listed in Table 1. We make this cut for all subhalos analyzed herein, and all galaxies naturally fall above this threshold. Using two different resolution runs of m11q (7070 M*⊙* & 880 M*⊙*), we found that the subhalo populations deviate around V_{\text{max}}$$\sim 4 km/s, which is below our cutoff in $V_{\text{max}}$. This is shown in Figure 8.

2.2 Our Sample: LMC-like centrals and dark matter cores

Our centrals are selected in the halo mass range $M_{\text{200m}}$ = 1.5 $-$ 3.4$\times 10^{11}$ M*⊙, where $M_{\text{200m}}$ refers to the mass measured within $r_{\text{200m}}$, defined as the radius at which the mean interior halo density equals 200 times the average matter density of the universe. The stellar masses of the centrals at $z=0$ are in the range $M_{\ast}$ = 0.34 $-$ 4.1$\times 10^{9}$ M⊙, compared to the 1.5$\times 10^{9}$ M⊙* of the LMC (McConnachie, 2012). These and other properties of our centrals are listed in Table 1. Our choice of ‘200m’ (as opposed to ‘200c’ which uses the critical density of the universe, $\rho_{\text{crit}}$) is motivated by its closer physical proximity to the ‘splashback’ radius, in which all (dark) matter that has passed through the core of the halo is enclosed (see Wetzel & Nagai, 2015, section 2, for a detailed explanation).

We analyze the highest resolution runs available for each system. There are three “low”-resolution ($m_{\text{bary}}=$ 7070 M*⊙) runs: m11c, m11d, m11e; one medium-resolution run ($m_{\text{bary}}=$ 2100 M⊙): m11c; and one high-resolution run ($m_{\text{bary}}=$ 880 M⊙*): m11q. For the first part of this analysis, we consider m11d, m11e, and m11q to be our primary focus because they have counterpart DMO simulations, allowing a comparative analysis of subhalo populations between DMO and baryonic versions of the same system. We examine all runs in Section 4, where we investigate the population of luminous companions around each host. All simulations were run with the core the FIRE-2 hydrodynamics and feedback models, while m11q, m11d, and m11e implemented an additional sub-grid model for the turbulent diffusion of metals in gas. By comparing m11q with and without this model, we have observed no difference in the mass and abundance of satellite galaxies, in agreement with Hopkins et al. (2018b) and Su et al. (2017), which quantified that including (or not including) metal diffusion physics did not change galaxy-wide properties. Phenomena sensitive to the physics of metal diffusion (such as metallicity gradients) are beyond the scope of this paper, so we proceed considering all runs equivalently.

To showcase our sample, the top panels in Fig. 1 show a visualization of the baryons (left) and the dark matter (right) within a region slightly larger than the virial radius of our highest resolution halo m11q. The dark matter component shows a large number of subhalos as expected within $\Lambda$CDM. On the left, stars (yellow) and gas (blue) have collapsed at the center of this LMC-mass halo to form the central dwarf, with $M_{\ast}$ $\sim$ 3$\times 10^{8}$ M*⊙*. However, several other satellite subhalos have also formed stars giving rise to the population of dwarf satellites. We highlight them with starred symbols and annotate their corresponding stellar masses. The dark-to-luminous mapping at these low masses is complex, with only a few subhalos forming stars and the rest remaining dark companions. Circles on the left panel indicate the fifteen subhalos of highest $V_{\text{max}}$ that have remained dark. The bottom row shows stars and gas for the rest of our sample and nicely illustrates the variations expected on the luminous companions of the LMC-mass halos.

Table 1 lists relevant properties of the LMC-analogs analyzed herein. Also listed, for comparison, is the Latte suite of isolated MW-mass hosts (see, e.g., Wetzel et al., 2016; Garrison-Kimmel et al., 2017b, hereafter GK17). See section 2.3 for more details on these simulations. The selection of all LMC-like hosts was blind with respect to satellite and subhalo population, aside from m11v, which was selected to have an ongoing merger at $z=0$, as can be seen by its two companions, each with M_{\ast}$$\gtrsim $M_{\ast}$(SMC) $\sim$ 4$\times 10^{8}$ M*⊙*. All were selected to be isolated from larger halos within $\sim$5 $r_{\text{200m}}$.

At $z=0$, our LMC-mass centrals all show a cored dark matter density profile (see left panel of Fig. 2) which contrasts with the denser and cuspy nature of the dark matter only runs (shown for halos m11q, m11d and m11e in dashed lines). This is better indicated in the bottom panel showing the ratio between the dark matter density in the dark matter only (DMO) run to the baryonic run, $\rho_{\text{dmo}}$ / $\rho_{\text{hydro}}$, as a function of radius. Our LMC-mass centrals show $3-10$ times lower densities in the inner regions when the effect of baryons is included, an effect that is much larger than that observed in more massive $M_{\text{200m}}$ $\sim 10^{12}$ M*⊙* hosts run with similar feedback555m11d and m11e are in the final stages of a merger at $z=0$, making their cores seem larger due to the poorly defined center of each halo. (dot-dashed purple line showing results from Latte, e.g. Wetzel et al., 2016, Lazar et al., in prep).

The formation of cores in our LMC-like galaxies responds to the rapid removal of self-gravitating gas at the center of the halo due mainly to supernova explosions, which changes temporarily the potential at the center of the halos and leads to the “heat up” of the orbits of dark matter particles (Navarro et al., 1996; Peñarrubia et al., 2010; Pontzen & Governato, 2012; Garrison-Kimmel et al., 2014). This is consistent with the starbursty nature of star formation reported in FIRE simulations (Oñorbe et al., 2015; Chan et al., 2015; El-Badry et al., 2016; Fitts et al., 2017) and also coincides with the regime where core formation due to the effect of baryons is expected to be maximal (Di Cintio et al., 2014; Chan et al., 2015; Tollet et al., 2016). However, El-Badry et al. (2017) showed that the degree of coring is subject to change on relatively short timescales as gas outflows and inflows change. This mass rearrangement also affects the circular velocity profiles, as shown by the right-hand panel of Fig. 2. Considering the effects of baryons, the maximum circular velocities of our halos are on the high end of agreement with that measured from HST proper motions of stars in the LMC (indicated by the orange shaded region van der Marel & Kallivayalil, 2014).

2.3 Our Sample: MW-mass centrals

To highlight the dependence on host mass of many phenomena explored herein, we compare to the Latte suite of seven zoom-in simulations of MW-mass hosts (m12b, m12c, m12f, m12i, m12m, m12r, and m12w; introduced in Wetzel et al. 2016) with both baryonic and dark matter only runs, all with a baryonic mass resolution of 7070 M*⊙. All MW-mass hosts were blindly selected, with the exception of m12r and m12w, which were chosen to have an LMC-mass satellite at $z=0$ in DMO (although, the LMC companion does not necessarily survive to $z=0$ in the baryonic runs, see Samuel et al., 2019). These simulations confidently resolve subhalos down to $V_{\text{max}}$ $\sim$ 4 km/s, and luminous satellite galaxies down to $M_{\ast}$ $\sim$ 1$\times 10^{5}$ M⊙* (galaxies with $M_{\ast}$ < 1$\times 10^{5}$ M*⊙* are mostly found in the higher resolution runs). Halo masses range from (1.1 - 1.7)$\times 10^{12}$ M*⊙, while stellar masses range from (1.5 - 10)$\times 10^{10}$ M⊙*. The mean stellar mass to halo mass ratio is $\sim$ 6$\times 10^{-2}$, compared to $\sim$ 8$\times 10^{-3}$ for the LMC-mass hosts. As shown in Fig. 2, the central densities of the MW-mass hosts are much greater than those of our LMC-mass hosts, likely due to the effects of baryonic contraction (Chan et al., 2015). This difference in host density is highly relevant to the discussion on subhalo abundances in the following section.

Any references to the names of various simulations are references specifically to the host/central galaxy of that simulation. The term ‘companion’ refers to a galaxy of nonzero stellar mass within $r_{\text{200m}}$ of its host, while ‘subhalo’ refers to the self-bound dark matter content of any object - luminous or dark - also within $r_{\text{200m}}$ of its host. ‘Satellite’ is a collective term for either. As it is a common phrase, we also use the term ‘virial radius’ to refer to $r_{\text{200m}}$, and ‘virial mass’ for $M_{\text{200m}}$. Note that we do not use the Bryan & Norman (1998) formulation for virial quantities.

3 Suppression of Dark matter subhalos in LMC-like hosts

The number of subhalos above a given $V_{\text{max}}$, referred to as the $V_{\text{max}}$ function, is a useful metric to evaluate the abundance and scale of substructure hosted by a central halo. In the last decade, it has been found that the $V_{\text{max}}$ function in MW-mass halos may be significantly suppressed when considering the increased tidal disruption due to the effect of the baryons in the central disk (D’Onghia et al., 2010; Kelley et al., 2019, GK17). The lower number of surviving subhalos results in a different prediction of the expected number of dwarfs around the MW and in particular for the inner regions where the suppression is maximal. This effect is critical for accurately predicting the radial distribution of satellites around the MW and M31 (Samuel et al., 2019). These and many other recent studies have primarily been concerned with massive host halos such as the MW. In this section we extend the scope of such inquiry to LMC-mass centrals, with a focus on the suppressing effect of the central baryonic galaxy, by comparing subhalo populations in dark matter only (DMO) simulations to subhalo populations of the same centrals simulated with hydrodynamics & baryonic physics.

Figure 3 shows the time-averaged subhalo $V_{\text{max}}$ function of three LMC-mass hosts (orange - m11q; blue - m11d; green - m11e) as a cumulative count of subhalos at a given $V_{\text{max}}$ and within different radial cuts from $r_{\text{200m}}$ to 0.2 $r_{\text{200m}}$ (left to right). Subhalos in each run are plotted from the highest $V_{\text{max}}$ present to the minimum converged $V_{\text{max}}$, as listed in Table 1. Time-averaging was computed over the most recent $\sim$1.3 Gyr, or $z\lesssim 0.1$ by sampling the $V_{\text{max}}$ function of every host at each successive snapshot (66 in total), then computing the mean number of subhalos at each $V_{\text{max}}$, including the $\rm 1\sigma$ deviation from the mean (shaded regions). Shading for the average depletion among the Latte suite of MW-mass simulations (bottom panels, purple dashed-dot lines) was calculated by adding the 1$\rm\sigma$ deviations in quadrature, for each run at each $V_{\text{max}}$. Within the virial radius of a $M_{\text{200m}}$ $\sim 10^{11}$ M*⊙* halo, on average one can expect order $10$ subhalos with $V_{\text{max}}$ $\geq 10$ km/s and order $10^{2}$ above $5$ km/s. We have explicitly checked that the radial distribution of those subhalos follows closely that of the dark matter in the host.

We additionally include the DMO subhalo populations of all hosts to illustrate how the additional baryonic potential and feedback effects change the distribution and total number of subhalos. This difference is quantified in the bottom panel where the number of subhalos at a given $V_{\text{max}}$ in DMO is divided by that number in Hydro. A number greater than one represents ‘suppression’ or ‘depletion’ of subhalos in the Hydro run, while a number less than one means there are more subhalos at that $V_{\text{max}}$ in Hydro than in DMO, or an ‘enhancement’ of substructure.

A close inspection of the bottom row in Fig. 3 shows that, on average, all subhalos with $V_{\text{max}}$ < 15 km/s are only slightly suppressed (by a factor of $\sim$ 1.3) in hydro runs compared to the DMO version, an effect that increases only mildly when looking into the inner regions (middle and right panels). While time-averaging is indeed implemented to mitigate the discreteness of the sample, there are not many (<10) satellites with $V_{\text{max}}$ > 15 km/s around LMC-like hosts and Poisson fluctuation may dominate. Take, for example, halo m11e (green). Although the middle panel seems to suggest an increase in the number of subhalos in the hydro run with respect to the DMO, the effect seems localized and it disappears for the inner regions ($r$ < 0.2 $r_{\text{200m}}$, right panel). We interpret this as a local fluctuation that results from poor numbers statistics: only 6 subhalos exist with $V_{\text{max}}$ = 10 km/s. We see that the overall trend in our LMC-mass hosts is a mild to negligible suppression of their subhalo populations by their central galaxies.

Also included in the bottom panel is the average suppression of subhalos for all MW-mass Latte hosts within the same fraction of their virial radii as our LMC-mass halos (see also GK17; Samuel et al., 2019). In all panels we find that the level of suppression in our LMC-mass halos is significantly smaller than in MW-mass halos. This is most evident in the inner regions $r$ < 0.2 $r_{\text{200m}}$, where subhalos are depleted by a factor $\sim 3-6$ in the m12’s compared to a maximum of $\lesssim$ 1.5 in our best resolved (though least cored) halo m11q. These results are mostly independent of the simulation method used. For example, using DMO runs with an added analytic disk D’Onghia et al. (2010) show that MW-mass environments suppress subhalos by a factor $\sim 3$ at $10^{7}$ M*⊙* (V_{\text{max}}$$=4-5 km/s) within 30 kpc of the center, compared to at most a factor $2$ found in our LMC-mass environments. The lower impact of baryons on the number of subhalos for LMC-mass hosts is consistent with previous arguments. Using hydrodynamical simulations and comparing to an analytic disk potential, GK17 showed that the primary mechanism for suppression of substructure is the enhanced tidal stripping of subhalos by the gravitational potential of the central galaxy. This is due to two effects in the scale of MW-mass hosts. First, the extra component added to the gravitational forces by the presence of the disk and, second, by the dark matter halo itself becoming more concentrated in the presence of the disk leading to a steeper gravitational potential well.

We note that we find no correlation between substructure depletion and resolution. For instance, Fig. 3 shows that “normal” resolution (m11d and m11e) can show either more or less substructure depletion than the high resolution run m11q. There is therefore not an obvious systematic effect with resolution. Instead, substructure depletion may depend on the specifics of the central baryonic mass and accretion history of each halo. This is consistent with the lack of dependence in the resolution tests for $N_{\text{dmo}}$/$N_{\text{hydro}}$ in MW-mass hosts as examined in GK17 and Samuel et al. (2019).

In the case of hosts within the $M_{\text{200m}}$ $=10^{11}$ M*⊙* regime, as analyzed here, the fraction of mass in the disk is much smaller than in MW-mass objects ($\sim$ 1% average for the LMC-analogs compared to $\sim$ 6% for the MW-analogs) partially explaining the lower suppression of substructure. Additionally, the dark matter halos of LMC-hosts in Hydro are predicted to be less dense compared to their DMO counterparts (see Fig. 2) due to the effect of stellar feedback; opposite to the trend found in MW-mass hosts, which are dominated by baryonic contraction. We conclude that substructure depletion is less significant in dwarfs centrals than expected in $\sim L_{*}$ galaxies and that hundreds of dark matter subhalos with $V_{\text{max}}$ $>5$ km/s are expected to be orbiting around field dwarf galaxies with mass comparable to the LMC.

3.1 Implications for Observational Subhalo Searches

Constraining the degree of subhalo suppression in regions near galaxies has strong implications on the viability of observational subhalo searches that are based on gaps in cold stellar streams. As noted by GK17, if all subhalos within 20 kpc of the MW are suppressed by the baryonic disk, it is less likely for surveys like Palomar-5 and GD-1 (Carlberg et al., 2012; Koposov et al., 2010) to detect any interaction with cold dark matter substructure (see also Chapman et al. (in prep), who detected a non-trivial infall rate of subhalos into the inner regions of FIRE MW-mass hosts, but that those subhalos are quickly destroyed). By examining the $z=0$ radial distribution of substructures with $V_{\text{max}}$ > 10 km/s as seen in Figure 10, we find that, on average, the MW-mass halos host no substructure within 30kpc (though individual halos may occasionally have a few subhalos in the range 20 kpc < d${}_{\text{sub}}$ < 30 kpc). On average, the LMC-mass halos are able to host substructure down to $\sim$20 kpc, while retaining roughly a factor of two more subhalos than MW-mass hosts up to approximately 45kpc. The total number of subhalos hosted by MW-mass halos only surpasses that of LMC-mass halos at $\sim$60-70 kpc. In addition, the overall reduced effect of host-subhalo tidal interactions make LMC-mass hosts somewhat cleaner systems to study substructure than their more strongly interacting MW-mass counterparts. We therefore argue that the study of cold streams around galaxies with $M_{\ast}$ $\sim 10^{9}$ M*⊙* may represent a more promising avenue to detect gaps associated to these dark subhalos.

Of course, the mass range of the substructure is also important. Previous work quantifying the sensitivity of cold stellar structures to disturbances by subhalos (e.g. Yoon et al. 2011) suggest that streams such as Palomar-5 are sensitive to DM substructures in the mass range 107 to 109 M*⊙. The range of median subhalo masses with $V_{\text{max}}$ present in the right panel of Figure 3 are 1.1$\times 10^{7}$ to 1.6$\times 10^{8}$ M⊙*, which, while on the lower end of the quoted sensitivity range, is still within the bounds. This agreement supports our prediction that stellar streams near Magellanic-like systems may have more success that those in the vicinity of the MW-mass objects.

Indirectly, gravitational lensing searches are another avenue to probe dark matter substructure. Line ratio anomalies (Chiba, 2002; Metcalf & Zhao, 2002) are strongly influenced by the degree of substructure predicted in the lens system. Novel methods using adaptive optics integral field spectroscopy to measure deviations in quasar narrow line emission (Nierenberg et al., 2014) may require the expected correction due to the baryonic effects. Our results are highly relevant to such searches, suggesting LMC-mass hosts are more likely to maintain a significant amount of substructure in the mass ranges to which such methods are sensitive.

4 Dwarf satellites in LMC-like hosts

A core prediction of $\Lambda$CDM is that halos and subhalos act as sites of galaxy formation. We therefore expect that some fraction of the surviving dark matter subhalos discussed above in Section 3 surrounding our LMC-mass hosts will host a luminous component consisting of stars and gas. The fact that even faint galaxies are much easier to detect than dark subhalos makes the luminous companions of LMC-mass systems a direct and testable prediction of the galaxy formation + $\Lambda$CDM model. In light of recent results by Gaia regarding proper motions of several MW dwarfs, and under the assumption that the availability of 6D information allows one to reconstruct the previous associations of the LMC to other fainter dwarfs (Sales et al., 2017; Kallivayalil et al., 2018), predictions are needed on the number and distribution of visible dwarf satellites expected to orbit around Magellanic dwarfs.

Fig. 4 highlights the mapping between the stellar mass content and the maximum circular velocities of the subhalos hosting luminous satellites with $M_{\ast}$ $\gtrsim$ 1$\times 10^{4}$ M*⊙* of our LMC-mass hosts. At $z=0$, the relation shows significant scatter, meaning that, at a given $V_{\text{max}}$, the corresponding stellar content may vary by several orders of magnitude. This effect is partially explained by the tidal stripping of subhalos within the LMC-mass host. This is shown by the gray open symbols indicating the same relation but using the peak value ever obtained for the maximum circular velocity ($V_{\text{peak}}$) of each satellite, roughly corresponding to its $V_{\text{max}}$ at infall. A similar increase in scatter in the $M_{\ast}$ $-$ $V_{\text{max}}$ relation has been observed for satellites in MW-mass hosts (e.g. Fattahi et al., 2013), however, scatter in the $M_{\ast}$ $-$ $V_{\text{peak}}$ relation is lower. This could have implications for abundance matching relations at small mass scales (see e.g. Garrison-Kimmel et al., 2017a). We return to tidal effects in sub-section 4.3. Additional factors contributing to the scatter is the stochasticity of the galaxy formation process near the low-mass end, as well as details of the assembly history of the halo and its inner dark matter density (Fitts et al., 2017).

We note that only a small fraction of the dark matter companions host a luminous dwarf companion (defined here as a galaxy within its host virial radius that has a stellar mass $M_{\ast}$ $\gtrsim 10^{4}$ M*⊙* in m11q). The bottom panel of Figure 4 shows the subhalo occupation fraction: the cumulative number of luminous $M_{\ast}$ $\gtrsim$ 1$\times 10^{4}$ M*⊙* satellites at a given $V_{\text{max}}$ divided by the cumulative number of dark + luminous subhalos at that $V_{\text{max}}$. We have stacked the subhalos and companions of all hosts to achieve a more complete sample across all ranges of $V_{\text{max}}$. To account for the different resolution in our runs, we only consider subhalos above the minimum $V_{\text{max}}$ thresholds introduced in Table 1 for each run.

The occupation fraction for galaxies with $M_{\ast}$ $\gtrsim$ 1$\times 10^{4}$ M*⊙* quickly decreases from about unity for $V_{\text{max}}$ $\geq 20$ km/s to only a few percent at $V_{\text{max}}$ $\sim 5$ km/s. For reference, our prediction is that about half of the subhalos with $V_{\text{max}}$ $=15$ km/s will host a luminous dwarf, the rest remaining dark or below the $M_{\ast}$ $=10^{4}$ M*⊙* resolved in our runs. These variations are commonly referred to as ‘stochasticity’ in the galaxy formation model, as dark matter subhalos of comparable mass may vary their stellar content by several orders of magnitude, including remaining totally dark, and predicting the exact form of the occupation fraction as a function of $V_{\text{max}}$ must consider many factors of the evolution of subhalos, such as their merger histories and accretion time, as well as external factors such as the onset and end of reionization and the form of the ionizing background.

4.1 New proper motions from Gaia

From the observational side, determining the dwarf satellites associated to the LMC prior to infall into the MW is not straightforward. Since tidal stripping due to the MW potential has already begun, the material once associated to the LMC in the past does not necessarily cluster around it today. However, because the LMC is inferred to be most likely in its first pericenter passage (Kallivayalil et al., 2013), cosmological simulations suggest that the stripped material may still retain its phase-space coherence, opening an avenue to disentangle previous associations (Sales et al., 2011, 2017; Deason et al., 2015; Jethwa et al., 2016). This coherence means that all subhalos within the LMC at infall are therefore expected to be distributed on the sky following the projection of positions and velocities of the LMC’s orbit. Note that final membership requires of a combination of position on the sky, galactocentric distance, and 3D velocity to be satisfied simultaneously.

This imposes strong constraints on the orbital poles expected for early companions of the LMC and can be used to single out possible associations. This criteria was used in Sales et al. (2011) to conclude that, with the exception of the SMC, no other classical dwarf with available proper motions at that time was consistent with an LMC association. With the arrival of new Gaia DR2 data, proper motions are now available for many dwarfs in the MW, including classical and ultrafaint dwarfs (we assume the cutoff for ultrafaints to be a stellar mass of $M_{\ast}$ $<10^{5}$ M*⊙*, as in Bullock & Boylan-Kolchin 2017). With the new 6D information, Kallivayalil et al. 2018 confirmed a likely association for 4 ultrafaints: Car2, Car3, Hor1 and Hyd1. The authors suggest follow up measurements for the confirmation of Dra2 and Hyd2, which have incomplete proper motion information. On the other hand, several of the classical dwarfs have updated proper motion measurements and their membership needs to be re-evaluated.

We have repeated the analysis in Sales et al. (2011, S11) but now using the Gaia DR2 proper motions presented in (Gaia Collaboration et al., 2018, H18). In particular, for classical dwarfs satisfying the galactocentric distance, position on the sky, and radial velocity constraints, S11 lists predictions for the orbital angular momentum expected in case of association. Following S11, we use a Cartesian coordinate system centered on the Milky Way, with $x$ in the sun-galactic center towards $l=0$°, $y$ towards $l=90$°in the direction of Galactic rotation, and $z$ coincident with the disk angular momentum, towards $b=90$°. For completeness, we list in Table 3 the $j_{x}$, $j_{y}$ and $j_{z}$ (all normalized to | $\vec{j}$ |) of all dwarfs included in H18 that were not part of the Kallivayalil et al. (2018) analysis. Errors are propagated from the quoted errors in H18 based on our calculation of $j_{x}$, $j_{y}$ and $j_{z}$. These calculated values are labelled ‘obs’ for each individual galaxy. However, of most relevance to this work are the dwarfs for which S11 had pre-determined possible association based on sky positions and radial velocity. For those cases, we also list the predicted angular momenta for the first pericentric passage ($t_{\text{1p}}$). Note that Sculptor and Leo II become a possible association only on a second-pericentric passage according to S11 ($t_{\text{2p}}$).

With the newest proper motions from Gaia, Carina and Fornax are also compatible with having been accreted as part of the LMC system (‘obs’ and $t_{1p}$ are consistent with each other within 1 $\rm\sigma$). With stellar masses of 3.8$\times 10^{5}$ M*⊙and 2.0$\times 10^{7}$ M⊙* , respectively, this newly confirmed pair of galaxies fills the classical dwarf scale (10${}^{5}\lesssim$ $M_{\ast}$ $\lesssim$ 107 M*⊙) in the satellite mass function of the LMC, which was previously populated only by ultrafaint dwarfs (Car2, Car3, Hor1, Hyd1 aside from the from the relatively bright SMC with $M_{\ast}$ $\sim$ 4.6$\times 10^{8}$ M⊙* ).

Pardy et al. (2019) suggest the possibility of LMC co-infall for Carina and Fornax by constraining their projected 2D orbital poles to within 30°of that of the LMC. We expand on this claim by using a more stringent criteria: matching angular momentum orientation with the one expected for and LMC-like debree at the same position on the sky of Carina and Fornax. Notice that the radial velocities for both dwarfs have been already found consistent with an LMC association in previous work (Sales et al., 2011, 2017). Also worth highlighting, the galactocentric distance of Carina is in good agreement with predictions of association whereas in the case of Fornax, the measured galactocentric distance ($\sim 140$ kpc) places it beyond the $\sim 100$ kpc preferred location of the debree (see for instance Fig.6 in Sales et al., 2011). The large distance of Fornax is more consistent with the scenario of a more massive infall halo mass for the LMC (whereas previous predictions were based on an LMC-analog with $\sim 10$ times lower mass), which would allow for a more extended distribution of the associated material. This caveat is an important one to bear in mind, and invites further investigation.

Note that Ursa Minor seems to meet the criteria set forth by S17, ($j_{x}<0$ and $|j_{x}|>>|j_{y}|,|j_{z}|$). However, its nearly perfect radial orbit ($j_{x}=-1$, $j_{y}\approx j_{z}\approx 0$) as well as its position in a completely distinct region of the sky than predicted for LMC debris (and where all currently known LMC satellites reside; see S17 Figure 1) suggest a non-Magellanic origin for Ursa Minor.

We summarize in Table 4 a complete list of MW dwarfs with their current understanding of association to the LMC. The rightmost columns labelled ‘possible’ (if follow up is needed) and ‘confirmed’ (if enough information exists to make the claim) with the relevant references. The previously confirmed satellites of the LMC (by S17 and K18) include Car2, Car3, Hor1, Hyd1 and the SMC. No label means that a given galaxy is unlikely to be associated with the LMC given the current data. Galaxies confirmed by our calculations using Gaia DR2 are labelled ‘this work’.

4.2 Simulated LMC Satellite Populations

With new observational context to the number of dwarf galaxies consistent with co-evolution and co-infall with the LMC, we can examine these results in a cosmological context. We provide this context by analyzing the satellite population of $\Lambda$CDM cosmological zoom-in simulations of isolated LMC-mass hosts. The left panel in Fig. 5 shows the stellar mass function of LMC satellites in FIRE (colored lines refer to the same simulations as previous figures, with the dashed lines representing an extrapolation to $M_{\ast}$ $\sim 10^{4}$ M*⊙* for the runs with resolution $m_{\text{bary}}$ = 7070 M*⊙*). In dark gray we show the observed stellar mass function of LMC satellties inferred from the kinematics of MW dwarfs from Gaia DR2 data, using starred symbols for the confirmed associations (SMC, Carina, Fornax, Hyd1, Car2, Hor1, and Car3) and in triangles including all ‘possible’ associations to the LMC, as determined by S17 and K18.

We find an overall good agreement between the inferred satellite population of the LMC and our simulated analogs. Our simulations predict between 1 and 5 classical satellites of the LMC, in agreement with the observational estimate of 3 for the LMC (SMC, Carina, and Fornax). There is an interesting mass dependence on the ability to predict relatively massive satellites for an LMC-like host. Only the two highest mass FIRE hosts (m11d with $M_{\text{200m}}$=2.8$\times 10^{11}$ M*⊙, and m11v with $M_{\text{200m}}$=2.9$\times 10^{11}$ M⊙) are able to reproduce the high-mass end of the LMC’s satellite mass function, in very close agreement with the halo mass estimates of the LMC ($\sim$3$\times 10^{11}$ M⊙) from other methods based on abundance matching (Behroozi et al., 2013b; Moster et al., 2013) and circular velocity measurements (van der Marel & Kallivayalil, 2014). In fact, the average halo mass of a LMC-SMC system in the EAGLE simulations is $\sim$3$\times 10^{11}$ M⊙* (Shao et al., 2018b; Cautun et al., 2019). The remaining three centrals with halo mass \sim 1.5$$\times 10^{11} M*⊙*, tend to have lower mass companions than the SMC. On the other hand, all runs have at least one satellite within a factor of two the stellar mass of Fornax, supporting its association to the LMC as suggested by the newly released Gaia kinematics.

One should keep in mind that the LMC-SMC association itself is rather unusual. Previous work on LMC-SMC selected pairs have showed them to be rare, though not impossible (Boylan-Kolchin et al., 2011a; Stierwalt et al., 2015). For example, Besla et al. (2018) used Illustris and SDSS to predict that the number of companions with M_{\ast}$$\sim 2$\times 10^{8}$ M*⊙* per LMC-mass dwarf is roughly 0.02 once projection effects have been taken into account. It is unclear how this figure changes with host mass, but following our results on the trend with virial mass, the likelihood of such a companion should increase if the LMC halo is on the massive end of the 1 $-$ 3 $\times 10^{11}$ M*⊙* range666m11v was specifically selected to be a multimerger, and thus cannot aid a discussion of the cosmological frequency of large satellites. In general, the number of classical dwarfs in our FIRE centrals is in good agreement with earlier predictions from semi-analytical calculations in Sales et al. (2013, light-blue dotted line) taken from the Millennium-II simulations (Boylan-Kolchin et al., 2009). However, the overall slope on the massive end of the FIRE runs is shallower than that in Sales et al. (2013), probably a result of a slightly different stellar mass - halo mass relation for low mass dwarfs in the semi-analytical catalog than in our hydrodynamical runs.

For the ultrafaint regime, our two highest resolution runs, m11q and m11c, predict $5-8$ companions with $M_{\ast}$ $\geq 10^{4}$ M*⊙* , which is in good agreement with the 5 inferred for the LMC from proper motion observations of MW dwarfs (Kallivayalil et al., 2018). In general, if extrapolating the medium resolution runs to $M_{\ast}$ $\geq 10^{4}$ M*⊙* (which, according to the highest resolution runs, would be a reasonable approach) the predicted number of satellites above this stellar mass may be as high as 11 for m11d. Although such predictions should be taken with caution given the low number statistics and numerical resolution, the assumption that some of the LMC ultrafaint companions still await discovery is a reasonable one. In particular some of the brighter ‘possible’ associations already identified, Hyd2 with $M_{\ast}$ $\sim$ 1.4$\times 10^{4}$ M*⊙* and Dra2 with $M_{\ast}$ $\sim$ 2.5$\times 10^{3}$ M*⊙*, might deserve a closer follow up to confirm or rule out their Magellanic origin.

Furthermore, our estimates for the number of ultra-faint companions can be regarded as lower limits, because of the implemented model for cosmic-ray heating in the ISM (in addition to the assumed UV background), which in these FIRE simulations induces too much heating in gas at early times ($z\gtrsim 10$) (Garrison-Kimmel et al., 2019). This effect has been explicitly tested with no observed impact on the population of classical dwarfs, however, the additional heating may in fact decrease the number of predicted ultrafaints compared to a method resulting in later reionization (for details, see Garrison-Kimmel et al., 2019, Section 3.3 and Appendix B). As such, our estimates are effectively a conservative lower limit and indicate that several yet-to-be-discovered companions to the LMC may yet exist.

Lastly, on the right panel of Fig. 5 we show the observed stellar mass function of MW satellites (black solid) along with satellites of M31 (black dot dashed). If the LMC brought along several of the dwarfs as estimated from the previous calculations (gray starred symbols), the MW satellite mass function must have looked rather different about 1 Gyr ago right before the infall of the LMC.777However, exact counts of luminous satellites are subject to change over this timescale as they are disrupted into streams, e.g. Sagittarius. This is shown as the gray dashed line, computed as the total MW satellites but subtracting the confirmed LMC associations. As usual, we define satellites as those within $r_{\text{200m}}$ of the host. In such case, the MW halo may have hosted a significantly lower number of dwarfs than now, although still in reasonable agreement with the predicted satellite population of Latte galaxies ($M_{\text{200m}}$ $\sim 10^{12}$ M*⊙, shaded in purple). The similarity between the high-mass end of the satellite mass function for the MW and the LMC argues once again for a rather massive pre-infall LMC halo, likely $\sim 3\times 10^{11}$ M⊙* and above, predicting several undiscovered ultra-faint dwarfs that were previously associated to the LMC.

4.3 The Predicted Dark Matter Content of LMC Satellites

Besides the number of dwarf galaxies expected around Magellanic-like systems, a further (arguably stronger) test for $\Lambda$CDM galaxy formation models is to reproduce the internal kinematics of the stars that are measured from observations. In the case of ultrafaint dwarfs, reported velocities from observations cover mostly the radius range $r<200$ pc; with at least half of the systems having measured velocities within 50 pc (Simon, 2019). Unfortunately, integrated quantities such as circular velocities are not yet converged in our simulations at such extreme small radii (typical gravitational softening $\epsilon_{\text{\tiny DM}}~{}\sim 20-40$ pc for all runs). We therefore analyze kinematic profiles of simulated classical dwarfs, but present a study of $V_{\text{max}}$ that can guide the conclusions in the ultrafaint regime.

Fig. 6 shows the circular velocity profiles of all simulated classical satellites ($M_{\ast}$ > $10^{5}$ M*⊙*) of our FIRE LMC-analogs. Lines are dashed below the dark matter convergence radius for each resolution as listed in Hopkins et al. (2018b). Individual curves are color coded according to the stellar mass content of each satellite, as indicated by the color bar. The simulation data hints to a correlation where the larger $M_{\ast}$ dwarfs will populate larger circular velocity subhalos, at least in the regime of classical dwarfs.

Observed dwarf spheroidals are also plotted as single points assuming $V_{\text{circ}}$($r_{1/2}$) = $\sqrt{3\langle\sigma_{\text{los}}\rangle}$ with $\sigma_{\rm los}$ the observed line-of-sight stellar velocity on the $y$-axis and the half-light radii $r_{1/2}$ on the $x$-axis (following Boylan-Kolchin et al. 2012). Sources for the observed $\sigma_{v,\ast}$ are listed in the legend of Fig. 6. Observed galaxies presumably associated with the LMC are shown in red while other MW dwarfs are plotted in gray. For completeness, we also include the full circular velocity profiles as a function of radius for Carina and Fornax from Wolf et al. (2010) and Strigari et al. (2010), respectively.

The predicted mass content for the brightest satellites of LMC-mass hosts, such as Carina and Fornax, are consistent with the observational values, supporting the possible association to the LMC inferred from their proper motions. It is important to highlight that the simulated curves in Fig. 6 are, if anything, lower limits to the true density of simulated dwarfs, as a smaller softening and increased resolution would result in higher inner velocities,888We have tested this explicitly, see Fig. 9. Also see Springel et al. (2008) for an evaluation of softening effects on $V_{\text{circ}}$. which would still accommodate the observed values. The good agreement with observations of Carina and Fornax is therefore encouraging and highly suggestive of a possible membership to the LMC group.

For the ultrafaint regime, we turn our analysis to maximum circular velocities, as we expect them to be well converged. Following Eq. 10 in Springel et al. (2008), the correction due to numerical effects is at most $\sim$0.1 km/s, or $\lesssim$1%. Interestingly, a closer inspection of our simulations reveals that tidal stripping is likely to play a major role on the present-day dark matter content of LMC companions, particularly in the ultrafaint regime. This is clearly seen in Fig. 7, which shows a comparison between peak maximum circular velocity achieved throughout a subhalo’s history ($V_{\rm peak}$) as a function of the circular velocity measured at present day (V_{\text{max}}$$(z=0)) for all subhalos in our sample of LMC-mass hosts in FIRE. Dark subhalos are indicated with gray dots, while large squares and triangles highlight the location of classical and ultrafaint simulated dwarfs, respectively. The symbols are color coded by lookback time since accretion (lower values are more recent infalls), and confirm that some of the surviving ultrafaint satellites in an LMC-system may have been accreted as early as 12 Gyr ago. The two regions of dark and light shaded gray indicate factors of 2 and 10 decrease in circular velocity.

As indicated by Fig. 7, the ultrafaint LMC companions are narrowly distributed at $V_{\text{peak}}$ $\sim$ 20 km/s at infall but show a large spread in $V_{\text{circ}}$ today, 6 $-$ 20 km/s. We find the present-day median $V_{\text{max}}$ for ultrafaint LMCs to be 14.4 km/s. As expected, the latest ones to infall remain close to the 1-to-1 line, as tidal disruption has not had sufficient time to affect their properties. We find that he amount of tidal stripping experienced is not dependent on stellar mass. However, consider $(i)$ the narrow range of $V_{\text{peak}}$ predicted for ultrafaints, and $(ii)$ their low stellar masses, at which feedback is not expected to affect the inner DM distribution. Both effects combined establish an interesting correlation in ultrafaints that might be used to assess their likelihood of association to the LMC. We emphasize that low central densities are a necessary rather than sufficient condition to determine an association to the LMC since an early infall onto the MW for ultrafaints will also induce tidal stripping and associated lower inner densities. Determining which host a satellite was first associated with is a task left to orbital phase-space analysis, as in sub-section 4.1.

Although tidal stripping proceeds mostly outside-in, subhalos that are affected by tides also register a drop in their inner dark matter content (e.g., Hayashi et al., 2003; Brooks & Zolotov, 2014). It is therefore expected that ultrafaint satellites of the LMC might show lower inner velocities than objects of similar mass that were accreted more recently into the MW. This prospect is interesting since Gaia DR2 data suggests that several of the ultrafaint galaxies are likely on their first infall onto the MW (Simon, 2018). Identifying ultrafaint dwarfs with orbital properties consistent with that of the LMC and that simultaneously show the lowest inner densities may help to constrain the most likely LMC group members999Note that the MW environment may also be conducive to tidal stripping such that low inner densities may be a telling but not sufficient condition for association to the LMC.. From that perspective, the likely associations for Hyd1 with $V_{\text{circ}}$ = 4.7 km/s and Car2 with $V_{\text{circ}}$ = 5.9 km/s, both within $r\leq 100$ pc, or the possibly associated Hyd2 with $V_{\text{circ}}$ = 6.2 km/s, seem favored compared to the more dense Car3, Hor1, or Dra2 with $V_{\text{circ}}$ $\geq 8$ km/s at $r\sim$ 50 pc (data from Simon, 2019). However, higher resolution simulations as well as more accurate observations are needed in order to make more definitive claims.

5 Summary and Conclusions

We used five $\Lambda$CDM cosmological zoom in simulations of LMC-mass hosts ($M_{\text{200m}}$ $=1-3\times 10^{11}$ M*⊙*) in FIRE to examine both the dark and luminous substructure of such galaxies, against which we compare to the Latte suite of seven simulations of MW-mass hosts. We summarize our primary findings here, and discuss our conclusions in the following paragraphs.

By comparing dark matter only to hydrodynamic baryonic simulations, we show that suppression of dark matter substructure is less strong for LMC-mass hosts than for MW-mass hosts, since this suppression is caused by tidal interactions with the central baryonic galaxy (Garrison-Kimmel et al., 2017b). We therefore expect that LMC-mass galaxies are promising laboratories for subhalo detection. 2. 2.

We calculate orbital angular momenta for 10 observed MW-dwarfs using Gaia DR2 data (Gaia Collaboration et al., 2018) and show that Fornax and Carina are highly consistent with predictions for tidal debris of a simulated LMC-MW infall scenario from Sales et al. (2017). This brings the inferred satellite mass function of the LMC up to seven members, including the SMC and the four ultrafaints identified in Kallivayalil et al. (2018). 3. 3.

We compare this to the simulated satellite mass function of our five LMC-mass hosts in FIRE and find excellent agreement on the bright end with our higher mass halos. Our simulations suggest that more ultrafaints are expected for a halo comparable to that of the LMC. In addition, the population of MW satellites pre-LMC infall remains consistent with simulated MW-mass hosts in FIRE. 4. 4.

We find that the tidal disruption of simulated LMC satellites, indicated by a reduction in their maximum circular velocities after infall, is an important effect, even at such host mass scales. We therefore expect that dwarfs associated with the LMC should have lowered densities, though this is not sufficient criteria for association by itself as tidal effects from the MW are also expected to affect nearby dwarfs.

We have shown that LMC-mass centrals suppress much less dark matter substructure compared to MW-mass centrals, with a particularly striking difference in suppression patterns for the inner 0.2 $r_{\text{200m}}$ of the halo. There, we observe a factor of $\sim$3.5 reduction (from dark matter only to baryonic runs) in the number of subhalos at $V_{\text{max}}$ = 10 km/s in MW mass hosts versus a factor of $\sim$1.4 for LMC mass hosts. This suppression of substructure is mainly due to tidal interactions with the central baryonic disk. We identify two features of the LMC-mass hosts in FIRE that explains this: (i) a much lower stellar mass to halo mass ratio than that of MW-mass objects and (ii) a significantly shallower inner density profile (larger core). The first implies that the gravitational potential of the galaxy is less significant compared to that of the halo. The second implies that the shape of the potential is much smoother than that seen in DMO versions of similar halos and that seen in halos hosting larger central galaxies. This combination provides a friendlier environment in the baryonic runs of cored $M_{\text{200m}}$ $\sim 10^{11}$ M*⊙* halos, leading to less suppression of subhalos due to tidal disruption on these scales compared to the harsher effects reported for $10^{12}$ M*⊙* hosts (D’Onghia et al., 2010; Garrison-Kimmel et al., 2017b; Samuel et al., 2019).

We compare the simulated satellite population of the LMC-mass hosts to the set of observed MW dwarfs that is consistent with an LMC co-infall scenario. We present revised calculations following Sales et al. (2011) and using updated Gaia DR2 proper motion for classical dwarfs ($M_{\ast}$ $\geq 10^{5}$ M*⊙) from Gaia Collaboration et al. (2018). We find that Carina and Fornax are now compatible with a common infall along with the LMC system, and are added to the previously suggested associations to the SMC and the group of ultrafaint candidates Car2, Car3, Hor1, and Hyd1 introduced in Kallivayalil et al. (2018). We find generally good agreement with the satellite population inferred for the real LMC. One of our halos, m11d with $M_{\text{200m}}$ $\sim 3\times 10^{11}$ M⊙*, is able to accurately predict the mass distribution of the three largest LMC satellites: the SMC, Fornax, and Carina, providing theoretical support to the claims of association. Furthermore, the predicted circular velocities for LMC satellites in this mass range is in good agreement with measurements of Carina and Fornax.

For fainter satellites, the LMC-mass systems in FIRE host comparable to slightly more ultrafaint companions than observed with the real LMC. On average, $7\pm 2$ satellites are predicted above $M_{\ast}$ $=10^{4}$M*⊙*, in reasonable agreement with the 5 above that mass limit presumed associated to the LMC (Car2, Hyd1, Carina, Fornax, and the SMC). We find that ultrafaint companions of the LMC are expected to have experienced significant tidal disruption within the LMC potential, as measured by the decrease on the subhalo maximum circular velocities since infall. As such, lower dark matter inner densities together with their orbital parameters (Kallivayalil et al., 2018), may help identify those ultra faint dwarfs that infalled onto the MW as part of the LMC group.

In summary, if the SMC, Carina, and Fornax are former satellites of the LMC, this may favor a relatively massive dark matter halo mass for the LMC prior to infall onto the MW, $M_{\text{200m}}$ $\sim 3\times 10^{11}$ M*⊙. This would push the expected ultrafaint dwarf numbers to the upper end of our predicted range, suggesting that some LMC associations still await discovery. According to our simulations, the missing dwarfs will lie roughly $30-80$ kpc from the LMC at infall and have $M_{\ast}$ $\sim 10^{4}$ M⊙* and $V_{\rm max}\approx 15$ km/s. This implies a relatively low central dark matter density that can be used as an additional membership criteria to discern from other, more recent individual infalls onto the MW. Based on their partially known orbital parameters and their low velocity estimates, the most promising candidate is Hyd2. Follow up observations are needed to confirm or dismiss the association.

In addition to LMC, this work presents the first observational and testable predictions using hydrodynamical cosmological simulations for the satellite mass function of a Magellanic-mass system down to the ultrafaint regime (semi-analytic modelling has been used for similar predictions, e.g. Dooley et al. 2017; Bose et al. 2018). Similar to how the study of satellite population of MW-mass hosts pushed forward our understanding of galaxy formation and cosmology in the past, the large predicted number of isolated dwarfs combined with upcoming deep surveys and wide field-of-view instruments such as WFIRST, as well as searches for satellites of LMC-mass hosts beyond the MW (e.g. MADCASH; Carlin et al., 2016), may turn the study of dwarf-dwarf systems into a valuable and essential test of the $\Lambda$CDM model.

Acknowledgments

The authors would like to thank Phil Hopkins for facilitating access to the FIRE runs and stimulating discussions of the simulations and science herein. We also thank Coral Wheeler, Shea Garrison-Kimmel, and the rest of the FIRE local universe collaboration for their thought-provoking insight and advice. We are grateful to the anonymous referee for a constructive report that helped improve this manuscript, as well as to Marius Cautun and Marcel Pawlowski for their helpful suggestions. LVS acknowledges support from NASA through the HST Programs AR-14582, AR-14583 and from the Hellman Foundation. AW was supported by NASA, through ATP grant 80NSSC18K1097 and HST grants GO-14734 and AR-15057 from STScI. MBK acknowledges support from NSF grant AST-1517226 and CAREER grant AST-1752913 and from NASA grants NNX17AG29G and HST-AR-14282, HST-AR-14554, HST-AR-15006, and HST-GO-14191 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. TKC was supported by NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. We ran simulations using 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 Additional Information

Here we include additional information relevant to the discussions above. We test the effects of resolution in our simulations by comparing the subhalos of m11q run at the ‘high’ and ‘normal’ resolutions as described in sub-section 2.2. We examine the subhalo $V_{\text{max}}$ functions as described in sub-section 2.1 in Figure 8 and the circular velocity profiles of subhalos with 8 km/s < $V_{\text{max}}$ < 12 km/s as described in sub-section 4.3 in Figure 9. We examine the radial distribution of subhalos in LMC-mass hosts and MW-mass hosts as described in sub-section 3.1 in Figure 10. We include our observational data set as described in sub-section 4.1 in Table 4, including stellar mass, galactocentric coordinates and radial velocities, distance to the LMC, as well as the current status of association to the LMC system.