ULTRAFAINT DWARF GALAXIES—THE LOWEST-MASS RELICS FROM BEFORE REIONIZATION

Our interest in galaxy systems with the lowest possible mass arises from a desire to understand the chemistry of the first stars (Frebel & Bromm 2012; Frebel et al. 2014, K13). Our assumption is that the surviving remnants of the first galaxies will have small dark halo masses and possibly low baryon mass fractions for the reasons outlined in Section 1. Many of these systems will have merged with larger galaxies long ago, but some are likely to be falling into larger systems today for the first time.

So what constitutes a first galaxy? This must retain a long-lived stellar system and importantly a confining dark matter halo.⁵ Whether low-mass stars formed alongside or immediately after the first stars is unknown. It is likely that the process was patchy, with the very first stars forming in rare density peaks, but with an increasing fraction of the universal volume becoming ionized as the first stellar generation took hold.

Here we make the simple assumption that the universe was enriched to a threshold metallicity of order [Fe/H] = −4 before reionization was complete (e.g., Bromm & Yoshida 2011). Some of this gas was accreted onto low-mass dark matter halos, and star formation began. Our threshold metallicity allows us to consider a "normal" IMF in order to describe how star formation proceeded at these early redshifts (Bromm et al. 2001; Schneider et al. 2002). Given the stochasticity of star formation at our low masses, there may be a significant amount of low-mass star formation before the first SN event. Some remnants of first galaxies may therefore contain stars formed from gas enriched only by Population III stars, thereby providing direct information on the first stellar yields.

However, in the systems relevant to our work (UFDs, very metal-poor DLAs), some of the detected elements (e.g., Sr) can only have arisen from SN ejecta. Thus, our minimum-mass models have the more stringent requirement that they can survive one or more SN events during their lifetime. In order to experience more than a single star formation episode, a very small galaxy must retain enough gas in the face of the energy outflow from a single star-forming episode for a second one to occur, perhaps at a much later time. A large burst of energy such as an SN could either heat the galaxy ISM and halt star formation or completely evacuate the star formation region of gas. This will cut off star formation until such time as the halo re-accretes the blown-out gas, assuming that this is allowed by the environment. Thus, a "first galaxy" is one that must retain sufficient gas to form stars in the wake of at least a single SN event.

This is not a contrived circumstance, and we show below (Section 6) that there may already be observations to support such a scenario. It is likely that there are situations in which only one or two SNe occur within a single star-forming "burst" (Karlsson et al. 2013; Ritter et al. 2014). Hurley-Keller et al. (1998) identify several starburst epochs in the Carina dwarf, with the first of these taking place of order ∼10 Gyr ago. At that time, the surface density of the star formation rate (${s}_{\star }$) within the core was roughly 10⁻³ ${M}_{\odot }$ yr⁻¹ kpc⁻² over the period of the burst, T₀. This timescale is now well established in starburst and dwarf galaxies such that ${T}_{o}\lesssim 10$ Myr (Tolstoy et al. 2009; Sharp & Bland-Hawthorn 2010).

To test the plausibility of a single SN event occurring in our models, we consider a simple model of star formation. We consider dark matter potentials with an Einasto profile over the virial mass range $\mathrm{log}[{M}_{\mathrm{vir}}({M}_{\odot })]$ = 5.5, 6.0, 6.5, and 7.0; we refer to these models as M55, M60, M65, and M70, respectively (see Figure 1). If we integrate over the Einasto profile, the total masses for these models in dark matter are $\mathrm{log}[{M}_{\mathrm{tot}}({M}_{\odot })]$ = 5.9, 6.4, 7.0, and 7.5. (These total masses do not include the small amount of baryons.) In Figure 2, we show how the low-mass halos are expected to grow in mass and radius with cosmic time. Halo growth is not treated in our simulations because the time frame of our analysis is short. Our results at a given halo mass at z = 10 correspond to more massive systems at later times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we show complete results for M65 and M70 and partial results for M55 and M60 because these models fail to retain sufficient baryons. The initial baryon fraction is assumed to be roughly 10% within the virial radius (e.g., Krumholz & Dekel 2012). We assume that star formation occurs where the gas is densest, i.e., within a radius enclosing 50% of the baryons (${r}_{0.5g}$) inside one scale radius r_s, the scale radius of the Einasto potential (see Section 3.1). Over this volume, we find that the gas density is within a factor of 6 of the maximum central value.

For the models (M55, M60, M65, M70), the enclosed gas masses within the half gas-mass radius (≲${r}_{{\rm{s}}}$/2) are ($2.8\times {10}^{3}$, 9.5 × 10³, 2.5 × 10⁴, 7.5 × 10⁴) ${M}_{\odot }$, and thus the average gas surface densities are (2.1, 2.6, 3.3, 4.2) ${M}_{\odot }$ pc⁻². In the clumpy gas models, there are order-of-magnitude variations in local gas density. From direct observation, we know that star formation takes place at these moderately low column densities over a range of metallicities down to [Fe/H] ≈ −1 (e.g., Bigiel et al. 2008). As these authors show, the star formation efficiency is observed to vary wildly below 9 ${M}_{\odot }$ pc⁻², with a median value of only a few percent. Little is known about how the star formation rate depends on metallicity, but it seems likely that it declines with [Fe/H] owing to the difficulty of forming H₂. We stress, however, that the inferred value of ${s}_{\star }$ in Carina is comparable to what is observed in the outer metal-poor regions of spiral galaxies (e.g., Bigiel et al. 2008), and this dwarf galaxy somehow managed to form long-lived stars at much lower metallicities.

For a star formation efficiency of 2% (Kennicutt 1989), the expected total stellar mass over the lifetime of the system is (56, 190, 670, 2400) ${M}_{\odot }$. If we assume a Kroupa mass function, roughly one-fifth of the mass in stars will end their lives as SNe. These estimates are susceptible to a lot of scatter in the low-mass/surface density limit (Bigiel et al. 2008). If we scale the projected star formation rate per unit area to the largest burst in Carina (∼10⁻³ ${M}_{\odot }$ yr⁻¹ kpc⁻²), the timescales required to form these stellar masses are (44, 54, 71, 86) Myr, for a half gas-mass radius of (20, 34, 55, 94) pc (see Section 3.1 for further details on the halo properties). If a typical burst lasts 10 Myr, as in Carina (Tolstoy et al. 2009), we would expect multiple bursts over the star-forming lifetime of the galaxy.

So can we consider an isolated SN event to be applicable to all of our models? It is conceivable that more than one event can occur within the timescale of a single burst T_o. The results given in Table 1 show that multiple SN events are rare in the M55 and M60 models and are relatively unusual in the M65 model, although they occur half the time in the M70 model. Multiple SNe in more gas-rich M70 models happen half the time, but in practice we can consider most of these eventualities to be "single SN events" because one of these is highly likely to be a low-mass event and physically separated from the other event. These events are likely to enrich their own immediate environs, although substantial mixing can occur between the events at late times (Ritter et al. 2014). While different choices of the size of the star-forming region or baryon fraction may alter these results, we would not expect this to alter the overall conclusion of a single SN event. We note that simulations of more massive halos will need to consider multiple events within the same burst.

Table 2.  Model Parameters

Cosmology
	z	h	${{\rm{\Omega }}}_{{\rm{\Lambda }}}$	${{\rm{\Omega }}}_{M}$
	10.0	0.7	0.7	0.3
α	${{\rm{\Delta }}}_{c}$	${\rho }_{c}$(g cm−3)	${{\rm{\Delta }}}_{u}$	$\langle {\rho }_{u}\rangle$(g cm−3)
0.18	177.5	3.68E-027	591.7	1.11E-027
Dark Matter
	rs	${\rho }_{{\rm{s}}}$	${x}_{\mathrm{vir}}$	${r}_{\mathrm{vir}}$
Model	(pc)	(g cm−3)		(pc)
M55	33.2	1.02E-023	5.98	198.5
M60	55.0	7.73E-024	5.29	291.3
M65	91.3	5.87E-024	4.69	427.6
M70	151.3	4.48E-024	4.15	627.7
	Ms	${M}_{\mathrm{vir}}$	M300	${M}_{\mathrm{tot}}$
Model	$({M}_{\odot })$	$({M}_{\odot })$	$({M}_{\odot })$	$({M}_{\odot })$
M55	5.64E+04	3.16E+05	4.10E+05	7.90E+05
M60	1.94E+05	1.00E+06	1.02E+06	2.73E+06
M65	6.73E+05	3.16E+06	2.36E+06	9.44E+06
M70	2.34E+06	1.00E+07	5.08E+06	3.29E+07
Baryons
	cs	${c}_{{\rm{s}}},\mathrm{iso}$	${T}_{\mathrm{gas}}$	${M}_{\mathrm{gas},\mathrm{vir}}$
Model	(km s−1)	(km s−1)	(K)	$({M}_{\odot })$
M55	2.70	1.91	5.33E+02	3.87E+04
M60	3.90	2.76	1.11E+03	1.19E+05
M65	5.63	3.98	2.32E+03	3.66E+05
M70	8.16	5.77	4.86E+03	1.15E+06
	${r}_{0.5\;{\rm{g}}}$	${\rho }_{\mathrm{gas},0}$	${n}_{{\rm{H}},0}$	${M}_{\mathrm{gas},{\rm{s}}}$
Model	$(\mathrm{pc})$	(g cm−3)	(cm−3)	$({M}_{\odot })$
M55	20.8	4.09E-23	17.2	5.63E+03
M60	34.5	3.11E-23	13.1	1.94E+04
M65	57.1	2.36E-23	9.91	6.73E+04
M70	94.8	1.79E-23	7.55	2.34E+05

Note. The cosmology is fixed by the halo formation redshift, z, the normalized hubble constant h = H_o/100 km s⁻¹ Mpc⁻¹, the dark energy density ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$, and the dark matter density ${{\rm{\Omega }}}_{M}$.

The Einasto halo properties (see text) are fixed by the concentration parameter, α, and the halo dark matter density contrast with the ambient medium, Δ, defined at redshift z. This density contrast can be expressed as ${{\rm{\Delta }}}_{c}$ averaged over the virial radius and normalized to the mean critical density ${\rho }_{c}$, or as ${{\rm{\Delta }}}_{u}$ normalized to the universal mass density $\langle {\rho }_{u}\rangle$. The halo potential is defined by its scale radius, r_s, and the mean density and mass within this radius, ${\rho }_{{\rm{s}}}$ and M_s. The halo virial mass is ${M}_{\mathrm{vir}}$ within the virial radius ${r}_{\mathrm{vir}}={x}_{\mathrm{vir}}{r}_{{\rm{s}}}$; M₃₀₀ is the mass within 300 pc, and M${}_{\mathrm{tot}}$ is the mass integrated to infinity. For the Baryons, c_s and ${c}_{{\rm{s}},\mathrm{iso}}$ are the adiabatic (physical) and isothermal sound speeds in the gas with mean temperature ${T}_{\mathrm{gas}}$. This gas has a central density (number density) of ${\rho }_{\mathrm{gas},0}$ (${n}_{H,0}$) with a virial (total) mass of ${M}_{\mathrm{gas},\mathrm{vir}}$ (${M}_{\mathrm{gas},\mathrm{tot}}$). Note that the initial baryon fraction ${f}_{{\rm{b}}}\approx {M}_{\mathrm{gas},{\text{}}\ {\rm{s}}}/{M}_{{\rm{s}}}\approx {M}_{\mathrm{gas},\mathrm{vir}}/{M}_{\mathrm{vir}}$ ≈ 10% and the radius containing half the initial gas within r_s is ${r}_{0.5\;g}$.

Download table as:  ASCIITypeset image

Massive stars evolve rapidly in the lead-up to the explosion, during which time the stellar wind and radiative luminosity can be substantial. This pre-SN phase can be an important factor in removing gas, particularly for halos with low baryon fractions. In the worst-case scenario, a star at the halo center can ionize an amount of gas, ${M}_{\mathrm{ion}}\sim {S}_{\star }{m}_{{\rm{p}}}/\alpha {n}_{{\rm{H}}}$ $\;{M}_{\odot }$, where α is the total recombination coefficient, ${m}_{{\rm{p}}}$ is the mass of the proton, and ${n}_{{\rm{H}}}$ is the hydrogen density. Using the Schaerer (2003) models at a metallicity of [Fe/H] = −4.0, the mean ionizing rate averaged over the Kroupa IMF is ${S}_{\star }\sim 4\times {10}^{48}$ photons s⁻¹ for the M55 model, $1\times {10}^{49}$ photons s⁻¹ for M60, and $3\times {10}^{49}$ photons s⁻¹ for M65. Taking the worst-case scenario where the most massive star falls near the center, the total ionized masses are $1.7\times {10}^{3}$ $\;{M}_{\odot }$ for M55, $6.1\times {10}^{3}$ $\;{M}_{\odot }$ for M60, $2.2\times {10}^{4}$ $\;{M}_{\odot }$ for M65, and $5.5\times {10}^{4}$ $\;{M}_{\odot }$ for M70. In each case, a large fraction of the neutral gas can be ionized. This is likely to have a significant impact in reducing subsequent star formation. But the more likely situation is not expected to be so severe. Because the total amount of star formation is small in our models, stochastic effects can dominate. The median radiation levels expected are more biased to B stars rather than rare O stars, and the ionized masses can be much less than 10%. We include the effects of preionization in our models below.

In our experiments, we look at a series of spherically symmetric gas distributions in hydrostatic equilibrium, detonate a single SN, and track the ability of the potential to retain significant amounts of gas over 25 Myr or more. The parameters of our models are shown in Table 2. For all of our models, we consider the period prior to the onset of cosmological reionization when stars were forming for the first time. At such early times, substantial amounts of gas have already begun to settle into dark matter halos up to and exceeding ${M}_{\mathrm{crit}}$ (Bullock et al. 2001; Ricotti & Gnedin 2005). This is consistent with the latest results from deep photometry for UFDs, where star formation appears to have been quenched by reionization in almost all cases (Brown et al. 2012, 2014; Webster et al. 2015a).

We consider both clumpy and smooth gas distributions within the low-mass halos. For the clumpy models, we make the important assumption that the confined gas has a fractal distribution before the SN explodes. There are two reasons for this assumption. First, the infalling gas becomes shocked and highly inhomogeneous close to the dark halo core (Abel et al. 2000). Second, a high-mass star that is destined to become an SN injects mechanical and ionizing energy into the surrounding gas prior to the explosion. This energy input drives turbulence in the surrounding gas (Elmegreen & Scalo 2004).

In looking for possible key factors in the retention, we consider the following scenarios.

• 1.  
  Abundances and cooling: We use chemical abundances that approximate primordial conditions starting at [Fe/H] = −4.0. This starting metallicity is conservative: Wise & Abel (2008) show that the initial metallicities of the first galaxies could have been as elevated as [Fe/H] ∼ −3.0 or higher. We compute models with cooling, and models $\left(\gamma =\frac{5}{3}\right)$ where we have artificially turned off microphysics affecting energy loss (e.g., cooling through metal and recombination lines). The cooling function is presented in Figure 3; the abundances are given in Appendix A.
• 2.  
  Clumpiness: We consider models with quasi-fractal cool star-forming gas, and comparison models where the gas is homogeneous and smooth. The key point here is that the degree of coupling between the SN energy and the surrounding medium will depend on whether the medium is smooth or structured.
• 3.  
  Central and off-center explosions: We consider a single 10⁵¹ erg SN explosion, located either at the geometrical center of the halo potential or off-center at a radius that encloses half the mass contained inside one scale radius of the potential. SNe that occur away from the center have an asymmetric environment.

This gives, in combination, eight scenarios for each of the halo masses considered. In each case, the mass and energy contained in the inner scale radius are integrated at 25 kyr intervals and plotted over the course of the simulation to obtain estimates for the masses retained at late times.

Zoom In Zoom Out Reset image size

In the Fyris Alpha models presented here, we use a fixed dark matter potential that dominates everywhere. Future models will include self-gravity in the gas. Therefore, we limit the initial baryonic mass inside a scale radius (r_s) to be 10% of the dark matter mass inside that radius. This results in an overall baryonic fraction out to the virial radius (${r}_{\mathrm{vir}}$) of ${f}_{{\rm{b}}}$ ≈ 12%, less than the canonical ${f}_{\mathrm{crit}}$ $\approx \;17$% global baryonic mass fraction (Planck Collaboration et al. 2014). Our choice of ${f}_{{\rm{b}}}$ is a conservative assumption because a lower gas budget generally assists the SN in accelerating and heating the gas. If the models retain gas with these lower mass fractions, they can be reasonably expected to do so with more gas in place.

In recent years, high-resolution cold dark matter (CDM) models have identified the need for a halo density profile that has an additional parameter beyond those used in the Navarro–Frenk–White (NFW) profile (Hayashi et al. 2004; Springel et al. 2008). The profiles seen in the latest simulations display local density derivatives that are power laws of radius, unlike the NFW profile that has asymptotic behavior at small radius. They invoke a generalized exponential function that has this property and note that it has been used in the past by Einasto (1965) in the context of galactic halo spatial denities.⁶

We use the Einasto function, a generalized exponential, to form a dimensionless potential defined by the parameters α, a halo concentration factor x_vir (often referred to as c), and the density contrast of the halo Δ (≈178) with respect to the universal background (Nichols & Bland-Hawthorn 2009). We stress that these parameters must be specified at a fixed redshift z because of their strong cosmic evolution (Power et al. 2003; Duffy et al. 2008).

We consider four halo mass models (M55, M60, M65, M70) defined in Section 2.2. The total halo masses are (3.0, 2.7, 2.5, 2.3) times the virial mass, respectively (e.g., Figure 4). The halo and gas properties are summarized in Table 3. Our normalization uses a density function value of 1/e at the scale radius, r_s. This gives a slightly different set of scaling constants, but recovers exactly the same physical quantities as the 1/4 scaling used in Nichols & Bland-Hawthorn (2009). This change was motivated by wanting to compare with dimensionless isothermal potentials used by Sutherland & Bicknell (2007), which also obtained a scaled density value of 1/e at the scale radius. (We refrain from using the popular "core" radius because, unlike the isothermal potentials, the generalized Einasto functions do not have a distinct core region; here we use r_s throughout.)

Zoom In Zoom Out Reset image size

Table 3.  Naming Convention for All Numerical Simulations

Models without Preionization
	Central Explosions				Off-center Explosions
	Smooth		Clumpy		Smooth		Clumpy
Halo Mass	Adiabatic	Cooling	Adiabatic	Cooling	Adiabatic	Cooling	Adiabatic	Cooling
${10}^{5.5}$ $\;{M}_{\odot }$	M55CSA	M55CSC	M55CCA	M55CCC	M55OSA	M55OSC	M55OCA	M55OCC
106 ${M}_{\odot }$	M60CSA	M60CSC	M60CCA	M60CCC	M60OSA	M60OSC	M60OCA	M60OCC
${10}^{6.5}$ $\;{M}_{\odot }$	M65CSA	M65CSC	M65CCA	M65CCC	M65OSA	M65OSC	M65OCA	M65OCC
107 ${M}_{\odot }$	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
Models with Preionization
	Central Explosions				Off-center Explosions
	Smooth		Clumpy		Smooth		Clumpy
Halo Mass	Adiabatic	Cooling	Adiabatic	Cooling	Adiabatic	Cooling	Adiabatic	Cooling
${10}^{5.5}$ $\;{M}_{\odot }$	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
106 ${M}_{\odot }$	⋯	⋯	⋯	M60CCH	⋯	⋯	⋯	M60OCH
${10}^{6.5}$ $\;{M}_{\odot }$	⋯	⋯	⋯	M65CCH	⋯	⋯	⋯	M65OCH
107 ${M}_{\odot }$	⋯	⋯	⋯	M70CCH	⋯	⋯	⋯	M70OCH

Download table as:  ASCIITypeset image

Our adopted halo parameters are based on those of Duffy et al. (2008) and Springel et al. (2008). As their fits to the simulated halos are measured over a relatively low redshift range, we extrapolate these properties back to a redshift of z = 10. At the time of writing, we became aware of new work (Correa et al. 2015) where the parameters are extended back to z = 20, but these appear to be consistent with our estimated values, particularly for the important M65 and M70 models.

In our simulations, we introduce a non-uniform medium to study the influence of inhomogeneity on the dynamical interaction between the SN ejecta and the surrounding gas. To assess the importance of inhomogeneity, we run both homogeneous (smooth) and inhomogeneous (clumpy) models and compare the outcomes. Following Sutherland & Bicknell (2007), to establish a non-uniform medium, we make use of an analogy with a turbulent medium. Rather than attempt to review this huge topic, we refer the reader to astrophysically oriented reviews as a starting point for background material (Elmegreen & Scalo 2004; Scalo & Elmegreen 2004).

We represent the non-uniform properties of the turbulent medium using three statistical characteristics: the variance, ${\sigma }_{F}^{2}$, the intermittency (by using lognormal distributions), and a self-similar power-law structure. We achieve both the property of intermittency and power-law structure simultaneously with an iterative scheme first developed by Lewis & Austin (2002); see Appendix B for more details. Hence, the initial distribution of the ISM that we employ should be regarded as a physically motivated generalization of a homogeneous model.

We further simplify our turbulent medium by neglecting the velocity structure, and we focus on choosing statistical parameters to describe the density alone. This is justified numerically by the relatively small turbulent velocities expected when compared to the very large velocities found in our global SN–ISM interaction. The velocities observed in typical warm ISM conditions fall in a range of transonic to mildly supersonic values, Mach 1–5 (e.g., Heiles et al. 2004). At temperatures at or below 10³ K, this corresponds to velocities $\ll$4 km s⁻¹. In our 0.1–1 kpc simulations over 10⁷ yr timescales, the resulting displacements amount to only a few cells at most and are insignificant compared to those of the energy bubble generated by the SN outburst, where velocities of many hundreds of kilometers per second occur. Consequently, we do not impose a turbulent velocity field on the cold ISM medium.

The modeling was carried out using the Fyris Alpha code (Bland-Hawthorn et al. 2007; Sutherland 2010). It solves ideal Euler hydrodynamic flows of gas with cooling and an equation of state appropriate to astrophysical gases, in the presence of a fixed gravitational field. The code uses a third-order semi-Lagrangian shock-capturing method that is robust and well suited to the extreme temperatures and Mach numbers encountered here.

The computational domain consists of a three-dimensional nested series of cubic meshes, as illustrated by the schematic in Figure 5. Each level is a factor of 3 smaller in physical domain and higher resolution than its outer containing level, covering the central third of the enclosing level. Each level is 216³ cells with an overall effective resolution of 1944³ cells.

Zoom In Zoom Out Reset image size

A relatively low resolution outer level, L0, covers the halo out to ±6 scale radii and provides an outer reservoir and boundary condition for the main level of interest, L1. The outer reservoir reduces or avoids completely any boundary artifacts on L1 due to gravity bringing new gas onto the main level when infall conditions apply, keeping level L1 as physically consistent as possible at all times. L0 contains the virial radius. Velocities at the outer boundary are extremely small at all times considered here, and the outer L0 boundary is essentially fixed.

The main level, L1, covers a region of just over ±2 scale radii and is the primary region where the movement and evolution of mass and energy are measured. Early epochs require higher resolution to resolve thin spherical shell radiative blast waves, so a third, inner injection level, L2, is added to resolve the early phases of the SNe. After 25 kyr, L2 is disabled and the remainder of the simulations are carried out on L0 and L1.

For each halo mass, the scale radius varies, and so the resolution is variable. The highest resolution, on level L2, is ${\rm{\Delta }}x\sim 0.25$, 0.37, and 0.62 pc per cell for the M55, M60, and M65 models, respectively. Thus, the effective mass resolution is $\sim {10}^{4}$ $\;{M}_{\odot }$/216³, or of order 10⁻³ ${M}_{\odot }$; in practice, cell densities can vary by more than 4 orders of magnitude across the grid.

The remaining key features of the calculations can be summarized as follows.

• 1.  
  A fixed gravity potential, interpolated locally by third-order PPM interpolation to match the interpolation of the hydrodynamical algorithm, taking conservative potential differences to compute mean accelerations.
• 2.  
  SNe simulated by depositing 10⁵¹ erg of energy as internal gas energy in a smoothed spherical region, as small as possible while retaining spherical geometry well enough to prevent gross non-spherical "pixel" errors, typically 6 cell radius on L2, corresponding to 1.5–3.7 pc. This then takes the form of a pulse of low-density, hot overpressure gas, which then expands, converting internal energy to kinetic energy and doing work against the potential. Energy is also lost via cooling.
• 3.  
  To allow for significant changes in the local metallicity in the models as the metal-rich ejecta propagate through the models, a generalized cooling model has been adopted that allows the local composition of each cell to be determined during the simulation. Instead of a single cooling function and uniform molecular weight, the cooling was separated into three components, representing hydrogen, helium, and heavier metals, and ionization was taken into account.

The metal-rich SN ejecta were followed using a scalar tracer variable, advected passively, representing the SN material, and hence the local helium fraction and metal fraction. The local enrichment of the gas was thus computed at each time step. Furthermore, the ionization fraction of the components, precomputed with the MAPPINGS IV code (Sutherland et al. 1993; Allen et al. 2008) as a function of temperature, allowed the equation of state to vary with the mean molecular weight, ensuring more accurate shock temperatures, and allowing for a wider range of cooling timescales caused by strongly varying composition. Figure 3 shows the resulting computed cooling functions for a range of encountered compositions on the grids.

We acknowledge that there are large uncertainties in the yields from essentially zero metallicity stellar SN models. The yield of metals, particularly oxygen, could have an impact on the outcome, lower-mass SNe giving less enrichment and cooling, while higher-mass stars could give more cooling. More detailed future work incorporating a range of SN progenitor masses will investigate the impact of different SN assumptions. By allowing for the ionization and composition to vary from cell to cell, we obtain more accurate species columns, for example, for neutral and ionized hydrogen.

Our focus is on the hotter atomic and ionized phases of high-pressure gas and the cooling therein, as that gas is volume filling and will have the primary effect on the global dynamic evolution if any. Gas in the molecular phase is likely to be confined to the centers of the densest cores, if present at all. As the outflows are all warm (>10⁴ K), they are expected to be atomic or ionized in any case. In the face of uncertainty in the formation of H₂ under primordial conditions and little or no information on CO cooling, we leave molecular cooling out. It would be essential to include it if we were tracking the (much longer) time required for the hot gas to re-cool back to star-forming temperatures $\ll 100$ K, but here we are focused on the initial halo gas retention problem. We defer the star formation problem to more substantial investigations.

That we neglect molecular cooling in these simulations can be considered a conservative approximation, but we expect it to have little effect on the global dynamics of gas retention. More than 99.9% of the internal energy is lost by the time the gas cools back below 10⁴ K from early temperatures of more than 10⁷ K, so only a tiny fraction of the SN energy is available to the molecular cooling in the present simulations.apj515565t3

In Table 4, we calculated the thermalization efficiency for a subset of the models as a check that our simulations correctly handle cooling. The thermalization efficiency was estimated by finding the ratio of the total energy at the time when gas is first leaving the grid to the total initial energy. Our efficiencies should be considered upper bounds, as it is likely that some cooling will occur at later times. From Thornton et al. (1998), we expect the thermalization efficiency to be ∼10%, which is consistent with the results of our models.

The preionization of the gas by the SN precursor plays a key role in the loss or survival of gas after the SN explosion, particularly in low-mass dark matter halos. We consider in detail the impact of an ionized H ii region that forms around the progenitor star prior to the explosion. The key parameters of this model are the stellar properties (T_eff, L) and the associated stellar atmospheres. These give the ionizing photon flux as a function of time up to the SN event. Transforming the grid coordinates to spherical coordinates around the star, radial summations to approximate the optical depth integrals are performed, resulting in thermal and ionization structures calculated from the MAPPINGS IV ionization code. These are then mapped back to Cartesian coordinates and used in the Fyris fluid computations. It was necessary to perform this remapping for every time step in the fluid dynamics simulation in order to track the evolution in opacity. This slowed the code by a factor of 5–10. We considered both clumpy and smooth external media with and without gas cooling.

We model the SN progenitor with a low-metallicity star of 25 ${M}_{\odot }$. Higher-mass stars produce stronger ionizing intensity and therefore assist gas loss by overheating the gas prior to the SN explosion. But since 80% of all SN events in a conventional IMF have lower mass and therefore weaker ionizing radiation fields, we consider our choice of the upper mass limit to be conservative. We use the evolutionary tracks of Meynet & Maeder (2002, hereafter MM02) for a rapidly rotating, low-metallicity massive star. The 25 ${M}_{\odot }$ track was missing from that work, and so we have interpolated the main-sequence phase of this model from the other stellar models. In Figure 6(A), we show the interpolated values of the luminosity and temperature for the 25 ${M}_{\odot }$ model at [Fe/H] ∼ −3.7 (MM02). These are the closest published tracks for our required starting metallicity of [Fe/H] = −4.0. We used the tabulated MS lifetime of 6.1 Myr; the interpolated temperature–luminosity track is shown in Figure 6(B).

Zoom In Zoom Out Reset image size

In Figure 7, we show the impact of the total ionizing budget during the pre-SN phase for our adopted stellar model and a perfect blackbody of the same surface temperature. The panels A, C, and E refer to the models with the progenitor at the center; panels B, D, and F are for the off-centered models. This gives some indication of the sensitivity to the stellar models. If the thin lines exceed the thick lines, then a Strömgren sphere is possible within the radius, and the stellar radiation would be incapable of extending beyond that, which is the initial case for all halos. However, as the gas responds to heating and evolves, the density is smoothed and falls on average owing to outflows. The quadratic dependence of recombinations means that the available recombinations falls rapidly, and once the photon curves exceed the recombination curves, the photoionized zones can grow rapidly and ionize more of the halo.

Zoom In Zoom Out Reset image size

In the upper panels, the M60 halo will be easily fully ionized after the first 2 Myr. The middle panels show that the M65 models are marginal if the photons are uniformly spread throughout the halo, and capable of ionizing the inner core region and possibly beyond before the SN occurs. The actual result will depend on the details of the radiative transfer, whereas this is an indicative integral view of the processes that does not take the geometry of the ionized volume into account. The lower panels for M70 indicate that the stars are insufficient to photoionize a significant volume of the halo.

These indications are broadly supported by the simulations: M60 is rapidly ionized, M65 is able to break out of the core radius before the SNe, and the ionized region is contained in the M70 models.

Stellar winds. To avoid the complexities of rapid post-main-sequence evolution, we assume that the star simply explodes as an SN at the end of the main sequence. Ekström et al. (2006) showed that mass loss during the main sequence in similar stars at very low metallicity may only amount to 1% of the initial mass, and Kudritzki (2002, 2005) showed that the stellar wind luminosities of low-metallicity ([Fe/H] ≤ −4.0) O stars are typically 10³⁴ erg s⁻¹ or less, several magnitudes below an equivalent solar-metallicity star. We combine these values to get a mean wind velocity, for a main-sequence mass loss of 1%:

$$\begin{eqnarray}&&{v}_{w}\sim 880\;{\left[\displaystyle \frac{{L}_{w}}{{10}^{34}\mathrm{erg}\;{{\rm{s}}}^{-1}}\right]}^{1/2}{\left[\displaystyle \frac{M}{25{M}_{\odot }}\right]}^{-1/2}{\left[\displaystyle \frac{{t}_{\mathrm{MS}}}{6.1\mathrm{Myr}}\right]}^{-1/2}\mathrm{km}\;{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 2 }$$

This wind is included in the pre-SN evolution, but the small mass flux and low ram pressure meant that it had negligible effect on the simulations. Once the H ii region pressurized the surrounding region, the stellar wind was unable to blow a wind bubble of any significant size.

The key parameter for these near-zero-metallicity O stars is their unusually high effective temperatures for a given mass. At 25 ${M}_{\odot }$, we estimate the initial temperature to be 48,000 K, compared to a similar-mass star at solar metallicity of only 40,000 K (see, e.g., Meynet & Maeder 2000 and more recently Georgy et al. 2012). This means that the star will be more effective as a source of ionizing photons than the equivalent star in the solar neighborhood. To convert the evolutionary tracks, ${T}_{\mathrm{eff}}(t)$ and L(t), to an ionizing photon–luminosity (see Figure 8), we used the ATLAS9 atmospheric grid (Castelli & Kurucz 2004), which uniquely contains a set of atmospheres for [Fe/H] = −4.0 with alpha-element abundance enhancements, perfectly matching the initial gas composition in our simulations. The interpolated stellar temperature, luminosity, and ionizing fluxes are shown in Figure 8. The derived ionizing spectrum is shown in Figure 9; for comparison, a photon flux assuming a blackbody instead of the stellar atmosphere is included.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our goal is to establish the low-mass limit of a spherical dark matter halo that can retain baryons after a single SN event. But first we consider the impact of the progenitor star on the surrounding gas prior to the SN explosion. The mechanical wind energy is completely dwarfed by the star's radiative phase for any progenitor of an SN.

We consider how the UV propagates in each of the four halos for both clumpy and smooth media that undergo cooling. This is a complex process, and so we include access to 3D movies that reveal how the ionization front develops. These can be found at: http://miocene.anu.edu.au/smallgalaxy. We show a time sequence of the full preionization phase and the subsequent SN explosive phase in two ways: (i) fixed reference frame; (ii) rotating reference frame. There are two other aids to visualization. First, we include half-tone figures that show critical time steps for the different models: M60 (Figures 12 and 13), M65 (Figures 14 and 15), and M70 (Figures 16 and 17). Second, we include color composite images in Figure 18 that directly compare the information embodied in Figures 12–17, i.e., the neutral and ionized components as the nebula develops before and after the SN event.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A summary of gas mass and energy retention after the preionization phase is given in Tables 5 and 6, respectively. The initial gas masses inside two cardinal radii, ${r}_{{\rm{s}}}$ and ${r}_{\mathrm{vir}}$, are included for each model. The M60, M65, and M70 models retain essentially all of the gas within ${r}_{\mathrm{vir}}$ for the cooling models, while the M55 case (not shown) is evacuated completely. For the adiabatic cases (not shown), the M60 halo is evacuated completely and the M65 halo appears to lose more than half the gas after 30 Myr. In this instance, as the gas reacts to the radiation field, it becomes photoheated, expands, smooths out, and drops in density as a result. The number of recombinations locally declines, and the radiation propagates into more of the surrounding gas.

Table 5.  Percentage of Initial Mass Retained inside r_s and ${r}_{\mathrm{vir}}$ for Cooling Models

Model	Inside rs		Inside ${r}_{\mathrm{vir}}$
Timea	Central	Off-Center	Central	Off-Center
M60H	${M}_{{\rm{g}},{\rm{s}}}$: 1.91E+04 ${M}_{\odot }$		${M}_{{\rm{g}},\mathrm{vir}}$: 1.22E+05 ${M}_{\odot }$
0 Myr	12.7	16.1	97.6	97.1
5 Myr	0.03	0.06	89.3	88.4
25 Myr	0.01	0.21	0.98	6.10
M65H	${M}_{{\rm{g}},{\rm{s}}}$: 6.61E+04 M${}_{\odot }$		${M}_{{\rm{g}},\mathrm{vir}}$: 3.74E+05 M${}_{\odot }$
0 Myr	50.3	78.8	100.0	99.4
5 Myr	2.42	41.7	99.7	97.9
25 Myr	13.8	33.9	81.5	73.3
M70H	${M}_{{\rm{g}},{\rm{s}}}$: 2.38E+05 M${}_{\odot }$		${M}_{{\rm{g}},\mathrm{vir}}$: 1.16E+06 M${}_{\odot }$
0 Myr	102.6b	98.2	100.0	100.0
5 Myr	106.2b	92.8	100.2b	100.1b
25 Myr	121.5b	107.0b	102.6b	101.2b

Notes. Bold: <50% retained. The time at 0 Myr is when the SN event takes place, at which point the preionization phase ends. Entries in bold fall below 50%.

^a0 Myr = t_SNe, = 6.1 Myr from start of simulation. ^bIncreased mass after some infall.

Download table as:  ASCIITypeset image

In Figures 12–17, we observe a number of important facts about the preionization phase. For the centered models, the Strömgren radius of the H ii region falls outside (M60), onto (M65), and inside (M70) the scale radius. From the zoomed figures, we observe how the photoheating smooths out the substructure. The effect is so strong that we observe a cellular structure as the ionization fronts of many expanding knots overlap, which is particularly clear in the color composite image (Figure 18). In the same figure, we see that broad channels appear where the radiation is able to escape to the outer halo over a substantial solid angle. Even before the SN event, the expanding Strömgren sphere has produced a diffuse neutral shell of gas around the hot bubble cavity. In the off-centered models, the expanding Strömgren sphere compresses the core region gas, and the radiation is blocked from propagating over at least 2π sr. This becomes important in the post-SN phase because substantial gas begins to re-accrete to the center along this axis.

Our complete set of simulations (see Table 3) comprises 24 cases for SN events without preionization with the massive progenitor, and 6 cases where the ionization of the SN progenitor is included. All 30 cases were run at L00 and L01, with a subset of the cooling models run at L02 to allow adequate sampling of the developing SN shell (see Figure 5). L02 is always found to be adequate to resolve the metal-enriched shell, which starts as a diffuse hot envelope that becomes thinner under expansion. The models do not suffer the overcooling through inadequate resolution identified by Dalla Vecchia & Schaye (2012).

We discuss in detail only the clumpy, cooling models for the centered and off-centered SN cases where the progenitor ionization is included. We present these results for three different halo masses (M60, M65, M70).⁷ Tables 5 and 6 include the total gas mass and energy retained within r = ${r}_{{\rm{s}}}$ and r = ${r}_{\mathrm{vir}}$ after periods of t = 5 Myr and t = 25 Myr. The halos that retain less than 50% of their original mass or energy have their retention fractions shown in boldface.

Each of the models (M60, M65, M70) is discussed in more detail below. In Figures 12–17, we present time-dependent 2D slices of the gas distribution before and after the SN explosion. These are not summed along the third dimension and therefore do not illustrate the full extent of the ionized gas (see Figure 18). Unless otherwise stated, all computations are performed over a 648³ grid with two levels of the hierarchy (L00, L01) illustrated in Figure 5 (see Section 3.4). In a subset of cases, it was necessary to investigate a three-level hierarchy (L00, L01, L02) over a 1944³ grid, which we clarify below.

It is worthwhile to pay particular attention to Figure 19. Here we present 1D plots of the evolution of the gas mass (A–F) and energy (G–L) retention fraction over much longer timescales than shown in the half-tone figures (up to 25 Myr). In each plot, we present two sets of curves for retention within $r\;\leqslant$ ${r}_{\mathrm{vir}}$ (upper) and $r\;\leqslant$ ${r}_{{\rm{s}}}$ (lower).

Zoom In Zoom Out Reset image size

The left panels (A, C, E, G, I, K) show all combinations—clumpy (solid curves) versus smooth (dashed curves), off-centered (thin curves) versus centered (thick curves), adiabatic (red curves) versus cooling (black curves). These plots are particularly illustrative and allow us to discuss the role of internal and kinetic energy in the developing explosion. For all models, the adiabatic cases (red curves) confirm earlier work in that there is little or no gas retention (Mac Low & Ferrara 1999). In the right panels (B, D, F, H, J, L), we show only the cooling clumpy models for the centered and off-centered cases; the gas retention is clearly more favorable for the larger halo masses.

One important insight emerges that explains why we have arranged the figures in this way. In each row, we compare two distinct halo masses (e.g., M55 versus M60 in the top row) to illustrate the impact of preionization due to the SN progenitor. The left panels (e.g., A) have no preionization, in contrast to the right panels (e.g., B), which include this effect. Broadly speaking, the results on retention are similar. The same holds true when we compare M60 models without preionization to M65 models with preionization, and M65 models without preionization to M70 models with preionization. These are important limiting cases when one considers enrichment in halos that do not have an important preionization phase, e.g., for star formation histories limited to low mass and/or SNe I.

The time-dependent behaviors of the M60 central and off-centered explosions are illustrated in Figures 12 and 13, respectively. The evolutionary phases of an explosive wind are well understood analytically (Weaver et al. 1977), although we see additional substructure due to the complexity of the initial medium. The SN ignites within a fully ionized, diffuse nebula confined by a neutral gas shell. A dense radiative shell quickly forms and drives a shock into the neutral shell. The momentum-driven shell does not stall and essentially evacuates all gas within the halo after 15 Myr. In the off-centered case, one hemisphere is initially evacuated within 15 Myr, with the opposite hemisphere losing most but not all of its gas. At late times, there is some re-accretion from the far side, leading to <10% gas retention within the core regions. The energy input from the preionization phase, in combination with the SN explosive energy, conspires to doom the M60 halos to invisibility.

Numerical artifacts. In order to examine the impact of artifacts, this case is treated with a three-level grid. No such anomalies are seen in the preionization phase. In the SN event, numerical artifacts evident in the first 600 kyr near the origin (Figure 20) are caused by the carbuncle instability (see Sutherland 2010). This effect is amplified by strong cooling in the early stages (<125 kyr) and leads to the "clover leaf" pattern at 6.3 Myr inside r_s (L01). Our algorithm stabilizes the instability when the radiative thin shell becomes only partially resolved (<3 cells thick) at later times, such that any structure that evolves in the interaction between the spherical shell and the Cartesian grid remains at all later stages. The instability is strongest where the radiative shock front becomes orthogonal to the underlying grid, in this case along the Cartesian axes, giving rise to the boxy structures with apexes along the x- and y-axes seen here at 16 Myr (L00). After the local density has dropped so that the blast wave is no longer radiative, the structure then simply expands with the now adiabatic blast. The additional L2 level is needed to better resolve the early blast, and the remaining visible non-spherical structure has radius errors of approximately 10%–15%.

Zoom In Zoom Out Reset image size

Gas retention and energetics. We now discuss the issue of gas retention in the M60 models in the context of the energy (internal+kinetic) of the disturbed gas. Table 6 shows the amount of energy retained at critical radii, r = r_s and r = r_vir, immediately before and after the SN explodes, and at 5 and 25 Myr after the SN event. We show density slices through our clumpy M60 models over this time frame in Figures 20 and 21.

Zoom In Zoom Out Reset image size

We recall that during the preionization phase, radiation leaks out from the center in the first few Myr, leading to a fully developed, low-density nebula by 6 Myr (see Section 3.5). In Figure 20, the central SN explodes into a warm low-density nebula, superheats the gas through shock heating, and removes essentially all of the gas from the dark matter halo out to ${r}_{\mathrm{vir}}$ after 25 Myr. This is borne out in Table 5.

For the off-centered case in Figure 21, the progenitor SN is placed at the half gas-mass radius to reflect the fact that most stars form off-center. During the preionization phase, radiation leaks out from the center in the first few Myr (see Section 3.5), and a conic "zone of avoidance" is set up diametrically opposite the progenitor. By 6 Myr, a low-density nebula extends across most of the core region. The SN explodes into an off-centered, warm low-density nebula, superheats the gas through shock heating, and removes much of the gas asymmetrically from the dark matter halo out to ${r}_{\mathrm{vir}}$ after 25 Myr. The early zone of avoidance leads to far-side gas re-accreting into the core region after the SN event has passed, but less than 10% of the original gas remains within ${r}_{\mathrm{vir}}$ on long timescales.

We now discuss gas retention in the M65 models within the context of the mass and energy of the disturbed gas. The time-dependent behaviors of the M65 central and off-centered explosions are illustrated in Figures 14 and 15, respectively, before and including the explosion, and in Figures 22 and 23, respectively, after the SN event. To aid comparison, the projected neutral and ionized gases in Figures 14 and 15 are combined into a color composite in Figure 18.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Tables 5 and 6 show the percentage of mass and energy retained at r_s and r_vir just before and after the SN explosion, and after delays of 5 and 25 Myr. We see that more gas mass is retained in the core in the off-centered case compared to the central case after 25 Myr post-SN. The tabulated energy retentions are percentages relative to the binding energy of the gas bound by the dark halo. A system loses energy through either cooling (e.g., radiative, adiabatic) or mass loss; conversely, a system gains energy through either heating (e.g., shocks, UV radiation) or accretion. In the off-centered model, we see the effects of mass loss through the virial radius at the same time that the core is accreting gas.

In Figure 14, the SN ignites at 6.1 Myr within a fully ionized, diffuse nebula confined by a neutral gas shell. A dense radiative shell quickly forms and drives a shock into the neutral shell. A lot of energy is radiated away in the time frame 0.2–2.5 Myr. The momentum-driven shell stalls at 1 Myr. At this point, the shock slows down and becomes non-radiative. The drop in temperature leads to a rapid increase in density. The neutral shell then becomes energy driven and more diffuse as it propagates outward at about the halo escape speed and ultimately stalls at the virial radius. Most of the gas just survives in the M65 progenitor ionization + SN model.

In Figure 15, the off-centered H ii region runs out down the halo density gradient rapidly on the side near the star, but struggles to ionize through the core to the opposite side. The geometry allows the ionization to affect a larger volume in the outer halo than indicated by the integrated analysis in Figure 7. The SN occurs again in a smoothed region and reaches a significant fraction of the core radius by 6.5 Myr.

In these models, no more than half the gas survives the SN explosion. The surviving gas is enriched with the specific abundance yields of the discrete SN event. Our highest-resolution simulations reveal why cooling is so effective in retaining gas compared to other factors. In the early stages, the super-hot metal-enriched SN ejecta exhibit strong cooling, leading to much of the explosive energy being lost. This remains true regardless of the assumed metallicity floor at the time of the explosion.

In Figures 22 and 23, we explore the post-SN phase over larger radial scales (L00) and on longer timescales. The UV radiation from the SN progenitor escapes along some sight lines over the full spatial extent of the models. In both cases, the development of the SN-driven bubble is well defined. A thick H i shell forms by 10 Myr ahead of the thin metal-enriched shell, and this material starts to re-accrete to the core by 25 Myr. In the off-centered case, the central column densities are higher because the shock front wraps around the core and compresses it.

We now consider the retention of gas and energy in the M70 models. The H i and H ii column densities before and just after the SN are shown in Figures 16 and 17, while Figures 24 and 25 show the H i density evolution after the SN. For ease of comparsion, Figures 16 and 17 are also shown on the color composite image of Figure 18. The M65 models walk the tightrope between survival and disruption, retaining half their baryons at best. Tables 5 and 6 show percentages of gas and energy retained just before and after the explosion, as well as after 5 and 25 Myr. The M70 models clearly survive a single SN event, losing very little gas, and by 25 Myr contain more gas than in their initial state as a result of accretion. As a result of cooling, the energies within both ${r}_{{\rm{s}}}$ and ${r}_{\mathrm{vir}}$ are much lower than the initial values at both 5 and 25 Myr, suggesting that the gas is not flowing outward.

Zoom In Zoom Out Reset image size

For the M70 models, for the higher total gas mass, more than a single event could explode at a given epoch. But this is a complication that we do not consider here because the combined impact of the two events will depend on their masses (one is likely to be in the low-mass limit) and locations within the halo. In Figure 16, the central SN ignites at 6.0 Myr. The photon flux is insufficient to ionize a significant part of the core volume, and the H ii region is stifled. In the off-centered case in Figure 17, the off-centered H ii region is able to reach the edge of the core, and the falling density slope to the outer halo ensures that the photoionization fronts accelerate outward on the side containing the star. The SN is able to just reach the scale radius by 6.5 Myr, but has only a small effect on the bulk of the halo.

The post-SN phase is shown in Figures 24 and 25. For the central SN model, the bubble at t = 2.5 Myr is smaller than in the M65 case and sweeps up denser gas. By 10 Myr the bubble collapses and dense gas accretes back to the center. In the off-centered case, the behavior is similar, although a small amount of radiation reaches the edge of the grid on the side nearest to the SN, while the gas on the opposite side of the halo is largely unaffected. As in the M65 off-centered case, the central densities are increased at 2.5 Myr when the shock compresses the core; however, the effect is weaker than in the M65 model because the shock is unable to wrap around the core.

Zoom In Zoom Out Reset image size

Our models consider dwarf galaxies to be isolated, but in practice the vast majority reside in the extended halo of a host galaxy, albeit beyond the virial radius. In the M65 and M70 models, there is clear evidence of the expelled gas beginning to fall back after 20 Myr. A proper consideration of this gas requires that we consider the environment of the halo, which is beyond the scope of the current work. We leave this to a later study.

As described in the preceding sections, the summary plots for the single SN event models in Figure 19 (supported by Figures 12–17) reveal the following.

• 1.  
  For all M55 models and M60 adiabatic models, essentially none of their original gas mass remains within a scale radius at 25 Myr. Thus, subsequent star formation does not occur.
• 2.  
  For all M60 cooling models, only a few percent of the original gas mass remains within a scale radius at 25 Myr. No significant star formation can occur.
• 3.  
  For the M65 models, all of the cooling scenarios retain 70%–80% of the original gas inside the virial radius at 25 Myr. This is sufficient for ongoing star formation. None of the adiabatic models retain more than 20% of the original gas.
• 4.  
  In the borderline M65 cooling cases, off-centered models retain 2–3 times more gas within a scale radius, but the SN metals mix less efficiently with this gas.
• 5.  
  For the M70 models, all of the gas survives an SN event and its radiative progenitor. We stress that we have not considered clustered multi-SN events, which can occur in these models.
• 6.  
  Cooling is the most important factor governing gas retention within the virial radius at late times.
• 7.  
  Clumpiness is far more important in the adiabatic case than in the cooling case for gas retention.

The least massive galaxies lie at the bottom of the CDM hierarchy. Our focus here is not their role as building blocks in galaxy evolution, but what they can tell us about the yields of the first stars prior to the reionization epoch, and the following generations during and immediately after reionization (e.g., Bland-Hawthorn et al. 2010; Norris et al. 2010; Frebel & Bromm 2012; Ritter et al. 2014). At least some of the first galaxies must have resided within low-mass halos. We anticipate that such systems, if they can be identified in the present-day universe, are likely to carry the most ancient chemical signatures because of their age and because of the relatively low number of enrichment events over cosmic time.

If we can identify such objects, it should be possible to obtain a relatively clean signature of the primordial yields from the earliest stellar generations. The star formation efficiency is greatly diminished in low-mass galaxies, thereby reducing the total number of enrichment events. There are likely to be instances where a large fraction of the metals are blown out of the galaxy, and this may obscure the signatures of earlier generations of stars. This can occur when the relatively rare SNe that go off in these galaxies are correlated in time. But we would still expect to find the signatures of a single enrichment event in at least some intrinsically low-mass galaxies.

We now discuss our models in light of two recent discoveries: UFDs (near field) identified in all-sky imaging surveys (Simon & Geha 2007; Kirby et al. 2008), and very metal-poor DLAs (far field) identified along quasar sight lines (Erni et al. 2006; Cooke et al. 2010), although a more extensive discussion is given elsewhere (Webster et al. 2014, 2015a, 2015b). We refer the reader to Figures 1 and 2, which show that slow evolution is expected between the formation redshift at z = 10, the epoch of the DLAs (z = 2–4), and present-day UFDs. The overall effect is additional dark matter that is mostly accreted to the outer envelope of the dwarf, rather than onto the central cusp (Figure 4).

In recent years, the Sloan Digital Sky Survey has doubled the number of known dwarf galaxies that orbit the Galaxy (for a review, see Tolstoy et al. 2009), particularly at the low-mass end. These newly discovered dark-matter-dominated galaxies are uniformly very old and cover a wide range in luminosity down to ∼10³ ${L}_{\odot }$. Unlike globular clusters, they show a wide range of eccentricities, which may indicate that some are undergoing tidal disruption (Martin et al. 2008). Objects less luminous than 10⁵ ${L}_{\odot }$ are referred to as UFDs.

There has been extensive discussion concerning the integrated masses of UFDs as observed today (Simon & Geha 2007; Martin et al. 2008; Strigari et al. 2008; Bovill & Ricotti 2009; Wolf et al. 2010; Frebel et al. 2014). Some studies quote a characteristic UFD mass within 300 pc, i.e., ${{\rm{M}}}_{300}\approx {10}^{7}$ $\;{M}_{\odot }$ (Strigari et al. 2008; Walker et al. 2009), although this is an arbitrary physical scale that is not easily related to well-known scaling relations. The characteristic virial mass ${M}_{\mathrm{vir}}$ today is highly uncertain, with estimates ranging from 10⁸ to 10⁹ ${M}_{\odot }$. For a dark halo described by the NFW profile, the mass within a fixed radius gives the mass within any other radius, assuming that this is an accurate description of the underlying gravitational potential. For our work, we need to know the virial mass when the UFD first formed. The halo parameters in cosmological simulations have a strong dependence with redshift, making any estimates rather uncertain (Power et al. 2003; Duffy et al. 2008). We have adopted ${M}_{\mathrm{vir}}$ ≲ 10⁷ ${M}_{\odot }$ because the anticipated total population of UFDs may be too large to support larger progenitor masses.

The simulated "observables" from our models are broadly consistent with the most metal-poor UFDs (Webster et al. 2014, 2015a). The data for these systems are summarized by Frebel & Bromm (2012, Figure 1). The metallicity range for our models is −1.0 < [Fe/H] < −4.0 with a mean [α/Fe] enhanced by a factor of a few with respect to the solar value depending on the IMF at the time of formation (Woosley & Weaver 1995). We have begun to consider the r- and s-process signatures in the gas, work that is to be presented in later papers. The internal kinematics of the baryons are less than 4 km s⁻¹ within the inner scale radius, consistent with observations (Walker et al. 2009).

Brown et al. (2014) find that the star formation history of UFDs is marked by one or two bursts of star formation before truncation by a global event. By construction, such a history is consistent with our models. UFDs have stellar populations with mean metallicities in the range −2.5 < [Fe/H] < −2.0 (Kirby et al. 2013). Within the framework of our models, this requires (i) 1–2 SN events in a low-mass halo, specifically the M65 event-centered models (Section 5.4); and (ii) 3–10 SN events in a starburst within a dark matter halo an order of magnitude more massive (e.g., M70). These cases can give rise to measurably distinct chemical histories, star formation histories, and star cluster populations. For example, [α/Fe] is enhanced in starbursts relative to quiescent star-forming regions. Direct evidence for this comes from X-ray observations of starburst winds (Martin et al. 2002).

At the time of publication, we note that Oñorbe et al. (2015) have presented a new model for UFDs where the baryons reside within dark matter halos that are orders of magnitude more massive than we consider (${M}_{\mathrm{vir}}$ ≈ $3\times {10}^{9}$ $\;{M}_{\odot }$). Given that all of the baryons are retained within our M70 models, for the same initial gas fraction, we cannot easily distinguish their model from ours in terms of the predicted chemical abundances.

In summary, if UFD masses are largely unchanged since their early formation, we would expect these halos to show significant scatter in [Fe/H] at low metallicity, and for some of the oldest stars to constrain the yields of early stellar generations (e.g., Frebel & Bromm 2012).

An interesting development is the recent discovery of very metal-poor DLAs along QSO sight lines (Pettini et al. 2008; Penprase et al. 2010). If these are protogalactic structures that have recently formed from the intergalactic medium (IGM), there is the prospect of identifying the chemical imprint of early generations of stars (e.g., Pettini et al. 2002). Interestingly, two of these systems with [Fe/H] ∼ −3 may bear the hallmarks of early stellar enrichment (Erni et al. 2006; Cooke et al. 2010, 2015).

First, Erni et al. (2006) identify a DLA toward the QSO Q0913+072 (z = 2.785) with an iron abundance characteristic of the IGM at that redshift. The C, N, O, Al, and Si abundances show an odd–even effect reminiscent of the most metal-poor stars in the Galactic halo. This pattern is created in models where the neutron flux is low (e.g., Heger & Woosley 2002), presumably owing to the low overall metal abundance. A more striking signature is the strong [N/H] depletion, which Pettini et al. (2002) have argued is further evidence for a system that has recently formed from the IGM. Erni et al. (2006) argue that the abundances appear to agree with 10–50 ${M}_{\odot }$ zero-metallicity Pop III models.

Second, Cooke et al. (2015) identify half a dozen DLA objects that resemble the properties of UFDs with mean [Fe/H] in the range −2.5 < [Fe/H] < −2.0 and line-of-sight dispersion velocities of 1 < ${\sigma }_{\mathrm{LOS}}$ < 5 km s⁻¹ corrected for thermal and turbulent broadening. These physical parameters are in line with our M65 and M70 models (Webster et al. 2015b).

Finally, we examined whether the line-of-sight properties of our models would qualify them as DLAs. For all very metal-poor DLAs, Cooke et al. (2015) infer the total mass of neutral gas to be ${10}^{4}\lesssim {M}_{\mathrm{gas}}\lesssim {10}^{7}$ $\;{M}_{\odot }$ within a linear scale of $\lesssim 30-1300$ pc, although their UFD candidates are at the low end of these ranges. There are no sight lines in the M55 models, and few sight lines in the M60 models, that reach the DLA threshold. But the M65 and M70 models fare better, with a significant fraction of sight lines within a radius of $\sim {r}_{{\rm{s}}}$ being classified as a DLA-type spectrum. Furthermore, the inferred gas metallicity of [Fe/H] ∼ −3 falls well within our predicted range. This gas is mostly self-shielded from any external radiation and therefore does not require an ionization correction.

We use a lognormal distribution to describe the single-point statistics of the density field of our non-uniform ISM. The lognormal distribution is a skewed continuous probability distribution. Unlike the normal distribution, it has a non-zero skewness, variable kurtosis, and in general the mode, median, and mean are unequal.

The lognormal distribution appears to be a nearly universal property of isothermal turbulent media in experimental, numerical, and analytical studies (e.g., Pumir 1994; Nordlund et al. 1999; Warhaft 2000). Moreover, it is encouraging that the lognormal distribution is the limiting distribution for the product of random increments, in the same way that the normal distribution plays that role for additive random increments. It is thus compatible, at least conceptually, with a generic cascading process consisting of repeated folding and stretching.

With a lognormal distribution, which is on average isotropic, the natural logarithm of the ISM density field is a Gaussian that has a mean m and variance s². The probability density function for the lognormal distribution of the mass density ρ is

$$\begin{eqnarray}&&P(\rho )=\displaystyle \frac{1}{s\sqrt{2\pi }\rho }\mathrm{exp}\left[\displaystyle \frac{-{(\mathrm{ln}\rho -m)}^{2}}{2{s}^{2}}\right].\end{eqnarray} \tag{ 5 }$$

The mean μ and variance ${\sigma }^{2}$ of the density are given by

$$\begin{eqnarray}&&\mu =\mathrm{exp}[m+{s}^{2}/2].\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\sigma }_{F}^{2}={\mu }^{2}\;(\mathrm{exp}[{s}^{2}]-1).\end{eqnarray} \tag{ 7 }$$

In these simulations we adopt $\mu =1.0$, ${\sigma }_{F}^{2}=5.0$ as our standard lognormal distribution. With these parameters, densities below the mean compose one-quarter of the mass and occupy three-quarters of the volume, and the mean is approximately 20 times the mode. See Sutherland & Bicknell (2007) for more details.

The two-point structure of a homogeneous turbulent medium may be described in Fourier space. We denote the Fourier transform of the density $\rho ({\boldsymbol{r}})$ by $F({\boldsymbol{k}})$ (where ${\boldsymbol{k}}$ is the wavenumber vector). The isotropic power spectrum D(k) is the integral over solid angle in Fourier space of the spectral density $F({\boldsymbol{k}}){F}^{*}({\boldsymbol{k}})$. In three dimensions,

$$\begin{eqnarray}&&D(k)=\int {k}^{2}F({\boldsymbol{k}}){F}^{*}({\boldsymbol{k}})\;d{\rm{\Omega }}.\end{eqnarray} \tag{ 8 }$$

Even if the spectral density is anisotropic, the angular integral averages the spectral density into a one-dimensional function of ${\boldsymbol{k}}$ only. For a power-law dependence on ${\boldsymbol{k}},$ $D({\boldsymbol{k}})\propto {{\boldsymbol{k}}}^{-\beta }$, and $\beta =5/3$, the spectrum is referred to as Kolmogorov turbulence.

In order to simultaneously achieve lognormal single-point statistics and a power-law self-similar structure, we have implemented the practical method developed by atmospheric scientists, for constructing two- and three-dimensional terrestrial cloud models, which are used in radiative transfer calculations (Lewis & Austin 2002).

The remaining choice in this procedure is to select the range of wave numbers over which to generate the fractal, in particular the minimum wavenumber ${{\boldsymbol{k}}}_{\mathrm{min}}$, which determines the largest structure scale in the resulting fractal with respect to the spatial grid. The main domain of interest on the intermediate level in the simulations covers at least two dark matter scale radii in each model. We choose ${{\boldsymbol{k}}}_{\mathrm{min}}=8$ over this domain, corresponding to approximately two cloud structures per scale radius. On the upper level, the domain covers three times the extent, and the fractal is extended in a ${{\boldsymbol{k}}}_{\mathrm{min}}=24$ model. The lowest level is initially populated with a simple resampled version of the intermediate-level density field. Each level of the simulations is sampled in 216³ cells, and so even the coarse ${{\boldsymbol{k}}}_{\mathrm{min}}=24$ upper level has well-resolved small clouds.

Here we consider purely spherical distributions, neglecting rotational formation of disk structures. The reasons are twofold. First, flattening due to rotation will provide a natural axis for the escape of the SN energy and potentially greatly reduce the efficiency of the SN in evacuating the region of potentially star-forming gas. We are looking for a conservative estimate to the lowest-mass halo that retains gas, and so we wish to retain spherical symmetry, in order to maximize the potential interaction between the host ISM and the SNe. Second, it is less clear that these very small halos include significant angular momentum, and so a non-rotating model is the simplest applicable one.

Taking the mean global temperature of the ISM to be $\tilde{T}$ and the line-of-sight turbulent velocity dispersion to be ${\sigma }_{t}$, the mean density of gas in the potential, Φ, is given by (Sutherland & Bicknell 2007)

$$\begin{eqnarray}&&\displaystyle \frac{\rho (r)}{\rho (0)}=\mathrm{exp}{\left[-\displaystyle \frac{{\sigma }_{D}^{2}}{{\sigma }_{g}^{2}}{\rm{\Phi }}(r)\right]}^{\beta }.\end{eqnarray} \tag{ 9 }$$

where we take e = 0, ${\sigma }_{g}^{2}={\sigma }_{t}^{2}+k\tilde{T}/\mu {\rm{m}}$, and ${\sigma }_{D}$ is the velocity dispersion of the dark matter, assuming also that $\tilde{T}$ is the mean ISM temperature and ${\sigma }_{t}$ is constant throughout the ISM.

This mean global density distribution is multiplied by the fractal generated above, giving a clumpy distribution with the same mean global properties, and which is stable enough over the timescale of the simulations to remain static in the absence of an SN explosion.