THE DiskMass SURVEY. III. STELLAR KINEMATICS VIA CROSS-CORRELATION

A cornerstone of the DiskMass Survey (DMS; Bershady et al. 2010a, hereafter Paper I) analysis consists of stellar kinematic measurements performed in the low-surface-brightness, low-velocity-dispersion regime of the disks of late-type galaxies via integral-field spectroscopy (IFS). As discussed in Bershady et al. (2010b, hereafter Paper II), the robust measurement of the line-of-sight (LOS) velocity dispersion (σ_LOS) is a critical first step when measuring the quantities of fundamental interest to our survey, which include: the shape of the stellar velocity ellipsoid (SVE); disk-mass surface densities, Σ_dyn; stellar mass-to-light ratios, ϒ^disk_*,λ; the dark-matter-halo density profile, ρ_DM(R); and the maximum disk mass fraction, $\mathcal {F}_{\ast,{\rm max}}^{\rm disk}$. This third installment of our series presents a detailed description of our general approach to extracting stellar kinematics from DMS galaxy-continuum IFS.

Although measurements of stellar kinematics via stellar-atmospheric absorption lines are prevalent in the literature, most studies typically target hot, spheroidal systems, such as elliptical/early-type galaxies or spiral-galaxy bulges, e.g., the SAURON Project (Emsellem et al. 2004; Falcón-Barroso et al. 2006). With the caveat that studies of elliptical galaxies are starting to push toward very low surface-brightness limits (e.g., Weijmans et al. 2009), observational constraints based on the desired spectral signal-to-noise (S/N) and spectral resolution in elliptical-galaxy studies can be less taxing than for similar observations of a typical spiral-galaxy disk: early-type galaxies have higher surface brightness—factors of order five in flux (∼2 mag arcsec⁻²)—and more substantial LOS stellar motions—factors of order 10 in velocity dispersion—such that sensitivity limits can be less stringent and spectral resolutions may be lower for similar quality kinematic measurements. Moreover, one expects spectra of early-type systems to be relatively simple compared to late-type-galaxy disks; the latter are typically more strongly effected by nebular emission lines and varied stellar populations (cf. Trager et al. 2000; Yi et al. 2005; Jeong et al. 2009).

Nonetheless, over the last ≳30 years (see van der Kruit & Freeman 1984), stellar kinematics have been measured in a number of late-type-galaxy disks. In particular, a series of work starting with Bottema et al. (1987) and summarized in Bottema (1993) remains one of the more exhaustive observational studies of this type. Some more recent examples include: MacArthur et al. (2009)—long-slit spectroscopy for eight Sbc–Scd galaxies out to ∼1.0 disk scale lengths; Noordermeer et al. (2008)—PPak IFS for four Sa–Sab galaxies out to ∼2.0 disk scale lengths; Ganda et al. (2006)—SAURON IFS for 18 Sb–Sd galaxies out to ∼1.5 disk scale lengths; and Kregel (2003)—long-slit spectroscopy for 17 edge-on Sb–Scd galaxies out to ∼2.0 disk scale lengths. All of these studies measure up to, and occasionally beyond, second moments in the stellar line-of-sight velocity distribution (LOSVD). The analysis techniques used by these studies vary: Bottema (1993) and Kregel (2003) focus on the cross-correlation (CC) method (Simkin 1974; Tonry & Davis 1979; Statler 1995) using single stellar templates; MacArthur et al. (2009) fit their galaxy spectra in wavelength space using a variant of the Movel algorithm with a library of stellar-population-synthesis templates; and the other studies mentioned have used the penalized-pixel-fitting (pPXF) method of Cappellari & Emsellem (2004), again combining stellar templates to provide an optimal match to the galaxy spectra in wavelength space. Many more methods exist; see reviews by Kuijken & Merrifield (1993) and de Bruyne et al. (2003). However, of all the implemented approaches, the CC method and methods that fit directly to the pixelated observations in wavelength space (hereafter referred to has "direct pixel fitting") are currently the most prevalent in the literature.

The CC method has recently been abandoned in many analyses due to complications in the error analysis (Rix & White 1992) and spectral masking (Cappellari & Emsellem 2004). This is despite prevalent lines of reasoning and evidence available in the literature demonstrating that a CC approach minimizes template mismatch and S/N limitations (Bottema 1988; Bender 1990; Statler 1995; Kregel et al. 2004), both of which are expected to be important for our survey. Thus, from its inception, we have preferred the CC method in our survey analysis over direct pixel fitting, and we pursued the CC method for the analysis of our pilot observations (Westfall 2009). In doing so, we were prompted to perform an in-depth study of the CC method from first principles due to lines of inquiry that we were unable to answer with readily available CC software. This study has become a critical aspect in our understanding of the robustness of our stellar kinematics. Moreover, our implementation of the CC method provides for proper accounting of spectral errors and enables (sometimes pervasive) use of masks to eliminate contaminant spectral features, and, thus, it represents a step forward in the applicability of the CC approach to stellar kinematics.

The primary concerns in any method used to measure stellar kinematics from galaxy-continuum spectra are (1) S/N limitations, (2) spectral masking, (3) template mismatch (including continuum differences), (4) appropriateness of the prescribed LOSVD, and (5) instrumental-broadening corrections. Some of these topics are specific to the observational setup, but most are ubiquitous for all stellar-kinematic analysis techniques. In this paper, our discussion focuses on items 1 and 2; items 3–5 are omitted for the following reasons.

Template mismatch is a rich subject, deserving of detailed discussions. The current modus operandi of most direct-pixel-fitting methods is to address template mismatch by creating an optimized template composed of multiple stellar spectra of different spectral type and luminosity class that is optimized by chi-square (χ²) minimization (as in, e.g., pPXF). Indeed, such approaches are physically motivated due to the composite stellar populations observed in any galaxy spectrum. We have begun to discuss template-mismatch issues in the context of our CC method in Paper II, taking advantage of the isolation of kinematic versus template-mismatch information in the CC function: We demonstrated that, for single giant-star templates ranging in spectral type from early G to early M, the measured velocity dispersions for an example elliptical galaxy (UGC 11356) vary by ∼15%, whereas they vary by ∼25% for an example spiral galaxy (UGC 6918), with only modest changes to χ² as assessed in wavelength space (see Figure 11 of Paper II). Moreover, we showed that the template with the minimum χ² does not necessarily represent the template with the least template mismatch as determined via the CC function. Thus, χ² optimization of templates may complicate the assessment of template mismatch; degeneracies in template matching in direct pixel fitting do not necessarily correspond to identically derived LOSVDs, and vice versa. That is not to say that single stellar templates are preferred in general or for our own CC method, but it does illustrate the complexity of investigating template mismatch both statistically and astrophysically. Our work in Paper II demonstrates that the CC function may provide a very effective means of diagnosing template mismatch, beyond what we have gleaned from the literature. Thus, further development regarding template mismatch using CC analysis will be central to forthcoming papers in our series. For now, Paper II provides our preliminary assessment of the power of the CC method in assessing template mismatch of single stellar templates. Herein, we demonstrate our handling of the central hurdle to this template-mismatch analysis, which concerns properly applied masking techniques, by laying out the formalism of our CC method and the basic performance in the absence of template mismatch.

Effectively any parameterized LOSVD can be accommodated within the formalism of the CC method as long as the spectroscopic data are able to provide robust measurements of all its defining parameters. For this paper, all our discussion focuses on results based on a Gaussian LOSVD. An extension of our CC method to include robust measurements of non-Gaussian parameters (such as the skew, h₃, and kurtosis, h₄, parameters in the Gauss–Hermite series of van der Marel & Franx 1993) will likely require a detailed and nuanced discussion, such that it is beyond the scope of this paper. Indeed, such higher-order moments of the LOSVD may be interesting for investigating the presence of thick disks or triaxial halos as evidenced by the stellar line profiles; however, in general, we expect the LOSVD for our nearly face-on disk galaxies to be relatively simple and nearly Gaussian. In his implementation of the Gauss–Hermite series in a CC method, Statler (1995) found that robust measurements of non-Gaussian moments were limited to spectra with S/N ≳ 30, which is near the upper limit for individual spectra in our survey. Thus, discussion of higher-order moments will likely be most (or only) relevant to the analysis of our stacked spectra (see Section 3.1 of Paper II); therefore, this discussion is deferred to later papers, if relevant to the analysis therein.

For the DMS, instrumental-broadening corrections are largely marginalized because our template stars are observed using the same spectrograph setup as used for the galaxy observations. Second-order effects may still exist given the wavelength dependence of the resolution (effecting resolution differences due to the relative Doppler shift of galaxy-spectrum and template-spectrum rest wavelengths), variation in the resolution due to different entrance angles of each fiber in the pseudo-slit, and errors in the reproduction of the spectrograph setup between observing allocations. Our derived kinematics include such corrections as formally described in Appendix A—this approach is applied to SparsePak (Bershady et al. 2004, 2005) and PPak (Verheijen et al. 2004; Kelz et al. 2006) IFS in K. Westfall et al. (2011, in preparation), hereafter Paper IV; however, instrumental-broadening corrections are more appropriately discussed on a case-by-case basis when measuring stellar kinematics for each galaxy data set, as opposed to the more general discussion appropriate to the present paper.

Thus, our discussion herein proceeds as follows: We discuss the CC method from first principles and address some primary concerns for its use in measuring stellar kinematics in Section 2. We also provide greater detail into the specific algorithm used by our CC software in Section 3. In Section 4, we present a large suite of Monte Carlo simulations used to test the accuracy and precision of our CC method. Given its prevalent use in the literature, Section 4 also presents a comparable set of simulations using the direct-pixel-fitting, pPXF code (Cappellari & Emsellem 2004) as a fiducial method by which to compare our CC approach. Finally, we provide a brief demonstration of our CC software as applied to SparsePak IFS of UGC 6918 in Section 5, which also includes a comparison with pPXF results. The majority of the discussion focuses on applications performed with the observed spectral resolution and range measured for SparsePak IFS near the Mg ib triplet; however, these results will be generally applicable. We summarize the primary findings of this study in Section 6.

Stellar kinematics from DMS spectra are extracted via a CC analysis (Simkin 1974; Tonry & Davis 1979) that largely follows and is mathematically identical to the prescription of Statler (1995). However, after an in-depth analysis of the systematics involved, we have developed a different approach that addresses the practicalities of spectroscopic observations. That is, while we have not altered the fundamental concepts behind the CC method, we approach CC analysis with a novel algorithm directed at working within the constraints of real data. Our algorithm involves enforcing a symmetry between the treatment of the galaxy spectrum and the convolved template used to fit it, just as in direct pixel fitting, which is not true of previous implementations. We emphasize that this algorithm is more accurate than previous usage of CC analysis to obtain stellar kinematics and is, therefore, the preferred method. The degree of inaccuracy of previous implementations depends on the spectral range and the equivalent width of spectral features near the edges of the observed spectra.

As the basis for our alterations to the nominal CC algorithm, we re-present the fundamentals of CC kinematic analysis and establish the basis of our CC approach in Section 2.1. We then explain the rational for our new CC approach as resulting from issues considered during the development of our CC software⁴ in Sections 2.2–2.4. Finally, in Section 2.5, we provide a concise statement of the fundamental tenants behind our CC approach.

2.1. Fundamentals

The optical continuum of a galaxy is largely constituted of emission from the stellar photospheres, plus associated stellar-atmospheric absorption lines, of a spatial ensemble of galactic stars. The ordered (bulk) motion of the stars in the galactic potential causes a systemic shift of the stellar-atmospheric absorption lines, while the random motions—due to the distribution of the stellar orbital eccentricities about the tangential speed—serve to broaden the absorption lines. Let I be the ideal template spectrum that identically represents the sum of all rest-frame emission from sources within the solid angle of an observed aperture, and let F be the true velocity distribution function (VDF) along the LOS of all emission sources within the same aperture. The true F could be a wavelength-dependent function allowing for dynamical sub-populations of varying spectral type as is seen in the Milky Way (e.g., Dehnen & Binney 1998); practically, however, such complex functions are not generally applied. Provided these definitions of I and F, the galaxy spectrum, G, can be written as

$$\begin{equation} G = I \otimes F, \end{equation} \tag{ 1 }$$

where ⊗ denotes convolution. The only element in Equation (1) that is truly available to the observer is G, and measurement of G is limited by observational technology. Therefore, in all practical applications of Equation (1) to observed data, I is replaced by T—a single or optimized composite stellar template—and F is replaced by a function B that is parameterized, at least to some degree. Assessing the ideality of T is difficult, but not intractable, as we have demonstrated in Paper II.

In direct pixel fitting, one solves Equation (1) by minimizing the difference, assessed by χ², between the two sides of the equation for an observed G and I = T by iteratively adjusting F = B. Contrarily, the CC approach to inverting Equation (1) is to perform a relative frequency component analysis between T and G using the CC function, X = G○T where ○ denotes correlation. The solution to the broadening function is then determined by finding B such that

$$\begin{eqnarray} G \circ T & \approx & (T \otimes B) \circ T \\ X & \approx & X_T. \nonumber \end{eqnarray} \tag{ 2 }$$

The accuracy of Equation (2) is dependent on the suitability of the parameterization of B and the degree to which T represents the sampled stellar population. The advantage of CC analysis is that the equivalent-width weighted average shape of all lines commonly found in G and T are represented via a single line profile in X, which is compared to the intrinsic line shape of the average of all lines in T. Therefore, by focusing on a correct representation of the peak of the CC function, the determination of B should be dominated by the line shape instead of the line depth, effectively separating spectral properties primarily dominated by kinematics (line shape) from those primarily dominated by the chosen template (relative line depth). In direct pixel fitting, the fitting algorithm is more equally subject to matching both the line shape and relative line depth.

As written, a fitting algorithm that minimizes the difference between the two sides of Equation (2) requires both a correlation and a convolution for every fit iteration. Therefore, use of Equation (2) with spectra containing few to many thousands of discrete samples can be computationally intensive. As such, literature definitions of the CC analysis method have employed the commutative properties of the convolution operator to write:

$$\begin{eqnarray} G \circ T & \approx & (T \circ T) \otimes B \\ X & \approx & A_B, \nonumber \end{eqnarray} \tag{ 3 }$$

where we use A to denote an auto-correlation (AC) function. Therefore, a fitting algorithm that minimizes the difference between the two sides of Equation (3) only requires one convolution for every fit iteration (just as in direct pixel fitting) because it is only necessary to calculate the AC of T once. Calculations are further expedited by performing convolutions and correlations in Fourier space.

The right-hand sides of Equations (2) and (3) are mathematically identical (X_T = A_B); however, as we demonstrate in the next section, this is not true for the practical application of these equations to real data. Observed spectra are both discretely sampled and truncated by our current detector technology, which has different effects on the practical calculation of X_T and A_B. These effects dictate our departure from the literature formulation of CC analysis, due to its inappropriate assumption that X_T = A_B, by following the procedural implications of Equation (2) as opposed to Equation (3). The motivation and relevance of this choice are developed in the next few sections.

2.2. Cross-correlation with "Detector-censored" Data

A fundamental consideration when using algorithms dependent on the Fourier analysis of spectra is the finite sampling of data, both in terms of the individual velocity channel width and the velocity window (the observed spectral range, OSR) available to the detector—see Section 2 of Brault & White (1971) for a detailed discussion of these concerns. These considerations are fundamental to our assertion in the previous section that the mathematically equivalent quantities X_T and A_B, as defined above, lose their equivalency when applied to real data. This section is primarily concerned with demonstrating this effect via an example illustration.

Let us identically (without noise) know I at all velocities and define F such that the velocity range over which F ≠ 0 is approximately ΔV_F centered at V = 0. We can then create a synthetic galaxy spectrum, I_F = I ⊗ F. Also, given our explicit knowledge of I and F, we can define X_I and A_F as the ideal counterparts of X_T and A_B, respectively, where the equivalency X_I = A_F must be true. Now, let us sample I_F and I over a limited velocity range, ΔV, centered at V = 0 to obtain (I_F)' and I', where ΔV ≫ ΔV_F. Using these spectra, we can create the detector-censored correlation functions $X^{\prime }_{I^{\prime }} = (I_F)^{\prime } \circ I^{\prime }$ and A'_F = (I'○I') ⊗ F. If the detector censoring is immaterial to the correlation functions, $X^{\prime }_{I^{\prime }}$ should be equivalent to A'_F; however, we show in Figure 1 that they are, in fact, not equivalent.

Zoom In Zoom Out Reset image size

Figure 1 follows the experiment outlined above using a K1 III template (HD 167042) as I and a Gaussian broadening function with a dispersion of 70 km s⁻¹ for F; we adopt ΔV_F ≈ 6σ_F = 420 km s⁻¹ as illustrated. We apply an observational detector censoring by limiting the velocity window (ΔV) of both I and I_F; we emphasize that the censoring is applied to I_F after performing the convolution. Although we set a smaller OSR than provided by our own real galaxy data, ΔV ≫ ΔV_F such that this idealized situation is now synonymous with a natural broadening of the template spectrum that is subsequently censored by "observation" and relevant to the practicalities of our galaxy observations. The subsequent calculation of $X^{\prime }_{I^{\prime }}$ and A'_F demonstrates significant differences due to the detector censoring. The way in which the detector censoring is applied to the data invalidates the commutation applied to Equation (2) to obtain Equation (3). One can also think of this in terms of the information content: the light gray regions in Figure 1 contain elements of I used to create I_F that remain in (I_F)', but these components are lost in the "detector censoring" such that they are not available to I', leading to a fundamental difference in the information content of $X^{\prime }_{I^{\prime }}$ and A'_F.

In Figure 1, we overplot A'_F and $X^{\prime }_{I^{\prime }}$, demonstrating the significant difference between the two, particularly toward positive velocities. The difference in the structure of $X^{\prime }_{I^{\prime }}$ at positive lags is mostly due to the red edge of the velocity window we have applied. This was done by design to demonstrate a strong example of the detector-censoring effect. Because the velocity window edge is placed very near the reddest line in the Mg i triplet, the detector censoring causes significant changes to the equivalent width of the line in (I_F)' relative to the expectation from I'. However, this example is not contrived; the Mg i lines can be, and indeed are, redshifted near the edge of our OSR for some of the galaxies observed in the DMS.

Thus, having established that Equations (2) and (3) are not equally good approximations when applied to real data, it is important to understand the significance of this difference on the measured stellar kinematics. We illustrate this in Figure 1(c) by defining a broadening function F' such that $A^{\prime }_{F^{\prime }}\approx X^{\prime }_{I^{\prime }}$, where the difference between F and F' is due to detector censoring. Given that the full width at half-maximum (FWHM) of $X^{\prime }_{I^{\prime }}$ is less than A'_F and that $X^{\prime }_{I^{\prime }} - A^{\prime }_F$ is asymmetric about V = 0, F' will have a small velocity shift and a systematically lower velocity dispersion when compared to F. To illustrate this, we produce two $A^{\prime }_{F^{\prime }}$ where we introduce a 10% ($\sigma _{F^{\prime }}=0.9\sigma _F$) and 20% ($\sigma _{F^{\prime }}=0.8\sigma _F$) systematic error in the velocity dispersion. The resulting functions suggest that the systematic error in σ will be of order 10%, with the systematic error in velocity being much less substantial.

It is important to note that, depending on the implemented procedure, this systematic error is not necessarily evident in simulations testing the recovery of input kinematics when using the nominal CC method. Without any detector censoring, the commutation applied to Equation (2) to obtain Equation (3) is perfectly valid, even if the spectra are not known at all wavelengths, because the information content is identical. Thus, tests of the CC method following Equation (3) should find no systematic error as long as the synthetic galaxy spectra are created using identically the same template spectrum used to create the fitted A_B function. Indeed, we found this to be true in our own preliminary implementations of the nominal CC method. The effect demonstrated in Figure 1 is not reproduced because of an incorrect assumption in these tests; the effect is a fundamental consideration when performing the CC method because spectra are detector-censored by the very act of observation. Therefore, fits to observed galaxy data should reflect this detector censoring in a similar way.

The detector censoring of data also effects the calculation of T_B due to convolution errors near the OSR boundaries of T. Indeed, if one does not properly account for these effects, systematic differences between G and T_B toward the OSR edges will effect systematic errors in the measured kinematics; this is also true of direct-pixel-fitting methods. Thus, the difference in the applicability of the approximations in Equations (2) and (3) are not automatically solved by preferring the procedural implications of the former over the latter. This preference is instead dictated by the ease with which one can mimic the effect of observation on the spectrum fitted to the observed galaxy spectrum and account for boundary effects on the convolution of the template with the broadening function. It is important to note that neither of these effects can be removed from the nominal CC analysis by applying a data-windowing (or apodization) function that downweights the contribution to the CC of pixels near the edges of the spectral window. The effects can be mitigated in this way, but not removed. Moreover, we make the point, as was implicit in our placement of the censoring window, that the significance of these effects is dependent on ΔV, ΔV_F, and the detailed structure of the spectra being considered. That is, one expects the significance of these effects to be correlated with the ratio ΔV_F/ΔV and the equivalent width within ΔV_F at either end of the OSR.

The solution we have taken to these issues is to treat the fitting function the same as is presented by our observing technology. That is, using Equation (2) we first censor T_B, creating (T_B)', in an identical way to how our observations censor G, along with eliminating regions near the OSR edges of both (T_B)' and G that suffer from boundary errors in the convolution. Thus, we implement a more symmetric comparison of the two spectra G and (T_B)'. In this sense, it is beneficial to observe template spectra with a larger rest OSR that accommodates the OSR of the galaxy spectra. As galaxy-continuum measurements are much more observationally expensive than template observations, it is preferable to do this instead of cutting out portions of the galaxy spectrum, if possible. However, such an approach complicates the detailed matching of, e.g., the resolution between the template and galaxy spectra given the different observational setups entailed. Until this is more practically accomplished, we instead further limit the OSR of the galaxy spectrum to accommodate the combination of the template OSR and the velocity width of the fitted convolution function.

2.3. Masking

The spectra of star-forming spiral galaxies (as observed by the DMS) have nebular emission lines. Moreover, observational artifacts are likely to exist such as poorly subtracted sky lines (especially with data taken in the region surrounding the Ca ii triplet), cosmic ray events that erroneously pass nominal removal algorithms, and telluric absorption features. None of these features (nebular emission lines nor observational artifacts) will be present in our stellar template spectra, at least not in the same sense as they are in the galaxy spectra. When fitting in direct pixel space, the effects of these spectral features on the results of the fitted kinematics are removed by masking pre-defined spectral regions from consideration in the fit-minimization statistic.

We wish to maintain the same masking abilities in our CC analysis. Although cast in the light of arguments concerning the edges of the OSR for G and T, the situation regarding lost information is only exacerbated by allowing for multiple mask regions interspersed throughout the OSR, unless all information available to G is also available to T. This is another key argument in our preference of a CC approach that first creates the detector-censored broadened template spectrum, (T_B)', in using Equation (2) to solve for B. By first creating T_B and then (T_B)' using the masks and OSR available to G, we avoid the complexities of performing equivalent operations through Equation (3) and account for the detector-censoring issues above. Thus, the process of masking spectral regions in our CC-fitting method becomes as uncomplicated as that used by direct-pixel-fitting methods.

Dalle Ore et al. (1991) have implemented a scheme that used multiple observations to replace spectral regions with artifacts (e.g., cosmic rays) by those without, and used regions of standard galaxy spectra to replace emission regions. Whereas analyses without masks or with a linear interpolation across masked regions were poor, Dalle Ore et al. (1991) report that such a splicing approach provided robust results. In our own analysis, we have implemented a scheme where we use the velocity offsets of galaxy spectra taken within a given radial range to remove sky lines and thereby create a relatively clean, azimuthally averaged spectrum (Paper II); however, for our SVE analysis, we wish to consider each fiber spectrum individually. Similar to the Dalle Ore et al. (1991) approach, we have considered a scheme where, instead of splicing in other galaxy spectra, we would replace undesired spectral regions with the appropriately convolved regions from the template spectrum. However, by making any non-intrinsic substitutions to the galaxy spectra, we effectively dilute the signal we are attempting to extract by adding (primarily) systematic error to our calculation. Instead, the scheme we apply in fitting the kinematics to each individual fiber spectrum involves defining the set of spectral regions to omit in the reference frame of either G or T, determining the equivalent regions in the opposite reference frame, and blanking the appropriate spectral regions in each reference frame.

Figure 2 demonstrates the use of our masking procedure in both stellar-continuum regions used in the DMS (Paper I) for the high-surface-brightness galaxy UGC 6918. One can imagine simply selecting the middle portion of the Mg i-region spectra to avoid contamination from the [O iii] and [N i] lines. However, such a scheme is impossible in the Ca ii region; the shear number of sky-line residuals obviates significant, uncontaminated portions of the spectrum. These two spectral regions also demonstrate the use of masks identified in the template reference frame (Mg i) and in the galaxy reference frame (Ca ii). Each set of masks are Doppler shifted to the opposite reference frame using the fitted velocity shift between the spectra (see the description in Section 3 of the iterative approach employed here). Additionally, the spectra are automatically masked such that the rest-frame spectral edges toward the blue and red are identical.

Zoom In Zoom Out Reset image size

When defining sky-line mask regions, the region masked should not be dependent on the fitted kinematics except to define the corresponding regions in the template frame. However, one could, for example, desire to mask the Mg i lines and define the mask in the template frame. Given the lines are not only redshifted but also broadened in the galaxy spectrum, we allow such masking regions to be broadened appropriately. This is demonstrated in Figure 2 for the masks applied to the [O iii] and [N i] lines. We expect a larger stellar velocity dispersion than gas velocity dispersion such that the resized masking regions should appropriately mask all of the nebular emission.

As shown in Figure 2, the application of the mask in the region surrounding Mg i is critical to a determination of the stellar kinematics. The comparison of the unmasked and masked CC functions shows this definitively. The large oscillations in the unmasked CC function toward negative velocity shifts are due to the correlation of the [O iii] line with the template stellar-absorption lines, which completely dominates any correlation signal from the stellar-absorption lines common to both the galaxy and template spectra. The difference in the masked and unmasked CC functions for the red spectra is less pronounced but still evident: a visual comparison shows that the poorly subtracted sky lines contribute significantly to the noise in the unmasked CC function. Therefore, while the masks have less effect in this spectral region, their application is no less critical to the robust determination of the stellar kinematics.

The application of our spectral masks forces pixels within the mask region to have zero flux. Given our preparation of the spectra before correlation, pixels with zero flux contribute nothing to the CC function. However, masking in this way can create undesired effects; we discuss these and other concerns in the next section.

2.4. Other Concerns Specific to Fast Fourier Transforms

2.5. Summary

By preferring Equation (2) over Equation (3), and by incorporating the fitted velocity into the determination of the mask transcription between the template and galaxy reference frames, the procedure followed by our fitting algorithm deserves some comment.

When fitting X by X_T, we do not consider the full velocity range available to the spectra in assessing the goodness of fit (χ²_X). As discussed above, we instead directly fit only the peak of the CC function, which is preferentially dominated by kinematic information than is its structure at large velocity lag. However, we find the calculation of the full velocity range of X and X_T useful for assessing template mismatch (Paper II). Fitting the CC function over the limited velocity range represented by the peak is also beneficial to the calculation of the error in X. Following Statler (1995), we define the covariance matrix in X to be

$$\begin{equation} (\delta X_{j,k})^2 = \sum _{n=1}^N T_{n-j} T_{n-k} (\delta G_n)^2, \end{equation} \tag{ 4 }$$

where N is the number of pixels in G and δG_n is the error in the measurement of the galaxy spectrum at pixel n. Note that this equation does not take into account any errors in T; instead, we follow Statler (1995) by designating errors in T as a contributor to template mismatch. Therefore, the calculation of

$$\begin{equation} \chi _X^2 = [X-X_T] [(\delta X)^2]^{-1} [X-X_T]^{\dagger } \end{equation} \tag{ 5 }$$

requires the inversion of the covariance matrix [(δX)²], once for every instance of the CC function, where † denotes the transpose operator.⁶ By only considering a window surrounding the CC peak, we can greatly reduce the computational expense of this calculation. There are varied approaches to the definition of this window in the literature; we define it by a multiple of the CC-peak FWHM and center it on the CC peak velocity. Although not presented in detail here, we find an optimal window size of approximately twice the FWHM of the CC peak based on minimizing both the systematic error in the kinematics due to the window and the error calculation time; this is the definition of the window we take throughout this document. Individual determinations of B minimize Equation (5) for a single calculation of X and [(δX)²]⁻¹ over this fitting window.

For all fittings discussed herein, we employ a Gaussian broadening function with the nominal three parameters: a normalization, a; a peak velocity, V; and a velocity dispersion, σ. The determination of B ostensibly allows for any value of σ; however, discrete convolutions are systematically in error when B is undersampled. We have run tests demonstrating that robust convolutions require that, when B is defined by a Gaussian, it must have a velocity dispersion that is no less than 0.85 pixels (in which case the Gaussian FWHM is no less than 2 pixels). Thus, whenever this (Nyquist) criterion is not met, we block-replicate the template spectrum—dividing each original pixel into a set of N pixels with the same value as the original pixel and a velocity width that is a factor of N smaller—to decrease the velocity width of each pixel. This avoids the undersampling of the Gaussian function and produces the correctly convolved spectrum when it is block-averaged—averaging the values of N pixels into a single pixel having a factor of N larger velocity width—back to the original sampling. We impose a maximum block-replication of 8 pixels for every one input pixel; Gaussian functions attempted during the fit minimization with σ below that allowed by this maximum block-replication factor are converted to delta functions.

Our fitting algorithm runs through three iteration tiers. The first iteration tier ensures the settling of the fit into a global—or, at least, a large local—minimum in χ² space. This is accomplished by starting the fit with an initial guess set of parameters that are then altered to produce the minimum χ²_X following the Levenburg–Marquardt minimization algorithm discussed in Press et al. (2007, Section 15.5.2). The derivatives of the fitting function X_T required by this minimization routine are determined via a finite-differencing method. Once the fitting algorithm minimizes χ²_X based on the initial guess parameters, we restart the minimization from a new set of guess parameters that are randomly distributed using a characteristic scale for each parameter (see Equation (10.5.1) of Press et al. 2007). For a Gaussian function, we choose characteristic scales of 1, 100 km s⁻¹, and 100 km s⁻¹ for a, V, and σ, respectively. The scale for a has been chosen based on input template and galaxy spectra that have had their mean value normalized to unity such that the normalization of B should also be of that order. The scales chosen for V and σ are based on the expectations of our survey. We find that our fitted kinematics are robust against the choices for these scales, meaning that the χ² landscape for our analysis is well behaved. This statement becomes less accurate for very low S/N spectra. First-tier iterations are typically performed five times and the best resulting χ²_X sets the final fitted parameters.

The second iteration tier adjusts the mask transcription function, which incorporates both the fitted V and σ, depending on the feature being masked (see Section 2.3). The adjustment of the mask requires recalculation of X and [(δX)²]⁻¹ and is, therefore, not included in the first-tier iterations due to the significance of the added computations. Each second-tier iteration updates the mask transcription function after completing a first iteration tier and then restarts the first iteration tier. Second-tier iterations continue until the fitted kinematics and the mask transcription parameters are marginally different (<0.1 pixel).

The third iteration tier incorporates an adjustment of the continuum between the template and the galaxy spectrum. It is often the case that the template spectrum will have a different continuum shape than that of the galaxy. In the spectral region near 5150 Å, for example, the flux contribution for UGC 6918 is divided roughly equally between A and K spectral types (Paper I). However, in terms of the line structure, the K1 III template provides the best fit to the galaxy spectrum (Paper II), meaning that the largest differential contribution of the A-star template is the continuum shape. Our "continuum-fitting" algorithm works to eliminate such relative continuum differences by fitting a low-order Legendre polynomial to the difference spectrum G − T_B. The purpose of this iteration tier is to force the continua of the two spectra to be minimally different such that B is, to first order, only affected by the high-frequency (line) fluctuations, as opposed to low-frequency (continuum) fluctuations. Again, because our goal is to measure B, we need only to "correct" the continuum of the template (or galaxy) spectrum in a relative sense. Each third-tier iteration uses the derived T_B from a completed second iteration tier to adjust the continuum and then restarts the second iteration tier. Third-tier iterations continue until the coefficients of the fitted polynomial are all marginally different from zero. Note, however, we employ a hard limit of ten iterations for both the second and third iteration tiers, such that fits to low-S/N spectra eventually converge or are considered non-convergent under the second- and/or third-tier criteria.

The first iteration tier is typical of previous implementations of CC analysis and, indeed, proper use of χ²-minimization algorithms; and the third iteration tier is an alteration to concepts previously accomplished by, e.g., applying a high-pass filter to the input spectra before using CC analysis to extract stellar kinematics (e.g., Statler 1995). The second iteration tier encompasses the primary novelty of our algorithm. It is also worth pointing out that, even though we only consider the CC peak in the first iteration tier, the formalism in Section 2.1 demonstrates that this peak considers all spectral information. Moreover, the third iteration tier considers the full wavelength range in adjusting the relative continuum shape between the template and galaxy spectra, which can have a significant effect on the kinematics. Finally, our assessments of template mismatch (Paper II) are critically dependent on considering the full velocity range in the CC function. Thus, no spectral information is discarded in our CC approach to extracting and assessing the stellar kinematics in DMS spectra.

As a by-product of the Levenburg–Marquardt algorithm, we use the covariance matrix to determine the formal parameter errors. These errors correlate well with the empirical errors generated from Monte Carlo simulations we describe in the next section. However, for this to hold true, one needs good estimates for the spectral errors. Our method for calculating spectral errors is provided in Westfall (2009).

We have performed a set of Monte Carlo (MC) simulations as an assessment of the ability of our software to robustly determine stellar kinematics in galaxies. For each simulation set, we define a value for 〈S/N〉—pixel⁻¹ with a pixel size of 7.6 km s⁻¹ and a spectral resolution element of ∼26 km s⁻¹—and a broadening function that are used to create a synthetic galaxy spectrum.⁷ Each realization of a synthetic galaxy spectrum is produced by convolving a K1 III template by a known broadening function and then adding normally distributed noise following a wavelength-dependent S/N function. Each simulation set consists of 50 such realizations. A previously seen example of a synthetic galaxy spectrum without noise is shown as spectrum I_F in Figure 1.

A wavelength-dependent S/N function is adopted in order to mimic the effects of the spectral vignetting function present in our galaxy spectra. To create this function, we average S/N(λ) across all fibers in the SparsePak IFS of UGC 6918 and fit it with a low-order polynomial. The effect of the vignetting function is thus shown to produce S/N ∼ 0.6 at the edges of the OSR with a maximum of S/N ∼ 1.3 at ∼5150 Å for a spectrum with a mean S/N of unity. For each MC simulation, we renormalize this functional representation of S/N(λ) such that each realization of a synthetic galaxy spectrum in the simulation has the appropriate 〈S/N〉 pixel⁻¹.

Our CC simulations explore all combinations of the following set of parameters:

$$\begin{eqnarray*} \langle {\rm S/N}\rangle & = & 1, 1.4, 2.0, 2.8, 4.0, 8.0, 16.0, 64.0 \nonumber \\ a_m & = & 1 \nonumber \\ V_m & = & 1000, 2500, 3600, 4160, 4640\ {\mbox{{km s$^{-1}$}}} \nonumber \\ \sigma _m & = & 5, 10, 20, 40, 80, 110\ {\mbox{{km s$^{-1}$}}}. \nonumber \end{eqnarray*}$$

Here, a_m, V_m, and σ_m represent the parameters for the Gaussian-model broadening function, B_m. The parameter space covered by these simulations is somewhat more limited than, but overall representative of, that which is available to the DMS. The range probed in 〈S/N〉 and σ_m reflect the ranges we expect for our galaxy data; although azimuthally averaged spectra may have greater S/N than represented. The chosen velocity shifts are selected to demonstrate any effects of the detector censoring for the minimum and first-quartile velocities of galaxies observed within the DMS (∼1000 and ∼2500 km s⁻¹), the most extreme velocity still containing all three of the Mg i triplet lines (∼3600 km s⁻¹), and those velocities retaining only two (∼4160 km s⁻¹) and one (∼4640 km s⁻¹) of those lines. Thus, we have performed 240 CC simulation sets for the analysis presented herein.

The limited representation of these simulations with respect to all the applications foreseen in the DMS deserves further discussion. Foremost among these limitations is that these simulations will provide no assessment of the effects of template mismatch either in the continuum or in the detailed line shapes/depths. Additional, more targeted simulations have been presented in Paper II as an exploration of the results for different templates and the Ca ii spectral region. For the present study, our simulations represent a scenario where we, in fact, have the ideal spectrum I—a practical improbability for our observed galaxy data. As such, these simulations represent a best-case-scenario in terms of the systematic errors we expect in B. Moreover, we only consider the Mg i region; however, we expect the results to be largely transferable between the two stellar-continuum wavelength regions and two IFS instruments (Paper I) used by the DMS. Finally, since our template observations are taken within the same observed wavelength region as the galaxy data, the true galaxy data actually have an extension of their rest-frame OSR toward the blue. In creating our synthetic galaxy spectra, we cannot reproduce this portion of the spectrum and simply supply a constant continuum value. Note, however, that we do limit the red end of the synthetic galaxy spectra to the observed wavelength range represented within our data. Having described these inaccuracies, it is worthwhile to note that such a scenario as provided by these simulations can be beneficial in establishing the systematic errors solely due to spectral noise. That is, using these simulations and an estimate of ∼4% for the template mismatch (Paper II), we can produce a transitionary S/N limit where the systematic errors due to spectral noise become comparable to those expected for template mismatch.

We present individual examples of synthetic G, fitted T_B, synthetic X, and fitted X_T for a few hand-picked simulations in Figure 3. These simulations are chosen as those providing reasonable results for B and, as such, are not necessarily "representative" of all the simulations with the same parameter set; we consider statistics of the ensemble of results below. Note that we differentiate between the "model" and "fitted" Gaussian parameters using the subscripts m and f, respectively. Figure 3 demonstrates that even at the lowest 〈S/N〉, one might expect to find a reasonable measurement of B. In fact, the systematic and random errors for the 〈S/N〉 = 1.0 example are comparable to the 〈S/N〉 = 16.0 example; the difficulty in the determination of B for the latter is the rather large σ_m and the censoring of the Mg i triplet by the velocity shift. Also note that, with a pixel size of 7.6 km s⁻¹ and a spectral resolution of ∼26 km s⁻¹, the CC algorithm is able to recover the 5 km s⁻¹ broadening to ∼14% even at 〈S/N〉 = 2.8; however, the random error in this case is >50%.

Zoom In Zoom Out Reset image size

Figures 4 and 5 provide statistical quantities for each of our 240 simulation sets. We calculate the mean and standard deviation (ε) of each Gaussian parameter after rejecting strong (3ε) outliers; only at 〈S/N〉 < 3 does the number of simulations used in the statistics, N_sim, become significantly different from the total number of simulations (50). Our assessments of the simulation sets demonstrate negligible influence of the velocity shift on the results. We expected this result specifically because of our approach to the CC method (Section 2.2), and we continue by largely omitting further discussion of velocity-offset effects.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 presents our assessment of the systematic errors in our kinematic measurements. As expected, the accuracy of the resulting fits is primarily a function of the spectral S/N. We find ≲10% systematic error in σ_f for all simulations when 〈S/N〉 ≳  4. Below this limit, the returned accuracy of σ_f becomes dependent on the input σ_m, as demonstrated in Figure 4 by the different line gray scale. In fact, the results exhibit an inverse correlation between the model and fitted values for both a and σ. There appears to be no such trend in V, at least consistent between all simulations; however, we find that the systematic differences (whether positive or negative) in V become more pronounced for models with higher σ_m. These trends make some intuitive sense in that, as the noise increases, the CC algorithm begins to become biased by the correlation between noise peaks. The effect is strongest for models with broader lines as the difference between their inherent widths and the width of a pixel (the natural correlation width of a noise peak) is larger than models with velocity dispersions of order the size of a pixel. We do not plot the random errors in Figure 4 for clarity, but we find that the random errors are substantially larger—factors of >2—than the systematic errors at all 〈S/N〉.

We plot the standard deviation, ε(x) for the variable x, in each simulation set against the input σ_m and 〈S/N〉 in Figure 5 as an assessment of our random errors. These data reinforce the point that the quality of the fitted kinematics is determined primarily by 〈S/N〉 and secondarily by σ_m. There is a strong correlation in both ε(V_f) and ε(σ_f)/σ_m with 〈S/N〉. Moreover, there is a direct correlation between the model value σ_m and ε(V_f). Finally, there appears to be a critical value where σ_m switches from a positive to negative correlation with ε(σ_f)/σ_m. For most of the simulations, these two quantities are anti-correlated; however, within the range 20 km s⁻¹ < σ_m < 40 km s⁻¹, ε(σ_f)/σ_m changes from an inverse to a direct correlation at small σ_m to large σ_m, respectively. This observation is strengthened by the fact that all simulations done with 〈S/N〉>1 demonstrate this effect; however, aside from simulations with σ_m = 5 km s⁻¹, the interdependence of ε(σ_f)/σ_m and σ_m is rather weak. As of now, the origin of this transition is unknown; it could be an observational effect due to sampling of σ_m by the spectral resolution or it could be due to the detailed line structure in the K1 III template. Although it may be possible to differentiate these effects with further analysis, it is not pertinent to our basic quantification of the random and systematic errors in our kinematic quantities as a function of 〈S/N〉—the fundamental goal of these MC simulations.

We provide here a list of general results drawn from Figures 4 and 5:

• 1.  
  Statistically aberrant fits for spectra that have no convergent kinematic solution are most probable when σ_m ≳ 80 km s⁻¹ and 〈S/N〉 ≲3; otherwise, we expect the probability of such spectra to be ≲10% at 〈S/N〉 ≳ 1.5.
• 2.  
  The quality of the kinematic measurements is largely independent of the velocity shift in the range V_m = 1000–5000 km s⁻¹, as sampled by our simulations.
• 3.  
  When 〈S/N〉 ≳ 2, systematic errors in V_f are always less than 15 km s⁻¹ and the mean systematic error in V_f of all simulations with the same 〈S/N〉 and σ_m are always less than ∼5 km s⁻¹.
• 4.  
  Systematic errors in σ_f are typically less than 20% for all simulations and ≲5% for means binned by input σ_m when 〈S/N〉 ≳ 2, except for the largest input velocity dispersion (σ_m = 110 km s⁻¹).
• 5.  
  Random errors in each fitted quantity (a_f, V_f, and σ_f) dominate over systematic errors by factors of two or more at all 〈S/N〉.
• 6.  
  We expect random errors of ≲10 km s⁻¹ in V_f when 〈S/N〉 ≳8.0 and σ_m ⩽ 110 km s⁻¹; this degree of error only occurs in spectra with σ_m ⩽ 20 km s⁻¹ at 〈S/N〉 ≲ 1.0.
• 7.  
  Random errors in σ_f can be as much as 100% when 〈S/N〉 = 2 and σ_m ≈ 5 km s⁻¹; however, this decreases to 20%–40% when σ_m ≳ 10 km s⁻¹.

Within the larger context of the DMS, we expect the majority of our spectra to have intrinsic values of 20 km s⁻¹ ≲σ_LOS ≲ 40 km s⁻¹; therefore, we expect random errors of ≲20% and, at most, 10% systematic error when 〈S/N〉 ≳ 3 pixel⁻¹ (the pixel size for SparsePak is 7.6 km s⁻¹ in the Mg i spectral region). We will not preclude ourselves from fitting kinematics at lower S/N given that formal errors in this regime tend to be larger than the systematic error; however, analysis of such results will be tempered by the expectations from these simulations. Moreover, the measured velocities should show only very modest systematic errors even at 〈S/N〉 ∼ 1.5.

There are only three galaxies observed for the DMS by SparsePak alone with velocities outside the range probed by these simulations (their respective V_sys values are 4856, 7152, and 9392 km s⁻¹). At the maximum of these redshifts, the rest wavelength overlap between a template star and the galaxy spectrum is ∼100 Å and the strong Mg i triplet will have been shifted out of the OSR. Thus, measurements of B for these galaxies will depend on the other (mostly Fe lines) present in the spectrum. We expect the S/N of these measurements to be lower due to the decreased total and mean equivalent width of these lines; however, there should be substantial kinematic information in these spectra, even given the limited overlap. Note that the bluest Mg i line is still within the OSR for PPak at redshifts of z ∼ 0.04, which is very near the upper limit of our survey. We will comment on the effects of limited rest wavelength overlap effects on our kinematics where appropriate in future papers; however, the simulations presented here predict this to be of little concern.

Below we analyze SparsePak IFS for UGC 6918 with 2 ≲ 〈S/N〉 ≲ 30 for all individual spectra that provide correlation peaks within a limited range of the expected systemic velocity. The fitted kinematics to these spectra demonstrate velocity dispersions within the range 10 km s⁻¹ ≲ σ_m ≲ 70 km s⁻¹. Therefore, Figure 4 predicts that the ≲4% mismatch with respect to the K1 III template (Paper II) should dominate the systematic error in our stellar kinematics. Moreover, Figure 5 predicts errors of ε(V_LOS) ∼ 1–2 km s⁻¹ and ε(σ_LOS) errors of ≲2% and ∼70% at, respectively, large and small intrinsic σ_LOS and 〈S/N〉. These estimates are determined by approximating that 〈S/N〉 ∝ σ_LOS and that σ_LOS is inversely proportional to the galaxy radius, R, as roughly appropriate for an exponentially declining dispersion. We test these predictions in Section 5.

4.1. Comparison with pPXF

We directly compare the internal test of our CC software (Detector-Censored Cross-Correlation, DC3) above with identical MC simulations performed using the pPXF software⁸ presented by Cappellari & Emsellem (2004). These MC simulations are necessary for providing a quantitative comparison between DC3 and pPXF and substantially augment the MC simulations performed by Cappellari & Emsellem (2004). Specifically, we expand the S/N regimes from 60 and 600 to eight values increasing logarithmically from 1 to 64, the velocity offsets from zero to a range of 1000–4640 km s⁻¹, and the velocity dispersion regimes from 48–360 km s⁻¹ (0.8–6.0 pixels) to 5–110 km s⁻¹ (0.7–14.6 pixels). To provide the most robust comparison between the two programs, we fit exactly the same synthetic spectra using pPXF as was done in the previous section using DC3. The mode in which we have run pPXF is designed to follow closely that of our CC approach in that we apply the following specifications: we adopt no penalization; we fit only the first and second moments (V and σ) of the LOSVD; we allow for a third-order polynomial fit to the continuum (where a constant is of zeroth order); we mask the same [O iii] and [N i] spectral regions as done in the CC method; and we fit only the K1 III template. Given that both our simulations above and the simulations we perform with pPXF do not allow for template mismatch, there is, ideally, no continuum differences between the synthetic galaxy spectra and the fitted template; however, we generally determine continuum adjustments in the nominal approach to fitting stellar kinematics for our real galaxy data. Therefore, the pertinent statistics in the recovery of the LOSVD should include the errors added by the uncertainty in the continuum.

Figure 6 compares the results of the MC simulations run using pPXF and DC3 in terms of both systematic and random errors. The overall conclusion given by this figure and the discussion below is that DC3 generally performs as well as pPXF, within the confines of these tests. Thus, in DC3, we have addressed the primary concerns with the CC method that have been invoked to explain the preference for direct-pixel-fitting methods. We note some features of our comparison between the DC3 and pPXF results in more detail below.

Zoom In Zoom Out Reset image size

We find differences in the systematic error returned from pPXF and DC3. For the pPXF simulations, we find the velocity systematic errors have a minimum range of 〈V_f〉 − V_m ∼ ±1 km s⁻¹ (∼0.1 pixels; Figure 6(a)) and the velocity dispersion has a ∼25% systematic error in 〈σ_f〉 for the σ_m = 5 km s⁻¹ case (Figure 6(f)). We emphasize that the synthetic galaxy spectra used here were created using DC3 sub-routines as explained above, such that these systematic effects do not necessarily represent failings of the pPXF software, but serve more as a warning when attempting to measure kinematics with undersampled broadening functions. That is, these systematic errors are likely due to numerical differences in the convolution algorithm performed by pPXF and that used by DC3 (Section 13.1 of Press et al. 2007); it is natural for DC3 to perform better given its use of exactly the same algorithm to fit the synthetic spectra as was used to generate them. Measurements of σ_f ≲ (0.85ΔV_p = 6.4) km s⁻¹  where ΔV_p is the pixel velocity width, become strongly biased toward larger velocity dispersion because the LOSVD is undersampled (i.e., below Nyquist). Both pPXF and DC3 employ a block-replication algorithm in an attempt to allow for measurements of σ_f below this limit (Section 3); however, differences in the implementation may be the cause for the systematic error seen in the pPXF results. The regimes in which these systematic differences exist are relatively uninteresting in terms of application to our DMS spectroscopy. Also, the trend in the systematic error in σ_f noted in the previous section for DC3—the bias of σ_f toward the size of a pixel at low 〈S/N〉—does not appear to be as strong for pPXF (compare panels (f) and (g) in Figure 6).

For both programs, the mean formal errors, 〈δV_f〉 and 〈δσ_f〉, calculated via the covariance matrix provide good measurements of the error when compared to the standard deviation in the returned parameters, ε(V_f) and ε(σ_f), at 〈S/N〉 ⩾ 4 (panels (c), (d), (h), and (i) in Figure 6). Below this 〈S/N〉, the pPXF calculations of 〈δV_f〉 and 〈δσ_f〉 at low σ_m can be substantial overestimates of the standard deviation derived from the MC, which is not reflected in the DC3 simulations. By comparison, DC3 provides formal errors that are good to better than 50% for the vast majority of simulations at all 〈S/N〉; the exceptions for V_f are in the highest σ_m cases where the formal errors tend to be under estimates when 〈S/N〉 < 4, and the exceptions for σ_f are in the highest and lowest σ_m cases when 〈S/N〉 ⩽ 2.

A comparison of the errors in the fitted parameters as returned by pPXF and DC3 is done by showing the ratio of the standard deviation in both V_f and σ_f in panels (e) and (j) of Figure 6. In this sense, we again find the performance of the two programs to be very comparable. In detail, the pPXF routine obtains, on average, a ∼5% lower error in the parameters; however, this difference is marginal when considering slight differences in the wavelength regions incorporated during the fitting procedure between the two programs. We have attempted to make the masked regions identical; however, they will inevitably be slightly different due to the details of the masking implementation of the two routines.

We have used DC3 to fit SparsePak IFS of UGC 6918 obtained during the pilot phase of the DMS. These data will be fully described and used in forthcoming work to measure disk and halo dynamical properties in this galaxy. Here, we fit the stellar kinematics using only a K1 III template as it minimizes template mismatch for this galaxy (Paper II), and is simple and easily reproducible. We emphasize that the use of a single stellar template here is not a requirement of DC3, but merely a convenience for clarity of this illustration.

Ignoring the results of our MC simulations in Section 4 for the moment, we attempt to fit any CC function with a peak within a few hundred km s⁻¹ of the systemic velocity of UGC 6918 (∼1109 km s⁻¹) regardless of the S/N. For all fits, a Gaussian function is used to parameterize B and a third-order Legendre polynomial is used to adjust the continuum difference between the broadened K1 III template, T_B, and the fitted galaxy spectrum in each fiber (see Section 3). Additionally, wavelengths surrounding the [O iii]λ5007 and [N i] (λ5198 and λ5200) nebular emission lines (as in Figure 2) are masked. The fitting window, where X is compared with X_T in a χ² sense, is centered on the CC peak and given a width of twice the CC-peak FWHM. Figures 7 and 8 present the results for two example spectra, one at high 〈S/N〉 and another at low 〈S/N〉. The low-〈S/N〉 spectrum shown here represents a spectrum near the low-〈S/N〉 limit where our fitting algorithm was able to converge to measurements of V_obs and σ_obs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 7, we find excellent agreement between X and X_T over the full velocity range. However, in comparing G and T_B, we note the applied [O iii] mask does not properly block the blue wing of the line profile. This artifact is easily removed by enlarging the mask of the [O iii] line; however, its presence provides a useful, albeit exaggerated, demonstration of the template-mismatch information available in the full CC function, outside the central peak. The effects of the unmasked remnant of the [O iii] are evident in the difference X − X_T; the [O iii] remnant is seen to correlate with the Mg i triplet at ∼−9000 km s⁻¹. Nonetheless, this artifact has a negligible effect on the fitted kinematics: we find the same results within the errors even after doubling the width of the mask. The predicted random errors of ε(V_obs) < 2 km s⁻¹ and ε(σ_obs) ∼ 2% from our MC simulations hold true for this spectrum.

It is difficult to assess any strong, systematic deviations between X and X_T or between G and T_B in the fit to the low-〈S/N〉 spectrum shown in Figure 8; the random noise overpowers these signatures. Even so, the correlation peak is still strong and obvious, and we are able to obtain a set of fitted kinematics. Our MC simulations demonstrate that, for the 〈S/N〉 and σ_obs of this spectrum, we should expect no systematic error in V_obs but a random error of 3–4 km s⁻¹; for σ_obs, we expect systematic errors of ∼10% with random errors of ∼60%. The measured random errors in the fitted kinematics are consistent with these expectations. Indeed, at ∼50%, the random errors are slightly better than we expect suggesting the predictive power of our simulations is good to ∼10%–15%, even in the presence of some, if minimal, template mismatch.

5.1. Comparison with pPXF

We have also fit our SparsePak UGC 6918 spectra using the pPXF software for comparison; the results for the high- and low-〈S/N〉 fibers discussed above are given in Table 1. For a fair and simple comparison, we use the K1 III star as the fitted template for both DC3 and pPXF. For clarity and reproducibility we have intentionally avoided using composite stellar templates, although Paper II outlines how such composites are constructed and applied in the context of DC3. While composites templates more realistically mimic the observed stellar populations, we also have shown in Paper II that we expect ⩽4% systematic error in our DC3 measurements of σ_f for UGC 6918 in this spectral region when using a single K1 III template. From our MC simulations, we expect no strong systematic difference between the results from DC3 and pPXF.

Table 1. Example Stellar Kinematic Fits

	〈S/N〉	$V_{f,{\rm \it DC3}}$	$V_{f,{\rm \it pPXF}}$		$\sigma _{f,{\rm \it DC3}}$	$\sigma _{f,{\rm \it pPXF}}$
Fiber	(pixel−1)	(km s−1)	(km s−1)	ΔVf/δ(ΔVf)	(km s−1)	(km s−1)	Δσf/δ(Δσf)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
52	27.8	1125.3 ± 1.5	1125.2 ± 1.4	0.02	68.3 ± 1.6	63.9 ± 1.6	1.95
55	1.7	1076.4 ± 3.9	1073.6 ± 3.9	0.52	11.2 ± 5.7	11.9 ± 5.9	−0.08

Notes. Columns are: (1) fiber number; (2) mean S/N per pixel; (3) fitted stellar velocity from DC3; (4) fitted stellar velocity from pPXF; (5) difference defined by $\Delta V_f = V_{f,{\rm \it DC3}} - V_{f,{\rm \it pPXF}}$ normalized by the error in the difference; (6) fitted stellar velocity dispersion from DC3; (7) fitted stellar velocity dispersion from pPXF; (8) difference defined by $\Delta \sigma _f = \sigma _{f, {\rm \it DC3}} - \sigma _{f, {\rm \it pPXF}}$ normalized by the error in the difference.

Download table as:  ASCIITypeset image

In Figure 9, we illustrate the difference in all kinematic measurements made by DC3 and pPXF via residual maps of $\Delta V_f = V_{f,{\rm \it DC3}} - V_{f,{\rm \it pPXF}}$ and Δσ_f. These differences increase with radius, as expected from the decrease in spectral S/N, and are incoherent in azimuth. The residuals are generally very small and track the kinematic errors, such that all velocity measurements are within the measurement errors determined by each program. For the velocity dispersion measurements, 82% and 94% of all spectra are within factors of, respectively, one and two times the measurement errors. However, as shown in Table 1, we find a marginally significant (∼2 standard deviations) difference between $\sigma _{f,{\rm \it DC3}}$ and $\sigma _{f,{\rm \it pPXF}}$ for the highest-〈S/N〉 fiber, albeit the difference is small in relative terms (6%). This difference may be due to different systematic effects of template mismatch on the two algorithms. Inclusion of our estimate template-mismatch error does indeed reconcile the two measurements.

Zoom In Zoom Out Reset image size

Thus, as with our comparison of the two algorithms using MC simulations, we find that a direct comparison of DC3 and pPXF using real data demonstrates that the two algorithms produce very comparable stellar kinematics.

Robust stellar kinematic measurements are crucial to the success of the DMS. Indeed, these measurements are the first step toward determining the disk dynamics and, e.g., mass surface density of each galaxy. Thus, the tools we use to perform these measurements deserve the extended discussion provided herein. Our method is generally applicable, and it is directly relevant to our kinematic analysis of the IFS in the DMS. We address complications in CC analysis that have been used as rationale for favoring direct-pixel-fitting methods (Rix & White 1992; Cappellari & Emsellem 2004), such that they should no longer be concerns. Our newly implemented CC software adjusts the nominal algorithm used in previous CC analyses to more properly determine stellar kinematics under the practical inevitability of detector-censored data. Moreover, as demonstrated in Paper II, our approach retains the primary attraction of the CC method, which is to largely consolidate kinematic and template-mismatch information toward different regions of the CC function.

The presentation of our CC approach follows both a mathematical discussion of the assumptions in the method and a practical application of these methods to real data. Our tests of the fundamental assumptions have isolated a ≲10% effect incurred when the broadened stellar template is not treated in a symmetric way with respect to the observed galaxy spectrum. This symmetric treatment is critical to our ability to apply spectral masks in conjunction with the CC approach. We have shown that such masks are necessary given the strong nebular emission lines, sky lines, and observational artifacts (such as cosmic rays) present in our data.

We extensively evaluate the performance—accuracy and precision—of our CC software using MC simulations. As expected, we find that the primary factor in the quality of the determined kinematics is the spectral S/N; the velocity dispersion, σ, is a strong secondary factor. With the caveat that we do not explore template mismatch, two principle results of our MC simulations are: (1) Random errors for spectra with 〈S/N〉 ≳ 3 and 20 km s⁻¹ ≲σ ≲ 40 km s⁻¹ are ≲20% in σ. (2) Systematic errors for spectra with 〈S/N〉 ≲ 2 are typically ≳10% in σ but are always ≲10% in σ when 〈S/N〉 ≳ 4. A detailed comparison with the pPXF software provided by Cappellari & Emsellem (2004) demonstrates comparable results.

Finally, we demonstrate a direct application of our CC software to SparsePak IFS of the stellar kinematics of UGC 6918, as well as a commensurate analysis using pPXF. As with our MC simulations, we find the two algorithms provide comparable kinematic measurements. Thus, we have been successful in addressing concerns over S/N and masking limitations in the CC method, effectively eliminating these concerns as reasons to favor direct-pixel-fitting methods for measurements of stellar kinematics. Our preference for DC3 is rooted the in the fundamental way in which CC methods isolate information pertaining to VDF and template mismatch. The latter is a topic explored further in future work.

It is a pleasure to acknowledge fruitful conversations concerning this work with D. Andersen, R. Swaters, and T. Martinsson. We also thank the referee for many useful suggestions that have increased the clarity of our presentation. K.B.W. and M.A.B. acknowledge support for this work from NSF grants AST-0307417, AST-0607516, and AST-1009491. K.B.W. also acknowledges supplementary support from the Wisconsin Space Grant Consortium and the AAS International Travel Grant. Portions of this work were produced as part of K.B.W.'s PhD thesis undertaken at the University of Wisconsin-Madison. K.B.W. is currently supported by NSF grant OISE-0754437. M.A.B. acknowledges support from the University of Wisconsin College of Letters & Sciences Ciriacko Faculty Fellowship. M.A.W.V. acknowledges travel support provided by the Leids Kerkhoven-Bosscha Fonds.

This appendix details the procedure for correcting velocity dispersion measurements derived via CC for differences in the instrumental broadening between the template and the galaxy spectrum. We ignore beam-smearing effects in order to simplify the following discussion.

From Paper II, we define δσ²_inst ≡ (σ^inst_G)² − (σ^inst_T)², independent of wavelength. However, the proper measurement of σ_LOS via the CC method requires σ^inst_T[λ_T] = σ^inst_G[λ_T(1 + z)] where the observed galaxy-spectrum wavelengths (λ_G) are identical to the Doppler-shifted template-spectrum wavelengths, λ_T(1 + z). That is, even if the template and galaxy observations are taken during the same night with the same spectrograph setup, one expects δσ²_inst ≠ 0 if there is a wavelength dependence of the instrumental dispersion and a non-zero Doppler shift between the galaxy and template. Thus, we not only have to account for gross differences between σ^inst_T and σ^inst_G due to, e.g., slight differences in the spectrograph setup, but also second-order effects incurred by the wavelength dependence of σ_inst.

We assume that the difference in instrumental broadening between G and T has had a negligible effect on the derived velocity between the template and the galaxy, which is a robust assumption for our expected velocity accuracy and small differences between σ^inst_T and σ^inst_G (at most a few km s⁻¹). In this case, only corrections to the measured velocity dispersion, σ_obs, are applied.

We define $\mathcal {T}$ to be the template with an instrumental resolution that ideally matches G. The CC of this template with a broadened version of itself is $\mathcal {X}_{\mathcal {T}} \equiv \mathcal {T}_B \circ \mathcal {T}$, where $\mathcal {T}_B = \mathcal {T} \otimes B$ and B is the correct LOS broadening function. Given a measurement of the Doppler shift (z = λ_G,j/λ_T,j − 1 for spectral feature j), one can create $\mathcal {T}$ by convolving T with a wavelength-dependent broadening function, B_δλ, with a Gaussian velocity dispersion of

$$\begin{equation} \sigma ^2_{\delta \lambda }[\lambda _T] =\big (\sigma ^{\rm inst}_G[\lambda _T (1+z)]\big)^2 - \big(\sigma ^{\rm inst}_T[\lambda _T]\big)^2 + \sigma ^2_{\rm off}, \end{equation} \tag{ A1 }$$

where σ_off is a constant put in place such that σ²_δλ ⩾ 0 at all λ_T. In principle, we need only force σ_off ⩾ 0. However, in practice, we require σ_δλ be large enough to avoid systematic errors in the discrete convolution due to sampling issues. Moreover, the CC method is predicated on a significant broadening between the template and galaxy spectra; therefore, the difference in resolution cannot be so large as to force σ_off to be larger than the astrophysical velocity dispersion in the galaxy spectrum. The latter is very rarely a problem for the DMS and always occurs at low S/N in our IFS; such occurrences will be noted where appropriate in the analysis of DMS stellar kinematics.

Thus, a comparison between $\mathcal {X}_{\mathcal {T}}$ and $\mathcal {X}=G\otimes \mathcal {T}$ requires a priori knowledge of the Doppler shift between G and $\mathcal {T}$. Instead of adding another iteration tier to our CC implementation, we make some assumptions that allow for mathematical and procedural shortcuts. The goal is then to generate a single correction, δσ²_inst, such that

$$\begin{equation} \sigma _{\rm LOS}= \big(\sigma _{\rm obs}^2 - \delta \sigma _{\rm inst}^2\big)^{1/2}, \end{equation} \tag{ A2 }$$

where B has velocity dispersion σ_LOS and the "observed" broadening function, B_obs, has velocity dispersion σ_obs. The measurement of σ_obs is that determined without correcting for differences in instrumental resolution; that is, B_obs is determined by minimizing G○T − (T ⊗ B_obs)○T as described in Section 3. Integrating the intensity-weighted effects of the instrumental dispersion differences on the CC function, we also define

$$\begin{equation} \sigma ^2_{\delta \lambda } = \delta \sigma _{\rm inst}^2 + \sigma ^2_{\rm off}. \end{equation} \tag{ A3 }$$

We calculate δσ²_inst by defining a fourth broadening function, $\mathcal {B}$, that minimizes

$$\begin{equation} (\mathcal {T} \otimes B_{\rm obs}) \circ T - (T \otimes \mathcal {B}) \circ T \end{equation} \tag{ A4 }$$

using our approach described in Section 3; the velocity offset of $\mathcal {B}$ is fixed at zero and we allow no continuum-fitting iterations. The minimization avoids complications from template mismatch and shot noise, as opposed to one incorporating G.

Using the definition $\mathcal {T}=T\otimes B_{\delta \lambda }$, we can rewrite the result of minimization of Equation (A4) to show

$$\begin{equation} \left(\left[T \otimes B_{\delta \lambda }\right] \otimes B_{\rm obs}\right) \circ T \approx (T \otimes \mathcal {B}) \circ T, \end{equation} \tag{ A5 }$$

or, after applying some commutations,

$$\begin{equation} A \otimes (B_{\delta \lambda } \otimes B_{\rm obs}) \approx A \otimes \mathcal {B}, \end{equation} \tag{ A6 }$$

where A ≡ T○T is the AC of the template spectrum. Such commutations are accurate given the symmetric treatment of the parenthetical terms on both sides of Equation (A5). Given that the convolution (or correlation) of two Gaussian functions produces a third Gaussian function with a dispersion equal to the quadrature sum of the two input Gaussian functions, Equation (A6) yields

$$\begin{equation} \sigma ^2_{\mathcal {B}} = \sigma ^2_{\delta \lambda } + \sigma _{\rm obs}^2. \end{equation} \tag{ A7 }$$

Substitution into Equation (A3) yields

$$\begin{equation} \delta \sigma _{\rm inst}^2 = \sigma ^2_{\mathcal {B}} - \sigma _{\rm obs}^2 - \sigma ^2_{\rm off}. \end{equation} \tag{ A8 }$$

The error in δσ²_inst is a propagation of the errors in σ_off (determined by σ^inst_G and σ^inst_T) and $\sigma _{\mathcal {B}}$ because σ_obs is taken to be a fixed quantity. For the high- and low-〈S/N〉 spectra analyzed in Section 5, we find δσ²_inst/(|δσ²_inst|)^1/2 ∼ 5 km s⁻¹.