The EDGE-CALIFA Survey: Molecular and Ionized Gas Kinematics in Nearby Galaxies

The EDGE-CALIFA survey (Bolatto et al. 2017) measured CO in 126 nearby galaxies with CARMA in the D and E configurations. Full details of the survey, data reduction, and masking techniques are discussed in Bolatto et al. (2017), but we present a brief overview here. The EDGE galaxies were selected from the CALIFA sample (discussed in the following section) based on their infrared (IR) brightness and are biased toward higher star formation rates (SFRs) (see Figure 6 of Bolatto et al. 2017). A pilot study of 177 galaxies was observed with the CARMA E-array. From this sample, 126 galaxies selected for CO brightness were re-observed in the D-array. These 126 galaxies with combined D and E array data constitute the main EDGE sample.¹⁵ The EDGE sample is the largest sample of galaxies with spatially resolved CO, at a typical angular resolution of 45 (corresponding to ∼1.5 kpc at the mean distance of the sample). Data cubes were produced with 20 km s⁻¹ velocity channels. At each pixel in the cube, a Gaussian is fit to the CO line. Velocity-integrated intensity, mean velocity, velocity dispersion, and associated error maps are created from the Gaussian fits. Pixels with velocities that differ from their nearest (non-blanked) neighbors by more than 40 km s⁻¹, generally caused by fitting failures in low signal-to-noise data, are replaced with the median value of the neighbors. This replacement is rare, occurring for only ∼0.5% of pixels in a given galaxy. Additional masking was applied to the CO maps where the Gaussian fitting introduced artifacts. This masking was based on a signal-to-noise ratio (S/N) cut using the integrated intensity and associated error map. Pixels with S/N < 1 were blanked in the velocity field. Average CO velocity dispersions are derived and are listed in Table 1. A beam smearing correction is applied and is discussed in Appendix B.

Table 1.  Parameters for the KSS

Name	Morph (Type)	log(${M}_{* }$)	log(SFR)	d	D25	CO ${V}_{\max }$	${\rm{\Delta }}{\text{}}V$	H i W50	H i W90	${{\rm{\Sigma }}}_{* }$	${\sigma }_{\mathrm{CO}}$	${\sigma }_{\mathrm{ADC}}$	${\sigma }_{{\rm{H}}\gamma }$	[S ii]/Hα	[N ii]/Hα
		(${M}_{\odot }$)	(${M}_{\odot }$ yr−1)	(Mpc)	(')	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(${M}_{\odot }$ pc−2)	(km s−1)	(km s−1)	(km s−1)
IC1199	Sbc (3.7)	10.8	0.2	68.3	1.2	199.0 ± 5.7	2.6 ± 4.3	141.8	229.0	189.6	11.1 ± 5.5	52.1 ± 26.5	34.6 ± 9.2	0.19 ± 0.01	0.39 ± 0.01
NGC2253	Sc (5.8)	10.8	0.5	51.2	1.4	174.2 ± 7.9	1.1 ± 7.7	211.0	243.9	205.6	9.2 ± 2.5	24.9 ± 37.2	28.0 ± 3.4	0.17 ± 0.01	0.37 ± 0.01
NGC2347	Sb (3.1)	11.0	0.5	63.7	1.6	285.6 ± 2.0	24.1 ± 0.4	267.7	289.2	336.5	10.1 ± 4.0	85.9 ± 5.0	36.6 ± 5.1	0.18 ± 0.01	0.38 ± 0.01
NGC2410	Sb (3.0)	11.0	0.5	67.5	2.2	227.2 ± 5.9	15.3 ± 5.1	247.3	267.3	198.6	17.6 ± 7.5	67.7 ± 16.5	38.2 ± 12.1	0.20 ± 0.01	0.47 ± 0.04
NGC3815	Sab (2.0)	10.5	0.0	53.6	1.4	185.8 ± 5.2	18.8 ± 3.6	176.5	223.5	215.4	15.1 ± 5.1	52.9 ± 8.0	26.9 ± 3.3	0.20 ± 0.02	0.38 ± 0.01
NGC4047	Sb (3.2)	10.9	0.6	49.1	1.5	216.3 ± 2.2	14.3 ± 0.5	220.3	243.2	296.6	9.5 ± 2.4	63.2 ± 9.5	29.5 ± 3.7	0.16 ± 0.01	0.35 ± 0.01
NGC4644	Sb (3.1)	10.7	0.1	71.6	1.5	186.1 ± 5.8	13.5 ± 2.8	171.0	214.7	132.6	16.6 ± 8.1	53.2 ± 10.9	31.3 ± 7.7	0.16 ± 0.01	0.40 ± 0.01
NGC4711	SBb (3.2)	10.6	0.1	58.8	1.2	142.8 ± 10.5	5.8 ± 0.9	167.3	185.2	121.1	10.0 ± 4.4	47.9 ± 37.2	32.9 ± 7.7	0.18 ± 0.01	0.36 ± 0.01
NGC5016	SABb (4.4)	10.5	−0.0	36.9	1.6	178.4 ± 0.8	17.7 ± 8.7	187.2	199.1	248.9	11.3 ± 5.6	56.4 ± 11.0	26.6 ± 3.8	0.16 ± 0.01	0.38 ± 0.01
NGC5480	Sc (5.0)	10.2	0.2	27.0	1.7	111.4 ± 6.6	4.8 ± 5.1	147.7	⋯	180.0	11.2 ± 3.5	37.1 ± 41.5	24.0 ± 2.9	0.20 ± 0.01	0.33 ± 0.01
NGC5520	Sb (3.1)	10.1	−0.1	26.7	1.6	162.3 ± 0.2	15.3 ± 1.4	158.1	170.1	202.3	13.4 ± 6.8	60.5 ± 0.4	26.5 ± 2.5	0.20 ± 0.01	0.39 ± 0.01
NGC5633	Sb (3.2)	10.4	0.2	33.4	1.1	187.5 ± 9.4	9.1 ± 2.9	200.6	⋯	396.0	11.3 ± 2.5	53.3 ± 23.4	25.1 ± 1.8	0.17 ± 0.01	0.37 ± 0.01
NGC5980	Sbc (4.4)	10.8	0.7	59.4	1.6	216.3 ± 4.4	8.2 ± 4.8	219.4	248.6	214.8	12.3 ± 2.3	52.5 ± 24.0	34.3 ± 3.8	0.17 ± 0.01	0.41 ± 0.01
UGC04132	Sbc (4.0)	10.9	1.0	75.4	1.2	238.5 ± 11.4	16.3 ± 8.0	255.6	291.5	206.5	15.8 ± 3.9	79.4 ± 40.0	33.1 ± 4.3	0.20 ± 0.01	0.40 ± 0.01
UGC05111	Sbc (4.0)	10.8	0.6	98.2	1.5	216.2 ± 11.9	9.1 ± 4.4	⋯	⋯	154.4	15.4 ± 3.6	52.2 ± 33.5	⋯	0.23 ± 0.02	0.41 ± 0.02
UGC09067	Sab (2.0)	11.0	0.7	114.5	0.8	211.7 ± 4.2	1.0 ± 2.1	212.7	⋯	110.7	14.1 ± 7.6	12.5 ± 40.7	29.0 ± 4.9	0.21 ± 0.01	0.35 ± 0.01
UGC10384	Sab (1.6)	10.3	0.7	71.8	1.2	187.4 ± 8.1	14.2 ± 7.0	180.7	200.9	74.5	18.9 ± 6.2	45.0 ± 12.9	35.4 ± 2.7	0.21 ± 0.01	0.35 ± 0.02

Note. The table lists relevant parameters for the KSS galaxies not already listed in Table 3. The morphology and types are from HyperLeda and are listed in Bolatto et al. (2017). Values for ${{\rm{M}}}_{* }$, SFR, distances (d), and the diameter of the 25th magnitude isophote (D₂₅) are taken from CALIFA and are also found in Bolatto et al. (2017). Errors on log(${M}_{* }$) and log(SFR) are ±0.1 in the corresponding units. CO ${V}_{\max }$ is the maximum CO rotation velocity, determined by the median of the CO RC at radii larger than twice the CO beam; the error is the standard deviation. Here, ${\rm{\Delta }}{\text{}}V$ is the median difference between CO and Hα RCs, as described in Section 4. H i W50 and W90 listed here are uncorrected for inclination. All values are taken from the GBT, except for NGC 5480, NGC 5633, and UGC 9067 (see Section 2.4). In this table, ${{\rm{\Sigma }}}_{* }$ are averages of stellar surface density radial profiles from (Utomo et al. 2017), ${\sigma }_{\mathrm{CO}}$ are the velocity dispersions estimated from the beam smearing corrected CO maps (Section 2.1), ${\sigma }_{\mathrm{ADC}}$ are the velocity dispersions estimated from the ADC (Section 5.5.2), and ${\sigma }_{{\rm{H}}\gamma }$ are the median velocity dispersions measured from the Hγ line; errors are the weighted standard deviations (Section 5.5.1). [S ii]/Hα and [N ii]/Hα are median intensity ratios, where the error is the standard error described in Section 5.6.

Download table as:  ASCIITypeset image

The CALIFA survey (Sánchez et al. 2012) observed 667 nearby (z = 0.005–0.03) galaxies. Full details of the CALIFA observations are presented in Walcher et al. (2014) and other CALIFA papers, but we present a brief overview for completeness. CALIFA used the PPAK IFU on the 3.5 m Calar Alto observatory with two spectral gratings. The low-resolution grating (V500) covered wavelengths from 3745 to 7500 Å with 6.0 Å (FWHM) spectral resolution, corresponding to an FWHM velocity resolution of 275 km s⁻¹ at Hα. The moderate-resolution grating (V1200) covered wavelengths from 3650 to 4840 Å with 2.3 Å (FWHM) spectral resolution, corresponding to an FWHM velocity resolution of 160 km s⁻¹ at Hγ (Sánchez et al. 2016a). The V500 grating includes many bright emission lines, including Hα, Hβ, Hγ, Hδ, the [N ii] doublet, and the [S ii] doublet. The V1200 grating contains many stellar absorption features used to derive the stellar kinematics as well as a few ionized gas emission lines, such as Hγ and Hδ. The typical spatial resolution of the CALIFA data are 25, corresponding to ∼0.8 kpc at the mean distance of the galaxies. The CALIFA galaxies were selected from the Sloan Digital Sky Survey (SDSS) DR7 to have angular isophotal diameters between 45'' and 792, in order to make the best use of the PPAK field of view. The upper redshift limit was set so that all targeted emission lines were observable for all galaxies in both spectral setups; the lower redshift limit was set so that the sample would not be dominated by dwarf galaxies. The data used for this study come from the final data release¹⁶ (Sánchez et al. 2016a) and data products come from  Pipe3D version 2.2 (Sánchez et al. 2016b, 2016c), provided in their final form by the CALIFA Collaboration.

The wavelength calibration of the data is detailed in Sánchez et al. (2012) and Appendix A.5 of Husemann et al. (2013) and is crucial to extract accurate line velocities. The wavelength calibration data are used to resample the spectra to a linear wavelength grid and to homogenize the spectral resolution across the band (6.0 Å FWHM for V500 and 2.3 Å FWHM for V1200). The calibration is done using HeHgCd lamp exposures before and after each pointing, using 16 lines for the V500 data and 11 lines for the V1200 data. The resulting accuracy of the wavelength calibration is ∼0.2–0.3 Å for the V500 data and ∼0.1–0.2 Å for the V1200 data. However, in our analysis of the V1200 data, we found errors in the wavelength calibration resulting from a bad line choice used to anchor the wavelength scale. This has been remedied in the current version of the data used here.

Once the data are calibrated,  Pipe3D fits and removes the stellar continuum, measures emission line fluxes, and produces two-dimensional data products for each emission line. Full details of  Pipe3D and its application to the CALIFA data can be found in Sánchez et al. (2016b, 2016c), and important details are reproduced here for completeness. The underlying stellar continuum is fit and subtracted to produce a continuum-subtracted or "emission line only" spectrum (Section 2 of Sánchez et al. 2016b). A Monte Carlo method is used to first determine the nonlinear stellar kinematic properties and dust attenuation at each pixel in the cube. Next, the results of this nonlinear fitting are fixed and the properties of the underlying stellar population are determined from a linear combination of simple stellar population (SSP) templates (see also Section 3.2 of Sánchez et al. 2016c). This model stellar spectrum is then subtracted from the CALIFA cube at each pixel to produce a continuum-subtracted cube. To determine the properties of the emission lines,  Pipe3D uses a nonparametric fitting routine optimized for weaker emission lines ("flux_elines") that extracts only the line flux intensity, velocity, velocity dispersion, and equivalent width (EW) (see Section 3.6 of Sánchez et al. (2016c) for full details). Each emission line of interest is fit using a moment analysis similar to optimal extraction. The line centroid is first guessed based on the rest-wavelength of the line, and a wavelength range is defined based on the input guess for the line FWHM. A set of 50 spectra in this range are generated using a Monte Carlo method and each is fit by a Gaussian. At each step in the Monte Carlo loop, the integrated flux of the line is determined by a weighted average, where the weights follow a Gaussian distribution centered on the observed line centroid and the input line FWHM. With the integrated flux fixed, the velocity of the line centroid is determined. The line fluxes are not corrected for extinction within  Pipe3D, so the desired extinction correction can be applied in the analysis. We do not apply an extinction correction, however, because it would have a minimal effect on the line centroid used here. Additional masking was also applied to the CALIFA velocity fields, using an S/N cut based on the integrated flux and error maps. Pixels with S/N <3.5 were blanked.

The linewidths of the V500 data (which covers Hα) are dominated by the instrumental linewidth (6.0 Å ≈ 275 km s⁻¹ at $\lambda ({\rm{H}}\alpha )=6562.68$ Å), and hence reliable velocity dispersions are not available for the V500 data. The instrumental linewidth can be removed from the V1200 data (2.3 Å ≈160 km s⁻¹ at $\lambda ({\rm{H}}\gamma )\,=4340.47$ Å). To determine the Hγ linewidth, we start with the continuum-subtracted cube and isolate the Hγ line. We fit the Hγ line at each spaxel using a Gaussian, where the linewidth is given by the width of the Gaussian fit. Pixels with S/N < 3 are blanked. We convert the resulting maps from wavelength to velocity using the relativistic convention, producing maps of the velocity dispersions for each galaxy. Independently,  Pipe3D does provide velocity dispersion maps derived from non-parametric fitting. The values in these maps, however, are frequently lower than the instrumental velocity dispersion over extended regions (a problem we do not find in our Gaussian fitting), and it is known that the pipeline systematically finds dispersions lower than obtained from Gaussian fitting (Section 3.6 of Sánchez et al. 2016c). We compare the velocity dispersions extracted from  Pipe3D and our Gaussian fitting to non-parametric fitting done with NEMO (Teuben 1995). In this fitting, we find the linewidth at each spaxel using the ccdmom mom=32 task. This finds the peak, locates the minima on either side of the peak, and takes a second moment over those channels. Velocity dispersions from this method agree much better with the Gaussian fitting results than with the  Pipe3D values, hence we adopt the Gaussian fitting results to determine the Hγ velocity dispersion. Before using these velocity dispersion maps in our analysis (Section 5.5.1), we remove the instrumental velocity dispersion, and model and remove the beam smearing effects (the latter is a small effect in the regions where we are interested in measuring the gas velocity dispersion). This procedure is discussed in Appendix B, and further caveats are discussed in Section 5.5.1. We regrid all CALIFA maps to the same grid as the corresponding EDGE map, using the Miriad task regrid (Sault et al. 1995).

CALIFA also derived effective radius (R_e) measurements for all EDGE galaxies as described in Sánchez et al. (2014). These values are listed in Table 4.

When comparing the velocity fields from the EDGE and CALIFA surveys, it is important to note that the velocities are derived using different velocity conventions: EDGE follows the radio convention, and CALIFA follows the optical convention. Because velocities in both surveys are referenced to zero, all velocities are converted to the relativistic velocity convention. In both the optical and radio conventions, the velocity scale is increasingly compressed at larger redshifts; typical systemic velocities in the EDGE-CALIFA sample are ∼4500 km s⁻¹. The relativistic convention does not suffer from this compression effect. Differences between these velocity conventions and conversions among them can be found in Appendix A. All velocities presented here are in the relativistic convention, unless otherwise noted.

In order to accurately compare the CO and ionized gas velocity fields, the EDGE and CALIFA data cubes were convolved to the same angular resolution. The convolution was done using the convol task in Miriad (Sault et al. 1995), which uses a Gaussian kernel. The EDGE beam was first circularized by convolving to a value 5% larger than the beam major axis. The CALIFA point spread functions are circular (Sánchez et al. 2016a). The EDGE and CALIFA cubes were convolved to a final 6'' resolution, corresponding to ∼2 kpc at the mean distance of the galaxies. Data products were reproduced as outlined in Sections 2.1 and 2.2. The CO and Hα velocity fields for NGC 2347 are shown in Figure 1. The rotation curves were derived as described in Section 3.1. There is excellent agreement between the native and convolved rotation curves for both CO and Hα, suggesting that, while it is best to match physical resolution, the convolution does not affect the results presented here.

Zoom In Zoom Out Reset image size

The EDGE collaboration obtained H i spectra for 112 EDGE galaxies from the Robert C. Byrd Green Bank Telescope (GBT) in the 2015B semester as part of GBT/15B-287 (PI: D. Utomo). We defer detailed discussion of these data to a future paper (T. Wong et al. 2018, in preparation) and here present only a brief overview. Observations were taken using the VEGAS spectrometer with a 100 MHz bandwidth, 3.1 kHz (0.65 km s⁻¹) spectral resolution, and a 3-σ sensitivity of 0.51 mJy. On-source integration time was 15 minutes for each galaxy. The GBT primary beam FWHM was 9', compared to the average EDGE D₂₅ = 16, so the galaxies are spatially unresolved. Data were reduced using standard parameters in the observatory-provided GBTIDL package. A first- or second-order baseline was fit to a range of line-free channels spanning 300–500 km s⁻¹ on either side of the signal range. The spectra were calibrated to a flux density scale, assuming a gain of 2 K Jy⁻¹ and a negligible coupling of the source size to the telescope beam. The widths containing 50% and 90% of the flux (W50 and W90 respectively) were derived from a Hanning smoothed spectrum to use as proxies for the maximum rotation velocity of the neutral atomic gas in these galaxies. These values are listed in Table 1, if available. The H i spectrum for NGC 2347 is shown in Figure 2, with the inclination-corrected W50 and W90 values marked.

Zoom In Zoom Out Reset image size

If H i data from GBT are not available, only W50 values were taken from Springob et al. (2005). Specifically, we use their W_C values, which are W50 corrected for the instrumental and redshift effects. For this work, these data are used for only three galaxies and come from either the Green Bank 300 ft telescope (NGC 5480 and NGC 5633) or Arecibo (line feed system, UGC 9067). These values are also listed in Table 1.

In this work, we use the H i rotation velocities as point of comparison to the CO and Hα rotation velocities. We convert the W50 and W90 values to rotation velocities where ${V}_{\mathrm{rot}}$=W $/(2\sin i)$, where i is the galaxy's inclination as listed in Table 3.

Rotation curves were determined for each galaxy using a tilted ring method (Rogstad et al. 1974; Begeman 1989), which has previously been applied to H i data (e.g Begeman 1989; Schoenmakers 1999; Fraternali et al. 2002; de Blok et al. 2008; Iorio et al. 2017), ionized gas data (e.g van de Ven & Fathi 2010; Di Teodoro et al. 2016), CO data (e.g Wong et al. 2004; Frank et al. 2016), and recently to $[{\rm{C}}\,{\rm{II}}]158\,\mu {\rm{m}}$ data in high-redshift galaxies (Jones et al. 2017). Galaxies were deprojected (position angles (PAs) and inclinations are listed in Table 3) and divided into circular annuli. The radius of each annulus was determined such that the width was at least half a beam. The center position, inclination (i), and PA are assumed to be the same for all annuli. The PA takes values between 0° and 360° and increases counterclockwise, where PA = 0 indicates that the approaching side is oriented due north. The rotation (${V}_{\mathrm{rot}}$), radial (${V}_{\mathrm{rad}}$), and systemic (${V}_{\mathrm{sys}}$) velocity components were determined in each ring through a first-order harmonic decomposition of the form

$$\begin{eqnarray}&&V(r)={V}_{\mathrm{rot}}(r)\cos \psi \sin i+{V}_{\mathrm{rad}}(r)\sin \psi \sin i+{\rm{\Delta }}{V}_{\mathrm{sys}}(r),\end{eqnarray} \tag{ 1 }$$

where r is the galactocentric radius and ψ is the azimuthal angle in the plane of the disk (Begeman 1989; Schoenmakers 1999). Before fitting, the central systemic velocity (${V}_{\mathrm{sys}}^{\mathrm{cen}}$) was subtracted from the entire map, such that the fitted systemic component is ${\rm{\Delta }}{V}_{\mathrm{sys}}(r)={V}_{\mathrm{sys}}(r)-{V}_{\mathrm{sys}}^{\mathrm{cen}}$.

The initial values for the PAs and inclinations were chosen from photometric fits to outer optical isophotes (Falcón-Barroso et al. 2017). If values were not available from this method, they were taken from the HyperLeda database (Makarov et al. 2014). Initial central systemic velocity values (${V}_{\mathrm{sys}}^{\mathrm{cen}}$) and center coordinates (R.A. and decl.) were taken from HyperLeda. The kinematic PAs were determined from the results of the ring fitting by minimizing ${V}_{\mathrm{rad}}$ at radii larger than twice the CO beam; an incorrect PA will produce a non-zero radial component. The ${V}_{\mathrm{sys}}^{\mathrm{cen}}$ values were refined by minimizing ${\rm{\Delta }}{V}_{\mathrm{sys}}$ at radii larger than twice the CO beam. The inclination is not as easily determined from kinematics; however, examining fits to individual annuli (rather than the rotation curve) can indicate whether the inclination is incorrect. Center offsets in R.A. and decl. (X_off, Y_off) were determined using a grid search method. At each point in the grid of X_off and Y_off values, a rotation curve was fit using that center. A constant was fit to the ${\rm{\Delta }}{V}_{\mathrm{sys}}$ component, and the combination of X_off and Y_off resulting in the best fit was selected as the center. The value of ${V}_{\mathrm{sys}}^{\mathrm{cen}}$ was then adjusted as necessary to again minimize ${\rm{\Delta }}{V}_{\mathrm{sys}}$. The sign of the offset is such that the correct center is (x_cen, y_cen) = (R.A.–X_off, decl.–Y_off). If a rotation curve could not be fit, either because there is little or no detected CO, or because the velocity field is very disturbed, the parameter values were unchanged from the initial values. The final values of the geometric parameters can be found in Table 3, including whether the PA, inclination, and ${V}_{\mathrm{sys}}^{\mathrm{cen}}$ values are derived from kinematics (this work), photometrically (Falcón-Barroso et al. 2017), or from HyperLeda.

Errors on the rotation curve were determined using a Monte Carlo method in which the geometrical parameters were drawn randomly from a uniform distribution. The center position was allowed to vary by 1'' in either direction because the average change in the CO (or Hα) center position from the original value is 07 over the whole EDGE sample. The inclination is varied by 2°, which is the average difference between the final and initial inclinations over the whole EDGE sample. The PA was also allowed to vary by 2°, which is the median difference between the final and initial PAs over the whole EDGE sample. This allows typical uncertainties in the kinematic parameters to be reflected in the rotation curves. The shaded error regions show the standard deviation of 1000 such rotation curves. The CO rotation curve showing ${V}_{\mathrm{rot}},{V}_{\mathrm{rad}},\mathrm{and}{\rm{\Delta }}{V}_{\mathrm{sys}}$ for one galaxy is shown in Figure 3(a).

Zoom In Zoom Out Reset image size

We note that this method to determine errors on the rotation curve differs from methods that use the differences between the approaching and receding sides of the galaxy, assuming that those differences are at the 2-σ level (e.g., Swaters 1999; de Blok et al. 2008). Typical uncertainties on the CO rotation velocity are ∼3–10 km s⁻¹ (1-σ). For the uncertainties stemming from the difference in rotation velocity between the approaching and receding sides to exceed the typical 1σ uncertainties in our Monte Carlo method, we find that rotation velocities of the approaching and receding sides would have to differ by 12–40 km s⁻¹. This seems unlikely, especially given that uncertainties on the rotation velocities derived from the differences between the approaching and receding sides presented by de Blok et al. (2008) are generally ∼10 km s⁻¹. Therefore, we conclude that the major sources of uncertainty in deriving our rotation curves are uncertainties in the geometric parameters.

Due to the beam size of the EDGE data, the observed velocities are affected by beam smearing, especially in the centers of the galaxies (Bosma 1978; Begeman 1987). Leung et al. (2018) analyzed the effect of beam smearing in the EDGE sample and found that it is only significant in the inner portions of the galaxy (≲0.5 R_e) where the velocity gradient is steep. For this study, we do not correct for beam smearing; instead, our analysis excludes points in the rotation curve within two beams from the center. The radius corresponding to twice the CO beam is referred to as ${R}_{2\mathrm{beam}}$ throughout. The excluded central region is shown in gray in Figure 3. Excluding the center of the galaxy also minimizes any effects from a bulge or an active galactic nucleus (AGN).

The CALIFA data were fit using the methods described in previous section, as well as the same PA and inclination as the CO listed in Table 3. In some cases, the Hα velocity contours are noticeably offset from the CO contours. CALIFA provides refinements to their astrometry in the headers of the data; however, these refinements are not large enough to account for some observed offsets. The CALIFA pipeline registers the R.A. and decl. for the center of the PPAK IFU to the corresponding center of the SDSS DR7 image (García-Benito et al. 2015). In DR2, 7% of the galaxies have registration offsets from SDSS >3'' (García-Benito et al. 2015). However, this registration process is known to fail in some cases. Indeed, in many of the galaxies for which we find offsets, this registration process has failed. Therefore, CALIFA centers were refit in the same way as the EDGE data, as described in Section 3.1. Because the V500 and V1200 data were taken on different days, the centers of the V500 and V1200 data were refit independently. For both the V500 and V1200 data, the average magnitude of the center offset is 09. The center offsets and ${V}_{\mathrm{sys}}^{\mathrm{cen}}$ values for the CALIFA data are presented in Table 4. The Hα rotation curve for NGC 2347 is shown in Figure 3(b).

Of the 126 EDGE galaxies, ≈100 have peak brightness temperatures ≥5σ (Bolatto et al. 2017). Reliable CO and Hα rotation curves could not, however, be derived for every detected EDGE galaxy. To best compare the CO and Hα rotation curves, a sub-sample of galaxies for which reliable CO and Hα rotation curves could be derived is used for the remainder of the analysis (the Kinematic Sub-Sample or KSS). A reliable rotation curve has small ${V}_{\mathrm{rad}}$ and ${\rm{\Delta }}{V}_{\mathrm{sys}}$ components at radii larger than ${R}_{2\mathrm{beam}}$ (as in Figure 3). In the centers of galaxies, there may be radial motions due to bars and other effects, but these should not affect the larger radii we consider here. Ensuring that both the CO and Hα have small ${V}_{\mathrm{rad}}$ components validates our assumptions that the CO and Hα have the same PA and inclination and that the PA and inclination do not change much over the disk (i.e., there are no twists or warps). In addition to the criteria on the rotation curves, there are four galaxies (NGC 4676A, NGC 6314, UGC 3973, and UGC 10205) for which the observed CO velocity width may not be fully contained in the band (Bolatto et al. 2017). One galaxy (UGC 10043) has a known Hα outflow (López-Cobá et al. 2017). These galaxies are also excluded from the subsample. Finally, we exclude galaxies with inclinations larger than 75°. At large inclinations, the line profiles can become skewed and a Gaussian fit to the line profiles is inappropriate and can lead to systematic biases in the mean velocities. We do, however, plan to analyze the highly inclined galaxies in a forthcoming paper (R. C. Levy et al. 2018, in preparation). Under these criteria, our sample size is reduced to 17 galaxies. Figure 18 shows CO and Hα velocity fields and rotation curves for all galaxies in the KSS. Specific notes on each galaxy in the KSS can be found in Appendix C. Table 1 lists global quantities for the KSS not listed in Table 3 or 4.

eDIG has been observed and discussed in the literature, so we highlight some results here for context. The importance of the warm ionized medium (WIM) as a significant fraction of the ISM in the Milky Way (MW) has been known for over four decades (e.g., Reynolds 1971; Reynolds et al. 1973; Kulkarni & Heiles 1987; Cox 1989; McKee 1990) and for over two decades in other galaxies (e.g., Dettmar 1990; Rand et al. 1990; Rand 1996; Hoopes et al. 1999; Rossa & Dettmar 2003a, 2003b). In particular, diffuse Hα can contribute >50% of the total Hα luminosity, with large variations (Haffner et al. 2009, and references therein). In the MW, half of the H ii is found more than 600 pc from the midplane (Reynolds 1993). How such a large fraction of diffuse gas can be ionized at these large scale heights is debated, but it is widely believed that leaky H ii regions containing O-star clusters can produce WIM-like conditions out to large distances from the cluster and the midplane (by taking advantage of chimneys and lines-of-sight with little neutral gas created by past feedback; see Reynolds et al. (2001) and Madsen et al. (2006)), although ionization sources with large penetration depths (such as cosmic rays) or re-accretion from the halo may also play a role.

Extraplanar H i exhibiting differential rotation has been detected in the MW (e.g., Levine et al. 2008) and in studies of individual galaxies (e.g., Swaters et al. 1997; Schaap et al. 2000; Chaves & Irwin 2001; Fraternali et al. 2002, 2005; Zschaechner et al. 2015; Zschaechner & Rand 2015; Vargas et al. 2017). Velocity gradients between the high-latitude gas and the dynamically cold midplane are generally −10–−30 km s⁻¹ kpc⁻¹, extending out to a few kpc, but there are large variations among and within individual galaxies.

Extraplanar Hα (i.e., eDIG) has also been found and studied in galaxies, primarily from photometry. In a recent study tracing the WIM in the spiral arms of the MW, Krishnarao et al. (2017) find an offset between the CO and Hα velocity centroids, although they do not interpret this offset as eDIG. Outside the MW, NGC 891 is the prototypical galaxy with bright eDIG extending up to 5.5 kpc from the midplane (Rand et al. 1990; Rand 1997) and a vertical velocity gradient of −15 km s⁻¹ kpc⁻¹ (Heald et al. 2006a), in agreement with measurements of its extraplanar H i (Swaters et al. 1997; Fraternali et al. 2005). Boettcher et al. (2016) find that the thermal and turbulent velocity dispersions (11 km s⁻¹ and 25 km s⁻¹, respectively) are insufficient to support eDIG in hydrostatic equilibrium with a 1 kpc scale height. NGC 5775 has observed H i loops and filaments with their rotation lagging the midplane (Lee et al. 2001), as well as Hα lags of −8 km s⁻¹ kpc⁻¹ detected up to 6–9 kpc from the midplane (Rand 2000; Heald et al. 2006a). NGC 2403 has eDIG that lags the midplane by 80 km s⁻¹, extending a few kpc above the midplane (Fraternali et al. 2004), in rough agreement with the lags observed in H i, which extend 1–3 kpc from the midplane (Fraternali et al. 2002). Finally, eDIG has been observed in the face-on galaxy M 83 with a lag relative to the midplane of 70 km s⁻¹ and a vertical scale height of 1 kpc (Boettcher et al. 2017). There is a range of eDIG velocity gradients and scale heights; moreover, H i and Hα vertical velocity gradients are not always similar or present (e.g., Zschaechner et al. 2015).

Apart from these case studies, there are several large photometric studies of eDIG independent of H i. Following the work of Rand (1996), Miller & Veilleux (2003a) and Rossa & Dettmar (2003a, 2003b) performed larger photometric surveys of nearby edge-on spiral galaxies. Rossa & Dettmar (2003a, 2003b) had a sample of 74 edge-on disks and found that 40.5% of the sample had eDIG extending 1–2 kpc from the midplane. In their sample of 17 galaxies, Miller & Veilleux (2003a) observe eDIG in all but one galaxy. Miller & Veilleux (2003b) did a spectroscopic follow-up study of nine edge-on galaxies with observed eDIG and found vertical velocity gradients ranging from −30 to −70 km s⁻¹ kpc⁻¹.

The stellar circular velocity curve, which accounts for stellar velocity dispersion, should agree with the CO rotation curve if CO is a dynamically cold tracer. Leung et al. (2018) test three different dynamical models of the galaxy's potential to determine stellar circular velocity curves and compare these to CO rotation curves of 54 EDGE galaxies. Overall, they find agreement between the CO rotation curves and the three models to within 10% at 1 R_e. We defer to Leung et al. (2018) for a complete discussion of the details of the stellar modeling. The agreement between the stellar dynamical modeling and the CO rotation curves verifies that CO is indeed a dynamically cold tracer, indicating that the Hα is exhibiting anomalous behavior rather than the CO. Moreover, our measured CO velocity dispersions (Table 1) are small (∼10 km s⁻¹), further indicating that the CO is dynamically cold.

It is possible that the observed difference between the CO and ionized gas rotation velocities could be produced by the inclination of the disk; however, we find no correlation between ${\rm{\Delta }}{\text{}}V$ and inclination, as shown in Figure 8. To determine a correlation, we calculate the Spearman rank correlation coefficient (r_s), which quantifies how well the relationship between the variables can be described by any monotonic function. Variables that are perfectly monotonically correlated will have r_s = ±1, assuming that there are no repeated values of either variable. The Spearman rank correlation coefficient does not, however, take the errors on the data into consideration. The errors on ${\rm{\Delta }}{\text{}}V$ are especially important here. Therefore, we use a Monte Carlo method to determine the correlation coefficient over 1000 samples drawn from a uniform random distribution within the error ranges on each point. The error reported on r_s is the standard deviation of all 1000 r_s values. As shown in Figure 8, the difference between the CO and Hα rotation velocities is not a result of the inclination of the galaxy (r_s = −0.01 ± 0.11).

Zoom In Zoom Out Reset image size

We note that, if the molecular and ionized gas disks had different inclinations, our assumption that they are the same could produce a ${\rm{\Delta }}{\text{}}V$. However, to produce only ${\rm{\Delta }}{\text{}}V\gt 0$ would require that all ionized gas disks are less highly inclined (with respect to us) than the molecular gas disks, which is extremely unlikely for a sample of 17 galaxies. We can, therefore, rule out that the inclination affects the results in this way.

We also plotted ${\rm{\Delta }}{\text{}}V$ against other global parameters of the galaxies, such as stellar mass (M_*), SFR, specific SFR (sSFR ≡ SFR/M_*), morphology, physical resolution, and CO V_max. There are no trends with any of these global parameters, but they are shown for completeness in Figure 17 of Appendix C. As mentioned in Section 2.3, the lack of trend with physical resolution (Figure 17(e)) justifies our choice to convolve to a common angular resolution rather than to a common physical resolution. These values and their sources are listed in Table 1. The lack of trends with these parameters agrees with Rossa & Dettmar (2003b), who also found no trends in the presence of eDIG with such global parameters.

In previous studies of eDIG, Rand (1996) and (Rossa & Dettmar 2003a) found a possible trend in the amount of eDIG with the SFR per unit area (Σ_SFR) as traced by the far-infrared (FIR) luminosity (L_FIR). The physical picture is that a minimum level of widespread star formation is needed to sustain a thick disk that covers the entire plane of the galaxy. Rossa & Dettmar (2003a) determine a threshold Σ_SFR above which they claim that an eDIG will be ubiquitous. This does not guarantee, however, that galaxies above this threshold will always have an eDIG or that galaxies below it cannot have eDIG. Rossa & Dettmar (2003a) define this threshold Σ_SFR as

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{FIR}}}{{D}_{25}^{2}}=(3.2\pm 0.5)\times {10}^{40}{\mathrm{erg}{\rm{s}}}^{-1}\,{\mathrm{kpc}}^{-2},\end{eqnarray} \tag{ 2 }$$

where D₂₅ is diameter of the 25th magnitude isophote. Catalán-Torrecilla et al. (2015) measure total IR (TIR, 8–1000 μm) luminosities (L_TIR) for 272 CALIFA galaxies, and L_TIR measurements are available for 15/17 KSS galaxies. The threshold defined by Rossa & Dettmar (2003a) can be converted to TIR by multiplying by 1.6 (Sanders & Mirabel 1996), such that

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{TIR}}}{{D}_{25}^{2}}=(5.1\pm 0.8)\times {10}^{40}\mathrm{erg}\ {{\rm{s}}}^{-1}\,{\mathrm{kpc}}^{-2}.\end{eqnarray} \tag{ 3 }$$

For the two galaxies without L_TIR measurements, we can estimate L_TIR from the SFR measured by CALIFA from extinction-corrected Hα, where

$$\begin{eqnarray}&&{L}_{\mathrm{TIR}}=\displaystyle \frac{1.6}{4.5\times {10}^{-44}}\left[\displaystyle \frac{\mathrm{SFR}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right]\,\mathrm{erg}\,{{\rm{s}}}^{-1}\end{eqnarray} \tag{ 4 }$$

(Kennicutt 1998) and the factor of 1.6 comes from converting from L_FIR to L_TIR (Sanders & Mirabel 1996). Using measurements of D₂₅ from HyperLeda (values are listed in Table 1), we compare the values of ${L}_{\mathrm{TIR}}/{D}_{25}^{2}$ for all KSS galaxies to the eDIG threshold (Equation (3)) in Figure 9. We find that ${94}_{-0}^{+6} \%$ of galaxies in the KSS have ${L}_{\mathrm{TIR}}/{D}_{25}^{2}$ greater than this threshold and should have eDIG based on this criterion.

Zoom In Zoom Out Reset image size

If many galaxies have a thick ionized disk, this could underestimate the dynamical mass of the galaxy derived from the ionized gas rotation velocity. Because we find that the ionized gas rotates more slowly than the molecular gas in galaxies with large Σ_SFR, this effect could be significant in local star-forming galaxies, and even more so at higher redshifts where there is more star formation occurring on average.

The trends between the eDIG and the SFR surface density discussed in Section 5.4 (e.g., Rand 1996; Rossa & Dettmar 2003a) are suggestive of star formation feedback playing an important role in forming the eDIG. In order for ionized gas to remain above or below the disk midplane in a long-lived configuration, it must have sufficient velocity dispersion (or at least a vertical bulk motion). This effectively acts as as an additional pressure term, allowing the gas to remain at larger scale heights. Therefore, we expect that galaxies with larger ionized gas velocity dispersions should have larger eDIG scale heights. Measuring the ionized gas velocity dispersion is, therefore, an important way to test these ideas.

5.5.1.  $H\gamma$ Velocity Dispersion Measurements

CALIFA data cannot, unfortunately, be used to accurately measure the Hα linewidth due to the low spectral resolution of the V500 grating employed (6.0 Å FWHM ≈ 275 km s⁻¹ at Hα). CALIFA observes Hγ with the moderate-resolution V1200 grating (2.3 Å FWHM ≈ 160 km s⁻¹ at Hγ), however, and those data can be used, in principle, to establish ionized gas velocity dispersions. Starting with the continuum-subtracted cubes, we fit the Hγ line with a Gaussian at each spaxel where the width of the Gaussian corresponds to the velocity dispersion, as discussed in Section 2.2. The measured linewidth is the convolution of the instrumental response with the actual gas velocity dispersion; with a small contribution from rotation smearing caused by the finite angular resolution. We apply a beam smearing correction that also accounts for the instrumental velocity dispersion. This method is described in detail in Appendix B. The accuracy of the resulting ionized gas velocity dispersion depends critically on our exact knowledge of the spectral resolution and response of the grating, because the instrumental velocity dispersion (σ_inst ≈ 68 km s⁻¹) is of the same order as the observed σ_Hγ before removal: in other words, the spectral resolution of the V1200 observation is marginal for the purposes of measuring the velocity dispersion in these galaxies, and our results should be considered tentative. Inspection of the maps suggests that it is likely that the beam smearing corrected σ_Hγ values reported here are lower limits to the real ionized gas velocity dispersion, and we caution against over-interpretation of these values. We do not correct our σ_Hγ for inclination, and hence assume that the velocity dispersion is isotropic. As discussed in Section 4, anisotropic velocity dispersions may be required to maintain ionized gas disks with vertical gradients in the rotation velocity (Marinacci et al. 2010). With these caveats in mind, we calculate the average beam smearing corrected Hγ velocity dispersion over the same region where ${\rm{\Delta }}{\text{}}V$ is calculated. Velocity dispersions range from ∼25–45 km s⁻¹ (Figure 10(a)). The scale height of the disk (h) corresponding to a given isotropic velocity dispersion (σ) is

$$\begin{eqnarray}&&h=\displaystyle \frac{{\sigma }^{2}}{\pi G{{\rm{\Sigma }}}_{* }(r)}\end{eqnarray} \tag{ 5 }$$

(van der Kruit 1988; Burkert et al. 2010). We use azimuthally averaged radial profiles of the stellar surface density (Σ_*) from Utomo et al. (2017) to find Σ_* over the same range of radii where ${\rm{\Delta }}{\text{}}V$ is calculated (listed in Table 1). The scale heights that correspond to the observed σ_Hγ are ∼0.1–1.3 kpc. Previous measurements of eDIG scale heights range from 1 to 2 kpc (Miller & Veilleux 2003b; Rossa & Dettmar 2003a; Fraternali et al. 2004) up to a few kpc above the disk (Rand 2000, 1997; Miller & Veilleux 2003a). Because our σ_Hγ are likely lower limits, the scale heights may indeed be larger than we report here.

Zoom In Zoom Out Reset image size

5.5.2. Velocity Dispersion Estimates from an Asymmetric Drift Correction (ADC)

It is possible to infer the velocity dispersion needed to produce the observed ${\rm{\Delta }}{\text{}}V$ using an ADC. Generally, an ADC is used to find V_circ given V_rot, σ, and Σ_*(r); however, because ${V}_{\mathrm{rot}}(\mathrm{CO})$ traces V_circ (Section 5.2), we can invert the ADC to find σ instead, with ${V}_{\mathrm{rot}}={V}_{\mathrm{rot}}$(Hα). If we assume that the velocity dispersion is isotropic (σ_r = σ_z = σ_ϕ), σ(r) =  constant, and Σ_*(r) = 2ρ(r, z)h(z) (Binney & Tremaine 2008), then

$$\begin{eqnarray}&&{\sigma }^{2}=\displaystyle \frac{{V}_{\mathrm{circ}}^{2}-{V}_{\mathrm{rot}}^{2}}{-d\mathrm{ln}{{\rm{\Sigma }}}_{* }/d\mathrm{ln}r}.\end{eqnarray} \tag{ 6 }$$

We use azimuthally averaged radial profiles for Σ_*(r) from Utomo et al. (2017) to find $d\mathrm{ln}{{\rm{\Sigma }}}_{* }/d\mathrm{ln}r$. We average V_circ, V_rot, and Σ_*(r) over the same radii as those where ${\rm{\Delta }}{\text{}}V$ is calculated; this excludes the central two beams (12''∼4 kpc) where beam smearing or a bulge can affect the rotation curve. Velocity dispersions from the ADC method (σ_ADC, Equation (6)) range from ∼15–85 km s⁻¹ in the KSS (Figure 10(b)). There is an apparent trend with ${\rm{\Delta }}{\text{}}V$ resulting from Equation (6). Using Equation (5), we find scale heights ranging from ∼0.1 to 2.0 kpc, again in rough agreement with previous measurements. For individual galaxies, the velocity dispersions and scale heights predicted from the ADC tend to be larger than those measured in the Hγ (Figure 11), but given the difficulty in the measurement, the agreement is reasonable.

Zoom In Zoom Out Reset image size

5.5.3. Trends between ΔV and the Velocity Dispersion

The velocity dispersion is not the only tracer of eDIG. The ratios of [N ii]λ6583/Hα ([N ii]/Hα) and [S ii]λ6717/Hα ([S ii]/Hα) increase with distance from the midplane and are used to probe the ionization conditions of the WIM (e.g., Miller & Veilleux 2003a, 2003b; Fraternali et al. 2004; Haffner et al. 2009). Although [S ii]/Hα varies only slightly with temperature, [N ii]/Hα is used to trace variations in the excitation temperature of the gas (Haffner et al. 2009). From observations of the MW and a few other galaxies, [S ii]/Hα = 0.11 ± 0.03 and [N ii]/Hα ∼ 0.25 in the midplane (Madsen 2004; Madsen et al. 2006), whereas [S ii]/Hα = 0.34 ± 0.13 and [N ii]/Hα ≳ 0.5 in the eDIG (Madsen 2004; Blanc et al. 2009). Observations of these ratios indicate that there must be additional heating sources aside from photoionization from leaky H ii regions to produce the WIM (Haffner et al. 2009, and references therein).

CALIFA provides Hα, [S ii], and [N ii] intensity maps for all galaxies. These were masked to cover the same radii as those where ${\rm{\Delta }}{\text{}}V$ is calculated. As shown in Figure 12, there are no trends between [S ii]/Hα or [N ii]/Hα with ${\rm{\Delta }}{\text{}}V$ (r_s = 0.03 ± 0.18 and r_s = − 0.02 ± 0.17, respectively). For both [S ii]/Hα and [N ii]/Hα, our ratios are all larger than for the plane of the MW, but only a few fall within the observed range for the eDIG. This is perhaps not unexpected, because emission from the plane is mixed with emission from the eDIG that would systematically lower the observed ratios. This is an encouraging hint that the observed ${\rm{\Delta }}{\text{}}V$ could be due to eDIG in a thick disk.

Zoom In Zoom Out Reset image size

To further investigate how the disk's geometry affects the observed rotation curve, we perform a suite of kinematic simulations using NEMO. Disks are given different scale heights and vertical rotation velocity distributions, as described in the follow subsections. The particle velocities are given by an input rotation curve that rises linearly from r = 0–1 units and is constant at V₀ = 200 km s⁻¹ from r = 1–6 units. The disk is then inclined and "observed" with a one-unit beam. The velocity is derived by fitting the peak of the line at each point in the simulated data cube. We then average ${V}_{\mathrm{rot}}$(r) between r = 1 and r = 5, in order to give $\bar{V}$ in the flat part of the rotation curve. The error on $\bar{V}$ (${\sigma }_{\bar{V}}$) is the standard deviation of ${V}_{\mathrm{rot}}$(r). A simulated ${\rm{\Delta }}{\text{}}V$ is computed by ${V}_{0}-\bar{V}$, which is analogous to the ${\rm{\Delta }}{\text{}}V$ defined previously (V₀ corresponds to ${V}_{\mathrm{rot}}$(CO) and $\bar{V}$ corresponds to the median ${V}_{\mathrm{rot}}({\rm{H}}\alpha )$ in the outer part of the galaxy). We test four disk configurations, which are described below. To convert the scale heights to physical units (i.e., kpc), we use the turn-over radius (R₀) of the rotation curve (one unit in the simulations) and find the average R₀ of the KSS Hα data. We fit the KSS Hα rotation curves with ${V}_{\mathrm{model}}={V}_{0}(1-{e}^{-R/{R}_{0}})$ (e.g., Boissier et al. 2003; Leroy et al. 2008) and fix V₀ to the value determined for that galaxy. The average R₀ over the sample is found, and the scaling is one unit = 1.77 ± 0.25 kpc. We note that the assumed beam in the simulations is one unit, which is nearly identical to the beam size of the observations (6'' = 1.73 kpc at the average distance of the KSS galaxies).

5.7.1. Thin Disks

First, we create a thin disk of particles with scale height h = 0. The simulated ${\rm{\Delta }}{\text{}}V$ as a function of inclination is shown in Figure 13(a). The resulting ${\rm{\Delta }}{\text{}}V$ values are all very small and are insufficient to explain the offsets seen in Figure 5. We recover the input rotation curve to within ∼2 km s⁻¹. This ∼2 km s⁻¹ offset is due to beam smearing. If a smaller beam is used, this offset disappears. Ergo, neither the inclination nor beam smearing of a thin disk can produce the observed ${\rm{\Delta }}{\text{}}V$.

Zoom In Zoom Out Reset image size

5.7.2. Thick Disks

5.7.3. Thick Disks with Vertical Rotation Velocity Gradients

We also test a thick disk with a vertical gradient in the rotation velocity. As mentioned in Section 5.1, vertical gradients in the rotation velocity have been observed in the extraplanar H i and eDIG of several galaxies. The physical rationale behind a vertical gradient in the rotation velocity is related to turbulent pressure support. Gas with larger velocity dispersions can be in pressure equilibrium at larger distances from the disk midplane. Material off the plane should have an orbit inclined with respect to the main disk enclosing the galactic center (like the stars); however, pressure support forces the gas to orbit parallel to the main disk. The gas farther off the plane rotates more slowly than gas closer to the plane, creating the vertical gradient in the rotation velocity. First, we impose a linear vertical rotation velocity gradient parameterized by η, where

$$\begin{eqnarray}&&{V}_{\mathrm{rot}}(z)={V}_{0}\left(1-\eta \displaystyle \frac{z}{h}\right).\end{eqnarray} \tag{ 7 }$$

As shown in Figure 13(c), there is a strong trend between η and ${\rm{\Delta }}{\text{}}V$. A linear vertical gradient in the rotation velocity can produce the observed values of ${\rm{\Delta }}{\text{}}V$ for η ≲ 0.3, meaning that at z = h, ${V}_{\mathrm{rot}}\gtrsim 0.7{V}_{0}$. For a given value of η, larger scale heights have less of an effect on ${\rm{\Delta }}{\text{}}V$, as seen from their shallower slope in Figure 13(c).

Next, we test a more realistic model for ${V}_{\mathrm{rot}}(z)$ than the linear vertical velocity gradient. The rotation velocity of material above the disk is governed by the potential. Here, we assume that the rotation velocity in the disk midplane is constant (V₀), hence the radial surface density profile, Σ(r), is described by a Mestel disk (Mestel 1963; Binney & Tremaine 2008). We also assume that the vertical density distribution is exponential. Therefore, the total density distribution of material that dominates the potential has the form

$$\begin{eqnarray}&&\rho (r,z)=\displaystyle \frac{{V}_{0}^{2}}{2\pi {{Grh}}_{p}}{e}^{-| z| /{h}_{p}},\end{eqnarray} \tag{ 8 }$$

where h_p is the vertical scale height of the material that dominates the potential. The potential of a thin Mestel disk is

$$\begin{eqnarray}&&\phi (r)=-{V}_{0}^{2}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{\max }}\right),\end{eqnarray} \tag{ 9 }$$

where r_max is the maximum extent of the disk (Binney & Tremaine 2008). Therefore, the total potential is

$$\begin{eqnarray}\phi (r,z) & = & \displaystyle \int \phi (r,z-z^{\prime} )\xi (z^{\prime} ){dz}^{\prime} \\ & = & -{V}_{0}^{2}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{\max }}\right)\displaystyle \frac{1}{{h}_{p}}\displaystyle \int {e}^{-| z^{\prime} | /{h}_{p}}{dz}^{\prime} \\ & = & {V}_{0}^{2}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{\max }}\right){e}^{-| z| /{h}_{p}}+c,\end{eqnarray} \tag{ 10 }$$

where the constant of integration (c) can be found by demanding that $\phi (r,z\to \infty )\to 0$. Therefore,

$$\begin{eqnarray}&&\phi (r,z)={V}_{0}^{2}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{\max }}\right){e}^{-| z| /{h}_{p}}.\end{eqnarray} \tag{ 11 }$$

The rotation velocity from a potential is given by ${V}_{\mathrm{rot}}^{2}=r\tfrac{\partial }{\partial r}\phi (r,z)$ (Binney & Tremaine 2008). Thus, from Equation (11),

$$\begin{eqnarray}&&{V}_{\mathrm{rot}}(r,z)={V}_{0}\sqrt{{e}^{-| z| /{h}_{p}}},\end{eqnarray} \tag{ 12 }$$

where the absolute value preserves the symmetry above and below the disk. The eDIG has its own density distribution and scale height (h_eDIG), which are largely independent of the potential and determined mostly by the star formation activity. Although large ratios of h_eDIG/h_p are physically unlikely, here we treat these two scale heights as independent quantities. We can find the eDIG density as a function of z where the vertical density distribution is described by the hydrostatic Spitzer solution (Spitzer 1942; Binney & Tremaine 2008; Burkert et al. 2010):

$$\begin{eqnarray}&&{\rho }_{\mathrm{eDIG}}(r,z)={\rho }_{0}{{\rm{sech}} }^{2}\left(\displaystyle \frac{z}{{h}_{\mathrm{eDIG}}}\right),\end{eqnarray} \tag{ 13 }$$

where ρ₀ is the density in the midplane. We repeat the NEMO simulations as before, but using ${V}_{\mathrm{rot}}$(r, z) given by Equation (12) (rather than Equation (7)) and an eDIG density distribution give by Equation (13) (rather than a Gaussian). The results of this more realistic model are shown in Figure 13(d). There is a strong trend between h_eDIG and ${\rm{\Delta }}{\text{}}V$. There is also a secondary trend between h_p and ${\rm{\Delta }}{\text{}}V$. With this model, it is possible to reproduce the observed range of ${\rm{\Delta }}{\text{}}V$ with h_eDIG ≲ 1.5 kpc. Using Equation (5) and the median Σ_* of the KSS at the radius where the measurements are done (≈200 M_⊙ pc⁻²), this implies a velocity dispersion ≲60 km s⁻¹, which agrees with the range of velocity dispersions inferred from the ADC and the measurements from Hγ.

In the previous subsections, we expounded our hypothesis that the observed difference between the molecular and ionized gas rotation velocities is due to eDIG in a thick disk with a vertical gradient in the rotation velocity. That is not, however, the only explanation. It is possible that we are instead measuring ionized gas velocities and velocity dispersions in the galactic bulge. To test this, we explore potential correlations between ${\rm{\Delta }}{\text{}}V$ and the measured bulge-to-disk (B/D) luminosity ratios (Méndez-Abreu et al. 2017). Here, we use only results from the r-band because it overlaps with Hα. This is shown in the top panel of Figure 14(a). To determine whether there are any underlying correlations, we follow the same methodology as was used for the velocity dispersions discussed in Section 5.5.3. There is one galaxy, with a large B/D ratio, that is a clear outlier from the rest; this galaxy is NGC 2347, which has the largest ${\rm{\Delta }}{\text{}}V$ in the subsample. It is excluded from the fitting of the linear regression. The bottom panel of Figure 14(a) shows the distance of each galaxy from the best-fit line. For the r-band, 33% of galaxies are consistent with the best-fit line within 1σ, and 67% are consistent within 3σ. It is possible that the bulge could be affecting the measured rotation velocities and contributing to the measured ${\rm{\Delta }}{\text{}}V$. This affect is likely minimal, however, as ${\rm{\Delta }}{\text{}}V$ is only measured at radii larger than ${R}_{2\mathrm{beam}}$ (12'', 4 kpc). Méndez-Abreu et al. (2017) also measure the effective bulge radius (R_bulge). We plot R_bulge versus ${\rm{\Delta }}{\text{}}V$ in Figure 14(b). The dashed vertical line marks ${R}_{2\mathrm{beam}}$, which is the smallest radius included when measuring ${\rm{\Delta }}{\text{}}V$. Bulge contamination is likely for NGC 2347, because ${R}_{\mathrm{bulge}}\approx {R}_{2\mathrm{beam}}$, but the other galaxies in the sample have much smaller bulges, so the likelihood of contamination from the bulge is lessened. The exact contributions of the bulge and eDIG to ${\rm{\Delta }}{\text{}}V$ are difficult to disentangle in detail, however, and it is likely that both contribute to the measured ${\rm{\Delta }}{\text{}}V$ at some level.

Zoom In Zoom Out Reset image size