The ALMaQUEST survey. VIII. What causes the velocity discrepancy between CO and H$\alpha$ rotation curves in galaxies?

We compare the CO(1-0) and H$\alpha$ kinematics in 34 nearby galaxies, selected from the ALMaQUEST and EDGE-CALIFA surveys. We use 3-D Barolo, a 3D tilted ring model, to derive the CO and H$\alpha$ rotation curves. Before comparing rotation curves in the 34 nearby galaxies, we found systematics between the MaNGA and the CALIFA data using eight MaNGA-CALIFA overlapping galaxies. We assume the rotation curves based on the MaNGA data are accurate and made the corresponding correction to the CALIFA data. Our result shows that $\sim$56% (19/34) of our galaxies present slower H$\alpha$ rotation curves compared to the CO rotation curves, with a median value of 6.5 km/s. The remaining galaxies (15/34) show consistent CO-H$\alpha$ rotation velocity within uncertainties. As a result, the Ha rotation may underestimate the total dynamical mass by ~6% for a circular velocity of 200 km/s (the median value in our sample). Furthermore, the velocity difference between the CO and H$\alpha$ rotation velocity is found to correlate with the velocity dispersion difference between CO and H$\alpha$, suggesting that the gas pressure plays a role in the velocity discrepancy. After incorporating the effect of pressure support due to the turbulent gas motion to our sample, the median value of the velocity differences decreases to 1.9 km/s, which in turn reduces the underestimation of dynamical mass to $\sim$2%. Finally, we also investigate the role that the extra-planar diffuse ionized gas (eDIG) plays in the CO-H$\alpha$ velocity discrepancy.


INTRODUCTION
The galaxy rotation curve is the relation between the mean orbital speed of stars or gas in a given galaxy and the radial distance from the galactic center. Practically, the rotation curves in galaxies can be measured by the redshifted emission or absorption lines such as HI 21 cm line (tracing the atomic gas), CO emission line (tracing the molecular gas), Hα line (tracing the ionized gas), and the stellar absorption lines (tracing the stars), as described in detail in Sofue & Rubin (2001). One of the important applications of rotation curves is the constraint on the dynamical mass distribution, which provides a strong evidence of the presence of dark matter (Sofue & Rubin 2001;de Blok et al. 2008). To derive the dynamical mass from rotation curves, one important assumption is that the rotating disk in the galaxy's radial equilibrium is balanced by the centrifugal force and the gravitational force. In other words, it assumes that the rotation velocity we measure is an optimal tracer of the circular velocity V c = R ∂Φ ∂R . This assumption requires the emitting material to be dynamically cold so that we can ignore the pressure term in the Virial theorem (Dalcanton & Stilp 2010;Iorio et al. 2017).
Several studies have shown that CO is a good tracer of the circular velocity (Davis et al. 2013;Leung et al. 2018;Lang et al. 2020). However, it seems not the case for the Hα emission line. Davis et al. (2013) compared the CO and ionized rotation curves in 24 early-type galaxies (ETG) from the AT LAS 3D survey. They find that 80% of their sample show slower ionized gas rotation velocities. The other 20% show consistent CO and ionized rotation within uncertainties. They attribute the systematic velocity discrepancy to the presence of arXiv:2206.05913v2 [astro-ph.GA] 14 Jun 2022 pressure support to the ionized gas. Similar studies have also been carried out for late-type galaxies. Levy et al. (2018) compare the CO and ionized rotation curves in 17 rotational-dominated late-type galaxies, selected from the EDGE-CALIFA survey (Sánchez et al. 2012;Bolatto et al. 2017). They also find that 75% of their sample show slower Hα rotation velocity than the molecular gas, while the other 25% show consistent CO and Hα rotations within errors. They attribute this to the presence of extra-planar diffuse ionized gas (eDIG). Since the rotation velocity of eDIG is slower than the ionized gas on the mid-plane (Levy et al. 2019;Bizyaev et al. 2017), the eDIG could reduce the observed mean rotation velocity, producing the systematic rotation velocity offset between CO and ionized gas. These studies suggest that the CO-Hα velocity discrepancy is common in nearby galaxies regardless of morphology types.
Understanding the reason behind such CO-Hα velocity discrepancy is important because the generation of mass models was conventionally done based on Hα rotation curves (Broeils & Courteau 1997;Östlin et al. 1999;de Blok et al. 2001;Kuzio de Naray et al. 2008;Epinat et al. 2009;Douglass et al. 2019) given that it is relatively easier to obtain a larger sample with Hα rotation measurement compared to other tracers. The lower velocity in Hα rotation implies that these mass models are potentially biased to a lower value. Hence, the comprehensive understanding of CO-Hα velocity discrepancy is important to quantify the potential mass bias based on Ha rotations.
The ALMA-MaNGA QUEnching and STar Formation (ALMaQUEST; Lin et al. 2020) and EDGE-CALIFA (Sánchez et al. 2012;Bolatto et al. 2017) surveys, which combined the CO(1-0) observations and integral field unit (IFU) spectroscopy data for 46 and 126 nearby galaxies respectively, provide us a good opportunity to compare CO and Hα rotation curves in a considerable amount of sample and also in a wide radial extension. Besides the high spatial and spectral resolution of CO(1-0) and Hα emission line provided by these datasets, the spatially matched resolution in CO(1-0) and Hα observations make it more ideal to compare the molecular and ionized gas rotation curves. In this study, with the selected 34 nearby galaxies from the combined ALMaQUEST and EDGE-CALIFA dataset, we aim to further investigate the CO-Hα velocity discrepancy with a larger sample size.
In Section 2, we describe the data used in this study. Section 3 and 4 present the method and tools we use to measure the CO and Hα rotation curves. The main results of this study and the interpretations of the CO-Hα velocity discrepancy are given in Section 5 and Section 6, respectively. Finally, we present the summary of this study in Section 7.

The MaNGA survey
Mapping Nearby Galaxies at APO (MaNGA; Bundy et al. 2015) is an integral field unit (IFU) survey as part of the SDSS-IV. With the technique of integral field spectroscopy, MaNGA provides the two-dimensional spatially resolved spectra of 10000 galaxies in the range 0.01<z<0.15. The spectral coverage of MaNGA is from 3600Å to 10300Å with spectral resolution R∼2000, which corresponds to σ inst ≈ 70km/s at the wavelength of Hα. The spatial resolution of MaNGA is 2 . 5, corresponding to the physical scale 1.5kpc at the average redshift (range from 1.1 to 5.9 kpc for the MaNGA sample included in this study). The spaxel size in the data products is 0.5 . Details about the MaNGA survey can be found in Bundy et al. (2015). The Hα data cubes for MaNGA galaxies are based on the SDSS Data Release 15 (DR15) analyzed through the MaNGA reduction pipeline (Law et al. 2016). The Hα datacubes used to measure the Hα rotation curves for MaNGA galaxies are taken from the MaNGA Data Analysis Pipeline (Westfall et al. 2019). The emission lines are extracted after subtracting the best-fit stellar continuum from the observed spectra. We only include the spaxels with S/N>3 for Hα flux in this study. The global SFR, stellar mass and the Σ * of the MaNGA galaxies used in this work are taken from the PIPE3D (Sánchez et al. 2016a,b) value-added catalog (Sánchez et al. 2018).

The CALIFA survey
The Calar Alto Legacy Integral Field Area (CALIFA) survey contains spatially resolved spectroscopy for 667 nearby galaxies (0.005<z<0.03) (Sánchez et al. 2012(Sánchez et al. , 2016c. CALIFA used PPAK IFU (Kelz et al. 2006) on the Calar Alto 3.5m telescope. The whole spectral coverage is from 3400Å to 7300Å, which is covered by two spectral gratings. The lower spectral resolution grating (V500) covers the range from 3745Å to 7300Å with spectral resolution R∼850, corresponding to σ inst ≈116km/s at Hα. The high-resolution grating (V1200) covers the range from 3400Å to 4750Å with spectral resolution R∼1650, i.e., σ inst ≈75km/s. The typical spatial resolution of CALIFA is 2 . 5, which corresponds to ∼0.8 kpc at the averaged redshift. More details about the selection procedure and statistical properties of galaxies in the CALIFA survey can be found in Walcher et al. (2014). The Hα data cubes for CAL-IFA galaxies are taken from the PIPE3D Dataproducts (Sánchez et al. 2016a,b). Similar to the processing of the MaNGA Hα data, The Hα emission lines are extracted after subtracting the best-fit stellar continuum from the observed spectra (Sánchez et al. 2016a,b), and only the spaxels with S/N>3 for Hα flux are included in this study.

The EDGE survey
Extragalactic Database for Galaxy Evolution (EDGE) survey (Bolatto et al. 2017) comprises 12 CO(J=1-0; 115.2712 GHz) measurement for 126 nearby galaxies with CARMA in the D and E configurations. The parent sample of EDGE-CALIFA galaxies is observed in E configurations, which consists of 177 galaxies selected based on their infrared brightness and biased to the higher star formation rate (SFR) in the CALIFA sample. The 126 galaxies are then selected and resampled in the D configuration and constitute the final EDGE-CALIFA data. More details about the EDGE-CALIFA survey can be found in Bolatto et al. (2017). The spatial resolution of EDGE data is ∼4 . 5, corresponding to the physical scale between 0.6 and 2.4 kpc for the EDGE sample included in this study. Since CALIFA and EDGE have different spatial resolutions, with EDGE having the worse resolution, we convolved the CALIFA data to match the EDGE resolution. We note that, in this case, the typical spatial resolution is similar to the spatial resolution of MaNGA data. The spectral resolution is σ inst ≈20 km/s in the data cube. We use the "dilated" mask to exclude noise (Section 3.1.1 in Bolatto et al. 2017). The mask is created by starting at 3.5σ (or greater) peaks in at least two consecutive channel maps in the data cube and then expands down to the surrounding 2σ contours. Finally, the additional 1 pixel is added to capture the low-level emission. To investigate the CO-Hα velocity discrepancy with a larger sample size, we include the sample used in Levy et al. (2018), which contains 17 star-forming, late-type galaxies. This set of galaxies are referred to as the EDGE-CALIFA Kinematics Sub-Sample (EDGE-CALIFA KSS hereafter).

The ALMaQUEST survey
The ALMaQUEST Survey contains 12 CO(J=1-0; 115.2712 GHz) MaNGA follow-up observation for 46 galaxies with ALMA in the C43-2 configuration. The field of view for each galaxy is ∼50 . with spatial resolution ∼2 . 5, which matches the MaNGA resolution.
The spectral resolution of 12 CO(1 − 0) for ALMaQUEST galaxies is σ inst ≈11 km/s. To investigate how gas drives the location and evolution of galaxies on the starforming main sequence, the ALMaQUEST survey covers a broad range of specific star formation rate, spanning from the green valley, main sequence, to the starburst regimes. More details in the ALMaQUEST survey can be found in Lin et al. (2020). In this study, we only include spaxels with S/N>3 for 12 CO(1-0).

METHOD
In this section, we introduce the general methodology for both the ALMaQUEST KSS and the EDGE-CALIFA KSS.

Tilted-ring model
Tilted ring model (Rogstad et al. 1974;van Albada et al. 1985;Begeman 1989) is a technique widely used to describe the kinematics of spiral galaxies. The main assumptions of this technique are the following: (i) The emitting material is settled in a thin disk, (ii) the rotation velocity of the material only depends on the distance from the center, and (iii) the motion of the material is dominated by circular motion. Under these assumptions, the tilted ring model breaks the rotating disk into a series of independent concentric rings with different radii, each of which has its geometric and kinematics parameters. Since the rotating disk is not necessarily face-on from the line of sight, the projection of each ring on the sky is usually parameterized as an elliptical ring. For a given projected ring, the projected velocity along the line of sight V los is given by: where V sys is the system velocity of the galaxies, V rot (R) is the deprojected rotation velocity of the ring with deprojected radius R, θ is the deprojected position angle (0 • is at the major axis of the projected elliptical rings), and i is the inclination between the line of sight and the rotating disk (0 • is face-on). The standard approach to apply the tilted ring model was to fit the 2D-velocity maps, but some drawbacks of such a 2D fitting algorithm exist. The most severe one is the beam smearing effect, which tends to flatten the gradient of the velocity field (van den Bosch et al. 2000;Swaters et al. 2009). In this work, we adopt an alternative approach by using a 3D-fitting software based on the concept of the tilted ring model called 3D BARALO (3DB hereafter; Di Teodoro & Fraternali 2015) to measure the kinematics and geometric parameters. d The radial averaged velocity dispersion derived from 3DB fit result. The region within one beam from the center is excluded from the calculation. For 8952-12701, which is lack of CO detection in the central region, We start the sampling from four beams away from the center. e The radial averaged velocity dispersion derived from 3DB fit result. The region within one beam from the center is excluded from the calculation. 3DB (Di Teodoro & Fraternali 2015) adopts a 3D tilted ring model that fits the 3-D datacube (two axes define sky coordinate and one axis defines spectrum direction) directly. In practice, 3DB fits the datacube ring by ring by building an artificial 3D model cube and finds the best fit through minimizing the residual calculated pixel by pixel between the datacube and the model cube. Before comparing the datacube and the model cube, the former need to be convolved to the same spatial and spectral resolution as the latter. This ensures the full control of instrumental effects, such as the beam smearing effect, which can significantly influence the fitting result at the region where the velocity gradient is large.

3D-Barolo
The residual function to be minimize in this study can be expressed as where n is the number of pixels, ∆r is the residual at each pixel, w(θ) is the weighting function of azimuthal angle. ∆r and the weighting function can be expressed as and where M and D are the flux values of the model and the data, respectively. θ is the position angle (θ = 0 • for the major axis). We note that there is no intensity weighting in the residual function. Hence, the faint pixels carry the same weights as bright ones.
In 3DB, each ring with deprojected radius R and width W is characterized by 8 parameters: Spatial coordinate of the centre (x 0 , y 0 ), V sys , V rot , velocity dispersion (σ v ), i, PA (the position angle between the North and the major axis, represented in degree from the North in  d The radial averaged velocity dispersion derived from 3DB fit result. The region within one beam from the center is excluded from the calculation. e The radial averaged velocity dispersion derived from 3DB fit result. The region within one beam from the center is excluded from the calculation. the counterclockwise direction), gas component surface density (Σ) and ring thickness along the axisymmetric axis (z d ).
In this work, we focus on the kinematics parameters V rot and σ v . To reduce the parameter space, some assumptions are made on other parameters before performing the fitting. First, we exclude the gas surface density from the fit by normalizing the flux in the model cube to that in the original data cube. Two different kinds of normalization are provided in 3DB: pixel-bypixel and azimuthally averaged. In the former case, the integrated flux of each pixel in the model cube is normalized to that in the data cube. In the second case, the model is normalized to the azimuthal-averaged flux in each ring. More details can be found in Di Teodoro & Fraternali (2015). In this study, we apply the pixelby-pixel approach since it allows us to accounts for the non-axisymmetric gas distribution. On the other hand, we fix the disk scale height to 100 pc. In Sec 6.2, we discuss the impact about this assumption. Finally, we assume a constant inclination, position angle, systematic velcity and kinematic center for each galaxy. 4. ANALYSIS

Fit rotation curves
To derive reliable rotation curves from the tilted ring model, it is crucial to set initial parameters properly. The initial dynamical center is the photometric center and the initial values for inclinations and PAs are derived from the photometric fit from NASA-Sloan Atlas (NSA) catalog (Blanton et al. 2011). The initial V sys is calculated as the midpoint between the velocities corresponding to the 20% of the peaks of the global line profile (Di Teodoro & Fraternali 2015). We assume that all rings in each galaxy share the same dynamical center, inclination, PA, and V sys . The initial values of rotation velocity and velocity dispersion are set to be 100 and 30 km/s, respectively.
We apply a two-step process to fit the rotation curves. In the first step, we use one unit of spatial resolution as the annulus width and fit the geometric parameters, rotation curves, and velocity dispersion at the same time. The purpose of the first step is to obtain reliable geometric parameters. In the second step, we fix the geometric parameters and fit the rotation velocity and velocity dispersion with half of the unit of the spatial resolution as the annulus width. These two steps ensure an accurate rotation curve measurement.
In the new version, we quantify the errors by 3DB's method for both the rotation velocity and the velocity dispersion. After finding the best model by minimizing the residual function, 3DB normalizes the value of the minimum of the residual function. Then, 3DB calculates a number of models by a Gaussian distribution oversampling. The center of the Gaussian distribution is the point that minimizes the residual function in the parameter space. Errors for each parameter are determined by the range where the residual increases 5 percent with respect to the minimum. i.e., the range where the residual function is smaller than 1.05. More details can be found in Di Teodoro & Fraternali (2015).
When measuring the rotation curves in galaxies, one of the important issues is the beam smearing effect, which causes the degeneracy between rotation velocity and velocity dispersion, and flattens the velocity gradient (e.g., van den Bosch et al. 2000;Swaters et al. 2009). The beam smearing effect becomes significant in the case where the velocity gradient is large, which is usually in the central region for a galaxy. Although 3DB allows better control of such instrumental effect compared to the 2D-fitting methods, some issues were found when measuring the central rotation curves, which will be discussed in more detail in Sec.5.1. The Hα rotation curves for our sample are fit using the same method under the assumption that CO and Hα are co-planar. In other words, CO and Hα share the same kinematics center, inclination, and PA.

The subsample selected for this study
To reliably measure rotation curves by the tilted ring model, there are practical upper and lower limits to the inclination of galaxies, as discussed in de . For high-inclination galaxies, it is difficult to derive rotation curves since the line-of-sight penetrates several layers with different distances from the center. For low-inclination galaxies, the decreased projected component of rotation velocity makes it difficult to derive reliable rotation curves. As adopted by de , the practical upper and lower lim-its of inclination are suggested to be 40 • and 75 • , respectively. We further exclude galaxies that show apparent bar structures, interacting galaxies, and galaxies with distorted morphology through visual inspection, as similarly done by Aquino-Ortíz et al. (2020). This selection makes sure that the gas motion is dominated by rotation and also validates the assumption of constant inclination and PA. After applying these selection criteria to the ALMaQUEST galaxies, 17 galaxies are selected and referred to as the ALMaQUEST Kinematic Sub-Sample (ALMaQUEST KSS hereafter). We note that the 17 galaxies from the EDGE-CALIFA KSS fit well to these selection criteria. The result of this research is based on the study on the combined 17 ALMaQUEST KSS and the 17 galaxies from the EDGE-CALIFA KSS.

The systematic discrepancy between CO and Hα rotation
Previous studies in the literature have found that Hα tends to rotate either consistently or slower than CO (Levy et al. 2018;Davis et al. 2013;Simon et al. 2005). In this study, we show that the CO-Hα velocity discrepancy is seen in both ALMaQUEST KSS and the EDGE-CALIFA KSS. The left panel in Figure 1 shows an example selected from the ALMaQUEST survey with consistent CO and Hα rotation. As a comparison, the right panel shows a galaxy from the EDGE-CALIFA survey which has an apparent velocity discrepancy between CO and Hα rotation. The same plots for all the galaxies included in this study can be found in Figure C.3 and Figure C.4. We note that, although the galaxies in the CALIFA survey are covered up to 2.5 effective radius (R e ; Walcher et al. 2014), which is larger than that of MaNGA galaxies (1.5R e ; Bundy et al. 2015), the radial range in which we compare the CO and Hα rotation curves are set to be between 0.3R e and 1.5R e for both the ALMaQUEST and the EDGE-CALIFA samples that are constrained by the coverage of CO. On the other hand, as we can see, there is a turnover feature shown in Hα rotation curves in the central region. Such a feature is commonly found in the Hα rotation curves among our samples. To test if this feature is artificial, we degrade the spectral resolution of our CO datacubes, which have original spectral resolution σ inst ≈11km/s, to the spectral resolution of MaNGA Hα data (σ inst ≈70km/s) by convolving them with a Gaussian kernel along the spectral axis. After fitting the convolved data by 3DB, we found such a feature also appears in the CO rotation curves. This phenomenon might be caused by the following two effects: (i) it is more difficult to derive reliable rotation velocity and velocity dispersion at the same time with a low spectral resolution; (ii) there are fewer data points in the most central region, which make it difficult to construct a reliable fit result. In this study, we exclude the rotation curves fit result within one resolution area from the center in the following analysis.
To quantify the velocity difference between CO and Hα rotation curves for a given galaxy, we first linearly interpolate the Hα rotation curve and resample the Hα rotation velocity at the same radii where we sample the CO rotation velocity. Then, for each galaxy, we calculate the velocity difference between CO and Hα rotation at each radius. We define ∆V for an individual galaxy as the variance weighted mean of the velocity differences V rot,CO − V rot,Hα . σ ∆V , the error of ∆V , is defined as the standard deviation of the velocity differences. In this study, we say that a galaxy has CO-Hα velocity discrepancy if ∆V >σ ∆V and has consistent CO-Hα rotation if ∆V ≤ σ ∆V .
As the spectral resolution differs between MaNGA and CALIFA datacubes, it is important to quantify whether there exists systematics or not in the measurement of the Hα rotation curves due to the difference in the instrumental characteristics. To conduct this test, we select in total 8 galaxies that are included in the main samples of both MaNGA and CALIFA surveys 1 . When comparing the rotation curves of the 8 MaNGA-CALIFA overlapping galaxies, we found systematics between the rotation curves fitted from the MaNGA Hα datacubes and that fitted from the CALIFA Hα datacubes. The latter tend to have slower rotation velocity. We adopt the definition of ∆V and use it to quantify the velocity differences between the MaNGA Hα and the CALIFA Hα rotation curves in these 8 overlapping galaxies (labeled as ∆V M C ). The median value of ∆V M C in these 8 galaxies is 9.8 km/s. Finally, we use bootstrap to test if this systematics is reliable, and a median value of ∆V = 8.5 with a standard deviation of 1.3 km/s is found. In the appendix B, we discuss the details of this test. This systematics causes ∆V in the EDGE-CALIFA KSS to become larger (or smaller in the ALMaQUEST KSS). To compensate for this systematics, we subtracted the ∆V in the CALIFA KSS by 8.5 km/s. All the discussions afterward have included this correction.
On the other hand, given the fact that the physical spatial resolutions of the ALMaQUEST KSS and the EDGE-CALIFA KSS are similar to each other, it is unlikely that the spatial resolution would raise an issue about the systematics between these two datasets. Figure 2 shows the relation between the ∆V and the physical spatial resolutions. As we can see, the physical spatial resolutions of the two datasets are comparable. Moreover, the lack of trend between the ∆V and the physical spatial resolutions are also found. The same results are found for the case which does not apply the MaNGA-CALIFA systematics correction to the EDGE-CALIFA KSS.
Based on the definition of ∆V , ∼56% (19/34) of the sample shows higher CO rotation velocity, while ∼44% (15/34) of the sample show consistent CO-Hα rotation. Some important parameters from 3DB fit result are listed in the table 1 and table 2. Figure 3 shows the normalized ∆V kernel density distribution (KDD) of our galaxies. In this figure, we use a normalized Gaussian kernel to represent a galaxy, which is centered at ∆V with the dispersion σ ∆V . After summing up all the Gaussians, the total distribution histogram is normalized to the unit area again. As we can see, the KDD shows that CO tends to have a higher rotation speed. The median ∆V of the ALMaQUEST+EDGE-CALIFA KSS galaxies is 6.5 km/s. In general, The distribution of ∆V in EDGE-CALIFA KSS (after correcting for the MaNGA-CALIFA systematics described in Appendix B) is similar in compared to the ALMaQUEST KSS. Figure 4 shows the relation between the ∆V and the radial mean of σ 2 Hα − σ 2 CO for each individual galaxy. The velocity dispersion differences are calculated within the region where ∆V is calculated. As we can see, there is a positive trend between these two parameters. We run the Spearman correlation analysis and find that the correlation coefficients give r s = 0.51 ± 0.19, 0.86 ± 0.22, and 0.71 ± 0.14 for the ALMaQUEST KSS, EDGE-CALIFA KSS, and the combined sample, respectively. Since the magnitude of the velocity dispersion is related to the magnitude of pressure gradient in the radial direction (Dalcanton & Stilp 2010;Levy et al. 2018;Davis et al. 2013), this gives us a hint that the radial pressure support may be responsible for the CO-Ha velocity discrepancy. We will discuss this scenario in the Sec.5.2.
We also study the dependence of the ∆V on other global parameters, such as star formation rate (SFR), stellar mass (M * ), and specific star formation rate (sSFR), morphology, inclination, CO V max , as well as on the Σ SF R and Σ M * within the area where the ∆V is calculated. These plots are shown in Figure C.5. Except for the Sersic index, we did not find any apparent trend between ∆V and other parameters. These plots suggest that (i) no correlation be found between the ∆V and the inclination; (ii) ∆V tends to be smaller in latetype galaxies (our galaxies are biased to the late-type, though).

The radial pressure support
Under the condition that gas motion is in equilibrium, the relation between the gravitational potential (Φ) and the V rot can be expressed by the following equation (Iorio et al. 2017): Where ρ is the volumetric density of the gas, and P is the gas pressure. Under the assumption that turbulence pressure dominates the pressure term, the pressure term can be expressed as ρσ 2 v (Dalcanton & Stilp 2010), where σ v is velocity dispersion. From eq. 5, we can find that the V rot is not a direct tracer of Φ if the pressure term is non-negligible. The asymmetric drift V A is defined as: , where V c = R ∂Φ ∂R is the circular speed. This could explain the CO-Hα velocity discrepancy if the V A term for Hα source is large. That is, there is larger radial pressure supporting the ionized gas that produce Hα emission. Under the assumption of a thin disk and isotropic velocity dispersion, V A can be expressed as Since 3DB can fit V rot , σ v and gas surface density simultaneously, we can calculate the V A for CO and Hα and determine the intrinsic circular velocity. Such a process is called asymmetric drift correction (ADC hereafter; Iorio et al. 2017). If this scenario can fully explain the "CO-Hα velocity discrepancy", we should expect

Notes about the asymmetric drift correction
From Eq.6 to Eq.7, one needs to assume a constant scale height and exponential distribution in the vertical direction, which is adopted in this study. Detail discussion about this assumption is in Sec 6.2. In addition, we also assume constant filling factors, path length, and temperature for the Hα and the X factor for the CO.
To obtain a smoother asymmetric drift correction, following Iorio et al. (2017), we use functional forms to describe the velocity dispersion radial profiles and the elements in the logarithm in Eq.7. The velocity dispersion is fitted by polynomials σ p (R, n p ) with the degree n p equal to or lower than 3. On the other hand, we fit σ 2 V Σ obs cosi by the the function (Bureau & Carignan 2002) Where f 0 is the normalization factor, and R c and R d are the characteristic radii. combining Eq.7 and Eq.9, we get No error of asymmetric drift correction is included in this study. The calculated V A for each galaxy is shown in Figure C.3 and Figure C.4. As expected, Hα tends to have higher V A due to its large velocity dispersion. In general, we did not find any radial trend of V A statistically. Figure 5 shows the KDD of ∆V from both the AL-MaQUEST KSS and the EDGE-CALIFA KSS after ADC, which is denoted as ∆V hereafter. Compared to Figure  3, we can see that ∆V is reduced in each of the samples after performing the ADC. The median value of ∆V becomes 1.9km/s, as opposed to the original value of ∆V ( = 6.5 km/s). This result suggests that the Hα emission source is sustained partially by radial pressure.

The presence of eDIG
Previous studies in the literature have shown that the rotation speed of extra-planar diffuse ionized gas (eDIG hereafter) decreases with increasing height above the mid-plane of the galaxy (Bizyaev et al. 2017;Levy et al. 2019). Such a vertical gradient in the rotation velocity could be responsible for the observed slower Hα rotation (Levy et al. 2018;Espinosa-Ponce et al. 2020;Sánchez 2020). Lacerda et al. (2018) suggests that EW(Hα) is a good proxy for diagnosing the ionization mechanism. Emissions with EW(Hα)>14Å trace the Figure 1. The left panel shows the rotation curves of one of our galaxies from the ALMaQUEST survey, 7815-12705. Blue one stands for Vrot(CO) and the red one is Vrot(Hα). ∆V , the variance weighted mean of the rotation difference between Vrot(CO) and Vrot(Hα), is calculated within the gray area and is presented as the green line. The purple line shows the velocity differences in each radius. The error of ∆V is the standard deviation of the velocity differences in sample points. The shadow areas present the errors of each quantity. The red bar at the lower-left corner is the beam size. The right panel shows the rotation curves in NGC5520, one of the galaxies in the EDGE-CALIFA survey which shows CO-Hα velocity discrepancy. We note that, in general, the radial range where we compare the CO and Hα rotation curves are similar between the ALMaQUEST samples and the EDGE-CALIFA samples, which covers from ∼ 0.3Re to ∼ 1.5Re. region where ionization is dominated by the HII region located in the mid-plane of galaxies, whereas emissions with EW(Hα)<3Å have recently been argued to primarily come from evolved, low mass stars (i.e., HOLMES) but the origin is still not well established. While the sources of the faint ionized gas have not been well established (e.g., Haffner et al. 2009), the EW(Hα) in between is likely the mixture of these two regions. Fol- lowing Haffner et al. (2009), we use the flux-weighted average of EW(Hα) in the region where we calculate the ∆V in each galaxy to divide our galaxies into two groups based on an EW(Hα)=14Å cut. The Hα intensity maps and the EW(α) maps are taken from the MaNGA DR15 PIPE3D value-added products (Sánchez Hα − σ 2 CO for the AL-MaQUEST KSS (red points) and the EDGE-CALIFA KSS (blue) samples. The error of the velocity dispersion difference is characterized by the standard deviation of the measurements in the annuli. The Spearman correlation coefficient (rs) shown in the upper-left corner is for the combined sample, which yields a strong correlation. The rs for the ALMaQUEST KSS and the EDGE-CALIFA KSS are 0.51 ± 0.19 and 0.86 ± 0.22, respectively.The positive correlation between ∆V and the velocity dispersion quadrature suggests that the gas pressure contributes to the CO-Hα velocity discrepancy. . We simply assume the presence of the eDIG is negligible in the high EW(Hα) group, whereas it might be non-negligible in the low EW(Hα) group.
We plot the ∆V density distribution for these two groups in Figure 6. Both groups show density distribution centered at the origin, but the high EW(Hα) group has a much smaller dispersion. We note that for the low EW(Hα) group, the ADC might not be accurate because of the following reasons: (i) the presence of eDIG "contaminates" the Hα velocity dispersion on the mid-plane; (ii) the assumption we made in Sec.5.2 (e.g. thin disk, constant scale height, etc) is no longer valid. Hence, the asymmetric correction alone for these galaxies may not be accurate, which leads to a larger dispersion. Therefore, this result suggests that the eDIG is likely to play a role in the CO-Hα velocity discrepancy, as suggested by Levy et al. (2018).

The MaNGA-CALIFA systematics
The results of this work assume that the rotation curves based on the MaNGA data are accurate and the MaNGA-CALIFA systematics is atrributed to the CAL-IFA data. Although the exact cause of the MaNGA-CALIFA systematics is not well known, this is a reasonable choice based on the fact that the MaNGA survey has higher spectral resolution, better control of instrumental dispersion and potentially lower systematics (Bundy et al. 2015;Law et al. 2021). For comparison, we discuss an opposite case, in which the offset is caused by the systematics in the MaNGA data only. We apply the correction of the MaNGA-CALIFA systematics to the AL-MaQUEST KSS. In this scenario, the median value of ∆V and the ∆V in the ALMaQUEST+EDGE-CALIFA KSS become greater, 15.0 km/s and 10.8 km/s, respectively as shown in the Figure 7. This result suggests that, while asymmetric drift correction again is able to further reduce the ∆ V, some other mechanism(s) is needed to fully account for the CO-Hα velocity discrepancy, which is beyond the scope of this work. 6. DISCUSSION

The robustness of the rotation curves and the ADC corrections
The spectral resolutions of both MaNGA IFS (σ inst ≈70 km/s) and EDGE-CALIFA survey (σ inst ≈116 km/s) of Hα are in general higher than the intrinsic astrophysical dispersion σ Hα . Hence, it is crucial to test if the relatively low spectral resolution of Hα influence the rotation curves and the ADC corrections. Hence, We convolve our ALMaQUEST KSS CO datacubes, which have original spectral resolution σ inst ≈ 11 km/s, to the  resolution of the CALIFA Hα observation (σ inst ≈ 116 km/s) with a Gaussian kernel and use 3DB to fit the CO datacube again. After this process, following the previous procedure, we compare the rotation curves derived from the original CO data to the rotation curves derived from the convolved CO data (using the same method described in Sec.5.1, denoted as ∆V CO ). The comparison of the rotation curves after ADC correction is denoted as ∆V CO . The result is shown in Figure 8. As we can see, both the ∆V CO and the ∆V CO density distribution show a Gaussian-like distribution centered at zero. By combining this result and the fact that the velocity dispersion of Hα is usually higher than that of CO, we conclude that: (i) The systematic velocity discrepancy between CO and Hα we found is not sensitive to their different spectral resolutions within the range considered in this work; (ii) Statistically, the ADC applied to the rotation curves presented in this work are reliable. Figure 8. We compare the Vrot(CO) fit from the original ALMaQUEST KSS datacubes, which have spectral resolution ∼11km/s, to the Vrot(CO) fit from convolved ALMaQUEST KSS datacube, which has spectral resolution ∼116km/s. The gray curve shows the kernel density distribution of ∆VCO before ADC, and the black one is for the one after ADC (∆V CO ).

Thin disk assumption
The discussion in this work so far relies on the thin disk assumption. With the presence of eDIG, the thin disk assumption is no longer valid. Hence, it is crucial to investigate whether thin disk assumption influences our result or not. The most straightforward test to do is use 3DB to fit our galaxies again with setting a larger scale height. We re-fit the ALMaQUEST KSS and the EDGE-CALIFA KSS again with a scale height of 1 kpc, and then compare the fit result to the ones with a scale height of 100 pc. The comparison is shown in Figure  9. As we can see, there is little impact on the rotation velocity. The same result is found on the velocity dispersion. Figure 10 show the scatter plots of spaxels which lie in the region where we calculate ∆V in the large EW (Hα) group. As we can see, most of the spaxels (90%) in these regions have EW (Hα) larger than 14Å, which are usually HII regions. Hence, the presence of DIG has little impact on the 3DB fit results for the high EW group galaxies. Therefore, we conclude that the thin disk model and a constant scale height are reasonable for these galaxies and the main conclusion of this work still holds.

Dynamical mass measurement
From the result presented above, one can find that the mass model based on ionized gas would potentially be biased low. For a galaxy with an intrinsic circular velocity of 200 km/s (the median value in our sample), the ∆V = 6.5 km/s (the median value of ∆V in our sample) would underestimate the total dynamical mass by ∼6%. In this study, we find that by applying ADC, the discrepancy between Ha and CO rotation curve can be further reduced to 1.9 km/s, corresponding to 2% underestimation in the dynamical mass measurement based on the Ha rotation curve. In the future, we will compare dynamical mass derived from the Hα rotation curves in our galaxies with other's studies, such as Zhu et al. (2018) and Aquino-Ortíz et al. (2020).

CONCLUSION
By combining the ALMaQUEST and the EDGE-CALIFA surveys, we analyze and compare the CO and Hα rotation curves in 34 rotational-dominated galaxies. Before combining the two datasets, we first compare the Hα rotation curves in 8 MaNGA-CALIFA overlapping galaxies. We find systematics between the Hα rotation curves measured from the MaNGA data and that from the CAL-IFA data, where the CALIFA rotation curves tend to have slower rotation velocities by 8.5 km/s. We compensate the systematics by adding the 8.5 km/s to the CALIFA dataset, and the conclusion is based on this premise. Our principal conclusions are: • 56 % of our galaxies show smaller Hα rotation velocity (10/17 in the ALMaQUEST KSS, and 9/17 in the EDGE-CALIFA KSS). The median value of the CO-Hα rotation velocity difference is 6.5 km/s (5.7 km/s for the ALMaQUEST samples, and 6.7 km/s for the EDGE-CALIFA samples). The remaining 44% of the sample shows consistency between CO and Hα rotations. For a galaxy with circular velocity 200km/s (the median value in our sample), ∆V =6.5 km/s would lead to ∼6% underestimation in the dynamical mass.
• The magnitude of velocity differences between CO and Hα rotation velocity, ∆V , correlates with the difference between CO and Hα velocity dispersion, which suggests the gas radial pressure gradient plays a key role to explain the CO-Hα velocity discrepancy. Under the assumption that the turbulence pressure dominates the radial pressure on the Hα emitting gas, we apply the asymmetric drift correction (ADC) to both CO and Hα rotation curves. After ADC, the median value of the CO-Hα rotation velocity difference reduces to , which assumes a Salpeter Initial Mass Function (IMF). When calculating the Σ SF R , only the star forming spaxels classified with the [SII] BPT diagnostic (Kewley et al. 2001(Kewley et al. , 2006 are included. Finally, we sum up the SFR spaxels within the region where we determine ∆V and divided it by the deprojected area the inclination correction applied.

B. THE TEST ON THE MANGA-CALIFA OVERLAPPING GALAXIES
This section aimed at investigating the systematic between the MaNGA and the CALIFA Hα data. We match the MaNGA and the CALIFA galaxies (Sánchez et al. 2016a,b), and 37 overlapping galaxies are found. To make the result applicable to the galaxies in our study, we apply the same selection criteria described in Sec 4.2, and finally, 8 galaxies are included in this test.
We use 3DB to fit the Hα rotation curves for these overlapping galaxies using both the MaNGA and the CALIFA datacubes. We follow the strategy described in Sec 4. For each galaxy, we use the same inclination and position for the MaNGA and the CALIFA fitting. The inclinations and the position angles are based on the MaNGA best fit.

B.1. The comparison of the rotation curves
After comparing the MaNGA and the CALIFA fitting result in these overlapping galaxies, we found a systematic in the rotation velocity. The best fit of the CALIDA data tends to have slower rotation velocity by 9.8 km/s, as shown in Figure B.1. We adopt the definition of ∆V and use it to quantify the difference between the MaNGA and the CALIFA rotation curves fit results in these overlapping galaxies, denoted as ∆V M C . To test the robustness of ∆V M C , we use bootstrap to re-sampling the velocity differences of the radial bins in all galaxies a number of times and calculate the median value of ∆V M C . Finally, the systematics of 8.5 km/s with a dispersion of 1.3 km/s is found in the median value of ∆V M C . This result suggests there is systematics between the MaNGA fit result and the CALIFA fit result based on our method. Given that MaNGA data has the better spectral resolution, we correct the ∆V s (i.e. the variance weighted mean of the velocity differences V rot,CO − V rot,Hα , as defined in Sec 5.1) in the EDGE-CALIFA KSS vaule by subtracting 8.5 km/s.

B.2. The comparison of the velocity dispersion radial profiles
The robustness of the velocity dispersion radial profiles are crucial for calculating the reliable asymmetric drift correction. Compared to H α , the measurement of the CO velocity dispersion should be more robust due to its higher spectral resolution. On the other hand, Law et al. (2021) found that the line spread function (LSF) of MaNGA Hα lines can be described well by a Gaussian, and the reliability of the astrophysical velocity dispersion of 20 km/s could be achieved when the S/N is sufficiently high enough. The low spectral resolution of the CALIFA data, however, raises the concern of the reliability of the Hα velocity dispersion measurement in the EDGE-CALIFA KSS. Under the condition that the spectral resolution is much lower than the astrophysical velocity dispersion, the LSF must be well known to get the reliable velocity dispersion. In this study, we adopt the instrumental dispersion of all CALIFA Hα data to be 116 km/s (Sánchez et al. 2016c). To validate this choice, we compare the Hα velocity dispersion radial profiles in these 8 overlapping galaxies measured between the MaNGA data and the CALIFA data.
In Fig B.2, the gray lines are the radial profile of the difference in the velocity dispersion between the MaNGA measurement and the CALIFA measurement for the 8 overlapping galaxies. The black vertical line is the typical uncertainty. As one can see, there is no apparent systematics in the velocity dispersion between these two measurements, suggesting that the typical instrumental dispersion we use for the CALIFA galaxies should be representative. On the other hand, the scatter of the velocity dispersion offset simply means the asymmetric drift correction measurement is subject to the uncertainty in the individual velocity dispersion measurement, which is not considered in this study. 3. The first column shows the CO velocity maps for the ALMaQUEST KSS. The color in these maps varies between ±300 km/s with negative velocity in blue. The interval of the contours is 50km/s. Both the axis values show the offset in arcsec from the center of the image. The black circles at the bottom left are the typical resolution of ALMaQUEST KSS data, which is 2 . 5. The second column shows the Hα velocity maps for the ALMaQUEST KSS with the features are identical to the first column. The third column shows the CO and Hα rotation curves. The x-axis is the galactocentric radius in the unit of effective radius (Re), and the y-axis is the velocity in km/s. Blue ones are CO rotation curves and red ones are Hα rotation curves. The purple lines connect the velocity differences between Vrot(CO) and Vrot(Hα) at each radial bin. ∆V, which is the variance weighted mean of the velocity differences, is calculated within the gray area. The error of ∆V is the standard deviation of the velocity differences in sample points. The green line shows the value of ∆V. The shadow areas present the errors of each quantity. The blue and red dash lines present the VA radial profiles of CO and Hα, respectively. The fourth column shows the CO and Hα rotation curves after asymmetric drift corrections (ADC). The features in the fourth column are identical to the third column. The fifth column shows the CO and Hα velocity dispersion radial profile. The x-axis is the galactocentric radius in the unit of Re and the y-axis is the velocity dispersion in km/s. The CO rotation curves are shown in blue and the Hα rotation curves are shown in red. The red points are ALMaQUEST KSS, and the black points are the EDGE-CALIFA KSS. The Sersic index comes from SDSS nsa catalog. For the ALMaQUEST KSS galaxies, the stellar mass and SFR are listed in Lin et al. (2020). For the EDGE-CALIFA KSS galaxies, the stellar mass and SFR comes from Bolatto et al. (2017). The Σ * and ΣSF R in this figure are calculated within the region where ∆V be calculated. Check Appendix A for more details.