Stellar Population Synthesis with Distinct Kinematics: Multi-Age Asymmetric Drift in SDSS-IV MaNGA Galaxies

We selected seven galaxies from the MaNGA survey to develop and test the algorithms in this study. The galaxies were chosen to have exceptionally regular gas and stellar kinematics (in terms of azimuthal symmetry, with little evidence for bar distortions), moderate inclinations ($20^{\circ}<i<50^{\circ}$), and clear evidence for asymmetric drift between their (single component) stellar velocities and ionized gas (as reckoned by the MaNGA data analysis pipeline data products; hereafter DAP). This list is by no means exhaustive of good candidates.

These galaxies were also chosen to have stellar masses and colors similar to the Milky Way (MW).¹¹1MW stellar mass is estimated to be $5.7\pm 1.3\times 10^{10}{\rm M_{\odot}}$ (Licquia et al., 2016), with an $i$-band absolute magnitude of ${\rm M}_{i}=-21.27\pm 0.37$ mag and $g-r$ color of $0.68\pm 0.06$ mag (Licquia et al., 2015). Photometric properties of these galaxies, sorted by stellar mass, are summarized in Table 1, where redshift, b/a (minor to major axis ratio), half-light radius, rest-frame $g-r$ and NUV-$i$ colors, SDSS $r$- and $i$-band absolute magnitude, stellar mass and Sersic index (n_S) are from the NASA-Sloan Atlas (NSA, Blanton et al., 2011).²²2http://www.nsatlas.org Except for n_S, all photometric quantities report their elliptical Petrosian aperture measurements. Magnitudes and colors have AB zero-points; absolute magnitudes and masses are scaled to H${}_{0}=100$ $\mathrm{km\ s^{-1}}$ Mpc^-1 throughout. A value of n${}_{S}=6$ is a hard upper limit in the Blanton et al. (2011) analysis. The NUV-$i$ color range is blueward of the red sequence (e.g., Figure 2 of Yan et al., 2016a) but in the redder half of the blue cloud. As we will see, the sample has peak rotation speeds between 200 to 300 $\mathrm{km\ s^{-1}}$, with both slow- and fast-rising rotation curves.

We refer to these galaxies by their plate-IFU (integral field unit) designation throughout the paper. Table 1 provides their corresponding, unique MaNGA ID (MID).

For the spectral analysis throughout this paper we use the “LOGCUBE” datacube from the Data Reduction Pipeline (DRP) of MaNGA (Law et al., 2016) from an internal release³³3MaNGA Product Launch (MPL) 5, v2_0_1. closely related to the versions in Data Release (DR)14 (Abolfathi et al., 2018). We spot-checked our analysis to verify that subsequent subtle changes in the data-processing through DR15 (Aguado et al., 2019) and the next internal data release⁴⁴4MPL 8, v2_5_3 do not alter the results in any significant way. Given the significant investment in MCMC analysis, we have retained the results from the older release.

MaNGA datacubes contain the reduced and combined spectra from the dithered IFU observations of a single galaxy. Spaxels have been resampled to a logarithmic bin of 70 $\mathrm{km\ s^{-1}}$ and to a common wavelength grid.⁵⁵5For reference, the median instrumental resolution (FWHM) varies between 250 $\mathrm{km\ s^{-1}}$ at the blue limit of 362 nm to 125 $\mathrm{km\ s^{-1}}$ at 950 nm, with a median FWHM value close to 150 $\mathrm{km\ s^{-1}}$ characteristic of the performance between 500 and 750 nm Law et al. (2016); Yan et al. (2016a). Data cubes also contain information on the mean and variance of the spectral line spread function (LSF). We use this LSF to match the resolution of our spectral templates to the data, spaxel-by-spaxel, before fitting. We measure the S/N of individual spaxels in MaNGA datacubes using the median of the spectra and inverse variance (IVAR in the datacube) within the wavelength range of 360-940 nm. In our analysis we ignore spaxels with S/N per Å less than 1.

Throughout this paper we conduct full-spectrum fitting of the observed spectra between 350 and 940 nm. The extensive wavelength coverage yields a wide range of stellar features including the full Balmer series dominated by hot young stars as well as strong stellar features from Ca H&K to the NIR triplet which are dominated by cool stars; together these features should contain the significant information on the the kinematics of the young and old populations.

Full-spectrum fitting uses pPXF. This code fits an observed spectrum with a set of spectral templates, convolved with a line-of-sight velocity distribution (LOSVD) in pixel-space. Throughout this study, we fit for the velocity and velocity dispersion moments of the LOSVD of all kinematic components. The code allows multiple kinematic components to be fit simultaneously; in our analysis we have leveraged this to fit the observed galaxy spectrum with a single component for the observed gas emission features and a single or multiple stellar kinematic components.

During our full-spectrum fitting of observed galaxy spectra, we fit for the following gas emission features; H$\alpha$, H$\beta$, H$\gamma$, H$\delta$, [OII]$\lambda 3727$, [OII]$\lambda 3729$, [OIII]$\lambda 4959$, [OIII]$\lambda 5007$, [OI]$\lambda 6300$, [OI]$\lambda 6364$, [NII]$\lambda 6548$, [NII]$\lambda 6584$, [SII]$\lambda 61716$ and [SII]$\lambda 6731$. The gas emission lines are modeled as Gaussian features and are treated as a single kinematic component independent of that measured for the stellar component(s). While we fix the relative flux ratio between the [OIII], [OI] and [NII] doublets, based on atomic physics, the fluxes of the gas emission lines are free parameters during the fitting process.

When fitting two stellar components we take advantage of the pPXF feature to constrain the relative weight of two kinematic components, referred to in the code as fraction but shortened to yfrac in this study. This parameter is a light-weighted quantity defined as the relative weight in the first kinematic component to that in all components:

where $w_{i}$ are the weights of the first component while $w_{j}$ are the weights of second component. In our analysis the first component will always correspond to the young stellar population and hence the use of the subscript ’y.’ We flux normalize the template spectra for our stellar components and hence, if the first (second) component represents the young (old) stellar population, yfrac is a relative indicator of star formation history. Though by default pPXF limits the user to constrain the relative weight of two kinematic components, in the course of our analysis we removed this limit to allow the development the algorithm in Section VI. This allowed us to constrain the relative weight of more than two kinematic components simultaneously during our full-spectrum fitting by providing a vector of relative fractional values (frac_n for the $n^{th}$ element), i.e.:

where $c$ is the number of kinematic components being fit, $\sum w_{n}$ is the sum of the weights of the $n^{th}$ component and $\sum_{k=1}^{c}\sum w_{k}$ is the sum of the weights of all components. Throughout our study, we allow the relative weight of the gas component to be free and unconstrained by yfrac.

The pPXF code can also minimize errors caused by imperfect spectral calibration, template mismatch, or scattered light by fitting multiplicative and/or additive Legendre polynomials or trigonometric series. In our analysis we simultaneously use additive and multiplicative polynomials of order 8 when fitting MaNGA galaxy spectra for their stellar kinematics with independent stellar population templates.⁶⁶6This order is based on experience gained from analysis of MaNGA data (Westfall et al. submitted) as well as analysis of SAURON data of galaxies in the Coma cluster (Shetty et al. submitted). When fitting mock spectra using the same stellar population templates we do not use any polynomials.

For our analysis, we adopt the SSPs of MIUSCAT (Vazdekis et al., 2012) as our template spectra for the stellar components. These models were developed using empirical stellar spectral libraries, particularly the MILES (Falcón-Barroso et al., 2011), Indo-US Library of Coudé Feed Stellar spectra (Valdes et al., 2004), and the CaT Stellar Library (Cenarro et al., 2001). These models were selected primarily due to their significant overlap in rest-frame wavelength coverage with MaNGA. The 346.5-946.9 nm span for MIUSCAT almost exactly matches the 360-1030 nm span for MaNGA at the typical redshift of the MaNGA sample. The spectral resolution (0.251 nm FWHM) is slightly higher at all wavelengths than the MaNGA data (cf.  Law et al., 2016; Yan et al., 2016a, respectively their Figures 18 and 20), which is important for properly convolving the data to the MaNGA instrumental resolution for kinematic analysis, as described in Westfall et al. (submitted). MIUSCAT SSPs used in our analysis cover an age range 0.03-14 Gyrs in 53 steps, metallicities ([M/H]) -0.66, -0.35, -0.25, 0.06, 0.15 and 0.26 (without $\alpha$-enhancement), and use the BaSTI isochrones (Pietrinferni et al., 2004) and a Kroupa initial mass function (IMF) (Kroupa, 2001). The metallicity range of these SSPs encompass the observed range of metallicities seen in large integral-field surveys sampling all galaxy types over a wide range of luminosity and radius (e.g., Sánchez-Blázquez et al., 2014; González Delgado et al., 2015; Zheng et al., 2017).

For the old stellar population we use SSPs older than 1.5 Gyrs as template spectra, and the remaining for the young stellar population modeling. This choice is motivated in Figure 1. Here we conduct a light-weighted full spectrum fit of each MIUSCAT SSP with a set of empirical stellar spectra representative of a broad range of stellar spectral types (see Appendix A) and sum the weights for the presented bins of stellar types. The plot demonstrates that except for the lowest metallicity models the relative weight in hot stars (O, B, and A type) falls below 10% above ages of 1.5 Gyrs, while the coolest stars begin to contribute substantially ($>$30%) after this time, namely the times-scale for the formation of the red-giant branch.

It is also the case that the intermediate stars (F,G) peak around 1 Gyr. While the hottest of these (F0 V) will still have significant Balmer absorption and Main Sequence (MS) anlifetimes of roughly 3 Gyr, for regions with on-going or recent star-foration the F-star contribution to the Balmer absorption equivalent width will be small compared contributions from hotter, more luminous MS stars. Nonetheless, based on this Figure, an argument can be made for an intermediate age range, say, between 0.5 and 3 Gyr, particularly for modeling galaxies or regions of galaxies that have not had recent star-formation. We take advantage of the multimodel distributions in Figure 1 to constrain the mix of stars used in the stellar library algorithms (Section VI.4). However, while an additional time-bin would greatly aid in determining the AVR, for simplicity of our development in this work we distinguish only between two age bins for measuring differential asymmetric drift signals.

In the MW solar cylinder from Aumer & Binney (2009) or from MW and M33 star clusters (Beasley et al., 2015) there is observed to be roughly a factor of 3 change in velocity dispersion between t=0 and 1.5 Gyr (assuming the birth population has the same dispersion as the molecular gas – something that is not clearly the case for the star-cluster measurements), and about a factor of 2 increase between 1.5 and 10 Gyr; in M31 (Dorman et al., 2015) the increase in velocity dispersion is roughly the same during these two time periods. Hence our choice of age bins appears a sensible starting point if the three large galaxies in the Local Group are representative of the larger population of intermediate-to-massive spiral disks.

We applied the MCMC algorithm to our sample of seven galaxies on a spaxel-by-spaxel basis. Figure 2 shows the measured, line-of-sight (projected) velocity and velocity dispersion maps for one galaxy for young and old stellar components in the second to fourth columns of the top row. (Maps for the remaining six galaxies in our sample are in Appendix D.) There is no smoothing applied to these maps; the color of each data point in the image corresponds to the kinematics of individual spaxels. Blank spaxels in these maps within the MaNGA IFU footprint are either regions masked in the DRP as foreground stars or, in the case of the MCMC only, represent spaxels where the MCMC did not converge within the maximum run-time permitted for jobs on the computer cluster to which we had access.⁸⁸8Job run-times were limited to 72 hours. Each job consisted of processing for several contiguous spaxels, leading sometimes to vertical sets of incompletion. For comparison we also show the derived velocity and velocity dispersion maps for gas and stars for a single stellar component. (In later sections these serve as the initial fitting step in multi-step alrgorithms.) These were measured directly from pPXF. Note the coherency of the velocity fields with the expected spider-diagram of isovels, as well as the smoothness of the kinematics on scales much larger than the beam size (2 arcsec fiber diameter). This smoothness and coherency persists for both single and two-component stellar kinematics. Keeping in mind that the kinematics are derived independently for each spaxel, this coherency indicates that our two-component stellar solutions from the MCMC algorithm are robust and measuring astrophysical kinematic signals. While there is much more scatter in the velocity dispersion maps, our approach is to avoid using these and higher kinematic moments.⁹⁹9NB: stellar velocity dispersion maps provided for the young and old components, denoted by a cross, are presented as qualitative metric to guide the development of the algorithm and are considered quantitatively unreliable much below the instrumental resolution until forthcoming improvements in the accuracy in the instrumental line-spread-function are concluded (D. Law, private communication). It is reassuring nonetheless to see that on average the younger component has lower dispersions than the older component as anticipated, and the old component has the qualitatively expected radial decline in value. These trends are seen for all galaxies in our sample.

While the asymmetric drift signal is readily apparent from Figure 2 to the practiced eye, this is more easily seen by extracting rotation curves using some fiducial geometry. Accordingly, Figure 3 shows the results of our MCMC algorithm in decoupling the deprojected tangential speeds of young and old stellar populations. In most cases and at most radii we see that both young and old stellar components lag in their tangential speed relative to the ionized gas (i.e., this is asymmetric drift), with the old component lag being larger, as one would expect from the presence of AVR. For comparison we show the deprojected tangential speeds for gas and stars for a single stellar component (again, derived directly from pPXF), from which the asymmetric drift signal also is readily apparent.

For Figure 3 we adopt geometric parameters derived from full, two-dimensional kinematic modeling of a monolithic inclined-disk. We use the method described in Westfall et al. (2011) and Andersen & Bershady (2013), applied simultaneously to the DAP gas and stellar kinematics. These parameters are summarized in Table 2. However, these geometric parameters are only used to define radial and azimuthal bins and deproject the kinematics; the individual fits to each spaxel are independent of these geometric parameters. Consequently the rotation curves in Figure 3 make no assumptions about kinematic or rotation-curve models except insofar as they share the geometries given in Table 2. Figure 2 and associated Figures in Appendix D are completely independent of any assumed geometric or rotation-curve model.

We focus on two fairly extreme cases: 8138-12704 and 8329-6103. Both galaxies have S/N above 10 Å^-1 even at large radii, and yfrac values between 0.1 to 0.2 (in the central regions) and 0.4 to 0.55 in the outer regions. The radial trend in yfrac is qualitatively what is to be expected for age gradients in spiral galaxies where the older population is centrally more concentrated.

The former galaxy 8138-12704 was chosen because of its exceptionally large asymmetric drift signal as measured for a single stellar kinematic component with respect to the ionized gas, yet clear presence of young and old stellar populations. At R = R_e, D_n4000 = 1.42 while $V_{\rm gas}-V_{\rm stars}=41.8$ $\mathrm{km\ s^{-1}}$. Of a sample of nearly 500 galaxies in MPL-5 selected to have regular kinematics and moderate inclinations (see Section  II.1), this galaxy is in the upper quartile in terms of deprojected gas rotation speed ($V_{\rm gas}=304.8$ $\mathrm{km\ s^{-1}}$). Since in general the asymmetric drift signal scales with velocity, it is useful to also look at the fractional difference between tangential speed of the stars and gas relative to the gas ($1-V_{\rm stars}/V_{\rm gas}$). For 8138-12704 this fraction is 13.7%, and it too is in the upper quartile.

In contrast 8329-6103 has an asymmetric drift signal that is in the lowest quartile for the sample of galaxies in both an absolute and relative sense: $V_{\rm gas}-V_{\rm stars}=11.8$ $\mathrm{km\ s^{-1}}$, and $1-V_{\rm stars}/V_{\rm gas}$ is only 6%. The galaxy has about 25% of the total stellar mass and 41% of the dynamical mass at R = R_e ($V_{\rm gas}=194.4$ $\mathrm{km\ s^{-1}}$) compared to 8138-12704, yet has a comparable D_n4000 of 1.36. Indeed, our analysis shows the two galaxies have yfrac values of 0.42 and 0.53 respectively. Compared to the MW (e.g., Licquia et al., 2016), 8329-6103 is almost a factor of two less luminous.

For the more massive galaxy, 8138-12704, the differential asymmetric drift signal is readily apparent at all radii sampled. To quantify this excellent kinematic separation of young and old stellar populations, Figure 4 shows the bivariate posterior probability distribution for velocity and velocity dispersion of two spaxels selected from near the major axis. One spaxel is on the approaching side at 0.3 R_e (3.5 arcsec radius) and the other is on the receding side at 1.1 R_e (12 arcsec radius). In both cases the peaks and 99.7% contours are well separated with the relative velocities and dispersions for the two populations. Moreover, they are differentiated as expected toward higher dispersions and lower speeds for the old population and vice versa for the younger population. While this is consistent with our priors, the presence of these distinct, age-dependent kinematics independently determined for all spaxels using an unconstrained model suggests there is ample information to extract a quantitative measure consistent with our qualitative astrophysical expectations.

For the less massive galaxy 8329-6103, Figure 3 shows the kinematics for the young and old populations are differentiated at smaller radii, as expected, but at larger radii the differentiation is less evident. This is reflected in Figure 4 as well. One spaxel is on the receding side at 0.6 R_e (4 arcsec radius) and the other is on the approaching side at 1.3 R_e (8.5 arcsec radius). For the inner spaxel, we see the anticipated velocity difference between young and old population, although the velocity dispersions are flipped in the opposite sense of what we might physically expect, i.e., the faster-rotating (young component) has a larger dispersion. However, the 67% probability contours are quite broad, and both components are at, or well below, the instrumental resolution. At the larger radius the young stellar component continues to appear to rotate more quickly, even though the probability contours substantially overlap. The velocity dispersions are well below the instrumental resolution in this limit.

These results are promising, indicating we should be able to derive separate kinematics for most of the MaNGA sample that have regular disk kinematics. Clearly Figure 3 shows some irregularities in the tangential speeds of young and old components at large radii, which correspond to low S/N. The behavior reflects an increasing frequency with decreasing signal-to-noise where the measured kinematics of the young stellar populations are dynamically hotter and rotate slower than the old stellar component. In Appendix B we demonstrate this behavior is due to degeneracies in disentangling the kinematics of the young and old stellar populations in this regime when using a global minimizer (i.e., MCMC) with a very broad range of template age and metallicity. Our subsequent algorithms that use a local minimizer with carefully constructed initial conditions largely eliminate this problem.

The simplest scheme for measuring the kinematics of two discrete stellar population components consists of giving pPXF the set of SSPs described in Section II. We again split the assignment of stellar population components at 1.5 Gyr regardless of metallicity. There are no constraints placed on yfrac; pPXF is free to assign weights to the templates as needed to minimize $\chi^{2}$. A single full-spectrum fit to the galaxy spectrum is made with these two (young and old) stellar kinematic components and, in the case of real galaxy spectra (as opposed to emission-line free mocks) a kinematically independent gas component. The initial velocities and velocity dispersions for the two stellar components are identical and are set at 0 and 210 $\mathrm{km\ s^{-1}}$respectively.

This simple algorithm performs well on our mock spectra (metric 1), yielding derived AD signals within 10% of the expected value 67% of the time even for AD signals as low as 5-10 $\mathrm{km\ s^{-1}}$(see Appendix C.1 and Table C1). While these results look promising, the test is highly idealized, absent noise and using the same templates in the mocks and the fitting templates. Consequently, we turn to apply this simple algorithm to real MaNGA spectra as a more definitive performance test. In applying the simple SSP algorithm to real galaxy data, we fit for the gas kinematic component simultaneously along with the stellar kinematic components and match the LSF of the SSP templates to that of the spaxels.¹⁰¹⁰10We use the PRESPECRES datacube extension providing each spaxels’ pre-pixelized (PRE) spectral resolution (SPECRES).

Application of the simple SSP algorithm to 8138-12704 (Section III) indicates the algorithm does not appear to reliably converge to the kinematics measured by the MCMC (metric 2). In Figure 2 we observe that for both the MCMC results and the simple SSP algorithm the younger component appears to have on average faster rotation than its older counterpart, as expected. However, the young velocities for the simple SSP algorithm are systematically lower (in amplitude) than derived from MCMC, while the older velocities are systematically higher; both young and old velocity fields are less smooth as well. Further, the $\sigma$ maps are noisy, particularly for the young component, and lack the expected structure of being azimuthally smooth with a radial decline. These qualitative conclusions for 8138-12704 are also found in the remaining six galaxies in our sample (Appendix D). This is quantified in the first three lines of Table 3.

Given that the more robust MCMC technique derived a solution for the kinematics different from our simple SSP algorithm, this suggests that the spectral fitting code (pPXF), which uses a Levenberg-Marquardt technique to find the local $\chi^{2}$ minima, is not able to converge to the global minima with real data.

To circumvent the pitfalls of a complex $\chi^{2}$ terrain we augment our simple SSP algorithm by adding an additional fitting step. This additional step allow us to improve initial conditions to position pPXF fairly close to where the global minimum is expected.

In a new first step, we conduct a full spectrum fit with only one stellar kinematic component and one gas emission-line component, using all SSPs (regardless of age) as templates for the single stellar component. This fit is essentially a stellar population synthesis of the spaxel and, accordingly, we use only multiplicative Legendre polynomials during this fit. Expectations based on the MW AVR are that the kinematics of the young stellar population should be close to that of the gas. Simultaneously the kinematics of the old stellar component should be similar to that from a single stellar component than the gas kinematics. Results from our MCMC analysis support these expectations. Hence, in the second fit we use the derived gas kinematics from our first fit as the initial conditions for the young component, and the derived stellar kinematics from our first fit as the initial conditions for the old component. During this second fit, we fix the kinematics of the gas to that derived from the first fit and use additive and multiplicative polynomials during the fit. We refer to this algorithm as ‘2-Step SSP’ in Table 3.

As seen in the Table this two-step SSP algorithm exhibits some modest quantitative improvement in the derived kinematics compared to the simple SSP algorithm. This demonstrates that the irregular $\chi^{2}$ space of the solution affects the reliability of the kinematic decomposition unless reasonable initial conditions can be placed on the kinematics. However, some observable differences persist in the velocities of the young component. The discrepancy suggests a lack of full convergence of the two-step SSP algorithms to the global minima of the $\chi^{2}$ space identified by the MCMC analysis. We find that although yfrac determined in the first step of the two-step SSP algorithm correlates well with that determined by the MCMC analysis, the relation has significant scatter; better initial conditions on yfrac could also be beneficial. Previously Katkov & Chilingarian (2012) used a Monte-Carlo technique to demonstrate that when conducting a spectral decomposition of co-spatial components an incorrect estimation of the relative light contribution of the components could bias the results.

To test if constraining the relative weights of the young and old components further improves results we provide a data-driven constraint on yfrac by using the template weights from the first full spectrum fit with a single stellar component. From this weight distribution yfrac is set in the second fitting step to the sum of the weights associated with the young SSP models ($\leq$ 1.5 Gyrs) over the sum of the weights of all SSPs. This constraint is imposed independently for each spaxel. We then add a third full spectrum fitting step to our SSP algorithm (i) where we update the initial conditions for the kinematics to that measured in the second step, but (ii) we remove the constraint on yfrac for the fit. The motivation for this step is to assume that the kinematic solution obtained in the second step is close enough to the global minima that the final step can free the fit to converge onto the global minima without constrains put in place on yfrac. We refer to this algorithm as ‘3-Step SSP’ in Table 3. It is immediately evident that across the board, all of the metrics improve significantly, and hence we adopt this as our final SSP algorithm. Figure 2 and Appendix D show the velocity and velocity dispersion fields for MCMC, simple and 3-step SSP algorithms.

The aim of using a multi-step approach to nudge pPXF toward a global minimum using iteratively improved initial conditions appears to be on target. In contrast, the initial two-step SSP algorithm had less refined initial conditions on the kinematics produced nosier kinematic maps for the two stellar components – compared both to this three-step algorithm and the MCMC results.

Our final SSP algorithm used to measure the two-component AD is a three-step process as follows:

1. 1.

   Use full spectral fitting with pPXF with one stellar and one gas (emission-line) kinematic component. The entire suite of SSP models are used as templates for the single stellar component. Since each SSP is normalized, the resulting spectral fitting weight distribution is light-weighted. The gas emission-line components assume Gaussian line profiles, and the kinematics for all lines are tied to a single value for velocity and dispersion per spaxel. During this fit, the initial conditions for the kinematic for all components are zero velocity with respect to the NSA redshift and velocity dispersion of 210 $\mathrm{km\ s^{-1}}$. All stellar kinematic fitting is done assuming a Gaussian LOSVD and multiplicative Legendre polynomials of order eight are used to match continuum shapes between observed and template spectra. The resulting kinematics and weight distribution of this fit provide informed constraints for the second fitting step.

2. 2.

   Use full spectral fitting with pPXF with two stellar kinematic components (young and old) and one gas (emission-line) kinematic component. The young kinematic component is restricted to use only SSP models with ages $\leq$ 1.5 Gyrs, while the remainder are used as templates for the old component; yfrac is constrained by the ratio of these weights derived in the first step. The initial kinematics of the young component is set to the derived gas kinematics from the first fit and the old component are set to those of the single stellar component. The kinematics of the gas component in this second fit is constrained to that derived by the previous fit, however the fitting code is free to optimize the relative weights of the different emission line features. Similarly, pPXF is allowed to optimize the relative weights of SSPs within young and old age bins, but the sum of these weights is constrained by yfrac from the first fit.

3. 3.

   Repeat the full spectral fitting in the same setup as the previous step with the initial kinematics of the stellar components updated to those derived in step 2, but with no constraint on yfrac. This full spectrum fitting provides the final measured kinematics for the two stellar components, and the weights in the two components provides the yfrac in the spaxel.

We use a representative subset of 51 stars from the empirical Indo-US Library of Coudé Feed Stellar Spectra (Valdes et al., 2004). As described in Appendix A these are well suited for our purposes of analyzing MaNGA spectra. Similar to what was done for the SSP algorithms, we begin by defining the templates for the two age components being fit.

In the context of SSPs (Section II.5) we define the young-component templates to be those containing significant Balmer features. For a stellar library, however this definition excludes stars cooler than late-F which we know are present in young stellar populations for any plausible present-day stellar initial mass function. This is quantified in Figure 7. Since the younger populations have both hot and cool stars, with metal lines contributed to the spectra by the latter, we select three ranges that always include the hottest stars but extended to progressively cooler limits in cases labeled A, B and C. Our aim is to include the smallest possible range of cool stars to minimize degeneracy with the older stellar populations, and hence simplify the fitting algorithm.

The choice of stellar templates for the old component is simpler since hot, Main Sequence stars are short lived and hence only present in young stellar populations, where they dominate the spectral continuum in the blue and visible portion of the spectrum. We therefore exclude hot stars (hotter than F) from the templates for the old component, effectively ignoring the possibility of extreme (very metal poor) blue Horizontal Branch stars from contributing significantly to the integrated light of a mixed stellar population. In one case (D) we restrict the older population to have only K and M type stars to understand if there is any sensitivity to the giant-branch effective temperature.

Table 4 summarizes the four distinct combinations of stellar templates we explore for the young and old stellar components. These template combinations place no restrictions on luminosity class.

As a starting point for the development of our algorithm, we test the following simple algorithms to identify pitfalls to be ameliorated by more complex techniques. These pitfalls highlight subtle but important aspects of fitting multiple kinematic components where velocities are offset by values comparable to the spectral resolution. There are no constraints on yfrac in any of these algorithms.

1. 1.

   Full Spectrum Fitting (Full): pPXF fits two kinematic components to the provided (mock or real) spectrum over the full wavelength range using the above mentioned template sets.

2. 2.

   Feature Fitting (Feature): Based on the knowledge that we expect to see hot stars only within young stellar populations, the kinematics of the young component are fit with pPXF using only wavelength regions of $\pm$1500 $\mathrm{km\ s^{-1}}$ about Balmer lines up to n=18 in the spectra while masking out the remainder, as shown in Figure 8. Because of the strong contamination from Ca H, we exclude H$\epsilon$; and because of the G band (CN), we only use $\pm$500 $\mathrm{km\ s^{-1}}$ for H$\gamma$. Given the $\pm$1500 $\mathrm{km\ s^{-1}}$ velocity width, the region contain H$\eta$ (n=9) and above form a continuous band from 385.36 to 367.21 nm. During this fit, two kinematic components are still applied in pPXF  despite the dominance of Balmer features in these regions for young populations, there remains an imprint of the kinematics of the older population in these regions, both in the Balmer lines as well as the metal lines within these windows. This fit is then followed by a fit using pPXF that excludes the Balmer regions. During this second fit, the velocity moments of the young component are kept fixed to the values derived from the first step.

Based on our tests using mock spectra (Appendix C.2) and Table C2), we find the performance of these two algorithms is comparable. This confirms expectations about where most of spectral information is stored for discriminating between young and old stellar populations. For simplicity and consistency with the SSP algorithms we adopt the Full algorithm.

Table C2 also reveals important performance differences between templates: Velocity systematics for the young component are minimized for template set C, while velocity systematics for the old component are minimized for template set B. Perhaps this is unsurprising given Figures 1 and 7. However it is concerning that neither template set optimizes both young and old component systematics. A further concern arises with template sets C and D. In these cases there are many occurrences where yfrac >0.85, i.e., the best-fitting solution yields too large a young-component contribution. This is due to the degeneracy in the contribution of cool stars in the young and old components; the pPXF likelihood minimization falls into a local minima where the stellar kinematics of the mock spectra is almost entirely modeled by the young component. An example of this is given in Appendix C.2, Figure C2. We to conclude that template set B is the best compromise. This leaves us with performance results, based on mocks, that are somewhat unsatisfactory.

Following the success of adding additional steps to improve the kinematic initial conditions with our SSP-based algorithms, we further developed the Full algorithm using the template set B by adding a step to provide initial conditions in an identical fashion as our two-step SSP algorithm. Foreshadowing further development, we refer to this algorithm as ‘2-Step 1-Bin SL’, where the ‘1-Bin’ designation refers to the grouping of the templates as a single unit for each of the stellar components.

An evaluation of the performance of this algorithm using mocks (Appendix C.2, Table C3) indicates no improvement, but we suspect this is due to systematics between the mocks generated by SSPs and the fits using the stellar library. We turn instead, and henceforth, to a comparison with the MCMC for real data. Velocity and velocity dispersion maps (Figure 2 and Appendix D) show marked qualitative improvement in the smoothness of the kinematic maps between our one-step and two-step stellar library algorithms. Table 5 quantifies this contrast, most pronounced for the young component velocities. There also are performance improvements also for the old-component velocities.

However, further inspection of the kinematic maps reveals velocities and velocity dispersions are quite noisy for the young component of galaxies with lower mass ($<2\times 10^{10}{\rm M}_{\odot}$) or bluer color ($({\rm NUV}-i)_{0}<3$). (These properties are correlated in general, and particularly so for our small sample). It is also the case that for $\sim$33% of the spaxels with $10<S/N<35$ both the simple and two-step 1-bin algorithms find very high values of yfrac compared to MCMC results. This is indicative of an incorrect allocation of cool stars to the young kinematic component in the stellar library algorithms, and suggests further algorithm improvement is desirable.

Since problems with our stellar library algorithms arise due to incorrect allocation of cool stars to the young kinematic component, we need an objective mechanism to appropriately assign weights for these stars in young and old components. The above stellar library algorithms attempt to steer the allocation of weights via a prescriptive template set (B). Additional constraint on the relative weights of the young and old component (yfrac) fails to improve algorithm performance with nearly all our quantitative metrics. Further improvement requires a better guided approach.

An astrophysically motivated mechanism to properly weight stellar contributions to young and old stellar populations is to assign priors based on what we would expect from stellar evolution. A limiting case is just to use SSPs but since the entire exercise here is to not restrict ourselves the specific weights produced by the mapping of theoretical isochrones to a template library, we aim to impose a less restrictive scheme that still is informed by the relative weights of different stellar types as we would expect from stellar evolution. As we show, it is possible to do this by coarsely categorizing stellar types in effective temperature while leaving pPXF the freedom to optimize template weights within these coarse (broad) categories.

Figure 7 shows the derived weights for our stellar library from fitting SSPs with pPXF; these results broadly match our expectations for stellar evolution and the impact of metallicity on, e.g., the temperature of the giant branch. The youngest populations, dominated by the hottest stars, indeed have most of their weight in B and A templates, but do contain significant weight in cooler stars. As the populations age (and the Main Sequence burns down while the asymptotic and then red giant branches are populated), the relative weight of the cooler stars gradually increases and the weight of the hotter stars decrease. The oldest SSPs are entirely dominated by the coldest stars, but the specific age where this occurs depends on metallicity because of the temperature of the stars on the horizontal and giant branches. Hence by constraining the relative weights of the templates used in the kinematic fitting such that they are at least consistent with that expected from stellar evolution at some granularity, we could disentangle the the relative weight of the cool stars in the two kinematic components in our fit. Just what granularity should be adopted is one question we answer.

To constrain the relative weights of the stellar templates in the two kinematic components we use pPXF to do a full spectrum fit of the target spectrum with SSPs for a single stellar kinematic component (plus a second component for gas, if the target is a real galaxy spectrum; see Section II). The resulting SSP weight distribution is translated into weights of our stellar templates for a young and old population via a reference look-up table containing the stellar weights derived from fitting our Indo-US stellar subset to each SSP (Figure 7). The weights are divided between young and old populations by summing over SSP weights younger or older than 1.5 Gyrs.

While this scheme provides reasonable estimates for individual stellar template weights, the detailed result can be strongly effected by noise in the galaxy spectra. Further, fixing the relative weights is almost equivalent to using the SSP spectra as templates themselves, which defeats the purpose of this algorithm. To relax these constraints we divided the stellar template into temperature bins. The purpose of this division is to constrain the relative weights only between but not within these bins for young and old populations during the full-spectrum fitting of the galaxy spectrum. The detailed distribution of weights for the individual stellar templates within each bin is a free parameter, tuned (by pPXF) to optimize the likelihood of the solution.

We divide templates into three temperature bins (each contains stars with a range of luminosity classes and metallicity) motivated by the fact that the young and old kinematic components both have cooler stellar types (F-M) as common templates, but the mix of intermediate-temperature stars (F-G) changes strongly with population age (see Figure 1), while only the young component has the hottest stellar types (O-A). The three bins are: Hot (which contains templates of stellar types O, B and A with $T_{\rm eff}$ roughly above 9000^∘ K), Intermediate (F, G with $T_{\rm eff}$ roughly between 5000^∘ and 9000^∘ K) and Cold (K and M roughly below 5000^∘ K). Since the unrealistic distribution of weights for the cool stellar types is causing the catastrophic failure in yfracwhen deriving the kinematics of the young and old kinematic components for our simple and 2-step 1-bin algorithms, this three-bin division fix the hot : intermediate : cool contributions in each kinematic component to reasonable values consistent with both the observed spectra and rough expectations from stellar evolution. Consequently we can dispense with the restrictions of the template sets A-D from Section VI.1. We refer to this as our 2-step 3-bin stellar library algorithm.

To implement constraints on the relative weights of the bins we utilize features in pPXF that (i) allow the relative weight to be fixed between multiple kinematic components, and (ii) grants the user the ability to tie together the kinematics of different components. Effectively we fit a set of template bins with tied kinematics for each stellar component. By fixing the relative weights of the bins of the young and old components we can avoid unrealistic template weights while giving freedom for template optimization. Nominally pPXF limits users to fixing the relative weight between only two kinematic components. For this study we modified the pPXF code to disable this limit, enabling us to split the templates for each component into as many bins as needed.

We find the addition of these multiple constraints on yfrac yield a significant qualitative improvements in the smoothness of the velocity and velocity dispersion maps for the lowest-mass and bluest galaxies in Appendix D. As seen in these Figures, this stellar library algorithm provides comparable qualitative performance to the SSP algorithm, and based on this we adopt this three-bin approach. Further comparison with the SSP algorithm is presented in Section VII.

The final algorithm (SL hereafter) can be summarized as a two-step process. The first step is identical to what is described in Section IV.3; it is a full-spectrum fit for a single stellar kinematic component using SSPs. This step is used to constrain the relative weights of empirical stellar spectra via a decomposition of the SSPs into relative weights for the empirical stellar spectra. Henceforth SSPs are not used. A second and final, two-stellar-component fit (also full-spectrum) uses the empirical stellar spectra with weights broadly constrained in three effective-temperature classes (Hot, Intermediate, and Cool). The second step uses a modified version of pPXF with fixed fractional contributions to the total fit for two stellar kinematic components, each with three different sub-populations. By definition, the three sub-populations for each kinematic component share the same kinematics. In detail:

1. 1.

   The relative weights for the Hot, Intermediate and Cool stellar bins for each (young and old) kinematic components are derived from the decomposition of the SSPs into individual stellar templates, as illustrated in Figure 1. The stellar template weights for the SSPs younger or older than 1.5 Gyr, respectively, are assigned to the young and old stellar components via the look-up table generated from the decomposition of each SSP (Figure 7).

2. 2.

   The initial kinematics of the young component are set to the derived gas kinematics from the initial fit, and the old component are set to those of the single stellar component.

3. 3.

   The kinematics of the gas component in this second fit is constrained to that derived by in the first fit. The fitting code is free to optimize the relative fluxes of the different emission lines.

4. 4.

   The fitting code is allowed to optimize the relative weights of the stellar templates within each temperature bin for each age component (young and old) but the relative sum of these weights is constrained by six yfrac values defined above. This flexibility minimizes the effect of template mismatch.

The velocity and velocity dispersion maps in Figure 2 and Appendix D show that our different algorithms display qualitatively similar kinematics. We focus further analysis on kinematics from spaxels within a $\pm 30^{\circ}$ wedge of each galaxies’ major axis, defined using the geometry in Table 2. Tangential velocities for each galaxy and tracer are corrected for their own systemic velocity determined by minimizing the difference in the approaching and receding components, excluding data at radii less than 2 arcsec. We find these systemic velocity values are consistent (at the $\sim$1 $\mathrm{km\ s^{-1}}$ level, on average) with values estimated from full two-dimensional kinematic modeling of a monolithic inclined-disk using the method described in (Westfall et al., 2011; Andersen & Bershady, 2013). There are variations between systemic velocities between tracers that are of order 3 to 5 $\mathrm{km\ s^{-1}}$. These differences are consistent with random errors, and they are negligible in our overall error budget.

Our two algorithms (SSP, SL) both use pPXF to minimize $\chi^{2}$. Because pPXF does not map the shape of the minima in $\chi^{2}$ space, these algorithms do not provide a direct estimate of the measurement uncertainty. In contrast the Markov-Chain Monte Carlo analysis we conducted for our sample galaxies do provide this information. From the MCMC chains we have marginalized the posterior probability distribution (PPD) for the kinematics (velocities, velocity dispersions) and relative light contribution of the young and old stellar populations of the galaxies. As mentioned in Appendix B, the shapes of the PPD for individual spaxels can sometimes be multi-modal, particularly at low S/N. We quantify the PPD width using a standard deviation (rather than, e.g., the median absolute deviation) specifically to include the impact of multi-modality.

We find that the uncertainty in the measured velocities of the two components correlates with the S/N of the spectra and yfrac. Figure 9 illustrates these trends by plotting the median width of the PPDs ($\tilde{w}_{\rm ppd}$) as a function of yfrac and S/N. The middle-panel shows that for the lowest values of yfrac the uncertainty in the measured velocity of young stellar component is the lowest for a given S/N. This is expected given that lower value of yfrac corresponds to lower contribution of the young stellar component to observed galaxy spectrum. Likewise, the right-panel illustrates that at a given S/N the reliability of the measured velocity of the old stellar component is higher when yfrac is lower, the relative contribution of the old stellar component is higher. Both panels illustrate the expected correlation between $\tilde{w}_{\rm ppd}$ and S/N. Since the systematic differences (seen in the next section) and the scatter between methods is smaller than the random errors exhibited in Figure 9, these values are well suited for characterizing the random errors in appplication of either the SSP or SL algorithms.

For reference, at a S/N (Å^-1) = 10 and yfrac = 0.4, the random errors in the young and old velocities are roughly 30 $\mathrm{km\ s^{-1}}$, while at S/N (Å^-1) = 30 the errors are between 5 and 10 $\mathrm{km\ s^{-1}}$ for young and old components respectively. On-line tables of quantities plotted in the middle and left panels of Figure 9 are available.

Figure 10 compares the median difference (spaxel by spaxel) between the derived velocities from our two algorithms binned in S/N and yfrac. This shows that the systematic differences are small in absolute value, but nonetheless non-zero. The SL algorithm’s (Section VI.5) measurement of the velocity of the young component is consistent with that of the SSP algorithm (Section IV.3) across S/N and yfrac to within $\pm$8$\mathrm{km\ s^{-1}}$, with an overall median difference of -5 $\mathrm{km\ s^{-1}}$ such that the SL algorithm tends to find the young component rotating faster. For the old component, the median difference is 7 $\mathrm{km\ s^{-1}}$ such that the SL algorithm tends to find the old component rotating slower.These systematic differences may reflect template mismatch or systematics in the final yfrac values of the two algorithms. We will further explore the cause of these systematics in future work, but it suffices to conclude here that the systematics are small and have little correlation with S/N or yfrac. As shown in the bottom panels of Figure 10, random errors almost always dominate over systematic differences between algorithms at the spaxel level.

As a summary measurement of asymmetric drift for two independent stellar components differentiated by age, we take the mean velocity from our two algorithms (SSP and SL), spaxel by spaxel, as our measure. For each spaxel we take $\tilde{w}_{\rm ppd}$ as the random error, and half the difference between the velocities from the two algorithms as the systematic error. When averaging the kinematics derived from, e.g., a radial bin of spaxels, we assume the random errors diminish with the usual quadrature averaging, while the systematic errors do not. Systematic errors in the binned data are taken as the median of the systematic error for individual spaxels in the bin.

Figure 11 shows the implementation of this summary scheme to a measure of asymmetric drift for young and old components relative to the ionized gas velocity for the seven galaxies in our sample. Line ratios (e.g., [OIII]$\lambda$5007/H$\beta$ versus [NII]/H$\alpha$) are consistent with HII-like regions for most galaxies at radii outside several arcsec; hence ionized gas velocities should be close to the circular speed of the potential. Results from Section VII.2 indicate that young-component velocities are unreliable for yfrac$<$0.15, and both young- and old-component velocities are unreliable when individual spaxels have S/N$<$10. These radial points have been indicated with open symbols. Referring to Figure 3, we see that most of the observed radial range for the galaxies in our sample are above these limits.

This final figure shows that we can indeed measure a distinct asymmetric drift for two stellar components for this sample of galaxies. Systematic errors dominate random errors in radial bins, but in most cases and at most radii the errors are substantially smaller than the difference between the asymmetric drift of young and old components. While the difference between ionized gas and young stellar component is small, it is non-zero ($6\pm 2$ $\mathrm{km\ s^{-1}}$ in the mean, with a median value of 8 $\mathrm{km\ s^{-1}}$). In all cases the asymmetric drift signals for young and old components straddle what is measured for a single stellar component, as one would expect. The asymmetric drift signal for the young component is generally small, with value of order 20 $\mathrm{km\ s^{-1}}$ or less, also as one might expect. The net effect of isolating the two population components by age is to elevate the asymmetric drift signal for the old component relative to the single-component asymmetric drift signal.

The results from Figure 11, combined with population age information will be exploited in future work to measure the AVR in MaNGA galaxies. However, it is already interesting to remark on the large range of AD signals observed in our sample of seven galaxies. There is also some indication of a trend with mass: The differential asymmetry drift signal for older populations is largest in the more massive galaxies, but this trend is far less regular (noting the small and large AD signals for 83320-9102 and 8482-3702, respectively) compared to the standard mass–rotation-speed scaling seen in Figure 3. While the addition of age estimates may transform the AD signal into a smoother trend in AVR with mass, the rather uniform $g-r$ colors of our sample suggests this expectation may be illusory. Given that there are at least three likely astrophysical paths to generating an AVR (see Section I), perhaps it is unsurprising to see a wide range in AD over such a modest range of stellar mass. Application of our methods to larger samples will provide statistics on the distribution of AD as a function of galaxy mass.

It is worth noting the limitations of the present methods and scope of analysis. We have selected galaxies with regular kinematics and with modest values of yfrac. One may well expect that interpretation of AD measurements will be problematic in cases where the kinematics are less regular or there are significant bar or oval distortions. It is imperative, therefore, to have an assessment of galaxy kinematic regularity, and to limit analysis to galaxies or radial zones within galaxies where the gas is on near-circular orbits and low-order distortions (e.g., m=1,2 modes) are weak or absent for both stars and gas.

We have also seen that our algorithms are limited in precision once the continuum S/N falls much below $\sim 20$ Å^-1. However, there is nothing preventing the methods described here to be applied in a more sophisticated analysis that fits multiple spaxels simultaneously in a parameterized model to increase the effective S/N; even a relatively simple tilted-ring approach would likely improve upon the S/N limitations of our current spaxel-by-spaxel scheme, e.g., at larger radii.

Nonetheless, the kinematics of the young component are not well measured for yfrac $<0.2$. This means that the AVR within the centers of galaxies with prominent, old bulge or pseudo-bulge populations, while typically observed at sufficient S/N, are not easily accessible given our current methods. While we have not probed yfrac $>0.8$ with our modest sample, we can expect that at least in the outskirts of later-type galaxies this condition will exist. We can also expect that measuring the kinematics for the old component will be difficult in this regime. Hence there is likely to be a sweet spot for our methods at intermediate radius, where both young and old disk components are comparably present. Such a sweet-spot is likely well-matched to a comparison with the MW solar neighborhood.

Finally, while we cannot offer a definitive recommendation between our SSP and SL algorithms, our preference is for the SL algorithm. Despite its greater complexity, the SL algorithm in principle should suffer less from template-mismatch systematics. Possibly supporting this statement is the observation that the SL algorithms’ velocity-dispersion maps for the young stellar component look smoother and more realistic (in amplitude) than for the SSP algorithm (Figure 2 and Appendix D).