ATOMIC CHEMISTRY IN TURBULENT ASTROPHYSICAL MEDIA. I. EFFECT OF ATOMIC COOLING

Our ionization network tracks the impact of 36 separate reactions acting on 24 species and six elements: hydrogen (H, H⁺ ), helium (He, He⁺ , He²⁺ ), carbon (C-C⁴⁺ ), oxygen (O-O⁷⁺ ), sodium (Na, Na⁺ ), magnesium (Mg, Mg⁺ , Mg²⁺ ), and free electrons (e⁻ ). For each species, we track several reactions, including radiative recombinations, dielectronic recombinations, and ionization due to electron impacts. A summary of the reactions considered as well as their source is given in Table 1.

Table 1.  List of the Reactions Considered and their Source

Number	Reaction	Source
0001	He+ + e− $\to$ He	1
0002	C+ + e− $\to$ C	5
0003	C2+ + e− $\to$ C+	4
0004	C3+ + e− $\to$ C2+	3
0005	C4+ + e− $\to$ C3+	2
0006	O+ + e− $\to$ O	7
0007	O2+ + e− $\to$ O+	6
0008	O3+ + e− $\to$ O2+	5
0009	O4+ + e− $\to$ O3+	4
0010	O5+ + e− $\to$ O4+	3
0011	O6+ + e− $\to$ O5+	2
0012	O7+ + e− $\to$ O6+	1
0013	Mg+ + e− $\to$ Mg	8
0014	Mg2+ + e− $\to$ Mg+	9
0015	Na+ + e− $\to$ Na	8
0016	H+ + e− $\to$ H	12
0017	H + e− $\to$ H+ + 2 e−	11
0018	He + e− $\to$ He+ + 2 e−	11
0019	He+ + e− $\to$ He2+ + 2 e−	11
0020	C + e− $\to$ C+ + 2 e−	11
0021	C+ + e− $\to$ C2+ + 2 e−	11
0022	C2+ + e− $\to$ C3+ + 2 e−	11
0023	C3+ + e− $\to$ C4+ + 2 e−	11
0024	O + e− $\to$ O+ + 2 e−	11
0025	O+ + e− $\to$ O2+ + 2 e−	11
0026	O2+ + e− $\to$ O3+ + 2 e−	11
0027	O3+ + e− $\to$ O4+ + 2 e−	11
0028	O4+ + e− $\to$ O5+ + 2 e−	11
0029	O5+ +e− $\to$ O6+ + 2 e−	11
0030	O6+ + e− $\to$ O7+ + 2 e−	11
0031	Mg + e− $\to$ Mg+ + 2 e−	11
0032	Mg+ + e− $\to$ Mg2+ + 2 e−	11
0033	Na + e− $\to$ Na+ + 2 e−	11
0034	He2+ + e− $\to$ He+	12
0035	He+ + H $\to$ He + H+	12
0036	He + H+ $\to$ He+ + H	12

Notes. Radiative recombination rates for all species are taken from Badnell (2006b). (1) Badnell (2006a), (2) Bautista & Badnell (2007), (3) Colgan et al. (2004), (4) Colgan et al. (2003), (5) Altun et al. (2004), (6) Zatsarinny et al. (2004b), (7) Mitnik & Badnell (2004), (8) Zatsarinny et al. (2004a), (9) Altun et al. (2006), (10) Verner & Ferland (1996), (11) Voronov (1997), and (12) Glover & Abel (2008).

Download table as:  ASCIITypeset image

Our implementation of this network follows the overall approach we have adopted in previous studies (Gray & Scannapieco 2010, 2013; Scannapieco et al. 2012). Each species is labeled by an index i, such that species i has Z_i protons, A_i nucleons, and a mass density of ρ_i. We then define a mass fraction of species i as X_i ≡ ρ_i/ρ, where $\rho ={{\sum }_{i}}{{\rho }_{i}}$, and define the molar mass fraction as Y_i ≡ X_i/A_i. For each species, a continuity equation for the molar mass fraction is then given as ${{\dot{Y}}_{i}}={{R}_{i}}$, where R_i is the total reaction rate due to all reactions.

Because of the complex ways that reaction rates depend on temperature and because of the large range of possible abundances, the resulting rate equations are often "stiff," i.e.,the change in timescales can be very different from one species to another. This requires implicit or semi-implicit methods to track all the relevant reactions throughout our simulations. To handle this, we use a variable-order Bader–Deuflhard method (e.g., Bader & Deuflhard 1983) which is well suited to problems with dimension >10 (e.g., Press et al. 1992; Bovino et al. 2013). The network presented here also has the advantage of being sparse. That is, the Jacobian, defined as ${{J}_{i,j}}\equiv \partial {{\dot{Y}}_{i}}/\partial {{Y}_{j}}$, is mostly comprised of zeros, with nonzero values falling on or near the diagonal. This allows us to use the MA28 sparse linear algebra package included with FLASH (Duff et al. 1986) to compute the matrix inverses required by the Bader–Deuflhard scheme, leading to a very fast and efficient implementation. In general, the chemistry package runs slightly faster than the hydrodynamics.

As the species evolve over a given time step, the temperature can change as the internal energy changes from ionizations and recombinations. Since the reaction rates are strong functions of temperature, the network can become unstable if too large a time step is used, especially in regions of strong cooling. To ensure testability but allow the simulation to evolve at the hydrodynamic time step, we subcycle the rate equations within each cell on a network time step, defined as:

$$\begin{eqnarray}&&{{t}_{{\rm net}}}\equiv \ {\rm min} \left( \alpha \frac{{{Y}_{i}}+0.1Y_{{\rm H}}^{+}}{{{{\dot{Y}}}_{i}}} \right),\end{eqnarray} \tag{ 1 }$$

where α is a constant that controls the maximum abundance change allowed, which we default to 0.1. Using the current species abundances, the instantaneous change in abundances ${{\dot{Y}}_{i}}$ is computed using the analytic ordinary differential equations at the current temperature. Note that a small fraction of $Y_{{\rm H}}^{+}$ is added as a buffer to prevent species with very small, but rapidly changing, abundances from causing prohibitively small time steps and excessive subcycling. In addition to this network subcycling, the Bader–Deuflhard method includes its own internal subcycling and time step controls.

At the beginning of each cycle, t_net is compared to the hydrodynamic time step. If the hydrodynamic time step is smaller, then no subcycling is done and the hydrodynamic time step is used to update the molar mass fractions. On the other hand, if the network time step is smaller, then the molar mass fractions are updated using the t_net, which is then recalculated and compared to the remaining hydrodynamic time step. This cycle continues until the network has advanced for the full hydrodynamic time step.

To test the implementation of our network, we emulated a series of equilibrium models by fixing the density and temperature in the simulation and running each case until the species reached equilibrium. We ran a series of seven such models that spanned the temperature range between 10⁴ and 10⁷ K, at a fixed density of 2.0 × 10⁻²⁰ g cm⁻³. All species were assumed to be neutral at the start of each simulation, and the final abundances were compared to equilibrium models from Cloudy (version 10.01) (Ferland et al. 1998), using the "coronal equilibrium" command that enforces only collisional ionizations.

We found that we matched the results from Cloudy very well over a wide range in temperature for all species. However, we did not match precisely for certain higher ionization states, as seen, for example, in the middle panel of Figure 1, but this is to be expected since our network does not follow every ionization state for some elements. We did find excellent agreement for those elements for which the ionization state changes by orders of magnitude with a modest increase in temperature. For the purposes of comparison, we summed together the higher ionization states that we do not follow and assigned their abundance to the highest state that we do follow. For example, we track the abundance of C⁴⁺ and not the two higher states. Therefore, the (dashed) navy line in the middle left panel of Figure 1 represents the summation of C⁴⁺ , C⁵⁺, and C⁶⁺.

Zoom In Zoom Out Reset image size

The second capability we implemented is individual cooling rates associated with each ion in the chemical network. To do this, we followed the procedure presented in Gnat & Ferland (2012) who compiled the ion-by-ion cooling efficiencies for several atomic species for between 10⁴ and 10⁸ K. These efficiencies include all the cooling processes considered in Cloudy, as described in Osterbrock & Ferland (2006). This includes collisional excitations followed by line emission, recombinations with ions, collisional ionizations, and thermal bremsstrahlung. Here we briefly outline this procedure, which we used to extend their work down to 5 × 10³ K for the ions in our network. Again, making use of Cloudy, we construct a grid of models that span 5 × 10³–10⁸ K in temperature for each ion i of element E. Each model is then composed of only hydrogen, electrons, and the ion under consideration. To ensure cooling from hydrogen is suppressed, the hydrogen density is set to 10¹⁰ times smaller than the ion density. Following Gnat & Ferland (2012), the elemental ion and electron density is set to ${{n}_{{{E}_{i}}}}$ = ${{n}_{{{e}^{-}}}}$ = 1.0 cm⁻³ and n_H = 10⁻¹⁰ cm⁻³. This produced a cooling rate for each ion as a function of temperature, and a comparison of these rates to those of Gnat & Ferland (2012) show identical results where the temperature range overlapped.

For those species in which we do not follow every ionization state in our simulation individually (such as C⁴⁺, which captures the joint abundance of C⁴⁺, C⁵⁺, and C⁶⁺) we formed a composite cooling curve as

$$\begin{eqnarray}&&{\Lambda }(T)=\frac{\mathop{\sum }\limits_{i}{{n}_{i}}(T){{{\Lambda }}_{i}}(T)}{\mathop{\sum }\limits_{i}{{n}_{i}}(T)},\end{eqnarray} \tag{ 2 }$$

where Λ_i(T) is the cooling rate of ion i and n_i(T) is the relative abundance of ion i in equilibrium at a temperature T. This composite is then used as the cooling rate for the highest ionization state followed. Similarly, we included cooling from nitrogen, neon, silicon, sulphur, calcium, and iron. Again, we assumed that the relative ionization abundances were determined by the collisional ionization equilibrium. This final cooling curve has the form

$$\begin{eqnarray}&&{{{\Lambda }}_{{\rm other}}}(T)=\mathop{\sum }\limits_{j}\frac{\mathop{\sum }\limits_{i}{{n}_{j,i}}(T){{{\Lambda }}_{j,i}}(T)}{\mathop{\sum }\limits_{i}{{n}_{j,i}}(T)},\end{eqnarray} \tag{ 3 }$$

where j loops over nitrogen, neon, silicon, sulphur, calcium, and iron.

With the addition of cooling, an important timescale is created, namely one that relates the total internal energy to energy loss rate per time as

$$\begin{eqnarray}&&{{t}_{{\rm cool}}}=\frac{{{\alpha }_{c}}{{E}_{i}}}{{{{\dot{E}}}_{i}}},\end{eqnarray} \tag{ 4 }$$

where E_i is the internal energy, ${{\dot{E}}_{i}}$ is the energy loss rate, and α_c is a constant between 0 and 1. As is the case with reaction rates, cooling rates are strong functions of temperature and species abundances, both of which can change drastically over a single chemistry time step. Therefore, we employ a method of subcycling the cooling on a cooling timescale within the chemistry routine. At the beginning of each cooling cycle, we calculate the cooling timescale using a given α_c(= 0.1). If this is shorter than the chemistry time step, we cycle over the cooling timescale, recalculating the temperature, cooling rate, and timescale at the end of each cycle. This proceeds until we have reached the chemistry time step.

Since we only consider two body interactions, the cooling time should scale linearly with density. To test this scaling, along with the overall dependence of our routines on temperature, we ran several tests in which the medium was initially completely ionized and the initial temperature was set to 1 × 10⁶ K, each of which had a different density between 5.0 × 10⁻²⁷ and 1.0 × 10⁻²⁵ g cm⁻³. These tests recovered the expected scaling with density, and also confirmed that cooling was always efficient at high temperatures and much less efficient at ≈10⁴ K, when most elements become neutral. Finally, we performed additional tests to ensure that the interpolation was smooth and always reproduced the table values under the correct conditions.

Our goal is to study the steady-state ionization structure of a gas that is being stirred with 1D velocity dispersions between 6 and 58 km s⁻¹ (corresponding to 3D velocity dispersions between 10 and 100 km s⁻¹), over a wide range of densities. In particular, we are interested in cases in which the heating from the turbulent stirring is balanced by atomic cooling.

The parameter space spanned by our simulations is greatly simplified by the dependencies of turbulent decay and cooling on density and length scale. For a given steady state, the average turbulent energy dissipation rate per unit volume, or consequently, the average heating rate per unit volume is

$$\begin{eqnarray}&&{\rm \bar{\ \Gamma\ }}=\bar{\rho }{{\sigma }^{3}}/L\ {\rm erg}\;{{{\rm s}}^{-1}}{\rm c}{{{\rm m}}^{-3}},\end{eqnarray} \tag{ 8 }$$

where σ is the velocity dispersion, $\bar{\rho }$ is the average density, and L is the driving scale of the turbulence. Conversely, the average cooling rate per unit volume is

$$\begin{eqnarray}&&{\rm \bar{\ \Lambda\ }}\propto {{\bar{\rho }}^{2}}\bar{\lambda }(T)/{{(\mu {{m}_{{\rm H}}})}^{2}}\ {\rm erg}\ {{{\rm s}}^{-1}}{\rm c}{{{\rm m}}^{-3}},\end{eqnarray} \tag{ 9 }$$

where μ is the average particle mass and $\bar{\lambda }(T)$ is the average cooling rate. Equating these two terms gives

$$\begin{eqnarray}&&\bar{\lambda }(T)\propto {{\sigma }^{3}}/(L\rho ),\end{eqnarray} \tag{ 10 }$$

or in other words that the overall thermal distribution of the gas is only dependent on σ and the column density. For the velocity dispersions we are interested in, for which cooling balances heating, this column density ranges between 10¹⁶ and 10²⁰ cm⁻².

Furthermore, because collisional ionization and recombination rates are also proportional to density squared, the ratio of the eddy turnover time defined here as ${{t}_{{\rm eddy}}}={{L}_{{\rm box}}}/2\sigma ,$ to the density-temperature averaged recombination time, t_rec, is also purely a function of σ and column density, with ${{t}_{{\rm rec}}}\geqslant {{t}_{{\rm eddy}}}$ for many species. This means that many of the species of interest experience a wide range of conditions within a recombination time, which, as we shall see below in more detail, can lead to abundance ratios that can not be predicted from local collisional equilibrium.

A summary of the simulation runs is given in Table 2, with the name of each run referring to its column density and 1D velocity dispersion. In the analysis that follows, we have chosen three representative models to describe in detail: N1E16_S6, N1E17_S20, and N1E19_S58, which span a range of densities and velocity dispersions, We will refer to these as Low, Medium, and High respectively, as they represent cases with progressively higher velocities and turbulent Mach numbers. Figure 2 shows the hydrodynamic and chemical evolution of N1E17_S20 in units of the eddy turnover timescale. The effect of the chemical kick is apparent at t_eddy ≈5 as the large jump in the fractional abundances of each species, e.g., ${{F}_{C}}={{F}_{{{C}_{i}}}}/{{\sum }_{i}}{{F}_{{{C}_{i}}}}$. After this adjustment, the species quickly find a global steady state that is distinct from instantaneous collisional ionization equilibrium. To ensure that the chemical kick did not introduce any bias in our model, we ran model N1E17_S20 without the kick. The statistical and atomic properties compared very well between these two models, although the run without the kick took over four times longer to reach a final steady-state.

Zoom In Zoom Out Reset image size

Table 2.  Summary of the Models

Name	$\bar{\rho }$	Column Density	${{\sigma }_{1{\rm D}}}$	T $_{{\rm MW}}$	TVW	MMW	MVW	$\bar{s}$	σs	sskew	skurt	σs,exp
N1E14_S6*	7e-31	1.1e+14	5.8	5.52	5.51	0.29	0.30	0.004	0.05	−0.60	1.01	0.15
N1E15_S6	7e-30	1.1e+15	5.8	2.55	2.54	0.29	0.30	0.005	0.05	−0.50	1.74	0.15
N1E16_S6 [Low]	7e-29	1.1e+16	5.8	1.49	1.48	0.51	0.51	0.005	0.13	−0.84	1.90	0.25
N1E17_S6	7e-28	1.1e+17	5.8	1.26	1.24	0.61	0.62	0.001	0.16	−0.79	1.48	0.30
N1E16_S12*	7e-29	1.1e+16	11.5	26.85	26.85	0.21	0.21	0.002	0.02	−0.54	1.95	0.10
N3E16_S12	2e-28	3.0e+16	11.5	7.01	6.99	0.44	0.45	0.006	0.08	−0.60	1.43	0.22
N1E17_S12	7e-28	1.1e+17	11.5	1.30	1.21	1.06	1.15	−0.061	0.46	−0.57	0.62	0.53
N3E17_S12	2e-27	3.0e+17	11.5	1.14	1.02	1.26	1.39	−0.119	0.61	−0.83	2.16	0.63
N7E16_S20*	5e-28	7.6e+16	20.2	34.05	33.99	0.36	0.37	0.006	0.06	−0.60	1.95	0.18
N1E17_S20 [Med]	1e-27	1.5e+17	20.2	5.54	5.31	0.91	0.96	−0.021	0.31	−0.68	1.16	0.46
N7E17_S20	5e-27	7.6e+17	20.2	1.01	0.82	2.39	2.95	−0.381	0.98	−0.19	−0.27	1.07
N1E18_S20	1e-26	1.5e+18	20.2	0.87	0.69	2.78	3.56	−0.497	1.11	−0.21	−0.11	1.20
N4E17_S35*	3e-27	4.6e+17	34.6	176.12	175.95	0.31	0.31	0.004	0.05	−0.49	1.34	0.16
N1E18_S35	1e-26	1.5e+18	34.6	1.12	0.89	4.13	5.25	−0.798	1.42	−0.33	0.11	1.44
N4E18_S35	3e-26	4.6e+18	34.6	0.90	0.73	4.63	6.49	−0.877	1.52	−0.37	0.16	1.56
N1E19_S35	1e-25	1.5e+19	34.6	0.83	0.64	5.32	7.65	−0.844	1.43	−0.17	−0.10	1.66
N3E18_S58*	2e-26	3.0e+18	57.7	67.55	64.53	0.91	0.96	−0.005	0.32	−0.54	0.73	0.46
N1E19_S58 [High]	7e-26	1.1e+19	57.7	1.06	1.09	7.98	9.83	−1.196	1.77	−0.32	0.07	1.80
N3E19_S58	2e-25	3.0e+19	57.7	0.86	0.71	8.79	13.55	−1.375	1.93	−0.32	−0.16	1.96
N1E20_S58	7e-25	1.1e+20	57.7	0.77	0.64	10.64	14.51	−1.441	2.01	−0.46	0.33	2.00
N1E16_S6_High	7e-29	1.1e+16	5.8	1.49	1.48	0.45	0.45	0.009	0.10	−0.44	1.12	0.22
N1E16_S6_C	7e-29	1.1e+16	5.8	1.50	1.49	0.49	0.49	0.007	0.12	−1.11	3.35	0.24
N1E17_S20_High	1e-27	1.5e+17	20.2	5.66	5.43	0.89	0.94	−0.021	0.32	−0.68	0.90	0.45
N1E17_S20_C	1e-27	1.5e+17	20.2	5.61	5.35	0.88	0.93	−0.028	0.34	−0.64	0.84	0.44
N1E19_S58_High	7e-26	1.1e+19	57.7	0.95	0.81	7.68	10.64	−1.173	1.74	−0.23	−0.22	1.84
N1E19_S58_C	7e-26	1.1e+19	57.7	1.07	1.14	7.36	10.41	−1.303	1.94	−0.57	0.66	1.83

Notes. Those models with the appended _C denote runs made with thermal conduction while those appended with _High denote high resolution models with twice the normal resolution. Those denoted with an asterisk denote models that result in thermal runaway. $\bar{\rho }$ is the mean density in units of g cm⁻³. Column density is the column density in units of cm⁻². ${{\sigma }_{1{\rm D}}}$ is the 1D velocity dispersion in units of km s⁻¹. T_MW and T_VW are the mass-weighted and volume-weighted temperatures in units of 10⁴ K. M_MW and M_VW are the mass-weighted and volume-weighted turbulent Mach numbers. $\bar{s}$ is the volumed averaged value of $s\equiv {\rm ln} \rho /\bar{\rho }.$ σ_s, s_skew and s_kurt are the rms, skewness and kurtosis excess of the volume-weighted probably density function of s. σ_{s, exp} is the expected variance from Equation (12).

Download table as:  ASCIITypeset image

Figure 3 shows the density, temperature, and the fraction of carbon in the commonly observed C³⁺ state on slices extracted from each of our representative models after they have reached a global steady state. Here we see that while density and temperature variations are small at low turbulent Mach numbers, the more highly turbulent simulations display a wide range of temperatures and densities, which are only weakly correlated with each other. Similarly, the distribution of the fraction of C³⁺, which has a recombination time that is comparable to t_eddy, displays features that sometimes follow the temperature distribution, but sometimes show substantial variations, indicating that F$^{{{{\rm C}}^{3+}}}$ is often very different than would be estimated from local collisional ionization equilibrium. Below we study each of these effects in turn.

Zoom In Zoom Out Reset image size

The top panel of Figure 4 shows the two-dimensional volume-weighted probability density functions, which quantify the fraction of the volume in the simulations located at various temperatures and densities. Consistent with the expectations from Figure 3, the gas in the Low simulation has a very small spread of temperatures and densities. In the Medium simulation, on the other hand, the temperatures and densities span roughly an order of magnitude. In this case, for which the (mass weighted) turbulent Mach number is M_MW = 0.91, T is roughly proportional to n⁻¹, as would occur exactly in a constant pressure medium, but occurs here with a considerable spread, consistent with the presence of transonic motions. Finally, in the High case, with a Mach number of M_MW = 8.0, the densities and temperatures span nearly six orders of magnitude and four orders of magnitude, respectively. Furthermore, these quantities are largely uncorrelated with each other, as would be expected in a medium filled with strong, supersonic shocks.

Zoom In Zoom Out Reset image size

While our simulations are the first to include full rate equations and cooling for a large number of atomic species, several studies of isothermal, supersonic turbulence have shown that the gas approximates a lognormal distribution (Vazquez-Semadeni 1994; Padoan et al. 1997; Klessen 2000; Ostriker et al. 2001; Li et al. 2003; Kritsuk et al. 2007; Federrath et al. 2008, 2010; Lemaster & Stone 2008; Schmidt et al. 2009; Glover et al. 2010; Collins et al. 2011; Padoan & Nordlund 2011; Price et al. 2011; Molina et al. 2012), defined as

$$\begin{eqnarray}&&{{p}_{s}}\ ds=\frac{1}{\sqrt{2\pi \sigma _{s}^{2}}}{\rm exp} \left[ -\frac{{{(s-{{s}_{0}})}^{2}}}{2\sigma _{s}^{2}} \right],\end{eqnarray} \tag{ 11 }$$

where s is the logarithmic density, $s\equiv {\rm ln} (\rho /\bar{\rho }).$

The mean logarithmic density in this case is related to the standard deviation as ${{s}_{0}}=-\sigma _{s}^{2}/2$. The density variation at a particular location is produced by the successive passage of shocks with Mach numbers independent of the local density, which gives a physical explanation for the lognormal density distribution. For an isothermal distribution, the variance of the logarithmic density corresponds to

$$\begin{eqnarray}&&\sigma _{s}^{2}={\rm ln} (1+{{b}^{2}}\;{{M}^{2}}),\end{eqnarray} \tag{ 12 }$$

where b is a constant that depends on the forcing that drives the turbulence. Federrath et al. (2008) showed that b = 1 for purely compressive, $\nabla \times {\boldsymbol{F}} =0,$ forcing, while b = 1/3 for purely solenoidal, $\nabla \cdot {\boldsymbol{F}} =0$, forcing.

The bottom, left panel of Figure 4 shows the probability density functions of the logarithmic density, which take on a Gaussian profile for runs Low, Med, and High. As expected, the width of the profile increases as the steady state Mach number increases, with Low having the smallest width and High having the largest. In Table 2 we show several important statistical quantities in terms of the logarithmic density for each run. In particular, we calculate the first four moments of the logarithmic density distribution: the mean, variance, skewness, and kurtosis.

The mean and variance have the usual definitions, while the skewness, $\left\langle {{(s-{{s}_{0}})}^{3}} \right\rangle /{{\sigma }^{3/2}}$, measures the symmetry of the distribution. A positive skewness denotes a distribution that has a long tail toward higher densities while a negative skewness represents a long tail toward lower densities. The kurtosis, $\left\langle {{(s-{{s}_{0}})}^{4}} \right\rangle /{{\sigma }^{2}}$ quantifies the peakedness of the distribution and measures the importance of the tails versus the peak. A larger kurtosis represents a flatter distribution while a lower value denotes a narrower distribution, with a Gaussian distribution having a kurtosis equal to 3. In Table 2 we give the kurtosis excess where this factor of 3 is subtracted.

Figure 5 shows the variance, skewness, and kurtosis excess as a function of Mach number for each model. The solid line in the top panel shows the expected σ_s–M relation from Equation (12), with b = 0.5. We find that for supersonic flows the models match very well with Equation (12) but find a significant mismatch for subsonic flows. A least-squares fit of Equation (12) to our points suggests a good fit to the $M\gtrsim 1$ points with b = 0.53. However, this does not match the expectation from Federrath et al. (2008) who suggest a value of b = 1/3 for purely solenoidal driving, for isothermal turbulent models. This difference may simply be due to the non-isothermal nature of our models. The local variation in temperature allows for a larger spread in density than in isothermal models.

Zoom In Zoom Out Reset image size

The right panel of Figure 4 show the 1D temperature PDFs for runs Low, Med, and High. As expected, the width of the PDF grows with the strength of the stirring, meaning that isothermality is not a good approximation for our high Mach number runs, which may help explain the discrepancy in b. Note also the large temperature spread causes a few cells in run High to reach the temperature floor, due to the extreme PdV work done by the stirring, but this is only a minor effect.

The top row of Figure 4 shows the two-dimensional PDFs of density and temperature for runs Low, Med, and High. For the Low case, the slight stirring has only produced very minor density contrasts, leaving a mostly uniform medium. However, there is a slight tail toward lower density and lower temperature. This leads to a highly concentrated phase diagram where essentially all the mass has the same temperature. This can also be seen in the left column of Figure 3 where we show slices of density, temperature, and the abundance of C³⁺ .

A similar trend is seen for run Med in the center panel of Figure 4 and the center column of Figure 3. Although, as mentioned below, the average Mach number of the flow is higher than Low; turbulence has not produced any strong density contrasts. In fact, the density is almost always within a factor of 3 of the mean. As in the case for Low, the density and temperature slices are largely uniform.

Finally, the High case is quite different than either Low or Med. Large density and temperature contrasts are seen, due to the much higher final Mach number (M = 8.0). This produces a very wide density-temperature PDF as shown in the right panel of Figure 4 with density and temperature spanning nearly six orders of magnitude and four orders of magnitude respectively. These contrasts are plainly seen in the right column of Figure 3.

The bottom row of Figure 3 shows slices of F$_{{{{\rm C}}^{3+}}}$. In general the density and temperature distributions are very similar. For example, the density and temperature slices for High show that the complex density structure is largely mirrored in the temperature. However, this is not the case for C³⁺ . In fact, there is only a weak relation between these quantities. As we discuss below, this highlights the need to follow the nonequilibrium atomic chemistry exactly when estimating atomic ionization states.

The large spread in both density and temperature in our simulations, particularly in the supersonic cases, naturally translates into a large spread in chemical species. In fact, the instantaneous chemical makeup in a given cell is a strong function of not only the current temperature and density but its thermodynamic history.

To get an idea of how the global steady state abundances depend on our full treatment of turbulence, cooling, and chemical reactions, we calculated approximate abundances using two approaches. For our first method, labeled F_Mean, we assumed the full medium was in collisional ionization equilibrium at the mass weighted average temperature, and we calculated the ionization states for each species using Cloudy. This method quantifies the errors involved with assuming constant temperature and density conditions in a region with significant turbulent motions. In our second estimate, labeled F_TurbEq, we used the final temperature PDF to calculate a mass-weighted average abundance for each species, again using Cloudy to calculate collisional ionization abundances for the full range of conditions. This method quantifies the errors involved with failing to track the nonequillibrium abundances accurately by including reaction rates within the simulations.

The left panel of Figure 6 shows the comparison between each of these estimates to the steady state abundances for Low, Med, and High. The name of each species is listed along the x-axis and the fractional species abundances are given along the y-axis. These values are also listed in Table 3. As shown above, for the Low simulation, with a mean temperature of 1.5 × 10⁴ K and a Mach number of 0.5, the spread in the temperature-density PDF is small. Hence, the differences between the global species abundances and those from our two estimates are mostly minor, with most of the true values matching the simpler calculations within a factor of ≈2 for most elements with significant abundances. Two notable exceptions to this agreement are C²⁺ for which the true abundance is ≈30 times greater than the F_Mean and F_TurbEq, values, and He⁺, which has a long recombination time, and a true abundance ≈20 times greater than expected from more simple estimates.

Zoom In Zoom Out Reset image size

In the Med case, for which the mass-weighted mean temperature is T_MW = 5.5 × 10⁴ and M = 0.92, the difference between the true values and simple estimates becomes more pronounced. Here the true abundances of species such as He, C, C⁺, C³⁺, O³⁺, and Mg, differ by a factor of ≈2 from our simple estimates, while He²⁺ and C⁴⁺ are more abundant by ≈5 than expected from either approximate approach. Note also that there is no clear trend for the more complicated F_TurbEq estimate to be a better predictor than the simple constant temperature approach.

Finally, in the High case, with ${{T}_{{\rm MW}}}=1.1\times {{10}^{4}}$ K and M = 8.0, the abundances of many species differ by several order of magnitude from simple expectations. For example, the abundance of H⁺ is over two orders of magnitude higher than it is in the constant temperature models, but close to the F_TurbEq estimate. The abundance of He, on the other hand, is very near the value expected from F_Mean and roughly twice the F_TurbEq estimate, while and He²⁺ is many orders of magnitude higher than expected by the F_Mean estimate, but still over 200 times less abundant than expected in the F_TurbEq estimate. Similarly, C²⁺, C³⁺, C⁴⁺, O²⁺, O³⁺, and O⁴⁺ which are almost completely absent in the F_Mean, are all present at significant levels, although less so than expected by F_TurbEq. On the other hand, Mg⁺ is present at a level much less than expected by the F_Mean estimate, but more than predicted by F_TurbEq. These results not only show that the presence of turbulence can introduce significant abundance differences when M ≈ 1, and change the abundances completely when M ≳ 1, but that these differences can only be predicted by full nonequilibrium calculations such as the ones carried out here.

To get an idea of the cause of the large discrepancy between carbon and oxygen in the High case, we compare the ionization and recombination timescales of these ions. This is shown in the right panel of Figure 6. We find that while the recombination time is short compared to the eddy turnover, the ionization timescale for these ions is long. Physically, this suggests that even though a parcel of gas is heated due to the crossing shocks, the ions do not have sufficient time to fully ionize. This explains the lower abundances found for F_TurbEq compared to F_Mean for C²⁺ through C⁴⁺ and O³⁺ through O⁷⁺ for case High.

In our simulations described above, energy transport is purely by advection. Yet even in a supersonic medium, the large difference in the thermal velocities of the electrons and the ions can lead to significant energy transport through conduction. In the case in which many collisions take place over the scale length for temperature variations, the thermal conductivity is given by

$$\begin{eqnarray}&&\kappa =\nu {{\left( \frac{{{m}_{p}}}{{{m}_{e}}} \right)}^{1/2}}=5.2\times {{10}^{-5}}T_{4}^{5/2}{{n}^{-1}}\qquad {\rm c}{{{\rm m}}^{2}}\ {{{\rm s}}^{-1}},\end{eqnarray} \tag{ 13 }$$

where ν is the plasma viscosity and T is the temperature in units of 10⁴ K (Braginskii 1958; Spitzer 1962).

Comparing this to the velocity dispersion and length scale of the turbulence gives ${{\sigma }_{1{\rm D}}}L/\kappa \propto {{\sigma }_{1{\rm D}}}nL$. Like the ratios of timescales discussed in Section 4.1, this is dependent only on the column density and turbulent velocity. In fact, the Re of the medium also only depends on these two quantities, although in practice the numerical viscosity in our simulations is much greater than the true physical value.

In the case in which the scale length for temperature variations is smaller than the collisional mean free path of the electrons, the conduction becomes saturated. The mean free path is ${{\lambda }_{i}}=1.3\times {{10}^{18}}\;{\rm c}{{{\rm m}}^{-2}}T_{7}^{2}/{{n}_{i,c}}$, where n_i,c is the ion density, the saturated flux is ${{F}_{{\rm sat}}}{{\alpha }_{{\rm ele}}}{{n}_{e}}\sqrt{{{(kT)}^{3}}/{{m}_{e}}},$ where α_ele is the electron conductivity flux-limiter coefficient, n_e is the electron number density, and m_e is the mass of the electron (Cowie & McKee 1977). From these scalings we can see that L/λ_i and ${{F}_{{\rm sat}}}/nkT(\sigma )$ are purely functions of column density and the temperature structure of the medium. This means that neither unsaturated or saturated conduction introduces a new parameter that must be spanned by our simulations.

To test the impact of electron conduction, we re-ran our representative simulation using the implicit SpitzerHighZ conductivity module within FLASH, a Larsen flux limiter, and an electron conductivity flux-limiter coefficient of α_ele = 0.2. The upper panel of Figure 7 compares the logarithmic density PDFs from our nominal models, with these conductions runs, which were roughly twice as expensive in terms of computer time. For run Low the PDF profile is essentially unchanged from the nominal model, in particular with regards to the mean and variance of s. Table 2 presents the statistical measures between our nominal runs and those with electron conduction. There is a slight increase in the skewness and kurtosis which produces a slightly longer distribution to lower densities and slightly more peaked distribution. Run Med shows no substantial change in the kurtosis, but conduction does produce a slightly less peaked distribution. Run High shows a similar trend as the other models, with a slight increase in both the skewness and kurtosis, but now with a noticeable difference in the extreme low-density tail of the log density PDF.

Zoom In Zoom Out Reset image size

The top right panel of Figure 7 shows the effect of electron conduction on the logarithmic temperature PDF. Like the density PDFs, the temperature PDFs are largely identical for Low and Med. Only for the High simulation is there a small decrease in the PDF at the highest temperatures at which conduction is most efficient, which is accompanied by an increase in the low-density end of the density PDF. Thus it appears that in the highest Mach number cases, conduction is able to operate quickly enough to move a noticeable amount of energy out of the hottest cells, moving some of this energy into low density regions which are able to expand a bit more due to the increase in pressure. However, these differences are small and suggest that electron conduction does not play a major role.

In order to simulate a large number of turbulent conditions for many eddy turnover times in the presence of chemistry at reasonable CPU cost, our simulations have adopted a relatively low resolution of 128³. In this case, two-point statistics such as the velocity and density power-spectra and structure functions are not expected to be well reproduced (e.g., Falkovich 1994; Kritsuk et al. 2007; Pan & Scannapieco 2010), and we do not attempt to analyze them here. On the other hand, the probability distribution function is a single-point quantity that is much easier to simulate at moderate resolution.

As a check of convergence, the bottom left row of Figure 7 compares the PDF obtained in our representative simulations with a second set of simulations with the same set of physical parameters, but a resolution of 256³. Again, we find that the resulting PDFs are nearly identical to the nominal runs and only in the High case is there a slight decrease in the high temperature tail and a slight dip in the low density tail. However, in terms of the mean and variance, the effect is minor, which suggests that our nominal model with N = 128 is sufficient to produce the key one-point statistical quantities and obtain estimates of the species abundances.

Finally, in Figure 8 we show the difference in the species abundances between our nominal runs and the comparison runs including electron conduction and an increase in resolution. The comparisons, which are also presented in Table 4, show that neither electron conduction nor an increase in resolution has a substantial effect on our abundance results.

Zoom In Zoom Out Reset image size

Table 4.  Effect of Resolution and Electron Conduction of the Final Elemental Abundances

Col. (cm−2)	1.1E16	1.1E16_High	1.1E16_C	1.5E17	1.5E17_High	1.5E17_C	1.1E19	1.1E19_High	1.1E19_C
${{\sigma }_{1{\rm D}}}$ (km s−1)	5.8	5.8	5.8	20.2	20.2	20.2	57.7	57.7	57.7
H	−0.23	−0.22	−0.15	−3.47	−3.50	−3.47	−0.37	−0.27	−0.29
H+	−0.39	−0.40	−0.52	0.00	0.00	0.00	−0.24	−0.34	−0.31
He	−0.00	−0.00	−0.00	−2.54	−2.59	−2.56	−0.03	−0.02	−0.05
He+	−3.26	−3.34	−3.71	−0.07	−0.07	−0.07	−1.14	−1.30	−0.95
He2+	−17.28	−17.39	−17.92	−0.86	−0.81	−0.84	−3.07	−3.14	−2.48
C	−0.67	−0.69	−0.60	−3.10	−3.14	−3.11	−0.37	−0.34	−0.36
C+	−0.11	−0.10	−0.13	−0.73	−0.75	−0.74	−0.25	−0.27	−0.26
C2+	−3.67	−3.60	−3.80	−0.10	−0.09	−0.09	−1.75	−2.26	−1.92
C3+	−14.19	−14.15	−18.30	−1.98	−1.91	−1.95	−3.31	−3.77	−3.20
C4+	−18.56	−18.41	−18.13	−3.47	−3.18	−3.63	−4.15	−4.26	−3.59
O	−0.28	−0.26	−0.17	−3.07	−3.10	−3.07	−0.37	−0.27	−0.30
O+	−0.33	−0.35	−0.50	−0.44	−0.46	−0.44	−0.25	−0.34	−0.31
O2+	−7.62	−8.05	−13.47	−0.21	−0.20	−0.21	−2.02	−2.38	−1.94
O3+	−17.97	−17.82	−17.52	−1.88	−1.85	−1.93	−3.20	−3.62	−3.02
O4+	−18.92	−18.73	−18.44	−5.92	−5.85	−6.01	−4.73	−4.81	−4.54
O5+	−18.72	−18.57	−18.29	−11.67	−11.25	−11.84	−6.21	−5.96	−6.09
O6+	−18.85	−18.66	−18.49	−13.14	−10.67	−16.48	−7.53	−7.05	−6.88
O7+	−19.68	−19.17	−19.82	−20.89	−19.84	−20.17	−15.17	−17.10	−17.03
Na	−4.26	−4.29	−4.29	−5.86	−5.85	−5.85	−2.10	−2.12	−2.22
Na+	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00
Mg	−1.40	−1.44	−1.36	−6.37	−6.39	−6.36	−0.65	−0.60	−0.65
Mg+	−0.18	−0.17	−0.12	−3.41	−3.44	−3.41	−0.31	−0.24	−0.24
Mg2+	−0.53	−0.54	−0.72	−0.00	−0.00	−0.00	−0.55	−0.76	−0.69

Notes. Each row gives a difference species while each column gives a different simulation. _High denotes models with twice the base resolution and those with _C denote models with electron conduction.

Download table as:  ASCIITypeset image

Having explored the evolution of three representative cases in detail, established convergence of abundances at a resolution at 128³, and shown that conduction has only a minor effect, we turn our attention to the full suite of models. Table 2 gives the overall parameters of each of these runs, which span a large range of columns, temperatures, and Mach numbers. In most cases, we found that heating from stirring and cooling from recombination are quickly balanced, as seen in the Med simulation presented in Figure 2. However, if the stirring is sufficiently strong that the average temperature of the model exceeds 10⁵ K these models will undergo thermal runaway. This can be explained by the fact that most elements have peaks in their cooling functions at ≈10⁵ K (e.g., Gnat & Ferland 2012; Oppenheimer & Schaye 2013). Therefore once the stirring passes this barrier, cooling can no longer balance turbulent heating and a meaningful steady state can no longer be achieved.

A summary of our results is given in Table 5, which shows the final abundances for all models that were able reach a steady state. Inspection of these abundances reveals many interesting trends affecting commonly observed species. For example, without any contribution from photoionization, He⁺ is often found at substantial levels, even when the average temperature is 10⁴ K, and in many of these cases the H⁺ fraction can exceed 50% (as seen when σ_1D = 35 km s⁻¹ and M_MW = 2.8, or when σ_1D = 58 km s⁻¹ and M_MW = 8.0).

Table 5.  Final Steady-state Abundances for Our Suite of Models

Col. (cm−2)	1.1E15	1.1E16	1.1E17	3.0E16	1.1E17	3.0E17	1.5E17	7.6E17	1.5E18	1.5E18	4.6E18	1.5E19	1.1E19	3.0E19	1.1E20
${{\sigma }_{1{\rm D}}}$ (km s−1)	5.8	5.8	5.8	11.5	11.5	11.5	20.2	20.2	20.2	34.6	34.6	34.6	57.7	57.7	57.7
TDW (104K)	2.55	1.49	1.26	7.01	1.30	1.14	5.54	1.01	0.87	1.12	0.90	0.83	1.06	0.86	0.77
MDW	0.29	0.51	0.61	0.44	1.06	1.26	0.91	2.39	2.78	4.13	4.63	5.32	7.98	8.79	10.6
H	−1.03	−0.23	−0.03	−3.96	−0.19	−0.10	−3.47	−0.17	−0.09	−0.42	−0.19	−0.09	−0.37	−0.15	−0.07
H+	−0.04	−0.39	−1.13	0.00	−0.45	−0.68	0.00	−0.48	−0.73	−0.21	−0.44	−0.73	−0.24	−0.53	−0.81
He	−0.00	−0.00	−0.00	−2.87	−0.17	−0.01	−2.54	−0.00	−0.00	−0.01	−0.00	−0.00	−0.03	−0.00	−0.00
He+	−2.20	−3.26	−4.58	−0.02	−0.49	−1.84	−0.07	−2.96	−3.68	−1.59	−2.72	−3.30	−1.14	−2.61	−2.94
He2+	−15.48	−17.28	−19.53	−1.27	−4.52	−7.37	−0.86	−9.96	−10.90	−5.22	−7.76	−9.17	−3.07	−5.66	−5.78
C	−1.51	−0.67	−0.21	−3.71	−0.67	−0.46	−3.10	−0.50	−0.28	−0.63	−0.34	−0.19	−0.37	−0.19	−0.12
C+	−0.02	−0.11	−0.41	−1.05	−0.11	−0.19	−0.73	−0.17	−0.33	−0.12	−0.27	−0.44	−0.25	−0.44	−0.61
C2+	−2.15	−3.67	−6.06	−0.05	−3.95	−4.54	−0.10	−4.02	−4.50	−2.46	−3.36	−4.03	−1.75	−3.30	−3.75
C3+	−10.10	−14.19	−20.35	−1.53	−12.45	−13.21	−1.98	−11.00	−10.54	−5.69	−7.42	−9.28	−3.31	−6.01	−5.65
C4+	−18.50	−18.56	−20.01	−3.96	−16.07	−21.04	−3.47	−14.78	−16.02	−7.24	−10.76	−13.39	−4.15	−6.89	−7.15
O	−1.25	−0.28	−0.05	−3.49	−0.50	−0.25	−3.07	−0.26	−0.10	−0.53	−0.20	−0.09	−0.37	−0.15	−0.08
O+	−0.03	−0.33	−0.99	−0.64	−0.17	−0.36	−0.44	−0.34	−0.70	−0.15	−0.43	−0.71	−0.25	−0.53	−0.79
O2+	−4.60	−7.62	−13.81	−0.13	−4.39	−7.06	−0.21	−5.74	−6.40	−3.17	−4.58	−5.47	−2.02	−4.20	−4.44
O3+	−18.07	−17.97	−17.73	−1.54	−14.15	−15.07	−1.88	−11.64	−11.91	−5.62	−8.00	−10.23	−3.20	−5.82	−5.70
O4+	−19.46	−18.92	−18.46	−5.51	−24.67	−25.59	−5.92	−20.29	−19.12	−8.90	−12.32	−16.24	−4.73	−7.68	−7.61
O5+	−19.54	−18.72	−18.22	−12.31	−26.39	−26.34	−11.67	−23.90	−26.03	−11.26	−16.89	−21.61	−6.21	−9.13	−10.90
O6+	−19.51	−18.85	−18.54	−18.62	−26.33	−26.33	−13.14	−26.34	−26.34	−14.61	−22.27	−26.33	−7.53	−10.96	−14.93
O7+	−21.43	−19.68	−19.64	−20.34	−26.54	−26.54	−20.89	−26.54	−26.55	−26.55	−26.56	−26.57	−15.17	−19.93	−23.88
Na	−5.31	−4.26	−3.87	−5.76	−3.88	−3.51	−5.86	−2.91	−2.45	−2.75	−2.19	−1.98	−2.10	−1.86	−1.57
Na+	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.00	−0.01	−0.01
Mg	−2.97	−1.40	−0.98	−7.23	−1.29	−1.04	−6.37	−0.88	−0.67	−0.98	−0.61	−0.48	−0.65	−0.44	−0.35
Mg+	−0.83	−0.18	−0.06	−3.84	−0.16	−0.09	−3.41	−0.13	−0.12	−0.25	−0.18	−0.20	−0.31	−0.25	−0.29
Mg2+	−0.07	−0.53	−1.57	−0.00	−0.60	−1.03	−0.00	−0.93	−1.42	−0.48	−1.00	−1.38	−0.55	−1.13	−1.46
${{{\rm H}}^{-*}}$	−8.39	−7.42	−7.39	−11.79	−7.36	−7.32	−11.13	−7.38	−7.36	−7.61	−7.42	−7.37	−7.58	−7.40	−7.40
${\rm H}_{2}^{+*}$	−7.08	−6.90	−7.51	−9.38	−6.99	−7.14	−9.07	−7.08	−7.20	−7.10	−7.11	−7.22	−7.18	−7.16	−7.30
${\rm H}_{2}^{*}$	−7.85	−6.24	−6.15	−14.08	−5.84	−3.90	−12.84	−3.16	−2.14	−2.86	−2.17	−1.84	−2.21	−1.86	−1.60

Notes. Abundance values are given as ${{{\rm log} }_{10}}({{F}_{i}}/{{F}_{T}})$. Here, ${{{\rm H}}^{-}}$, ${\rm H}_{2}^{+}$, ${{{\rm H}}_{2}}$ are computed according to Equations (14)–(16) and are thus approximate.

Download table as:  ASCIITypeset image

On the other hand, C³⁺ mostly appears at substantial levels when the mean temperature is significantly higher than 10⁴ K, such as in the ${{\sigma }_{1D}}=11.5$ km s⁻¹, ${{T}_{{\rm MW}}}=7.0\times {{10}^{4}}$ K, and ${{\sigma }_{1D}}=20$ km s⁻¹, T_MW = 5.5 × 10⁴ K runs. Note however, that it is also seen in the high-Mach number σ_1D = 58, km s⁻¹ run with T_MW = 1.1 × 10⁴ K, even though in CIE it is only expected when $T\geqslant 5\times {{10}^{4}}$ K as shown in Figure 1. Similarly, O³⁺, which has a slightly lower ionization potential, is most abundant in runs with elevated temperatures or high Mach numbers, although at even higher fractions than seen in C³⁺.

Interestingly, although its ionization potential is only 5.1 eV, substantial levels of neutral sodium are found in most of our simulations, particularly those with higher Mach numbers. Inspections plots of the physical distribution of this ion show that it is mostly found in under-dense regions, which drop to relatively low temperatures through adiabatic cooling in expansions. Finally, Mg⁺ which recombines at 7.6 eV and is ionized at 15 eV is extremely abundant in all our simulations, with the exception of the highest temperature ${{\sigma }_{1{\rm D}}}=11.5$ km s⁻¹ ${{T}_{{\rm MW}}}=7.0\times {{10}^{4}}$ K, and ${{\sigma }_{1D}}=20$ km s⁻¹ T_MW = 5.5 × 10⁴ K cases. In fact, in the absence of photoionization, the existence of significant, cospatial Mg⁺, C³⁺, and O²⁺ is a clear indicator of a broad temperature-density PDF, caused by supersonic turbulence.

At highest Mach numbers, when the dispersions in temperature and density are large, regions in the simulation are able to reach conditions in which the fraction of molecular hydrogen may not be negligible. Because parcels of gas move between regions with significantly different conditions on the timescale in which molecules are formed, calculating the H₂ fraction exactly would require solving a compete set of rate equations for H₂ formation within our simulations themselves. While this is beyond the scope of this work, we can nevertheless make a rough estimate of the molecular fractions calculating cell by cell abundances of H⁻, ${\rm H}_{2}^{+}$, and H₂ in the approximate case in which the local conditions last long enough for these additional species to come into equilibrium with the species tracked exactly by our simulations. In this case, in each cell we can compute

$$\begin{eqnarray}&&\begin{array}{l} {{F}_{{{{\rm H}}^{-}}}}\;= \\ \left( \begin{array}{ccccccccccccccc} \left( \begin{array}{ccccccccccccccc} {{k}_{1}}{{F}_{{\rm H}}}{{F}_{{{{\rm e}}^{-}}}}+{{k}_{23}}{{F}_{{\rm H}_{2}^{+}}}{{F}_{{{{\rm e}}^{-}}}} \\ \end{array} \right)/ \\ \left( \begin{array}{ccccccccccccccc} ({{k}_{2}}+{{k}_{15}}){{F}_{{\rm H}}}+({{k}_{5}}+{{k}_{16}}){{F}_{{{{\rm H}}^{+}}}}+{{k}_{14}}{{F}_{{{{\rm e}}^{-}}}} \\ \quad +\;({{k}_{21}}+{{k}_{22}}){{F}_{{\rm H}_{2}^{+}}}+{{k}_{28}}{{F}_{{\rm H}{{{\rm e}}^{+}}}}+{{k}_{29}}{{F}_{{\rm He}}} \\ \end{array} \right) \\ \end{array} \right) \\ \end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\begin{array}{l} {{F}_{{\rm H}_{2}^{+}}}\;= \\ \left( \begin{array}{ccccccccccccccc} \left( \begin{array}{ccccccccccccccc} {{k}_{3}}{{F}_{{\rm H}}}{{F}_{{{{\rm H}}^{+}}}}+{{k}_{7}}{{F}_{{{{\rm H}}_{2}}}}{{F}_{{{{\rm H}}^{+}}}} \\ \quad +\;{{k}_{16}}{{F}_{{{{\rm H}}^{+}}}}{{F}_{{{{\rm H}}^{-}}}}+{{k}_{25}}{{F}_{{{{\rm H}}_{2}}}}{{F}_{{\rm H}{{{\rm e}}^{+}}}} \\ \end{array} \right)/ \\ \left( {{k}_{4}}{{F}_{{\rm H}}}+{{k}_{6}}{{F}_{{{{\rm e}}^{-}}}}+({{k}_{21}}+{{k}_{22}}){{F}_{{{{\rm H}}^{-}}}} \right) \\ \end{array} \right) \\ \end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\begin{array}{l} {{F}_{{{{\rm H}}_{2}}}}= \\ \left( \begin{array}{ccccccccccccccc} \left( \begin{array}{ccccccccccccccc} {{k}_{2}}{{F}_{{{{\rm H}}^{-}}}}{{F}_{{\rm H}}}+{{k}_{4}}{{F}_{{\rm H}}}{{F}_{{\rm H}_{2}^{+}}} \\ \quad +\;{{k}_{21}}{{F}_{{{{\rm H}}^{-}}}}{{F}_{{\rm H}_{2}^{+}}}+{{k}_{R}}{{F}_{{\rm H}}}{{F}_{{\rm H}}} \\ \end{array} \right)/ \\ \left( \begin{array}{ccccccccccccccc} {{k}_{7}}{{F}_{{{{\rm H}}^{+}}}}+({{k}_{8}}+{{k}_{23}}){{F}_{{{{\rm e}}^{-}}}}+{{k}_{9}}{{F}_{{\rm H}}} \\ \quad +\;{{k}_{11}}{{F}_{{\rm He}}}+({{k}_{24}}+{{k}_{25}}){{F}_{{\rm H}{{{\rm e}}^{+}}}} \\ \end{array} \right) \\ \end{array} \right) \\ \end{array}\end{eqnarray} \tag{ 16 }$$

where the temperature dependent reaction rates, k₁, k₂, etc are taken from the compilation presented in appendix A of Glover & Abel (2008). Furthermore, to estimate the formation of H₂ from neutral–neutral surface reactions on dust we follow Glover & Jappsen (2007), Equation (36) which assumes that every collision between a metal atom and a grain results in a reaction. In that case, the reaction rate per unit volume for an atomic species i is given by ${{k}_{R}}=3.0\times {{10}^{-18}}{{T}^{1/2}}$ $/\left( 1.0+4\times {{10}^{-2}}{{T}^{1/2}} \right.$ $+\;\left. 2\times {{10}^{-3}}T+8.0\times {{10}^{-6}}{{T}^{2}} \right)$. The resulting fractions are shown in Table 5, denoted by asterisks because they are only approximate measurements. While the majority of hydrogen in our simulations remains atomic in all cases, these estimates indicate that in the highest Mach number runs, ≈1% of the gas could be in the form of H₂.

In addition to Table 5 we make available online the data files from each run.⁵ Each file presents a variety of information, including the true abundances and the simple abundance estimates presented in Section 4.3. We also give Doppler parameters for each species, i, defined as,

$$\begin{eqnarray}&&b_{i}^{2}=\sigma _{i,1{\rm D}}^{2}+2{{k}_{b}}{{T}_{{\rm MW}}}/{{A}_{i}}{{m}_{{\rm H}}},\end{eqnarray} \tag{ 17 }$$

where ${{\sigma }_{{\rm i},1{\rm D}}}$ is the 1D velocity dispersion, k_b is the Boltzmann's constant, and A_i is the ion atomic mass. While the masses of hydrogen and helium are small enough that the temperature term makes a substantial contribution, the Doppler parameters of the heavier elements are very close to σ_i,1D in general, providing a good measurement of the local velocity dispersion.