The MUSE Atlas of Discs (MAD): Ionized gas kinematic maps and an application to diffuse ionized gas Free

ABSTRACT

We have obtained data for 41 star forming galaxies in the MUSE Atlas of Discs (MAD) survey with VLT/MUSE. These data allow us, at high resolution of a few 100 pc, to extract ionized gas kinematics (V, σ) of the centres of nearby star forming galaxies spanning 3  dex in stellar mass. This paper outlines the methodology for measuring the ionized gas kinematics, which we will use in subsequent papers of this survey. We also show how the maps can be used to study the kinematics of diffuse ionized gas for galaxies of various inclinations and masses. Using two different methods to identify the diffuse ionized gas, we measure rotation velocities of this gas for a subsample of six galaxies. We find that the diffuse ionized gas rotates on average slower than the star forming gas with lags of 0–10 km s⁻¹ while also having higher velocity dispersion. The magnitude of these lags is on average 5 km s⁻¹ lower than observed velocity lags between ionized and molecular gas. Using Jeans models to interpret the lags in rotation velocity and the increase in velocity dispersion we show that most of the diffuse ionized gas kinematics are consistent with its emission originating from a somewhat thicker layer than the star forming gas, with a scale height that is lower than that of the stellar disc.

1 INTRODUCTION

The global motions of stars and gas in galaxies are determined primarily by the mass distribution inside the galaxy. In particular, the analysis of rotation curves in the outer parts of galaxies based on neutral hydrogen observation has led to the firm establishment of dark matter in galaxies (e.g. Rubin, Ford & Thonnard 1980; Bosma 1981; van Albada et al. 1985, and many others). Although the 21 cm line has been one of the most widely used lines for this, it is also possible to trace rotation curves with ionized (e.g. Mathewson, Ford & Buchhorn 1992; Garrido et al. 2002; Erroz-Ferrer et al. 2016) or molecular gas (e.g. Sofue 1996; Sofue et al. 1997). As these emission lines are produced by gas in different physical states, they are often found in different parts of the galaxy.

Much can be learned by comparison of kinematic tracers with each other. By comparing the velocity dispersion of low surface brightness CO, Caldú-Primo et al. (2013) and Mogotsi et al. (2016) provide evidence for a faint, diffuse, higher dispersion CO component in nearby spiral galaxies that appears similar in dispersion (and therefore in thickness) as the neutral hydrogen. The kinematics of different tracers can be different. Molecular gas has on average a lower velocity dispersion than atomic gas. Gas dynamical tracers have in general a lower velocity dispersion than stars. For stars in spiral galaxies, the velocity dispersion is in turn a function of the age of the population, with older stars showing higher dispersion and a bigger vertical extent (Wielen 1977; Carlberg et al. 1985).

The velocities of ionized gas should be very close to those of the molecular gas, as the stars that are responsible for the ionization have formed very recently from molecular gas and have a low asymmetric drift (e.g. Quirk et al. 2018). However, molecular and ionized gas do not always agree with each other. Davis et al. (2013) study early-type galaxies and find that the ionized gas, contrary to the molecular gas, does not necessarily trace the circular velocity of the galaxies, mostly because of a different distribution of the molecular gas. Recently Levy et al. (2018) compared ¹²CO (J=1–0) rotation curves with H α rotation curves derived from CALIFA data (Sánchez et al. 2012). They found that the H α gas usually, but not always, shows velocity lags, with median values between 0 and 25 km s⁻¹, with respect to the CO gas, which they attribute to the presence of extraplanar diffuse ionized gas.

First identified in the Milky Way (Reynolds 1984, but already suggested by Hoyle & Ellis 1963), a major part of the H α emission from star forming galaxies comes from diffuse ionized gas (DIG) instead of directly from H ii regions (Reynolds 1990; Walterbos & Braun 1994; Zurita, Rozas & Beckman 2001). The excitation mechanisms of this diffuse gas are not completely understood, but have been attributed to UV radiation from leaky H ii regions (Hoyle & Ellis 1963; Reynolds 1984), cosmic rays (e.g. Dahlem, Dettmar & Hummel 1994; Vandenbroucke et al. 2018), turbulent mixing layers (Binette et al. 2009; Haffner et al. 2009) and evolved stars (e.g. Kaplan et al. 2016; Zhang et al. 2017). All of these lead to different line ratios in optical strong emission lines (Haffner, Reynolds & Tufte 1999; Hoopes & Walterbos 2003; Madsen, Reynolds & Haffner 2006). The layer of ionized gas in the Milky Way, the Reynolds layer, extends to beyond a kpc above the galactic disc (Reynolds 1989). Also in external galaxies, such diffuse gas layers have been found (e.g. Rossa & Dettmar 2003), extending on average 1–2 kpc above the galactic plane.

The study of the kinematics of this DIG is interesting as it potentially could reveal the origin and ionization source of the DIG.

The kinematics of ionized gas have traditionally been studied with long slits (e.g. Rubin et al. 1980; Mathewson et al. 1992), as well as Fabry–Perot (FP) interferometers (Tully 1974; Ryder et al. 1998; Hernandez et al. 2008). Recently, 2D mapping of ionized gas has become possible by the developments of integral field units such as SAURON (Bacon et al. 2001), SparsePak (Bershady et al. 2004), PMAS (Roth et al. 2005), and VIRUS-P (Hill et al. 2008) to study ionized gas kinematics in late-type galaxies (Ganda et al. 2006; Bershady et al. 2010; García-Lorenzo et al. 2015).

Each of these methods for studying ionized gas has its advantages and disadvantages. The FP interferometers provide higher spectral resolution but at the cost of wavelength baseline. Some of the mentioned IFUs have lower spectral resolution and/or lower spatial resolution than the FP interferometers, but provide a longer wavelength range. The Multi-Unit Spectroscopic Explorer (MUSE) instrument (Bacon et al. 2010) on the Very Large Telescope (VLT) falls in between the properties of the mentioned instruments. Although its spectral resolution is not as high as that of FP interferometers, its spatial resolution is superb, and its large spectral range at intermediate resolution, combined with a 1 arcmin × 1 arcmin field of view, makes this instrument also suitable to study more diffuse gas that has surface brightness comparable to the stellar continuum emission. This is a particular niche for MUSE that has been difficult to observe with previous instruments.

In this paper we compare the kinematics of the DIG with that of the star forming gas for a sample of star forming galaxies with various inclinations. Previous studies of the kinematics of the diffuse ionized gas have mainly focused on extraplanar gas in (almost) edge-on galaxies such as NGC 891 (Heald et al. 2006b), NGC 4302 (Heald et al. 2007), NGC 5775 (Rand 2000; Tüllmann et al. 2000; Heald et al. 2006a) and NGC 2403 (Fraternali, Oosterloo & Sancisi 2004), and recently on survey-scale with 67 edge-on galaxies in the Mapping Nearby Galaxies at APO survey (MaNGA; Bizyaev et al. 2017). There seems to be a consensus that the rotational velocity of the gas above the plane is lower than that of the gas in the mid-plane, although this gas may show evidence for non-circular motions (Fraternali et al. 2004). There exist no measurements of the kinematics of diffuse gas in relatively face-on galaxies with the exception of the work of Boettcher, Gallagher & Zweibel (2017), who also find lags in the rotation velocity of DIG in M 83.

For the extraplanar gas, models have been developed to explain this slower rotation. Collins, Benjamin & Rand (2002) use ballistic models to explain this gas. The clouds, which are launched from the mid-plane of the galaxy disc, show a decrease in rotation signal as they travel away from the plane. The clouds also move radially outwards due to a combination of conservation of angular momentum and lower gravitational force. Barnabè et al. (2006) (see also Benjamin 2002) proposed hydrostatic models to predict the rotation velocity of gas outside the galactic mid-plane. These models also explain the velocities of extraplanar diffuse gas, although they require high temperatures (10⁴ ≲ T ≲ 10⁶ K) for the gas, which might be more appropriate for galactic coronae than for the cold gas that has been observed. Fraternali & Binney (2006) and Marinacci et al. (2010) develop models in which there is a continuous launching of clouds from the disc, which eventually lose velocity due to a drag with a hot corona.

This paper presents some of the first results of the MUSE Atlas of Discs (MAD) Survey. This survey is mapping out the inner parts of a sample of 45 (mostly) nearby star forming galaxies with MUSE. The goal of the survey is to understand the formation and evolution of disc galaxies through studies of the kinematics of gas and stars in these galaxies, together with their star formation histories to understand how disc galaxies have formed and evolved. The description of the survey will be presented in Carollo et al. (in preparation, henceforth Paper 1), as well as in Erroz-Ferrer et al. (2019, henceforth Paper 2) which describes the resolved metallicity and star formation properties in the galaxies. With the high spatial resolution of MUSE, we check the scenario proposed by Levy et al. (2018) that H α rotation curves may be lagging the ¹²CO (J=1–0) rotation curves because of the presence of extraplanar DIG. We do this by separating the star forming gas from the diffuse ionized gas and analysing the kinematics separately. We also outline the procedures used for deriving the kinematics, which we will use in future papers. This paper is structured as follows. In Section 2, we briefly recapitulate the sample selection and data reduction. In Section 3 we present the procedures for the derivation of the kinematics and present the kinematic maps. Section 4 shows the methods to identify DIG and the rotational velocity measurements. The rotation velocities of DIG and star forming (SF) gas are then presented in Section 5, followed by a discussion in Section 6.

2 OBSERVATIONS AND DATA REDUCTION

The MAD survey has observed a sample of 45 nearby star forming galaxies with VLT/MUSE as part of a Large Guaranteed Time Program and provides us with high spatial and spectral resolution observations of the centres of these galaxies. Although the sample selection and data reduction will be described in detail in paper I (Carollo et al. in preparation), we briefly summarize the selection criteria and reduction procedures here as well.

Besides observability constraints (visible for at least one consecutive hour per night from VLT with airmass <1.5), galaxies were selected to be bright (M_B < 13) and on the star forming main sequence (Brinchmann et al. 2004; Daddi et al. 2007; Noeske et al. 2007). In order to make optimal use of the field of view of MUSE, a size cut was applied to ensure that galaxies were not significantly smaller than the field of view and that our observations would extend to at least 0.75 effective radius. To ensure a uniform coverage of spectral features, the host galaxy redshifts were limited to the range z = 0.002–0.012. Target galaxies were additionally selected to have good-quality optical Hubble Space Telescope (HST) imaging in at least one red passband (F606W, F814W, or equivalent) available from either HST/WFPC2, HST/ACSWFC, or HST/WFC3, moderate inclinations (0.3 < ϵ < 0.95), and limited foreground extinction. The final hand-picked sample contains 45 galaxies. In this paper we discuss 41 sample galaxies; we note that the remaining four galaxies were chosen as a reference sample containing merging galaxies. The data reduction and analysis of these remaining galaxies is more complicated, as two of these galaxies are at cosmological distances, and the other two are large mosaics. As the kinematics of these particular four galaxies have little relevance for this paper, we postpone their presentation to a future paper. The 41 sample galaxies are listed in Table 1 and shown in Fig. 1. All galaxies were observed as part of the MUSE GTO observations, except for NGC 337, which was observed during the commissioning of MUSE, and NGC 1097, which was observed by program 097.B-0640 (PI Gadotti).

We reduced the data with the MUSE Instrument Pipeline (Weilbacher, Streicher & Palsa 2016, version 1.2 to 2.2, depending on when the galaxy was observed). We used the pipeline to perform basic steps like bias subtraction, flat fielding, and wavelength calibration. Each galaxy was observed with multiple (at least 3) exposures, which were interlaced with sky integrations. The total on-target exposure time was always 1 h per galaxy. We experimented with different sky observation strategies for the first few galaxies, until we decided on a Object–Sky–Object observation pattern with 2 min sky observations. To minimize the influence of bad pixels, the on-target observations were dithered with a small few-arcsecond dither pattern. We aligned the individual exposures by generating a narrow-band H α image for each cube and using the image registration task tweakreg from the DrizzlePac package (Avila et al. 2015) to measure shifts between different cubes. We then used the pipeline to drizzle the individual exposures to a common reference cube for each galaxy.

We used zap (Soto et al. 2016) to subtract the sky from the individual exposures. zap performs a principal component analysis (PCA) on the separate sky observation, and encaptures the variation in the line spread function and the Poisson noise of bright sky lines in this way. After subtraction of the average sky determined on the sky frame from a target observation, zap also subtracts residuals in the sky subtraction by fitting the PCA components to the spectra. Our galaxies were processed with zap versions 1.0 and 2.0, depending on the observation date of the galaxy. We inspected the sky subtracted by zap to ensure that no galaxy light was subtracted.

The sky subtracted cubes were then combined to a single cube for each galaxy by taking the median value. We then corrected this cube for Milky Way foreground extinction using the values of Schlafly & Finkbeiner (2011), and masked foreground stars based on a by-eye identification on HST images and spectral identification in the MUSE cubes. We excluded the centre of NGC 3783 from the analysis as the broad Balmer lines from the type 1 AGN dominate most of the spectrum in this region.

3 DERIVATION OF THE KINEMATICS

In order to derive accurate velocities and dispersions from gas emission lines, we first remove the stellar continuum by modelling this continuum with synthesized stellar templates. As the signal-to-noise of the continuum is not high enough to obtain good fits to the stellar continuum, we bin the data using the voronoi package of Cappellari & Copin (2003). To avoid stellar absorption lines and sky emission lines we use a 100 Å wide area around 5700 Å to determine the signal-to-noise. We then bin the spectra to achieve a S/N of 50 per Å.

We fit PEGASE-HR templates (Le Borgne et al. 2004) to the stellar continuum using pPXF (Cappellari 2017). The spectral resolution of this library (R = 10 000) is higher than the MUSE resolution and we therefore convolve the libraries with a wavelength-dependent line spread function (LSF). We determine the MUSE LSF using both arc and sky lines, and after ensuring that our findings are similar to the one used by Guérou et al. (2017) we adopt their parametrization of the LSF. Before fitting with pPXF we mask all spectral pixels within a 400 km s⁻¹ window of a strong emission line (Table 2) as well as the region around the NaD absorption lines. During the fit we allow for an additive polynomial up to order 4.

For analysing the ionized gas emission lines we subtract the best-fitting continuum pixel-by-pixel. This is justified as long as the stellar populations do not change radically within one Voronoi bin. For each spaxel we scale the best-fitting spectrum to the spectrum of the spaxel using the median value of each spectrum in the fitted range. We do not see evidence for any systematics in the continuum subtracted cubes. Further discussion on the reliability and possible influence of the continuum subtraction can be found in Section 6.1.1.

We then proceed to analyse the emission lines. Even though the emission lines are often much brighter than the continuum, it is in many cases still necessary to bin pixels together to get a sufficiently high S/N for analysing the emission lines kinematics. In order to estimate the signal-to-noise of the H α line without knowing its actual velocity, flux or dispersion we take the height of the H α line divided by the average dispersion of the surrounding continuum as a measure for the signal-to-noise. This is usually a lower limit to the significance with which the line is detected, but can sometimes lead to an overestimate of the S/N. After we have estimated the S/N per pixel, we tesselate the galaxy to reach S/N of 10 per bin using the aforementioned voronoi code. For the most massive galaxies in the sample (e.g. NGC 4030, NGC 3521), this leads to bins with essentially the size of 1 spaxel.

We fit a Gaussian line for each of the emission lines to the spectra binned using the H α tesselation. We note that some authors (e.g. Boettcher et al. 2017) use combinations of broad and narrow Gaussians, but given the 2.5 Å spectral resolution of MUSE such a decomposition is not always warranted by our data. However, we will discuss the use of such a decomposition in Appendix  C. We broaden each line by the LSF at that wavelength and a velocity dispersion.

Different lines can trace gas with different physical conditions and therefore do not need to necessarily have the same kinematics. Here, we divide the lines in two groups, one containing the two visible Balmer lines, and the other group the other lines from Table 2. Lines inside each group share common kinematics (V, σ) but not common fluxes. Both groups of lines are fit simultaneously to the continuum-subtracted binned spectra using the Levenberg–Marquardt algorithm. Unless specified otherwise, the kinematics used in the remainder of the paper are based on the kinematics of the two Balmer lines.

We present the velocity and velocity dispersion of the H α and H β lines in Figs 2 and 3. The velocity maps show often irregular structure. NGC 3256 is a merging system (Vorontsov-Velyaminov 1959) with two nuclei separated ∼5 arcsec in projection (Lira et al. 2008) for which the observed velocity field is very irregular. NGC 337 is an asymmetric galaxy with an off-centred bar, which too has been suggested to be a merger (Sandage & Bedke 1994).

We note that the velocity dispersions that we measure for the ionized gas are often below the instrumental resolution. At the wavelength of H α, which, out of the two Balmer lines in our spectral range, is the dominant line for the determination of the kinematics because of its generally higher S/N, the instrumental resolution is σ ≈ 47 km s⁻¹. Velocity dispersions of ionized gas in disc galaxies are often between 10 and 50 km s⁻¹ (e.g. Epinat et al. 2008; Erroz-Ferrer et al. 2015). It is therefore a priori not clear how reliable the determined velocity dispersions are. In Appendix  A we perform a set of idealized simulations to see if we can recover the intrinsic dispersion of the gas. Our conclusions are that we see a bias for dispersions below 25 km s⁻¹, and that, although maybe the absolute values of the dispersions are biased, the relative values are still robust. We note that the typical Voronoi bin size is much smaller than the scale over which the rotation velocity varies on order of the dispersion value, and that therefore beam smearing effects are not important.

4 KINEMATICS OF THE DIG

4.1 Identification of the DIG

A critical step in the analysis is the identification of DIG. In edge-on galaxies, the identification is usually assumed to be the extraplanar gas. In less inclined galaxies, it is customary to either identify the star forming regions by the line ratios which differ from those for diffuse ionized gas, or by identifying peaks in the spatial distribution of emission line maps such as H α.

In this paper, we use two methods to identify the DIG in our maps. The first method we use was developed by Blanc et al. (2009) and has since been used also by other authors to distinguish between DIG and star forming gas (e.g. Kaplan et al. 2016; Kreckel et al. 2016). This method assumes that the surface brightness of H α is composed of a contribution from DIG and of H ii regions. In regions where the surface brightness is below a cut-off, f₀, to be discussed below, it is assumed that the DIG is completely dominant and the DIG fraction C_DIG = 1, while at higher surface brightnesses we can write: C_DIG = f₀/f_H α. The surface brightness f₀, which determines this fraction, is based on measurements of the [S ii]/H α ratio throughout a galaxy, which is known to be much higher for diffuse ionized gas than star forming regions (Madsen et al. 2006). For our data we first correct the [S ii] and H α lines for extinction by using the Balmer decrement, assuming an intrinsic decrement of 2.86. Then we sum the fluxes of the two [S ii] 6716 and 6731 Å lines together. We fit Formula 8 of Blanc et al. (2009) to the data using the MCMC code emcee (Foreman-Mackey et al. 2013), where we define the likelihood as the sum of the squared difference in the ratio divided by the squared errors on the ratio. Following Blanc et al., we adopt an intrinsic value of [S ii]/H α = 0.35 for DIG and [S ii]/H α = 0.11 for H ii regions, both of which we scale with a free parameter Z to account for metallicity differences between the MAD galaxies and the Milky Way (we refer the reader to the Appendix for results based on fits in which the metallicity was not free but adopted from Paper 2). To make the fit robust against outliers (see Fig. 4), we exclude the 5 per cent worst points from the likelihood.

In Fig. 4 we show for the galaxy NGC 4030 the [S ii]/H α ratio as a function of H α surface brightness, together with the fit shown in the solid red line. From this fit we decide how to divide the spaxels in the datacube into those where, statistically, H ii regions dominate (f_H α > 2f₀) and where DIG dominates (f_H α < 2f₀). The values of f₀ for the galaxies in the subsample are given in Table 3. The value of f₀ found for NGC 7162 is much lower than for the other galaxies. It is unclear if this is because the Blanc criterion breaks down for lower spatial resolution or if the intrinsic properties of this galaxy are simply different.

We also include a second identification method based on finding peaks in the H α emission, similar to what was done in Weilbacher et al. (2018). We run the publicly available astrodendro code.¹ Given an intensity map, this code splits an image into leafs, branches, and trunks (the leafs being the star forming regions in our maps). We run astrodendro with min_delta = 6.0 × 10⁻¹⁹ erg s⁻¹/cm⁻² and min_npix = 10. However, in order to minimize overlap between DIG and star forming regions we require for the dendrogram method that the DIG is at least 4 pixels away from the star forming regions.

In Fig. 5 we show for NGC 4030 which pixels are identified as SF gas, DIG, and neither for both methods. Similar figures can be found in Appendix F (online version only) for the five other galaxies in the kinematic subsample.

In Appendix F (online version only) we also compare the light and pixel fractions of DIG and SF gas. We find that the Blanc et al. method gives a relatively low fraction of pixel and light in DIG. The dendrogram method points at on average 30 per cent of the luminosity and 60 per cent of the pixels being in the form of DIG.

4.2 DIG velocity subsample

One of our goals is to derive rotation velocities for DIG and compare these with the velocities of SF gas. In order to carry this analysis out, we only look at galaxies for which the gas appears to be on circular orbits. For several galaxies, we see very complicated velocity structure close to their centres, for example the S-shapes seen in NGC 4941 and several other galaxies. In order to avoid very detailed modelling of the gas flows in these centres, we focus, for the velocity studies in this paper, on six galaxies without obvious kinematic twists and no large variation of systemic velocity of the gas with radius. The kinematic parameters of these six galaxies are summarized in Table 3. Although we tried to avoid galaxies with bars, NGC 7162 might be weakly barred (Buta et al. 2015).

4.3 Kinematic analysis

In order to derive the kinematics we need to know the location of the kinematic centre, the kinematic position angle, and the flattening (1 − b/a, with a and b the length of the major and minor axes) of the galaxy. Two of those quantities, the kinematic centre and the flattening, are often not easy to determine. The position of the kinematic centre is difficult to determine because the motion of the gas is often irregular close to the centre. The flattening is also difficult to infer from the kinematics, particularly for galaxy centres such as those that we are looking at.

We use kinemetry (Krajnović et al. 2006) to fit radially the best-fitting ellipse with free flattening and kinematic position angle. For the kinematic centre, we chose to adopt the photometric centre of the galaxy, which we determine on the white light image of the MUSE cube.

For the kinematic position angle, we use the median value of the position angle in the outermost 15 arcsec. The inclination is difficult to determine from our data; the radial profiles derived with kinemetry show significant scatter. We therefore use inclination values based on the axis ratios from the S4G photometry (Salo et al. 2015). These values are given in Table 3. We estimate the inclination i through cos (i) = q, with q = b/a. This assumes an infinitely thin disc; assuming an intrinsic thickness of q₀ = 0.2 leads to an increase in inclination of four degrees for the most flattened galaxy.

We measure the rotation velocity of the diffuse and star forming gas in elliptical annuli defined by the position angle and flattening of the galaxy. Each annulus is 3 arcsec wide, as a compromise between radially blurring the rotation signal and having sufficient bins to infer the rotation velocity for each of the components. We then calculate the angle of each bin along the ellipse and determine the amplitude (and uncertainty on the amplitude) of the sinus that best fits the velocities along the ellipse using lmfit (Newville et al. 2014), taking into account the uncertainties on the bin velocities. We do this for both the DIG bins and the star forming bins. For the velocity dispersions, we calculate the maximum likelihood estimate of the dispersion, without any angular dependence.

In order to compare with lower resolution data, we re-bin the six data continuum subtracted data cubes of the kinematic subsample to a sampling of 0|${^{\prime\prime}_{.}}$|6, which we subsequently smooth to a resolution of 2 arcsec. Contrary to before, we do not distinguish between DIG and SF gas for these reduced-resolution data, but instead measure the velocities and dispersions in a luminosity-weighted way.

5 RESULTS

In Figs 6 and 7 we show the rotation velocities and dispersions of star forming and diffuse ionized gas. For both DIG identification methods, the DIG seems to be, on average, lagging behind the star forming gas. This is expected if we are indeed observing extraplanar gas.

The difference between the velocity dispersion of diffuse ionized gas and star forming gas is more pronounced. (Note that the light-weighted low-resolution dispersion values can be lower than either measurement of the Voronoi binned gas as the convolution can spread out a single narrow peak over many pixels). We also note that the SF gas does not always agree with the rotation curves based on stellar kinematics. Given the simplistic derivation of the latter, one does not expect these curves to agree in detail; however, the SF gas rotating slower in most cases might point at the necessity of an asymmetric drift correction also for this component. The measured velocity dispersion, which is already corrected for instrumental broadening, consists of a thermal component and a gravitational component. As the thermal broadening is typically of the order of 10 km s⁻¹, most of the broadening must be due to gravitational broadening.

In order to quantify the difference in rotation velocity between the two components, we follow Levy et al. (2018) and take the median difference between the rotational velocities of the DIG and star forming gas as well as the median error on these differences. We note that taking a uncertainty-weighted average gives similar values. In Fig. 8 we show these median differences in rotational velocity between the diffuse ionized gas and the star forming gas for both DIG identification methods. For the dendrogram method, the velocity difference is consistently higher than zero, meaning that star forming gas rotates faster than DIG. For the surface brightness cut method, two galaxies (NGC 7162 and NGC 4592) do not show any noticeable difference. In Fig. 9, we compare the lags found with the dendrogram method to the lags observed between H α and CO measurements by Levy et al. Our difference measurements are on average 5 km s⁻¹ lower than between H α and CO. A Kolmogorov–Smirnov test suggests that the velocity differences in this work and those of Levy et al. are not drawn from the same sample, although with a relatively low significance level of α = 0.05. We note that the average differences between the integrated low-resolution measurement (black dash-dotted curves in Fig. 6) and the star forming gas, which are the two quantities most equivalent with the H α and CO measurements of Levy et al., are even smaller. Although it is possible that for the most-inclined galaxies the H α kinematics are affected by extinction and therefore miss the fastest rotating gas near the mid-plane, we see that (although low in number) also the less inclined galaxies in our sample show on average lower values.

Since the velocity fields are quite irregular, it is difficult to carry out the same analysis for the full sample. We therefore decided to only look at differences in dispersion between DIG and SF gas. To minimize the influence of bulges, bars, and AGN, we only look at radii outside 14 arcsec. We exclude NGC 5643, IC 2560, NGC 4941, NGC 3081, NGC 4593, and NGC 3393, NGC 1566, NGC 1097 which are known AGN, and for which the BPT diagram (see Paper 2) shows that most of our gas observations are dominated by the AGN.

In Fig. 10 we show how the median difference in velocity dispersion depends on the stellar mass and the axial ratio. We note that on differences in velocity dispersion are consistent with being zero or lower than zero, meaning that the SF has a equal or lower dispersion than the DIG. The uncertainty weighted sample average is Δσ = −6.9 ± 2.7 km s⁻¹.

If the gas we identified as DIG as indeed coming from extraplanar gas, one might expect a correlation between the mass of the galaxy and the dispersion of the dIG, since more massive galaxies have more massive discs (and higher surface mass density). From Fig. 10 we note that there is also no evidence for a trend with stellar mass; once removing the point with the highest difference in σ, the merger remnant NGC 3256, the dependence on mass looks flat.

We perform a similar check by looking at the differences in dispersion versus the axial ratio of the host galaxies. Assuming that all the galaxies are axially symmetric (which is not the case for NGC 3256), this is a proxy for the inclination of the host galaxy. Some DIG models (Marinacci et al. 2010) allow for anisotropic velocity dispersions. We do not see any evidence from our data that the dispersions of the DIG depend on the inclination of the galaxy.

6 DISCUSSION

6.1 Possible influences on the measurement

6.1.1 Continuum subtraction

We have shown in the previous sections how DIG can be seen to rotate slower than SF gas. We also showed that the strength of this rotational lag signal is dependent on the used method to identify DIG. Both methods however compare higher S/N data with lower S/N data. Although we do not believe that the Gaussian line fitting would lead to lower velocities for low S/N data, it is possible that during the fitting of the stellar continuum template mismatch alters the velocity of the H α line. As the stars generally rotate slightly slower and have a somewhat higher dispersion, underestimation of H α absorption in the stellar continuum could induce fake line emission at that location.

To show that this is not the case, we repeat the analysis for the six galaxies with DIG identification by H α surface brightness, but instead of fitting the velocity of the Balmer lines, we fit the velocity of the two [N ii] lines at 6548 and 6583 Å. These two lines are much less sensitive to continuum subtraction problems, as they are on opposite sides of H α, while they are still bright enough to accurately determine velocities. In Fig. 11 we show the rotation velocities of star forming and diffuse gas as determined from Balmer and the [N ii] lines for NGC 1084. We have checked the other five galaxies and find the same behaviour there. The [N ii] kinematics are different for SF gas and DIG, with the SF gas rotating faster than the DIG. This confirms our suspicion that the measurement is not affected by the continuum subtraction. For the gas identified as SF gas, the rotational velocities are lower for the [N ii] lines than for the Balmer lines. We think that this is because the relative contribution of the DIG to the observed flux is higher for [N ii], as the ratio of [N ii]/H α is known to enhanced for DIG (e.g. Haffner et al. 1999; Hoopes & Walterbos 2003; Madsen et al. 2006).

6.1.2 Extinction correction

We correct the H α fluxes for dust attenuation using the Balmer decrement, assuming an intrinsic ratio of 2.86 and a Fitzpatrick (1999) dust law, before identifying the DIG. It is possible that the light of a heavily attenuated H ii region mixed with the light of the DIG with different kinematics, will lead to an observed region with SF-like brightness but DIG-like kinematics. We have therefore checked if the results are significantly different if we would not correct for extinction. For this, we have re-calculated the f₀ values in the same way as in Section 4.1, but without extinction correction, and re-measured the kinematics for the new divisions between DIG an SF gas. Although we see some small quantitative differences, we do not see any qualitative differences in the results.

6.2 Most of the gas we identify as DIG is close to the mid-plane

In this section we explore the possible vertical distribution of the DIG in the six galaxies. We argue that both from a velocity point of view, as well as from a dispersion point of view, our data suggest a majority of the DIG emission that we are seeing is close to the mid-plane. For this, we assume that the Jeans equation holds, and that the gas and the potential follow a vertically exponential distribution.

1. For the two DIG identification methods used in this paper, the velocity lags are modest. Using the kinematic separation technique, we find velocity differences that are slightly higher than the other two methods. We use equation (B7) (see Appendix) as a measure to identify the scale height of the DIG normalized by the scale height of the potential h_z. We assume that the velocity of the star forming gas is equal to the rotation velocity in the mid-plane, and that the scale height of the potential is close to the scale height of the stars. A 10 per cent lag in velocity would then lead to a scale height that is 20 per cent of that of the stars (note that for a self-consistent potential this leads to the stars rotating at |$\frac{2}{3} V_0$|⁠, with V₀ the circular velocity in the mid-plane). Using a measurement for the lag in the velocity may actually in some instances be a more solid indicator of the spatial distribution of the gas than the velocity dispersion. For the latter, it is not always known what the contribution of non-gravitational motions is, whereas for a velocity measurement such components would most likely average out. However, we do acknowledge that there is still considerable uncertainty in interpreting these lags, as there is no guarantee that the gas is on circular orbits. We convert the velocity lags to scale heights for the velocity lags found by separating DIG and SF gas with the dendogram method. With the exception of two peaks in NGC 4592 and NGC 3521, the inferred scale heights are very moderate and always much lower than the scale height of the stars.

2. The dispersions can provide an independent measurement of the vertical scale height (equation B5). Using the median dispersion ratios between DIG and SF gas from the dendogram method, we find that dispersions of the DIG are 16 per cent (NGC 4030) to 35 per cent (NGC 1084) higher than that of the SF gas. At face value, this would lead to scale heights of the DIG that are similar fractions higher, and therefore only marginally higher than the scale height of the SF gas. However, it is not known what the contribution of non-gravitational motions is to the velocity dispersion of the gas, and therefore this fraction can be much higher. Using an estimate for the stellar scale height and stellar surface mass density in equation (B5), we can obtain an upper limit on the scale height. From the pPXF fits to the stellar continuum, we can obtain a crude estimate of the stellar surface mass density. The details of this procedure are outlined in Paper 2. The templates we use are based on a Kroupa IMF. We note that using a Salpeter IMF would lead to a higher surface mass density and therefore a lower scale height of the DIG.

In order to infer the vertical distribution of the stars, we make use of the relation between scale length and scale height of discs given by equation (1) in Bershady et al. (2010). The stellar scale heights for all six galaxies vary between 178 pc (for NGC 4790) and 371 pc (for NGC 4030). Solving equation (B5) with these values, and using the radial dispersion profiles for the DIG as identified from the dendrogram method, we find upper limits to the scale heights of the DIG. These are shown in Fig. 12. We see that the upper limits increase with radius. This is likely a consequence of using the stellar surface mass density as the dominant part of the potential. It is expected that the potential at larger radii will be more dominated by dark matter. As the scale heights for the velocity lags are based on a Mestel disc with a constant rotation curve, the scale heights based on velocity do not show radially rising profiles. Most important is however that under the assumption of an exponential vertical distribution almost all the upper limits on the DIG scale heights are consistent with the DIG originating from a layer that has a lower scale height than the stars.

We note that there is an alternative way of identifying DIG, which consists of performing double Gaussian component fits to emission lines. In Appendix  C, we explore this method and see that our conclusions of most of the DIG has a scale height lower than that of the stars is robust.

7 SUMMARY

In this paper we have presented the methodology for the kinematic extraction and the presentation of the kinematics maps of the first 41 galaxies of the MAD Survey with VLT/MUSE. We have outlined the methodology for the measurement of the gas kinematics, which we will use in future papers.

We use two methods, one based on the ratio of [S ii]/halpha (Blanc et al. 2009), the other on finding peaks in the H α maps, to identify diffuse ionized gas. For a subsample of six regularly rotating galaxies it is possible to measure differences in rotational velocity between ionized and star forming gas. We find median velocity lags of 0–10 km s⁻¹. These measurements are dependent on the method used to identify DIG.

Using the dendrogram method we measure differences in ionized gas velocity dispersions for DIG and SF gas for the rest of the sample after removing galaxies affected by an AGN. Although we caution against possible biases in the velocity dispersion due to the limited spectral resolution of our data, we do find a consistently higher velocity dispersion for DIG than for star forming gas. We do not find any dependence for this value on mass or inclination.

Interpreting the measured velocity lags and dispersions differences with Jeans models, we show that the observed diffuse gas in our data is consistently located close to the mid-plane of the galaxy, with scale heights that are below those of the stars. This limits the models for the origin of DIG gas to processes such as expanding shells, leaky H ii regions, and possibly ionizing radiation from older stars.

SUPPORTING INFORMATION

Appendix F. Identification of Dig.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 60.A-9100(C), 095.B-0532(A), 096.B-0309(A), 097.B-0165(A), 098.B-0551(A), 099.B-0242(B), 100.B-0116(A). We thank the ESO staff for their assistance during the observations. JB acknowledges support by FCT/MCTES through national funds by grant UID/FIS/04434/2019 and through Investigador FCT Contract No. IF/01654/2014/CP1215/CT0003. MdB thanks Kyriakos Flouris for sharing his reduction script and Sebastian Kamann for comments on the draft. This research made use of astropy, a community-developed core python package for Astronomy (Astropy Collaboration, 2013). This research made use of astrodendro, a python package to compute dendrograms of Astronomical data (http://www.dendrograms.org/). 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.

Footnotes

REFERENCES

APPENDIX A:

SIMULATIONS TO CHECK THE RECOVERY OF THE VELOCITY DISPERSION

For our analysis, we measure the position and the width of the Balmer lines and the [N ii] and [S ii] lines. However, most lines have widths that are well below the width of the line spread function of MUSE. This is generally not a problem to estimate the centroid of the line, however, it becomes increasingly difficult to estimate the dispersion correctly for low values of the dispersion. In order to assess how reliable the dispersion measurements are, we perform a small suit of simulations.

We focus exclusively on the Balmer lines, since we do not use the dispersion of the other lines. We generate H α and H β in a ratio of 3:1 on an spectrum that is oversampled 10 times. We assume idealized conditions to test the kinematics and therefore do not simulate the stellar continuum in the mock observations.

For each velocity dispersion we simulate 1000 spectra. We assume the emission line spectrum is redshifted by a velocity of 1600 km s⁻¹, which is typical of galaxies in the sample, and we additionally add a random velocity drawn from a Gaussian distribution with a width of 10 km s⁻¹. The width of the Gaussian is obtained by taking the squared sum of the line spread function and the dispersion. It is known that the LSF varies spatially. In order to capture the spatial variation in our simulations, we use a reduced cube of arc lines. As our data are the average of four different rotator angles, we use the same rotation angles to generate four arc cubes, which we than stack with small random Gaussian offsets (σ = 1 arcsec). We find that the spatial variation of the FWHM of the LSF at 6595 Å is 0.09 Å. At the wavelength of H α this corresponds to an uncertainty of σ = 2 km s⁻¹. This variation is higher than the one measured by Guérou et al. (2017), likely because they average many more observations. We include this uncertainty to our simulations by adding this as a Gaussian random variable to the LSF.

After generating the two emission lines, we rebin the spectrum and add Gaussian noise in such a way that the signal-to-noise is the same as what we would measure from the Voronoi binning of the cube on the H α line.

We then fit the lines in the same way as we fit the data, except that we skip the continuum subtraction part. These simulations are thus very much idealized and present an upper limit for how well we understand the measured velocities and dispersions.

Fig. A1 shows the measured velocity dispersion versus the input velocity dispersion and shows that we can recover unbiased dispersions down to about 35 km s⁻¹. The width of the error bars shows the 16th and 84th percentiles of the recovered dispersions. In the right-hand panel we show the bias in the recovery of the dispersion as a function of the input dispersion. The error bars in the figure show the 68 per cent confidence range and are normalized by the median uncertainty as determined from the fit to the simulated data at that input dispersion. The recovered uncertainties are close to 1 for high velocity dispersions, but become higher towards lower dispersions. This is consistent with the added uncertainty in the LSF.

APPENDIX B:

JEANS MODELS FOR EXTRAPLANAR GAS

We assume that the density distribution of the galaxy is dominated by the stellar mass in the plane of the galaxy, and that the vertical distribution has an exponential profile. We can then write the physical density at radius R and elevation above the disc z as |$\rho (R,z) = \frac{\Sigma (R)}{2 h_z} \exp \left(-z/h_z\right)$|⁠, where h_z is the scale height of the stars, and Σ(R) the surface mass density.

APPENDIX C:

TWO-COMPONENT FITS

For the kinematic subsample for which we determine the rotation velocity of the DIG (Section 4.2), we also fit double Gaussian component line profiles to the continuum subtracted spectra. To do this, we first tesselate these data to a high S/N of 50 at the location of H α. We assume that each line is composed of a narrow and a broad component, and that in this case the kinematics of each component of the H α line and lines from [N ii] and [S ii] are the same: for each voronoi bin the narrow components of all emission lines have the same velocity V₀ and velocity dispersion σ₀ and similar for the broad components (V₁, σ₁). The one-component fits were carried out with a damped χ² fitting method. For the two-component fits we chose to use an MCMC sampler (emcee; Foreman-Mackey et al. 2013), as such a sampler is more likely to determine the global minimum.

We acknowledge that such a decomposition is difficult at the spectral resolution of MUSE. We have removed the fits for which we did not see any convergence in the velocity dispersion of the second component, by removing all fits with δσ > 60 km s⁻¹. Although these decompositions are uncertain (not only are velocities and dispersions between different components degenerate, but the fits are also sensitive to small residuals from the sky subtraction), they do provide an interpretation of the data that is not necessarily less valid than a single component fit. The kinematics of these decompositions of the six galaxies are shown in Fig. C1.

It is interesting that in three of the six galaxies (NGC 4790, NGC 4592, and NGC 7162), we do not find evidence for gas at very high dispersions. It is possible that such a component is present but that our data do not have the spectral resolution or sensitivity to identify this. NGC 1084 shows higher dispersions in the second component in regions with little star formation. NGC 3521 shows a biconical structure in the second component, which resembles an outflow. The velocity maps of NGC 3521 and NGC 4030 clearly show lower rotation velocities for the second component and also a less pinched velocity field.

In Fig. C2 we show a quantitative comparison between the velocities and velocity dispersions of the narrow and broad components and the stellar rotation curves and stellar velocity dispersions. We derive the profiles by fitting the mean radial velocity/velocity dispersion profile of each component. As in particular for the broad component there are many outliers, we clip all values which are more than 5σ (standard error) away from the mean value. We then fit the values again (i.e. we do this only for one iteration). The narrow component is rotating with a velocity close to the circular velocity in the mid-plane, while the broad component is rotating slower. The dispersion values of the broad component are (per definition) higher than the narrow component, but, except for NGC 1084, do not seem to exceed the values found for the stars. Although differences between SF gas and DIG are more pronounced, these tentative results strengthen our conclusion that majority of the DIG that we identify in our data is coming from gas close to the mid-plane.

We additionally check for the presence of high-dispersion gas by stacking the continuum subtracted cubes. We use linear interpolation to convert the Voronoi binned spectra (of Section 3) to rest-frame wavelength using the H α velocity. We then carry out a similar two-component fit to the stacked spectrum of each galaxy. For the second component, we find dispersion values as summarized in Table C1. We note that these values are slightly higher than those found from the radial fits in Fig. C2. As the lags in velocity between the two component fits are sometimes of the order of 50 km s⁻¹, the stacking will artificially broaden either component by the same order of magnitude. Another reason for the discrepancy is the inclusion of the central 7 arcsec in the data when we stack the cube. When we exclude these central 7 arcsec we find values that are only slightly higher than the values in Fig. C2. Fits with three kinematic components did not converge.

APPENDIX D:

STELLAR ROTATION CURVES

We derive an estimate of the circular velocity in the mid-plane of the subsample of six galaxies by fitting Jeans Anisotropic Models (JAM) by Cappellari (2008) to the second moment of stellar kinematics. This second moment was derived by using |$v^2_{\mathrm{ RMS}} = \sigma ^2 + V^2_{\mathrm{ los}}$|⁠. The free parameters of these models are a multi-Gaussian expansion (MGE) of the mass distribution in the galaxy, an MGE expansion of the light distribution, a value for the orbital anisotropy β_z, and the inclination.

We fix the inclination to the value used in this paper based on the flattening of the galaxy light. We use the structural parameter fits of Salo et al. (2015) to derive the luminous distribution of the tracer population, by converting each morphological component into an MGE using the mge_fit_1d code (Cappellari 2014). These different luminous components were then scaled by an M/L ratio (per morphological component) and combined to form the MGE of the mass distribution. We allowed β_z to vary between −1.0 and 0.5. The JAM models allow for a convolution with a PSF, for which we used a single Gaussian with an FWHM of 0|${^{\prime\prime}_{.}}$|7, which is typical for the sample.

We use the best-fitting JAM models to obtain a mass distribution of the galaxy. Using the same deprojection as in the JAM models, we calculate the circular velocity in the z = 0 plane and multiply this by sin (i) to correct for the inclination. Given the crude modelling (e.g. the assumptions on the light profile, no proper dark matter halo), this velocity is more indicative than accurate.

APPENDIX E:

RESULTS FOR CONSTRAINED FITS TO f₀

In the text we followed Blanc et al. to determine the contribution of emission of the DIG to the observed H α surface brightness. We assumed that there was a value, f₀, below which the surface brightness of H α is completely dominated by the DIG. We determined a galaxy wide value of f₀ by fitting the observed line ratios of [S ii]/H α as a function of H α surface brightness, assuming a free scaling, assumed to be a metallicity, between the intrinsic values of [S ii]/H α for DIG and H ii regions in the Milky Way, and those in the galaxy we fitted. Here we perform a similar fit, but instead do not use a free parameter for the metallicity. We use instead the values determined in Paper 2 with the M13 method and assume a solar abundance of 12 + log  O/H = 8.69.

The f₀ values that we find for these more constrained fits are higher than the values found in the text using metallicity as a free parameter. The fits to the [S ii]/H α line ratios, with the metallicity fixed to the Paper 2 value, do not look as good as those with one more free parameter used in the main text of this Paper. For completeness, we do include the results of a repeated analysis with these somewhat higher values for f₀ in Figs E1 and E2.