The Production of Cold Gas Within Galaxy Outflows

TL;DR

This paper uses 3D simulations to study how hot gas in galaxy outflows cools and produces observable cold gas, revealing the dominant cooling processes and their impact on emission properties.

Contribution

It introduces a detailed simulation framework including various cooling mechanisms and constructs a steady-state model linking outflow properties to observed emissions in ULIRGs.

Findings

Compton cooling is the most efficient cooling process.

Atomic/ionic cooling enhances density inhomogeneities.

Model luminosities match observed emission-line properties in ULIRGs.

Abstract

I present a suite of three-dimensional simulations of the evolution of initially-hot material ejected by starburst-driven galaxy outflows. The simulations are conducted in a comoving frame that moves with the material, tracking atomic/ionic cooling, Compton cooling, and dust cooling and destruction. Compton cooling is most efficient of these processes, while the main role of atomic/ionic cooling is to enhance density inhomogeneities. Dust, on the other hand, has little effect on the outflow evolution, and is rapidly destroyed in all the simulations except the case with the smallest mass flux. I use the results to construct a simple steady-state model of the observed UV/optical emission from each outflow. The velocity profiles in this case are dominated by geometric effects, and the overall luminosities are extremely strong functions of the properties of the host system, as observed in ultra-luminous infrared galaxies (ULIRGs). Furthermore the luminosities and maximum velocities in several models are consistent with emission-line observations of ULIRGs, although the velocities are significantly greater than observed in absorption-line studies. It may be that absorption line observations of galaxy outflows probe entrained cold material at small radii, while emission-line observations probe cold material condensing from the initially hot medium at larger distances.

Figures

Tables

Table 1. Table 1 Simulation Parameters

Name	SFR	${\dot{M}}_{\mathrm{hot}}$	$ϵ$	$v_{\star }$	$v_{\mathrm{final}}$	$T_{\star }$	$n_{\star }$	Atomic	Compton	Dust
	$M_{\odot }$ yr-1	$M_{\odot }$ yr-1		km s-1	km s-1	${10}^6$ K	cm-3	Cooling	Cooling	Cooling
M82_10_1	10	10	1.0	500	1000	18	4.7	N	N	N
M82_10_1A	10	10	1.0	500	970	18	4.7	Y	N	N
M82_10_1AC	10	10	1.0	500	920	18	4.7	Y	Y	N
M82_10_1ACD	10	10	1.0	500	900	18	4.7	Y	Y	Y
M82_5_1ACD	10	5	1.0	710	1320	36	1.7	Y	Y	Y
M82_2_1ACD	10	2	1.0	1120	2120	90	0.42	Y	Y	Y
M82_1_1ACD	10	1	1.0	1580	3020	180	0.15	Y	N	Y
M82_5_0.5ACD	10	5	0.5	500	910	18	2.4	Y	Y	Y
M82_2_0.5ACD	10	2	0.5	790	1480	45	0.6	Y	Y	Y
ULIRG_50_1.0ACD	100	50	1	710	1220	36	17	Y	Y	Y
ULIRG_20_1.0ACD	100	20	1	1120	1950	90	4.2	Y	Y	Y
ULIRG_50_0.5ACD	100	50	0.5	500	860	18	24	Y	Y	Y
ULIRG_20_0.5ACD	100	20	0.5	790	1360	45	6.0	Y	Y	Y

Equations

$$\frac{1}{r^2}\frac{d}{dr}(\rho v{r}^2)={\dot{q}}_{\mathrm{m}},$$

$$\rho v\frac{dv}{dr}=-\frac{dP}{dr}-{\dot{q}}_{\mathrm{m}}v,$$

$$\frac{1}{r^2}\frac{d}{dr}\left[\rho v{r}^2\left(\frac{v^2}{2}+\frac{\gamma }{\gamma -1}\frac{P}{\rho }\right)\right]={\dot{q}}_{\mathrm{e}},$$

$${\dot{q}}_{\mathrm{m}}$$

$${\dot{q}}_{\mathrm{e}}$$

$$c_{s,\mathrm{hot}}(r)\approx c_{s,\mathrm{hot},0}=0.82{\left(\frac{ϵ\dot{E}}{{\dot{M}}_{\mathrm{hot}}}\right)}^{1/2},$$

$$\frac{d{E}_k}{dt}$$

$$\frac{dI}{dt}$$

$$\rho ={\rho }_{\star }(v/v_{\star }{)}^{-2}(r/r_{\star }{)}^{-1},$$

$$P_j\equiv c\int {\sigma }_j(\nu )n_{\nu }d\nu ,$$

$$P_j\equiv {\tilde{P}}_j\frac{SFR}{M_{\odot }{\mathrm{yr}}^{-1}}{\left(\frac{r}{300\mathrm{pc}}\right)}^{-2},$$

$$H_C=\frac{{\sigma }_{T}F}{m_e{c}^2}(⟨h\nu ⟩-4k_{B}T)n_e={\Gamma }_C\Delta ϵn_e,$$

$$m_g=4\pi a^3{\delta }_g/3.$$

$$H_{\mathrm{coll}}=\{\begin{array}{ll}5.38\times 10^{-10}{n}_e{a}^2{T}^{3/2}& x^{\ast }>4.5\\ 1.47\times 10^{-3}{n}_e{a}^{2.41}{T}^{0.88}& x^{\ast }>1.5\\ 6.48\times 10^6{n}_e{a}^3& x^{\ast }\le 1.5,\end{array}$$

$$G_0=1650\frac{\dot{M}}{M_{\odot }{\mathrm{yr}}^{-1}}{\left(\frac{R}{300\mathrm{pc}}\right)}^{-2}.$$

$$C_d^{\mathrm{rec}}\approx 4.65\times 10^{-30}\frac{\mathrm{erg}}{\mathrm{s}{\mathrm{cm}}^3}\left(\frac{\mathcal{D}}{{\mathcal{D}}_{\odot }}\right)T^{0.94}{\left(\frac{G_0\sqrt{T}}{n}\right)}^{\delta }n{n}_e,$$

$$\frac{da}{dt}=-An{T}_6^{-1/4}{e}^{-B{T}_6^{-1/2}},$$

$$\frac{d{m}_g}{dt}=4\pi a^2{\delta }_g\frac{da}{dt},$$

$$n(\mathbf{x})=n_{\star }+\frac{3n_{\star }}{8}\left[\frac{r_4+0.5r_2+0.25r_1}{1+0.5/8^{1/2}+0.25/8}\right],$$

$$\frac{r}{r_{\star }}$$

$$v_{\mathrm{clump}}=x_{\mathrm{clump}}\frac{da}{dt}\approx x_{\mathrm{clump}}\frac{2}{3}\frac{v}{r}a<c_s,$$

$$x_{\mathrm{clump},\mathrm{com}}\text{lower}2.15277pt\frac{\text{buildrel}<}{\sim }\frac{3r}{2M{a}^2}\approx r_{\star }\left[\frac{1}{M}\frac{3}{2^{7/3}}\frac{r_{\star }^{1/3}}{r^{1/3}}\right]\approx \frac{r_{\star }^2}{4r},$$

$${\left(\frac{dI}{dt}\right)}_{\mathrm{adiabatic}}=I(\gamma -1)\frac{d\ln n}{dt}.$$

$$P(k)\equiv ⟨\tilde{n}(k){\tilde{n}}^{\ast }(k)⟩,$$

$$\frac{⟨n^2⟩}{{⟨n⟩}^2}-1=\frac{4\pi }{{⟨n⟩}^2}\int {\left(\frac{k}{2\pi }\right)}^{3}P(k)\mathrm{d}\ln k.$$

$$\frac{d{L}_i}{d{v}_z}$$

Full text

1. Introduction

Many of the observed features of galaxies are thought to be caused by outflows. The material they transport helps to determine the number density and baryonic content of dwarf galaxies (e.g. Cole et al., 2000; Scannapieco et al., 2001, 2002; Benson et al., 2003); the metals they remove helps to determine the relationship between the interstellar medium metallicity and galaxy stellar mass, (e.g. Tremonti et al., 2004; Kewley & Ellison, 2008; Mannucci et al., 2010); and the enriched material they disperse helps to determine the metallicity history of the intergalactic medium (e.g. Ferrara et al., 2005; Martin et al., 2010; Peeples et al., 2014; Turner et al., 2015). Yet despite their importance, the dynamical and microphysical processes that control galaxy outflows remain extremely poorly understood.

From a theoretical perspective, the most important outstanding issue is that the highly efficient cooling within the interstellar medium (ISM) makes it impossible to model supernovae by simply adding thermal energy to the medium, unless both the Sedov and shell formation stages are sufficiently well resolved (Kim & Ostriker, 2015). At the same time, currently-available computational resources do not allow for direct modeling of supernovae within a galaxy-scale simulation (e.g. de Avillez & Breitschwerdt, 2004; Cen & Chisari, 2011; Barai et al., 2013; Creasey et al., 2015; Walch et al., 2015). The connection to observations is also complicated by the possibility that radiation pressure on dust (Thompson et al., 2005; Murray et al., 2011; Hopkins et al., 2011), non-thermal pressure from cosmic rays (Socrates et al., 2008; Uhlig et al., 2012; Booth et al., 2013; Salem & Bryan, 2014), and energy input from gravitationally driven motions (Sur et al., 2016) may also contribute to driving winds.

From an observational perspective, the most important issue is that the multiphase nature of galaxy outflows causes us to only have a partial view of the material they eject. In fact the outflowing material spans an enormous range of temperatures: from $10^{7}-10^{8}$K plasma observed in X-rays (e.g. Martin, 1999; Strickland & Heckman, 2007, 2009), to $\approx 10^{4}$K ions observed at optical and near UV wavelengths (e.g. Pettini et al., 2001; Tremonti et al., 2007; Martin et al., 2012; Soto et al., 2012), to $10-10^{3}$ K molecular gas observed at radio wavelengths (e.g. Walter et al., 2002; Sturm et al., 2011; Bolatto et al., 2013). Furthermore, the easiest phase to interpret, the X-ray emitting gas, is detectable only in the nearby Universe. The colder phases, on the other hand, can be easily observed out to high redshifts, and their kinematics have been carefully studied in many cases (e.g. Schwartz & Martin, 2004; Rupke et al., 2005; Erb et al., 2012; Weiner et al., 2009; Kornei et al., 2013; Arrigoni Battaia et al., 2015). On the other hand, not only is the connection between the motion of this material and the X-ray emitting gas highly uncertain, but its overall origin remains unknown.

Supernovae are not evenly distributed in time and space, and so complete thermalization of the outflowing medium may not be achieved, in particular when the star formation rate is small and supernova events not so frequent. Incomplete thermalization would cause the outflow to be born as a collection of hot gas surrounding cold clouds that would then be driven out of the host galaxy by ram pressure acceleration from hotter material (e.g. Veilleux et al., 2005). However, this model runs into serious difficulties in explaining the observations, both because: (i) shocks and conduction from the exterior medium tend to compress the cloud perpendicular to the direction of the flow, greatly reducing the momentum flux it receives; and (ii) instabilities and evaporation lead to rapid cloud disruption (Klein et al., 1994; Mac Low & Zahnle, 1994; Orlando et al., 2006, 2008; Scannapieco & Brüggen, 2015; Brüggen & Scannapieco, 2016). Together these imply that the lifetimes of the clumps are likely to be much shorter that the time required to accelerate them to the observed speeds.

A second possibility is that thermalization of the outflowing medium is efficient and the cold clouds are instead formed later, directly from the cooling of high temperature ($10^{7}-10^{8}$K) material, already moving at high radial velocities. This was first proposed by Wang (1995a, b), who calculated steady-state one-dimensional profiles of outflowing winds including cooling, parameterized as a power law. A more detailed one-dimensional model was then developed by Efstathiou (2000), which incorporated the infall of cooling gas from a halo, and the outflow of hot gas from the interstellar medium. Silich et al. (2003, 2004), Tenorio-Tagle et al. (2007), and Wünsch et al. (2011) presented steady-state one and two-dimensional hydrodynamic calculations of galaxy outflows including radiative cooling, accounting for chemical evolution, and comparing their results with observations of several nearby galaxies. Martin et al. (2015) studied simple models of cooling winds in the context of the Ly$\alpha$ line profiles of ultra-luminous infrared galaxies (ULIRGs). Most recently, Thompson et al. (2016) derived one-dimensional models of wind evolution over a large parameter space, including a radiative cooling function that captured the key role of metal-line cooling at temperatures below $10^{7}$ K (see also Zhang et al., 2015). They also provided scalings for the radius at which gas cools efficiently as a function density, column density, emission measure, radiative efficiency, and cool gas velocity, highlighting the dependence of the results on the mass entrained by the wind and the host galaxy star formation rate surface density.

Here I present three-dimensional simulations of the evolution of cooling gas in galaxy outflows that include several physical processes beyond those considered in previous studies. To maximize the resolution of the calculations, I adopt a comoving frame that moves radially outward with the outflowing material, capturing both the global cooling of the flow, and the formation of clumps within it. The cooling routines make use of the code capabilities developed in Ferrara & Scannapieco (2016) to study quasar-driven outflows, including: (i) radiative line cooling that fully accounts for the impact of the ionizing background and the additional heating due to phoionization, (ii) cooling due to the Compton scattering between free elections and photons from the host galaxy, which can exceed metal-line cooling by over an order of magnitude at small radii, and (iii) cooling from dust grains, calculated by self-consistently tracking the destruction of dust by sputtering by the hot outflow material.

The structure of this paper is as follows. In §2, I describe the simulation methods, focusing on the approach used to work in a comoving frame, the modeling of atomic and ionic cooling and Compton cooling, and the modeling of dust cooling and destruction. In §3, I present the simulation results, describing first the impact of the various physical processes on a fiducial case, and then contrasting this evolution with outflows with varying energy input rates, mass loading factors, and star formation rates. In this section, I also examine the size distribution of the cold clumps produced in the simulations, and in §4, I use these results to assess the observability of the cold gas in each of the cases studied. Conclusions are given in §5.

2. Methods

To study the formation of cold clouds in starburst-driven outflows, I carried out a set of simulations using FLASH (version 4.2), a multidimensional hydrodynamic code (Fryxell et al., 2000). The simulations used an unsplit hydrodynamic solver with a predictor-corrector type formulation (Lee, 2013), and they made use of the shock-detect flag which lowered the prefactor in the Courant-Friedrichs-Lewy timestep condition from 0.4 to 0.25 in the presence of strong shocks. The simulations included atomic / ionic cooling and heating process, Compton cooling, and dust cooling and destruction using the methods described in detail in Ferrara & Scannapieco (2016), modified to handle the current problem as described below. Unlike the simulations in Ferrara & Scannapieco (2016) however, I also modified the FLASH code follow a cubic periodic volume in a comoving frame that allows the medium to expand, conserving mass, as it moves outward from the galaxy. Here I describe each of the physical processes included in the calculations, beginning with this new approach.

2.1. Comoving Frame

We are interested in gas condensing within a thermalized hot wind expanding rapidly from a starbursting galaxy. For most galaxies, this hot medium moves many scale heights per Myr, while a typical starburst episode lasts for many Myrs (e.g. Greggio et al., 1998; Förster Schreiber et al., 2003). Thus, the distribution is expected to be well-approximated by an equilibrium configuration, and this seems to be born out observationally in most objects in which reliable X-ray analyses can be made (e.g. Heckman et al., 1990, 1995; Ott et al., 2005; Strickland & Heckman, 2007; Yukita et al., 2012). In this case, assuming a conical expansion with a fixed opening angle, the equations of mass, momentum, and energy conservation are:

[TABLE]

[TABLE]

[TABLE]

where $\rho$, $v$, $P,$ and $\gamma$ are the density, radial velocity, pressure, and the ratio of specific heats of the outflow and the mass and energy input rates are

[TABLE]

respectively, where $r_{\star}$ is the driving radius of the flow, $c_{s,0}$ is the sound speed of the hot medium at $r=0,$ $\Lambda$ is the radiative cooling function, and $n=\rho/\mu m_{p}$ is the total number density of the medium, with $\mu$ the average particle mass in units of the proton mass, $m_{p}$. For reference, in the case of M82, $r_{\star}\approx 300$ pc (McKeith et al., 1995; Strickland & Heckman, 2009). Note that here I have ignored radiative cooling within the driving radius, including this in the definition of $c_{\rm s,0}$.

The solution to these equations is given by the Chevalier & Clegg (1985) model, which can be used to gain a good understanding of the hot wind conditions as a function of radius. For a $\gamma=5/3$ gas, the sound speed when $r<r_{\star}$ is approximately constant at

[TABLE]

where $\dot{E}$ is the total energy input from supernovae per unit time, $\epsilon$ is the fraction of this energy that is deposited into the hot medium, $\dot{M}_{\rm hot}$ is the total mass ejected from the galaxy per unit time.

We can solve eqs. (1)-(3) outside the driving radius by simplifying the mass and momentum equations, using the fact that $\rho vr^{2}$ is a constant at these radii. Defining the kinetic and internal energy per unit mass as $E_{k}=v^{2}/2$, and $I=(\gamma-1)c_{s}^{2}$, this gives

[TABLE]

where $\dot{Q}\equiv-\Lambda\,n_{b}^{2}/\rho$ is the rate of radiative energy loss per unit mass, $n_{b}$ is the number density of baryons (protons + neutrons) and $\Lambda$ is the cooling - heating rate as discussed in more detail below. These equations can be used to track $v$ and $r$ as a function of time from the edge of the driving region for any choice of $\dot{Q}.$ The mean density of the medium as a function of time is then given by

[TABLE]

where $v_{\star}$ and $\rho_{\star}$ are the velocity and density of the medium at the driving radius. This approach allows the code to work in a comoving volume that travels outward with the medium, conserving mass. I have modified the cosmology routines in the FLASH code to allow for such a calculation, integrating eqs. (6) and (7) forward using a fourth-order Runge-Kutta method, and computing the expansion factor $a(t)=\left[\rho(t)/\rho_{0}\right]^{-1/3}.$ Note that unlike usual cosmological calculations in which $a(t)=1$ at the end of the calculation, $a(t)=1$ at the onset of the calculation when $r=r_{\star},$ and is greater than 1 at all subsequent times.

2.2. Atomic, Ionic, and Compton Cooling and Heating

In order to track atomic/ionic cooling processes in the presence of significant photoionization from the central starburst, I followed the approach in Ferrara & Scannapieco (2016) and made use of the cooling/heating functions developed in Gnedin & Hollon (2012) (see also Sazonov et al., 2005). This requires a calculation of

[TABLE]

where $n_{\nu}$ is the radiation field expressed as a number density of photons at a frequency $\nu$ and $\sigma_{j}$ is the dissociation/photoionization cross section for $j\in$ $\{$H2, HI, HeI, CIV$\}$. For the starburst spectrum, I took a case in which the galaxy has been forming solar metallicity stars at a constant rate as modeled by Schaerer (2003). Defining

[TABLE]

where SFR is the star formation rate. Making use of the analytic fits provided in Verner et al. (1996), this gives $\tilde{P}_{j}$ of $1.03\times 10^{-7}$ s*-1*, $5.28\times 10^{-8}$ s*-1*, $2.05\times 10^{-8}$ s*-1* and 0 for Lyman-Werner dissociation, and HI, HeI, and CVI ionization, respectively.

In the left panel of Figure 1, I show the cooling and heating rate normalized by the mean baryonic number density, $n_{B}^{2},$ computed for a range of number densities from $0.01$ to $100$ cm*-3* for a medium at a distance of 300 parsecs from a starburst with $\dot{M}_{\rm hot}=10M_{\odot}\,{\rm yr}^{-1}$. At the very highest densities, recombination times are short and the cooling rate above $10^{4}$ K is similar to the case without a background, although heating by photoionization is important below $10^{4}$K. At 1 cm*-3*, heating is boosted only slightly, but the lowest temperature peak of the cooling function is significantly lessened, due primarily to hydrogen remaining much more ionized than would be the case in the absence of a background. Finally, at the lowest densities, the mass fractions of neutral hydrogen, neutral helium, and other species with low-ionization potentials are extremely small, leading to a large decrease in the cooling rate at temperatures below $\approx 3\times 10^{5}$K. Together these effects shift the point at which cooling and heating are balanced from $\approx 10^{4}$K for gas with \rho/\mu m_{p}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};} 10 cm*-3* to $\approx 3\times 10^{4}$K for gas with \rho/\mu m_{p}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};} 0.1 cm${}^{-3}.$

In addition to atomic/ionic cooling as tabulated in these models, the presence of radiation from the starburst also leads to additional gas cooling and heating through Compton scattering with free elections. To estimate this contribution, I modified the overall cooling/heating rate of the gas by a factor of

[TABLE]

where $F=L/(4\pi r^{2})$, $\langle h\nu\rangle=L^{-1}\int h\nu L_{\nu}d\nu$ = 9.14 eV such that $T_{C}\equiv\langle h\nu\rangle/4k_{B}=2.65\times 10^{4}$K, $L=5.55\times 10^{43}$ erg s*-1* SFR/ ($M_{\odot}$ yr*-1*) is the luminosity per unit star formation rate (SFR). Including this effect shifts the cooling and heating curves from their Gnedin & Hollon (2012) values to those shown in the right panels of Figure 1, for a case with star formation rate of 10 $M_{\odot}$ yr*-1* and a distance of 300 pc. Because the energy losses from Compton scattering are proportional to $n_{e}$ rather than density squared, they have the largest impact at the lowest densities, and because these losses are proportional to $T-T_{C}$, these have the largest effect when this gas is the hottest. Thus Compton cooling is larger than atomic/ionic cooling for all temperatures above $T\approx 10^{6}$K at densities lower than \rho/\mu m_{p}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};} 10 cm${}^{-3}\,(\dot{M}_{\rm hot}/10M_{\odot}{\rm yr}^{-1})(r/300{\rm pc})^{-2}.$

2.3. Dust Cooling and Destruction

As in Ferrara & Scannapieco (2016), the simulations include the impact of dust, which is able to provide cooling through two mechanisms. In the first mechanism, ion-grain and electron-grain collisions transfer thermal energy to the dust, which is subsequently radiated away in the infrared. The cooling rate due to this mechanism can be written as $C_{d}^{\rm rad}=n_{g}H_{\rm coll}$ (Draine, 1981) and Dwek & Werner (1981) where $n_{g}={\cal D}_{0}n(\mu m_{p}/m_{g,0})$ is the number density of dust particles in the gas, with $m_{g,0}$ the initial grain mass, calculated by approximating the grains as sphere of radius $a$ and constant density $\delta_{g},$

[TABLE]

The heating rate deposited by particles into the grain (in ergs) can be parametrized as

[TABLE]

where $x^{*}=1.26\times 10^{11}a^{2/3}/T$ (Montier & Giard 2004).

The second mechanism acts when a strong UV field is present, leading to positively-charged grains. In this case the thermal energy of electrons ejected from the grain surface is less than the energy of electrons recombining, leading to a net cooling rate that depends on the flux in the UV (“Habing”) band 6-13.6 eV, $G_{0}=F/F_{H}$, where the source flux is normalized to the Habing flux $=1.6\times 10^{-3}\,{\rm erg\,{s^{-1}}}$ cm*-2* (Habing, 1968), such that

[TABLE]

Then, the cooling rate due to recombinations on dust for a solar metallicity gas with a Milky Way dust-to-gas ratio is

[TABLE]

with $\delta=0.74/T^{0.068}$.

Dust is included as a separate fluid, following its evolution using a passive tracer field that travels with the flow (Ferrara & Scannapieco, 2016). In each cell, ${\cal D}={\cal D}_{0}$ initially and dust is modeled as a single-grain population, whose pre-shock radius is $\langle a\rangle=0.1\mu$m. As the grains encounter the hot medium, they suffer a gradual erosion due to thermal sputtering, described by a simple fit to the Dwek et al. (1996) results (see also Draine & Salpeter, 1979):

[TABLE]

where $n$ and $T_{6}$ are the density and temperature (in units of $10^{6}K$), and $(A,B)=(1.2\times 10^{-5}\mu$m yr${}^{-1},3.85)$. Finally, the rate of decrease of dust mass of a single grain is given by

[TABLE]

where the grain density $\delta_{g}=3$ g $\,{\rm{cm^{-3}}}$.

3. Results

Having implemented these code modifications, I carried out a set of simulations that span a large range of star formation rates (SFR), mass outflow rates ($\dot{M}_{\rm hot}$), and energy input rates $\dot{E},$ as summarized in Table 1. All cases adoped a driving region with a radius of $r_{\star}=300$ pc, a value consistent with the nearby starburst M82 (McKeith et al., 1995; Strickland & Heckman, 2009) and also used in Zhang et al. (2015). The initial position of the simulation region was taken to be at $r_{\star}$ and the simulations tracked material within a 5123 cell comoving box, $150$ parsecs on a side, such that the comoving resolution was $0.29$ parsec. All simulations assumed that 1051 ergs of energy were released per 100 solar masses of stars formed, and that a fraction $\epsilon$ of the energy went into the wind. As described in §2.1, these values then determined the initial radial velocity and temperature of the material, $v_{\star}$ and $T_{\star}$, and the subsequent evolution of the flow determines the final value $v_{\rm final}$, which is 2 $v_{\star}$ in the absence of cooling. Each of these values are also shown in Table 1. All simulations also assumed a total opening angle of outflowing material of $\Omega=\pi$ steradians as observed in large statistical samples of low-redshift starbursts (Martin et al., 2012), and this angle in turn determined the initial density of the material, as also given in the table.

As they are approximating an initially thermalized medium, the simulations were run with a uniform initial temperature, but a fluctuating initial density field given by a white noise spectrum down to the scale of $1.2$ comoving parsecs (4 cells), and the amplitude of the fluctuations dropping by 50% at $0.58$ parsecs and by 25% at $0.29$ parsecs. The fluctuations were normalized to have an rms of $n_{\star}/8$ on 4 parsec scales, such that

[TABLE]

where are $r_{1}$, $r_{2}$, and $r_{4}$ are random numbers selected from a uniform distribution between -1 and 1, which vary on scales of 1, 2, and 4 cells respectively.

The simulations were of two types, a set with a SFR of 10 $M_{\odot}$ yr*-1*, similar to that of the nearby starburst M82 (Förster Schreiber et al., 2003) and second set with a SFR of 100 $M_{\odot}$ yr*-1* similar to what is observed in ULIRGs (Daddi et al., 2005; Piqueras López et al., 2012). For each type of simulation, I varied the mass and energy input into the wind, and in the M82 case I also varied the physical assumptions, helping to disentangle the role that the various physical processes described in §2 have on the evolution of the outflow.

3.1. Global Evolution and the Impact of Dust and Cooling Processes

In the first simulation, denoted as M82_10_1, I took $\dot{M}_{\rm hot}$ = 10 $M_{\odot}$ yr*-1* and $\epsilon=1.0$ (such that $\dot{E}=10^{50}$ ergs / yr), I neglected all radiative cooling processes, and assumed that the gas contained no dust. Although the outflow is a steady-state structure in this model, we can nevertheless plot the properties of the material as a function of time since passing though the driving radius. The top panels of Figure 2 show the evolution of the radial distance from the galaxy, radial velocity, sound speed, and average temperature in this case. Because $\dot{Q}=0,$ an exact solution to eqs. (1) - (3) can be obtained as described in Chevalier & Clegg (1985). In the region outside of the driving radius, this is

[TABLE]

where the Mach number, $M\equiv v/c_{s},$ and the second equality holds when $\gamma=5/3.$ The lower left panel of Figure 2 shows the fractional difference between this expression and the Mach number in the simulations. Here we see that the numerical and analytical results agree closely, differing by less than $0.2\%$ as the gas moves from as distance of 300 pc to 30 kpc from the host.

The central panel of the lower row of Figure 2 shows the fraction of the gas in the simulation with temperatures between $10^{4}$ K and $3\times 10^{4}$ K, a rough measure of the gas that will be visible though measurements of low-ionization state ions. As cooling within the simulation is uniform and caused by adiabatic expansion, this leads to 100% of the simulation moving into the cold state at a time set by the initial conditions and eq. (19).

Finally, the right panel of Figure 2 shows the evolution of the average density in the simulation, the rms fluctuations about the mean density, $(\left<n^{2}\right>-\left<n\right>^{2})^{1/2}$, and the rms cold gas density, $\left<n_{\rm cold}^{2}\right>^{1/2}$. While the overall average density is determined directly by $v$ and $r$ as $n\propto(vr^{2})^{-1},$ the rms density fluctuations give a measure of the inhomogeneity of the medium as a function of time, which must be solved for numerically. The evolution of the rms fluctuations shows some early rearrangements, followed by a graduate decline in $(\left<n^{2}\right>\left<n\right>^{-2}-1)^{1/2}$ as pressure differences smooth inhomogeneities as described in more detail below. The decrease in $(\left<n^{2}\right>\left<n\right>^{-2}-1)^{1/2}$ continues for roughly $\approx 10^{6}$ years, until the simulation enters the phase of high Mach number expansion, after which inhomogeneities stay largely constant.

This freeze-out of structure at late times can be understood from the fact that pressure differences can only smooth out structures whose sound speed is greater than their rate of expansion. The overall expansion factor is $a(t)=[n(t)/n_{\star}]^{-1/3}=[vr^{2}/(v_{\star}r_{\star}^{2})]^{1/3},$ such that at large radii, $a(t)\approx(2r/r_{\star})^{2/3}.$ This means, at large radii, the condition that the expansion rate of a clump of physical size $x_{\rm clump}$ be less than the sound speed of the medium gives

[TABLE]

which in comoving coordinates, $x_{\rm clump,com}\equiv x_{\rm clump}\,a^{-1}$ gives

[TABLE]

where in the final expression I have also assumed that cooling is minimal such that eq. (19) holds approximately and $a\approx 0.6M.$ Thus the drop in the sound speed locks in density fluctuation into the medium when $r/r_{\star}$ becomes large, an issue we will return to in §3.3.

Finally, as the density fluctuations never become large in this simulation, the rms cold gas density $\left<n_{\rm cold}^{2}\right>^{1/2}$ is very close to the average density during the period in which the cold fraction is one. In this adiabatic simulation, this occurs when the mean density has dropped significantly, meaning that emission from low-ionization state ions is likely to be negligible in this model, as discussed in more detail in §3.4.

Figure 3 examines the impact of additional cooling processes on the evolution of an M82 type outflow with $\dot{M}_{\rm hot}$ = 10 $M_{\odot}$ yr*-1* and $\dot{E}=10^{50}$ ergs / yr. The figure shows results for three different cases: case M82_10_1A including atomic/ionic cooling assuming solar metallicity, case M82_10_1AC, including atomic/ionic and Compton cooling, and case M82_10_1ACD, including atomic/ionic cooling, Compton cooling, and dust cooling and destruction assuming an initial Milky Way dust-to-gas ratio. For comparison, in these plots I also compute the rate of cooling due to adiabatic expansion:

[TABLE]

In the M82_10_1A, case atomic/ionic cooling overtakes the adiabatic cooling rate at $\approx 10^{6}$ years, when the temperature has dropped to $\approx 10^{6}$K and emission processes are very efficient. At this point, the sound speed and temperature drop rapidly, causing a decrease in pressure that leads to a radial acceleration that is larger than in the adiabatic case. On the other hand, this same pressure decrease means that at very large radii the radial velocity is always less than in the adiabatic case, leading to a final velocity of 970 km s*-1* as opposed to 1000 km s${}^{-1}.$

A more dramatic effect of atomic/ionic cooling is its impact on the evolution of density fluctuations. Because cooling occurs first in the densest regions, the pressure within them drops significantly with respect to their surroundings, leading to a rapid increase in $(\left<n^{2}\right>\left<n\right>^{-2}-1)^{1/2}$. Thus radiative cooling not only causes the gas to reach $10^{4}\,-\ 3\times 10^{4}$K temperatures at higher densities, but also amplifies densities inhomogeneities, increasing the observability of the medium.

Moving to the M82_10_1AC case that includes Compton cooling, we uncover even larger changes with respect to the adiabatic case. In fact, for the first 300,000 years of evolution, cooling by this process is over an order of magnitude more efficient than atomic/ionic cooling, and almost as efficient as adiabatic cooling itself. This causes the temperature to be significantly lower at very early times and leads to a short phase in which $v$ exceeds the adiabatic value, but it ultimately causes the final velocity to be 920 km s*-1* as compared to 1000 km s*-1*. Because Compton cooling depends on the scattering of electrons with photons, the energy loss per unit mass is independent of density, and unlike atomic/ionic cooling, it does not lead to the growth of inhomogeneities. On the other hand, it allows the gas to cool into the $10^{4}\,-\ 3\times 10^{4}$K range at smaller radii and at higher densities.

Finally, the lower two rows of Figure 3 show the results from the M82_10_1ACD simulation that includes dust cooling and destruction. While dust cooling is extremely efficient at the high densities that occur when $r\approx r_{\star}$ in this model, these high densities also lead to rapid dust sputtering from collisions between dust particles and the $10^{7}$K medium. Thus, at least for this choice of parameters, dust makes a difference only at the very earliest times, having an overall impact that is similar to slightly reducing $\epsilon,$ the fraction of the available supernova energy added into the hot medium.

3.2. Effect of Varying Outflow Properties

Next we turn our attention to the relationship between the initial wind properties and the subsequent evolution of the outflowing gas. In Figure 4, as in the previous section, we consider three M82-like cases in which $r_{\star}$ = 300 parsec, the SFR = 10 $M_{\odot}$ yr${}^{-1},$ and $\dot{E}=10^{50}$ ergs / yr, but now we take $\dot{M}_{\rm hot}$ = 5 $M_{\odot}$ yr*-1* (case M82_5_1AC with $v_{\star}=710$ km s*-1*), $\dot{M}_{\rm hot}$ = 2 $M_{\odot}$ yr*-1* (case M82_2_1AC with $v_{\star}=1120$ km s*-1*) and $\dot{M}_{\rm hot}$ = 1 $M_{\odot}$ yr*-1* (case M82_1_1AC with $v_{\star}=1580$ km s*-1*). Note that the last of these cases corresponds to pure supernova ejecta with no energy losses and no additional entrained material. All models include atomic/ionic cooling, and cooling by Compton scattering and dust emission. Decreasing the mass outflow rate while leaving the energy input rate unchanged has several important effects. First it leads to higher temperatures and pressures, which in turn boost the radial velocity of the flow. This means that the density of the medium decreases due both to lower mass fluxes, and also due to the material is moving outward more rapidly. As atomic/ionic cooling and dust cooling per unit mass are proportional to $n$, this means that lower $\dot{M}_{\rm hot}$ cases not only start with more internal energy per particle, but this energy is more difficult to radiate away.

On the other hand, as the Compton cooling rate per particle is proportional to the gas temperature and independent of the gas density this process increases in importance in lower $\dot{M}_{\rm hot}$ cases. Thus in the M82_5_1ACD case, Compton cooling dominates over atomic/ionic cooling for roughly the first $10^{6}$ years of evolution, as opposed to $\approx 3\times 10^{5}$ years in the M82_10_1ACD case. While dust cooling is rapid, it is limited to very early times, much as was found in the M82_10_1ACD case. All together, these processes are less important than adiabatic expansion. Thus cooling has only a moderate impact in this simulation, decreasing the final velocity of the flow by $\approx 100$ km/s, and slightly increasing the inhomogeneity of the medium.

Radiative cooling is even less important in the $\dot{M}_{\rm hot}$ = 2 $M_{\odot}$ yr*-1* case shown in the central panels of Figure 4. In this case atomic/ionic cooling is always at least an order of magnitude smaller than adiabatic cooling, causing it to have little impact on the evolution of the medium or the evolution of the density structure. Furthermore, while Compton cooling is significantly higher in the initially hot medium, it still falls far short of the adiabatic cooling rate. While dust cooling is also much smaller than adiabatic cooling, dust sputtering proceeds slowly, due to the lower densities, which leads to appreciable dust densities up to $\approx 10^{5}$ years after the material has moved out of the driving region. This is late enough that the material has undergone significant acceleration, and is well into the supersonic regime. In this case, no stage of rapid cooling is seen whatsoever, and the difference between the final velocity (2120 km/s) and the final velocity in the adiabatic case (2240 km/s) is less than 5%.

In the M82_1_1ACD case with $\dot{M}_{\rm hot}$ = 1 $M_{\odot}$ yr${}^{-1},$ Compton, atomic/ionic, and dust cooling are to too week to lead to observable effects, and the reduction in the final velocity is again only a few percent. Interestingly however, gas densities in this model have dropped to the point that sputtering is no longer able to destroy the dust before it leaves the galaxy. This means that this case is able to deposit dust into the surrounding circumgalactic medium directly through the hot medium, without requiring the any cold clouds to be entrained by the hot gas. This is the hottest and most rapidly-moving of all the simulated cases, but it is also the only one in which dust is not completely destroyed within the hot wind. This suggests that, perhaps counterintuatively, starbursting galaxies whose winds entrain the least amount of material may nevertheless be the most efficient at depositing the dust observed at distances of 20 kpc to several Mpc outside of galaxies (Ménard et al., 2010; Ménard & Fukugita, 2012).

Table 1 also include describes models including energy losses within the driving region, in which the fraction of supernovae energy deposited in the outflowing medium was taken to be 50%, $\epsilon=0.5.$ In the M82_5_0.5ACD case with $\dot{M}_{\rm hot}$ = 5 $M_{\odot}$ yr${}^{-1},$ this reduces the initial velocity to $v_{\star}=500$ km s${}^{-1},$ from 710 km s*-1* and in the M82_2_0.5ACD case with $\dot{M}_{\rm hot}$ = 2 $M_{\odot}$ yr${}^{-1},$ it reduces the initial velocity to $v_{\star}=790$ km s${}^{-1},$ from 1120 km s${}^{-1}.$ It also reduces the initial temperature of the winds by a factor of $2$ with respect to their $\epsilon=1$ values. Such slower, colder winds, are naturally much more susceptible to cooling. In the $\dot{M}_{\rm hot}$ = 5 $M_{\odot}$ yr*-1* case, this leads to an $\approx 9\%$ reduction in the final velocity and the rapid production of heterogeneous, cold gas within $\approx 1$ kpc of the galaxy. In the 2 $M_{\odot}$ case, cooling processes are more modest, such that the final velocity is reduced by $\approx 6\%$, and no phase of rapid structure growth occurs. In both cases, the densities are high enough that dust is rapidly destroyed by sputtering.

Figure 5 shows the results of simulations of a more extreme case of a ULIRG type galaxy with a SFR of $\dot{M}$ = 100 $M_{\odot}$ yr*-1* and $\epsilon=1$ such that the energy input is $10^{53}$ ergs yr*-1*. Such galaxies have similar sizes as M82 (e.g. Soifer et al., 2000), and so for ease of comparison I take $r_{\star}=300$ pc. I also consider mass outflow rates of 50 $M_{\odot}$ yr*-1* and 20 $M_{\odot}$ yr*-1* such that $M_{\odot}$ yr*-1*/SFR is 0.5 and 0.2 as in the M82 cases above. Thus the initial conditions in the ULIRG_50_1ACD case are the same as in the M82_5_1ACD case except the number density of the gas and the photon flux from the galaxy are both shifted upward by a factor of 10. This in turn increases the importance of atomic/ionic, Compton, and dust cooling, boosting them well above adiabatic cooling. The result is widespread cooling and growth of inhomogeneities within $10^{5}$ years, at a distance just outside the driving radius. Thus, unlike the M82_5_1ACD case, the ULIRG_50_1ACD case is so radiatively efficient that it is $10^{4}$K throughout the full evolution outside the driving region. In the 20 $M_{\odot}$ yr*-1* case, cooling processes are slightly reduced, but again the gas cooling occurs very near to the galaxy. Here $10^{4}$K temperatures are reached within $3\times 10^{5}$ years and $r\approx 600$ pc.

Table 1 also describes ULIRG models with outflow rates of 50 $M_{\odot}$ yr*-1* and 20 $M_{\odot}$ yr*-1* in which the efficiency of energy input is 50% ($\epsilon=0.5$). The higher densities and lower temperatures in these models boost the cooling efficiencies even further such that the final velocities are much lower than their ($2v_{\star}$) adiabatic values and $10^{4}$K temperatures are reached within $10^{5}$ years, corresponding to radii that are barely beyond the driving radius. Furthermore, although the star-formation in ULIRGs is dust-enshrouded, the high densities and temperatures in the outflow mean that dust is destroyed almost immediately as it leaves the galaxies. This suggests that the most likely state of outflowing material from the dustiest, most rapidly star-forming galaxies in the Universe is cold and dust-free.

3.3. Growth of Structure

Figure 6 shows the density distribution on slices taken from the M82 type simulations with various cooling processes included, at representative times of $\approx 3\times 10^{4},$ $\approx 3\times 10^{5}$, and $\approx 3\times 10^{6}$ years. Each panel is labeled with its corresponding values of the outflow radius and the mean temperature. Between $3\times 10^{4}$ and $3\times 10^{5}$ years, cooling is dominated by adiabatic expansion and by Compton cooling (in the simulations in which it is included). Because both these processes lead to uniform cooling per unit mass, they do not cause pressure differences, and thus the overall changes in the structures between these two times is minor.

When the medium cools below $\approx 10^{6}$K however, the efficiency of atomic/ionic cooling increases dramatically. Unlike expansion and Compton cooling, this process is density-dependent, leading to the most rapid cooling per unit mass in the overdense regions. As these regions cool, the higher pressure gas around them compresses them further, leading to a sharp increase in the magnitudes of density fluctuations, as evidenced in the increase in the rms density fluctuations seen in Figure 3. Thus by $3\times 10^{6}$ years, the density fluctuations have increased from their values at $3\times 10^{5}$ years by roughly a factor of 10. This rapid growth of fluctuations is soon halted, however, as the exterior medium has a chance to cool and pressure contrasts become much more modest. At this time the average temperature in the medium is $\lesssim 10^{4}$K, which corresponds to a sound speed of $\approx 10$ km/s. As the medium is moving $\approx 1000$ km/s, its expansion is rapid with respect to these motions. From eq. (21) this means that all fluctuations with comoving sizes above $\approx 0.6r_{\star}M^{-1}(r/r_{\star})^{-1/3}\approx 2\,{\rm pc}(r/r_{\star})^{-1/3}$ are subsequently frozen into the medium.

Figure 7 shows the density distribution on slices taken from M82 type simulations with mass outflow rates smaller than the fiducial value of 10 $M_{\odot}.$ Note that in this figure, the density fluctuations are shown on the same scale in all panels. As in the 10 $M_{\odot}$ case, above there are only minor changes in the density distribution in the period between $3\times 10^{4}$ and $3\times 10^{5}$ years, when adiabatic and Compton cooling are dominant. At later times, however, the evolution of the medium is strongly dependent on the mass outflow rate. The M82_5_1ACD case with $\dot{M}_{\rm hot}$ = 5 $M_{\odot}$ yr*-1* shows a significant increase in the magnitude of density fluctuations at $3\times 10^{6}$ years, which is caused by density-dependent atomic/ionic cooling. However, because the mean density in the simulation is less than in the $\dot{M}_{\rm hot}$ = 10 $M_{\odot}$ yr*-1* case, the overall impact of this cooling is also less, increasing density fluctuations by a factor of $\approx 3$, rather than a factor of 10 as seen above.

In the 2 $M_{\odot}$ yr*-1* case, the impact of cooling from radiative processes in decreased even further, such that by the time the simulation moves into the $\lesssim 10^{6}$ temperature regime, the outflow has expanded to the point that density fluctuations are already largely frozen into the medium. Thus there is no particular increase in the rms fluctuations at $3\times 10^{6}$ years seen in this model in Figure 7, which is also consistent with the overall evolution of $(\left<n^{2}\right>-\left<n\right>^{2})^{1/2},$ as shown in Figure 4. Similarly, in the 1 $M_{\odot}$ yr*-1* case, cooling from atomic/ionic processes is even less important, and there is no phase in which density contrasts are increased.

To quantify these effects further, I computed the density power spectra from the simulations as

[TABLE]

where $n({\bf k})\equiv\int d^{3}x\,n({\bf x})\,\exp(-i{\bf k\cdot x})$ and the wave number $k$ is equal to $2\pi$ divided by the wavelength of the perturbation. In this case, the total density fluctuations in units of the mean density are

[TABLE]

Figure 8 shows the evolution of the structure in the simulations normalized as $4\pi(k/2\pi)^{3}P(k)\left<n\right>^{-2},$ meaning that the maximum values of these curves will be roughly equal to the variance of the density fluctuations in units of the mean density.

In the adiabatic M82_10_1 case, which appears in the upper left panel of this figure, there are no processes in place to enhance density inhomogeneities, leading to a steady decline in the power spectrum over time. This is most dramatic in the interval from $t=0$ to $t\approx 3\times 10^{4}$ in which pressure differences in the initial conditions act to rapidly smooth structure. While this smoothing occurs on all scales, it is most significant at the shortest distances, corresponding to $2\pi/k\lesssim 2$ comoving parsecs on this plot. Note also that the simplified initial conditions as given by eq. (18) cause some aliasing in the $3\times 10^{4}$ power spectrum, leading to a saw-tooth pattern that would not be expected in a less-idealized configuration. The pressure smoothing continues at later times, strongly suppressing structures below $\approx 3$ comoving parsecs by $3\times 10^{5}$ years, but having a much more subtle impact at later times, as the medium expands to the point that the remaining homogeneities become frozen into the medium.

In the M82_10_1A case in which atomic and ionic cooling are included, the early evolution is almost identical to that seen in the M82_10_1 case. Between $3\times 10^{5}$ years and $3\times 10^{6}$ years, however, as the outflow cools adiabatically to below $\approx 10^{6}$K, density-dependent radiative cooling becomes dominant. This leads to the rapid growth of structure, boosting the power spectrum at all scales, and increasing the overall rms value of the fluctuations by a factor $\approx 10$, as also seen in Figure 6.

The inclusion of Compton and dust cooling, as shown in cases M82_10_1AC and M82_10_1ACD has only a minor impact of the overall evolution of the power spectra, but the impact of changing the mass outflow rate is much more significant. Thus the curves from the M82_5_1AC and M82_2_1ACD cases illustrate quantitatively how decreasing $\dot{M}_{\rm hot}$ from 10 $M_{\odot}$ yr*-1* to 5 $M_{\odot}$ yr*-1* leads to a muted period of structure growth through density-dependent cooling, while decreasing it to 2 $M_{\odot}$ yr*-1* leads to an overall evolution that is extremely similar to the adiabatic M82_10_1A case.

Similarly, when the fraction of the supernova energy that is deposited into the hot medium is reduced to $\epsilon=0.5,$ not only is the radial velocity of the outflow reduced, but the density of the medium goes up, increasing the efficiency of radiative cooling. In the M82_5_0.5ACD case with 5 $M_{\odot}$ yr*-1* and $\epsilon=0.5$, this causes the phase of structure growth from density-dependent cooling to become much more significant, leading to an overall history similar to the M82_10_1ACD case. On the other hand, when the mass outflow rate is 2 $M_{\odot}$ yr${}^{-1},$ reducing $\epsilon$ from 1 to 0.5 does not increase the importance of radiative cooling to the point that it leads to significant structure growth. Instead, the large difference in structure growth between these two similar models means that the observability of the material cooling from outflows is likely to be a strong function of the host galaxy properties, hinting at the wide range of Lyman-$\alpha$ luminosities seen in rapidly star formation galaxies (Martin et al., 2015).

The last four panels of Figure 8 show the evolution of the density power spectra in the ULRIG cases. Unlike the M82 runs, structure grows most rapidly in the ULIRG simulations between $3\times 10^{4}$ and $3\times 10^{5}$ years, because the high-mass fluxes in these models lead to rapid density-dependent cooling much earlier in the outflow evolution. On the other hand, in the ULIRG simulations, this early cooling also causes the underdense regions to quickly drop to the $\approx 2\times 10^{4}$K point at which cooling and photo-heating are balanced (see Figure 1). Since this now happens are small radii, these fluctuations are not yet frozen into the medium. Thus pressure difference begin to act to reverse the growth of structure, leading to density distributions at $3\times 10^{6}$ years that are somewhat more homogeneous than in the M82 cases.

4. Observability

In Figure 9, I show the radial distribution of $\left<n_{\rm cold}^{2}\right>^{1/2},$ the rms number density of material between $10^{4}$K and $3\times 10^{4}$K, a measure of the emission expected in low ionization state metal lines. Note that although each simulation tracks a single comoving region 150 parsec on a side as it moves outward, this is a piece of an overall steady-state outflow that will continue for as long as the mass and energy input are sustained. Thus at any given time, the entire radial distribution depicted in these figures would be expected to contribute to observed emission. Moreover, as the models only include photoheating from the source galaxy and not the additional contribution from the metagalactic background or the flux from nearby galaxies, at large distances the density of gas with observable low-ionization state ions may be somewhat greater than shown here. This means the value in Figure 9 should be thought of as a conservative lower limit on the rms number density of observable cold material, especially at large distances.

The upper left panel of this figure shows the results from the fiducial M82 simulations with mass outflow rates of $10M_{\odot}$ yr${}^{-1},$ and various cooling processes included. Here we see that it is Compton cooling that has the largest impact on the radius at which gas cools efficiently. This is due both to the suppression of atomic/ionic cooling and the enhancement of Compton cooling at short distances from the starburst. On the other hand, the presence of dust has only a small impact on the cold gas properties, because it is quickly destroyed at the high densities found in this case.

Moving to M82 models with the same energy flux ($\epsilon=1$) but lower mass fluxes, the impact of atomic/ionic cooling and Compton cooling is dramatically reduced. This is due to the fact that the outflows are launched at higher temperatures and at higher velocities, such that these processes have more internal energy per particle to remove and less time to do so. Thus these cases only cool below $3\times 10^{4}$K at distances beyond $\approx 2$ kpc and at rms densities below $\approx 10^{-2}$cm${}^{-3},$ meaning that their emission is likely to be very difficult to observe. When the energy flux is reduced by a factor of two ($\epsilon=0.5$), cooling processes are substantial in the $\dot{M}_{\rm hot}=5M_{\odot}$ yr*-1* case, but minimal in the $\dot{M}_{\rm hot}=2M_{\odot}$ yr*-1* and $\dot{M}_{\rm hot}=1M_{\odot}$ yr*-1* cases, illustrating the strong sensitivity of the wind to initial conditions.

In the ULIRG case in which the star formation rate, mass outflow rate, and energy input rate are all boosted by a factor of ten, rapid cooling is a generic feature of all models. Thus in the $\dot{M}_{\rm hot}=50M_{\odot}$ yr*-1* cases $\left<n_{\rm cold}^{2}\right>^{1/2}$ approaches $\approx 10$ cm*-3* at distances just outside the driving radius, while in the $\dot{M}_{\rm hot}=20M_{\odot}$ yr*-1* cases, $\left<n_{\rm cold}^{2}\right>^{1/2}$ approaches $\approx 1$ cm*-3* at distances of $\approx 500-800$ parsecs. Note also that in the $\dot{M}_{\rm hot}=50M_{\odot}$ yr${}^{-1},$ $\epsilon=0.5$ case the outflow shows a double peaked profile, due to the clumping of structure that occurs as $T$ drops rapidly through the unstable regime.

As a more direct measure of observability I also computed the approximate luminosity per line-of-sight sight velocity for a simple model in which observations are taken straight along the axis of the outflow, which encompasses $\pi$ steradians, $\pi/2$ on each side of the galaxy. In this case, the luminosity in line $i$ per line-of-sight velocity is given by

[TABLE]

where $\Lambda_{i},$ which has units of erg s*-1* cm3 is the average power per unit volume emitted in a line of interest by a medium with an average cold gas density squared of $\left<n_{\rm cold}^{2}\right>.$ Here a positive $v_{z}$ is defined as a velocity towards the observer, corresponding to a blueshift of the emission line, and typical values of $\Lambda_{i},$ are on the order of 10*-23* to 10*-22* ergs s1 cm3 (Boehringer & Hensler, 1989; Sutherland & Dopita, 1993; Wiersma et al., 2009; Gnedin & Hollon, 2012).

In Figure 10, I plot $\Lambda_{-22,i}^{-1}\frac{dL_{i}}{dv_{z}}$ from each of the simulations, the luminosity normalized by $\Lambda_{-22,i}\equiv\Lambda_{i}/(10^{-22}\,{\rm ergs}\,{\rm s}^{-1}{\rm cm}^{3}).$ As in Figure 9, the normalized luminosities shown here represent lower limits because none of these models includes a metagalatic background. The upper left panel of this figure shows the results from the fiducial M82 simulations with mass outflow rates of $10M_{\odot}$ yr*-1* and various cooling processes included. Here we see that while the inclusion of additional cooling processes has a large impact on the radial distance at which gas cools, it has a much smaller impact on the luminosity as a function of velocity. This is due to the fact that although the additional processes cause the gas to cool at higher densities, they also cause the volume of the outflow in the $10^{4}$ to $3\times 10^{4}$K range to decrease significantly. Note also that the range of line-of-sight velocities over which outflowing gas is observed has very little to do with its radial distance, and the profiles are instead dominated by geometric effects.

In M82 models with $\epsilon=1$ and lower mass fluxes, the relative inefficiency of atomic/ionic cooling and Compton cooling cause the cold gas to only be found at low densities and radial velocities approaching the maximum possible value of $2v_{\star}.$ These effects cause the overall normalized luminosity to decrease and the line-of-sight velocity to increase rapidly as a function of $M_{\odot}$. Thus the $5M_{\odot}$ yr*-1* case is roughly 3 times fainter than the $10M_{\odot}$ yr*-1* case, with a velocity offset that is 1.5 time as large, while the $2M_{\odot}$ yr*-1* and $1M_{\odot}$ yr*-1* cases are 100 and 1000 times fainter than the fiducial case, respectively, with velocity offsets that are 2 and 3 times as large. Furthermore, while decreasing the energy input ($\epsilon=0.5$) boosts the luminosities and decreases the distances somewhat, the very strong trend with mass flux means that the low $M_{\odot}$ outflows are likely to remain undetectable. Note that the strong dependence of the total luminosity of the outflow to the properties of the host is consistent with the large dynamic range of Ly$\alpha$ luminosities observed in outflowing galaxies (Martin et al., 2015).

It also means that moving to higher mass fluxes, as in the ULIRG cases, yields overall luminosities that are much higher than in the M82 cases. Furthermore, the strong cooling in these models, driven both by the higher densities as well as the higher background fluxes, means that the cold material is found at velocities well below $2v_{\star},$ the maximum velocity without cooling included. Thus the $\approx 1000$ km/s maximum velocities found in the ($\epsilon=0.5$) cases are consistent with maximum Ly$\alpha$ and [OIII] $\lambda$5007 blueshifts observed in ULIRGs (Martin et al., 2015) although the velocities are slightly higher than typically found, while the luminosities are slightly lower, suggesting that models with $\epsilon<0.5$ may provide even better matches to the observations.

On the other hand, the velocities in the models are significantly greater than observed in absorption line measurements of Na I (e.g. Martin, 2005, 2006; Rupke et al., 2005) and OH (e.g. Veilleux et al., 2013). One possibility is that the molecular and NaI absorption arise from cold clouds that are accelerated and shredded at small distances by the hot wind. This entrainment and heating of this material would then lead to additional radiative losses within the driving regime, such that smaller $\epsilon$ values would be favored, leading to strong Ly$\alpha$ and [OIII] emission as the hot wind itself cooled as it expanded into the circumgalactic medium. Further modeling and comparisons with ULIRG observations are needed to properly assess this scenario.

5. Conclusions

Starburst-driven galaxy outflows play a key role in the history of structure formation and are well observed at a variety of redshifts. However, our ability to interpret these observations is strongly hampered by the fact that the $\approx 10^{4}$K material that is easiest to observe is also the most difficult to model theoretically. While it is often assumed that this material arises from cold clouds accelerated by the hot wind, detailed models of cold-cloud hot-wind interactions show that instabilities and evaporation are likely to shred clouds well before they are significantly accelerated. A second possibility is that the cold material is formed directly from the cooling of hot wind material, cooling it into the observable range after it is already moving at high radial velocities.

To test this scenario, I have carried out a suite of fully-three dimensional simulations of initially-hot material ejected by starburst-driven galaxy outflows. The simulations are conducted in a comoving frame that moves with the material, and, building on the cooling routines presented in Ferrara & Scannapieco (2016), they simulations include atomic/ionic cooling that accounts for the ionizing background, as well as Compton cooling, and cooling from dust grains. I find that atomic/ionic cooling is significantly reduced due to photoionization, particularly at small radii around galaxies with higher star formation rates. Instead, it is Compton scattering with photons from the host galaxy that has the largest impact, matching or exceeding temperature losses due to all other mechanisms including adiabatic expansion. Compton cooling becomes particularly important when the SFR is large, such that the galaxy is luminous, and when the mass flux is large, such that the outflowing material spends a long time near the galaxy.

The most important role of atomic/ionic cooling, on the other hand, is enhancing inhomogeneities in the medium. While it is subdominant at high temperatures, these processes overtake the adiabatic and Compton cooling rates once the overall temperature has dropped to below $\approx 10^{6}$K. Because collisional cooling occurs most rapidly in the densest regions, the pressure within them drops significantly with respect to their surroundings, leading to further increases in density. This rapidly amplifies densities inhomogeneities by roughly an order of magnitude until the medium as a whole reaches $\approx 10^{4}$K and the sharp decrease in sound speed freezes out further growth of structure.

On the other hand, while dust is an extremely efficient coolant, it has little effect on outflow evolution, as it is rapidly destroyed in the extreme conditions found in most galaxy outflows. In fact the only case that dust is preserved throughout the full outflow evolution is the fastest moving, lowest density case I considered. This suggests that small starbursts with modest mass outflow rates, rather than large dusty starbursts, may be the most efficient at polluting the circumgalactic medium with dust.

The simulation results can be used to construct a simple steady-state model of the observed UV/optical emission from each outflow, adopting a simple approximation that the outflow is observed along its central axis. In this case, the radial distance at which the outflow cools has only a minor impact on the observations, as models that cool more quickly have higher cold gas densities but also much smaller cold-gas volumes. In fact, the spread in the of line-of-sight velocities of the cold material has very little to do with its radial distance, and the profiles are instead dominated by geometric effects.

On the other hand, the total luminosities and maximum velocities are strong functions of the properties of the host, which is consistent with the large dynamic range of Ly$\alpha$ and [OIII] luminosities observed in outflowing ULIRGs. In fact, the luminosities and maximum velocities found in the models with significant energy losses ($\epsilon=0.5$) are consistent with such observations, although the velocities are greater than observed in NaI and OH absorption line studies. It may be that absorption lines probe entrained cold material, while emission lines probe clumps of cold material condensing out the hot wind, a possibility that deserves further theoretical study.

Acknowledgments

I would like to thank Andrea Ferrara, Crystal Martin, and Sylvain Veilleux for insightful discussions that greatly improved this paper. This work was supported by the National Science Foundation under grant AST14-07835 and PHY11-25915 and NASA theory grant NNX15AK82G. I thank the Simons Foundation and the organizers of the workshop Galactic Winds: Beyond Phenomenology (J. Kollmeier and A. Benson) and the Kavli Institute for Theoretical Physics and the organizers of the program The Cold Universe (P. Caselli, A. Ferrara, M. Ouchi, R. Schneider, and J. Tan), for discussions that helped to inspire and refine the work carried out here. I would also like to thank the Texas Advanced Computing Center (TACC) at The University of Texas at Austin, and the Extreme Science and Engineering Discovery Environment (XSEDE) for providing HPC resources via grant TG-AST130021 that have contributed to the results reported within this paper.