THE ORIGIN AND KINEMATICS OF COLD GAS IN GALACTIC WINDS: INSIGHT FROM NUMERICAL SIMULATIONS

ULIRGs are starburst galaxies with infrared luminosity >10¹² L_☉, and are usually found in major mergers and interacting galaxies (Sanders et al. 1988). They are believed to go through starburst phases twice, when the gas in a galaxy with a prograde orbital geometry is tidally disturbed during the first encounter with another galaxy and when both galaxies meet again and finally merge (Mihos & Hernquist 1996; Murphy et al. 2001; Li et al. 2004). We choose to model a molecular gas-rich spiral galaxy on its first encounter with another galaxy of similar mass.

We set up a molecular disk with M_g = 10¹⁰ M_☉ in a dark matter halo with M_halo = 5 × 10¹² M_☉. CO observations of ULIRGs at both first and second passages show the presence of molecular disks with M_g = (0.4–1.5) × 10¹⁰ M_☉ (Sanders et al. 1988; Solomon et al. 1997), which is in the range found for gas-rich spiral galaxies. However, the emission originates in regions a few hundred parsecs in radius, yielding surface densities of ∼(0.5–1) × 10⁴ M_☉ pc⁻², within which the molecular mass dominates the dynamical mass (Sanders et al. 1988; Solomon et al. 1997). At these high surface densities, molecular hydrogen will dominate (Blitz & Rosolowsky 2006), as it can form within a few million years in turbulent regions with densities over 100 cm⁻³ (Glover & Mac Low 2007). The density of H₂ traced by CO emission is ∼500 cm⁻³, comparable to the envelope of giant molecular clouds, while a region of much higher density in ULIRGs is traced by HCN emission, ∼10⁵ cm⁻³, comparable to star-forming cloud cores (Solomon et al. 1992).

We assume that the entire interstellar medium (ISM) is a scaled-up version of a normal galactic disk with the ambient densities a factor of ∼100 higher, making even the intercloud medium a molecular region. Thus we assume that the surface density distribution of the molecular disk is exponential, with Σ(R) = Σ₀ exp(−R/R_d) where R_d is a scale radius (see also Σ(R) of Arp 220 by Scoville et al. 1997).

We choose to model a disk with a central surface density of Σ₀ = 10⁴ M_☉ pc⁻², with a disk scale radius R_d = 0.7 kpc. The disk is in hydrostatic equilibrium with a Navarro et al. (1997; hereafter NFW) halo potential, and a disk potential based on the thin disk approximation (Toomre 1963), since M_dyn ≈ M_g. The NFW potential is

$$\begin{equation} \Phi (x) = \frac{G M_{\rm halo}}{R_v} \frac{\ln (1+cx)/x}{F(c)}, \end{equation} \tag{ 1 }$$

where we set the virial radius R_v = 326 kpc, x = r/R_v, c is a halo concentration factor, set to c = 5, appropriate for a large halo (Jimenez et al. 2003), and F(c) = ln(1 + c) − c/(1 + c).

The velocity dispersion of the molecular gas is observed to be 90 km s⁻¹ in Arp 220, which appears to be at the end of the merging process (Scoville et al. 1997). Such a high velocity dispersion yields a scale height of 15 pc in its disk (Scoville et al. 1997). Molecular clouds with such high dispersion will get destroyed by colliding with other clouds, so the cooling time behind the shocks must be shorter than the destruction time interval. We assume the gas is supported by turbulence with a similarly high velocity dispersion c_s = 55 km s⁻¹. It is lower than 90 km s⁻¹ because gravity from our disk gas and halo cannot confine the gas with a higher velocity dispersion. As a comparison, we also model a disk with a higher surface density of Σ₀ = 5 × 10⁴ M_☉ pc⁻² with R_d = 0.17 kpc and c_s = 90 km s⁻¹, based on the central surface density observed in Arp 220. This is the highest Σ₀ of all observed ULIRGs.

Figure 1 shows the vertical density distributions of both disks. The exponential scale heights are rather small, 7 and 2 pc, but the gas within them is very dense, ∼5000 and 10⁵ cm⁻³, respectively. The gas density is still ∼500 cm⁻³ at Z ≳ 4R_d. At higher altitudes, where the gas is less dense, the gas is physically atomic or even ionized. We do not take that state change into account in our model, though. When the number density drops to n = 10⁻² cm⁻³ we set the gas density in the halo constant as it is no longer dynamically important on the length and timescales treated in our model.

Zoom In Zoom Out Reset image size

We assume a single starburst that occurs at the center of the disk, and that all the kinetic energy of the starburst supernovae is released in a central wind of constant mechanical luminosity. In reality, the discrete energy inputs from supernovae generate blastwaves that become subsonic in the hot interior of the bubble first produced by stellar winds, and hence can be treated as a continuous mechanical luminosity in the study of bubble dynamics (Mac Low & McCray 1988). These assumptions mean that a single superbubble forms, evolving to produce a bipolar outflow of gas.

Figure 2 shows the evolution of mechanical luminosity L_mech as a function of time for an instantaneous starburst with 10⁹ M_☉ of gas turning into stars and for continuous starbursts with star formation rates of 100 and 500 M_☉ yr⁻¹, based on the Starburst 99 model (Leitherer et al. 1999). The Starburst 99 model uses a power-law initial mass function with exponent α = 2.35 between low-mass and high-mass cutoff masses of M_low = 1 M_☉ and M_up = 100 M_☉ with solar metallicity. Star formation rates in ULIRGs are estimated to be ≳100 M_☉ yr⁻¹ based on far-infrared luminosities and the assumption of continuous star formation (see Table 1 of Martin 2005; note that the star formation rates given there correspond to a low mass cutoff of 0.1 M_☉ and must be divided by a factor of 2.55 before comparison to the population synthesis models).

Zoom In Zoom Out Reset image size

The amount of mechanical power supplied per unit stellar mass depends on the star formation history. About ∼40 million supernovae, for example, will be produced by an instantaneous starburst of M_* = 10⁹ M_☉. For a continuous starburst, L_mech increases until the death rate of massive stars catches up to their birth rate, after about 40 Myr. The power rises particularly rapidly over the first few Myr, the period modeled by our simulation. Figure 2(a) shows the evolution for a continuous starbursts with 500 M_☉ yr⁻¹ and 100 M_☉ yr⁻¹. For our ULIRG models, we use constant mechanical luminosity winds with values L_mech = 10⁴³, 10⁴², and 10⁴¹ erg s⁻¹. The highest of these corresponds to the mechanical luminosity expected from stellar winds during the first 2 Myr of an instantaneous starburst with M_* = 10⁹ M_☉. The subsequent supernovae will result in a far higher mechanical luminosity, but that is likely to be vented out of the galactic disk through the hole opened by the initial blowout. L_mech = 10⁴³ erg s⁻¹ also corresponds to a model with a constant star formation rate (SFR) of 15.9 M_☉ yr⁻¹. This correspondence assumes that the birth and death rates of massive stars are in equilibrium, which is achieved about 40 Myr after the burst begins. In this scenario, the burst would been ongoing through its initial stages before our simulation starts, but gas displaced by the initial feedback had been replaced by inflows, and prior feedback energy was largely radiated. Our simulation is not such a good representation of this scenario.

We note it is an oversimplification to assume a single starburst at the disk center for modeling a galactic outflow. Bipolar outflows are seen in some starburst galaxies such as M82 (Strickland & Stevens 2000; Strickland et al. 2004) for example, although some ULIRG winds appear to require starburst regions extended to ≳1 kpc to launch the cool outflow (Martin 2006). Three-dimensional, hydrodynamic simulations of dwarf starbursts have shown that extended, multiple energy sources, as well as a single central energy source, form a bipolar outflow (Fragile et al. 2004). These simulations also demonstrated that the main effect of multiple sources was to reduce the fraction of metals and energy ejected from the galaxy from almost unity to around 50%. Our assumption thus represents a reasonably strong lower bound to the amount of kinetic energy that will be deposited in the observed cold gas. We therefore start with this assumption and do not expect the results to differ much from those expected with more extended star-forming regions.

In addition, we neglect the effects of UV radiation on molecular hydrogen in the disk. The UV radiation from massive stars may photodissociate some of molecular hydrogen outside star-forming cores to atomic hydrogen. However, the assumed turbulent pressure with c_s = 55 km s⁻¹ is 14 times greater than the increased thermal pressure by photodissociation. We thus safely neglect the effects of UV radiation on the disk gas structure.

We define the model with Σ₀ = 10⁴ M_☉ pc⁻² and L_mech = 10⁴³ erg s⁻¹ as our fiducial model (U1/X1). We list the parameters for all the other runs in Table 1.

Table 1. Parameters for Starburst Models

Modela	log10Σ0b (M☉ pc−2)	log10Lmechc (erg s−1)	$\dot{M_{\rm in}}$d (M☉ yr−1)	RSNe (pc)	Tfloorf (K)	Resolution (pc)
X1	4	43	17	50	104	0.2
X1-0	4	43	17	25g	104	0.1
X1-2	4	43	17	50	104	0.4h
X1-4	4	43	17	50	104	0.8h
U1	4	43	17	50	102	0.2
U1-A	4	43	1.7	50	102	0.2
U1-B	4	43	49	50	102	0.2
U1-C	4	43	120	50	102	0.2
X2	4	42	17	50	102	0.2
X3	4	41	17	50	102	0.2
S1	4	43	17	25	102	0.2
V1	4.7	43	17	50	102	0.2

Notes. ^aAll models run with ZEUS. ^bCentral surface density. ^cMechanical luminosity. ^dMass loading rate. ^eSize of source region where supernova energy is injected. ^fThe minimum temperature floor for cooling. ^gThe level of noise on the surface of the source region is kept same with that of our standard model, X1 by setting the number of cells covering the source region the same. ^hRatioed grids are used (the resolution is 0.2 pc within the source regions).

Download table as:  ASCIITypeset image

In the molecular disk of a ULIRG with a very small scale height, a bubble with L_mech ⩾ 10⁴¹ erg s⁻¹ quickly blows out of the disk at t ≪ 1 Myr. We can only simulate the evolution of the bubble up to t ∼ 0.3–1 Myr before it leaves our grid, which is rather small because of the cost of high resolution. We argue in Section 5 that cooled, swept-up shells, which we identify as Na i-absorbing gas, can acquire maximum velocities primarily determined by when they blow out, accelerate, and fragment. The question remains whether these shell fragments and clumps survive as the hot interior wind streams through them.

Each individual shell fragment after blowout remains unstable to smaller scale R–T instabilities while being further accelerated by the wind, requiring extremely high resolution to fully resolve. Resolution is not as critical for previously published numerical studies on starburst winds that explored feedback parameters. For example, reasonable assumptions about cooling losses indicate that moderate luminosity starbursts do not remove a significant fraction of the galaxy's gas from the halos (Mac Low & Ferrara 1999), and the mass loss rates would only decrease further with more cooling. A large fraction of the heavy elements do escape from the halos (Mac Low & Ferrara 1999; Fujita et al. 2004), and this result does not depend on resolution below ∼ a few tens of parsecs, so long as at least a few R–T modes are excited in the swept-up shells to let the metal-enriched gas escape. Therefore, a major question that must be answered to understand observations of cold gas at high velocity is whether clumps of dense gas survive in the lower density wind (see Klein et al. 1994; Marcolini et al. 2005), where they will fragment because of hydrodynamic instabilities, such as R–T (strictly speaking, Richtmyer–Meshkov; see Richtmyer 1960; Meshkov 1969) and Kelvin–Helmholtz.

Weaver et al. (1977) argued that the density in the hot interior of bubbles in uniform gas is dominated by conductive evaporation from the dense shell. However, our model does not explicitly include thermal conduction, or material ablated off of high-density molecular clouds associated with the central starburst. Instead, we add additional mass to our central luminosity source to account for this process, multiplying the mass input rate by a mass-loading factor ξ. The amount of thermally evaporated mass in a bubble expanding into a uniform medium is proportional to L^27/35_mechρ^−2/35. The internal temperature is

$$\begin{equation} T_b(t)=(\gamma -1)\frac{\mu }{k_{B}} \left[\frac{5}{11}\int _{0}^{t}L_{\rm mech}(t^{\prime }) dt^{\prime } \right] \bigg/ \left[\int _0^t\xi M_{SN}(t^{\prime }) dt^{\prime } \right], \end{equation} \tag{ 2 }$$

where L_mech(t) and the mass of supernova ejecta M_SN(t) as a function of time are taken from Starburst 99 model, the adiabatic index γ = 5/3, mean mass per particle μ = 14/22m_H, and k_B is the Boltzmann factor. Weaver et al. (1977) showed that 5/11 of the total input mechanical energy goes into the hot, pressurized region bound by the inner and outer shock fronts (see their Equation (14)). This was applied to superbubbles by McCray & Kafatos (1987). We show T_b(t) in Figure 3(a).

Zoom In Zoom Out Reset image size

We can estimate the terminal velocity of the wind driven by such a bubble by equating the bubble interior energy with the kinetic energy of the mass-loaded wind, so that

$$\begin{equation} v^2_{\rm wind}(t)=2\left[\frac{5}{11}\int _0^t L_{\rm mech}(t^{\prime }) dt^{\prime }\right]\bigg/\left[\int _0^t \xi M_{SN}(t^{\prime })dt^{\prime }\right], \end{equation} \tag{ 3 }$$

which we plot in Figure 3(b). Since the mechanical luminosity L_mech(t) and the mass of supernova ejecta M_SN(t) are both linearly proportional to the amount of gas converted to stars, T_b and v_wind are the same for all strengths of starburst with the same mass-loading rate ξ. We take a fiducial ξ value of about 8 corresponding to a mass-loading rate $\dot{M_{\rm in}}=17 (L_{\rm mech}/10^{43})\;{M}_{\odot }\;{\rm yr}^{-1}$, but run test simulations with ξ varying between ∼1 and 15 to explore the effects of mass loading on the shell kinematics. Comparisons to observations suggest ξ ≈ 10 (Suchkov et al. 1996—see M82 comparison; Martin et al. 2002). Throughout this paper, we designate as the wind the hot interior gas freely streaming outward after blowout of the shell.

4.1.1. Mechanical Luminosities

Our fiducial model (X1) has mechanical luminosity L_mech = 10⁴³ erg s⁻¹. This model corresponds to the first 2 Myr of a starburst in which 10⁹ M_☉ of gas turns into stars instantaneously. We compare this to models with lower mechanical luminosities L_mech = 10⁴² and 10⁴¹ erg s⁻¹ (models X2 and X3, respectively) in the same molecular disk. These mechanical luminosities correspond to instantaneous bursts of 10⁸ and 10⁷ M_☉. Figure 4 shows the density distribution of our fiducial model in its right panel. This may be compared to Figure 5, which shows the density distributions of the two models with lower L_mech at t ≈ 0.49 Myr and 0.85 Myr, respectively. The swept-up shells fragment due to R–T instability into multiple clumps and shells.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Secondary Kelvin–Helmholtz instabilities ablate the sides of these fragments as the hot gas streams through them. Look at the clumps, for instance, at (R,Z)≈(40, 65), (30,110), and (15–25,120) pc in X1. Note also that the swept-up shells in the horizontal direction are also R–T unstable, because our disk gas is stratified in the radial direction, too, due to its exponential surface density profile.

In fact, the degree of fragmentation is larger in X2 and more so in X3 because $\dot{M_{\rm in}}\propto L_{\rm mech}$ and so the density of interior gas is lower. In particular, most shell fragments in X3 are already falling back to the disk. A mechanical power of 10⁴² erg s⁻¹ is just too small in such a dense environment to accelerate the bulk of the shells to the disk's escape velocity.

4.1.2. Surface Density

4.1.3. Mass Loading

Mass loading from thermal conduction and molecular clouds determines the density of the bubble interior and wind. Figure 6 shows the density distributions of our fiducial ULIRG model with T_floor = 10² K with mass-loading rates of 1.7, 17, 49, and 120 M_☉ yr⁻¹ (models U1-A, U1, U1-B, and U1-C). These mass-loading rates correspond to bubble interior temperatures, T_b = 1.7 × 10⁸, 2.2 × 10⁷, 7.5 × 10⁶, and 3.2 × 10⁶ K and wind terminal velocities expected when all the thermal energy is converted to kinetic energy, v_wind = 2700, 1000, 600, and 250 km s⁻¹, respectively.

Zoom In Zoom Out Reset image size

Figure 6 shows the bubbles just before they leave the grid, at t = 0.23, 0.27, 0.35, and 0.41 Myr, respectively. Since the input mechanical luminosity is the same in all these models the bubbles initially grow at about the same rate, driven by the thermal pressure of the hot interior gas. However, a bubble with a lower mass-loading rate and higher wind terminal velocity expands faster into the halo once the swept-up shells fragment and the hot gas blows out between the fragments. In addition, the blowout occurs earlier and in more places with a lower mass-loading rate because the density of the bubble interior gas is lower. A higher density contrast between the hot interior and the swept-up shells promotes shell fragmentation by R–T instability, the growth of which is proportional to the density contrast. In the least dense model U1-A, all the swept-up shells quickly fragment into fingers and filaments and the dense clumps at their edges are subject to secondary Kelvin–Helmholtz instabilities.

As the hot bubble interior becomes transonic after blowout, it accelerates the shells and shell fragments by ram pressure rather than thermal pressure. This low-density wind thus can accelerate the shells to higher velocity after blowout. We quantitatively study the effects of wind ram pressure in the next section (Section 5.4). The observed X-ray temperature T_X is ∼0.67 keV =7.7 × 10⁶ K in all kinds of starburst galaxies from dwarfs to ULIRGs (Martin 1999; Heckman et al. 2001; Huo et al. 2004; Grimes et al. 2005). This corresponds to ξ ≈ 10 in Equation (3). However, the X-ray emission is proportional to n² so it is biased toward high-density regions such as the interface between the hot interior gas and the shells and their fragments. At this interface, conductive evaporation and turbulent ablation raise the density. The recent observations of diffuse hard X-ray emission in starburst galaxies suggest the existence of a very hot (log T ≳ 7.5) metal-bearing gas (e.g., Strickland et al. 2004; Strickland & Heckman 2007). Then the bulk of the hot wind may still be very hot ∼10⁸ K, the temperature which Strickland & Stevens (2000) modeled for M82 with ξ = 1.

4.2.1. Grid Resolution

Figure 7 shows the density distributions of our fiducial model with grid resolution varying from 0.1 to 0.8 pc. We chose T_floor = 10⁴ K for this resolution study, because shell densities are not too far from what we expect analytically with this high minimum temperature. The growth of R–T instability is significantly enhanced in the highest resolution run, X1-0, and suppressed in the lower resolution runs, X1-2 and X1-4. In particular, all the outermost shells seen in X1 are further fragmented by R–T instability in X1-0 with the resolution increased by a factor of only 2. The positions of outer shock fronts in the horizontal direction agree very well among the four simulations, since the shells there are not subject to severe hydrodynamic instabilities.

Zoom In Zoom Out Reset image size

Figure 8 shows the density profiles at the outer shock fronts before any fragmentation occurs for the runs in our standard resolution study. The bubble in our highest resolution run X1-0 grows a little slower in the beginning, because we chose the size of the source region in X1-0 to be half of that in X1 in order not to overproduce noise at the contact discontinuity. We will show below that this noise feeds R–T instability and must be maintained the same in order to study the effects of resolution on shell fragmentation alone. Thus the density profile of X1-0 at t = 0.06 Myr is not directly comparable to those of other runs at t = 0.05 Myr, but Figure 8 still demonstrates the trend in resolving the peak shell density as a function of resolution.

Zoom In Zoom Out Reset image size

Figure 8 shows that the shell density is progressively better resolved as the resolution increases. However, the cooling times at the shock front are typically of order 10² yr (see Section 4.2.2), so a shock with cooling floor equal to the background temperature of 10⁴ K may be treated as isothermal. At the time displayed, the Mach number of the outer shock is ${\cal M}=19$. The shell density we expect from an isothermal shock propagating into a background density ρ_bg is $\rho _{bg} {\cal M}^2 = 6.1\times 10^{-19}\;{\rm g\;cm}^{-3}$. Figure 8 shows that the shell density ρ_shell is still unresolved by a factor of ∼4 even with our highest grid resolution of 0.1 pc.

The difference we see in the development of R–T instability develops because of two factors. First, the postshock density of the swept-up shells is not fully resolved. Increasing density contrast drives faster R–T growth. Second, the linear R–T instability grows more quickly at smaller wavelengths, but the nonlinear development moves to increasing larger wavelengths as competition between growing modes becomes important (Youngs 1984). At least 10–25 zones is required to resolve the smallest modes, though (e.g., Mac Low & Zahnle 1994), so grid resolution matters critically for the initial development and transition to nonlinearity. However, we believe that we reached the point where resolution effects are no more important than physics we have not included such as magnetic fields, nonequilibrium cooling, thermal conduction, and photoionization, as well as the assumption of azimuthal symmetry (see Section 3).

For example, strong magnetic fields B ∼ 20 μG observed in the Antennae merging galaxy (Chyży & Beck 2004) can potentially inhibit the formation of cold, dense shells or suppress their fragmentation. Our study does show that the degree of R–T fragmentation is important to reproduce the observed wide range of cooled shell fragments. Recall, however, that the shell density in our simulations can be more than an order of magnitude underresolved. Thus the degree of fragmentation will not be significantly overestimated in our simulations unless the magnetic fields are strong enough to reduce the shell density by more than that.

The density ρ_sh expected behind an isothermal shock running into a magnetic field with Alfvenic Mach number $\cal {M_A}$ and having $1 \ll {\cal M}_A \ll {\cal M}$ is $\rho _{sh} = \rho _0 \sqrt{2}{\cal M}_A$ (Draine & McKee 1993). Before the effects of magnetization become important in limiting R–T instability, the ratio between the magnetized and unmagnetized postshock densities must be of order $2^{1/2} {\cal M}_A / {\cal M}^2 < 10^{-2}$. With ρ ∼ 10⁻²⁴ g cm⁻¹ and v_s = 500–800 km s⁻¹, the shell density expected with B ∼ 20 μG becomes three order of magnitude lower than that without B. The effect of magnetic field thus becomes substantial only when a bubble grows to the high-Z, low-density part of the disk, above Z>100 pc. However, most of the fragmentation occurs within Z ∼ 100 pc especially with the highest-resolution run, Thus our results on the wide absorption profiles are robust. If anything, the suppression of fragmentation by magnetic pressure at high latitude will allow the hot wind to keep accelerating the outermost shells. Magnetic pressure might thus even increase the amount of cool gas with high terminal velocity.

We also note that resolving the shell density profile is not important to following the overall dynamical evolution of bubbles driven by the thermal energy of the hot, pressurized regions (Castor et al. 1975; Weaver et al. 1977), but is important to follow the details of shell fragmentation due to R–T instability and the fate of dense fragments and clumps by shocks and following hydrodynamic instabilities. We will show below that the resolution of 0.2 pc is still not sufficient to properly model hydrodynamic instabilities acting on shells and clumps, but is just sufficient for the purpose of demonstrating a wide range of velocities in shell fragments caused by R–T instability.

4.2.2. Cooling Temperature Floor

We show the density distributions of our fiducial model with different cooling curve cutoffs of T_floor = 10² K (U1) and 10⁴ K (X1) at the time of blowout in Figure 4. Figure 4 shows that the temperature floor we choose for our cooling function appears to have a negligible influence on the evolution of bubbles in our models.

However, closer examination reveals that the density of shocked shells is approximately the same in both simulations despite the difference in temperature floor. For example, at t = 0.05 Myr before any fragmentation occurs, the shell density in both simulations is ∼8 × 10⁻²⁰ g cm⁻³ in the vertical direction where the background disk density is ρ_bg = 2 × 10⁻²¹ g cm⁻³. This is because the density peak in the simulations is limited by resolution in these models, not by the strength of cooling. The shock velocity in the vertical direction at t = 0.05 Myr is ∼230 km s⁻¹, or Mach number ${\cal M} = 15$ in 10⁴ K background gas. The immediate postshock temperature is then T = 7.6 × 10⁵ K. This shocked gas quickly cools to or below 10⁴ K because the exponential cooling time (e.g., Mac Low & McCray 1988) is very short,

$$\begin{equation} \tau _{\rm cool} = 3kT/4n_{bg}\Lambda \approx 64\;{\rm yr}, \end{equation} \tag{ 4 }$$

with the mean mass per particle μ = 14/22m_H for ionized gas, and Λ(T) = 4.1 × 10⁻²³ erg cm³ s⁻¹ from the MacDonald & Bailey (1981) cooling curve.

The shell density expected from an isothermal shock will then be $\rho _{\rm shell}=\rho _{bg} {\cal M}^2\approx 5\times 10^{-19}\;{\rm g\;cm}^{-3}$ with Mach number ${\cal M}=15$ if T_floor = 10⁴ K. If T_floor = 10² K, the shell density will reach values even greater than the isothermal value. However, since the shell density is far from being resolved even with 0.1 pc resolution in our simulations, the cooling floor has a negligible influence on our models.

4.2.3. Effects of Source Region on Shell Fragmentation

We use a limited grid size in order to maintain high resolution, so we simulate the evolution of bubbles only up to the time of blowout. We showed in the previous section that the bubbles blow out very early, at t ≪ 1 Myr, because of the small scale height of our molecular disk. For comparison, ultraluminous starbursts have ages of up to 50 Myr (Murphy et al. 2001). To extrapolate our computational results to later times, we use the ballistic approximation (Zahnle & Mac Low 1995; Fujita et al. 2004) that after blowout, shell fragments travel on radial ballistic orbits in the gravitational potential of the galaxy, with no further accelerations by gas pressure gradients.

The strength of this approximation depends on the behavior of the wind, and conditions in the region into which the wind penetrates. This approximation was successfully used by Zahnle & Mac Low (1995) to follow the ejecta of a typical Shoemaker-Levy 9 fragment falling back to Jupiter's atmosphere, and was found to give results quantitatively consistent with the observations of the size and shape of the infrared bright spots. That case does differ from the starburst case in having neither an ongoing wind nor a complex density distribution above the site of the explosion. In the starburst case, continuing supernova explosions at late times drive an ongoing wind (see Figure 2), while the mergers that drive most starbursts yield a complex geometry above the site of the blowout. However, even if the stronger winds expected from supernovas at late times in an instantaneous starburst do accelerate the fragments further, the ballistic approximation will still yield a lower limit to their velocities. Since the question of most interest is whether cold gas can be accelerated to such high velocities, a lower limit is already useful. The complex geometry of tidal tails may also not be of great concern, as models suggest that they only cover a small fraction of the solid angle visible from the nucleus, and so are unlikely to be dynamically dominant.

Under the ballistic approximation, the equation of motion for a shell fragment at a distance r from the galactic nucleus is

$$\begin{equation} v(r) = \lbrace v^2(r_b)+2[\Phi (r_b)-\Phi (r)]\rbrace ^{1/2}, \end{equation} \tag{ 5 }$$

where Φ(r) is the total halo and disk potential, and v(r_b) = v_b is the shell velocity at blowout at a position of r_b.

In Figure 9, we plot v(r) for several initial velocities starting at a fixed initial position, r_b = 200 pc, just above our model disk. By setting r_b = 200 pc, we can ignore the disk potential which is negligible above the disk compared to the halo potential. Thus we set Φ = Φ_halo in Equation (5) and solve the equation analytically using Equation (1). The upper limit to the radius reached by a shell fragment at time t is given by the linear approximation r ⩽ v_bt, because the shell will only decelerate in the potential. We show actual radii for different initial velocities at t = 10 Myr by vertical ticks in the left panel of Figure 9. Although Figure 9 is plotted as a function of radial velocity, we find similar results if total velocity is used instead, because shell fragments are traveling in nearly radial directions. This figure shows that the linear approximation is quite good for v_b>100 km s⁻¹, so we use it to follow the cold gas for times of order 10 Myr.

Zoom In Zoom Out Reset image size

For times t ≫ 100 Myr, the linear approximation fails and most gas falls back to the center after reaching heights of a few hundred kpc. Only gas with v_b ≳ 1000 km s⁻¹ can completely escape the potential of the disk and halo (the halo alone has an escape velocity v_esc = 800 km s⁻¹ at radii less than 0.01R_v). We show below that very little gas actually escapes the halo, with ≲0.1% of the total shell mass accelerated above v_b ≳ 1000 km s⁻¹ in our fiducial model. Our results confirm the result from previous simulations of smaller galaxies that the loss of ISM mass is inefficient (e.g., Mac Low & Ferrara 1999; D'Ercole & Brighenti 1999). However, significant mass is circulated over scales of 10 kpc, presenting a significant absorption cross section as suggested by Martin (2006).

Only the coolest gas can be observed in Na i absorption. We chose a temperature cutoff of T < 5 × 10⁴ K to trace this gas because cooling remains numerically limited even with our highest resolution grid. Gas below this cutoff is so dense that the cooling time (Equation (4)) is very short, τ_cool < 0.01 Myr, so it is physically expected to reach temperatures where Na i is present. As an example, the least dense shell fragment in Figure 10 has ρ = 3.0 × 10⁻²⁴ g cm⁻³ and T = 3.3 × 10⁴ K. Equation (4) gives an exponential cooling time τ_cool ≈ 0.007 Myr, using Λ(T) ≈ 4.2 × 10⁻²³ erg cm³ s⁻¹ from the MacDonald & Bailey (1981) cooling curve.

Zoom In Zoom Out Reset image size

Figures 10 and 11 show the gas density, temperature, and radial velocity in models X1 and X1-0 along a line of sight through the central continuum source at an angle θ = 19° from the vertical axis. Note that our temperature cutoff picks up only the densest and the coolest parts of the shells excluding numerically diffused interfaces between the shells and the hot gas.

Zoom In Zoom Out Reset image size

We compare the characteristics of this temperature-selected gas to the measured properties of Na i absorption in ultraluminous starbursts. As argued above in Section 2.2, our assumption of a central point source of wind luminosity likely represents a lower limit to the amount of cold gas entrained. Fragile et al. (2004) found that, in galaxies where a central energy source ejects nearly all its kinetic energy in the hot gas (Mac Low & Ferrara 1999), distributed supernovas deposit as much as half of the energy in the cold gas. However, much of this energy is radiated promptly rather than converted to kinetic energy, so determining the significance of distributed supernovas requires further work.

5.2.1. Cool Gas Column Densities

Figure 7 shows that the lines of sight plotted in Figures 10 and 11 go through multiple distinct shells, adding to a total column density of cool gas N_H ∼ 7 × 10²¹ cm⁻². The inferred column densities in ultraluminous starbursts are a bit lower, exceeding 10²¹ cm⁻² in only one out of four ultraluminous starbursts (Martin 2005). Since the model column densities fluctuate by as much as an order of magnitude between neighboring sightlines, we plot the column density distribution as a function of angle in Figure 12. We find that most sightlines with N_H>5 × 10²² cm⁻² lie within 30° of the galactic disk where the shells are usually not subject to hydrodynamic instabilities. The typical value is N_H ∼ 10²² cm⁻², similar to the largest column density estimated from observations of ultraluminous starbursts (Martin 2005, 2006; Rupke et al. 2002).

Zoom In Zoom Out Reset image size

5.2.2. Geometric Dilution

5.2.3. Resolution Effects

A second reason we overestimate the cool gas column density is that the fragmentation of the shells is limited by numerical resolution, as well as our assumption of azimuthal symmetry. However, as discussed above in Section 4.2.1, the overall kinematics described in Section 5.1 should remain similar unless shell fragments are entirely mixed with hot gas and destroyed. To distinguish real physical effects from numerical artifacts, we study the kinematics of the cold gas along a θ = 19° line of sight in our fiducial simulation X1 and the high-resolution simulation X1-0.

The first two fragments from the center are parts of a filament remaining from the initial fragmentation of the swept-up supershell due to R–T instability at t ∼ 0.1 Myr. These high-density peaks are traveling with similar velocities v ∼ 400–500 km s⁻¹ in both simulations and dominate the column densities of cool gas in the studied sightlines. The outermost shell in X1 keeps sweeping up high-altitude disk gas, while in X1-0 this shell further fragments by R–T instability to the last two fragments in Figure 11. The coherence of the shell in X1 clearly occurs because of the lower resolution in that run. However, the shell in X1 and the outermost fragment of the two in X1-0 travel with similarly high velocity of ∼800–900 km s⁻¹, though they contribute very little to the total column density.

The hot wind overruns and shocks the first fragment in X1 and the first two fragments in X1-0. As a result, secondary R–T and K–H instabilities act on them, removing gas from the cold cloud and mixing it into the hot wind. In both cases, these clumps of gas are resolved by about ∼10 cells, demonstrating that this is a minimum size below which fragments tend to survive because the secondary instabilities cannot be resolved. Mac Low & Zahnle (1994) show that fragments must be resolved by at least 25 zones to resolve the secondary instabilities. Fragments smaller than that remain as artificially massive clouds in our simulations, overpredicting the cool gas column densities.

Klein et al. (1994) suggest that these clumps of gas should be destroyed, via hydrodynamic instabilities, over ∼10 shock crossing timescales

$$\begin{eqnarray} t_{cc} &=& r_c(v_w-v_c)\sqrt{\frac{\rho _{c}}{\rho _{w}}} =(0.04\;{\rm Myr}) \left(\frac{r_c}{2\;{\rm pc}}\right)\nonumber\\ &&\times\left(\frac{v_w-v_c}{500\;{\rm km\;s}^{-1}}\right) \left(\frac{\rho _c/\rho _w}{100}\right)^{1/2}, \end{eqnarray} \tag{ 6 }$$

where r_c is the radius of the clump, ρ_w and ρ_c are the density, and v_w and v_c are the velocity of the wind and the clump, and the scaling parameters hold for typical 10 zone clumps in our model. In reality, t_cc may be substantially longer because the density ratios ρ_c/ρ_w are probably underestimated: shell densities, and thus clump densities, are underresolved, while the wind density may be overestimated by our mass-loading scheme. We need ρ_c/ρ_w ≈ 10⁴ for these clumps to survive for more than ∼5 Myr, which may just be reachable in a rapidly diverging wind.

Moreover, recent numerical studies clumps may be stabilized against Kelvin–Helmholtz and R–T instabilities, reducing mass loss, by either thermal conduction (Marcolini et al. 2005; Vieser & Hensler 2007a), or even weak magnetic fields (Mac Low et al. 1994; Shin et al. 2008). Since we are not able to address the fate of shell fragments further in this study, we assume the bulk of the cool mass remains cool for the duration of starburst to be observed as Na i-absorbing gas. Larger, denser clumps are more likely to survive, in reality.

The bottom left panels of Figures 10 and 11 show the distribution of column density of cool gas N_H as a function of its radial velocity v_r at the end of the simulations. So long as the linear version of the ballistic approximation is correct, the velocity spread remains constant. Cool gas is seen over a wide range of velocities >450 km s⁻¹ at both resolutions. Shells that fragment earlier have been accelerated less, and so their fragments travel more slowly.

We now explore whether absorption from multiple shell fragments can explain the large absorption line widths measured for ultraluminous starbursts by considering the effects of resolution and viewing angle. Figure 13 shows the distribution of velocity widths Δv seen in cool gas with T < 5 × 10⁴ K in runs X1-0, X1, X1-2, and X1-4, with grid size increasing from 0.1 to 0.8 pc, along sightlines spaced by 1° from the vertical axis. We define Δv as the difference between the maximum and the minimum velocities seen in the cool gas along a given sightline—that is, full width at zero. This resolution study demonstrates that the fraction of sightlines with large velocity spread increases as the resolution improves up to 0.2 pc, but then appears to converge.

Zoom In Zoom Out Reset image size

Comparison to Figure 7 shows that the largest velocity widths are seen in sightlines intersecting the largest number of fragmented shells. Better resolving postshock densities in the swept-up, cooled, shells leads to quicker shell fragmentation and reformation, producing a larger number of shells, and thus the wider range of velocities in the fragments.

Figure 13 also shows that most sightlines in X1 and X1-0 are above the average observed line width 〈v〉 = 320 ± 120 km s⁻¹ (dashed lines). Observers measure the average of the sightlines toward continuum sources subtending about an arcsecond, and this corresponds to a length of 1.84 h_0.7 kpc at redshift 0.1. Assuming that the observers are likely to view ULIRGs at t = 5–10 Myr after the onset of starburst, we can suggest all the shell fragments within ∼20°–40° at a given sightline in our simulations will contribute to the absorption profiles. Hence, our models suggest that observers will measure a large line width regardless of viewing angle.

To quantify this, we compute the average velocity width of all sightlines in an axisymmetric cone within θ as

$$\begin{equation} \overline{\Delta v}(\le \theta) = \frac{\int _0^{\theta } \Delta v(\alpha) \sin (\alpha) d\alpha }{\int _0^{\theta } \sin (\alpha) d\alpha }. \end{equation} \tag{ 7 }$$

The value of $\overline{\Delta v}$ is lowered by sightlines near the disk midplane where the bubbles are not blowing out and sightlines intersecting holes made by blowout with nearly zero Δv. The average velocity widths are clearly smaller for lower resolution runs, however. We find $\overline{\Delta v}(\theta \le 30^{\circ })$ is 220 km s⁻¹ for run X1-0, and 210 km s⁻¹ for X1, but only 110 km s⁻¹ for X1-4.

Figure 14 shows the distribution as a function of angle with degree spacing of mass-weighted average velocities of cool gas v_av (solid lines) for our standard resolution study (runs X1-0 through X1-4). The average shell velocities plotted may be misleading, since shell mass differs significantly at each sightline. Therefore, we also plot a mass-weighted average shell velocity across a 10° arc centered on each sightline v_av,10. This ought to be the quantity most directly comparable to observations of the average velocity.

Zoom In Zoom Out Reset image size

On sightlines with multiple fragments, the average mass-weighted velocity reflects the velocity of the more massive fragments. These are fragments of the initial swept up shell, which follow ballistic orbits with little acceleration after shell fragmentation. Thus, average mass-weighted velocity tends to lie substantially below peak velocity.

The observable quantity v_av,10 shows a clear converging trend with resolution. The converged value toward the pole appears to be under 400 km s⁻¹, with correspondingly lower values at other angles, as shown in Figure 14. Figure 14 suggests that cool gas in shells and their fragments will be observed to be traveling with v_av,10 ≈ 200–350 km s⁻¹ at angles to the pole θ < 60°, the angle where blowout occurs in all models.

Our model is consistent with observations that constrain ULIRG wind geometry. Figure 14 shows that absorption at velocities v>100 km s⁻¹ is detected at all angles greater than 10° from the disk plane. It follows that 98 % of all random, radial sightlines would exhibit such outflows. This result is consistent with 15 of 18 observed objects showing such outflows in Martin (2005), and a similar fraction seen by Rupke et al. (2005).

Figure 15 shows the distributions of v_av and v_av,10 for models X2 and X3 run with decreasing mechanical luminosity at the same resolution as model X1. The mass-weighted average velocities of shells should depend on the mechanical luminosity L_mech driving the bubble by thermal pressure. In a spherical bubble, the dependence would be L^1/5_mech. This dependence is actually slightly stronger in our models of blowout. We find v_av,10 dropping from 400 km s⁻¹ to ∼100 km s⁻¹ moving from model X1 to X3, with each step having an order of magnitude lower mechanical luminosity. This gives an empirical relation closer to L^0.3_mech. Note the high spike in v_av,10 of X3 around ∼15° is biased by a small amount of cool mass present in its vicinity. By way of comparison with observations, Figure 6 of Martin (2005) shows velocities reaching v = 700 km s⁻¹ for an SFR of 1000 M_☉ yr⁻¹. At an SFR of 1 M_☉ yr⁻¹, the sparse data shown in that Figure suggest v = 30–40 km s⁻¹, while our empirical relation would suggest a value v ≃ 90 km s⁻¹. More data at low SFRs and a broader range of models will be required to establish whether there is truly a discrepancy with the model results.

Zoom In Zoom Out Reset image size

The mass-weighted average velocities v_av,10 in X1 and X2 agree with the observed shell velocity, but v_av,10 in X3 is well below the observed value. In disks as massive as these, starbursts with mechanical luminosity lower than 10⁴² erg s⁻¹ produce the lower end of the outflow velocity range observed in ULIRGs, which average v_s,obs = 330 ± 100 km s⁻¹ in the Martin (2005) sample and 170 km s⁻¹ in the column-density weighted average velocities from the Rupke et al. (2005) sample. Column-density weighting with velocity can be confounded by variation in covering factor with velocity, as has been shown for AGN outflows (e.g., Arav et al. 1999, 2005). We cannot model covering factors well with our two-dimensional simulations, so we postpone consideration of the question of column-weighted velocity to future work. We do note that overpredicting acceleration is much harder to do than underpredicting it, though, so we think our model is likely to be robust. The velocity spread, which is our most important result from the simulations, is similar to that found in both ULIRG studies.

5.5.1. Resolution

Figure 14 also shows the distribution of terminal velocity of cool gas v_term at each sightline. We define the maximum velocity at blowout as the terminal velocity a shell will ever acquire, following the ballistic approximation (Section 5.1). The maximum velocity of shell fragments is very high, v_term>500 km s⁻¹ at the angle θ < 60° where blowout occurs in all the runs. This high-velocity cool gas is found in the fragments of the outermost shells in our highest resolution model X1-0. The hot wind continues to accelerate a piece of shell, sweeping up the ambient gas, until it fragments further. As discussed above, shell fragmentation is very sensitive to resolution. Thus it is important to test the convergence of the mass of high-velocity gas.

Figure 16(a) shows the mass distribution of cool gas as a function of velocity in X1-0, X1, and X1-2 at t = 0.27 Myr. The fraction of cool gas traveling at high velocities is low. For example, the mass traveling with v ⩾ 500 km s⁻¹ is less than a few percent of the total shell mass, and the mass traveling faster than the observed terminal velocity of 750 km s⁻¹ is <0.1%.

Zoom In Zoom Out Reset image size

However, these results are best understood as upper limits. The amounts of cool gas in the high-velocity end are progressively smaller as the resolution increases. Roughly factor of 2 decreases in the mass of cool gas with velocities above 500 km s⁻¹ occur between runs with factor of 2 improvements in linear resolution. This suggests that we have not yet converged on the actual mass of high-velocity gas, although we have set good upper limits. We think this high-velocity gas is not likely to go away entirely, even if we run a simulation with a higher resolution and with additional physics. However, we cannot fully quantify its amount with our simulations.

The convergence properties of the peak velocity for cold gas v_term are somewhat better. The highest resolution models show <20% variations, suggesting that the general result is reasonably robust. The observed maximum value of 750 km s⁻¹ is consistent with the models at angles θ < 60°.

5.5.2. Mass Loading

Figure 17 shows the distributions of v_av, v_av,10, and v_term for simulations with different mass-loading rates and so different wind terminal velocities: U1-A, U1, U1-B, and U1-C. Cool gas has higher terminal velocities in bubbles with hotter winds that themselves have higher velocities v_wind. The hottest wind, in model U1-A, has v_wind = 2700 km s⁻¹, while the coldest wind, in U1-C, has v_wind = 250 km s⁻¹. Figure 16(b) shows the mass distribution of cool gas as a function of velocity for the four runs at the time of blowout. The cool gas with high v_term carries only a small amount of mass (see Section 5.5.1). For example, the fraction of cool gas with v>500 km s⁻¹ is ≲2% of the total shell mass, and the fraction with v>750 km s⁻¹ is ≲0.5%. The bulk of swept-up and cooled gas is driven to v ∼ 400 km s⁻¹ by thermal pressure of hot sonic gas, but a small fraction of it seems to be accelerated to higher velocity by the ram pressure of the same hot gas as it accelerates to supersonic velocities during blowout.

Zoom In Zoom Out Reset image size

5.5.3. Mechanical Luminosity

5.5.4. Summary

To directly compare to observed profiles, we generate Na i λ5890 absorption line profiles along sightlines through our simulations, as well as generalizing this procedure to fully model the observed doublet Na i λ, λ5890, 5896. To generate the profiles, we begin with the line intensity

$$\begin{equation} I_{\nu }=I_{\nu }(0) \exp ^{-\tau _{\nu }}=I_{\nu }(0) \prod _{i=1}^N \exp ^{-\tau _i} \end{equation} \tag{ 8 }$$

where the optical depth through cell i at frequency ν is τ_i, the sightline intersects N cells, and the background continuum is I_ν(0). We normalize the continuum by setting I_ν(0) = 1 and compute the profile as a function of the macroscopic velocity of Na i-absorbing gas, v. We compute I_v in the simulations by setting v = c(ν − ν₀)/ν₀, and computing the optical depth contributed by each cell in each of 1000 velocity bins. The optical depth in each cell

$$\begin{equation} \tau _{i}(v)=N_{{\rm Na}\,\hbox{\scriptsize {\sc i}}}\, s \, \lambda _{{\rm Na}\,\hbox{\scriptsize {\sc i}}} P(\Delta v), \end{equation} \tag{ 9 }$$

(Spitzer 1978) where $N_{{\rm Na}\,\hbox{\scriptsize {\sc i}},i}$ is the column density of Na i in a cell i, and the absorption cross section s_ν integrated over frequency ν for hν ≫ kT is s ≈ 2.654 × 10⁻²f₅₈₉₀ with the oscillator strength f₅₈₉₀ = 0.6. The Maxwellian velocity distribution function, P(Δv), is

$$\begin{equation} P(\Delta v)=\frac{1}{\sqrt{\pi }b}\exp {-(\Delta v^2/b^2)}, \end{equation} \tag{ 10 }$$

with Δv = v − v_i and $b=\sqrt{2k_B T_i/\mu m_H}$ for thermal broadening. Our simulations do not track chemical abundances, ionization state, or dust depletion, so we do not directly predict the Na i column. For the purpose of illustration, we compute the total Na i column from the total H i gas column using the conversion used by Martin (2005) to estimate $N_{{\rm H}\,\hbox{\scriptsize {\sc i}}}$ from their $N_{{\rm Na}\,\hbox{\scriptsize {\sc i}}}$ measurements, $N_{{\rm Na}\,\hbox{\scriptsize {\sc i}}} = 1.122\times 10^{-6}N_{{\rm H}\,\hbox{\scriptsize {\sc i}}}$.

We generated Na i 5890 line profiles along lines of sight through models X1 and X1-0 to show the effect of numerical resolution on the line profiles in Figure 18. Each sightline is described by its inclination from the polar axis of the simulation. The profile shapes reflect the different structure in these simulations. We can, for example, compare the absorption profiles in the middle left panel in Figure 18 with the density and velocity distributions of shells in Figures 10 and 11. The sightlines at θ = 19° intersect three and four shell fragments in X1 and X1-0, respectively. Each of these fragments generates a Δv ≳ 50–100 km s⁻¹ absorption line. A few of these lines are optically thick with the line profile becoming completely black at the line center. As the wind expands, column densities will drop linearly with radius due to geometric dilution. The complex of lines is spread out over ∼400 km s⁻¹.

Zoom In Zoom Out Reset image size

Figure 18 shows that the line profiles from the high and low resolution models present very similar structure for the most part. One minor exception is the minimum velocity. From 10° to 20°, the line profiles show a sharp cutoff at a velocity ∼500 km s⁻¹ in X1, but ∼3–40 km s⁻¹ in X1-0. This cutoff reflects the velocity of the first R–T fragments that form from the accelerating swept-up shells, which is lower in the high-resolution simulation. The overall lineshape does appear reasonably well converged.

Observers see the average over parallel sightlines subtending a few kiloparsecs of the disk. They must also contend with the instrumental response function, and the doublet nature of the Na i line, which has components at 5890 Å and 5896 Å, with optical depths differing by a factor of ∼2. To compare to the actual observed profiles, we generate the Na i 5890/5896 doublet in the frame of Na i 5890 by assuming the ratio of equivalent widths is only 1.3, typical of the observed lines in ULIRGs (rather than the optically thin limit of two).

We present average line profiles of the doublet over a 20° wide set of radial sightlines spaced at degree intervals. Although this is not exactly what is measured, it is comparable because these rays will subtend about 1–2 kpc of the shell when the bubble is ∼10 Myr old, so they intersect the same region of the shell, albeit at slightly different angles. We are actually doing the analysis at an earlier time, when the bubble is a factor of 100 smaller, and correspondingly higher column density, though, so we also must apply a geometric dilution factor to the column densities.

Figure 19 compares our simulated doublet line profiles, with and without geometric dilution of the column densities by a factor of 100, to five typical ULIRG spectra from Martin (2005) that have an instrumental broadening with FWHM ≃ 65 km ⁻¹. We convolve in a Gaussian of that width to simulate the broadening. The widths of the absorption profiles are very similar between our model at θ ⩽ 50° where blowout occurs and the observations. The undiluted absorption in our model much exceeds that observed. Our model overpredicts the cool gas column for two reasons (see Section 5.2). First, because we analyzed the models at a very early time in the starburst. At later times, the spherical expansion of the wind and the cold gas embedded within it dilutes their column densities. Second, the fragmentation of the shell is incompletely resolved, so some material remains cool that in reality would have been mixed into the hot wind. We account for the first factor by directly reducing the column densities, producing the diluted profiles shown in the figure, which reproduce the observed intensities well.

Zoom In Zoom Out Reset image size

The ULIRG spectra show absorption beginning from the systemic velocity, and extending to high velocities. Although this zero velocity absorption is absent in Figure 18, it is seen in the simulated doublet profiles shown in Figure 19, as it comes primarily from contributions from the λ5896 line of the doublet blueshifted into zero velocity of the λ5890 line (note our model does not include the contribution from stellar absorption).

5.6.1. Long Slit Observations

Measured Doppler shifts across ULIRGs show the cool outflow is extended spatially; and some outflows present a significant velocity gradient over kpc scales (Martin 2006). Our simulation ends long before the wind has reached such large scales and does not include the rotation of the galactic disk. The simulation does demonstrate, however, the amplitude of the velocity gradient that arises across the minor axis due to projection effects.

The wind is launched without any net angular momentum in our two-dimensional simulation. We examined whether, in the absence of rotation, the position affected the velocity much. Figure 20 shows terminal and mass-weighted average velocity along slits oriented at 30° from the major and minor axes of the disk, using parallel lines of sight through a three-dimensional rotation of the two-dimensional simulation. The absorption properties along the major axis of the galaxy will be symmetric about its center. An asymmetric velocity gradient is produced along the minor axis because one side of the outflow cone is directed more along the sightline. Substantial variations in velocity width and average velocity are seen at ∼40 pc and ±30 pc. These occur where sightlines intersect massive, slow-moving shells at θ>45° (see Figure 14) as well as light, fast-moving blowout fragments at θ < 45°. On the other hand, very high terminal velocities are seen on most sightlines since they intersect the fast blowout components.

Zoom In Zoom Out Reset image size

The simulated sightlines encounter both fast-moving and slow-moving gas at various latitudes, unlike the radial sightlines that we studied in previous sections. Our simulations viewed along nonradial sightlines naturally account for a wide velocity range of cool gas starting at v ≈ 0 and a high terminal velocity.

Finally, we crudely try to separate low-ionization, Na i-absorbing gas and ionized Hα-emitting gas in the cooled shell fragments in our simulations. Our aim is to compare the line widths of both components. Heckman et al. (2000) observed no correlation between the absorption and emission line widths for these lines, and drew the conclusion that only one of the lines could originate in the swept-up shells.

To model the Hα emission, we must decide whether the emitting gas is primarily photoionized by the central starburst or shock ionized in the wind. Shock ionization dominates at least some starburst-driven winds, such as NGC 1482 (Veilleux & Rupke 2002). Some observed ULIRGs also show extended shock-excited nebulae (Monreal-Ibero et al. 2006). However, the dynamics observed by Heckman et al. (2000) seem unlikely to come from shock ionized gas, since they see no correlation with the motions of the cold gas that presumably trace the shock. Therefore we assume photoionization by the central starburst, choosing an ionizing photon luminosity Q and a density distribution from a time the galaxy is likely to be observed. Figure 2 shows the ionizing photon luminosities as a function of time for three starburst cases that we consider. Since the propagation of ionization fronts critically depends on both the densities and the positions of shells, we cannot extrapolate the density distribution at blowout to t ≈ 10 Myr using the ballistic approximation. We can, however, still vary the strength of photon luminosity over a few orders of magnitude to examine the effect of attenuating the photon flux by 1/r² on our existing simulations. For example, a shell will experience about three orders of magnitude less photons per unit area after it travels with 500 km s⁻¹ for 10 Myr from r_b = 200 pc. For our purpose of demonstrating the lack of correlation between the line width of Na i-absorbing and Hα-emitting gas, this crude method is sufficient.

To solve for the transfer of ionizing radiation across our grids, we use the photon-conserving radiative transfer code developed by Abel et al. (1999) in the same manner as it was used in Fujita et al. (2003). It computes the propagation of ionization fronts around a point source; in our study this is a central starburst source. This code is a post-processing step that operates on a given density distribution at a given time step in our simulation. However, the ionization fronts propagate sufficiently rapidly to adjust almost instantaneously to a changing density distribution.

Figure 21 shows the column density as a function of velocity for Na i-absorbing gas (diamonds) and Hα-emitting gas (triangles) in run U1 at θ = 13° with ionizing luminosities of Q = 10⁵³, 10⁵⁴, 1.9 × 10⁵⁴, and 10⁵⁵ photons s⁻¹. With Q ⩽ 10⁵⁴ photons s⁻¹, the line widths are very small for Hα-emitting gas, Δv_Hα < 100 km s⁻¹, but large for Na i-absorbing gas, $\Delta v_{{\rm Na}\,\hbox{\scriptsize {\sc i}}}=480\;{\rm km\;s}^{-1}$. The densest shell at r = 110 pc traps the photons rather effectively. As Q is increased, that first shell is ionized, but the second shell then traps the photons. The bottom left panel shows that the Na i line width is now very small, $\Delta v_{{\rm Na}\,\hbox{\scriptsize {\sc i}}}=38\;{\rm km\;s}^{-1}$, while Δv_Hα = 330 km s⁻¹. With Q ≳ 2 × 10⁵⁴, all the shells are ionized. In reality, the transition from large $\Delta v_{{\rm Na}\,\hbox{\scriptsize {\sc i}}}$ to $\Delta v_{{\rm Na}\,\hbox{\scriptsize {\sc i}}}\approx 0$ should not be this abrupt because many combinations of fragments and clumps are possible within the observers' 10° × 10° field of view. As we note, we are far from presenting realistic distributions of Na i absorbing and Hα-emitting gas. However, Figure 21 still demonstrates that the line widths of Na i-absorbing and Hα-emitting gas can easily be uncorrelated although they both originate in the swept-up shells.

Zoom In Zoom Out Reset image size

This study clearly is not the final word on this subject. Rather it is a proof of the concept that a starburst wind can accelerate neutral gas up to high velocities without special circumstances, and that the resulting flow configuration can reproduce many of the observable properties of ULIRG winds. We here recap the approximations we have made in order to treat this problem, and add some discussion of how they might be lifted.

The biggest approximation of our study is that we analyze the kinematics of shell fragments at t ≪ 1 Myr in our small grids and extrapolate the results by the ballistic approximation including geometric dilution for comparison with the observations at t ≫ 1 Myr. Much of the total mechanical energy from even an instantaneous burst of star formation has yet to be deposited at an age of 2 Myr when our simulation ends. These dense clumps may further fragment and ultimately be destroyed by R–T and Kelvin–Helmholtz instability as the high-velocity flow of gas from within the bubble streams through them. Conversely, they may be dense enough to be accelerated further prior to their destruction by the ongoing starburst wind.

As the grid resolution of our models improves to 0.1 pc, we do begin to resolve the details of hydrodynamic instabilities. Marcolini et al. (2005) used 0.1 pc resolution for their simulations of the dynamical shredding of radiatively cooling clouds, and seem to have reached adequate resolution. Thus, in future work, continuing to evolve our bubbles in bigger grids with the same resolution may well improve our results.

We have made four other substantial physical approximations as well. First, the assumption of aziumthal symmetry further suppresses R–T instability. Mac Low et al. (1989) pointed out that the typical spike and bubble structure of the R–T instability is limited to rings in an axisymmetric blowout. The detailed distribution of the fragments will certainly be different in three dimensions, but finer fragments will probably have a broader velocity range, and thus are unlikely to change our qualitative result. The behavior of individual fragments will also not be qualitatively changed by the dimensionality, although two-dimensional fragments split into larger pieces initially (Stone & Norman 1992b; Korycansky et al. 2002).

Second, we neglected magnetic fields in this study. We showed that magnetic fields reduce the peak shell density, and thus the fragmentation (Section 4.2.1), but that they don't act to do so before fragmentation is essentially complete in our model, so that this effect is secondary. Magnetic fields can, on the other hand, help preserve individual fragments from further breakup (Stone & Norman 1992b; Mac Low et al. 1994; Shin et al. 2008), possibly even allowing their further acceleration in the continuing starburst wind.

Third, thermal conduction is not explicitly treated in our model. This acts on parsec scales, so numerical diffusion and turbulent mixing will still dominate over thermal conduction under most circumstances at the resolutions that we consider. It is worth noting, though, that studies of individual clump fragmentation have found that thermal conduction can have stabilizing effects (Marcolini et al. 2005; Vieser & Hensler 2007b). Fourth, the dynamical effects of photoionization have also been neglected. This could heat at least low column-density shells and fragments up to 10⁴ K, reducing their density contrast with the wind and enhancing their tendency to fragment.

We also assume a single starburst at the disk center for modeling a galactic outflow. In realistic galaxies, we expect multiple star clusters to be scattered around the disks (e.g., Vacca 1996; Martin 1998). A more complicated structure of cool shells and their fragments is expected with a realistic distribution of star formation as shown by models of supernova-driven turbulence in our own galaxy (e.g., Avillez 2000; Avillez & Berry 2001; Joung & Mac Low 2006). Such a distribution was modeled with a fractal density distribution by Cooper et al. (2008). Fragile et al. (2004) showed that in dwarf galaxies, distributed supernova explosions resulted in increased transfer of energy to the cold gas compared to the centrally concentrated energy injection assumed here and by Mac Low & Ferrara (1999), so our approximation likely gives a lower limit to the kinetic energy of the cold gas. This supports our result of a wide velocity range in Na i-absorbing gas. The next obvious step is to study the effects of realistic star formation on the shell kinematics in three-dimensional, AMR simulations.

Another related issue is that our models assume a low-density, uniform, gas distribution above the site of blowout. However, ULIRGs are formed in galaxy mergers. The actual situation near the nucleus of merging galaxies is likely to be more chaotic. We rely on the idea that the large-scale tidal tails and other merger structures are generally going to be well removed from the wind generation region and will not cover a large fraction of the sky. Ultimately, to remove this limitation, full models of merging galaxies at high resolution will be needed, but this remains some years in the future.

We study the origin and the kinematics of cool gas that produces Na i absorption lines in galactic winds, using hydrodynamic simulations of the blowout of starburst-driven superbubbles from the molecular disk of ULIRGs. The bubbles sweep up the dense disk gas, which quickly cools to form dense, thin shells. The cooled shells fragment by R–T instability as they accelerate in the stratified atmosphere of the disk. This blowout occurs very early in our models, at t ≪ 1 Myr in models with L_mech ⩾ 10⁴¹ erg s⁻¹. The dense fragments left by the R–T instability lag behind the low-density, high-velocity wind. These fragments carry most of the mass that is swept up by the bubbles, and should have the highest column of Na i. The results of our numerical convergence study suggest that superbubble blowout, combined with subsequent geometrical dilution, can reproduce not just qualitative but quantitative properties of the observed lines.

As a result of R–T fragmentation, multiple shell fragments and clumps travel at different velocities. A sightline going through them reproduces the observed broad line width of the Na i absorption profiles: 〈v〉 = 320 ± 120 km s⁻¹. This result does not appear to depend strongly on physical parameters such as mass-loading rate, mechanical luminosity, or surface density. However, this result requires sufficient numerical resolution to follow secondary fragmentation of the shell, so that any line of sight runs over multiple cold gas fragments. We find that a resolution of at least ∼0.2 pc is required to produce this effect, if mass-loading rates remain moderate ξ ≲ 10. The suggestion by Heckman et al. (2000) that swept-up shells will not show a wide range in velocity is based on a hydrodynamic simulation of M82 with a resolution of 4.9 pc (Strickland & Stevens 2000), more than an order of magnitude worse.

The mass-weighted average velocity of cool gas in galactic outflows is found to be ∼400–500 km s⁻¹ among all the models with L_mech = 10⁴³ erg s⁻¹, but ≲200 km s⁻¹ in models with L_mech = 10⁴² and 10⁴¹ erg s⁻¹. The bulk of shells and their fragments are accelerated by the thermal pressure of the bubble interior gas, which is proportional to L_mech. No other parameters influence the results. The mass-weighted average velocities in our high-resolution simulation with L_mech = 10^42–43 erg s⁻¹ agree with the observed value v_s,obs = 330 ± 100 km s⁻¹. A mechanical luminosity of L_mech = 10⁴³ erg s⁻¹ corresponds to an instantaneous burst of M_* = 10⁹ M_☉ or a continuous SFR=500 M_☉ yr⁻¹ at t < 1 Myr after the onset of starburst.

As the swept-up shells fragment by R–T instability, the bubble interior gas blows out and becomes a low-density supersonic wind, as its thermal energy is converted to kinetic energy by expansion. The ram pressure of this wind continues to accelerate entrained cold fragments, and to sweep up high-altitude disk gas, producing small amounts of cold gas with velocities of the order of the terminal velocity of the wind. The correlation between wind terminal velocities and terminal velocities of cool gas in simulations with different mass-loading rates supports this picture. If the wind velocity reaches v_wind ≳ 1000 km s⁻¹, a significant fraction of sightlines through a simulation encounters a terminal velocity close to or greater than the observed average terminal velocity of Na i absorption profiles v_t,obs = 750 km s⁻¹. The mass of cool gas with such high velocity is less than a few percent of that of the total cool gas, however. These outermost shells and fragments with high velocity do fragment further in our highest-resolution simulation, but their fragments are traveling with similarly high velocity. Although the mass in the high-velocity tail decreased as the resolution is increased, we believe it is not likely to go away entirely in a higher-resolution simulation.

The clumps and fragments seen in the simulations may be observed as Na i-absorbing gas or Hα-emitting gas, depending on the amount of ionizing photons produced from the central starburst source. To study this, we used a ray-tracing method to model the location of the ionization front in our models as a function of time, and so to trace the gas emitting in Hα. By varying the photon luminosity, we showed that the velocity range of the ionized and neutral components do not show any correlation with each other. Thus the lack of correlation in the observed line widths of Na i absorbing gas and Hα-emitting gas does not rule out the swept-up shells as the origin of both components.

The hot wind in our model is purely energy-driven. By construction the interior bubble can never become momentum-driven because radiative cooling is turned off in the bubble interior in our models. Future work must determine whether this mechanism reproduces the empirical correlation observed between the maximum outflow velocity and the escape velocity of the host galaxy (Martin 2005). It should also predict the scaling of the mass-loss efficiency with galaxy mass, a key input in cosmological models (e.g., Oppenheimer & Davé 2006, 2008)