© ESO, 2016

1. Introduction

Complex feedback phenomena, such as galactic winds (GWs), play a major role in the evolution of galaxies and the surrounding intergalactic medium (IGM) in which they are embedded. In particular, GWs self-regulate the growth of both the stellar and black hole masses in galaxies (e.g., Veilleux et al. 2005 for a review) and they are also invoked for the color transformation of spiral galaxies into passive ellipticals caused by gas exhaustion (e.g., Weiner 2009). The GW feedback has been discussed by a number of authors (e.g., Veilleux et al. 2005; Rupke et al. 2005c) and is now part of many galaxy formation models (e.g., Dutton & van den Bosch 2009; Fabian 2012; Hopkins 2015).

Starburst-driven GWs are generated by the radiation and mechanical energy liberated by massive stars and the sequential explosions of SNe (Chevalier & Clegg 1985; Hopkins et al. 2012); and they are common in a range of galaxy types: dwarfs (e.g., Schwartz & Martin 2004), starburst (e.g., Chen et al. 2010) and luminous infrared galaxies (LIRGs; L_IR = L_{(8−1000 μm)} = 10¹¹−10¹² L_⊙) and ultraluminous infrared galaxies (ULIRGs, L_IR ≥ 10¹²; e.g., Martin 2005; Rupke et al. 2005b,c; Arribas et al. 2014). Active galactic nuclei (AGN) are invoked as the power source for the most extreme (i.e., fastest and powerful) outflows observed in quasars (e.g., Cicone et al. 2014; Sturm et al. 2011; Villar Martín et al. 2014).

Over the past 20 years, the observational multiwavelength data on outflows is steadily increasing both in quality and quantity. Outflows are a multiphase phenomena and, therefore, their different phases (hot, warm, neutral, molecular) can be studied via different tracers (e.g., FeXXV at 6.7 keV, Hα, NaD, CO(1−0) at 115.271 GHz) in different wavelength ranges (X-rays, optical, and radio; Veilleux et al. 2005).

GWs are observed at any redshift. During the peak of star formation (z ~ 2−3) warm-ionized outflows are very prominent and ubiquitous (Maiolino et al. 2012; Cano-Díaz et al. 2012; Newman et al. 2012; Förster Schreiber et al. 2014).

In the local Universe, ionized (e.g., Arribas et al. 2014), neutral (e.g., Rupke et al. 2005a,b,c), and molecular (e.g., Cicone et al. 2014) GWs have been found in galaxies undergoing intense and spatially concentrated star formation (i.e., starburst galaxies) in which the star formation rate (SFR) per unit area exceeds Σ_SFR ≥ 10^-1M_⊙ yr^-1 kpc^-2 (Heckman 2002).

Of particular interest is to investigate the amount of gas entrained in winds since its evacuation may have a significant impact on the gas reservoir available for star formation. A method to detect the neutral gas in GWs consists in tracing NaD absorption against regions with a bright stellar continuum.

In this context, [U]LIRGs, which are mostly powered by star formation, are ideal objects to investigate the neutral phase of SNe-driven outflows although they have a more important AGN contribution with increasing luminosity (Veilleux et al. 2009). In particular, [U]LIRGs offer the opportunity to study the feedback phenomenon in environments similar to those observed at high-z, but with a much higher signal-to-noise (S/N) and spatial resolution. Indeed, local [U]LIRGs share some basic structural (Arribas et al. 2012) and kinematical (Bellocchi et al. 2013) properties with distant star-forming galaxies, although they may differ in terms of gas content (which increases with increasing redshift; Tacconi et al. 2010).

Most of the previous observational works on outflows in local starbursts and [U]LIRGs have been based on long-slit data, which often provide a limited knowledge of their geometry and kinematics. Integral field spectroscopy (IFS) provides a better characterization of the morphology and the velocity structure of GWs (e.g., Colina et al. 1999; Jiménez-Vicente et al. 2007; Westmoquette et al. 2011; Bellocchi et al. 2016).

In this paper, we use IFS-data taken with the VIMOS at the Very Large Telescope (VLT) to significantly expand in number previous samples of spatially resolved neutral outflows (e.g., Rupke & Veilleux 2013, 2015), especially in the less studied LIRG luminosity range. In particular, we analyze the NaD absorption doublet in 30 LIRGs and 8 ULIRGs (51 individual galaxies). Our main goal is to characterize the kinematics and structure of the neutral outflows and investigate how they are related to the host properties, for instance, the infrared luminosity (i.e., star formation), and dynamical mass.

The paper is organized as follows. In Sect. 2 we briefly describe the sample, observations, data reduction, and line fitting. In Sect. 3 we present the data analysis, including details about the continuum modeling, and the generation of integrated spectra and spectral maps. In Sect. 4, the properties of the neutral gas kinematics are commented and discussed with special focus on the characteristics of GWs and their relation with other galaxy properties. Finally, the main results and conclusions are summarized in Sect. 5. In Appendix A, the kinematic maps, comments on the individual objects are presented. Throughout the paper, we assume H₀ = 70 km s^-1 Mpc^-1 and standard Ω_m = 0.3, Ω_Λ = 0.7 cosmology.

2. Sample, observations, data reduction, and line fitting

2.1. Sample and observations

The sample contains 51 individual galaxies mainly drawn from the IRAS Revised Bright Galaxy Sample (RBGS; Sanders et al. 2003), which have been observed via IFS by Arribas et al. (2008) (see Table 1). This sample contains a large number (i.e., 43/51) of sources in the less-studied LIRG luminosity range (L_IR = 2.9 × 10¹¹L_⊙, on average), while a smaller number of objects (i.e., 8/51) are classified as ULIRGs with L_IR = 1.6 × 10¹²L_⊙ on average. The mean redshift of the LIRGs and ULIRGs subsamples are 0.024 and 0.069, respectively (see Table 1).

The sample is not complete either in luminosity or in distance (Rodríguez-Zaurín et al. 2011; Bellocchi et al. 2013), however, it is representative of the [U]LIRGs population. Indeed, this sample covers a variety of morphological types along the merging process (i.e., rotating disks, interacting systems, and mergers) and encompasses objects with different nuclear ionization types (i.e., HII, Seyfert and LINER). See Table 1, for details.

Nearly half of the objects in the sample (i.e., 22/51) show an excess of 24 μm, hard X-ray emissions, or they have optical line ratios indicating the presence of AGNs (Table 1). However, the AGN contribution to the infrared luminosity or to the ionized gas kinematics is substantial in the cases of the ULIRGs IRAS F05189-2524, F23128-5919 (S) and for the LIRG F07027-6011 (N) (Table 1).

Throughout the paper, we made use of the SFR derived by Rodríguez-Zaurín et al. (2011) inferred from L_IR, following Kennicutt (1998; with a Salpeter initial mass function); we also made use of the dynamical masses (M_dyn) and angle of inclination of the galaxies derived by Bellocchi et al. (2013).

The observations were carried out with the integral field unit (IFU) of VIMOS at the 8.2-m telescope of the ESO-VLT (Le Fèvre et al. 2003). In brief, the IFS data of all the targets were acquired using the high resolution orange grating, which provide an intermediate spectral resolution of 3470 (dispersion of 0.62 Å pix^-1) in a wavelength region from 5250 to 7400 Å. The field of view (FoV) in this configuration is 27′′ × 27′′ with a spaxel scale of 0.67′′ per fiber. The 1600 spectra obtained are organized in a 40 × 40 fiber array, which constitute one single pointing. A square of 4 pointing dithered pattern for each target was used, providing an effective FoV of 29.5′′ × 29.5′′. For more details about the observations, see Arribas et al. (2008).

We excluded from analysis the galaxies 08355-4944 and F21130-4443 from the [U]LIRGs sample described above since they show very weak NaD absorption preventing any robust study. In addition, nine objects (i.e., F06035-7102; F06076-2139 (N, S); F06259-4780 (S); F06295-1735; F08520-6850; F12596-1529; F17138-1017; F22491-1808) have a large percentage of low S/N spectra in individual spaxels and were not suitable for the spatially resolved analysis. In summary, of the total sample of 51 individual objects observed, the integrated and spatially resolved properties have been studied for 49 and 40 galaxies, respectively.

2.2. Data reduction

We reduced VIMOS raw data with the pipeline provided by ESO via ESOREX (version 3.6.1 and 3.6.5), which allows us to perform the following steps: sky and bias subtraction, flat-field correction, spectra tracing and extraction, correction of fiber and pixel transmission, flux and wavelength calibration flux. We also used a set of customized IDL and IRAF scripts. Once the four quadrants were individually reduced for the four dithered positions, they were combined into a single data-cube made of 44 × 44 spaxels (i.e., 1936 spectra). We checked the width of the instrumental profile and wavelength calibration using the [O I] λ6300.3 Å sky line. The average values for the full width half maximum (FWHM) and central wavelength for the whole sample were (6300.29 ± 0.07) Å and (1.80 ± 0.07) Å, respectively. For each spectrum, we corrected for the effect of instrumental dispersion by subtracting it in quadrature from the observed line dispersion, i.e., σ_line= . A more detailed description of the data reduction process is given in Monreal-Ibero et al. (2010).

2.3. Line fitting

The observed NaD-absorption features are modeled both in the IFS and integrated spectra with two Gaussian profiles (i.e., single kinematic component) as in Cazzoli et al. (2014; also Davis et al. 2012). This approach has the limitation that it is only correct in the cases of either uniform covering factor and low optical depth or in the case of covering factor that varies with velocity as a Gaussian (independently to the optical depth). Absorption line fitting methods that consider both optical depth and covering factor are presented in Rupke et al. (2005a), and references therein. Although the observed NaD-absorption features show generally well-resolved lines and not saturated profiles (i.e., flat-bottomed), we stress that the validity of our fitting approach is limited to the cases mentioned above. As in many cases the line ratio within the two doublet lines indicates an optically thick regime, it is implicity assumed that the covering factor varies with velocity as a Gaussian in these cases.

We found some spectra that have strongly asymmetric profiles and their modeling requires two Gaussian pairs (i.e., a single component fitting is certainly an oversimplification). After visual inspecting all the spectra, we decided to apply a two-components line fitting for the spectra of four objects (F07027-6011 (S), F10409-4556, F11506-3851, and F12116-5615) similarly as in Cazzoli et al. (2014). In a few other cases, there is evidence for the presence of multiple kinematic components in individual spectra (e.g., F04315-0840, F10257-4310, F13229-2934). However, the limited S/N does not allow us a robust kinematic decomposition into two different components. Therefore, a single kinematic component fitting of the NaD absorption doublet was preferred in these cases. Figure 1 shows examples of the two different Gaussian fits for the galaxy IRAS F10409-4556. The existence of two kinematical components is expected as the NaD absorption feature may take part in the galaxy ordered rotation and be entrained in winds with blueshifted velocities. Therefore, we identify these components according mainly to their spatial distribution and velocities.

We also fit the NaD line profile in the integrated spectra (after the subtraction of the stellar contribution, Sect. 3.1) with one or two components following the approach described above. The majority (21/28) of these integrated spectra show a prominent nebular emission line HeI λ5876 Å, which was included in the line fitting (i.e., modeled selectively with one or two Gaussian components and then subtracted). This strategy has also been applied in a few cases for the spatially resolved data.

3. Data analysis

3.1. Integrated spectra: generation and modeling of the continuum and line profile

From the lFS data cubes, we generated a spatially integrated spectrum per galaxy via the S/N optimization method proposed by Rosales-Ortega et al. (2012), using the pingsoft tool (Rosales-Ortega 2011). Before integration, the spectra were corrected from the general stellar velocity pattern. For this step, we used the narrow component of the ionized gas velocity field, which describes the systemic behavior (Bellocchi et al. 2013). This method significantly improves the modeling of weak features of the stellar continuum compared to other techniques (e.g., S/N cutoff; see Rosales-Ortega et al. 2012 for a detailed discussion).

We applied a penalized PiXel fitting analysis (pPXF; Cappellari & Emsellem 2004) for the recovery of the shape of the stellar continuum. As in Cazzoli et al. (2014), we used the Indo-U.S. stellar library (Valdes et al. 2004) to produce a model of the stellar spectra that matches the observed line-free continuum (any interstellar medium (ISM) features, including NaD, are ignored in the fitting). The result of this approach is a model that in general reproduces the continuum shape well (the residuals are typically <10%) except few cases. The integrated and model spectra are shown for each galaxy in Appendix A.

We checked all the stellar continuum modeling by visual inspection of the residuals especially in the line-free continuum regions and nearby the NaD feature (e.g., we checked if any potential NaD-HeI blending may affect the continuum model). We then flagged each object with the parameter Q, indicating the quality of stellar modeling (Table 2). From the total sample of 49 objects, 31 were flagged with good (Q = 2) or very good (Q = 3) modeling. For these objects, we generated a pure ISM-spectra by subtracting from the observed NaD profile the model of the stellar spectra obtained via the pPXF method. This model was also used to fix the stellar (systemic) zero velocity of the spectra. Then, we applied the Gaussian fitting described in Sect. 2.3 to the spectra to derive the neutral gas kinematics. Figure 2 shows the purely ISM NaD absorption and its modeling, while the kinematics properties are summarized in Table 2 for each galaxy. For only 3 out of 31 (i.e., F01159-4443 (N), F21453-3511, F23128-5919 (N)) was not possible to model the purely ISM NaD line profile since the stellar subtraction resulted in an almost undetectable neutral gas NaD doublet (the stellar contamination is >95%, Table 2).

In ~61% (i.e., 17/28) of the cases a single kinematic component (a Gaussian pair) already gives a good fit, suggesting that if a second component exists in these galaxies, it is weak. Only 11 out of 28 require two kinematic components; we find that a two-Gaussian component model per line led to a remarkably good fit of the NaD absorption, significantly reducing the residuals with respect to one-Gaussian fits (Fig. 2).

We tested the feasibility of the stellar continuum modeling on a spaxel-by-spaxel basis. Unfortunately, the spectra in individual spaxels, in general, lack the S/N required making the determination of the stellar and ISM contributions in 2D highly inaccurate. Therefore, we used a different (and more simplistic) approach to evaluate the stellar contamination in the NaD maps for the 2D analysis, as discussed in Sect. 4.1.

3.2. Two-dimensional maps

The 2D analysis is based on the spectral maps (i.e., velocity field (V), velocity dispersion (σ), and equivalent width (EW)) of the different kinematic components, which were generated after the line fitting procedure (Sect. 2.3). To produce the maps themselves, we used a set of IDL procedures (i.e., jmaplot) developed by Maíz-Apellániz (2004). The maps are shown, for each galaxy (for a total of 40 objects), in Appendix A. As a reference we also show the continuum image and Hα map of the systemic (narrow) component generated from the VIMOS-IFU data cube (Bellocchi et al. 2013). The integrated spectrum (Sect. 3.1) is also shown for each object . The maps of the LIRG F11506-3851 are already discussed in detail in Cazzoli et al. (2014) and, therefore, they are not included in Appendix A, but the results are considered in the overall statistics and figures.

4. Results and discussion

4.1. Contamination from stellar NaD absorption

In order to measure the stellar contribution to the NaD doublet (i.e., stellar contamination), we computed the ratio between the EW(NaD) in the stellar model obtained with pPXF (Sect. 3.1) and the total EW(NaD) observed in the integrated spectra. This ratio is listed as percentages for each individual galaxy in Table 2 and a histogram of its distribution is presented in Fig. 3.

We divided the sample into three groups: i) objects for which the NaD absorption is dominated by the ISM (i.e., stellar contribution <35%), ii), medium-contaminated objects (i.e., stellar contribution between 35% and 65%), and stellar dominated objects (i.e., stellar contribution >65%). Considering only the 31 objects with good stellar continuum modeling (Fig. 3), the stellar contribution is dominant for 23% of the cases, while the NaD absorption is dominated by the ISM for 42% of the cases.

As mentioned above, the stellar continuum fitting method could not be applied on a spaxel-by-spaxel basis. Therefore, to evaluate the origin of the NaD absorption in the maps, we compare the values of the sodium equivalent width, EW(NaD), with that of a purely stellar NaD doublet. In fact, we consider that the NaD feature is mainly originated in the ISM when EW(NaD) > 1.3 Å. We established this (conservative) threshold considering the EW(NaD) of the stellar models obtained with pPXF, which are shown in Fig. 4 as a function of the stellar NaD fraction. In fact, independent of the fraction of NaD originated in the stars, the EW of the stellar models is on average 0.9 Å. Therefore, our choice of EW(NaD) > 1.3 Å to identify a NaD feature dominated by the ISM is rather conservative as it is above the stellar contribution in 90% of the cases (or 100% within 1σ). Therefore, the generally strong NaD absorption (with EW(NaD) > 1.3 Å) seen in many regions of the [U]LIRGs (Appendix A) can be robustly associated with the cold neutral ISM. As for comparison, starburst galaxies with EW(NaD) > 0.8 Å are considered strong ISM-NaD absorbers by Chen et al. (2010).

4.2. Galactic winds in one dimension (integrated spectra)

When identifying GWs in the integrated ISM spectra (Sect. 4.1) only components with central line velocities blueshifted with respect to the systemic by more than 60 km s^-1 are considered to be outflowing (Col. 9 in Table 2). This cutoff, which is on the order of two spectral elements of ~0.6 Å each, is similar to those assumed in previous studies (e.g., 50 km s^-1 in Rupke et al. 2005b). Table 1 lists the classification based on this criterion for each individual object in the survey.

Neutral gaseous outflows are found in 19 out of 28¹ [U]LIRGs. The incidence of winds in our local [U]LIRGs, i.e., 68%, is similar to that found in other samples of infrared-luminous galaxies (e.g., Heckman et al. 2000; Rupke et al. 2005b). In five objects the NaD absorption is found at systemic velocity (i.e., −60 km s^-1<V< + 60 km s^-1), while in one case (F06259-4780 (C)) is redshifted.

We found that, generally, the outflows have velocities, V, in the range from 65 km s^-1 to 260 km s^-1 (on average ~165 km s^-1). In Fig. 5 (left) we show the neutral gas outflow velocities obtained from the spatially integrated analysis as a function of the SFR. A regression of the type V ∝ SFRⁿ to this mainly LIRGs sample excluding objects with evidence of strong AGNs yields n = 0.15 ± 0.06. This dependency is in rather good agreement with the results of Rupke et al. (2005c; n = 0.21 ± 0.04) for their sample of (mostly) ULIRGs. These results are also fairly consistent with those of Martin (2005; n = 0.35), who considered a sample of ULIRGs and three dwarf galaxies of low SFRs. As shown in Fig. 5, the different LIRG and ULIRG samples complement each other well, sampling the SFR rather homogeneously over nearly 2 orders of magnitude. The general trend defined with these 1D (i.e., integrated and long-slit) data is consistent with an index n in the range 0.1−0.2. The effects of inclination and the presence of AGNs make it difficult to constrain this value further with those samples.

The comparison of the kinematic properties of the neutral outflows with their ionized counterparts (Arribas et al. 2014) for the same objects, is shown in Fig. 6.

On the one hand, we find that the outflow velocities (measured at the center of the line) are significantly higher for neutral than for ionized outflows in all cases. On the other hand, the neutral gas velocities seem on average slightly smaller than the commonly used ionized gas maximum outflow velocities, which are defined as V_max = | −ΔV | + FWHM/2 (Westmoquette et al. 2011; Genzel et al. 2011).

Our measurements seem to contradict previous single aperture studies for which ionized and neutral winds are correlated in Seyferts but not in starburst (Rupke et al. 2005b). On a spatially resolved basis, only weak evidence of a correlation between the velocities of the neutral and ionized gas outflow phases (within a given galaxy) are found by Rupke & Veilleux (2013) for a sample of six nearby merger [U]LIRGs (including three obscured QSOs). However, all these findings are not fully comparable as they differ in terms of type of object and ionized gas tracers. For the velocity dispersions, both neutral and ionized outflows show similar median values of 115 and 138 km s^-1, respectively, if the strongest AGNs are excluded.

4.3. Galactic winds in two dimensions: Identification and detection rate

We carried out the search of GWs in 2D by inspecting the NaD spectral maps for regions with, simultaneously:

• 1.

  significant blueshifted velocities (≳60 km s^-1) with respect to the systemic, which cannot be explained by rotation. For this identification, we consider the velocity field of the narrow component of Hα as reference, for which Bellocchi et al. (2013) have shown follows systemic behavior (Appendix A);

• 2.

  EW(NaD) > 1.3 Å, to guarantee that the NaD feature is dominated by ISM absorption; and

• 3.

  a relatively broad kinematic component (σ> 90 km s^-1).

This procedure identifies 22 objects out of 40 (55%) with outflows (see Table 1), implying a detection rate that is slightly lower than that obtained from the analysis of the integrated spectra (i.e., 71%). This is likely a consequence of adopting a stricter criterion for the identification of outflows in the spectral maps. However, the agreement in the identification of a GW and/or rotation between these two procedures is high. From the 25 objects for which it was possible to make a kinematical classification via both, spatially resolved and spatially integrated spectra, 18 have outflow signatures and four rotation or systemic signatures, according to the two methods. Therefore, we find inconsistency only in three cases. In two cases (F12115-4656 and F12043-3140 (S), Figs. A.27 and A.26), this is because the putative outflows only cover a small faint region of the maps and, as a result, are undetected in the integrated spectra. In the other case (F07160-6215, Fig. A.14), the neutral gas velocity field has a complex structure and the wind partially overlaps the approaching side of the rotation pattern, making the 2D identification difficult.

Previous works have shown that the ability of detecting neutral GWs with absorption tracers depends on the galaxy inclination angle. Specifically, Heckman et al. (2000) found a probability of ~70% of detecting outflowing gas in absorption in starburst galaxies with an inclination less than 60°. Similar trends are also seen in more recent studies of the relation between EWs(NaD) and the galaxy inclination angle (e.g., Chen et al. 2010). As seen in Fig. 7, in our sample GWs are detected in galaxies observed in a wide range of inclination angles (i.e., 10°−80°). Although we also find a large percentage of GWs (~78%) detected in galaxies at low inclination angles (<60°), this percentage is similar to the total number of galaxies with low inclination whether a GW is detected or not.

In addition to the 22 objects with neutral gas GW detection, 10 targets show spider-like NaD velocity field in at least one kinematic component. This indicates disk rotation and this subsample is discussed in Sect. 4.8. Nine sources (eight LIRGs and one ULIRG) lack any clear outflows or rotation signatures either because of a lack of blueshifted velocities or because the putative disk kinematic is very irregular; these sources are excluded for the following discussion. The maps and integrated spectra of these nine sources, however, are shown in Appendix A. Table 1 (Col. 9) summarizes the different kinematic patterns found in the maps.

4.4. Spatially resolved neutral GWs properties: two-dimensional kinematic and geometry

We derived the (typical) wind velocity (V_w) as the deprojected median value² over the region identified as GW from the spectral maps. The inclination-corrected outflow velocities of the neutral gas entrained in GWs are in the range ~80−700 km s^-1, which is similar to the velocities found previously for starburst driven winds in local star-forming galaxies and [U]LIRGs (Heckman et al. 2000; Martin 2005; Rupke et al. 2005c; Chen et al. 2010). An exception to the general behavior is the outflow in the ULIRG F06206-6315 with velocities that exceed 1000 km s^-1, which host an AGN (Table 1).

The outflow velocity inferred from the maps also has a clear correlation with the global SFR (Fig. 5, right). In particular, on the one hand, V ∝ SFRⁿ with n = 0.40 ± 0.07 for all cases, and n = 0.30 ± 0.05 excluding the objects with an AGN. This is in rather good agreement with previous results, including those obtained from the integrated spectra (see Fig. 5, left). On the other hand, the wind velocity does not seem to correlate with the dynamical mass of its host (Fig. 8).

Results from a previous IFS survey of local major-merger [U]LIRGs (Rupke & Veilleux 2013) indicate that in such systems the typical neutral gas outflow velocity dispersion is high, up to 1000 km s^-1 (FWHM) in the case of MRK231. In our sample of (mainly) LIRGs, the typical velocity dispersions of the observed winds are lower (σ ~ 95−190 km s^-1; FWHM ~ 230−460 km s^-1). However, these values are significantly higher than the thermal velocity dispersion of the warm neutral gas (i.e., 8 km^-1, Caldú-Primo et al. 2013). This indicates that a wide range in neutral gas velocity is integrated along the line of sight or that the winds are turbulent and associated with shocks (as seen in the ULIRG F10565+2448 by Shih & Rupke 2010). In the standard GW scenario, (e.g., Heckman et al. 2000) broadening effects are consistent with outflowing gas in interaction (e.g., via shocks) with the surrounding material and with the presence of turbulent mixing layer on cloud surface.

In our sample, outflows are in many cases consistent with being along the minor axis of the ionized gas rotation (Appendix A). Thanks to the present IFS data, we are also able to infer the morphology outflow in about half of the sample (see figures in Appendix A). These outflows appear to be extra-planar, conical (the projected area is triangular in shape), and extended on kpc scales, which is consistent with the expectations from the standard GW model (Heckman et al. 2000). With the present data set, it is unpractical to develop a customized and more detailed model invoking superbubbles or different geometries (which is beyond the scope of this paper). We, therefore, assume a simple outflow model where the GWs are emerging perpendicular to the disk in all cases.

The measured outflow extension values (Table 3) are in good agreement with those reported previously for neutral gas GWs in [U]LIRGs (e.g., Veilleux et al. 2005). However, because of the lack of a bright continuum at large radii, the outer regions of the outflows may not have been detected with absorption line techniques (e.g., via NaD), and thus the quoted extensions should be considered lower limits.

In those outflows with a cone-like morphology, we identify the cone apex (near the galaxy nucleus) and boundaries in the spectral maps, and then we measure the wind 3D opening angles (i.e., C_Ω; see Table 3). We find C_Ω values are in the range between 0.1 and 0.6, in agreement with those (indirectly) estimated from the wind detection rate by Rupke et al. (2005c) and Veilleux et al. (2005) in local starburst and LIRGs (i.e., C_Ω ~ 0.4). These are, however, lower than those inferred by these authors in local ULIRGs (C_Ω ~ 0.8). Our results do not support a correlation between C_Ω and SFR (or the V_w/L_IR, Fig. 9), although we note that the uncertainties associated with the determination of C_Ω are large.

While in half of the cases it is not possible to infer the morphology of the outflows, it is relatively straightforward to estimate the wind projected area for all 22 GWs detected (Table 3). On average, the projected area of the neutral winds is 5 kpc². A spatially resolved broad Hα emission is seen in many objects of our sample (Bellocchi et al. 2013), but its presence shows no obvious correlation with the observed neutral GWs. Specifically, the neutral outflows are more extended with respect to the Hα broad component in 15 of the cases. Considering these cases, the areas covered by the outflow and the broad Hα component generally only overlap partially (12 cases).

4.5. Galactic winds feedback in 2D

While it is relatively straightforward to demonstrate that a GW is present, it is more difficult to robustly calculate the rates at which mass and energy are transported out by the wind (i.e., its feedback effects). We adopted a free wind (FW) model to quantify the neutral wind feedback in the form of outflowing mass, outflow mass rate, and loading factor. Details of the model are given in Heckman et al. (2000), Rupke et al. (2002, 2005c).

Briefly, this model consists of a starburst surrounded by thin shells of a free-flowing wind, a shocked wind (i.e., GW ionized phase), and entraining clouds of neutral ISM. One drawback of this model is that a thin shell could be broken leading to the formation of filaments. Indeed, a starburst-driven wind, in its expansion through the ISM, entrains clouds that are denser than the ambient ISM. A dense cloud embedded in a subsonic flow experiences pressure differences along the surface. This may lead to Rayleigh-Taylor instabilities and the fragmentation of clouds with the consequent formation of filaments that are not accounted for by the FW model (Heckman et al. 2000). Although this scenario is possible, the spatial resolution of the VIMOS data does not allow us to probe possible filaments and, therefore, the shell is a natural and simple interpretation of our observations of GWs. Other approaches are, however, possible (e.g., single radius FWs or superbubbles; Rupke & Veilleux 2013).

To estimate how much of the neutral gas is brought out into the outflow, knowledge of the column density of the wind (N_H) is required. To obtain spatially resolved column densities, we used a method that relates N_H to the strength of the NaD absorption (i.e., its equivalent width) via reddening (E_{B − V}). Specifically, first, Turatto et al. (2003) via a light-curve fitting of low resolution spectroscopy of SNe found that the EW of the NaD correlates with the reddening well (see also Veilleux et al. 1995). Second, Bohlin et al. (1978) found that the E_{B − V} follows a well-defined linear relation with N_H, analyzing a large sample of galactic stars. Combining these relations³, we found the following dependence of N_H, with EW(NaD): (1)The median column density in our sample is 4.2 × 10²¹ cm^-2 (see Table 3). This approach, which has also been followed in Davis et al. (2012) and Cazzoli et al. (2014), is applied here on a spaxel-by-spaxel basis. However, this ignores issues such as the spatial variation of the ionization state of NaD (Murray et al. 2007), metal depletion on dust grains, line saturation effects (e.g., the square-root part of the curve of growth), and variations in gas-to-dust ratio (Wilson et al. 2008).

Alternatively, column densities may be computed following Hamann et al. (1997). By using spatially resolved and spatially integrated data, we found values in rather good agreement with the current estimate and previous works (e.g., 1.6 × 10²¹ cm^-2 by Rupke et al. 2005c). Specifically, we found average values of 4.6 × 10²¹ cm^-2 and 2.0 × 10²¹ cm^-2 for 2D and 1D cases, respectively.

Additionally, we also estimated N_H with the EW-ratios (R_NaD) of the (1D) integrated spectra (after removing the stellar contribution) via the sodium curve of growth (Draine 2011). For this, we only considered optically thin (R_NaD ≥ 1.5, Rupke et al. 2005a) outflowing components (14/24 cases), obtaining a median value of 1.1 × 10²¹ cm^-2. This is smaller than the results obtained from Eq. (1) quoted above and the typical value found by Rupke et al. (2005c).

For the present paper, we consider the values obtained from Eq. (1) which are in good agreement with those values from the method described in Hamann et al. (1997). A more detailed estimation of N_H (e.g., from more complex line-profile modeling as in Rupke et al. 2005a) is beyond the aim of this work.

Our IFS observations probe the wind shape in half of the cases (as discussed in Sect. 4.4) allowing us to constrain the radial extent and variation of the velocity of the outflow with position. We customized the FW model, allowing each spatial element of the wind to have its own velocity and distance from the wind’s origin. We used the column density and velocity measured in each spatial resolution element (k) to derive the wind mass (M_w) and outflow rate (), following Rupke et al. (2005c), as These equations describe the mass and mass outflow rate for a GW flowing into a solid angle C_Ω,k with a cloud covering factor C_f. We consider the observed GWs as a series of thin shells, each one located at radius R_w,k with the corresponding inclination-corrected velocity V_w,k and column density (N_H,k) inferred, as previously discussed. The extent R, median velocity V and average column density N_H, as well as the mass and mass outflow rate, which are derived from our observations assuming the FW model, are listed in Table 3. The parameter C_f is assumed to be 0.37 as in Rupke et al. (2005c).

As mentioned in the previous section, we cannot infer the wind morphology in ten cases; this is mainly because of projection effects (Table 3). In these cases, to calculate the wind mass and mass rate, we assumed a fiducial radius of 3 kpc and a wind opening angle of 0.4, which are the characteristic values seen for the subsample for which we can derive the actual wind morphology. We also considered the median values of V, and N_H distributions measured over each GW region.

The wind mass estimates range between 4 × 10⁷ M_⊙ to 7.5 × 10⁸ M_⊙ with a mean value of 1.6 × 10⁸M_⊙. These values are roughly in agreement (lower by a factor of 3 when comparing averages) with those reported by Rupke et al. (2005c) obtained for neutral outflows in their local LIRGs sample (i.e., 6.3 × 10⁸ M_⊙ on average).

Previous long-slit studies of Rupke et al. (2005b,c) and Martin (2006), on the basis of the FW wind model, found that the mass of neutral outflows in [U]LIRGs can be up to 10% of the dynamical mass (estimated via CO measurement). For our sample, the wind mass is typically ~3% of the dynamical mass of the galaxy. These results indicate that neutral outflows in [U]LIRGs may carry away a significant amount of gas mass that would otherwise be available for further star formation (Sato et al. 2009).

However, in general, we found that the mass outflow rate (i.e., ) is not larger than the global SFR (Table 3). Specifically, the outflow loading factor (i.e., η = ) is η< 1 in nearly all objects (Table 3), indicating that the mass-loss rates (in the neutral phase) are relatively unimportant for quenching the SFR. However, an estimation of the mass-loss rates, which take into account other wind phases (e.g., molecular) may lead to an increase of the significance of the feedback.

These results are in rough agreement with previous measurements of η in nearby [U]LIRGs (Rupke et al. 2005c) and empirical models (e.g., Zahid et al. 2012) for which η ≤ 1. Values of η above the unity, which indicates a strong mass loading of the outflow by the galaxy ISM, have been found in active galaxies (Veilleux et al. 2005).

In the left panel of Fig. 10, we show the loading factor as a function of the infrared luminosity, which are in negative correlation. Excluding those galaxies for which we were unable to estimate the wind morphology in detail, the Paerson coefficient, r_PC hereafter, is ~–0.6. A similar trend was found by Rupke et al. (2005c), although these authors used the K-band magnitude as a tracer of the infrared luminosity. Nevertheless, considering the same subsample, our analysis also supports the absence of any dependence of η on dynamical masses of the host for the neutral phase (r_PC< 0.1); this is seen in the right panel of Fig. 10, where we present the dynamical mass versus the mass loading factors. This is in contrast to the inverse dependency on η and the dynamical mass found for the ionized phase of GWs in 32 LIRGs without AGNs (Arribas et al. 2014). We note, however, that the present sample (12 objects) is relatively small and it includes 6 galaxies hosting weak AGNs.

In summary, the present results indicate that the neutral gas loading factors are small (<1) and, therefore, the feedback effects are not expected to be large if we only take the neutral gas phase of the outflows into consideration. For the objects with a well-defined morphology, the loading factors show a trend with the infrared luminosity, but they seem uncorrelated with the host mass. The relative small number statistics (i.e., 12 objects, including 6 weak AGNs) prevent us from considering the latter trend as a firm conclusion.

4.6. The gas cycle: IGM and ISM pollution via outflows

The GWs that originated in the disk of spiral galaxies, in principle, could have the ability to eject large percentage of the cold gas reservoir available for star formation into the IGM (Veilleux et al. 2005). We compare the wind velocity to the escape velocity of each galaxy to investigate the fate of the outflowing material.

Following a simple recipe (Arribas et al. 2014 and reference therein) for a galaxy with dynamical mass (M_dyn), we estimate the escape velocity (v_esc) at r = 3 kpc (considering an isothermal sphere truncated at R_max as the gravitational model for the host galaxy) as (4)where G is the gravitational constant (i.e., 4.3 × 10^-3 pc M_⊙ (km s^-1)²) and R_max/r = 10. This approach allows us to take into account both the dependence of the dynamical mass estimation to rotation and dispersion motions and assumes that the halo drag is negligible. We refer to Bellocchi et al. (2013) for a detailed description of the dynamical mass estimation.

For the two galaxies (F05189-2524 and F07027-6011 (N)) for which the dynamical mass estimation was not possible as a result of AGN contamination (Bellocchi et al. 2013), we calculate v_esc, as in Veilleux et al. (2005) as follows: (5)Similar to the previous case, this simpler approach assumes a truncated isothermal gravitational potential and no halo drag. We consider that the Hα velocity amplitude, which is defined as the half of the observed peak-to-peak velocity corrected for the inclination of the galaxy, (Bellocchi et al. 2013) is a good proxy for v_rot and R_max/r = 10 (as in the previous case).

Table 4 lists the escape velocities and outflow velocity in units of the escape velocity i.e., V_w/v_esc. This ratio ranges from 0.2 to 4.2 (Fig. 11) with a median (average) value of 0.9 (1.1).

There are ten objects for which V_w/v_esc> 1 and, therefore, a significant amount of their outflowing gas could pollute the IGM (Fig. 11, right). The most extreme case is F06206-6315, where an AGN is likely playing a role in boosting the velocity of the outflowing gas to exceed the escape velocity (see also Sect. 4.7). Excluding this object, the sample indicates a (weak) correlation in the sense that V_w/v_esc is higher for the less massive galaxies as found for the ionized outflows (Arribas et al. 2014). This result is consistent with the prediction that the less massive galaxies have the most severe impact on the IGM pollution.

In the majority of the cases, we find that the median wind velocities do not exceed the escape velocities (i.e., V_w/v_esc< 1, Table 4), indicating that most of the cold outflowing material that we probe with our NaD measurement is likely falling back to the disk. This circulation of the gas could in principle favor the redistribution of metals, modifying the abundance gradients in the galaxy disk (Spitoni et al. 2010).

This fountain scenario is consistent with the presence of clouds of neutral gas, either above the main disk in disk-like layers (in slow rotation) or sparsely distributed (i.e., raining back showing, in projection, an irregular velocity field). In the latter case, the neutral gas is still raining back in the form of high velocity clouds (HVCs; Spitoni et al. 2013), which are gravitationally bound to the host galaxy, but not virialized.

However, V_w refers here to the median velocity in the outflow and, therefore, even for values of V_w/v_esc< 1, part of the gas entrained in the wind could in principle escape.

The efficiency of the observed winds of polluting the IGM, is in any case, a lower limit. Specifically, by using absorption lines as tracers of the outflowing gas we are limited to observing gas at a distance smaller than the projected size of the stellar disk. Tenuous and hot material (in absence of radiative cooling) or molecular gas, may be found at larger radii (up to 10 Kpc in M82, Lehnert et al. 1999) and is more likely to escape.

In Fig. 11 we consider the ratio V_w/v_esc as a function of dynamical mass and dynamical ratio of the ionized disk (i.e., V/σ, from Bellocchi et al. 2013) in order to study the factors controlling the recycling of gas and metals. As expected, gas and metals are more likely to escape to the IGM (higher V_w/v_esc) in galaxies with small potential wells.

In addition, we note that there is a tendency for V_w/v_esc to increase in turbulent and thicker ionized gas disk (i.e., higher disk height and lower dynamical ratio, Fig. 11, right). In such disks, the gas turbulence may help the wind to escape, thereby enhancing its vertical growth. This is in agreement with a scenario in which the wind fluid follows the path of least resistance (Cooper et al. 2008).

4.7. Wind engine and energetics

The outflow kinematic properties are more extreme in galaxies that host a powerful AGN (Veilleux et al. 2005). Therefore, the outflow velocity can be used as a discriminant of starburst-driven versus AGN-driven winds as suggested by Rupke et al. (2005) and Rupke & Veilleux (2013). The influence of AGNs on the velocity (and power) of outflows has also been seen in ionized (Arribas et al. 2014) and molecular (Cicone et al. 2014; Sturm et al. 2011) wind phases.

As mentioned in Sect. 4.3, the ULIRG F06206-6315 likely hosts an AGN-driven GW, according to the observed outflow kinematics (Table 3). In addition to this extreme case, there are ten other objects where evidence for the presence of an AGN has been found according to X-ray and infrared observations or optical emission line ratios and kinematics (Table 1). Despite this, from the observed neutral wind kinematics in these objects (Table 3) it is not obvious that the AGN is significantly boosting the outflows; this is probably because the AGN is not powerful enough.

To further address the issue about the wind driver, the energy (E_w) and power (P_w) of the winds are calculated, using the shell formalism and the FW model outlined in Sect. 4.5, as follows: where σ is the velocity dispersion and b, the Doppler width (defined as: ; Rupke et al. 2005c) for each spaxel. As in Sect. 4.5, for the ten cases for which we cannot estimate the wind morphology, we measured V, N_H, and b on spaxel-by-spaxel basis and used their median values (calculated over the area where the wind is observed). In addition, we assumed a wind radius and opening angle of 3 kpc and 0.4, respectively.

Excluding the case of F06206-6315, the estimated wind energies are in the range from 3 × 10⁵⁵ erg to 1.2 × 10⁵⁷ erg, while the powers are in the range between 6 × 10⁴⁰ erg s^-1 and 1.4 × 10⁴³ erg s^-1. The median energy and power are 4 (2) × 10⁵⁶ and 14 (6) × 10⁴¹ erg s^-1, respectively, excluding (including) the winds because it was not possible to determine their morphology.

For comparison, Rupke et al. (2005c) found that the typical wind energy and power in local LIRGs are 5 × 10⁵⁶ erg and 4 × 10⁴¹ erg s^-1, respectively. The average wind energy in the present work is generally consistent within errors with that of the above mentioned work, although our estimation of the wind power is generally larger (by a factor of about three).

To investigate the wind origin, in Fig. 12, we compared the wind power of the neutral outflows and the kinetic power of the starburst associated with SNe (i.e., KP). We calculated the latter with the recipe of Veilleux et al. (2005): KP ~ 7 × 10⁴¹ SFR/M_⊙ yr^-1. In nearly all of the cases, the wind power is within 1% and 20% of the kinetic power supplied by the starburst (about ~5% on average). This is in agreement with a scenario in which these winds are originated in dense core of powerful nuclear starbursts (Veilleux et al. 2005). However, our result is in partial contrast with hydrodynamical simulations, which often assume wind thermalization efficiencies in the range 10−100% (e.g., Strickland & Stevens 2000) and direct measurements of the thermalization efficiency. For example, for the nearby starburst M82, Strickland & Heckman (2009) found that medium-to-high thermalization efficiencies (>30%) are required in hydrodynamical models to match the set of observational constraints derived from hard X-ray observations.

We see only two cases, for which an unlikely thermalization efficiency of ~90−100% could be required. For the ULIRG F06206-6315, the presence of an AGN suggests that its outflow is rather driven by the energy liberated by the accretion of gas into the black hole. For the LIRG F06592-6313, an AGN-driven outflow could be a reasonable possibility, although there is not evidence for this outflow at optical wavelengths. However, we note that the large uncertainties associated with P_w also make these cases consistent with lower thermalization efficiencies. For the other galaxies, the AGN contribution is in principle not required to explain the observed energetics.

4.8. Thick disks in slow rotation

The velocity fields of the neutral gas observed in 11 [U]LIRGs have the spider diagram pattern characteristic of a rotating disk (figures in Appendix A; Cazzoli et al. 2014 for the case of the LIRG F11506-3851). In three galaxies (F07027-6011 (S), F10409-4556, and F11506-3851), the rotation is detected in one kinematic component, but an additional outflowing component is also detected (Figs. A.13 and A.23, and Cazzoli et al. 2014). However, while in case of an ideal rotating disk the velocity dispersion map should be centrally peaked, the observed neutral gas velocity dispersion maps are generally rather irregular.

The observed NaD feature in these 11 galaxies shows a wide range of stellar contribution (Table 2), though in most of the objects (8/11) the stellar absorption is not dominant (i.e., <50%). For the three galaxies where the NaD absorption is dominated by the stellar contribution (F09437+0317 (S), F12043-3140 (S) and F12115-4656), the rotation curves may be difficult to interpret in terms of neutral gas motions. However, it is important to consider that the stellar contribution is computed in the integrated light. Therefore, in some cases even if the stellar absorption in the integrated spectrum is substantial, the 2D distribution may help to identify regions where the NaD absorption is dominated by neutral gas or the stellar component, as in F11056-3851 (Cazzoli et al. 2014). Without the knowledge of the stellar properties and kinematics in these galaxies, we cannot investigate this further.

From our spectral maps, we extracted the velocity and velocity dispersion values for both the neutral and ionized (i.e., Hα narrow disk-like component) ISM phases in a ~10 pseudo-slit along the NaD major kinematic axis. The corresponding position-velocity diagrams (PV diagrams) are plotted in the upper panels of Fig. 13, and show a variety of shapes for the rotation curves of both ionized and neutral gas components in the inner regions (typically within 3 to 8 kpc). Similarly, in the lower panels of that figure, the velocity dispersion radial profiles are shown.

We found that for 9 out of 11 galaxies, the rotation axes of the neutral and ionized gas disks are fairly aligned, with offsets smaller than 15° (Table 5), although in two cases (F10409-4556 and F09437+0317 (S)) the kinematic axes of the neutral gas disks are poorly constrained. The neutral gas is found in slower rotation than the ionized gas in the majority of the cases (i.e., 8/11), while the neutral gas seems to corotate with the ionized gas in only two cases. In addition, in F06259-4780 (C) the neutral gas disk is observed in counter-rotation with respect to the ionized gas disk. In this context, it is interesting to note that a counter-rotating neutral gas component was also found in the inner few kpc of M82 (Westmoquette et al. 2012).

The radial profiles of the velocity dispersion are either flat (e.g., F18093-5744 (S)) or with considerably large deviations from what is expected for a thin rotating disk at large radii. These deviations are particularly evident in the case of F22132-3705 (with values in the range of 110−200 km s^-1).

In Fig. 14, we show a comparison of the ionized and neutral V/σ ratios of the ISM disks. We found that the neutral disks are dynamically hotter (i.e., low V/σ values, Table 5) than the ionized disks and, therefore, they are likely to be thicker (i.e., larger h_z). This confirms that the ionized gas resides in regions of high density close to the innermost regions of the disk, while the neutral gas is located further out.

Therefore, to estimate the h_z of the disks we follow two different approaches, which are outlined in Cresci et al. (2009) and Binney & Tremaine (2008). Specifically, for the ionized disks we used the approximation proper for thin disks, for which h_z can be computed as (8)For the neutral disks, however, we used an approximation more suitable for thick disks that is written as (9)In both cases, we considered as ΔV the (inclination-corrected) semiamplitude of the velocity field, as σ the mean velocity dispersion across the galaxy disk (excluding the nuclear regions) and a nominal radius (R) of 2 kpc. We chose this radius since it is the typical distance at which the rotation curve starts to flatten (Fig. 13). The results are summarized in Table 4. We find that the neutral gas disks are thicker by a factor up to 46 (~15 on average) than the ionized gas disks.

When the same approach (formula) is used for the neutral and ionized disks, the results indicate that the neutral disks are still significantly thicker (i.e., by a factor 8 on average) than the ionized disks.

5. Summary and conclusions

We have studied the properties of neutral gas outflows in a sample of 51 local [U]LIRGs (z ≤ 0.09) on the basis of VLT/VIMOS IFS of the NaD feature. For the analysis, we followed two approaches. First, for each galaxy we combined the spectra in the data cube to obtain a high S/N spatially integrated spectrum. Second, for a subsample of [U]LIRGs, we analyzed the NaD spectra on a spaxel-by-spaxel basis to trace the spatially resolved 2D structure of the neutral gas. The main conclusions can be summarized as follows:

• 1.

  Stellar and ISM contributions to the NaD feature. We evaluated the contribution of old stars to the NaD absorption via a stellar continuum modeling (pPXF analysis) of the spatially integrated spectra. The fraction of stellar contribution ranges from about 20% to 100%, but in nearly half of the cases for which high quality modeling was possible, the NaD is mainly originated in the ISM (i.e., stellar fraction <35%). This analysis also allows us to derive a threshold of EW(NaD) = 1.3 Å, above which the NaD absorption in [U]LIRGs is likely dominated by the ISM. Since the results from the stellar continuum modelling on a spaxel-by-spaxel basis could be uncertain, we consider that threshold to identify regions dominated by ISM absorption in the maps.

• 2.

  Outflow kinematics from the integrated spectra.  For the objects with a reliable stellar modeling, we generate a purely ISM NaD absorption spectrum after subtracting the stellar model. In 22 objects we measure blueshifted NaD profiles, which indicate typical neutral gas outflow velocities in the range 65−260 km s^-1. Excluding the galaxies with powerful AGNs, the neutral outflow velocity shows a dependency with the SFR of the type V ∝ SFR^0.15, which is in fair agreement with previous results. The neutral outflow (central) velocities are significantly higher than those for the ionized gas, but they are in rather good agreement with the ionized gas maximum outflow velocities that we considered. This suggests that the wind entrains and accelerates the cold ambient gas, which is likely located at relatively large distances from the regions where the ionized gas resides. The velocity dispersions for both neutral and ionized outflows show similar median values (if the strongest AGNs are excluded).

• 3.

  2D Neutral outflows. Detection, morphology, and kinematics. In the neutral gas velocity fields of 22 out of 40 targets, we found clear signatures of GWs. These neutral winds are conical in shape in 12 out of 22 objects. For the remaining objects, we were not able to constrain the morphology mainly owing to projection effects. We generally observe collimated outflows with C_Ω ~ 0.4. The inclination -corrected outflow velocities of the neutral gas entrained in GWs is in the range ~80−700 km s^-1, except for the ULIRG F06206-6315 for which the typical wind velocity exceeds 1000 km s^-1. The typical velocity dispersions are in the range ~95−190 km s^-1 (i.e., 230−460 km s^-1 in FWHM), indicating either that a wide range of velocities are integrated along the line of sight or that these winds are turbulent. The V-SFR relation inferred for the 2D analysis is similar to that obtained from the integrated spectra. No clear correlation is seem between the outflow velocity and the dynamical mass.

• 4.

  GWs Feedback. Slowed star formation and gas recycling.  Based on a simple FW model, we found that the wind mass estimates range from 0.4 × 10⁸ M_⊙ to 7.5 × 10⁸ M_⊙ (1.6 × 10⁸ M_⊙, on average), reaching up to ~3% of the dynamical mass of the host. The mass rates are only ~0.2−0.4 times the corresponding SFR in most cases indicating that, generally, mass losses are small for slowing down significantly the star formation. The derived mass loading factors () correlate better with the starburst infrared luminosity (r_PC ~ −0.6) than with the host galaxy mass (r_PC< 0.1). The comparison of the median wind velocity and host escape velocity indicates that, in the majority of the cases, most of the outflowing neutral gas rains back into the galaxy disk. We found that on average V_w/v_esc is higher in less massive galaxies, confirming that the galaxy mass has a primary role in shaping the recycling of gas and metals. We also found a tendency in the sense that V_w/v_esc is higher in more turbulent and thicker ionized gas disks.

• 5.

  GWs power source. The comparison between the wind power and kinetic power of the starburst associated with SNe indicates that the starburst could be the only main driver of the outflows in nearly all the [U]LIRGs galaxies, as wind power is generally lower than 20% of the kinetic power supplied by the starburst. In the case of F06206-6315, the outflow is likely AGN driven, as also indicated by its kinematics.

• 6.

  Kinematics of disks. A significant number of [U]LIRGs (11/40) show spider diagrams, such as neutral gas velocity fields, plus rather irregular velocity dispersion maps in one kinematic component. The comparison between the rotation curves of the ionized and neutral disks indicates that the neutral disk lags compared to ionized disk in the majority of the cases (i.e., 8/11). While in two cases the neutral gas seem to nearly corotate with the ionized gas, we also find a case (IRAS F06259-4780 (C)) in which the neutral gas disk counter-rotates with respect to the ionized disk. Our kinematics measurements indicate that nearly all the neutral gas disks are dynamically hotter and thicker (by a factor up to 46, 15 on average) with respect to the ionized disks.

---

Acknowledgments

This work was funded by the Marie Curie Initial Training Network ELIXIR of the European Commission under contract PITN- GA-2008-214227 and the grants AYA2010-21161-C02-01 and AYA2012-32295 by the Spanish Ministry of Science and Innovation (MICINN). S.C. gratefully acknowledges the logistic and financial support provided by the Cavendish Laboratory (Cambridge, UK). We are grateful to E. Bellocchi for kindly providing us the Hα maps shown in the Appendix. We also appreciate some information and help provided by M. Pereira-Santaella. This work is based on observations carried out at the European Southern Observatory, Paranal (Chile), Programs 076.B-0479(A), 078.B-0072(A), and 081.B-0108(A). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Comments, maps, and integrated spectra of individual sources

The first part of this appendix is devoted to comment briefly on the spatially resolved and/or the integrated characteristics for each galaxy. In the second part, we present the neutral gas equivalent width (EW), velocity field (V), and velocity dispersion (σ) maps for those objects for which we were able to obtain spectral maps (Sect. 2.1). The VIMOS data for the LIRG IRAS 11506-3851 are discussed in detail by Cazzoli et al. (2014) thus the comments and correspondent maps are omitted here.

F01159-4443 (N). The neutral gas kinematics of this galaxy shows deprojected blueshifited velocities (~130 km s^-1) that are larger than those of the systemic component of the ionized gas at the same location. Therefore, the observed neutral gas velocities cannot be explained with rotation and are interpreted as outflowing.

F01159-4443 (S). This galaxy is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD because of the low quality of stellar modeling and lack of any clear GW or disk signature, respectively.

F01341-3735 (N). We found a prominent blue wing in the pure-ISM NaD line profile in the spatially integrated spectrum of this LIRG. In the spectral maps, we identified an outflowing neutral gas component with high velocity and velocity dispersion.

F01341-3735 (S). This galaxy was only studied via the integrated spectrum, where the NaD doublet was modeled with a single kinematic component found at the systemic velocity.

F04315-0840. The neutral gas spatial distribution and kinematics of this post-coalescence late-merger starburst LIRG is very different from that of the ionized gas (traced via the Hα narrow component). The neutral gas is mainly seen in the main body of the object, while it is almost absent in the spiral arm. A neutral wind originates from the nucleus toward the southeast with velocity dispersions values higher than 90 km s^-1. The blue velocities of the Hα broad component (Bellocchi et al. 2013) spatially overlaps the region where the neutral wind is detected, indicating the presence of a multiphase GW. Interestingly, the velocities of both wind phases are similar at ~200−350 km s^-1. In the integrated spectrum, we found a blueshifted kinematic component at V ~ 190km s^-1 in agreement with our spatially resolved results.

F05189-2524. This merger ULIRG has a very compact nucleus and hosts an AGN that powers the high velocity outflow seen via NaD. The morphology of the outflow could not be constrained well owing to projection effects. A noticeable blue wing is seen in the NaD ISM-absorption line profile in the integrated spectrum, along with the NaD resonant emission line (which is not a common feature of GWs in nearby galaxies) that is studied in detail by (Rupke & Veilleux 2015).

F06035-7102. This galaxy is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD because of the low quality of stellar modeling and lack of any clear GW or disk signature, respectively.

F06076-2139 (N) and (S). We only study these galaxies via their spatially integrated spectra. In both spectra, the NaD line profile were modeled with one broad and blueshifted kinematical component. However, the complex NaD line profile observed in the spectrum of the northern galaxy likely require more than one (or two) kinematical component(s).

F06206-6315. This ULIRG has two nuclei, which are clearly seen in the near-IR images (Rodríguez-Zaurín et al. 2011) but not in the optical continuum image. The tidal tail, starting in the north and bending toward the southeast, is almost undetected via NaD. While we did not analyze the low S/N integrated spectrum, in 2D, we observed part of the neutral gas as entrained in an AGN-driven outflow. The comparison with the galaxy escape velocity and wind velocities indicates that the neutral gas entrained in the wind is able to escape. A molecular gas outflows, traced by the infrared OH lines, was found by Spoon et al. (2013) at higher observed (terminal) velocity (~800 km s^-1).

F06259-4780 (N). In this galaxy, we detected extended neutral gas that is outflowing nearly everywhere in contrast to what is seen for the ionized gas (narrow component). In addition, neutral and ionized ISM phases also have remarkably different spatial distribution. We identify the GW according to the most blueshifted velocity, the highest velocity dispersion values, and also its direction. Indeed, the outflow is oriented along the minor axis of the rotation seen weakly in the ionized velocity field.

F06259-4780 (C). The 2D-kinematic of the neutral gas in the central galaxy of the system (the southern in Fig. A.8) is classified as rotating disk. This gas is observed in counter-rotation with respect to the ionized disk. This makes this object unique in our sample since, in the nearly all cases, we found slow neutral gas rotation with respect to the ionized gas.

F06259-4780 (S). The southern galaxy of this system is not analyzed in this work. Difficulties in modeling the observed weak NaD absorption in the spatially integrated and spatially resolved data of this galaxy are partially due to its edge-on orientation.

F06295-1735. This LIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD. The quality of the stellar modeling is low and a large number of spectra in individual spaxels were not suitable for the analysis.

F06592-6313. The integrated spectrum of this LIRG shows a strong NaD absorption in which the stellar contribution is rather low (25%). We modeled the purely ISM NaD line profile with one blueshifted kinematic component. In the spectral maps, toward the south, we also observe strongly blueshifted absorption doublets (with velocities up to 500−600 km s^-1), which we interpret as a signature of GW. The broad component of the Hα emission line (Bellocchi et al. 2013), seen with velocities up to 500 km s^-1, is rather oriented with the wind well (PA ~ 130°) and partially overlaps the area covered by the neutral wind. However, the overall morphology of the broad Hα component is poorly constrained and it is unclear how it is related to NaD.

F07027-601 (N). The velocity map of this galaxy shows outflowing neutral gas with velocities up to ~300 km s^-1 extended within 1.5 kpc of the nucleus. Despite evidence of an AGN, found by Arribas et al. (2014) according to its optical spectrum, we found no strong evidence that the outflow we detect is boosted by the AGN. The morphology of the outflow could not be constrained well owing to projection effects. Excluding the outflowing gas, the NaD is seen with an irregular spatial distribution, velocity field, and velocity dispersion map. This gas in not-ordered motion could be gas falling back onto the disk in HVCs.

F07027-601 (S). This southern galaxy is an important case to study because its NaD IFS spectra, when fitted with multiple components, provide the possibility of disentangling the ordinary rotation and outflow. One kinematic component was found with rather regular rotation pattern, but lagging with respect to the kinematics of the ionized gas (i.e., Hα-narrow component) with a small positional misalignment. The other (outflowing) kinematic component is evident within ~1.3 kpc of the nucleus toward the northeast. This neutral gas component is broad (σ ~ 140 km s^-1, on average) and has (observed) velocities within 100−300 km s^-1. The morphology of the outflow could not be constrained well owing to projection effects.

F07160-6215. The velocity field of the neutral gas as a disk-like pattern in the main body of this dusty edge-on LIRG as seen for the ionized gas (traced via Hα  and Brγ λ2.1655 μm; Bellocchi et al. 2013 and Piqueras López et al. 2012, respectively). The velocity dispersion is almost flat, excluding the region toward the northeast, where a broad and blueshifted NaD component was found and interpreted as a signature of GW. However, no ionized wind-phase counterpart is seen in VIMOS and SINFONI IFS data (Bellocchi et al. 2013; Piqueras López et al. 2012).

08355-4944. This LIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD. The very weak NaD absorption, seen in both the integrated and the spatially resolved spectra, prevents any robust study.

08424-3130 (N) and (S). These spiral galaxies in interaction are only partially observed in our VIMOS data. For this reason, we did not study these galaxies in detail even though in the northern galaxy, we detect neutral gas with blueshifted velocities which are higher than those of the ionized gas (i.e., Hα-narrow component) at the same location, which can be interpreted as outflowing.

F08520-6850. This LIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS. The quality of the stellar modeling is low and a large number of spectra in individual spaxels were not suitable for the analysis.

09022-3615. This ULIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD. The quality of the stellar modeling is low and the spectral maps lack of any clear GW or disk signature.

F09437+0317 (N). Two pointings (i.e., northeast (NE) and northwest (NW); Rodríguez-Zaurín et al. 2011) sample the northern galaxy. However the NW is not useful because of the low S/N of the spectra of the data cube. In the spectral maps obtained from the NE data cube, we observe that the neutral gas disk kinematic shares the same major kinematic axis with the ionized gas rotation, although the neutral gas is found in slower rotation compared to the warm ionized gas (Bellocchi et al. 2013).

F09437+0317 (S). In the southern galaxy (observed with one VIMOS pointing) we observe a shared major kinematic axis and a high kinematical coupling within the neutral and ionized ISM phases. The neutral gas rotation is only slightly slower (of about 40 km s^-1, on average) than that measured for the ionized gas at the same distance.

F10015-0614. The NaD-clouds spatial distribution only partially reproduce the spiral structure of this LIRG. No hints of rotation pattern are found in the neutral gas velocity field and velocity dispersion map. On the contrary, we found blueshifted velocities up to 200 km s^-1 along the kinematic minor axis of the ionized gas rotation. We did not consider as part of the outflow the NaD-clouds toward the north and northeast because of their overlap with the approaching side of the ionized gas rotation pattern. The morphology of the outflow could not be constrained well because of projection effects.

F10038-3338. Our data for this LIRG, show deep, broad, and blueshifted NaD line profiles in the integrated spectrum and spectral maps (toward the northeast). We interpret this feature as a GW.

F10257-4339. This merger LIRG is the nearest object in our VIMOS-IFS survey. The neutral outflow is observed out to distances of 1.4 kpc. We did not consider the NaD-clouds (with negative velocities) toward the north and northwest as part of the outflow, since we observe an irregular and asymmetric NaD line profile in the IFS spectra. We tested two-component modeling, but this modeling gives ambiguous results.

F10409-4556. The observed NaD feature in this LIRG shows generally complex absorption line profiles, best fitted with two absorption doublets. We considered the most negative velocities as part of the outflow, which have nearly no spatial overlap with the approaching side of the ionized gas rotation. Toward the north, this fast kinematic component associated with a GW dominates the NaD profile. However, the morphology of this outflow could not be constrained well because of the ambiguity in the identification of the region where it is originated. The outflow is also seen via the integrated spectrum.

F10567-4310. The NaD absorption seen in the integrated spectrum and spectra of the IFS cube of this LIRG galaxy is significantly deep and mainly interstellar in origin. Indeed, our stellar continuum modeling indicates a rather low stellar contamination to the NaD doublet (28%). The neutral gas spatial distribution is very compact and the neutral gas velocity field is dominated by a starburst driven GW with outflowing velocity ranging 100 to 300 km s^-1. The morphology of the outflow could not be constrained well due to projection effects.

F11255-4120. The observed NaD integrated line-profile in this LIRG show the lowest stellar contamination of the full sample (i.e., 18%). The purely-ISM line profile was modeled with two kinematical components both having blueshifted velocity. Despite of that, the observed NaD line profile in our IFS-data is modeled with one kinematic component well. Our NaD IFS-maps show evidence of a GW, in the inner region (within R> 1.2 kpc), that consist in: large blueshfited velocities (up to 300 km s^-1) and high turbulence (σ (NaD) up to 280 km s^-1). The morphology of the outflow could not be constrained well owing to projection effects. Weak evidence for an ionized wind counterpart was found (Bellocchi et al. 2013).

F12043-3140 (N). This galaxy is not analyzed in this work. The stellar continuum modeling of the integrated spectrum is rather uncertain. Our spectral maps indicate a generally blueshifted velocity pattern and rather low values for the velocity dispersion and EWs. Therefore, the evidence for claiming a GW detection are weak.

F12043-3140 (S). The absorption due to stars contributes significantly (77%) to the NaD line profile, so the rotation pattern seen in this VIMOS-spectral map might be tracing the stellar kinematics. This result is also supported by the low depth of the absorption found in the EW map. However, additional observations are needed to confirm this interpretation. The 2D neutral gas disk rotation is slower with respect to that of the ionized gas and the major kinematic axes are not aligned well. The high velocity dispersion values seen toward the north are possibly due to additional turbulence this outer region.

Along the southern minor axis there is a region with high blueshifted velocities, and high dispersions that may indicate the presence of an outflow. However, the values of high velocity dispersion are found in just a few spaxels, likely due to low S/N, as the general region has low values. Additionally, in the region of the putative wind, the average value of the residual maps (i.e., V_NaD − V_Hα) is slightly lower then the limit we consider. Therefore, we did not include this case from the sample with detected outflows.

F12115-4656. The velocity field is fairly regular with the kinematic center in positional agreement with the continuum and Hα flux peaks. The neutral gas disk is ~45 times thicker than that of the neutral gas. This LIRG was already studied in great detail by Arribas et al. (2008). We do not excluded the presence of a GW oriented along the kinematic minor axis, but it is not the main feature either in the spectral maps or in the integrated spectrum (where the stellar fraction to the NaD is rather high, 73%).

F12116-5615. The NaD line profile seen in this LIRG (both in the integrated spectrum and in the spectra of individual spaxels) was fitted with two kinematic components. One kinematical component was found dominated by not-ordered motions and the other by a neutral GW via IFS. The wind is seen almost face-on preventing us from measuring its extent and opening angle. The neutral gas entrained in the wind is likely swept out polluting the IGM. The GW is likely multiphase, since the Hα broad component is seen as outflowing (Bellocchi et al. 2013).

F12596-1529. This LIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD. The quality of the stellar modeling is low and a large number of spectra in individual spaxels were not suitable for the analysis.

F13001-2339. The ionized gas velocity field and velocity dispersion map (traced via the Hα narrow component) have an irregular pattern and the neutral gas is seen mainly as outflowing. The outflow spatially overlap the region where the ionized gas velocity field has positive velocities. The morphology of the outflow could not be constrained well because of projection effects and the lack of any rotation pattern of reference.

F13229-2934. This LIRG hosts a neutral gas wind detected with both the spatially integrated and spatially resolved spectra. As for the galaxy F04315-0840, an ionized wind counterpart was observed (Bellocchi et al. 2013). This has roughly the same velocities observed in the present neutral gas velocity map. Even though this galaxy hosts an AGN, the comparison of the wind power and kinetic power of the starburst indicate that the multiphase GW is likely driven by SNe.

F14544-4255 (E). We found that the NaD originated in stars gives a modest contribution (41%) to the integrated line profile for this galaxy, but the modeling is rather uncertain. The neutral gas velocity field is disk-like and the kinematic center seems to be in positional agreement with the continuum peak. The rotation seen for the neutral gas lags compared to that of the ionized gas. For both ISM phases we did not observe a central peak in the velocity dispersion maps (see also Bellocchi et al. 2013); this may indicate that the contribution of the bulge component is negligible.

F14544-4255 (W). The western galaxy has not been classified either in 1D or in 2D. In the spectral maps, the EWs(NaD) are generally low and we found no clear rotation pattern or evidence indicating a GW.

F17138-1017. This LIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD. The quality of the stellar modeling is low and a large number of spectra in individual spaxels were not suitable for the analysis.

F18093-5744 (N). The purely ISM NaD profile seen in the integrated spectrum of this galaxy is modeled with one blueshifted kinematic component well. However, in the spectral maps the putative wind (likely seen toward the northeast) has both low velocities and EWs(NaD). This characteristic makes its detection is rather uncertain.

F18093-5744 (C). This galaxy is not analyzed in this work since it was not possible to make a kinematical classification via both the spatially-resolved and spatially-integrated spectroscopy of the NaD absorption.

F18093-5744 (S). The spider pattern seen in the neutral gas velocity field of the galaxy is fairly regular. The rotation is seen for a disk three times thicker and slower with respect to that of the ionized gas. The stellar contamination, estimated with our stellar continuum modeling of the integrated spectrum, is modest (41%).

F21130-4446. This ULIRG is not analyzed via NaD in this work. The NaD absorption is very weak both in the integrated and spatially resolved spectra.

F21453-3511. In this galaxy, a GW at rather slow velocity is seen emerging from the nucleus. The outflow is also seen via a blue wing in the NaD integrated line profile (despite its strong stellar contamination, i.e., 97%). The morphology of the outflow could not be constrained well because of a patchy distribution of the NaD outflowing clouds and projection effects. Only a weak indication of the presence of an ionized wind was found (Bellocchi et al. 2013). Excluding the region where the wind is detected, the neutral gas velocity field is rather irregular and roughly follows the kinematic of the Hα disk but with strong deviations (up to 150 km s^-1).

F22132-3705. The disk-like velocity field is fairly regular with the kinematic center in positional agreement with the continuum and Hα flux peaks. The neutral gas disk is ~20 times thicker than that of the neutral gas. The stellar fraction to the NaD is modest (47%). The high values found in the velocity dispersion map toward the northeast (particularly evident in the velocity dispersion radial profile) can be explained by either addition turbulence in these outer regions or GW. The evidence for claiming a GW detection however is weak.

F22491-1808. This ULIRG is not analyzed in this work either with the spatially integrated spectrum or via spatially resolved IFS of NaD. The quality of the stellar modeling is low and a

large number of spectra in individual spaxels were not suitable for the analysis.

IRAS-F23128-5919 (N) and (S). Despite the high integrated stellar contamination, the NaD velocity field of the two galaxies does not show any ordered rotation pattern typically seen, e.g., the stars. Most of the cold gas associated with this system is outflowing reaching velocities up to 500−600 km s^-1. The fast neutral outflowing gas has a conical shape (with the apex in the southern galaxy) as also seen via the Paα and Hα line profiles (Piqueras López et al. 2012; Bellocchi et al. 2013). For this system, the capability for the wind of polluting the IGM is doubtful because of the inability to estimate the contributions to the velocity field of each individual galaxy reliably (Bellocchi et al. 2013), which leads to additional uncertainties for the calculation of v_esc.

All Tables

All Figures