Parametric Recovery of Line‐of‐Sight Velocity Distributions from Absorption‐Line Spectra of Galaxies via Penalized Likelihood

In an ideal situation in which the LOSVD is perfectly sampled by the data, the (h₃,...,h_M) parameters of the LOSVD (eq. [4]) are essentially uncorrelated to (V,σ). Therefore, one expects the best‐fitting parameters to change little, irrespective of whether (V,σ) are fitted first or together with the other parameters. To the lowest order, the difference between the two approaches is given by eq. (10) of van der Marel & Franx (1993).

An advantage of fitting (V,σ) first is that the Gauss‐Hermite parameters that one obtains coincide with the true Gauss‐Hermite moments integrated over the LOSVD (see van der Marel & Franx 1993 for a detailed discussion). This also means that the value of the parameters (h₃,..., h_M) does not depend on the adopted number M of terms in the expansion. This precise relation between the LOSVD and the extracted kinematical parameters can be very useful for the dynamical modeling (e.g., Rix et al. 1997).

Moreover, when only a single template is adopted (K = 1 in eq. [3]), for each pair of (V,σ) values the best‐fitting value of all the remaining parameters (w₁,h₃,...,h_M,b₀,...,b_L) can be determined linearly. These ideas led van der Marel (1994) to design an efficient method to fit the parameters of the LOSVD in pixel space, where (V,σ) are fitted first and the remaining parameters are linearly expanded at the best‐fitting (V,σ) solution.

For extracting the stellar kinematics of the 72 galaxies of the SAURON survey (de Zeeuw et al. 2002), we needed to operate in pixel space to be able to deal with the contamination due to emission‐line gas present in most of the observed objects (Emsellem et al. 2004). The speed of the van der Marel (1994) pixel‐fitting algorithm was an attractive feature, given the large number (∼200,000) of independent spectra from the survey and the need to compute accurate errors by Monte Carlo simulations. However, the observed LOSVDs are not always well sampled by the SAURON pixel scale of dv ≈ 60 km s ^-1. The Nyquist critical frequency for h₄ is ∼σ/2. This means that undersampling becomes a problem for the derivation of the first four Gauss‐Hermite moments of the LOSVD when the observed velocity dispersion is of the order of ∼2 pixels (for SAURON, 120 km s ^-1).

To test the accuracy of the recovery of the Gauss‐Hermite moments using this method when undersampling becomes significant, we first created a realistic model spectrum by fitting a (logarithmically rebinned) library of Vazdekis (1999) galaxy model spectra to the average SAURON spectrum of the barred lenticular galaxy NGC 3384. This spectrum was oversampled by smoothly interpolating it on a 30 times finer spectral grid and was subsequently convolved with the LOSVD of Figure 1. It was integrated over the SAURON pixels, in the wavelength range 4800–5380 Å, and Poissonian noise was finally added to represent the simulated galaxy spectrum G(x). The template T(x) was obtained by integrating the original oversampled spectrum over the SAURON pixels. In Figure 2 we show an example of such a spectrum, with noise added at S/N = 60.

Zoom In Zoom Out Reset image size

We subsequently tried to recover the true Gauss‐Hermite moments by fitting (V,σ) first and expanding (h₃,...,h_M) linearly at the best‐fitting (V,σ) location. To check our results we computed the true moments (h ˆ₃, h ˆ₄), by analytically integrating over the LOSVD (eq. [8] of van der Marel & Franx 1993) using (V,σ) of the best‐fitting Gaussian computed in § 3. We found h ˆ₃ = -0.144 and h ˆ₄ = 0.029, which, as expected, are very close to the best‐fitting parameters computed before.

The result of this test, for both S/N = 60 and 600, is shown in Figure 3. Here the shape of the input LOSVD was held fixed, but the velocity scale was varied so that the σ_in of the best‐fitting Gauss‐Hermite series varied in the range 48–360 km s ^-1, corresponding to 0.8–6 SAURON pixels, for 1000 different Monte Carlo realizations. By construction, for σ_in = 114.8 km s ^-1 the LOSVD reduces to the one presented in Figure 1. The recovery of the true moments is strongly biased when σ_in≲240 km s ^-1 (4 pixels). This can be understood by the fact that when the LOSVD is not very well sampled, an asymmetry in the profile can be compensated during the Gaussian fit by a small V‐shift, while a symmetric deviation from a Gaussian is nearly equivalent to a change of σ. The bias in the recovery of the input parameters remains unchanged even for a noise‐free model spectrum, although of course the scatter in the values disappears. And no improvement is observed if the template is oversampled before convolving it with a well‐sampled LOSVD in equation (3). These results are representative of what we found with different realistic input LOSVDs.

Zoom In Zoom Out Reset image size

What was seen in the simulation is also found on real SAURON data. The approach of fitting (V,σ) first becomes insensitive to any deviation from a Gaussian LOSVD when the observed velocity dispersion is low (≲120 km s ^-1, or 2 pixels). As a practical example this method gives the misleading impression that the LOSVD in the nucleus of a low‐dispersion object like NGC 3384, containing a fast‐rotating disklike structure, is consistent with being essentially symmetric, as, e.g., the central LOSVD of the high‐dispersion nonrotating giant elliptical galaxy M87. This apparent similarity of the overall shape of LOSVD in the two different galaxies is an artifact of the extraction method, and we will show that it is still possible to recover, with a different method, the strong asymmetry of the LOSVD from the SAURON data of NGC 3384, even at low dispersion.

The bias present in the recovery of the Gauss‐Hermite moments of an undersampled LOSVD in the previous section was mainly due to the fact that (V,σ) were not fitted simultaneously with the other LOSVD parameters. The fit should then become generally unbiased when all (V,σ,h₃,...,h_M) parameters are varied simultaneously to optimize the fit. In this way one should recover the best‐fitting parameters of the Gauss‐Hermite series, which do not precisely coincide with the true Gauss‐Hermite moments of the LOSVD. By definition this method will provide the lowest χ² in the fit, and for this reason it was also originally chosen by van der Marel & Franx (1993).

In Figure 4 we repeated the experiment of the previous section while fitting all the above nonlinear parameters together. As expected the solution is now almost unbiased and the measurements tend to be spread around the true input values for a wider σ_in range. However, when σ_in≲120 km s ^-1 (2 pixels) and the LOSVD is correspondingly not well sampled, the spectrum does not contain enough information to constrain all the free parameters, and the scatter in the recovered values increases dramatically.

Zoom In Zoom Out Reset image size

Moreover, at low σ_in the σ and h₄ values are not unbiased anymore, not being symmetrically distributed around the known input values. The reason for this asymmetry can be understood by looking at the shape of the χ² contours, which measure the agreement between the input spectrum G(x) and its best‐fitting model G_mod(x). In Figure 5 we plot the Δχ² = χ² - χ²_min contours for fits obtained by convolving G(x) with the LOSVD of Figure 1, while keeping V and h₃ fixed to the known best‐fitting values (§ 3) and varying σ and h₄ of G_mod(x). The narrow curved valley of nearly constant χ², for decreasing σ and increasing h₄ from the location of the minimum χ²_min, is due to the strong correlation between these two parameters, which is caused by the undersampling, and produces the effects seen in Figure 4.

Zoom In Zoom Out Reset image size

When the errors in the measured parameters are large, one has to make a trade‐off decision of scatter versus bias. Usually, when the data are unable to tightly constrain a large number of parameters, one prefers to retain in the fit only the parameters that are required to significantly decrease the χ². For this reason, while fitting a parametric LOSVD, it is common practice to reduce the fit to (V,σ) at low S/N: although the fit may be biased by the lack of freedom in the model, an additional flexibility would not lead to meaningful measurements, and structures and trends in the (V,σ) values can still be detected because of the decreased scatter.

Another problem with the large scatter that is present when fitting all parameters together is specific to the two‐dimensional kinematical measurements obtained with an integral‐field spectrograph. When presenting kinematics maps there is no easy way to also attach errors to the displayed values, and it becomes difficult to estimate from the maps when some structure is real and when it is due to noise. For this reason it would be preferable to display on the kinematics maps only the features that are statistically significant.

The LOSVD of galaxies is generally well reproduced by a Gaussian (e.g., Bender, Saglia, & Gerhard 1994). For this reason, it would be preferrable to use a technique to fit the LOSVD in which the solution is free to reproduce the details of the actual profile when the S/N is high, but the solution tends to a Gaussian shape when the S/N is low. This was done for the case in which the LOSVD is described by a nonparametric function by using the maximum penalized likelihood formalism (e.g., Merritt 1997), and we refer to that paper for details. Here we apply the formalism to the case in which the LOSVD is parametrically expanded as a Gauss‐Hermite series. We show that this leads to a much simpler and faster implementation than for the general case.

The idea is to fit the parameters (V,σ,h₃,...,h_M) simultaneously as in § 3.2, but to add an adjustable penalty term to the χ² to bias the solution towards a Gaussian shape when the higher moments are unconstrained by the data so that the penalized χ² becomes

A natural form for the penalty function is given by the integrated square deviation of the line profile (v) from its best‐fitting Gaussian (v):

This penalty does not suppress noisy solutions, which are already excluded by the use of a low‐order parametric expansion for (v). It was shown by van der Marel & Franx (1993) that in the case in which (v) has the form of equation (4), then equation (7) is well approximated by

In principle one could then define the penalty as = ² and optimize χ²_p (eq. [6]). In practice, however, this is not desirable for two reasons:

• 1.  
  As explained in § 2, it is computationally much more efficient to minimize the residuals r_n (eq. [2]) than to explicitly compute the χ².
• 2.  
  One needs a way to automatically adjust the penalty factor α according to the χ² of the observed fit.

We found that a simple and effective solution to these problems consists of using the following perturbed residuals as input to the nonlinear least‐squares optimizer:

where the variance is defined as

The qualitative interpretation of this formula is that a deviation of the LOSVD from a Gaussian shape will be accepted as an improvement of the fit only if it is able to correspondingly decrease the scatter σ(r) by an amount related to . To quantify, one can compute the objective function of the fit

The sum of the residuals in the second term is zero by construction, because of the fact that the weights are optimized for given (v) (§ 2; this is required for this form of perturbation to work). Considering the definition of variance (eq. [10]), one can finally write

which is of the desired form of equation (6) with α = λ²χ² automatically scaled according to the χ² of the fit. In practice σ(r) in equation (9) is computed using a robust biweight estimator (Hoaglin, Mosteller, & Tukey 1983) so that equation (12) is only valid as an approximation. One can see from this formula that with λ = 1, a deviation = 10% of the LOSVD from a Gaussian (e.g., h₃ = 0.1 and h_m = 0) requires a corresponding decrease in the scatter σ(r) of the unperturbed residuals by more than (1 + ²)^1/2 = 0.5% to be accepted by the optimization routine as an improvement of the fit.

A statistical interpretation of equation (12) can be obtained by considering that for a good fit, χ²∼N. For a given λ a deviation from a Gaussian needs to decrease χ² by more than Δχ²∼Nλ²² to be accepted, and this variation can be associated with a specific confidence level (see, e.g., Press et al. 1992, § 15.6) at which a Gaussian shape is excluded. Although this formula may serve as a guideline, we believe it is safer to test a choice of λ using simulations such as the ones presented.

A test of the application of equation (9) to the same model spectra as § 3.2 is shown in Figure 6. For the same S/N = 60 as before, using λ = 0.7, the measurements are now essentially unbiased when σ_in≳120 km s ^-1 (2 pixels). For lower input dispersion the LOSVD converges toward a Gaussian and, as expected, the bias becomes similar to Figure 3. However, the crucial difference with the previous case consists of the fact that for higher S/N = 600 the solution tends to converge to the known best‐fitting values, because the penalty scales with the χ².

Zoom In Zoom Out Reset image size

We mention here that alternative forms to equation (9) are possible, with the requirement that the objective function has the form of equation (6). One possibility is to include a multiplicative perturbation of the residuals:

In this case the objective function becomes

and if one neglects the small term containing ⁴, this becomes again of the form of equation (6), with α = 2λχ². We verified that in practice this alternative form of perturbing the residuals produces virtually the same results as using equation (9) if the same α parameter is adopted. We prefer the perturbation of equation (9), given its robustness against outliers, because of the use of the biweight in the computation of σ(r). The form of equation (13) may be useful in implementations in which the average residuals are not zero for a given (v).

In Figure 7 we compare the three different approaches to the recovery of the LOSVD by using them to extract the stellar kinematics from SAURON (Bacon et al. 2001) integral‐field spectroscopic observations of the galaxy NGC 3384 (see de Zeeuw et al. 2002 for a description of the observations), which has a rather low velocity dispersion over the whole observed field. The data cube was spatially binned to a minimum S/N = 60 using the Voronoi two‐dimensional (2D) binning method of Cappellari & Copin (2003). The same differences that we observed using synthetic‐model spectra can also be seen from real data:

• 1.  
  The Gauss‐Hermite parameters are everywhere significantly suppressed when fitting V and σ first (Fig. 7, top panels).
• 2.  
  The simultaneous fit of all parameters (Fig. 7, middle panels) tends to produce a velocity dispersion map that is noisier and has smaller values, while h₄ tends to be strongly positive. The h₃ values fluctuate in the outer parts between large positive and negative values.
• 3.  
  The use of a penalty function (Fig. 7, bottom panels) overcomes the problems of the two previous approaches. Only the statistically significant non‐Gaussian features are preserved; otherwise, the solution is smoothly reduced to a Gaussian.

We also applied the three methods to the extraction of the kinematics of all the 48 early‐type galaxies of the SAURON survey (de Zeeuw et al. 2002). As expected from the simulations, the three techniques gave very different results only for galaxies like NGC 3384, which have a low‐velocity dispersion. The kinematics extracted from SAURON observations of galaxies with dispersion ≳180 km s ^-1 over the whole observed field provided very similar result with the three methods (although by construction not identical).

Zoom In Zoom Out Reset image size

The availability of fast computers makes Monte Carlo simulations the preferred way to estimate measurement errors for LOSVD extraction methods. This consists of repeating the full measurement process for a large number of different realization of the data, obtained by adding noise to the original spectra.

The method discussed in this paper, like many other available methods for extracting the LOSVD, makes use of some form of filtering or penalization of the solution, thus biasing the results to suppress the noise. It may be worth emphasizing that when computing Monte Carlo errors, one has to correct for the bias in the measurements introduced by the noise‐suppression mechanism. Failure to correct for the bias in the kinematics can result in a significant underestimate of the confidence intervals of parameters determined by fitting dynamical models to the data.

One way to obtain proper errors from penalized methods consists of correcting the error estimates by increasing the percentile intervals by an amount given by the bias (e.g., Beers, Flynn, & Gebhardt 1990). In real measurements one does not know the actual bias (otherwise the measurements would have been corrected for it), so an estimate of the maximum expected bias can be used instead. This requires one to perform Monte Carlo simulations of the kinematic extraction to determine the maximum bias at different S/N levels and σ_in values.

A simpler order‐of‐magnitude estimate of the errors on the parameters may be obtained by noting that Figure 4 provides a more realistic representation of the errors in the measurements than Figure 6, at low σ_in values. In the case of the penalized pixel‐fitting (pPXF) routine described in the previous section, one can then obtain the errors by setting the parameter λ in equation (9) to a very small value before running the Monte Carlo simulations. This simple approach will generally provide a conservative estimate of the actual errors.

To summarize the discussion of the previous sections, the suggested algorithm for recovering the LOSVD parametrized as a Gauss‐Hermite series is as follows:

• 1.  
  Start with an initial guess for the parameters (V,σ), while setting the initial Gauss‐Hermite parameters h₃,...,h_M = 0.
• 2.  
  Solve the subproblem of equation (3) for the weights (w₁,...,w_K,b₀,...,b_L).
• 3.  
  Compute the residuals r_n from the fit using equation (2).
• 4.  
  Perturb the residuals as in equation (9) to get r^'_n.
• 5.  
  Feed the perturbed residuals r^'_n into a nonlinear least‐squares optimization routine and iterate the procedure from step (2) to fit for the parameters (V,σ,h₃,...,h_M).

An IDL routine implementing the algorithm described in this work is available online.²

We thank the SAURON team for making the data available, and in particular Richard McDermid for fruitful discussions and testing of the method, and Tim de Zeeuw, Jesus Falcón‐Barroso, Davor Krajnović, and Glenn van de Ven for commenting on the draft. We are grateful to the referee Roeland van der Marel for useful comments that improved the presentation of this work. M. C. acknowledges support from a VENI grant awarded by the Netherlands Organization for Scientific Research (NWO).