format: text

text:
# Gas Kinematics and the Black Hole Mass at the Center of the Radio Galaxy NGC 4335Based on observations with the NASA/ESA Hubble Space Telescope obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

## 1 Introduction

The Hubble Space Telescope has made it possible to measure the masses of black holes (BHs) in the centers of many nearby active and quiescent galaxies using stellar and/or gaseous kinematics (for reviews, see e.g., Richstone et al. 1998; van der Marel 1999; Ho 1999; de Zeeuw 2001; Kormendy & Gebhardt 2001). To date, BH masses have been measured in about 40 galaxies, both spirals and ellipticals, and this number continues to increase. The BH masses correlate loosely with host spheroid luminosity (Kormendy & Richstone 1995) and more tightly with inner stellar velocity dispersion (Gebhardt et al. 2000; Ferrarese & Merritt 2000).

Central emission-line gas is detected in virtually all nearby radio galaxies, defined here as galaxies which harbor kpc-scale radio-jets. The equivalent widths of the gas emission lines are generally much larger than those of the absorption lines in the integrated stellar light, so that measurement of the kinematics of nuclear emission-line gas is an efficient way to determine the central gravitational potential and BH mass (e.g., Harms et al. 1994; Ferrarese, Ford & Jaffe 1996; Macchetto et al. 1997; van der Marel & van den Bosch 1998; Ferrarese & Ford 1999; Verdoes Kleijn et al. 2000, hereafter VK00; Sarzi et al. 2001; Barth et al. 2001). A drawback of this approach is that the gas kinematics might be affected by non-gravitational motions. Nevertheless, using the gas kinematics to determine accurate BH masses is particularly interesting for radio galaxies because determining the BH mass and the properties of the gas disk can advance our understanding of radio-jet formation and evolution. In turn, the radio-jet offers an extra diagnostic of the BH accretion and immediate BH surroundings. Current questions include: what is the lower-limit to BH masses that can form kpc-scale jets? Is this lower-limit higher or lower than typical BH masses in spirals which (almost) never show kpc-scale jets? Is there a correlation between BH mass and jet properties, such as total power or jet velocity? For these and other reasons we are performing a systematic study of a complete sample of nearby radio galaxies with Fanaroff & Riley (1974) Type I radio jets, the ‘UGC FR-I sample’ (Verdoes Kleijn et al. 1999, hereafter VK99; Xu et al. 2000) using observations at multiple wavelengths. In particular, we have performed a STIS spectroscopic survey of the inner gas distributions to measure the kinematics and the physical state of the gas (Noel-Storr et al. 2002).

Here we concentrate on one galaxy from our sample, NGC 4335, which is a relatively unknown isolated giant elliptical ($`M_B=20.7^m`$; Paturel et al. 1997) at 66 Mpc. The gas appears embedded in a well-defined dust disk (diameter $`750\mathrm{pc}`$) and gas kinematics can be traced sufficiently far out along the three slits to allow detailed gas dynamical modeling. In addition we perform stellar dynamical modeling for NGC 4335 using the stellar kinematics derived from a WHT/ISIS long-slit observation.

The paper layout is as follows. Sections 2 and 3 present HST/WFPC2 broad- and narrowband imaging, HST/STIS gas emission-line spectroscopy and ground-based WHT/ISIS stellar absorption-line spectroscopy, including the basic data reduction and derivation of the gaseous and stellar kinematics. Section 4 describes the modeling of the gas disk flux distribution, the derivation of the stellar mass distribution and the fits to the observed gaseous kinematics to estimate the BH mass. Section 5 describes two-integral modeling of the WHT stellar kinematics to determine the stellar mass-to-light ratio and to constrain the BH mass independently. Section 6 discusses the implications of the BH mass measurements in NGC 4335 for our understanding of BH demography, and for the techniques used to measure black hole masses.

We adopt $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ throughout this paper. This does not directly influence the data-model comparison for any of our models, but does set the length, mass and luminosity scales of the models in physical units. Specifically, distances, lengths and masses scale as $`H_0^1`$, while mass-to-light ratios scale as $`H_0`$.

## 2 Imaging

### 2.1 WFPC2 Setup and Data Reduction

We imaged NGC 4335 in the context of HST program GO-6673. We used the WFPC2 instrument (described in, e.g., Trauger et al. 1994; Biretta et al. 1996) on March 15, 1997 to obtain images of NGC 4335 in the F555W ($`V`$) and F814W ($`I`$) filters and linear ramp filter FR680N which includes H$`\alpha `$+\[N II\]. The observing log is presented in Table 1. The ‘Linear Ramp Filters’ (LRFs) of the WFPC2 have a central wavelength that varies as a function of position on the detector. The LRF FR680N image (WF CCD) was used as ‘on-band’ image, with the galaxy position chosen so as to center the filter transmission on the H$`\alpha `$+\[N II\] emission lines. A combination of the broad-band images F555W and F814W (PC CCDs) was chosen as ‘off-band’ image and covers primarily stellar continuum. This combination of filters ensures that the effect of dust extinction is the same in the on- and off-band image. The final H$`\alpha `$+\[N II\] emission image was obtained by subtracting the stellar continuum, as determined by the off-band image, from the on-band image. The quality of the H$`\alpha `$+\[N II\] image is set by the Poisson noise of input images and the uncertainty in the alignment and scaling of the on-band image and off-band image. The stellar surface brightness is steeply rising at the center of NGC 4335. The resulting emission image depends quite sensitively on precise sub-pixel alignment and scaling. This introduces a $`15\%`$ flux error which is estimated by varying the alignment and scaling within a reasonable range. The $`V`$-band and H$`\alpha `$+\[N II\] images are shown in Figure 1. Extended H$`\alpha `$+\[N II\] emission is present with a roughly elliptical morphology, elongated in the same direction as the dust disk. As described in Section 4.1, the emission-line gas is consistent with being in a thin disk, but the uncertainties involved in constructing the image are too large to prove this conclusively. We cannot rule out a more spheroidal distribution, especially towards the nucleus.

These data were presented previously in VK99, who also give a more detailed description of the data reduction.

### 2.2 The Dust Disk

The central dust detected in NGC 4335 extends to more than 2 kpc from the nucleus. In the outer regions ‘arms’ of dust gradually align with the well-defined central dust disk which appears as roughly half an ellipse on the west side of the nucleus (see Fig. 1). This is commonly seen in galaxies with inclined nuclear dust disks (e.g., NGC 7052, van der Marel & van den Bosch 1998; NGC 6251, Ferrarese & Ford 1999; 3C 449, Martel et al. 1999). On the east side no dust obscuration is evident in either the $`V`$\- or $`I`$-band image and its $`VI`$ color is typical for dust-free giant elliptical galaxies. We conclude that the eastern half of the dust disk produces negligible obscuration in the optical.

We fitted by eye an elliptical contour to the outline of the dust disk in the $`V`$-band image, and measured the position angle of the major axis of the dust disk to be PA=$`156^{}`$. This is well-aligned with the galaxy major axis just outside the dust disk (cf. Fig. 1 and VK99). If we assume that the dust disk is thin and intrinsically circular, then the inclination is $`i=66^{}\pm 7^{}`$ (VK99).

The total mass of the kpc-scale dust distribution estimated using the observed color excess is $`5.8\times 10^5\mathrm{M}_{}`$ if the dust resides in the mid-plane of the galaxy (VK99). Following the prescription by Goudfrooij & de Jong (1994), the IRAS $`60\mu `$ and $`100\mu `$ flux measurements for NGC 4335 from Knapp et al. (1989) indicate a dust mass $`2.2\times 10^6\mathrm{M}_{}`$. This is not inconsistent with the aforementioned result, given that the IRAS measurements cover a much larger area than that occupied by the nuclear dust disk studied with HST. Furthermore, dust mass estimates using color excess are not sensitive to a smoothly distributed dust component, if present (cf. e.g., Tran et al. 2001).

## 3 Spectroscopy

### 3.1 HST/STIS Gas Emission-line Spectroscopy

In the context of GO-program 8236, we used STIS (see Kimble et al. 1998) on September 22, 1999 to obtain long slit spectra of NGC 4335 during two orbits of the HST telescope. We obtained spectra along three adjacent slit positions using a $`0.2^{\prime \prime }`$ wide slit (52X0.2). The lay-out of the slits is shown in Figure 1. The G750M grating was used in combination with two pixel spatial and spectral on-chip rebinning (in order to obtain enough signal to noise while retaining sufficient spatial and spectral resolution) yielding spectra covering the wavelength range from 6436 Å to 7100 Å over 511 pixels. Wavelength calibration spectra of the internal arc lamp were obtained at the beginning of the first orbit and at the conclusion of the second orbit.

The slits were aligned along PA=$`32.44^{}`$. This is almost parallel ($`\mathrm{\Delta }`$PA=$`8.44^{}`$) to the galaxy major axis (PA=$`156^{}`$). We denote the central slit by C, and the adjacent slits on the east- and west-side by E and W, respectively. A log of the observations is provided in Table 2. Target acquisition uncertainties and other possible systematic effects could cause the slit positions on the galaxy to differ slightly from those commanded to the telescope. We determined the actual slit positions from the data themselves, using the STIS continuum and emission-line fluxes and the WFPC2 images. We compared the continuum counts as function of position along the slit for the three STIS slits with the counts in the WFPC2 $`V`$ and $`I`$ images at these positions. The STIS continuum was determined blue-ward of H$`\alpha `$+\[N II\] and red-ward of \[S II\] respectively. This analysis indicates that the positioning errors for the STIS slits with respect to the WFPC2 image are less than $`0.05^{\prime \prime }`$ in both the direction along and perpendicular to the slit. Another estimate of the positioning of the slits with respect to the galaxy nucleus was obtained by interpolation of the central few H$`\alpha `$+\[N II\] fluxes assuming the maximum flux occurs at the galaxy nucleus. This analysis implies that the central slit was properly centered on the nucleus to an accuracy of $`0.01^{\prime \prime }`$. The small slews performed by the telescope between observations to dither the slits and observe the adjacent slits have a nominal positioning error of $`0.003^{\prime \prime }`$.

Most of the necessary data reduction steps were performed by the HST/STIS calibration pipeline (CALSTIS version 2.3), including flat-fielding, hot pixel and cosmic-ray removal, absolute sensitivity calibration and wavelength calibration which takes into account the Earth’s motion. To facilitate hot pixel and cosmic-ray removal two exposures were taken at each slit position dithered by 2 pixels in the direction along the slit. The uncertainty in the relative wavelength accuracy is $`0.1`$Å or $`5\mathrm{km}\mathrm{s}^1`$ (Diaz-Miller 2001).

### 3.2 Gas Kinematics and Fluxes

The spectra show several emission lines, of which the H$`\alpha `$+\[N II\] complex at 6548, 6563, 6583 Å has sufficiently high signal-to-noise ratio $`S/N`$ for a kinematical analysis. The \[S II\] doublet at 6716 and 6731 Å can be fitted only in the central $`0.2^{\prime \prime }`$ in the central slit and at a lower signal-to-noise ratio. The results for these few points are consistent with the H$`\alpha `$+\[N II\] kinematics described below. To quantify the H$`\alpha `$+\[N II\] gas kinematics we fitted the spectra assuming that each emission line is a Gaussian on top of a flat continuum. This yields for each line the total flux, the mean velocity $`V`$ and the velocity dispersion $`\sigma `$. We fitted the \[N II\] doublet under the assumption that the individual lines have the same $`V`$ and $`\sigma `$ and that the ratio of their fluxes equals the ratio of their transition probabilities (i.e., 3; see references in Mendoza 1983). The H$`\alpha `$+\[N II\] complex is influenced by blending of the lines, and we made the additional assumptions that H$`\alpha `$ and the \[N II\] doublet have the same kinematics.

We did the Gaussian fitting using software described in van der Marel & van den Bosch (1998). Prior to the emission-line fit the spectra were rebinned along the slit where necessary to obtain sufficient $`S/N`$, and hence sufficiently small errors on the gas kinematics (Table 3), to discriminate between the black hole masses in dynamical models discussed in Section 4.3 and later. The same gaussian fitting was performed on modeled emission-line spectra which resulted from the dynamical models. Differences between the velocity profiles (VP) and the single Gaussian fits are discussed in Section 4.3.

The Gaussian fits to the emission lines are shown in Figure 2. The resulting kinematics are shown in Figure 3 and tabulated in Table 3. The systemic velocity used in Figure 3 and Table 3 was estimated from the HST data themselves, by including it as a free parameter in the dynamical models described below (see Section 4.3). This yields $`v_{\mathrm{sys}}=4620\pm 13\mathrm{km}\mathrm{s}^1`$. This result is consistent with $`v_{\mathrm{sys}}=4615\pm 22\mathrm{km}\mathrm{s}^1`$ by Huchra et al. (1983) and $`v_{\mathrm{sys}}=4595\pm 47\mathrm{km}\mathrm{s}^1`$ from the RC3 catalog (de Vaucouleurs et al. 1991).

The mean-velocity profiles are smooth functions of position and are indicative of a rotating circular disk of emission gas. The rotating disk picture is further supported by the fact that the mean velocities on the minor axis are close to $`0\mathrm{km}\mathrm{s}^1`$ in all three slits. The velocity dispersion profiles in all three slits show irregularities. We verified that this result does not depend sensitively on chosen region to fit the continuum and masked deviant pixels during the fit. The velocity dispersion peaks at $`240\mathrm{km}\mathrm{s}^1`$ at the nucleus in the central slit C. In slit W the dispersion peaks near the nucleus. However, in slit E the dispersion actually dips at the position closest to the nucleus. It seems therefore that the kinematics of NGC 4335 are consistent with a rotating gas disk, but with additional turbulent motion.

We also performed single Gaussian fits to H$`\alpha `$ and the \[N II\] doublet independently. The two resulting flux profiles have very similar shapes. The mean velocities for the two sets of emission lines are very similar as well. The median difference between the mean velocities of the two components is $`20\mathrm{km}\mathrm{s}^1`$. This is slightly larger than the typical formal errors given in Table 3. These formal errors are derived by propagating the spectral flux errors. We will assume $`20\mathrm{km}\mathrm{s}^1`$ as a more realistic error estimate for the H$`\alpha `$+\[N II\] mean velocities. The H$`\alpha `$ and \[N II\] velocity dispersion profiles show very similar behavior, but H$`\alpha `$ has systematically lower velocity dispersions by about $`2070\mathrm{km}\mathrm{s}^1`$. The differences in dispersion could be caused by the fact that (i) the H$`\alpha `$ and \[N II\] have intrinsically different kinematics or (ii) the H$`\alpha `$ kinematic measurements are affected by stellar H$`\alpha `$ absorption. A similar effect and of the same magnitude is seen for the gas disk in the S0 galaxy NGC 3245 (Barth et al. 2001).

### 3.3 WHT/ISIS Stellar Absorption-line Spectroscopy

A long-slit stellar absorption spectrum at PA$`=148^{}`$ (i.e., within $`10^{}`$ from the major axis of NGC 4335) was obtained in service mode with the ISIS (blue arm) spectrograph on the 4.2m William Herschel Telescope on La Palma on June 13, 2001. A slit-width of $`1.0^{\prime \prime }`$ was used with the R300B grating and the EEV12 CCD detector. The dispersion is 0.86 Å/pix and the spatial scale is $`0.19^{\prime \prime }`$/pix. The observed wavelength range was $``$3430–6900 Å. The observations included two 1800s exposures of NGC 4335, a 60s exposure of the B2IVp star BD+33 2642 to facilitate flux calibration, and standard calibration exposures. The instrumental velocity dispersion $`\sigma _{\mathrm{instr}}75\mathrm{km}\mathrm{s}^1`$ was determined from sky lines and the emission-line widths in the calibration lamp spectra. Both the star and galaxy spectra were observed with a seeing of FWHM $``$$`1.2^{\prime \prime }`$.

The data reduction was performed in IRAF using the CCDPROC (Version Dec93) and TWODSPEC packages. The basic data reduction steps include bias subtraction, flat-fielding with dome-flats and cosmic-ray removal. The contribution from dark current is negligible. The wavelength calibration was done with contemporary CuAr lamp spectra taken directly before and after the observations of NGC 4335. Either 31 or 32 arc-lamp emission lines were used to map the wavelength as a function of CCD coordinate. The final residuals per fitted emission line are $`0.3`$ Å (or $`18\mathrm{km}\mathrm{s}^1`$ at 5000 Å). We subsequently performed the wavelength calibration and spatial rectification to align the wavelength direction with the CCD columns in one step. The flux calibration took into account telescope throughput and atmospheric extinction. The Galactic extinction towards NGC 4335 is $`A_B=0.063`$ (Schlegel et al. 1998) and hence negligible. As a cross-check for the flux calibration we verified that the broad-band flux observed through a $`1^{\prime \prime }\times 1^{\prime \prime }`$ central aperture in the WHT/ISIS spectrum and our previously observed HST $`V`$-band image (VK99) agreed to within $`25\%`$, correcting for differences in filter passbands and point spread function (PSF). This level of agreement is satisfactory given the uncertainties in WHT/ISIS seeing conditions and extinction corrections. Finally, after verifying that the final calibrated spectra of stars and galaxies did not have any remaining large-scale patterns, we logarithmically rebinned in the wavelength direction to perform the kinematical analysis.

### 3.4 Stellar Kinematics

We used the wavelength range $`5200`$–5700Å in the NGC 4335 absorption spectrum to derive the stellar kinematics. This range includes the red-shifted Mg $`b`$, Fe1 and Fe2 absorption features. To obtain the stellar kinematics we used the pixel space fitting method of van der Marel (1994), which compares the NGC 4335 spectrum to a stellar ‘template’ spectrum. A suitable template star (spectral type K0V) was obtained during a different observing run with WHT/ISIS in identical setup. To check for template mismatch we also used stellar templates (spectral types K1III and K2) described in van der Marel & Franx (1993). These data have an instrumental velocity dispersion $`\sigma _{\mathrm{instr}}60\mathrm{km}\mathrm{s}^1`$, i.e., very similar to the $`\sigma _{\mathrm{instr}}`$ of the NGC 4335 spectrum. The inferred stellar kinematics agree within the errors. The observed stellar kinematics as a function of radius are listed in Table 4 (and also shown in Figure 12 which is discussed in Section 5). Stars and gas rotate in the same direction in the core of NGC 4335. The galaxy spectrum integrated along the slit over the central $`25^{\prime \prime }`$ yields an aperture velocity dispersion $`\sigma _0=282\pm 15\mathrm{km}\mathrm{s}^1`$. The $`25^{\prime \prime }`$ length corresponds to two effective radii as determined by fitting the surface brightness distribution shown in Figure 6. The listed formal errors include photon noise, uncertainties in the continuum fitting and the variation in derived $`\sigma _0`$ observed by using different stellar templates. We are not aware of previous determinations of the central velocity dispersion for NGC 4335, but note that the central value $`\sigma _0=300\pm 60\mathrm{km}\mathrm{s}^1`$ inferred from the Faber & Jackson (1976) relation agrees well with the observations.

## 4 Modeling and Interpretation of the H$`\alpha `$+\[N II\] kinematics

### 4.1 Flux Distribution

To model the H$`\alpha `$+\[N II\] gas kinematics we need a description of the intrinsic (i.e., deconvolved and de-inclined) emission-line flux distribution. We use fluxes from both the WFPC2 H$`\alpha `$+\[N II\] image and STIS spectra to constrain the intrinsic distribution. Figure 4 compares the actual flux distribution observed through the STIS slits with the expected flux distribution through the STIS slits as inferred from the WFPC2 image. The agreement for slits C and W is good considering the uncertainties in the absolute flux calibrations (especially of the LRF filter). For slit E there is a clear discrepancy, but this does not affect the inferred BH mass (see Section 4.4). We model the (face-on) intrinsic flux profile as a double exponential,

$$F(R)=F_1\mathrm{exp}(R/R_1)+F_2\mathrm{exp}(R/R_2),$$
(1)
and assume that the disk is infinitesimally thin and viewed at an inclination $`i=66^{}`$, (cf. Section 2.2). The total flux contributed by each of the two exponential components is $`I_i=2\pi R_i^2F_i`$ ($`i=1,2`$), and the overall total flux is $`I_{\mathrm{tot}}=I_1+I_2`$.

The best-fitting parameters of the model flux distribution were determined by a simultaneous fit to the WFPC2 and STIS H$`\alpha `$+\[N II\] fluxes. We use a $`0.8^{\prime \prime }\times 0.8^{\prime \prime }`$ square region of the H$`\alpha `$+\[N II\] image centered on the nucleus. The flux errors take into account both Poissonian errors and a $`15\%`$ error due to uncertainties in the alignment and subtraction of the on- and off-band image (see Section 2.1). From the STIS observations we use the H$`\alpha `$+\[N II\] flux of the Gaussian emission-line fits (Section 3.2). The fluxes from Gaussian fits agree within the errors with fluxes estimated by directly measuring the spectral lines. The STIS flux errors take into account readnoise and Poissonian errors. Typically, the STIS fluxes have much smaller relative errors than the WFPC2 fluxes. We perform an iterative fit of the double exponential to the H$`\alpha `$+\[N II\] flux data, taking into account flux errors and the necessary convolutions with the appropriate PSF, pixel size and aperture size for each setup. The diamonds in Figure 5 show the predictions of the model that best fits all available data simultaneously. This model has parameters $`R_1=0.049^{\prime \prime }`$, $`R_2=0.78^{\prime \prime }`$ and $`I_1/I_{\mathrm{tot}}=0.31`$. The absolute calibration gives $`I_{\mathrm{tot}}=7.7\times 10^{13}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ for H$`\alpha `$+\[N II\]. The lowest panel in Figure 5 shows the intrinsic (i.e., deconvolved) flux distribution as a function of radius.

The thin disk flux model provides a decent fit to the observed STIS fluxes, apart from a slight under-prediction of the central fluxes in slit W. The WFPC2 measurements are also reasonably well-fitted, although the observed flux profile on the major axis suggests a narrower central flux peak. A second fit was made as well in which we included only the WFPC2 data. This resulted in a considerably narrower intrinsic profile (Figure 5). However, the two model emission-flux distributions result in BH masses which agree to within $`25\%`$ as discussed in Section 4.4.

### 4.2 The Stellar Luminosity Density

For the purpose of dynamical modeling we parameterize the three-dimensional stellar luminosity density $`j`$. We assume that $`j`$ is oblate axisymmetric, that the isoluminosity spheroids have constant flattening $`q`$ as a function of radius, and that $`j`$ can be parameterized as

$$j(R,z)=j_0(m/a)^\alpha [1+(m/a)^2]^\beta ,m^2R^2+z^2q^2.$$
(2)
Here $`(R,z)`$ are the usual cylindrical coordinates, and $`\alpha `$, $`\beta `$, $`a`$ and $`j_0`$ are free parameters. When viewed at inclination angle $`i`$, the projected intensity contours are aligned concentric ellipses with axial ratio $`q^{}`$, with $`q^2\mathrm{cos}^2i+q^2\mathrm{sin}^2i`$. The projected intensity for the luminosity density $`j`$ is evaluated numerically.

We assume $`i=66^{}`$ as discussed in Section 2.2. We take $`q^{}=0.81`$ within a radius of $`3^{\prime \prime }`$ based on the isophotal analysis of VK99 (their Figure 2). They determined that the isophotes are very close to elliptical and that $`q^{}`$ varies from $`0.81`$ to $`0.76`$ between major-axis radius $`R3^{\prime \prime }`$ and $`R7^{\prime \prime }`$. The isophotal major-axis PA=$`156^{}`$ is nearly constant, varying by less than $`<5^{}`$ over the same radius range. A model with constant $`q`$ and PA is expected to be adequate in view of the main uncertainties for the BH mass from our dynamical modeling as discussed in Section 6. For example, we verified that the variation in ellipticity changes the derived circular velocities by $`<6\%`$, which is small compared to our velocity errors.

We determine the mean flux profile over three rows (corresponding to a width $`0.14^{\prime \prime }`$) along the unobscured eastern galaxy semi-minor axis (cf. Section 2.2). We use the $`I`$-band fluxes as opposed to $`V`$-band fluxes because the flux contribution of emission-lines is negligible in this wavelength region. The observed surface brightness profile is shown in Figure 6. The projected intensity profile of the best-fit model using the fluxes on the semi-minor axis between $`r=0.0^{\prime \prime }`$ and $`r=10^{\prime \prime }`$ is shown by the solid curve. This model has $`\alpha =0.36`$, $`\beta =1.02`$, $`a=0.35^{\prime \prime }`$ and $`j_0=1.6\times 10^2\mathrm{L}_{}\mathrm{pc}^3`$.

The surface mass density of the dust disk, which is proportional to the optical extinction (cf. VK99; and references therein), does not show a strong increase towards the center but varies by just about a factor of two over the disk (cf. Figure 2 in VK00). Taking a Galactic value of the gas-to-dust ratio of 100 (e.g., Bohlin, Savage & Drake 1978), the total mass of the disk, which extends over a region $`r2.5^{\prime \prime }`$ is $`5.8\times 10^7\mathrm{M}_{}`$. The ratio of enclosed disk mass to stellar mass is at most $`6\times 10^4`$ in the range $`m=[10\mathrm{pc},1000\mathrm{pc}]`$ (i.e., $`m=[0.03^{\prime \prime },3^{\prime \prime }]`$. Hence we can neglect the gravitational potential exerted by the gas and dust disk.

There is one complication in estimating the stellar mass profile. The central $`0.23^{\prime \prime }`$ of the galaxy has a color $`VI1.35`$ on both the eastern and western side of the galaxy nucleus. This color is typical of dust-free giant ellipticals. It implies that the dust either does not extend all the way to the center (i.e., has a ring structure) or becomes optically thick towards the center. In the latter case it would hide approximately half of the stellar mass in this region as the dust is very close to the center. We discuss this scenario in Section 4.4 and find that it requires an unusually low stellar mass-to-light ratio to explain the observed gas kinematics.

### 4.3 Dynamical Models for a Thin Rotating Gas Disk

Our thin-disk models for the gas kinematics are similar to those employed in van der Marel & van den Bosch (1998) and VK00. The galaxy model is axisymmetric, with the stellar luminosity density $`j(R,z)`$ chosen as in Section 4.2 to fit the available surface photometry. The stellar mass density $`\rho (R,z)`$ follows from the luminosity density upon the assumption of a constant mass-to-light ratio $`\mathrm{{\rm Y}}`$. $`\mathrm{{\rm Y}}`$ is included as a free parameter in the dynamical model. We assume that the gas is in circular motion in an infinitesimally thin disk in the equatorial plane of the galaxy, and has the circularly symmetric flux distribution $`F(R)`$ given in Section 4.1. Thus the position angle of the gas disk is set by the position angle of the galaxy which is known to an accuracy of $`1^{}`$. The inclination $`i`$ of the gas disk is assumed to be identical to the dust disk inclination. Below, we will verify PA and $`i`$ also from the two-dimensional gas kinematics independently. The circular velocity $`V_\mathrm{c}(R)`$ is calculated from the combined gravitational potential of the stars and a central BH of mass $`M_{}`$. The line-of-sight velocity profile of the gas at position $`(x,y)`$ on the sky is a Gaussian with mean $`V_\mathrm{c}(R)\mathrm{sin}i\frac{x}{R}`$ and dispersion $`\sigma _{\mathrm{gas}}(R)`$, where $`R=\sqrt{x^2+(y/\mathrm{cos}i)^2}`$ is the radius in the disk. The velocity dispersion of the gas is assumed to be isotropic, with contributions from thermal and non-thermal motions: $`\sigma _{\mathrm{gas}}^2=\sigma _{\mathrm{th}}^2+\sigma _{\mathrm{turb}}^2`$. For practical reasons, we refer to the non-thermal contribution as ‘turbulent’, although its origin is not certain (see Section 4.5). We parameterize $`\sigma _{\mathrm{turb}}`$ through:

$$\sigma _{\mathrm{turb}}(R)=\sigma _{\mathrm{t0}}+[\sigma _{\mathrm{t1}}\mathrm{exp}(R/R_\mathrm{t})].$$
(3)
This functional form is meant merely to provide a fit to the observed dispersion and is not based on an underlying physical mechanism for the turbulence. The temperature of the gas is not an important parameter: the thermal dispersion for $`T10^4\mathrm{K}`$ is $`\sigma _{\mathrm{th}}10\mathrm{km}\mathrm{s}^1`$, and is negligible with respect to $`\sigma _{\mathrm{turb}}`$ for all plausible models.

The predicted VP for any given observation is obtained through flux-weighted convolution of the intrinsic VPs with the PSF of the observation and the size of the aperture. The STIS PSF is represented as a sum of 5 Gaussians (see Table 5). We ignore the effect of ‘pixel-bleeding’ for the STIS PSF, since this changes the kinematics predicted by the model only by a few $`\mathrm{km}\mathrm{s}^1`$ (Barth et al. 2001). The convolutions are described by the semi-analytical kernels given in Appendix A of van der Marel (1997), and were performed numerically using Gauss-Legendre integration. A Gaussian is fit to each predicted VP for comparison to the observed $`V`$ and $`\sigma `$. We define a $`\chi ^2`$ quantity that measures the quality of the fit to the kinematical data, and the best-fitting model was found by minimizing $`\chi ^2`$ using a ‘downhill simplex’ minimization routine (Press et al. 1992).

As mentioned, the observed VP and the VP predicted by the dynamical modeling are fitted with a single Gaussian. Thus we have to verify that the deviations from a Gaussian are the same in observations and model. Figure 2 shows that in some spectra, especially near the galaxy nucleus, the emission lines tend to have smaller peaks and broader wings than the Gaussian fit. The same trend is indeed seen in the modeled spectra. We verified that, apart from this difference, there are no systematic differences as a function of radius between the VP and the Gaussian fit. Furthermore, the modeled spectra allow us to determine the difference between the moments of the full VP and the Gaussian fit qualitatively, also for blended H$`\alpha `$+\[N II\] emission lines, as is the case in many spectra. The representation by single Gaussian fits turns out to be also quantitatively adequate in the sense that the differences between true and gaussian $`v`$ and $`\sigma `$ of the VPs are typically within the quoted errors for the observed Gaussian moments.

### 4.4 A Thin Disk Fit to the Gas Mean Velocities

We first analyze the H$`\alpha `$+\[N II\] gas disk mean velocities and defer discussion of the velocity dispersions to Section 4.5. We start out with models with no BH in which the gas disk has the orientation as determined from the photometry. Thus the gas disk is assumed to have the same inclination $`i=66^{}`$ and position angle PA=$`156^{}`$ as the dust disk and galaxy (see Section 2.2). This corresponds to an intrinsic axial ratio of $`q=0.77`$ for our spheroid model. The STIS long-slit observations were performed at a small angle $`\theta =8.44^{}`$ to the major axis (cf. Section 3.1).

First we vary the $`I`$-band mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ and find that a model with $`\mathrm{{\rm Y}}_I=3`$ in solar units fits the observed gas kinematics best (see Figure 7). We will refer to this model as the ‘standard model’. There are some differences between the model and observations. For slit E, the model predicts a steeper mean velocity gradient than observed. For slit W, the observed mean velocities around $`x0.5^{\prime \prime }`$ are systematically larger than predicted by the model. The total and reduced $`\chi ^2`$ of this model with one free parameter are $`89`$ and $`3`$, respectively. However, almost $`40\%`$ of the total $`\chi ^2`$ is contributed by just the two points at $`x=0.406^{\prime \prime }`$ in slit W and at $`x=0.356^{\prime \prime }`$ in slit E. These two large deviations might be caused by local deviations from the model, for instance local turbulence. In fact, at the latter location the velocity dispersion shows an off-nuclear peak. Discarding these two points, the reduced $`\chi ^2=2`$. Thus given the idealization of the model, a rotating disk model with localized turbulence is not too unlikely.

We can estimate $`\mathrm{{\rm Y}}_I`$ independently from existing empirical correlations between $`\mathrm{{\rm Y}}`$ and galaxy optical luminosity. For NGC 4335, $`m_B=13.5^m`$ (Paturel et al. 1997) and we assume $`BV=1.0`$, $`VI=1.3`$ and $`RI=0.45`$, which is typical for bright ellipticals and consistent with our WFPC2 observations (cf. VK99). We derive $`\mathrm{{\rm Y}}_I=2.5`$ using the relation between $`\mathrm{{\rm Y}}_R`$ and absolute $`B`$ magnitude determined by van der Marel (1991). The value $`\mathrm{{\rm Y}}_I=3.5`$ is found using the correlation between $`\mathrm{{\rm Y}}_V`$ and absolute $`V`$ magnitude as determined by Magorrian et al. (1998). The two values are consistent because both correlations show scatter of a factor $`1.5`$ in the observed mass-to-light values at given host luminosity. The value $`\mathrm{{\rm Y}}_I=3.0`$ derived from the gas kinematics falls nicely in the expected range $`2.5<\mathrm{{\rm Y}}<3.5`$.

The inclination of the galaxy and its embedded gas disk are assumed to be similar to the larger scale dust disk. Which $`i`$ is preferred by the gas disk kinematics alone? To address this question, we derive the stellar mass model and intrinsic emission-line flux model assuming $`i=40^{}`$ ($`q=0.4`$) and $`i=80^{}`$ ($`q=0.8`$). The profile shape of our stellar mass model (Section 4.2) does not depend on inclination but the central stellar mass density changes, i.e., it increases for decreasing inclination. The changes in the intrinsic emission-line flux profile are negligible. For a given $`\mathrm{{\rm Y}}_I`$, the increase in central stellar density increases the predicted velocities, while the decrease in inclination decreases the line-of-sight component of the predicted velocities. Thus these two effects compete. As Figure 8 shows, compared to the fit for $`i=66^{}`$ the fits are poorer for $`i=40^{}`$ and $`i=80^{}`$ and the best fit $`\mathrm{{\rm Y}}_I45`$ is higher than expected from the empirical correlations with host luminosity.

The PA of the gas disk is fixed in the model to be identical to that of the galaxy major axis, which is determined to within $`1^{}`$ by the WFPC2 isophotal analysis. We have varied the PAd of the gas disk major axis by $`\pm 20^{}`$, and find that the circular disk model fits best for PA$`{}_{\mathrm{d}}{}^{}=152.5^{}`$, i.e. $`3.5^{}`$ smaller than the PA of the galaxy major axis. The $`\chi ^2`$ increases by $`>40`$ for PAd which are $`\pm 10^{}`$ away from PA$`{}_{\mathrm{d}}{}^{}=152.5^{}`$. Thus the PA of the gas disk as preferred by the gas kinematics agrees well with the PA derived from the WFPC2 stellar isophotes.

In summary, we see that the kinematics of the gas disk support the disk inclination and PA and the mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ implied by independent methods. The slope of the central velocity measurements for each slit are well fitted by this model with no BH. This a posteriori validates our inclusion of the central velocity gradient to determine $`\mathrm{{\rm Y}}_I`$ because this gradient is evidently not affected significantly by the central BH mass. It also implies that the data provide an upper limit to the BH mass.

We now determine this BH mass upper limit. The models indicate that only the gradient between the central three points closest in absolute radial distance from the nucleus (i.e., at $`|x|<0.15`$ in slit C) are significantly affected by a BH with mass $`M_{}10^8\mathrm{M}_{}`$. We will therefore use these mean velocities to compute the $`\chi ^2`$. Since the model parameters other than $`M_{}`$ are optimized to fit the whole velocity profile, the $`\chi ^2`$ from the central three mean velocities is expected to follow approximately a $`\chi ^2`$ probability distribution with $`N_{\mathrm{df}}=31=2`$ degrees of freedom (3 velocity measurements and 1 free parameter, $`M_{}`$) and hence an expectation value $`\chi ^2=2`$. In agreement with this the standard model has $`\chi ^2=1.8`$. Figure 9 shows the increase in $`\chi ^2`$ as a function of increasing BH mass. From the $`\mathrm{\Delta }\chi ^2`$ for the one-parameter model, a $`M_{}>10^8\mathrm{M}_{}`$ is formally ruled out at the 3 $`\sigma `$ level.

A more conservative upper-limit can be obtained by assuming the gas disk is maximally face-on. This will decrease the inferred intrinsic axis ratio of the galaxy (see Section 4.2). Since no giant ellipticals flatter than E6 are observed, we take $`q=0.4`$ as the minimum acceptable axis ratio for NGC 4335. This implies a galaxy inclination $`i=40^{}`$. Figure 9 shows that a $`M_{}>1.8\times 10^8\mathrm{M}_{}`$ is ruled out at more than the $`3\sigma `$ level for $`i=40^{}`$. Such a model cannot fit the complete velocity curves as well as the standard model for any $`\mathrm{{\rm Y}}_I`$ (Figure 8). Moreover, the best fit is provided by $`\mathrm{{\rm Y}}_I5`$ for $`i=40^{}`$ which is at least $`40\%`$ larger than expected, as discussed above.

Section 4.1 shows that the WFPC2 flux profile suggests a narrower flux profile than the STIS data. We performed the gas dynamical modeling also with the double exponential fit to the WFPC2 data only. This results in a $`25\%`$ higher upper-limit to the inferred BH mass. The other notable discrepancy between the WFPC2 and STIS fluxes in slit E does not influence the BH mass measurement directly.

In Section 4.2 we discussed the possibility that there might be a central dust distribution inside $`r0.23^{\prime \prime }`$ which is completely opaque. This will alter the assumed stellar mass density, because inside $`r=0.23^{\prime \prime }`$ approximately half of the galaxy light will be obscured. We determine the stellar mass density assuming the intrinsic stellar flux inside $`r=0.23^{\prime \prime }`$ is two times the observed flux and assume that the obscuration by dust is negligible outside $`r=1^{\prime \prime }`$. The fitted model has $`\alpha =0.0`$, $`\beta =1.4`$, $`a=0.31^{\prime \prime }`$ and $`j_0=6.2\times 10^2\mathrm{L}_{}\mathrm{pc}^3`$. In this case a mass-to-light ratio $`\mathrm{{\rm Y}}_I1.5`$ is needed to obtain a good fit to the galaxy-dominated kinematics. This is an unusually low $`\mathrm{{\rm Y}}_I`$ and lies well outside the observed scatter of observed absolute magnitude vs. $`\mathrm{{\rm Y}}_I`$ correlations (cf. Section 4.4). Moreover, the emission-line flux is strongly peaked towards the center which suggests we have a direct view of the nucleus without significant obscuration by dust. We conclude that NGC 4335 probably does not have an opaque central dust distribution. If however it is present, then our upper limit on the BH mass is an overestimate.

In conclusion, there is good agreement between the disk inclination $`i`$ and PA as determined from gas kinematics and dust morphology and between the mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ from gas kinematics and from independent methods. These agreements form support for the thin disk model, which indicates $`M_{}<10^8\mathrm{M}_{}`$. The quantities $`\mathrm{{\rm Y}}`$, $`i`$ and PA are mainly constrained by the gas kinematics at $`r>0.25^{\prime \prime }`$, which is outside the BH radius of influence $`r_{\mathrm{BH}}`$ (with $`r_{\mathrm{BH}}=GM_{}/\sigma ^2`$, and $`\sigma `$ the central stellar velocity dispersion). The upper limit to the BH mass is determined from the velocity gradient at $`r<0.15^{\prime \prime }`$. The inferred BH mass is thus based on the assumption that the thin disk model holds for $`r<0.15^{\prime \prime }`$. The next Section considers if this assumption is supported by the gas velocity dispersions.

### 4.5 The Gas Velocity Dispersions

The velocity dispersion is only roughly fitted by a model of the form (3) with $`\sigma _{\mathrm{t0}}=50.0\mathrm{km}\mathrm{s}^1`$, $`\sigma _{\mathrm{t1}}=200\mathrm{km}\mathrm{s}^1`$ and $`R_t=0.2^{\prime \prime }`$ (Figure 10). The observed velocity dispersion varies between $`\sigma _{\mathrm{gas}}50\mathrm{km}\mathrm{s}^1`$ and $`\sigma _{\mathrm{gas}}250\mathrm{km}\mathrm{s}^1`$ (Figure 10). The value of $`50\mathrm{km}\mathrm{s}^1`$ corresponds to the instrumental and thermal line broadening. In the absence of a BH, the induced line broadening due to differential rotation over the apertures increases the predicted dispersion to at most $`70\mathrm{km}\mathrm{s}^1`$. The presence of a BH with mass $`10^8\mathrm{M}_{}`$ would only affect the central aperture in each slit, increasing the predicted dispersions to $`100\mathrm{km}\mathrm{s}^1`$. Thus for many apertures, the gas shows a larger velocity dispersion than predicted by the dynamical model of a thin rotating disk. Also a double Gaussian fit to each emission-line, to fit a broad and a narrow component of the emission-line, yields a narrow-line component with a velocity dispersion in excess of that expected from the dynamical model (e.g., $`\sigma _{\mathrm{gas}}190\mathrm{km}\mathrm{s}^1`$ at the nucleus). Thus, the ‘excess’ velocity dispersion cannot be ascribed easily to the presence of a second kinematic component in addition to the thin rotating disk. An excess velocity dispersion increasing towards the nucleus, i.e., similar to seen here, is often observed in other gas disks (e.g., van der Marel & van den Bosch 1998; VK00; Barth et al. 2001), with M 87 being a notable exception (Macchetto et al. 1997). Significant excess dispersion is only absent at $`r0.5^{\prime \prime }`$. These points are well-fitted by our model which is strong evidence that (i) the model of a thin disk in rotation applies at these radii and (ii) the model parameters that depend sensitively on these points, i.e., $`\mathrm{{\rm Y}}_I`$, PA, and $`i`$, are trustworthy regardless of the effects of excess velocity dispersion observed closer to the nucleus. What causes the excess velocity dispersion? Is it gravitational or non- gravitational in origin?

A gravitational origin would imply that the gas rotates at speeds lower than the circular velocity. Only if the excess in velocity dispersion is modest, $`\sigma /v1`$, do approximate formulae for this effect of ‘asymmetric drift’ (e.g., Binney & Tremaine 1987) allow the circular velocity to be derived from the observed mean velocity (e.g., Barth et al. 2001). This regime does not apply to the observed kinematics in the central regions if interpreted as asymmetric drift. Unfortunately, the BH mass depends exclusively on these nuclear apertures. It is therefore useful to determine the BH mass in the most extreme case that all dispersion in the central spectrum is gravitational. One option is that the gas forms a spheroidal distribution of collisionless cloudlets. To first approximation we can model the gas assuming a spherical isotropic model and solve the Jeans equations in similar fashion as done in VK00. A black hole mass $`M_{}6\times 10^8\mathrm{M}_{}`$ is then required to produce the observed $`\sigma =240\mathrm{km}\mathrm{s}^1`$. The model is simplistic and deviations from axisymmetry and isotropy can change the result by factors of a few (de Bruijne, van der Marel, & de Zeeuw 1996). Moreover the model assumes a constant number density to flux density ratio for the emission-line flux. Different ionization mechanisms (e.g., central point source or shocks) quite likely produce quite different relations between number density and emission-line emissivity as a function of radius. A second option is to assume that asymmetric drift is present except at the very center. The central velocity dispersion is then due entirely to differential rotation over the aperture. Also in this case, a $`M_{}6\times 10^8\mathrm{M}_{}`$ is required to produce the central $`\sigma =240\mathrm{km}\mathrm{s}^1`$ for $`i=66^{}`$ and $`\mathrm{{\rm Y}}_I=3`$. Thus if the thin disk model breaks down within the central pixel we might have significantly under-predicted the BH mass.

The thin disk model could be a good approximation if the excess velocity dispersion is caused by localized turbulence, which leaves the gas on average in circular rotation. The highly irregular behavior of the velocity dispersion points indeed to a non-gravitational origin. The velocity dispersion of the gas peaks towards the nucleus in slit C and W, but dips towards the nucleus for slit E. Interestingly, the location of irregular peaks in the velocity dispersion profile seems to coincide with deviations of the observed mean velocities compared to the model, especially in slit E. The relative flux of \[N II\]6583 and H$`\alpha `$ can shed light on the excitation mechanism. Figure 11 shows that the flux ratio \[N II\]6583 / H$`\alpha `$ roughly peaks towards the nucleus in all slits. This is qualitatively in agreement with the excess dispersion being due to turbulent shocks as they tend to enhance the emission-line ratio (Dopita et al. 1997). Interestingly, the maximum line ratio (\[N II\]6583 / H$`\alpha `$ $`5`$) in fact occurs at the location of the dip in the velocity dispersion in slit E (and a possible dip in slit W). No special feature is detected in the WFPC2 image (Figure 1) at these locations. However, the spine of the straight radio-jet observed at arcsecond-scale resolution with the VLA (PA=$`79^{}`$; Wrobel, Machalski & Condon 2002; Xu et al. 2000) passes right through the most conspicuous dip in slit E (between $`x=0.04^{\prime \prime }`$ and $`x=0.12^{\prime \prime }`$) if it continues straight down to the sub-arcsecond scales of interest here. Unfortunately, no radio jet emission is detected at VLBA scales (Xu et al. 2000). Surprisingly, the velocity gradient in slit E is regular and only slightly shallower than predicted by our standard model. Thus it seems that the irregular excess in dispersion does not significantly alter the mean rotational velocity.

The irregularity of the velocity dispersion profile and its correspondence to changes in the \[N II\]6583 / H$`\alpha `$ flux ratio favor a non-gravitational origin for the excess velocity dispersion. The irregularities do not correspond to irregularities in the gas mean velocities. Instead, the mean velocities suggest normal rotation. These observations are consistent with a gas disk for which the bulk of the gas orbits at the circular velocity but locally posesses additional kinetic energy in the form of turbulent motions. Testing the viability of such a hydrodynamical model is very interesting but beyond the scope of this paper. This leaves room for an alternative scenario, in which the excess velocity dispersion is predominantly gravitational (despite its irregular behavior) implying $`M_{}6\times 10^8\mathrm{M}_{}`$. The WFPC2 imaging and kinematics cannot rule out such a spherical gas distribution at the very center (Section 2.1).

## 5 Stellar Dynamical Models

The WHT/ISIS stellar kinematics can place constraints on the mass-to-light ratio $`\mathrm{{\rm Y}}_I`$ and the BH mass independently from the gas kinematics. We in turn discuss the constraint on $`\mathrm{{\rm Y}}_I`$ and BH mass. A fully general stellar dynamical modeling of an axisymmetric galaxy requires ‘three-integral’ models in which the phase-space distribution function $`f(E,L_z,I)`$ depends on the orbital energy $`E`$, the angular momentum component parallel to the symmetry axis, $`L_z`$, and a third, non-classical integral (e.g., Binney & Tremaine 1987). However, the inferred mass-to-light ratio from ground-based observations is generally insensitive to assumptions about the structure and dynamics of the galaxy. Models with different inclination (van der Marel 1991) or anisotropy (van der Marel 1994; 1999) yield the same value of $`\mathrm{{\rm Y}}`$ to within $`20`$%.

Thus, to derive $`\mathrm{{\rm Y}}_I`$, we constructed axisymmetric two-integral models for NGC 4335 in which the phase-space distribution function $`f(E,L_z)`$ depends only on the two classical integrals of motion. We use the same galaxy mass model as for the gas kinematics (i.e., $`\rho =\mathrm{{\rm Y}}_Ij`$ with $`j(R,z)`$ as derived in Section 4.2) and we also include a BH mass. The modeling procedure is the same as used (and described in more detail) by VK00 and similar to applications of models to large samples by, e.g., van der Marel (1991) and Magorrian et al. (1998). Figure 12 shows the observed and predicted kinematics as a function of major axis distance. A value of $`\mathrm{{\rm Y}}_I3.5`$ in solar units provides a good fit to the stellar mean velocity and velocity dispersion profiles at $`r>2^{\prime \prime }`$. This is $`20\%`$ higher than inferred from the gas. Given the $`20\%`$ uncertainty for $`\mathrm{{\rm Y}}`$, the stellar kinematics out to a radius of $`5^{\prime \prime }`$ are consistent with the mass-to-light ratio inferred by the gas velocity profile in the inner $`1^{\prime \prime }`$.

The two-integral models can only reproduce the observed peak in the stellar velocity dispersion in the inner $`2^{\prime \prime }`$ by invoking a BH of mass $`M_{}3\times 10^9\mathrm{M}_{}`$. If this BH mass is correct, the HST gas kinematics are resolving the radius of influence $`r_{\mathrm{BH}}2.1^{\prime \prime }`$. HOwever, these same gas kinematics firmly rule out such a massive BH. Instead the gas models predict a BH mass which is at least between 5 and 30 times smaller. This is very surprising, because the two gas modeling variants discussed in Section 4.5 ascribe all nuclear gas kinetic energy to gravitation. The gas is assumed to be in either a disk or a spherical cloudlet distribution. These two models are expected to bracket the true nuclear gas distribution. However, there are reasons to mistrust the high BH mass from stellar dynamical modeling. If radial anisotropy in the stellar motions is present, the BH mass is overestimated. As an example, Magorrian et al. (1998) determined the BH masses in a sample of 36 galaxies with similar modeling as performed here (i.e., two-integral modeling with ground-based kinematics and HST-based photometry). They found that for the 15 most spherical galaxies, a spherical model with radial anisotropy remove the need for a BH in all cases (see also van der Marel 1999). It turns out that this degeneracy can be lifted if stellar kinematics at HST resolution is available (i.e., resolving the sphere of influence of the BH) and the complete parameter-space for axisymmetric models is considered with three-integral modeling. Early-type galaxies do indeed harbor BHs but the masses inferred from two-integral and three-integral axisymmetric modeling differ systematically. We compare nine E/S0 galaxies in Magorrian et al (1998) which also have three-integral modeling (as listed in Kormendy & Gebhardt 2001). The average ratio of two-integral BH mass over three-integral BH mass is 3.7 with a minimum of 0.41 and a maximum of 10. Thus it is not too unexpected that, if we were to have stellar kinematics at HST resolution as well, three-integral modeling would indicate a BH mass in better agreement with the value $`M_{}6.0\times 10^8\mathrm{M}_{}`$ as derived from the gravitational modeling of the gas velocity dispersion. By contrast, the factor $`>30`$ difference in BH mass between the resolved thin disk modeling and the stellar dynamical modeling seems a far stretch. Unfortunately, we cannot expect to obtain a stronger constraint on BH mass with more general axisymmetric three integral modeling with the current data because we (i) lack stellar kinematics at HST resolution and hence do not resolve the $`r_{\mathrm{BH}}0.4^{\prime \prime }`$ of a $`M_{}=6\times 10^8\mathrm{M}_{}`$ and (ii) have ground-based stellar kinematics only along the major axis and with reliable measurements of only first and second order moments.

## 6 Discussion and Conclusions

Before analyzing the implications of the inferred black hole mass we first discuss the two main differences between our gas disk dynamical modeling and the modeling used for NGC 3245 presented by Barth et al. (2001) which constitutes the current state-of-the-art. First, instead of the observed emission-line flux distribution, we have used a double exponential fit to the flux distribution. For NGC 3245 the use of an exponential fit changes the inferred $`M_{}`$ by $`30\%`$. Moreover, it accounted less well for the wiggles in the velocity profile. In our case the deviations of the exponential model from the observed photometric and spectroscopic emission-line fluxes are estimated to have a $`25\%`$ effect on BH mass (see Section 4.4). Second, we do not incorporate the velocity shifts due to asymmetric illumination of the slit by the gas disk in the dispersion direction. In other galaxies this could be important to derive the BH mass as they produce surface brightness ’caustics’ (Maciejewski & Binney 2001). These caustics are not observed in NGC 4335. STIS observations of stars by one of us show that the velocity shift amounts to at most $`30\mathrm{km}\mathrm{s}^1`$ for a $`0.2^{\prime \prime }`$ wide slit in the extreme case of a star at the edge of the slit. Only the (central) apertures in the adjacent slits have an asymmetric flux gradient in the dispersion direction. The kinematics from these apertures are not taken into account in deriving the upper limit to $`M_{}`$. Hence we conclude that the velocity shift has negligible effect on the derived $`M_{}`$.

Black hole mass $`M_{}`$ correlates in general rather tightly with observed velocity dispersion $`\sigma `$ in the inner region of galaxies (Gebhardt et al. 2000; Ferrarese & Merritt 2000). Recently, Tremaine et al. (2002) completed a detailed analysis of the correlation and we will use their results in what follows. They find a best fit correlation of the form $`\mathrm{log}M_{}=\alpha +\beta \mathrm{log}(\sigma /\sigma _0)`$ with $`\alpha =8.13\pm 0.06`$ and $`\beta =4.02\pm 0.32`$ for $`\sigma _0=200\mathrm{km}\mathrm{s}^1`$. The observed spread in the correlation indicates that the intrinsic dispersion in $`\mathrm{log}M_{}`$ is 0.3 dex (perhaps smaller if observational errors are underestimated). Measuring the central $`\sigma `$ in similar fashion as in Tremaine et al. (2002) gives $`\sigma =282\mathrm{km}\mathrm{s}^1`$ for NGC 4335 (cf. Section 3.4) which corresponds to a predicted $`M_{}=5.4\times 10^8\mathrm{M}_{}`$. The $`3\sigma `$ upper-limit $`M_{}<10^8\mathrm{M}_{}`$ from the thin disk modeling of the gas mean velocities (Section 4.4) falls well below the relation, even when including the reported intrinsic dispersion in $`M_{}`$. The residuals are $`0.73`$ dex and $`0.43`$ dex for the upper-limit of $`M_{}=1.0\times 10^8\mathrm{M}_{}`$ (best-fit model) and $`M_{}=2.0\times 10^8\mathrm{M}_{}`$ (for a maximally face-on gas disk), respectively (see Figure 13). The residuals are even larger if we use the best-fit relation as determined by Ferrarese (2002) who infers a larger $`\beta =4.58`$. The $`M_{}3\times 10^9\mathrm{M}_{}`$ derived from stellar dynamics corresponds to an equally large, but positive residual of $`0.7`$ dex. Finally, the gravitational modeling of the gas velocity dispersions yields $`M_{}6\times 10^8\mathrm{M}_{}`$ corresponding to a residual of $`0.05`$ dex, well within the reported intrinsic dispersion of the correlation.

Which BH mass are we to trust? The analysis strongly supports a thin rotating gas disk model at $`r0.5^{\prime \prime }`$. The model provides a reasonable fit to the gas mean velocities and dispersions, the dust disk morphology and gas disk kinematics indicate the same inclination and PA, and the $`\mathrm{{\rm Y}}`$ from stellar and gas kinematics is consistent. However, there are doubts for the validity of the model at $`r0.5^{\prime \prime }`$. The WFPC2 and STIS fluxes are consistent with a thin disk surface brightness profile but the signal-to-noise of the WFPC2 image is too low to rule out a more spherical distribution. An excess velocity dispersion with an irregular profile is observed at $`r0.5^{\prime \prime }`$. Only an ad hoc explanation of localized random motion exists for this excess. It is not clear if such quasi-stationary turbulence in a thin disk is physically viable (e.g., Wada , Meurer & Norman 2002; and references therein). If the excess velocity dispersion is (partly) due to gravitational motion around the BH, then the thin disk model underestimates the true BH mass. Ascribing all gas kinetic energy, including the excess velocity dispersion, as counterbalancing the gravitational potential in a simple manner of isotropically moving collisionless cloudlets yields a BH mass which agrees well with the $`M_{}\sigma `$ correlation (Tremaine et al. 2002). Ascribing the nuclear velocity dispersion to rotational motion from a spatially unresolved disk yields an equally good agreement (Section 4.5). As discussed in more detail in Section 5, it is not too unreasonable to assume that the stellar dynamical modeling overestimates the BH mass by a factor of $`5`$. This would be due to radial anisotropy (often observed in bright ellitpicals such as NGC 4335) which the two-integral modeling does not take into account. In conclusion, the gas spheroidal model infers a BH mass in accordance with the empirical $`M_{}\sigma `$ relation, but remains very simplistic. If doubts about the validity of the thin disk modeling for $`r0.5^{\prime \prime }`$ were proven true, its inferred BH mass is expected to be an underestimate, driving the expected true BH mass to higher values in better agreement with the $`M_{}\sigma `$ relation. Similarly, the expected corrections for the stellar dynamical model bring its predicted BH mass in better agreement with the relation.

What are the implications if in reality NGC 4335 indeed harbors a BH with $`M_{}6.0\times 10^8\mathrm{M}_{}`$? Our results then suggest that gas dynamical modeling assuming thin rotating disks cannot be used to derive accurate black hole masses, even when all of the following are true: (i) the gas mean velocities very clearly suggest rotation; (ii) the surrounding dust disk appears regular; (iii) the inclination and position angle of the inner gas disk are well-constrained due to the use of three adjacent slits and these angles are both indicated by the gas disk kinematics and independent methods; and (iv) the derived $`\mathrm{{\rm Y}}_I`$ from gas and stars agree, ruling out asymmetric drift at $`r0.5^{\prime \prime }`$. This then would cast doubt on other $`M_{}`$ values determined from gas kinematics using similar models. Figure 13 shows that these measurements have a large influence on the upper-end of the $`M_{}`$ detections and hence on the slope of the correlation. Moreover, to date all measurements of $`M_{}`$ in nearby radio galaxies are based on gas kinematics. However, the fact that the inferred BH masses using this method follow the best-fit relation more closely than the stellar dynamical measurements argues against this worry of the validity of the gas dynamical modeling in general. The $`\chi ^2`$ per degree of freedom as used in deriving the best-fit relation (cf. Equation 3 in Tremaine et al. 2002) is 0.27 for the gas dynamical measurements while for the stellar dynamical measurements it is 1.30 (1.17 if one excludes the Milky Way BH mass measurement which is obtained from a stellar dynamical model which is completely different from those used for external galaxies). Moreover, the average residual of the gas dynamical BH mass measurements is positive, counter to what is expected if the gas mean velocities systematically underestimate the circular velocity and hence the inferred BH mass.

Independent determinations of $`M_{}`$ in galaxies from gaseous and stellar kinematics observed at high spatial resolution are crucial to address the worries about the gas dynamical modeling. This has been performed for IC 1459 (Verdoes Kleijn et al. 2000; Cappellari et al. 2002). Verdoes Kleijn et al. derive from the gaseous kinematics at 6 FOS pointings a black hole mass ranging from $`M_{}=1.5\times 10^8\mathrm{M}_{}`$ (thin disk model) to $`6\times 10^8\mathrm{M}_{}`$ (isotropic spheroidal model). However, Cappellari et al. (2002) infer $`M_{}=(2.6\pm 1.1)\times 10^9\mathrm{M}_{}`$ based on three-integral modeling of combined ground-based and HST/STIS stellar kinematics. Unfortunately, the HST stellar kinematics cannot be measured accurately inside the sphere of influence for this BH mass. Cappellari et al. also present a more complete view of the gas kinematics from a recent STIS long-slit observation. This indicates that while the data are consistent with the earlier FOS measurements, and the inferred BH mass is similar (modeled with independently developed software), the gas mean velocities are rather perturbed in this particular case. Moreover, IC1459 has an irregular dust distribution. Thus, the gas and dust properties in IC 1459 are quite different from the photometrically and kinematically well-behaved central disk in NGC 4335. An independent $`M_{}`$ measurement based on the stellar kinematics for NGC 4335 at HST resolution, and for other galaxies with similar gas disk kinematics, would test the validity of the gas dynamical modeling in well-behaved gas disks as opposed to irregular gas disks.

It is similarly crucial to understand the origin of the excess velocity dispersion commonly seen in nuclear gas disks. Presently, it is not clear under which circumstances the thin gas disks become locally turbulent (for example, gravitational or MHD instabilities) and if they remain globally stable (e.g., Wada, Meurer & Norman 2002; and references therein). A better idea for the origin of the excess is needed to improve on the highly idealized collisionless spherical model discussed in this paper. In this respect, it is also useful to determine the ionization mechanism of the emission-line gas in nuclear disks. If it is shocks we should be wary about the modeling with unperturbed, infinitely thin rotating disks. Furthermore, the disks of gas and dust might perhaps interact with ambient hot X-ray gas which is often present in the centers of bright ellipticals (e.g., Gunn 1979). The ultimate goal for gas dynamical modeling is thus to explain self-consistently the density, ionization state and dynamical state of the central gas distributions. This requires high S/N two-dimensional photometry, kinematics and line ratios to determine the dynamics, ionization state, density and temperature as a function of disk radius. This will provide not only accurate BH masses but also a vast improvement in our understanding of the fueling of BHs.

Kiloparsec-scale radio jets are only seen in early-type galaxies but never in spirals (with one possible exception known to the authors; Ledlow et al. 2001). For instance, the UGC FR-I galaxies all have Hubble types E-S0 (cf. VK99). Prior to this study, only BH masses $`M_{}2\times 10^8\mathrm{M}_{}`$ have been reported in nearby radio galaxies (i.e., NGC 4261, NGC 4374, M87, NGC 5128, NGC 6251 and NGC 7052, see Tremaine et al. 2002 for references). Interestingly, none of the seven galaxies with a Hubble type later than S0 in the the sample compiled by Tremaine et al. (2002) have $`M_{}>1\times 10^8\mathrm{M}_{}`$. We are aware of only one galaxy with a Hubble type later than S0 and a $`M_{}>2\times 10^8\mathrm{M}_{}`$, which is the Sombrero galaxy (M104, Hubble type Sa) with $`M_{}=1.0\times 10^9\mathrm{M}_{}`$, but this measurement is based on models less general than three-integral stellar dynamical models. By contrast, all BH mass upper limits and detections in a sample of 16 mostly early-type disk galaxies (Hubble type S0-Sb; Sarzi et al. 2001) are below $`2\times 10^8\mathrm{M}_{}`$. These results suggest that the differences in BH mass might be the underlying factor for this host preference of radio-jets. However, if our determination of a $`3\sigma `$ upper limit of $`M_{}<10^8\mathrm{M}_{}`$ on the BH mass of NGC 4335 is correct, then NGC 4335 would illustrate that BHs with $`M_{}<10^8\mathrm{M}_{}`$ are also capable of producing FR-I radio jets. This would argue against BH mass being the (only) parameter underlying the host morphology preference of radio galaxies.

It is a pleasure to thank Marcella Carollo and Chris O’Dea for stimulating advice, and to thank Chris Benn and Andy Longmore for carrying out our WHT service observations and providing helpful suggestions during the analysis. We would like the referee for his/her comments which helped to improve the paper significantly. Support for proposal #8236 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. This paper made use of the LEDA database: http://leda.univ-lyon1.fr. G.V.K. is grateful to the Leids Kerkhoven-Bosscha Foundation for financial support for a visit to STScI in February 2001.


text:
# NLO corrections to the photon impact factorTalk given at the 9th International High-Energy Physics Conference in Quantum ChromoDynamics (QCD ’02), Montpellier, France, 2-9th July 2002.

## 1 Introduction

Understanding the total cross section for the scattering of two highly virtual photons, having virtualities $`Q_1^2`$ and $`Q_2^2`$ at large centre-of-mass energy $`s`$ ($`sQ_1^2,Q_2^2`$) should be in reach of perturbative QCD. The process $`\gamma ^{}\gamma ^{}\mathrm{hadrons}`$ is considered to be an excellent testing ground for the applicability of perturbative QCD in the Regge limit . If the energy is high enough to validate Regge asymptotics but not too high in order to suppress unitarity corrections we expect the $`\gamma ^{}\gamma ^{}`$ cross section to be described by the BFKL  equation.

To leading logarithmic accuracy (LLA) the predicted cross section, based on the BFKL equation, rises too quickly with increasing $`s`$. The situation is very different at NLO. The calculation of NLO corrections to the BFKL Kernel was initiated in  and finally completed in . The corrections were first seen to be very large. However, their size is under control when additional collinear logarithms are taken into account  or when the kinematical conditions are forced to avoid these extra logarithms (rapidity vetoes) . The NLO corrections tend to lower the power rise of cross sections to values that seem to be compatible with the data .  seem to confirm this although the $`\gamma ^{}\gamma ^{}`$-cross section is considered using only NLO corrections to the kernel, ignoring those to the photon impact factor. These studies, however, can at best be viewed as an estimate of higher order corrections since they do not take care of higher order corrections to the coupling between external particles, virtual photons in the $`\gamma ^{}\gamma ^{}`$ case, and the NLO BFKL ladder. In order to make reliable predictions, being consistent to NLO, the NLO corrections to the coupling of virtual photons to the exchanged BFKL ladder, described by the impact factor, has to be taken into account. These corrections are currently under study  and the status of this work is reviewed in this contribution.

Besides $`\gamma ^{}\gamma ^{}`$ scattering perturbative QCD at high energies may be studied in any situation where a large rapidity gap between targets is observed and where a hard scale is involved in observing that final state (cf. Fig. 1). Prominent examples are the observation of forward jets in $`\gamma ^{}p`$ collisions at HERA  or the production of Mueller-Navelet jets  in hadron-hadron collisions. The coupling of the BFKL ladder to the relevant jet production vertex at NLO has been finished recently .

## 2 The calculational program

We focus our discussion on the example of $`\gamma ^{}\gamma ^{}`$-scattering which may serve as the canonical example for a scattering process in perturbative QCD at very high energy and is sometimes referred to as the gold-plated process to test BFKL dynamics.

At high energy we anticipate Regge factorization and as a result the total cross section for $`\gamma ^{}\gamma ^{}`$-scattering is written as a convolution (cf. Fig. 2)

$$\sigma _{\gamma ^{}\gamma ^{}}(s)=\mathrm{\Phi }_\gamma ^{}𝒢_\omega \mathrm{\Phi }_\gamma ^{}$$
(1)
where $`𝒢_\omega (𝒓^2,𝒓_{}^{}{}_{}{}^{2},s_0)`$ is the Green’s function for the exchange of two reggeized gluons, projected into the colour singlet state, obtained as a solution of the (NLO) BFKL equation. $`\mathrm{\Phi }_\gamma ^{}`$ is the impact factor for virtual photons under discussion. At leading order, this impact factor (Fig. 3) is calculated from cut quark box diagrams: the virtual photon splits into a $`q\overline{q}`$-pair and the reggeized gluons from the $`t`$-channel couple to the $`q\overline{q}`$ pair in all possible ways.

At NLO we have contributions from virtual corrections to the leading order diagrams as well as contributions with an additional gluon in the intermediate state (cf. Fig. 4). We express the impact factor in terms of the particle-Reggeon vertices $`\mathrm{\Gamma }_{\gamma ^{}q\overline{q}}^{(0)}`$, $`\mathrm{\Gamma }_{\gamma ^{}q\overline{q}}^{(1)}`$ and $`\mathrm{\Gamma }_{\gamma ^{}q\overline{q}g}^{(0)}`$, denoting the coupling of a virtual photon to Reggeon via a $`q\overline{q}`$-pair at LO, NLO and with an additionally emitted gluon, respectively. These vertices result from perturbative QCD amplitudes with colour octet exchange, taken in the high energy limit and using Regge factorization. Expanding the impact factor in powers of the strong coupling, $`\mathrm{\Phi }_\gamma ^{}=g^2\mathrm{\Phi }_\gamma ^{}^{(0)}+g^4\mathrm{\Phi }_\gamma ^{}^{(1)}`$ we would naively write the NLO impact factor in terms of these vertices as

$`\mathrm{\Phi }_\gamma ^{}^{(1)}`$ $`={\displaystyle \frac{dM_{q\overline{q}}^2}{2\pi }𝑑\varphi _{q\overline{q}}\mathrm{\hspace{0.17em}2}\mathrm{Re}\mathrm{\Gamma }_{\gamma ^{}q\overline{q}}^{(1)}\mathrm{\Gamma }_{\gamma ^{}q\overline{q}}^{(0)}}`$
$`+{\displaystyle \frac{dM_{q\overline{q}g}^2}{2\pi }𝑑\varphi _{q\overline{q}g}\left|\mathrm{\Gamma }_{\gamma ^{}q\overline{q}g}^{(0)}\right|^2}.`$ (2)
$`dM_i^2`$ and $`d\varphi _i`$ denote the invariant mass and the phase space of the respective intermediate states $`i=q\overline{q},q\overline{q}g`$.

The virtual corrections $`\mathrm{\Gamma }_{\gamma ^{}q\overline{q}}^{(1)}`$ have been calculated in . We have expressed all loop integrals in analytic form as an expansion in $`ϵ=(4D)/2`$, quite in contrast to  where all integrals are kept as they are. For the real corrections, we considered the square of the particle-reggeon vertex $`|\mathrm{\Gamma }_{\gamma ^{}q\overline{q}g}^{(0)}|^2`$, which has been calculated in  for longitudinally polarised virtual photons. In  we recently completed the real corrections by adding the contributions from transversely polarised photons. The real corrections have been considered in  as well, but not in a very suitable form for the task of finally evaluating the impact factor.

However, calculating the amplitudes as they are does not quite complete the task. The individual contributions are still infrared divergent and have to be combined in order to get the expected cancellation that has been shown previously . At the same time, one has to consider the subtraction of leading logarithmic terms. These are present in the virtual corrections and proportional to the well-known LLA gluon trajectory function. In the real contribution to the impact factor they arise as the additional gluon is emitted with a large rapidity separation to the $`q\overline{q}`$-pair. Both of these LLA-terms are individually infrared divergent as well as the emitted gluon becomes soft.

In  we consider the subtraction of infrared divergences and leading log terms in combination. We have extracted the infrared divergent contributions from real and virtual corrections and defined suitable subtraction terms. The difference of our results and the respective subtraction terms is finite upon integration over the gluon phase space. Re-adding the subtracted terms with the integrations over the gluon phase space performed, explicitly allows us to exhibit the infrared divergences as poles in $`ϵ`$ and cancel them successfully against those from the virtual corrections. The final result for the NLO impact factor reads

$$\mathrm{\Phi }_\gamma ^{}^{(1)}=\mathrm{\Phi }_\gamma ^{}^{(1,\mathrm{v})}|^{\mathrm{f}in}+C_A\mathrm{\Phi }_\gamma ^{}^{(1,\mathrm{r})}|_{C_A}^{\mathrm{f}in}+C_F\mathrm{\Phi }_\gamma ^{}^{(1,\mathrm{r})}|_{C_F}^{\mathrm{f}in}$$

$$\frac{2\mathrm{\Phi }_\gamma ^{}^{(0)}}{(4\pi )^2}\left\{\beta _0\mathrm{ln}\frac{𝒓^2}{\mu ^2}+C_F\mathrm{ln}(𝒓^2)\right\}$$

$$+|\mathrm{\Gamma }_{\gamma ^{}q\overline{q}}^{(0)}|^2\{C_A[\mathrm{ln}^2\alpha (1\alpha )s_0\mathrm{ln}^2M^2]$$

$$+2C_F[83\mathrm{ln}\alpha (1\alpha )\mathrm{\Lambda }^2+\mathrm{ln}^2M^2$$

$$+\mathrm{ln}^2\frac{\alpha }{1\alpha }]\}\frac{dM^2}{2\pi }\frac{d\varphi _{q\overline{q}}}{(4\pi )^2}.$$
(3)
Therein, $`𝒓^2`$ is the transverse momentum of the reggeized gluon in the $`t`$ channel, $`\alpha `$ is the light cone momentum fraction of the outgoing quark and the invariant mass of the $`q\overline{q}`$-pair is denoted with $`M^2`$. In the first line we have the finite remainders from the calculation from real (r) and virtual (v) corrections. The second line arises upon renormalization and includes a term proportional to $`\beta _0`$, the only term depending on the renormalization scale $`\mu `$. The last three lines are the finite remainders of the subtraction terms. The dependence on $`\mathrm{\Lambda }`$, characterising the cone in which the emitted gluon is considered to be collinear, will cancel against a similar term, implicit in the first line.

The subtraction of the leading logarithmic terms induces a scale $`s_0`$ which can be translated into a rapidity cutoff beyond which the emitted gluon will belong to the leading logarithmic term. However, since the particular choice of this scale is arbitrary, the NLL impact factors depend on it. This dependence was irrelevant at LLA, since a change in the scale

$$\mathrm{ln}\frac{s}{s_0}=\mathrm{ln}\frac{s}{s_1}+\mathrm{ln}\frac{s_1}{s_0}=\mathrm{ln}\frac{s}{s_1}+\mathrm{NLLA}$$
(4)
is of higher order w.r.t. the LLA. These NLLA terms are now taken care of and may be phenomenologically important.

## 3 Outlook and Conclusions

Besides the above discussion of the NLO impact factor, our calculations have the potential to give further insight into the photon wave function picture. This picture, in conjunction with the saturation model has been applied successfully to the description of both deep-inelastic and diffractive scattering cross sections at HERA, e.g. . First steps in this direction have been done in , showing that an extension of the current picture to a higher $`q\overline{q}g`$ Fock-state of the virtual photon is in principle possible. Further steps in this direction include a consistent treatment of infrared divergences in configuration space and remain to be done.

In order to complete the calculation of the impact factor we have to calculate the phase space integrals over the remaining infrared finite terms, defined in . We will express the phase space integrals in terms of a set of standard integrals. For first phenomenological applications this might best be done numerically.

With these results we will be able to calculate the $`\gamma ^{}\gamma ^{}`$ cross section to NLL accuracy. In combination with the NLO jet vertex  an interesting variety of phenomenological applications of the NLO BFKL equation is now possible. Numerical results for the cross sections for $`\gamma ^{}\gamma ^{}`$-scattering, production of forward jets at HERA and Mueller-Navelet jets at the Tevatron and the LHC are important goals of future work and remain to be done.