© ESO 2021

1. Introduction

Galactic-scale outflows are expected to play a major role in regulating star formation (Hopkins et al. 2012; Hirschmann et al. 2013; Chisholm et al. 2017), mediating the coevolution of galaxies and their central supermassive black holes (Fabian 2012; King & Pounds 2015), and setting the observed masses of nearby early-type galaxies. There are two main types of outflows: generated either by the active galactic nuclei (AGNs) in galaxy centers or through powerful winds driven by intense star formation. The galaxy baryonic mass distribution is affected by feedback at low mass (stellar feedback) and at high mass (AGN feedback). In either case, feedback can potentially cause the transfer of large quantities of matter and energy throughout the galaxy interstellar medium (ISM), leading to the enrichment of the ISM (Oppenheimer et al. 2010), and possibly to the suppression and quenching of star formation (Krumholz et al. 2017; Su et al. 2017).

Kiloparsec-scale feedback driven by star formation within the ISM is usually described by two physical phenomena: ionization and mechanically driven outflows. Radiation-pressured stellar winds originate from the intense radiation of young OB stars, which is strong enough to photo-ionize the ISM gas, mainly in the initial phase of the feedback process, while the outflow velocity builds up (Chevalier 1997; Veilleux et al. 2005; Girichidis et al. 2020). Mechanically driven outflows then dominate after enough mass is dragged by the feedback, resulting in high-velocity shocks (Sharp & Bland-Hawthorn 2010). Supernova events are usually related with the latter type of feedback, displaying characteristic shock-driven ionized gas spectra. Studies of star-formation-driven outflows using integral field spectroscopy (IFS) observations only observe the domination of one of these processes since they usually happen on smaller scales than the spatial resolution and are probably mixed in the resulting observed spectra (e.g., Ho et al. 2014; Rodríguez del Pino et al. 2019; D’Agostino et al. 2019).

Expected to dominate the feedback effect in the lower-mass end, star formation regulates itself by outflowing and heating the neighboring gas, which is more abundant in regions with higher star-formation rates (SFRs). The higher SFR, however, correlates with the mass outflow rate, creating a regulation cycle that prevents star-forming regions from growing in stellar mass. This relation is given by the equation Ṁ_out = η_SFR, where Ṁ_out is the mass outflow rate and η is the mass loading factor (Schaye et al. 2010; Vogelsberger et al. 2013; Somerville & Davé 2015).

From this relation, we expect star-formation-driven outflows to be more powerful at redshifts z ∼ 2, where the SFR presents its cosmological peak (Madau & Dickinson 2014). Indeed, both the SFR and outflow velocities seem to be higher in galaxies at higher redshifts than compared to nearby galaxies (Sugahara et al. 2017; Davies et al. 2019). Outflows in high-z galaxies, however, are more difficult to measure, usually demanding methodologies that result in loss of the spatial resolution, such as the stacking of large sections of galaxies (Herrera-Camus et al. 2021) or even the summing of the integrated spectra of multiple galaxies (Ginolfi et al. 2020).

In this sense, IFS studies of low-redshift galaxies present the advantage of providing more detailed information regarding the outflow component for even not-so-extreme star-formation activity (e.g., Rodríguez del Pino et al. 2019; Zaragoza-Cardiel et al. 2020). With the emergence of large IFS surveys such as CALIFA (Sánchez et al. 2012), SAMI (Croom et al. 2012), and MaNGA (Bundy et al. 2015), in addition to powerful instruments mounted on large light-collecting telescopes such as VLT/MUSE (Bacon et al. 2010), GTC/MEGARA (Carrasco et al. 2018), and Gemini/GMOS (Allington-Smith et al. 2002), statistical studies can now be performed in spatial resolutions on sub-kiloparsec scales. This allows us to investigate feedback properties, such as ionization mechanisms, outflow extent, and kinematics, within a large number of galaxies on local scales.

Aside from providing the benefit of higher resolution probing smaller physical scales, low-redshift studies also allow us to trace the bulk of the feedback phenomena in normal star-forming systems instead of those with extreme starburst activity seen at higher redshift. In fact, starburst galaxies represent only 1% of low-redshift star-forming galaxies, accounting for only 3−6% of the stellar production (Bergvall et al. 2016). Even when we move to higher redshifts, such as ∼2, starburst galaxies are only 2% of the galaxies forming stars and contribute to 10% of the SFR density (Rodighiero et al. 2011). This indicates that star-forming galaxies that experience feedback processes in their cosmic evolution are usually located within (or close to) the star-formation main sequence (Brinchmann et al. 2004; Peng et al. 2010).

Characterization of star-formation outflows is a difficult task due to several aspects. In addition to the SFR dependence with redshift, outflow velocities are also dependent on other parameters, such as stellar mass. Furthermore, works in the literature use different methodologies in order to trace the outflow components, such as non-parametric (e.g., Heckman et al. 2015), one-component (e.g., Du et al. 2016), and two-component (e.g., Rubin et al. 2014) measurements and the use of line ratio maps (e.g., Veilleux & Rupke 2002) in neutral (e.g., Cazzoli et al. 2014; Concas et al. 2019), molecular (e.g., Leaman et al. 2019; Salak et al. 2020), and ionized (e.g., Arribas et al. 2014; Roberts-Borsani et al. 2020) gas phases. Comparing observational results with one another and with theoretical predictions is not simple, and more information regarding all these aspects of the outflowing gas is needed to better understand the role of stellar feedback in galaxy evolution.

In this work we aim to probe the stellar feedback characteristics in a sample of low-redshift star-forming galaxies using spatially resolved MUSE data. The paper is structured as follows: In Sect. 2 we present our data, how we selected our sample, and our methodology. Section 3 presents the results we obtain by measuring the outflowing gas and the characterization of its parameters. We estimate the outflow properties and energetics in Sect. 4. In Sect. 5 we discuss our results in comparison with others from the literature using the estimated outflow energetics. Finally, we conclude in Sect. 6. Throughout this paper we use a cosmology with H₀ = 70.5 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.73, and Ω_M = 0.27 and a Chabrier (2003) initial mass function (IMF).

2. Data and analysis

In order to probe the properties of star formation and outflows, we selected a sample of 15 star-forming galaxies at 0.03 <  z <  0.2, which were originally part of the Valparaíso ALMA/APEX Line Emission Survey (VALES; Villanueva et al. 2017; Hughes et al. 2017) sample of 91 galaxies at 0.02 <  z <  0.35. The final sample was selected following several criteria, described below.

2.1. The VALES sample

The VALES sample has been obtained from the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010), which provides far-IR observations of extragalactic sky with the PACS and SPIRE cameras in the 100, 160, 250, 350, and 500 μm bands (Ibar et al. 2015). The first 67 galaxies of the VALES sample were selected from the three equatorial fields covered by H-ATLAS (totaling ∼160 deg²), satisfying the following criteria: (1) flux of S_160 μm >  150 mJy, close to where typical star-forming galaxies have their flux peaks; (2) sources with no S_160 μm >  160 mJy (3σ) detections within 2 arcmin; (3) having an unambiguous identification in the Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009), with RELIABILITY > 0.8 (Bourne et al. 2016; 4) Petrosian SDSS radii smaller than 15″; (5) high-quality spectra from the Galaxy and Mass Assembly survey (GAMA; Z_QUAL ≥ 3 Liske et al. 2015; Driver et al. 2016). ALMA CO(1−0) line observations (as part of the VALES collaboration) and GAMA far-UV (FUV) to far-IR (FIR)/submillimeter photometry are also available for these galaxies (Villanueva et al. 2017). A total of 24 additional starburst galaxies complemented the initial VALES sample observed with the Atacama Pathfinder Experiment APEX/SEPIA Band 5 targeting the CO(J = 2−1) emission line at z = 0.1−0.2 (Cheng et al. 2018), resulting in the final VALES sample of 91 galaxies.

2.2. GAMA data

GAMA multiwavelength data are ideal for the analysis of the spectral energy distributions (SEDs) of our sample. In this work we make use of 19 broadband photometric data sets from the ultraviolet to the far-IR included in the GAMA Panchromatic Data Release (Driver et al. 2016) from several facilities: FUV and near-UV (NUV) bands from GALaxy Evolution eXplorer (GALEX); u, g, r, i and z optical bands from SDSS; Y, J and H near-IR bands from the Visible and Infrared Telescope (VISTA); W1, W2, W3 and W4 mid-IR bands from the Wide-field Infrared Survey Explorer (WISE); and 100, 160, 250, 350 and 500 μm far-IR bands from H-ATLAS, resulting in SED photometry between 0.15 and 500 μm available for our sample.

2.3. MUSE data

Optical IFS data were obtained with the MUSE instrument (Multi Unit Spectroscopic Explorer; Bacon et al. 2010), mounted on the Very Large Telescope (VLT) for 18 galaxies of the original 91 from the VALES sample as part of the observing program 0100.B-0764(A) (PI Hughes). These galaxies were chosen because of the possibility of tracing the CO(1–0) emission line with ALMA at higher resolution. Observations were performed between November 2017 and February 2018 in dark time. The wide-field mode was used with a field-of-view of ∼1′×1′ and a 0.2″ spatial sampling. The wavelength range covered in each data cube is ∼4750−9350 Å with 1.25 Å spectral sampling and a line spread function (LSF) with a full width at half maximum of ∼2.5 Å around Hα (∼48.5 km s⁻¹). The spectral resolution ranges from 1750 in the blue domain to 3750 in the red (Bacon et al. 2017).

Each observation was split into five Detector Integration Time (DIT) exposures with offset 90° angles. We measured seeing-limited spatial resolutions of 0.58−0.94 arcsec. These represent physical scales of 0.4−3.3 kpc at the galaxies, with a median of 1.1 kpc. The per-pixel signal-to-noise ratio (S/N) at 5500 Å was found to be well above 5 in the inner disks of the galaxies. Peak S/N values range from 10−373 with a median of 41. The quality of the data ensures that we can extract reliable properties for our science goals.

The data reduction pipeline (v2.4.1 Weilbacher et al. 2016) was used to process the raw data into science-ready data cubes. This includes bias and dark-frame removal, flat-field corrections, wavelength and flux calibration. For the astrometric calibration and telluric line subtraction we used the default configuration files of the pipeline. Sky exposures near each target were used to correct for the sky continuum. Finally, the five DITs were aligned using in-field point sources and combined into a 3D data cube. We assessed the quality of the reduced data by comparing nuclear spectra with those from the GAMA database and found good correspondence in spectral shape and absolute flux.

2.4. SED modeling with CIGALE

To reproduce the SED profiles of our sample with GAMA data, we made use of the CIGALE code (Code Investigating GALaxy Emission¹; Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019). CIGALE computes the star-formation history (SFH) of a galaxy and combines it to single stellar population models, estimates the emission from ionized gas and applies flux attenuation due to dust. Dust emission based on energy balance in the mid- and far-IR is also taken into account. This is added to possible AGN and radio synchrotron emission to model the multiwavelength data given as an input to the code. As a result, CIGALE gives many physical parameters obtained from the fit, such as SFHs, stellar mass, dust luminosities and many others. One can inspect the quality of the fit through evaluation of the best χ² and reduced χ² as , where N is the number of bands fitted. Aside from the best fit obtained through χ² minimization, CIGALE also performs a Bayesian estimate of the physical properties through a likelihood-weighted mean, computed for all grid models, and a 1σ uncertainty represented by the peak width of the probability distribution function. The values used from the SED fits in this paper were obtained from the Bayesian approach, unless stated otherwise, as more reliable estimates should be obtained (see Noll et al. 2009, for a more detailed discussion). Mock galaxy catalogs can also be created to estimate the robustness of the fit, where CIGALE computes the flux densities from the fit added with instrumental noise, and then reanalyzes this mock data to obtain the physical parameters and compares with the results of the original fit. The main goal of the SED fitting in this work is to obtain SFHs and other spatially integrated parameters for our sample and use them to compare with possible outflowing scenarios in the sample. The top left panel of Fig. 1 displays an example of CIGALE fit for one galaxy of our sample.

Fig. 1. Example of spectral fittings for one galaxy (HATLASJ083832). Top left: CIGALE SED fit, where points represent the broadband magnitudes from the GAMA survey and the red, orange, blue, and green lines represent the dust spectra, stellar spectra, stellar spectra attenuated by the dust, and emission-line continuum spectra, respectively. Blue points represent the UV and optical bands, and the red points represent the IR bands. Top right: STARLIGHT stellar continuum fit, with the black and red lines representing observed and modeled spectra, respectively. Bottom left: single-Gaussian fit of the Hα+[N II] emission lines, where the points and the black line represent the observed spectrum, and the red line represents the modeled lines. A residual spectrum is shown in the bottom part of the panel. Bottom right: same as the previous panel, but showing the double-Gaussian fit, where both components are represented by blue lines. The Gaussian and STARLIGHT fits are from the same spaxel, located in the central region of the galaxy.

We used the following GAMA broadbands in the SED fit: GALEX FUV and NUV; SDSS u, g, r, i, and z; VISTA VIRCAM Y, J, H, and Ks; WISE1, 2, 3, and 4; and Herschel PACS100, PACS160, SPIRE 250, 250, and 500 μm. We also included the Hβ and Hα integrated fluxes to the fit, obtained from the single component fit with MUSE data. As we are more interested in the SFH parameters of the fit, we decided to first fit the dust contribution of the SED distribution, which is a source of uncertainty in the fit, since the IR bands present considerably larger uncertainties in comparison to the other bands. We then fitted the IR bands (from WISE3 to longer wavelengths) isolated just with the Draine et al. (2014) models to fit the dust contribution with four free parameters.

We then used the resulting dust luminosities to fit the other bands, including the stellar (Bruzual & Charlot 2003) and nebular contributions, the dust attenuation (modified attenuation law from Charlot & Fall 2000, with two free parameters), and the SFH function. We assumed a delayed exponential SFH with two constant burst/quenching terms, as described by the equation:

(1)

where r₁ and r₂ are the constant burst/quenching terms for two recent time periods t₁ and t₂. Assuming we can constrain the recent SFR, since we include Hβ and Hα emission lines together with the UV bands in the fitting procedure, we fixed t₁ = 100 Myr and t₂ = 10 Myr ago (in lookback time). Aside from r₁ and r₂, we also left τ as a free parameter. A summary of the parameters used in the fit is described in Table 1 and the resulting SED fit for each galaxy is shown in Appendix A. We also show the reduced chi-square and residuals of the UV and optical part of the SED fit for the VALES sample fit in comparison with the mock models in Fig. 2. As shown, the best-fit results for the VALES sample are comparable with the mock models results, indicating quality sample fits.

Fig. 2. Quality of the SED modeling of the lower-wavelength photometric broadbands (from FUV to WISE2), with the addition of the Hα and Hβ integrated fluxes from MUSE data. Total reduced chi-square and median residuals (in percentage of the observed magnitudes from GAMA and MUSE), for both the estimated value for the VALES sample and the exact value for the mock sample, are displayed in the left and right panels, respectively.

2.5. Stellar population synthesis with STARLIGHT

In order to obtain ionized emission-line gas properties, we must first remove the stellar continuum and absorption contribution to the IFS spectra. This was performed using the STARLIGHT spectral synthesis code (Cid Fernandes et al. 2005). STARLIGHT fits the observed spectrum with a combination of single stellar populations, with different ages and metallicities, assuming a single dust screen with R_V = 3.1, following Cardelli et al. (1989) dust law. Before we applied the code to the MUSE data cubes, we corrected the spectra for redshift and Galactic extinction using the PYCASSO software (Cid Fernandes et al. 2013; de Amorim et al. 2017). We used a base of 45 simple stellar-population models from Bruzual & Charlot (2003) with a Chabrier (2003) IMF and STELIB stellar library (Le Borgne et al. 2003) and “Padova (1994)” evolutionary tracks (Alongi et al. 1993; Bressan et al. 1993; Fagotto et al. 1994a,b,c; Girardi et al. 1996). The base of templates present 15 age bins between 1 Myr and 13 Gyr and 3 metallicity bins, Z = 0.004, 0.02 and 0.05. The top right panel of Fig. 1 shows a STARLIGHT fit example for one spaxel of one of the galaxies in the sample.

2.6. Emission-line measurements

Emission-line parameters such as integrated flux, centroid velocity and velocity dispersion were obtained by fitting Gaussian profiles to the stellar-continuum free residual spectra resulted from the STARLIGHT fittings to the MUSE data cubes. A fitting code developed by our group², which makes use of the LMFIT package (Newville et al. 2014), a nonlinear least-square minimization solver based on the Levenberg-Marquardt algorithm from SCIPY (Virtanen et al. 2020), was used to fit the emission lines. The main emission lines (Hβ, [O III]λλ4959,5007, [O I]λλ6300,6366, Hα, [N II]λλ6548,6584 and [S II]λλ6717,6731) were fitted. We constrain the line central wavelength and width to be the same for the same ionic species ([N II] and Hα lines would also share these parameters) and the [O III]λ5007/4959, [O I]λ6300/6366, and [N II]λ6584/6548 ratios are fixed to the values predicted by atomic physics. Before we applied the emission-line fittings, we employed the Voronoi tessellation technique (Cappellari & Copin 2003) with a target S/N ratio of 10 at 6400–6500 Å rest-frame spectral region.

We fitted the data cubes with two strategies: (i) one Gaussian profile for each emission line; and (ii) initially one Gaussian for each emission line but then selecting spaxels where a second component would be necessary. In the latter case, a second, broader, component would be added, and the selected spaxels would be refitted. The criterion selected to flag a spaxel for the double component fit is for two consecutive spectral pixels, within the Gaussian profiles (not within the fitted continuum), to be above or below the 5σ range, where σ is the noise surrounding the fitted emission line. Finally, we created map masks for “good” spatial pixels using the following criteria, which were applied to the single component Hα fit: (i) flux uncertainties should be below 20% of the measured integrated flux; (ii) peak velocity and velocity dispersion uncertainties should be below 10 km s⁻¹. Uncertainties were obtained through LMFIT fitting, which explores the parameter space, estimating confidence intervals for variable parameters. Regions that would not fulfill the above criteria were masked out. We also apply the line fitting accounting for the MUSE LSF such as , using the instrumental dispersion values described in Sect. 2.3. No corrections regarding the point-spread function (PSF) was taken into account when fitting the emission lines. Figure 1 shows an example of the two emission-line fitting strategies.

We inspected the possibility of the second component arising due to beam smearing due to the spatial resolution limitations. We follow the methodology applied by Gallagher et al. (2019), comparing the second component velocity dispersion with the stellar continuum velocity dispersion we obtain with STARLIGHT. If beam smearing effect would to be considerable, the broadening measured by the second component would also be present in the stellar velocity dispersion. We find that σ_2nd is usually considerably higher (∼2×) than σ_* for the Voronoi regions we trace the second component, with median values of σ_2nd = 110 km s⁻¹ and σ_* = 47 km s⁻¹, and σ_2nd exceeding σ_* ∼92% the 1σ level. This indicate that beam smearing has no effect in the broad results of the paper.

We also note that, for the SED fitting described in Sect. 2.4, we used the integrated Hα and Hβ fluxes summing the first and second components from the double component emission-line fit.

2.7. Final sample

From the 18 galaxies observed with MUSE, we removed two galaxies (HATLASJ090750 and HATLASJ084630) due to them being interacting galaxies and one galaxy (HATLASJ085835) due to it being an active galaxy. These galaxies were removed due to our emission-line fitting strategy. The first two clearly present more than one narrow component, probably because of a superposition of different galaxies along the line of sight³. The third galaxy presents clear AGN type I broad hydrogen profiles, thus affecting the measurement of star-forming emission. We do not rule out any further galaxies from being AGNs: In fact, there are a couple of galaxies that suggest the presence of AGNs through diagnostic diagrams (HATLASJ085828 and HATLASJ090532). Regions related to AGN ionization in the BPT diagram (Baldwin et al. 1981) were removed from our analysis.

The final sample consists of 15 galaxies, distributed along the M_*-SFR diagram as shown in Fig. 3, spanning the stellar mass interval of 10.2 < log(M_*/M_⊙) < 11.3 (median log (M_*/M_⊙) = 10.6) and 0.2 < log SFR (M_⊙ yr⁻¹) < 1.8 (median log SFR (M_⊙ yr⁻¹) = 0.6). We note that some galaxies in our sample present SFRs of up to 10 times the star-formation main sequence. The sample is, however, generally closer to the main sequence than to starburst values. Further details about each galaxy of the sample are described in Appendix A.

Fig. 3. Galaxy sample distribution in the M*-SFR diagram (blue dots). The median of the distribution is represented by the yellow star. Values were obtained from the CIGALE fit of the sample. Dashed blue and red lines represent the star-forming main sequence and retired-galaxy relations obtained by Cano-Díaz et al. (2019). Dotted lines indicate SFR offsets by the displayed multipliers.

In comparison, the values reported from Villanueva et al. (2017) for the first 67 VALES galaxies are 0.1 < log SFR (M_⊙ yr⁻¹) < 1.9 (median log SFR (M_⊙ yr⁻¹) = 1.2) and 10.7 < log(M_*/M_⊙) < 11.3 (median log (M_*/M_⊙) = 10.7). These parameters were obtained using MAGPHYS SED fitting code (da Cunha et al. 2008). The sample we are probing, therefore, has similar values with respect to the VALES mother sample.

2.8. Derived spatially resolved maps

Several distribution maps are derived from the methodology described above to aid our analysis. We estimate the SFR surface density (Σ_SFR) using the dust-corrected Hα flux (extinction law from Cardelli et al. 1989) measured summing the first and the second component from the double component emission-line fit. We used the relation log SFR = log L_Hα + k (Kennicutt & Evans 2012) with the k coefficients from Boquien et al. (2021). This coefficient, which we also mapped its spatial distribution, depends on the assumed IMF (Chabrier) and the gas metallicity. Metallicity maps were constructed using the Marino et al. (2013) relation using the O3N2 calibrator. Specific SFRs (sSFR) of Σ_SFR/Σ_M* were also obtained using the stellar-mass surface density that resulted from the STARLIGHT fit. We constructed Σ_SFR, Σ_M*, and sSFR maps for the galaxy sample, which are shown in the individual figures in Appendix A. A summary of the galaxy sample physical properties is described in Table 2.

3. Results

3.1. Single-component fit

We begin investigating signatures of outflowing gas by looking into the single-component emission-line fit of the spatially resolved MUSE data. Assuming that the presence of a secondary component would drive the single-component fit to present a somewhat higher velocity dispersion and a displaced peak velocity in relation to what should be expected if this second component were not present (as we can observe in the example shown in the bottom panels of Fig. 1), we can expect the kinematic properties of the single-component fit to indicate the presence of the secondary component by somehow measuring this velocity shift and the increase in velocity dispersion.

To probe this scenario we fitted a disk model to the gas velocity field. Following Bertola et al. (1991), we assumed a spherical potential, with the gas following only circular orbits. The observed radial velocity at the position (R,Ψ) is given by the equation

(2)

where v_sys is the systemic velocity (which is close to zero since we have redshift-corrected the spectra), A is the centroid velocity amplitude, R and Ψ are the radial and angular coordinates of a given pixel in the plane of the sky, Ψ₀ is the position angle of the line of nodes, c₀ is a concentration parameter (constraining the radius at which the centroid velocity reaches 70 per cent of the amplitude A) and θ is the disk inclination (θ = 0° for a face-on disk). Finally, the parameter p measures the slope of the rotation curve after reaching the maximum amplitude. This parameter was fixed to p = 1, which corresponds to an asymptotically flat rotation curve at large radii. The best fit was obtained also using the LMFIT package with a Levenberg-Marquardt least-squares minimization. We used the optical continuum image of each galaxy as weight in order to better reproduce the observed velocity maps. The resulting velocity model maps for each galaxy of the sample are shown in Appendix A, as well as the residual velocities (v_res = v_obs − v_model).

Our main interest in this model is to retrieve the residual velocities. These are representative of the centroid velocity deviation due to the presence of a second component in the line profile. Figure 4 shows the relation between the velocity dispersion and the Σ_SFR with the residual velocities for the entire sample. Although some scatter is observed, no correlation between the two parameters is observed with the residual velocities, indicating that the hypothesis raised in the beginning of this section does not hold. Figures of these relations for individual galaxies, shown in Appendix A, show that this scenario could be true for some galaxies of our sample. This could be the case of HATLASJ083832, HATLASJ085356, and HATLASJ090949, for example, where higher residual velocities are observed for higher Σ_SFR and velocity dispersion. The Σ_SFR and σ values where v_res increases, however, do not seem to be the same between these galaxies. This, together with the fact that we may not have the presence of outflows in a good part of the sample, may “wash out” the outflow signature in Fig. 4. We explore this scenario further in Sect. 3.3, where we compare these slopes with the SFR obtained through the SED fitting.

Fig. 4. Relation between the velocity dispersion σ (left) and SFR density ΣSFR (right) with the residual velocities obtained from the rotation model fitted to the Hα velocity fields. Each point represents a Voronoi region. The top panels show the distribution for the whole sample, while the middle and bottom panels are separated into redshifted and blueshifted residuals, respectively. Square symbols represent median values, with the dashed lines representing 25 and 75 percentiles of the distribution. No clear correlation between the two parameters is observed with vres.

3.2. Double-component fit

For regions where a single component did not result in a satisfactory fit, we added a second component, following the guidelines mentioned in Sect. 2.6, considering it to be related with outflowing gas due to star-forming activity. Because it is difficult to obtain strong constraints on its parameters, we decided to constrain the second component to (i) have lower integrated flux than the first component, so F_2nd/F_1st <  1; (ii) its velocity dispersion needs to be higher than 1.2 and lower than 3.0 times the first-component velocity dispersion. These values are empirical, and were chosen after testing several velocity dispersion conditions; (iii) its centroid velocity should not differ by more than 200 km s⁻¹ (in absolute values) than the one fitted for the first component. We also fixed the Hα/Hβ and [N II]/Hα ratios between the first and second components, so both ratios are the same for both components. The Voronoi regions where we fitted a second component represent a small portion of some galaxies in the sample, as can be seen in the flux maps of the second component in the figures shown in Appendix A, for each galaxy in our sample. In fact, we detect the second component in 13 of the 15 galaxies.

Figure 5 shows the relation between the velocity dispersion of the second component σ_2nd and the integrated Hα flux ratio F_2nd/F_1st. The left panel includes regions where the emission-line fit results in some parameter boundary (e.g., F_2nd/F_1st = 1), which we call “boundary regions”. The right panel is the same, removing such regions. From now on, we only consider non-boundary regions, as we consider these to be satisfactorily constrained in our measurement methodology. We can observe in Fig. 5 that the bulk of the regions present a 80 km s⁻¹  <  σ_2nd <  150 km s⁻¹ and a F_2nd/F_1st ∼ 0.2. The distribution is scattered for lower velocity dispersion, but narrows for values of σ_2nd >  130 km s⁻¹, where the F_2nd contribution is lower than half of the F_1st. Few points reach values σ_2nd >  200 km s⁻¹, with 0.2 <  F_2nd/F_1st <  0.4. Although there is a suggestion that the F_2nd/F_1st decreases with σ_2nd in the σ_2nd >  200 km s⁻¹ region, the region is barely populated, and we would need more points to make conclusions for this region.

Fig. 5. Relation between the velocity dispersion of the second component and the flux ratio between the second and the first components, from the double-component fit. Left panel: all the Voronoi regions where a double-component fit was performed. These include what we call the “boundary regions” (explained in the text), and the number of such regions, discriminated by criteria, is shown in the top right. Right panel: distribution without boundary regions.

The left panel of Fig. 6 shows the relation between Σ_SFR and F_2nd/F_1st flux ratio. The F_2nd/F_1st ratio displays scattered distribution up to log Σ_SFR ∼ 2.0 M_⊙ yr⁻¹ kpc⁻², but concentrated in values F_2nd/F_1st <  0.5. F_2nd/F_1st generally decreases with increasing Σ_SFR in this interval. This could be due to that, although the S/N ratio is high overall for all the points in this plot⁴, the S/N ratio is probably lower for low Σ_SFR values, resulting not only in the scatter observed in the F_2nd/F_1st ratio but also restricting us to trace small F_2nd contributions, and thus not being able to populate low F_2nd/F_1st ratio for low Σ_SFR values. This is indicated by the apparent concentration of larger Voronoi regions in the left lower corner of the plot, coming from regions where more binning was needed to reach a satisfactory S/N ratio. The trend shifts with increasing F_2nd/F_1st ratio for log Σ_SFR >  2.0   M_⊙ yr⁻¹ kpc⁻², where the flux ratio is much less scattered, as shown by the 1-σ interval lines, indicating a higher flux contribution from the second component in high-SFR regions.

Fig. 6. Relations between SFR surface density and sSFRs and the second-over-first-component flux ratio. Black points represent Voronoi regions, where the point size is proportional to the inverse of the Voronoi covered area, so small points represent regions where more spaxels had to be binned to reach a satisfactory S/N. Strong blue lines represent a running median with a 20-point window, with weak blue lines representing the 1σ of the distribution.

The right panel of Fig. 6 shows the relation between sSFR and F_2nd/F_1st flux ratio, where we observe again a decrease in sSFR and F_2nd/F_1st with a shift in very high values (log sSFR >  0.3 Gyr⁻¹), although in this case there are fewer points in this trend shift. Again, this could be just dominated by the lower S/N ratio. In this case, there is no clear shift as shown in the left panel (or at least not as clear), where the flux ratio seems to increase after a certain SFR value. Visually inspecting the regions where we trace the second component, we notice these are usually located in galaxy centers, where stellar mass is increased, which could explain the higher F_2nd contribution in lower sSFR values. However, since it is hard to disentangle the dependence factor in this relation, further analysis would be needed to address this relation.

In Fig. 7 we display the relation between the Σ_SFR and the second-component velocity dispersion, σ_2nd. Although the σ_2nd are not high, concentrating between 80 <  σ_2nd <  120 km s⁻¹, we can observe a correlation between the two. This is expected if we consider that the second component is indeed tracing an outflowing gas due to increased star-formation activity. The low σ values, however, seem to indicate that these outflows are not strong and of low velocities.

Fig. 7. Second-component velocity dispersion as a function of the SFR surface density. Voronoi regions are represented by black points, with point size being anticorrelated with the Voronoi region area. The black line and blue shaded region show the best-fit relation with a 3σ uncertainty, and the relation is displayed in the top-left corner. Pearson coefficients are also displayed, indicating a significant positive correlation.

We note that, due to our measurement strategy of fitting a second component in regions where the residuals above the 5σ threshold, we expect these regions to present higher signal-to-noise ratios. Thus the fact that we are probing low velocity (and luminosity, as we will see in the following sections) outflows are a result of the deeper observations in this work, when compared to other studies using galaxy surveys (e.g., MaNGA).

3.3. SED fit

After inspecting the relations in the spatially resolved MUSE data, now we take a look at the results obtained with the integrated broadband SED photometry fitted with CIGALE. In order to probe whether the outflow activity affects the recent bursting or quenching of star formation, we compare the ratio of the SFRs at 10 and 100 Myr obtained from the SED fitting to the physical properties associated with outflowing ionized gas. In Fig. 8 we compare this ratio with the slopes we observe for each galaxy between Σ_SFR and the single-component σ with the residual velocities (Appendix A). As indicated in Sect. 3.1, these properties do not seem to show outflow signature for the entire sample, but could for individual galaxies. However, this does not seem to be the case. We do not observe any correlation between the SFR_10 Myr/SFR_100 Myr ratio and the Σ_SFR and v_res slope, even though some of the galaxies present steep slopes between v_res and Σ_SFR (e.g., HATLASJ083832). Although a weak anticorrelation is observed between the σ × v_res slope and the SFR ratio, the significance is too small to take this result into consideration.

Fig. 8. Relation between the slopes of log ΣSFR (left panel) and the one-component fit σ (right) versus residual velocities and the SFR10 Myr/SFR100 Myr ratio obtained with the SED fit. The slopes are obtained through the fits shown for each individual galaxy in Appendix A. Pearson coefficients are given in the top-right corners; they show no correlation in the first relation and only weak anticorrelation with little significance in the second. The cyan line shows the SFR10 Myr/SFR100 Myr = 1 ratio.

This result, in addition to the lack of evidence of outflowing gas in Fig. 4 as discussed in Sect. 3.1, indicate that, if gas outflows are present, the simple disk kinematics model we are using here is not detailed enough to account for the non-circular motions these outflows could generate. This is also a clear indication that these outflows would not be strong, since no major kinematics deviation from the galaxies rotating disk is observed.

Figure 9 displays the dependence of the SFR_10 Myr/SFR_100 Myr ratio on the integrated Hα luminosity ratio between the second and the first components, for the case of the two-component fit. We observe a moderate anticorrelation with a small significance, but higher than the relation observed in the right panel of Fig. 8. The significance is still not good enough to consider this anticorrelation to be true, and a larger sample would be needed to confirm this trend. If this is indeed a true anticorrelation, we could be probing the outflow impact timescales in SFR, with the feedback being effective in inhibiting star formation between 10 and 100 Myr. However, even though we could be tracing the outflow signature in the single- and double-component emission-line measurements for some galaxies, the results do not show a clear indication of impact in the recent SFH between 10 and 100 Myr at global scales.

Fig. 9. Dependence of the SFR10 Myr/SFR100 Myr ratio with the spatially integrated Hα luminosity ratio between the second and the first components. Pearson coefficients are given in the top-right corner; they show a moderate anticorrelation with a small significance.

We note that for five of the galaxies in our sample, companions can be observed within the MUSE field of view (FoV), and these have different GAMA IDs. This means that the GAMA data used to do the SED modeling do not integrate the entire MUSE FoV emission, which could lead to discrepancies when comparing the emission-line and SED results, the main topic of this section. Thus, for Figs. 8 and 9 we estimated the slopes and the luminosity ratios masking out the companion galaxies. These companion galaxies, however, have negligible contribution to the v_res slopes and L(Hα)_2nd/1st ratio, and the masking does not change the results obtained here. These galaxies are HATLASJ084217, HATLASJ085406, HATLASJ090949, HATLASJ114244, and HATLASJ142128.

3.4. Gas excitation

In order to investigate the gas excitation of the outflowing gas, we analyze the line ratios in the Voronoi regions where we measure a second component in the emission-line spectra, using MUSE resolved data. Figure 10 shows the distribution of these regions in the BPT diagrams (Baldwin et al. 1981) and their dependence on the second-component velocity dispersion σ_2nd. We also display the relation between the [N II]/Hα, [O III]/Hβ and [S II]/Hα line ratios with σ_2nd. It is important to note that since we do not let line ratios vary between the two components in the emission-line measurements, we decide to use the line ratio values from the single-component fit, which are better constrained. We observe that the regions where the outflow seems to be more intense (σ_2nd ≳ 200 km s⁻¹) display low [O III]/Hβ ratios and are close to the composite region in the BPT diagram. These values fall quite far from what should be expected by shock ionization models shown by the line grids in the figure, obtained using the MAPPINGS V code and the 3MdB database (Sutherland & Dopita 2017; Alarie & Morisset 2019). In fact, the regions with the highest σ_2nd appear to be the most distant from the regions covered by the shock models.

Fig. 10. BPT diagrams (top panels) and the dependence of the second-component velocity dispersion σ2nd on the [N II]/Hα, [O III]/Hβ, and [S II]/Hα line ratios of the Voronoi regions where we do measure a second component related to outflowing gas (bottom panels). The BPT diagrams are colored with its σ2nd. Dashed lines display the theoretical (Kewley et al. 2001) and empirical (Kauffmann et al. 2003) regions, separating pure star formation and star formation combined with AGN ionization, respectively. Red, black, and green grids display MAPPINGS V shock models for LMC (Z = 0.0071), Dopita_2005 (Z = 0.0065), and solar metallicities (Z = 0.0183), respectively. Grids are for electron density ne = 1 cm−3, with magnetic fields 0.0001 > B < 10 μG and shock velocities 100 <  vshock <  400 km s−1 (shock velocities usually vary along the x axis of the BPT diagram). Some shock velocity values are indicated, as is the general direction of increasing magnetic field intensity. Blue lines in the bottom panels represent the running median with the 1σ distribution of the Voronoi regions. Dashed green lines represent solar metallicity shock models with B = 1 μG and 0.01 <  ne <  1000 cm−3, identified by the labels. Y-axis values for the dashed green lines are vshock.

The dependence of σ_2nd on the [N II]/Hα ratio would indicate that shocks are driving the more turbulent motions. However, this is not in agreement with the [S II]/Hα, a ratio that is more sensitive to shocks, which shows a decrease in regions displaying higher σ_2nd, where it should also increase. This is illustrated by the dashed lines, which represent how the shock models vary with shock velocities in relation to these line ratios, for several electron densities. The observed [O III]/Hβ shows a similar behavior in comparison to the shock models (it decreases with increasing shock velocities), but with a considerable shift in line ratio values (observed [O III]/Hβ is ∼2 lower than what expected from the models), and this difference causes the observed values to not populate the parameter space of the shock models.

It is important to note that in this comparison we used simple shock models, without the precursor component. Adding this contribution further increases the difference observed in the [O III]/Hβ ratio.

We also compared the observed line ratios with photoionization models. We obtained emission-line intensities from the 3MdB database for CLOUDY models (Ferland et al. 2013), which is described in detail in Morisset et al. (2015). In summary, we used the “CALIFA” table for H II regions in the 3MdB database, which consists in a CLOUDY model using the simple stellar populations templates used by Cid Fernandes et al. (2013) in the analysis of CALIFA observations. These templates comprise four bins of metallicity (Z= 0.0037, 0.0076, 0.019 and 0.0315) and 39 different stellar ages (between 10⁶ and 1.4 × 10¹⁰ years). The templates are then used to compute the SED passed to CLOUDY to run the photoionization models. Aside from the stellar age and metallicity, four other parameters are considered in the 3MdB database: the nebula geometry (“form factor” f_S = 0.03 and 3.0 for a full sphere and a hollow spherical bubble, respectively); the nitrogen over oxygen abundance ratio (log N/O = − 0.5, −0.25, 0, 0.25 and 0.5); the ionization parameter (log(U) from −4 to −1.5 in steps of 0.25); and the Hβ luminosity fraction of the nebula (Hβ fraction between ∼20% to ∼100% in steps of ∼20%).

For this comparison we used a full sphere geometry (f_S = 0.03), a total nebula Hβ luminosity (Hβ fraction ∼100%) and a log N/O = 0.25. Grids were created using the values of log(U) between −4.0 and −2.0 (in steps of 0.5) and stellar ages of 1, 3, 4, 5.6, and 8.9 × 10⁶ years. Figure 11 shows the BPT diagrams with the data values shown in Fig. 10 now compared with the photoionization models described above, with metallicity values of Z = 0.0076 (upper panels) and Z = 0.019 (bottom panels).

Fig. 11. BPT diagrams as shown in Fig. 10, here in comparison with the photoionization grids obtained from the 3MdB database using CLOUDY models. The top and bottom panels show the model results for Z = 0.0076 and Z = 0.019, respectively. For comparison, the sample median metallicity is Z = 0.0086. Blue lines represent different stellar ages of 1, 3, 4, 5.6, and 8.9 × 106 years as ionization sources and are roughly oriented from top right to bottom left as age increases, although not linearly. Red lines represent variation in the photoionization parameter, ranging from log(U) between −4.0 and −2.0 (in steps of 0.5), roughly oriented from top left to bottom right as U decreases. Some stellar age values are displayed in the top-left panel, along with the increasing U orientation.

The observed line ratios match more closely the grid models for the Z = 0.0076 than the ones with Z = 0.019. This agrees with the measured median metallicity of the sample, Z = 0.0086, being closer to the lower-metallicity models. In general, the observed line ratios are better reproduced by the photoionization models than the shock models displayed in Fig. 10. The regions presenting high σ_2nd, observed in the bottom-right corner of the H II classification in the [N II]/Hα BPT diagram, are closely approached by the photoionization from very young ages (∼1 Myr) and low photoionization parameters (log(U)  ∼ −4.0). Slightly older populations (∼5 Myr) with also slightly higher photoionization parameters (log(U)  ∼ −3.5) appear to better reproduce the outflow ratios observed in the [S II]/Hα and [O I]/Hα BPT diagrams, which present more unfolded model grids and thus could better differentiate the photoionization parameters.

It is worth noting that stellar ages older than 5 Myr do not produce a strong [O III] emission line due to its high ionization energy level, thus the grids fall to very low [O III]/Hβ ratios (below the limits shown here). [O III]/Hβ would increase again for older populations (t >  0.5 Gyr), due to hot evolved stars, which are not representative of our sample and thus not shown in Fig. 11.

To summarize, we conclude that the observed line ratios related to the outflowing regions (regions where we measure an emission-line second component) are better represented by photoionization models than compared to shock models, with a subsolar metallicity. We further discuss how this compare to other results from the literature in Sect. 5.2.

4. Outflow parameters

The second component fitted to the emission-line profiles was detected in 13 of the 15 galaxies in our sample, as displayed in the individual maps for each galaxy in Appendix A. Assuming that this second component is related to outflowing gas in the star-forming regions in our galaxy sample, we can use this component to estimate the outflow parameters.

4.1. Outflow velocity

We first derive the velocity shift between the two components and the outflow velocity, considering also the velocity dispersion of the second component. The velocity shift Δv = v_1st − v_2nd as a function of the SFR surface density is shown in the left panel of Fig. 12. Although we allowed for the second component to assume both positive and negative velocities in relation to the first component, we observe that the outflow velocities are mainly toward blueshifted values, which is expected, since the outflow redshifted component is generally most susceptible to the host-galaxy extinction. The values are concentrated between 0 and −100 km s⁻¹, with a large scatter in the lower Σ_SFR values, reaching ∼ − 40 km s⁻¹ around Σ_SFR ∼ 1.0 M_⊙ yr⁻¹ kpc⁻², where the distribution concentrates, and remaining almost constant toward higher Σ_SFR values.

Fig. 12. Velocity shift between the two emission-line components and outflow velocity as a function of SFR surface density for the Voronoi regions, shown as black points. Blue lines in the left panel display the running median with a 1σ distribution. The black line and blue shaded region in the right panel represent the best-fit relation (displayed in the top-left corner, along with Pearson coefficients) with a 3σ uncertainty.

We adopt the Genzel et al. (2011) definition to estimate the maximum outflow velocity, with v_out ∼ |Δv − 2σ_2nd|. The right panel of Fig. 12 shows its dependence on Σ_SFR, where a correlation can be observed. Outflow velocities are usually low, mainly ranging between 200 to 300 km s⁻¹, with v_out ∼ 250 km s⁻¹ around Σ_SFR ∼ 1.0 M_⊙ yr⁻¹ kpc⁻², the peak of the distribution. A branch of higher outflow velocities (400 <  v_out <  550 km s⁻¹) appears in the higher Σ_SFR end, indicating the presence of stronger outflows in the sample, but this region is weakly populated. With a 3σ uncertainty, the estimated correlation is described as

(3)

4.2. Outflow extent and electron density

The estimation of the outflow extent, or outflow radius, and the outflow geometry in general, is perhaps one of the parameters leading to the higher uncertainties in outflow characterization. Here we have assumed that we do not spatially resolve individual H II regions⁵, and thus we also do not resolve the outflow spatially. However, we do detect the emission-line second component in regions smaller than the spatial resolution.

In order to deal with this issue, we only estimated the total outflow extent within each galaxy in the sample. This was performed by summing up the spaxels area in kpc², and assuming a projected circular geometry. Thus the resulting outflow projected area can be represented by , where R_out is the radius of the entire outflow area in the galaxy. The estimated R_out values are displayed in the top-left panel of Fig. 13, as well as its relation with the second-component Hα luminosity. A weak positive relation is observed with a large scatter. Outflow extents vary between 0.5 and 20 kpc, with an average of 6.8 kpc and a 1σ distribution of 5.4 kpc.

Fig. 13. Outflow properties of the analyzed sample. All panels are in function of the outflow SFR surface density, except the top left. Top left: outflow extent (definition described in the text) as a function of the second-component Hα luminosity. Top right: electron density. Middle left: mass outflow rate. Middle right: mass loading factor. Bottom left: ratio between the outflow kinetic power and the total bolometric luminosity. Bottom right: ratio between the momentum outflow rate and the SFR photon momentum. We display the Pearson coefficients for five of the displayed relations, and the best-fit model with a 3σ uncertainty for the panels is shown as a black line with blue shaded regions.

We also estimated the outflow electron density, as it is necessary to further estimate the outflow energetics. The [S II] emission lines are quite weak within the Voronoi regions, resulting in unconstrained estimations of the electron density. We thus estimated the electron density by integrating the fluxes of both lines in the regions where we detect the outflowing component, and used the ratio from these integrated fluxes in the calculation. The electron densities were obtained by using the PYNEB package (Luridiana et al. 2015) assuming a nebular temperature of 10 000 K, and are shown as a function of the SFR surface density in the top-right panel of Fig. 13. Density values vary between 10 and 730 cm⁻³, with a mean value of cm⁻³.

It is worth noting that the SFR surface density values displayed in Fig. 13 are obtained within the integrated outflowing regions in each galaxy, divided by the outflowing area (calculated as explained above). Thus, these values should not be directly compared to the SFR surface density values displayed in the previously shown figures, which take the SFR within the Voronoi regions into account.

4.3. Mass outflow rate and loading factor

After estimating outflow properties such as velocity, density, and extent, we can now derive the mass outflow rate and its relation with the SFR. Following the assumption that the outflow does not vary radially and that the second component is indeed related to the outflowing gas (either by the expanding shocked gas or by photoionization), we can derive the mass outflow rate as (Newman et al. 2012):

(4)

where m_H is the hydrogen atomic mass, γ_Hα is the Hα emissivity at gas temperature of T_e = 10⁴ K (γ_Hα = 3.56 × 10⁻²⁵ erg cm³ s⁻¹), n_e is the outflow electron density, v_out and R_out are the outflow velocity and radial extent, respectively, and the L(Hα)_2nd is the Hα luminosity of the second component, related to the outflowing gas.

We used the R_out and n_e values described in the previous section, and displayed in the top panels of Fig. 13. For v_out, we used a median value for each galaxy and the L(Hα)_2nd values used in this calculation are also shown in the top-left panel of Fig. 13. The relation between the derived mass outflow rate and the SFR surface density is shown in the middle-left panel of Fig. 13, and displays a strong correlation (r = 0.75) with a p value indicating less than 0.5% chance of the relation originating at random. Mass outflow rate values vary between 1.2 × 10⁻³  M_⊙ yr⁻¹ and 1.9 M_⊙ yr⁻¹, with a median of 0.1 M_⊙ yr⁻¹. The best-fit relation obtained is given by

(5)

considering a 3σ uncertainty.

The middle-right panel of Fig. 13 shows the relation between the SFR surface density and the mass loading factor, given by η = Ṁ_out/SFR. Although there seems to be a positive correlation, there is a considerable scatter (p-value = 0.22) indicating that a larger galaxy sample would be necessary to confirm this trend. We further discuss these results in Sect. 5.3.

4.4. Outflow energetics

Finally, we use the derived outflow properties to obtain the outflow kinetic power and momentum rate. The outflow kinetic power is given by , and its dependence with Σ SFR is shown in the bottom-left panel of Fig. 13, displayed as the ratio with the bolometric luminosity assuming that L_bol ∼ 10¹⁰ × SFR L_⊙ (Kennicutt 1998). The outflow momentum rate is given by , and is displayed in the bottom-right panel of Fig. 13 as a ratio to the photon momentum rate derived from the bolometric luminosity.

As expected from the relation between the Ṁ_out and Σ_SFR both and ṗ_out also relate with SFR surface density, although with a weaker correlation and an increased scatter. The best-fit relations are

(6)

(7)

where c is the light speed.

Outflow kinetic power accounts for very little, with values varying between 7.2 × 10⁻⁶% and 0.02% of the bolometric luminosity, with a median value of only 5.2 × 10⁻⁴%. The outflow, however, represents a median value of 1.3% of the momentum driven by the bolometric luminosity alone (varying between 0.02% and 30%), indicating that although the outflow does not contribute much to the energy deposition in the ISM, its contribution is not as negligible to the gas kinematics.

5. Discussion

5.1. Feedback impact on SFR

The results obtained with the SED fitting reported in Sect. 3.3 indicate no clear feedback impact in the recent SFH probed by our sample. We observe in Figs. 8 and 9 that the ratio between the SFRs at 10 Myr and at 100 Myr in entire galaxies does not change significantly either with single- or with double-component properties. These results indicate that even though we can trace outflowing gas signatures, the feedback is not strong enough to impact the recent star-forming activity.

We explore these results further by investigating not only the relation between the Hα luminosity ratio and the SFR ratio, but with SFRs at 10 and 100 Myrs and their difference, displayed in Fig. 14. We can observe that the outflow has increased luminosity contribution in higher-SFR galaxies in both timescales. The bottom panel of Fig. 14 also indicates a decrease in SFR between 100 and 10 Myr (up to ∼100 M_⊙ yr⁻¹) for some galaxies presenting higher outflow contribution. This relation exhibits a much higher significance (p value = 0.02) when compared to the other relations shown in Figs. 8 and 9 and follows the trend suggested in Fig. 9, where the SFR at 10 Myr is decreased in relation to the SFR at 100 Myr. However, as indicated in Fig. 14 (top panel), the SFR decrease at 10 Myr is minor, not significant enough to change the overall trend of increasing SFR with increasing outflow luminosity.

Fig. 14. Relation between the Hα luminosity ratio between the second and first components and the SFR measured at 10 (blue) and 100 Myr (red points, top panel) and the difference in SFR between these two epochs (bottom panel). Points at the same x-axis value in the top panel represent the same galaxy at different timescales. The Pearson coefficients in the bottom panel indicate a moderate anticorrelation with a higher significance than observed in previous comparisons made in Figs. 8 and 9.

To summarize, although we can observe outflow Hα luminosity contributions of up to ∼20% of the rotating ionized gas observed in the higher-SFR galaxies, we only observe a clear change in recent SFR when comparing absolute values between 10 and 100 Myr. Circumstantial evidence, which could be confirmed with a larger sample, gives similar implications when comparing the outflow luminosity contribution with the SFR ratio and when looking into the individual slopes between the velocity residuals and the velocity dispersion. These results point toward a scenario where the feedback is responsible for suppressing some of the star formation between 100 and 10 Myr timescales.

5.2. Outflow ionization mechanism

As discussed in Sect. 3.4, the comparison of the observed line ratios favors a scenario where the outflowing gas is driven by photoionization, while shock models cannot reproduce the data distribution in the BPT diagrams. These results are somewhat surprising because outflows in star-forming galaxies have been linked to shock excitation in the literature. Rodríguez del Pino et al. (2019) have reported, in their study of nearby (z <  0.015) MaNGA star-forming galaxies, a correlation between the [S II]/Hα line ratio, a known tracer of shocks, and the velocity dispersion of both narrow and broad emission line components, which is not observed in our data (see Fig. 10). López-Cobá et al. (2017), comparing the outflows in the inclined disk galaxy UGC 10043 (z ∼ 0.007, SFR = 0.36 M_⊙ yr⁻¹, log(M_*/M_⊙) = 9.79) to the MAPPINGS III photoionization and shock models, in a very similar fashion as shown in this work, have concluded that shocks do describe the observed line ratios better. Other studies relating star-formation-driven local universe feedback showing signatures of strong shocks include Sharp & Bland-Hawthorn (2010), Ho et al. (2014), and Avery et al. (2021), among others.

Although shock-induced feedback seems to be a majority in the reports found in the literature, some do find evidence that photoionization is the main outflow ionization mechanism in star-forming galaxies. Davies et al. (2019) have reported low [N II]/Hα ratio values for both narrow and broad components in z ∼ 2.3  star-forming sample (9.4 < log(M_*/M_⊙) < 10.8, median log(M_*/M_⊙) = 10.1; main-sequence offset of −0.8 < log(SFR/SFR_MS)< 1.0, median log(SFR/SFR_MS) = 0.3), indicating low contribution of shocks. They also noted the results from Sharp & Bland-Hawthorn (2010), which can give some helpful insights into our results. Studying a sample of 5 nearby starburst galaxies (z <  0.03), Sharp & Bland-Hawthorn (2010) modeled the time evolution of ionizing and mechanical luminosities that arise from star-formation activity, mainly separating the photoionization and the wind-driven shocks mechanisms, for an instantaneous and a continuous star-forming scenarios (see their Fig. 29). For the starburst scenario the photoionization from young stars dominates up to ∼10 Myr, when then shocks from winds created by the radiation pressure increase in luminosity contribution, dominating the luminosity emission up until ∼100 Myr, when both mechanisms fall to negligible values. On the other hand, the continuous star-formation scenario indicates a domination of the photoionization mechanism over wind-driven shocks in all time evolution in general.

As we observe that the photoionization mechanism dominates, the scenario proposed by Sharp & Bland-Hawthorn (2010) is in agreement either with the early stages of the starburst model or with any timescale of the continuous star-formation model. As mentioned in the previous section, we do not see any major changes between the SFRs at 10 and 100 Myr (see Fig. 14), indicating that the galaxies in our sample are better reproduced by the continuous star-formation scenario. Indeed, since these are galaxies lying close to the star-forming main sequence, this seems to be a reasonable scenario. In this case, the Sharp & Bland-Hawthorn (2010) results do not help us to further constrain the outflow timescale, as the photoionization model dominates at all times.

5.3. Outflow properties

As further studies are performed in star-forming feedback, better we constrain the gas outflow properties. However, some properties are harder to constrain than others, which sometimes drives differences between studies when used differently. These are specially the case of the outflow electron density and radial extent. In this section we try to compare these properties with other results in the literature.

As discussed in Sect. 4.2, the outflow extent used in our work is not conventional since we integrate the entire region where we find outflowing gas for each galaxy. The outflow radius for each galaxy is displayed in Fig. 13, and it has a mean value of R_out = 6.8 kpc. While studying a large sample of low-redshift (z <  0.2) luminous infrared galaxies (LIRGs), Arribas et al. (2014) reported that outflows are generally quite compact, usually ∼0.7 kpc in radius, although it can reach up to ∼2 kpc for some galaxies. Davies et al. (2019) used the mean value between the PSF-resolved spatial resolution (0.7 kpc) and maximum radius extent observed in their sample (2.6 kpc), R_out = 1.7 kpc. Assuming biconical geometry outflows in a MaNGA galaxy subsample, Bizyaev et al. (2019) derived that the cone radius is typically ∼1 kpc in the disk plane. These values are all smaller than those that we used, which is expected given the different methodologies.

Regarding the outflow electron density, changes should not be as dramatic, as it is an intensive property, in principle characteristic of the outflowing gas regardless of the area it covers. We find density values between 10 and 730 cm⁻³, with a mean value of cm⁻³, while the papers mentioned above find (or use) 70 cm⁻³ (Bizyaev et al. 2019), cm⁻³ (Davies et al. 2019) and 315 cm⁻³ (Arribas et al. 2014). Higher values from Arribas et al. (2014) and Davies et al. (2019) are probably due to the higher-SFR nature of their galaxy sample, while Bizyaev et al. (2019) electron density could be affected by their outflow geometry assumption, lowering its value in comparison to ours. Regardless of the origins of these differences, they are not large.

Figure 15 displays the difference in the mass outflow rate dependence with the SFR surface density using the values mentioned for these studies. The curves present roughly the same inclination, indicating that the mass loading factor does not vary with these different property measurements, aside from a constant represented by the shift in mass outflow rate for a given SFR surface density. With this comparison we conclude that regardless of the radial extent and electron density measurement strategy we apply, the main results regarding outflow scaling relations are not affected materially.

Fig. 15. Mass outflow rate as a function of the SFR surface density using outflow radial extent and electron density values from different articles. As in Fig. 13, shaded regions display the 3σ distributions, and the equations in the bottom-right corner describe the best-fit models.

5.4. Comparison with other observational results from the literature

The results we obtained for the outflow properties are different than other studies regarding the mass outflow rate and the mass loading factor, usually indicating that we are tracing weaker outflows, but maybe expected for our galaxy sample, which presents moderate star-formation activity.

The slope between the mass outflow rate and Σ_SFR we obtain in this work (0.9 ± 0.3) is considerably steeper than the reported in MaNGA outflows in nearby galaxies by Avery et al. (2021) (0.26 ± 0.08). In their case, they do not separate the AGN contribution in this relation, which could account for this difference since AGN emission is higher in the central parts that usually contribute less to Σ_SFR than disk regions. Those authors found no correlation between the Σ_SFR and η, not far from our results, where we find a quite weak positive correlation, but with a very scattered distribution.

The traced outflows in our sample are in general quite weak in comparison with other studies, not only represented by the very low contribution to the bolometric luminosity (∼5.2 × 10⁻⁴%), but also in the mass loading factor, with values of 0.01 <  η <  1.6 with a median of 0.12. This is lower than reported by Rodríguez del Pino et al. (2019) for their star-forming MaNGA sample (median value η ∼ 0.25) and more distant to the values reported for galaxies with higher SFRs (η ∼ 0.3–0.5, Arribas et al. 2014; Davies et al. 2019). Avery et al. (2021) report values of outflow momentum rate ∼10% of total momentum output, including AGNs (although the AGN luminosity contribution is usually less than 20%), while we observe ∼0.7%, even though their η values seem to be similar to what we obtain (see their Fig. 11). In another study using MaNGA data, tracing outflows using NaD neutral gas, Roberts-Borsani et al. (2020) obtain η values also somewhat higher than we do (−1 <  log η <  1, with a median probably η ≳ 1). Zaragoza-Cardiel et al. (2020), using a star-formation self-regulator model in a large sample of galaxy disks with MUSE data, find 1 <  η <  5 for 7.5 <  log Σ_* <  9 M_⊙ kpc⁻².

As discussed in Sect. 5.3, the outflow properties can change the estimated contribution, also affecting the mass loading factor. If we use the outflow extent and the electron density from Arribas et al. (2014) and Davies et al. (2019), however, we obtain median values of η = 0.22 and η = 0.08, respectively, still lower than the reported values. This indicates that the luminosity contribution of the outflows is indeed quite low for the SFR regime probed by our sample. Thus the discrepancy between our results and the others from the literature could be explained by either the deeper observations (in comparison to MaNGA data, for example) or the lower SFRs (e.g., in comparison to ultra-LIRGs and LIRGs).

5.5. Comparison with theoretical results

Although our mass loading factor values do not exactly match the values reported from other observational studies, it seems to be in agreement with theoretical studies. FIRE simulations presented by Muratov et al. (2015), modeling small-scale feedback up to z = 0, indicate that mass loading factors for disk progenitors decrease to η ≪ 1 at low redshift. Using a small-box simulation resolving supernova remnants in the ISM, Li & Bryan (2020) report values of η ≳ 0.08 for their breakout wind ISM condition, where hot outflows are considered, higher than the partial breakout and cold outflow models (0.01 <  η <  0.08 and η <  0.01, respectively). Schneider et al. (2020) also report a 0.1 <  η <  0.2 for the outflow hot gas phase, depending on the outflow radial extent.

Usually feedback theoretical studies prefer to characterize the outflow properties, such as the mass loading factor, along with the stellar mass instead of the SFR, as we have used throughout this work. Figure 16 displays the dependence of the mass loading factor with the stellar-mass surface density. η values are the same as displayed in Fig. 13, while Σ_M* was obtained integrating the stellar mass derived from STARLIGHT fit within the same outflowing radius determined in Sect. 4.2. We observe a clearer correlation between these parameters, in contrast with the comparison between η and Σ_SFR. The obtained relation is given by

(8)

Fig. 16. Relation between the stellar-mass surface density and the mass loading factor. As in the other figures, the best-fit model with a 3σ uncertainty is represented by the black and the blue shaded region and is described by the equation given with the Pearson coefficients.

The correlation between the mass loading factor and the SFR or stellar mass has been debated in recent studies. Although observational relations have been observed to be positive in several scenarios (e.g., Arribas et al. 2014; Davies et al. 2019; Avery et al. 2021, including this work; although Zaragoza-Cardiel et al. (2020) do find an anticorrelation), models seem to indicate otherwise (e.g., Creasey et al. 2013; Muratov et al. 2015; Li et al. 2017). Several aspects of the feedback characterization could be driving this discrepancy, such as the different gas phase measurements, physical properties (as discussed in Sect. 5.3), and the galaxy scales probed. Further constraints in this topic seem to be necessary to understand what is the reason for this slope uncertainty.

6. Conclusions

In this work we have used detailed MUSE and GAMA observational data of 15 low-redshift (z ∼ 0.15) star-forming galaxies included in the VALES sample in order to characterize the ionized gas feedback due to star-formation activity.

We initially measured the emission-line spectra with a single component in order to investigate the kinematical deviation due to feedback from the galaxy rotation, which we have modeled. We do not observe evidence of an increase in residual velocities v_res when compared to the velocity dispersion σ and SFR surface density Σ_SFR for the entire sample, although we do observe this for some individual galaxies. We then applied a double-component fit for galaxy regions where the initial measurements were not satisfactory in order to trace a broader component, interpreted to be related with outflowing gas. The positive relation between the second-component velocity dispersion σ_2nd and Σ_SFR strengthens this interpretation.

We used the results obtained from the ionized gas emission lines to compare with the SED modeling of the GAMA photometric broadbands. Modeling the SFHs of the sample, we could conclude that the SFR at 10 Myr is somewhat suppressed in galaxies with higher Hα luminosity contribution from outflows, in comparison with the SFR at 100 Myr ago. The effect of the feedback at these timescales is also observed in comparison to the σ versus v_res slope for individual galaxies obtained from the single-component emission-line fit, although with less significant statistics.

When comparing the gas excitation with shocks and photoionization models, we see that the observed line ratios are better reproduced by the ionization from young stellar population radiation, in contrast with the bulk of recent reports on feedback analysis that suggest that shocks are most common in stellar feedback. This difference could be explained by the broadly constant recent star formation (< 100 Myr), while a “burstier” activity could result in shock-dominated ionization.

The derived outflow properties, such as the mass outflow rate (∼0.1 M_⊙ yr⁻¹), the outflow kinetic power (∼5.2 × 10⁻⁴%L_bol), and the mass loading factor (∼0.12), indicate that the outflows we observe are weak and probably have little effect on the galaxies’ evolution. Even though stellar feedback seems to be a common phenomenon within star-forming galaxies, as we have detected outflowing regions in 13 out of 15 galaxies, our results indicate that feedback in low-redshift star-forming galaxies is not effective in large gas removal and dramatic star-formation quenching, resulting in minor (and apparently local) effects within these galaxies.

Acknowledgments

G. S. C. acknowledges the support from the Comití Mixto ESO-Chile, the DGI at University of Antofagasta, and from CONICYT/ANID FONDECYT project No. 3190561. G. S. C. would also like to thank Laurent Chemin for insightful discussions about the rotation model results in this work. T. M. H. acknowledges the support from the Chinese Academy of Science (CAS) and the National Commission for Scientific and Technological Research of Chile (CONICYT) through a CAS-CONICYT Joint Postdoctoral Fellowship administered by the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. M. B. acknowledges partial support from FONDECYT regular grant 1170618 and 1211000. E. I. acknowledges partial support from FONDECYT through grant N° 1171710. R. L. acknowledges support from CATA, BASAL grant AFB-170002. Y. Q. X. acknowledges support from NSFC-12025303, 11890693, 11421303, the CAS Frontier Science Key Research Program (QYZDJ-SSW-SLH006), the K.. C. Wong Education Foundation, and the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A06.

References

Appendix A: Parameters of individual galaxies

Here we present relevant information for each individual galaxy, shown in Figs. A.1-A.15, including distribution maps, physical parameters and SED and BPT plots. The caption in Fig. A.1 details these information and applies to the other figures.

Fig. A.1. Individual data of HATLASJ083601. Listed from top left to bottom right: (first row) Hα emission-line integrated fluxes of the one-component fit and the first and second component of the two-component-fit maps; (second row) maps of the stellar-mass surface density, SFR surface density in logarithmic scale, sSFR surface density in logarithmic scale, and nebular metallicity; (third row) maps of the Hα peak velocity, rotation model, residuals between Hα peak velocity and the rotation model, and velocity dispersion; (fourth row) maps of the [N II]/Hα and [O III]/Hβ ratios, BPT coloring weighted by the Hα flux (using an HSV (hue, saturation, value) color scheme), and BPT diagram; (fifth row) maps of the electron density and optical extinction, and the relations between velocity dispersion and SFR surface density with the rotation model residual velocities; (sixth row) SED fit from CIGALE and SFHs obtained from this fit. A list of physical parameters is also given in the top-right corner of the figure. Aside from the identified Hα flux maps, all maps are retrieved from emission-line measurements. Stellar masses are obtained from the STARLIGHT stellar continuum fit. Nebular metallicity was estimated using the O3N2 calibrator with the relation presented by Marino et al. (2013). The BPT diagram (Baldwin et al. 1981) also displays in dashed lines the commonly used theoretical (Kewley et al. 2001) and empirical (Kauffmann et al. 2003) lines that aim to separate pure star formation and star formation combined with AGN ionization, respectively. Each point in the BPT diagram represents a Voronoi region. Electron densities are obtained by using the PYNEB package (Luridiana et al. 2015) and the [S II]λ6717/6731 ratio, assuming a nebular temperature of 10 000 K. Plots of the dependences of residual velocities with velocity dispersion and SFR surface density are separated into all-Voronoi regions, only regions presenting redshifted residuals, and regions presenting blueshifted residuals, from top to bottom, respectively. Green lines represent a simple fit between the parameters. Finally, the SED fit shows the UV/optical and IR bands fitted in blue and red dots. Blue, orange, green, and red lines represent stellar attenuated emission, stellar unattenuated emission, nebular emission, and dust emission, respectively.

Fig. A.2. Same as Fig. A.1, but for galaxy HATLASJ083832.

Fig. A.3. Same as Fig. A.1, but for galaxy HATLASJ084217.

Fig. A.4. Same as Fig. A.1, but for galaxy HATLASJ085111.

Fig. A.5. Same as Fig. A.1, but for galaxy HATLASJ085346.

Fig. A.6. Same as Fig. A.1, but for galaxy HATLASJ085356.

Fig. A.7. Same as Fig. A.1, but for galaxy HATLASJ085406.

Fig. A.8. Same as Fig. A.1, but for galaxy HATLASJ085450.

Fig. A.9. Same as Fig. A.1, but for galaxy HATLASJ085828.

Fig. A.10. Same as Fig. A.1, but for galaxy HATLASJ085957.

Fig. A.11. Same as Fig. A.1, but for galaxy HATLASJ090532.

Fig. A.12. Same as Fig. A.1, but for galaxy HATLASJ090949.

Fig. A.13. Same as Fig. A.1, but for galaxy HATLASJ091157.

Fig. A.14. Same as Fig. A.1, but for galaxy HATLASJ114244.

Fig. A.15. Same as Fig. A.1, but for galaxy HATLASJ142128.

All Tables

All Figures

Fig. 1. Example of spectral fittings for one galaxy (HATLASJ083832). Top left: CIGALE SED fit, where points represent the broadband magnitudes from the GAMA survey and the red, orange, blue, and green lines represent the dust spectra, stellar spectra, stellar spectra attenuated by the dust, and emission-line continuum spectra, respectively. Blue points represent the UV and optical bands, and the red points represent the IR bands. Top right: STARLIGHT stellar continuum fit, with the black and red lines representing observed and modeled spectra, respectively. Bottom left: single-Gaussian fit of the Hα+[N II] emission lines, where the points and the black line represent the observed spectrum, and the red line represents the modeled lines. A residual spectrum is shown in the bottom part of the panel. Bottom right: same as the previous panel, but showing the double-Gaussian fit, where both components are represented by blue lines. The Gaussian and STARLIGHT fits are from the same spaxel, located in the central region of the galaxy.

In the text

Fig. 2. Quality of the SED modeling of the lower-wavelength photometric broadbands (from FUV to WISE2), with the addition of the Hα and Hβ integrated fluxes from MUSE data. Total reduced chi-square and median residuals (in percentage of the observed magnitudes from GAMA and MUSE), for both the estimated value for the VALES sample and the exact value for the mock sample, are displayed in the left and right panels, respectively.

In the text

	Fig. 3. Galaxy sample distribution in the M*-SFR diagram (blue dots). The median of the distribution is represented by the yellow star. Values were obtained from the CIGALE fit of the sample. Dashed blue and red lines represent the star-forming main sequence and retired-galaxy relations obtained by Cano-Díaz et al. (2019). Dotted lines indicate SFR offsets by the displayed multipliers.
In the text

Fig. 4. Relation between the velocity dispersion σ (left) and SFR density ΣSFR (right) with the residual velocities obtained from the rotation model fitted to the Hα velocity fields. Each point represents a Voronoi region. The top panels show the distribution for the whole sample, while the middle and bottom panels are separated into redshifted and blueshifted residuals, respectively. Square symbols represent median values, with the dashed lines representing 25 and 75 percentiles of the distribution. No clear correlation between the two parameters is observed with vres.

In the text

Fig. 5. Relation between the velocity dispersion of the second component and the flux ratio between the second and the first components, from the double-component fit. Left panel: all the Voronoi regions where a double-component fit was performed. These include what we call the “boundary regions” (explained in the text), and the number of such regions, discriminated by criteria, is shown in the top right. Right panel: distribution without boundary regions.

In the text

Fig. 6. Relations between SFR surface density and sSFRs and the second-over-first-component flux ratio. Black points represent Voronoi regions, where the point size is proportional to the inverse of the Voronoi covered area, so small points represent regions where more spaxels had to be binned to reach a satisfactory S/N. Strong blue lines represent a running median with a 20-point window, with weak blue lines representing the 1σ of the distribution.

In the text

Fig. 7. Second-component velocity dispersion as a function of the SFR surface density. Voronoi regions are represented by black points, with point size being anticorrelated with the Voronoi region area. The black line and blue shaded region show the best-fit relation with a 3σ uncertainty, and the relation is displayed in the top-left corner. Pearson coefficients are also displayed, indicating a significant positive correlation.

In the text

Fig. 8. Relation between the slopes of log ΣSFR (left panel) and the one-component fit σ (right) versus residual velocities and the SFR10 Myr/SFR100 Myr ratio obtained with the SED fit. The slopes are obtained through the fits shown for each individual galaxy in Appendix A. Pearson coefficients are given in the top-right corners; they show no correlation in the first relation and only weak anticorrelation with little significance in the second. The cyan line shows the SFR10 Myr/SFR100 Myr = 1 ratio.

In the text

	Fig. 9. Dependence of the SFR10 Myr/SFR100 Myr ratio with the spatially integrated Hα luminosity ratio between the second and the first components. Pearson coefficients are given in the top-right corner; they show a moderate anticorrelation with a small significance.
In the text

Fig. 10. BPT diagrams (top panels) and the dependence of the second-component velocity dispersion σ2nd on the [N II]/Hα, [O III]/Hβ, and [S II]/Hα line ratios of the Voronoi regions where we do measure a second component related to outflowing gas (bottom panels). The BPT diagrams are colored with its σ2nd. Dashed lines display the theoretical (Kewley et al. 2001) and empirical (Kauffmann et al. 2003) regions, separating pure star formation and star formation combined with AGN ionization, respectively. Red, black, and green grids display MAPPINGS V shock models for LMC (Z = 0.0071), Dopita_2005 (Z = 0.0065), and solar metallicities (Z = 0.0183), respectively. Grids are for electron density ne = 1 cm−3, with magnetic fields 0.0001 > B < 10 μG and shock velocities 100 <  vshock <  400 km s−1 (shock velocities usually vary along the x axis of the BPT diagram). Some shock velocity values are indicated, as is the general direction of increasing magnetic field intensity. Blue lines in the bottom panels represent the running median with the 1σ distribution of the Voronoi regions. Dashed green lines represent solar metallicity shock models with B = 1 μG and 0.01 <  ne <  1000 cm−3, identified by the labels. Y-axis values for the dashed green lines are vshock.

In the text

Fig. 11. BPT diagrams as shown in Fig. 10, here in comparison with the photoionization grids obtained from the 3MdB database using CLOUDY models. The top and bottom panels show the model results for Z = 0.0076 and Z = 0.019, respectively. For comparison, the sample median metallicity is Z = 0.0086. Blue lines represent different stellar ages of 1, 3, 4, 5.6, and 8.9 × 106 years as ionization sources and are roughly oriented from top right to bottom left as age increases, although not linearly. Red lines represent variation in the photoionization parameter, ranging from log(U) between −4.0 and −2.0 (in steps of 0.5), roughly oriented from top left to bottom right as U decreases. Some stellar age values are displayed in the top-left panel, along with the increasing U orientation.

In the text

Fig. 12. Velocity shift between the two emission-line components and outflow velocity as a function of SFR surface density for the Voronoi regions, shown as black points. Blue lines in the left panel display the running median with a 1σ distribution. The black line and blue shaded region in the right panel represent the best-fit relation (displayed in the top-left corner, along with Pearson coefficients) with a 3σ uncertainty.

In the text

Fig. 13. Outflow properties of the analyzed sample. All panels are in function of the outflow SFR surface density, except the top left. Top left: outflow extent (definition described in the text) as a function of the second-component Hα luminosity. Top right: electron density. Middle left: mass outflow rate. Middle right: mass loading factor. Bottom left: ratio between the outflow kinetic power and the total bolometric luminosity. Bottom right: ratio between the momentum outflow rate and the SFR photon momentum. We display the Pearson coefficients for five of the displayed relations, and the best-fit model with a 3σ uncertainty for the panels is shown as a black line with blue shaded regions.

In the text

Fig. 14. Relation between the Hα luminosity ratio between the second and first components and the SFR measured at 10 (blue) and 100 Myr (red points, top panel) and the difference in SFR between these two epochs (bottom panel). Points at the same x-axis value in the top panel represent the same galaxy at different timescales. The Pearson coefficients in the bottom panel indicate a moderate anticorrelation with a higher significance than observed in previous comparisons made in Figs. 8 and 9.

In the text

	Fig. 15. Mass outflow rate as a function of the SFR surface density using outflow radial extent and electron density values from different articles. As in Fig. 13, shaded regions display the 3σ distributions, and the equations in the bottom-right corner describe the best-fit models.
In the text

	Fig. 16. Relation between the stellar-mass surface density and the mass loading factor. As in the other figures, the best-fit model with a 3σ uncertainty is represented by the black and the blue shaded region and is described by the equation given with the Pearson coefficients.
In the text

Fig. A.1. Individual data of HATLASJ083601. Listed from top left to bottom right: (first row) Hα emission-line integrated fluxes of the one-component fit and the first and second component of the two-component-fit maps; (second row) maps of the stellar-mass surface density, SFR surface density in logarithmic scale, sSFR surface density in logarithmic scale, and nebular metallicity; (third row) maps of the Hα peak velocity, rotation model, residuals between Hα peak velocity and the rotation model, and velocity dispersion; (fourth row) maps of the [N II]/Hα and [O III]/Hβ ratios, BPT coloring weighted by the Hα flux (using an HSV (hue, saturation, value) color scheme), and BPT diagram; (fifth row) maps of the electron density and optical extinction, and the relations between velocity dispersion and SFR surface density with the rotation model residual velocities; (sixth row) SED fit from CIGALE and SFHs obtained from this fit. A list of physical parameters is also given in the top-right corner of the figure. Aside from the identified Hα flux maps, all maps are retrieved from emission-line measurements. Stellar masses are obtained from the STARLIGHT stellar continuum fit. Nebular metallicity was estimated using the O3N2 calibrator with the relation presented by Marino et al. (2013). The BPT diagram (Baldwin et al. 1981) also displays in dashed lines the commonly used theoretical (Kewley et al. 2001) and empirical (Kauffmann et al. 2003) lines that aim to separate pure star formation and star formation combined with AGN ionization, respectively. Each point in the BPT diagram represents a Voronoi region. Electron densities are obtained by using the PYNEB package (Luridiana et al. 2015) and the [S II]λ6717/6731 ratio, assuming a nebular temperature of 10 000 K. Plots of the dependences of residual velocities with velocity dispersion and SFR surface density are separated into all-Voronoi regions, only regions presenting redshifted residuals, and regions presenting blueshifted residuals, from top to bottom, respectively. Green lines represent a simple fit between the parameters. Finally, the SED fit shows the UV/optical and IR bands fitted in blue and red dots. Blue, orange, green, and red lines represent stellar attenuated emission, stellar unattenuated emission, nebular emission, and dust emission, respectively.

In the text

	Fig. A.2. Same as Fig. A.1, but for galaxy HATLASJ083832.
In the text

	Fig. A.3. Same as Fig. A.1, but for galaxy HATLASJ084217.
In the text

	Fig. A.4. Same as Fig. A.1, but for galaxy HATLASJ085111.
In the text

	Fig. A.5. Same as Fig. A.1, but for galaxy HATLASJ085346.
In the text

	Fig. A.6. Same as Fig. A.1, but for galaxy HATLASJ085356.
In the text

	Fig. A.7. Same as Fig. A.1, but for galaxy HATLASJ085406.
In the text

	Fig. A.8. Same as Fig. A.1, but for galaxy HATLASJ085450.
In the text

	Fig. A.9. Same as Fig. A.1, but for galaxy HATLASJ085828.
In the text

	Fig. A.10. Same as Fig. A.1, but for galaxy HATLASJ085957.
In the text

	Fig. A.11. Same as Fig. A.1, but for galaxy HATLASJ090532.
In the text

	Fig. A.12. Same as Fig. A.1, but for galaxy HATLASJ090949.
In the text

	Fig. A.13. Same as Fig. A.1, but for galaxy HATLASJ091157.
In the text

	Fig. A.14. Same as Fig. A.1, but for galaxy HATLASJ114244.
In the text

	Fig. A.15. Same as Fig. A.1, but for galaxy HATLASJ142128.
In the text