THE BLACK HOLE MASS IN THE BRIGHTEST CLUSTER GALAXY NGC 6086

We use a combination of R-band (0.6 μm) and I-band (0.8 μm) photometry to constrain the stellar mass profile of NGC 6086. For radii out to 10'', we adopt the high-resolution surface brightness profile presented in Laine et al. (2003), obtained with WFPC2 on HST. This surface brightness profile has been corrected for the WFPC2 PSF by applying the Lucy–Richardson deconvolution method (Richardson 1972; Lucy 1974); specific details of the implementation are described in Laine et al. (2003).

At larger radii out to 86'', we use R-band data from T. R. Lauer, M. Postman, & M. A. Strauss (1995, private communication), obtained with the 2.1 m telescope at Kitt Peak National Observatory (KPNO). The KPNO data have a field of view (FOV) of 52 × 52, which enables accurate sky subtraction. To create a single surface brightness profile, we assessed the individual profiles from WFPC2 and KPNO data at overlapping radii between 5'' and 10''. We measured the average R − I color for these radii and added it to the WFPC2 profile. The two profiles were then stitched together such that their respective weights varied linearly with radius between 5'' and 10'': the WFPC2 data contribute 100% to the combined profile for r ⩽ 5'' and the KPNO data contribute 100% for r ⩾ 10''. PSF deconvolution was not necessary for the KPNO data, as they contribute to the combined surface brightness profile at radii well beyond the seeing full width at half-maximum (FWHM). Our translation of the WFPC2 profile to the R band assumes no R − I color gradient; Lauer et al. (2005) find a median color gradient, $\frac{\Delta \left(V-I \right)}{\Delta \rm log \left(\it r \right)}$, of −0.03 mag for BCGs and other core-profile galaxies.

At radii beyond 1'', isophotes of NGC 6086 all have major-axis position angles within 5° of true north, with an average apparent axis ratio of 0.7. We adopt 0° east of north as the major-axis position angle of NGC 6086, and we assume edge-on inclination. We deprojected the surface brightness using the procedure of Gebhardt et al. (1996), which assumes spheroidal isodensity contours. The resulting major- and minor-axis luminosity density profiles are shown in Figure 1.

Zoom In Zoom Out Reset image size

We performed integral-field spectroscopic observations of NGC 6086 with OH-Suppressing Infra-Red Imaging Spectrograph (OSIRIS; Larkin et al. 2006) on the 10 m W. M. Keck II telescope and Gemini Multi-Object Spectrograph North (GMOS-N; Allington-Smith et al. 2002; Hook et al. 2004) on the 8 m Gemini Observatory North telescope. The instrument properties and our observations are summarized in Table 1. Our observations with OSIRIS used the W. M. Keck Observatory LGS–AO system (Wizinowich et al. 2006; van Dam et al. 2006); the inner component of the resulting H-band (1.6 μm) PSF has an FWHM value of ≈01. The GMOS data were collected under excellent seeing conditions; images of point sources from the Gemini North Acquisition Camera⁶ indicate an I-band FWHM of 04.

In Figure 2, we display the reduced mosaic of NGC 6086 from OSIRIS, summed over all spectral channels. Usable data from OSIRIS and GMOS extend to radii of 084 and 49, respectively. For radii out to 30'', we use major-axis kinematics from CBH99, obtained with the Intermediate Dispersion Spectrograph (IDS)⁷ on the 2.5 m Isaac Newton Telescope.

Zoom In Zoom Out Reset image size

2.2.1. OSIRIS

OSIRIS is a near-infrared (NIR), IFS built for use with the Keck AO system. It features two-dimensional spatial sampling at four scales between 002 and 01. We observed NGC 6086 with the 005 spatial scale, which provided adequate signal to noise and placed several pixels within the radius of influence. To minimize noise in individual spectra, we used the broad H-band filter, which covered several Δν = 3, ¹²CO band heads at observed wavelengths from 1.54 μm to 1.71 μm (at z ≈ 0.032). We chose to detect H-band features instead of the more prominent ν = 0–2 ¹²CO band head in the K band, which suffered from higher thermal background at an observed wavelength of 2.37 μm.

We recorded nine science exposures of the galaxy center and five sky exposures of a blank field 50'' away, for total integration times of 2.25 hr and 1.25 hr, respectively. Our dithers repeated an "object–sky–object" sequence, such that every science frame was immediately preceded or followed by a sky frame. We also recorded spectra of nine spectral template stars, using the same filter and spatial scale as for NGC 6086. To measure telluric absorption, we recorded spectra of several A0V stars, covering a range of airmasses similar to those for NGC 6086 and template stars.

We used version 2.2 of the OSIRIS data reduction pipeline⁸ to subtract sky frames, correct detector artifacts, perform spatial flat-fielding, calibrate wavelengths, generate data cubes with two spatial dimensions (x,y) and one spectral dimension (λ), and construct a mosaic of NGC 6086 from multiple data cubes. The pipeline uses an archived calibration file to perform spectral extraction of the raw spectra across the detector and assemble a data cube; the calibration file was generated by illuminating individual columns of the OSIRIS lenslet array with a white light source. We used custom routines to remove additional bad pixels from detector images, extract one-dimensional stellar spectra from three-dimensional data cubes, and calibrate galaxy and template spectra for telluric absorption. Although one-dimensional stellar spectra from OSIRIS comprise an average over many spatial pixels, spatial variations in instrumental resolution are negligible relative to the velocity broadening in NGC 6086: (Δσ_inst)² ∼ 5 × 10⁻³ σ².

Contamination from telluric OH emission presents a severe challenge for observing faint, extended objects with OSIRIS. The small FOV (08 × 32 for broadband observations at 005 per spatial pixel) does not allow for in-field sky subtraction, and subtracting consecutive science and sky frames only provides partial correction, as the relative flux from different vibrational transitions in OH varies on timescales of a few minutes. After subtracting a sky frame from each science frame, we are forced to discard the spectral channels with strong residual signals from OH, which compose approximately 15% of our spectral range. In Figure 3, we illustrate representative spectra from OSIRIS and distinguish kinematic fitting regions from residual telluric features. At both ends of the H-band spectrum, atmospheric water vapor acts as an additional contaminant. We have restricted our kinematic analysis to observed wavelengths between 1.48 and 1.73 μm.

Zoom In Zoom Out Reset image size

The second challenge for studying the centers of galaxies with OSIRIS is accurate determination of the PSF. We must construct an average PSF for mosaicked data from several hours of observations over multiple nights, during which seeing conditions and the quality of AO correction can change significantly. To estimate the PSF, we recorded a one-time sequence of exposures of the LGS–AO tip/tilt star for NGC 6086, using the OSIRIS spectrograph with the same filter and spatial scale settings as for our science frames. Data cubes were then collapsed along the spectral dimension to produce images of the star. In Appendix A, we discuss different methods for estimating the PSF, and how PSF uncertainty influences our modeling results.

2.2.2. GMOS

GMOS-N is a multi-purpose spectrograph on Gemini North. GMOS includes an IFS mode, in which hexagonal lenslets divide the focal plane and fibers map the two-dimensional field to a one-dimensional slit configuration. A second set of lenslets samples a field ∼60'' from the science target, allowing for simultaneous sky subtraction. We observed the center of NGC 6086 in IFS mode with the detector's CaT filter, detecting the infrared Ca ii triplet at observed wavelengths between 0.87 and 0.90 μm. A representative spectrum from the center of NGC 6086 is shown in Figure 4. GMOS data were reduced using version 1.4 of the Gemini IRAF software package.⁹ This standard pipeline subtracts bias and overscan signals, removes cosmic rays, mosaics data from three CCDs, extracts spectra, corrects throughput variations across fibers and within individual spectra, calibrates wavelengths using arc lamp exposures, computes an average sky spectrum, and performs sky subtraction. We stored individual spectra from each GMOS exposure, along with their spatial positions relative to the center of NGC 6086, for eventual spatial binning.

Zoom In Zoom Out Reset image size

With seeing-limited spatial resolution, GMOS poorly resolves the black hole sphere of influence. Nonetheless, kinematics derived from GMOS provide a good complement to those from OSIRIS. The Ca ii triplet region in GMOS spectra has a more clearly defined continuum than H-band spectra, and with less telluric contamination, as is evident from comparing Figures 3 and 4. As a result, LOSVDs extracted from GMOS spectra have lower systematic errors than LOSVDs extracted from OSIRIS spectra. Additionally, the angular region yielding high signal-to-noise spectra from GMOS is four times larger than that for OSIRIS.

In order to attain a sufficient signal-to-noise ratio (S/N) for effective kinematic extraction, we perform spatial binning on our two-dimensional grids of spectra from OSIRIS and GMOS. For the mean-normalized galaxy spectrum Y, mean-normalized stellar template T, and LOSVD $\mathcal {L}$ from the best fit over N_c spectral channels, we define

$$\begin{equation} {\rm S}/{\rm N} \equiv \left(\sum _{i = 1}^{N_c} \left[ \, Y_i - (T *\mathcal {L})_i \, \right]^2 / \, N_c \right)^{-1/2}. \end{equation} \tag{ 1 }$$

At the very center of NGC 6086, we spatially bin spectra until S/N>20 is achieved. This requires binning 2 × 2 spatial pixels from OSIRIS; consequently, our kinematic data have a central spatial resolution of 01, similar to the PSF FWHM. At the center of the GMOS mosaic, we bin seven hexagonal pixels, corresponding to an approximate diameter of 055. The remaining spectra from each data set are grouped to match angular and radial bins defined within the orbit models and to maintain S/N between 25 and 40. Our resulting binning schemes for both OSIRIS and GMOS use only two angular bins on each of the positive (north) and negative (south) sides of the major axis. The angular bins span 0°–369 and 369–90° from the major axis. Axisymmetric models perform LOSVD fitting in one quadrant of the projected galaxy. Symmetry about the major axis is enforced by co-adding spectra from the positive and negative (east and west) sides of the minor axis, before LOSVD extraction. LOSVDs extracted from the negative (south) side of the major axis are inverted before being input to the models. We define systemic velocity relative to the template star separately for OSIRIS and GMOS data.

Additionally, a spectral binning factor is necessary to smooth over channel-to-channel noise in spectra of NGC 6086. Our final kinematic extraction uses smoothing factors of 30 and 12 spectral pixels for OSIRIS and GMOS spectra, respectively. These values are chosen by comparing the best-fit LOSVDs from a large range of smoothing factors and identifying the minimum factor above which LOSVDs in each data set are stable between −500 km s⁻¹ and 500 km s⁻¹. Our smoothing values are consistent with the range of optimal values determined by Nowak et al. (2008) for near-infrared spectra with S/N ∼ 25–50.

H-band spectra from OSIRIS contain several absorption features that are potentially useful for kinematic extraction, but some are compromised by incompletely subtracted telluric OH emission lines, which are masked from the fit. Three broad features are relatively insensitive to the narrow OH lines: the ν = 3–6 ¹²CO band head at 1.6189 μm rest, the ν = 4–7 ¹²CO band head at 1.6401 μm rest, and the Δν = 2 band of OH between 1.537 and 1.545 μm rest. Additionally, the Δν = 2 OH band between 1.526 and 1.529 μm rest does not intersect any strong telluric emission features. We have verified that these four spectral features offer a robust comparison between stellar and galaxy spectra by repeating the fits with a large range of spectral smoothing factors. When we add other features to the fit, the root-mean-squared (rms) residual (essentially S/N⁻¹) becomes unstable to small changes in spectral smoothing.

To extract LOSVDs from GMOS spectra, we analyzed the λ8498 and λ8542 lines of Ca ii. The third line in the well-known calcium triplet, λ8662, is compromised by a flat-field artifact and discarded from kinematic analysis. The thick black lines in Figures 3 and 4 indicate the spectral channels used in our final extraction of LOSVDs from OSIRIS and GMOS, respectively.

Uncertainties for each LOSVD are determined by 100 Monte Carlo trials. In each trial, random noise is added to the galaxy spectrum, according to the rms value of the original fit, and the LOSVD fitting process is repeated. At each velocity bin, the uncertainties $\sigma _{\mathcal {L}+}$ and $\sigma _{\mathcal {L}-}$ are computed from the distribution of trial LOSVD values. We then adjust the uncertainties in the wings of each LOSVD, so that $\mathcal {L} - \sigma _{\mathcal {L}-} = 0$.

Stellar template mismatch can be a major source of systematic error in determining LOSVDs (e.g., CBH99; Silge & Gebhardt 2003; Emsellem et al. 2004). To address this issue, we have directly observed a diverse set of late-type template stars. Our nine templates from OSIRIS are giant, supergiant, and dwarf stars with spectral types from G8 to M4. We have found the most appropriate template for OSIRIS spectra of NGC 6086 to be HD 110964, an M4III star. In Appendix A, we describe our method for choosing the optimal template star, and how fitting LOSVDs with different template stars influences measurements of M_• and M_⋆/L_R. We fit GMOS spectra of NGC 6086 with a single G9III template star, HD 73710, which we observed with GMOS. The calcium triplet region is less sensitive to template mismatch than other optical- and near-infrared regions used to measure kinematics (Barth 2002).

Our integral-field observations uncover complex kinematic structures within the central 3.1 kpc (49) of NGC 6086. In Figure 6, we display two-dimensional maps of kinematic moments from OSIRIS and GMOS data; we have computed v_rad, σ, h₃, and h₄ by fitting a fourth-order Gauss–Hermite polynomial to each non-parametric LOSVD. For data with modest S/N, non-parametric LOSVDs must be used with caution, as noise may falsely introduce strong non-Gaussian components to the fit. In Figure 7, we compare two estimates of v_rad and σ from OSIRIS spectra. One estimate is obtained by fitting Gaussian profiles to non-parametric LOSVDs (Figures 7(a) and (b)), and the other is obtained by forcing a Gaussian LOSVD to fit the original spectra (Figures 7(c) and (d)). In every spatial region, v_rad and σ from the two fitting options are consistent within errors, and so we can trust the non-parametric LOSVDs. GMOS data show similar agreement between the two estimates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The stellar velocity dispersions measured by OSIRIS and GMOS each peak within 250 pc of the galaxy center, but not at the central spatial bin. Central decreases in velocity dispersion have been observed in several other early-type galaxies with known black holes (e.g., van der Marel 1994; Pinkney et al. 2003; Gebhardt et al. 2007; Nowak et al. 2008). Possible physical explanations include an unresolved stellar disk or a localized population of young stars. No dust features are present in photometry of NGC 6086, nor is there any evidence of an active galactic nucleus. The radial velocities are highly disturbed within the central 200 pc, which are only resolved by OSIRIS: at maximum, Δv_rad = 194 ± 52 km s⁻¹. Gebhardt et al. (2007) found similar patterns in v_rad and σ in the central 100 pc of NGC 1399, which were reproduced by models with a high prevalence of tangential orbits. Likewise, our best-fitting model of NGC 6086 is tangentially biased in the central 200 pc; the average value of σ_r/σ_t is 0.55. However, the two-dimensional structure of v_rad is not consistent with a resolved stellar disk. Axisymmetric modeling of elliptical galaxies by Gebhardt et al. (2003), Shapiro et al. (2006), and Shen & Gebhardt (2010) suggests that tangential bias is common within the black hole sphere of influence.

Figures 6(c), (d), (g), and (h) illustrate the two-dimensional behavior of the third- and fourth-order Gauss–Hermite moments, h₃ and h₄. Within errors, our measurements are largely consistent with h₃ = 0, while h₄ is significantly negative, corresponding to LOSVDs with "boxy" shapes and truncated wings.

We use major-axis kinematics from CBH99 to constrain stellar orbit models at radii out to 18.9 kpc (293), several times the extent of our IFS data. To incorporate these data into our models, we have adopted higher uncertainties than the values quoted in CBH99; our treatment attempts to account for additional systematic errors, which are described by CBH99 but excluded from their published measurements for NGC 6086. In Figure 8, we compare kinematic moments from CBH99 to the moments derived from OSIRIS and GMOS, selecting the spatial bins along the galaxy's major axis. We invert the sign of v_rad and h₃ for bins on the southern half of the galaxy. Values of v_rad, σ, and h₃ measured from the three data sets largely agree, although the OSIRIS data yield somewhat smaller values of σ. At radii between 06 and 49, GMOS spectra from the southern half of the galaxy yield significantly lower values of σ than spectra from the northern half. The average asymmetry is 35 km s⁻¹; our median error for individual GMOS measurements is 16 km s⁻¹. The asymmetry is not seen in long-slit data, which more consistently agree with GMOS along the north side of the major axis. In spatial bins corresponding to the minor axis, v_rad and σ behave similarly to the major-axis trends depicted in Figure 8 and agree with minor-axis kinematics from Loubser et al. (2008). Beyond the central 200 pc, we find no convincing signs of kinematically distinct stellar populations dominating galaxy spectra at 0.5, 0.9, and 1.6 μm.

Zoom In Zoom Out Reset image size

The most significant discrepancy in the major-axis kinematics is between the negative values of h₄ derived from IFS data and the positive values of h₄ measured by CBH99 (Figure 8(d)). Given the uniformity of the positive and negative values over many radial bins, we attribute this discrepancy to systematic errors in at least one set of measurements. As noted by CBH99 (and references therein), stellar template mismatch can bias h₄; however, this study and CBH99 both perform careful analysis with multiple stellar templates (9 and 28 stars, respectively). Nowak et al. (2008) demonstrated that large spectral smoothing factors can also bias h₄ to negative values. Still, a smoothing factor >100 would be necessary to produce the full discrepancy between our values and those from CBH99. Another source of error could be adjustments to the equivalent widths (EWs) of absorption features in galaxy and template spectra; we address this issue in Appendix A. Regardless of the cause, including discrepant data in stellar orbit models can influence the best-fit solutions. We discuss the effects on measurements of M_• and M_⋆/L_R in Section 4.3.2.

We generate stellar orbit models of NGC 6086, using the static potential method introduced by Schwarzschild (1979). We use the axisymmetric modeling algorithm described in detail in Gebhardt et al. (2000b, 2003); Thomas et al. (2004, 2005); and Siopis et al. (2009). Here we provide a summary of the procedure. Similar models are presented in Richstone & Tremaine (1984), Rix et al. (1997), Cretton et al. (1999), and Valluri et al. (2004).

We assume that the central region of NGC 6086 consists of three mass components—stars, a central black hole, and an extended dark matter halo—described by the radial density profile

$$\begin{equation} \rho (r) = \frac{M_\star }{L_R} \nu (r) + M_\bullet \delta (r) + \rho _{\rm halo}(r) \,. \end{equation} \tag{ 2 }$$

The stellar distribution is assumed to follow the observed (deprojected) luminosity density ν(r) (see Figure 1) with a constant stellar mass-to-light ratio M_⋆/L_R. For the dark matter halo, we compare two density profiles: the commonly used Navarro–Frenk–White (NFW) form (Navarro et al. 1996) and a logarithmic (LOG) profile:¹⁰

$$\begin{equation} \rho _{\rm halo}(r) = \frac{v_c^2}{4\pi G} \frac{3r_c^2 + r^2}{\big(r_c^2 + r^2 \big) ^2} \,. \end{equation} \tag{ 3 }$$

The free parameters in the LOG profile are the asymptotic circular speed v_c and the core radius r_c, within which the density is approximately constant. The enclosed halo mass for this profile,

$$\begin{equation} M_{\rm halo}(<r) = \frac{v_c^2 \, r}{G} \left(1 - \frac{r_c^2}{r_c^2 + r^2} \right) \;, \end{equation} \tag{ 4 }$$

is predominantly set by v_c. The difference between the NFW and LOG profiles is greatest at small radii, where the NFW profile yields higher densities, ρ_halo ∝ r⁻¹. However, each profile is greatly exceeded by the stellar mass density in the inner regions of NGC 6086. As described below, we have compared LOG and NFW profiles in a subset of models and find no significant differences in the best-fit values of M_•.

For a given set of input parameters M_•, M_⋆/L_R, and ρ_halo, we compute a continuous, static gravitational potential from Equation (2). Azimuthal symmetry about the z-axis (corresponding to the projected minor axis) is imposed, as well as symmetry about the equatorial plane (z = 0). We then generate stellar orbits by propagating test particles through the potential. Orbits are tracked in a finely spaced polar grid, (r, θ), where θ is the polar angle from the z-axis. Our models of NGC 6086 use 96 radial and 20 polar bins per quadrant. Each orbit is sampled at a random set of azimuthal angles, ϕ.

The initial phase space coordinates of test particles are chosen to sample thoroughly three integrals of motion: energy E, angular momentum component L_z, and non-classical integral I₃. Computational noise and finite propagation steps introduce noise into test-particle trajectories; this is mitigated by allowing each particle to complete 200 circuits of the potential and then determining its average orbit. Orbits that escape the potential are not included in subsequent fitting. For a given potential, our model of NGC 6086 includes approximately 16, 000–19, 000 bound orbits. Identical counterparts with the opposite sign of L_z raise the total to 32,000–38,000 orbits. Each orbit in the model is assigned a scalar weight; initially, all bound orbits are given equal weights.

The set of best-fit orbital weights is determined by comparing projected LOSVDs from the orbits to the observed LOSVDs for the galaxy. Our models use non-parametric LOSVDs, defined in 15 velocity bins between −1000 and 1000 km s⁻¹. Each observed LOSVD spatially maps to a linear combination of bins within the model, according to the spatial boundaries of the corresponding spectrum and to the instrument-specific PSF. A corresponding model LOSVD is computed from the projected velocity distributions of individual orbits in each spatial bin, the appropriate combination of spatial bins, and the orbital weights. Only the orbital weights are varied to determine the best-fit solution.

The best-fit solution is determined by the method of maximum entropy, as in Richstone & Tremaine (1988). This method maximizes the function f ≡ S − αχ², where

$$\begin{equation} \chi ^2 = \sum _i^{N_b} \sum _j \frac{[ \mathcal {L}_{i, \rm data}(v_j) - \mathcal {L}_{i, \rm model}(v_j)] ^2 }{\sigma ^2_i (v_j)} \end{equation} \tag{ 5 }$$

and

$$\begin{equation} S = - \sum _k w_k \, {\rm ln} \left(\frac{w_k}{V_k} \right). \end{equation} \tag{ 6 }$$

Here, $\mathcal {L}_{i,\rm data}$ and $\mathcal {L}_{i,\rm model}$ are LOSVDs in each of the i = 1,  ...,  N_b spatial bins, σ²_i(v_j) is the squared uncertainty in $\mathcal {L}_{i,\rm data}$ at velocity bin v_j, w_k is the orbital weight for the kth orbit, and V_k is the phase volume of the kth orbit. The parameter α is initially small so as to distribute orbital weights broadly over phase space and is increased over successive iterations so that the final optimization steps exclusively minimize χ². A further constraint for all solutions is that the summed spatial distribution of all weighted orbits must match the observed luminosity density profile.

Our combination of IFS and long-slit data within 30'' is not sufficient to measure directly the dark matter halo profile or enclosed mass. We therefore consider different halo masses and profiles in our analysis. For a given dark matter halo profile, we generate a set of 32,000–38,000 stellar orbits for an input M_• and M_⋆/L_R and obtain χ² for the best-fit orbital weights using the method described in the previous subsection. This process is repeated over a finely sampled grid in M_• and M_⋆/L_R. We have completed all trials for NGC 6086 on supercomputers at the Texas Advanced Computing Center (TACC), totaling ∼10, 000 CPU hours. Our results are summarized in Table 2.

In Figure 9, we illustrate how χ² in the orbit models varies with M_• and M_⋆/L_R. Three dark matter halo masses are shown: no dark matter (left panel), an intermediate LOG halo with v_c = 300 km s⁻¹ (middle), and a maximal LOG halo with v_c = 500 km s⁻¹ (right). The latter is chosen to approximate the measured line-of-sight velocity dispersion of 302 km s⁻¹ for NGC 6086's host cluster (Zabludoff et al. 1993), which corresponds to a full three-dimensional velocity dispersion of 523 km s⁻¹. We set the core radius in Equation (3) to be r_c = 8.0 kpc, reflecting the value of 8.2 kpc determined by Thomas et al. (2007) for NGC 4889, the Coma BCG. Our dynamical models are constrained within a radius of 18.9 kpc, corresponding to outermost radius of 293 for long-slit data in CBH99. We also have run models with a single NFW dark matter profile, constructed to contain the same enclosed mass within 18.9 kpc as our most massive LOG halo. In cosmological N-body simulations, the NFW scaling parameters c and r_s are correlated according to the relationship

$$\begin{equation} r_s^3 = \left(\frac{3 \times 10^{13} \, M_\odot }{200 \frac{4\pi }{3} \rho _{\rm crit} \, c^3} \right) 10^{\frac{1}{0.15} \left(1.05 - \rm log_{10} \it c \right)} \end{equation} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

(Navarro et al. 1996; Rix et al. 1997),¹¹ where ρ_crit = 3H²₀/8πG. Combining this relationship with our enclosed mass constraint, we obtain c = 9.6 and r_s = 94.0 kpc.

Dark matter is ubiquitous in galaxies, thus motivating its inclusion in stellar orbit models. Furthermore, models with dark matter produce better fits to our full set of kinematics: when the dark matter component is removed, χ²_min increases by ∼100 (Table 2). In Appendix B, we compare our full sets of observed LOSVDs to the best-fitting models with and without dark matter. The largest discrepancies between the two models occur at radii beyond 15'', where we only have a few data points from CBH99. Without thorough radial coverage or multiple long-slit position angles, we cannot fully untangle degeneracies between M_⋆/L_R, v_c, and r_c (or M_⋆/L_R, c, and r_s in the case of an NFW profile).

In Figure 10, we display the dark matter fraction as a function of radius, for each of the halos described above. In each case, we use the best-fit values of M_• and M_⋆/L_R, described below, to compute the total enclosed mass. Using our surface brightness profile from HST/KPNO, we compute an effective radius, R_eff, of 317 (20.4 kpc), defined as the semimajor axis of the elliptical isophote containing half of the total luminosity. We assume a total luminosity of 1.82 × 10¹¹ L_☉,R from M_V = −23.11 in Lauer et al. (2007) and V − R = 0.64. Within R_eff, dark matter composes 44%–74% of the total mass (for the intermediate-mass LOG halo and the NFW halo, respectively). This range agrees with dynamical models of other cD galaxies with LOG and NFW halos: Thomas et al. (2007) found ∼50%–75% dark matter within R_eff for NGC 4889 and NGC 4874, and Gebhardt & Thomas (2009) found ∼40% dark matter within R_eff for M87. Within the 49 outer radius of GMOS data, the maximum dark matter fraction in our models is 20%.

Zoom In Zoom Out Reset image size

In Figure 11, we illustrate the variation of M_• and M_⋆/L_R with enclosed halo mass. We compute the best-fit values for M_• and M_⋆/L_R by integrating the two-dimensional likelihood function from each χ² surface; we describe this method and our determination of errors in Section 4.3. Figures 9 and 11 show that the best-fit values of M_• and M_⋆/L_R are substantially influenced by the presence of dark matter in the stellar orbit models. This occurs because our innermost kinematics sample the nucleus of the galaxy, where orbits are dominated by the enclosed mass of stars and the central black hole, whereas enclosed stellar and dark halo masses are both important at larger radii. The significant presence of dark matter at large radii drives the best-fit models to lower values of M_⋆/L_R. In turn, the decreased stellar mass requires a higher black hole mass to reproduce the kinematics in the nucleus. This trend was initially demonstrated by Gebhardt & Thomas (2009) for M87. Negative covariance between M_• and M_⋆/L_R is also visible in χ² contours for individual dark matter halo models (Figure 9). In Section 5, we compare our best-fit values of M_• from different dark matter halo models to the predictions from the M_•–σ and M_•–L relationships.

Zoom In Zoom Out Reset image size

Using an NFW profile yields a 7% decrease in the best-fit value of M_⋆/L_R, relative to the LOG profile with the same enclosed mass at 18.9 kpc. This is because the majority of our kinematic measurements occur at r ⩽ 5 kpc, where the more centrally concentrated NFW profile yields higher dark matter densities. In the central 100 pc, the stellar core of NGC 6086 varies nearly as r⁻¹ in luminosity density, mimicking the slope of the NFW profile. Near the black hole, the lower stellar mass density balances the higher density in dark matter, and so the best-fit values of M_• are identical for the NFW and LOG profiles. We find χ²_min, LOG − χ²_min, NFW = 0.8, indicating no significant difference in the goodness of fit. Thomas et al. (2005, 2007) and Gebhardt & Thomas (2009) have found similar difficulties in distinguishing between NFW and LOG profiles.

A second way to address the influences of dark matter on M_⋆/L_R and M_• is to fit the orbit models only at radii where dark matter composes a small fraction of the enclosed mass. For NGC 6086, we have run two trials in which we only fit LOSVDs from OSIRIS and GMOS: one trial with the maximum-mass LOG dark matter halo described above and one trial with no dark matter. In both of these trials, the best-fit values of M_• and M_⋆/L_R agree with our results from fitting IFS and long-slit data with the maximum-mass LOG halo (see Table 2). This agreement provides strong evidence that M_• ∼ 3 × 10⁹ M_☉, regardless of our insensitivity to the exact structure of the dark matter halo in NGC 6086. In contrast, forcing models without dark matter to fit long-slit data biases the best-fit black hole mass to a substantially lower value, ∼6 × 10⁸ M_☉, and increases the best-fit stellar mass-to-light ratio to 6.7 M_☉ L⁻¹_☉. This mass-to-light ratio is highly inconsistent with stellar population estimates, as we discuss in Section 5. Even though excluding long-slit data results in consistency between models with and without dark matter, these models are not as thoroughly constrained, and we obtain slightly larger confidence intervals in M_• and M_⋆/L_R for each trial.

Our data from OSIRIS, GMOS, and CBH99 play complementary roles in constraining the gravitational potential of NGC 6086. In Figure 12, we compare model results for one dark matter halo (LOG; v_c = 500 km s⁻¹), using data only from GMOS and CBH99, versus only from OSIRIS and CBH99. The GMOS data are sufficient to detect a black hole, in part because of excellent seeing. Yet the strong diagonal contours in the left panel of Figure 12 indicate that the black hole mass derived from GMOS is degenerate with the enclosed stellar mass. LOSVDs from OSIRIS have large statistical errors and by themselves cannot place strong constraints on the black hole mass. However, the OSIRIS data help separate the respective influences of the stars and the black hole. Using GMOS and OSIRIS data together reduces covariance between M_• and M_⋆/L_R and lowers the statistical uncertainties of both quantities. Long-slit data from CBH99 confirm the presence of dark matter and tighten constraints on M_⋆/L_R and M_• for individual dark matter halo models.

Zoom In Zoom Out Reset image size

It is not clear how to interpret the significant increase in the best-fit value of M_•, from 1.9^+1.5_−1.1 × 10⁹ M_☉ with OSIRIS and long-slit data only, to 7⁺³₋₃ × 10⁹ M_☉ with GMOS and long-slit data only. Using central σ = 329 km s⁻¹ from Loubser et al. (2008), and M_• = 3.6 × 10⁹ M_☉ from our combined-data trial, we compute r_inf = 022. In this case, GMOS marginally resolves the sphere of influence, with a seeing FWHM ∼2r_inf. At small radii, LOSVDs from GMOS yield slightly higher velocity dispersions than overlapping LOSVDs from OSIRIS (Figure 8(b)), which could contribute to the increase in M_•. With no obvious way to assess independently the accuracy of data from OSIRIS versus GMOS, we favor including both sets of LOSVDs. The corresponding black hole mass is 3.6^+1.7_−1.1 ×  10⁹ M_☉, which lies between the two partial-data values and has the narrowest confidence interval. The confidence interval in M_⋆/L_R is also minimized by including all data, though the best-fit value of 4.6^+0.3_−0.7 does not change significantly upon exclusion of OSIRIS or GMOS data.

None of the χ² surfaces in Figures 9 and 12 are completely smooth. In particular, the models without any dark matter show large variations in χ² over small changes in M_• and M_⋆/L_R (Figure 9). We suspect that these variations arise from numerical noise in propagating test particles through different potentials: each small change in the potential may send a given test particle through a different set of spatial regions. The cumulative effect is that each model creates a different set of test-particle LOSVDs, and χ² can change abruptly in spite of the freedom to adjust orbital weights. When the models include a constant dark matter component, the relative changes in the potential are smaller, and the noise in χ² is less pronounced. Still, the χ² surface for each dark matter halo exhibits a noise floor at the level of Δχ² ∼ 1. In Section 4.3.1, we describe how this noise influences our measurements of confidence intervals.

The figure of merit for evaluating confidence intervals in M_• and M_⋆/L_R is Δχ² ≡ χ² − χ²_min, where χ²_min is the lowest output value among all models. For NGC 6086, we determined confidence intervals by integrating the relative likelihood function, $P \propto e^{-\frac{1}{2}\Delta \chi ^2}$. Although Δχ² is a better statistical indicator than χ² per degree of freedom (e.g., van der Marel et al. 1998; Gebhardt et al. 2003), the latter is useful for crudely indicating the level of agreement between the data and the model with the best fit. For each model, the number of degrees of freedom, N_dof, depends on the number of observed LOSVDs and the number of velocity bins evaluated per LOSVD. Because we used a spectral smoothing factor in determining LOSVDs for OSIRIS and GMOS, the velocity bins are not entirely independent. We estimate that each LOSVD from OSIRIS or GMOS has 1 degree of freedom per 2 velocity bins, whereas long-slit data from CBH99 has 1 degree of freedom per velocity bin. For all experiments with NGC 6086, we find χ²_min/N_dof between 0.7 and 1.8, indicating reasonable agreement.

4.3.1. Confidence Intervals

We determine confidence intervals by using Δχ² as an empirical measure of relative likelihood between models with different M_• and M_⋆/L_R and by numerically integrating this likelihood with respect to M_• and M_⋆/L_R. In contrast to the majority of previous studies (e.g., Gebhardt et al. 2000b, 2003, 2007; Nowak et al. 2007, 2008; Gebhardt & Thomas 2009; Gültekin et al. 2009b; Siopis et al. 2009; Shen & Gebhardt 2010; cf. van der Marel et al. 1998), we do not use fixed values of Δχ² to define confidence intervals. The fixed Δχ² method is appropriate only if the orbit models cleanly sample a well-defined likelihood function of M_• and M_⋆/L_R; other studies typically assume a two-dimensional Gaussian likelihood function. Our models of NGC 6086 produce noisy χ² contours (Figure 9), in which case the fixed Δχ² method is sensitive to noise in individual models. This effect is especially pronounced for models without a dark matter halo.

We define likelihood P such that two models with χ²₁ and χ²₂ have relative likelihood:

$$\begin{equation} \frac{P_1}{P_2} = e^{-\frac{1}{2} \left(\chi ^2_1 - \chi ^2_2 \right)}. \end{equation} \tag{ 8 }$$

This form of P is valid, provided that χ² is measured from independent, Gaussian-distributed data points (Cowan 1998); in our case, these are the observed LOSVDs. To evaluate likelihood with respect to a single variable (i.e., x ≡ M_•), we marginalize the two-dimensional surface with respect to the other variable (y≡ M_⋆/L_R), such that

$$\begin{equation} P \left(x \right) \propto \sum _{y_{\rm min}}^{y_{\rm max}} e^{-\frac{1}{2} \chi ^2 \left(x, y \right)} \; \delta y {,} \end{equation} \tag{ 9 }$$

where δy is the interval between sampled values of y. Confidence intervals in x are determined by evaluating the cumulative distribution:

$$\begin{equation} C(x) = \frac{ \int _{x_{\rm min}}^x P (x^{\prime }) dx^{\prime } }{\int _{x_{\rm min}}^{x_{\rm max}} P (x^{\prime }) dx^{\prime }}. \end{equation} \tag{ 10 }$$

In practice, we define [x_min, x_max] and [y_min, y_max] by expanding our range of models until the marginalized likelihood functions P(M_•) and P(M_⋆/L_R) are nearly zero at the minimum and maximum modeled values of M_• and M_⋆/L_R. The physical limit M_• = 0 is included in all trials. For confidence level k, we define confidence limits at $C = \frac{1}{2} \left(1 \pm k \right)$. For example, the 68% confidence interval comprises all x for which 0.16 ⩽ C ⩽ 0.84. For each trial, we define the best-fit values of M_• and M_⋆/L_R as the median values from P(M_•) and P(M_⋆/L_R), corresponding to $C = \frac{1}{2}$. In Figure 13, we show P(M_•) and P(M_⋆/L_R) for each dark halo setting, along with cumulative distributions, median values, and confidence intervals. To estimate precise confidence limits, we linearly interpolate C between discretely sampled values of M_• and M_⋆/L_R.

Zoom In Zoom Out Reset image size

Our empirical treatment yields wider 68% confidence intervals than those derived from fixed Δχ² and a Gaussian likelihood function. Our intervals for confidence levels ⩾90% typically fall near those derived from fixed Δχ². By construction, our confidence intervals do not include M_• = 0. For the maximum-mass LOG dark matter halo in NGC 6086, the marginalized likelihood corresponding to M_• = 0, P(M_• = 0) is 0.06% of the maximum marginalized likelihood value. For a Gaussian likelihood function, this likelihood ratio would indicate >99.98% confidence for a black hole detection. For models without dark matter, P(M_• = 0) is 23.2% of the maximum value, and the detection falls to 91% confidence.

An alternative way to determine confidence intervals from noisy χ² is the method of van der Marel et al. (1998), in which random noise is added to the LOSVDs output by each orbit model, the χ² surface is re-computed, and confidence limits are determined using fixed Δχ². This process is repeated in a Monte Carlo fashion, and the extrema of the confidence limits from all trials are adopted. This treatment assumes that numeric noise in orbit models produces fluctuations about an intrinsically Gaussian likelihood distribution. The global likelihood function of our models, however, is visibly non-Gaussian (see Figure 13, top, in particular).

4.3.2. Systematic Errors

The confidence intervals measured from χ² surfaces account for statistical errors in the observed LOSVDs and random noise within the stellar orbit models. Systematic errors must be addressed separately. We have directly tested several systematic effects.

Our largest systematic error arises from discrepancies between LOSVD shapes derived from IFS versus long-slit data, as indicated by the parameter h₄ (Figure 8(d)). To test how this discrepancy biases M_• and M_⋆/L_R, we fit Gaussian profiles to the original LOSVDs from CBH99, constructing an alternative set of LOSVDs with h₃ = 0 and h₄ = 0. Using these new LOSVDs in combination with the OSIRIS and GMOS data, we repeated the stellar orbit models for the maximum dark matter halo and found M_• = 4.5^+1.4_−1.1 × 10⁹ M_☉ and M_⋆/L_R = 4.2^+0.3_−0.3 M_☉ L_☉⁻¹. This is a 22% increase in M_• and a 9% decrease in M_⋆/L_R, relative to the corresponding trial with original LOSVDs from CBH99. The direction of the bias indicates that more enclosed mass is required to produce LOSVDs with extended wings (h₄ > 0) at ∼3–19 kpc. Since the additional mass is obtained by increasing M_⋆/L_R, the innermost LOSVDs drive the best fit toward a smaller value of M_•.

Two additional sources of systematic error are the uncertainties of the optimal stellar template for fitting LOSVDs and the average PSF for OSIRIS. We find that both effects yield small errors. Trials with three PSFs yield a standard deviation of 8% in M_• and 1% in M_⋆/L_R. Trials with two stellar templates differ by 5% in M_• (corresponding to 3.6% standard deviation) and 2% in M_⋆/L_R (1.6% standard deviation). We describe these tests in detail in Appendix A.

We use the following prescription to compute the total error in M_• for a given dark halo:

$$\begin{equation} \sigma _{+, \, \rm tot} = \big(\sigma _{+, \, \chi ^2}^2 + \sigma _{\rm PSF}^2 + \sigma _{\rm temp}^2 + 2\delta _{h_4}^2 \big) ^ {\frac{1}{2}} \end{equation} \tag{ 11 }$$

and

$$\begin{equation} \sigma _{-, \, \rm tot} = \big(\sigma _{-, \, \chi ^2}^2 + \sigma _{\rm PSF}^2 + \sigma _{\rm temp}^2 \big) ^ {\frac{1}{2}} \;. \end{equation} \tag{ 12 }$$

Here, σ_+, tot and σ_−, tot are the upper and lower portions of the 68% confidence interval, including all errors; $\sigma _{+, \, \chi ^2}$ and $\sigma _{-, \, \chi ^2}$ are the contributions from integrating the empirical likelihood function as in Section 4.3.1 and represent statistical errors; σ_PSF is the standard deviation in best-fit M_• from trails with different OSIRIS PSFs; σ_temp is the standard deviation from trials with different template stars; and $\delta _{h_4}$ is the difference in best-fit M_• from using re-fit versus original LOSVDs from CBH99. As excluding h₃ and h₄ from these LOSVDs introduces a bias toward higher M_•, we assign this effect solely to σ_+, tot, with a magnitude of $2 \times \frac{1}{\sqrt{2}}\delta _{h_4}$. Our equations defining σ_+, tot and σ_−, tot for M_⋆/L_R are similar to Equations (11) and (12), except we apply $\delta _{h_4}$ entirely to σ_−, tot. To apply our results to different dark matter halo settings, we define the systematic errors as percentages, such that $\delta _{h_4}$, σ_PSF, and σ_temp scale with M_• and M_⋆/L_R. We list the 68% confidence intervals, including all errors, inside parentheses in Table 2. The confidence intervals outside parentheses in Table 2 only include $\sigma ^2_{\pm, \, \chi ^2}$.

By adding in quadrature, we have assumed zero correlation between different sources of systematic error; this is the most conservative approximation. $\delta _{h_4}$ and σ_temp are likely correlated (see, e.g., CBH99), but the contributions from σ_temp are small, and we have not run extensive tests to measure covariance between stellar templates and overall trends in h₄. There is no obvious reason to suspect covariance between other terms.

4.3.3. Other Potential Sources of Error