Be it therefore resolved: Cosmological Simulations of Dwarf Galaxies with Extreme Resolution

TL;DR

This study uses ultra-high-resolution cosmological simulations to explore the formation, structure, and star formation histories of dwarf galaxies, revealing insights into their internal structures, metallicities, and observability.

Contribution

It presents the first ultra-high-resolution simulations resolving individual supernova remnants, providing new insights into dwarf galaxy formation and internal structures.

Findings

Ultra-faint dwarfs have early star formation truncated by reionization.

Classical dwarfs continue star formation to low redshift.

Simulated dwarf properties match observed sizes and dynamics.

Abstract

We study a suite of extremely high-resolution cosmological FIRE simulations of dwarf galaxies ($M_{\rm halo} \lesssim 10^{10}$$M_{\odot}$), run to $z=0$ with $30 M_{\odot}$ resolution, sufficient (for the first time) to resolve the internal structure of individual supernovae remnants within the cooling radius. Every halo with $M_{\rm halo} \gtrsim 10^{8.6} M_{\odot}$ is populated by a resolved {\em stellar} galaxy, suggesting very low-mass dwarfs may be ubiquitous in the field. Our ultra-faint dwarfs (UFDs; $M_{\ast}<10^{5}\,M_{\odot}$) have their star formation truncated early ($z\gtrsim2$), likely by reionization, while classical dwarfs ($M_{\ast}>10^{5} M_{\odot}$) continue forming stars to $z<0.5$. The systems have bursty star formation (SF) histories, forming most of their stars in periods of elevated SF strongly clustered in both space and time. This allows our dwarf with $M_{\ast}/M_{\rm halo} > 10^{-4}$ to form a dark matter core $>200$pc, while lower-mass UFDs exhibit cusps down to $\lesssim100$pc, as expected from energetic arguments. Our dwarfs with $M_{\ast}>10^{4}\,M_{\odot}$ have half-mass radii ($R_{\rm 1/2}$) in agreement with Local Group (LG) dwarfs; dynamical mass vs. $R_{1/2}$ and the degree of rotational support also resemble observations. The lowest-mass UFDs are below surface brightness limits of current surveys but are potentially visible in next-generation surveys (e.g. LSST). The stellar metallicities are lower than in LG dwarfs; this may reflect pre-enrichment of the LG by the massive hosts or Pop-III stars. Consistency with lower resolution studies implies that our simulations are numerically robust (for a given physical model).

Tables

Table 1. Table 1: Properties of dwarfs in the suite. Each row lists a different resolved, central or satellite galaxy at z = 0 𝑧 0 z=0 . Columns give: (1) M halo subscript 𝑀 halo M_{\rm halo} : halo mass. (2) R vir subscript 𝑅 vir R_{\rm vir} : virial radius. (3) V max subscript 𝑉 max V_{\rm max} : maximum circular velocity. (4) M ∗ subscript 𝑀 ∗ M_{\ast} : bound stellar mass after removal of satellites and contamination. (5) R max subscript 𝑅 max R_{\rm max} : radial extent of stars (maximum radius of any bound star), as determined from visual inspection. (6) R 1 / 2 subscript 𝑅 1 2 R_{1/2} : mean projected (2D) half-stellar-mass radius. (7) v rot / σ subscript 𝑣 rot 𝜎 v_{\rm rot}/\sigma : Ratio of the stellar velocity shear v rot subscript 𝑣 rot v_{\rm rot} to dispersion σ 𝜎 \sigma . (8) g soft subscript 𝑔 soft g_{\rm soft} : typical minimum gas gravitational+hydrodynamic force softening reached in star-forming gas (this is adaptive). (9) DM soft subscript DM soft {\rm DM}_{\rm soft} : dark matter force softening (held constant). (10) n crit subscript 𝑛 crit n_{\rm crit} : minimum gas density required for star formation (in addition to self-shielding, Jeans instability, and self-gravity).

	$M_{\mathrm{halo}}\left({10}^9{M}_{\odot }\right)$	$r_{\mathrm{vir}}\left(\mathrm{kpc}\right)$	$V_{\max }\left({\mathrm{kms}}^{-1}\right)$	$M_{\star }\left({10}^3{M}_{\odot }\right)$	$r_{\max }^{\ast }\left(\mathrm{kpc}\right)$	$R_{1/2}\left(\mathrm{pc}\right)$	$v_{\mathrm{rot}}/\sigma$	$g_{\mathrm{soft}}\left(\mathrm{pc}\right)$	${\mathrm{DM}}_{\mathrm{soft}}\left(\mathrm{pc}\right)$	$n_{\mathrm{crit}}\left({\mathrm{cm}}^{-3}\right)$
	Resolved Central Galaxies in Our High-Resolution ($30M_{\odot }$) Simulations ($>100$ star particles in AHF-Identified Halo)
${\mathrm{m10q}}_{30}$	7.7	51	34	5200	7.7	720	0.17	0.40	14	1000
${\mathrm{m10q}}_{30}\mathrm{Sat}$	0.34	6.3	16	1.2	1.2	560	0.85	0.40	14	1000
${\mathrm{m10v}}_{30}$	9.0	54	30	330	8.2	330	0.45	0.10	14	1e5
${\mathrm{m10v}}_{30}B$	3.2	37	24	41	5.6	280	0.20	0.10	14	1e5
${\mathrm{m10v}}_{30}C$	1.1	26	16	2.9	3.9	540	0.39	0.10	14	1e5
${\mathrm{m10v}}_{30}D$	0.75	24	16	3.7	4.8	860	0.37	0.10	14	1e5
${\mathrm{m10v}}_{30}\mathrm{Sat}$	0.52	20	13	2.0	2.8	430	0.29	0.10	14	1e5
${\mathrm{m10v}}_{30}F$	0.44	20	12	1.9	2.5	560	0.66	0.10	14	1e5
${\mathrm{m10v}}_{30}G$	0.27	17	12	2.3	2.5	570	1.1	0.10	14	1e5
${\mathrm{m09}}_{30}$	2.5	35	22	12	4.0	200	0.56	0.10	14	1e5
${\mathrm{m09}}_{30}B$	0.67	23	15	1.8	3.4	620	0.79	0.10	14	1e5
	Resolved Central Galaxies in Our Standard-Resolution ($250M_{\odot }$) Simulations ($>100$ star particles in AHF-Identified Halo)
${\mathrm{m10q}}_{250}$	7.5	51	34	2700	7.7	550	0.19	1.0	29	1000
${\mathrm{m10v}}_{250}$	8.4	53	30	300	8.0	310	0.31	1.0	29	1000
${\mathrm{m10v}}_{250}B$	2.7	37	24	66	5.5	350	0.14	1.0	29	1000
${\mathrm{m09}}_{250}$	2.5	36	22	27	5.3	420	0.14	1.0	29	1000

Full text

††pagerange: Be it therefore resolved: Cosmological Simulations of Dwarf Galaxies with Extreme Resolution–References††pubyear: 2018

1 Introduction

Although the currently favored cosmological paradigm – $\Lambda$ Cold Dark Matter Theory ($\Lambda$CDM) – has been widely successful in predicting the counts, clustering, colors, morphologies, and evolution of galaxies on large scales, as well as a variety of cosmological observables (Eisenstein et al., 2005; Viel et al., 2008; Reid et al., 2010; Komatsu et al., 2011), several challenges have arisen to this model in recent years, most of them occurring at the smallest scales – those of dwarf galaxies ($M_{\rm\star}\lesssim 10^{9}M_{\odot}$; see e.g. Bullock & Boylan-Kolchin 2017 for a review.) Among these small-scale challenges, perhaps best known is the Missing Satellites Problem (MSP; Klypin et al. 1999; Moore et al. 1999; Bullock 2010): counts of galaxies predicted from a naive assignment of stellar mass to dark-matter only simulations of Milky Way (MW)-mass galaxies drastically over-predicts the actual number of currently observed dwarf galaxies around the MW.

The severity of the missing satellites problem is sensitive to the low-mass edge of galaxy formation: any halo mass threshold below which galaxy formation cannot proceed will result in firm predictions for the abundance and distribution of low-mass galaxies around the Milky Way. In fact, Graus et al. (2018) suggest that the radial distribution of Milky Way satellites requires galaxy formation to persist in dark matter halos with virial temperatures below the atomic cooling limit, potentially presenting an issue in the opposite sense of the classic missing satellites problem.

Whether there is a well-defined halo mass scale at which galaxy formation ceases to operate, and the precise location of this low-mass cutoff, remain unknown. Any low-mass cutoff almost certainly would be affected by the timing of reionization – both the onset and the end – as well as by the overall flux of ionizing photons, the spectrum of the radiation, the proximity to more massive structure and the self-shielding ability of the gas itself (Efstathiou, 1992; Dijkstra et al., 2004; Hoeft et al., 2006; Weinmann et al., 2007; Oñorbe et al., 2015). Likewise, the properties of the lowest-mass galaxies that do manage to form should reflect the imprint of the cosmic reionizing background. So-called “fossils” of reionization, as first proposed by Ricotti & Gnedin (2005), are galaxies that managed to form some stars before having their star formation shut down by reionization (Bullock, Kravtsov & Weinberg, 2000). Observations of six UFDs around the Milky Way show that they have stellar ages that are indistinguishable from ancient globular clusters, with the entirety of their star formation occurring before $z\sim 2$ (Brown et al., 2014). This suggests that ultra-faint satellites of the Milky Way are indeed fossils of the reionization era.

However, the fact that all of the UFDs considered in Brown et al. (2014) were satellites, rather than isolated galaxies, makes it more difficult to distinguish this effect from other quenching mechanisms, such as ram-pressure stripping (Gunn & Gott, 1972). Wheeler et al. (2015, hereafter W15) used simulations run with the first generation Feedback in Realistic Environments (FIRE; Hopkins et al., 2014)111http://fire.northwestern.edu. FIRE uses the pressure-entropy version of SPH hydrodynamics. and with baryonic particle masses $m^{\rm p}_{\rm bar}=250\rm\,M_{\odot}$ to show that both isolated and satellite dwarfs in the UFD mass range had uniformly ancient stellar populations, suggesting that they were indeed reionization fossils. Additionally, Rodriguez Wimberly et al. (2018) removed halos from dark-matter-only simulations (Garrison-Kimmel et al., 2014) that would have been destroyed by the galactic disk according to their pericentric distance (Garrison-Kimmel et al., 2017) and showed that there is a vanishingly small probability that all of the observed ultra-faints fell into the Milky Way by $z=2$, making environmental quenching an unlikely explanation for their ancient stellar ages.

Another well known and long-standing small-scale challenge to $\Lambda$CDM is the Core-Cusp Controversy (CCC; Flores & Primack, 1994; Moore, 1994; de Blok, 2010), in which the predicted density of dark matter halos as measured from dark-matter-only simulations suggests the presence of a central density “cusp”222Note that Baushev & Pilipenko (2018) argue that the predictions of cusps themselves may be numerical artifacts. ($\rho\sim r^{-1}$; Navarro, Frenk & White 1997), while observations of some dwarf galaxies suggest that the actual profile shape can flatten at the center, into a shallower-density “core” ($\rho\sim\rm const$; Salucci & Burkert 2000; van den Bosch & Dalcanton 2000; de Blok & Bosma 2002; Oh et al. 2008).

This issue has since been shown to be largely a result of comparing dark-matter-only simulations to observations, and that a dark matter core can be created in galaxies through repeated fluctuations in the gravitational potential via regular expulsions of the galactic gas supply from bursty star formation and its resulting feedback (Mashchenko, Couchman & Wadsley, 2006; Governato et al., 2012; Pontzen & Governato, 2012; Chan et al., 2015; Oñorbe et al., 2015; Tollet et al., 2016). However, the ability of supernova feedback to reduce the central density of the dark matter halo is limited by the competing effects of the halo potential and the total mass in stars formed (Garrison-Kimmel et al., 2013). Galaxy formation is highly inefficient in the regime of UFDs ($M_{\rm\star}\lesssim 10^{5}$), meaning the supernova energy input per unit binding energy of the dark matter halo is much lower than in higher mass systems. Di Cintio et al. (2014a) showed that galaxies with $M_{\rm\star}/M_{\rm{vir}}\lesssim 10^{-4}$ fail to produce enough star formation to significantly alter the inner halo density cusp to a core. Indeed, most cosmological $\Lambda$CDM simulations to date fail to form cores in UFDs (Munshi et al., 2013; Fitts et al., 2017, although see Read, Agertz & Collins 2016 for an idealized study). Several UFDs have been found to host globular clusters near their centers, leading some authors to argue that their existence and the lack of dynamical friction implied is evidence for cores in these objects (Goerdt et al., 2006; Amorisco, 2017; Caldwell et al., 2017), but there remains little direct observational evidence for cores in UFDs. If stronger evidence arises, it may require new physics in the dark sector (see, e.g., Bullock & Boylan-Kolchin 2017; Tulin & Yu 2018; Buckley & Peter 2018 for recent reviews).

Our understanding of the severity of these challenges to $\Lambda$CDM is complicated by computational difficulties in dealing with low-mass galaxies. Most cosmological hydrodynamic simulations of dwarfs run to $z=0$ – including smoothed particle hydrodynamics (SPH), moving-mesh or adaptive-mesh-refinement (AMR) methods – have a quasi-Lagrangian mass resolution333Mass resolution is set by a combination of, e.g., the numerical cell/particle masses and the physical minimum gas mass at which a self-gravitating structure can be identified, which are similar in our study but do not have to be in general (see Hopkins et al. 2018b). ranging from $\sim 250-10^{4}\rm\,M_{\odot}$. This means that in UFDs, where the mass of the stellar content can be as low as only a few 100s of solar masses, these simulations may be unable to resolve the galaxies; if they do, these “galaxies” often have very few, and sometimes even a single, star particles (Sawala et al., 2014; Wheeler et al., 2015; Munshi et al., 2017). This problem becomes more and more challenging at lower masses, because $M_{\rm\star}/M_{\rm{vir}}$ drops rapidly (meaning a larger and larger number of resolution elements are required to represent the small number of stars that should be present at lower $M_{\rm{vir}}$).

As more and more UFDs are discovered and scrutinized at increasing levels of detail, it becomes imperative for cosmological hydrodynamic simulations to continue to push to higher resolution to allow for comparisons with observations. The stellar particle count must not only be large enough to ensure that the galaxy actually forms (i.e. is not a random grouping of particles), but it also must be sufficient to accurately estimate, e.g., the half-light radius for the galaxies, as well as rotation or other higher-order properties. Cosmological simulations must be able to reliably resolve UFDs to make firm predictions for the next generation of telescopes.

In this paper, we introduce a new set of high-resolution cosmological hydrodynamic zoom-in simulations (GIZMO/FIRE-2) of isolated dwarf galaxies with baryonic particle masses of $m^{\rm p}_{\rm bar}=30\rm\,M_{\odot}$ – the highest resolution ever run to $z=0$. This new generation of $m^{\rm p}_{\rm bar}\sim 10M_{\odot}$ simulations marks a transition point between simulations that treat star formation within a single stellar population in the aggregate and simulations that model the collapse and fragmentation of a molecular cloud into individual stars, and allows us to probe smaller physical scales than previously possible in cosmological simulations. This in turn enables the comparison between a larger set of “resolved” simulated UFDs and an ever-increasing set of observations at the low-mass end of galaxy formation. We introduce the suite in Section 2, give an overview of results and compare to observations in Section 3, including the star formation histories (SFHs, Section 3.2), halo structure (Section 3.3), kinematics (Section 3.4), and chemical abundances (Section 3.5). We conclude in Section 4.

2 Simulations

2.1 Resolution and Motivation

Our highest resolution suite consists of cosmological zoom-in simulations of Lagrangian volumes surrounding three “primary” isolated dwarf galaxy halos ($\rm m10q_{30}$, $\rm m10v_{30}$, $\rm m09_{30}$), each with baryonic particle mass of $m^{\rm p}_{\rm bar}=30\rm\,M_{\odot}$ ($m_{\rm dm}$ is larger by the universal baryon fraction), and $z=0$ virial masses444We define virial overdensity with the spherical top hat approximation of Bryan & Norman (1998). $\sim 2-10\times 10^{9}\rm\,M_{\odot}$ (see Table 1 for a full list of simulation properties).

Before going forward, we stress that the mass resolution achieved here, $\sim 30\,M_{\odot}$, is not simply an incremental improvement. It reaches a critical physical scale where the cooling radius of a SN remnant (approximately the radius enclosing $\sim 3000\,\rm\,M_{\odot}$, with only very weak residual dependence on metallicity or gas density) is resolved with $\gtrsim 100$ elements, which is essential for capturing the basic dynamics. Furthermore, many independent studies (e.g., Lapi, Cavaliere & Menci 2005; Kim & Ostriker 2015; Martizzi, Faucher-Giguère & Quataert 2015; Walch & Naab 2015; Hu 2018; Hopkins et al. 2018a) have shown that with mass resolution of $\lesssim 100\,M_{\odot}$, predictions for SN feedback become nearly independent of the detailed numerical implementation: whether one simply “dumps” $\sim 10^{51}\,{\rm erg}$ into surrounding gas in thermal or energy, or applies a more sophisticated injection model, the asymptotic behavior of the blastwave will converge to the same behavior at this resolution. As a result, Hopkins et al. (2018a) showed the predictions of galaxy-formation simulations become vastly less-sensitive to the sub-grid model for feedback. We will therefore compare various properties at our high resolution and a resolution $\sim 8\times$ poorer. In Hopkins et al. (2018b), the primary galaxies $\rm m10q_{30}$ and $\rm m10v_{30}$ are also studied at even lower resolution.

2.2 Numerical Methods

All details of the numerical methods and initial conditions555The ICs used here are publicly available at http://www.tapir.caltech.edu/~phopkins/publicICs, are presented in Hopkins et al. (2018b), where lower-resolution versions of these volumes were studied extensively. We briefly summarize essential elements here.

The simulations are run with the GIZMO (Hopkins, 2014)666A public version of GIZMO is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html code using the updated FIRE-2 implementation of star formation and stellar feedback from Hopkins et al. (2018b). FIRE-2 uses the “meshless finite mass” (MFM) Lagrangian finite-volume Godunov method for the hydrodynamics, accounts for gas heating from a variety of processes including the UV background777We adopt the UV background from the December 2011 update of Faucher-Giguère et al. (2009) (available here: http://galaxies.northwestern.edu/uvb/), which reionizes the universe rapidly around $z\sim 10$ and completes reionization by $z=6$, and was designed to produce a reionization optical depth consistent with WMAP-7 and local sources and cooling from $T=10-10^{10}\,K$, and models star formation via a sink-particle method in gas that is locally self-gravitating, Jeans-unstable, self-shielding/molecular, and exceeds a critical density $n_{\rm crit}$ (see table 1). We adopt a standard, flat $\Lambda$CDM cosmology with $h\approx 0.70$, $\Omega_{\rm m}=1-\Omega_{\Lambda}\approx 0.27$, and $\Omega_{\rm b}\approx 0.045$ (consistent with current constraints; see Planck Collaboration et al. 2014).888For the sake of comparison with other work, $\rm m09_{30}$ adopts slightly different cosmological parameters than the other two simulations. These differences are at the $\sim 1\%$ level, and matter far less than stochastic run-to-run variance.

Once stars form, stellar feedback from SNe (Ia & II), stellar mass-loss (O/B and AGB) and radiation (photo-electric and photo-ionization heating, and radiation pressure, with a five-band radiation transport algorithm; see Hopkins et al. 2018a) are included following a STARBURST99 stellar evolution model (Leitherer et al., 1999) based on the star particle age and metallicity (adopting a Kroupa 2002 IMF). In lower-resolution simulations, the feedback quantities (e.g. SNe rates, stellar spectra) can simply be IMF-averaged. At the high resolution here, however, this produces un-physical outcomes: for example, it effectively assigns $\sim 1/3$ of an O-star to each new star particle, rather than having $\sim 1/3$ of star particles in young stellar populations contain an O-star. We therefore follow Ma et al. (2015); Su et al. (2018b); Grudić & Hopkins (2018) to sample the IMF discretely. Note that while Su et al. (2018b) show in a variety of tests that this has little effect on galaxy properties when compared to run-to-run stochastic variation, Applebaum et al. (2018) show that stochastic sampling of the IMF makes feedback burstier, stronger, and quenches star formation earlier in small dwarf galaxies.

2.3 Analysis

We use the Amiga Halo Finder (Knollmann & Knebe, 2009, AHF;) to identify gravitationally bound members of each halo in post-processing. We initially require $>100$ bound star particles to consider a galaxy resolved. Because AHF can incorrectly assign central stellar halo stars to subhalos, we visually inspected each simulation to determine whether all star particles within the virial radius of each halo are physically associated with the central galaxy or with a satellite, in some cases adjusting the radial extent, $r^{*}_{\rm{max}}$, and stellar mass of the galaxy accordingly. All figures reflect these adjustments. Table 1 summarizes the properties of all dwarfs that meet this criterion. Note that a less conservative $\sim 15$-particle cut yields $>100$ additional galaxies; these will be studied in future work. Fig. 1 shows the DM and stellar distribution in the primary galaxies of each of our zoom-in simulations as well as in the satellites of $\rm m10q_{30}$ and $\rm m10v_{30}$.

3 Results

3.1 Halo and Stellar Masses

Fig. 2 shows the DM subhalo and satellite stellar mass functions of the three “primary” galaxies. Fig. 3 compares the stellar mass-halo mass relation of the simulations for all galaxies in the high resolution region we consider well-resolved to observational constraints implied by abundance-matching.

First, we note that while stochastic run-to-run variations (e.g. the magnitude of the largest starbursts at high redshift; see El-Badry, Weisz & Quataert 2017; Su et al. 2018b, a; Keller et al. 2019) can lead to changes as large as a factor $\sim 2$ in $M_{\ast}$, there does not appear to be any systematic dependence on resolution in these properties for the $\sim 4$ galaxies that are well-resolved at both resolution levels. We have also run volumes m10v and m10q at even lower resolution; at resolution $(30,\,250,\,2\times 10^{3},\,1.6\times 10^{4},\,1.3\times 10^{5})\,M_{\sun}$, the primary galaxy in m10v has stellar mass $M_{\ast}=(3.3,\,3.0,\,1.5,\,4.8,\,2.6)\times 10^{5}\,M_{\sun}$, and the primary in m10q has $M_{\ast}=(5.2,\,2.7,\,3.0,\,2.0,1.3)\times 10^{6}\,M_{\sun}$. This suggests the sub-grid algorithm that couples SNe mass, energy and momentum to surrounding gas particles properly handles the transition between explicitly resolved and unresolved SNe remnants, as it is specifically designed to do (see Hopkins et al. 2018a for extensive tests and discussion). That, in turn, is encouraging for the robustness of previous (lower-resolution) predictions from FIRE simulations.

Our simulated galaxies broadly sample the scatter about extrapolated abundance matching relationships that are tuned to reproduce the stellar mass function of the Local Group, suggesting that the underlying prescriptions in our runs should also reproduce observations when applied to LG-like environments at this resolution. However we stress this comparison is not rigorous: almost all of the simulated galaxies are below the detection limits (in surface brightness or stellar mass) of the observations used to calibrate the abundance-matching relations; moreover, the observations are of LG satellites, not isolated systems. For a more rigorous comparison of classical dwarf mass functions, see Garrison-Kimmel et al. (2018).

Our simulations suggest, at least absent a massive host such as the MW, UFDs form in all DM halos with $M_{\rm{vir}}^{z=0}\gtrsim 5\times 10^{8}\rm\,M_{\odot}$ (equivalently, $V_{\rm max}^{z=0}>13~{}{\rm km\,s}^{-1}$). This is just for galaxies with $\gtrsim 100$ star particles: we form galaxies with at least 15 star particles in all halos with $V_{\rm max}^{z=0}>10~{}{\rm km\,s}^{-1}$ and at least 1 star particle in all halos with $V_{\rm max}^{z=0}>7.5~{}{\rm km\,s}^{-1}$. To get a (very crude) sense of the implications for LG dwarf populations, we can compare to the number of $V_{\rm infall}>13\,{\rm km\,s}^{-1}$ halos+subhalos (where $V_{\rm infall}$ is $V_{\rm max}$ when the satellite crosses MW’s virial radius) in DM-only simulations of 12 LG analogues from Garrison-Kimmel et al. (2014). This naively predicts $\sim 400-610$ such MW satellites with $M_{\rm\star}\gtrsim 10^{3}\rm\,M_{\odot}$, $\sim 180-380$ “isolated” dwarfs above this mass in the Local Field999The Local Field is the immediate environment of the MW and M31 outside of each massive galaxy’s virial radii. (within $1\rm\,Mpc$), and up to $140$ satellites of other field dwarfs in the same volume, where the range spans one standard deviation from the mean for the sample.

However, the mere presence of a MW-mass galaxy at the center of a halo has been shown to completely destroy $\sim 1/3-1/2$ of all of its subhalos (Garrison-Kimmel et al., 2017), so the number of satellites may be significantly less than these rough estimates. Using a higher $M_{\rm\star}\gtrsim 10^{4}\,\rm\,M_{\odot}$ ($V_{\rm max}\gtrsim 20~{}{\rm km\,s}^{-1}$) threshold gives results qualitatively consistent with the abundance-matching inferences for these populations in Dooley et al. (2017). Interestingly, the fact that we do not see a sharp cut-off in halo mass for UFDs is also consistent with the recent study by Graus et al. (2018), which argued the radial distribution of MW satellites requires halos with $V_{\rm infall}\gtrsim 6-10~{}{\rm km\,s}^{-1}$ be populated with UFDs. We will study the less well-resolved, but more abundant, lower-mass UFD population in future work.

3.2 Star Formation Histories

Figs. 4-5 show the archaeological SFH for each galaxy (the distribution of formation times of all $z=0$ stars), in cumulative and differential form. In every case, the more massive galaxies have SF down to $z<0.5$, while all galaxies with $M_{\rm\star}<10^{5}\rm\,M_{\odot}$ had SF cease well before $z=2$. The least massive UFDs, with $M_{\rm\star}<10^{4}\rm\,M_{\odot}$ in halos with $M_{\rm{halo}}\lesssim 10^{9}\rm\,M_{\odot}$, have the most ancient populations, with SF ceasing before $z=10$ (the midpoint of reionization, as modeled here).101010The one exception appears to be a single star particle from $\rm m10v_{30}\,Sat$, likely a contaminant from its host $\rm m10v_{30}$. This is consistent with the SFHs of LG dwarfs (Brown et al., 2014). Since our galaxies do not have a massive host, it is clear that the UFD quenching was driven by reionization, not environmental effects. High-mass UFDs ($M_{\rm\star}\sim 10^{4-5}\,M_{\sun}$) have their accretion cut off by reionization at $z\sim 10-6$, but form stars for another $\sim 1$ Gyr (until $z\sim 2-4$). This residual SF comes from gas that was accreted pre-reionization and is self-shielding to the UV background, so continues to form stars until it is exhausted or blown out by SNe (see also Oñorbe et al. 2015).

We note that the primary m10q galaxy appears to “self-quench,” i.e. quench without external influence, and exhaust its cold gas before $z=0$. Although this does not occur in the specific lower-resolution run here (which has on-going SF to $z=0$), in several previous studies (see e.g. Su et al., 2018b) we have shown that this occurs semi-stochastically in lower-resolution runs (as a particularly large burst of SNe concurrently can eject the remaining small amount of cold gas in the galaxy and shut down SF). Roughly, comparing iterations of m10q with intentionally small perturbations to the ICs or run-time parameters in Su et al. (2018b), Hopkins et al. (2018b) and the ensemble of dwarfs studied in Fitts et al. (2017), we find this occurs with order-unity probability. Without a statistical sample, we cannot say if it is more likely in high-resolution runs, but it does indicate that such gas expulsion is not purely an artifact of lower resolution. Whether this is consistent with observations requires a larger sample: Geha et al. 2012 argue the non-star forming fraction in the field for galaxies with $10^{7}\rm\,M_{\odot}<M_{\rm\star}<10^{9}\rm\,M_{\odot}$ is zero, but several obviously quenched counter-examples with $M_{\rm\star}\la 10^{7}\,\rm\,M_{\odot}$ exist out to $>1$ Mpc from the LG (Karachentsev et al., 2001, 2014, 2015; Makarov et al., 2012; Makarova et al., 2017; Cole et al., 2014), and Fillingham et al. (2018) argue that the quenched fraction in the LG and field dwarfs requires some internal quenching (i.e. un-connected to a massive host) must occur at these masses.

Fig. 5 shows the SFR averaged in $10$ and $200$ Myr windows, analogous to timescales over which H$\alpha$ and UV continuum measurements are sensitive; we find the SFRs are bursty, as found previously at lower resolution (Muratov et al., 2015; Sparre et al., 2017; Faucher-Giguère, 2018). SFRs for the other simulated UFDs are similar to $\rm m09_{30}$ shown. The results are qualitatively similar at low/high resolution, although the periods of low SFR in the low-resolution runs drop to [math] while remaining finite at high resolution – this is simply an artifact of the discrete nature of our star particles (at $\dot{M}_{\ast}\sim 10^{-6}\,\rm\,M_{\odot}\,{\rm yr^{-1}}$, it requires $\sim 250\,$Myr to form a single star particle at $\sim 250\,M_{\sun}$ resolution). Quantitatively, whether the burstiness changes with resolution depends on our definition. The unit of SF (a star particle) is larger at low resolution, so sampling variations decrease at higher resolution. Based on the metric used in Hopkins et al. (2014), namely the dispersion in $\log(\langle\dot{M}_{\ast}(\Delta t_{1})\rangle/\langle\dot{M}_{\ast}(\Delta t_{0})\rangle)$ (where $\Delta t_{1}\sim 10\,$Myr and $\Delta t_{0}\sim$ Gyr), we actually find the burstiness slightly increases in higher-resolution runs. However, quantifying burstiness in terms of the fraction of the total $M_{\ast}$ formed in periods where $\langle\dot{M}_{\ast}(\Delta t_{1})\rangle>1.5\,\langle\dot{M}_{\ast}(\Delta t_{0})\rangle$ (Sparre et al., 2017), we find slightly decreasing burstiness at higher resolution (for $\Delta t_{1}\sim 10\,$Myr; for much longer averaging times, the effect vanishes). In all cases, the resolution dependence is small ($\sim 10\%$ for order-of-magnitude change in resolution).

3.3 Cusps and Cores

Fig. 6 shows the DM profiles of the primary galaxies vs. resolution. It is well-known that convergence in the halo mass profiles of $N$-body calculations is almost entirely determined by mass resolution (see e.g. Power et al., 2003), so the simulations here provide an important convergence test. Using the most conservative definition of a converged radius from Power et al. (2003), our ultra-high resolution DM profiles should be converged down to $\sim 60-100\,$pc; using the more-aggressive criterion from Hopkins et al. 2018b gives $\sim 30-40\,$pc. Inside of this Power radius, profiles will tend to flatten for numerical reasons.

However, episodes of strongly-clustered star formation that cycle between dense GMCs and explosive outflow can produce physical cored DM profiles (Pontzen & Governato, 2012; Di Cintio et al., 2014b; Chan et al., 2015; Oñorbe et al., 2015). Hopkins et al. (2018b) considered a “full physics” resolution study of the DM profiles in m10v and m10q here at lower resolution (from $250$ to $10^{5}\,\rm\,M_{\odot}$), and argued that convergence was dominated by the baryonic effects, not traditional $N$-body considerations, especially in the more massive m10q. We confirm this here: at higher resolution, the core in m10q is more pronounced, exactly as expected given its larger stellar mass in that particular run. Given the weak dependence of burstiness on resolution, it is not surprising that the cusp/core behavior also remains robust.

Most importantly, as predicted by the much lower-resolution simulations referenced above, all our UFDs (which all have $M_{\rm\star}/M_{\rm{vir}}<10^{-4}$), including those not shown, exhibit cusps down to at least $50-100\,$pc. For those UFDs that exhibit $\sim 100\rm\,pc$ cores, the deviation from NFW occurs just outside of the converged radius, and is likely in part due to imperfect centering on the dark matter. Furthermore, the half-light radii of the UFDs are all substantially larger than their numerical convergence radii, further indicating that measurable dynamics of such UFDs (which are sensitive to the total mass within the deprojected half-light radius; see below and Wolf et al. 2010) would point to dark matter cusps. However, it is important to note that the lack of cores predicted in UFDs is resolution-limited, and allows for the presence of small $\lesssim 100\rm\,pc$ cores to be detected observationally without posing a challenge to the model.

3.4 Kinematics and Galactic Structure

3.4.1 Sizes and Surface Brightness Distributions

Our result that every isolated dark matter halo with $M_{\rm{vir}}>4.4\times 10^{8}\rm\,M_{\odot}$ forms a UFD suggests that these objects may be ubiquitous in the field. However, these low-mass halos have shallow potential gravitational wells, causing the galaxies that form within them to have larger effective radii and extremely low surface brightnesses (Kaufmann, Wheeler & Bullock, 2007; Bullock et al., 2010; Bovill & Ricotti, 2011a, b, W15). This means that, despite their abundance, they may be very difficult to detect. W15 calculated values for $R_{1/2}$, the 2D projected half-mass radius, for their lowest mass UFDs ($M_{\rm\star}\lesssim 10^{4}\rm\,M_{\odot}$) that ranged from $200-500\rm\,pc$, making them undetectable with current surveys such as the Sloan Digital Sky Survey (SDSS, York et al., 2000). However, because the galaxies studied had only tens of star particles, and the size of galaxies is highly sensitive to resolution, there was a possibility that the $R_{1/2}$ of the UFDs was not resolved.

The left panel of Fig. 7 shows the mass-size relation for our much better-resolved dwarfs alongside data for classical MW dwarfs (McConnachie 2012, open black circles), year 1 and 2 UFD candidates galaxies from the Dark Energy Survey (DES, Drlica-Wagner et al., 2015; Bechtol et al., 2015), Hydra II (Martin et al., 2015), Virgo I (Homma et al., 2016), and Cetus III (Homma et al., 2018); for the last two galaxies, we used their stellar mass estimates under the assumption of a Kroupa IMF in order to be consistent with the simulations. Also shown in the figure is a surface brightness detection limit of $\rm\mu_{V}=30~{}mag~{}arcsec^{-2}$ for solar absolute magnitude $\rm M_{\odot,V}=4.83$ assuming a stellar mass-to-light ratio of $M_{\rm\star}/L\approx 1\,\rm\,M_{\odot}/L_{\odot}$ (so this corresponds to a physical, bolometric $0.036\,L_{\sun}\,{\rm pc}^{-2}$). For, e.g., a Plummer profile with central surface brightness $\Sigma_{\rm peak}=L/\pi\,R_{1/2}^{2}$, this corresponds to the surface brightness detection limit for SDSS. We also compare the improved limit $\rm 32.5~{}mag~{}arcsec^{-2}$ which is anticipated for upcoming surveys such as the co-added Large Synoptic Survey Telescope (LSST).

Our more massive galaxies agree well with observed systems in Fig. 7, while the lowest mass UFDs have even lower surface brightnesses than those of W15 in the same mass range. Every UFD with $M_{\rm\star}<10^{4}\rm\,M_{\odot}$ in our sample has $R_{1/2}>400$, lying close to or above the likely LSST detection limit, which suggests that they may well go undetected for some time. The sizes of these are much larger than our force softening or the mass-resolution-based Power et al. (2003) convergence radius discussed above ($<30\,$pc), so this is likely robust. Moreover, where we do see (small) changes in size with resolution, the galaxies essentially move along close-to-constant surface brightness tracks.

Interestingly, there does not seem to be any tight correlation in the simulated galaxies between $R_{1/2}$ and $M_{\rm\star}$, at $M_{\ast}\ll 10^{7}\,M_{\sun}$. This is yet another indication that galaxy formation prescriptions with simplistic recipes for determining galaxy size, or the assumption that $R_{1/2}\propto r_{\rm{vir}}$ (as in e.g. Mo, Mao & White, 1998; Kravtsov, 2013), fail for the lowest mass galaxies (see also Kaufmann, Wheeler & Bullock 2007). The recent detection of ultra-diffuse galaxies (UDGs) likewise suggests that even more massive dwarfs have wildly varying effective radii (several dex range) for galaxies at a single $M_{\ast}$ (van Dokkum et al., 2015). The apparent size-mass relation in the observations (Fig. 7) is clearly an effect of the surface brightness limit (where most of the galaxies pile up, as noted by the authors of the survey studies cited above). Crucially, it is likely that most UFDs lie at lower surface brightness (larger radii at fixed stellar mass) than we currently can detect. The discovery of a new Milky Way satellite in the Gaia Data Release 2 with $\rm\mu_{V}=32.3~{}mag~{}arcsec^{-2}$ – 100 times more diffuse than most UDGs – is likely the first indication of what is to come in the extremely low surface brightness sky (Torrealba et al., 2018).

It is noteworthy that our simulations do not reproduce those UFDs with half-light radii in the $100\rm\,pc$, given how well they agree with many other properties of observed dwarfs. To check whether this is also a selection effect, we restrict those dwarfs in our sample that have average surface brightnesses $\rm>30~{}mag~{}arcsec^{-2}$ to an annulus that encloses only the central region of the galaxy with at least that surface brightness. Doing so reduces both $R_{1/2}$ as well as $M_{\rm\star}$, which yields $R_{1/2}\sim 100\rm\,pc$ ($R_{1/2}<200\rm\,pc$) for $\rm m09_{30}$ ($\rm m10v_{30}\,B$) – bringing both into the region occupied by the DES UFDs. This suggests that observations may only be sensitive to the “bright” core of more massive objects and raises an intriguing question: do the DES UFDs all have diffuse (and relatively massive) outer halos that are currently invisible to us?

3.4.2 Dynamical Masses and Sizes

Although our lowest-mass simulated galaxies show no overlap with observations in the stellar mass – size plane, there is significant overlap with observations in the dynamical mass – size plane. The right panel of Fig. 7 shows $R_{1/2}$ vs.$M_{1/2}$, where $M_{1/2}\simeq 4\,\langle\sigma_{\rm los}^{2}\rangle\,R_{1/2}/G$ $M_{1/2}$ is an estimate of the total (dynamical) mass from $\sigma_{\rm los}$, the line-of-sight stellar velocity dispersion (see, e.g., Wolf et al., 2010). We compare observed systems from Collins et al. (2013) with $R_{1/2}$ taken from Kirby et al. (2014) where available, and McConnachie (2012) otherwise. The simulations agree well with the data at the radii/masses both sample (see also Campbell et al. 2017; González-Samaniego et al. 2017; Errani, Peñarrubia & Walker 2018), and with the true mass within the 3D half-stellar-mass-radius (endpoints of colored bars).

The overlap between all simulated and observed galaxies in the right-hand panel of Fig. 7 means that the very diffuse galaxies in the left-hand panel, which would be invisible to current observational surveys, do not inhabit dark halos that are fundamentally different from those hosting analogs of observed systems. The truly remarkable aspect of these galaxies is their large half-light radii; the underlying mass distributions, which are strongly dominated by dark matter at all radii, look just like standard predictions of CDM. For $r\ll r_{\rm s}$, where $r_{\rm s}$ is the NFW scale radius, we expect that $M_{\rm enclosed}(<\!r)\propto r^{2}$; this is precisely the correlation seen in Figure 7.

3.4.3 Rotation

Observations of both dSphs and dIrrs in the Local Volume suggest that most dwarf galaxies have kinematics dominated by random motions (Wheeler et al., 2017, hereafter W17). This can be seen in Fig. 8 for an observational sample of thirty dSphs and dwarf ellipticals and ten isolated dwarfs in the LG taken from W17. Following their prescription to calculate $v_{\rm{rot}}/\sigma$ for our simulations, where $v_{\rm{rot}}$ is the rotation across its axis and $\sigma$ is the underlying (constant) line-of-sight velocity dispersion, we perform a Bayesian analysis on the positions and velocities of each simulated galaxy along independent lines of sight. We explore two models for $v_{\rm{rot}}$: a flat model that assumes constant rotation, $v_{\rm{rot}}(R)=v_{o}$, and a radially varying pseudo-isothermal sphere, $v_{\rm{rot}}=v_{o}\sqrt{1-R_{o}/R\arctan(R/R_{o})}$, where $R$ is the distance from the rotation axis on the plane of the sky and $v_{o}$ and $R_{o}$ are the rotation velocity and rotation radial scale parameters respectively (see W17 for details).

Fig. 8 shows $v_{\rm{rot}}/\sigma$ vs ellipticity for all simulations (with three independent viewing angles for each), as well as data from W17. The solid line is the locus occupied by oblate isotropic rotators (galaxies flattened primarily by rotation; see Binney 1978). Objects below the line are generally understood to be pressure-supported. Seven of the ten isolated dwarfs (black stars) show no clear signs of rotation; the same is true for nearly all of the LG satellites (grey downward triangles). In addition to posterior probabilities for $v_{\rm{rot}}$ and $\sigma$, the model calculates the Bayesian evidence for both rotation vs non-rotation and flat model vs. radially varying models. Objects that prefer a flat model are shown as open symbols for the isolated observed dwarfs and the simulations. None of the observed satellites prefer a flat model. Fig. 8 also shows $v_{\rm{rot}}/\sigma$ vs. $M_{\ast}$.

Our simulated galaxies show a remarkably good overlap with the observed data in both of these planes. With higher resolution and the newer FIRE-2 code, we find marginally more rotation than the FIRE-1 simulations from W17, but the difference is small. Although all three of the more massive dwarfs do have Bayesian evidence for rotation, none have high $v_{\rm{rot}}/\sigma$. At high masses ($\gtrsim 10^{6}\,M_{\sun}$), our sample has just one galaxy (m10q) that has low $v_{\rm{rot}}/\sigma$. Since most observed systems are also non-rotating at this mass, this is expected, but large $v_{\rm{rot}}/\sigma\sim 2$ begins to appear at these masses; it therefore would be interesting to explore this mass regime with better statistics.

3.5 Chemical Abundances

Kirby et al. (2013a) used stellar metallicities to demonstrate that the relationship between stellar mass and stellar metallicity extends, unbroken, down to Milky Way and M31 dwarf spheroidals, irregulars, and even UFDs down to $M_{\rm\star}\sim 10^{3.5}\rm\,M_{\odot}$. This is a striking result, as the sample includes both satellites and isolated dwarfs, meaning the relationship is unaffected by infall into a more massive host. Fig. 9 explores the extremely low-mass end of the (stellar) mass - (stellar) metallicity relation (MZR), comparing observations from Kirby et al. (2013a) and Vargas, Geha & Tollerud (2014) for MW+Andromeda satellites. There is an obvious and intriguing discrepancy between observed and simulated galaxies, which increases to nearly 2 dex towards lower masses.

We emphasize that at higher masses, $\gg 10^{7}\,M_{\odot}$, previous studies with the same FIRE physics have extensively compared stellar and gas-phase MZRs at both $z=0$ and higher redshifts and found remarkably good agreement between observations and simulations in $\rm[Fe/H]$ as well as other species ([O/H], [Mg/H], [Z/H], etc.; see Ma et al. 2017a, b; Escala et al. 2018; Hopkins et al. 2018b; El-Badry, Weisz & Quataert 2017; Wetzel et al. 2016). The discrepancy reported here appears to be specific to low-mass dwarfs. We also emphasize that the discrepancy is not unique to $\rm[Fe/H]$, as [Mg/H] and other species show a similar (albeit slightly weaker) offset, nor to the method by which $\rm[Fe/H]$ is weighted (e.g., light-or-mass weighting or taking $\langle$[Fe/H]$\rangle$ vs. [$\langle$Fe$\rangle$/$\langle$H$\rangle$] give a qualitatively similar result).

What could cause this? Recall, our simulations adopt a very simple IMF-sampling and yield model: even though SNe are discrete, the yields are IMF-averaged, and most yields are independent of the progenitor metallicity (see Hopkins et al. 2018b for details); they therefore reflect yields at solar metallicity, where they are best understood both theoretically and observationally. Yield predictions for extremely metal-poor stars in the literature differ by as much as $\sim 1$ dex (Woosley & Weaver, 1995; François et al., 2004), so a strong dependence on progenitor metallicity could explain the offset seen here. It is also possible that the discrepancy has its origin in differential re-incorporation of metals (vs. ejection in galactic winds): our dwarfs produce sufficient metals to lie on the MZR, if they retained them all and re-incorporated them into new stars in a closed-box fashion. But this is difficult to reconcile with the strong outflows required to explain their very low stellar masses (Kirby, Martin & Finlator, 2011).

Perhaps more importantly, our simulations include no explicit model for hyper metal-poor or metal-free Population-III stars (this is the reason why they are initialized with a metallicity “floor” of $\rm[Fe/H]=-4$). Theoretical calculations suggest a single massive Pop-III star exploding in a pair-instability supernova could produce $\gtrsim 100\,M_{\sun}$ in heavy elements (Kozyreva, Yoon & Langer, 2014) – sufficient to enrich $\gtrsim 3\times 10^{6}\,M_{\sun}$ worth of gas (comparable to the mass of our most massive galaxy studied here) to the minimum metallicity $\rm[Fe/H]\sim-2.5$ observed. Indeed, a number of recent studies have argued that an early Pop-III phase should pre-enrich almost all star-forming galaxies in halos more massive than $10^{8}\,M_{\sun}$ to $\rm[Fe/H]\sim-2$, even by $z\gtrsim 15$ (Chen, 2015; Jaacks et al., 2018). Since these would leave no other relics, they would not change our other predictions, but presumably they would leave unique, testable abundance patterns in the observed dwarfs.

Another interesting possibility is that pre-enrichment occurs not via Pop-III stars but environmentally (from the more massive MW/Andromeda host), as the observed dwarfs are all LG satellites. Given typical progenitor masses of MW+Andromeda-mass systems at $z\gtrsim 6$, if just $\sim 1\%$ of the metals produced by SNe at these times escapes to $\sim 100\,$kpc physical radii, then a $\sim 1\,$Mpc co-moving volume would be polluted to $\rm[Fe/H]\sim-2$. To test this, we consider the simulations (using the identical code and physics) of LG-mass, MW+Andromeda-like pairs, presented in Garrison-Kimmel et al. (2018). While the resolution of these simulations ($\sim 3000\,M_{\sun}$) is extremely high for such massive halos, it is not sufficient to resolve the lowest-mass UFDs. However, we can directly compare $\rm[Fe/H]$ values for galaxies at $M_{\rm\star}>10^{5}\rm\,M_{\odot}$. We find that there is still a $0.5-1$ dex offset between the simulated satellite galaxies within $300\rm\,kpc$ of the galactic center and the observations, suggesting that pre-enrichment from a massive neighbor is not enough to relieve the discrepancy for massive dwarfs.

To get a sense for how the presence of a host galaxy may affect lower-mass satellites, we consider that all of the galaxies with $M_{\rm\star}<10^{5}\rm\,M_{\odot}$ from Kirby et al. (2013b) are within $150\rm\,kpc$ of the Galactic center. If we measure the $\rm[Fe/H]$ values for the gas enclosed within the radius occupied by $90\%$ of those stars traced back to higher redshift, while excluding the inner $10\%$ to eliminate the galactic disk, we determine that $\sim 20\%$ ($1\%$) of the gas has been enriched to $\rm[Fe/H]\sim-2.5$ by $z=4$ ($z=10$). This means that, for the more massive UFDs that are still forming stars to $z\sim 2$, pre-enrichment from nearby massive host could explain the discrepancy. Additionally, preliminary examination of an extremely-high resolution ($\sim 900\,M_{\sun}$) version of the MW-mass halo studied in Wetzel et al. (2016), run to $z\sim 4$, shows several satellites pre-enriched to $\rm[Fe/H]\sim-2.5$ by $z=5$. However, for the lowest mass UFDs that have completed star formation by $z=10$, this mechanism alone is insufficient to explain the discrepancy.

4 Conclusions

We have run a suite of hydrodynamic cosmological zoom-in simulations with explicit treatment of star formation and stellar feedback from SNe, stellar mass-loss, and radiation to $z=0$. These simulations have baryonic mass resolution of $30\,M_{\sun}$ (maximum resolved cold gas densities $\gtrsim 10^{5-6}\,M_{\sun}$, spatial scales $\sim 0.1-0.4$ pc, timescales $\sim 10-100$ yr), allowing us to probe smaller physical scales than previously possible in cosmological simulations. They resolve ultra-faint dwarf galaxies ($M_{\rm\star}<10^{5}\rm\,M_{\odot}$) with $\gtrsim 100-1000$ star particles and, for the first time, the internal structure of of individual SNe remnants within the cooling radius. Using these simulations, we have shown that:

• •

“Well-resolved” galaxies ($\gtrsim 100$ star particles) exist in all DM halos with $M_{\rm{halo}}>5\times 10^{8}\rm\,M_{\odot}$ or $V_{\rm max}>13~{}{\rm km\,s}^{-1}$. This suggests these are ubiquitous in the field. Using DM-only LG simulations, we estimate there are anywhere from $\sim 180-380$ such UFDs ($M_{\rm\star}\sim 10^{3-5}\rm\,M_{\odot}$) “isolated” in the Local Field, up to another $\sim 140$ as satellites of other Local Field dwarfs, and $\sim 400-610$ within the virial radii of the MW and Andromeda. However, we also predict most UFDs have very low surface brightness, below the detection capabilities of current surveys. It may be more efficient to search for UFDs as satellites of isolated LG field dwarfs: in the simulations, these have properties that may render them detectable by fully co-added LSST data.

• •

All UFDs have uniformly ancient stellar populations: they complete star formation by $z>2$, even if they are isolated field centrals, owing to reionization quenching. More massive field dwarfs continue SF to $z<0.5$.

• •

Dwarfs in this mass range have bursty and strongly clustered SF in both time and space, which produces violent outflows. Despite this, only dwarfs with $M_{\rm\star}/M_{\rm{vir}}\gtrsim 10^{-4}$ produce appreciable cores in their DM profiles; all the UFDs are below this threshold and exhibit DM cusps down to $\sim 100\,$pc or smaller (equivalent to $\sim 0.2\,R_{1/2}$). This is expected from simple energetic arguments: the UFDs simply produce insufficient stellar feedback to strongly perturb their DM potential. Neither the burstiness nor cusp/core threshold is strongly sensitive to resolution.

• •

Properties that can be measured in lower-resolution simulations (e.g., masses, sizes, cusp/core profiles, SFHs of the more massive galaxies here) appear robust to resolution, given the FIRE-2 numerical method adopted here. This is particularly important since many studies have argued the evolution of SNe blastwaves, bubble overlap, and galactic outflows becomes much less strongly-sensitive to sub-grid numerical implementation choices at a mass resolution $<100\,M_{\sun}$ (owing to the ability to resolve individual SN blast waves self-consistently).

• •

Where the surface brightness is high enough to be detected, the simulations agree well with the location of observed galaxies in size-stellar mass-dynamical mass space. While there is a physical correlation between stellar size ($R_{1/2}$) and dynamical mass, driven by sampling a larger fraction of the halo mass at larger radii, there is no intrinsic tight correlation between stellar mass and size at these UFD masses. Rather, the observed relation reflects the surface-brightness limits of current surveys.

• •

The simulations have $v_{\rm{rot}}/\sigma$ and ellipticity values consistent with most observed LG isolated and observed dwarfs. Both simulated and observed systems are primarily non-rotating, even in star-forming, isolated field galaxies. Hints of discrepancies between rotation in observed and simulated FIRE dwarfs suggested in El-Badry et al. (2017, 2018) appear to manifest only at higher masses ($M_{\ast}\sim 10^{7-8}\,M_{\sun}$).

• •

Our dwarfs appear to under-predict the metallicities of observed LG satellites, with the discrepancy growing below $M_{\rm\star}\ll 10^{7}\,M_{\sun}$ to just under 2 dex for UFDs. We argue that some of this may owe to pre-enrichment in the LG environment by the massive host, suggesting the observed LG dwarfs may not be universally representative. It may also reflect our neglect of any treatment of Pop-III stars (a single one of which could, in theory, produce more than enough metals to account for the discrepancy), or more generally the effects of more detailed progenitor mass and metallicity-dependent yields.

Overall, our study confirms many results previously found in lower-resolution simulations in the $M_{\rm\star}\sim 10^{3}-10^{6}M_{\rm\star}$ mass scale. Given that these simulations are beginning to resolve the evolution of individual supernovae, this suggests that the sub-grid approximations used in the lower-resolution simulations may be robust, at least insofar as their consequences for galaxy properties. The major uncertainties at the resolution level simulated here may no longer lie in sub-grid treatments of the collective effects of stellar feedback (which can be explicitly resolved) but rather in the fact that different individual stars have different evolutionary histories and feedback properties. One hopes that this may signal a turning point in simulations of dwarfs where we move from cosmological simulations with effective or sub-grid ISM physics models to those that directly model collapse and fragmentation of molecular clouds into individual stars, akin to previous studies of first stars or the rich studies of individual patches of the ISM (which actually often contain more mass than a UFD).

Acknowledgments

CW is supported by the Lee A. DuBridge Postdoctoral Scholarship in Astrophysics. Support for PFH and SGK was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342, and NASA grants NNX15AT06G, JPL 1589742, 17-ATP17-0214. ABP acknowledges generous support from the George P. and Cynthia Woods Institute for Fundamental Physics and Astronomy at Texas A&M University. SGK acknowledges additional support by NASA through Einstein Postdoctoral Fellowship grant number PF5-160136 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. 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 the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522, by NASA through grants NNX15AB22G and 17-ATP17-0067, and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592.