HYDRODYNAMICAL COUPLING OF MASS AND MOMENTUM IN MULTIPHASE GALACTIC WINDS

We seek to model the impact of a hot, supernova-driven outflow on cooler material. Given a small set of assumptions, including spherical symmetry and negligible radiative cooling, the hot phase of the wind can be modeled analytically at distances close to the plane of a galaxy. All of our simulations use an analytic model of the hot wind as a constant background state with properties set using the adiabatic wind model of Chevalier & Clegg (1985, hereafter CC85). The CC85 model envisions a hot wind driven by central energy and mass input from stellar feedback processes. Three input parameters determine the solutions to the model: the energy input as a function of time, $\dot{E}$, the mass input as a function of time, $\dot{M}$, and the size of the driving region within which energy and mass are injected, R_*. By $\dot{M}$ we mean the total mass injection rate to the wind due to supernovae and mass loading, not the star formation rate. With these parameters, a solution to the set of spherical hydrodynamic equations can be found that smoothly transitions from subsonic within the driving region to supersonic at further radii. The solutions cross the sonic point at $r={R}_{* }$.

We choose the input parameters for our version of the CC85 model according to the fits derived by Strickland & Heckman (2009) using Chandra X-ray observations of the nearby starburst galaxy M82. In particular, we set $\dot{E}={10}^{42}$ erg ${{\rm{s}}}^{-1}$, $\dot{M}=2\,{M}_{\odot }$, ${\mathrm{yr}}^{-1}$, and ${R}_{* }=300$ pc. In interpreting their results, we have made additional assumptions about the wind mass-loading factor β and the supernova thermalization fraction α, for which they give a range of correlated values. Here we are using the Strickland & Heckman (2009) interpretation of β, meaning the fraction of total mass injected into the wind as compared to the mass injected by supernovae and stellar winds, $\beta =\dot{M}/{\dot{M}}_{\mathrm{SN}+\mathrm{SW}}$. Likewise, our definition of α corresponds to their and refers to the fraction of the supernova energy that is deposited in low-density gas and does not suffer large radiative losses before being incorporated into the wind. We take $\beta =1.42$ and $\alpha =0.33$, values near the middle of the acceptable range of fits reported by Strickland & Heckman (2009). These choices give a central temperature of ${T}_{{\rm{c}}}\gtrsim {10}^{7.5}$ K in the driving region, consistent with estimates made from the X-ray emission (see Strickland & Heckman 2009, Table 2). The resulting values for number density, velocity, temperature, and Mach number as a function of radius for this model are displayed in Figure 1.

Zoom In Zoom Out Reset image size

For the background flow in our multiphase wind simulations, we use the physical parameters of the CC85 wind model shown in Figure 1 at r = 1 kpc. These values are

$$\begin{eqnarray*}{n}_{\mathrm{wind}} & = & 5.2626\times {10}^{-3}\,{\mathrm{cm}}^{-3},\\ {v}_{\mathrm{wind}} & = & 1.1962\times {10}^{3}\,\mathrm{km}\,{{\rm{s}}}^{-1}=1.2225\,\mathrm{pc}\,{\mathrm{kyr}}^{-1},\\ {P}_{\mathrm{wind}}/k & = & 1.9881\times {10}^{4}\,{\mathrm{cm}}^{-3}\,{\rm{K}},\end{eqnarray*}$$

where k is the Boltzmann constant. At r = 1 kpc, the wind pressure corresponds to a temperature of ${T}_{\mathrm{wind}}=3.7778\ \times {10}^{6}$ K, as displayed in Figure 1. We choose a radius of 1 kpc to set the background wind properties in our simulations for several reasons. First, we wish to capture the mass loading of the wind outside the driving region, which restricts us to hot wind properties at radii $r\gt 300$ pc. Second, we will model the interactions between cool and hot material in a simulation volume with a physical length of 160 pc. Given that the wind properties in our simulations remain approximately constant across this volume and the hot wind density, temperature, and Mach number change most rapidly just outside the driving region (see Figure 1), we favored an initial radius of $r\sim 1\,\mathrm{kpc}$ over smaller radii. Further, the best observations of a multiphase wind come from M82, where cool material clearly resides at $r\gt 1\,\mathrm{kpc}$ above the disk (Leroy et al. 2015). The fits derived by Strickland & Heckman (2009) suggest that the hot wind should not yet have suffered serious radiative losses at $r\sim 1$ kpc, which would invalidate the CC85 model (Zhang et al. 2015; Thompson et al. 2016).

To capture the multiphase nature of galactic winds, our simulations also include a cool component representing interstellar or circumgalactic material. As with previous studies of galactic winds, we model this second component as dense clouds initially at rest with respect to the hot wind (e.g., Scannapieco & Brüggen 2015; Brüggen & Scannapieco 2016). Our aim is to extend previous studies by examining the detailed momentum and energy coupling between different wind phases, so we consider both idealized spherical clouds and more realistic turbulent clouds with a distribution of interior densities set by turbulent processes. Observations have revealed cool gas in outflows at a variety of densities and temperatures (see references in Section 1). By varying the median density of the cold gas, $\tilde{n}$, and its interior density structure, we are able to model a range of properties for the dense component of the wind. Both the median number density and the cloud morphology affect the integrated surface density of the cold gas, ${{\rm{\Sigma }}}_{\mathrm{cl}}={M}_{\mathrm{cl}}/\pi {R}_{\mathrm{cl}}^{2}$, where M_cl is the total mass of the cloud and R_cl is the cloud radius. The surface density influences how momentum transfers from the hot wind to the cold component, as discussed below.

As the hot wind destroys the clouds, cool material will both heat and rarify. This material will accelerate as momentum transfers from the hot wind, and enter the outflowing wind with a range of velocities. To adequately track all of this material, we perform simulations using boxes with long aspect ratios. The simulation boxes feature transverse dimensions equal to 8R_cl and a long dimension along the wind direction of 32R_cl. Figure 2 displays an example of an entire box, showing a density projection of the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ turbulent cloud at time t = 0. The clouds are initially centered 2R_cl from the left boundary of the box to capture the resulting bow shock. The long dimension of the box enables the simulations to follow the bulk of the cloud as it accelerates and track material stripped from the main cloud body.

Zoom In Zoom Out Reset image size

2.2.1. Cloud and Sphere Models

Each cloud initially sits at rest and in thermal pressure equilibrium with the surrounding hot wind. The clouds in our study are not dense enough to be gravitationally confined, allowing us to neglect self-gravity. The density of the cloud material therefore determines its initial temperature. We initialize the spherical clouds with constant interior temperature and initialize the interior temperatures of regions within the turbulent clouds along an appropriate isobar. We list the temperatures and median and mean densities of the cloud initial conditions in Table 1. We list the full range of temperatures for the turbulent clouds—their median temperature matches the spheres. For both the spheres and turbulent clouds we taper the densities at the edge, starting at a radius of 4.5 pc, such that the cloud density smoothly transitions from the median to the wind density. The density taper has the form

$$\begin{eqnarray}&&n{(r)}_{\mathrm{cl}}=\tilde{n}\,\exp [\mathrm{ln}({n}_{\mathrm{wind}}/\tilde{n})/({R}_{\mathrm{cl}}-4.5)]| r-4.5| ,\end{eqnarray} \tag{ 1 }$$

where r is the distance from the center of the cloud and $\tilde{n}$ is the median cloud density. ${R}_{\mathrm{cl}}=5$ pc for all of our clouds. To create the turbulent clouds, we excise a region from a Mach 5 isothermal turbulence simulation (Robertson & Goldreich 2012) and scale the density linearly to match the desired median. Figure 3 shows a histogram of the initial gas densities for the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ turbulent cloud simulation.

Zoom In Zoom Out Reset image size

Table 1.  Summary of Simulations Performed, and Physical Parameters of Interest for the Cloud in Each Simulation

Model	${\tilde{n}}_{\mathrm{cl}}$ a	${\bar{n}}_{\mathrm{cl}}$ b	${{\rm{\Sigma }}}_{\mathrm{cl}}$ c	NHd	Tcle	Mclf	χg	tcch	tcooli
Name	(cm−3)	(cm−3)	(M⊙ pc−2)	(cm−2)	(K)	(M ${}_{\odot }$)		(kyr)	(kyr)
S01	0.1	0.086	0.013	2.87 × 1018	1.988 × 105	1.05	1.9 × 101	17.8	26.0
S05	0.5	0.45	0.064	1.43 × 1019	3.976 × 104	5.04	9.5 × 101	39.8	3.96
S1	1.0	0.88	0.13	2.85 × 1019	1.988 × 104	10.0	1.9 × 102	56.4	1.30
T01	0.1	0.16	0.022	2.05 × 1019	2.2 × 103	1.73	1.9 × 101	17.8	26.0
T05	0.5	0.83	0.11	1.02 × 1020	4.4 × 102	8.61	9.5 × 101	39.8	3.96
T1	1.0	1.7	0.24	2.23 × 1020	2.0 × 102	18.8	1.9 × 102	56.4	1.30

Notes. Simulations with a spherical cloud begin with an "S"; those with a turbulent cloud begin with a "T." All simulations listed in Table 1 are run in a volume with a numerical resolution of $2048\times 512\times 512$ cells, corresponding to a physical box size of $160\times 40\times 40$ pc.

^aMedian density of the cloud, calculated including all material with a density greater than 1/10 ñ. ^bMean density of the cloud, calculated including all material with a density greater than 1/10 $\tilde{n}$. ^cSurface density of the cloud, calculated as ${M}_{\mathrm{cl}}/(\pi {R}_{\mathrm{cl}}^{2})$. ^dMaximum initial column density sight line through the cloud (excluding the hot wind). ^eInitial temperature of the cloud material, in pressure equilibrium with the wind. The median temperature is provided for the spherical clouds, and the temperature of the highest-density regions is provided for the turbulent clouds. ^fInitial mass of the cloud, calculated including all material with a density greater than 1/10 $\tilde{n}$. ^gInitial density contrast of the cloud with the background wind, calculated using $\chi ={\tilde{n}}_{\mathrm{cl}}/{n}_{\mathrm{wind}}$. ^hCloud crushing time, calculated using Equation (2) with the median density and a cloud radius of 5.0 pc. ⁱCloud cooling time, calculated using Equation (3) with the median cloud density and temperature.

Download table as:  ASCIITypeset image

2.2.2. Cloud Crushing Time

To interpret our results in relation to previous studies of cloud–shock interactions, we make use of the "cloud crushing time" t_cc initially devised by Klein et al. (1994). When the hot wind first interacts with the cool cloud, a shock drives into the overdense gas (and a reverse shock reflects into the oncoming wind). The cloud crushing time estimates how long the initial shock takes to pass through and compress the cloud. This timescale can be calculated in terms of known quantities by relating the initial density contrast of the cloud to the wind $\chi ={n}_{\mathrm{cl}}/{n}_{\mathrm{wind}}$ along with the radius of the cloud R_cl and the wind velocity v_wind as

$$\begin{eqnarray}&&{t}_{\mathrm{cc}}={R}_{\mathrm{cl}}{\chi }^{\tfrac{1}{2}}/{v}_{\mathrm{wind}}.\end{eqnarray} \tag{ 2 }$$

We use the cloud crushing time as an evolutionary timescale for the clouds in our simulations and report t_cc for each of our simulations in Table 1. The derivation leading to Equation (2) assumes that the pre-shocked gas in the cloud evolves adiabatically, allowing for an estimate of the shock speed within the cloud from the density contrast and wind speed. In a radiatively cooling cloud the pre-shock conditions change as the cloud cools, decreasing the sound speed and slowing the shock. The cloud crushing times listed in Table 1 therefore represent imperfect estimates of the cloud compression timescale, and our simulations show that they underestimate the duration of this phase. In our simulations, maximum compression is typically reached at ∼2t_cc.

2.2.3. Cloud Cooling Time

The cooling time provides another important evolutionary timescale for the clouds we simulate. If the cloud cooling time greatly exceeds the cloud crushing time, we expect the clouds to evolve similarly to the well-studied adiabatic case (Klein et al. 1994; Xu & Stone 1995; Poludnenko et al. 2002; Nakamura et al. 2006). If the cloud crushing time grows much longer than the cooling time, the radiative loss of energy may affect the cloud properties substantially before the initial shock can disrupt it. We can estimate the cooling time of the clouds in our simulation as their thermal energy divided by the rate at which energy is lost owing to radiative cooling as

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{3{n}_{\mathrm{cl}}{{kT}}_{\mathrm{cl}}}{2{\rm{\Lambda }}({T}_{\mathrm{cl}})},\end{eqnarray} \tag{ 3 }$$

where k is the Boltzmann constant and ${\rm{\Lambda }}(T)$ is the cooling rate in erg ${{\rm{s}}}^{-1}$ ${\mathrm{cm}}^{-3}$, evaluated at the cloud temperature (see Appendix C for information on our calculation of cooling rates). In theory, we might like to use T_cl and n_cl post-shock to calculate the cooling time. Given the complex nature of the cooling function, these quantities often prove impossible to calculate analytically (Creasey et al. 2011). Additionally, the lognormal gas distribution in our turbulent clouds means that n and T vary for different parts of the cloud. Instead, we calculate the cloud cooling times using the median pre-shock density and temperature to compare most directly with the cloud crushing time. These cooling times appear in Table 1. Our simulations sample the transition from clouds dominated by adiabatic evolution to those in the radiatively cooling regime. As discussed in Section 4, our simulations reproduce this transition and thereby justify our assumptions used to compute the cooling time.

In Figures 4 and 5 we illustrate the general evolution of the $\tilde{n}=1\,{\mathrm{cm}}^{-3}$ turbulent and spherical cloud simulations (models T1 and S1; Table 1). These figures show the column density of gas in the box (including the hot wind), in both a y-axis projection with the hot wind entering the box from the left and an x-axis projection in the direction of the wind. The figures show snapshots of the simulations at five representative times—the initial conditions, and after 2t_cc, 5t_cc, 10t_cc, and 20t_cc.

Zoom In Zoom Out Reset image size

In both simulations, the maximum average cloud density is reached at t = 2t_cc. The maximum average density is $\bar{n}=5.2$ ${\mathrm{cm}}^{-3}$ for T1 and $\bar{n}=6.8$ ${\mathrm{cm}}^{-3}$ for S1. (We calculate the mean density using material with density greater than 1/10 the initial median density, as in Table 1.) Despite the similarity of this compression timescale, Figures 4 and 5 show a drastic difference in cloud morphology at t = 2t_cc. While the initial shock has propagated at different speeds through regions of different density for the turbulent cloud, the shock has very effectively compressed the sphere into a single flat pancake. As a result, the evolution of the two types of clouds diverges strongly at later times. Low-density material is quickly accelerated between $t=5\mbox{--}20$t_cc in the turbulent cloud, leaving only a few high-density cloudlets at late times. We quantify this rapid mass loss in Figure 6, which shows the normalized cloud mass as a function of cloud crushing time. The mass is calculated using material at or above 1/3 the initial median density. By t = 20t_cc, nearly 80% of the original mass of the turbulent cloud has been mixed into the hot wind and fallen below the density threshold of n = 0.33 ${\mathrm{cm}}^{-3}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The loss of material proceeds more rapidly for turbulent clouds than for spheres, even with the cloud crushing time calculated as a function of the median density. As can be seen in the third panel of Figure 5, after the initial pancake stage, the spherical cloud compresses into a single core, with a small surface area and high column density in the wind direction. This morphology is a result of the original shock moving in the wind direction combined with the compression from shocks on the sides of the cloud, which push material toward the center. This high-density core presents a smaller surface area for ablation of material than the many small high-density regions in the turbulent cloud (compare the third panels of Figures 4 and 5). The single bow shock at the front of the spherical cloud protects the core, while also resulting in a pressure gradient through the cloud. Cloud elongation in the direction of the wind gradually breaks up the original core into smaller fragments, each of which has an individual bow shock and loses material. However, as demonstrated in Figure 6, the overall process is much slower for spheres than for turbulent clouds. Only ∼40% of the original spherical cloud material has been lost by t = 20t_cc.

Figures 4 and 5 show the evolution of clouds with a median density of $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$. We have also run four high-resolution simulations with lower median densities. Models T05 and S05, the turbulent and spherical clouds with a median density of $\tilde{n}=0.5$, show a mass and morphology evolution similar to their higher-density counterparts. Figure 6 shows that the mass loss as a function of cloud crushing time is nearly identical for the $\tilde{n}=0.5$ ${\mathrm{cm}}^{-3}$ and $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ clouds out to 25t_cc. However, as the original median density continues to decrease, the evolution begins to follow a qualitatively different path. The original shock that passes through the cloud does not compress the gas enough for it to reach densities where it can cool efficiently. The warm cloud gas quickly gets rarified and accelerated with the hot wind, causing the mass loss of the clouds to proceed much more rapidly at lower original median densities. This difference is clearly visible in the $\tilde{n}=0.1$ evolutionary tracks in Figure 6.

In model S01, almost none of the gas reaches the requisite density to cool. In fact, only a small ring of gas within the cloud reaches the threshold density. This ring is visible in Figure 7, which shows surface density projections of the $\tilde{n}=0.1$ sphere simulation at 2t_cc. The ring structure is caused by the shock moving in the wind direction colliding with the oblique shocks coming in from the sides of the cloud. While the densities in this collision do not get large enough to result in a single compact core as seen in the higher-density models, this small amount of dense gas is enough to cause the extended tail seen in the mass evolution track of Figure 6. However, most of the cloud is destroyed in less than 5t_cc, consistent with adiabatic simulations of cloud–shock interactions (Schneider & Robertson 2015).

Zoom In Zoom Out Reset image size

The evolution of the $\tilde{n}=0.1\,{\mathrm{cm}}^{-3}$ turbulent cloud (model T01) is less dramatically different as compared to its higher-density counterparts. Although the median cloud density is below that required for effective cooling, the lognormal density distribution of initial cloud material spans a range up to $n\approx 10$ ${\mathrm{cm}}^{-3}$. As a result, parts of the cloud that were initially above the median density are still able to cool effectively after being shocked. As Figure 6 shows, the post-shock mass loss for the $\tilde{n}=0.1$ turbulent cloud proceeds faster than the $\tilde{n}=1$ or $\tilde{n}=0.5$ clouds, but more slowly than the $\tilde{n}=0.1$ sphere. We expect that the speed of mass loss for turbulent clouds would continue to increase as the initial median density is lowered, until eventually the evolution proceeded on an adiabatic track (shown by the dotted black line in Figure 6).

Figure 8 shows density–temperature phase diagrams for the $\tilde{n}=1$ turbulent cloud at four different times in its evolution—the initial conditions, and t = 5t_cc, 10t_cc, and 20t_cc. Each bin is colored according to the total mass it contains, to better focus attention on the locations with the majority of the cloud material. Throughout the simulation, most of the mass is in the hot wind, visible as the red region at the upper left of the distribution in each panel (and as a single bin in the initial conditions). At t = 0, the cloud's mass is distributed across a range of densities and temperatures, with most of the mass in regions at or above the median density of $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$. Each panel also contains an isobar showing the original pressure of the gas (recall that the cloud's pressure is matched to the wind in the initial conditions), as well as a dashed line that shows the equilibrium location between heating and cooling in our simulations.

Zoom In Zoom Out Reset image size

After the cloud has been shocked, the density–temperature phase diagram takes on a characteristic shape. At early times, most of the gas is at higher pressure than the initial conditions. As the simulation proceeds, the cloud gas slowly evolves back toward thermal pressure equilibrium with the incoming wind. Cloud material quickly fills out the entire range of densities between the original cloud density and the wind. A large amount of mass has been compressed to high densities, much of it over an order of magnitude higher than the initial densities. The high-density material cools effectively, with much of the cloud mass located just above the equilibrium cooling temperature, at the lower right in each panel. There also exists a mass concentration around $\mathrm{log}(n)=0.5$ and $\sim 2\times {10}^{4}$ K. This buildup reflects the shape of the cooling curve, with maximally efficient cooling around 10⁵ K and a steep falloff around 10⁴ K. These features remain throughout the subsequent evolution, as the mass in high-density bins is gradually reduced.

The phase diagrams for simulations T01 and S01 yield an estimate of the density threshold required for efficient cloud cooling in our simulations. Figure 9 shows the mass-weighted density–temperature phase diagrams for both simulations at t = 2t_cc. The turbulent cloud, displayed in the left panel, has a large amount of mass in the cooling sweet spot just above ${\mathrm{log}}_{10}(n)=0$. In contrast, most of the shocked sphere material never crosses this density threshold, and as a result, the gas remains at high temperatures and low densities, where it quickly mixes with the hot wind. Animations of the phase diagrams show clearly the path that the cool gas takes in the turbulent cloud simulation. Cloud gas originally above the median density shocks to densities just above ${\mathrm{log}}_{10}(n)=0$, after which the gas begins to cool rapidly down to 10⁴ K. Only a small fraction of the sphere material reaches densities of ${\mathrm{log}}_{10}(n)=0$, and as a result, most of the mass of the cloud is lost at early times.

Zoom In Zoom Out Reset image size

While cool clouds have been observed in galactic outflows at a variety of distances and velocities, the primary physical mechanism responsible for fast-moving cool gas continues to be debated. In this section, we investigate the velocity of the cool gas in our simulations and show that in the purely hydrodynamic case, hot winds are unlikely to accelerate cool gas to the high velocities seen in many outflows.

Figure 10 shows mass-weighted density–velocity histograms for model T1. The figure shows two representative snapshots, at t = 5t_cc and 15t_cc, with the velocity in the wind direction plotted on the y-axis. As in the density–temperature diagrams, the hot wind appears as a mass concentration at the upper left in each panel, with a small spread around the initial wind density and velocity. While cloud material has clearly acquired a range of densities, only low-density material travels at high speeds. Only 3% of the cloud mass above the original median density of $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ is moving faster than 200 km ${{\rm{s}}}^{-1}$ by 15t_cc. The average velocity of the dense material ($n\gt 1$ ${\mathrm{cm}}^{-3}$) is only ${v}_{x}\approx 120$ km ${{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

One of the goals of this work was to investigate how different cloud density structures change the effectiveness of the hot wind in accelerating high-density material. This question is addressed in Figure 11, which compares the density–velocity histograms for the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ turbulent cloud at t = 10t_cc to the sphere. This figure indicates that accelerating high-density material proves more difficult in the turbulent cloud case. This difference results from the different initial column densities of the cloud material, as we describe below.

Zoom In Zoom Out Reset image size

In a simple model, the ram pressure of the hot wind will accelerate cold cloud material (we neglect other forces in the following analysis). The ram pressure of the hot wind scales as

$$\begin{eqnarray}&&{P}_{\mathrm{ram}}={n}_{\mathrm{wind}}{v}_{\mathrm{wind}}^{2},\end{eqnarray} \tag{ 4 }$$

which gives a constant ram pressure per unit mass of ${P}_{\mathrm{ram}}/{m}_{{\rm{H}}}=7.2\times {10}^{13}$ ${\mathrm{cm}}^{-1}\,{{\rm{s}}}^{-2}$ for our background wind conditions. The associated acceleration of the cloud, g_cl, will then be

$$\begin{eqnarray}&&{g}_{\mathrm{cl}}=\displaystyle \frac{{P}_{\mathrm{ram}}}{{m}_{{\rm{H}}}}\displaystyle \frac{{m}_{{\rm{H}}}}{{{\rm{\Sigma }}}_{\mathrm{cl}}}=\displaystyle \frac{{P}_{\mathrm{ram}}}{{m}_{{\rm{H}}}}\displaystyle \frac{1}{{N}_{{\rm{H}}}},\end{eqnarray} \tag{ 5 }$$

where ${{\rm{\Sigma }}}_{\mathrm{cl}}$ is the average surface density of the cloud, and ${N}_{{\rm{H}}}={n}_{\mathrm{cl}}L$ is the column density of a particular region of the cloud. If the cloud is a constant-density sphere, then L can be approximated as twice the cloud radius. For the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ spherical cloud with a radius of $r\approx 5$ pc, the column density is approximately ${N}_{{\rm{H}}}\approx 3\times {10}^{19}\,{\mathrm{cm}}^{-2}$ across the whole cloud, and the resulting acceleration is

$$\begin{eqnarray}&&{g}_{\mathrm{cl}}\sim 2.3\times {10}^{-6}\,\mathrm{cm}\,{{\rm{s}}}^{-2}{\left[\left(\displaystyle \frac{{n}_{\mathrm{cl}}}{1{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{L}{10\mathrm{pc}}\right)\right]}^{-1},\end{eqnarray} \tag{ 6 }$$

or

$$\begin{eqnarray}&&{g}_{\mathrm{cl}}\sim 0.75\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{kyr}}^{-1}{\left[\left(\displaystyle \frac{{n}_{\mathrm{cl}}}{1{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{L}{10\mathrm{pc}}\right)\right]}^{-1}.\end{eqnarray} \tag{ 7 }$$

The cloud crushing time for the $\tilde{n}=1$ clouds is t_cc $\approx \,56.4$ kyr. The acceleration in Equation (7) gives a cloud velocity of ${v}_{\mathrm{cl}}\sim 85\,\mathrm{km}\,{{\rm{s}}}^{-1}$ after 2t_cc (113 kyr), which is roughly consistent with our results for the velocities of the densest gas after 2t_cc, shown in the right panel of Figure 11. In fact, the average velocity of gas denser than n = 500 ${\mathrm{cm}}^{-3}$ at 2t_cc for the $\tilde{n}=1$ sphere is $77\,\mathrm{km}\,{{\rm{s}}}^{-1}$. After the cloud has been crushed, Equation (7) is no longer an adequate estimate of the cloud acceleration, because the column densities have increased by over an order of magnitude (compare the first and second panels in Figure 5).

In contrast, at t = 0 the densest sight lines through the $\tilde{n}=1\,{\mathrm{cm}}^{-3}$ turbulent cloud have column densities as high as ${N}_{{\rm{H}}}=2\times {10}^{20}\,{\mathrm{cm}}^{-2}$, almost an order of magnitude larger than the sphere (see the first panel in Figure 4). At these column densities, the cloud acceleration is only ${g}_{\mathrm{cl}}\approx 0.11\,\mathrm{km}\,{{\rm{s}}}^{-2}\,{\mathrm{kyr}}^{-1}$, yielding a dense gas velocity of ${v}_{\mathrm{cl}}\sim 12\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at t = 2t_cc. As a result, the regions of high column density are accelerated less in the turbulent cloud, and the high-density gas at t = 2t_cc is traveling considerably slower than in the spherical cloud case. The average velocity of gas with $n\gt 500$ ${\mathrm{cm}}^{-3}$ is only $v\approx 16\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at t = 2t_cc.

In both the turbulent cloud and sphere simulations, much less acceleration occurs at late times than predicted by Equation (7) (see right panel of Figure 10). As mentioned previously, this minimal acceleration can be explained by the increased column densities that result from the cloud crushing. In both the turbulent and spherical clouds, compression toward the center caused by the initial shock results in column densities at t = 2t_cc that are an order of magnitude higher than the original values (compare the first and second panels in both Figures 4 and 5). Thus, we conclude that entrainment of cool material to high velocities in a hot wind remains very inefficient under these circumstances, and the problem only compounds for clouds with more realistic shapes and density distributions.

These simulations also show that the average surface density of the cloud, ${{\rm{\Sigma }}}_{\mathrm{cl}}={M}_{\mathrm{cl}}/\pi {R}_{\mathrm{cl}}^{2}$, as given in Table 1 fails to adequately predict how much dense gas will be accelerated. This fact can be demonstrated by comparing the dense gas acceleration for the $\tilde{n}=0.5$ ${\mathrm{cm}}^{-3}$ turbulent cloud (not shown) and the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ sphere. The $\tilde{n}=0.5$ ${\mathrm{cm}}^{-3}$ turbulent cloud originally has a lower average surface density than the $\tilde{n}=1$ sphere—${{\rm{\Sigma }}}_{\mathrm{cl}}=0.11$ and ${{\rm{\Sigma }}}_{\mathrm{cl}}=0.13\,{M}_{\odot }$ ${\mathrm{pc}}^{-2}$, respectively. However, the spatial coherence of the densest regions within the turbulent cloud leads to individual column densities several times higher than the average column density of the sphere, around ${N}_{{\rm{H}}}\approx {10}^{20}\,{\mathrm{cm}}^{-2}$ as compared to ${N}_{{\rm{H}}}=3\times {10}^{19}\,{\mathrm{cm}}^{-2}$. As a result, the average velocity of the densest material in the $\tilde{n}=0.5$ ${\mathrm{cm}}^{-3}$ turbulent cloud simulation is ${v}_{x}\approx 30$ km ${{\rm{s}}}^{-1}$, less than half the speed of the dense gas in the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ spherical cloud simulation.

Having established that in our simulations momentum does not transfer efficiently enough to accelerate the dense phase to the wind velocity, we would like to better quantify the momentum gained by other phases of the outflow. Figure 6 showed the evolution of cloud mass for each of our high-resolution simulations as an integrated quantity, the total sum of material above a given density threshold. To better understand the way that gas evolves within the cloud, we now divide this evolution into multiple density bins. Figure 12 shows this binned mass evolution plot for the $\tilde{n}=1$ cloud. The black line matches the one displayed in Figure 6, comprising all the material above a density threshold of 1/3 $\tilde{n}$ but no longer normalized by the initial cloud mass. Colored lines in Figure 12 show the evolution of total cloud mass in four density bins: low, from $0.02\lt n\lt 0.2$ ${\mathrm{cm}}^{-3};$ medium, from $0.2\lt n\lt 2.0$ ${\mathrm{cm}}^{-3};$ high, from $2.0\lt n\lt 20$ ${\mathrm{cm}}^{-3}$, and very high, $n\gt 20$ ${\mathrm{cm}}^{-3}$. The other three panels show one-dimensional (1D) histograms of the mass as a function of density. Integrating any of the colored regions in the histogram yields the value represented as a circle on the line of the same color in the top left panel.

Zoom In Zoom Out Reset image size

Much of the evolution in Figure 12 takes place early on, with the most drastic shift occurring before 2t_cc as the cloud is crushed. Initially, much of the cloud mass (12 ${M}_{\odot }$) is in the high-density bin, $2.0\lt n\lt 20$ ${\mathrm{cm}}^{-3}$, with a significant amount (5.9 ${M}_{\odot }$) also in the medium-density bin that surrounds the initial median density of n = 1 ${\mathrm{cm}}^{-3}$. Only 0.64 ${M}_{\odot }$ is at $n\gt 20.0$ ${\mathrm{cm}}^{-3}$ initially. As the cloud is crushed, almost all of the material increases by about an order of magnitude in density, such that at 2t_cc the medium-density bin has only 0.85 ${M}_{\odot }$, the high-density bin has 6.0 ${M}_{\odot }$, and 11.6 ${M}_{\odot }$ is in the very high density bin. These values are highlighted by circles at 2t_cc in the first panel of Figure 12. Very little mass is in the low-density bin at this time, and no cloud mass has been lost.

After 2t_cc, the evolution proceeds more gradually, with material slowly moving from higher-density bins to lower-density ones. At 5t_cc mass is concentrated in the same locations evident in the density–temperature phase diagram—there is a significant amount of mass (10.0 ${M}_{\odot }$) at very high densities, as well as a bump around ${\mathrm{log}}_{10}(n)=0.5$ as a result of the shape of the cooling curve. These features are still present at 10t_cc, though the mass in the highest-density bin has been substantially reduced. At all times in our simulation after 2t_cc, the majority of the cloud mass is in the highest-density bin. Eventually we expect this mass to shift to lower-density bins as the last of the high-density gas is destroyed.

In the same way that we have integrated the total cloud mass, we can integrate total cloud momentum in the simulations. Figure 13 shows how the total momentum is divided among the same four density bins. The top left panel shows the evolution of the total momentum in the wind direction, computed as

$$\begin{eqnarray}&&{p}_{\mathrm{tot}}=\displaystyle \sum _{n={n}_{\mathrm{low}}}^{n={n}_{\mathrm{high}}}{M}_{i}{v}_{x,i}\end{eqnarray} \tag{ 8 }$$

for each cell, i, in the simulation, where M_i is the integrated mass in that cell and v_x is the direction of the wind. The sum is taken over the same density ranges shown in Figure 12. In the other three panels, we show histograms of the momentum as a function of density at t = 2t_cc, 5t_cc, and 10t_cc. Integrals of the histograms are again shown as circles on the colored lines in the top left panel. At t = 2t_cc the momentum is distributed in a similar way to the mass, with most of the momentum in the two highest-density bins. At t = 5t_cc, however, the distribution of total momentum has shifted. The three lower-density bins have continued to gain momentum, but the highest-density bin actually loses some, despite the fact that the highest-density bin still contains the majority of the cloud mass. At this point, the total momentum in the medium- and high-density bins ($0.2\,{\mathrm{cm}}^{-3}\lt n\lt 20\,{\mathrm{cm}}^{-3}$) is ${p}_{\mathrm{tot}}=495\,{M}_{\odot }$ km ${{\rm{s}}}^{-1}$, 50% more than is in the highest-density bin ($n\gt 20\,{\mathrm{cm}}^{-3}$), ${p}_{\mathrm{tot}}=315\,{M}_{\odot }$ km ${{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

This shift indicates that momentum is transferred efficiently from the hot wind to lower-density gas. Gas with densities between $1\,{\mathrm{cm}}^{-3}\lt n\lt 10\,{\mathrm{cm}}^{-3}$ is in a relatively long-lived phase, giving it more time to gain momentum from the hot wind. The total momentum in the warm phase likely also increases as gas from higher-density bins moves to lower density. Gas in the highest-density bin with the most momentum may be the most likely to decrease in density. Without tracer particles, this shift is difficult to quantify in our simulations, but the highest-density bin clearly loses mass at every time $t\gt 2$ t_cc and at least some high-density material with significant momentum shifts to the lower-density bins. Overall, the distribution of total momentum across the different gas phases in the cloud is remarkably equal. The density–velocity histograms shown in Figure 10 also indicate this equality, showing that lower-density material moves more quickly than higher-density material. At very late times, $t\gt 17$ t_cc, the highest-density bin again has the most total momentum—this reflects the lack of mass in the lower bins by this time.

We can also compare the distribution of momentum between the turbulent and spherical clouds by looking at the integrated momentum in the wind direction, a "momentum surface density":

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{p}=\displaystyle \sum _{i=0}^{i={N}_{x}}{n}_{i}{v}_{x,i},\end{eqnarray} \tag{ 9 }$$

where N_x is the total number of cells in the wind direction, n_i is the number density of the gas in cell i, and ${v}_{x,i}$ is the velocity in the wind direction in cell i. We show projections of this quantity for the $\tilde{n}=1$ turbulent cloud and sphere at t = 2t_cc in Figure 14. As one would expect given the differences in the density–velocity phase diagrams, the momentum surface density is much higher for the sphere, reaching a maximum of $2.4\times {10}^{23}$ km ${{\rm{s}}}^{-1}$ ${\mathrm{cm}}^{-2}$ in the highest column. The momentum surface density of the background wind is $3.1\times {10}^{21}$ km ${{\rm{s}}}^{-1}$ ${\mathrm{cm}}^{-2}$. In contrast, the momentum has a larger spatial extent and is spread over a wider range of gas column densities in the turbulent cloud, with the maximum column of $4.3\times {10}^{22}$ km ${{\rm{s}}}^{-1}$ ${\mathrm{cm}}^{-2}$. Thus, even though the average momentum surface density at this time is similar between the two simulations, the momentum is clearly much more concentrated in the high-density gas in the spherical cloud simulation. This result is consistent with the total momentum being shared across every density bin in Figure 13.

Zoom In Zoom Out Reset image size

Cooper et al. (2009) presented the only previous study that has investigated the destruction of clouds with a density distribution that is not symmetric along multiple axes. In that work, the authors compared the destruction of a radiatively cooling fractal cloud to that of a radiatively cooling spherical cloud. The background hot wind properties in Cooper et al. (2009) resembled ours. Our results agree well with theirs, in that they found that fractal clouds were more susceptible to fast destruction than spheres. However, their computational volume was too small to follow the evolution of the gas for many cloud crushing times, and their resolution was relatively poor, and thus the fate of the small dense clumps that result from the destruction of an inhomogeneous cloud remained unclear. In our work, we find that while these small cores can survive for tens of cloud crushing times as a result of efficient cooling, they are very difficult to accelerate due to their high column density. Additionally, higher resolution tends to amplify these effects. As a result, the cloudlets in our simulations are gradually destroyed over the course of $t\sim 20$t_cc as the dense gas gets eaten away by the hot wind.

In the past several years, several studies have investigated the ability of hot winds to entrain cool gas and carry it large distances from a galaxy. Scannapieco & Brüggen (2015) studied cloud–wind interactions using adaptive mesh refinement simulations with a maximum resolution equivalent to the fixed resolution of our simulations. The primary purpose of their work was to explore the effects that different background wind parameters had on cool cloud lifetimes and acceleration. Using spherical clouds of $T\approx {10}^{4}$ K, the authors derived a scaling for cloud destruction time that only depended on the Mach number of the hot wind, with different density contrasts accounted for in the cloud crushing time. Here, we compare our results to that scaling at t₅₀, defined as the time when 50% of the cloud material is at or above 1/3 the initial cloud density. Scannapieco & Brüggen (2015) find

$$\begin{eqnarray}&&{t}_{50}=4{t}_{\mathrm{cc}}\sqrt{1+{M}_{\mathrm{wind}}},\end{eqnarray} \tag{ 10 }$$

which for our background wind would correspond to ${t}_{50}=10$t_cc. Looking at Figure 6, we see that this actually fits the turbulent clouds quite well, but our spherical clouds have a much longer lifetime.

At first glance, this result may seem contradictory. However, the longer spherical cloud lifetimes may be explained by the different treatment of cooling in our simulations. As explained in Appendix C, we allow gas in our simulations to cool to temperatures as low as T = 10 K using cooling rates calculated for solar-metallicity gas. We also include the effects of a photoionizing background. The resulting heating keeps low-density gas ($n\lt 1$ ${\mathrm{cm}}^{-3}$) warm, but is much less effective at higher densities. In contrast, Scannapieco & Brüggen (2015) simulated larger-scale clouds with the assumption of complete ionization and therefore only allowed cooling above $T={10}^{4}$ K. The ability of gas to cool to temperatures well below $T={10}^{4}$ K in our simulations results in greater cloud compression. The resulting smaller surface area decreases the efficiency with which material ablates and correspondingly increases the cloud lifetime.

We find further evidence for this explanation by performing simulations with a different mean molecular weight μ. All the simulations presented in this paper used $\mu =1$ when converting from mass density to number density,

$$\begin{eqnarray}&&n=\rho /\mu {m}_{{\rm{p}}},\end{eqnarray} \tag{ 11 }$$

where m_p is the mass of a proton. However, we also performed a low-resolution turbulent cloud simulation with $\mu =0.6$, the value appropriate for fully ionized solar-metallicity gas. In comparing the lifetime of the cloud in this simulation with the same cloud in the standard simulation, we find that the $\mu =0.6$ cloud is destroyed more slowly, though the effect is small. This slight difference can be explained by the higher number densities for a given mass density in the $\mu =0.6$ simulation. These higher number densities lead to more efficient cooling, which in turn leads to higher average densities at early times. The resulting high-density clumps of gas prove more difficult for the hot wind to destroy than in the equivalent $\mu =1$ model. Hence, the $\mu =0.6$ cloud lives longer.

Scannapieco & Brüggen (2015) also find higher final velocities for their clouds than we do, even when comparing only spherical clouds. For clouds with similar background wind conditions they find velocities $v\sim 200\mbox{--}300$ km s⁻¹ at t = 15t_cc (see their Figure 7), while we never see velocities above 200 km s⁻¹. Again, this difference is likely a result of lower-temperature gas in our simulations. Because our clouds can compress to higher densities in the initial stages of the wind interaction, their column densities increase more and acceleration is more difficult. This effect is compounded in the turbulent cloud case with higher initial column densities.

Other studies of cool gas in the context of galactic winds have investigated the effects of additional physics in cloud–wind interaction simulations. Brüggen & Scannapieco (2016) perform a similar study to Scannapieco & Brüggen (2015) but incorporate the effects of conduction. They find that conduction can result in complete evaporation of the cloud at early times in cases where the column density of cloud material is below $N\simeq 1.5\times {10}^{18}$ ${\mathrm{cm}}^{-2}$. The clouds we simulate do not start with column densities this low, and we therefore do not expect conduction to result in their rapid, complete destruction. We note that the initial conditions of our hot wind most resemble their Mach 6.4 run that shows the least amount of difference in cloud evolution between the conduction and nonconduction models.

Brüggen & Scannapieco (2016) demonstrate that when clouds in their simulations do survive, conduction tends to decrease the cross section of the cloud presented to the hot wind. The decreased cross section provides less surface area for ram pressure to act on the dense gas, which in turn decreases the efficiency of cloud acceleration. This effect is qualitatively similar to the impact of inhomogeneous column densities for turbulent clouds described in Section 6.1, though the origin is completely different. We suggest that the inhomogeneous cloud structure may work in concert with conduction to further increase the early destruction of low column density material and decrease the efficiency with which high column material accelerates.

Magnetic fields may also play a role in cool cloud evolution. A number of cloud–shock studies investigated the effects of planar magnetic fields in spherical clouds, with inconclusive results regarding whether the presence of field lines increased or decreased the time until cloud destruction (Gregori et al. 2000; Fragile et al. 2005; Shin et al. 2008). More recently, McCourt et al. (2015) performed wind–cloud simulations testing the effects of tangled magnetic fields incorporated within a spherical cloud. They showed that the presence of the magnetic field drastically increased the lifetime of the resulting cloudlets and increased their acceleration, allowing them to reach the hot wind speed without destruction.

The McCourt et al. (2015) simulations help motivate a future extension of our high-resolution wind–cloud simulations with realistic initial conditions to include MHD. The simulations by McCourt et al. (2015) use a resolution of 32 cells/R_cl. Scannapieco & Brüggen (2015) demonstrated that spherical cloud simulations with resolution less than 64 cells per cloud radius tend to prolong cloud lifetimes. Our simulations of turbulent clouds show a trend toward faster mass loss with increasing resolution that has not yet converged at 128 cells/R_cl. In addition, we have shown that spherical clouds have longer lifetimes than those with more realistic internal structure. While the results of McCourt et al. (2015) suggest that including a tangled magnetic field would likely prolong the lifetime of the dense components of the turbulent clouds we simulate, the combined effects of resolution and density structure make a quantitative prediction difficult. Incorporating magnetic fields in wind simulations with turbulent clouds is an avenue that we aim to pursue in future work.

In addition to potentially important additional physics, parameters such as the metallicity of the cloud gas and nature of the background wind in our simulations will affect the quantitative results we have presented. In this work we have focused exclusively on the effects of cloud structure and surface density. The primary effect of changing the metallicity of the gas would be to change the cooling rates. As noted in Section 8.2, tests with a different value of μ indicate that more rapid cooling leads to longer cloud survival, and our comparison with the simulations of Scannapieco & Brüggen (2015) indicates that less efficient cooling reduces cloud survival time. However, this is not an order-of-magnitude effect, and we therefore do not expect a change in metallicity to drastically alter our results.

On the other hand, as Figure 1 shows, the state of the background wind in the Chevalier & Clegg (1985) model changes rapidly with radius as the wind escapes the galaxy. Given our choice of a single background wind state, we do not regard our simulations as completely generic with respect to galactic winds, though we do expect the general result of more rapid destruction of turbulent clouds to hold. Scannapieco & Brüggen (2015) demonstrated a scaling with background wind Mach number that indicates that clouds survive longer with increasing Mach number (see their Equation (22)). In the future, we would like to investigate the relationship between turbulent cloud lifetime and the background wind parameters, sampling a variety of distances from the galaxy that would inform the initial conditions for clouds at each distance.

As a final note, we consider the final fate of the dense gas in our simulations. Our simulations do not include gravity, and if the simulations continued running indefinitely, eventually all of the cool material would mix into the outflowing hot wind. However, we can estimate the effect of the host galaxy's gravitational potential. Using an analysis similar to that in Section 6.1, we can compare the expected acceleration of dense cloud regions owing to the wind's ram pressure, P_ram, to their expected deceleration owing to gravity.

We can estimate the ram pressure acceleration as a function of column density, ${N}_{{\rm{H}}}={n}_{\mathrm{cl}}L$, via

$$\begin{eqnarray*}&&{g}_{\mathrm{ram}}\sim 2.3\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{kyr}}^{-1}{\left(\displaystyle \frac{{N}_{{\rm{H}}}}{{10}^{19}{\mathrm{cm}}^{-2}}\right)}^{-1}.\end{eqnarray*}$$

(This expression is equivalent to Equation (7).) Similarly, we can estimate the gravitational acceleration as a function of column density. At 1 kpc, the gravitational acceleration of M82 is

$$\begin{eqnarray}{g}_{\mathrm{grav}} & \sim & 1.4\times {10}^{-7}\,\mathrm{cm}\,{{\rm{s}}}^{-2}\left(\displaystyle \frac{{M}_{{\rm{M}}82}}{{10}^{10}\,{M}_{\odot }}\right){\left(\displaystyle \frac{R}{1\mathrm{kpc}}\right)}^{-2},\\ {g}_{\mathrm{grav}} & \sim & 0.044\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{kyr}}^{-1}.\end{eqnarray} \tag{ 12 }$$

Given these estimates, we would expect cloud regions with column densities greater than ${N}_{{\rm{H}}}\approx 5\times {10}^{20}$ ${\mathrm{cm}}^{-2}$ to begin to fall back toward the galaxy. None of the clouds in our simulations initially have column densities quite this high—the densest sight lines in the $\tilde{n}=1$ ${\mathrm{cm}}^{-3}$ turbulent cloud reach ${N}_{{\rm{H}}}\approx 3\times {10}^{20}$ ${\mathrm{cm}}^{-2}$. Nonetheless, the acceleration due to ram pressure and deceleration due to gravity for these column densities are very similar, so in the presence of gravity the densest gas in our turbulent cloud simulations would likely be accelerated very little or possibly fall back toward the central galaxy.

We implement radiative cooling in Cholla using an operator-split approach. After the hydrodynamic quantities have been updated for a given time step according to the ideal hydrodynamic equations, we add a source term to the internal energy to account for losses (or gains) as a result of radiative cooling (or heating) of the gas, such that

$$\begin{eqnarray}&&{e}^{n+1}={\tilde{e}}^{n}+\dot{e}{\rm{\Delta }}{t}_{\mathrm{ad}},\end{eqnarray} \tag{ 15 }$$

where $\dot{e}$ is the rate of change of the internal energy and ${\rm{\Delta }}{t}_{\mathrm{ad}}$ is the hydrodynamic time step. Here, ${\tilde{e}}^{n}$ represents the updated internal energy after the hydrodynamic time step. The internal energy is either calculated from the gas pressure assuming an ideal gas equation of state or tracked directly via the dual-energy formalism. Because the adiabatic time step can be large compared to the radiative cooling time $e/\dot{e}$, the rate of change in internal energy is often a nonlinear function of the temperature over the course of a single ${\rm{\Delta }}{t}_{\mathrm{ad}}$. Thus, we employ the common approach of subcycling the radiative cooling steps (e.g., Smith et al. 2008; Gray & Scannapieco 2010), calculating a new $\dot{e}$ many times over the course of a single adiabatic time step so as to limit the change in internal energy for a radiative substep, ${\rm{\Delta }}{t}_{\mathrm{rad}}$, to less than 5% of the current internal energy:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}e}{e}\lt 0.05.\end{eqnarray} \tag{ 16 }$$

In practice, this update means that at the end of each hydrodynamic time step, we calculate the temperature for a given cell's number density and internal energy according to the ideal gas law

$$\begin{eqnarray}&&T=\displaystyle \frac{P}{{nk}},\end{eqnarray} \tag{ 17 }$$

where n is the number density of the gas in ${\mathrm{cm}}^{-3}$ calculated assuming $n=\rho /(\mu {m}_{{\rm{p}}})$, m_p is the mass of a proton, and k is Boltzmann's constant. We have assumed a mean molecular weight of $\mu =1$ for all calculations and have verified that this does not influence our results compared with values of $\mu =0.6$, as described in Section 8.

We next look up the tabulated net cooling rate associated with this density and temperature using a bilinear interpolation (described below) and calculate the resulting change in temperature over the full adiabatic time step given the cooling rate, ${\rm{\Delta }}T$ = $\dot{T}$ ${\rm{\Delta }}{t}_{\mathrm{ad}}$. If the change in temperature is greater than 5%, we shrink the radiative substep such that ${\rm{\Delta }}{t}_{\mathrm{rad}}$ results in a temperature change of 5%. We then update the temperature with this smaller time step and repeat the process until we have synchronized with the full adiabatic time step.

We tabulate the cooling and heating rates per unit volume, Λ and Γ (measured in erg ${{\rm{s}}}^{-1}$ ${\mathrm{cm}}^{-3}$), using the Cloudy code, version 13.03 (Ferland et al. 2013). For the calculations presented in this work, we assume that the gas is in photoionization equilibrium, subject to the z = 0 cosmic microwave background and the HM05 UV/X-ray background, as described in Hazy, the Cloudy documentation. The gas metallicity is assumed to be solar, using the GASS10 abundances in Cloudy (Grevesse et al. 2010). Examples of the absolute cooling rates for gas of several different densities are shown in Figure 20.

Zoom In Zoom Out Reset image size

Our cooling and heating rates are tabulated on a grid, with points calculated at equally spaced logarithmic intervals of 0.1 between 10 K $\lt T\lt {10}^{9}$ K and 10⁻⁶ ${\mathrm{cm}}^{-3}\lt n\lt {10}^{6}$ ${\mathrm{cm}}^{-3}$. To perform the interpolation necessary to calculate Λ and Γ at an arbitrary n and T, we take advantage of a novel feature of GPUs—texture memory. This special memory space allows a user to copy a 1D, 2D, or 3D array, or "texture," onto the GPU, and then use built-in GPU functions to quickly retrieve arbitrary values from that texture using bilinear interpolation. We have tested this method against a similar CPU-based method using GSL bilinear interpolation functions and find that the GPU texture memory approach typically speeds up the radiative cooling calculation by orders of magnitude as compared to a CPU function. In fact, when implemented using the operator-splitting approach described above for simulations run with ∼128³ cells/GPU, the time spent calculating radiative cooling is completely negligible compared to the hydrodynamic calculations.