The In Situ Origins of Dwarf Stellar Outskirts in FIRE-2

Due to the low masses and small sizes of dwarf galaxies, studying their detailed structure requires high-resolution simulations. Moreover, their shallow potential wells make the structure of dwarf galaxies relatively more sensitive to the details of the implementation of feedback from star formation, such as stellar winds and supernovae (SNe), than the structure of more massive galaxies. Thus, while realistic feedback prescriptions are crucial to understanding the structural properties of low-mass galaxies, such galaxies are also among the most sensitive tests of these same feedback prescriptions (e.g., Brooks & Zolotov 2014; Hu et al. 2016; Wheeler et al. 2017; Hu 2019; Wheeler et al. 2019; Dashyan & Dubois 2020; Smith et al. 2021). Furthermore, cosmological simulations are necessary to properly capture the significant effect of reionization on these small galaxies (e.g., Bovill & Ricotti 2009; Boylan-Kolchin et al. 2015; Fitts et al. 2017; Weisz & Boylan-Kolchin 2017; Graus et al. 2019). Cosmological simulations also supply the variety of environments (e.g., Garrison-Kimmel et al. 2019; Jahn et al. 2019) and the assembly histories (e.g., Fitts et al. 2017) needed to study the dwarf population, especially given the relatively large impact that the interactions even with purely dark subhalos can have on the star formation histories (SFH) of dwarfs (Starkenburg et al. 2016).

We use simulations from the FIRE project, specifically the suite of cosmological-baryonic zoom-in simulations run with the FIRE-2 feedback models (Hopkins et al. 2018) along with some variations described below. The FIRE project contains a suite of simulated dwarf galaxies that are uniquely well suited for the structural study at hand—as we have argued above, assessing the origin of dwarf stellar structure formation requires simulations that are both cosmological in nature and high resolution. The baryonic particle mass of the dwarf simulations at hand range between 2100 and 7100 M_⊙, with a minimum gravitational force softening of ∼2 pc. Though there are a relatively low number of initial conditions (we study nine in this work), the resolution of the FIRE-2 simulations is critical for robust studies of the stellar structure of low-mass galaxies (Ludlow et al. 2021). Indeed, it has been shown that the stellar structure of FIRE-2 dwarfs at even lower masses remains robust to the resolution effects (Wheeler et al. 2019). Thus, we expect the structure of these relatively more massive galaxies to be robust to resolution effects due to the high resolution of the FIRE-2 simulations—moreover, as a small forward reference, we note that we find no evidence for a dependence between the resolution and the dwarf properties that we measure in this work.

All simulations use the Meshless-Finite-Mass mode of the GIZMO ¹² gravity+magnetohydrodynamic code (Hopkins 2015), which provides adaptive spatial resolution, conservation of mass, energy, momentum, and excellent shock-capturing and conservation of angular momentum, reproducing the advantages of both smoothed-particle hydrodynamics and Eulerian adaptive mesh refinement schemes. 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).

Table 1 lists the simulations included in this study. All are simulations of isolated galaxies (the main halo being simulated is the highest-mass object in the zoomed box). The names of the simulations follow a convention: (1) in the prefix "mXX," the number XX denotes the order of magnitude of the halo mass of the most massive halo in the zoomed region, in solar masses; (2) the letter following the prefix is a label representing the simulation's unique initial conditions; (3) the parenthetical after the resolution indicator is a shorthand reference to the physics included in the run (described further below).

We select the z = 0 snapshots of the set of isolated dwarf galaxy simulations with stellar masses 10⁸ ≲ M_⋆/M_⊙ < 10¹⁰, along with a set of three higher-mass galaxies (M_⋆ > 10¹⁰ M_⊙) that we use as a reference sample. We will show in this work that the stellar halos of the high-mass reference set are built largely from accreted particles, as expected for L_⋆ galaxies, and use this high-mass sample as a benchmark for the classical stellar halo assembly.

We use simulations run with three different variations of the FIRE physics engine, denoted in the "Run Type" column of Table 1 as (Hydro+, no MD), (Hydro+, MD), (MHD+), and (CR+) to explore the sensitivity of our conclusions to small changes in the baryonic physics implemented in the simulations. Runs marked (Hydro+, no MD) and (Hydro+, MD) use the core FIRE-2 physics described in full in Hopkins et al. (2018) with and without metal diffusion (MD), respectively. Runs marked (Hydro+, MD) in Table 1 use the physics described above, but also include the numerical implementation of turbulent metal diffusion described in Escala et al. (2018). Including a model for turbulent metal diffusion allows neighboring gas cells to exchange gas as it is enriched, as we presume happens in nature. This has the main effect of narrowing the distribution of abundances of different metals in the gas and stars, bringing the width of the abundance distributions in line with observations from dwarf galaxies without changing the structural properties of the galaxies themselves (Escala et al. 2018). Since we do not expect metal diffusion to affect the properties studied in this work, both (Hydro+) designations may be treated interchangeably, as separate realizations of similar dwarfs from the same initial conditions that differ only in the stochastic variation induced by a different pattern of stellar feedback—we propagate the presence of metal diffusion in the naming scheme for completeness.

In brief, radiative heating and cooling is treated from 10–10¹⁰ K, including free–free, photoionization/recombination, Compton, photoelectric and dust collisional, cosmic ray, molecular, and metal-line and fine-structure processes (following each of 11 tracked species independently), and accounting for photoheating both by a UV background (Faucher-Giguère et al. 2009) and local sources, as well as self-shielding. In the UV background model used for the simulations in this work, reionization occurs at z ∼ 10, significantly earlier than the current empirical constraints (see, e.g., Planck Collaboration et al. 2020). The early reionization onset mainly affects the SFH of ultra-faint galaxies (up to around M_* ∼ 10⁵ M_⊙; Wheeler et al. 2019) with far less stellar mass than the simulated galaxies studied here.

Star formation occurs only in gas identified as self-gravitating according to the Hopkins et al. (2013) criterion, which is also molecular and self-shielding (following Krumholz & Gnedin 2011), Jeans unstable, and exceeds a minimum density threshold ${n}_{\min }=1000\,{\mathrm{cm}}^{-3}$. Once a star particle forms, the simulations explicitly follow several different stellar feedback mechanisms, including (1) local and long-range momentum flux from radiation pressure (in the initial UV/optical single-scattering, and reradiated light in the IR), (2) energy, momentum, mass, and metal injection from SNe (Types Ia and II) and stellar mass loss (both OB and AGB), and (3) photoionization and photoelectric heating. Every star particle is treated as a single stellar population with known mass, age, and metallicity from which all feedback event rates, luminosities, energies, mass-loss rates, and other relevant quantities are tabulated directly from stellar evolution models (starburst99; Leitherer et al. 1999), assuming a Kroupa (2001) initial mass function (IMF).

Runs marked (MHD+) again use the same physics as the (Hydro+) runs but solve the equations of ideal magnetohydrodynamics (MHD) as described and tested in Hopkins & Raives (2016) and Hopkins (2016), with anisotropic Spitzer–Braginskii conduction and viscosity as described in Hopkins (2017), Su et al. (2017), and Hopkins et al. (2020a). Runs marked (CR+) include all of the physics implemented in the (MHD+) runs, with the addition of the magnetohydrodynamic treatment of cosmic rays described in Chan et al. (2019), Hopkins et al. (2020a), and Ji et al. (2020). Magnetic fields, conduction, and viscosity have only small effects on bulk galaxy properties at any mass scale (Hopkins et al. 2020a), but cosmic rays can can suppress SF and total stellar masses by factors ∼2–4 for the galaxies with M_vir ≳ 10¹¹ M_⊙ at late times (z ≲ 1–2) by building up hot gas in the galaxy halo that suppresses the galactic fountain (Hopkins et al. 2020a). This only occurs in halos massive enough to confine the hot (T ∼ 10⁶ K) gas, and the set of dwarfs that we consider here are right on the boundary of this effect, which is more pronounced in the more massive systems used for comparison (last five rows of Table 1). However, we do not see any systematic difference in the shapes of the galaxies produced in the (CR+) runs relative to the (MHD+) and (Hydro+) runs.

These simulations are uniquely well suited to study the connection between the formation of dwarf galaxies and their structure. The combination of physics implementations in the (Hydro+) runs has already been shown to reproduce the mass–size relation of observed galaxies across ≳ 5 orders of magnitude in stellar mass (El-Badry et al. 2016; Chan et al. 2018). The sample includes a large range of different assembly histories (9 in total for the dwarf sample) that give rise to a wide variety of present-day galaxy morphologies, from thin disks to highly diffuse structures (El-Badry et al. 2018b). With the FIRE-2 feedback model, the resolved ISM produces SFH that are generically bursty for dwarf galaxies (e.g., Sparre et al. 2017; Faucher-Giguère 2018; Flores Velázquez et al. 2021), an effect that can be suppressed in lower-resolution simulations that do not resolve the multiphase ISM. On the other hand, the burstiness of the star formation in these simulations is slightly enhanced relative to the models that explicitly include local radiative feedback, which has the effect of warming the gas and disrupting some giant molecular clouds (GMCs) before collapse, leading to a smoother SFH in massive dwarfs (Hopkins et al. 2020b).

Because these simulations resolve the individual large GMCs (Benincasa et al. 2020) where the star formation occurs, the star formation has a very high spatial granularity across the galaxy; clustering of star-forming clouds is thus naturally taken into account down to ∼pc scales. The granularity of the resulting stellar feedback in both space and time is also very high, since the mass of the individual star particles (2100–7100 M_⊙ per particle) is sufficiently small to support the assumption that each represents a single-age, single-metallicity stellar population while still being large enough to fully sample the IMF at the high-mass end (Sanderson et al. 2020). The inclusion of turbulent metal diffusion additionally reproduces the observed dependence of the width of abundance spreads on stellar mass in dwarfs (Escala et al. 2018). Comparing the runs with additional physics to the (Hydro+) runs for identical initial conditions allows us to both confirm that these variations do not strongly influence the structure of dwarfs, and set a bound on the degree of scatter due to stochastic SN feedback for a fixed assembly history.

We parameterize the three-dimensional shapes of the sample using the ratio of the semiprincipal axis diameters A, B, and C where A ≥ B ≥ C. The ratios of these axes, B/A and C/A, give a quantitative description of the intrinsic shape of each galaxy; in Figure 1, we show three extreme cases in both B/A versus C/A space at left and as wireframe renderings at right. In particular, we show a prolate galaxy in green (C/A ∼ B/A ≪ 1), a disky galaxy in blue (C/A ≪ B/A ∼ 1), and a spheroidal galaxy in red (C/A ∼ B/A ∼ 1). For context, we also show the 1σ contours of the observed dwarfs sample of Kado-Fong et al. (2020) in blue (measurements at 1R_eff) and red (measurements at 4R_eff).

Zoom In Zoom Out Reset image size

Making the structural measurements for age-separated stellar components of low-mass halos requires simulations of sufficiently high resolution. Even at the mass resolution of the FIRE-2 simulations studied in this work, we use cumulative (i.e., the shape of the stars within a given radius), rather than differential (the shape of the stars at a given radius), intrinsic shape measurements in order to overcome the instability due to low stellar densities in the outskirts of the dwarfs; though we note that the radial density profile remains stable. This is a qualitatively different approach than what is used in integrated light measurements of intrinsic shape, where the galaxy intrinsic shape distribution is inferred at fixed radius (Padilla & Strauss 2008; van der Wel et al. 2014; Kado-Fong et al. 2020). The need for a cumulative measure of the intrinsic shape will also lead us to measure the intrinsic shapes as a function of stellar age instead of galactocentric radius, comparing in particular the young (ages < 2 Gyr) and old (ages > 12 Gyr) stellar populations (see Section 3.1). Based on observations of nearby dwarfs, this age separation is expected to correlate with the observational differential shapes (see, e.g., Aparicio & Tikhonov 2000; Aparicio et al. 2000; Zaritsky et al. 2000; Hidalgo et al. 2003; Demers et al. 2006; Bernard et al. 2007; Stinson et al. 2009; Strader et al. 2012; Nidever et al. 2019; Pucha et al. 2019). However, we note that the age-separated populations that we use in this work are not the exact equivalent to the radially separated stellar populations of Kado-Fong et al. (2020). The simulated and observed measurements also differ in that the observed measurements of the galaxy-projected shapes are used to infer the intrinsic shape distribution of the galaxy population (whereas the simulated values are direct intrinsic shape measurements of individual galaxies), and that the simulated measurements are mass-weighted in bins of stellar age (whereas the observed measurements are light-weighted at fixed radii). Thus, the emphasis of this work will be on the change in galaxy structure within the simulations and observations as a function of galaxy properties, rather than a direct quantitative comparison between the observed and simulated measurements. However, in Appendix B, we also demonstrate that these cumulative measurements well trace the mock measurements made using the observational inference machinery of Kado-Fong et al. (2020).

As mentioned above, unlike in observations, in the simulations, we are able to directly measure the intrinsic shape of the galaxies in our sample. To do so, we compute the reduced moment of inertia tensor I from the set of star particles at distances less than ${R}_{\max }=50\,\mathrm{kpc}$ for the simulated galaxies with M_halo ∼ 10¹¹ M_⊙ and ${R}_{\max }=70\,\mathrm{kpc}$ for the galaxy with M_halo ∼ 10¹² M_⊙. The reduced moment of inertia tensor I is given by

$$\begin{eqnarray}{\boldsymbol{I}}=\left[\begin{array}{ccc}{I}_{{xx}} & {I}_{{xy}} & {I}_{{xz}}\\ {I}_{{yx}} & {I}_{{yy}} & {I}_{{yz}}\\ {I}_{{zx}} & {I}_{{zy}} & {I}_{{zz}}\end{array}\right]\end{eqnarray} \tag{ 1 }$$

where each element is computed from the star particles (of mass m_k ) as

$$\begin{eqnarray}&&{I}_{{ij}}=\displaystyle \frac{{\sum }_{{r}_{k}\lt {R}_{\max }}{m}_{k}{q}_{i,k}{q}_{j,k}/{r}_{k}}{{\sum }_{{r}_{n}\lt {R}_{\max }}{m}_{n}},\end{eqnarray} \tag{ 2 }$$

for i, j ∈ (x, y, z) and where ${{\boldsymbol{r}}}_{{\boldsymbol{k}}}=({q}_{x}^{k},{q}_{y}^{k},{q}_{z}^{k})$ is the distance from the principal halo center. The principle axes calculated from the reduced moment of inertia may overweight the star particles at large distances relative to the unreduced moment of inertia tensor; for our sample, however, we find a median fractional difference in the principal axis ratios of 5%, and conclude that our choice of moment of inertia tensor does not significantly affect our results. The eigenvalues λ₁, λ₂, and λ₃ of this matrix are the inverse squares of the semiprincipal axis diameters A, B, and C of the Poinsot ellipsoid that corresponds to the moment of inertia tensor:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{B}{A} & = & \sqrt{\displaystyle \frac{{\lambda }_{1}}{{\lambda }_{2}}}\\ \displaystyle \frac{C}{A} & = & \sqrt{\displaystyle \frac{{\lambda }_{1}}{{\lambda }_{3}}}.\end{array}\end{eqnarray} \tag{ 3 }$$

This method is, in principle, sensitive to the choice of R_max and emphasizes the contribution of particles at large r (see Equation (2)). We test the effect of our choice of intrinsic shape measurement method as follows. First, though we choose a larger R_max for the higher-mass halo, we would derive the same shape parameters B/A and C/A to within Δ(B/A) ∼ Δ(C/A) ∼ 0.01 at ${R}_{\max }=50\,\mathrm{kpc}$. Furthermore, although we adopt a fixed value of ${R}_{\max }=50\,\mathrm{kpc}$ for the dwarf galaxies considered in this work, we find that varying R_max does not qualitatively affect our results down to ${R}_{\max }\sim 10\,\mathrm{kpc}$. Finally, we fit two-dimensional Sérsic profiles to projections of each galaxy along its principal axes, and find no evidence for a systematic offset between the Sérsic fits and moment of inertia-derived axis ratios. Thus, although the moment of inertia method potentially overemphasizes the star particles at large distances, we find that our results are not impacted by this bias in practice.

Observational works indicate that dwarfs in this mass range should host a central thick disk (see, e.g., van der Wel et al. 2014) with a round extended stellar component (Kado-Fong et al. 2020) that is populated by intermediate/old stars (see Stinson et al. 2009, and citations therein). If a disk-halo system is indeed present in the FIRE-2 dwarfs, we should be able to detect the same transition from disk to halo when computing the intrinsic shape in bins of stellar age—we demonstrate that this transition is indeed detectable in the mock Hyper Suprime-Cam (HSC) data in Section 3.3. We visualize an example of such a divide in Figure 2, which shows the projected density of the dwarf galaxy m11h (MHD+) in young stars (left column, blue), old stars (middle column, red), and neutral hydrogen (right column).

Zoom In Zoom Out Reset image size

We present the main results of this section in Figure 3. In each panel, we show the intrinsic axis ratios of the simulated galaxies, B/A and C/A, where C ≤ B ≤ A. At left, for reference, we show the well-defined stellar disk and halo systems of the high-mass (M_⊙ > 10¹⁰ M_⋆) reference sample (m11f, m11g, and m12i, see Table 1). The intrinsic shapes of young star particles (youngest 25th percentile) and old star particles (oldest 25th percentile) are shown in blue and red, respectively. In this panel, the young stars clearly are assembled in a well-formed disk (C < B ∼ A), and the old stars occupy a round stellar halo (A ∼ B ∼ C). Similarly, in the middle panel, we show the dwarf galaxies in our sample that have a clear disk-halo system. These galaxies are characterized by substantial increases in C/A as a function of age along with relatively little evolution in B/A compared to C/A over the same comparison, and were categorized by visual inspection of the galaxies' stellar components. Although the oldest (ages greater than 12 Gyr) stars may not dominate the light in the outskirts of the FIRE-2 dwarfs, we see the same trend toward a more spheroidal shape when comparing the young (ages less than 2 Gyr) and intermediate (ages between 2 and 12 Gyr) star particles.

Zoom In Zoom Out Reset image size

We note that the young stellar disks in these dwarfs tend to be thicker (higher C/A) than the disks in the higher-mass galaxies—this phenomenon is in good agreement with previous observational findings (Padilla & Strauss 2008; Sánchez-Janssen et al. 2010; Kado-Fong et al. 2020). Both the high-mass and dwarf disks are largely axisymmetric, maintaining a B/A ∼ 1 in both their young and old stellar populations. Though the old stellar population is significantly rounder than the young stellar disks of these dwarfs, the dwarf stellar halos are not as round (meaning C/A approaching unity) as those of the old stellar population in the high-mass comparison sample. This is again consistent with the observations of Kado-Fong et al. (2020), who find that dwarfs have an average minor-to-major axis ratio 〈C/A〉 ∼ 0.5 at 4R_eff, significantly flatter than the outer stellar halo of the Milky Way (MW; C/A ≳ 0.8, Das & Binney 2016; Iorio et al. 2018). This may be due to a difference in the formation mechanism of dwarf stellar halos and massive galaxy stellar halos, which we will explore in later sections.

The evolution of C/A as a function of age is also markedly different between the disky dwarfs and high-mass reference sample. In Figure 4, we show the change in C/A from the minimum C/A for each disky dwarf (purple) and high-mass galaxy (gray). Because C/A increases monotonically as a function of age for both of these groups, this is equivalent to showing the change in C/A from the C/A of the youngest (t < 2 Gyr) star particles. The x–z projections as a function of age for one such disky dwarf (top) and one high-mass galaxy (bottom) are shown at right. The figure serves to illustrate two main points: first, the overall change in C/A is systematically smaller for the disky dwarf sample; second, the high-mass galaxies approach their maximum C/A values more rapidly as a function of age, with a sharp increase in C/A at 5 ≲ t_age ≲ 8 Gyr and relatively little change beyond that. This can also be seen qualitatively in the stark change in shape between the [4, 6] Gyr and [6, 8] Gyr surface densities of the high-mass galaxy. The dwarf galaxies, meanwhile, increase in C/A gradually as a function of age.

Zoom In Zoom Out Reset image size

This difference in the three-dimensional shape as a function of age between the massive and dwarf disk galaxies can be understand when one considers the origins of the star particles that constitute these age-binned groups. In Figure 5, we show the histories of a random subset of star particles in two massive disks (top two rows, gray panels) and two disky dwarfs (bottom two rows, purple panels). In the left panels, we show the particle distance from the galaxy center as a function of time, while on the right we show the distribution of the maximum distance attained by these particles in three age bins (ages < 2 Gyr in blue, between 6 and 8 Gyr in orange, and >12 Gyr in red). The particle tracks on the left are colored by age, as are the histograms at right. We find in particular that the old, spheroidal components of the massive galaxies (the red tracks and histograms show particles with ages >12 Gyr) correspond to accreted components on orbits characterized by large apocenters that persist down to z = 0. The young (ages < 2 Gyr) components are formed in situ close to the host center, and the intermediate (ages between 6 and 8 Gyr) components represent a transition between the ex situ and in situ-dominated stellar populations.

Zoom In Zoom Out Reset image size

The disky dwarfs, meanwhile, are dominated by stars formed in situ (though some accreted particles are present) at all ages—in the right panels, we show that, although the majority of old stars in the massive disks attain a maximum distance greater than the radius containing 90% of the star particles at z = 0 (R_⋆,90(z = 0)), the star particles in dwarf disks are characterized by maximum distances less than R_⋆,90(z = 0) at all ages. That is, the difference in the evolution of three-dimensional shapes as a function of age between the massive and dwarf disks is driven by a difference in the origin of the stars that populate each age group. The spheroidal old component of the massive disks is driven by an accreted population that falls in on relatively circular orbits compared to the flatter, in situ population that dominates the old stellar component of the disky dwarfs. We will further examine the origin and shape evolution of the disky dwarfs in Section 3.2 by examining the role of accreted star particles in building the stellar outskirts of the full sample (Figure 6).

Zoom In Zoom Out Reset image size

Finally, in the right panel of Figure 3, we show the intrinsic shapes of the nondisky dwarfs in our sample. Although we refer to these galaxies as "nondisky" dwarfs, this morphological class refers to galaxies that lack a young stellar disk and do not show a monotonic increase in C/A as a function of stellar age. These dwarfs do not have a well-defined disk-halo structure, and appear to be significantly different in intrinsic shape to the population of observed galaxies from Kado-Fong et al. (2020). ¹³ The lack of a well-defined disk in particular may point to an overly vigorous or bursty SFH, preventing the formation of a rotationally supported disk at low redshift. El-Badry et al. (2018b) found in particular that the ability of FIRE-2 galaxies to form a gaseous disk and maintain a fairly quiescent (less bursty) SFH is linked to the accretion of high angular momentum gas at low redshifts, wherein the low-mass FIRE-2 galaxies struggle to build up stores of this high angular momentum gas due to both efficient gas removal via star formation feedback and inefficient cooling of high angular gas in the circumgalactic medium.

While the star formation feedback and/or cooling prescription may lead to a lack of disk-halo structure in some dwarfs, those same prescriptions are capable of creating the dwarf stellar disk-halo systems that are remarkably similar in structure to the observed population. Thus, in this work, we will examine the two facets of the sample separately, examining first the dwarf stellar halo (in the simulated galaxies that do form disks) and then the origin of the nondisky dwarfs.

As a final note, we find that there is no change in the morphological classification of the galaxy as a function of the physics implementation for the majority of the dwarfs in the sample. That is, it appears that the z = 0 stellar intrinsic structure is largely unaffected by the additional physics implemented in the (MHD+) and (CR+) runs as described in Section 2.1. However, we do find that two galaxies, m11i and m11v, are disky in their (MHD+) runs and nondisky in their (CR+) or (Hydro+) runs. We find that the change in morphological classification in m11i is due to its proximity to the boundary between our disky and nondisky classifications—indeed, in the middle panel of Figure 3, the disky m11i (MHD+) galaxy is shown by the pair of points at the lowest B/A value of the dwarf sample. Thus, we conclude that the change in morphological class of m11i is due to the rigid nature of the binary classification that we use in this work. Though the disky and nondisky dwarfs are generally well-separated morphological classes, there do exist edge cases such as m11i in which the classification can change as a result of relatively small perturbations in intrinsic axes space. m11v does undergo a significant shape change between the MHD+ (B/A = 0.90 and C/A = 0.46 for t < 2 Gyr stars) and CR+ (B/A = 0.25 and C/A = 0.13) runs. However, m11v is also experiencing an ongoing interaction that is the likely cause of this shape variation.

In Figure 6, we show data for a selection of galaxies the maximum distance of the star particles from the host halo center since z ∼ 6.2 as a function of the present-day distance. This figure acts to visualize the contributions to the young (top, blue) and old (bottom, red) stellar populations by in situ and ex situ components: in situ star particles flare out from the 1:1 line, while accreted star particles occupy the horizontal bands that correspond to the maximum distance of the progenitor system over the time considered. In the set of panels at left, we show the z = 0 distances versus maximum distance attained from host center for a disky dwarf (m11h (MHD+)) and nondwarf (m11f (MHD+)), while at right we show the fraction of particles that attain a displacement of greater than 15 kpc (i.e., those particles that lie above the dotted line in the left panels) for the full set of disky dwarfs and massive galaxies. The right panel shows this metric for both the old stars (filled points) and for all particles with a z = 0 host distance of greater than 10% R_vir.

As expected, the young stellar populations in all the galaxies shown in Figure 6 are dominated by in situ contributions. However, the difference in how the dwarfs and the MW-mass galaxy have assembled their old stellar populations is apparent. Whereas the old star particles in the dwarfs (and therefore their stellar halos) are dominated by in situ star particles, the old stellar population in the MW-mass galaxy is completely dominated by accreted star particles. We find no significant difference between the accreted halo fraction of the nondisky dwarfs and disky dwarfs; though the sample size is too small to straightforwardly extrapolate this behavior to the general dwarf population.

The diagnosis of the halo star particle origin that we use in this work is not meant to produce exact ex situ fractions, but is instead aimed to identify bulk changes in the halo formation mechanism. We test our ΔR > 15 comoving kpc threshold against a cut in formation distance, which has been used in prior work with the FIRE-2 simulations to track ex situ contributions (Sanderson et al. 2018) and has been shown to produce the accreted fractions in agreement with subhalo tracking from merger trees for MW-like galaxies in FIRE-2 (Necib et al. 2019). We compare these two ex situ flagging methods for the disky dwarf m11h (Hydro+, MD) and the high-mass galaxy m12i (Hydro+, no MD), and find that our ΔR cut produces accreted fractions consistent to within 10% of those based on a cut on the formation distance. We also consider the effect that our choice of a threshold cut at ΔR = 15 kpc plays on our estimated ex situ fractions, and find that a cut of ΔR = 30 kpc (ΔR = 7 kpc) changes the estimated ex situ fractions by no more than 12% (5.7%), and does not affect the trend with stellar mass shown in Figure 6. We also find that this trend holds when we consider all stars in the outskirts, regardless of age (R(z = 0) > 0.1R_vir (open points in Figure 6). Indeed, the divergent behavior seen in Figure 6 between the dwarf and massive galaxy stellar halos is consistent with the supposition that the stellar halos of massive galaxies are dominated by ex situ stars, while the dwarf stellar halos are largely in situ structures.

The thick and in situ nature of the dwarf stellar outskirts in FIRE-2 also draws a clear structural similarity to the thick disks in more massive galaxies. Though the origin, existence, and nature of thick disks in the MW and beyond are still matters of significant debate (see, e.g., Yoachim & Dalcanton 2006; Bovy et al. 2012; Belokurov et al. 2020; Agertz et al. 2021; Park et al. 2021), it is straightforward to make an internal comparison between the dwarf stellar outskirts and massive thick disks within the context of the FIRE-2 simulations. It has been demonstrated that the thickened shape of the thick disks in the more massive FIRE-2 galaxies is caused by star formation characterized by larger vertical scale heights during bursty phases of star formation at lookback times of t_lb ≳ 5 Gyr (Yu et al. 2021). It is thus of interest to examine whether the shape of the FIRE-2 dwarf stellar populations are set by stellar migration or by a change in the star formation configuration as a function of cosmic time. To address this, we show in Figure 7 the difference in the scale height of the star particles in bins of age at z = 0 and at z = z_form, where z_form is the formation redshift (i.e., the lookback time equivalent to the z = 0 age of the star particle population), of star particles in a disky dwarf galaxy (m11h) and a massive galaxy (m12i) as a function of star particle age. We compute the scale heights from an exponential fit to the vertical stellar surface density in age bins spaced by 1 Gyr with widths of 2 Gyr, and exclude the star particles with ages greater than 10 Gyr for the massive galaxy in this figure, as the massive galaxies are expected to be dominated by ex situ stars at this age (see Figure 6), and the majority of the massive FIRE-2 thick disks are assembled at times more recent than 10 Gyr (Yu et al. 2021).

Zoom In Zoom Out Reset image size

We first consider the massive FIRE-2 galaxy; Figure 7 demonstrates that the star particles in this galaxy show a negligible change in vertical scale height between z = z_form and z = 0. This lack of scale height evolution is consistent with the results of Yu et al. (2021), who find that the thick disks in these massive galaxies are formed as thick disks from birth, rather than from stellar displacement after formation. If, however, we turn our attention to the dwarf disk, we see that the z = 0 vertical scale height is consistently higher than the vertical scale height at formation. This evolution is instead consistent with the picture of El-Badry et al. (2016), where the stellar outskirts are indeed formed from the displacement of stars after formation. This implies that the dwarf stellar outskirts in FIRE-2 are distinct in formation mechanism from both the mechanisms that create the massive stellar halos and massive stellar thick disks in the FIRE-2 simulations.

As discussed in Section 2.2, the way in which we measure the intrinsic shape of the FIRE-2 dwarfs is significantly different than the measurement methods applied to the observations in Kado-Fong et al. (2020). Thus, it is of interest to quantify the extent to which the stellar halos around the FIRE-2 dwarfs would be observable in HSC Subaru Strategic Program (HSC-SSP), the imaging survey used by Kado-Fong et al. (2020). This test will help to understand the extent to which the FIRE-2 dwarf stellar halos and observed dwarf stellar outskirts track the same physical structure.

In order to quantify the observability of the FIRE-2 dwarf stellar halos, we generate mock observations of the sample in this work using the Flexible Population Synthesis package (Conroy et al. 2009; Conroy & Gunn 2010) to generate the simple stellar populations over a grid of stellar ages (between 10⁵ and 10^10.3 yr) and metallicities (in the range $-4\lt {\mathrm{log}}_{10}(Z/{Z}_{\odot })\lt 1$). We then compute the stellar mass-to-light ratio in the HSC i band for each model spectrum. Each star particle in the simulated dwarfs is then assigned an HSC i-band mass-to-light ratio based on its age and metallicity via linear interpolation of the model grid; this process allows us to compute the surface brightness and light-weighted age maps for each galaxy in our sample. Following Sanderson et al. (2018), we construct these maps via a simple binning of star particles with a box size of 0.25 kpc on a side.

We show these surface brightness and light-weighted age maps in Figure 8 for a set of galaxies that span the stellar masses and morphologies of the sample at hand. Each panel shows the isophote at the nominal HSC i-band surface brightness limit for galaxy profiles (μ_i = 28.5 mag arcsec⁻²) as a lime-green contour. In order to best quantify the ability of HSC to detect the old, round stellar outskirts of the FIRE-2 galaxies, we project each galaxy along the y-axis. We also compute an effective radius in the HSC i band using this projection via a two-dimensional, single-component Sérsic fit; in each panel, we note the maximum radial distance of the μ_i = 28.5 mag arcsec⁻² isophote in terms of the HSC i-band effective radius, which we will denote as $\max {R}_{\mathrm{HSC},\mathrm{obs}}$.

Zoom In Zoom Out Reset image size

Naturally, we find that the value of $\max {R}_{\mathrm{HSC},\mathrm{obs}}$ increases as a function of stellar mass; adopting the stellar mass bins used in Kado-Fong et al. (2020), we find that the FIRE-2 galaxies with stellar masses of M_⋆ < 10^8.5 M_⊙ have a mean (standard deviation of) $\max {R}_{\mathrm{HSC},\mathrm{obs}}$ of 2.0 (0.2) R_eff,HSC−i, while the galaxies with stellar masses of 10^8.5 < M_⋆ < 10⁹ M_⊙ have a mean (standard deviation of) $\max {R}_{\mathrm{HSC},\mathrm{obs}}$ of 3.6 (0.8) R_eff,HSC−i, and the galaxies with stellar masses of 10⁹ < M_⋆ < 10^9.6 M_⊙ have a mean (standard deviation of) $\max {R}_{\mathrm{HSC},\mathrm{obs}}$ of 4.1 (1.7) R_eff,HSC−i. These values are consistent with the spatial limit of the HSC imaging for the dwarfs in Kado-Fong et al. (2020).

As shown in Figure 8, in the cases where a young stellar disk is formed (i.e., excluding the nondisky dwarfs due to the lack of an analogous population in the observations), the mock HSC observations extend beyond the young stellar disk into the rounder and older stellar outskirts. This finding supports our initial assertion in Section 2.2 that the age-separated intrinsic shape measurements in the simulated galaxies and the radially separated intrinsic shapes inferred from the observed dwarf population are tracing the same physical structure.

Although diagnosing, in detail, the formation path of nondisky dwarfs is not the main aim of this paper, it is informative to compare their global properties to those of the disky FIRE-2 dwarfs and to observed dwarf galaxies.

It is known that the dwarfs in FIRE-2 are, on average, somewhat lower in their H i gas fractions as compared to observed dwarfs (El-Badry et al. 2018a). Here we expand upon this trend to demonstrate that the gas content of the FIRE-2 dwarfs is also linked to their three-dimensional structure. To do so, we compare the H i gas content ¹⁴ of the FIRE-2 simulations to that of the observed sample of Bradford et al. (2015). We compute the H i and stellar masses within 200 kpc of the galaxy center to derive the H i gas fractions of the FIRE-2 galaxies (though we note that a choice of 50 kpc from galaxy center would change the H i gas fractions by a few percent at most, and generally much less). Drawn from the Sloan Digital Sky Survey spectroscopically confirmed targets with Arecibo Legacy Fast ALFA 21 cm observations, this sample consists of isolated galaxies with stellar masses 10⁷ ≲ M_⋆/M_⊙ ≲ 10¹¹. We compare this isolated H i sample to the FIRE-2 galaxies in Figure 9; the observed galaxies are shown in gray scatter, while the simulated galaxies are shown by the crosses (colored by morphology, as indicated). Notably, while the disky dwarfs lie on the main relation found by Bradford et al. (2015), the nondisky dwarfs are systematically more gas-poor at fixed stellar mass than the observed galaxies.

Zoom In Zoom Out Reset image size

Next, we contrast the SFH of the nondisky and disky dwarfs. In Figure 10, we show at left the fraction of stellar mass built in starbursts over the last 2 Gyr versus the galaxy's stellar mass. We define a starburst as the snapshots for which

$$\begin{eqnarray}&&\displaystyle \frac{\langle \mathrm{SFR}{\rangle }_{10\,\mathrm{Myr}}}{\langle \mathrm{SFR}{\rangle }_{1\,\mathrm{Gyr}}}\gt 1.5.\end{eqnarray} \tag{ 4 }$$

We find that the FIRE-2 dwarfs that succeed in forming disks are those that either are particularly gas rich (as seen in Figure 10; for easy comparison, we outline the galaxies where the H i gas fraction exceed 0.7 with a dashed circle) or have particularly nonbursty recent SFH for their stellar mass. To visualize the difference in SFH, we show the SFH for a nondisky dwarf at the top (green), and a disky dwarf with a low burst fraction at the bottom (purple). Though both are highly bursty at high redshift, the disky dwarf becomes significantly less bursty at low redshift.

Zoom In Zoom Out Reset image size

As the FIRE-2 simulation suite is a set of cosmological zooms, we will first contextualize our results via comparison to similar studies done on cosmological and noncosmological simulations.

Because of the high resolution necessary to study the stellar outskirts of low-mass halos, there are relatively few cosmological simulations that make similar structural measurements on dwarf galaxies to which we can compare. Pillepich et al. (2019) use TNG50, the highest-resolution volume of the IllustrisTNG simulations (Pillepich et al. 2018), to present the mass-weighted intrinsic shape measurements made in an ellipsoidal aperture at twice the half-mass–radius, 2R_1/2,mass, with a thickness of 0.4R_1/2,mass. TNG50 reaches a baryonic mass resolution of 8.5 × 10⁴ M_⊙ and a DM mass resolution of 4.5 × 10⁵ M_⊙ (about an order of magnitude lower resolution in both baryons and dark matter than the galaxies studied in this work) and a minimum gravitational force softening of 74 pc for the gas (compared to ∼2 pc for the galaxies in this work). TNG50 contains approximately ∼5500 dwarf galaxies with stellar masses in the range 10⁸ < M_⋆/M_⊙ < 10¹⁰. The intrinsic shape of the dwarf galaxies in TNG50 shows a strong dependence on stellar mass, with the 10⁹ < M_⋆/M_⊙ < 10¹⁰ galaxies characterized by thick disks and the 10⁸ < M_⋆/M_⊙ < 10⁹ dwarfs significantly more spheroidal. We see no such trend with mass in the FIRE-2 simulation suite. It is of interest that the high-mass TNG50 dwarfs form thick disks, even in the old stellar population. This is in contrast to the FIRE-2 dwarfs that succeed in forming thick disks; the FIRE-2 stellar disks are dominated by young stars, while the old stellar populations that dominate the stellar mass budget maintain largely spheroidal shapes at all radii (for cumulative measurements). Observations indicate that high-mass dwarfs should generically host stellar disks, as they do in TNG50 (Kado-Fong et al. 2020), but also find morphological differences as a function of age in resolved star studies, as is seen in FIRE-2 (Aparicio & Tikhonov 2000; Aparicio et al. 2000; Zaritsky et al. 2000; Hidalgo et al. 2003; Demers et al. 2006; Bernard et al. 2007; Stinson et al. 2009; Strader et al. 2012; Nidever et al. 2019; Pucha et al. 2019).

It is also informative to compare to the higher-resolution, noncosmological simulations of dwarf galaxies. In particular, Smith et al. (2021) examine the effect of photoionization and photoelectric heating from young stars on the efficacy of SN feedback and the resultant impact on the SFH and disk structure. They find that the inclusion of photoionization significantly dampens the burstiness of the simulated dwarf SFH, and that the runs that include only SNe feedback result in superbubbles that significantly disrupt the gaseous disk. It is important to note that the galaxy studied in Smith et al. (2021) is initialized with both a stellar and gas disk, and is significantly lower in stellar mass (M_⋆ ∼ 10⁷ M_⊙) than the galaxies studied in this work; it is observationally uncertain whether galaxies in this mass range are generally disky (Kado-Fong et al. 2020; Carlsten et al. 2021). However, their finding that overvigorous star formation feedback can disrupt the disk is in line with our finding that it is the bursty and gas-poor FIRE-2 dwarfs that do not form young stellar disks.

Finally, it has been noted that the dark matter particles may induce numerical heating in the stellar component of simulated galaxies. We do not expect this numerical heating to significantly affect our results—the resolution of the FIRE-2 dwarfs surpasses the dark matter particle mass limit estimated by Ludlow et al. (2021) at which numerical heating accounts for no more than 10% that of the virial velocity (m_DM ≲ 4.5 × 10⁴ M_⊙; the maximum particle mass of the FIRE-2 dwarfs is m_DM = 3.8 ×10⁴ M_⊙). Indeed, for the massive FIRE-2 galaxy m12i, Ludlow et al. (2021) estimate that the change in scale height attributable to numerical heating is only 70 parsec, which is much smaller than the change in scale height observed in the FIRE dwarfs (see Figure 7). Additionally, a study of the somewhat smaller (M_halo ∼ 10¹⁰ M_⊙) FIRE-2 galaxies run at 30M_⊙ resolution showed that the structural measurements in dwarfs at fiducial FIRE-2 resolution are robust to resolution effects (Wheeler et al. 2019).

Due to the steepening stellar-to-halo mass relation at low masses, it has long been thought unlikely for dwarf galaxies to assemble a significant stellar halo population via the accretion of satellite galaxies. Thus, the observed presence of extended old (Lin & Faber 1983; Minniti & Zijlstra 1996; Minniti et al. 1999; Aparicio & Tikhonov 2000; Aparicio et al. 2000; Zaritsky et al. 2000; Hidalgo et al. 2003; Demers et al. 2006; Bernard et al. 2007; Stinson et al. 2009; Strader et al. 2012; Nidever et al. 2019; Pucha et al. 2019) and round (Kado-Fong et al. 2020) stellar outskirts around dwarf galaxies points to a different mode of stellar halo formation at low masses.

In this work, we have shown that the FIRE-2 dwarf galaxies are capable of producing an extended old and round stellar population in conjunction with a young stellar disk, and we find that the stellar outskirts of the FIRE-2 dwarfs reproduce several key qualities of observed high-mass dwarf stellar halos. However, FIRE-2 also produces a large contingent of nondisky dwarfs whose structural and gas properties are at odds with the observed population. Nevertheless, it is informative to consider the mechanism through which the high-resolution and cosmological dwarfs simulated in FIRE-2 form disk-halo systems, and to consider the plausibility of such a mechanism in the context of our knowledge of the observed dwarf population.

First, in a marked contrast to the ex situ stellar halos of massive galaxies, these dwarf stellar halo populations are dominated by the in situ stars that occupy an increasingly spheroidal component as a function of age. This behavior explains the apparent ubiquity of round stellar outskirts around high-mass dwarfs (M_⋆ > 10^8.5 M_⊙), as feedback should be a generic process in all galaxies that proceeds regardless of accretion history or environment. The existence of an old, in situ population in the outskirts of the FIRE-2 dwarfs is not surprising, as it has been shown that the radial migration of stars due to baryon-driven potential fluctuations operates most strongly in FIRE-2 dwarfs in this mass range (El-Badry et al. 2016; Graus et al. 2019). This feedback-driven migration has previously been suggested as a formation mechanism for in situ stellar halos (Stinson et al. 2009; Maxwell et al. 2012). Star formation feedback-driven size fluctuations have also recently been found via correlations between size and star formation rate in both the FIRE-2 dwarfs and a sample of nearby dwarfs from the Local Volume Legacy Survey (Emami et al. 2021).

Furthermore, the old stellar population of the FIRE-2 dwarfs, while largely spheroidal compared to the young stellar disk populations, is significantly more flattened (a lower minor-to-major intrinsic axis ratio C/A) than are the halos of the higher-mass FIRE-2 galaxies. This is again due to their in situ formation mechanism, wherein the disk gradually thickens into the halo as the stars age. These preferentially flatter dwarf stellar halos are in good agreement with observations of dwarf stellar outskirts, which find that dwarf stellar outskirts at 4 effective radii are preferentially flatter than the MW at the same distance. Indeed, the mean C/A of the young stellar population of disky dwarfs in FIRE (〈C/A〉 = 0.21 ± 0.03, where 〈〉 denotes the sample mean (see Table 2), is similar to the maximum a posteriori estimate of the mean C/A of the massive dwarf population inferred by Kado-Fong et al. (2020) (〈C/A〉 = 0.30 ± 0.01) for dwarfs of stellar mass 10⁹ < M_⋆/M_⊙ < 10^9.6 and 〈C/A〉 = 0.32 ± 0.02 for dwarfs of stellar mass 10^8.5 < M_⋆/M_⊙ < 10⁹ (see Table 1 of Kado-Fong et al. 2020). However, we stress that, because of the differences in the intrinsic shape measurement method (differential versus cumulative) and the exact stellar populations being measured (i-band light-weighted versus age-separated), we are not emphasizing here the apparent agreement in the numerical value of 〈C/A〉 between the observed and simulated dwarfs. Rather, we point to the difference in 〈C/A〉 of the simulated (observed) dwarfs and simulated (observed) high-mass galaxies as the key observational constraint on the difference in the assembly mode of their stellar halos.

Table 2. Mean Intrinsic Shape Ratios of the FIRE-2 Dwarf Sample

	$\langle \tfrac{B}{A}{\rangle }_{\mathrm{young}}$	$\langle \tfrac{C}{A}{\rangle }_{\mathrm{young}}$	〈N〉young	$\langle \tfrac{B}{A}{\rangle }_{\mathrm{old}}$	$\langle \tfrac{C}{A}{\rangle }_{\mathrm{old}}$	〈N〉old
Nondwarfs	0.95 ± 0.02	0.13 ± 0.01	9.6 × 105	0.93 ± 0.02	0.69 ± 0.04	9.2 × 105
Disky Dwarfs	0.82 ± 0.03	0.21 ± 0.03	6.3 × 104	0.89 ± 0.04	0.50 ± 0.04	8.6 × 104
Nondisky Dwarfs	0.55 ± 0.05	0.39 ± 0.04	9.9 × 104	0.56 ± 0.02	0.42 ± 0.02	9.6 × 104

Note. Mean intrinsic shape measurements (denoted by 〈〉) and the associated error on the mean for the sample considered in this work, divided by morphology. 〈N〉 refers to the number of star particles used in the shape computation for the old (ages greater than 12 Gyr) and young (ages less than 2 Gyr) stellar populations.

Download table as:  ASCIITypeset image

Before continuing on from this point, we stress again that the FIRE-2 simulation suite consists of a relatively small number of initial conditions, and that the isolated dwarfs in the simulation suite are biased toward the conditions that are more isolated than the typical field dwarf. Due to the small number of initial conditions and the selection of these initial conditions, it is key to remember that the FIRE-2 dwarfs should not be treated as a dwarf population study. Furthermore, of these relatively small number of initial conditions, the FIRE-2 suite produces a large number of nondisky dwarfs that are inconsistent with average observed dwarf galaxies in both their stellar structure and gas content. We thus emphasize that our goal is not to validate the ability of the FIRE-2 simulations to reproduce the overall dwarf population, but to demonstrate that the in situ mode of dwarf stellar halo creation that has been suggested from observational results can indeed be naturally produced by contemporary cosmological hydrodynamical simulations.

Although the FIRE-2 simulations succeed in producing dwarfs with a disk-halo structure that well reflects the observations, it must be noted that these disky dwarfs are not the dominant mode of dwarf stellar structure in the FIRE-2 simulations. Indeed, most of the dwarfs in the sample have significantly lower values of B/A than the observed dwarfs, and do not host a clear young stellar disk. To better understand this discrepancy, we suggest three possibilities below.

First, the production of the nondisky dwarfs at z = 0 could stem from the implemented physics of the FIRE-2 simulations. That the nondisky dwarfs tend to be gas-poor with a bursty recent SFH (see Figure 10) suggests that the star formation feedback prescription in FIRE-2 may be preventing the formation of both stellar and gaseous disks in these nondisky dwarfs (see also El-Badry et al. 2018b). It is also notable that the observed high-redshift dwarfs lack a well-defined disk; Zhang et al. (2019) found that high-z, high-mass dwarfs are characterized by triaxial prolate ellipsoids at high redshift (z ≳ 1) and transition to thick-disk structures at low-z (see also Ceverino et al. 2015, for a similar transition in simulated dwarfs). It is thus possible that the FIRE-2 dwarfs are failing to undergo this transition at the appropriate cosmic time. At the same time, the existence of the disky dwarfs, which also host in situ stellar halos, indicates that it is possible to form disk-halo structures with this feedback prescription. Interestingly, El-Badry et al. (2018a) also showed that the low-mass FIRE-2 galaxies tend to be overly dispersion-supported due to inefficient accretion of high angular momentum gas, with the same trend that the rotation-supported gas disks occupy the galaxies with less bursty recent SFH. This interplay between the gaseous structure, stellar structure, and recent SFH of the FIRE-2 dwarfs points to a link between the way in which star formation proceeds in these dwarfs and their ability to form a young stellar disk. Indeed, the idea that star formation feedback can significantly influence dwarf structure is well established in regards to both the baryonic morphology (see, e.g., Governato et al. 2007; Smith et al. 2019, 2021) and in the dark matter halo (the core-cusp problem; see, e.g., Peñarrubia et al. 2012; Pontzen & Governato 2012). Indeed, as the most easily observable leg of this set, the stellar structure of galaxies in the dwarf mass regime could prove to be a powerful tool to constrain the feedback prescriptions in conjunction with requirements for cored/cuspy dark matter profiles.

Second, it is possible that the observations are missing a significant population of nondisky dwarfs due primarily to surface brightness sensitivity limits. We find that the nondisky dwarfs are not systematically lower in stellar surface density (within the stellar half-light radius) than the disky dwarfs. However, there is evidence that the FIRE-2 dwarfs with more bursty star formation also have lower surface brightnesses (Chan et al. 2018). It is thus plausible that the observed samples preferentially include disky dwarfs.

Explaining the lack of observational analogs to nondisky dwarfs in FIRE-2 as being due to observational incompleteness would also necessitate the existence of a substantial and likely dominant population of low surface brightness dwarfs at 10⁸ ≲ M_⋆/M_⊙ ≲ 10¹⁰. A forthcoming probe of the distance and stellar mass distribution of a sample of low surface brightness galaxies will help to directly address this question (Greco et al., in preparation). However, the observational studies of the intrinsic shapes of low surface brightness galaxies (Kado-Fong et al. 2021) and ultra diffuse galaxies in clusters (Burkert 2017; Rong et al. 2020) indicate that, although LSBGs are relatively round, they do not show significantly lower values of B/A than spectroscopic dwarf samples, as is the case for the nondisky dwarfs. We thus find it unlikely that the nondisky dwarfs represent a dominant low surface brightness population that is undetected in current-generation surveys. Similarly, because the nondisky galaxies lack a disk at all stellar age bins considered, we find it unlikely that the difference in the observed and simulated intrinsic shape distributions originates from a difference between light-weighted differential and age-binned, mass-weighted cumulative shapes. We note that some of the nondisky dwarfs are consistent with previous inferred shapes for cluster ultra diffuse galaxies (Rong et al. 2020). ¹⁵ It is of interest, but beyond the scope of this work, to ask whether these galaxies differ from the general set of nondisky dwarfs.

A last possibility is that, because the FIRE-2 suite represents a relatively small number of dwarf galaxies, the nondisky dwarfs may be a smaller proportion of the total dwarf population than is reflected in the FIRE-2 suite. Part of this discrepancy could be due to bias in the dwarfs chosen for the suite. Galaxies in FIRE-2 zooms are selected to be isolated, and to have small Lagrangian regions, and so it is possible that there is some systematic bias in the sample that causes them to be more nondisky. However, the observational results do not indicate that the nonisolated dwarfs are more nondisky (Kado-Fong et al. 2020); further investigation into the behavior of the FIRE-2 zoom suite as a function of environment would be enlightening, but is outside the scope of this work. To place a limit on the sensitivity of the observational methods of Kado-Fong et al. (2020) to a morphologically distinct nondisky population, we construct a toy model wherein we simulate a mock galaxy population composed of a disky and nondisky population and rerun our observational analysis, as detailed in Appendix A. We find that the hypothetical presence of a small minority (≲10%) of galaxies with nondisky shapes would not affect the observationally inferred intrinsic shape parameters. However, the fraction of nondisky galaxies in FIRE-2 greatly exceeds this 10% threshold. Thus, we find that the overrepresentation of nondisky dwarfs in the simulated sample is unlikely to be explained by observational uncertainties or by Poisson fluctuations (N_non−disky = 10 ± 3).