FORMATION OF COMPACT STELLAR CLUSTERS BY HIGH-REDSHIFT GALAXY OUTFLOWS. I. NON-EQUILIBRIUM COOLANT FORMATION

As the minihalos we are considering in this paper are made up of primordial gas, their chemical makeup is highly restricted, with contributions from only hydrogen, helium, and low levels of deuterium. Yet even these three isotopes can exist in a variety of ionization states and molecules and are thus associated with a substantial network of chemical reactions that must be tracked throughout our simulations.

2.1.1. Implementation

The chemical network that was implemented into FLASH is outlined by Glover & Abel (2008, hereafter GA08). Throughout our simulations we track three states of atomic hydrogen (H, H⁺, and H⁻) and atomic deuterium (D, D⁺, and D⁻), three states of atomic helium (He, He⁺, and He⁺⁺), two states of molecular hydrogen (H₂ and H⁺₂) and molecular hydrogen deuteride (HD and HD⁺), and electrons (e⁻). For simplification, any reaction that involved molecular deuterium (D₂) and all three-body reactions were neglected. As stated in GA08, the very small amount of D₂ and D⁺₂ produced makes any cooling by these molecules irrelevant, while three-body reactions only become important at n ≳ 10⁸ cm⁻³ (e.g., Palla et al. 1983), many orders of magnitude denser than the conditions considered here. With these constraints, a total of 84 reactions were used out of the 115 described in GA08.

Photodissociation rates due to an external radiation field were also included as given in Glover & Savin (2009). These rates are calculated assuming a T_eff = 10⁵ K blackbody source and their strength is quantified by the flux at the Lyman limit, J(ν_α) = 10⁻²¹J₂₁ erg s⁻¹ cm⁻²Hz⁻¹ sr⁻¹. Note that once H₂ and HD are produced in sufficient quantities, some molecules are self-shielded from the background radiation. However, for simplicity, we considered only the case where there was no self-shielding, and thus our results place an upper limit on the effect of a dissociating background. This process adds an additional 7 reactions for a total of 91 reactions in the chemical network.

In reactions that involve free electrons recombining with ions, there are two possible choices for the reaction rate, depending on the overall optical depth of the cloud to ionizing radiation. In the optically thin case (case A; Osterbrock 1989) ionizing photons emitted during recombination are lost to the system, while in the optically thick case (case B), ionizing photons are reabsorbed by neighboring neutral atoms, which have the effect of lowering the recombination rates by essentially not allowing recombination to the ground state. There are three reactions (H⁺ + e⁻ → H + γ, He⁺ + e⁻ → He + γ, and He⁺⁺ + e⁻ → He⁺ + γ) where this is a concern, and as we shall see below, in all cases our clouds are optically thick, such that case B recombination rates are appropriate.

The binding energy from each species is also important in the total energy budget and on the evolution of the gas. In FLASH these are defined such that all of the neutral species (H, D, e⁻, and H) have binding energies equal to zero. As the gas is heated and the atomic species begin to ionize, the endothermic reactions remove the binding energy between the nucleus and electron(s) from the internal energy of the gas. For H and D this requires 13.6 eV, while the ionized states for helium (He⁺ and He⁺⁺) have ionization potentials of 24.5 eV and 79.0 eV, respectively. H⁻ and D⁻ are only weakly bound and have similar binding energies of 0.75 eV. Finally, H₂ and HD have binding energies of 4.4 eV, and H⁺₂ and HD⁺ have binding energies of 10.9 eV, somewhat lower than the atomic species.

To describe the evolution of our 14 species, we enumerate them with an index i such that each has Z_i protons and A_i nucleons, following the structure and syntax from Timmes (1999). Next, we consider a gas with a total mass density ρ and temperature T, and denote the number and mass densities of the ith isotope as n_i and ρ_i, respectively. For each species we also define a mass fraction

$$\begin{equation} X_i \equiv \rho _i/ \rho = n_i A_i/(\rho N_A), \end{equation} \tag{ 1 }$$

where N_A is Avogadro's number, and we define the molar abundance of the ith species as

$$\begin{equation} Y_i \equiv X_i/A_i = n_i/(\rho N_A), \end{equation} \tag{ 2 }$$

where conservation of mass is given by ∑^N_iX_i = 1. Each of the 14 species can then be cast as a continuity equation in the form

$$\begin{equation} \dot{Y}_i \equiv \frac{dY_i}{dt} = \dot{R_i}, \end{equation} \tag{ 3 }$$

where $\dot{R_i}$ is the total reaction rate due to all the binary reactions of the form i + j → k + l, defined as

$$\begin{equation} \dot{R_i} \equiv \sum _{j,k} Y_l Y_k \lambda _{{\rm kj}}(l) - Y_i Y_j \lambda _{{\rm jk}}(i), \end{equation} \tag{ 4 }$$

where λ_kj and λ_jk are the creation and destruction chemical reaction rates for a given species. If the species in question is affected by UV background radiation, the continuity equation takes the following form:

$$\begin{equation} \dot{R_i} \equiv \sum _{j,k} Y_l Y_k \lambda _{{\rm kj}}(l) - Y_i Y_j \lambda _{{\rm jk}}(i) - Y_i J(\nu _{\alpha }), \end{equation} \tag{ 5 }$$

where the last term accounts for the amount of these species that are destroyed by the background radiation, J(ν_α). Throughout our simulations, changes in the number of free electrons are not calculated directly, but rather at the end of each cycle we use charge conservation to calculate their molar fraction, as

$$\begin{equation} Y_{\rm elec} = Y_{{\rm H}^+} + Y_{{\rm D}^+} + Y_{{\rm HD}^+} + Y_{{\rm H}_2^+} + Y_{{\rm He}^+} + 2 Y_{{\rm He}^{++}} - Y_{{\rm H}^-} - Y_{{\rm D}^-}. \end{equation} \tag{ 6 }$$

Because of the often complex ways that the chemical reaction rates depend on temperature and the intrinsic order of magnitude spread in the rates, the resulting equations are "stiff," meaning that the ratio of the minimum and maximum eigenvalue of the Jacobian matrix, $J_{i,j} = \partial \dot{Y}_i / \partial Y_j$, is large and imaginary. This means that implicit or semi-implicit methods are necessary to efficiently follow their evolution. To address this problem, we arrange the molar fractions of the 13 species, excluding e⁻, into a vector Y, and solve the resulting system of equations using fourth-order accurate Kaps–Rentrop, or Rosenbrock method (Kaps & Rentrop 1979). In this method, the network is advanced over a time step h via

$$\begin{equation} {\bf Y}^{n+1} = {\bf Y}^{n} + \sum ^{4}_{i=1} b_i \bDelta _i, \end{equation} \tag{ 7 }$$

where the $\bDelta _i$ vectors are found by successively solving the four matrix equations

$$\begin{equation} (\hat{1}/\gamma h - \bar{J}) \cdot \bDelta _1 = f({\bf Y}^n), \\ \end{equation} \tag{ 8 }$$

$$\begin{equation} (\hat{1}/\gamma h - \bar{J}) \cdot \bDelta _2 = f({\bf Y}^n+a_{21}\bDelta _1) + c_{21}\bDelta _1/h, \\ \end{equation} \tag{ 9 }$$

$$\begin{equation} (\hat{1}/\gamma h - \bar{J}) \cdot \bDelta _3 = f({\bf Y}^n+a_{31}\bDelta _1 + a_{32}\bDelta _2) + (c_{31}\bDelta _1+c_{32}\bDelta _2)/h,\\ \end{equation} \tag{ 10 }$$

and

$$\begin{eqnarray} (\hat{1}/\gamma h - \bar{J}) \cdot \bDelta _4 &=& f({\bf Y}^n+a_{41}\bDelta _1 + a_{42}\bDelta _2+a_{43}\bDelta _3)\nonumber \\ &&+ (c_{41}\bDelta _1+c_{42}\bDelta _2+c_{43}\bDelta _3)/h. \end{eqnarray} \tag{ 11 }$$

Here, b_i, γ, a_ij, and c_ij are fixed constants of the method, $f({\bf Y}) \equiv \dot{\bf Y},$ $\hat{1}$ is the identity matrix, and $\bar{J}$ is the Jacobian matrix. Note that the four matrix equations represent a staged set of linear equations and that the four right-hand sides are not known in advance. At each step, an error estimate is given for the difference between the third- and fourth-order solutions. For comparison we also carried out tests, using a multi-order Bader–Deuflhard method (Bader & Deuflhard 1983). However, in the end, the Rosenbrock method was chosen over this method because of its efficiency and speed.

As the species evolve, the temperature of the gas changes from the release of internal energy from recombinations or the loss of internal energy from ionizations and dissociations. These changes can in turn affect the reaction rates. Thus, to ensure the stability of the chemistry routine while at the same time allowing the simulation to proceed at the hydrodynamic time step, we developed a method of cycling over multiple Kaps–Rentrop time steps within a single hydrodynamic time step. Here, we estimated an initial chemical time step of each species as

$$\begin{equation} \tau _{{\rm chem},i} = \alpha _{\rm chem} \frac{{Y}_i+0.1Y_{H^+}}{\dot{Y_i}}, \end{equation} \tag{ 12 }$$

where α_chem is a constant determined at runtime that controls the desired fractional change of the fastest evolving species. The changes in molar abundances, $\dot{Y}_i$, were calculated from the ordinary differential equations that make up the chemical network, and the molar fractions of each species Y_i are given by the current values. In both the tests and simulations, we chose a value of α_chem = 0.5.

Note that we offset the subcycling time step by adding a small fraction of the ionized hydrogen abundance to Equation (12). This is because there are conditions where a species is very low in abundance but changing very quickly, for example, rapid ionization of atomic species, which will cause the subcycling to run away with extremely small time steps. In regions in which most species are neutral, this has little effect since the chemical time step is likely longer than the hydrodynamical one, and in regions in which abundances are rapidly changing, then this extra term buffers against very small times steps. It also prevents rapid changes in internal energy as energy is removed as atomic species are ionized and gained as they recombine.

Once calculated, these species time steps are compared to each other and the smallest time step, associated with the fastest evolving species, is chosen as the subcycle time step. If this is longer than the hydrodynamic time step, the hydrodynamic time step is used instead and no additional subcycling is done. If subcycling is required, the species time step is subtracted from the total hydrodynamic time step and the network is then updated over the chemical time step. The species time steps are recalculated after each subcycle and compared to the remaining hydrodynamic time step. This is repeated until a full hydrodynamic time step is completed, as schematically shown in Figure 1.

Zoom In Zoom Out Reset image size

In cases in which the gas is extremely hot or cold, the chemical make-up can be determined directly from the temperature, avoiding the need for matrix inversions. If the temperature is above 10⁵ K then all atomic species become ionized and all molecular species are dissociated, and the network can be bypassed. If the temperature is between 2.0 × 10⁴ and 1.0 × 10⁵ K, then we "prime" the solutions and ionize 5% of available neutral hydrogen, 5% of neutral helium (4.5% into singly ionized helium and 0.5% into doubly ionized helium), before entering the iterative solver, to help accelerate the routine toward the correct solution. Finally, if the temperature is less than 50 K, then all species are kept the same, and no reactions are calculated. This is done because cooling and chemistry rates become unimportant at such low temperatures. It also has the benefit of speeding up the simulation slightly as very little time is spent in either the cooling or chemistry routines. In all other cases, the full network is evolved without alterations.

2.1.2. Chemistry Tests

To test our implemented chemical network, we carried out a series of runs in which initially dissociated and ionized gas was held at constant temperature and density for 10¹⁶ s. A small initial time step (t₀≈ 10⁶ s) was used and allowed to increase up a maximum time step of 10¹² s. Models were run with total hydrogen number densities varying from 0.01 cm⁻³ to 100 cm⁻³, and temperatures ranging between 10² K and 10⁴ K, and no external radiation. In each case, the results were compared to the results of a different implementation of the same chemical network within the ENZO code (S. C. O. Glover 2009, private communication, G09 hereafter), yielding the molar fractions shown in Figure 2. In this figure, the three columns correspond to runs with different temperatures, the curves correspond to runs with different densities, and the rows correspond to the evolution of different species.

Zoom In Zoom Out Reset image size

The match between our tests and the numerical results from G09 is excellent. In all cases and at all temperatures, the curves closely track each other, in most cases leading to curves that are indistinguishable. Although the abundances of several species change by many orders of magnitude throughout the runs, the two methods track each other within 10% in all cases except for H₂ at 10⁴ K, which is unimportant as a coolant at this temperature but nevertheless consistent within a factor of 1.5 at all times. Furthermore, this agreement between methods is also seen for deuterium species, which are not shown in this figure as they follow H exactly, maintaining a $\frac{1}{6000}$ ratio between both species at all times.

At T = 100 K, all ionized species quickly recombine with the free electrons to form neutral atoms. However, even during this relatively quick transition from ionized to neutral, H⁺ and H⁻ ions (not shown) persist for long enough to catalyze the formation of substantial amounts of molecular gas, leading to final H₂ molar fractions of ≈10⁻⁴. At T = 1000 K, the evolution is very similar to the T = 100 K case, although the species do not reach equilibrium as quickly, leading to even higher levels of H₂ formation. Finally, at T = 10⁴ K, it takes even longer for the ionized hydrogen to recombine, but in this case, less molecular species are formed, as collisional dissociation of H₂ and HD are more prevalent, limiting the maximum amount of these species.

Also apparent in these plot is the dependence of the species evolution on the density of the gas. Chemical reactions are fundamentally collisional processes whose rates are quadratic in number density. Thus, as we are not considering three-body interactions, the timescale associated with chemistry should decrease linearly with the density. This is seen for all temperatures and species shown in Figure 2, as in every case each line is separated from its neighbor by a factor of 10 in time, exactly corresponding to the density shift between cases.

2.1.3. Effect of the Background Radiation

Background radiation with photon energies between 11.2 and 13.6 eV can excite and dissociate molecular hydrogen. In the absence of other coolants, this can have drastic effect on the evolution of the cloud. Two extremes are immediately apparent, a strong background case in which any H₂ or HD formed is quickly dissociated, and a background-free case in which no molecules are photodissociated. A simple test was constructed to study the effect of the background and determine a fiducial value for J₂₁. J₂₁ was varied between 0 and 1 at five different values. For each value of J₂₁ the number density was varied between n = 10⁻¹ and 1.0 cm⁻³. Each test was run at a constant temperature and constant density with evolving chemistry and no cooling. The results are given in Figure 3. From this, we determine that only background levels at or above J₂₁ = 0.1 give an appreciable difference in the abundance of H₂ and HD over a megayear timescale, which as we shall see below, is the timescale of shock–minihalo interactions. At the same time, J₂₁ = 0.1 provides a reasonable upper limit to the level of background expected before reionization (e.g., Ciardi & Ferrara 2005). Therefore, we use this value as a fiducial value in the simulations with a background.

Zoom In Zoom Out Reset image size

The second major process added to the code was radiative cooling, which was divided into two temperature regimes. At temperatures ⩾10⁴ K, cooling results mostly from atomic lines of H and He, with bremsstrahlung radiation also becoming important at temperatures above 10⁷ K. Below 10⁴ K, on the other hand, the net cooling rate is determined by molecular line cooling from H₂ and HD, which, as it is an asymmetric molecule, can radiate much more efficiently than H₂, and thus can be almost as important although it is much less abundant. Cooling from H₂ operates down to T ⩽ 200 K and to number densities n>10⁴ cm⁻³ (Glover & Abel 2008; Galli & Palla 1998), while HD which can cool the gas to slightly lower temperatures and to higher number densities (Bromm et al. 2002). As we are restricting ourselves to primordial gas in this study at any given temperature the overall cooling rate, Λ_Total, is the combination from both regimes:

$$\begin{equation} \Lambda _{\rm Total} = \Lambda _{\rm Atomic} + \Lambda _{\rm Molecular}. \end{equation} \tag{ 13 }$$

Each cooling rate has the form

$$\begin{equation} \Lambda _{i,j} = n_{i} n_{j} \lambda _{i,j}, \end{equation} \tag{ 14 }$$

where Λ_i,j is the energy loss per volume due to species i and j, n_i and n_j are the number densities of each species, and λ_i,j is the cooling rate in erg cm^-3 s⁻¹. Cooling rates for the collisional excitation between H₂ and H, H₂, H⁺, and e⁻ and between H⁺₂ and H or e⁻ are taken from GA08. The cooling rate for the collisional excitation between HD and H is taken from Lipovka et al. (2005). Finally, cooling rates from hydrogen and helium atomic lines are calculated using CLOUDY (Ferland et al. 1998). In calculating these rates, we followed the procedure described in Smith et al. (2008) and used the "coronal equilibrium" command which considers only collisional ionization. The cooling curve was calculated assuming case B recombination for the recombination lines of hydrogen and helium, as discussed further in Section 3.1.

Any cooling routine contains a natural timescale that relates the total internal energy to the energy loss per time:

$$\begin{equation} \tau _{\rm cool} = \frac{\alpha _{\rm cool} \times E_i}{\dot{s}}, \end{equation} \tag{ 15 }$$

where α_cool is a constant between 0 and 1, in all cases set at 0.1, E_i is the internal energy, and $\dot{s}$ is the energy loss per time. Cooling rates are very dependent on temperature and species abundances, and these quantities can change rapidly over a single chemical time step.

A method of subcycling over cooling time steps was developed to ensure that the correct cooling rates are used. An initial cooling timescale is calculated assuming α_cool = 0.1 using Equation (15) which is then compared to the chemical time step. If τ_cool is smaller than the fraction of the chemistry time step then that fraction of energy is subtracted from the internal energy and temperature. The cooling rate and cooling time step is recalculated with the updated temperature. This continues until the chemistry time step is reached. This is schematically given in Figure 4.

Zoom In Zoom Out Reset image size

2.2.1. Cooling Tests

As a test of our cooling routines, we reproduced the example curves given in Prieto et al. (2009). In this work, the authors present the effects of H₂ and HD cooling in a primordial gas. The gas begins at an initial temperature of T = 500 K with initial number densities, relative to hydrogen: $n_{{\rm H}^+} = 10^{-4}; n_{{\rm H}^-} = n_{{\rm H}_2^+} = 10^{-12}; n_{{\rm H}_2} = 10^{-3}; n_{\rm D} = 10^{-5}; n_{{\rm D}^+} = 10^{-9}; n_{{\rm HD}} = 10^{-6}; n_{{\rm HD}^+} = 10^{-18}; n_{{\rm He}^+} = n_{{\rm He}^{++}} = 0.0$1 and with initial hydrogen and helium densities of ρ_H = 0.75 × ρ_tot and ρ_He = 0.24 × ρ_tot, where ρ_tot is the total baryonic matter density.

Three models were run with total number densities of n_tot = 1.0, 10.0, and 100.0 cm⁻³. Cooling was tracked for 10⁸ yr with chemistry evolving simultaneously. The results of this calculation are shown in Figure 5, which indicates good agreement with Prieto et al. (2009). It should be noted that temperature evolution in this plot has a linear dependence of the number density of the gas. For example, a gas with 10 times the number densities of another gas will cool 10 times quicker. This is again because most of the cooling is coming from the collisions between two species, in this case H₂ or HD and H.

Zoom In Zoom Out Reset image size

As mentioned above, HD can be more important than H₂ for gas cooling at higher densities and colder temperatures. To determine whether or not HD cooling is important in this simulation, we apply the cooling test to two different scenarios. First, we use the same initial abundances as described above and second, using primordial abundances with a small fraction (0.01%) of each atomic species ionized. Each test was run twice, once with deuterium and once without. The results of these tests are given in Figure 6. At high number densities, HD cooling does not have a perceivable effect. At intermediate temperatures, HD cooling is important for a gas with initial abundances from the cooling tests. Finally, at low temperatures, HD cooling is very important in both cases.

Zoom In Zoom Out Reset image size

A simple model is used for the gas and dark matter of the protocluster whose collapse redshift of z_c = 10 (a cosmic age of ≈0.5 Gyr) is taken to be just before the epoch of reionization, and whose total mass of M_c = 3.0 × 10⁶ M_☉ is taken to be on the large end of minihalos formed at this redshift. The gas is assumed to have a primordial composition of 76% neutral atomic hydrogen and 24% neutral atomic helium by mass. Initially, the cluster has a mean density that is enhanced by a factor Δ = 178 (e.g., Eke et al. 1998) above the background, ρ_c = ΔΩ₀(1 + z_c)³ρ_crit = 6.54 × 10⁻²⁵g cm^-3. In this case, the cloud's virial radius is R_c = 0.393 kpc and virial velocity is v_c = 6.55 km s^-1. We assume that the radial profile is given by Navarro et al. (1997):

$$\begin{equation} \rho (R) = \frac{\Omega _0\rho _c}{cx(1+cx)^2} \frac{c^2}{3F(c)} \ \rm g\,cm^{-3}, \end{equation} \tag{ 16 }$$

where c is the halo concentration factor, x = R/R_c, and $F(t) \equiv \ln (1+t)-\frac{t}{1+t}$. We assume that as the gas collapses inside the dark matter halo, it is shock-heated to its virial temperature, T_c = 1650 K, and develops a density distribution of isothermal matter in the CDM potential well:

$$\begin{equation} \rho _{\rm gas}(R) = \rho _0e^{ - (\frac{ v_{\rm esc}^2(0)-v_{\rm esc}^2(R) }{ v_c^2 }) } \rm {\rm g}\,cm^{-3}, \end{equation} \tag{ 17 }$$

where the escape velocity as a function of radius is given by v(xR_vir) = 2v²_c[F(cx) + cx(1 + cx)⁻¹][xF(c)]⁻¹. From Madau et al. (2001), we take a typical value of the halo concentration to be c = 4.8, although some observations suggest that high-redshift halos may be less concentrated than expected from this estimate (Bullock et al. 2001). With this value of c, we can compute the central density as

$$\begin{eqnarray} \rho _0 &=& \frac{(178/3)c^3\Omega _b e^A\left(1+z\right)^3}{ \int _0^c \left(1+t\right)^{A/t} t^2 dt} \rho _{\rm crit} = 39215 \ \Omega _b \rho _{\rm crit} \left(1+z_c\right)^3 \nonumber\\ &=& 2.16\times 10^{-23}\, {\rm g \;cm^{-3}}, \end{eqnarray} \tag{ 18 }$$

where A ≡ 2c/F(c) = 10.3 and t = cx.

To determine which case to use in our chemistry routine the optical depth for H⁺, He⁺, and He⁺⁺ recombination was calculated from this profile:

$$\begin{equation} \tau _{\nu }(r) = \int ^{r}_{r_0} \sigma _{\nu } n(r^{\prime }) dr^{\prime }, \end{equation} \tag{ 19 }$$

where σ_ν is the cross section of interaction and n(r') is the number density. For hydrogen-like atoms, the cross section is

$$\begin{equation} \sigma _{\nu } = \frac{7.91 \times 10^{-18}}{Z^2} \left(\frac{\nu _1}{\nu }\right)^3 \ {\rm g} \rm \ \mathop {\rm cm}\nolimits ^2, \end{equation} \tag{ 20 }$$

where hν₁ = 13.6 Z² eV and Z is the nuclear charge. Calculating the optical depth from Equations (17)–(20) yields the results given in Table 1. We assume for the case of He⁺⁺ recombination that the surrounding helium is singly ionized. In all cases the optical depth is much greater than 1 (see Section 3) and therefore we use case B rates for all recombination reactions in our fiducial simulations.

As the minihalo in our simulation is made up of both gas and collisionless dark matter, a two-part gravity scheme was required. First, the Ricker (2008) multigrid Poisson solver was used to calculate the gravitational potential due to the gas component. Second, the acceleration due to the total matter was calculated using Equation (17) above. The general equation for the gravitational acceleration of an ideal gas is given by $a_{\rm Grav} = \frac{k}{\rho (R) m_p} \big[ T(R) \frac{\partial \rho (R)}{\partial R} + \rho (R) \frac{\partial T(R)}{\partial R} \big] \ \rm cm\,s^{-2}.$ To account for the dark matter halo's contribution to gravity, we calculated the gravitational acceleration of the total matter and subtracted its contribution from the baryonic matter in the initial configuration using the above equation with constant temperature. Finally, we added the gravitational contribution from the self-gravity of the gas. The total acceleration is simply given by

$$\begin{equation} a_{\rm Tot} = a_{\rm M,0} - a_{\rm gas,0} + a_{{\rm SG}}, \end{equation} \tag{ 21 }$$

where a_M is the acceleration from the total initial mass density as given by Equation (17), a_gas,0 is the contribution from the baryonic matter in the initial configuration, and a_SG is the contribution from the self-gravity calculated from the Poisson solver. Initially, when the minihalo gas is in hydrostatic balance with its surroundings, these last two terms will cancel each other and the cloud will remain unchanged. When the cloud is disrupted and cooling takes effect, the self-gravity will cause the cloud to collapse.

The above equations are correct up to the viral radius of the cloud. To ensure a smooth density transition from the cloud, we simply keep the gas outside of the cloud be gravitationally bound to the cloud and solve for the expected density profile. The acceleration here is then

$$\begin{equation} a_{\rm Grav} =-\frac{G\!M(R > R_c)}{R^2}, \end{equation} \tag{ 22 }$$

and the density is

$$\begin{equation} \rho (R > R_c) = \rho (R_c)e^{({\frac{R_0}{R} - \frac{R_0}{R_c}})}, \end{equation} \tag{ 23 }$$

with R₀ = GM_cm_p/k_bT. Here, G is the gravitational constant, m_p is the mass of a proton, k_b is Boltzmann's constant, and, as above, M_c and T_c are the mass and temperature of the cloud, respectively. As R → ∞, the density goes down to a small fraction of ρ(R_c). A test of our gravity routine showed that the cloud was able to maintain hydrostatic balance for many dynamical times in the absence of an impinging galaxy outflow.

A Sedov–Taylor solution is used to estimate the properties of the galactic outflow. The initial input energy is taken to be E = E₅₅(erg), where E₅₅ is the energy of the SNe driving the wind in units of 10⁵⁵ erg, and the wind efficiency is derived from the amount of kinetic energy from the SNe that is channeled into the outflow. The shock expands into a gas that is δ times greater than the background at a redshift of z_c. As in Scannapieco et al. (2004), we assume that the cloud is a distance R_s = 3.6 kpc, using fiducial values: z_c = 10, E₅₅ = 10, M₆ = 3, δ = 44, and = 0.3.

With these values, the velocity of the blast front is v_s = 415 km s⁻¹, when it reaches the minihalo, and the resulting temperature of the fully ionized post-shock medium is T = 2.4 × 10⁶ K. By the time the shock has covered the separation distance, R_s, it will have entrained a mass

$$\begin{equation} M_{s, \rm Total} = 4.4 \times 10^{7}\, M_{\odot }, \end{equation} \tag{ 24 }$$

with a surface density of

$$\begin{equation} \sigma _{s} = 2.6 \times 10^{5}\, M_{\odot }\,\rm kpc^{-2}. \end{equation} \tag{ 25 }$$

Note that while the above equations are the solutions that come from a simple spherical blast wave, the wind in the simulation is still well approximated by a plane wave solution, because the size of cloud is much smaller than the distance between the SN and the cloud.

To model this wind, a time-dependent boundary condition is imposed at the leftmost boundary of the simulation volume. The expected lifetime of the shock is given by σ_s = v_post ρ_post t_s, where v_post is the post-shock velocity of the blast wave, ρ_post is the post-shock density, and σ_s is the surface density of the entrained material. Solving for t_s and putting in the appropriate values, the expected shock life time is t_s = 2.5 Myr. After this time the shock begins to taper off with the density decreasing and temperature increasing and keeping the pressure constant. This is done to prevent the excessive refinement that a sharp cutoff would cause. The density falls off as

$$\begin{equation} \frac{\rho (t)}{\rho (0)} = 0.01 + 0.99 e^{-\tau _s/1.5} \end{equation} \tag{ 26 }$$

and the temperature rises as

$$\begin{equation} \frac{T(0)}{T(t)} = 0.01 + 0.99 e^{-\tau _s/1.5} \, \end{equation} \tag{ 27 }$$

where τ_s is defined as $\frac{t - 1 \,{\rm Myr}}{1 \,{\rm Myr}}$. This also prevents the hydrodynamic (Courant) time step from becoming extremely short behind the shock, in order to maintain pressure equilibrium in an extremely rarified medium.

Note that in this initial study, the dark matter and gas distributions have been somewhat idealized, and more complicated geometries could be used to model these components in greater detail. For example, a triaxial instead of a spherical distribution could be assumed for the dark matter halo, inhomogeneities could be added to the minihalo gas, and the shock could be assumed to impact the minihalo off-axis. While each of these possibilities would be qualitatively interesting and would naturally alter the final outcome of the halo, they are nevertheless beyond the scope of this study.

In this simulation, several distinct stages of evolution are identified during the interaction between the cloud and the outflow, as shown in Figures 7 and 8. Initially, the cloud is in hydrostatic equilibrium as the shock enters the simulation domain. If it were not supported by pressure, the cloud would collapse on the free-fall time which, using the average cloud density, is

$$\begin{equation} t_{\rm ff} = \sqrt{\frac{3 \pi }{32 G \rho }} \approx 100 \rm \ Myr. \end{equation} \tag{ 28 }$$

As the cloud is initially in hydrostatic balance, the initial sound crossing time is similar to the free-fall time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the shock contacts and surrounds the cloud, it heats and begins to ionize the gas. The shock completely envelops the cloud on a characteristic "intercloud" crossing timescale, defined by Klein et al. (1994) as

$$\begin{equation} t_{\rm ic} = \frac{2 R_c}{v_{s}} \approx 4.5 \rm \ Myr. \end{equation} \tag{ 29 }$$

As the cloud is enveloped, the shock moves through the outer regions fastest and ionizes this gas first. This in turn promotes rapid molecule formation as the gas cools and recombines incompletely, leaving H⁺ and H⁻ to catalyze the formation of H₂ and HD. Interestingly, because the shock slows down as it moves through denser material, the gas behind the center of the halo remains undisturbed until the enveloping shocks meet along the axis at the back of the halo. The leads to a "hollow" H₂ distribution at 6.6 Myr, in which the molecular coolants are confined to a shell surrounding the undisturbed, purely atomic gas.

After the enveloping shocks collide at the back of the cloud, a strong reflected shock is formed that moves away from the rear of the cloud and back through the halo material. Without cooling, this reflected shock would eventually lead to cloud disruption (Klein et al. 1994). However in our case, the shock has the opposite effect. It moves through the cloud, and the gas is briefly ionized, but then quickly cools and recombines, forming H₂ and HD throughout the cloud. This can be seen in the upper row of Figure 8, which shows the conditions at ≈8 Myr.

At this point the cloud is denser, smaller, and full of new coolants. Using the conditions from the center of the cloud 8 Myr after the start of the simulation, we calculate new timescales. Now the free-fall time is 21 Myr and the sound crossing time is ≈27 Myr. The cloud is cold and dense enough to start collapsing.

The timescale for the formation of H₂, given in GA08, is

$$\begin{equation} t_{\rm H_2} = \frac{X_{H_2}}{k_1X_{e}n} \rm (s), \end{equation} \tag{ 30 }$$

where $X_{{\rm H}_2}$ is the mass fraction of H₂, X_e is the mass fraction of electrons, k₁ is the reaction rate for the formation of H⁻ (H + e⁻ → H⁻ + γ), and n is the total number density (≈1 cm⁻³ at 8 Myr). Initially, as the shock begins to impact the cloud, this timescale is very short, on the order of 0.1 Myr to get a final abundance of ≈10⁻⁵. As the abundance of H₂ increases and the abundance of electrons decrease, this timescale quickly increases. Although as the cloud collapses the density increases which lowers this timescale.

The H₂ cooling timescale, given by Klein et al. (1994), is

$$\begin{equation} t_{\rm cool} = \frac{1.5 n k T}{n_{{\rm H_2}} n_{{\rm H}} \Lambda _{{\rm H,H_2}}} \rm \ (s), \end{equation} \tag{ 31 }$$

where $n_{{\rm H}_2}$ and n_H are the number densities of H₂ and H, respectively, and $\Lambda _{{\rm H,H_2}}$ is the cooling rate between H and H₂. At 8 Myr the H₂ cooling time in most of the cloud is only 0.2 Myr, meaning that pressure support drops dramatically after this time. Any expansion due to shock heating is halted as the gas is quickly cooled by H₂ and HD as they form. Furthermore, as the cloud collapses, the chemistry and cooling timescales decrease, rapidly accelerating the collapse.

The final state of the cloud in our simulation is a thin cylinder stretching from the center of the dark matter halo to several times the initial virial radius. The temperature of this gas is 100–200 K, much colder than the initial virial temperature. The gas is also much denser than the initial minihalo, reaching values of up to 10⁻²¹  cm⁻³ or n ≈ 10³ cm⁻³, in the center of the cloud, and even this density is probably only a lower limit set by the resolution of our simulation. On the other hand, the cloud is quite extended along the x-axis, with substantial differences in velocity along the cylinder. Thus, it is continually stretched and fragments until the end of the simulation at 14.7 Myr (row 3 in Figure 8).

Figure 9 shows rendered density contours of the major stages of evolution of the cloud from t = 0 through the end of the simulation. The first panel shows the initial configuration, with the cloud in hydrostatic equilibrium, and the shock front entering from the left side of the simulation volume. The next panel shows the cloud after being impacted by the shock, highlighting the density enhancement in the outer shell of the minihalo gas. The bottom left panel shows the cloud as it begins to cool and collapse, at a time at which the reverse shock has already passed though the cloud and coolants are found throughout the shocked, recombined material. Finally, the last panel shows the distribution at the end of the simulation. The cloud has now been stretched over a large distance and much of its mass has been accelerated to above the escape velocity, moving outside of the dark matter halo. The dense knots of this material in this figure are tightly gravitationally bound, have number densities approaching 10³ cm⁻³, and are destined to form extremely compact stellar clusters.

Zoom In Zoom Out Reset image size

While the collision happens on order of the shock crossing time of the halo, the final distribution of the clumps evolves on the longer timescale defined by R_vir/v_c ≈ 100 Myr. To study the final state of the stretched and collapsed distribution without continuing the simulation out to such extremely long times, we divided the x-axis into 100 evenly spaced bins between x = 0 kpc and x = 1.4 kpc. We then calculated the mass of each bin by summing up the gas from each cell from the FLASH simulation in a cylinder with a 24 pc radius and length of the bin. Similarly, we calculated the initial velocity of each bin by adding the momentum from each cell within this cylinder and dividing by the total mass in each bin.

We evolved this distribution forward in time using a simple numerical model, which assumed that motions were purely along the x-axis and pressure was negligible at late times. In this case, acceleration could be calculated directly from the gravity between each pair of particles and from the potential of the dark matter halo. Furthermore, if any given particle moved past the particle in front of it, we merged them together, adding their masses and calculating a new velocity from momentum conservation.

Evolving the distribution in this way for an additional 200 Myr past the end of the simulation yielded the results shown in Figure 10. As the stretched cloud continues to move outward, particles begin to attract each other, and eventually merge together to create larger clumps. By 100 Myr most of the particles have merged, after which their motions are purely ballistic. This can be seen in the top panel as the lines for the late times overlap each other and in the middle panel as the velocity profiles overlap.

Zoom In Zoom Out Reset image size

At the final time of 200 Myr after the end of the numerical simulation, three small, stable clumps with masses of 5.0 × 10⁴ M_☉, 4.0 × 10⁴ M_☉, and 3 × 10⁴ M_☉ can be seen in the top panel of Figure 10. Each of these new peaks is located far outside of the original dark matter halo.

At temperatures above 10⁴ K, the primary source of cooling is atomic lines from hydrogen and helium. Although we have shown that for the primordial cloud, case B rates should be used for both the chemical network and cooling functions, Figure 11 shows a comparison between our fiducial run, HBN, and a run in which reaction and cooling rates are taken for case A recombination (HAN). The high temperature H–He cooling curve is taken from Wiersma et al. (2009). The upper panels show density contours and the bottom show contours of H₂ abundance.

Zoom In Zoom Out Reset image size

As expected, the case B simulation produces greater molecular coolant abundance at similar overall densities. This difference can be seen in the lower two panels of the first column of Figure 11 at 6.63 Myr. To remain in pressure support as the cloud becomes denser from the shock, the cloud must get hotter. However, because case A cools slightly faster, this support is quickly removed and the cloud takes on a more extended shape as evident in the first two columns of Figure 11. Although, by 14 Myr the abundance of molecular coolants are very similar between HBN and HAN with each containing $X_{\rm H_2}$ ≈10^−2.5.

In both cases, the fate of the minihalo gas is the same. Atomic cooling occurs sufficiently rapidly to sap the shock of its energy and drop the post-shock temperature to ≈10⁴ K, and non-equilibrium processes step in to provide molecular coolants below 10⁴ K. The gas is then able to collapse and form into a long dense filament within which clumps are formed. In fact, the only substantial differences between the runs are the details of the distributions of clumps, which is somewhat more extended in the case A run as compared to the case B run.

A more uncertain aspect of our simulation is our assumption of a negligible dissociating background. In fact, the presence of at least a low level of dissociating background is necessary in order for the minihalo not to collapse and form stars on its own, cooling by H₂ and HD left over from recombination. To set an upper limit on the impact of such a background we modify the rates in our chemical network to approximate a relatively large dissociating background of J₂₁ = 0.1, as discussed in Section 2.1. Furthermore, as these rates are modified for all reactions throughout the simulation, this background is taken to affect even the densest regions of the cluster. This is equivalent to assuming that the cloud is optically thin to 11.2–13.6 eV photons at all times during the simulation.

Figure 12 shows the comparison between the run including this background (HBY) and the fiducial run HBN. As expected, the abundance of H₂ is reduced in the case with the UV background, peaking at about ≈10⁻⁴ instead of ≈10⁻² in the run without a background. Interestingly, this difference persists even after a few megayears into the simulation, and the abundance of the HPY run remains stable at about ≈10⁻⁴.

Zoom In Zoom Out Reset image size

However, this value is more than sufficient to cool the gas to the same temperature as in HBN. Even with this lower mass fraction, the cooling time is smaller than the dynamical time, and the evolution of the cloud remains essentially unchanged. The cloud collapses and is stretched into the same configuration as found without a background. Dense clouds are again found between 0.2 and 0.4 kpc, at 0.55 kpc, and 0.9 kpc and the density and temperature of each of these clouds is comparable to those found in the fiducial without a background. By neglecting any molecular self-shielding, this represents a worst case scenario for H₂, yet shock–minihalo interactions continue to make compact stellar clusters.

Finally, we considered the impact of the maximum refinement level on our results. Figure 13 compares the fiducial run with six levels of refinement and an effective resolution of 4.55 pc to a lower resolution run (LBN) with five maximum levels of refinement and an effective resolution of 9.1 pc. Here, density is shown in the upper two rows in each pair while the mass fraction of H₂ is shown in the bottom two rows.

Zoom In Zoom Out Reset image size

Only a slight dependence on the formation for H₂ is found with resolution. This is most apparent in the second column, corresponding to t = 7.7 Myr, which shows that the shocks are slightly broadened in the lower resolution case. Since both chemical reactions and cooling go as n², this smearing out has the effect of slightly decreasing H₂ formation and cooling in the lower resolution run. However, enough H₂ is produced in both cases for the cloud to collapse efficiently, and evolve in the same manner up until late times, when the difference in H₂ abundance is small.

Furthermore, we also conducted similar resolution studies using case A recombination, and also modifying chemistry and cooling to account for the presence of a dissociating UV background. Again comparisons between the high-resolution runs (HAN and HBY) with the low-resolution runs (LAN and LBY) uncovered only weak differences with resolution. Compact stellar clusters were formed in all cases.

If indeed shock–minihalo interactions lead to massive clusters of stars, the longest-lived stars in these objects will remain observable down to low redshift, providing a direct observational constraint on our results. Furthermore, the clusters formed in our simulations are extremely compact, and thus unlikely to be tidally disrupted as they eventually become gravitationally bound to the even larger structures that form over cosmological time. As discussed in Scannapieco et al. (2004), the most natural candidate for these old and dense stellar clusters is the population of metal-poor globular clusters, associated with the halos of galaxies (e.g., Zinn 1985; Ashman & Bird 1993; Larsen et al. 2001; Strader et al. 2005; Brodie & Strader 2006).

There are several detailed properties of the clusters in our simulations that support this association. Halo globular clusters, like the more metal-rich (disk) globular clusters, are extremely compact, with typical half-light radii of ≈3 pc (e.g., Jordán et al. 2005). Unlike higher-metallicity globular clusters, however, which may have mostly formed at intermediate redshifts (Elmegreen 2010) the age of all metal-poor globular clusters is between 10 and 13 Gyr, placing their formation within or before reionization, during the epoch in which minihalos had not yet been evaporated by ionization fronts.

While globular clusters exits at a range of masses, their mass distribution is well described by a Gaussian in log₁₀(M_gc) with a 0.5 dispersion, and mean at 10⁵ solar masses, precisely spanning the masses of the clumps formed in our simulations. Although the lower mass limit of globular clusters may be set by destruction through mechanical evaporation (e.g., Spitzer & Thuan 1972) and shocking that occurs as the cluster passes through the host galaxies (e.g., Ostriker et al. 1972), the upper mass limit appears to be an intrinsic property of the initial populations (e.g., Fall & Rees 1985; Peng & Weisheit 1991; Elmegreen 2010). Furthermore, this upper limit of ≈10⁶ M_☉ in stars corresponds roughly with the 10⁴ K virial temperature above which atomic cooling becomes effective in high-redshift virialized clouds of gas and dark matter (Fall & Rees 1985). No ≳10⁶ M_☉ stellar clusters can result from shock–minihalo interactions, because no minihalos exist with gas masses ≳10⁶ M_☉.

The second clue as to the origin of globular clusters comes from observations of stars being tidally stripped from these objects. If, like galaxies, globular clusters are contained within individual dark matter halos, the stripping of stars would be highly suppressed, due to the increased gravitational potential. Surprisingly, no such evidence of dark matter halos is seen (Moore 1995), suggesting that these clusters formed through a mechanism markedly different than star formation within galaxies.

Again, the triggered star formation seen in our shock–minihalo simulations provides a natural explanation of this key observed property. While gas is initially gathered together by a collapsed dark matter potential, the lack of effective coolants prevents this collapse from continuing to the point that stars are formed at the center of the potential. Instead the momentum of the impinging galaxy wind is such that it is able to accelerate the gas above the escape velocity of the halo—even as it induces cooling and compresses the gas to the point that it remains gravitationally bound on its own. The trigger for star formation and the mechanism that removes the dark mater halo are one and the same, resulting in a population of stellar clusters that form at the moment they break free from their dark matter enclosures.