ORBIT-BASED DYNAMICAL MODELS OF THE SOMBRERO GALAXY (NGC 4594)

John R. Jardel; Karl Gebhardt; Juntai Shen; David B. Fisher; John Kormendy; Jeffry Kinzler; Tod R. Lauer; Douglas Richstone; K. Gültekin

doi:10.1088/0004-637X/739/1/21

1. INTRODUCTION

Most galaxies are thought to host supermassive black holes (SMBHs) at their centers. The masses of these SMBHs have been observed to correlate with several properties of their host elliptical galaxies and of the classical bulge components of their host disk galaxies. For example, M_• correlates with galaxy/bulge mass (Dressler 1989; Magorrian et al. 1998; Laor 2001; McLure & Dunlop 2002; Marconi & Hunt 2003; Häring & Rix 2004), luminosity (the M–L relation) (Kormendy 1993; Kormendy & Richstone 1995; Kormendy & Gebhardt 2001; Gültekin et al. 2009b), velocity dispersion (the M–σ relation) (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Gültekin et al. 2009b), and globular cluster (GC) content (Burkert & Tremaine 2010; Harris & Harris 2011). These and other similar correlations suggest that galaxy formation and black hole (BH) growth are fundamentally linked. To better understand this interplay, accurate BH masses are needed.

One challenge that limits the accuracy is the determination of the host galaxy's inclination. Projection effects are difficult to model and cause loss of information, leading to systematic uncertainties. Therefore, SMBHs in galaxies whose inclination is confidently known have the best chance of being accurately and robustly measured. Another issue that limits the accuracy is the effect that a dark matter (DM) halo has on the determination of SMBH mass. Gebhardt & Thomas (2009) show that orbit-based models can underestimate BH mass when DM is not considered in the modeling; however, Schulze & Gebhardt (2011) find the effect is small when the BH's sphere of influence is well resolved.

NGC 4594 (M104, or the Sombrero Galaxy) is a nearly edge-on Sa type spiral with a prominent stellar disk and large, classical bulge (Kormendy & Kennicutt 2004). The shape of this disk indicates that it (and thus the entire galaxy) is inclined at an angle very close to 90°. Throughout this paper we assume a distance to NGC 4594 of 9.8 Mpc, calculated from surface brightness fluctuations (Tonry et al. 2001). Unless otherwise stated, all distance-dependent quantities are scaled to this value. Tonry et al. (2001) use a value of H₀ = 74 km s⁻¹ Mpc⁻¹ in their distance determinations; however, we compare our M_• and L_V to Gültekin et al. (2009b) who adopt H₀ = 70 in their work. We therefore scale the Gültekin et al. (2009b) distances down by 6%. BH mass scales as M_• ∝ D and luminosity as L_V∝D⁻²; these quantities are adjusted accordingly.

NGC 4594 was one of the first galaxies in which a BH was detected and it has a long history of study. Kormendy (1988, hereafter K88) first found evidence for a massive BH of M_• = 5.4^+11.8_{− 3.7} × 10⁸ M_☉ using only ground-based observations. With isotropic Jeans models, Emsellem et al. (1994b) measured a BH of mass M_• = (5.4 ± 0.5) × 10⁸ M_☉. Later, Kormendy et al. (1996, hereafter K96) used high-resolution kinematics from the Faint Object Spectrograph (FOS) on the Hubble Space Telescope (HST)—the same data set we include—to measure logM_• = 8.8 ± 0.5 M_☉. This corresponds to a mass of 5.8^+12.4_{− 4.0} × 10⁸ M_☉. With isotropic models, Magorrian et al. (1998) obtained a value of M_• = 6.9^+0.2_{− 0.1} × 10⁸ M_☉. These values for M_• all lie toward the high-mass end of the M–σ and M–L relations. Massive SMBH measurements are frequently being revised, and we expect the confidently known inclination of NGC 4594 to lead to one of the more secure measurements of a high-mass SMBH.

We present new Gemini spectroscopy of the major axis as well as Spectroscopic Areal Unit for Research on Optical Nebula (SAURON) integral field kinematics covering the central region of the galaxy. We also use high-resolution HST/FOS kinematics of the nucleus and kinematics derived from GCs at large radii. We combine these kinematic data sets with HST and ground-based photometry to run axisymmetric orbit-based models. These models allow us to measure the BH mass, stellar mass-to-light ratio (M/L), and DM halo of NGC 4594. In addition, we recover information about the internal orbit structure of the galaxy.

2. DATA REDUCTION AND ANALYSIS

Dynamical modeling requires as input the three-dimensional luminosity density distribution ν(r) as well as the line-of-sight velocity distribution (LOSVD) at many locations in the galaxy. We use HST and ground-based images for the photometry. Our kinematics include high-resolution HST/FOS spectra, long-slit spectra from the Gemini Near-Infrared Spectrograph (GNIRS) on Gemini, SAURON integral field kinematics, and individual velocities of GCs. We discuss each in turn.

2.1. Photometry

In order to cover a large enough dynamical range to have leverage on both the central SMBH and dark halo, we use surface brightness profiles from HST and ground-based images. The stellar disk of NGC 4594 dominates at intermediate radii on the major axis causing the isophotes in this region to be substantially flattened. This abrupt change in ellipticity introduces an additional challenge to the deprojection. Our standard technique is to assume that the surfaces of constant luminosity density ν are coaxial, similar spheroids (Gebhardt et al. 1996). Clearly the presence of a disk invalidates this assumption, so we decompose the surface brightness into bulge and disk components, deprojecting each separately so that our assumption holds for each component. Afterward, we re-combine the deprojected profile of each component ν_bulge + ν_disk and input the total ν(r) into our modeling program.

The bulge–disk decomposition fits directly to a projected image. We construct a model disk by considering a Sérsic (1968) profile

$\begin{equation} \mu (R) = \mu _0 \mathrm{ exp}[-(R/R_0)^{1/n}], \end{equation} \tag{ 1 }$

where μ₀ is the surface brightness at R = R₀ and n is the Sérsic index. For n = 1, the profile is an exponential. For inclinations other than 90°, the projection of our disk model is an ellipse. By specifying the inclination of the disk i, the axial ratio b/a of the ellipse is given by b/a = cos i for a thin disk.

We construct many disk models by varying μ₀, R₀, i, and n (keeping n close to 1). Each model is then subtracted from the image until the residual brightness distribution has elliptical isophotes. The remaining light is assigned to the bulge component. A one-dimensional major axis bulge profile is produced by averaging the bulge light in elliptical, annular isophotes. Hence, we are left with an analytic disk model and a non-parametric bulge model. We identify the best bulge and disk models as those that minimize the rms residuals of the model-subtracted image.

In addition to the obvious main disk, NGC 4594 hosts a well-studied nuclear disk (Burkhead 1986, K88, K96) at small radii. We fit the nuclear disk in the HST image and the main disk in the ground-based image. Because we fit directly to the images, dust lanes and object masking become important. We keep a bad pixel list which instructs our code to ignore trouble spots. Dust lanes are selected by eye, while SExtractor (Bertin & Arnouts 1996) is used to identify foreground stars, background galaxies, and GCs.

2.1.1. HST Image

To probe the nuclear region, we use a point-spread function (PSF)-deconvolved HST Wide Field Planetary Camera 2 (WFPC2) image (GO-5512; PI: Faber). This image is presented in K96 and provides an excellent view of the central region of the galaxy. Centered on the PC1 camera, the image is taken in the F547M filter and has a scale of 0 farcs 0455 pixel⁻¹ of the central 34'' × 34'' of the galaxy. The PSF deconvolution uses the Lucy–Richardson algorithm (Richardson 1972; Lucy 1974) for 40 iterations and is well tested on WFPC2 images (Lauer et al. 1998). The best-fit parameters from the bulge–disk decomposition are listed in Table 1.

Table 1. Summary of Disk Parameter Fit

Disk	μ₀ (mag arcsec⁻²)	R₀ (arcsec)	n	i
Outer 1	18.8	66.8	1.0	80
Outer 2	16.7	40.1	1.0	80
Nuclear	20.4	4.1	1.1	83

Download table as: ASCII Typeset image

NGC 4594 is also thought to have weak LINER emission (Bendo et al. 2006) and there is a point source in the HST image. Furthermore, heavy dust absorption also makes the determination of the central bulge surface brightness profile difficult for R lesssim 0 farcs 17. To deal with these issues, we extrapolate the bulge surface brightness μ_bulge(R) inward to R = 0 farcs 02 with a constant slope fit to the region near R = 0 farcs 17. Figure 1 shows the result of this extrapolation as well as the fits of the other components in the ground-based image.

**Figure 1.** Radial plot of all resulting components from our bulge–disk decomposition along the major axis. Dotted lines are the disk profiles with parameters from Table 1. The solid black line is the total surface brightness (bulge + disks). Solid colored lines are the bulge profiles, the red line is the result from fits to the *HST* image, and the blue line is from the ground-based image. Diamonds indicate the raw *HST* bulge profile before we apply our dust correction and point-source removal. The gap between the *HST* and ground-based bulge profiles is interpolated over before deprojection. Plotted in green is the globular cluster surface brightness profile, arbitrarily scaled to match the stellar surface brightness at its innermost point.
Download figure:
Standard image High-resolution image

2.1.2. Ground-based Image

We obtained a wide-field, I-band image from the Prime Focus Camera on the McDonald 0.8 m telescope. This instrument provides a large unvignetted field of view (45 × 45 arcmin²) and a single CCD detector. Therefore, we can more robustly carry out sky subtraction and accurately constrain the faint isophotes. The image is corrected for bias, flat field, and illumination using standard routines in IRAF.

In our fit to the ground-based image, we ignore the central 20'' due to overexposure and contamination from the nuclear disk. We attempt fits in the region 20''–900'' with only one stellar disk given by Equation (1); however, these produce unacceptable residuals. Instead of modifying Equation (1), we add a second disk (in addition to the nuclear disk fit only in the HST image). This approach is similar to the Multi-Gaussian Expansion technique used to model the light distribution of bulges and ellipticals (Emsellem et al. 1994a). A summary of the fits of all components is given in Table 1 and plotted in Figure 1.

2.1.3. Globular Cluster Profile

GCs are essentially bright test particles that allow us to probe the potential at radii where the stellar light is faint. They have been used in orbit-based models of other galaxies (Gebhardt & Thomas 2009; Shen & Gebhardt 2010; Murphy et al. 2011). To include them in our models, we use the GC number density profile (Rhode & Zepf 2004) as an analog to the stellar density. The number density profile is converted to a surface brightness profile by arbitrarily adjusting the zero point to match the stellar profile in log space.

The green line in Figure 1 shows that the slope of the GC surface brightness profile is different from that of the stars. We run models using both the measured luminosity density distribution of the GCs and assuming that of the stars. We find a significant preference for the measured GC profile.

2.1.4. Bulge Profile

Our bulge–disk decomposition returns a non-parametric form of the bulge profile. It is not necessary to have a parameterized bulge profile for our dynamical models; however, we fit a Sérsic profile to ground-based bulge model using Equation (1). The bulge is well fit by a Sérsic function, with the rms residuals equal to 0.08 mag arcsec⁻². We measure μ₀ = 13.5 mag arcsec⁻², R₀ = 0 farcs 1, and n = 3.7.

We can convert the central surface brightness μ₀ and radius R₀ parameters to the more familiar "effective" parameters μ_e and R_e. The effective radius R_e is given by R_e = (b_n)ⁿR₀ and the effective surface brightness μ_e = μ₀ + 2.5log (e) b_n (MacArthur et al. 2003). The factor b_n depends on n; an expansion for b_n can be found in MacArthur et al. (2003). Applying these conversions, we obtain μ_e = 21.3 mag arcsec⁻² and R_e = 156 farcs 2.

We obtain a simpler estimate for the half-light radius of the bulge from integration of the surface brightness profile; no fitting functions are required. We estimate R_e = 117'' ± 12''. The integrated magnitudes are calculated for component x by L_x = 2π∫I_x(r)r dr. This does not take into account the ellipticity of each component, so we scale the luminosity by L^true_x ≈ (1 − epsilon _x)L_x, where epsilon _x is the ellipticity of each component, assumed to be constant with radius. The bulge profile is known to become rapidly circular for r gsim 100'' (Burkhead 1986) so our procedure almost certainly underestimates M_bulge and bulge-to-total ratio (B/T). These numbers are computed as a sanity check only and do not affect the dynamical models.

The absolute integrated magnitudes are M_disk = −21.4 and M_bulge = −22.5 in F814W, corrected for Galactic extinction along the line of sight (Schlegel et al. 1998). Using the HST calibration package SYNPHOT (described in detail below), we convert these F814W magnitudes to V-band Vega magnitudes. We obtain M_{V, bulge} = −22.1 and M_{V, disk} = −21.0. These structural parameters lie exactly on the fundamental plane for bulges and ellipticals as presented in Kormendy et al. (2009). Our integrated magnitudes translate to a bulge-to-total ratio B/T = 0.73 with the nuclear disk contributing 1% of the total light. This value of B/T is lower than previous measurements—Kormendy et al. (2011) report B/T = 0.925 ± 0.013. Our B/T, however, is in good agreement with a recent measurement by Gadotti (2011, model BD). Regardless of the value of B/T, our dynamical models are unaffected, because we add all the bulge and disk light together again after the deprojection.

We do not explore the possibility of fitting an exponential stellar halo in addition to a bulge and disk as Gadotti (2011) do. Our bulge–disk decomposition produces a non-parametric bulge profile which could in principle be a combination of a Sérsic bulge plus exponential halo. However, this resulting profile is well fit by a Sérsic function with n significantly larger than 1. We therefore do not agree with the claim made by Gadotti (2011) that the bulge of NGC 4594 is actually an exponential stellar halo.

2.1.5. Deprojection

We combine the HST and ground-based bulge profiles by zero pointing both to F814W. We calculate the F547M photometric zero point for the HST image from the SYNPHOT package in IRAF. Spectral template fitting (Section 2.2.1) shows that in the central region of the galaxy gsim 85% of the light comes from K6III stars. We therefore convert the F547M zero point to F814W with SYNPHOT using the Bruzual Atlas⁶ template for a K6III star.

Before deprojection, we extrapolate the one-dimensional profiles μ(R) with a constant slope to R = 1800''. The three disk profiles are then combined and deprojected via Abel inversion in the manner described in Gebhardt et al. (1996). We assume an inclination of i = 90°. The inclinations of the combined disk components imply an ellipticity of e = 0.83. Our composite bulge profile is deprojected in a similar fashion, assuming a constant ellipticity of 0.25 (Burkhead 1986). We then add ν_disk(r, θ) + ν_bulge(r, θ) to obtain the total luminosity density distribution ν(r, θ) input to our models.

The GC luminosity density profile is obtained via a similar deprojection, but with the additional assumption of spherical symmetry. The normalization of the GC light profile is irrelevant, as our models fit only to the slope of the profile.

2.2. Kinematics

Kinematics for NGC 4594 come from four sources. The first uses near-IR data from Gemini/GNIRS long-slit observations along the major axis. These data were taken under good seeing conditions (around 0 farcs 5) and have high signal-to-noise ratio. The second set comes from the FOS on HST using the square aperture of 0 farcs 21 × 0 farcs 21 and is published in K96. The third set of data is from the SAURON instrument (Emsellem et al. 2004). The SAURON data for NGC 4594 have not been published previously. Individual velocities from GCs are our fourth source of kinematics. These data are published in Bridges et al. (2007). We describe each data set in detail.

2.2.1. Gemini Kinematics

We use GNIRS (Elias et al. 2006) on the Gemini South Telescope to measure near-IR spectra of NGC 4594. The data were taken on 2005 January 17. We placed the 150'' × 0 farcs 30 slit along the major axis with the galaxy nucleus centered within the slit. We use a spatial pixel size of 0 farcs 15. With the 32l l/mm grating in third order, we obtain a wavelength coverage of 19800–26200 Å at 6.4 Å pixel⁻¹. Using sky lines, we measure a resolving power around 1700 or an instrumental dispersion of 75 km s⁻¹. The total on-target exposure is 24 minutes, taken in 12 × 2 minute individual exposures. Sky frames of equal exposure are taken throughout.

From both setup images and images of telluric standards, we measure an FWHM in the spatial direction of 0 farcs 5, assuming a Gaussian distribution. We use this PSF for the dynamical models.

We use a custom pipeline to reduce the GNIRS data; however, it produces very similar results to the Gemini GNIRS reduction package. The pipeline includes dark subtraction, wavelength calibration for the individual exposures, sky subtraction, registration, and summing.

There is adequate signal to extract kinematics out to a radius of 45''. Figure 2 shows an example spectrum, where we plot the data in black and the template convolved with the best-fit LOSVD. The velocity templates come from the GNIRS spectral library (Winge et al. 2009), where we select stars with a range of types from G dwarf to late giant. The kinematic extraction program performs a simultaneous fit to the LOSVD and relative weights of the templates. This procedure is described in Gebhardt et al. (2000) and Pinkney et al. (2003). We present these data in the form of Gauss–Hermite moments in Table 2.

Table 2. Gemini Kinematics

R (arcsec)	V ( km s⁻¹)	ΔV ( km s⁻¹)	σ ( km s⁻¹)	Δσ ( km s⁻¹)	h₃	Δh₃	h₄	Δh₄
0.00	19	14	253	16	−0.087	0.033	0.023	0.046
0.15	−36	11	257	12	−0.008	0.042	−0.016	0.038
0.30	−75	12	249	8	−0.053	0.048	−0.017	0.037
0.52	−112	11	234	8	0.020	0.037	−0.047	0.035
0.82	−144	12	221	9	0.060	0.038	−0.051	0.032
1.20	−172	9	202	9	0.065	0.039	0.003	0.031
1.73	−190	8	185	9	0.107	0.031	0.018	0.031
2.40	−208	7	184	9	0.089	0.032	0.024	0.031
3.30	−232	8	175	9	0.140	0.029	0.044	0.027
4.57	−236	9	171	7	0.173	0.032	0.023	0.027
6.45	−235	7	178	9	0.165	0.033	0.058	0.025
8.77	−189	11	203	11	0.028	0.039	−0.007	0.036
11.40	−171	10	192	10	0.071	0.033	−0.009	0.033
14.32	−187	16	185	15	0.227	0.042	0.102	0.060
17.70	−201	13	229	16	−0.023	0.042	0.020	0.050
22.20	−228	12	198	14	0.125	0.040	0.052	0.047
28.58	−235	11	194	14	0.122	0.042	0.036	0.051
36.08	−285	6	141	10	0.046	0.041	0.032	0.034
44.40	−277	7	149	9	0.094	0.042	0.101	0.039
−0.15	45	12	240	15	−0.026	0.034	0.002	0.041
−0.30	112	11	243	15	−0.052	0.036	−0.016	0.038
−0.52	130	10	224	15	−0.080	0.043	0.014	0.037
−0.82	162	8	212	10	−0.068	0.049	0.019	0.030
−1.20	176	7	206	11	−0.066	0.047	0.065	0.031
−1.73	205	7	190	12	−0.087	0.043	0.022	0.037
−2.40	226	7	173	12	−0.101	0.037	0.053	0.039
−3.30	244	9	187	11	−0.113	0.039	0.037	0.037
−4.65	246	8	184	12	−0.130	0.035	0.025	0.040
−6.15	237	9	220	14	−0.133	0.042	0.094	0.035
−8.40	221	9	211	13	−0.191	0.051	0.079	0.036
−10.95	162	9	209	14	−0.093	0.048	0.028	0.043
−14.25	178	9	198	14	−0.092	0.059	0.046	0.042
−19.20	204	10	198	17	−0.166	0.060	0.094	0.048
−24.52	207	11	193	17	−0.045	0.057	0.029	0.052
−29.33	241	11	182	16	−0.149	0.054	0.092	0.050
−36.15	271	11	180	12	−0.084	0.046	0.053	0.042
−45.15	274	10	141	9	−0.053	0.041	−0.029	0.028

Notes. Kinematics along the major axis of NGC 4594. Gauss–Hermite moments were derived from the LOSVDs that are the input to the dynamical models.

Download table as: ASCII Typeset image

Figure 3 shows the kinematics derived from our analysis of the Gemini spectra. Between 1'' and 5'', V rises and σ drops. This is the result of the nuclear disk which becomes important at this radial range (K88). Beyond 10'', we see similar behavior in V and σ, it is caused by the main stellar disk.

2.2.2. HST/FOS Kinematics

K96 present HST/FOS kinematics of the nuclear region of NGC 4594. The FOS has a 0 farcs 21 × 0 farcs 21 aperture. There are three pointings with accurately known positions for NGC 4594 (GO-5512; PI: Faber). The dynamical models include the exact placement and aperture size for the FOS pointing (K96), and use the HST PSF (Gebhardt et al. 2000).

2.2.3. SAURON Kinematics

We also include SAURON integral field kinematics. The SAURON data are from a single pointing exposing on the central region, taken in the low-resolution setting of the instrument (Bacon et al. 2001). In addition to V and σ, the SAURON data also include the higher order Gauss–Hermite moments h₃ and h₄. Details of the data reduction and analysis can be found in Bacon et al. (2001) and Emsellem et al. (2004).

Our modeling code fits to the entire LOSVD rather than its moments, so we reconstruct LOSVDs from the Gauss–Hermite moments. We create 100 Monte Carlo realizations of a non-parametric LOSVD from the uncertainties in the Gauss–Hermite parameters of each SAURON bin (Gebhardt & Thomas 2009). The 1433 reconstructed SAURON LOSVDs are spatially sampled more finely than our modeling bins. We therefore average the SAURON data to match our binning by weighting according to the uncertainties in the LOSVDs.

We reconstruct Gauss–Hermite moments from the combined SAURON LOSVDs for plotting purposes only. Figure 3 shows these moments near the major and minor axes. The major axis V for the SAURON data is significantly lower than that measured for the Gemini data. The reason for this is that SAURON data are binned to match the gridding of our model bins. Near the major axis, the bins range in polar angle from θ = 0° to 11°. These bins are described by a single LOSVD constructed by averaging individual LOSVDs which sample the region at smaller spatial scales. Thus, the average LOSVD contains contributions from LOSVDs as much as θ = 11° above the major axis.

2.2.4. Globular Cluster Kinematics

At large radii, we use individual GC velocities published in Bridges et al. (2007) to derive LOSVDs. The data contain positions and radial velocities for 108 GCs in NGC 4594. We discard the innermost 14 GCs as there are too few GCs inside R lesssim 130'' to reconstruct an LOSVD in the inner parts of the galaxy. Assuming axisymmetry, the positions of the GCs are folded about the minor and major axes. In order to preserve rotation, we flip the sign of the velocity for all GCs that are folded about the minor axis. The GCs are then divided into annular bins extending from θ = 0° to 90° at radii of 131'', 214'', 350'', 574'', and 941'' with roughly 20 GCs per bin.

Within each spatial bin, we calculate the LOSVD from the discrete GC velocities by using an adaptive kernel density estimate adapted from Silverman (1986) and explained in Gebhardt et al. (1996). Each LOSVD contains 15 velocity bins. The velocity bins are highly correlated for the GCs, and there are likely only a few degrees of freedom per LOSVD. The 1σ uncertainties in the LOSVDs are estimated through bootstrap resamplings of the data (Gebhardt et al. 1996; Gebhardt & Thomas 2009).

We compute Gauss–Hermite moments from the GC LOSVDs—again for plotting purposes only—and show these in Figure 3. Uncertainties are calculated by fitting moments to each resampling of the LOSVD during the bootstrap. The GC kinematics resemble the minor axis stellar kinematics in many of the panels. For example, their dispersions appear to be an extrapolation of the minor axis velocity dispersions. There is slight rotation (small h₃) and possible evidence of radial anisotropy (positive h₄).

3. DYNAMICAL MODELS

The dynamical models rely on the orbit superposition technique first developed by Schwarzschild (1979). We assume axisymmetry and match the luminosity density profile and kinematics of the galaxy to those reconstructed from an orbit library. The library is populated with orbits carefully chosen to sample E, L_z, and the third, non-classical integral I₃.

The code used in this paper is described in Gebhardt et al. (2000, 2003), Thomas et al. (2004, 2005), and Siopis et al. (2009). Similar axisymmetric codes are presented in Rix et al. (1997), van der Marel et al. (1998), Cretton et al. (1999), and Valluri et al. (2004). van den Bosch et al. (2008) present a fully triaxial Schwarzschild code. The basic outline of our code is as follows: (1) convert the luminosity density distribution ν(r) into the stellar density ρ(r) via an assumed stellar mass-to-light ratio M/L_I. (2) Add to this density the contribution from a BH of mass M_• and a DM halo with density profile ρ_DM(r). (3) Calculate the potential Φ associated with this density distribution and integrate a large number of orbits (typically ~ 20,000) over many dynamical times. (4) Assign a weight w_i to each orbit and determine the w_i values by minimizing the χ² difference between the observed kinematics and luminosity density of the galaxy and those resulting from the PSF-convolved orbit library, subject also to the constraint of maximum entropy.

We maximize the entropy-like quantity $\hat{S} \equiv S-\alpha \chi ^2$ , where S is the Boltzmann entropy and α controls the relative weight of S or χ². For small values of α, reproducing the observed kinematics becomes unimportant, and the models act to only maximize entropy. As α increases, maximizing entropy becomes less important and more weight is given to matching the observations. In practice, we start with a small value of α and gradually increase it until χ² asymptotes. The interested reader may see Siopis et al. (2009) or Shen & Gebhardt (2010) for more details.

Our model grid consists of 19 radial and 5 azimuthal bins covering a radial range of 0 farcs 03–1800'' spaced logarithmically. Additionally, we use 15 velocity bins to describe our LOSVDs. We incorporate the effects of seeing by convolving the light distribution for each orbit with a model PSF before comparing with data (Gebhardt et al. 2000). We approximate the PSF as Gaussian with an FWHM of either 0 farcs 94, 0 farcs 5, or 0 farcs 09 depending on whether the data are from SAURON, Gemini, or HST observations, respectively. The convolution extends to a radius of 10 × FWHM.

We run over 8500 models with different values of the model parameters M/L_I, M_•, and ρ_DM. We use Δχ² statistics to determine the best-fit parameter values and their uncertainties. Models whose values of χ² are within Δχ² = 1 of the minimum for a given model parameter (marginalized over the others) define the 1σ or 68% confidence band of that parameter.

3.1. Model Assumptions

Our fiducial density profile is a combination of stellar mass, DM, and a central SMBH

$\begin{equation} \rho (r,\theta) = \frac{M}{L} \nu (r,\theta) + \rho _{{\rm DM}}(r,\theta) + M_{\bullet } \delta (r), \end{equation} \tag{ 2 }$

where M/L_I is the stellar mass-to-light ratio, assumed constant with radius and δ(r) is the Dirac delta function. The angle θ is the angle above the major axis. While ρ_DM can in principle be a function of θ, we do not consider flattened models. We assume a spherically symmetric, logarithmic halo of the form

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

This profile is cored for radii r lesssim r_c and produces a flat rotation curve with circular speed V_c for r r_c. It has two free parameters, r_c and V_c, which are varied in the fitting process. Including M/L_I and M_•, this brings the total number of model parameters to four.

Recently, Gebhardt & Thomas (2009) have shown that the inclusion of a DM halo can significantly affect modeled BH masses. To test for this, we run a smaller suite of models without a dark halo.

4. RESULTS

Our best-fit values for the four model parameters are $M/L_I=3.4 \pm 0.05 \, \frac{M_{\odot }}{L_{\odot }}$ , M_• = (6.6 ± 0.4) × 10⁸ M_☉, V_c = 376 ± 12 km s⁻¹, and r_c = 4.7 ± 0.6 kpc. Figure 4 shows the χ² minima around each of the model parameters. Each dot represents a single model and the solid curve is a smoothed fit to the minimum. The points of the solid curve at Δχ² = 1 above the minimum determine the 1σ confidence limits on the parameters. All four model parameters have well behaved χ² curves with sharp, well-defined minima. This allows robust determination of the model parameters with small 1σ uncertainties.

The uncertainties we present are derived strictly from Δχ² statistics. The 1σ error bars on quoted parameters correspond to models within Δχ² = 1 of the minimum value. Systematic effects are likely to contribute in addition to this quoted uncertainty. While in general for other galaxies one of the biggest sources of systematic uncertainty is inclination, for NGC 4594 inclination uncertainties are unimportant. Other sources of uncertainty may include effects due to non-axisymmetries, but these are likely small or zero since only the most massive ellipticals are thought to be significantly triaxial (Binney 1978; Kormendy & Illingworth 1982; Tremblay & Merritt 1996). For more on systematic uncertainties, the reader is referred to Gebhardt et al. (2003) and Gültekin et al. (2009a, 2009b).

Figure 5 shows correlations among the four model parameters. Plotted are the different projections of the four-dimensional parameter space; every small dot corresponds to a model run. Red dots are models that lie within Δχ² = 4 of the minimum, and large black dots are within Δχ² = 1. There appears to be a slight correlation between M/L_I and M_•—much less severe than in M87 (Gebhardt & Thomas 2009). Not surprisingly, the high resolution of our HST kinematics is able to break the degeneracy between M_• and M/L_I. We discuss this further below. The dark halo parameters do not show any obvious correlation, indicating the GC and stellar kinematics were able to break the degeneracy usually observed between these two parameters.

Our best-fit model has (unreduced) χ² = 582.6. It is non-trivial to calculate the number of degrees of freedom ν_DOF. Roughly, ν_DOF = N_LOSVD × N_bin; however, there are complicated correlations between velocity bins (Gebhardt et al. 2003). With this crude estimate for ν_DOF, our best-fit model has reduced χ²_ν = 0.6.

We compare the modeled value of our stellar mass-to-light ratio with that obtained from evolutionary population synthesis models (Maraston 1998, 2005). We adopt values of 10 Gyr and 0.1 for the stellar age and metallicity of NGC 4594 (Sánchez-Blázquez et al. 2006) and use these to derive the predicted I-band M/L_I from the Maraston models. For a Salpeter initial mass function (IMF) with stellar masses drawn from the range 0.1–100 M_☉, this analysis yields M/L_I = 3.99. If instead the stars obey a Kroupa IMF drawn from the same range, then M/L_I = 2.58. We multiply these by a factor of 1.096 corresponding to A_I = 0.099 to correct for Galactic extinction along the line of sight (Schlegel et al. 1998) to obtain M/L_I = 2.83 and M/L_I = 4.37. Our dynamically derived stellar M/L_I = 3.4 ± 0.05 falls nicely between these two values.

In Figure 6, we plot the total mass-to-light ratio as a function of radius for our best-fit model with 1σ uncertainties (gray region). The red cross-hatched region represents the range in stellar M/L_I from stellar population models described above. Total M/L_I rises near the center of the galaxy due to the contribution of the SMBH. As we go out in radius, the stars become more important to the total mass over roughly the range 5''–50''. Here the total M/L_I approaches both our dynamically determined M/L_I and the range derived from stellar population models. Past 50'' M/L_I once again rises due to the importance of the dark halo.

Figure 7 plots the enclosed mass of each component as well as the total mass of the galaxy. At our innermost bin, the total mass is almost two orders of magnitude greater than the stellar mass, meaning we are probing the BH's sphere of influence quite well. The green line plotted is the mass profile Kormendy & Westpfahl (1989) derived from their gas rotation curve. It agrees well with the total mass distribution derived here.

4.1. Models without Dark Matter

We run 189 models with no dark halo. In these models, we exclude the GC data and use only the stellar kinematics. Correspondingly, the number of degrees of freedom impacting the unreduced χ² are proportionately fewer. We measure a BH mass of M_• = (6.6 ± 0.3) × 10⁸ M_☉ and stellar mass-to-light ratio of M/L_I = 3.7 ± 0.05. The minimum unreduced χ² = 628, proving models without a dark halo are a worse fit.

We do not see the dramatic change that Gebhardt & Thomas (2009) see in M87 where the inclusion of a DM halo causes their determination of M_• to double. Instead, our results mirror those of Shen & Gebhardt (2010) in NGC 4649 where the inclusion of a DM halo does not significantly change the modeled M_•. The likely explanation for this behavior is the inclusion of high-resolution HST kinematics in both NGC 4649 and NGC 4594. Schulze & Gebhardt (2011) find the same effect for a larger sample of galaxies. Whenever the data have high enough resolution to resolve the BH's radius of influence R_inf ~ GM_•/σ², DM has no significant effect on the determination of M_•. For NGC 4594 we measure R_inf sime 57 pc sime 1 farcs 2. We use HST/FOS kinematics whose central pointing has a PSF of 0 farcs 09 sime 0.08 R_inf. Additionally, the light profile of NGC 4594 is more centrally concentrated than that of M87. These factors combine to allow a more accurate determination of M_•, removing the freedom that the models have to trade mass between M_• and M/L_I. This is evidenced by the lack of correlation among M_• and M/L_I in Figure 5.

4.2. Orbit Structure

Having already determined the orbital weights that provide the best fit to the data, we reconstruct the internal unprojected moments of the distribution function. We perform this analysis on our best-fit model and the models that define the 68% confidence region (over all combinations of the four model parameters), yielding internal moments at each grid cell.

We define the tangential velocity dispersion to be $\sigma _t \equiv \scriptstyle\sqrt{\frac{1}{2}(\sigma _{\phi }^2+\sigma _{\theta }^2)}$ , where σ is actually the second moment, containing contributions from both streaming and random motion in the phgr direction. Figure 8 shows the radial run of the ratio σ_r/σ_t. The second moment of the distribution function is tangentially biased where the disk is important (gray region) as expected but is mostly isotropic elsewhere. The red region plots σ_r/σ_t for stars near the minor axis, showing almost perfect isotropy. The green region indicates that at large radii, GC kinematics show significant radial anisotropy. We discuss the implications of this in Section 5.2 below.

**Figure 8.** Radial run of the ratio of the radial to tangential components of the velocity dispersion tensor. Shaded ares represent 68% confidence regions, with gray indicating stars near the major axis, red meaning stars near the minor axis, and green representing GCs averaged over all angles.
Download figure:
Standard image High-resolution image

5. DISCUSSION

5.1. Black Hole–Bulge Correlations

We discuss the position of NGC 4594 on the M–σ and M–L relations, as defined by Gültekin et al. (2009b, hereafter G09) and compare our values of the correlation parameters to previous measurements. We calculate the effective velocity dispersion σ_e similar to G09

$\begin{equation} \sigma _e^2 \equiv \frac{{\int }^{R_e}_{R_{\mathrm{inf}}} (V^2(R) + \sigma ^2(R)) I(R) dR}{{\int }^{R_e}_{R_{\mathrm{inf}}} I(R) {\it dR}}, \end{equation} \tag{ 4 }$

where V(R) is the rotational velocity and I(R) is the surface brightness profile. This makes σ_e essentially the surface-brightness-weighted second moment. Instead of integrating from the center of the galaxy (R = 0) as G09 did, we integrate from R > R_inf to ensure that we do not bias σ_e with the high dispersion near the BH. Our outermost kinematic data point is at R = 45'', thus we have a gap in kinematic coverage between 45'' < R < R_e = 114''. The velocity dispersion near the end of our long-slit data is dropping sharply; however, V may still contribute to the integral for R > 45''. To investigate this, we use the gas rotation curve presented in Kormendy & Westpfahl (1989) which extends well beyond R_e. Truncating the integral at R = 45'' gives σ_e = 292 km s⁻¹ while using the extended rotation curve yields σ_e = 297 km s⁻¹.

The problem with this definition of σ_e is that it includes a contribution from the rotation of the disk. It has been shown that BH mass does not correlate with disk properties (Kormendy et al. 2011) so this is not ideal. However, to compare with G09 we must be consistent in our calculation of σ_e. We therefore quote this value of σ_e when we compare to the M–σ relation determined by G09. As we expect BH mass to track bulge quantities, disk contribution to σ_e is likely to add a source of intrinsic scatter to spiral galaxies in the M–σ relation. In fact, spiral galaxies are observed to have larger scatter about M–σ than ellipticals of similar σ_e.

We also discuss some possible alternatives to σ_e where we attempt to remove the disk contribution. One option is to remove V²(R) from Equation (4) altogether. Bulges are known to rotate, however (Kormendy & Illingworth 1982), and this will likely underestimate σ_e. This crude calculation gives σ_e = 200 km s⁻¹.

Another option is to assume some degree of bulge rotation a priori. If we assume NGC 4594 rotates isotropically (Kormendy & Illingworth 1982), then its flattening determines its position on the V/σ– epsilon diagram (Binney 1978). Kormendy (1982) shows that the relation

$\begin{equation} \frac{V}{\sigma } \approx \sqrt{\frac{\epsilon }{1-\epsilon }} \end{equation} \tag{ 5 }$

approximates the isotropic rotator line to roughly 1% accuracy. We use our value of the bulge ellipticity epsilon = 0.25 in Equation (5) and assume this value of V/σ applies globally to the entire bulge. We then use our measured dispersion profile σ(R) to determine the bulge velocity V_bulge(R). Using these quantities, we determine σ_e = 230 km s⁻¹. We compare this to the kinematics listed in Kormendy & Illingworth (1982). These data include long-slit spectra taken at a position angle parallel to the major axis, but 30'', 40'', and 50'' above it. From these data, it is apparent that the bulge σ off the major axis is roughly constant at ~220 km s⁻¹. The rotation velocity rises from 0 to 100 km s⁻¹ at large radii. We estimate the luminosity-weighted mean V ~ 50 km s⁻¹. Adding this in quadrature to the constant bulge σ = 220 km s⁻¹ gives σ_e ≈ 226 km s⁻¹. This estimate does not contain any rotation from the disk and is consistent with our determination of σ_e = 230 km s⁻¹ obtained by assuming a constant V/σ.

Our BH mass M_• = (6.6 ± 0.4) × 10⁸ M_☉ agrees nicely with that of G09 which uses the K88 value. In fact, when corrected for their various distance determinations, most values of M_• in the literature agree quite well (K88; Emsellem et al. 1994b, K96; Magorrian et al. 1998) despite the many modeling techniques and data sets used. This is likely due to the high degree of isotropy as evidenced in Figure 8, deduced from the V/σ– epsilon diagram (Kormendy & Illingworth 1982), and noted in K88.

Figure 9 plots the position of NGC 4594 on the G09 M–σ and M–L relations. We plot each determination of σ_e in the left-hand panel. Straightforward application of Equation (5) leads to a value of σ_e that falls directly on the G09 M–σ line (blue asterisk). Next closest is the method of calculating σ_e by assuming a value of V/σ (green square). This point lies 0.44 dex above the G09 line; however, this is still within the estimated scatter. The calculation of σ_e that ignored all rotation is, not surprisingly, farthest from the G09 line. Calculation of the relevant quantities for comparison with the M–L relation is straightforward, and we plot our value of M_• and L_V (green square) along with that from G09 in the right-hand panel.

5.2. Globular Clusters

As demonstrated in Section 4.2, we find significant radial anisotropy in the GCs. It is interesting that the stellar kinematics at smaller radii do not show this feature. This difference in orbital properties combined with the difference in their light profiles might suggest the GCs and stars are two distinct populations of tracer particles. This could also indicate the two populations have different formation scenarios.

Unfortunately, there is no radius in the galaxy where we have simultaneous coverage of both stellar and GC kinematics. Thus, we are unable to test whether the stellar orbits become more radial in the ~50'' between where the stellar kinematics run out and the GCs begin. However, since the light profiles of both populations are significantly different (Figure 1) there is no reason to assume they should share similar orbit properties.

Figure 8 shows the GCs in NGC 4594 are radially anisotropic (σ_r/σ_t > 1) over roughly the radial range 100''–1000'' (approximately 1–10 R_e). Previous studies of the GC systems of galaxies have found their velocity ellipsoids to be isotropic (Côté et al. 2001, 2003). However, these studies used spherical Jeans modeling instead of the more general axisymmetric Schwarzschild code we use.

Rhode & Zepf (2004) determine with high confidence that the color distribution of the GC system in NGC 4594 is bimodal. This may indicate different subpopulations of GCs with different orbital properties that formed at different epochs in the galaxy's history. In our analysis, we make no distinction between red and blue subpopulations. We use the light profile and kinematics of all available GCs, regardless of color. However, since we use different sources for our kinematics and photometry data, there is the possibility that each source draws from a different GC subpopulation.

5.3. Dark Halo

The parameters V_c and r_c of our model dark halo imply a central DM density of ρ_c = 0.35 ± 0.1 M_☉pc⁻³ Using an improved Jeans modeling technique, Tempel & Tenjes (2006) model NGC 4594 and find a DM halo with central density ρ_c = 0.033 M_☉pc⁻³, 10 times lower than our value. They, however, measure a larger stellar M/L_I_V = 7.1 ± 1.4 in the bulge.

We plot the fraction of enclosed mass that is DM as a function of half-light radius R_e in Figure 10. At 1 R_e there is already a roughly 50–50 mix of stars and DM. Inside of R_e the DM still contributes a non-negligible fraction to the total mass content.

sstarf — **Figure 10.** Fraction M_DM/(M + M_DM) of enclosed mass that is dark matter as a function of radius.
Download figure:
Standard image High-resolution image

In a study measuring DM properties in 1.7 × 10⁵ local (z < 0.33) early-type galaxies from the Sloan Digital Sky Survey, Grillo (2010) find a correlation between the fraction of DM within R_e and the logarithmic value of R_e. With our measured value of R_e, this correlation predicts a DM fraction at R_e of 0.68. Our value of 0.52 is smaller, but still within their 68% confidence limit.

Thomas et al. (2009) derive scaling relations for halo parameters based on observations of early-type galaxies in the Coma cluster. These relations are constructed for similar galaxies using the same halo parameterization and modeling code used in this paper. This makes comparison to our parameters straightforward. We compare to the observed relations between halo parameters r_c, V_c, and ρ_c and total blue luminosity L_B. Our value of V_c falls directly on the V_c–L_B relation; however, our measured r_c is smaller by roughly an order of magnitude. Since ρ_c∝V²_c/r²_c, the discrepancy in r_c causes our measurement of ρ_c to be high when compared to the Thomas et al. (2009) ρ_c–L_B relation. Scatter in this relation is large, however, and the environment of NGC 4594 is different from that of the Coma galaxies.

Kormendy & Freeman (2004; 2011, in preparation) also derive scaling laws for similar parameters in galaxies of later Hubble type (Sc-Im). We measure a much higher density and much smaller core radius than these relations imply at the L_B of NGC 4594. We interpret this as the result of severe compression of the halo by the gravity of the baryons (Blumenthal et al. 1986). Such an effect is expected in early-type galaxies with massive bulges.

We thank Eric Emsellem for providing reduced SAURON data and helpful comments. This work would not be feasible without the excellent resources of the Texas Advanced Computing Center (TACC). K.G. acknowledges support from NSF-0908639. D.R. is grateful for hospitality and support from the Institute for Advanced Study in the form of a Corning Glass Works Foundation Fellowship.

ORBIT-BASED DYNAMICAL MODELS OF THE SOMBRERO GALAXY (NGC 4594)

ABSTRACT