A STELLAR DYNAMICAL MEASUREMENT OF THE BLACK HOLE MASS IN THE MASER GALAXY NGC 4258

Here, we describe the acquisition and reduction of the photometric data, which are needed mainly to determine the stellar mass distribution and the kinematic-tracer profile, in the way described in Sections 5.1 and 5.4.

A journal of the imaging observations, from which the surface photometry was derived, is shown in Table 1. A wide variety of data are available in different colors, both from HST and the ground. For reasons explained in Section 4.1, we produced three luminosity profiles. (1) The near-infrared (NIR) profile was created by combining ground-based and HST J-, H-, and K-band images, as described in Section 2.1.1. (2) The V-band profile was generated from images in the V-band (HST) and in the R-band (ground based), as described in Section 2.1.2. (3) The I-band profile was created from the HST I-band image, as described in Section 2.1.2. Although this image was saturated in the nucleus, it was useful to consult for the luminosity weighting (normalization) of the line-of-sight velocity distributions (LOSVDs), as it identifies the kinematic tracer population (Section 5.4).

2.1.1. Reduction of the J, H, and K-Band Data

The Two Micron All-Sky Survey (2MASS) J-, H-, and K-band images (Figure 1, top row) were kindly provided to us by Tom Jarrett (Jarrett et al. 2003), and extend out to R ≈ 86 from the center¹⁴. For the morphological study of the large-scale properties of the galaxy (Section 3), the signal-to-noise ratio (S/N) was boosted by constructing a composite J + H + K image, which was then smoothed with a Gaussian point-spread function (PSF) having full width at half-maximum (FWHM) of 7'' for r ⩽ 60'' and FWHM = 11'' for r > 60''. The 2MASS K-band image was reduced in a manner similar to that of the J + H + K image. It proved to have a sufficiently high S/N to give a reasonably adequate measurement of the outer disk. A comparison of the K-band and J + H + K-band photometry provides an estimate of the uncertainties; these are dominated by the effects of patchy star formation and other irregularities, not by photon statistics. Therefore, we decided to discard the J + H + K-band profile as a tracer of stellar mass in the dynamical modeling, and to use the K-band profile instead, as being less affected by these irregularities.

Zoom In Zoom Out Reset image size

The Kitt Peak National Observatory (KPNO) 2.1-m telescope K-band image extends out to R ≈ 3' from the center (Figure 1, bottom row). Comparison of the 2MASS K-band profile with the KPNO K-band profile, which has a higher spatial resolution, showed that the 2MASS profile is affected by smoothing only at r < 32'', so the 2MASS profile was truncated inside this radius. Thus, the KPNO profile was used for 2'' < r < 32'' and the 2MASS profile at r ⩾ 32''.

Chary et al. (2000) obtained HST NICMOS images of NGC 4258 in bandpasses that approximate J, H, and K. Eric Becklin and Ranga Chary kindly made the reduced images available to us (Figure 2), along with the corresponding PSF images. Details of the reductions can be found in Chary et al. (2000). The HST WFPC2 V-band image (Section 2.1.2) shows extensive dust at r ≳ 2'' on the SW side of the center and a few thinner dust patches on the NE side. There is also a scatter of faint sources in the bulge, apparently associated with some ongoing star formation. In the NICMOS K-band image, these faint sources are not visible, and the dust patches on the NE side are essentially invisible as well. However, the dust on the SW side, although less conspicuous in K than in V, remains a substantial problem that we had to correct. After considering a number of alternatives, we chose the strategy of creating a symmetrized image of the galaxy, whereby the southwest (SW) side was replaced with the northeast (NE) side flipped (not rotated) across the major axis. As a check, the profiles from both the symmetrized and the unsymmetrized K-band images were compared and found to agree extremely well at r ≲ 2''. Therefore, this symmetrization is not a factor in the V − K color gradient discussed in Section 4.1. We applied the same procedure to the J- and H-band images, as well.

Zoom In Zoom Out Reset image size

The nucleus of NGC 4258 emits broad permitted lines (Ho et al. 1997), and is classified as a Type 1.9 Seyfert. The AGN at the center is very red (Chary et al. 2000) and, in the K band, it appears as a bright, unresolved point source that clobbers the central brightness distribution of the stars and needs to be subtracted. Unfortunately, the PSF provided by Chary et al. (2000) was determined using the TinyTim routine and not from actual images taken at the same time as the galaxy observations. It proves to be a mediocre match to the AGN at r < 02, in the sense that there is a residual diffraction pattern no matter how the PSF is scaled. We applied the following procedure to improve the quality of the PSF subtraction (see Figure 3). Using VISTA (Lauer et al. 1983), we scaled the TinyTim PSF image by various factors and measured the resulting surface brightness profiles in the residual image using the unsymmetrized galaxy image. We also examined the residual images to see which ones had the small-scale structure, illustrated in the top-left image of Figure 3, optimally subtracted. Most of the weight in the choice of the final PSF scale factor came from how well the first diffraction ring of the AGN was subtracted. The scale factor that we finally adopted in this way yields ellipticity and position angle profiles that are both smooth and well behaved, which is reassuring. This is not the case with scale factors that differ by more than ±10% from the adopted value. However, ultimately no value of the scale factor allows both a convincing subtraction of the first diffraction ring and a smooth profile at r < 0 2 that follows the power law seen at a larger radii. This is probably a result of the imperfect match of the TinyTim-generated PSF to the K-band data. It is unlikely that a real galaxy feature is involved, considering that no such feature is detectable in the higher resolution V-band image in the radial range 0 05 < r < 0 2 (see Section 2.1.2).

Zoom In Zoom Out Reset image size

We also applied 40 iterations of Lucy deconvolution (Lucy 1974; Richardson 1972) to the HST NICMOS J-, H-, and K-band images using the PSFs from Chary et al. (2000). Because the PSFs are not ideal, we applied the deconvolution to images with the PSF of the AGN subtracted as explained above. Still, the deconvolved K-band image showed bad "ringing," so we truncated it inside r = 0 46. The J-band profile appeared better behaved, at least partly since the AGN is much dimmer at this bluer wavelength, although there was also some sign of a residual AGN at r < 0 1. The J-, H-, and K-band deconvolved profiles agree very well at r ≳ 0 5, where the AGN is not a problem, and the ellipticity profiles are all well behaved.

Based on the preceding analysis, we decided (1) to use the original (not deconvolved) K-band photometry for r ⩾ 2 8 (KPNO data at 2 8 ⩽ r < 32'', and 2MASS data at 32'' ⩽ r ⩽ 560''), (2) to average the deconvolved NICMOS K- and J-band photometry for 0 2 ⩽ r < 2 8, and (3) to use deconvolved J-band data for 0 038 ⩽ r < 0 2 in order to avoid the AGN as much as possible. The resulting profile is shown in Figure 4. However, even in J, the points at r < 0 11 are probably still affected by the AGN. We elaborate on this point in Section 4.2.

Zoom In Zoom Out Reset image size

2.1.2. Reduction of the V- and I-Band Data

We used the MDM 1.3 m telescope to obtain an image of the galaxy in the R band (Figure 5). The image was reduced and dust-clipped in the usual way. The luminosity profile was then extracted and matched to the V-band HST profile, discussed below, in the range 5'' ≲ r ≲ 7'', using a V-band zero point. There were no significant color, ellipticity, or position angle (PA) gradients between the V-band and R-band profiles over the radial range in which they were matched. Henceforth, we will refer to this composite profile, made by stitching together the HST V-band profile with the R-band profile shifted by a constant offset, as the V-band profile. It extends out to r ≈ 150''.

Zoom In Zoom Out Reset image size

We obtained HST WFPC2 images of the NGC 4258 nucleus through the F547M filter under program GO–8591. The F547M filter was selected to exclude any bright emission lines associated with the AGN; it is roughly equivalent to a narrow V-band filter, and is referred to as such in the remainder of this paper. The galaxy was centered in the PC1 chip and dithered in a 2 × 2 square raster of ∼0.5 pixel steps over four separate 400 s exposures to maximize the spatial resolution of the complete data set. The four dithered images were combined into a single Nyquist-sampled "super image" using the algorithm of Lauer (1999a), which provides for the optimal combination of undersampled images without any associated degradation of the resolution. This final image has a pixel scale of 0 0228, twice as fine as that of the original images. However, the Nyquist image remains modulated by the HST PSF; thus, it was deconvolved prior to further analysis using 80 iterations of the Lucy (1974); Richardson (1972) algorithm. The PSF was estimated using the TinyTim package, but incorporates the WFPC2 pixel-response function recovered by Lauer (1999b). The final deconvolved V-band image of NGC 4258 is shown in Figure 6.

Zoom In Zoom Out Reset image size

The V-band surface brightness profile of NGC 4258, shown in Figure 4, was estimated from the F547M Nyquist image using a combination of the high-spatial resolution algorithm of Lauer (1985) and the least-squares estimator of Lauer (1986). The patchy dust seen throughout the bulge of NGC 4258 complicates measurement of the brightness profile; as is visible in Figure 6, there is a compact dust feature just slightly offset from the nuclear source that was especially problematic. For all dust features other than the central patch, pixels affected by dust could simply be excluded from the least-squares profile estimator. In practice, this was done iteratively by comparing the image to a model reconstructed from the profile to isolate faint dust absorption that may have been missed during the initial inspection of the image. Unfortunately, pixels could not be masked out of the high-resolution profile algorithm, which is used in preference to the least-squares estimator for r < 0 5, so the nuclear dust patch was first filled in with an estimated correction provided by the lower resolution least-squares algorithm. This initial correction was then refined by replacing it with an improved estimate provided by a model reconstructed by the high-resolution profile.

Figure 6 shows the central portion of the V-band image divided by a model reconstructed from the final brightness profile. Apart from the nuclear dust patch, the residuals are flat, showing that the brightness profile extracted is faithful to the central structure of NGC 4258. The nuclear dust patch, itself, is clearly compact, affecting only a single quadrant of the nucleus, which again allowed it to be isolated from the profile estimation. As a "sanity check" on the brightness profile, we compared the profile to an intensity trace taken along the bulge semimajor axis on the side of the nucleus opposite the nuclear dust patch. The two measures agreed extremely well, showing that the final profile is not strongly affected by any residual dust features or the pixels excluded from the analysis.

In addition to the HST F547M data, we also analyzed two archival PC1 F791W (I-band) images obtained by Cecil et al. (2000) under program GO–6563. Unfortunately, both images were saturated at r < 0 09 from the nucleus, and thus provide information only at somewhat larger radii. The images were combined and deconvolved; however, without a full accounting of the light contributed by the nucleus, the I-band profile should be regarded with caution for r < 0 5.

Using the F547M and the F791W images, we created a F547M/F791W color map of the nuclear region of the galaxy, approximating V − I (Figure 7). The color map shows more clearly the distribution of dust in the center, and will be further discussed in Section 4.1. A larger scale view of the central region of the galaxy is shown in Figure 8, which is a composite color image created by combining the F547M and the F791W images along with an archival F300W (U-band) image, also obtained by Cecil et al. (2000).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We used the Ca ii triplet absorption line (8498 Å, 8542 Å, 8662 Å; hereafter CaT) to determine the stellar LOSVDs at the center of NGC 4258, and at various positions along the major axis (out to R = 18 18) as well as along the minor axis (out to R = 11 69). The spectra were obtained from the ground as well as with the HST Space Telescope Imaging Spectrograph (STIS). Spectrograph configurations are shown in Table 2.

In particular, we made 24 cosmic ray split exposures of NGC 4258 with STIS to obtain high-quality nucleus-centered CaT spectra along the major axis of the galaxy, out to R = 0 80. The STIS slit was placed at the PA = 140° instead of directly on the major axis (PA = 150°) because of guide-star constraints. Our reduction procedure closely followed that of Pinkney et al. (2003). The procedure bypasses the pipeline (Calstis) reduction so that we can employ "self-darks," i.e., dark frames that are constructed by combining our dithered data sets with the galaxy light rejected. Our method works best if all of the exposures are taken within a short period of time because of the rapidly changing dark pattern on the STIS CCD. This was not quite the case here, since the exposures were taken on 2001 March 12 and 21. Nevertheless, a single self-dark worked well.

There is one small difference between our STIS-reduction procedure and that described by Pinkney et al. (2003): the initial dark subtraction was done with an archival, "weekly" dark rather than by a sum of the galaxy frames with a region around the galaxy masked out. This created a smoother dark on which to begin further iterations using the same "self-dark" strategy.

We also used the Modspec spectrograph on the MDM 2.4 m telescope to obtain ground-based CaT spectra along both the major and minor axes at larger radii. The configuration is again shown in Table 2. The seeing varied between 1'' and 15 (FWHM). Kinematic parameters were extracted out to R = 18 18 and R = 11 69 along the major and minor axes, respectively. Again, the reduction procedure closely followed that described by Pinkney et al. (2003).

Figures 9 and 10 show spectra extracted at several radial positions from the STIS and the ground-based data, respectively. The emission-line feature visible near 8620 Å (rest frame) in the nuclear STIS spectrum in Figure 9 is most probably the Fe ii 8618 Å line, and had to be removed before the LOSVD fitting could be done. Consequently, there was a concern that the central (or central few) LOSVDs would be less reliable, but it proved easy to interpolate under the line (see Figure 9), so in the end this was not a problem.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We decided to discard the central ground-based LOSVD along the major axis, because of concerns about template mismatch and AGN contamination. In retrospect, the central ground-based LOSVD on both the major and minor axes should have been dropped since the STIS data had comparable statistical quality and much better understood PSFs. In fact, we used the minor axis central LOSVD, except in one experiment, where its exclusion had no effect (see Section 5.3).

In a dust-free axisymmetric galaxy that also exhibits symmetry about the equatorial plane, such as assumed by our modeling method, kinematic quantities also manifest symmetries about the center. This entitles us to "symmetrize," with respect to the center of the galaxy, the LOSVDs along the major and minor axes, for both STIS and Modspec. This process helps alleviate potential biases in the extraction of the velocity profile. Gebhardt et al. (2003) provide more details on how this was done. Unfortunately, NGC 4258 contains substantial quantities of dust, especially near the dynamically important central region (Figures 6, 7), which calls into question the validity of this assumption. Potential implications are discussed in Section 6.

We used the maximum penalized likelihood method (MPL), as described by Pinkney et al. (2003), to derive LOSVDs from all the CaT spectra. A nonparametric form of the LOSVDs, rather than Gauss–Hermite moments, was used for the dynamical modeling. Each LOSVD was represented by 13 equally spaced velocity bins, with the uncertainty in each bin determined from Monte Carlo simulations.

Figure 11 illustrates V, σ, H3, and H4, the fitting parameters from the truncated Gauss–Hermite series expansions of the LOSVDs (Gerhard et al. 1993; van der Marel & Franx 1993), along the major and minor axes of the galaxy, both from STIS and from the ground. The same parameters in tabular form can be found in Table 3. We stress again that a nonparametric form of the LOSVDs was actually used for the dynamical modeling. We anticipate our best-fitting model discussed in Section 5.1 by illustrating its projected Gauss–Hermite moments here.

Zoom In Zoom Out Reset image size

Color gradients in the outer disk of NGC 4258 are small, even between the R and K bands (Figure 4), implying that the M/L ratio in the disk is not seriously affected by star formation or dust. However, inside ∼2'' from the nucleus, the V-band profile "peels off" from the NIR profile and becomes flatter. This color gradient can be seen in Figure 4, but becomes more apparent in Figure 12, which shows the inner 10'' of the surface brightness profile. Also, the color map shown in Figure 7 clearly identifies the color gradient as a "red," diffuse area surrounding the nucleus out to r ≈ 2'' (R ≈ 1''–2''). This is an issue for the dynamical modeling, because it takes place in the dynamically important central region, which affects the entire model: we must determine which, if either, of the two profiles (V-band or NIR) traces stellar mass.

Zoom In Zoom Out Reset image size

Could the color gradient be caused by variations in the stellar population? The strongest evidence against variations in the stellar population is the absence of any marked color gradient between the J and K profiles, which should be present if the central NIR brightening were caused by either a metallicity gradient or red supergiant stars. Furthermore, although some traces of scattered star formation can be seen on the V-band image (Figure 6) as well as on the color map (Figure 7), they are absent in the NIR images (Figure 2), even after accounting for the different FWHM of the PSFs.

More information can be obtained from the equivalent width (EW) profiles of the three CaT line components, which we computed from the STIS spectra extending out to 1 5 from the center (Figure 13). A decreasing EW profile would be evidence of a young stellar population or of a brightening I-band nonthermal continuum. A drop is seen in the central 0 3, and is discussed in the following subsection. Here, we are concerned with possible trends beyond that radius. For r > 0 3, there is a hint of a slow increase in the EW away from the center, but the data are also consistent with a flat profile. However, the data points at r > 0 8 suffer from the low S/N of the originating spectra, and we have no data as far out as ∼2'', where the V-band profile begins to peels off.

Zoom In Zoom Out Reset image size

The other possible explanation is the presence of diffuse dust, which would produce a small gradient in J − K but a larger one in V − J, just as observed. The observed smoothness in the V − J color gradient would then imply a similarly smooth increase in the projected density of diffuse dust toward the center.

In the end, we were unable to determine unequivocally whether the origin of the color gradient is a variation in the stellar population or the presence of a diffuse dust component. We decided to construct models for both of these possibilities, and investigate the uncertainties induced on the black hole mass (Section 5.4).

Dynamical modeling only strongly constrains the total mass, M_tot(r₀) = M_• + M_*(r₀), inside a small radius r₀ around the center, where M_*(r₀) is the enclosed stellar mass, and r₀ is comparable to the spatial resolution of the nuclear spectrum (∼01). Therefore, there is a degeneracy between M_• and M_*(r₀), which can be resolved only if we know M_*(r₀) separately, i.e., only if we know the luminosity profile and the stellar M/L ratio inside r₀. The latter is assumed to be constant throughout the galaxy and is predominantly constrained by the dynamical model at large radii (Section 5.1). In essence, M_• = M_tot(r₀) − M_*(r₀) corresponds to the "excess" dynamical mass inside r₀ after subtraction of the mass due to stars. Therefore, we must again determine which, if either, of the two luminosity profiles (V band or NIR) traces stellar mass, M_*.

The V-band profile, which is well determined down to r ≃ 0 011, manifests a progressively more shallow slope for r ≲ 0 15, interrupted by a peak within 004 of the nucleus. Although the cause of the peak is not entirely clear, it is unlikely that it is due to light from an old stellar population. The CaT EW profile (Figure 13) shows a ∼30% decline within 02 of the nucleus, which strongly suggests contamination by a continuum source, possibly OB stars or a nonthermal AGN.

Notwithstanding the uncertain origin of the peak, it seems unlikely that it traces stellar mass. Therefore, we replaced the innermost three V-band data points in a way that extrapolates inward the profile at 0 04 ≲ r ≲ 0 15 (Figure 12; Section 5.4). However, we also created models with the central peak to investigate its effect on M_•. Not surprisingly, the effect was negligible because the light in the peak corresponds to a very small mass (Section 5.4).

In Section 2.1.1, we saw that, at r ≲ 0 11, the J-band profile begins to show signs of contamination by AGN light, and the PSF of the AGN could not be subtracted reliably. We also saw earlier that the central depression in the CaT EW profile begins at r ≃ 0 2, and implies a brightening of the CaT continuum, which lies in the I band, inside that radius. The AGN is red, so it seems likely that the J-band profile is even more affected by the AGN continuum radiation than the I-band profile. While luminous AGNs normally have a blue featureless continuum, low-luminosity AGNs such as NGC 4258 tend to have emission spectra that peak in the mid-infrared (Ho 1999), probably as a consequence of their low accretion rates (Ho 2008). Taking all these considerations into account, we decided that the J-band profile inward of r ≃ 0 2 cannot be trusted to trace mass (i.e., the old stellar population). Accordingly, inside r ≃ 0 2, we replaced it with a profile that follows the V-band profile without the central peak. As explained above, this profile is presumably not affected by the AGN. The related dynamical models are discussed in Section 5.4.

A detailed account of our dynamical modeling method is provided in the Appendix. Here, we present a short description with emphasis on NGC 4258.

The first step involves a deprojection of the surface photometry to obtain the three-dimensional luminosity density. This is done under the assumption that emissivity is constant on similar, coaxial spheroids. The deprojection is performed via a nonparametric Abel inversion (see Gebhardt et al. 1996) assuming a value for the inclination angle (i = 72° in the base model). The gravitational potential can then be recovered by specifying a central point mass (black hole) M_• and the stellar M/L ratio, ϒ, assumed constant throughout the galaxy. The dynamical models for NGC 4258 extend out to r = 100''.

The dynamical models are constructed for a grid of (M_•, ϒ) values. Each model consists of a superposition of representative orbits, appropriately weighted to match the observed kinematics, subject to constraints that force the stellar density distribution to match the luminous density distribution of the galaxy. Each model is computed using N_o ≈ 7000 orbits. Orbits are integrated for about 100 crossing times.

Self-consistency is enforced by requiring that the model satisfies exactly the photometric constraints on a polar grid of 15 radial shells as described in the Appendix, each subdivided into four bins in a polar angle. The grid used to solve Poisson's equation is 4 times finer than the grid used to bin the kinematic results. For most orbits, energy is conserved to much better than 1%.

For each model, the orbit weights that best fit the observed kinematics are determined by minimizing

$$\begin{equation} \chi ^2 = \sum _k \left[ l_{k,o} - \sum m_{ik} w_i \right]\Big/ {\sigma _k^2}, \end{equation} \tag{ 1 }$$

where l_k,o is the light at a specific position and velocity (the index k runs over both position and velocity) with uncertainty σ_k, and m_ik is the mass deposited by the i^th orbit weighted by w_i in the k^th projected position and velocity bin; see Section 2.2). We use a maximum entropy technique as described in the Appendix to minimize χ² using non-negative orbit occupation numbers.

For each data set, we determine χ²(M_•, ϒ). The minimum determines the best (most probable) combination of M_• and ϒ for the galaxy. The Δχ² = 1 contour on the (M_•, ϒ) plane, where Δχ² ≡ χ² − χ²_min determines the nominal "1σ" uncertainty ranges for M_• and ϒ.

The process of converting observational quantities to model input entails assumptions and systematic uncertainties that have to be taken into account when determining the uncertainty of the "best" M_•, which may be larger than the nominal "1σ" confidence band on the (M_•, ϒ) plane. Therefore, the previous steps must be repeated for a number of plausible parameter choices (other than M_• and ϒ) in order to determine the overall uncertainty. We cannot afford to make an exhaustive exploration of the parameter space because of the time required to create the models. Instead, we opt to construct a base model corresponding to our "best guess" for the parameter values and investigate the effect on M_• of varying each parameter independently.

The results from the dynamical modeling are summarized in Table 4, and the corresponding χ² contour maps are shown in Figure 14. The black hole mass for the base model, which uses the V-band profile without the central peak for both the mass and the kinematic tracer profiles, is M_• = 3.31^+0.22_−0.17 × 10⁷ M_☉, with a total (not per degree of freedom) χ²_min ≈ 290. The internal velocity moments of the base model are illustrated in Figure 15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Dynamical Models

Model	Massa Profile	Tracera Profile	Kinem.b Data set	ic (°)	d	M•e (107 M☉)	ϒe	Sections
Base	Vcorr	Vcorr	S+G	72	0.35	3.31+0.22−0.17	3.6+0.1−0.1	5.1
D1	Vcorr	Vcorr	S+G	72	0.45	3.48+0.42−0.38	3.7+0.2−0.2	5.2.1
D2	Vcorr	Vcorr	S+G	62	0.35	3.62+0.38−0.49	3.8+0.1−0.2	5.2.2
K1	Vcorr	Vcorr	S1 + G	72	0.35	3.25+0.21−0.13	3.6+0.1−0.1	5.3
K2	Vcorr	Vcorr	S3 + G	72	0.35	2.20+0.54−0.31	3.6+0.1−0.1	5.3
K3	Vcorr	Vcorr	G	72	0.35	1.03+1.00−0.28	3.7+0.1−0.2	5.3
K4	Vcorr	Vcorr	S+G1	72	0.35	3.33+0.21−0.19	3.7+0.1−0.1	5.3
P1	NIRcorr	Vcorr	S+G	72	0.35	2.49+0.09−0.10	3.5+0.1−0.1	5.4
P2	Vcorr	NIRcorr	S+G	72	0.35	3.33+0.18−0.17	3.6+0.1−0.1	5.4
P3	NIRcorr	NIRcorr	S+G	72	0.35	2.43+0.21−0.30	3.5+0.1−0.1	5.4
P4	V	V	S+G	72	0.35	3.25+0.22−0.14	3.6+0.1−0.1	5.4

Notes. ^aMass and tracer profiles named as in Figure 12 (see Section 5.4). ^bS = All STIS kinematic data listed in Table 3. S₁ = As in S but without the central kinematic datapoint.S₃ = As in S but without the three central kinematic datapoints. G = All ground-based (Modspec) kinematic data listed in Table 3. G₁ = As in G but without the central minor-axis datapoint. ^cDisk inclination (see Section 5.2.2). Edge-on corresponds to i = 90°. ^dIsophote ellipticity, assumed constant throughout the galaxy (see Section 5.2.1). ^eBlack hole mass and stellar M/L ratio (ϒ), assumed constant throughout the galaxy (see Section 5.1). Uncertainty values correspond to the nominal "1σ" one-dimensional uncertainties from the contour maps (Figure 14). The M/L ratio calibration is in the V band.

Download table as:  ASCIITypeset image

The total number of parameters is 27 (LOSVDs) times 13 (velocity bins), which is 351. Typically, the reduced χ² is 0.3–0.5 in our models (Gebhardt et al. 2003), which is less than the expected value of unity because of covariance in the LOSVD components. For NGC 4258, the reduced χ² is near unity, but most of the contribution to χ² comes from the minor axis. The velocity centroid on the minor axis varies much more than the stated uncertainties (see Figure 11), whereas an axisymmetric model would have zero velocity along the minor axis. Regardless of whether this contribution is due to underestimated uncertainties or real nonaxisymmetric motion, it is the main cause for the inflated χ². Along the major axis, the reduced χ² is about 0.4, which is typical for the orbit-based models.

5.2.1. Isophotal Ellipticity ()

5.2.2. Disk Inclination Angle (i)

In Section 2.2, we mentioned that the presence of an emission line near the CaT in the nuclear STIS spectrum may have affected the reliability of the LOSVD extraction. In fact, the next two STIS spectra (at R = 0 05 and 0''10) show (weaker) signatures of the same emission line, as well. We also mentioned that we have some concerns about the quality of the central LOSVD along the minor axis in the ground-based data. In this section, we remove the suspect LOSVDs from the χ² fit, and examine the effect on M_•.

In models K1 and K2, we isolate the uncertainty introduced to M_• by the potentially unreliable LOSVD extraction. Model K1 is identical to the base model, except the nuclear STIS LOSVD is not included as a constraint. In model K2, the three innermost STIS LOSVDs are not included. In model K1, the effect on M_• is very small. This is reassuring, since it is the nuclear spectrum that would have been affected the most. Model K2 yields a significantly lower M_• = 2.20^+0.54_−0.31 × 10⁷ M_☉ but with a much higher uncertainty, which overlaps with the base model at the ∼1.5σ level.

An interesting question is whether HST spectroscopy is required to find evidence for a black hole. Using ground-based kinematic data alone should, of course, yield a black hole mass consistent with the mass obtained by also including HST data, albeit with larger uncertainties, depending primarily on the angular size of the radius of influence of the black hole on the sky. We address this question for NGC 4258 with model K3, which is identical to the base model, except only ground-based kinematic data are used. It turns out that the stellar M/L ratio can be constrained quite well, but the black hole mass, M_• = 1.03^+1.00_−0.28 × 10⁷ M_☉, is very uncertain. This is not surprising, considering that the radius of influence for the black hole in the base model is GM_•/σ²_e ≈ 0 32, where σ_e = 105 km s⁻¹ is the velocity dispersion in the main body of the bulge. This is substantially smaller than the spatial resolution of the ground-based kinematic data, for which FWHM = 1''–15. Therefore, HST kinematic data are critical for a well constrained estimate of the black hole mass in this galaxy. HST photometric data are also critical for galaxies where the luminosity profile near the center cannot be reliably extrapolated from the profile farther out, as is also the case with this galaxy (unfortunately, it is still not possible to determine unambiguously the inner mass profile, as discussed in Section 5.4).

We also constructed model K4, which uses all kinematic data (STIS+Modspec) except the central ground-based LOSVD along the minor axis (the central ground-based LOSVD along the major axis was not included in any model). We find that M_• is essentially identical to that of the base model.

In Section 4, we saw that stellar population gradients and/or dust extinction in the innermost ∼2'' of the galaxy, as well as light from the unresolved AGN, can induce variations in the measured stellar M/L ratio near the center of the galaxy. These have to be taken into account to determine which, if any, of the available luminosity profiles traces the mass–density profile, and which, if any, corresponds to the luminosity profile of the kinematic tracer population. The latter profile is needed for the luminosity weighting (normalization) of the LOSVDs. Unfortunately, we were unable to make these determinations convincingly, for the reasons explained in Section 4.

Consequently, we decided to identify three possible luminosity profiles, and to make models that use them as either mass or kinematic tracer profiles. The expectation is that these three profiles correspond to limiting cases, and that the true profiles lie somewhere in between. Then, the spread in M_• that we obtain from these limiting cases would be a good indication of the M_• uncertainty caused by our inability to determine the true profiles.

The profiles, which are depicted in Figure 12, are the following.

• 1.  
  The V profile. This is identical to the V-band profile as measured.
• 2.  
  The V_corr profile. This is identical to the V profile, with the three innermost points replaced by a shallower core, as explained in Section 4.2.
• 3.  
  The NIR_corr profile. For r > 0 2, this is identical to the spliced J, J + K, and K profiles in the top panel of Figure 4. Inward of 02, that profile is replaced by the one that parallels (shows no color gradient with respect to) the V_corr profile, for the reasons mentioned in Section 4.2.

If the observed color gradient is mostly or entirely due to the presence of diffuse dust or metallicity, then the NIR_corr profile is a better tracer of mass than the V or the V_corr profiles. However, the V-band profile has higher intrinsic resolution and a dimmer AGN, and the WFPC2 is better understood photometrically than NICMOS. For these reasons, we try all three versions as stellar mass profiles.

The kinematic tracer population was observed in the CaT wavelength, which lies in the I band. It would, thus, be natural to use the I-band profile as the kinematic tracer profile. However, we saw in Section 2.1.2 that the I-band image is saturated in the central 009, and thus cannot be deconvolved. The NIR profile shares the same slope with the I-band profile at larger radii (where I is not saturated), and could arguably be used instead. Unfortunately, it is also ill determined near the center due to the AGN and reduced resolution (Section 4.1). Consequently, we decided to again test both NIR_corr and V_corr as the kinematic tracer profile, and examine the induced uncertainty in M_•. For the base model, we adopted V_corr because of its higher resolution and smaller AGN contamination. This is somewhat arbitrary, but the choice of the kinematic tracer profile turns out to have a very small effect on M_•, as we will see below.

These scenarios were tested in models P1 through P4. Models P1–P3 substitute NIR_corr for V_corr in the mass and/or kinematic tracer profiles, and show that the choice of the stellar mass profile matters much more than the choice of the tracer profile. When V_corr is used for the mass profile (base model and P2), M_• = 3.31^+0.22_−0.17 × 10⁷ M_☉ or M_• = 3.33^+0.18_−0.17 × 10⁷ M_☉, depending on whether V_corr or NIR_corr is used, respectively, for the tracer profile. These numbers become 2.49^+0.09_−0.10 × 10⁷M_☉ and 2.43^+0.21_−0.30 × 10⁷ M_☉, respectively, when NIR_corr is used for the mass profile (models P1 and P3). The relative insensitivity to the choice of the tracer profile is not surprising, considering that, e.g., for a spherical system, the tracer profile ν(r) affects mass M(r) within the radius r only via dln ν/dln r, which by construction is nearly identical for both V_corr and NIR_corr near the center of the galaxy.

Model P4 substitutes V for V_corr in the mass profile. It yields M_• = 3.25^+0.22_−0.14 × 10⁷ M_☉, very similar to the base model, reflecting the fact that the only difference between P4 and the base model is the small luminosity excess within 004 of the nucleus, well below the spatial resolution of the kinematic data (∼01). The stellar M/L ratio is constrained at larger radii, so that the only difference in the stellar mass profile, M_*(r), between P4 and the base model is due to this small luminosity excess. From the M_• − M_*(r₀) degeneracy, discussed in Section 4.2, the only option for the modeling algorithm is to reduce M_• by a (small) amount equal to the difference in stellar mass inside r₀ ≈ 0 1 between P4 and the base model.

This point is illustrated in Figure 16, which shows the stellar mass profile, M_*(r), and the total mass profile, M_• + M_*(r), for the base model and for models P1 and P4. The mass profiles of these models correspond to the three profile choices (V_corr, NIR_corr, and V, respectively) discussed above. It is clear that dynamic modeling assigns black hole masses to P4 and to the base model such that they both contain the same total mass inside r ≃ 0 04 (since both models are constrained by the same kinematic data and the same V-band mass profile beyond r ≃ 0 04). This is not the case with model P1 which, although again constrained by the same kinematic data set, is characterized by a different mass profile (NIR_corr) outside of ∼0 1, the spatial resolution of the kinematic data.

Zoom In Zoom Out Reset image size

We construct the model on a grid in coordinates r, η, where r is the radius and η is the latitude. We restrict ourselves to potentials (and mass distributions) symmetric about the equator. The grid is evenly spaced in ν = sin (η) from ν = 0 to 1, and in r in even intervals in k defined by

$$\begin{equation} k = {1 \over a} \log \left(1 + {a \over b } r \right). \end{equation} \tag{ A1 }$$

Bins in (r, η) are bounded by the grid points. Below, we also use standard cylindrical coordinates (ϖ, z). It is essential to explore the phase space of orbits with a resolution at least as fine as the resolution used to construct the models.

Orbits in axisymmetric potentials always have two isolating integrals, the energy E, and the z-axis angular momentum J_z. Regardless of the possible presence of a third integral, conservation of E and J_z confine the orbit within a region of (ϖ, z) space defined by a boundary where the orbits' velocities are zero:

$$\begin{equation} \Phi (\varpi,z) + {J_z^2 \over 2 \ \varpi ^2} < E. \end{equation} \tag{ A2 }$$

The allowed volume intersects the equator only over a limited range in radius with an upper and a lower limits. We choose (E, J_z) pairs so that there is an orbit with upper and lower limits in the middle of each radial bin. We associate an area with each grid point (ΔE × ΔJ_z) by halving the distance to adjacent grid points. For each specified (E, J_z), the section of allowed phase space (we ignore ϕ because of the conservation of J_z) that lies on the equator is a two-dimensional surface with coordinates ϖ and v_ϖ (a surface of section). A regular orbit is defined by a curve on this surface. We systematically tile this surface in the manner invented by Thomas et al. (2004), launching an orbit from the center of each bin and assigning a phase-space volume to each orbit by

$$\begin{equation} \Delta \Omega = \Delta E \, \Delta J_z \int T (\varpi, v_\varpi) d\varpi \times d v_\varpi, \end{equation} \tag{ A3 }$$

based on the work of Binney et al. (1985), where T is the time between successive crossings of the equatorial plane.

The procedure above gives a starting position for each orbit in the equatorial plane of the model with a well defined ϖ, v_ϖ, E, and J_z. We use E to determine v_z, and we integrate the motion of the orbit in the (ϖ, z) plane. Assuming that all orbits cross the equator, it is a complete survey of orbits. We follow the orbits of stars in spherical coordinates to take advantage of the near sphericity of the mass distribution in the region where a black hole may dominate the gravitational field. Using r as the radial coordinate, η as the equatorial angle, and ϕ as the axial angle, the equations of motion are

$$\begin{equation} \frac{dv_r}{dt} = {v_\eta ^2 + v_\phi ^2 \over r} - {\partial \Phi \over \partial r}, \end{equation} \tag{ A4 }$$

$$\begin{equation} \frac{dv_\eta }{dt} = -{ v_r v_\eta \over r} - \tan \eta \, {v_\phi ^2 \over r} - {1 \over r} {\partial \Phi \over \partial \eta }, \; {\rm and} \end{equation} \tag{ A5 }$$

$$\begin{equation} v_\phi = {l_z \over r \cos \eta }. \end{equation} \tag{ A6 }$$

The first two equations are integrated using a standard fourth-order Runge–Kutta method with variable time step, and v_ϕ is recovered at each time step from Equation (A6). We obtained good energy conservation by letting dt = r/v_circ(r), where v_circ is the velocity of a circular orbit at that radius, with between 0.01 and 0.1. Orbits are followed for about 100 crossing times. We terminate the calculation of an orbit if its energy error exceeds 1%. Typical energy errors are 0.1%. The potential is set equal to zero at ∞.

The accelerations for the orbit integrations are determined by decomposing the assumed galaxy density distribution in spherical harmonics. Specializing immediately to axisymmetry gives

$$\begin{equation} \Phi ({\bf x}) = - 4 \pi G \sum _{l=0}^{l_{\rm max}} \left(A_l r^{-(l+1)} + B_l r^l \right) P_l (\nu), \end{equation} \tag{ A7 }$$

with

$$\begin{equation} A_l = \int _0^r \left[ \int _0^1 \rho (s,\nu) P_l(\nu) d \nu \right] s^l s^2 ds, \end{equation} \tag{ A8 }$$

$$\begin{equation} B_l = \int _r^\infty \left[ \int _0^1 \rho (s,\nu) P_l(\nu) d \nu \right] s^{-(l+1)} s^2 ds, \end{equation} \tag{ A9 }$$

and where ν = sin η. The gravitational acceleration on the orbit is determined by differentiating Equation (A7) with respect to r and η.

We tested our expansion in a Kuzmin–Kutusov model (Dejonghe & de Zeeuw 1988) with a/c = 0.5, which corresponds to a very flat galaxy (about E5). This model has analytical density, potential, and accelerations. Truncating the expansion at P₆ accurately reproduces the exact accelerations to 1% or better everywhere in the galaxy, so we set l_max = 6 in Equation (A7).

The additional acceleration from a mass point is straightforward to add to the radial component of the accelerations above. We store and model the galaxy in the angular range 0° ⩽ η ⩽ 90°. We store the three lowest order internal moments in each bin for each orbit (mass m, mv, mvv) for later use.

Simultaneously, the projected line-of-sight velocities are binned in a three-dimensional grid in projected radius R, angle β (measured from the major axis of the galaxy on the plane perpendicular to the line of sight), and velocity v. The spacing in R and β is defined as above. Since the orbits are being followed in r and η, ϕ is not defined. At each time step, in order to project the orbit onto the line of sight, we choose 100 values of ϕ, evenly spaced from 0 to π.

The next step is to find the superposition of orbits consistent with the stellar mass distribution that best matches the observed LOSVDs. Richstone & Tremaine (1988) provide an efficient method to match the mass distribution that can be augmented to minimize χ². Their method maximizes an objective function

$$\begin{equation} F = - \sum _i C_i(x_i), \end{equation} \tag{ A10 }$$

while satisfying a set of constraints

$$\begin{equation} Y_j = \sum _i m_{ij} x_i. \end{equation} \tag{ A11 }$$

We choose to maximize the objective function

$$\begin{equation} S - \alpha \, \chi ^2 = - \sum _{i=1}^{N_{\rm orb}} w_i \log \left({w_i \over \Delta \Omega _i} \right) - \alpha \sum _{k=1}^{N_d} {(l_k - l_{k,o})^2 \over \sigma _k^2}, \end{equation} \tag{ A12 }$$

where w_i is the weight of the i^th orbit and l_k is the light in the k^th bin in the observational phase space composed of the position on the sky and the line of sight velocity. The more familiar LOSVD is the distribution of light in line of sight velocity at a single position on the sky. Thus, the index k identifies a small range in both projected velocity and projected position. The variable α is an adjustable parameter discussed further below. There are N_orb orbits in the orbit library, and k varies from 1 to N_d, where N_d is the number of positions at which the LOSVD has been measured times the number of velocity intervals at each measurement. The first sum in Equation (A12) is an entropy, and the second is a χ² that measures the goodness of a fit of the model LOSVDs to the observed LOSVDs. The sum of the orbits must also match the deprojected (three-dimensional) light distribution in the N_bin bins:

$$\begin{equation} M_j = \sum _{i=1}^{N_{\rm orb}} m_{ij} w_i \ {\rm for} \ j = 1, N_{\rm bin}, \end{equation} \tag{ A13 }$$

where m_ij is the mass contribution of the i^th orbit to the j^th bin.

In order to put Equation (A12) in the form of Equations (A10) and (A11), we make the following substitutions and identifications. First, we define x_i = w_i for i = 1, N_orb, and add an additional set of variables that measure the error in the light projected into each LOSVD bin $x_{N_{\rm orb} + k} = l_k - l_{k,o}$ for k = 1, N_d. Since the light in each bin of the model LOSVDs is given by the contribution of each orbit to that projected velocity and position there is an additional set of constraint equations of the form the form

$$\begin{equation} l_{k,o} = \sum _{i = 1}^{N_{\rm orb}} m_{ij} w_i - x_j \ {\rm for} \ j = N_{\rm bin} + 1, N_{\rm bin} + N_d, \end{equation} \tag{ A14 }$$

where k = j − N_orb, and m_ij is the contribution of the i^th orbit to the k^th component of the set of LOSVDs.

Then, we have

$$\begin{eqnarray} C_i = & x_i \log \left({x_i / \Delta \Omega _i} \right)&\quad {\rm for }\quad i \le N_{\rm orb} \nonumber \\ C_i = & \alpha x_i^2 / \sigma ^2_{i - N_{\rm orb}} &\quad {\rm for }\quad N_{\rm orb} < i \le N_{\rm orb} + N_d\quad \end{eqnarray} \tag{ A15 }$$

defining C_i in Equation (A10), and Equations (A13) and (A14) correspond to Equation (A11). Note that for i > N_orb, −C_i is maximized at x_i = 0. The first set of C_i (for i ⩽ N_orb) effectively keeps the orbit's weight positive. The second set (for i > N_orb) minimizes the components of χ². We then use the method described by Richstone & Tremaine (1988) to iteratively solve the equations.

An important feature of this problem is the choice of the parameter α. Although we have not devised a method to specify appropriate values of this parameter in advance, if the quantity χ² reflects a satisfactory accounting of the estimates of errors in the recovery of the LOSVDs of the target galaxy, then the choice of α and interpretation of the results are reasonably straightforward. A small α will generally produce a solution. Increasing α slowly decreases χ², but not beyond some limit. In practice, we start with α = 0 and increase it slowly until the change in χ² is less than 0.02 in a single iteration.

In this case, we specify a dynamical model with a gravitational potential and a distribution function (DF) for the stars (the starlight need not, and does not, determine the mass in this model). The gravitational potential is

$$\begin{equation} \Phi (r) = v_{\rm circ}^2 \log (r) - {G {\it M}_\bullet \over r}. \end{equation} \tag{ A16 }$$

This potential includes a point mass and a flat circular velocity curve (with velocity v_circ) at large r. We investigated Michie-like DFs of the form

$$\begin{equation} f = A \exp \left\lbrace - \left[ {E + J_z^2\big/\big(2 \varpi _a^2\big) \over \sigma ^2} \right] \right\rbrace J_z^{2N} \sqcap (E_1, E, E_2), \end{equation} \tag{ A17 }$$

where A is normalization, E and J_z are energy and angular momentum, wϖ_a (an anisotropy distance), N, and σ (a velocity dispersionlike quantity) are model parameters, and E₁ and E₂ are upper and lower limits of the energy. The function ⊓ is defined by

$$\begin{equation} \sqcap (E_1,E, E_2) = \left\lbrace\begin{array}{@{}ll} 1, & \mbox{if } E_1 \le E \le E_2 \\ 0, & \mbox{otherwise}. \end{array} \right. \end{equation} \tag{ A18 }$$

The DF can be rewritten as

$$\begin{eqnarray} f &=& \pi A \varpi ^{2N} \exp (- \Phi /\sigma ^2) \sqcap (v_1, v_\phi, v_2) v_\phi ^{2N}\nonumber\\ && \times\,\exp \big(-v_\phi ^2\big/(2 \tilde{\sigma }^2)\big) \exp \big(-t^2/(2 \sigma ^2)\big) d t^2 dv_\phi,\nonumber\\ \end{eqnarray} \tag{ A19 }$$

where t² = v²_ϖ + v²_z, $\tilde{\sigma }^2 = \sigma ^2 (\varpi _a^2/(\varpi _a^2 + \varpi ^2))$ and $v_i = \sqrt{2(E_i - \Phi)}$. This function has a number of interesting properties: N > 0 tends to create v_ϕ enhanced anisotropy and a flattened density distribution; $\tilde{\varpi }_a < \infty$ tends to depress v_ϕ and creates a prolate density distribution. The cutoffs in energy avert a divergence in the luminous matter distribution near the black hole.

The density is

$$\begin{eqnarray} \rho &=& 2 \pi \sigma ^2 A \varpi ^{2N} \exp (- \Phi /\sigma ^2) \int _{v_1}^{v_2} v_\phi ^{2N} \exp \big(-v_\phi ^2/(2 \tilde{\sigma }^2)\big) \nonumber\\ && \times\,\big[\exp \big(-t_1^2/(2 \sigma ^2\big)\big) - \exp \big(-t_2^2/(2 \sigma ^2)\big)\big] dv_\phi,\qquad \end{eqnarray} \tag{ A20 }$$

where t₂ = v²₂ − v²_ϕ and t₁ = max(v²₁ − v²_ϕ, 0).

Note that the distribution functions with different choices of parameters can be added to each other (in the same potential), and similarly the densities and density-weighted moments are additive. A Monte Carlo realization of the model is constructed by sampling each density distribution and choosing the velocities randomly from the velocity distribution. Note that we have the additional freedom of choosing an arbitrary fraction of the ϕ velocities in the positive sense, and hence setting the spin of any model between a maximum amount and zero. The total number of random variates chosen sets the normalization.

For the test here, we adopted a potential with v_circ = 220 kms⁻¹ and M_• = 1.126 × 10⁸ M_☉. We added two DFs. The first had σ = 160 kms⁻¹, ϖ_a = 600 , N = 0, E₁ = Φ(10 ), E₂ = Φ(1000 ), with equal numbers of the ϕ velocities chosen in the positive and negative sense. The second DF had σ = 120 kms⁻¹, w_a = 200 , N = 2, E₁ = 10 , E₂ = 1000  and 3/4 of the ϕ velocities were chosen positive. Each DF was assigned 1/2 × 10⁹ random points, so they had equal mass. The second DF, featuring a moderately spinning disky structure, dominates the model near the center of the galaxy, but declines rapidly. The first DF, a mildly prolate structure with no net spin dominates further out. No stars are found beyond a radius of 1000  from the center of the model. The 10⁹ random points each are converted to projected positions and velocities (we assumed the galaxy was edge on), and binned into "observations" mimicing our combined HST and ground-based data, and also producing data mimicing data that would be obtained with an integral field unit. For the example, below we used the IFU data with a resolution of 10 at the center, and lower resolution data extending to a distance of 800  from the center. At each location, the projected velocity variates were binned at a resolution of 20 kms⁻¹. Gaussian noise was added to the data in a series of 20 realizations. These data sets were then fed into the modeling program, which produced (in each case) a χ² map in v_circ and M_•.

Marginalizing over each variable, in turn, produced the determinations of M_• and v_circ illustrated in Figure 17. Each data set had LOSVDs at 32 positions with 15 velocity bins, a product of 480. For each realization of this set, the best-fit model had a χ² near 500. Since the velocity bins are uncorrelated in these data sets, this is about the expected χ². For the presentation in Figure 17, the plots of χ² against black hole mass and v_circ were arbitrarily shifted to a minimum of 500. The error bars judged from the individual χ² profiles are about 20% smaller than the error (above) determined from the spread in minima. The mean values of M_• and v_circ from these experiments were (1.15 ± 0.03) × 10⁸ M_☉ and 222 ± 5 kms⁻¹, in excellent agreement with the target values.

Zoom In Zoom Out Reset image size

An earlier version of the code suffered from a units conversion error that led to an underestimate of black hole masses (and M/L). All black hole masses obtained with this program published in or prior to 2003 should be multiplied by 1.09.