text:
# The counterrotating core and the black hole mass of IC 1459Based on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the AURA, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal no. 7352.

## 1 Introduction

The E3 giant elliptical galaxy IC 1459 ($`M_V22.3`$ \[LEDA, see Paturel et al. 1997\], $`D=30.3\pm 4.0`$ Mpc \[Ferrarese & Merritt 2000\], $`R_e40\stackrel{}{\mathrm{.}}6`$ \[Burkert 1993\]) is a member of a loose group consisting mainly of spirals. It shows clear signs of past interaction, evidenced by stellar “spiral arms” detected on deep photographs (Malin, 1985) and shells at large radii (Forbes & Reitzel, 1995). The nuclear region harbors a fast counterrotating stellar component (Franx & Illingworth, 1988), which cannot be identified as a separate structure in photometric observations. A blue LINER source (unresolved at HST resolution), surrounded by patchy dust absorption and by a gaseous disk, is found at the very center (Forbes, Franx, & Illingworth, 1995; Carollo et al., 1997; Verdoes Kleijn et al., 2000, hereafter VK00).

VK00 measured the mass of the central supermassive black hole in IC 1459 by analyzing kinematical observations of the nuclear gas disk, obtained with the Faint Object Spectrograph (FOS) on board the Hubble Space Telescope (HST). They found a central black hole mass of $`M_{}=(26)\times 10^8M_{}`$ (scaling their values to our adopted distance111The choice of the distance $`D`$ does not influence our conclusions but sets the scale of our models in physical units. Specifically, lengths and masses scale as $`D`$, while mass-to-light ratios scale as $`D^1`$. ). They also calculated isotropic two-integral axisymmetric dynamical models for ground-based major-axis stellar kinematics, resulting in a significantly larger mass estimate ($`M_{}=[46]\times 10^9M_{}`$). For comparison, the $`M_{}`$-$`\sigma `$ relation (Ferrarese & Merritt, 2000; Gebhardt et al., 2000) as given by Tremaine et al. (2002) predicts a black hole mass of $`(1.1\pm 0.3)\times 10^9M_{}`$.

In this paper, we construct axisymmetric three-integral models for IC 1459, using Schwarzschild’s orbit-superposition method (Rix et al., 1997; van der Marel et al., 1998; Cretton et al., 1999). The models are constrained by HST/STIS stellar kinematics measurements along the major axis and deep, ground-based spectroscopic observations along five different position angles of the slit. We show that this data-set allows us to derive information on the nature of the kinematically decoupled component (KDC). We also re-address the value for the mass of the nuclear BH in this galaxy.

This paper is organized as follows. We discuss the spectroscopic and photometric data in Section 2 and present the adopted dynamical model in Section 3. We discuss our results in Section 4 and summarize our conclusions in Section 5.

## 2 Data

### 2.1 HST spectroscopy

An HST spectrum of IC 1459 was obtained with STIS (GO 7352, PI: C.M. Carollo) along the galaxy major axis (as determined from the HST central isophotes), using the G430L grating with the 52″$`\times `$0$`\stackrel{}{\mathrm{.}}`$1 slit. Four exposures, slightly shifted in the spatial direction, were obtained for a total of 9812 seconds. The availability of four shifted images allowed for a good removal of the numerous hot and bad pixels, together with cosmic rays, during the coaddition.

The extraction of the stellar kinematics was performed using the pixel-space Line-Of-Sight Velocity Distribution (LOSVD) fitting method of van der Marel (1994b). In pixel space, the masking of the prominent gas emission lines during the fit becomes straightforward. Specifically, a K3 III template star, observed with the 0$`\stackrel{}{\mathrm{.}}`$2 STIS slit (GO 7388, PI: D. Richstone; no acceptable stellar template was available for the 0$`\stackrel{}{\mathrm{.}}`$1 slit), was convolved with a parametrized LOSVD. The parameters of the LOSVD were optimized to the lowest $`\chi ^2`$ by direct comparison with the observed spectrum. The galaxy spectrum was rebinned spatially before fitting, to an average S/N$``$20 per Å. The error images produced by the STScI pipeline were used to weight the measurements and estimate the errors in the fit. Our IDL222http://www.rsinc.com fitting procedure differs from that of van der Marel (1994b) in the following aspects:

1. We used the elegant DIRECT deterministic global optimization algorithm (Jones, Perttunen, & Stuckman, 1993) to perform the $`\chi ^2`$ optimization for the nonlinear variables $`V`$ and $`\sigma `$. This method is *guaranteed* to find the global minimum of a (possibly non-smooth) nonlinear function inside a finite box, even in the presence of multiple minima. It can be effectively used in many non-smooth real-world optimization problems with a small number of variables ($`N10`$).
2. We performed a $`3\sigma `$ clipped fit to make the fit stable against outliers (e.g. residual bad pixels, unmasked gas emission lines or template mismatch). In practice, once the global minimum is found, the galaxy pixels which deviate more than $`3\sigma `$ from the best fitting template spectrum are added to the list of masked pixels. The global optimization is repeated from scratch on the new set of pixels. This process is iterated until the set of bad pixels does not change anymore.
3. The errors on the $`V`$ and $`\sigma `$ fitting parameters are determined from the confidence levels for a $`\chi ^2`$ distribution with two degrees of freedom.

The resulting dispersion measurements $`\sigma _{\mathrm{fit}}`$ were corrected by approximating the line-spread-function (LSF) by a Gaussian, for both the galaxy and the template star. The template was observed with an effective instrumental dispersion $`\sigma _{\mathrm{LSF},t}100`$ km s-1 (Leitherer et al., 2001), while the galaxy was observed with $`\sigma _{\mathrm{LSF},g}167`$ km s-1. The effective instrumental dispersion for the template star is smaller than that of the galaxy, even though the star was observed with a larger slit. The reason for this is that the STIS PSF falls completely within the 0$`\stackrel{}{\mathrm{.}}`$2 slit that was used to observe the template, so that the template LSF is determined essentially by the PSF size. The corrected dispersion $`\sigma `$ was obtained from

$$\sigma ^2=\sigma _{\mathrm{fit}}^2\sigma _{\mathrm{LSF},g}^2+\sigma _{\mathrm{LSF},t}^2.$$
(1)
This correction is of the same order of magnitude as the measurement errors. No kinematic measurements could be extracted from the two central rows (0$`\stackrel{}{\mathrm{.}}`$05 from the galaxy center) of the STIS spectrum due to the sharp decrease of the absorption line S/N of the stellar component, caused by the rise of a featureless non-thermal continuum from the central source. The non-thermal nature of the continuum is indicated by the sharp decrease of the line-strength parameter $`\gamma `$ towards the center.

The measured kinematical profiles are presented in Table 1 and in Figure 8.

### 2.2 Ground-based spectroscopy

We also analyzed a set of ground-based spectroscopic observations taken with the CTIO 4m telescope in November 1988. These observations were made along five different position angles with a 1$`\stackrel{}{\mathrm{.}}`$5 wide slit, using a CCD with a 0$`\stackrel{}{\mathrm{.}}`$73 pixel size. The seeing PSF was $`1\stackrel{}{\mathrm{.}}5`$ FWHM during the observations. The major axis, two $`\pm 45^{}`$ intermediate axes and two position angles close to the minor axis were observed. The location of the different slit positions is shown in Figure 1. The major axis data were presented previously in van der Marel & Franx (1993). The mean velocity $`V`$, the velocity dispersion $`\sigma `$ and the higher order Gauss-Hermite moments ($`h_3`$$`h_6`$) were extracted from the high signal-to-noise spectra by the method described in that paper. The data are presented in Table 2.

### 2.3 HST imaging: correction of dust absorption

For our dynamical modeling, we need to derive an $`I`$-band surface density model from the available WFPC2/F814W (460 s) Hubble Space Telescope image of IC 1459. A complication is that this galaxy contains gas and dust in its central regions, which has to be corrected for in order to obtain an accurate determination of the central stellar surface brightness. We derive this correction in a manner that is similar to Carollo et al. (1997, hereafter C97) and assume

1. the intrinsic color of the galaxy pixels that are affected by dust is the same as that of the surrounding pixels;
2. the dust is a screen in front of the stellar emission;
3. the Galactic extinction law (Cardelli, Clayton, & Mathis, 1989) applies.

In this process, we use the available WFPC2/F555W (1000 s) frame in addition to the WFPC2/F814W image.

In practice, C97 corrected the pixel values by assuming that all pixels on the same isophote have the same intrinsic color. We add the assumption that the intrinsic galaxy color varies smoothly as a function of radius. Specifically here we assume the color varies linearly with the logarithm of the radius. This assumption is justified by Figure 2. The pixel correction then involves the following steps:

1. We measure the average ellipticity and position angle on the original galaxy image.
2. We construct a calibrated color map (in our case $`VI`$).
3. We perform a straight line fit to the pixel colors as a function of $`\mathrm{log}m`$, where $`m`$ represents the semimajor axis length (see Figure 2). The fit is done by minimizing the absolute deviation to reduce the effect of large systematic pixel value deviations due to dust.
4. We compute an $`E(VI)`$ map by calculating, for every pixel, the color predicted from the previously fitted relation and subsequently subtracting this from the pixel value. The result is shown in Figure 3.
5. We correct the pixels above a given $`E(VI)`$ threshold (in general as a function of $`m`$), using the standard Galactic extinction curve.

This ad hoc method corrects the major effects of patchy dust absorption and can also be used to measure dust corrected color gradients (compare Figure 2 with Figure 1o in C97).

## 3 Dynamical model

### 3.1 The mass model with MGE

IC 1459 is not a perfectly axisymmetric body, evidenced by a small $`5^{}`$ isophotal major axis twist in the range 25″–100″ (Franx et al. 1989; C97). Nevertheless, assuming axisymmetry seems reasonable to study its nuclear dynamics, given that our kinematical observations are all located inside the region that can be well reproduced by assuming axisymmetry.

The first step in the dynamical modeling is to obtain a parametrization of the stellar surface density. For this we adopted the Multi-Gaussian Expansion (MGE) method (Monnet, Bacon, & Emsellem, 1992; Emsellem, Monnet, & Bacon, 1994), which allows for non-elliptical isophotes and radial ellipticity variations.

We obtained an MGE fit to the WFPC2/F814W dust-corrected images of IC 1459 using the method and software developed by Cappellari (2002). The MGE fit was performed using both the PC1 (at full resolution) and the entire WFPC2 mosaic simultaneously. The PSF was parametrized by fitting a circular MGE model, of the form $`PSF(R^{})=_{k=1}^MG_k\mathrm{exp}[R^2/(2\sigma _k^2)]/(2\pi \sigma _k^2)`$, to a PSF computed with TinyTim (Krist & Hook, 2001) for the center of the PC1 CCD. The numerical values of the relative weights $`G_k`$ (normalized such that $`_{k=1}^MG_k=1`$), and of the dispersions $`\sigma _k^{}`$ are given in Table 3.

Figure 4 shows a comparison between the observed photometry and the MGE model along a number of profiles, while Table 4 gives the corresponding numerical values of the analytically deconvolved MGE parametrization of the galaxy surface brightness

$$\mathrm{\Sigma }(x^{},y^{})=\underset{j=1}{\overset{N}{}}\frac{L_j}{2\pi \sigma _j^2q_j^{}}\mathrm{exp}\left[\frac{1}{2\sigma _j^2}\left(x^2+\frac{y^2}{q_j^2}\right)\right].$$
(2)
Also tabulated are the distance-independent central surface brightness of each Gaussian $`I_j^{}=L_j/(2\pi \sigma _j^2q_j^{})`$. The central Gaussian (no. 1) represents the contribution of the blue unresolved central spike, which is of non-thermal origin and will not be included in the dynamical model calculation. The comparison between the isophotes of the convolved model and the actual image is shown in Figure 5. We confirm the finding of C97 that there is no photometric evidence for a nuclear disk in the dust corrected images.

### 3.2 Three-integral models

We model IC 1459 by using an axisymmetric three-integral implementation of Schwarzschild’s (1979) orbit superposition method, developed by Rix et al. (1997); van der Marel et al. (1998); Cretton et al. (1999) and Richstone et al. (2002). The method is described extensively in the above references and we will not repeat the details here. Very briefly, the stellar density is derived by deprojecting (see Appendix A) an MGE parametrization of the observed surface density (equation ), assuming axisymmetry and a value for the stellar mass-to-light ratio $`\mathrm{{\rm Y}}`$. Stellar orbits are integrated in the potential that is generated by both the stellar density and a central supermassive BH. Each orbit is projected onto the space of the observables, taking into account PSF convolution and aperture binning. Finally, the distribution of orbital weights that best fits both the luminosity density distribution of the model and the observed kinematics is determined by solving a Non-Negative Least-Squares (NNLS, Lawson & Hanson 1974) problem.

We modified the software implementation by van der Marel et al. (1998, hereafter M98) such that it can deal with the MGE parametrization of the surface brightness derived in the previous paragraph, in a way that is similar to Cretton & van den Bosch (1999). The differences between our software and the implementation used by M98 are:

1. Simplified calculation of the potential. This is possible because of the MGE formalism and the details described in Appendix A;
2. Substitution of the numerical integration of the luminosity density along the line of sight with the analytic MGE projection;
3. Simplified and more accurate (due to the smaller number of numerical integrations involved) evaluation of the density integration on a polar or Cartesian grid (Appendix B).

### 3.3 Orbit library

To include a representative set of orbits in the orbit library, we sample them on a grid that covers the full extent of the integral space $`(E,L_z,I_3)`$, where $`E`$ is the energy, $`L_z`$ is the angular momentum along the $`z`$ symmetry axis and $`I_3`$ is a non-classical integral. This is achieved by sampling the energy logarithmically in the radius of the circular orbit with the given energy $`E(R_c)`$, from 0$`\stackrel{}{\mathrm{.}}`$05 to 300″. This grid includes 99% of the total luminous mass of the galaxy. The quantity $`\eta L_z/L_{\mathrm{max}}`$ is sampled linearly on an open grid in the interval $`[0,1]`$. $`I_3`$ is sampled linearly in the angle $`w`$ on an open grid in the interval $`[0,w_{\mathrm{th}}]`$, where $`w_{\mathrm{th}}`$ is the angle $`w`$ for the “thin tube” orbit at the given $`(E,L_z)`$ (see Figure 5 in M98 for details). After some testing we chose a $`20\times 14\times 7`$ grid in integral space (including both positive and negative $`L_z`$), but similar results (with longer computation times) were obtained with a $`25\times 26\times 13`$ grid. This choice of parameters is similar to that used in the modeling of other nearby galaxies (e.g. van der Marel et al. 1998; Verolme et al. 2002).

## 4 Results

### 4.1 Stellar kinematics

Our model has three free parameters: the galaxy inclination $`i`$, the BH mass $`M_{}`$ and the stellar mass-to-light ratio $`\mathrm{{\rm Y}}`$. Inclinations smaller than $`i55^{}`$ can be excluded based on the observed axial ratio of IC 1459, since these values correspond to intrinsic stellar densities that are too flat for an elliptical galaxy.

We fitted dynamical models to the $`V`$, $`\sigma `$ and $`h_3`$$`h_6`$ observations along all five position angles and to the STIS/G430L $`V`$ and $`\sigma `$ measurements. We varied the parameters $`\mathrm{{\rm Y}}`$ and $`M_{}`$ for three different values of the inclination $`i=60^{}`$, $`i=75^{}`$ and $`i=90^{}`$. The overall best fit was found for an inclination $`i=90^{}`$ (corresponding to edge-on viewing). Contours of equal $`\chi ^2`$, using the confidence levels for a $`\chi ^2`$-distribution with three degrees of freedom, are shown in the upper panel of Figure 6. The best-fit parameters are a mass-to-light ratio of $`\mathrm{{\rm Y}}=3.1\pm 0.4`$ (in the $`I`$-band) and a central black hole mass of $`M_{}=(2.6\pm 1.1)\times 10^9M_{}`$ (3$`\sigma `$ level). The inclination is not strongly constrained by our data and our results do not depend significantly on the adopted value. In particular, the best-fitting black hole mass at $`i=60^{}`$ is $`M_{}=(2.5\pm 1.2)\times 10^9M_{}`$. The best-fit model has $`\chi ^2867`$ with $`N=594`$ kinematical constraints. The fact that $`\chi ^2>N`$ even for the best model, is due to some small systematic errors, which are in particular visible for the even Gauss-Hermite moments $`h_4`$ and $`h_6`$. The excellent match (Figure 7, 8) to the observed kinematical constraints along all position angles suggests that the decoupled core in IC 1459 can be well described as an axisymmetric component. This allows us to study its internal dynamics, shape and mass distribution.

The solid lines in Figure 7 and 8 were calculated by fitting to only the photometry and kinematics. The kinematical constraints are distributed over 95 apertures (not all independent), each of which provides 6 Gauss-Hermite moments and 12 STIS apertures with $`V`$ and $`\sigma `$ measurements. Selfconsistency is ensured by fitting to the mass in the meridional (120 constraints) and projected plane (140 constraints), as well as to the mass fractions contained in the kinematical apertures. The models therefore have to satisfy a total of 961 constraints, which are to some extent correlated. Since our library contains $`20\times 14\times 7=1960`$ orbits whose orbital weights are the unknowns that have to be determined, there is no unique solution that defines the internal dynamics of this galaxy (in the approximation of this model). In fact, the NNLS matrix equation that is solved is numerically rather ill-conditioned, giving rise to a distribution function composed of sharp isolated peaks. Such DFs are unphysical. A smoother and better (see Cretton et al. 1999; Verolme & de Zeeuw 2002) solution can be obtained by using regularization. The regularization scheme that is used here minimizes, up to a certain degree, the differences in weights between neighboring orbits. Although these additional regularization constraints result in a slightly poorer fit to the data, a considerable amount of smoothing can be forced upon the solution before systematic deviations in the fit become apparent.

Figure 9 shows the integral space \[the weights on the $`(E,L_z,I_3)`$ grid\] of our best-fitting model, for a selected set of energies (corresponding to the radii effectively constrained by the observations). The top row shows the integral space for the unsmoothed solution which, as a consequence, appears very noisy. A small amount of regularization ($`\mathrm{\Delta }=4`$; see M98 for details) was added as a constraint in the middle row. While the fit to the data remains essentially unchanged, some coherent structures clearly emerge in the solution space. The bulk of the galaxy rotation can be recognized as a large bright (high weight) region at low and negative $`\eta L_z/L_{\mathrm{max}}`$ and high $`w`$ on the range $``$1″–13″. At radii smaller than $``$0$`\stackrel{}{\mathrm{.}}`$8, the regularized solution contains a well-separated group of nearly circular orbits ($`\eta 1`$), which counterrotates with respect to the main galaxy body. These orbits appear necessary to explain the observed rapid counterrotation in the kinematics. The contribution of the counterrotating component seems to decrease inside $`R0\stackrel{}{\mathrm{.}}2`$. This is probably not a real effect: there is not enough data inside this radius to constrain the phase space in detail and its appearance is dictated by the smoothness constraint.

The structure of integral space remains qualitatively the same at high regularization ($`\mathrm{\Delta }=1`$, bottom row). This latter model is still able to produce a reasonable fit to all constraints (see dashed line in Figure 7, 8). Although our conclusions do not depend strongly on the regularization parameter $`\mathrm{\Delta }`$, in the following we adopt $`\mathrm{\Delta }=4`$ for our plots. Both Cretton et al. (1999) and Verolme et al. (2002) have carried out a number of tests of the Schwazrschild code, and established which value of $`\mathrm{\Delta }`$ gives the best reconstruction of a known distribution function. The result is a $`\mathrm{\Delta }4`$. Detailed comparison with independent results obtained by Gebhardt with a maximum entropy technique (which starts at the opposite end with a fully smooth distribution function, and then makes it less and less so) shows that again an amount of smoothing equivalent to our $`\mathrm{\Delta }4`$ is optimal.

Both the best fit values and the details of the internal kinematic structure are in excellent agreement with our previous results (Cappellari, 2000; Cappellari et al., 2001), which started from a triple power-law parametrization of the stellar surface density instead of the MGE expansion and were constrained by ground-based data only. This comparison shows that the best-fitting parameters do not depend on the details of the adopted stellar density parametrization. However, our best-fitting black hole mass is slightly larger than the value $`(1.1\pm 0.3)\times 10^9M_{}`$ that is predicted by the $`M_{}`$-$`\sigma `$ relation (Ferrarese & Merritt, 2000; Gebhardt et al., 2000) as given by Tremaine et al. (2002).

The contribution of a dark halo becomes evident in the observed stellar kinematics only well beyond an effective radius, and hence can be safely ignored (Carollo et al. 1995; Gerhard et al. 2001). Furthermore, we follow standard practice and assume that the stellar mass-to-light ratio does not vary with radius. To determine whether our results depend on this assumption, we also ran models in which $`\mathrm{{\rm Y}}`$ follows the change of the $`VI`$ color presented in Fig. 2. The color gradient of IC 1459 is $`d(VI)/d(\mathrm{log}R)0.084`$. Using models from Charlot, Worthey & Bressan (1996), the color variation can be translated into a fractional decrease of $`\mathrm{{\rm Y}}`$ by $`15`$ per cent per decade in radius, more or less independent of whether the color gradient is due to age or to metallicity variations. The results obtained from models with a non-constant $`\mathrm{{\rm Y}}`$ do not show any significant difference in the best fitting BH mass (bottom panel of Figure 6) or in the internal dynamical structure, compared to the models with constant $`\mathrm{{\rm Y}}`$.

Given that the counterrotating component appears so well isolated from the rest of the orbital distribution (see also Figure 10), we can make a quantitative estimate of the mass involved. This is not possible by any other means, due to the fact that the counterrotating component is not apparent from the photometry and does not appear as a double component in the LOSVD. By summing the solution weights of the orbits that contribute to the counterrotation, we find that this component carries $`0.5`$% of the luminous galaxy mass which is $`3\times 10^9M_{}`$ (determined from the Gaussian components of Table 4, using the best fitting mass-to-light ratio).

Figure 11 shows the $`V/\sigma `$ field in the meridional plane of the model. This was calculated by adding the velocity moments of the orbits in the best fitting solution with low regularization of Figure 9. The flattened appearance of the counterrotating stellar component is clearly visible. Figure 12 presents the velocity moments in polar coordinates. We see the nuclear kinematics shows mild tangential anisotropy, which compares well with the corresponding measurements obtained for other galaxies (Gebhardt et al. 2002). However, the underlying structure is more complex as the distribution function is in fact the sum of the bulk of the galaxy and the counter rotating component, which is not evident in the second moments.

Cretton & van den Bosch (1999) showed that the Gauss-Hermite (GH) parametrization for the LOSVD, which is used to constrain the kinematics, can produce artificial counterrotations in some special cases. This is essentially due to the fact that the model is trying to reproduce a small number of GH moments and not the actual observed LOSVD and when two functions coincide on a small set of GH moments this does not guarantee that they are the same function. This phenomenon can be easily recognized after the fit, however, by comparing the LOSVDs predicted by the model and the observed ones. Before interpreting the internal dynamics of IC 1459, we check that our best fit dynamical model does not suffer from this shortcoming.

Figure 13 shows the LOSVD of the model with low regularization ($`\mathrm{\Delta }=4`$), obtained by direct summation on the velocity histograms of the best fitting orbits, together with the observed LOSVD. The latter was reconstructed from the Gauss-Hermite parametrization introduced by van der Marel & Franx (1993), using the mean velocity, the velocity dispersion and the Gauss-Hermite moments that are plotted in Figure 7 and 8. The Hermite polynomials are tabulated in the Appendix C. The plot only represents the most extreme observed LOSVD, but similar results are obtained at the other ground-based apertures. It is apparent from this comparison that the model provides an excellent fit to the observed LOSVD, at least well within the errors suggested by the observed differences in the LOSVD at opposite sides of the galaxy, that are visible in Figure 13. We can thus safely affirm that we are *not* in the situation of Figure 10 of Cretton & van den Bosch (1999) and that our model reproduces the data correctly.

### 4.2 The STIS gas kinematics

The very central STIS/G430L spectrum of IC 1459 covers the wavelength range between 3500 and 5800 Å. It is dominated by an unresolved source of featureless, non-stellar continuum and by prominent and broad ($`\sigma 500`$ km s-1) forbidden gas emission lines of \[O II$`\lambda `$3727, \[Ne III$`\lambda `$3869, \[S II$`\lambda \lambda `$4068,4076, \[O III$`\lambda `$4363, \[O III$`\lambda \lambda `$4959,5007, \[N I$`\lambda `$5200 and the Balmer series emission-lines H$`\delta `$, H$`\gamma `$ and H$`\beta `$. A detailed analysis of the physical properties of this spectrum will be presented elsewhere. In this Section we focus on the gas kinematics.

Of all emission lines, \[O II$`\lambda `$3727 can be traced reliably to the largest distance from the center ($``$1″). Single Gaussians provide a good fit to the emission lines and are therefore used to determine the mean velocity and velocity dispersion at each row of the spectrum (Figure 14, Table 5). The reality of the velocity curve that was determined in this manner can be verified in Figure 15. The line intensity was normalized by dividing each column of the spectrum by the maximum intensity of the \[O II\] line in that column. The mean velocity curves for the other species agree within the errors with the \[O II\] line in the regions of overlap. The galaxy nucleus is assumed to be at the peak in the continuum flux. The most salient features are a central steep velocity gradient in the direction opposite to the nuclear stellar rotation and a clear asymmetry of the velocity curve with respect to the nucleus. The velocity curve reaches a peak velocity difference of $`\pm 200`$ km s-1, then drops rapidly to zero and even changes sign at $`0\stackrel{}{\mathrm{.}}5`$ on one side of the nucleus.

The asymmetry in the rotation curve with respect to the nucleus prevents a good fit to the complete curve using an axisymmetric gas dynamical model. However, the velocity curve inside $`R0\stackrel{}{\mathrm{.}}3`$ shows a clear sign of rotation and one might follow standard practice of assuming a gas disk. VK00 obtained 6 FOS emission-line spectra in this region which together with ground-based CTIO emission-line spectra suggested a thin disk in Keplerian rotation. Under this assumption, they determined a BH mass of $`M_{}=1.5\times 10^8M_{}`$ from a combined best-fit to the H$`\beta `$ and H$`\alpha `$+\[N II\] gas kinematics (scaling their values to our adopted distance). The rotation velocities that VK00 obtained for H$`\alpha `$+\[N II\] with the smallest FOS aperture are compared to the STIS results in Figure 14. The data are in broad agreement, but the STIS data clearly provide a more complete and detailed view of the gas kinematics within the central arcsecond. In view of this we have performed a new analysis of the gas kinematics.

We constructed a model for the STIS kinematics using the IDL software we developed in Bertola et al. (1998). The model assumes the gas moves on circular orbits in a thin disk in the galaxy symmetry plane, in the combined potential of the stars and a possible central BH. Our model can deal with a general gas surface brightness distribution (e.g. an observed narrow band image) and takes into account pixel binning, PSF and slit effects to generate a two-dimensional model spectrum with the same pixel scale as the observations. Similar to the observations, the mean model velocity was determined by fitting a Gaussian to each row of the model spectrum. As the STIS observations use the narrow 0$`\stackrel{}{\mathrm{.}}`$1 slit, whose width is comparable to the PSF FWHM, we neglect the velocity offset correction (van der Marel, de Zeeuw, & Rix, 1997; Maciejewski & Binney, 2001; Barth et al., 2001). A detailed comparison between our gas modeling software and the completely independent one used by VK00 shows agreement to better than 5 km s-1 everywhere, for both $`V`$ and $`\sigma `$.

For our STIS model we used the $`i=60^{}`$ inclination, the parametric approximation of the disk surface brightness and the ‘turbulent’ velocity dispersion parameters as derived by VK00. The stellar potential was determined by deprojecting the axisymmetric MGE fit to the galaxy surface brightness that was used for the stellar kinematic modeling and multiplying this with the best fitting stellar mass-to-light ratio $`\mathrm{{\rm Y}}`$. This $`\mathrm{{\rm Y}}`$ equals the one used by VK00 (taking into account the differences in assumed distance).

In Figure 14, we compare the gas velocities with models with several different BH masses. The main result is that none of the models properly fits the gas kinematics. This is not entirely surprising, given the somewhat disturbed nature of the gas kinematics, as discussed above. Nonetheless, it is striking to note that the model with $`M_{}=2.6\times 10^9M_{}`$, the best-fit value inferred from the stellar kinematics, does not even come close to fitting the gas kinematics. It predicts a central rotation curve gradient that is steeper than what is observed and a a rotation curve amplitude that is larger than what is observed. The BH mass must be lowered by a factor of $`7`$ to provide a reasonable fit in the central region $`R0\stackrel{}{\mathrm{.}}3`$. The combined STIS, FOS and CTIO data in this region are best fit by a model with $`M_{}3.5\times 10^8M_{}`$. This is somewhat larger than the value inferred by VK00, but is otherwise in reasonable agreement. A model with no BH predicts a rotation curve gradient that is considerably more shallow than what is observed.

### 4.3 The BH mass discrepancy

Discussions of BH demography have recently focused on the scatter in the relation between BH mass and velocity dispersion. Stellar and gas kinematical BH mass determinations are consistent with the same relation and suggest a scatter of $`0.3`$ dex (see e.g. the compilation by Tremaine et al., 2002). It is therefore surprising that our analyses of the stellar and gaseous kinematics in one and the same galaxy have yielded BH masses that differ by a factor $`7`$, i.e., $`0.9`$ dex. To understand the origin of this discrepancy we proceed by discussing the potential systematic problems that may have plagued our analysis.

#### 4.3.1 Possible problems with the interpretation of the gas kinematics

There are several possible reasons why our analysis of the gas kinematics may have yielded an incorrect estimate of the BH mass.

First, it is possible that the underlying assumptions of our model are correct, but that we have used an incorrect value for the inclination of the inner gas disk ($`i=60^{}`$) or the angle between its projected major axis and the slit ($`\theta =0^{}`$). These angles are derived from images of the emission-line gas disk at arcsec scales and larger (see VK00 for discussion). However, the ellipticity of the gas disk isophotes in the range 0$`\stackrel{}{\mathrm{.}}`$25–1$`\stackrel{}{\mathrm{.}}`$0 is $`ϵ=`$0.17–0.37. VK00 interpreted this as thickening of the gas disk. Alternatively, by allowing the disk to be actually more face-on, the ellipticity is consistent with inclinations between $`10^{}`$ and $`22^{}`$. By constructing models for various $`(i,\theta )`$, we find that a model with $`i20^{}`$, $`\theta 0^{}`$ and $`M_{}=2.6\times 10^9M_{}`$ indeed provides a reasonable and much better fit to the STIS mean gas velocities at $`R0.3^{\prime \prime }`$ than a model with $`M_{}=3.5\times 10^8M_{}`$ (Figure 14). However, this interpretation implies that the gas disk is warped (e.g. Statler 2001) and that it cannot reside in the equatorial plane of the galaxy for $`R1^{\prime \prime }`$ (the galaxy inclination cannot be lower than $`i55^{}`$: see Section 4.1). Although the spherical BH potential dominates the galaxy potential for $`R0\stackrel{}{\mathrm{.}}65`$, the severe warping likely invalidates the assumption of a thin circular gas disk in equilibrium.

Second, it is possible that the underlying assumptions of our model are incorrect, and that the gas does not move on circular orbits in an infinitely thin disk. In our models we have assumed that the observed velocity dispersion of the gas is ‘turbulent’ and does not contribute to its hydrostatic support. However, it may be that the gas is better modeled as individual clouds that move ballistically. In this case the velocity dispersion would contribute to the hydrostatic support and the gas would rotate slower than the circular velocity (asymmetric drift). This would have caused our models to underpredict the BH mass. If $`\sigma /V`$ is small, then the upward correction to the BH mass is small and it can be estimated using approximate equations for asymmetric drift (e.g. Barth et al., 2001). However, for IC 1459 the value of $`\sigma /V`$ of order unity, and to properly estimate the BH mass one would have to construct axisymmetric collisionless models for the gas kinematics. This is outside the scope of the present paper. However, one limiting case that can be calculated fairly easily is to assume that the gas can be modeled through the Jeans equations as an isotropic spherical distribution of collisionless cloudlets. This may not be entirely unreasonable, given that the isophotes of the gas disk become rounder towards the center. We find that the BH mass must then be $`M_{}1.0\times 10^9M_{}`$ to reproduce the observed central velocity dispersion of the gas (this supersedes the somewhat smaller value that we quoted for this scenario in VK00).

Third, it is possible that our models for the gas kinematics of IC 1459 are flawed at an even more basic level. The gas may not be in equilibrium. The fact that IC 1459 has a counter-rotating core, possibly indicative of a recent accretion event, may be of relevance in this context. Or alternatively, the gas kinematics (mean velocities and velocity dispersions) may not be dominated by gravitational forces. For galaxies in which the central emission-line gas displays a LINER-type spectrum the gas excitation mechanism might be shock-ionization (e.g. Dopita et al., 1997) as opposed to photo-ionization (as inferred for Seyferts and Quasars). IC 1459 indeed shows LINER-type emission and hence shocks may influence the gas kinematics.

IC 1459 was originally identified by us as a good candidate for thin-disk modeling because of the regular, smooth and prominent rotation evident in ground-based CTIO spectra. With FOS, spectra could only be obtained for a few individual apertures (VK00). This yielded insufficient information to test the hypothesis that the gas is indeed rotating in a thin disk. The new STIS spectra provide a more complete view of the gas disk kinematics and clearly show that thin disk models are oversimplified. The gas does not rotate at all at $`R1^{\prime \prime }`$, indicating that the true state of the gas must be considerably more complicated than what we have been able to represent in any of our models. In hindsight, IC 1459 is a poor galaxy for which to attempt a gas kinematical determination of the central BH mass. The values that we have inferred range from $`3.5\times 10^8M_{}`$ (thin gas disk) to $`1.0\times 10^9M_{}`$ (spherical distribution of cloudlets), but we cannot attach strong confidence to either of these numbers.

#### 4.3.2 Possible problems with the interpretation of the stellar kinematics

The highest spatial resolution stellar kinematical data presented here are based on a STIS spectrum obtained in 4 orbits of HST time. While this is a significant exposure time, the S/N of the spectrum is relatively low (due mostly to the very narrow slit employed). This has several consequences. The low S/N, the relatively low spectral resolution, and the presence of prominent gas emission lines filling important absorption features, prevent the reliable extraction of higher-order Gauss–Hermite moments. As a result, we can only obtain the mean velocity $`V`$ and velocity dispersion $`\sigma `$ of the best-fit Gaussian, with error bars that are considerably larger than those for the ground-based kinematics. Furthermore, the unresolved source of featureless non-thermal continuum that dominates the central spectrum of IC 1459 prevents the measurement of the stellar kinematics inside 0$`\stackrel{}{\mathrm{.}}`$1. The latter is particularly unfortunate, because the influence of a BH is always expected to be largest close to the galaxy center.

The $`V`$ and $`\sigma `$ profiles derived from the STIS spectrum (for $`R`$0$`\stackrel{}{\mathrm{.}}`$1) do not show any obvious evidence for the presence of a BH (Figure 8). While the central velocity gradient is steeper than the one measured from the ground, there is no strong central increase in $`\sigma `$. Also, the value of $`\sigma `$ at 0$`\stackrel{}{\mathrm{.}}`$1 is no larger than the central dispersion observed from the ground. This might suggest a low black hole mass, and indeed, a dynamical model with the value, derived from the gas kinematics, of $`M_{}=3.5\times 10^8M_{}`$, is consistent with the STIS data. However, our model with the much larger mass of $`M_{}=2.6\times 10^9M_{}`$ (inferred from our analysis of all the stellar kinematical data) is also consistent. The reason for this surprising fact can be understood by considering the LOSVD predicted by the model with $`M_{}=2.6\times 10^9M_{}`$ for the STIS apertures closest to the center ($`R=\pm 0\stackrel{}{\mathrm{.}}10`$; Fig. 16). It is strongly non-Gaussian and has very broad wings, caused by the high-velocity stars near the central BH. Qualitatively similar profiles, close to the BH, were predicted with analytic models (e.g. van der Marel, 1994a) and have been observed with HST since (e.g. Joseph et al., 2001). Figure 17 shows that the true centered moment $`\sigma _{\mathrm{true}}`$ displays a large increase inside the BH sphere of influence ($`R_{\mathrm{BH}}0\stackrel{}{\mathrm{.}}65`$). This rise is absent in $`\sigma `$ because the best-fitting Gaussian is quite insensitive to the wings, and as a result the model can fit the flat $`\sigma `$ profile even with a large black hole mass.

As a further check of the consistency of our $`M_{}=2.6\times 10^9M_{}`$ dynamical model with the STIS observations, we convolved a well-sampled K3 III stellar template with the predicted model LOSVD at the different apertures along the STIS slit and rebinned this to the STIS pixel scale. The spectrum thus obtained was processed through the same procedure, described in Section 2.1, to extract the $`V`$ and $`\sigma `$. The resulting values agree to within the errors with the model profiles presented in Figure 8.

For the above reasons, our black hole mass estimate is determined primarily by the ground-based data. Our axisymmetric models take into account the limited spatial resolution of the data, and they formally rule out the BH mass $`M_{}=3.5\times 10^8M_{}`$ that was suggested by the gas kinematics. This comes about because the central $`\sigma `$ measured from the ground is larger than the models can handle with a black hole mass of less than $`10^9M_{}`$, while also fitting the higher order moments. However, Fig. 18 shows just how small this effect is. The effect is formally quite significant, but this assumes that there are no uncertainties involved in the data-model comparison other than Gaussian random errors. This is obviously an oversimplification. On the one hand, the underlying assumptions of the models may not be fully correct. For example, IC 1459 could be significantly triaxial, the stars may not yet have settled into an equilibrium configuration, or line-of-sight projections may be partially compromised by dust. On the other hand, there could be systemic problems with the data, e.g., due to continuum subtraction or template mismatching. Such effects can easily compromise the data at the few percent level (van der Marel et al., 1994c). Since the very central region is affected by a continuum component that is not present at larger radii, it is quite possible that these errors could depend on radius. While there are no obvious indications that any of these issues may be affecting our analysis, it does mean that our BH mass determination via the stellar kinematics should be treated with caution.

The accuracy of stellar dynamical BH mass determinations hinges critically on the availability of full LOSVD measurements inside the BH sphere of influence. Unfortunately, we have not been able to secure such a measurements with HST due to the complications introduced by the modest $`S/N`$, the presence of a non-thermal nuclear point source and of strong gas emission lines. Our BH mass measurement $`M_{}=(2.6\pm 1.1)\times 10^9M_{}`$ therefore rests primarily on $`1.5^{\prime \prime }`$ resolution ground-based spectroscopy. This makes the accuracy of our measurement rather modest, despite the use of three-integral models and the well-defined formal errorbar.

## 5 Discussion and conclusions

We studied the dynamics of the elliptical galaxy IC 1459 using self-consistent axisymmetric dynamical models based on Schwarzschild’s orbit superposition method. The intrinsic density is obtained by deprojecting the observed HST/WFPC2 surface density using an MGE parametrization, after correction for the effects of patchy dust absorption. The axisymmetric model provides an excellent fit to all HST and ground-based kinematical constraints along five slit positions, while also matching the observed HST photometry. Moreover we verified that the model velocity histograms indeed reproduce the observed LOSVDs and not only its moments.

The good match to the observations obtained with an axisymmetric model provides a test that the counterrotating core in IC 1459 can be accurately described as an axisymmetric stellar component. We explored the $`\chi ^2`$ map as a function of the three free parameters of the problem $`(i,\mathrm{{\rm Y}},M_{})`$ and subsequently analyzed the best fit solution by adding linear regularization constraints in the integral space. The observations can only be reproduced if $`0.5`$% of the galaxy stellar mass ($`3\times 10^9M_{}`$) counterrotates in the radial range $`R\mathrm{\hspace{0.17em}2}\mathrm{}`$ on nearly circular orbits in a flattened disk. This counterrotating component is found to be very well separated from the bulk of the stars in integral space. The mass estimation of the counterrotating stellar component is only possible by using dynamical modeling, since the decoupled component does not show in the surface brightness or color maps and cannot be isolated from the observed line profiles.

These results could be explained by a scenario in which the stellar counterrotating core in IC 1459 has an external origin, as evidenced by its clear isolation in integral space from the bulk of the stars in the galaxy. For example, the core may have formed by acquisition of a gas-rich component that settled into a disk in the galaxy potential and later turned into stars. In this scenario it is likely that the acquisition must have happened long ago. The current stellar core counterrotates not only with respect to the outer stellar body but also with respect to the nuclear gas-disk of IC 1459, which must have formed later. Numerical simulations of gas accretion by a BH that leads to disk formation and BH fueling were made by e.g. Bekki (2000) However, other formation scenarios are also possible. For example, Holley-Bockelmann & Richstone (2000) showed that, in the presence of a central BH, a rapidly rotating counterrotating core can also be produced by acquisition of a dense stellar satellite.

IC 1459 has a rotating gas component near its center. This offers the opportunity to estimate the BH mass independently from either gaseous or stellar kinematics. For both we have available high spatial resolution long-slit kinematics from STIS, as well as ground-based kinematics along multiple position angles. Unfortunately, we have found that neither the gaseous nor the stellar kinematics yield a very stringent determination of the BH mass. The main problem for the gas kinematical modeling is that there is evidence that the gas motions are disturbed, possibly due to non-gravitational forces acting on the gas. The main problem for the stellar kinematical modeling is that we were unable to obtain an HST measurement of the stellar LOSVD inside the BH sphere of influence, due to the modest $`S/N`$ of our STIS spectrum and due to the presence of the non-thermal nuclear point source of IC 1459 and prominent gas emission lines. Consequently, the stellar dynamically inferred BH mass hinges almost exclusively on the low spatial-resolution ground-based data. These complications may explain why we find rather discrepant BH masses with the different methods. The gas kinematics suggests that $`M_{}=3.5\times 10^8M_{}`$ if the gas is assumed to rotate at the circular velocity in a thin disk. Alternatively, the central velocity dispersion of the gas implies that $`M_{}=1.0\times 10^9M_{}`$ if the gas is modeled as an isotropic spherical distribution of ballistic cloudlets. The stellar kinematics suggest that $`M_{}=(2.6\pm 1.1)\times 10^9M_{}`$. These different estimates bracket the value $`M_{}=(1.1\pm 0.3)\times 10^9M_{}`$ predicted by the $`M_{}`$-$`\sigma `$ relation (Ferrarese & Merritt, 2000; Gebhardt et al., 2000) as given by Tremaine et al. (2002).

To our knowledge, IC 1459 is only the second galaxy (Pinkney et al. 2000 modeled NGC 4697) for which a BH mass determination has been attempted with HST observations of both gaseous and stellar kinematics. Unfortunately, as it turned out, IC 1459 may not been the best galaxy to attempt such a comparison. However, it certainly is important to perform such comparisons. They provide insight into the reliability of BH mass determinations, which is necessary to properly understand the correlation between the BH mass correlation and other global galaxy parameters (Tremaine et al. 2002 and the refs therein). This is particularly relevant if one wants to understand the scatter in such relations and their implications for our understanding of the galaxy formation process.

With the availability of 2D spectroscopic detectors (e.g. SAURON on WHT, VIMOS on VLT, GMOS on Gemini) it is possible to obtain fully two-dimensional distributions of LOSVDs and line-strengths of kinematically decoupled cores (e.g. Davies et al., 2001; de Zeeuw et al., 2002). These observations will provide a much larger number of kinematical constraints than can be obtained from a small number of slit positions, which results in dynamical models that better constrain the model parameters (e.g. Verolme et al., 2002), especially if they resolve the radius of influence of the BH. In turn, this will improve the reliability of internal dynamical studies such as that presented in this paper.

Support for proposal 7352 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This work was initiated during a visit of MC to Leiden Observatory. We thank Karl Gebhardt for useful comments on an earlier draft of this manuscript.

## Appendix A MGE potential calculation

In this Appendix, we discuss the techniques that were used to efficiently and accurately evaluate the MGE potential. The MGE surface brightness defined in equation (2) can be deprojected analytically (Monnet et al., 1992) to derive an axisymmetric intrinsic luminosity density consistent with the observations. Using the equations in Cappellari (2002) the intrinsic luminosity density can be written as

$$\rho (R,z)=\underset{j=1}{\overset{N}{}}\frac{L_j}{(\sigma _j\sqrt{2\pi })^3q_j}\mathrm{exp}\left[\frac{1}{2\sigma _j^2}\left(R^2+\frac{z^2}{q_j^2}\right)\right],$$
(A1)
where $`N`$ is the number of Gaussians required in the expansion, $`L_j`$ is the Gaussian total luminosity, $`\sigma _j`$ the corresponding dispersion and the intrinsic axial ratio $`q_j=\sqrt{q_j^2\mathrm{cos}^2i}/\mathrm{sin}i`$, where $`i`$ is the galaxy inclination.

The general form for the potential that corresponds to an MGE density profile is

$$\mathrm{\Psi }(R,z)=G\mathrm{{\rm Y}}\sqrt{2/\pi }\underset{j=1}{\overset{N}{}}\frac{L_j𝒬_j(R,z)}{\sigma _j}$$
(A2)
(Emsellem et al., 1994), with

$$𝒬_j(R,z)=_0^1\frac{\mathrm{exp}\left\{\frac{T^2}{2\sigma _j^2}\left[R^2+\frac{z^2}{1(1q_j^2)T^2}\right]\right\}}{\sqrt{1(1q_j^2)T^2}}dT.$$
(A3)
Since the Gaussians of an MGE parametrization sample many orders of magnitude in radius, it is worthwhile to check if it is necessary to do the actual integration in equation (A3) or whether a limiting expression can be used instead. We desire a maximum relative error of $`ϵ\mathrm{\hspace{0.17em}10}^4`$ for the potential. This is at least two orders of magnitude more accurate than what can generally be achieved for photometric measurements.

An approximate expression for the potential at large radii can be found using a multipole expansion of the potential (e.g. Binney & Tremaine, 1987). Terminating the expansion at the first term, $`𝒬_j(R,z)`$ in equation (A3) equals

$$𝒬_j(R,z)=\frac{\sigma _j\sqrt{\frac{\pi }{2}}}{\sqrt{R^2+z^2}},$$
(A4)
and the potential reduces to the potential of a point mass with mass $`=\mathrm{{\rm Y}}L_j`$. Even in the thin disk limit, a maximum relative error $`ϵ\mathrm{\hspace{0.17em}10}^4`$ is obtained outside the sphere $`100\sigma _j\sqrt{R^2+z^2}`$. In the case of positive Gaussians this is a conservative estimate for the total relative error in the potential, since it applies only to one Gaussian, while in practice not all Gaussians in the expansion will be in the least favorable condition at the same time. This error limitation is generally not true when the MGE is composed by both positive and negative Gaussians, due to cancellation effects.

An approximate expression for the potential at small radii can be obtained by calculating a Taylor expansion of the $`\mathrm{exp}(x)`$ term in equation (A3) around $`x=0`$. The first term in this expansion equals

$$𝒬_j(R,z)=\frac{\mathrm{arcsin}\left(\sqrt{1q_j^2}\right)}{\sqrt{1q_j^2}}.$$
(A5)
This corresponds to replacing the potential by its central analytic value $`\mathrm{\Psi }(0,0)`$. Assuming a lower limit on the axial ratio $`0.1<q_j`$ a maximum relative error $`ϵ\mathrm{\hspace{0.17em}10}^4`$ is obtained inside the sphere $`\sqrt{R^2+z^2}\sigma _j/171`$.

An example of the usefulness of this approximation is given by the present application. Our orbit grid spans the radial range $`0\stackrel{}{\mathrm{.}}05<R<300\mathrm{}`$; a comparison with Table 4 therefore shows that these approximations are used often and only a small number of integrations has to be carried out. A related advantage of this is that the regions where the approximations are used are precisely those where the integrand is badly behaved and the numerical integration may become inaccurate.

Similar considerations are valid for the two components of the acceleration in the meridional plane, which are required by the orbit integration.

When the integral in equation (A3) has to be evaluated, we use a globally adaptive, recursive Gauss-Kronrod integration scheme (Favati, Lotti, & Romani, 1991). The use of such an algorithm is very beneficial at small and large distances from the galaxy center (but before the switch to the limiting regimes), where the integrand of equation (A3) shows a sharp peak close to $`T=1`$ or to $`T=0`$ respectively. In these cases the integrand is evaluated more densely around the peaks.

## Appendix B Integration of a Gaussian on a grid

Axisymmetric Schwarzschild modeling is aimed at finding the linear combination of orbits that reproduces the observed kinematics at various apertures, the intrinsic density $`\rho (R,z)`$ on the meridional plane and, to limit the effect of discretization, also the surface density $`\mathrm{\Sigma }(x,y)`$ on the projected plane. In general, the comparison with the density requires a double or triple integration. However, these computations can be simplified for an MGE parametrization. We discuss the different cases below. We always assume that the Gaussians have a maximum value of one, which means that intrinsic quantities have to be multiplied by $`L_j/[(\sigma _j\sqrt{2\pi })^3q_j]`$ and projected quantities by $`L_j/(2\pi \sigma _j^2q_j^{})`$ (see Cappellari, 2002, for details). In the following equations, the index has been dropped from the $`\sigma _j`$ and $`q_j`$ variables and $`q`$ has to be substituted by $`q^{}`$ when dealing with projected plane quantities.

### B.1 Polar grid in the meridional plane

We integrate an axisymmetric three dimensional Gaussian defined by $`\rho (R,z)`$ in equation (A1) on the “torus” $`D`$ of Figure 19 around the galaxy symmetry axis, which is delimited by the intervals $`[r_0,r_1]`$ and $`[\theta _0,\theta _1]`$ in spherical coordinates ($`\theta =0`$ on the symmetry axis). The integrated mass is given in spherical coordinates by

$$_D\rho =\mathrm{\hspace{0.17em}2}\pi _{\theta _0}^{\theta _1}_{r_0}^{r_1}e^{(r/c(\theta ))^2}r^2\mathrm{sin}\theta \mathrm{d}r\mathrm{d}\theta $$
(B1)
where $`c(\theta )=2q\sigma /\sqrt{1+q^2+(1q^2)\mathrm{cos}2\theta }`$. The innermost integral can be performed analytically and the integral becomes

$$\pi _{\theta _0}^{\theta _1}\left\{\frac{c(\theta )\sqrt{\pi }}{2}\left[\mathrm{erf}\left(\frac{r_1}{c(\theta )}\right)\mathrm{erf}\left(\frac{r_0}{c(\theta )}\right)\right]+r_0e^{(r_0/c(\theta ))^2}r_1e^{(r_1/c(\theta ))^2}\right\}c(\theta )^2\mathrm{sin}\theta \mathrm{d}\theta $$
(B2)
### B.2 Polar grid in the projected plane

In the case of the integration of a two dimensional Gaussian on a polar cell $`D`$, delimited by the intervals $`[r_0,r_1]`$ and $`[\theta _0,\theta _1]`$ in polar coordinates, we have to compute the integral

$$_D\mathrm{\Sigma }=_{\theta _0}^{\theta _1}_{r_0}^{r_1}e^{(r/c(\theta ))^2}rdrd\theta ,$$
(B3)
where $`\theta =0`$ on the Gaussian minor axis, which gives

$$_{\theta _0}^{\theta _1}\frac{c(\theta )^2}{2}\left[e^{(r_0/c(\theta ))^2}e^{(r_1/c(\theta ))^2}\right]d\theta .$$
(B4)
### B.3 Cartesian grid on the projected plane

We compute the integral of a two-dimensional Gaussian on a Cartesian grid cell $`D`$ delimited by the intervals $`[x_0,x_1]`$ and $`[y_0,y_1]`$. Assuming the cell is parallel to the Gaussian major axis, the integral is given by

$$_D\mathrm{\Sigma }=\frac{\pi q\sigma ^2}{2}\left[\mathrm{erf}\left(\frac{x_0}{\sigma \sqrt{2}}\right)\mathrm{erf}\left(\frac{x_1}{\sigma \sqrt{2}}\right)\right]\left[\mathrm{erf}\left(\frac{y_0}{q\sigma \sqrt{2}}\right)\mathrm{erf}\left(\frac{y_1}{q\sigma \sqrt{2}}\right)\right].$$
(B5)
If the Gaussian major axis makes an angle $`\theta `$ with the $`x`$ axis, the double integral on the grid cannot be evaluated explicitly, but can still be expressed as a single integral in one of the coordinates as follows

$$\frac{q\sigma \sqrt{\pi }}{p}_{x_0}^{x_1}[\mathrm{erf}(\frac{f_0(x)}{2pq\sigma })\mathrm{erf}\left(\frac{f_1(x)}{2pq\sigma }\right)]e^{[x/(p\sigma )]^2}dx,$$
(B6)
with $`p=\sqrt{1+q^2+(1q^2)\mathrm{cos}2\theta }`$, and we defined

$$f_k(x)=(1q^2)x\mathrm{sin}2\theta p^2y_k\mathrm{for}k=0,1.$$
(B7)
## Appendix C Hermite polynomials

We tabulate the first seven Hermite orthogonal polynomials, following the normalization of van der Marel & Franx (1993). They are written in Horner form for efficient and stable evaluation. The first five are also given in equation (A5) of that paper.

$`H_0(y)`$ $`=`$ $`1,H_1(y)=\sqrt{2}y,H_2(y)={\displaystyle \frac{2y^21}{\sqrt{2}}},H_3(y)={\displaystyle \frac{y\left(2y^23\right)}{\sqrt{3}}},`$
$`H_4(y)`$ $`=`$ $`{\displaystyle \frac{y^2\left(4y^212\right)+3}{2\sqrt{6}}},H_5(y)={\displaystyle \frac{y\left[y^2\left(4y^220\right)+15\right]}{2\sqrt{15}}},`$
$`H_6(y)`$ $`=`$ $`{\displaystyle \frac{y^2\left[y^2\left(8y^260\right)+90\right]15}{12\sqrt{5}}}.`$ (C1)