SIMULATIONS OF HIGH-VELOCITY CLOUDS. I. HYDRODYNAMICS AND HIGH-VELOCITY HIGH IONS

Kyujin Kwak; David B. Henley; Robin L. Shelton

doi:10.1088/0004-637X/739/1/30

1. INTRODUCTION

The halo of the Milky Way contains clouds of neutral hydrogen (H i) with |v_LSR| ≳ 90 km s⁻¹, known as high-velocity clouds (HVCs; Muller et al. 1963; Wakker & van Woerden 1997). Although some HVCs are relatively close to the Galactic disk (e.g., |z| ≲ 4 kpc for Complex M; Danly et al. 1993; Keenan et al. 1995), other HVCs are known to be in the upper Galactic halo, where the Galaxy interacts with its surroundings (e.g., d = 10 ± 2.5 kpc for Complex C; Thom et al. 2008). The distant HVCs are thought to be gas stripped off satellite galaxies (e.g., Gardiner & Noguchi 1996; Putman et al. 2004), extragalactic gas falling into the Galaxy (Oort 1966; Blitz et al. 1999), or gas left over from the formation of the Galaxy (Maller & Bullock 2004).

Observations of interstellar absorption of light from active galactic nuclei show that the halo also contains highly ionized gas (including high ions, such as C iv, N v, and O vi) moving with velocities comparable to those of the H i HVCs (Sembach et al. 2003; Fox et al. 2004, 2005, 2006; Collins et al. 2007). These high ions trace gas with T ∼ (1–3) × 10⁵ K, and enable us to probe the connection between the hottest Galactic gas (T ≳ 10⁶ K) and cooler phases of the interstellar medium (ISM). These ions also help us investigate the interaction between HVCs and the ambient halo gas, and, by extension, the properties of the clouds and the halo gas. Some highly ionized HVCs are seen in the same directions as H i HVCs, but there are also highly ionized HVCs without corresponding high-velocity H i (Sembach et al. 2003). The sky covering fraction of O vi HVCs (∼60%; Sembach et al. 2003) is larger than that of H i HVCs (∼37%; Lockman et al. 2002), suggesting the possibility that in the halo, warm-hot (a few times 10⁵ K) high-velocity gas traced by high ions could be more common than cool high-velocity gas traced by H i.

High ions in HVCs could hypothetically result from several physical processes, including mixing of cool neutral HVC gas with hot (T ≳ 10⁶ K) highly ionized ambient gas, thermal conduction between cool HVC gas and hot ambient gas, radiative cooling of hot high-velocity gas, and shock heating. For example, when a cool cloud is embedded in hot gas, thermal conduction can form a layer of intermediate-temperature gas (T ∼ a few times 10⁵ K) in which high ions are abundant (e.g., Borkowski et al. 1990, for stationary clouds). Also, simple radiative cooling of hot gas will result in it passing through the temperature range optimal for high ions. Collisional ionization equilibrium (CIE) and non-equilibrium ionization (NEI) calculations of this cooling, for stationary gas at least, have been carried out by many authors (see Gnat & Sternberg 2007, and references therein). High ions could also be produced in the hot gas behind the shocks that form as HVCs travel supersonically through the ambient medium. The observed ratios of high ions' column densities or line intensities can be compared with the results of these models to determine how high ions are produced; e.g., see Figure 7 in KS10. Although this particular comparison was made for low-velocity ions, one could test models for high-velocity ions in a similar way.

Some of these hypothetical processes require that hot gas is plentiful in the halo. Observations of the diffuse soft X-ray background (Burrows & Mendenhall 1991; Snowden et al. 1991, 1998; Kuntz & Snowden 2000; Smith et al. 2007; Galeazzi et al. 2007; Henley & Shelton 2008; Lei et al. 2009; Yoshino et al. 2009; Henley et al. 2010) and of X-ray absorption lines, such as O vii and O viii (e.g., Yao & Wang 2005, 2007; Fang et al. 2006; Bregman & Lloyd-Davies 2007; Yao et al. 2009) show that hot gas with T ∼ (1–3) × 10⁶ K is indeed present in the halo. Although the distance and origin of this hot gas are uncertain at present, it is likely that some of this hot gas exists in the same region as the H i HVCs. If so, then it is plausible that high ions are produced in the turbulent mixing layers (TMLs) between the cool HVCs and the hot ambient ISM. This turbulent mixing of cool and hot gas arises from the Kelvin–Helmholtz or shear instability induced by the velocity difference across the interface between the two types of gas. The TML model was first suggested by Begelman & Fabian (1990) and later developed by Slavin et al. (1993). Column densities of various high ions (including C iv, N v, and O vi) have been predicted for the TML model using a variety of techniques: analytical calculations (Slavin et al. 1993), three-dimensional (3D) magnetohydrodynamic (MHD) simulations (Esquivel et al. 2006), and two-dimensional (2D) hydrodynamic simulations incorporating NEI calculations (Kwak & Shelton 2010, hereafter KS10). Note that these studies all used a plane-parallel geometry, rather than a HVC-like geometry. TML column density predictions have been compared with observations of HVCs (Sembach et al. 2003; Fox et al. 2004; Collins et al. 2007; KS10), and should be of value for studying other astrophysical situations in which cool or warm material slides past hot material.

In this paper, we present the results of simulations that trace the interaction between an H i HVC and the hot ambient medium in a more realistic geometry than previous TML studies, that is, a spherical cloud falling through the hot ISM. As in our previous study (KS10), we trace the ionization states of carbon, nitrogen, and oxygen with non-equilibrium calculations which allow us to estimate the amounts of interesting ions (C iv, N v, and O vi) more accurately than CIE calculations. The new simulations enable us to investigate how the cold gas ablates from a spherically shaped cloud due to shear instabilities, how it mixes with the hot ambient gas, and where the high ions that are produced in the process of ablation and mixing reside. In addition, we can examine the velocities of the high ions: there are high-velocity high ions that move with similar velocities as the H i HVCs and low-velocity high ions that move with similar velocities as the ISM. Using the results of our simulations, we will answer the following questions: (1) How many high ions are produced by H i HVCs traveling through hot ambient ISM?, (2) How many high-velocity high ions result?, (3) How do the ratios of high ions and H i compare with observations?, and (4) Are HVCs likely to reach the Galactic disk intact? In this paper, we will concentrate on the high-velocity ions, deferring the discussion of low-velocity ions to Paper II (D. B. Henley et al. 2011, in preparation).

This paper is organized as follows. In the next section, we briefly summarize our numerical methods and present the parameters for our suite of seven simulational models. In Section 3, we present the results of our simulations. Specifically, Section 3.1 describes in detail the hydrodynamical evolution of one of our model clouds, Section 3.2 describes the evolution of the amount of neutral material, Section 3.3 describes the effect of the different model parameters on the cloud evolution, Section 3.4 discusses whether or not HVCs can reach the Galactic disk, and Section 3.5 describes the production of high ions in our simulations. In Section 4, we compare the column densities and column density ratios predicted by our simulations with observed Complex C values. In Section 5, we discuss the effects that the various assumptions in our simulations have on our results. We summarize our results in Section 6.

2. NUMERICAL METHODS AND MODEL PARAMETERS

We carried out our simulations using the same code as KS10, namely, FLASH version 2.5 (Fryxell et al. 2000). We used the FLASH NEI module to track the ionization evolution of carbon, nitrogen, and oxygen.¹ The simulations include radiative cooling, although the cooling curve was calculated assuming CIE. See KS10 for more details of the code, including some discussion of CIE versus NEI cooling rates. As in KS10, we used the abundances from Allen (1973), which are the default abundances in FLASH. Our results can be rescaled to give results for different abundance tables, assuming that the radiative cooling curve does not change significantly with the assumed abundances (an assumption which may not always be valid).

In KS10 the simulations were run in 2D Cartesian coordinates to study TMLs in a plane-parallel geometry. Here, we use 2D cylindrical coordinates to study initially spherical clouds as they fall through the hot ISM. As in KS10, we use 2D simulations in order to minimize the unaffordably large computing resources (particularly memory) needed to track the ionization of the three elements in our NEI simulations. In addition, instead of the cloud moving relative to a stationary ISM, in our simulations the ambient medium moves relative to an initially stationary cloud. In this way, we can trace the cloud evolution for a long time without requiring a large computational domain, reducing the amount of memory and CPU time needed to run the simulations. However, although the cloud is initially stationary in the computational domain, all velocities stated in this paper are in the observer's frame, in which the ISM is stationary. In the computational grid's reference frame, the ISM is initially moving in the +z direction with velocity |v_{z, cl}|, where −|v_{z, cl}| is the initial velocity of the cloud in the observer's frame. Therefore, the z-velocity of gas in a given grid cell in the observer's frame is v_z = v_{z, sim} − |v_{z, cl}|, where v_{z, sim} is the z-velocity of the gas in that cell obtained directly from the simulation. Note that, although we present results for HVCs with negative velocities (i.e., moving toward us), our results are equally applicable to HVCs with positive velocities if the signs on the velocities are changed.

The parameters of our various model clouds are presented in Table 1 and Figure 1. We ran seven different models in which we varied the cloud's initial radius (Table 1, Column 2), radial profile (Figure 1), z-velocity, v_{z, cl} (measured in the observer's frame; Table 1, Column 3), and number density expressed in terms of the hydrogen number density, where n_H/n_He = 10 (Table 1, Column 4). The ambient number density is 10³ times smaller than the initial cloud number density (Table 1, Column 7). All cloud material with a density ⩾5 times the ambient density initially moves at v_{z, cl} relative to the ambient medium, while the less dense outskirts of the cloud are initially set to the ISM's velocity. As in KS10, we required that the cloud and the ISM were initially in pressure balance. We achieved this by varying the temperature from T_cl = 10³ K at the cloud center to T_ISM = 10⁶ K in the ambient medium, such that the pressure was constant.

Table 1. Model Parameters

Model	Radius	v_{z, cl}	n_{H, cl}	M_init,T	M_init,v	n_{H, ISM}
	(pc)	(km s⁻¹)	(cm⁻³)	(M_☉)	(M_☉)	(cm⁻³)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
A	20	−100	0.1	120	130	1.0 × 10⁻⁴
B	150	−100	0.1	4.9 × 10⁴	5.1 × 10⁴	1.0 × 10⁻⁴
C	150	−150	0.1	4.9 × 10⁴	5.1 × 10⁴	1.0 × 10⁻⁴
D	150	−300	0.1	4.9 × 10⁴	5.1 × 10⁴	1.0 × 10⁻⁴
E	150	−150	0.1	4.9 × 10⁴	4.9 × 10⁴	1.0 × 10⁻⁴
F	300	−100	0.1	4.0 × 10⁵	4.2 × 10⁵	1.0 × 10⁻⁴
G	150	−100	0.01	4.9 × 10³	5.1 × 10³	1.0 × 10⁻⁵

Notes. Column 1: model identifiers. Column 2: approximate radius of the model cloud (except for Model E). See Figure 1 for the detailed density profile of each model cloud. Model E is a uniform-density cloud with an exact radius. Column 3: initial velocity of the cloud along the z-direction measured in the observer's frame. Column 4: initial hydrogen number density of the cloud at its center. Column 5: initial mass of cloud having a temperature T < 10⁴ K. Column 6: initial mass of cloud moving with v_{z, cl}. Note that all cloud material with a hydrogen number density greater than 5n_{H, ISM}, where n_{H, ISM} is the hydrogen number density of the ISM (Column 7), moves initially at speed v_{z, cl} relative to the ISM. All models have the same cloud temperatures (T_cl = 10³ K at cloud centers) and ISM temperatures (T_ISM = 10⁶ K).

Download table as: ASCII Typeset image

Because the clouds generally do not have sharp edges, the cloud mass is not well defined. For this reason, we have calculated two masses for each cloud: the mass of material with T < 10⁴ K (Table 1, Column 5), and the mass of material initially moving at v_{z, cl} relative to the ambient medium (Table 1, Column 6). These two masses are generally similar to each other. The cloud in Model E has a sharp edge, and so the two masses are identical in this model.

By running Models B, C, and D, we can see the effects of varying the cloud's initial velocity (Section 3.3.1), because these models otherwise have the same initial parameters. Apart from the radial density profile, Models C and E are identical, and so these two models show the effect of the cloud's initial density profile (Section 3.3.2). The cloud and ISM densities in Model G are an order of magnitude smaller than those in the other models, and so comparing Models B and G reveals the effect of varying the cloud's density (Section 3.3.2). Finally, Models A, B, and F reveal the effect of varying the cloud size (Section 3.3.3).

Although the number densities of HVCs and of the upper halo are not well constrained, our chosen model parameters are consistent with existing observational and theoretical constraints. The H i column densities on sightlines passing through the centers of our model clouds range from 9.3 × 10¹⁸ cm⁻² (Model G) to 1.9 × 10²⁰ cm⁻² (Model F), and are consistent with measured column densities of H i HVCs (e.g., Hulsbosch & Wakker 1988). The density of the hot gas in the upper halo is more uncertain, with observational estimates including n_H < 6.3 × 10⁻⁴ cm⁻³ from pulsar dispersion measures (Gaensler et al. 2008) and n_H = 7.5 × 10⁻⁴ cm⁻³ from O vii column density measurements (assuming a uniform spherical halo of radius 20 kpc; Bregman & Lloyd-Davies 2007).² The halo density has also been constrained by combining observations with theoretical models. Peek et al. (2007) obtained an upper halo density of n_H ∼ 2 × 10⁻⁴ cm⁻³ by considering the drag force on an HVC complex observed in H i. Grcevich & Putman (2009) estimated that the number density in the upper halo is >(2–3) × 10⁻⁴ cm⁻³, assuming that the Milky Way's satellite galaxies lost their gas through ram-pressure stripping as they passed through the upper halo. The halo density may not decrease dramatically with distance—Weiner & Williams (1996) estimated that the density near the Magellanic Stream (d ∼ 50 kpc) is n_H ∼ 1 × 10⁻⁴ cm⁻³, assuming that the observed Hα emission arises from an interaction between the Magellanic Stream and the ambient gas. Note that most of these estimated halo densities are slightly larger than the halo density used in most of our simulations (n_H = 1 × 10⁻⁴ cm⁻³ for Models A through F). However, Sembach et al. (2003) found that their observations of O vi HVCs favored a low-density extended hot halo (n_H ≲ 10⁻⁴ to 10⁻⁵ cm⁻³, r ≳ 70 kpc). The lower halo density in Model G (n_H = 10⁻⁵ cm⁻³) was chosen to investigate the lower range of halo densities, as well as a lower cloud density.

3. RESULTS

3.1. Hydrodynamical Evolution of High-velocity Clouds

We chose Model B as our representative model as it exhibits many of the physical processes also seen in the other models. In Section 3.3, below, we will discuss the differences and similarities between the various models. Here, we discuss Model B, whose evolution is shown in Figure 2. Figure 2 plots, from top to bottom, the hydrogen number density, the temperature, v_z (in the observer's frame), and the ion fractions³ of C iv, N v, and O vi. Note that at t = 0 Myr the ion fractions for all three high ions are low almost everywhere on the grid. This is because the cloud is too cold and the ambient medium too hot for these ions.

Figure 2 shows that the initially spherical cloud deforms during the simulation. This is due to the Bernoulli effect. As the ISM flows around the cloud, the speed of the ISM relative to the cloud is greater along the sides of the cloud (r ∼ 150, z ∼ 0 pc) than near the top or the bottom of the cloud (r ∼ 0, z ∼ ±150 pc). As a result, according to Bernoulli's equation, the pressure will be greater immediately in front of and behind the cloud than at the edge. Examining the pressures in the output from the hydrodynamical code confirms this expectation.

The cloud also deforms due to Kelvin–Helmholtz or shear instabilities (e.g., Chandrasekhar 1961), instigated by the velocity difference between the cloud and the ISM. Because the speed of the ISM relative to the cloud is largest at the edge of the cloud, the instabilities grow most rapidly here (Chandrasekhar 1961, Section 101), and this part of the cloud is pulled outward (see the hydrogen number density and v_z plots in Figure 2). As the ISM flows around the initially spherical cloud, a vortex develops behind the cloud. The flow in this vortex becomes more complicated as the shear instabilities start to ablate material from the edge of the cloud. This material impedes the flow of the ISM into the vortex behind the cloud, while new vortices form around the ablating material. This complicated flow causes ISM material some way behind the cloud to move in toward the r = 0 axis; this material then flows down (i.e., in the −z-direction) and then out again. As this material flows outward along the back of the cloud, it helps to stretch the cloud out in the horizontal direction.

The gas ablated from the cool, dense cloud (T_cl = 10³ K, n_{H, cl} = 0.1 cm⁻³) mixes with the hot, tenuous ISM (T_ISM = 10⁶ K, n_{H, ISM} = 10⁻⁴ cm⁻³), creating mixed gas of intermediate temperature. Eventually, the mixed gas reaches a temperature of a few times 10⁵ K, which is optimal for radiative cooling through line emission. As the mixed gas flows back from the edge of the cloud, the flow splits: some of the mixed gas is drawn into the vortex behind the cloud, while some flows further back from the cloud, continuing to mix with the hot ISM. The temperature of the mixed gas flowing away from the cloud continues to increase above a few times 10⁵ K, as the continued mixing raises the temperature more rapidly than the gas can radiatively cool. In contrast, the mixed gas flowing into the vortex behind the cloud cools, due to both mixing with cooler gas and radiative cooling. In particular, this gas cools more efficiently after it reaches the region near the r = 0 axis, resulting in the accumulation of gas with T ∼ 10⁴ K along this axis in the temperature plots of Figure 2. We traced the fractions of original cloud material and ISM material contained in the mixed gas and found that a significant fraction of the cool gas along the r = 0 axis at later times was initially hot ISM, indicating that this gas is mixed gas that has undergone radiative cooling.

The plots of the high ion fractions (the last three rows in Figure 2) show that the fractions of these high ions are higher in the mixed gas than in the initially cool cloud gas or in the initially hot ISM gas. As the cool and hot gas mix, the ions that we are interested in are produced both by ionization during the heating of the initially cool gas, and by recombination during the cooling of the initially hot gas. The fraction of Li-like ions for each element depends somewhat on the temperature of the gas; for example, in the hotter gas the O vi fraction is higher than the C iv fraction. The NEI ionization levels in Figure 2 are different from those expected from CIE, because changes in the ionization levels lag behind the changes in temperature that are due to mixing or radiative cooling. This is similar to what we found in our previous study of TMLs (KS10).

It is possible that our use of 2D cylindrical geometry resulted in an overestimate of the amount of cool gas that accumulates along the r = 0 axis. In such a geometry, material that flows toward the r = 0 axis tends to stick to the axis, because of the reflecting boundary condition there. The 2D cylindrical geometry is also responsible for the protuberance at the front of the cloud seen at later times, for example, the feature indicated by the arrow in the far right density plot in our Figure 2 (Vieser & Hensler 2007b). However, although the high ions contained in the cooled gas that accumulates along the symmetry axis give rise to large ion column densities along this sightline (see Section 3.5), this column of material does not contribute much to the total number of high ions on the grid.

3.2. Evolution of the Number of H i Atoms

The default FLASH NEI module assumes that hydrogen is fully ionized, and so we were unable to trace the ionization evolution of hydrogen. We therefore assumed that the hydrogen in gas with T < 10⁴ K is entirely neutral (H i), while the hydrogen in hotter gas is fully ionized. Our initial clouds "lose" their H i content via the physical processes discussed above, namely, ablation of material from the cloud and its subsequent temperature increase due to mixing with the hot ISM. (Note, however, that H i can also be replenished when hot gas cools below 10⁴ K, and such gains can count against the losses.) Here, we consider the loss of H i from the clouds due only to heating, and the loss of H i due to heating and/or ablation. In the former case, we consider all H i atoms, regardless of their velocity. In the latter case, we consider only H i atoms moving with HVC-like velocities: v_z ⩽ −80 km s⁻¹ for Models A, B, F, and G (in which the initial velocity of the cloud was v_{z, cl} = −100 km s⁻¹), and v_z ⩽ −100 km s⁻¹ for Models C, D, and E (in which v_{z, cl} was −150 or −300 km s⁻¹). These definitions were chosen to correspond with observational analyses, in which one can distinguish between high- and low-velocity H i.

To investigate the loss of H i from the clouds, for each model we calculate the ratio, β(t), of the number of neutral hydrogen atoms lost since the beginning of the simulations to the initial number of neutral hydrogen atoms, i.e.,

$\begin{equation} \beta (t) \equiv \frac{\mathcal {N}_\mathrm{H\,\mathsc{i},init}- \mathcal {N}_\mathrm{H\,\mathsc{i}}(t)}{\mathcal {N}_\mathrm{H\,\mathsc{i},init}}, \end{equation} \tag{ 1 }$

where $\mathcal {N}_\mathrm{H\,\mathsc{i},init}$ and $\mathcal {N}_\mathrm{H\,\mathsc{i}}(t)$ are the initial and current numbers of H i atoms, respectively. As noted above, we have two different ways of counting $\mathcal {N}_\mathrm{H\,\mathsc{i}}(t)$ , leading to two different values of β(t): β_All(t) for H i at all velocities and β_HVC(t) for HVC-like H i. Note that the initial number of H i atoms, $\mathcal {N}_\mathrm{H\,\mathsc{i},init}$ , is the same in both cases, as all the H i is high velocity at the start of the simulation. The cool gas that accumulates along the r = 0 axis after ablation, mixing, and cooling has low velocities; it is thus included in β_All(t) but not β_HVC(t). We did not include material that escapes from the top of the computational domain. This means that the true value of $\mathcal {N}_\mathrm{H\,\mathsc{i}}(t)$ is larger than the value obtained from the computational domain, and so our estimates of β(t) are upper limits. However, the material that escapes from the top of the domain has low densities and low velocities, and so the amounts of HVC-like material are not significantly affected.

Figure 3 shows β(t) for each of our seven models. The black lines show β_All(t), and the colored lines show β_HVC(t). In all cases, β generally increases with time, indicating that H i is lost throughout the simulation. For Models A through E plus G, the black and colored lines are similar to each other (they are almost identical for Model A), indicating that most of the ablated material is "hot" (i.e., above 10⁴ K). However, for Model F the colored line is clearly above the black line, indicating that some of the ablated material is not yet ionized, or that some of the ablated gas has been heated and subsequently cooled. We find that the latter explanation is the more important: the amount of radiatively cooled ablated gas that accumulates along the r = 0 axis is larger in Model F than in the other models.

Here, we note a number of features from Figure 3. (1) β is relatively insensitive to the cloud's initial velocity, over a wide range of velocities (100–300 km s⁻¹; compare Models B, C, and D). (2) β is also relatively insensitive to the cloud's initial density profile (compare Models C and E) and the cloud and ISM's initial densities (compare Models B and G). (3) A smaller cloud loses its H i content, as a fraction of its initial mass, faster than a larger cloud (compare Models A, B, and F). With these various trends in mind, we will discuss the effects of the different model parameters in more detail in the following section.

3.3. Differences between the Models—The Effects of Different Model Parameters

Figure 4 compares the evolution of Models A, C, D, E, F, and G by plotting the temperature on a logarithmic scale at various times (the evolution of Model B is shown in Figure 2). Note that for Models C through G, the four panels correspond to similar stages in these clouds' evolution, while the first Model A panel corresponds to a similar stage of evolution as the final panels for the other models (see Section 3.3.3). The computational domain for Model D has a larger height than those for Models B, C, E, and G, because of the higher initial velocity of the cloud. In simulations with higher initial cloud velocities (Models C, D, and E), the location of the cloud shifts upward further than in the simulations with lower cloud velocities (Models B, F, and G). This upward shift is because, in our simulations, the cloud is initially stationary while the ISM flows upward, pushing the cloud upward (see Section 2). However, despite the low velocity in Model A, the cloud in this simulation is still shifted upward significantly, due to its low inertia relative to the ISM.

3.3.1. The Effect of the Cloud Velocity (Models B, C, and D)

The initial conditions for the cloud and the ISM are the same in Models B, C, and D, apart from the initial velocity of the cloud: −100, −150, and −300 km s⁻¹ in the observer's frame, respectively. Because the sound speed of the ISM (T = 10⁶ K) is ∼150 km s⁻¹, these velocities correspond to the subsonic, transonic, and supersonic regimes, respectively. We would therefore expect a bow shock to develop in Model D (supersonic case); Figure 4 shows that a bow shock does indeed develop at early times in this model, and persists until the end of the simulation.

Apart from the formation of a bow shock in Model D, Models C and D both evolve with similar hydrodynamical processes that were seen in Model B: Bernoulli's effect, ablation of the cloud's material due to shear instabilities, mixing of the ablated gas with the hot ISM, and cooling of the mixed gas. However, the shear instabilities grow more rapidly the larger the velocity difference between the cloud and the ISM (Chandrasekhar 1961, Section 101). The faster-growing instabilities clearly affect the evolution of the cloud: a faster cloud disrupts more violently than a slower cloud. The temperature plots show that the amount of cool cloud material (shown in white) is smaller for a faster cloud at a given time: compare the third, fifth, seventh, and ninth temperature panels in Figure 2 (Model B) with the Model C and D panels in Figure 4 (these sets of panels correspond to the same times: t = 30, 60, 90, and 120 Myr). These differences in the severity of the cloud disruption due to shear instabilities lead to different morphologies for the clouds with different velocities.

Although different cloud velocities lead to different cloud morphologies, we find that the faster-growing instabilities do not necessarily lead to material being ablated and/or ionized at larger rates for the faster clouds. Figure 3 shows that the rates at which the Model B, C, and D clouds lose their H i to ablation and/or ionization, offset by material that has cooled, are rather similar. Note that in Model D, the bow shock that forms in front of the cloud helps protect the cloud from ablation, because the shocked ISM has a lower velocity than the unshocked ISM. Note also that in Model D, the mixing of the ablated gas with the hot ISM is also different from the mixing in Models B and C: the fast-moving ISM constrains the ablated material so that it remains close to the cloud. This gas mixes and cools as it moves along the cloud periphery. The mixed gas is so closely constrained to the edge of the cloud that it cannot be clearly seen in the temperature plots (Figure 4), but high ions are abundant along the periphery of the cloud where the mixed gas exists (see Figure 5).

**Figure 5.** O vi ion fraction for Model D at t = 30 Myr (the same time as the first Model D panel in Figure 4).
Download figure:
Standard image High-resolution image

3.3.2. The Effect of the Cloud Density Profile (Models C and E) and the Cloud Density (Models B and G)

Models C and E have the same initial conditions, apart from the initial density profile of the cloud: in Model C, the cloud density decreases smoothly at the edge of the cloud until it equals the ambient density, whereas in Model E the cloud has a uniform density and a sharp edge (see Figure 1). Models E and C allow us to compare the results for a uniform, sharp-edged density profile with those for a more realistic density profile. Note that similar comparisons have been made in previous works (Vieser & Hensler 2007b; Heitsch & Putman 2009).

The first two panels for Model C in Figure 4 show that a tail of low-temperature gas propagates upward immediately behind the cloud, whereas Model E has not developed such a tail at early times. These tails are composed of the low-density gas that was initially at the edge of cloud in Model C, due to the smooth transition between the cloud density and the ISM density at the cloud edge. Similar tails are also seen in the other models at early times, except in Model E. The low-density gas at the edge of the Model C cloud is ablated via shear instabilities more quickly than the high-density gas at the edge of the Model E cloud. This ablated material accumulates near the r = 0 axis at early times after mixing with the ISM and radiatively cooling, and moves upward along this axis, forming the tail. Cool gas does eventually start accumulating along the r = 0 axis in Model E, but at later times (see the third Model E panel in Figure 4, corresponding to t = 90 Myr).

Models B and G have the same initial conditions, except for the cloud and ISM densities being an order of magnitude lower in Model G. The Model G cloud generally evolves in a similar fashion to the Model B cloud; this is because the ratio of densities across the cloud–ISM interface are the same in both models, and so the growth of perturbations due to the shear instability should occur at the same rate (Chandrasekhar 1961, Section 101). However, the lower density in Model G reduces the cooling rates, so the mixed gas cools more slowly. As a result, very little of the ablated gas that accumulates along the r = 0 axis in Model G is cool.

As mentioned at the end of Section 3.2, the detailed differences in the hydrodynamical evolution of the cloud between Models C and E as well as between Models B and G do not affect β significantly (Figure 3), indicating that the rate at which the cloud loses its H i does not strongly depend on either the initial density profile or the initial density of the cloud, provided that the density contrast between the cloud and the ISM is the same.

3.3.3. The Effect of the Cloud Size (Models A, B, and F)

Models A, B, and F have the same initial conditions, apart from the size of the cloud; in these three models, the radii of the initial clouds are ≈20, ≈150, and ≈300 pc, respectively (see Figure 1). As these clouds all have the same initial velocity (−100 km s⁻¹), they should all be subject to the same physical processes.

As noted at the end of Section 3.2, a smaller cloud loses its H i content more rapidly than a larger cloud. In order to better understand this trend, we consider the simple case of a uniform spherical cloud losing mass at a rate proportional to its surface area, i.e., dM(t)/dt∝ − 4πr²(t). In this case, the radius decreases linearly with time, r(t) = r₀ − kt, where r₀ is the initial radius of the cloud and k is the rate at which the radius decreases. Because the cloud in this simple model is uniform and spherical, $\mathcal {N}_\mathrm{H\,\mathsc{i}}(t) \propto r^3(t)$ , and so β(t) = 1 − (1 − kt/r₀)³ (from Equation (1)).

The smooth thick light gray curves in Figure 3 show β(t) calculated according to this simple model with k = 0.1 pc Myr⁻¹ and r₀ = 20, 150, and 300 pc (top to bottom). This simple model represents the loss of H i from the cloud due to all processes, and so should be compared with the colored (HVC-like H i) curves; as previously noted, these curves represent H i lost to ablation and to ionization. Note that this model is not an accurate physical model—Figures 2 and 4 show that the cloud does not remain spherical during its evolution. As a result, we did not fit the thick light gray curves in Figure 3 to the corresponding curves derived from the hydrodynamical simulations. Nevertheless, this simple model can provide some insight into the relative behavior of Models A, B, and F.

The curve derived from this simple model for r₀ = 20 pc is in reasonably good agreement with the Model A curve. For Models B and F, the simple model overestimates the mass loss at earlier times, and for Model F it slightly underestimates the mass loss at later times. Despite these shortcomings, this simple model indicates that a major reason that a smaller cloud loses its mass more rapidly is because it has a larger surface area relative to its mass.

If the clouds evolve according to this simple model, with the same value of k for all clouds, then the Model A cloud will be at a similar phase in its evolution (i.e., at the same value of r(t)/r₀) at t = 16 Myr as the Model B cloud at t = 120 Myr and as the Model F cloud at t = 240 Myr. We find that the morphology of the cloud in all three models varies in a similar fashion according to this timescale, i.e., the Model A, B, and F clouds have similar shapes at t = 16, t = 120, and t = 240 Myr, respectively. To see this, compare the first Model A panel in Figure 4 (t = 15 Myr), the final temperature panel in Figure 2 (Model B at t = 120 Myr), and the final Model F panel in Figure 4 (t = 240 Myr). Therefore, the Model A cloud after t = 16 Myr (second through fourth Model A panels in Figure 4) is at a later phase in its evolution than the Model B cloud at t = 120 Myr and the Model F cloud at t = 240 Myr.

At a given time t, although a smaller cloud will have lost more of its H i mass relative to its initial mass than a larger cloud, the above simple model predicts that the larger cloud will have lost more mass. Our simulations bear out this prediction—the total mass lost from the Model F cloud is larger than the mass lost from the other, smaller model clouds. Furthermore, if the ablated material were to mix and cool at the same rate in all models, then the fraction of the ablated gas that is cool would be the same in all models, regardless of the cloud's initial size. This would result in the ratio of β_All to β_HVC being the same in all models. However, Model F, which has the largest cloud, yields larger ratios than the other models, i.e., the difference between the black and colored lines for Model F in Figure 3 is larger than for the other models. This probably indicates that more radiatively cooling takes place in the ablated material in Model F than in the models with smaller clouds.

3.4. The Fate of High-velocity Clouds

Heitsch & Putman (2009) modeled HVCs of various sizes (initial H i mass, $M_\mathrm{H\,\mathsc{i},init}= 10^{3.1}$ –10^4.6 M_☉), and found that HVCs with $M_\mathrm{H\,\mathsc{i},init}< 10^{4.5} \;M_{\odot }$ will lose all their H i after traveling <10 kpc. Thus, smaller HVCs are unlikely to reach the disk as neutral hydrogen; nor will larger HVC complexes, if they are in fact composed of small cloudlets. However, it is possible that the HVC material could still reach the disk in the form of warm ionized material (Heitsch & Putman 2009; Shull et al. 2009; Bland-Hawthorn 2009).

Our suite of models includes one cloud (Model F) with $M_\mathrm{H\,\mathsc{i},init}= 10^{5.5} \;M_{\odot }$ ,⁴ which is larger than the masses of the clouds simulated by Heitsch & Putman (2009). The Model F cloud has β_HVC ≈ 0.3 at the end of the simulation (t = 240 Myr), i.e., ≈70% of the cloud's initial high-velocity H i mass remains at T < 10⁴ K at this time, although it is possible that some of this high-velocity material has broken off the main cloud. The cloud will travel much further before completely dissipating—the value of β_HVC from the end of the Model A simulation and the scaling discussed in Section 3.3.3 imply that ∼30% of the Model F cloud's initial high-velocity H i will remain below 10⁴ K at t ∼ 900 Myr.

As well as our simulating a more massive cloud, there are other differences between our and Heitsch & Putman's simulations: we used a 2D geometry and a cooling curve calculated with solar abundances, whereas Heitsch & Putman carried out 3D simulations with 1/10 solar abundances. Both of these differences would tend to stabilize our clouds against disruption (Section 5), leading to longer-lived clouds in our simulations. However, comparing similar-sized clouds in our and Heitsch & Putman's simulations indicates that the cloud lifetimes agree within a factor of ∼3–5. Therefore, even taking this into account, our Model F simulations indicate that very large clouds (≳10^5.5 M_☉) will live for at least a few hundred megayears, and travel a few tens of kiloparsecs (assuming a speed of ∼100 km s⁻¹). Such large cloud masses are not implausible: e.g., masses of ≳2 × 10⁶ M_☉, ∼10⁷ M_☉, and ∼10⁷ M_☉ has been measured for the Smith Cloud (Nichols & Bland-Hawthorn 2009), Complex C (Wakker et al. 2007; Thom et al. 2008), and Complex H (Lockman 2003), respectively (although the complexes are not single clouds). The largest HVCs may therefore survive as far as the Galactic disk at least partially as neutral hydrogen.

3.5. High Ions

In our simulations, the gas that ablates from the cold cloud and mixes with the hot ISM is rich in high ions (C iv, N v, and O vi). In this section, we investigate the properties of these high ions in more detail. First, we estimate the quantities of high ions that are produced via ablation and ionization of the cloud material, and then we discuss the variation of column density with radius for each high ion; these column densities can be directly compared with observations.

Figure 6 shows the masses of the high ions in our simulational domains as functions of time. In each case, we calculate two values: the total mass of a given ion in the simulational domain, regardless of velocity (black lines), and the mass of a given ion moving with HVC-like velocities (colored lines). When calculating the total mass of a given ion, integrated over all velocities, we subtract off the mass of that ion that was contained in the hot ambient ISM at t = 0. Our values do not include the ions that have escaped from the domain. This means that the total masses of the high ions are actually lower limits. However, as with the amount of HVC-like H i (Section 3.2), the masses of the HVC-like ions are not significantly affected by our neglecting material that has flowed off the domain, because such material has a low density and a low velocity in the observer's frame. Note from Figure 6 that the masses of the high ions with HVC-like velocities are smaller than the corresponding total masses integrated over all velocities. This is because the gas that ablates from the cloud decelerates as it mixes with the ambient gas. Eventually, the mixed gas slows to halo-like velocities, but still contains high ions.

In Models B through G, the masses of the high ions generally increase with time, although there are fluctuations on small timescales. Because the high ions are abundant in the material ablated from the cloud, the time evolution of their masses is similar to that of the lost H i in two distinct ways. First, similar amounts of high ions are produced regardless of the cloud's initial velocity (compare Models B, C, and D; cf. Section 3.3.1) or density profile (compare Models C and E; cf. Section 3.3.2). Typically, the amounts of high ions produced in Models B through E agree within a factor of two. Second, more high ions are produced from the larger cloud than from the smaller cloud (compare Models A, B, and F; cf. Section 3.3.3). The density of the Model G cloud is 1/10 that of the Model B cloud, and so the masses of the high ions produced in Model G are commensurately lower. Unlike the other models, Model A shows a decrease in the mass of high ions at later times. This occurs because Model A is sampling a much later phase of the cloud evolution (see Section 3.3.3), such that by the end of the Model A simulation, the cloud is mostly destroyed and so provides fresh cool gas for mixing with the ISM more slowly.

To investigate the spatial distributions of the high ions, we calculated column densities for sightlines running vertically through the Model B domain. Figure 7 shows H i, C iv, N v, and O vi column densities as functions of radius for t = 30, 60, 90, and 120 Myr. The footprints of the high ion distributions increase with time, especially for the ions with slow velocities. The column densities of the high ions generally increase with time while the H i column density decreases, because the cool H i-rich gas is ionized and converted to high ion-bearing gas. Most of the material containing high ions with HVC-like velocities was recently ablated from the cloud and ionized to this level, which is why the footprints of the HVC-like high ions (gray lines in Figure 7) are similar to the footprint of the H i cloud.

**Figure 7.** Column density as a function of radius for H i (top left), C iv (top right), N v (bottom left), and O vi (bottom right), from Model B. In each plot, the panels correspond to t = 30, 60, 90, and 120 Myr (top to bottom). The black lines are the column densities integrated over all velocities, and the gray lines for the column densities for material with HVC-like velocities. Note that for H i (top left), these two lines are almost identical.
Download figure:
Standard image High-resolution image

The column densities of H i calculated throughout the Model B simulation are above the current 21 cm detection limit (a few times 10¹⁸ cm⁻²) over most of the cloud, so if the Model B cloud were a real cloud, it could be identified as an HVC throughout the whole simulation period. We find that simulations with similar clouds (i.e., Models C, D, and E) have similar H i column densities to Model B. The H i column densities for larger (Model F), smaller (Model A), or less dense (Model G) clouds than the Model B cloud vary according to their size or density. Even the smaller column densities (from Models A and G) are above the current detection limit over much of the clouds up to the ends of the simulations.

At the edges of the clouds, however, the H i can become undetectable, while the high ions remain detectable. For example, in Model B at t = 120 Myr, the H i column density falls below the 21 cm detection limit beyond r ≈ 300 pc, but there are substantial amounts of O vi (up to ∼10^13.5 cm⁻²) out to larger radii (see Figure 7). Also, it should be noted that the column densities plotted in Figure 7 are for vertical sightlines through our model domains. The high-velocity high ions reside mainly behind the main body of the cloud, in material that has ablated from the cloud, mixed with the ambient medium, and fallen behind the cloud. Therefore, diagonal sightlines through the model domain could intersect significant column densities of high ions, while missing most of the H i. Our simulation results could therefore partially explain the presence of highly ionized high-velocity gas on sightlines without corresponding high-velocity H i (Sembach et al. 2003). Our results offer only a partial explanation because they account only for sightlines with high-velocity high ions that have high-velocity H i nearby; some of the high-velocity O vi detections in Sembach et al. (2003) are on sightlines without any high-velocity H i within several degrees.

We noted above that similar amounts of high ions are produced in Models C, D, and E as in Model B (Figure 6). The top four panels in Figure 8 show that the O vi column densities from Models C, D, and E are similar to those from Model B at t = 120 Myr, although the footprints vary from model to model. In contrast, the Model G O vi column densities (fifth panel in Figure 8) are about an order of magnitude smaller than those from Model B, which is commensurate with the density difference between these two models.

As was noted in Section 3.3.3, Model A is mostly in a later evolutionary phase than the other models. From Figure 9, we see that the O vi column densities in this model peak around t = 30 Myr and subsequently decrease (see also the bottom panel in Figure 6). When we compare Models A, B, and F at similar stages of their evolution (t = 15, 120, and 240 Myr, respectively, see Section 3.3.3), we find that the O vi column densities with HVC-like velocities from Models A and F are factors of ∼5 smaller and ∼1.5 larger than those from Model B, respectively (∼2 × 10¹³, ∼1 × 10¹⁴, and ∼1.5 × 10¹⁴ cm⁻² for Models A, B, and F, respectively). The ratio of the O vi column densities from Models A, B, and F (approximately 0.2:1:1.5) is similar to but not equal to the ratio of initial radii, 0.13:1:2.

In the following section, we compare our column density predictions with observations of Complex C. In particular, we will consider the ratios of the high ions to H i and to each other. Such ratios predicted by different models are frequently compared with the observed ratios in order to identify the physical process(es) by which the high ions are produced.

4. COMPARISON WITH OBSERVATIONS OF COMPLEX C

Among the several HVC complexes, Complex C has been studied most extensively, perhaps because of its location in the northern hemisphere and its large size on the sky (∼1800 deg²; Wakker & van Woerden 1991). Its distance is reasonably well constrained (10 ± 2.5 kpc; Thom et al. 2008), and it is of rather low metallicity (∼0.13 solar; Collins et al. 2007). Observations of high-velocity high ions along several sightlines passing through the Complex C region have been reported by Sembach et al. (2003), Fox et al. (2004), and Collins et al. (2007). In addition, Complex C is thought to be colliding with the Milky Way's thick disk (Tripp et al. 2003). In this section, we will compare our model predictions with observations of Complex C in several different ways. First, in Section 4.1, we compare the ion column densities, N, predicted by our models with those measured for Complex C (Sembach et al. 2003; Fox et al. 2004; Collins et al. 2007). Sembach et al. (2003) looked for but did not find a correlation between the Complex C O vi and H i column densities. We discuss this lack of correlation in terms of our model predictions in Section 4.2. Finally, in Section 4.3 we calculate the ratios of column densities averaged over the whole of the model cloud, and we compare these ratios with observations in Sections 4.4 and 4.5.

4.1. Ion Column Densities

Sembach et al. (2003) searched for high-velocity O vi in the survey of high-latitude FUSE observations (Wakker et al. 2003). They detected a total of 84 high-velocity O vi absorption features on 59 out of the 102 sightlines in the survey, spread over the high-latitude sky. The median O vi column density of these 84 detected features is 9.3 × 10¹³ cm⁻², although there is significant sightline-to-sightline variation (between ∼10¹³ and ∼3 × 10¹⁴ cm⁻²; Sembach et al. 2003). For the nine sightlines toward Complex C with detections, the median high-velocity O vi column density is 7.6 × 10¹³ cm⁻² (range: (4.7–16.6) × 10¹³ cm⁻²; Sembach et al. 2003). In Model B, the column densities of HVC-like O vi approach the median of the observed column densities at some times for some impact parameters, but are generally lower than the observed median. Furthermore, along some sightlines toward Complex C, there are measurements of C iv and N v column densities (Fox et al. 2004; Collins et al. 2007). The observed high-velocity C iv column density ranges from <6.9 × 10¹² to 6.5 × 10¹³ cm⁻² while that of N v ranges from <6.9 × 10¹² to <3.1 × 10¹³ cm⁻² (Collins et al. 2007). For certain choices of age and impact parameter, our Model B results match the column densities at the higher end of the observed ranges (e.g., r ∼ 150 pc at t = 120 Myr for C iv and r ∼ 0 pc at t = 120 Myr for N v) and at the lower end (e.g., r ∼ 150 pc at t = 60 Myr for C iv and r ≳ 100 pc at t = 60 Myr for N v).

4.2. N(O vi) versus N(H i)

The filled stars in Figure 10 show $N(\mbox{O\,{\sc {vi}}})$ against $N(\mbox{H\,{\sc {i}}})$ for nine sightlines through Complex C (Sembach et al. 2003, Table 8; we followed Sembach et al. and did not include the H i components in brackets in their table). As noted by Sembach et al., there is no correlation between $N(\mbox{O\,{\sc {vi}}})$ and $N(\mbox{H\,{\sc {i}}})$ ; $N(\mbox{O\,{\sc {vi}}})$ varies by ≈0.5 dex over Complex C while $N(\mbox{H\,{\sc {i}}})$ varies by ≈1.5 dex. The other symbols in Figure 10 show the predictions from our various models, extracted for multiple vertical sightlines through our model domains at several different times. For any given model at a given time, the predictions for different sightlines exhibit a similar trend to the observations, i.e., $N(\mbox{O\,{\sc {vi}}})$ is approximately constant over a large range of $N(\mbox{H\,{\sc {i}}})$ . The trend in the predictions is because, in our simulations, the H i column density decreases significantly toward larger radii, whereas the O vi column densities remain more or less constant out to larger radii (see Figure 7). Not only do the predictions exhibit the same trend as the observations, but also some of the predicted values are in good agreement with the observations. However, not all of the observations can be matched by a single model at a specific time. For example, five of the nine observed data points (the four with the largest values of $N(\mbox{O\,{\sc {vi}}})$ and the one with $N(\mbox{H\,{\sc {i}}})\approx 10^{19.0} \;\mbox{cm}^{-2}$ , $N(\mbox{O\,{\sc {vi}}})\approx 10^{13.8} \;\mbox{cm}^{-2}$ ) are similar to values from Model C at t = 120 Myr (cyan squares) and from Model F at t = 240 Myr (orange diagonal crosses). The remaining data points are similar to predicted values at other times and/or from other models.

The predicted values from Models A and G (small and low-density clouds, respectively) have H i and O vi column densities that are smaller than the observed values for Complex C, regardless of simulation time and sightline. However, it is possible that Complex C is composed of many small or low-density clouds, each of which individually resembles a Model A or G cloud. If this is the case, the total column densities would be obtained by multiplying the Model A or G predictions (which are for a single cloud) by the number of clouds along the line of sight, resulting in a shift toward larger values of $N(\mbox{H\,{\sc {i}}})$ and $N(\mbox{O\,{\sc {vi}}})$ . This means that the predictions for multiple Model A or G clouds would belong upward and to the right of the Model A or G points in Figure 10, and therefore could become consistent with the observations. Note that, in order to carry out this shift, we must assume that each cloud contributes the same amount of O vi and H i at the same period in time.

The consistency between the observed and predicted values of $N(\mbox{O\,{\sc {vi}}})$ against $N(\mbox{H\,{\sc {i}}})$ supports the idea that Complex C has interacted with the hot ISM to produce O vi. In particular, this consistency confirms that the loss of cloud material due to shear instabilities and its subsequent mixing with the hot ISM, subject to radiative cooling, is a viable physical process for the production of O vi. However, the timescale for Complex C's interaction with the hot ISM is uncertain, as it is dependent on the size and density of the constituent clouds—a set of smaller clouds could produce similar values of $N(\mbox{O\,{\sc {vi}}})$ and $N(\mbox{H\,{\sc {i}}})$ to a single larger cloud that has been evolving for a longer time.

In the following subsections, we will look at the evolution of the column density ratios for each model, and compare them with the Complex C observations. This will enable us to place a rough constraint on the timescale for Complex C's interaction with the hot ISM.

4.3. Column Density Ratios

For comparison with the observed column density ratios, we calculated ratios of column densities that had been averaged over the model clouds. We considered only sightlines for which the column density was above some cutoff, N_cut, so that material not significantly affected by the cloud–ISM interaction is not included in the averages. For a given ion, the average column density, $\bar{N}(t)$ , at time t is given by

$\begin{equation} \bar{N}(t) = \frac{ \int N(t,r) \alpha (t,r) r\,dr }{\int \alpha (t,r) r\,dr }, \end{equation} \tag{ 2 }$

where N(t, r) is the column density for a sightline at radius r, and

$\begin{equation} \alpha = \left\lbrace \begin{array}{@{}rl@{}}1 & \mbox{if $N > N_\mathrm{cut}$;} \\ \ms 0 & \mbox{otherwise.} \end{array} \right. \end{equation} \tag{ 3 }$

For H i, N_cut = 10¹⁶ cm⁻² for Models B through F and 10¹⁵ cm⁻² for Models A and G. For the high ions, N_cut = 10¹¹ cm⁻² for Models B through F and 10¹⁰ cm⁻² for Models A and G. We found that the choices of N_cut did not significantly affect the average column densities and the resulting ratios.⁵ We did not subtract off a background column density due to the hot ambient medium: because the high ion fractions are small in this medium, the background column density is negligible.

The evolution of the column density ratios calculated above is presented in Figure 11 (ion to H i) and Figure 12 (ion to ion). Because we will compare the models with Complex C observations, we plot only the ratios for HVC-like material. These curves were calculated using solar abundances (Allen 1973), but subsolar abundances have been reported for Complex C (Fox et al. 2004; Collins et al. 2007). We can adjust our predicted ratios by shifting the curves in Figures 11 and 12 according to the ratios of the Complex C abundances to solar abundances. The black vertical lines in Figures 11 and 12 indicate the amounts by which the curves should be shifted downward to go from solar abundances to the abundances measured along the PG 1259+593 sightline (Fox et al. 2004; see Table 8 in KS10 for the relevant multiplication factors). However, it is important to note that the abundances measured along the PG 1259+593 sightline are the lowest among the Complex C measurements (Collins et al. 2007), and so the shifts indicated in Figures 11 and 12 would be the most extreme cases. It is also important to note that the Complex C abundance measurements use low ions, and thus may only be sampling cloud material, rather than the turbulently mixed material in which the high ions are produced. The turbulently mixed material includes a contribution from the ambient medium, which may have a different metallicity from the cloud. Note that nitrogen is more depleted along the PG 1259+593 sightline than carbon and oxygen.

Figures 11 and 12 also show the observed ratios for Complex C from Fox et al. (2004) and Collins et al. (2007). In some cases, these two studies give slightly different column densities for the same sightline, due to differences in their measurement methods (e.g., the determination of the continuum and the velocity range used for the column density integration). In such cases, we plot the measurements from both studies. In the following two subsections, we discuss the evolution of the ion-to-H i and ion-to-ion ratios that we have plotted, and compare these ratios with the Complex C observations.

4.4. Ion-to-H i Column Density Ratios

The predicted ratios of high ions to H i generally increase over time, although large fluctuations occur in some models at earlier times (Figure 11). The ratios are largest for the smallest cloud (Model A) and smallest for the largest cloud (Model F) at a given time, and so an observed ratio can correspond to a smaller cloud that has evolved for a short period of time, or a larger cloud that has evolved for a longer period of time (cf. Section 4.2). Either interpretation would be possible without knowing Complex C's 3D structure (i.e., whether Complex C is composed of many small clumps, fewer larger clumps, or a few large clumps mingled with many small clumps).

The observed ratios vary by more than an order of magnitude from sightline to sightline. The variations in the predicted ratios, due to differences in the age or size of a cloud, are as large as the observed variations. The range of ratios predicted with solar abundances covers the observed range. However, when we shift the predicted ratios according to the subsolar Complex C abundances, the models can only match the lower observed ratios.

Although there may be a few ways to interpret the observed ratios in terms of our models, we provide the following interpretation as a possible example. The observed $N(\mbox{N\,{\sc {v}}})/N(\mbox{H\,{\sc {i}}})$ ratios are almost all upper limits, and do not help to constrain the models. The $N(\mbox{C\,{\sc {iv}}})/N(\mbox{H\,{\sc {i}}})$ and $N(\mbox{O\,{\sc {vi}}})/N(\mbox{H\,{\sc {i}}})$ ratios observed toward PG 1259+593 are similar to those predicted by Model A at t ∼ few × 10 Myr or by Models B through E at t ∼ 100 Myr, after the abundance shifts indicated by the vertical bars in Figure 11 have been applied. Nearly all the other sightlines have larger observed $N(\mbox{C\,{\sc {iv}}})/N(\mbox{H\,{\sc {i}}})$ and $N(\mbox{O\,{\sc {vi}}})/N(\mbox{H\,{\sc {i}}})$ ratios than the PG 1259+593 sightline. As the predicted ratios rise with time, the model predictions (after the abundance shift) will match these other sightlines at even later times. As noted above, the abundance shifts indicated in Figure 11 are the most extreme cases, as the PG 1259+593 sightline has the lowest abundances of the Complex C sightlines. However, even the solar-abundance predictions from some models (e.g., B through E) do not match the ratios observed along some sightlines until t ∼ 100 Myr. These results suggest that the models require at least ∼100 Myr of evolution to match the observed ratios, if the radii of the composite clouds were initially ≳150 pc (Models B through E). As the model clouds' initial velocities were 100–300 km s⁻¹ (see Table 1), this implies that Complex C has moved ≳10–30 kpc from its initial location to its current location, suggesting that Complex C likely originated as extragalactic material. Note, however, that the measured metallicity of Complex C is too high for a purely extragalactic origin (Tripp et al. 2003; Collins et al. 2003) and requires a metal-enrichment process. In preliminary calculations, we have found that mixing with the Galactic gas noticeably enhances the metallicity of the H i HVC gas, but a full discussion is beyond the scope of our current study.

As noted above, there is significant sightline-to-sightline variation in the observed ratios. Sembach et al. (2003) reported a weak trend for the ratio of O vi to H i to increase toward lower Galactic longitude for nine sightlines through Complex C. Complex C is oriented diagonally on the Galactic coordinate grid, running from low latitudes and longitudes to high latitudes and longitudes (e.g., Wakker & van Woerden 1991), with the higher-latitude, higher-longitude region of the complex being higher above the disk (Thom et al. 2008). Hence, the trend reported by Sembach et al. (2003) means that the ratio of O vi to H i increases toward lower Galactic latitude and lower height above the disk. The ratio of C iv to H i also increases toward lower latitudes, although there are fewer data points for this trend (the N v measurements are mostly upper limits).

This variation of $N(\mbox{O\,{\sc {vi}}})/N(\mbox{H\,{\sc {i}}})$ with longitude was discussed by Tripp et al. (2003), who suggested that the lower-longitude, lower-latitude region of Complex C is closer to the disk, where the ambient medium is denser. They suggested that this greater density leads to a more vigorous interaction between the ambient medium and the HVC, leading to greater production of O vi and greater ionization of H i. Our simulations are unsuitable for directly examining this suggestion, as we have not investigated the effect of an increasing ambient density on high ion production. An increasing ambient density will lead to both greater ram-pressure stripping of material, and faster-growing shear instabilities, because of the lower density contrast between the cloud and ambient medium (Chandrasekhar 1961, Section 101). Both of these effects should lead to greater production of high ions.

An alternative explanation for the observed $N(\mbox{O\,{\sc {vi}}})/N(\mbox{H\,{\sc {i}}})$ trend is that the region of Complex C nearer the disk is at a later stage in its evolution, having passed through more ambient medium, and so has a higher high-ion-to-H i ratio (see Figure 11). However, the range of heights spanned by Complex C (∼2 kpc, using data from Thom et al. 2008) divided by a velocity of ∼100 km s⁻¹ corresponds to a timescale of only ∼20 Myr. From Figure 11, we can see that this timescale is insufficient to explain the range of observed ratios, which vary by more than an order of magnitude from sightline to sightline. The trend could also be explained by the region of Complex C toward lower longitudes and latitudes being initially composed of smaller clouds than the higher-longitude, higher-latitude region, which give rise to larger high-ion-to-H i ratios (see Figure 11). However, there is no obvious reason for the size of the clouds of which Complex C was composed when it began interacting with the ISM to vary systematically with position.

Another possible explanation for the observed $N(\mbox{O\,{\sc {vi}}}){/} N(\mbox{H\,{\sc {i}}})$ trend is that the turbulently mixed material toward lower longitudes and latitudes has a higher metallicity. Such a difference in metallicity could arise from the mixing of metal-poor Complex C gas with the metal-rich ISM. If the metallicity of the halo decreases with height and/or Galactocentric distance,⁶ then the lower-longitude, lower-latitude region of Complex C will be interacting with higher-metallicity gas than the higher-longitude, higher-latitude region. As a result, the turbulently mixed material toward lower longitudes and latitudes would have higher abundances and thus higher ratios of high ions to H i. If this speculation is correct, then variation in elemental abundances with latitude may be apparent from other ions. O i, for example, does not show an anticorrelation between abundance and latitude (using data from Collins et al. 2007), but these abundance measurements are probably probing undisturbed cloud material rather than turbulently mixed material. More measurements of elemental (as opposed to ionic) abundances for Complex C are needed to test whether or not the abundances vary with latitude.

4.5. Ion-to-Ion Column Density Ratios

Unlike the ion-to-H i ratios, the observed ratios of C iv and N v to O vi do not vary greatly from sightline to sightline, nor do the predicted ratios vary systematically with cloud size or time. The observed C iv-to-O vi ratios are in good agreement with the predicted values; the agreement is even better when the predicted ratios are shifted according to the Complex C abundances. For N v to O vi, all but the Mrk 279 measurement from Fox et al. (2004, filled gray square) are upper limits, and so do not strongly constrain the models. If the ratio of gas phase nitrogen to oxygen is equal to the solar ratio, then our results are in good agreement with the observed column density ratio. However, if the ratio of gas phase nitrogen to oxygen is about 1/10 solar, as it is on the PG 1259+593 sightline (Fox et al. 2004), then our results are overabundant in N v compared to O vi by about a factor of 10. Note that PG 1259+593 is more than 10 deg higher in Galactic latitude than Mrk 279, so it is not implausible that the relative abundances along these sightlines could differ.

Another way to look at the ratios of C iv and N v to O vi is to plot them in $N(\mbox{N\,{\sc {v}}})/N(\mbox{O\,{\sc {vi}}})$ – $N(\mbox{C\,{\sc {iv}}})/N(\mbox{O\,{\sc {vi}}})$ space (Figure 13). Although the ratios for ions moving with HVC-like velocities predicted by the simulations in this paper overlap the ion ratios predicted by our previous numerical models of TMLs (KS10), the ratios from some our models are lower than KS10's, and better approach the ratios observed along three sightlines toward Complex C. However, as mentioned earlier, if our results are rescaled according to the Complex C abundances measured along the PG 1259+593 sightline (Fox et al. 2004), our predicted data points are shifted by −0.97 and −0.35 along the x- and y-axes, respectively. After such a shift, some of our models' predictions would match the observed ratios (note that the observed $N(\mbox{N\,{\sc {v}}})/N(\mbox{O\,{\sc {vi}}})$ ratios are upper limits, so our model predictions would be consistent with the observations, even after a shift of ≈1 dex to the left). Figure 13 confirms that ablation of cloud material followed by mixing with the hot ambient gas is a viable physical process for the production of high-velocity high ions. Note that some other models shown in Figure 13 are in agreement with the observed ion ratios toward Complex C. However, these models do not distinguish between high- and low-velocity ions; such a distinction is needed for an accurate comparison with observations of HVCs.

5. NEGLECTED PHYSICAL PROCESSES AND THEIR EFFECTS

Our simulations used a 2D geometry because the memory requirements for tracking the ionization evolution of carbon, nitrogen, and oxygen in three dimensions are prohibitive. This assumed geometry prevented us from modeling realistic magnetic field configurations, and so magnetohydrodynamic effects were ignored. We also neglected thermal conduction, external heating, such as photoelectric heating by Galactic or extragalactic UV radiation, and photoionization.

KS10 discuss the effect of assuming a 2D geometry on turbulent mixing, compared with a 3D geometry, as well as the effect of a magnetic field. Generally, turbulent mixing is more severe in 3D simulations than in 2D, which would lead to greater production of high ions, while a magnetic field is known to suppress turbulence (see references in KS10). It is unclear whether the combination of a magnetic field and a 3D geometry would result in more or fewer high ions, compared with our 2D simulations. In addition, our previous 3D MHD simulations showed that the dynamical effect of the magnetic field depends strongly on the orientation of the field (Kwak et al. 2009). As a result, the formation and size of the TMLs, and thus the production of high ions, would also be affected by the orientation of the magnetic field. However, it should be noted that, for a plane-parallel geometry, the 3D MHD TML simulations of Esquivel et al. (2006) predicted similar high ion column density ratios as the 2D NEI simulations of KS10. Finally, it should be noted that the magnetic field strength is likely to be small in the upper halo, which is the region of interest of this study.

Our simulations did not include thermal conduction. Vieser & Hensler (2007a, 2007b) studied numerically the effect of thermal conduction on the evaporation of a cool cloud embedded in a hot ambient medium. They found that thermal conduction suppresses shear instabilities by smoothing out the steep temperature and density gradients at the interface between the hot and cool gas (Vieser & Hensler 2007b). However, in our simulations, except for Model E, there are not steep gradients at the initial interface between the cool cloud and the hot ambient medium (see Figure 1). Also, in our simulations, heat diffusion due to thermal conduction is unlikely to be important relative to heat diffusion by turbulence (Esquivel et al. 2006; KS10). For the mixed gas in our simulations (T ∼ 10⁵ K, n ∼ 10⁻² cm⁻²), the Spitzer thermal diffusion coefficient⁷ is κ_Sp ∼ 10²⁴ cm² s⁻¹. In contrast, the turbulent diffusion coefficient $\kappa _\mathrm{turb}\sim \mbox{$\frac{1}{3}$}v_\mathrm{turb}L_\mathrm{inj}$ (Cho et al. 2003) is ≳10²⁵ cm² s⁻¹, where the turbulent velocity is v_turb ∼ 100 km s⁻¹ and the energy injection scale is L_inj ≳ 1 pc (note that L_inj corresponds to the size of the largest eddies; Cho et al. 2003). These diffusion coefficients mean that turbulent diffusion of heat would dominate over diffusion by thermal conduction, even if thermal conduction were included in our simulations.

We included radiative cooling in our simulations (albeit under the assumption of CIE), but not photoelectric heating due to UV irradiation of dust grains. Including such heating would, in principle, reduce the effect of radiative cooling. However, in practice, the photoelectric heating rate is unlikely to be important for our simulations of HVCs in the upper halo, because of the decrease in the UV radiation field and the dust density with height. For example, the photoelectric heating rate for the upper halo (|z| ≳ 3 kpc) used by Joung & Mac Low (2006) in their ISM models would be ≲10⁻³² erg cm⁻³ s⁻¹ in our mixed gas, whereas the radiative cooling rate is ∼10⁻²⁵ erg cm⁻³ s⁻¹ (using the 1993 version of the Raymond & Smith 1977 code).

A more important cooling-related effect is that the lower metal abundances seen on some sightlines through Complex C (Fox et al. 2004; Collins et al. 2007) would reduce the cooling rate, as the cooling is dominated by metal line emission. Radiative cooling tends to stabilize a cool cloud against disruption as it moves through hot ambient gas (Vietri et al. 1997), and so reduced cooling rates would tend to accelerate the disruption of the cloud and so may lead to the production of more high ions than are currently predicted by our simulations. We defer investigation of this effect to a later study.

Finally, we did not consider radiative transfer, and the possibility that some of the ion populations may be affected by photoionization. Including photoionization is beyond the scope of this study. However, we note that, if it is important, photoionization by escaping Galactic radiation would tend to increase the amount of C iv relative to O vi.

6. SUMMARY

We have presented the results of hydrodynamical simulations of HVCs that self-consistently trace the NEI evolution of carbon, nitrogen, and oxygen. This is the first time such simulations have been carried out with an initially spherical cloud moving through hot, tenuous ambient gas. We have concentrated particularly on the production of high ions (C iv, N v, and O vi) in the TMLs that form between the cold clouds and the hot ambient medium. In contrast to the previous plane-parallel TML models (e.g., KS10), the new simulations enable us to investigate how the cold gas ablates from an initially spherical cloud and subsequently mixes with the hot ambient gas.

Material is constantly being ablated from the clouds (Section 3.1). This initially cold material mixes with the hot ambient medium, and becomes ionized, so the clouds are constantly losing their H i content (Section 3.2). However, despite this loss of material, we find that HVCs can survive for long times: at least a few hundred megayears for a cloud with initial mass ∼4 × 10⁵ M_☉. Note that clouds that are currently ∼10 kpc from the disk (e.g., Complex C with z ∼ 8 kpc; Thom et al. 2008) and traveling toward the disk at ∼100 km s⁻¹ would require only ∼100 Myr to reach the disk. It is therefore possible that the more massive regions of HVC complexes, if not composed of many small cloudlets, may reach the disk at least partially as neutral hydrogen (Section 3.4).

High ions are produced in the mixed gas, both by ionization of the initially cool cloud material and recombinations in the initially hot ambient material (Sections 3.1 and 3.5). The high ions initially appear while the mixed gas is still traveling at HVC-like velocities. The material stripped from the clouds slows while remaining abundant in high ions, and ultimately ends up with halo-like velocities (Paper II), ceasing to be recognizable as HVC material.

The simulations may also explain why some high-velocity highly ionized gas is seen on sightlines that lack high-velocity H i (Sembach et al. 2003), at least for sightlines that have high-velocity H i nearby (Section 3.5). At the edges of the clouds, the high-velocity H i column density can fall below the 21 cm detection limit ( $N(\mbox{H\,{\sc {i}}})\lesssim 10^{18} \;\mbox{cm}^{-2}$ ), while the column densities of high-velocity high ions remain detectable (e.g., $N(\mbox{O\,{\sc {vi}}})\ \gtrsim 10^{13} \;\mbox{cm}^{-2}$ ). In addition, the high-velocity high ions exist mainly in material that has ablated from the main body of the cloud and mixed with the ambient gas. This material continues to travel at HVC-like velocities, although it has a tendency to fall behind the cloud. As a result, sightlines that cross the ablated material but not the cloud (such sightlines would cross our domain diagonally) would intersect large numbers of high ions but small numbers of neutral hydrogen atoms.

We investigated a suite of models with a range of model parameters. We found that the cloud's initial velocity does not affect the rate of H i loss (Section 3.3.1) or high ion production (Section 3.5). The rate of H i loss also does not strongly depend on the initial density of the cloud or its density profile (Section 3.3.2). However, the cloud's initial size does affect the rate of H i loss—a smaller cloud loses its H i mass relative to its initial mass more rapidly than a larger cloud. This difference is at least partially due to the fact that a smaller cloud has a larger surface-area-to-volume ratio (Section 3.3.3).

Significant column densities of high ions are produced in our simulations. In some models, at some times, these column densities are large enough to be observed, and are in reasonable agreement with the ion column densities observed for Complex C (Section 4.1). As the ion column densities generally rise throughout our simulations, it is likely that the predicted column densities would be in better agreement with the observed values if our simulations were continued beyond their current end times. Our models also compare well with the Complex C observations in terms of the trend between $N(\mbox{O\,{\sc {vi}}})$ and $N(\mbox{H\,{\sc {i}}})$ (Section 4.2), the ion-to-H i ratios (Section 4.4), and the ion-to-ion ratios (Section 4.5). We therefore conclude that the ablation of cloud material and its subsequent turbulent mixing with the hot ambient medium is a viable mechanism for the production of high-velocity high ions.

The column density ratios of high ions to H i depend on the initial cloud size. A smaller cloud reaches a given ion-to-H i ratio at an earlier time than a larger cloud. If Complex C were initially composed of clouds with r ≳ 150 pc (which cannot be ruled out, given the lack of knowledge of Complex C's 3D structure), then the observed ion-to-H i ratios would imply an age for Complex C of ≳100 Myr. This in turn implies that Complex C has traveled ≳10–30 kpc from its initial location, suggesting that Complex C formed from extragalactic material. Furthermore, if parts of Complex C are composed of clouds as large as our Model F cloud (M ∼ 4 × 10⁵ M_☉), these clouds may reach the disk partially in the form of neutral hydrogen, as noted above. HVCs could therefore provide new material to the disk for star formation. We also noted that there is a trend in Complex C for the O vi-to-H i ratio to increase as one gets closer to the disk, and speculated on possible reasons for this trend in Section 4.4.

In this paper, we have concentrated particularly on high-velocity ions in the Galactic halo. However, low-velocity high ions are also observed in the halo, and, as noted above, low-velocity ions are produced abundantly in our simulations. These low-velocity ions will be discussed in Paper II.

We thank Orly Gnat for providing us with the data from Gnat et al. (2010) plotted in our Figure 13. We also thank the anonymous referee whose comments have helped improve this paper. The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. The simulations were performed at the Research Computing Center (RCC) of the University of Georgia. This work was supported by NASA grant NNX09AD13G, awarded through the Astrophysics Theory and Fundamental Physics Program. D.B.H. acknowledges funding from NASA grant NNX08AJ47G, awarded through the Astrophysics Data Analysis Program.

SIMULATIONS OF HIGH-VELOCITY CLOUDS. I. HYDRODYNAMICS AND HIGH-VELOCITY HIGH IONS

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL METHODS AND MODEL PARAMETERS

3. RESULTS

3.1. Hydrodynamical Evolution of High-velocity Clouds

3.2. Evolution of the Number of H i Atoms