THE M87 BLACK HOLE MASS FROM GAS-DYNAMICAL MODELS OF SPACE TELESCOPE IMAGING SPECTROGRAPH OBSERVATIONS^*

Supermassive black holes are thought to play a fundamental role in the growth and evolution of galaxies (Silk & Rees 1998; Di Matteo et al. 2005). This idea is supported by the empirical relationships between the mass of the black hole, whose gravitational influence is limited to a very small region at the center of the galaxy, and the large-scale properties of the host galaxy, such as the bulge stellar velocity dispersion (σ_⋆; e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000a; Tremaine et al. 2002; Gültekin et al. 2009) and the luminosity (L_bul; e.g., Dressler 1989; Kormendy & Richstone 1995; Marconi & Hunt 2003; Sani et al. 2011). Other interpretations have also been put forth to explain the black hole mass–bulge correlations. For example, the relations may simply be the result of random mergers without the need for additional physical processes like active galactic nucleus (AGN) feedback (Peng 2007; Jahnke & Macciò 2011). Our understanding of the interplay between black holes and their host galaxies is highly incomplete. Thus, accurately characterizing the shape and cosmic scatter of the relations is crucial for distinguishing between the various theoretical interpretations, as well as for determining the black hole mass function and space density of supermassive black holes.

The black hole mass–host galaxy correlations are composed of 87 black hole mass (M_BH) measurements, which are most often made by modeling the motions of stars or nuclear gas disks (Kormendy & Ho 2013). In spite of the impressive number of observations, the local black hole mass census is still incomplete, notably for low- and high-mass black holes, and there remain major open questions. At the high-mass end, the slope, intrinsic scatter, and possibly the functional form of the M_BH–host galaxy relations are not well constrained (Gültekin et al. 2009; Graham & Scott 2013; McConnell & Ma 2013). Properly quantifying the scatter at the high-mass end of the M_BH correlations may be especially helpful as non-causal interpretations of the scaling relations predict a decrease in the scatter at higher black hole and galaxy masses (Jahnke & Macciò 2011). Furthermore, the current M_BH–σ_⋆ and M_BH–L_bul relationships make strongly divergent predictions for the black hole masses in the most luminous and the highest-dispersion galaxies (Lauer et al. 2007; Bernardi et al. 2007), leading to questions about which of the two correlations is more fundamental. At least one of the scaling relations must be wrong at the high-mass end, but currently it is not known which, or how to reconcile the difference. In addition, there have been very few consistency checks between the stellar and gas-dynamical methods, therefore it is unclear how much of the intrinsic scatter is due to inconsistencies between stellar and gas-dynamical measurements, or if there a systematic difference between the masses derived using the two methods.

Only direct comparisons between the stellar and gas-dynamical methods within the same galaxy can address some of these questions, but currently such tests have only been attempted on a very small number of galaxies. Stellar-dynamical black hole masses have been made for both IC 1459 and NGC 3379, but the ionized gas turned out to be highly disturbed and a useful gas-dynamical M_BH measurement could not be obtained (Verdoes Kleijn et al. 2000; Gebhardt et al. 2000b; Cappellari et al. 2002; Shapiro et al. 2006; van den Bosch & de Zeeuw 2010). In NGC 4151, the gas-dynamical determination (Hicks & Malkan 2008) is generally consistent with the stellar-dynamical mass, but the authors labeled their stellar-dynamical result as tentative because they were unable to find single best fitting model (Onken et al. 2007). Also, the gas-dynamical measurement for the black hole in M81 (Devereux et al. 2003) agrees with a previous, but preliminary, stellar-dynamical determination (Bower et al. 2000). This leaves six galaxies (NGC 3227, NGC 3998, NGC 4258, NGC 4335, M87, Cen A) for which the two types of black hole mass measurements can be meaningfully compared. The stellar and gas M_BH measurements for NGC 3227, NGC 4258, and Cen A are consistent (Davies et al. 2006; Pastorini et al. 2007; Neumayer et al. 2007; Hicks & Malkan 2008; Siopis et al. 2009; Cappellari et al. 2009), while the stellar-dynamical mass exceeds the gas-dynamical determination by a factor of ∼2–5 for the remaining three galaxies (Macchetto et al. 1997; Verdoes Kleijn et al. 2002; de Francesco et al. 2006; Gebhardt et al. 2011; Walsh et al. 2012).

In particular, the supermassive black hole in the giant elliptical galaxy M87 has been the subject of numerous mass determinations, beginning with the pioneering work of Sargent et al. (1978) and Young et al. (1978), who suggested the presence of a ∼(3–5) × 10⁹ M_☉ black hole based on the surface brightness and stellar velocity dispersion profiles. There were concerns associated with these and other early stellar-dynamical black hole mass measurements (see the review by Kormendy & Richstone 1995, and references therein), but a much improved stellar-dynamical measurement was recently obtained by Gebhardt & Thomas (2009) and Gebhardt et al. (2011) using three-integral, axisymmetric Schwarzschild models that also include the mass distribution of the dark halo. More specifically, Gebhardt et al. (2011) applied orbit-based models to high angular resolution spectroscopy obtained with Gemini's Near-Infrared Integral Field Spectrometer assisted by adaptive optics and additional large-scale data (Emsellem et al. 2004; Murphy et al. 2011), finding M_BH = (6.6 ± 0.4) × 10⁹ M_☉.

In addition to the stellar-dynamical M_BH measurements, there are also two gas-dynamical measurements for the black hole in M87. In fact, M87 was the first gas-dynamical target with the Hubble Space Telescope (HST). Using spectra acquired with the Faint Object Spectrograph (FOS), Harms et al. (1994) modeled the emission-line kinematics as a rotating disk to infer a central mass of (2.4 ± 0.7) × 10⁹ M_☉. A few years later, Macchetto et al. (1997) used spectroscopy from the Faint Object Camera (FOC) and rotating, thin-disk models that, for the first time, incorporated the important effects of the propagation of light through the telescope and spectrograph optics to determine M_BH = (3.2 ± 0.9) × 10⁹ M_☉.

The work of Harms et al. (1994) and Macchetto et al. (1997) were important milestones for HST and the field of supermassive black hole detection. However, neither study was able to fully map out the kinematic structure of the gas disk. Harms et al. (1994) relied on just six spectra obtained with the FOS and Macchetto et al. (1997) used FOC long-slit observations at three, non-contiguous, parallel positions that only covered a fraction of the disk. Also, the disk inclination has been a source of uncertainty in the M87 gas-dynamical black hole mass. Harms et al. (1994) adopted a value of i = 42° ± 5°, which was measured by Ford et al. (1994) from HST Hα narrowband images, while Macchetto et al. (1997) allowed the inclination to be a free parameter in their dynamical models, but found a wide range of acceptable values with 47° < i < 65°. Perhaps more importantly, later investigations of circumnuclear gas disks revealed the prevalence of an intrinsic velocity dispersion (e.g., van der Marel & van den Bosch 1998; Verdoes Kleijn et al. 2000; Dalla Bontà et al. 2009). Although the origin of the intrinsic velocity dispersion is not known, it may be dynamically important, contributing pressure support to the disk. Models that treat such systems as simple, thin, dynamically cold disks will underestimate M_BH (Barth et al. 2001; Neumayer et al. 2007; Walsh et al. 2010). Only by carrying out detailed modeling of the emission-line profiles, including the contributions to the line widths from rotational and instrumental broadening, will it be possible to determine if there is significant intrinsic velocity dispersion in the M87 disk and its subsequent effect on the gas-dynamical M_BH measurement.

Given that M87 is located at the extreme high-mass end of the M_BH–host galaxy relations and is a very important test case for a gas/stars comparison study, the factor of ∼2 discrepancy between the best current stellar and gas-dynamical masses is troubling. Furthermore, the mass of the black hole in M87 is a crucial parameter for making inferences about the size of the innermost stable circular orbit of the black hole and the black spin from very long baseline interferometry observations (Doeleman et al. 2012). In this paper, we revisit the gas-dynamical black hole mass measurement. We use new HST data acquired with the Space Telescope Imaging Spectrograph (STIS) and carry out more comprehensive gas-dynamical modeling than has previously been attempted for this galaxy. We describe the observations in Section 2, the emission-line kinematics in Section 3, the modeling procedures in Sections 4 and 5, the results and error analysis in Section 6, and conclude in Section 7. Throughout the paper, we assume a distance to M87 of D = 17.9 Mpc in order to allow for a direct comparison to the stellar-dynamical M_BH measurement by Gebhardt et al. (2011).

We obtained new HST STIS observations of M87 on 2011 June 7. The observations have a wavelength scale of 0.554 Å pixel⁻¹ and a spatial scale of 00507 pixel⁻¹. We used the G750M grating, centered on 6581 Å, to provide coverage of the Hα spectral region, which includes the [O i] λλ6300, 6364, Hα, [N ii] λλ6548, 6583, and [S ii] λλ6716, 6731 emission lines. However, the [O i] and [S ii] lines were weak, except at locations very close to the nucleus, and we focus solely on the Hα and [N ii] lines for the remainder of the paper. We used the E1 aperture position to place the nucleus near the readout end of the CCD, to avoid charge transfer efficiency losses. The STIS 52x0.1 slit was placed at five parallel positions without any space between adjacent positions, and was oriented at a position angle (P.A.) of 51° (within 15° of the nuclear gas disk's major axis as determined by Macchetto et al. 1997). Between 3 and 5 dithered exposures were acquired at each slit position to aid in the removal of cosmic rays, leading to exposure times that ranged from 1521 to 2911 s per slit position.

The data from each subexposure were reduced by trimming the overscan region, subtracting the bias and dark files, and applying flat-field corrections. After this initial processing, we used LA-COSMIC (van Dokkum 2001) to remove cosmic rays and hot pixels. We completed the reduction of each subexposure by wavelength and flux calibrating the data and rectifying for geometric distortions. The reductions were carried out using individual IRAF⁵ tasks within the standard Space Telescope Science Institute (STScI) pipeline, with the exception of the cosmic-ray cleaning step. The reduced, geometrically rectified, subexposures were aligned and combined using the IRAF tasks IMSHIFT and IMCOMBINE in order to produce the final two-dimensional (2D) STIS image at each of the five slit positions.

Furthermore, from the HST archive we retrieved Wide Field Planetary Camera 2 (WFPC2) F658N, F547M, and F702W images, which were originally made under programs GO-5122 and GO-5476 and had the M87 nucleus placed on the PC detector. The images in each filter were taken as a sequence of two individual exposures with total exposure times of 2700 s, 800 s, and 280 s for the F658N, F547M, and F702W images, respectively. The subexposures were combined and the cosmic rays were rejected using the IMCOMBINE and LA-COSMIC IRAF tasks. We then created a continuum-subtracted Hα+[N ii] image by rotating the F658N and F702W images to a common orientation, and subtracting a scaled version of the F702W image from the F658N image. The scale factor was chosen such that a background region in the continuum-subtracted image would have a mean flux as close to zero as possible. In Figure 1, we present the emission-line image with the location of the STIS slits overlaid, as well as the F547M WFPC2 image. A compact central source and an extended disk structure with several spiral-like wisps is seen in the continuum-subtracted image, as was previously found by Ford et al. (1994).

Zoom In Zoom Out Reset image size

From each of the 2D STIS images, we extracted spectra from individual rows out to ∼05 from the slit center, which is as far as the emission lines were detectable. At locations where the emission lines became weak, we binned together three STIS rows in order to improve the signal-to-noise ratio (S/N). Prior to measuring the emission-line kinematics, we subtracted the starlight continuum on a row-by-row basis by fitting a line between wavelengths 6295–6862 Å excluding the regions near the emission lines, then subtracting this fit from the spectrum. Due to the low S/N in the continuum and the small wavelength range, a straight line is a sufficient description of the continuum near Hα and we do not attempt a more complicated subtraction. Below we detail our spectral fitting procedure and describe the observed velocity fields.

3.1. Spectral Fitting Method

We simultaneously fit three single Gaussians to the continuum-subtracted Hα and [N ii] emission lines using the MPFIT library in IDL (Markwardt 2009), which employs a Levenberg–Marquardt least-squares minimization technique. We then performed a Monte Carlo, adding random Gaussian noise to the spectrum based upon the model residuals, and refit the new spectrum to determine the velocity, velocity dispersion, and flux for each of the three emission lines during each iteration. From the resulting distributions, we were able to determine whether the [N ii] λ6548 and λ6583 velocities were consistent within 2σ. If the velocities were consistent, we calculated the mean and the 2σ errors from the combined distribution of the [N ii] velocities, which was taken to be the final velocity and error for that spectral row. The same process was adopted for the [N ii] velocity dispersions and for 3 × the flux of the [N ii] λ6548 line and the flux of the [N ii] λ6583 line. If instead the kinematics of the [N ii] lines were not consistent, we fit progressively more constrained models (such as tying together the [N ii] velocity widths) until the [N ii] kinematics were consistent. With our approach, we are assuming that the [N ii] velocities should be equal, the [N ii] velocity dispersions should be the same, the fluxes of the [N ii] lines should maintain a 3:1 ratio (as dictated by atomic physics), and that any differences that are seen are due to fitting error. We chose not to fit Gaussian models that enforce these constraints from the start because we found doing so produced very small errors on the kinematics—typically we measured velocity errors of only a few km s⁻¹ for the good S/N spectra (while the STIS pixel scale is ∼25 km s⁻¹).

We also tested whether the spectral fits could be significantly improved with the addition of a broad Hα component described by a single Gaussian. Including a broad component either caused the χ² per degree of freedom ($\chi ^2_\nu$) to decrease by less than 15%, produced a broad component with zero flux, or resulted in unreasonably small fluxes and/or line widths for the Hα and [N ii] narrow components. We cannot fully rule out the existence of a compact broad-line region, but the data do not require the presence of a physically distinct broad-line component, similar to the giant elliptical galaxy M84 (Walsh et al. 2010). In both these galaxies, the very broad central line widths are consistent with coming from the inner portion of the rotating disk. Thus, the kinematics presented here come from single Gaussian fits to the Hα and [N ii] lines only without the inclusion of a model for a broad component.

While the above fitting method worked well for many of the observed spectra, the spectra extracted near the center of slit positions 2–4 were complex with severely blended emission lines. In addition, some of these rows showed artifacts, in the form of a periodic fluctuation in the flux, from the geometric rectification step of the data reduction (Dressel et al. 2007). Moreover, spectra extracted near the center of slit position 5 exhibited blended [N ii] λ6548 and Hα emission lines and a noisy [N ii] λ6583 line that may be double-peaked. It was therefore difficult to decompose the spectra into Hα and [N ii] components with accurate and unique mean velocities, velocity dispersions, and fluxes, even when applying constraints to the three single Gaussians that were being fit. Consequently, we consider the measurements from the central rows of slits 2–5 to be unreliable and we do not use them to constrain our gas-dynamical models. The innermost kinematic measurements provide relatively little useful information when constraining the disk models because the line profiles near the nucleus are dominated by rotational broadening and blurring by the point-spread function (PSF). For completeness, however, we incorporate into our final error budget the effect on the black hole mass when the measurements deemed as unreliable are used to constrain the dynamical model as well. In Figure 2, we show example fits to spectra extracted from slit positions 3 and 5, including rows where the kinematic measurements are considered unreliable.

Zoom In Zoom Out Reset image size

3.2. Observed Velocity Fields

The observed velocity fields are modeled assuming that the gas participates in circular rotation in a thin disk. We provide a brief overview of the modeling here, but refer the reader to Barth et al. (2001) and Walsh et al. (2010) for further details. At each radius (r) in the disk, the circular velocity is calculated relative to the galaxy's systemic velocity ($\upsilon _\mathrm{sys}$) based on the enclosed mass, which depends on the black hole mass, the stellar mass profile, and the mass-to-light ratio (ϒ). The model velocity field is created on a pixel grid that is oversampled relative to the size of a STIS pixel, and is projected onto the plane of the sky for a given value of the disk inclination angle (i).

The intrinsic line-of-sight velocity profiles are assumed to be Gaussian before passing through the telescope optics and are generated on a velocity grid with a bin size that matches the STIS scale. The intrinsic line profiles are centered on the projected line-of-sight velocity at each point on the model grid and weighted by the emission-line flux distribution. In addition, the widths of the intrinsic line profiles are set based upon contributions from the thermal velocity dispersion of the gas and the velocity dispersion due to turbulent motion internal to the gas disk.

Next, the model velocity field is "observed" in a manner that matches the STIS observations. This synthetic observation includes accounting for the blurring by the telescope PSF by convolving a Tiny Tim (Krist & Hook 2004) PSF model with each velocity slice of the line profile grid. Also, we propagate the line profiles through the finite width of the STIS slit by applying shifts to the velocity depending on the location along the slit width that the photon enters (Maciejewski & Binney 2001). After rebinning to the STIS pixel size, we are left with a model 2D spectral image similar to the observed STIS data. At this stage, we convolve the model 2D STIS image with the Tiny Tim CCD charge diffusion kernel in order to account for the charge that is spread between non-subsampled, neighboring pixels.

Finally, we extract spectra on a row-by-row basis from the model 2D image and fit a single Gaussian to the emission line, analogous to the measurements that were made from the data. We thus are able to measure the model velocity, velocity dispersion, and flux as a function of position along the slit. The best model is taken to be the one that most closely matches the observations in the χ² sense. We use the downhill simplex algorithm by Press et al. (1992) to minimize χ² and determine the best-fit parameters. We separately optimize the fits to the observed line widths and fluxes before determining the black hole mass from the model fit to the observed radial velocities. When fitting to the observed velocities, the model parameters are M_BH, ϒ, i, $\upsilon _\mathrm{sys}$, θ, x_offset, and y_offset. Here, θ is the angle between the slit and the projected major axis of the disk, and x_offset and y_offset are the offsets between the center slit and the nucleus.

Gas-dynamical models for M87 were calculated on a pixel grid oversampled by a factor of s = 6 relative to the STIS pixel size. We found that preliminary models with a subsampling factor of s = 6 were able to best reproduce the observed velocities and emission-line fluxes, and could be run in a reasonable amount of time. Thus, we selected s = 6 for the final models. We also incorporated the effect on M_BH of our subsampling factor choice into the error budget, which will be summarized in Section 6.1. In addition, we used a PSF generated with Tiny Tim for a monochromic filter passband at 6600 Å. The model is an s = 6, 04×04 portion of the full STIS PSF. Such a model excludes part of the extended PSF wings, with 12% of the total flux being omitted. As will be shown in Section 6.1, the black hole mass is not strongly dependent on the size of the PSF model, but because the computation time increases rapidly with increasing PSF size, we chose to use a 04-diameter PSF. Finally, we fit disk models over the entire radial range covered by the observed kinematics. With a fitting radius of 05, not only is there an adequate number of measurements to constrain the dynamical models, but also the disk model continues to provide a good description of the observations out to this radius. As discussed previously, we exclude the uncertain kinematics measured from the central rows during the modeling.

Below we provide additional details about the M87 gas-dynamical models, including the stellar contribution to the circular velocities, the parameterization of the intrinsic emission-line flux distribution, the components that make up the width of the line profiles, and our approach for allowing the turbulent velocity dispersion to be dynamically important.

5.1. Stellar Mass Model

We determined the stellar contribution to the gravitational potential using the V-band surface brightness profile of M87 from Kormendy et al. (2009). The profile combines measurements from ground-based images with high-resolution HST WFPC1/PC data, which were originally presented in Lauer et al. (1992) and have already been deconvolved by the telescope PSF. In addition to the surface brightness profile, Kormendy et al. (2009) present the ellipticity () of the isophotes as a function of major axis radius. At radii within 180'', the isophotes are nearly circular, with < 0.2. Therefore, we deprojected the surface brightness profile assuming that the intrinsic density distribution of the central regions could be modeled as the sum of spherically symmetric components with Gaussian profiles. A similar approach has been used to derive the luminous density distribution for past gas-dynamical black hole mass measurements (e.g., Sarzi et al. 2001; Barth et al. 2001; Coccato et al. 2006; Walsh et al. 2010). To carry out the deprojection, we used the Multi-Gaussian Expansion (MGE) software developed by Cappellari (2002), which is based on the methods described in Emsellem et al. (1994). In addition, the M87 nucleus is classified as a Type 2 Low Ionization Nuclear Emission-line Region (Ho et al. 1997). We attribute the innermost compact Gaussian component of the MGE model, with a dispersion of 3.5 pc, to light from the AGN, and exclude this component from the stellar mass distribution.

The final product of the multi-Gaussian deprojection of the surface brightness profile is the stellar contribution to the circular velocity in the galaxy, which is shown in Figure 3. The figure also provides the circular velocities due to the presence of a 3.0 × 10⁹ M_☉ black hole, which is a plausible lower limit on the black hole mass from our analysis. Clearly, the black hole's sphere of influence is well resolved by our STIS observations, and is further supported by the nearly Keplerian shape of the multi-slit velocity curves. In other words, the black hole dominates the galaxy's gravitational potential and the stellar contribution is negligible over the radial range probed by the STIS kinematics. Consequently, we found that our dynamical models are insensitive to the stellar mass-to-light ratio, and we chose to fix ϒ = 4 (V-band, solar units) during the model fitting. This value was selected using the B−V color for the galaxy from the Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al. 1991) and the relations presented in Bell et al. (2003). Such a stellar mass-to-light ratio is typically found at the centers of giant elliptical galaxies (Cappellari et al. 2006). We also tested the effect on M_BH when another fixed value for ϒ is adopted, as well as when ϒ is allowed to be a free parameter. The results will be described in Section 6.1.

Zoom In Zoom Out Reset image size

5.2. Emission-line Flux

5.3. Line Widths

5.4. Asymmetric Drift Correction

Currently, there is no clear consensus for how to interpret the physical nature of the intrinsic velocity dispersion, hence we explore two scenarios that should bracket the range of possible masses for the black hole in M87. In one case, we assume that the intrinsic velocity dispersion does not affect the circular velocity, which could result if the intrinsic velocity dispersion is due to local microturbulence but the bulk motion of the gas remains in circular motion (van der Marel & van den Bosch 1998). In the second case, we assume that the internal velocity dispersion is due to turbulent motion that supplies pressure support to the disk, causing the observed mean rotational speed ($\upsilon _\phi$) to be smaller than the local circular velocity ($\upsilon _c$) for a given black hole mass (e.g., Coccato et al. 2006; Neumayer et al. 2007; Walsh et al. 2010). For M87, we estimate an upper bound to the black hole mass by applying an asymmetric drift correction following Barth et al. (2001). Under the assumption that the gas motions are isotropic in the r and z (cylindrical) coordinates, the correction is

$$\begin{equation} \upsilon _c^2\ -\ \upsilon _\phi ^2\ =\ \sigma _r^2 \left [ -\frac{d \ln \nu }{d \ln r}\ -\ \frac{d \ln \sigma _r^2}{d \ln r}\ -\ \left (1\ -\ \frac{\sigma _\phi ^2}{\sigma _r^2} \right) \right ] \;. \end{equation} \tag{ 1 }$$

We take the intrinsic emission-line flux distribution to be a proxy for the number density of tracer particles, or clouds in the gas disk, ν. Since the asymmetric drift correction depends on an azimuthally averaged value for the radial derivative of ν, we refit the observed fluxes with an analytic function composed of two Gaussians (similar to the flux model presented in Section 5.2), but required that the isophotes have the same center, P.A., and axis ratio. For a given intrinsic radial velocity dispersion, σ_r, we computed σ_ϕ and the projected velocity dispersion as outlined by Barth et al. (2001), before adding the projected velocity dispersion in quadrature to σ_th and propagating the line profiles through the telescope and instrument. Although we initially attempted to parameterize σ_r as a constant+exponential function, we found that the intrinsic radial velocity dispersion was best described as constant. By calculating disk models with an asymmetric drift correction, and fitting to the observed line widths, we found a best-fit value of σ_r = 191 km s⁻¹.

After optimizing the model fits to the observed emission-line fluxes and velocity dispersions, we fit disk models to the observed radial velocities. The results returned by the χ²-minimization routine seemed to strongly depend on the initial guesses for M_BH and i in particular. Therefore, to ensure a global minimum was found, we used a two-step process. First, we ran a coarse grid of models with fixed M_BH and i values while the other parameters (θ, $\upsilon _\mathrm{sys}$, x_offset, y_offset) were allowed to vary. Second, we used the parameter values of the top five models (with the lowest χ²) from the grid search as the starting points to do five runs in which all the parameters were allowed to vary. By doing these five runs, we were able to distinguish between local and global minima.

The results of the best-fit model without an asymmetric drift correction are summarized in Table 2 and the model predictions are shown alongside the data in Figures 4 and 5. We measure a black hole mass of M_BH = 3.5 × 10⁹ M_☉. The other model parameters, including the disk inclination angle, are well constrained. During the numerous preliminary model runs, the parameters varied at most by 7°, 9 km s⁻¹, 2°, 001, and 004 for i, $\upsilon _\mathrm{sys}$, θ, x_offset, and y_offset, respectively. The model has χ² = 66, found by comparing the model predictions to the observed velocities. Given that 62 velocity measurements were used to constrain six model parameters, the χ² per degree of freedom is $\chi ^2_\nu = 1.2$. As can be seen in Figure 4, the model is able to reproduce the shape of the observed velocity curves and the emission-line flux profiles very well. At some locations in the disk, an intrinsic velocity dispersion is needed to match the data, although this velocity dispersion is relatively small given the large velocities found in M87. Moreover, our model for the intrinsic velocity dispersion still severely underestimates the velocity dispersions located −010 and −015 away from the center of slits 3 and 4. We will examine this discrepancy and the possible impact on the inferred black hole mass in Section 6.1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Disk Model Parameters

MBH	i	$\upsilon _\mathrm{sys}$	θ	xoffset	yoffset
(M☉)	(°)	(km s−1)	(°)	('')	('')
$(3.5_{-0.7}^{+0.9}) \times 10^9$	42$_{-7}^{+5}$	1335$_{-9}^{+3}$	6 ± 2	−0.02 ± 0.01	$0.03_{-0.04}^{+0.02}$

Notes. These results give the best-fitting parameters for models calculated without an asymmetric drift correction. The errors given for the black hole mass are 1σ uncertainties, and the errors provided for the other parameters correspond to the range of values found from all the preliminary models that were run. During the modeling, ϒ was frozen at 4 (V-band solar units). A relative angle of θ = 6° corresponds to a disk major axis position angle 45° east of north.

Download table as:  ASCIITypeset image

6.1. Uncertainty in the Black Hole Mass

We measured the statistical errors on M_BH by holding M_BH fixed while allowing the other parameters to vary such that χ² was minimized. During this process, we froze M_BH to values between (1.0–9.7) × 10⁹ M_☉ in steps of 0.1 × 10⁹ M_☉. We searched for the range of M_BH values that caused the minimum χ² to increase by 1.0 and 9.0 to estimate the 1σ and 3σ model fitting uncertainties. Figure 6 presents the results of these disk models, specifically the χ² variation with M_BH after marginalizing over the other parameters. Ultimately, we determine that the range of black hole masses is (3.4–4.0) × 10⁹ M_☉ (1σ uncertainties) and (3.0–6.8) × 10⁹ M_☉ (3σ uncertainties).

Zoom In Zoom Out Reset image size

These M_BH errors account for the random noise in the gas dynamical models and the uncertainties associated with the parameters i, $\upsilon _\mathrm{sys}$, θ, x_offset, and y_offset. However, there are additional sources of uncertainty that are not incorporated into the statistical errors. Below, we explore these sources of uncertainty and their impact on M_BH, in order to arrive at the final range of possible black hole masses for M87.

Subsampling factor. We ran gas-dynamical models on pixel grids that were subsampled by a factor of s = 2–10, in intervals of s = 2. We found that the black hole mass varied between (3.5–3.9) × 10⁹ M_☉ with a root-mean-square (rms) scatter of 2.3 × 10⁸ M_☉, which is 7% of the best-fit black hole mass.

PSF size. When calculating dynamical models using a Tiny Tim PSF with diameters of 03, 04, 05, and 06, the rms scatter in the black hole mass was 1.2 × 10⁸ M_☉. The PSF size has a minimal impact on the inferred black hole mass, as the scatter is just 3% of the best-fit mass.

Stellar mass-to-light ratio. The stellar contribution to the galaxy's gravitational potential is negligible over the region probed by the STIS data (∼40 pc). As a result, our models were unable to constrain the stellar mass-to-light ratio, and preliminary model runs returned a wide range of values for ϒ. We chose to fix ϒ = 4 (V-band, solar) during the final model run, but also tested freezing ϒ to 9, which is the stellar mass-to-light ratio found for M87 by Murphy et al. (2011) and adopted by Gebhardt et al. (2011), and allowing the parameter to float. The test produced in an rms scatter in the black hole mass of 1.3 × 10⁸ M_☉, or 4% of the best-fit M_BH.

Emission-line flux model. The intrinsic emission-line flux distribution was modeled as the sum of two Gaussian components, whose isophotes had different centers, P.A.s, and axis ratios. An additional four analytic forms of varying complexity were also able to reproduce the observed emission-line fluxes: two exponentials, one Gaussian+one exponential, three Gaussians, three exponentials, and two Gaussians+one exponential. These analytic functions also represented elliptical isophotes with different centers, P.A.s, and axis ratios. Although the quality of the fit to the observed velocities varied significantly depending on the intrinsic flux distribution used (χ² ranged from 66 to 130), there was a small effect on the black hole mass. The disk models with different parameterizations of the emission-line flux distribution returned best-fit masses between 3.1 × 10⁹ M_☉ and 4.2 × 10⁹ M_☉, with an rms scatter of 5.0 × 10⁸ M_☉, or 7% of the best-fit black hole mass.

Differences between the [N ii] and Hα kinematics. Gas-dynamical models were run using the [N ii] kinematics because [N ii] λ6583 was the strongest emission line in our spectra. If gas-dynamical models are instead fit to the Hα kinematics, the black hole mass increases by 11% to 3.9 × 10⁹ M_☉.

Including all kinematic measurements. We had difficulty fitting the spectra extracted near the center of slits 2–5. We consider the kinematics and emission-line fluxes measured from these rows unreliable and did not use them to constrain our final disk model. When these velocities, velocity dispersions, and emission-line fluxes are instead included in the fit, the best-fit black hole mass decreased by 11% to 3.1 × 10⁹ M_☉.

Asymmetric drift correction. By assigning a dynamical origin to the intrinsic velocity dispersion and calculating an asymmetric drift correction as discussed in Section 5.4, the black hole mass increased by 6% to 3.7 × 10⁹ M_☉. Such a minor change to M_BH is not surprising considering the small amount of intrinsic velocity dispersion relative to the large velocities that are seen in M87. Compared to the dynamically cold, thin-disk model, the model with an asymmetric drift correction was a worse fit to the observed velocities, with a χ² = 71.

The predictions from the model with asymmetric drift correction, however, continue to be significantly smaller (by ∼140–240 km s⁻¹) than the observed dispersions located −010 and −015 away from the center of slits 3 and 4. In order to establish the maximum possible increase in the black hole mass, we searched for a new parameterization of the intrinsic velocity dispersion that would match these few points, at the expense of being able to reproduce the observed velocity dispersions at the other locations. Since the distribution of the observed line widths is asymmetric, and the velocity dispersions at negative Y Offsets are generally larger than those measured at positive Y Offsets, we modified the input kinematics by simply replacing the observed velocity dispersions from the positive Y Offset side with the values from the negative Y Offset side. We then fit gas-dynamical models with an asymmetric drift correction first to the modified line widths and then to the velocities. We determined that the intrinsic velocity dispersion was best characterized as a constant+exponential function, with σ₀ = 202 km s⁻¹, σ₁ = 2040 km s⁻¹, and r₀ = 4.4 pc, and that the black hole mass increased by 20% to M_BH = 4.2 × 10⁹ M_☉. As can be seen in Figure 7, the new model is able to sufficiently match the high velocity dispersions observed at Y Offsets of −010 and −015 from slits 3 and 4.

Zoom In Zoom Out Reset image size

The formal model fitting error and the additional sources of uncertainty were added in quadrature to derive the final range of black hole masses. Thus, we find the M87 black hole mass is (2.8–4.4) × 10⁹ M_☉ (1σ uncertainties) and (2.7–6.9) × 10⁹ M_☉ (3σ uncertainties) with a best-fit value of 3.5 × 10⁹ M_☉. We propagated the 6% increase in the black hole mass resulting from the asymmetric drift correction into the upper uncertainty of M_BH only. Even when using an extreme characterization of the intrinsic velocity dispersion to better match the observed line widths at Y Offsets of −010 and −015 in slits 3 and 4, the black hole mass does not increase beyond the final 1σ uncertainties.

Using newly acquired STIS observations of M87 with the 01-wide aperture at five adjacent parallel positions, we have measured the [N ii] velocity, velocity dispersion, and flux well within the black hole's dynamical sphere of influence. Through gas-dynamical modeling of the STIS data, we have determined that the mass of the black hole is $M_\mathrm{BH} = (3.5_{-0.7}^{+0.9}) \times 10^9\, M_\odot$. This black hole mass is similar to the expectations from the M_BH–bulge relations. With σ_⋆ = 324 km s⁻¹ (Gebhardt et al. 2011) and a V-band luminosity calculated from M_V = −22.71 (Lauer et al. 2007), the most recent calibrations of the M_BH–σ_⋆ relation suggest that M_BH ∼ 3 × 10⁹ M_☉ (Graham & Scott 2013; McConnell & Ma 2013) and the M_BH–L_bul correlation predicts M_BH = 2 × 10⁹ M_☉ (McConnell & Ma 2013).

Unlike the past gas-dynamical studies of Harms et al. (1994) and Macchetto et al. (1997), we have mapped out the complete kinematic structure of the emission-line disk. We have constrained the gas disk inclination angle, which has previously been a source of uncertainty, to be between 35° < i < 47° with a best-fit value of i = 42°. In addition, our models included a detailed treatment of the effects of the telescope and instrument. For the first time, we have found that a relatively small amount of additional velocity dispersion internal to the gas disk is needed to reproduced the observed line widths. Macchetto et al. (1997) did not find evidence for such an intrinsic velocity dispersion in the M87 disk, but the instrumental resolution of their FOC observations was low (FWHM ∼ 430 km s⁻¹). Thus, the order-of-magnitude improvement in spectral resolution of STIS compared to the older FOC data has given us the sensitivity to detect an intrinsic velocity dispersion. We considered the possibility that the intrinsic velocity dispersion provides dynamical support to the gas disk, and determined that inferred mass increases by just 6%. We incorporated this effect, along with a number of other sources of uncertainty, into the error budget.

Before comparing our mass measurement to the prior mass determinations, there are some additional systematics to keep in mind when interpreting the results. Our thin-disk model matches the overall shape of the observed velocity curves very well and returns $\chi ^2_\nu = 1.2$; however, the data are not randomly scattered about the best-fit model at all locations. Instead the model systematically deviates from the observations in certain regions, at Y Offset = 025–04 at slit 3, for example. Consequently, there may be some regions that depart from circular rotation, whose small-scale velocity structure cannot be matched with any disk model. This is a common problem for gas-dynamical models in general, and makes estimating the uncertainties in the black hole mass challenging. Past studies have opted to rescale the error bars on the kinematics, such that the best-fit model would have a $\chi ^2_\nu \approx 1$, prior to measuring the statistical errors on M_BH (e.g., Barth et al. 2001; Atkinson et al. 2005; de Francesco et al. 2006; Dalla Bontà et al. 2009). For this work, we do not adjust the kinematic errors because our best-fit model already has $\chi ^2_\nu \approx 1$. Despite the reasonable $\chi ^2_\nu$ of our best-fit model, we note that the model is not a perfect description of the gas kinematics everywhere in the disk.

Furthermore, we treat the (small amount of) intrinsic velocity as dynamically important by applying an asymmetric drift correction. The correction is applicable in the limit of collisionless particles and for σ_r/v_c ≪ 1. Here, we applied the correction to gas clouds, which are not collisionless, and to a system that has σ_r/v_c < 0.3 over the radial extent of the STIS kinematics. Also, the asymmetric drift correction depends upon the number density of clouds in the disk, ν(r), which is unknown. While we assumed that the intrinsic emission-line flux distribution can be used as a proxy for ν(r), this is at best a rough approximation. With these considerations in mind, our model with an asymmetric drift correction only provides estimation, albeit a useful one, of the dynamical influence of the intrinsic velocity dispersion.

Our mass measurement is consistent with the two previous gas-dynamical measurements of the black hole in M87. The Harms et al. (1994) and Macchetto et al. (1997) masses are M_BH = (2.9 ± 0.8) × 10⁹ M_☉ and M_BH = (3.8 ± 1.1) × 10⁹ M_☉, after scaling to our adopted distance. Of the remaining disk model parameters, our best-fit inclination angle is the same as the one used by Harms et al. (1994) and is consistent with the lower bound of acceptable values (i = 47°–65°) found by Macchetto et al. (1997), while our best-fit systemic velocity ($\upsilon _\mathrm{sys} = 1335$ km s⁻¹) is larger than the ones used by Harms et al. (1994; $\upsilon _\mathrm{sys} = 1309$ km s⁻¹) and determined by Macchetto et al. (1997; $\upsilon _\mathrm{sys} = 1290$ km s⁻¹). For comparison, the recession velocities measured from optical lines given by the NASA Extragalactic Database span a wide range of values, from 800 to 1747 km s⁻¹. We established that the projected major axis of the gas disk is 45° east of north, which is roughly perpendicular to the jet and similar to P.A. of 38° measured by Macchetto et al. (1997).

The most recent M87 stellar-dynamical black hole measurement of M_BH = (6.6 ± 0.4) × 10⁹ M_☉ (Gebhardt et al. 2011) is a factor of two larger than our mass measurement, and there is a 2σ discrepancy between our gas-dynamical mass and the Gebhardt et al. (2011) stellar-dynamical mass when accounting for the uncertainties associated with each measurement. Recent stellar-dynamical work has shown the importance in some galaxies of including a dark halo, incorporating orbital libraries that better sample the phase-space occupied by tangential orbits, using triaxial geometries, and accounting for a spatially varying stellar mass-to-light ratio (Gebhardt & Thomas 2009; Shen & Gebhardt 2010; van den Bosch & de Zeeuw 2010; Schulze & Gebhardt 2011; McConnell et al. 2013). While Gebhardt et al. (2011) address the first two potential sources of systematic error, they do not examine the latter two sources. In particular, they do not fit for the galaxy's intrinsic shape and orientation, and instead use axisymmetric orbit-based models that assume an edge-on inclination. Therefore, it is still possible for the stellar-dynamical black hole mass measurement to change again if more general triaxial models are applied, although the magnitude and direction of the effect on the inferred black hole mass is unknown a priori. Also, Gebhardt et al. (2011) construct stellar-dynamical models under the assumption of a constant stellar mass-to-light ratio, however, as shown by McConnell et al. (2013), spatial gradients in the mass-to-light ratio can have a non-negligible effect on M_BH. In contrast, our dynamical models rely on gas in circular orbits at small radii and are insensitive to the overall normalization of the stellar mass-to-light ratio (Section 6.1), making our gas-dynamical measurement immune to this uncertainty.

Our gas-dynamical measurement and the stellar-dynamical determination of Gebhardt et al. (2011) utilize high-quality, high-resolution observations and apply some of the most up-to-date modeling techniques, yet there continues to be a discrepancy in the M87 black hole mass. This highlights the need for carrying out more such cross-checks between the two methods. Currently, M87 is one of six consistency tests that have yielded a meaningful comparison between the gas and stellar-dynamical techniques, and in three cases (this study included) the stellar-dynamical mass exceeds the gas-dynamical determination by a factor of ∼2–5. Since the stellar and gas-dynamical techniques are completely independent methods with different systematic effects, only by applying both techniques to the same object can conclusions be made regarding the consistency of the methods, the subsequent effects on the M_BH scaling relations, and the magnitude and distribution of the intrinsic scatter of the M_BH–host galaxy relations.

J.L.W. has been supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under Award No. 1102845. Support for program 12162 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. Research by A.J.B. has additionally been supported by NSF grant AST-1108835. Some of the data presented in this paper were obtained from the Multimission Archive at the Space Telescope Science Institute (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. We dedicate this paper to the memory of Wallace L. W. Sargent, whose early work on M87 has had a profound impact on the field of black hole searches. A.J.B. and L.C.H. have also benefited enormously from many years of interaction and fruitful scientific collaboration with Wal.