THE BLACK HOLE MASS IN M87 FROM GEMINI/NIFS ADAPTIVE OPTICS OBSERVATIONS

We observed the central region of M87 in queue mode using the NIFS (McGregor et al. 2002) on the Gemini Telescope. We used AO corrections with the laser guide star system, ALTAIR (Boccas et al. 2006). The AGN in M87 provided the low-order corrections for tip–tilt and focus in a manner similar to Krajnović et al. (2009). An important feature in M87's center is the nearly point-like knot, located about 1'' off the nucleus along the outbound jet, named HST-1 (Harris et al. 2003; Perlman et al. 2003; Cheung et al. 2007; Madrid 2009), which allows us to monitor the telescope's point spread function (PSF). The PSF serves as an important input to the kinematic modeling. The data were taken over five nights in 2008 April and May with 23 dithered positions of 10 minute exposure each on M87. The telluric standards HIP 59174 and HIP 61138 were observed to monitor and correct for atmospheric absorption. We used the K_G5605 grating that provides wavelength coverage from 2.00 to 2.43 μm, with a spectral resolution of 5290 over a field of 3'' × 3'' sampled at 004 north–south and 0103 east–west across the image slices.

We used the Gemini NIFS package of IRAF (Tody 1993) scripts (developed mainly by T. Beck) for the majority of the reductions. This package provides the flattening, registration (spatially and spectrally), hot pixel identification, and sky subtraction. This package does not yet handle the error frames appropriately, so we also passed our original uncertainty frames (dark current, read noise, and Poisson noise) through the same reductions as the science data as a first modification. Our second modification involves interpolating over the telluric Brackett γ line and dividing the telluric standard by a 10⁴ K Planck function; this modification to the telluric correction retains a proper relative spectral shape in the science data. We also skipped the final NIFS script step where the data are resampled to equal x- and y-steps since we found that this interpolation enhanced the residual structure from our PSF fits.

We generally took sky exposures before and after the on-target frames, although there were some on-target frames that only had one sky exposure. We used the sky nearest in time for each on-target exposure. The sky subtraction usually worked well judging from the inspection of residuals under sky lines and tests with subtraction between different pairs of sky exposures. However, it is clear that some atmospheric emission-line variability occurred between our sky nods and caused uncertainty beyond our direct, statistical noise. There are techniques to bundle atmospheric emission lines into common transition families and fit a series of scalings between science and sky frames to minimize the residuals (Oliva & Origlia 1992; Rousselot et al. 2000; Davies 2007), but we have not employed them here. The public CO line lists extend only to 2.27 μm and therefore do not cover the CO bandheads through which we measure all our kinematics. However, since we have so many sets of exposures, with each set at a different dither position, the problematic sky regions are mitigated. Furthermore, for the spectral extractions we mask out wavelengths near the CO bandheads that have large variations between sets of sky exposures above the thermal background. The final product from the NIFS reduction package is a wavelength-calibrated spectrum for each integral field unit (IFU) element at each of the 23 different dither positions. We next find the relative position for each exposure, the PSF, and remove the AGN and jet continua if present.

A crucial step for dynamical modeling with AO data is to determine the PSF and, in the case of M87, to remove non-stellar features from the spectra. For galaxies with shallow light profiles like M87, determination of the PSF is particularly important, since one has to know how much light from the outer regions is in the central spatial elements. Fortunately, for M87, we are in the excellent situation of having the point source (HST-1) within the field that can be used as a measure of the PSF. However, the central AGN is so bright that it significantly contaminates the stellar spectra within the central few spatial elements.

There are a few techniques to estimate the PSF with AO data outlined in Davies et al. (2004), which we considered. Davies et al. (2006) present observations where the PSF is measured from Brackett γ in an unresolved AGN. Since M87 has Hα in the central regions, it should have some Brackett γ; however, the redshifted line falls, to within one pixel, on one of the brightest sky lines at 21798 Å. The residuals from the sky line are strong enough to not detect Brackett γ emission. The only observational handles we have on the AGN and jet flux contributions to any particular pixel are the spatial brightness variations and the change in the spectral slope. The AGN and jet are intrinsically much redder than the stellar populations. We wish to use this information without making assumptions about the stellar surface density. Thus, another goal is to study the stellar profile within the region dominated by the AGN; using the integrated light profile does not allow one to do this, whereas that information is contained in the spectra.

As opposed to measuring the light model (PSF and AGN fraction) from the reconstructed IFU data, we could use the stellar light profile as measured from HST. There are three reasons for not forcing the light profile from a previously measured HST image. First, the jet has knots that are moving on short timescales, and we cannot be sure that the AGN and jet fractional light and the spatial position will be the same. Second, we desire to use the spectral information to attempt to measure the underlying stellar component into the center. Third, the light profile and AGN fraction depend on the specific color. Data with K-band filters (Corbin et al. 2002) show color profiles that are flat with radius, although the analysis cannot easily go below 1'' radii.

We make the assumptions of a constant stellar population in the center, with a spatially constant spectral slope and CO equivalent width (EW) only altered by the relative AGN/jet contamination. These assumptions form a closed constraint on the AGN/jet components and the PSF. Similar to Lauer et al. (1992) with visible light HST data and the approach of CLEAN algorithms (Högbom 1974), we model the inner AGN jet as a set of point sources. Additionally, we fit PSFs to the telluric standards as a verification of our in situ PSF models. We find in both that the sum of an inner, AO-corrected, non-circular function and an outer, natural seeing function fit the data well without spatially coherent residuals. The two-component PSF is common with these types of data, but circular PSFs are commonly assumed (Neumayer et al. 2007).

It is important to get a robust spatial model for the AGN, jet, and stellar light since the kinematics in this spectral region are sensitive to the EW of the CO lines. Silge et al. (2005) show that a mismatch between the EWs of the velocity templates with the data can bias the velocity dispersion either high or low by up to 30%. Given that the additional continuum of the AGN will dilute the EWs, it is important to use as much information as possible to constrain the relative contribution.

Thus, (1) we assume that the stellar population (and therefore color and spectral slope) do not vary with radius near the center, and (2) we treat the AGN and jet as a set of point sources with unknown brightnesses with a flat continuum. We use the following operational definition of the CO EW:

$$\begin{equation} \rm EW_{CO}=\Delta \lambda \displaystyle \sum _{\lambda =2.29\,\mu \rm m}^{2.42\,\mu \rm m} \frac{D(\lambda)- a_{{\rm tot}}\times \lambda ^{\alpha _{{\rm tot}}}}{a_{s}\times \lambda ^{\alpha _{s}}}, \end{equation} \tag{ 1 }$$

where D(λ) are the counts in each pixel, a is the fitted zero point for the continuum (from a power-law fit), α is the fitted power-law exponent for the continuum, and a_s and α_s are the zero point and spectral index for the stellar light, respectively. Note that under assumption (1), this EW should not vary with position.

We begin the analysis of each science frame by considering all pixels outside of 09 from the center of the AGN and 03 from the center of HST-1; in this way the fit to the stellar model uses a relatively pure stellar signal. We make a least-squares power-law fit from 2.1 to 2.27 μm to each pixel and find a robust estimate for EW_CO and α_s from a bi-weight calculation (Beers et al. 1990). The estimate of any pixel's stellar continuum strength, a_s, can then be found with a least-squares power-law minimization for a_tot and α_tot followed by the application of Equation (1). The EW_CO and α_s over all frames are 220 ± 26 Å and −3.11 ± 0.43. We make continuum and EW maps for all pixels, subtract off the non-stellar continuum, and normalize by the stellar continuum extrapolated through the CO bandheads. This procedure produces the reduced spectra which are later used for extraction of the kinematics.

2.2.1. PSF Model

The PSF determination requires further analysis. We smooth the stellar continuum intensities using a 02 × 02 boxcar. We assume that the PSF has the form of an anisotropic Gaussian plus a power law (in the form of a Moffatt function), given by

$$\begin{eqnarray} {\rm PSF}(x,y) &= N(x,y) + M(x,y), \nonumber\\ N(x,y) &= \frac{a_1}{2\pi a_2^2 a_4} \exp\left(-\left(x_c^2+\frac{y_c^2}{a_4^2}\right)\bigg/{2a_2^2}\right) \nonumber \\ M(x,y) &= (1-a_1)(a_6-1)/\pi a_5 \left(1+\frac{x^2+y^2}{a_5^2}\right)^{a_6} \\ x_c &=x\cos (a_3)+y\sin (a_3) \nonumber \\ y_c &=y\cos (a_3)-x\sin (a_3), \nonumber \end{eqnarray} \tag{ 2 }$$

where PSF(x, y) is the value of the PSF at a given position x and y. N(x, y) is the Gaussian model for the inner PSF, with normalization a₁, width a₂, position angle a₃, and axis ratio a₄. M(x, y) is the Moffatt model for the outer PSF, with normalization (1 − a₁), width a₅, and exponent a₆. In this model, the PSF is assumed constant over the NIFS 2'' field; this approximation is very accurate for a laser guide star system at a wavelength of 2 μm.

Due to the large number of parameters needed to solve for the PSF (shape parameters for the PSF, multiple components for the AGN/jet, and stellar profile parameters), these fitting procedures involve minimizing a complicated function with local minima. We resort to simulated annealing as a global minimization tool (Press et al. 1992) with temperature schedules that reduce by 30% over each iteration and terminate with 10⁻⁴ fractional convergence. We first perform a simultaneous fit to the stellar-continuum-subtracted data with three central AGN/jet point source components, a fourth point source component for the HST-1 clump, and PSF parameters given by Equation (2).

We subsample each PSF evaluation over each pixel by five E-W and three N-S since the individual IFU elements do not properly Nyquist sample the PSF. This fit determines reasonable locations and strengths for all source components, but it improperly lets the central AGN/jet drive the PSF fit, and we know, in fact, that this feature is resolved given the multiple components. We therefore re-minimize the PSF terms and the HST-1 terms with data within 12 of the preliminary HST-1 position, but outside of 05 of the AGN center. The isolated clump, HST-1, then delivers a clean PSF determination. Finally, we re-minimize all source components, but hold the PSF terms fixed across the whole map. A representative decomposition example is shown in Figure 1 for one of the 23 data sets, where we show the raw IFU data, the stellar model, AGN/jet model, and residuals. The bottom right panel shows the residuals plotted on a linear scale. There remains some structure in the residuals, but the maximum residual is less than 1% of the measured value. We have tried a variety of PSF models and additional point source models, and find no improvement. Given the tightness in the HST-1 position determination, we register all frames off this cleanly isolated feature. From the median of all frames, we estimate AGN and HST-1 spectral indices of −0.67 ± 0.52 and −1.9 ± 2.6, respectively.

Zoom In Zoom Out Reset image size

This computation also delivers the range of PSFs for the 23 data sets. Most of the PSFs are close to the diffraction limit of 006 for the inner component, and the full range of the FWHM for the inner component is 0.055–0.19. The fraction of light in the inner component ranges from 0.14 to 0.45. The fraction of light in the central component is indicative of the Strehl ratio; however, in practice, the Strehl ratio is hard to measure given uncertainties in the PSF model (Gebhardt et al. 2000c). We make a two-dimensional image of each of the 23 analytic individual PSFs. From these 23 images, we then average to make the two-dimensional array which we use for the PSF in the dynamical modeling.

Figure 2 plots the flux, ellipticity, and position angle versus radius for the combined PSF. The values are reported in Table 1. This PSF has the inner Gaussian FWHM of around 008 with a fraction of the light in the central component of a₁ = 0.38. We use the analysis of the PSF as measured from the reconstructed IFU data, and we have also compared this PSF to that measured from the telluric standards. We find a similar FWHM for the inner and outer components, but the fraction of light in the inner component changes. For the PSF measured from the tellurics, the amount of light in the central component ranges from 60% to 70%, which is 1.5 to a factor of two larger than that determined from the science frames. We argue that using the telluric PSF is too optimistic for multiple reasons. First, the telluric star is used as the reference star for the PSF and tip/tilt corrections (natural guide star mode), whereas the M87 data use the laser as the reference. Second, the M87 data use the nucleus for the tip/tilt corrections, which is more extended than the telluric star. Third, the M87 frames come from long exposures of 10 minutes compare to exposure times of only a few seconds for the telluric. Thus, the M87 PSF is expected to be more extended. In any case, we also run a full set of dynamical models using the telluric PSF and find insignificant differences. It is encouraging to see that both PSFs give similar results; we attribute this robustness to the well-resolved kinematics in the region influenced by the BH.

Zoom In Zoom Out Reset image size

Table 1. PSF Parameters

R	Flux	Ellipticity	P.A.
''			°
0.00	1.00	0.169	170.4
0.02	0.96	0.169	170.4
0.04	0.86	0.166	170.4
0.06	0.72	0.161	170.4
0.08	0.56	0.154	172.7
0.10	0.42	0.142	170.9
0.12	0.31	0.125	167.5
0.14	0.22	0.105	161.6
0.16	0.17	0.087	151.3
0.18	0.13	0.086	139.0
0.20	0.11	0.084	128.8
0.22	0.090	0.083	122.1
0.24	0.076	0.083	117.9
0.26	0.065	0.074	114.9
0.28	0.056	0.061	112.7
0.32	0.043	0.035	109.9
0.36	0.034	0.018	107.7
0.40	0.027	0.009	105.9
0.44	0.021	0.004	104.3
0.48	0.017	0.002	103.0
0.52	0.014	0.001	101.7
0.56	0.011	0.000	100.6
0.60	0.010	0.000	99.6
0.64	0.008	0.000	98.7
0.68	0.006	0.000	97.7
0.72	0.005	0.000	96.6
0.76	0.005	0.000	96.1
0.80	0.004	0.000	94.3
0.86	0.003	0.000	92.2
0.94	0.002	0.000	90
1.00	0.001	0.000	90

Download table as:  ASCIITypeset image

2.2.2. Central Stellar Distribution and Offset With AGN

To align the data with the kinematic axis at large radius from the fit of Kormendy et al. (2009), we take a position angle of −25° (E of N) for the major axis. Figure 3 shows the radial and angular bins used in the modeling, following the same binning of Gebhardt & Thomas (2009). Figure 3 only plots the spatial region for the NIFS data (the model goes out to 2000''), with a gray-scale image from one of the 23 reconstructed IFU images. Spectra on opposite sides of the major axis are combined. We generally have 10 spectra at a given radius, since we have five angular bins on each side of the minor axis. The two central radial bins are not used due to AGN contamination (discussed below), and the next two outer radial bins require a sum over all angles in order to obtain adequate signal-to-noise ratio (S/N). The total number of spectra from the NIFS data that are used for dynamical modeling is then 40. Table 2 provides the locations for these 40 bins. We also include the signal to noise per pixel for each bin.

Zoom In Zoom Out Reset image size

Figure 4 shows the spectra for three different spatial regions. The top spectrum comes from radii 008 < R < 018, where we have summed over all angular bins; the middle is from 018 < R < 03, and the bottom is from a radius of 06. The wavelength range shown here is the region used for the kinematic extractions. There are no significant absorption lines on either side, although the blueward region is used to help estimate the relative contribution of AGNs and stars (as discussed previously). The spectrum in the middle panel of Figure 4 represents the innermost point used in the dynamical models. The black line shows the data, and the dashed black lines are regions not used in the kinematic extractions due to large variations in the night sky. In this wavelength region, the residuals due to sky subtraction tend to be positive, even though the subtracted sky frames use the same exposure time as the on-target frames. We have tried different sky subtraction levels and find little difference in the kinematic extractions because we excise the regions with sky lines. In the other wavelength regions, the sky residuals average to zero. The red line in the plots is the fitted LOSVD convolved with the template. The template comes from a library of 10 stars observed with GNIRS, as reported in the GNIRS template library (Winge et al. 2009). We select stellar types from G dwarfs to M supergiants. We rely on the GNIRS template library as opposed to the NIFS template library, since the GNIRS library contains a larger range of spectral types. The template library is an important consideration for this wavelength region (Silge et al. 2005).

Zoom In Zoom Out Reset image size

The kinematic extraction simultaneously fits a non-parametric LOSVD and template weights for the individual stars. The template composition is allowed to vary spatially. The technique is described in Gebhardt et al. (2000a) and Pinkney et al. (2003). The LOSVD is defined in 15 velocity bins of 260 km s⁻¹. There is a smoothing parameter applied to the LOSVD, but given the high S/N for most of the spectra, the smoothing has little effect on the extractions; thus, there is only a modest correlation between adjacent velocity bins. We use Monte Carlo simulations to determine the uncertainties in the LOSVD. The S/N of each spectrum determines the noise to be used in the Monte Carlo simulations; from 1000 realizations of each spectrum, we generate an average LOSVD and the 68% uncertainty.

The dynamical modeling uses the non-parametric LOSVD directly. However, it is sometimes convenient to express the LOSVD in a parameterized form as Gauss–Hermite moments, to show the radial run of the kinematics, and to compare the data with the models. Table 2 shows the first four Gauss–Hermite moments for the NIFS data. Figure 5 plots the velocity dispersion versus radius, where the dispersion is measured from a Gauss–Hermite fit to the LOSVDs. Figure 5 plots all of the data at each radius, and there are between 1 and 10 angular bins at each radii; thus, there are multiple points at nearly all radii in the figure. There is no rotation seen at a significant level in the NIFS data.

Zoom In Zoom Out Reset image size

We input 107 LOSVDs in the dynamical models. These LOSVDs come from 40 spatial bins from the NIFS data, 25 from the SAURON data, and 42 from the large radial data of Murphy et al. (2011). The data in Murphy et al. come from the IFU VIRUS-P, where we have nearly complete angular coverage. The S/N of those data is very high (50–100 per resolution element). The solid line in Figure 5 plots the velocity dispersion from the best-fit dynamical model. The model generates LOSVDs, and their dispersions come from Gauss–Hermite fits to those LOSVDs. For the dynamical model dispersions, we average along angles at a given radius for clarity. In Figure 5, we plot both the NIFS and VIRUS-P dispersions, which have different PSFs. The model is convolved to each of the PSFs, and the plotted dispersions include the convolution.

Data within the central 018 are excluded. Within R < 008, the number of individual NIFS spatial pixels is only 50, whereas the number of spatial pixels for bins used in the models ranges from 250 to 3000. Given the shallow surface brightness profile for M87 and the low number of NIFS pixels, the signal from the central stars is low, and contamination from the AGN is high. We have tried a variety of models for the AGN, PSF, and stellar light profile; in all cases, we find that little information is contained in the central spectrum. We do not further discuss this spectrum.

The spectrum coming from the radial region 008 < R < 018 has higher S/N, but still low enough that the kinematic extraction is highly uncertain. Figure 4 plots this spectrum in the top panel. It has lower stellar S/N compared to any spectra we use, and is further compromised by the uncertainty in the AGN subtraction. However, we still attempt a kinematic extraction. The red line in Figure 4 is the best-fit convolved template from the region one radial bin further in radius, from 018 < R < 03. This region has a velocity dispersion of 480 km s⁻¹. There are wavelength regions, for example at 2.31 μm <λ< 2.35 μm, where the model and data are offset. We do not attribute this offset to poor LOSVD modeling, but instead to poor AGN and stellar light discrimination. We use this region as an example and extracted LOSVDs including and excluding various regions. The range in velocity dispersion over all tests is from 350 to 620 km s⁻¹. The uncertainties on the dispersion for each individual extraction, coming from Monte Carlo simulations, are much smaller than this range, indicating that we are dominated by systematics as opposed to noise. For these reasons, we exclude this spectrum. We note that our best-fit dynamical model predicts a velocity dispersion of 451 km s⁻¹ at this location (the solid line in Figure 5 shows the predicted dispersion from the model). This value is in the middle of our range of dispersions using the various extractions. For the central radial bin (R< 008), the model prediction is 430 km s⁻¹.

In Gültekin et al. (2009), we assign M87 a velocity dispersion of σ_e = 375 km s⁻¹, which in turn comes from the analysis of Gebhardt et al. (2000b). In Figure 5, however, we see that this value is reached only at r < 2'', a location that is clearly within the inwardly rising portion of the dispersion profile associated with the BH's kinematic influence; this value unlikely represents the M87 galaxy overall.

The velocity-dispersion parameter used in the BH–σ relation is the effective velocity dispersion, σ_e, which in Gebhardt et al. (2000b) is defined as $\sigma _e^2 = \int _0^{R_e} I(r)V^2(r) dr/\!\!\int _0^{R_e} I(r) dr$, along the major axis, where I is the surface brightness, V is the projected second moment of the velocity distribution, and R_e is the half-light radius, which for M87 is 100'' as reported in Lauer et al. (2007) and Kormendy et al. (2009). This operational definition of σ_e appears to provide a good correlation with BH mass, but there are many different ways in which one can integrate the kinematics in order to provide one number for the galaxy. With this definition σ_e = 352 km s⁻¹, using the kinematics and surface brightness profile presented in this paper. The previous value of 375 km s⁻¹ results from using the older kinematic and light profiles. It is clear that σ_e contains a substantial contribution from the light inside where the BH influences the kinematics. If instead we exclude radii within this region (defined as r_s = GM_BH/σ² and equal to 21 for our models) from the integral that determines σ_e, we find σ_e = 324 km s⁻¹, about 8% smaller. We choose 324 km s⁻¹ as our best estimate of σ_e, with a range from 312 km s⁻¹ to 352 km s⁻¹.

Churazov et al. (2010) show that there exists a radial "sweet spot", where the velocity dispersion at that radius is robustly related to the circular velocity. By providing a dispersion value that is indicative of the galaxy as a whole, this estimate may correlate well with the BH. Based on the kinematics from Murphy et al. (2011), the dispersion value of the "sweet spot" for M87 is 312 ± 10 km s⁻¹. Furthermore, Cappellari et al. (2007) measure a value of 306 km s⁻¹ by integrating the two-dimensional data within a radius of 30''. There are a variety of ways to represent a velocity dispersion for a galaxy, and until there is a physically motivated model it is not obvious which measure is optimal. Thus, in order to be consistent with uses of σ_e for other galaxies and the current incarnation of the BH–σ correlation, one should use σ_e as reported above (324 km s⁻¹), but other correlations should be studied.

We note that contamination of σ_e by the light from stars within BH's kinematic influence is likely to be less important for most other galaxies, since M87 is both close and has an unusually massive BH. At the same time, it may be prudent that this issue be examined for all galaxies in the context of refining the BH–σ relation overall.

Gültekin et al. (2009) present two BH–σ relations, one for all galaxies, and one for elliptical galaxies alone. Using σ_e = 324⁺²⁸₋₁₂ km s⁻¹ gives log(M_BH) = 9.0^+0.4_−0.2 in the first case and log(M_BH) = 9.1^+0.4_−0.2 in the second. Likewise, evaluating the Gültekin et al. BH–L relation with M_V = −22.71 (Lauer et al. 2007) gives log(M_BH) = 9.0 ± 0.2. Both relations thus give M_BH = 1 × 10⁹ M_☉, in contrast to our determination of M_BH = 6.6 × 10⁹ M_☉. Thus, our measurement differs from the predictions of this BH–σ and BH–L relations by 0.82 dex. However, the intrinsic scatter in these relations is estimated by Gültekin et al. to be 0.44 (BH–σ, all galaxies), 0.31 (BH–σ, ellipticals only), and 0.38 (BH–L, ellipticals only). Adding this scatter in quadrature gives estimates of log(M_BH = 9.0^+0.6_−0.5, 9.1^+0.5_−0.4, 9.0 ± 0.4), respectively. Thus, our measured value of log(M_BH = 9.82 ± 0.03) differs from the predictions by 1.4–2σ. Given the present uncertain state of knowledge of the high-mass ends of both the BH–σ and BH–L relations, we do remark on the larger significance of M87 deviation from the relations, except to say that it does highlight the need to improve the sample of BH mass determinations from the most massive galaxies.

Our best-fit BH mass is (6.6 ± 0.4) × 10⁹ M_☉. Sargent et al. (1978) report a BH mass of 6 × 10⁹ M_☉ (after scaling to our assumed distance), which is within 1σ of our reported value. Their model is based on lower spatial resolution data (about 15), assumes that the velocity distribution is isotropic, and does not include a dark halo. It is impressive that after three decades of improvement in data quality, modeling, and understanding, there is essentially no change in the measured BH mass. Part of the reason for the robustness of the Sargent et al. result is that the radial influence on the projected kinematics from the BH extends to nearly 10'' (see Figure 5), so the influence of the BH was clearly visible in their kinematic data. They also use isotropic models, whereas we run axisymmetric models with no restrictions on the anisotropy. To study the effect of the assumption of isotropy, we fit isotropic models to the kinematic data presented in this paper. The comparison between the projected dispersions of the isotropic models and the data is poor, with an increase in χ² by over a factor of two. The poor fit makes it difficult to assign a best-fit mass and the range of equally poor-fitting models has BH masses that range from 6 × 10⁹ to 8 × 10⁹ M_☉, consistent with the models of Section 3, which show significant tangential anisotropy (Figure 7). Thus, in M87, the assumption of isotropy does not have a significant effect on the measurement of the BH mass, although isotropic models provide a poor fit to the data. Sargent et al. also do not include a dark halo, which has been shown to cause the BH to be underestimated. Their velocity dispersions at large radii are lower than ours (245 compared to 300 km s⁻¹), which is most likely because their template library was incomplete and their spectra had lower S/N. The lower dispersion causes the assumed M/L of the stars to be lower, an error of the opposite sign to the error caused by neglect of the dark halo. Thus, the impressive agreement between our value and that of Sargent et al. (1978) appears to be due in part to the competing effects of observational errors (dispersions too small, which makes the stellar M/L too low and the BH mass too large) and oversimplified models (no dark halo or velocity anisotropy, both of which make the BH mass too small). Another often-quoted BH mass determination from stellar kinematics comes from Magorrian et al. (1998) who report a value of 4.2 × 10⁹ M_☉ (for our distance). The likely reason for the difference is that they do not include a dark halo and thus overestimate the stellar M/L.

The BH mass reported here is nearly the same as that reported in Gebhardt & Thomas (2009), within 4%. There is very little kinematic data in common between the two studies. The kinematic data in Gebhardt & Thomas come from older long-slit data at a spatial resolution of 10 (van der Marel 1994), while in this paper we use two-dimensional coverage at a spatial resolution of 008. We further use ground-based data from Murphy et al. (2011) that have excellent S/N and radial extent. There is some data from SAURON (Emsellem et al. 2004) in common between the two studies, but this provides only 10% of the LOSVDs used in the models. Thus, the dynamical models from the two studies use nearly independent kinematic data sets, and give approximately the same answer.

The uncertainties on the BH mass from these two studies are similar even though the data presented here are superior in many ways: the previous uncertainty is 0.5 × 10⁹, whereas the uncertainty with the NIFS data is 0.4 × 10⁹. In order to keep the BH mass measures independent, the models presented in this paper do not include the SAURON data inside of 25. The similarity in the BH mass uncertainty is then due primarily to the fact that the two sets of data have similar accuracy on the kinematics in the central 25. Combining all NIFS data, the accuracy on the velocity dispersion is 0.2% (1 km s⁻¹). Combining all SAURON data within 25 provides the same accuracy. Thus, as long as one has a reliable PSF and no systematic differences in the kinematic extractions, then it is expected that the uncertainty on the BH mass is similar using either data set. We have run a subset of models including both the NIFS data and all SAURON data; in this case, the uncertainty on the BH mass decreases to 0.25 × 10⁹ (with no change in the best-fit mass). We report and utilize the result using only the NIFS data within 25 in order to (1) provide as independent result as realistically possible and (2) control potential systematic differences in the kinematic extractions. Murphy et al. (2011) find a difference in the velocity dispersion of the SAURON data at large radii compared to their measurements, which they attribute to template issues. While we do not find an offset in the dispersion values in the central region, we desire to maintain the independence. The major difference, however, is that there is no degeneracy with the stellar M/L using the NIFS data, whereas the degeneracy is very strong otherwise. Thus, the systematic uncertainty from the M/L profile is effectively removed with the AO data, making the result on the BH mass and orbital structure much more robust.

For M87, the AO data have removed the systematics due to the M/L profile, but the systematics due to the extraction of the kinematics remain important. These systematics include continuum placement, template mismatch, and removal of the AGN contribution. The first two are general and the latter is specific to M87. Getting any of these controlled to better than 1% of the velocity dispersion will be very difficult.

Corrected to our distance, the BH masses reported from gas kinematics are (2.9 ± 0.8) × 10⁹ M_☉ in Harms et al. (1994) and (3.8 ± 1.1) × 10⁹ M_☉ in Macchetto et al. (1997). As discussed in Gebhardt & Thomas (2009) the mass reported here is in conflict with these by about 2σ. Possible reasons for the differences are discussed in Gebhardt & Thomas, with the most likely reason being uncertainty in the inclination of the gas disk. Macchetto et al. assume a value of 51° based on the gas kinematics. Harms et al. assume a value of 42° based on the imaging of the gas emission. The reported difference provides a measure of the systematic uncertainty in the inclination (i.e., whether the gas kinematics or the gas distribution are more affected by non-gravitational forces). Applying this 9° difference in the inclination changes the Macchetto et al. BH mass from (3.8 ± 1.1) × 10⁹ to (5.4 ± 1.3) × 10⁹ M_☉, which would lead to an insignificant difference of 0.6σ between our result and theirs. Of course, the analysis is more complicated than this simple application since one would need to re-model the gas kinematics with a different inclination. A proper treatment would be to include the gas kinematics with the stellar dynamical models. Our focus in this paper is on the stellar kinematics, and we do not attempt to merge the gas kinematic analysis.

While the kinematics obtained from the AO study produce effectively the same BH mass and its uncertainty from kinematics taken in native image quality, the robustness of the measures is greatly strengthened. For example, the BH mass is not dependent on the assumption of constant M/L. Trying to generalize this result to other galaxies with BH mass determinations is difficult since the measure of the BH mass depends on many aspects. There are two observational extremes that we highlight as examples. The first is having a measure of BH mass that comes from observations that resolve well the kinematic influence of the BH. In the most extreme case, high S/N spectra could potentially see the high-velocity wings in the LOSVD due to the BH (as discussed in van der Marel 1994). The second example would be to allow the poorer resolution of the BH, but provide a very accurate measure of the mass-to-light profile. In this paper, we rely on the first strategy; Gebhardt & Thomas rely on the second. That the two strategies give consistent results, at least for M87, suggests that both may be reliable.

Other studies have reported robust measures of the BH mass from ground-based studies that only poorly resolve the BH's kinematic influence. Shapiro et al. (2006) measure a BH for NGC 3379 from SAURON data that is consistent with that measured from HST data using both stars (Gebhardt et al. 2000a) and gas kinematics. Kormendy (2004) summarizes the history of BH mass measures for many galaxies and finds that, in general, the differences are within the reported uncertainties. If one has sufficient signal to noise and two-dimensional coverage (e.g., SAURON or VIRUS-P), then it should be possible to measure a BH mass robustly. Thus, it is not necessarily required to resolve the region influenced by the BH.

Being able to use data that do not resolve well the BH's influence on the kinematics allows us to study BHs that are either distant or of low mass. Both of these regimes are important for understanding the physical nature of the BH correlations with the host galaxy. For example, McConnell et al. (2011) measure a BH mass in NGC 6086, which is 133 Mpc distant. The kinematic influence of the BH is barely resolved, and the degeneracy between the BH mass and M/L profile is strong. However, as demonstrated for M87, as long as one properly characterizes the mass profile at large radii, high signal-to-noise data can measure the BH mass accurately.

It is possible that systematic uncertainties bias the current crop of BH correlations. One obvious consequence could be that without accounting for the effect of systematic uncertainties, the measured intrinsic scatter would increase. Gültekin et al. (2009) measure the scatter of 0.44 dex for the full sample of galaxies with measured BH masses and 0.31 dex for ellipticals. Once systematic effects are understood and included, the intrinsic scatter may decrease. Other consequences include increasing the mass density of BHs, if BH masses are all underestimated, and changing the slope or curvature of any correlation. Schulze & Gebhardt (2011) re-analyze the set of 12 galaxies from Gebhardt et al. (2003) including a dark halo. They find an increase of 50% in the BH mass, due primarily to improved dynamical modeling (more complete orbit sampling) and partly to including a dark halo. The increase correlates with BH mass. It is important to re-evaluate all BH mass estimates.

The key to understanding all of these effects comes from high spatial resolution data. Data from HST (mainly from Space Telescope Imaging Spectrograph, STIS) are generally regarded as providing the most significant results for BH mass studies. The small and stable PSF is a central aspect for the robustness of the data from HST. Future uses of STIS will play an important role for quantifying BH masses. The main obstacle for HST though is that it is a relatively small mirror and requires substantial observing time. For example, in order to measure the BH mass in M87 at the same accuracy presented here would require nearly 100 orbits. While this amount of time could be justified for a small number of objects, going to a much larger sample using HST is difficult. Fortunately, AO observations are in a mature stage where they can provide much larger samples.