STARBURST-DRIVEN GALACTIC WINDS: FILAMENT FORMATION AND EMISSION PROCESSES

The simulations were performed using the Piecewise Parabolic Method with a Lagrangian Remap (PPMLR) code utilized in Paper I. PPMLR is a multidimensional hydrodynamics code based on the method described by Colella & Woodward (1984) and has been extensively modified (e.g., Sutherland et al. 2003a, 2003b) from the original VH-1 code (Blondin 1995). Thermal cooling, based on output from the MAPPINGS III code (Sutherland & Dopita 1993) has been implemented. This allows for a realistic evolution of a radiatively cooling gas (Sutherland et al. 2003b; Saxton et al. 2005). The simulations discussed in this paper are three-dimensional and utilize Cartesian (x, y, z) coordinates. They were performed on the SGI Altix computer operated by the Australian Partnership for Advanced Computing.

In order to study the wind/cloud interaction in sufficient detail, whilst still retaining the ability to follow the evolution and survival of the cloud over a period of approximately 1 million years, the simulations cover a physical spatial range of 50 × 50 × 150 pc, with the cloud centered on the origin. As our intent is to compare these simulations to the formation and survival of the clouds found in Paper I, we choose initial conditions for the wind based upon the results of those simulations (Table 1). In order to understand our choice, we briefly recount the formation and evolution of the starburst-driven winds which we simulated in Paper I.

Table 1. Initial Conditions

Parameter	Symbol	Value
Wind temperature (K)	Tw	5 × 106
Wind density (cm−3)	nw	0.1
Wind velocity (km s−1)	vw	1200
Wind mach number	$\mathcal {M}_{\rm w}$	4.6
Average cloud temperature (K)	Tc	5 × 103
Cloud velocity (km s−1)	vc	0
Cloud radius (pc)	rc	5
Fractal cloud volume (pc3)	$V_{\rm c\_frc}$	1491
Spherical cloud volume (pc3)	$V_{\rm c\_sph}$	523

Download table as:  ASCIITypeset image

• 1.  
  A series of small bubbles of hot (T ≳ 10⁷ K) gas forms in the starburst region. As these bubbles expand, they merge and follow the path of least resistance out of the disk of the galaxy, i.e., the tenuous hot gas surrounding the denser disk clouds.
• 2.  
  As the bubble breaks out of the disk, it begins sweeping up the surrounding halo gas entering the "snow-plow" phase of its evolution. The structure of the wind in the phase is characterized by five different zones: (a) the injection zone, (b) a supersonic free wind, (c) a region of hot, shocked turbulent gas, (d) a cooler "shell" of swept-up halo gas, and (e) the undisturbed ambient gas.
• 3.  
  Clouds of disk gas inside and surrounding the central injection zone are broken up by the freely expanding wind; fragments of disk gas are accelerated into the outflow by the ram pressure of the wind.

Since it is the interaction of the clouds with the freely expanding wind that results in their fragmentation, we select the temperature, density, and velocity of our supersonic wind to be that of the inner free-wind region, namely T_w = 5 × 10⁶ K, n_w = 0.1 cm⁻³, and v_w = 1200 km s⁻¹, respectively.

Our simulations make use of two cloud geometries: a fractal shaped and a spherical shaped cloud. The fractal cloud was chosen in order to allow us to compare the breakup of the clouds in this study to the results of our global model in Paper I and has the same form as the clouds in the inhomogeneous disk used in the global study. The use of the fractal cloud also allows us to investigate the effects of inhomogeneities in the clouds' density distribution on its evolution. A spherical cloud was also modeled in order to allow us to better understand the importance of the assumed cloud geometry. To create the fractal cloud, we first created a 1024 × 1024 × 1024 sized fractal cube using the method described in Sutherland & Bicknell (2007) and Paper I. A single cloud was isolated and extracted from this cube using a blob coloring technique where each cell is examined, and any discontinuous group of nonzero cells given a unique label. An appropriate cloud was then selected and placed at cell number (256, 256, 256) of a 512 × 512 × 1536 sized grid, i.e., the origin of the simulation. The smaller resolution grids used in our fractal simulations were created by downsizing this larger grid (see Table 2). These arrays represent the density of the cloud.

The initial grid was setup by first setting the density (n_w) and pressure (P_w = n_wkT_w/μm_p) to be that of the hot wind. The fractal cloud was then created by adding the density array representing the cloud to the density of the halo gas. To create the spherical cloud, a high-density spherical region of radius r_c = 5 pc was centered on the origin of the computational grid. The radial profile of the density of the spherical cloud is described by an exponential with a scaling radius of 3, in order to mimic the tapered density distribution in the fractal cloud's core. A similar tapered density distribution was considered by Nakamura et al. (2006) for an adiabatic cloud. The boundary condition on the inner z-axis was set to have a fixed inflow with the same properties as that of the hot wind (e.g., T_w = 5 × 10⁶ K, n_w = 0.1 cm⁻³, and v_w = 1200 km s⁻¹). All other boundaries were set to be inflowing/outflowing.

Figure 1 shows the initial density distribution of both the fractal and spherical clouds. The average number density of the fractal cloud is set to be n_c = 63 cm⁻³ and has a total mass of M_c = 1387 M_☉, occupying a volume of $V_{\rm c\_frc} = 1491 \,\rm pc^{3}$. The spherical cloud is setup to occupy a similar volume to the dense core of the fractal cloud, having a radius of r_c = 5 pc and occupying a volume of $V_{\rm c\_sph} = 523 \,\rm pc^{3}$. The average density of the spherical cloud is n_c = 91cm⁻³ and has a total mass of M_c = 523 M_☉. The lower average density of the fractal cloud is due to the large volume of less dense (n_c = 30 cm⁻³) gas that surrounds the cloud core (Figure 1, left panel). In order to understand the effect of the cloud's assumed initial density on its evolution, a simulation in which the density of the fractal cloud was doubled (n_c = 126 cm⁻³) was also performed. In all simulations, the temperature and velocity of each cloud was set to be T_c = 5 × 10³ K and v_c = 0 km s⁻¹, respectively.

Zoom In Zoom Out Reset image size

While these simulations were designed to be applicable to the problem of starburst-driven winds and therefore have densities and temperatures typically found in a such an environment, they are also applicable to other astronomical phenomena that involve the interaction of a supersonic wind with a cloud of gas (e.g., supernova remnants). When cooling is included, the PPMLR code has a one parameter scaling, which is discussed in Sutherland & Bicknell (2007). In general, the density scale of the simulations is inversely proportional to the spatial scale, so that within reason, these simulations can be adapted to problems on both larger and small scales.

The interaction of a supersonic wind with a cloud of dense gas is a problem of some complexity. Whilst there are many factors that could affect a cloud's evolution and survival, such as thermal conduction (Marcolini et al. 2005; Orlando et al. 2005), magnetic fields (Fragile et al. 2005), and photoevaporation (Tenorio-Tagle et al. 2006; see Section 6), here we focus on the importance of radiative cooling and the effect on the clouds initial structure. We have also performed a comprehensive resolution study in order to ascertain the degree as to which the assumed resolution affected our global simulations in Paper I.

We adopt the following naming convention for our simulations: an r or an a indicates whether the simulation includes radiative cooling or is adiabatic, while an f or an s indicates if the geometry of the cloud is fractal or spherical, respectively. The numerical value indicates the number of cells in the x-plane of the computational grid. For example, the simulation rf384 includes radiative cooling, incorporates a fractal cloud and utilizes a computational grid of size 384 × 384 × 1152 cells. An exception to the naming convention is rfd384, which is identical to rf384, but whose cloud is twice as dense.

In total, eleven simulations were performed, with the purpose of each falling within the four different categories outlined below.

• 1.  
  Resolution study. Models rf064, rf096, rf128, rf192, rf256, rf384, and rf512 form our resolution study. These seven simulations all include radiative cooling and utilize a fractal cloud. The resolution of each simulation is given in Table 2 and ranges from 0.79 to 0.10 pc per cell width.
• 2.  
  Cloud structure. Model rs384 has a resolution of 0.13 pc per cell width, includes radiative cooling, and utilizes a spherical cloud. By comparison to rf384, this model is designed to investigate the effect of the shape of the cloud on its breakup and survival.
• 3.  
  Radiative cooling. Models af384 and as384 have the same properties as rf384 and rs384, respectively, but are adiabatic in nature. Both models are designed to help understand the importance of radiative cooling on the evolution and survival of a cloud.
• 4.  
  Cloud density. Model rfd384 is identical to rf384, but has a larger cloud density and mass of n_c = 126 cm⁻³ and M_c = 2770 M_☉, respectively. This model is designed to investigate the effect of the clouds initial density on its evolution.

A summary of the parameters used in each simulation is given in Table 2.

For each simulation, we record the density, temperature, pressure, velocity, emissivity, and a cloud gas tracer in each cell at intervals of 0.01 Myr. With the exception of rf512, each simulation is followed until the time at which the cloud flows off computational grid. In the case of rf512, a computational error occurred while the simulation was in progress. Unfortunately, this simulation is too computationally expensive to rerun, and we are therefore only able to follow the evolution of rf512 to a time t = 0.7 Myr. As such, our resolution study only considers the first 0.7 Myr of the evolution.

Based on their two-dimensional adiabatic simulations, Klein et al. (1994) suggest that a resolution of 120 cells per cloud radius is necessary in order to sufficiently capture the hydrodynamics of the wind/cloud interaction. While this proposed benchmark is easily achieved for a two-dimensional study, it becomes computationally demanding for a fully three-dimensional simulation that includes more complicated physics (e.g., radiative cooling, thermal conduction). In order to overcome this difficulty, other authors have performed two-dimensional axisymmetric simulations (Orlando et al. 2008), assumed symmetry in the solution (Gregori et al. 1999, 2000; Melioli et al. 2005; Orlando et al. 2005), and neglected the more complex physics in their three-dimensional calculations (Orlando et al. 2005). We do not assume symmetry as we show in Section 3 that, even in the idealized case of a spherical cloud, the solution becomes asymmetric as the cloud is broken up and accelerated into the turbulent flow. As the PPMLR code used in this work and Paper I utilizes a uniform grid, the computational resources required to model the entire cloud at the 120 cells per cloud radius resolution suggested by Klein et al. (1994) would be excessive. We are forced to limit the main simulations in this study to a resolution of 38 cells per cloud radius (we consider a resolution of 50 cells per cloud radius in our resolution study). Nevertheless, it has been shown that the global properties of the interaction, as well as the averaged characteristics of cloud ablation process, can be well described at resolutions below this criterion (see, for example, Gregori et al. 2000; Poludnenko et al. 2002; Melioli et al. 2005).

3.1.1. Spherical Cloud

We start with the simple case of the interaction of a spherical radiative cloud (rs384) with a supersonic wind, before describing the evolution of the fractal cloud. In order to illustrate the hydrodynamics of the wind/cloud interaction, throughout this section, we will use slices of the density, temperature, pressure, and velocity through the central y = 0 plane. Figure 2 shows the evolution of the spherical cloud at six different epochs from 0.1 to 1.1 Myr at 0.2 million year intervals. Each panel represents the logarithm of the number density (cm⁻³). The top panels of Figure 3 show the logarithm of the temperature (K), the middle panels show the logarithm of the pressure (cm s⁻²), and the bottom panels show the magnitude of the velocity (km s⁻¹). The evolution at two different epochs is given: 0.35 Myr (left panels) and 0.75 Myr (right panels).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The spherical cloud is initially at rest and is immersed in a hot (T = 10⁶ K), supersonic (V = 1200 km s⁻¹) wind. A bow shock is immediately formed upstream of the cloud. At 0.10 Myr (Figure 2, upper left panel), the "front" of the cloud has been exposed to the high pressure of the shock, while gas is ablated from the "back" of the cloud, a result of the Kelvin–Helmholtz instability and the strong rarefaction that is formed. A high-density shock (n = 1000 cm⁻³) begins to propagate through the cloud, reflecting of the back wall at approximately 0.30 Myr (Figure 2, middle left panel). The Kelvin–Helmholtz instability continues to work, stripping material from the edge of the cloud (e.g., 0.5–1.1 Myr). This material is funneled to approximately 5 pc behind the cloud, where it combines and condenses forming a tail of dense (n ∼ 10 cm⁻³), cool (T ∼ 10⁴ K) cloudlets, which are entrained into the hot turbulent flow downstream from the main cloud (Figure 3, upper left panel).

As the cloud evolves (Figures 2, right panels), the tail of cool gas continues to grow, becoming thicker as the cloud elongates. The cloud material appears as a filament of cool, low-velocity (v < 400 km s⁻¹) gas, immersed inside a region of hot (T > 10⁷ K), turbulent gas with velocity v ∼ 400–1000 km s⁻¹. Small cloudlets continue to be broken off the main cloud, as the Kelvin–Helmholtz instability acts to shed its outer later. However, at 0.75 Myr only 12% of the mass of material that remains on the computational grid is found mixed into the hot wind. The bulk of the mass is still found in the cloud's elongated core and tail, and despite the driving of radiative shocks into the cloud(s) there is little increase in the temperature of the cloud material as it radiates heat, remaining cool, and cohesive as it leaves the computational grid.

3.1.2. Fractal Cloud

Figure 4 shows the evolution of the radiative fractal cloud in model rf384 from 0.1 to 1.1 Myr through density slices at intervals of 0.2 Myr. Figure 5 is identical to Figure 4, but shows the evolution of a cloud (rfd384) with twice the density and mass than rf384. As with the radiative spherical cloud, a bow shock immediately forms upstream of the cloud. However, a significant effect of the inhomogeneous structure of the cloud is the formation of a shock off each ridge on the cloud's surface that is exposed to the wind. This has the effect of creating a "web" of interacting shocks.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By 0.3 Myr (Figure 4, middle left panel), the density at the front of the cloud has jumped to n ∼ 1000 cm⁻³, but fallen to approximately n = 10 cm⁻³ at the rear of the cloud. Despite the high pressure exerted on the cloud, very little heating occurs, with the cloud maintaining temperatures of T ∼ 10³–10⁴ K (Figure 6, left upper two panels). As the cloud evolves (Figure 4, right panels), small cloudlets are broken off the main cloud by the Kelvin–Helmholtz instability. These cloudlets form a filamentary structure downstream of the main cloud and have velocities in the range of 0–400 km s⁻¹, somewhat slower than the velocity of the surrounding stream (v > 800 km s⁻¹) (Figure 6, bottom panels).

Zoom In Zoom Out Reset image size

The cloud is exposed to the high temperature and velocity of the wind at all points along the front surface of the cloud. The cloud breaks up fastest in regions where the density is lowest and the radius of the cloud is at its smallest. The cloud is continually eroded by the Kelvin–Helmholtz instability, the fragments of which are immersed in a low pressure, turbulent gas. At 0.75 Myr, the percentage of the mass of cloud material remaining on the computational grid that has mixed into the hot wind is ∼25%. The bulk of the original cloud mass is found in the remaining core remnant and the stream of small, dense ∼1 pc cloudlets. If these cloudlets are exposed to the wind they produce their own high-pressure bow shock upstream of their position. However, the majority of the cloudlets are sheltered from the wind by other fragments broken off from the main cloud. These cloudlets survive to leave the computational grid.

The higher density of the cloud in model rfd384 (Figure 5) results in the cloud retaining its structural integrity far longer than the lower density cloud in rf384. While the evolution of the higher density cloud is overall very similar to that in rf384, the breakup of the cloud via the Kelvin–Helmholtz instability is slower: at 1.1 Myr (bottom right panel), a significant bulk of the cloud material still remains on the computational grid and displays a similar morphology to the cloud in rf384 at 0.7 Myr (Figure 4, top right panel). In contrast to rf384, at 0.75 Myr only 11% of the clouds mass remaining on the computational grid is found mixed into the hot wind. However, it is likely that this amount will increase as the cloud is further eroded by the Kelvin–Helmholtz instability.

As observed by Xu & Stone (1995), the initial shape of the cloud has an effect on its subsequent evolution, with our fractal cloud fragmenting faster than the spherical cloud. While this is due in part to the lower average density of the fractal cloud, even the high-density fractal cloud, whose average density and total mass is greater than that of the spherical cloud (see Table 2), has a greater degree of fragmentation and is less cohesive when it leaves the computational grid. This is a result of the inhomogeneous nature of the fractal cloud's initial density distribution and the larger cross section it presents to the incident wind (Figure 1, left panel). The cloud first fragments along regions where the wind finds paths of least resistance, i.e., regions of low density. As a result, the fractal cloud breaks up into multiple core fragments. As the wind finds no regions of least resistance in the spherical cloud, it is able to retain a single cohesive structure for a longer period of time. Thus, not only is the initial geometry of the cloud important in determining its evolution, the distribution of the cloud's density determines how quickly the cloud begins to fragment. More homogeneous density distributions would result in less initial fragmentation.

3.2.1. Adiabatic Case

In order to understand the degree as to which the inclusion of radiative cooling affects the survival of a cloud, we performed two simulations in which cooling was neglected: as384 and af384 for the spherical and fractal case, respectively. Figure 7 shows the density of the adiabatic clouds in as384 (left) and af384 (right) at 0.3, 0.6, 0.9 Myr epochs. It can be immediately seen that in the absence of radiative cooling there is a far greater degree of mixing of cloud material with the surrounding stream, with 40% and 59% of the fractal and spherical clouds' masses respectively on the computational grid found mixed into the hot wind at 0.75 Myr.

Zoom In Zoom Out Reset image size

The destruction of the adiabatic fractal cloud occurs faster and is more complete than the adiabatic spherical cloud. At 0.9 Myr, the cloud has been almost completely destroyed (Figure 7, bottom right panel), with the cloud material mixed into the hot gas. The initial interaction with the wind is the same as in the radiative case, with a bow shock forming upstream of the cloud. However, the cloud gas quickly begins to heat to temperatures of the order of T ∼ 10⁶ K (Figure 8, upper left panel) and the cloud expands. Again the cloud breaks up first in regions where the density is low and the cloud radius is at its smallest. The high pressure exerted on the cloud (Figure 8, middle left panel) and the Kelvin–Helmholtz instability act to strip material from the cloud. While this material is able to survive in the radiative model, here the cloudlets broken off the main cloud are quickly heated and destroyed. The bulk velocity of the gas downstream of the main cloud is lower than that found in the radiative case (Figure 8, lower panels).

Zoom In Zoom Out Reset image size

In the case of the adiabatic spherical cloud, the shedding of the cloud's outer layer by the Kelvin–Helmholtz instability, also seen in the radiative model, is greatly enhanced. By 0.9 Myr (Figure 7, bottom left panel), only the cloud core remains and is subsequently destroyed by the Kelvin–Helmholtz instability as the simulation progresses. As with the fractal cloud, the bulk velocity of the gas downstream of the main cloud is lower than that of the radiative spherical cloud at the same time (Figure 9, lower right panel). The evolution of a spherical adiabatic cloud is discussed in more detail below.

Zoom In Zoom Out Reset image size

3.2.2. Cloud Survival

One of the significant effects of the inclusion of radiative cooling is the longer lifetime of the impacted cloud. While this effect has been observed in the past by other authors (e.g., Mellema et al. 2002; Melioli et al. 2005), our study is of higher resolution and does not assume any symmetry, making it is a useful exercise to directly compare the simple case of the evolution of the spherical cloud in both our adiabatic (as384) and radiative (rs384) models in order to determine the mechanism behind the radiative cloud's survival. The initial interaction of the wind and cloud is shown via density slices in Figure 10 at 0.05, 0.20, and 0.35 Myr epochs in both adiabatic (left) and radiative (right) models. It can be clearly seen that while the evolution begins almost identically with a bow shock forming upstream of the cloud, in the adiabatic case cloud material immediately starts being ablated from the back of the cloud.

Zoom In Zoom Out Reset image size

As in the radiative model, a shock propagates through the adiabatic cloud. However, the initial density increase observed in the cloud is not as extreme (e.g., n < 1000 cm⁻³). The cloud material is heated to temperatures of the order of T > 10⁵ K and the cloud expands. The shock travels though the cloud, reflecting off the back surface at approximately 0.27 Myr, somewhat faster than in the radiative cloud. The adiabatic cloud expands transversely as it is accelerated downstream. However, this transverse expansion is suppressed in the radiative cloud as a result of the lower degree of heating of the cloud gas (Figure 3, upper panels). In both cases, the Kelvin–Helmholtz instability acts to strip material from the edges of the cloud forming a tail of material downstream of the cloud position (Figure 10, middle and lower panels). In the adiabatic model, this tail is geometrically thick with density and temperature n ∼ 1 cm⁻³ and T ∼ 10⁶ K, respectively. In contrast, the tail formed in the radiative model is geometrically thin with density n ∼ 10 cm⁻³ and temperature T ∼ 10⁴ K.

In the adiabatic model, the internal cloud shock reflects again off the front the cloud at approximately 0.35 Myr. At this time, the transverse expansion of the cloud persists and the Kelvin–Helmholtz instability continues to strip material from the clouds exterior into the downstream tail of gas. (Figure 10, lower left panel). The transverse expansion (Figure 9, middle panels) results in a higher rate of acceleration in the adiabatic model. As a consequence, the cloud has a lower relative Mach number relative to the stream (Figure 9, lower left panel) than the radiative cloud (Figure 3, lower left panel). The growth rate of the Kelvin–Helmholtz instability is lower for higher Mach numbers and its effect is strongly diminished for the radiative cloud. The adiabatic cloud is more easily disrupted and destroyed. This can be dramatically seen in the middle left panel of Figure 7 where the Kelvin–Helmholtz instability has stripped the entire outer layer of the adiabatic spherical cloud (see also Figure 3 of Orlando et al. 2005).

Since a radiative cloud is broken up via the Kelvin–Helmholtz instability into a filamentary structure of small ∼1 pc sized clouds, the survival of these small clouds is of interest. We now compare the cloud crushing time (t_crush) and the Kelvin–Helmholtz timescale (t_KH) to the cooling time (t_cool) of a cloud with radius R_c = 1 pc, density ρ_c = 10 cm⁻³, temperature T_c = 10⁴ K, and velocity v_c = 200 km s⁻¹. The cloud crushing time of such a cloud is t_crush ≈ R_c/v_sh ≈ (ρ_c/ρ_w)R_c/v_w = 3 × 10¹² s, where the density and velocity of the wind is ρ_w = 0.1 cm⁻³ and v_w = 1000 km s⁻¹, respectively. The Kelvin–Helmholtz timescale is t_KH = R_c(ρ_c + ρ_w)/(v_c − v_w)(ρ_cρ_w)^1/2 = 3 × 10¹¹ s. The cooling time for a 1 pc cloud in our simulations is of the order of 10¹⁰ s, somewhat shorter than the cloud crushing time and the Kelvin–Helmholtz timescale, suggesting that the cloudlets may remain sufficiently stable to ablation and survive to later times. In addition, we note that the self-gravity of the clouds may cause them to collapse, becoming more difficult to disrupt (Mellema et al. 2002).

In order to test the dependence of our results on the numerical resolution of the code, we have performed seven simulations of the radiative fractal cloud interacting with a supersonic wind at increasing resolutions from 0.78–0.10 pc per cell width (see Table 2). We calculated the flux of mass through a surface at z = 75 pc over the first 0.7 million years of each simulation in our resolution study. Figure 11 shows this mass flux for each simulation as a function of time (left panel), as well as the total mass flux integrated over the first 0.7 Myr of the simulation as a function of resolution (right panel). The mass flux over the first 0.4 Myr of the evolution is similar at all resolutions. This is a result of the well resolved hot stream of gas passing through the flux surface. After this point, the cloud material begins to pass though the surface and the mass flux starts to vary with resolution.

Zoom In Zoom Out Reset image size

The initial drop in the mass flux seen in the left hand panel of Figure 11 is due to the rarefaction that passes though the flux surface (see Figure 2, top left panel). From approximately 0.2 to 0.4 Myr, the density of the stream increases as mass is ablated from the rear of the cloud resulting in a similar flux at all resolutions. During this time, the mass flux gradually increases as the n = 0.1 cm⁻³ mass loaded stream of gas passes through the flux surface. At approximately 0.4 Myr, the mass flux begins to vary rapidly, dramatically increasing in each simulation. The large variation in the mass flux is due to the turbulent nature of the gas passing though the flux surface and the dense cloudlets immersed within this gas. The cloudlets pass through the surface at different times resulting in fluctuations in the mass flux.

There is a general trend of increasing mass flux with the resolution of the simulation. This trend can be explained by the increase in fragmentation of the cloud with increasing numerical resolution. Figure 12 shows volume renderings of the projected density at 0.7 Myr in models rf064, rf128, rf256 and rf512, which have resolutions of 0.78, 0.39, 0.20, and 0.10 pc per cell width, respectively. The increase in the fragmentation of the cloud with resolution can clearly be seen, and will be discussed in more detail in Section 4.2. At the highest resolution attempted in this study (0.1 pc), the filaments resemble a "foam" of cloudlets, while at low resolution the cloud has been broken up into only a few large fragments. As a consequence, the cross section of dense material that passes through the flux surface at any given time after 0.4 Myr is larger in the higher resolution simulations.

Zoom In Zoom Out Reset image size

The total integrated mass flux passing through the surface at z = 0.75 pc over the first 0.7 Myr also increases with numerical resolution (Figure 11, right panel). Again, this is caused by the larger degree of fragmentation at high resolution, resulting in more ablation of cloud material from the back of the cloud. Between our highest and lowest resolution simulations the difference in mass flux is approximately 10%. This discrepancy is likely to increase at higher resolutions, although it is possible that convergence may occur at extremely high-resolution simulations (>0.1 pc) that utilize an adaptive mesh. For the 2 pc per cell width resolution of our global simulations in Paper I, this error is increased to approximately 20%. We will discuss the impact of the numerical resolution on the results of Paper I in Section 5.

The most significant effect of increasing the numerical resolution of the simulation is the increase in fragmentation of the cloud. The increase in cloud fragmentation can clearly be seen in Figure 12, with the cloud broken into only a few large fragments in the lowest resolution simulations, but hundreds of fragments at higher resolution. Nevertheless, filamentary structure, where the concentration of cloudlets is higher, can still be made out at high resolution. These "filaments" are located in similar positions to the filaments in the low-resolution simulations.

We are able to calculate the properties of each cloudlet by using an algorithm which allows us to pick out and select fragments. Note that we impose a minimum mass of M_c = 10⁻³ M_☉ for a fragment to be selected. Figure 13 shows the number of cloudlets produced as a function of numerical resolution. There is a general trend from the number of cloudlets produced in the interaction to increase as a power law with increasing resolution. Even with larger computational resources and an adaptive mesh, this trend is likely to continue ad infinitum. The dependence of the degree of fragmentation on the numerical resolution has been observed by other authors in both two-dimensional (Klein et al. 1994) and three-dimensional (Stone & Norman 1992) simulations of a spherical cloud. This effect can be explained by the growth rate of Kelvin–Helmholtz instability, which is faster at smaller wavelengths. As the resolution is increased, this instability is increasingly resolved and more fragmentation of the cloud is observed. This is illustrated in Figure 14, where the Kelvin–Helmholtz instability can clearly be seen to increasingly fragment the cloud at higher numerical resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 15 gives velocity (left) and mass (center) histograms of the cloudlets, as well as the total mass of the cloudlets as function of the cloudlet mass (right), at 0.7 Myr for the four resolutions shown in Figure 12. The massive fragments (M ≳ 10²–10³ M_☉) present at all resolutions are the remnants of the cloud core. Since the number of cloudlets increases with resolution, the highest resolution simulations are comprised of numerous low mass cloudlets. In general, there is a trend toward increasing smaller mass fragments at all resolutions. However, in all cases, the majority of mass of the cloudlet system is found in the massive fragments (>40 %) rather than the smaller fragments. The velocity of the cloudlets does not vary significantly with resolution. Despite the large number of lower mass fragments present in the high-resolution simulations, they still fall within the velocity range of v_c = 150–400 km s⁻¹. The bulk of the cloudlets at all resolutions have a velocity in the range of v_c = 180–220 km s⁻¹. This is likely to increase as the evolution progresses and the cloudlets are further accelerated by the wind.

Zoom In Zoom Out Reset image size

Unfortunately, we are unable to fully resolve the interaction of a radiative cloud with a supersonic wind at this time. This limits our ability to draw any reliable conclusions regarding the small-scale evolution of the wind/cloud interaction. A significant increase in numerical resolution and/or the use of an adaptive mesh would be required in order for convergence to possibly occur. Nevertheless, there is a clear trend in the large-scale evolution of the wind/cloud interaction at all resolutions considered in this study, allowing us to draw some physical conclusions. For example, the soft X-ray luminosity is sufficiently resolved and is almost constant with numerical resolution (see Section 5.2). The effect of radiation in keeping the cloud cool, suppressing the transverse expansion, and minimizing the effect of the Kelvin–Helmholtz instability, is also not affected by the resolution of the simulation. In all cases, the cloud is not immediately destroyed and mixed into the hot wind as seen in adiabatic models, but is instead broken up into numerous small cloudlets. The major effect of increasing the resolution is the increased fragmentation of the cloud. However, we still see the same qualitative structure at high resolution, with the cloud breaking up to form a filamentary structure that becomes finer and finer as more detail is resolved.

At optical wavelengths (such as Hα), starburst-driven winds appear as spectacular filamentary systems extending several kpc along the minor axis of the host galaxy, e.g., M82 (Shopbell & Bland-Hawthorn 1998), NGC 3079 (Veilleux et al. 1994), NGC 1569 (Westmoquette et al. 2008). While it has long been proposed that this filamentary material was expelled from the central region of the galaxy (Lynds & Sandage 1963; Bland & Tully 1988), until now the mechanism behind the formation of the filaments has not been completely understood. In Paper I, we proposed that the filaments were formed via clouds of disk gas that are brokenup and accelerated into the outflow by the ram pressure of the wind. In order for this mechanism to be viable, the cloud fragments need to survive and remain sufficiently cool to emit at Hα temperatures. The simulations presented in this paper allow us to address this important issue.

As in Paper I, we define the Hα emitting gas to be cloud material with temperatures in the range of T = 5 × 10³ − 3 × 10⁴ K. We note that as photoionization is not included in our model, the Hα emission discussed in this paper arises solely from shock ionization. However, photoionization is known to play a role in the ionization of the filaments in many winds. For example, the filaments in M82's wind are known to be photoionized at low distances above and below the galactic plane, with shock ionization becoming dominant at large distances. An investigation into the effects of photoionization on the cloudlets is warranted, but is beyond the scope of this study. Figure 16 shows a three-dimensional volume rendering of the density of the Hα emitting gas at 0.5 Myr in model rf384 (online: movie of evolution of the Hα emitting gas in rf384 over the first 1.37 Myr). It can be immediately seen in Figure 16 that the Hα emitting material corresponds to the dense cloud material. Thus, the survival mechanism for a cloud proposed in Section 3.2.2 can be invoked to explain the filaments observed in starburst-driven winds.

Zoom In Zoom Out Reset image size

As discussed in Section 3, the cloud is broken up via the Kelvin–Helmholtz instability, with the fragments subsequently entrained into the outflow forming a filamentary structure. Figure 17 (left panel) gives the emission weighted histogram of the z-velocity along the filament at 0.7 Myr in model rf384. The majority of the Hα emission has velocities in the range of v ∼ 0 − 30 km s⁻¹, with the velocity increasing with distance along the z-axis. This is consistent with the Hα gas reaching higher velocities at larger distances above the galaxy plane, which was observed in the global simulations in Paper I. Note that this result represents only the velocity dispersion at the base of a single filament early in its evolution, with the velocity likely to increase as the cloudlets, which form the filament, are further accelerated in the direction of the flow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As in our global simulations, we again find that it is the ram pressure of the wind that accelerates the clouds. The main difference between the filaments formed in Paper I and the filaments found here is the number of cloudlets that comprise the filament; a direct result of the higher numerical resolution of the simulations in this work. If we were able to increase the resolution of the global simulations, we would find similar fine structure in the filaments that are formed. Figure 18 (left) shows the mass flux of the Hα emitting gas through a surface at z = 75 pc at each resolution considered in our study. As expected, there is no Hα emitting gas passing though the surface until approximately 0.3 Myr, as the dense cloud material has yet to encounter the flux surface. As in Figure 11, the difference in the mass flux at each resolution is a result of the increasing fragmentation of the Hα emitting clouds at high resolution. The right-hand panel of Figure 18 shows the total integrated mass flux, over the first 0.7 Myr, of the Hα emitting gas as a function of resolution. We see an increase in the total mass passing through the flux surface, again caused by the increase in fragmentation as the Kelvin–Helmholtz instability is further resolved. Even at the high resolution of these simulations the Hα mass flux is yet to converge.

In Section 3.2.2, we discussed how a cloud could survive the interaction with a hot supersonic wind. The ability of a cloud to radiate heat is crucial for the clouds survival, allowing it to remain stable to ablation and emit at Hα temperatures. Without this ability, the cloud quickly heats above T = 10⁶ K, expands and becomes susceptible to the Kelvin–Helmholtz instability. While a radiative cloud is still disrupted, the small sized fragments that are broken off the main cloud have cooling times faster than the cloud crushing time and the Kelvin–Helmholtz growth rate, and thus, may possibly survive. In the adiabatic case, the fragments get heated and destroyed, with the cloud material quickly becoming mixed into the hot wind. If the fragments survive, they are drawn out into strands and form a filament downstream of the original cloud position, reminiscent of the filaments seen in starburst-driven winds.

In Paper I, we proposed four mechanisms that could give rise to the soft X-ray emission that is observed to be spatially correlated to the Hα emitting filaments. This correlation has been observed by Chandra in many starburst-driven winds (e.g., Cecil et al. 2002; Strickland et al. 2004a, 2004b). Here, we summarize the proposed mechanisms.

• 1.  
  The mass-loaded wind. This is the largest contributor to the soft X-ray emission in the global simulations. As mass is ablated from the clouds it is mixed into the surrounding turbulent gas, creating a region of hot (T ≳ 10⁶ K) rapidly cooling gas that emits strongly at soft X-ray energies.
• 2.  
  The intermediate temperature interface between the hot wind and cool filaments. Gas at the boundary between the hot and cool gas mixes to produce a thin region of intermediate density and temperature. Like the mass-loaded wind, this mixed gas is a strong emitter of soft X-rays.
• 3.  
  Bow shocks. Gas is heated to X-ray temperatures as a bow shock is formed upstream of dense clouds accelerated into the flow.
• 4.  
  Colliding bow shocks. When two bow shocks interact, the gas is further shock heated to X-ray temperatures.

The first two processes involve the mixing of hot and cold gas and could be the result of numerical diffusion in the simulations and therefore not physical. Our resolution study allows us to examine the effect of increasing resolution on the soft X-ray emission and determine the realism of the above possible emission processes.

As in Paper I, we infer the X-ray luminosity in the soft (0.5–2.0 keV) energy band using broadband cooling fractions obtained from the MAPPINGS IIIr code (Sutherland & Dopita 1993). Figure 19 shows the soft X-ray emissivity in models rf128 (left) and rf512 (right) at 0.3, 0.5, and 0.7 Myr epochs. The strongest X-ray emitter in both models is the bow shock that immediately forms upstream of the cloud as it interacts with the wind. Regions where bow shocks are interacting result in the highest X-ray emissivities. Apart from a few marginally bright tails coming off some of the cloudlets, particularly in model rf128 (see Figure 19, bottom left panels), we see no evidence that mass loading of the wind by ablation from the cloud is a significant contributor to the soft X-ray luminosity. We also see very little evidence that the intermediate temperature interface plays a significant role. There are a few bright regions upstream of cloudlets that have been exposed to the wind. However, it is likely that this enhanced emission is related to the X-ray emitting bow shocks that have formed around each cloudlet.

Zoom In Zoom Out Reset image size

The main difference between the simulations shown in Figure 19 at high (right panel) and low (left panel) resolution is in the structure of the main bow shock and the emission from the cloudlets. While we observe some structure in the low-resolution simulations, as the resolution is increased we see clear regions where colliding bow shocks lead to a significant increase in the X-ray emissivity. This is more evident at later times when the main cloud has fragmented and there are many X-ray emitting bow shocks upstream of the resulting cloudlets (Figure 19, bottom right panels). It might be expected then that the increase in fragmentation at higher resolutions would lead to an increase in the X-ray luminosity as more bow shocks are formed. However, the majority of the cloudlets formed are sheltered from the impacting wind by the main cloud and do not form a bow shock. Thus, they are not seen at soft X-ray energies and do not contribute to the X-ray luminosity (Figure 19, right panels).

Indeed, the soft X-ray luminosity hardly varies with resolution. Figure 20 shows the soft X-ray luminosity as a function of resolution for our radiative fractal cloud at 0.7 Myr. At all resolutions, the luminosity is of the order of L_x ∼ 10³⁶ erg s⁻¹. This amount varies only negligibly from low to high resolution. This is a result of the strongest X-ray source being the main bow shock upstream of the original gas cloud, which is present at all resolutions. As discussed above, we see little evidence that the mass-loaded component of the wind, nor the intermediate temperature interface between the hot and cold gas plays a large role in the soft X-ray emission. While high resolution (<0.5 pc/cell) global simulations are needed in order to confirm this result, we conclude that bow shocks and their interaction are the main source of the soft X-ray emission in starburst-driven winds.

Zoom In Zoom Out Reset image size

The importance of radiative cooling in the formation and survival of the filaments was discussed in Section 3.2.2. For our proposal of the formation of a filament, via the breakup and acceleration of cool disk gas into the wind, to be a viable mechanism, cooling must be present in the outflow. Observations of the O vi emission line could be used to detect this cooling in the filaments. Heckman et al. (2001) report the detection of O vi emission in the dwarf starburst galaxy NGC 1705, which they associate with cooling in the outflow from gas at temperatures of T ≳ 3 × 10⁵ K. They propose turbulent mixing layers (e.g., Slavin et al. 1993) as a possible origin for this emission.

In order to determine where the O vi emission may arise in our simulations, we have produced a "map" of the predicted O vi emission in models rf384, rs384, af384, and as384 at 0.7 Myr (Figure 21). We assume that the O vi emission in our simulations falls within the temperature range of T = (1 × 10⁵) − (4 × 10⁵) K. When cooling is neglected (bottom two panels), O vi emission is only observed around the surviving cloud core. This emission is the result of the high degree of mixing of the hot and cold gas seen in the adiabatic simulations. In the radiative models (top two panels), O vi emission is observed throughout the flow and is closely aligned to the filamentary gas. This emission is caused by the mixing of hot and cold gas in the vicinity of each cool cloudlet that comprises the filament.

Zoom In Zoom Out Reset image size

The distribution of the O vi emission is most significant in the radiative fractal model (rf384; top panel), where the structure of the filament can still be clearly seen. The emission weighted histogram of the z-velocity for this model at 0.7 Myr is shown in the right-hand panel of Figure 17. The velocity dispersion of the O vi gas is similar to that of the Hα emitting filaments, falling in the range of v ∼ 0–40 km s⁻¹. While our investigation of the O vi emission is only preliminary, further detection of O vi kinematics similar to those proposed in this work would lend support to the premise of cooling in the filaments of starburst winds.

The effects of thermal conduction on the wind/cloud interaction has been investigated over the last decade by a number of authors (e.g. Vieser & Hensler 2000; Hensler & Vieser 2002; Marcolini et al. 2005; Orlando et al. 2005; Recchi & Hensler 2007). The immediate effects on the evolution of a radiative cloud has been described in detail by Marcolini et al. (2005) who modeled the interaction of a cloud of radius R_c = 15 pc and temperature T_c = 10⁴ K with a hot wind of temperature T_w = 10⁶ K. Here, we summarize the evolution of a radiative cloud in their two-dimensional simulation that includes the effects of thermal conduction.

• 1.  
  In the early phases of the cloud's evolution, thermal conduction results in a converging shock forming around the cloud. The cloud shrinks in size until the central pressure in the cloud's core is high enough to halt the collapse.
• 2.  
  The cloud begins to re-expand and over time becomes elongated in the downstream direction.
• 3.  
  As the cloud evolves, it is compressed in a nonuniform manner due to the shape of the cloud and the complexity of the flow.
• 4.  
  After 1 Myr, the cloud forms a filamentary structure along the symmetry axis of the simulation.

The key difference to the evolution of a cloud in a simulation that neglects thermal conduction is that the cloud does not suffer hydrodynamical instabilities and fragment. This is because thermal conduction smooths the density and velocity gradients across the surface of the cloud, inhibiting the growth of instabilities (also see Vieser & Hensler 2000). However, Marcolini et al. (2005) note that there is also significant mass loss from the cloud through evaporation in their thermal conduction models.

Orlando et al. (2005) also considered the effects of thermal conduction on a cloud's evolution by investigating the survival of a cloud of radius R_c = 1 pc being overrun by Mach 30 and 50 shocks. They performed both three-dimensional simulations that ignored radiative cooling and thermal conduction and two-dimensional simulations where they were included. Similar to Marcolini et al. (2005), they found that thermal conduction inhibited the growth of hydrodynamical instabilities. In their simulations, the structure of a cloud consists of a cool dense core surrounded by a "corona" of gas, dominated by thermal conduction, that gradually evaporates as the cloud evolves. The cloud does not fragment and there is significant mass loss through evaporation. In the case of the slower Mach 30 shock, the dense core is unaffected by heat conduction and still fragments into multiple cloudlets. They also investigated the X-ray emission that would arise from their simulations (Orlando et al. 2006). They found that for the slower shock, the soft X-ray emission arose in the thermally conducting corona, while the cold core would likely emit at optical wavelengths. In the case of the faster Mach 50 shock, they suggest that there would be no optical component.

In this paper, we find that the soft X-ray emission arises primarily from the main bow shock upstream of the original cloud. The optically emitting filamentary gas originates from the cool cloudlets that are the remnants of this cloud. As the simulations of Marcolini et al. (2005) are the closest in initial conditions to our study, we use them as a basis for our supposition on the effect of thermal conduction on our results. Most significantly, thermal conduction would likely act to inhibit the growth of the Kelvin–Helmholtz instability, which is responsible for the high degree of cloud fragmentation seen in our simulations. While the initial evolution of the cloud would be similar to that discussed in Section 3, the breakup of the cloud would likely be totally or partially suppressed. The initial structure of the cloud (e.g., fractal or spherical) would play a role in determining the degree of any fragmentation that occurs. A more cohesive Hα emitting filamentary structure would be formed. In their study, Marcolini et al. (2005) reported the bow shock upstream of the initial cloud to be the main source of the soft X-ray emission in their nonthermal conducting simulations. When thermal conduction was considered, they found an increase in the soft X-ray emission in the region between the cloud surface and the bow shock. It is likely that our simulations would show a similar result. Note that this emission is still complimentary to the filamentary Hα emission, as seen by observations of starburst winds.

We conclude that thermal conduction may help the survival time of a radiative cloud overrun by a supersonic wind by stabilizing it against the destructive effect of the Kelvin–Helmholtz instability. While there will be less mixing of cloud material into the wind by hydrodynamical effects, there will still be significant mass loss from the cloud via evaporation. However, as discussed below, the presence of a magnetic field may counteract any benefits to cloud stability gained through thermal conduction.

Magnetic fields have been shown by a number of authors to have a varying effect on the wind/cloud interaction. A variety of different circumstances have been taken into consideration, such as the strength and orientation of the magnetic field. Early models were largely two-dimensional and ignored the effects of radiative cooling (e.g., Mac Low et al. 1994), while later models have been three-dimensional (e.g., Gregori et al. 1999, 2000; Shin et al. 2008), and/or included radiative cooling and thermal conduction (e.g., Fragile et al. 2005; Orlando et al. 2008). These studies have shown that the orientation of the magnetic field (i.e., parallel, perpendicular, or oblique to the flow) has a significant effect on the evolution of the cloud. For example, Shin et al. (2008) found that in the case of a strong oblique magnetic field, the cloud was actually pushed out of the central y = 0 plane as it evolved. It has also been suggested that magnetic fields can help suppress the hydrodynamical instabilities that shred the cloud in a nonmagnetized simulation (Mac Low et al. 1994; Fragile et al. 2005), with less fragmentation occurring with increasing field strength. On the other hand, Gregori et al. (1999) show that in a three-dimensional model, the growth of hydrodynamical instabilities can actually be accelerated by the presence of a magnetic field.

It is widely thought that the presence of a magnetic field will act to suppress thermal conduction (e.g., Chandran & Cowley 1998; Narayan & Medvedev 2001). Marcolini et al. (2005) suggest that the coefficient of thermal conductivity κ may be significantly reduced below the Spitzer level in the presence of a magnetic field. To investigate how a reduced degree of conductivity may effect their results, they performed a number of simulations, where κ was below the standard Spitzer level κ_sp = 6.1 × 10⁻⁷T^5/2 erg s K⁻¹. They found that at sub-Spitzer levels, thermal conduction is not as efficient at suppressing the hydrodynamical instabilities that fragment the cloud. More recently, Orlando et al. (2008) ran a series of two-dimensional axisymmetric simulations that included radiative cooling, thermal conduction, and magnetic fields. They showed that while thermal conduction was not completely suppressed in the presence of a magnetic field, regardless of the orientation of the field, the cloud would be broken up and destroyed by hydrodynamical instabilities. However, they also found that this effect was lessened in the case of a strong magnetic field.

While clearly the interplay between thermal conduction and magnetic fields will influence the evolution and survival of the filaments in starburst winds, without detailed three-dimensional simulations (with initial conditions appropriate to the problem) that consider radiative cooling, thermal conduction, and the complex nature of magnetic fields, it is difficult to hypothesize what effect the magnetic field would have on a filament's formation and survival. Nevertheless, a filamentary structure is still likely to be formed by the wind/cloud interaction, with the most significant difference being in the number of cloudlets that comprise the filament. We do not anticipate that there would be a large effect to the Hα, soft X-ray, and O vi emission processes discussed in Section 5.