THE COMPLEX STRUCTURE OF THE MULTI-PHASE GALACTIC WIND IN A STARBURST MERGER

Galactic winds, driven by supernova-heated gas (Chevalier & Clegg 1985) or active galactic nuclei (AGNs; Veilleux et al. 2005), have been frequently observed in starburst galaxies (e.g., Lehnert & Heckman 1995). They play a significant role in enriching the intergalactic medium with metals and suppressing the star formation and growth of super massive black holes (Veilleux et al. 2005, and references therein). They are an indispensable element in theoretical models of galaxy evolution, in which they provide the feedback mechanism that enables the models to reproduce galaxy scaling relations consistent with those observed (e.g., Dutton & van den Bosch 2009). Winds have also been invoked to explain a number of issues related to contemporary cosmology (e.g., Veilleux & Rupke 2004), such as the discrepancies between the observed and predicted galaxy luminosity functions (e.g., Springel & Hernquist 2003).

Ultraluminous infrared galaxies (ULIRGs, defined to have L_IR ⩾ 10¹² L_☉) are the late stages of a merger of two gas-rich galaxies, which may then form a QSO (Sanders et al. 1988). Though rare at low redshift, their co-moving density increases rapidly at high redshift (Le Floc'h et al. 2005). They are hosts to intense starbursts and may dominate the star formation in the universe at z ⩾ 1 (Le Floc'h et al. 2004; Caputi et al. 2007). Detailed absorption-line studies of a large number of ULIRGs have shown that winds are ubiquitous in these systems (Heckman et al 2000; Rupke et al. 2002, 2005a, 2005b, 2005c; Martin 2005, 2006) and may play a crucial role in their evolution.

Our current knowledge of ULIRG winds, especially their geometry, is very limited. Though previous studies of ULIRGs have yielded valuable insights (including velocity distributions, detection rates, and rough mass estimates), they give little or no information on the spatial variation of wind properties (though see Martin 2006 for a one-dimensional study). Previous work on estimating properties of ULIRG winds assumed spherical symmetry of the wind with uniform velocity and density (Rupke et al. 2005b). Since ULIRG winds occur in merging systems that have more complicated dynamics than isolated disk galaxies, interaction between two galaxies can possibly affect the outflow. For example, Rupke et al. (2005b) found that winds are detected twice as often in ULIRGs as in LIRGs, which suggests that the ULIRG winds may have a wider opening angle compared to winds from less luminous galaxies. If this is true, then a spherically symmetric wind model may be too simplistic to represent ULIRG winds.

To gain a better understanding of ULIRG winds, we obtained two-dimensional spectra of F10565+2448 with the integral field unit (IFU) on Gemini North. Outflows in our target have been studied with one-dimensional spectra. Martin (2006) has studied the velocity variations in this galaxy via long-slit spectra, and Rupke et al. (2005b) estimated the mass, energy, and momentum in the wind. Using the IFU data we can: (1) measure how extended the outflow is, (2) determine whether the basic model for a disk wind (Veilleux et al. 2005, and references therein) applies to the late stage of galactic merger, (3) map the velocity structures of different phases of the outflow, and (4) refine estimates of the mass and energetics of the neutral phase, in which most of the mass likely resides.

F10565+2448 is one of the nearest and brightest starburst ULIRGs selected from the survey of Rupke et al. (2005a), with z = 0.0431 (Downes & Solomon 1998) and L_{IR(8–1000 μm)} = 1.1 × 10¹² L_☉ (Sanders et al. 2003). The CO redshift measured by Downes & Solomon (1998) is consistent with the optical redshift measured by Rupke et al. (2005a). It is classified as an H ii galaxy by its spectral type (Veilleux et al. 1995). However, in a more recent study, Yuan et al. (2010) classified it as a starburst-AGN composite galaxy, and the AGN contribution to infrared luminosity is estimated to be 10%–20% (Veilleux et al. 2009b; Yuan et al. 2010). The outflow's absorption lines' structure in this galaxy shows spatial variation in long-slit spectra (Martin 2006). A Hubble Space Telescope (HST) image shows that our target has three nuclei, and the one we observed is by far the most luminous. Figure 1 shows the HST image and a green box indicating the nucleus we observed.

Zoom In Zoom Out Reset image size

We observed this object with the Gemini Multi-Object Spectrograph (GMOS) IFU in two-slit mode using the B600 disperser centered on 6300 Å, which give us a wavelength range of 5600–6980 Å with spectral resolution of 4.26 unbinned pixels (roughly 1.75 Å). Each image covers a science field of 50 × 70, which at this redshift is about 4.2 × 5.9 kpc, with 1000 lenslets. There were two dithering positions with offsets of 15 from the center in the east–west direction. With eight exposures of 1800 s each, a total of four hours worth of data was collected. Seven frames were taken on the night of 2007 February 12 (UT) and one was taken on 2007 July 4 (UT). The seeing was variable within the range of 10–13.

The data were reduced using the Gemini IRAF package (version 1.8) following the standard steps: overscan subtraction, bias subtraction, flat-fielding, wavelength calibration, sky subtraction, flux calibration, and spectra extraction. The reduced two-dimensional images were assembled into a three-dimensional data cube with coordinates (x,  y,  λ) using the GFCUBE routine. For each data cube, differential atmospheric refraction was corrected by shifting the image slices at each wavelength to keep the centroid of the nucleus constant. After centering and combining all the cubes, and then trimming off the uneven edges, the total science field that was used for analysis is 78 × 66 or 6.5 × 5.5 kpc. The original resolution of the images was 02 pixel⁻¹. To get better signal-to-noise ratio (S/N) we binned the data cube to 06 pixel⁻¹, and since this is smaller than the seeing we did not lose significant spatial information. The final combined data cube was sliced into 143 one-dimensional spectra that can be used by our absorption-line fitting routine.

In our final data cube, there appears to be a resolution offset between the two sides of the field of view (FOV). After working with the Gemini Helpdesk, we learned that there is a known resolution difference between the two pseudo-slits that is probably flexure dependent. Since the reason for the offset is still unknown, there is no standard procedure to fix the problem. After inspecting the reduced data, we determined that we can make a meaningful analysis without any correction. This resolution difference in the wavelength direction can be up to 10 km s⁻¹ which could possibly affect our spectral fitting, however, we have checked the fitted parameters and found no clear systematic effects between the two pseudo slits.

Methods for analyzing superwinds using spectroscopy utilize both emission and absorption lines. However, emission lines can be contaminated by light from star-forming regions and they have an observational degeneracy between infall and outflows. Absorption lines have two advantages in the case of ULIRGs. First, emission lines in ULIRGs are dominated by rest-frame starburst component (e.g., Veilleux et al. 1995; Rupke et al. 2005b), while the absorption-line rest-frame component is typically swamped by the outflow component (Rupke et al. 2005b). This is at least partly due to the different density dependence of emission and absorption-line strengths. Infall and outflow sources can also be easily distinguished for absorption lines, since an absorbing cloud has to be between the observer and the light source. In emission lines, on the other hand, blueshifted gas could be inflowing or outflowing, depending on the (unknown) geometry. Rupke et al. (2005a) developed software for detailed analysis of the outflows via fitting the Na i D λλ5890, 5896 absorption lines. Because of its low ionization energy, Na i D probes the neutral gas, which carries significant amounts of mass and energy and is easily detectable in these galaxies. We used this software to analyze the properties of the wind in F10565+2448.

We use the fitting routine, NAFIT, developed by Rupke et al. (2005a) to analyze the Na i D feature. This routine fits multiple velocity components to the Na i D absorption lines by varying the covering factor, optical depth, line width (or more specifically the Doppler parameter b = $\sqrt{2} \sigma$ = FWHM/$[2 \sqrt{\ln \,2}]$), and central wavelength of each component and minimizing χ² using a Levenberg–Marquardt optimizer.

For each line, assuming a continuum level of unity, the intensity of the absorption line is given by

$$\begin{equation} I(\lambda) = 1 - C_{f}(\lambda) + C_{f}(\lambda)e^{-\tau (\lambda)}, \end{equation} \tag{ 1 }$$

where C_f is the covering factor and τ is the optical depth. Assuming a Maxwellian velocity distribution, τ can be expressed as a Gaussian:

$$\begin{equation} \tau (\lambda) = \tau _{c}e^{-(\lambda - \lambda _{c})/(\lambda _{c}b/c)^{2}}, \end{equation} \tag{ 2 }$$

where τ_c and λ_c are the central optical depth and central wavelength in the line and b is the width of the line. As the Na i D feature is a doublet, and the ratio of the optical depth of the two lines is determined by atomic physics, the central optical depth of the stronger line is always twice that of the weaker line, i.e., τ_2,c = 2 × τ_1,c. At any particular wavelength, τ(λ) = τ₁(λ) + τ₂(λ) where τ_i(λ) is given by Equation (2) for each line. We also assume an identical covering factor for both lines in a doublet. With this assumption, both optical depth and covering factor can be constrained simultaneously. For a detailed discussion on accounting for blended lines and overlapping absorbers, see Rupke et al. (2005a).

We wrote a wrapper to this fitting software to fit all of our spectra automatically. Starting at the center of the galaxy, where signal to noise is the highest, we step outward in a spiral pattern, using the output parameters of the previous fit as the initial guess for the next spectrum. Figure 2 is a signal-to-noise map of the continuum near the Na i D line, in units of S/N per unit spectral resolution element. The S/N is as high as 50 at the center of the nucleus and declines to as low as 2 or 3 at the edges of the FOV. As shown in Figure 3, at the continuum peak, the absorption line is best fitted by two velocity components. (Note: each velocity component contains a pair of doublet lines). The one-component fit can reproduce the basic shape of the line but leaves significant residuals. A two-component fit lowered the residuals significantly, and three or more components overconstrains the fit. As we move outward to regions with lower signal to noise, we find spectra that can be fit by one component equally well as by two. However, switching between one- and two-component fits within the same data cube introduces artificial patterns to our velocity maps. Therefore, we first perform two-component fits for all spectra; we then visually inspect the fits and remove any components that have negligible contribution to the equivalent width of the line. Examples of significant versus insignificant components are shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some spectra show the He i λ5876 emission line which appears immediately blueward of the Na i absorption lines. We include this emission line in our fit, where necessary, to remove its effect on our absorption-line fitting. For most of the spectra, however, the He i emission is weak enough that it blends in with the noise and can be safely ignored. Figure 5 shows examples of He i line fitting.

Zoom In Zoom Out Reset image size

All parameters are allowed to vary freely during the fit. Upper and lower limits of all parameters are the same for all spectra except for the central wavelength which is ±5 Å of the previous fit. The upper limit of the optical depth τ is set to 5 since our fits are insensitive to optical depths beyond that. The covering factor is constrained to be 0–1. To ensure that we are not extracting velocity dispersion widths smaller than the spectral resolution, the lower limit of the Doppler parameter is set to the spectroscopic resolution of 75 km s⁻¹ determined from the line width of the calibration lamp lines. The fitted velocity widths are typically well above this lower limit. It should also be noted that Rupke et al. (2005a) have shown that not de-convolving the instrumental resolution before the fitting does not change the result significantly. The upper limit of the velocity width is set to 500 km s⁻¹ and no spectra came close to exceeding this limit.

We have checked our Monte Carlo simulation outputs for possibilities of fitting degeneracies. Since this step essentially fits slightly modified versions of the same spectrum 1000 times, degeneracies for any parameters should show up as a bimodal or more complicated distribution in the fitting results. We find that almost all of the fitted parameters show a roughly Gaussian distribution. Exceptions occur mostly for the lower S/N spectra around the edges that are later discarded anyway due to their low quality.

We fit five emission lines ([O i] 6300 Å and 6364 Å, Hα, and [N ii] 6548 Å and 6583 Å) using a suite of routines that simultaneously estimates the stellar continuum and fits all emission lines at once, with the widths and redshifts constrained to be the same for each emission line. The stellar continuum, minus emission lines, was estimated using the stellar synthesis models of González-Delgado et al. (2005). The emission lines were fit with two Gaussians, one narrow and one broad; the broad component was discarded near the edges of the data, where the signal to noise was too low. The emission-line fitting routines are described further in S. N. Rupke et al. (2010, in preparation).

Figure 1 is the HST/ACS two-color image of our target, with F814W shown in red and F435W shown in green. The green box on the lower right nucleus indicates the GMOS FOV, and the continuum image of our data is shown in the lower left. The continuum flux map of our target is integrated over 6400–6500 Å. As mentioned in Section 2, the pixels are binned 3 × 3 to get better signal to noise.

The Hα flux map is shown in the top two panels of Figure 6. We have applied a global Balmer-decrement-determined extinction correction to the fluxes, E(B − V) = 1.75 (Veilleux et al. 1995), which corresponds to A_V = 1.75 × 3.1 = 5.43. For both components, the Hα flux peaks at the same pixel as the continuum flux map. The bottom two panels show the [N ii]/Hα flux ratio map for both components. For the narrow component, the line ratio is fairly uniform and symmetric around the center and has an average value of 0.5, consistent with star formation (Kewley et al. 2006). Away from the nucleus, however, the ratio increases. The broad component has a much wider range of line ratios and a peak of 1.8 appears in the south and northwest.

Zoom In Zoom Out Reset image size

Using the Hα flux from the emission-line fitting, we can estimate a star formation rate (SFR). Integrating over the Hα flux map, we obtained a total flux of F(Hα) = 4.08 ×10⁻¹² erg s⁻¹ cm⁻². Using Ned Wright's Javascript Cosmology Calculator (Wright 2006), assuming H_o = 71 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_vac = 0.7, we get a luminosity distance of D_L = 188 Mpc and therefore L(Hα) = 1.7 ×10⁴³ erg s⁻¹. The SFR can be estimated using the following equations (Kennicutt 1998):

$$\begin{equation} {\rm SFR} (M_{\odot }\, \rm yr^{-1}) = 7.9 \times 10^{-42} L(\rm H\alpha) (\rm erg\, s^{-1}) \end{equation} \tag{ 3 }$$

$$\begin{equation} {\rm SFR} (M_{\odot }\, \rm yr^{-1}) = 4.5 \times 10^{-44} L_{{\rm IR}(8\hbox{--}1000\,\mu \mathrm{m})} (\rm erg\, s^{-1}). \end{equation} \tag{ 4 }$$

The SFR estimated from L(Hα) is about 136 M_☉ yr⁻¹, and 190 M_☉ yr⁻¹ from L_IR. The SFR estimated from Hα flux is slightly smaller as expected since the extinction from Veilleux et al. (1995) is a lower limit.

Figure 7 shows the velocity map of the Na absorption lines and emission lines, respectively. (Note: since the signal to noise is very low at the edges of the FOV, we did an additional 3 × 3 binning at the four corners so we can constrain the parameter better. This only applies to the absorption-line maps and not the emission-line maps. This additional binning shows up as four large pixels, which are nine times larger than the rest of the pixels, at the four corners of all the absorption-line velocity and velocity dispersion maps.) Velocities are calculated with respect to a fixed redshfit of z = 0.0431. The narrow emission-line map, which we interpret as systemic, is blueshifted in the east and redshifted in the west and clearly shows rotation (see Section 5 for more discussion). The position angle (P.A.) of the major axis according to Downes & Solomon (1998), 100° east of north, is plotted as red lines on the absorption lines and narrow emission-line velocity maps as a reference. Most of the northern half of the broad line emission map has values very close to systemic, but the southern half shows significant blue shifts. The absorption lines are separated into higher-outflow-velocity (blue) and lower-velocity (red) components. Both components are generally blueshifted with respect to systemic and their velocity gradients are east–west oriented, roughly aligned with the major axis.

Zoom In Zoom Out Reset image size

Figure 8 shows the velocity dispersion, or full width at half-maximum (FWHM), maps of the emission and absorption lines, respectively. Both narrow and broad emission lines have a high velocity dispersion region in the southeast, roughly in the same region as the peak of [N ii]/Hα. The red absorption component has a typical FWHM of roughly 200 km s⁻¹, while the bluer component has a higher typical FWHM at about 300 km s⁻¹. The systemic emission-line FWHM increases away from the nucleus from about 150 to 350 km s⁻¹.

Zoom In Zoom Out Reset image size

NAFIT calculates the error in each fitted parameter using Monte Carlo simulations described in detail in Rupke et al. (2005a). Errors in the velocities are 5–10 km s⁻¹ close to the center and remain ⩽ 30 km s⁻¹ for most of the pixels, although for some pixels close to the edge the error can be > 50 km s⁻¹. FWHM errors start from ⩽10 km s⁻¹ in the center, remain below 50 km s⁻¹ for most pixels, and increase up to >80 km s⁻¹ around the edges.

Other than the random errors due to noise, the systematic effect between the two pseudo slits mentioned in Section 2 may also cause errors. The two pseudo slits are oriented north–south with respect to each other. Judging from the images of the fitted parameters, we can see that there is no clear systematic offset in the north–south direction (the direction of the systematic offset). The velocity gradient is oriented in the east–west direction, and from the arc lamp measurements we find that the resolution offset is small (about 10 km s⁻¹) compared with the velocity gradient observed. Therefore, we can safely conclude that any errors caused by the systematic offset have a negligible effect on our analysis.

Our target is known to be dusty and we can test whether the neutral gas absorption corresponds to the dust features. A map of the Na i D line equivalent width is shown in the right panel of Figure 9. The equivalent width shown here is the total equivalent width of both velocity components. A map of the dust features is shown in the Sloan Digital Sky Survey (SDSS) g/i ratio map in left panel of Figure 9. We are using g/i ratio as a proxy for E(B − V) and assuming that the stellar population does not vary too much across the FOV. A lower g/i ratio indicates a redder color and corresponds to the easily visible dust filaments in the high-resolution HST image (Figure 1). The SDSS g and i images are rotated and resampled to match the P.A. and resolution of our data. A pixel-to-pixel comparison plot of g/i ratio and equivalent width images is shown in Figure 10. The Na i D line strength correlates well with the dust reddening.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compared our IFU data with the freshly analyzed long-slit observations of our target along a P.A. of 53°(Rupke et al. 2005a). The long-slit data probe further from the nucleus and show the presence of a neutral (ionized) wind up to 7–8 (10) kpc from the nucleus, as well as regions of large line broadening and shock excitation that extend 10 kpc from the nucleus. This is consistent with similar results from Martin (2006) along a P.A. of 109°. Quantities calculated from Rupke et al. (2005a) long-slit data, including outflow velocities, velocity dispersions, Hα fluxes, and [N ii]/Hα ratios, agree with our measurements from the IFU data.

5.1. The Outflow Velocity Structure

5.2. An Extended Multi-phase Outflow

5.3. Outflow Parameters

We can refine the estimation of mass and kinetic energy in the neutral wind using a wind model that depends on the physical parameter outputs of our fitting. We choose a spherical thin shell instantaneous outflow similar to the one used by Rupke et al. (2005b). The difference is that Rupke et al. assumed spherically symmetric velocity and density, but we have information on the spatial variations of these quantities. An alternative to the thin shell model is an outflow with finite thickness and velocity that depends on radius. In Veilleux et al. (1994) and Martin & Bouché (2009), the thick shell outflow has V ∝ Rⁿ, where R is the (spherical) radial distance from the nucleus and V is the outflow velocity. Because we observe the galaxy in a face-on orientation, we cannot easily probe variations in velocity with R, so we assume a single, characteristic value for R but let V vary across the outflowing shell according to the velocity output of NAFIT.

We have de-projected the velocities from observed to purely radial velocities. The instantaneous mass outflow rate, average outflow rate, and total mass are

$$\begin{equation} dM/dt_{\rm thin}^{\rm inst} = \mu m_{p} R^{2} \sum \frac{N(\theta,\phi)}{dr} v(\theta,\phi) {\rm sin}(\theta) d\theta d\phi \end{equation} \tag{ 5 }$$

$$\begin{equation} dM/dt_{\rm thin}^{\rm avg} = \mu m_{p} R \sum N(\theta,\phi) v(\theta,\phi) {\rm sin}(\theta) d\theta d\phi \end{equation} \tag{ 6 }$$

$$\begin{equation} M_{\rm thin} = \mu m_{p} R^{2} \sum N(\theta,\phi) {\rm sin}(\theta) d\theta d\phi, \end{equation} \tag{ 7 }$$

where R is the radius of the thin shell, theta and phi are the usual spherical angular coordinates, dr is the thickness of the shell (dr ≪ R), μ is the mean molecular weight, m_p is the mass of a proton, and N is the column density of hydrogen. For all our calculations, we choose R = 5 kpc to be consistent with Rupke et al. (2005b) and dr = 0.5kpc ≪ R. Using the mass outflow rate we can calculate the energy and momenta outflow rates (see Rupke et al. 2005b), the values are summarized in Table 1. The column density of hydrogen, N(H), can be estimated from N(Na i) given a metallicity and ionization fraction:

$$\begin{equation} N(\rm H) = N(\rm Na\,\mathsc{i})(1 - y)^{-1}10^{-(a+b)}, \end{equation} \tag{ 8 }$$

where $y = (1 - \rm Na\,\mathsc{i}/\rm Na) = 0.9$, $a = \log [N(\rm Na)/N(\rm H)] = -5.69$ which we set to solar value, and $b = \log [N(\rm Na)/N(\rm H)_{\rm total}] - \log [N(\rm Na)/N(\rm H)_{\rm gas}] = -0.95$ is the assumed depletion.

Table 1. Outflow Characteristics

Quantity	Total	Outflow Rate (Instantaneous)	Outflow Rate (Average)
Mass	9.5	3.1	2.1
Energy	57.1	43.3	42.3
Momentum	49.8	36.0	35.0

Notes. Mass, energy, momentum, and the outflow rates in logarithmic values. Column 1: the log of total quantity of mass, energy, and momentum; Column 2: the log of instantaneous outflow rate; Column 3: the log of average outflow rate. The units are the following: mass—M_☉, mass outflow rates—M_☉ yr⁻¹, energy—erg, energy outflow rates—erg s⁻¹, momentum—dyne s, and momentum outflow rates—dyne.

Download table as:  ASCIITypeset image

The values that we calculated for mass, energy, momentum, and their outflow rates are within a factor of 3–4 of the values in Rupke et al. (2005b). The new mass, energy, and momentum are smaller than what is predicted from long-slit data, while the outflow rates are larger. With the new spatial information on velocity and column densities, the new values should be more accurate. The thin shell geometry that we assumed is, however, uncertain. Theory and simulation predicts a thin shell of warm gas behind the forward shock of the wind (e.g., Castor et al. 1975; Suchkov et al. 1994; Strickland & Stevens 2000), but it is unclear whether it is accompanied by neutral gas. Since we have limited information on the radial depth of the wind, we cannot rule out a thick outflowing shell in which neutral clouds are continuously entrained in the wind. Switching from a thin shell to thick shell model will decrease our measured mass. Furthermore, we are only considering the part of the wind where we have reliable spectral fitting; the actual wind probably extends further than what we observed and carries more mass.

In some optically thick cases, our fitting routine finds only a lower limit for τ for one of the velocity components, which is set to 5 in our fits. This is applicable to ∼9% of the spectra used for this section's calculations and may cause us to underestimate the mass. We check the extent of this effect by comparing our values to those calculated by using only one velocity component fit throughout. When only one component is used for fitting, only <2% of the fitted τ used for calculation are at the lower limit. The estimated mass and energy are slightly smaller than but within a factor of two of the values presented in Table 1. From the present data alone, we cannot rule out the possibility that optically thick lines cause a significant underestimate of the mass in the wind. However, we believe it is more likely that our mass estimates are not biased in this way because: (1) the saturated components in the current data arise primarily at low signal to noise; (2) the one-component fits, which have more sensitivity to the overall shape of the lines in the low S/N region, did not give a drastically different mass estimate. Higher S/N data will likely resolve this issue and better constrain the optical depths of each component.

Since our observation only trace the Na i D clouds, it is possible that we underestimated the total mass of the wind by not considering the other phases. Rupke et al. (2005b) compared ratio of cold, neutral gas outflow rates to the expected mass injection rated from the hot gas created by supernovae. The ratios are found to be close to unity on average but it can vary from 0.001 to 10 for individual galaxies. On the other hand, the theory suggests that the entrained clouds dominate the mass budget of the wind (e.g., Strickland & Stevens 2000), and we suspect that the Na i D clouds are the entrained gas (Rupke et al. 2005b). The conversion factor used to "convert" the neutral Na column density to that of the total Na has also been shown to be uncertain (Murray et al. 2007). However, the total mass in the cold gas cannot be more than an order of magnitude larger than our calculated value, otherwise the mass ratio between cold and hot gas will be too large (Strickland & Stevens 2000; Strickland & Heckman 2009).

Within the above constraints, the present data put on firmer footing the coarse mass outflow rates measured by Rupke et al. (2005b). It is clear that these outflows are massive and energetic. The neutral outflow rate in this particular system is comparable to the SFR. However, this particular system has a higher outflow rate than the average value in ULIRGs (Rupke et al. 2005b). Further spatially resolved and multi-phase observations are necessary to understand the exact relationship between mass outflow rate and SFR in galactic superwinds.

We have presented Gemini integral-field data of a galactic wind in the starburst ULIRG F10565+2448 to study the structure of the wind in detail. We fitted the Na i D absorption lines and the Hα and [N ii] emission lines across the 7 × 6 kpc FOV using multiple velocity components. The systemic rotation velocity determined from the narrow emission lines is consistent with the molecular gas disk rotation. The absorption-line gas velocities shows asymmetric outflows with velocity gradients along the major axis opposite to that of the systemic rotation. We see widespread evidence of shock ionization through collisional excitation of the gas and large velocity dispersions, with particularly strong shocks observed in the southeast corresponding to the region with lower absorption-line outflow velocities. This supports the notion of gas deceleration in the east, which should be accompanied by shocking. Our results strengthen the hypothesis that ULIRG winds may differ in structure from those in the more quiescent galaxies.

We observed both neutral and ionized outflowing gas. The correlation between the outflow velocities of these two phases is weak, suggesting a complex relationship between the two. We found a spatial correlation between the dusty region and the high Na i D equivalent width region, confirming spatial coincidence of dust and neutral entrained outflow material. The presence of dust in the outflow raises the possibility that winds may be responsible for removing dust around obscured AGN and revealing a bright quasar.

We calculated the mass, energy, and momentum outflow rate assuming a spherical thin-shell outflow. The numbers that we obtained are in general agreement with those in Rupke et al. (2005a). However, by studying the structure of winds in other starburst mergers, we should be able to improve on the current results. First, by studying edge-on disks we will be able to better constrain the opening angles of these winds. Second, by studying a sample of winds, we will place better constraints on the relationship between the ionized and neutral phases of the wind and whether asymmetric velocity structures are common. Finally, and perhaps most importantly, we must compare observations to detailed simulations of mergers that include realistic wind prescriptions—such simulations do not yet exist but are necessary for understanding detailed data sets like the present one.

The authors thank Jabran Zahid for sharing routines that were the framework for our stellar and emission-line fitting software. This paper used data from the NASA/ESA Hubble Space Telescope and obtained these data from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA), and the Canadian Astronomy Data Centre (CADC/NRC/CSA). This paper used data from the Sloan Digital Sky Survey (SDSS). Funding for the Sloan Digital Sky Survey (SDSS) has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. D.S.N.R. is supported by a University of Hawaii start-up grant to Lisa Kewley and by NSF CAREER grant AST07-48559. This work is based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership.