Resolving the Formation of Protogalaxies. I. Virialization

TL;DR

This paper uses advanced simulations to analyze how baryonic matter in early protogalaxies virializes, revealing the roles of shocks, turbulence, and cooling in galaxy formation.

Contribution

It provides detailed insights into the virialization process of baryons in early protogalaxies, including the effects of cooling and turbulence, using high-resolution cosmological simulations.

Findings

Virial shocks occur near the virial radius in non-cooling halos.

Cooling causes the virial shock to shrink and filaments to penetrate deeper.

Turbulence in cooling halos is supersonic and driven by mergers.

Abstract

(Abridged) Galaxies form in hierarchically assembling dark matter halos. With cosmological three dimensional adaptive mesh refinement simulations, we explore in detail the virialization of baryons in the concordance cosmology, including optically thin primordial gas cooling. We focus on early protogalaxies with virial temperatures of 10^4 K and their progenitors. Without cooling, virial heating occurs in shocks close to the virial radius for material falling in from voids. Material in dense filaments penetrates deeper to about half that radius. With cooling the virial shock position shrinks and also the filaments reach scales as small as a third the virial radius. The temperatures in protogalaxies found in adiabatic simulations decrease by a factor of two from the center and show flat entropy cores. In cooling halos the gas reaches virial equilibrium with the dark matter potential through its turbulent velocities. We observe turbulent Mach numbers ranging from one to three in the cooling cases. This turbulence is driven by the large scale merging and interestingly remains supersonic in the centers of these early galaxies even in the absence of any feedback processes. The virial theorem is shown to approximately hold over 3 orders of magnitude in length scale with the turbulent pressure prevailing over the thermal energy. The turbulent velocity distributions are Maxwellian and by far dominate the small rotation velocities associated with the total angular momentum of the galaxies. Decomposing the velocity field using the Cauchy-Stokes theorem, we show that ample amounts of vorticity are present around shocks even at the very centers of these objects.

Figures

Tables

Table 1. Table 1 Simulation Parameters

Name	l	zend	Npart	Ngrid	Ncell	Cooling model
	[Mpc]
A0	1.0	15.87	2.22 $\times$ 107	30230	9.31 $\times$ 107 (4533)	Adiabatic
A6	1.0	15.87	2.22 $\times$ 107	40486	1.20 $\times$ 108 (4943)	H,He
A9	1.0	18.74	2.22 $\times$ 107	45919	1.21 $\times$ 108 (4953)	H,He,H2
B0	1.5	16.80	1.26 $\times$ 107	23227	6.47 $\times$ 107 (4023)	Adiabatic
B6	1.5	16.80	1.26 $\times$ 107	21409	6.51 $\times$ 107 (4023)	H,He
B9	1.5	23.07	1.26 $\times$ 107	20525	5.59 $\times$ 107 (3823)	H,He,H2

Table 2. Table 2 Halo Properties

Name	zcoll	Mtot	${\rho }_c$	$T_{\mathrm{vir}}$aaVirial temperatures are calculated with $\mu$ = 1.22 in all cases.	Tc	$⟨T⟩$
		[$M_{\odot }$]	[cm-3]	[K]	[K]	[K]
H2 Induced Collapse
A0	18.74	9.8 $\times$ 106	8.1	9200	10000	5700
A6	18.74	9.8 $\times$ 106	17	9200	7700	5500
A9	18.74	9.8 $\times$ 106	1.6 $\times$ 106	9200	590	5000
B0	23.07	6.2 $\times$ 106	13	8300	12000	6000
B6	23.07	6.2 $\times$ 106	15	8300	8900	5600
B9	23.07	6.7 $\times$ 106	3.0 $\times$ 106	8700	580	4200
Ly$\alpha$ Induced Collapse
A0	15.87	3.6 $\times$ 107	4.9	19000	17000	12000
A6	15.87	3.6 $\times$ 107	1.8 $\times$ 106	19000	8700	7300
B0	16.80	3.5 $\times$ 107	3.8	19000	31000	11000
B6	16.80	3.6 $\times$ 107	4.0 $\times$ 106	20000	9000	7500

Equations

$$\rho \frac{D{v}_i}{Dt}=-\frac{\partial }{\partial x_i}p+\rho g_i$$

$$\frac{1}{2}\frac{D^{2}I}{D{t}^2}=2\mathcal{T}+\mathcal{V}+3(\gamma -1)\mathcal{E}-\int p\mathbf{r}\cdot d\mathbf{S},$$

$$\mathcal{V}=-\frac{1}{2}G{\int }_V\frac{\rho (\mathbf{x})\rho ({\mathbf{x}}^{\prime })}{|\mathbf{x}-{\mathbf{x}}^{\prime }|}d\mathbf{x}d{\mathbf{x}}^{\prime },$$

$$r_{200}={\left[\frac{GM}{100{\Omega }_{\mathrm{CDM}}(z)H^2(z)}\right]}^{1/3},$$

$$t_{\mathrm{cool}}(r_{\mathrm{cool}})\equiv H(z{)}^{-1}$$

$$V_c=\sqrt{\frac{GM}{r_{200}}}\text{and}{T}_{vir}=\frac{\mu m_p{V}_c^2}{2k},$$

$${\mathcal{M}}_{turb}=\frac{v_{rms}}{c_s};c_s=\sqrt{\frac{\gamma kT}{\mu m_h}}$$

$$\beta =\frac{3(\gamma -1)\mathcal{E}+2\mathcal{T}}{E_s-\mathcal{V}}-1,$$

$$\frac{d{j}^2}{dr}>0\text{for rotational stability,}$$

$${\mathbf{u}}^{\prime }=\mathbf{u}+\frac{1}{2}\omega \times \mathbf{h}+\frac{1}{3}(\nabla \cdot \mathbf{u})\mathbf{h}+\frac{1}{2}\mathcal{D}\cdot \mathbf{h}$$

$$\dot{\rho }+\nabla (\rho \mathbf{u})=0,$$

$$\frac{1}{\rho }\frac{D\rho }{Dt}=-\nabla \cdot \mathbf{u}.$$

$$\frac{\partial \omega }{\partial t}+\nabla \times (\omega \times \mathbf{u})=\frac{1}{{\rho }^2}\nabla \rho \times \nabla p,$$

Full text

1. INTRODUCTION

The process of virialization is clearly fundamental to all scales of galaxy formation. Lynden-Bell (1967) demonstrated that violent relaxation occurs during the virialization of a dissipationless system, but does the equivalent occur for the baryonic matter? If it does, how it achieves virial equilibrium should be inherently different because of hydrodynamical effects and radiative cooling. Additionally, this would create a Maxwellian velocity distribution for the baryons as well. Turbulent velocities would exceed rotational ones. This would be at odds with the standard galaxy formation theories, which generally assume smooth solid body rotating gaseous distributions (e.g. Crampin & Hoyle, 1964; Fall & Efstathiou, 1980; Mo et al., 1998). The first occurrence of widespread star formation can be regarded as the commencement of galaxy formation, and its feedback on its host will affect all subsequent star formation. It is crucial to model the initial stage of galaxy formation accurately. Differences in initial configurations of a collapsing halo may manifest itself in different types of central luminous objects, whether it be a stellar disk (Fall & Efstathiou, 1980; Mo et al., 1998), a starburst (see §4 in Kennicutt, 1998, for a review), or a massive black hole (Bromm & Loeb, 2003; Volonteri et al., 2005; Spaans & Silk, 2006; Begelman et al., 2006). These differences may result from varying merger histories and the ensuing virial heating or turbulence generation of the new cosmological halo.

For galaxy clusters, cosmological virialization has been studied extensively (Norman & Bryan, 1999; Nagai & Kravtsov, 2003a; Nagai et al., 2003b; Schuecker et al., 2004; Dolag et al., 2005). It is customary to connect the velocity dispersion to a temperature through the virial theorem for a collisionless system, where the potential energy equals twice the kinetic energy. However X-ray observations and such cosmological simulations of galaxy clusters have indicated that turbulent energies are comparable to thermal energies. Central turbulent pressure decreases the density, but the temperature is largely unchanged. This leads to an increased entropy and a flatter entropy radial profile (Dolag et al., 2005) that is in better agreement with X-ray observations (e.g. Ponman et al., 1999). Simulations of merger dynamics suggest that turbulence is mostly generated in Kelvin-Helmholtz instabilities between bulk flows and virialized gas during minor mergers (e.g. Ricker & Sarazin, 2001; Takizawa, 2005). Alternatively turbulence can be generated by conduction (Kim & Narayan, 2003; Dolag et al., 2004) or acoustic transport of energy (Norman & Bryan, 1999; Cen, 2005).

In standard galaxy formation models (Rees & Ostriker, 1977; Silk, 1977; White & Rees, 1978; Blumenthal et al., 1984; White & Frenk, 1991; Mo et al., 1998), gas shock-heats to the virial temperature as it falls into DM halos. These models succeed with considerable accuracy in matching various observables, such as star formation histories, galaxy luminosity functions, and the Tully-Fisher relationship (White & Frenk, 1991; Lacey & Silk, 1991; Cole et al., 1994, 2000). Galaxy formation models depend on the virial temperature, most notably through the cooling function that controls star formation rates and their associated feedback mechanisms. Atomic hydrogen and helium radiative cooling is efficient in halos with masses between $10^{8}$ and $10^{12}{M_{\odot}}$ as cooling times can be less than the dynamical time of the system, a condition that galaxy clusters do not satisfy. This strong cooling suggests that in galaxies thermal energy may be less important for virialization than turbulent kinetic energy. Motivated by the results of galaxy cluster turbulence, we investigate this potentially important role of turbulence and radiative cooling in galaxy formation, using a series of high resolution numerical simulations of protogalactic halos in this work.

We consider kinetic energy and pressure forces in our virial analysis of protogalactic halos (cf. Shapiro et al., 1999; Iliev & Shapiro, 2001). This allows us to investigate the equilibrium throughout the entire halo and determine the importance of each energy component in the virial theorem. Kinetic energy can be decomposed into radial and azimuthal motions along with turbulence, which affects the collapse of gas clouds primarily in three ways. First as seen in galactic molecular clouds, turbulence plays an integral part in current theories of star formation as the density enhancements provide a favorable environment for star formation (Larson, 1979, 1981; Myers, 1999; Goldman, 2000). Second if the turbulence is supersonic, gas dissipates kinetic energy through radiative cooling, which aids the gaseous collapse (Rees & Ostriker, 1977). Conversely turbulent pressure adds an additional force for the collapsing object to overcome and can delay the collapse into a luminous object. Last, turbulence provides an excellent channel for angular momentum transport as the halo settles into rotational equilibrium to satisfy Rayleigh’s inviscid rotational stability argument (Rayleigh, 1920; Chandrasekhar, 1961) in which the specific angular momentum must increase with radius. Cosmological hydrodynamic simulations have just begun to investigate angular momentum transport within turbulent collapsing objects, and turbulence seems to play a large role in segregating low (high) angular momentum gas to small (large) radii (Norman & Bryan, 1999; Abel et al., 2002; Yoshida et al., 2006b).

We study idealized cases of structure formation where stellar feedback is ignored because it provides a convenient problem to focus on the interplay between cosmological merging, hydrodynamics, and cooling physics during the assembly of early halos. Some of the discussed physical principles should, however, be applicable to galaxies of all masses. These simulations provide the simplest scenario to which we can incrementally consider further additional physics, such as H2 and HD cooling physics (Saslaw & Zipoy, 1967; Palla et al., 1983; Flower et al., 2000), primordial stellar feedback (Whalen et al., 2004; Kitayama et al., 2004; Alvarez et al., 2006; Yoshida et al., 2006a; Abel et al., 2007), metal enrichment from primordial stars (Heger & Woosley, 2002; Tumlinson, 2006), AGN feedback (Springel et al., 2005; Kuhlen & Madau, 2005), and “normal” metal-enriched star formation (see Larson, 2003, for a review).

We present a suite of adaptive mesh refinement simulations that are described in §2. Then we analyze the local virial equilibrium and shocks in halos in §3. There we also differentiate between infall through voids and filaments and its associated virialization. We discuss the situations in which virial heating and turbulence occur. Next in §4, we decompose the velocity distribution in principle axes to explore virialization in both the DM and baryonic components. Furthermore we decompose velocities into shear and compressible flows to study turbulent flows in the virialized gas. We discuss the implications of these results on star and galaxy formation in §5. Finally we summarize in the last section.

2. THE SIMULATIONS

To investigate protogalactic ($T_{\rm{vir}}$ $>10^{4}$ K) halo virialization in the early universe, we utilize an Eulerian structured, adaptive mesh refinement (AMR), cosmological hydrodynamical code, Enzo111See http://http://lca.ucsd.edu/portal/software/enzo (Bryan & Norman, 1997, 1999; O’Shea et al., 2004). Enzo solves the hydrodynamical equations using the second order accurate piecewise parabolic method (Woodward & Colella, 1984; Bryan et al., 1995), while a Riemann solver ensures accurate shock capturing with minimal viscosity. Additionally Enzo uses an adaptive particle-mesh $n$-body method to calculate the dynamics of the collisionless dark matter particles (Couchman, 1991). Regions of the simulation grid are refined by two when one or more of the following conditions are met: (1) Baryon density is greater than 3 times $\Omega_{b}\rho_{0}N^{l(1+\phi)}$, (2) DM density is greater than 3 times $\Omega_{\rm{CDM}}\rho_{0}N^{l(1+\phi)}$, and (3) the local Jeans length is less than 4 cell widths. Here $N=2$ is the refinement factor; $l$ is the AMR refinement level; $\phi=-0.3$ causes more frequent refinement with increasing AMR levels, i.e. super-Lagrangian behavior; $\rho_{0}=3H_{0}^{2}/8\pi G$ is the critical density; and the Jeans length, $L_{J}=\sqrt{15kT/4\pi\rho G\mu m_{H}}$, where $H_{0}$, $k$, T, $\rho$, $\mu$, and $m_{H}$ are the Hubble constant, Boltzmann constant, temperature, gas density, mean molecular weight in units of the proton mass, and hydrogen mass, respectively. The Jeans length refinement ensures that we meet the Truelove criterion, which requires the Jeans length to be resolved by at least 4 cells on each axis (Truelove et al., 1997).

We conduct the simulations within the concordance $\Lambda$CDM model with WMAP first year parameters (WMAP1) of $h$ = 0.72, $\Omega_{\Lambda}$ = 0.73, $\Omega_{M}$ = 0.27, $\Omega_{b}$ = 0.024$h^{-2}$, and a primordial scale invariant ($n$ = 1) power spectrum with $\sigma_{8}$ = 0.9 (Spergel et al., 2003). $h$ is the Hubble parameter in units of 100 km s*-1* Mpc*-1*. $\Omega_{\Lambda}$, $\Omega_{M}$, and $\Omega_{b}$ are the fractions of critical energy density of vacuum energy, total matter, and baryons, respectively. Last $\sigma_{8}$ is the rms of the density fluctuations inside a sphere of radius 8$h^{-1}$ Mpc.

Using the WMAP1 parameters versus the significantly different WMAP third year parameters (WMAP3; Spergel et al., 2006) have no effect on the evolution of individual halos as are considered here. However these changes play an important role in statistical properties. For example, halos with mass $10^{6}{M_{\odot}}$ at redshift 20 correspond to $2.8\sigma$ peaks with the WMAP1 but are $3.5\sigma$ peaks for WMAP3. The $\Omega_{M}$/$\Omega_{b}$ ratio also only changed from 6.03 to 5.70 in WMAP3.

We also have verified that there is nothing atypical about the mass accretion rate histories of the objects we study. The mass accretion history of these objects exhibit smooth growth during minor mergers and accretion and dramatic increases when a major merger occurs. This behavior is consistent with typical halo assemblies in extended Press-Schechter calculations (Bond et al., 1991; Bower, 1991; Lacey & Cole, 1993, 1994; van den Bosch, 2002) and cosmological numerical simulations (e.g. De Lucia et al., 2004; Gao et al., 2005). The mass accretion histories in our simulations are well described by the fitting function of van den Bosch with $M_{0}=3\times 10^{7}{M_{\odot}}$, $z_{f}=17$, and $\nu=12.5$. We also compare our data against the mass accretion histories of Gao et al., who tested their data against an extended Press-Schechter calculation of the growth history of the halos. We find no major discrepancies between the two histories.

The initial conditions of this simulation are well-established by the primordial temperature fluctuations in the cosmic microwave background (CMB) and big bang nucleosynthesis (BBN) (Hu & Dodelson, 2002; Burles et al., 2001, and references therein).

We perform two realizations with different box sizes and random phases. In the first simulation (simulation A), we set up a cosmological box with 1 comoving Mpc on a side, periodic boundary conditions, and a 1283 top grid with three nested child grids of twice finer resolution each. The other simulation is similar but with a box side of 1.5 comoving Mpc (simulation B). We provide a summary of the simulation parameters in Table 1. These volumes are adequate to study halos of interest because the comoving number density of $>$104 K halos at $z=10$ is $\sim$6 Mpc*-3* according to an ellipsoidal variant of Press-Schechter formalism (Sheth & Tormen, 2002). We use the COSMICS package to calculate the initial conditions at $z$ = 129 (119)†††To simplify the discussion, simulation A will always be quoted first with the value from simulation B in parentheses. (Bertschinger, 1995, 2001), which calculates the linear evolution of matter fluctuations. We first run a dark matter simulation to $z=10$ and locate the DM halos using the HOP algorithm (Eisenstein & Hut, 1998). We identify the first dark matter halo in the simulation with $T_{\rm{vir}}$ $>$ 104 K and generate three levels of refined, nested initial conditions with a refinement factor of two, centered around the Lagrangian volume of the halo of interest. The nested grids that contain finer grids have 8 cells between its boundary and its child grid. The finest grid has an equivalent resolution of a 10243 unigrid. This resolution results in a DM particle mass of 30 (101) ${M_{\odot}}$ and an initial gas resolution of 6.2 (21) ${M_{\odot}}$.

Enzo employs a non-equilibrium chemistry model (Abel et al., 1997; Anninos et al., 1997). We conduct three simulations for each realization with (i) the adiabatic equation of state with an adiabatic index $\gamma=5/3$, (ii) a six species chemistry model (H, H+, He, He+, He++, e-), and (iii) a nine species chemistry model that adds H2, H${}_{2}^{\rm+}$, and H- to the six species model. In the nine species model, we use the molecular hydrogen cooling rates from Galli & Palla (1998). These models are differentiated in the text by denoting 0, 6, and 9, respectively, after the simulation name (e.g. simulation B0). Compton cooling and heating of free electrons by the CMB and radiative losses from atomic and molecular cooling are also computed in the optically thin limit.

To restrict the analysis to protogalactic halos in the H2 models, we suppress H2 formation in halos that cannot undergo Ly$\alpha$ cooling by reducing the residual electron fraction to $10^{-12}$ instead of a typical value of $\sim 10^{-4}$ only at the initial redshift (Shapiro et al., 1994). This mimics an extreme case where all H2 is dissociated by an extremely large radiation background, and the halo can only collapse and form stars when free electrons from ionized hydrogen can catalyze H2 formation.

We end the simulations with non-equilibrium cooling when the gas begins to rapidly cool and collapse. We choose a final resolution limit of $\sim$3000 (4000) proper AU, corresponding to a refinement level of 15. We end the adiabatic simulations at the same redshift. In a later paper, we will address the collapse of these halos to much smaller scales.

3. VIRIAL ANALYSIS

The equation of motion for an inviscid gas in tensor notation reads:

[TABLE]

where $D/Dt=\partial/\partial t+v_{j}\partial/\partial x_{j}$ is the total derivative. Here $v$ is velocity; $p$ is pressure; $\rho$ is density; and $g=\nabla\Phi$ where $\Phi$ is the gravitational potential. From this Chandrasekhar & Fermi (1953) derived the general virial theorem for a region contained within a surface ${\bf S}$, in scalar form,

[TABLE]

where

[TABLE]

$\mathcal{T}=\frac{1}{2}\int\rho{\bf v}^{2}d{\bf x}$, $I=\int\rho{\bf x}^{2}d{\bf x}$, and $\mathcal{E}=\int\varepsilon\;d{\bf x}$ denote the gravitational potential energy, the trace of the kinetic energy and inertia tensor, and the total internal thermal energy, respectively. The surface term $E_{s}$ (the last term in eq. [1]) is often negligible in the outer regions of the halo. The system is not necessarily in virial equilibrium if $\ddot{I}=0$, but the time-averaged quantity is zero when the entire system is in virial equilibrium. A system is expanding or contracting whether $\ddot{I}$ is positive or negative, respectively, based on energy arguments. Ballesteros-Paredes (2006) gives counterexamples to this simple interpretation. However, in the cases presented here, spherically averaged radial velocities are always negative.

We define the halo as the material contained in a sphere with a radius $r_{200}$ enclosing an average DM overdensity of 200 and as such relates to mass by

[TABLE]

where $M$ is the mass of the halo, $\Omega_{\rm{CDM}}(z)$ is evaluated at a redshift $z$, and $H$ is the Hubble parameter at $z$. The region where the cooling time is shorter than a Hubble time is denoted as the cooling radius $r_{\rm{cool}}$,

[TABLE]

(White & Frenk, 1991). Mass and radius define a circular velocity and virial temperature, which are

[TABLE]

for a singular isothermal sphere (see Bryan & Norman, 1998, with $\beta$ = 1 and $\Delta_{c}$ = 200). 0.59 and 1.22 in units of the proton mass are appropriate values for $\mu$ for the fully ionized and completely neutral states of a primordial hydrogen and helium mixture of gas, respectively. We use $\mu$ = 1.22 throughout this paper. We note that Iliev & Shapiro (2001) considered non-singular, truncated isothermal spheres, and the resulting virial temperature is $\sim$15% lower than the one calculated in equation (5).

For $\gamma=5/3$, $T_{\rm{vir}}$ is the temperature at which an ideal adiabatic gas reaches virial equilibrium with the specified potential. Please note that for an isothermal gas where $\gamma$ is close to unity virial equilibrium is established between the turbulent energies and the gravitational potential as the $3(\gamma-1)\mathcal{E}$ term in equation (1) goes to zero.

3.1. Local Analysis

We evaluate the terms of equation (1) with respect to radius (i.e. the volume contained in a radius $r$). Figure 1 illustrates the radial structure of (a) the turbulent Mach number,

[TABLE]

(b) turbulent and (c) thermal energies per $m_{h}$, and (d) a “virialization” parameter333This is a modified version of the $\beta$ used in Shaw et al. (2006) to account for kinetic energies and so that it does not diverge when $\mathcal{V}\rightarrow 0$ in the center. It still has the same behavior of $\beta\rightarrow 0$ as $\ddot{I}\rightarrow 0$.

[TABLE]

of the adiabatic and radiative cooling simulations when the cooling halo collapses. Here $v_{rms}$ is the three-dimensional rms velocity and is assessed using the gas velocities relative to the mean gas velocity of each spherical shell. In the top row of Figure 1, we also plot the Mach number, using $v_{rms}$ with respect to the mean velocity of the gas within $r_{200}$. In the six and nine species simulation, this occurs at zc,H2 = 18.7 (23.4) and zc,Lyα = 15.9 (16.8), respectively. The radial profiles are centered on the densest point in the simulation with the collapsing halo. Several properties of the most massive halo in each simulation are detailed in Table 2. The sections of the Table compare the halo in the adiabatic, Ly$\alpha$, and H2 simulations.

Figure 2 shows the mass-weighted radial profiles of gas mass enclosed and gas density at z = zc,Lyα in the adiabatic and cooling cases. Both realizations are remarkably similar. Halos in the adiabatic case have a central core with a radius $\sim$50 pc and gas density of $\sim$5.0 (3.5) cm*-3*. Core densities in simulation A are slightly higher than simulation B, which has larger thermal and turbulent pressures (see Figure 1). With radiative cooling, gas infalls rapidly as it cools and undergoes a self-similar collapse with $\rho\propto r^{-12/5}$.

3.1.1 Virial Radius

We define the virial radius $r_{\rm{vir}}$ when $\beta=0$ and $d\beta/dr<0$. When we radially average the Ly$\alpha$ halo, $r_{\rm{vir}}$ = 419 (787) pc in the adiabatic simulations where the corresponding $r_{200}$ value is 615 (576) pc.

A well defined shock exists between the interface between voids and the halo. This material shock-heats to $T_{\rm{vir}}$ and virializes at a radius comparable to $r_{200}$ in the adiabatic cases. When we include radiative cooling, this radius decreases everywhere around the halo-void virial shock and is low as $r_{200}$/2. In contrast to the voids, the filamentary gas shock-heats at an even smaller radius. Dekel & Birnboim (2006) also studied the stability of cold inflows within a hot virialized medium and found similar results. Figures 2 and 3 illustrate these changes. At $r_{200}$, densities in the adiabatic case begin to increase more rapidly than the cooling case as material accretes at the virial shock. No significant increase in $d\rho/dr$ is seen in the cooling case, indicative of a self-similar collapse.

3.1.2 Adiabatic Model

We start with the discussion of the adiabatic model as it is the simplest case and later compare the calculations with radiative cooling to this model. Virialization should transfer potential energy to kinetic energy that dissipates in shocks to thermal energy, which is the implication of the dissipationless virial theorem. The solid lines in Figure 1 represent the energies in adiabatic models. The physics illustrated in this Figure are as follows:

1. Thermal energy— The gas shock-heats to the virial temperature at the virial shock. Virial heating continues with decreasing radius as the surface term becomes significant in the interior. The resulting central temperature of the halo is 10000 (12000) K, which is 1.2 (1.5) $T_{\rm{vir}}$, at a redshift of zc,H2 when the H2 model collapses. At the time (z = zc,Lyα) of collapse caused by Ly$\alpha$ cooling, the central temperature is 17000 (31000) K, corresponding to 0.9 (1.6) $T_{\rm{vir}}$.

2. Kinetic energy— It increases along with the thermal energy during virialization. The gas is generally turbulent, appearing as a velocity dispersion with a bulk radial inflow. At $r_{200}$, the kinetic energy is equivalent to the thermal energy, $\mathcal{T}/\mathcal{E}\sim 1$. This ratio steadily drops toward the center, where $\mathcal{T}/\mathcal{E}\sim 1/3$. This decrease in kinetic energy is apparent in all the calculations except simulation B at zc,Lyα, increasing by a factor of two in the center.

3. Turbulent Mach number— At $r_{200}$, the turbulent Mach number $\mathcal{M}_{turb}$ is maximal and varies from 1–3 in all adiabatic simulations. $\mathcal{M}_{turb}$ decreases to subsonic values $\sim$0.15 but never below in the interior. Note that $\mathcal{M}_{turb}$ does not increase as the turbulent energy towards the center in simulation B because of the also growing sound speed there.

4. Virialization parameter— Virial equilibrium is quantified by the virialization parameter $\beta$, where the collapse is retarding or accelerating when it is negative or positive, respectively. At zc,H2 and zc,Lyα and in both simulations, $\beta$ is within 20% of being virialized ($\beta=0$). For comparison purposes, this corresponds to a halo having 80% of the required velocity dispersion for virialization in the dissipationless case. At $r_{\rm{vir}}$, $\beta$ is nearly zero which defines the virialized object. Within $r_{\rm{vir}}$, the values decrease to values around –0.1 but stays $\,\buildrel<\over{\sim}\,$ 0.

Characteristics of turbulence in our adiabatic models are similar to ones found in galaxy cluster simulations (Norman & Bryan, 1999; Dolag et al., 2005). Both groups find that turbulence provides $\sim$5–30% of the total pressure, i.e. $\mathcal{T}/(\mathcal{T}+\mathcal{E}$), in the cluster cores. Our protogalactic halos have $\sim$25% of the pressure in the turbulent form. Also the galaxy clusters in Norman & Bryan (1999) have comparable Mach numbers of $\sim$1.6 at $r_{\rm{vir}}$, $\sim$0.5 at $r_{\rm{vir}}$/3, and $\sim$0.3 in the core. These similarities suggest that virial turbulence is generated over a large range of mass scales.

3.1.3 Ly$\alpha$ Cooling Model

At zc,Lyα, halos in calculations with the H+He cooling model start to rapidly collapse. The dashed lines in the third and fourth columns of Figure 1 illustrate the energies of this model.

1. Thermal energy— Compared to the adiabatic models, the gas can radiatively cool through Ly$\alpha$ emission to T $\sim$ 8000 K within $r_{\rm{cool}}$ $\sim$ $r_{\rm{vir}}$. The entire halo is isothermal at this equilibrium temperature. Below this temperature, the cooling function of pristine gas drops by several orders of magnitude, and the gas can no longer cool efficiently. The thermal energy is $\sim$65% lower than the adiabatic case.

2. Kinetic energy— In response to the lesser thermal energy, the system tends toward virial equilibrium by increasing kinetic (turbulent) energy. The gravitational potential and surface terms do not appreciably change with the inclusion of radiative cooling. Turbulent energy within $r_{\rm{vir}}$ increases as much as a factor of 5 when compared to the adiabatic case.

3. Turbulent Mach number— The changes in thermal and kinetic energies equate to a increase of $\mathcal{M}_{turb}$ by a factor of 2–3 to values up to 1.5. The turbulence is supersonic in all cases at the virial shock, but when we include radiative cooling, this trait emanates inward as the halo begins to rapidly cool. When the central core becomes gravitationally unstable, the entire halo is supersonically turbulent.

4. Virialization parameter— The increased kinetic energies compensate for the loss in thermal energy and the halo remains in a similar virial state. This is apparent in the remarkably similar radial characteristics of $\beta$ in the adiabatic and H+He models of simulation A.

3.1.4 H2 Cooling Model

The collapses caused by H2 cooling at z = zc,H2 have very similar dynamics as the halos described in the previous section. The dashed lines in the first and second columns of Figure 1 illustrate the energies of this model.

1. Thermal energy— H2 cooling is efficient down to 300 K, so gas can depose a much larger fraction of its thermal energy. Inside $r_{\rm{cool}}$ $\sim$ 0.32 (0.19) $r_{\rm{vir}}$, thermal energies are only 5% of the values in the adiabatic models.

2. Kinetic energy— The turbulent energies must increase as in the Ly$\alpha$ case, and they increase by 93% (44%) on average inside $r_{\rm{cool}}$.

3. Turbulent Mach number— Similarly, $\mathcal{M}_{turb}$ increases up to a factor of 10 to become supersonic at values up to 3 throughout the halo. They are somewhat larger than the Ly$\alpha$ cases since H2 can cool to significantly lower temperatures than the virial temperature.

4. Virialization parameter— The virial equilibrium of the halos are also similar to the other models. $\beta$ smoothly transitions from nearly equilibrium at $r_{\rm{vir}}$ to an increased radial infall with $\beta$ = –0.2 at 70 pc. Then it increases to 0.4 inside 10 pc, which corresponds to the gas decelerating from the rapid infall as it encounters the central molecular cloud.

3.1.5 Model Summary

Baryons are close to virial equilibrium over three orders of magnitude in length scale by gaining both thermal and kinetic energies independent of cooling physics. Central temperatures of the adiabatic simulations are up to twice the nominal virial temperature. Similar to galaxy cluster studies, turbulence in the adiabatic model contributes $\sim$25% to the energy budget with Mach numbers $\sim$ 0.3 in the center. In cooling cases, atomic and molecular cooling inhibit virialization through heating, therefore the object must virialize by gaining kinetic energy up to five times the energy seen in the adiabatic models. This translates into the flow becoming supersonically turbulent with Mach numbers ranging from one to three.

3.2. Variations in the virial shock

Using the adiabatic invariant, $K=T/n^{2/3}$, which we label “entropy”, allows us to differentiate between gas accreting from voids and filaments. As a precaution, we note that $K$ is not an invariant when $\gamma$ varies; however, this is not the case in our simulations in which we permit molecular hydrogen cooling. Here molecular fractions remain low, $<10^{-3}$, and $\gamma\approx 5/3$ even in the densest regions. The top row of Figure 3 depicts the variance of $K$ with respect to radius in two-dimensional histograms, where the intensity of each pixel represents the mass having the corresponding $K$ and $r$. The middle and bottom rows display two-dimensional slices of $K$ and density, respectively, through the densest point in the halo. The virialized gas from the voids has low density and does not significantly contribute to the mass averaged radial profiles. Figure 3 illustrates this gas at $r\sim r_{200}$ and K\lower 3.00003pt\hbox{,\buildrel>\over{\sim},}10^{4.5} K cm2. The gas in filaments has lower entropy than the rest of the halo at $r>$ 150 pc and K\lower 3.00003pt\hbox{,\buildrel<\over{\sim},}10^{3.3} K cm2. In equation (1), the pressure in the surface term is the constant at a given radius. The accreting, denser, unshocked gas in filaments has lower temperatures than the more diffuse accreting gas. The gas remains cool until it shocks and mixes well inside $r_{\rm{vir}}$ and as small as \sim$$r_{\rm{vir}}/4 in the most massive filaments. Similar characteristics of cold accretion flows have been noted and discussed by Nagai & Kravtsov (2003a), Kereš et al. (2005), and Dekel & Birnboim (2006).

Entropy in the exterior of the halo differ little between adiabatic and cooling runs outside of $r_{\rm{cool}}$. But as the gas falls within $r_{\rm{cool}}$, it cools and condenses, which gives a lower entropy, and the $r$-$K$ histograms and entropy slices display this clearly. Another significant difference in the cooling simulations is the contraction of the virial shock by a factor of 1/3 when compared to adiabatic runs. This is caused by the contraction of the cooling gas. Here the cold filaments penetrate to even smaller radii. This is also evident in the radial density profiles of Figure 2.

3.3. Virial Heating and Turbulence

In order for a system to remain in virial equilibrium as it grows in mass, additional gravitational energy is balanced through two possible mechanisms: heating of the gas ($\dot{\mathcal{E}}>0$) and increasing the kinetic energy of the gas ($\dot{\mathcal{T}}>0$). We differentiate between two main cases of virialization by comparing the cooling time, $t_{\rm{cool}}$ = kT/n$\Lambda$, of the system to the heating time, $t_{\rm{heat}}$ = $T_{\rm{vir}}$/$\dot{T}_{\rm{vir}}\sim\frac{3}{2}M_{\rm{vir}}/\dot{M}_{\rm{vir}}$ in the case of rapid mass accretion. Birnboim & Dekel (2003), Dekel & Birnboim (2006), and Wang & Abel (2007) find that radiative cooling rates are greater than heating rates from virialization for halos with masses below $10^{12}{M_{\odot}}$.

1. Thermalization ($t_{\rm{cool}}$ $>$ $t_{\rm{heat}}$)— When no efficient radiative cooling mechanisms (e.g. H2, Ly$\alpha$, He I) exist, the system virializes by injecting energy into heat $\mathcal{E}$ and partly into kinetic energy $\mathcal{T}$. In the process, the halo becomes pressure supported and virialized. Traditional galaxy formation scenarios only consider this thermalization while neglecting the kinetic energy term of equation (1). However it is important to regard kinetic energy, even in adiabatic models, as the gas violently relaxes. Turbulence velocities are similar to the velocity dispersion of the system and contributes notably to the overall energy budget as seen in adiabatic cases in Figure 1.

2. Turbulence generation ($t_{\rm{cool}}$ $<$ $t_{\rm{heat}}$)— When a cooling mechanism becomes efficient, the system now dispenses its thermal energy and loses pressure support within $r_{\rm{cool}}$. The gas will cool to a minimum equilibrium temperature. As the cooling halo collapses and radial velocities increase, the gas still lacks enough kinetic and thermal energy to match the gravitational energy and surface term in equation (1). The gas becomes more turbulent in order to virialize. We see this in the second row of Figure 1, where turbulent energies are significantly increased well inside the halo in the cooling models as compared to the adiabatic calculations.

Through virial analyses, we have shown that turbulent energies are comparable, if not dominant, to thermal energies in galaxy formation. In the next Section, we further investigate the significance and nature of the turbulence through velocity distributions and decompositions in order to study any small-scale anisotropies in the internal flows.

4. VELOCITY DISTRIBUTIONS

In CDM cosmogony, collisionless dark matter dominates the gravitational potential and oscillates as it stabilizes. Lynden-Bell (1967) showed how a collisionless system undergoes violent relaxation if embedded within a rapidly time-varying potential. Individual mass elements do not conserve energy during violent relaxation, only the entire system conserves energy. This behavior randomizes the energies of the mass elements, and statistical mechanics makes the resulting energy (velocity) distribution to tend to Maxwellian. Furthermore, the system “forgets” its original configuration during virialization or the incorporation of a lesser halo. Later studies have inferred two baryonic scenarios of virialization. First, violent relaxation and the accompanying phase mixing also applies to the gaseous component (van den Bosch et al., 2002; Sharma & Steinmetz, 2005). In the other case, the gaseous component dissipates all turbulent motions and finally rigidly rotates as a solid body with a velocity appropriate to the overall spin parameter (e.g. Loeb & Rasio, 1994; Mo et al., 1998; Bromm & Loeb, 2003).

Figure 4 shows the velocity distributions of the dark matter and baryonic components of the halo at zc,H2 in the left column and zc,Lyα in the right column. It also overplots the simulations with adiabatic and radiative cooling. We plot the radial and tangential velocity distribution on the top and bottom rows, respectively. The velocities are taken with respect to the bulk velocity of the halo. We also transform the velocity components to align the z-axis and total angular momentum vector of the DM halo.

The radial velocity distributions at z = zc,H2 are approximately Maxwellian in both dark matter and gas with a skew toward infall. The infall distributions are shifted by $\sim$1 km s*-1* in the cooling case when compared to adiabatic. However at z = zc,Lyα, the effects of Ly$\alpha$ cooling become more prevalent in the halo when compared to H2 cooling, shifting the radial velocity distribution by $\sim$5 km s*-1* that is caused by faster infall. These distributions have two components that represent virialized gas and infalling gas in filaments. We further discuss this in the next Section.

The tangential velocity distributions are nearly Maxwellian in all cases except for the dark matter in simulation A at z = zc,Lyα (right panels in Figure 4a). This deviation from Maxwellian arises from two major mergers that occur between 25 and 85 Myr (z = 17–21) before the final collapse. The residual substructure from the major merger causes three distinct populations with Gaussian distributions centered at –0.2, +13.6, and –6.7 km s*-1* with $\sigma$ = 11.6, 4.2 and 3.6 km s*-1*, respectively. These distributions clearly do not resemble a solid body rotator, whose velocity distribution would contain all positive velocities. In other words the turbulent velocities exceed the typical rotational speeds.

Distributions in dark matter are broader than the gas in both simulations and collapse redshifts as expected because for the gas we only give the bulk velocities and do not add the microscopic dispersion (cf. van den Bosch et al., 2002; Sharma & Steinmetz, 2005).

4.1. Halo and Filament Contrasts

The dark matter velocity distributions are typical of a virialized system with the majority of the matter having a Maxwellian distribution with a dispersion corresponding to the main halo (Boylan-Kolchin & Ma, 2004; Dieman et al., 2004; Kazantzidis et al., 2004). Substructure appears as smaller, superposed Gaussians, which are stripped of its outer material as it orbits the parent halo. Dynamical friction acts on the substructure and decreases its pericenter over successive orbits, and the subhalo is gradually assimilated in the halo.

The filaments penetrate deep into the halo and provide mostly radial infall inside $r_{200}$. They do not experience a virial shock at $r_{200}$, and this contrast is apparent in the radial velocity distributions. When we restrict our analysis scope to the filaments (i.e. $r>150$ pc and $K<10^{3.3}$ K cm*-2*), the radial velocity distribution (Figure 4) is skewed toward infall, centered at –15 km/s, which is approximately the circular velocity of the halo. The rest of the gas outside of this region in $r-K$ space has already been virialized, shock heated, and roughly exhibits a Maxwellian distribution, centered at zero, with its associated substructures. Hence the mass in filaments dominate the radial velocity distributions at negative values in Figure 4.

4.2. Turbulence

Radial inflows can create turbulence in the halo. Filaments provide an influx material with distinct angular momentum. This gas virializes in the presence of an already turbulent medium that has a relatively high specific angular momentum at $r>r_{200}/4$. The Rayleigh inviscid instability criterion requires

[TABLE]

where $j$ is the specific angular momentum. If this is not satisfied, the system will become unstable to turbulence. The onset of turbulence can be delayed if viscosity were large enough so that Reynolds numbers are below the order of 102 or 103. However there are many modes of instability if equation (8) is not met, and even a gas with low Reynolds number will eventually become fully turbulent (Shu, 1992).

Velocity dispersions can characterize the general magnitude of a turbulent medium, but its local nature is better detailed by applying the Cauchy-Stokes decomposition,

[TABLE]

that decomposes the velocity field into bulk motion ${\bf u}$, vorticity ${\bf\omega}=\nabla\times{\bf u}$, expansion and contraction $\nabla\cdot{\bf u}$, and a distortion $\mathcal{D}$ without change in volume. Here the vector ${\bf h}$ describes the separation between gas parcels at position ${\bf x}$ and ${\bf x}^{\prime}$, and ${\bf u}^{\prime}$ is the velocity at ${\bf x}^{\prime}$.

We relate $\nabla\cdot{\bf u}$, which is plotted in Figure 5, to a convergence timescale through the continuity equation,

[TABLE]

that can be rewritten in terms of the total derivative $D/Dt$ as

[TABLE]

$D\rho/Dt$ describes density changes along the fluid flow lines, and the $1/\rho$ factor converts this change into an inverse timescale on which local densities e-fold. We denote $-(\nabla\cdot{\bf u})^{-1}$ as the “Lagrangian convergence timescale” (LCT). Converging flows ($\nabla\cdot{\bf u}<0$) are ubiquitous within the halo. On large scales, the smallest LCTs on the order of 20 kyr exist at the virial shock, adjacent to both the filaments and voids. In the cooling models, these timescales are also small in turbulent shocks well within $r_{200}$. Analogous to the dynamical time, typical shocked LCTs decrease toward the center as density increases.

The local magnitude and nature of turbulence is further illustrated by the vorticity ${\bf\omega}$, whose component perpendicular to the slice, $\omega_{y}$, is shown in Figure 6. The local rotation period, $4\pi/|\omega|$ is also helpful to quantify and visualize the nature of the flow. The vorticity equation reads

[TABLE]

where the source term is non-zero when the density and pressure gradients are not aligned, i.e. baroclinic vorticity generation. This occurs at and near shocks throughout the halo, regardless of radiative cooling. In the adiabatic models, vorticity exists even at modest but sufficient resolution in the large pressure-supported cores (see $\mathcal{M}_{turb}$ for $r<100$ pc in Figure 1) and generates a turbulent medium with $\mathcal{M}_{turb}$ = 0.3. In the cooling models, this large-scale vorticity is still present but increases in the collapsing core. As shocks become abundant in the center, we do not see any dampening of kinetic energy. Perhaps this mechanism maintains turbulent motions during virialization, even in the presence of dissipative shocks. In Figure 6, adjacent, antiparallel fluid flows, i.e. a sign change in $\omega_{y}$, are ubiquitous, which visually demonstrates that turbulence exists throughout the halo. The length scale of these eddies decrease with increasing density as with the LCTs.

Hence we believe significant turbulence generated during virialization should be present in all cosmological halos. The cooling efficiency of the gas, the total halo mass, and partly the merger history determines the magnitude of turbulence. We discuss some implications of virial turbulence in the following section.

5. DISCUSSION

We have investigated the virialization of early pre-galactic cosmological halos in this paper with a suite of AMR simulations with varying chemistry and cooling models and collapse epochs. When analyzing the local virial equilibrium of the halo, we do not assume that it is in equilibrium but explicitly calculate all of the relevant terms in the virial theorem. In both adiabatic and radiative cooling cases, we find that the kinetic (turbulent) energy is comparable, if not dominant, with the thermal energy. Turbulence appreciates as radiative cooling becomes efficient because thermal energy alone cannot bring the system into virial equilibrium. In this case, the gas attempts to virialize by increasing and maintaining its kinetic energy.

Besides violent relaxation, at least two other hydrodynamic processes will augment virial turbulence. The first occurs when radial inflow interacts with the virialized gas. Due to the Rayleigh criterion, the high angular momentum gas creates an instability when it is deposited by filaments at small radii. The second happens when minor and major mergers create Kelvin-Helmholtz instabilities and drives additional turbulence (e.g. Ricker & Sarazin, 2001; Takizawa, 2005). Our results show that this turbulence is acting to achieve close to virial equilibrium at all stages during assembly and collapse.

Virial turbulence may be most important in halos which can cool rapidly when compared to virial heating from mass accretion. Interestingly all halo masses below $\sim$10${}^{12}{M_{\odot}}$ that can cool by Ly$\alpha$ emission satisfy this condition (Birnboim & Dekel, 2003; Dekel & Birnboim, 2006; Wang & Abel, 2007).

Turbulence appears to mix angular momentum efficiently so that it redistributes to a radially increasing function, and thus only the lowest specific angular momentum material sinks to the center. This segregation allows a collapse to proceed as it were self-similar and basically non-rotating. Similar results have been reported in cosmological simulations of collapses of the first stars (Abel et al., 2000, 2002; Yoshida et al., 2003, 2006b; O’Shea & Norman, 2007). Here and in protogalactic collapses, the turbulent velocities become supersonic. One would expect even higher Mach numbers in larger potential wells that still have $t_{\rm{cool}}$ $<$ $t_{\rm{dyn}}$.

The inclusion of the surface term allows us to study the virial equilibrium in the halo’s interior where the gravitational potential is not influential. Here our simulations show thermal and kinetic energies balancing the surface term and potential energy to achieve virial equilibrium. Before cooling is efficient, gas virially heats and its temperature can exceed the traditional virial temperature within $r_{200}$/2 as seen in our adiabatic simulations. The consequences of this additional heating is substantial because halos can collapse and form stars before the virial temperature reaches the critical temperature, such as $\sim$10,000K for Ly$\alpha$ cooling.

The timescale at which the center collapse occurs is crucial to the type of object that forms there. If the collapse occurs faster than a Kelvin-Helmholtz time for a massive star ($\sim$ 300 kyr), a black hole might form from the lowest angular momentum gas. Conversely if the collapse is delayed by turbulent pressure, star formation could occur in the density enhancements created by turbulent shocks. The ensuing radiative feedback may create outflows and thus slow further infall and possibly prevent the formation of a central black hole. The nature of the first galaxies poses an important question in the high-redshift structure formation, and to address this problem we must consider their progenitors – the first stars.

5.1. Pop III Feedback

Numerical simulations have shown that the first, metal-free (Pop III) stars form in isolation in its host halo. They are believed to have stellar masses $\sim$100 $M_{\odot}$ (Abel et al., 2002; Omukai & Palla, 2003; Tan & McKee, 2004; Yoshida et al., 2006b) and produce \sim$$10^{50} photons s*-1* that can ionize hydrogen and dissociate H2 (Schaerer, 2002). One-dimensional radiative hydrodynamical calculations (Whalen et al., 2004; Kitayama et al., 2004) and recently three-dimensional radiative hydrodynamical AMR (Abel et al., 2007) and SPH (Yoshida et al., 2006a) simulations found that pressure forces from the radiatively heated gas drive a $\sim$30 km s*-1* shock outwards and expels the majority of the gas in the host halo. Additionally the star ionizes the surrounding few kpc of the intergalactic medium (IGM).

Pop III stellar feedback invalidates some of our assumptions in the calculations presented here, but the general aspect of kinetic energy being dominant should hold in the presence of these feedback processes. In a later paper, we will expand our simulations to include radiative feedback from primordial stars (cf. Yoshida et al., 2006a; Abel et al., 2007) and the metal enrichment from pair instability supernovae (Barkat et al., 1967; Bond et al., 1984; Heger & Woosley, 2002) of the IGM and subsequent star formation.

6. SUMMARY

We have investigated the process of virialization in pre-galactic gas clouds in two cosmology AMR realizations. Our virial analyses included the kinetic (turbulent) energies and surface pressures of the baryons in the system. The significance of each energy component of the gas varies with the effectiveness of the radiative cooling, which we quantify by performing each realization with adiabatic, hydrogen and helium, and H2 cooling models. We highlight the following main results of this study as:

1. Inside $r_{\rm{cool}}$, gas cannot virialize alone through heating but must gain kinetic energy. It is up to a factor of five greater than thermal energy throughout the protogalactic halos. This manifests itself in a faster bulk inflow and supersonic turbulent motions.

2. In the radiative cooling models, supersonic turbulence ($\mathcal{M}$ = 1–3) leads to additional cooling within turbulent shocks. We expect turbulence in larger galaxies, up to 10${}^{12}{M_{\odot}}$, to be even more supersonic.

3. Baryonic velocity distributions are Maxwellian that shows violent relaxation occurs for gas as well as dark matter. Turbulent velocities exceed typical rotational speeds, and these halos are only poorly modeled as solid body rotators.

4. Virial shocks between the void-halo interface occur between $r_{200}$/2 and $r_{200}$. Dense, cold flows in filaments do not shock-heat until well within $r_{200}$ and as small as $r_{200}$/4.

5. Turbulence generated during virialization mixes angular momentum so that it redistributes to a radially increasing function (the Rayleigh criterion).

After the halo virializes, its central part will undergo turbulent collapse, such as in primordial star formation and galactic molecular clouds. These collapses should be ubiquitous in early structure formation as turbulence can be generated through virialization, merging, and angular momentum segregation. We conclude that turbulence plays a key role in virialization and galaxy formation.

This work was supported by NSF CAREER award AST-0239709 from the National Science Foundation. We appreciate Marcelo Alvarez and an anonymous referee providing very constructive feedback on this paper. We thank Michael Norman, Ralph Pudritz, Darren Reed, Tom Theuns, and Peng Wang for helpful discussions. We also thank Paul Kunz for his invaluable help with our analysis software. We are grateful for the continuous support from the computational team at SLAC. We performed these calculations on 16 processors of a SGI Altix 3700 Bx2 at KIPAC at Stanford University.