The Black Hole Mass of NGC 4151 from Stellar Dynamical Modeling

High-spatial resolution Hubble Space Telescope (HST) imaging of NGC 4151 had been previously acquired with the Wide Field Camera 3 (WFC3) through the F547M (medium-$V$; GO-11661, PI: M. Bentz) and F160W ($H$; GO-13765, PI: B. Peterson) filters.

Details of the reduction and processing of the F547M images are described in Bentz & Manne-Nicholas (2018), but in brief, four images through the F547M filter were acquired in a single orbit on 2010 July 3 with a two-point dither pattern to cover the central gap between the detectors. At each point in the dither pattern, a short exposure and a long exposure were obtained, allowing for correction of saturation in the nucleus due to the AGN while also providing good sensitivity to the extended host galaxy. The total exposure time was 2310.0 s, and the drizzled image covers a field of view of $2\farcm 7\times 2\farcm 7$ at a pixel scale of 0$\farcs$04.

Six individual F160W images were acquired in pairs with equal exposure times and a two-point dither pattern over the course of three orbits on 2015 December 7, 2015 December 28, and 2016 January 5. The total exposure time was 3317.6 s, and the drizzled image covers a field of view of $2\farcm 3\times 2\farcm 1$ at a pixel scale of 0$\farcs$1283.

The drizzled F547M and F160W images were each fit separately with Galfit  (Peng et al., 2002, 2010) to construct two-dimensional surface brightness profiles (see Bentz et al. 2013; Bentz & Manne-Nicholas 2018 for more details). The surface brightness models allowed the PSF of the central AGN and the background sky to be removed from each image, leaving just the host galaxy for further analysis.

To compare the surface brightness profiles of the galaxy in the two filters, we rotated the F160W image to match the orientation of the F547M image, blurred and rebinned the F547M image to match the PSF profile and spatial scale of the F160W image, and then cropped the two images so they were both centered on the AGN and included a common field of view. Comparing the positions of the few resolved sources in common between the two images, we found that this process aligned the two images to an accuracy of $\sim 0.5$ pix at the final spatial scale. We then used the IRAF task ellipse to fit elliptical isophotes to both images and measure their one-dimensional surface brightness profiles in a consistent manner. Finally, we applied the most recent calibration of the Vegamag zeropoints, corrected for Galactic extinction based on the Schlafly & Finkbeiner (2011) recalibration of the Schlegel et al. (1998) dust map of the Milky Way, and applied small color corrections to account for the differences between a $V-$band filter and F547M, and an $H-$band filter and F160W. Figure 2 (top) displays the one-dimensional surface brightness profiles in the inner regions of NGC 4151, with $V$ shown in blue and $H$ shown in red. The $V-H$ color is displayed in the bottom panel. The pixels inside a radius of $\sim 0\farcs 5$ (equivalent to a radius of about 4 pixels, and denoted by the vertical line in the bottom panel) are affected by residuals from the subtraction of the AGN PSF, but at radii outside this region the $V-H$ color of the galaxy is quite flat across the inner $\sim 1\farcm 0$.

The surface brightness of NGC 4151 in both the WFC3 F547M and F160W images were examined for their appropriateness in constraining the stellar luminosity density within the dynamical models. The F547M image shows clear regions of dust extinction and excess stellar flux, while the F160W image shows a smoother galaxy profile but at the expense of a lower intrinsic spatial resolution than the F547M image.

We first tried a Multi-Gaussian Expansion (MGE) analysis (Emsellem et al., 1994; Cappellari, 2002) using the AGN- and sky-subtracted images of the galaxy discussed above. MGE surface-brightness decomposition begins by finding the center then binning the field in angle and radius. Photometry is measured for each point on the grid, and then nested two-dimensional Gaussian profiles with variable axis ratios and central intensities are fit to this two-dimensional intensity map. During this process, we required that the PAs of any elongation in the Gaussians aligned with the kinematic PA determined previously (this restriction is mandated by our choice of axisymmetric modeling strategy). Although the observed velocity field shows twisting and misalignment with the photometric PA, indicative of the presence of a bar, we defer the study of bar models to the future.

The parameters of the best-fit nested Gaussian profiles are described by the intensity in counts, the width in pixels ($\sigma$), and axis ratio ($q$). These parameters can then be converted into surface density in $L_{\odot}$ pc^-2 and $\sigma$ in arcsec. For both the F547M and the F160W images, the MGE solutions provided a surface brightness profile that included a lot of bumps and other structures, rather than a smooth increase in surface brightness from larger radii into the galaxy center.

We then examined the potential of using the surface brightness components that had been fit to the images with Galfit in the process of subtracting the AGN and sky. These surface brightness components included an exponential disk, Sérsic bulge, and an elongated Sérsic component to represent the weak bar, with the central positions of all the galaxy components tied together, and the position angles of the central galaxy components (not including the disk) set to the kinematic position angle of the inner galaxy. While smoother than the MGE profile, a significant amount of structure was still visible in the surface brightness profile for both the F547M and the F160W image.

Finally, we also examined the results of fitting a Nuker profile to the images using Galfit, to promote smoothness in the surface brightness profile while still retaining some flexibility to account for the complicated morphology. Figure 3 demonstrates that the Nuker fit to the F547M image is very close to the combination of several Sérsic profiles, but is smoother. However, both these profiles overpredict the amount of stellar light in the central $0\farcs 1$ – the region so dominated by the AGN emission that any residual stellar light is poorly constrained. A comparison with the MGE model from Onken et al. (2014), suitably scaled in normalization, shows that the latter has a central core with the radius of the innermost Gaussian component $\sim 0\farcs 11$, while the Nuker model has a cusp. As will be shown later in Section 5, a cuspy profile results in a best-fit black hole mass of zero, which is clearly unphysical for this bright AGN. We therefore constructed another model described by the Zhao (1996) profile, which closely follows the Nuker model except the innermost $0\farcs 1$ – the region dominated by the AGN emission. We note that both Nuker and Zhao profiles are equally flexible and can represent cuspy or cored systems for different choices of parameters; here we used the former parametrization for the cuspy profile and the latter for the cored one. These two extreme cases likely encompass the true range of possible density profiles for this galaxy. The 3d luminosity density for our adopted Zhao profile is

with $\gamma=0$, $\beta=2.47$, $\alpha=0.65$, $r_{0}=0\farcs 113$, and $j_{0}=4.08\times 10^{6}\,L_{\odot}\,\mathrm{kpc}^{-3}$. We note that although $\gamma=0$ formally corresponds to a cored profile, the zero logarithmic slope is achieved only asymptotically, and in practice the density remains weakly cuspy at the resolution limit, unlike the MGE parametrization, which is much more obviously cored within $0\farcs 1$.

Determination of how mass traces light is integral to understanding the gravitational potential of the galaxy for modeling. Based on OSUBSGS $H-$band imaging combined with SDSS $g-$ and $i-$band imaging, Onken et al. (2014) produced color-color maps of NGC 4151 and found very flat colors across the galaxy bulge with average values of $g-i=1.1\pm 0.1$ mag and $i-H=2.4\pm 0.2$ mag. This is in line with our findings based on HST F547M and F160W imaging as described in Section 2.2, where we also find a very flat color across the galaxy bulge, with $V-H=2.9\pm 0.1$ mag.

Onken et al. (2014) adopted the models of Zibetti et al. (2009) to predict $\Upsilon_{H}\simeq 0.4\pm 0.2$  M${}_{\odot}/$L_⊙. Comparison of many different prescriptions for estimating $\Upsilon$ by Roediger & Courteau (2015), however, finds that Zibetti et al. (2009) consistently predicts the lowest values at these galaxy colors. Furthermore, when compared to the known $\Upsilon$ for resolved stars in M31 (Telford et al., 2020), the prescriptions of Zibetti et al. (2009) seem to be biased low at these colors.

Based on the analysis by Roediger & Courteau (2015), we therefore adopt $\Upsilon_{H}\simeq 0.7\pm 0.1$ M${}_{\odot}/$L_⊙ and $\Upsilon_{V}\simeq 2.9\pm 0.3$ M${}_{\odot}/$L_⊙ as our best estimates for the $H-$ and $V-$band stellar mass-to-light ratios in the nucleus of NGC 4151.

Observations of HI 21 cm emission from the large-scale disk of NGC 4151 led Davies (1973) to report an inclination of $26^{\circ}\pm 8^{\circ}$. Simkin (1975) examined the ratio of the minor to major axis using density tracings of wide-field photographic plate images and found the outermost isophotes of the galaxy to suggest an inclination of $21^{\circ}\pm 5^{\circ}$. The nearly circular shape of the galaxy disk in the HI 21 cm maps presented by Bosma et al. (1977) also support a face-on orientation.

This is at odds with what has generally been determined from optical imaging alone. For example, surface brightness decomposition of optical ground-based imaging with a modest FOV generally provides a minor-to-major axis ratio of $\sim 0.7$ for the galaxy disk (e.g., Bentz et al. 2009), suggesting an inclination of 46^∘. However, as argued by Davies (1973), this elongated structure is likely to be more representative of the weak galaxy bar than the disk, although they also note that there are difficulties with this interpretation as well.

Thus the inclination of NGC 4151 is somewhat uncertain (Simkin 1975 go so far as to describe NGC 4151 as “pathological”!), but is likely to be in the range of $20-45^{\circ}$.

For a nearly face-on orientation, it is very difficult to estimate the intrinsic flattening $q$ from the photometry (even a very flat disk would still appear nearly round in projection). Therefore, our strategy is to explore the parameter space of the inclination angle $i$ and intrinsic flattening $q$ independently. Specifically, we first take the 1d surface brightness profile $I(R)$ and deproject it under the assumption of spherical symmetry, obtaining the luminosity density profile

We then assume that the intrinsic luminosity density is actually axisymmetric and ellipsoidally stratified, with the axis ratio $q\equiv z/R$ being a free parameter in the model, where $q=1$ is spherical. If such a density profile is projected face-on, the result would be identical to $I(R)$, while for a nonzero inclination angle $i$, the projected minor axis is $q^{\prime}=\sqrt{1-(1-q^{2})\,\sin^{2}i}$  times smaller than the projected major axis. To preserve the angularly-averaged surface brightness profile, we multiply the scale radius of the 3d profile by $s\equiv 1/\sqrt{q^{\prime}}$, so that the geometric mean of the projected major and minor axes is the same as in the spherical model. For instance, with an inclination $i=30^{\circ}$ and flattening $q=0.25$, the projected major and minor axes are $1.07$ and $0.935$ times the original scale radius of the one-dimensional surface brightness profile. We believe that the extra flexibility allowed by independent variation of flattening and inclination (with additional constraints coming from kinematics) is more important than an accurate reproduction of the density profile per se.

When performing SD modeling, the distance to the galaxy is a large source of uncertainty in determining $M_{\mathrm{BH}}$. The scale or radii of the galactic features being probed and the value of $M_{\mathrm{BH}}$ scale linearly with the assumed distance. Based on the information available at the time, Onken et al. (2014) adopted a recessional velocity distance measurement of 13.9 Mpc from Pedlar et al. (1992), cautioning that this distance was highly uncertain but that no better estimates were currently available.

Indeed, the only other information available was a Tully-Fisher distance to NGC 4151 that was reported by Tully et al. (2009) as 3.9 Mpc and is highly-discrepant from the group-averaged distance to NGC 4151 of $11.2\pm 1.1$ Mpc. Bentz et al. (2013) recalculated the group-averaged distance to NGC 4151, excluding NGC 4151 itself, and reported a distance of $16.6\pm 3.3$ Mpc, albeit based on only three galaxies.

Since then, additional studies have worked to pin down the distance more accurately. Hönig et al. (2014) measured the physical size of the inner radius of the dusty torus using broad-band RM and the angular size of the torus with near-infrared interferometry, reporting a dust-parallax distance of $19.0^{+2.4}_{-2.6}$ Mpc to NGC 4151. Tsvetkov et al. (2019) reported the discovery of a Type IIP supernova in NGC 4151. Adopting a standardizable candle approach, they find a distance of $16.6\pm 1.1$ Mpc, but they also report a distance of $20.0\pm 1.6$ Mpc if they instead adopt an expanding photosphere analysis.

Most recently, a Cepheid distance to NGC 4151 has been obtained with HST WFC3 optical and near-infrared imaging (GO-13765, PI: B. Peterson), resulting in a distance of 15.8$\pm$0.4  Mpc (Yuan et al., 2020), which we have adopted for our modeling. At this distance, 1″ corresponds to 77 pc.

To characterize the PSF of our NIFS datacube, we took a slice at a wavelength corresponding to an AGN broad emission line and subtracted off a slice at a wavelength corresponding to the stellar and AGN continuum, in effect creating narrowband images with an on-band image and an off-band image. The residual image thus contains an image of the unresolved AGN point source at a wavelength corresponding to a broad emission line, allowing for characterization of the PSF.

Using Galfit, we modeled the PSF image with multiple two-dimensional Gaussian components, each having common centers but different widths ($\sigma$) and flux contributions. The final three-component model with circular Gaussian components is shown in Table 2. This three-component model accurately matches the wings and the core of the PSF. If only two components were used, the fit was not as good but the values agreed more closely with what was found by Onken et al. (2014).

For the SAURON datacube, the PSF was characterized by Dumas et al. (2007) as a single Gaussian component with width $\sigma=2\farcs 0$.

The kinematic position angle was derived using the method of Krajnović et al. (2006). The kinematic measurements derived from an initial pPXF run were smoothed and bi-symmetrized before determining the best-fit systemic velocity along with the kinematic position angle of the system as measured from the y-axis of the data cube. This angle represents the line of maximum galaxy rotation, perpendicular to the axis of rotation, and for our kinematics is 37.2^∘ counterclockwise from the y-axis.

The photometric major axis of the galaxy is somewhat uncertain. At intermediate radii ($\gtrsim 30\arcsec$), the galaxy is visibly elongated in a direction almost perpendicular to the kinematic major axis (see Figure 1), while at smaller radii (within the footprint of the SAURON dataset, shown by a green rectangle), its elongation is still offset by some $30^{\circ}$ from the kinematic major axis. Furthermore, 21 cm imaging of the galaxy shows a nearly circular gas disk extending several arcmin beyond the high surface brightness stellar features, with a position angle of $26^{\circ}$ (Davies, 1973). The intrinsic galaxy shape is thus non-axisymmetric and radius-dependent. However, in the present study we do not attempt to model this complex morphology and instead assume an intrinsically axisymmetric shape. The use of axisymmetric models requires that the kinematic and photometric major axes are aligned, which motivates the constraint on the position angle during the photometric fit.

Next, we employed the adaptive binning method of Cappellari & Copin (2003). This binning procedure is based on Voronoi (1908) tessellations and assigns adjacent spaxels into bins in a scheme that approximates constant signal-to-noise (S/N) across the FOV. Bins that are similar in shape to the spaxels are prioritized, i.e. not overly irregularly shaped, even at the expense of some irregularity in the S/N. This Voronoi binning method preserves the high spatial resolution at the center of the galaxy where the S/N is highest while adaptively binning the spaxels at the edges of the field to increase their collective S/N.

The data and noise cubes were first cropped, ensuring that the center of the FOV was defined by the centroid of the AGN position, and that spaxels at the edges of the cube with very low counts due to dithering between exposures were rejected. This left us with 49$\times$49 spaxels, with the odd number of spaxels in each dimension ensuring that the central spaxel contained the AGN. We then produced both a slice of the data cube (the signal) and a slice of the noise cube that well-represented the cubes as a whole, focusing on the region between the 16200 Å  CO (6-3) bandhead and the strong [FeII] emission line at 16440 Å. These signal and noise slices were then used to assign the individual spaxels to bins.

We identified the Voronoi bin pattern for only one-quarter of the FOV, bounded by the kinematic major and minor axes, and then we replicated the bin pattern for the other three quadrants of the FOV. This scheme ensures that each bin has symmetric analogs in the other quadrants of the FOV, while the central spaxel containing the AGN was assigned to its own single bin. With this approach, we divided the FOV into 165 total bins.

Once the Voronoi bin assignments were determined, the spaxels that were assigned to each particular bin were combined for both the data cube and the noise cube. For all the spaxels assigned to a particular bin, the spectra from the data cube were co-added, channel by channel, and the spectra from the noise cube were added in quadrature. Voronoi binning the NGC 4151 integral field spectroscopy is an improvement of the analysis over that of Onken et al. (2014), where the cubes were rectified with a constant spatial sampling of 0$\farcs$2$\times$0$\farcs$2 across the FOV, thus degrading the spatial resolution in the crucial region near the black hole sphere of influence.

The NIFS spectra cover the wavelengths $\sim 15100-17700$ Å. Notable features within the wavelength range include forbidden Fe emission on the blue side, including the strongest line at 16440 Å, and hydrogen Brackett 11 and 10 lines on the red side at 16811 Å and 17367 Å, respectively. We focused on the wavelength range $\sim 15200-16400$ Å, which included the $\sim$16200 Å  CO (6-3) bandhead and blueward absorption. The spectra were fit within the wavelength ranges $\sim 15150-15330$, $15450-15500$, $15530-16000$, $16060-16130$, and $16180-16350$ Å, the same as Onken et al. (2014). These ranges satisfactorily masked the strong emission lines while focusing the analysis on the strongest expected stellar absorption features.

As LOSVDs can be described by Gauss-Hermite polynomials, when the noise of the spectra is high or the data are undersampled with a low velocity resolution, the proper use of pPXF is to bias or penalize higher order Gauss-Hermite terms toward zero, defaulting them to more Gaussian shapes. In the original study, Onken et al. (2014) adopted no penalty. Based on the details of the NIFS integral field spectroscopy, however, we adopted a BIAS of 0.3. We also examined the largest changes in kinematic measurements between BIAS = 0.3 and BIAS = 0 to determine which bins produce unreliable kinematic measurements that might need to be masked during the modeling process.

While convolving the velocity template spectra to match the observations, pPXF provides the option to include Legendre polynomials (Whittaker & Watson 1920) to improve the fit to the continuum. The additive polynomials (DEGREE) can mediate some of the effects of template mismatch (where poor fits of kinematics are obtained due to templates that incompletely represent the observations) and imperfect sky subtraction in the data reduction, while the multiplicative polynomials (MDEGREE) can aid with slight imperfections in the flux calibration and reduce or remove the need for reddening correction. We found a good balance between the weights of the various polynomial components and the goodness of fits to the galaxy spectra when adopting 2nd order additive and 2nd order multiplicative Legendre polynomials, identical to what was found by Onken et al. (2014). These low-order polynomials improved the overall shape of the best fits to the spectra without introducing localized fluctuations that might have inhibited our ability to constrain the stellar kinematics.

We adopted the point-symmetric two-sided pPXF fitting of our spectra, in which two identically-shaped bins symmetric across the center ($x,y$ and $-x,-y$) are fed to pPXF and fitted simultaneously. Point-symmetric fitting improves the S/N of the data and provides symmetric measurements that utilize the full information contained in each of the two spectra. This is another improvement in the analysis over that of Onken et al. (2014), where the bins were fit independently and then the output measurements ($V,\sigma,h_{3},h_{4}$) were four-fold symmetrized (also known as bisymmetrization, described below) instead and their errors from individual bins were added in quadrature and divided by two. We opted not to perform such bisymmetrization, in which spectra from four bins at positions ($\pm x,\pm y$) are fitted simultaneously, for two reasons. First, bisymmetrized kinematic maps hide any non-axisymmetric features (although we currently explored only axisymmetric models, in the future we plan to extend the analysis to a more general geometry). Second, in the point-symmetrized case we have essentially two independent variants of kinematic maps (one in the opposing pair of quadrants $x>0,y>0$ and $x<0,y<0$, the other – in the second pair $x>0,y<0$ and $x<0,y>0$; here $x$ and $y$ are coordinates aligned with the kinematic major axis). Even though for most models we fitted both halves of the dataset simultaneously (and hence the model maps, which are bisymmetric by construction, lie halfway between the two variants of the observed maps), we can fit the two variants separately to quantify the systematic differences in the model parameters arising from the asymmetries in the data – this is explored later in Section 5.2 and in Figure 8.

In our tests, we found that the G star was only used in 20% of the best-fit spectra, and even when it was used its contribution had a low weight, and so we omitted it as a template. The M star contributed 98% of the time and the K star contributed 97% of the time, and so we included both the K and M templates in our final fit. This is a slight difference from the method of Onken et al. (2014), who used only the M star template in the final kinematic fits.

With all of these adopted parameters, we determined the final stellar kinematics for the central region of NGC 4151 as shown in Figure 4. We were not able to reliably recover the stellar kinematic signature at the center because of the very bright AGN, so for the NIFS data we omit these bins by masking the central $3\times 3$ spaxels. This left 156 bins in our final kinematic maps, which are shown in Figure 4

In summary, the stellar kinematic maps presented here are based on the following improvements over the analysis procedures of Onken et al. (2014): (1) an improved telluric absorption correction; (2) the inclusion of noise information carried through the full reduction pipeline (rather than assuming a constant noise of 2%); (3) the adoption of Voronoi binning to preserve the native spatial resolution near the galaxy nucleus and maintain a constant S/N across the field; (4) the adoption of a penalty term in the implementation of pPXF; and (5) fitting the co-added symmetrized spectra in Voronoi bins with pPXF instead of symmetrizing the pPXF measurements in 0.2″bins.

We also reanalyzed the SAURON data of NGC 4151. The data cube was provided to us having already been Voronoi binned, with 482 total bins that were not symmetric across the FOV. We followed many of the same procedures as the original analysis carried out by Dumas et al. (2007), but we updated a few details, including the use of the MILES stellar template library of Sánchez-Blázquez et al. (2006) with the updates and corrections of Falcón-Barroso et al. (2011). The library contains 985 empirical stellar spectra covering $3525-7500$ Å with a well-constrained spectral resolution of 2.51 Å (Beifiori et al., 2011). We focused on a subset of 148 stars to consider as templates, selecting those stars that were in common with both the Indo-US and Elodie libraries (Valdes et al., 2004; Prugniel & Soubiran, 2001; Soubiran et al., 1998; Katz et al., 1998) and spanned spectral types F0 to K7.

During the fitting with pPXF, we included the entire spectral range of $4825-5380$ Å  except for small regions around the strong AGN emission lines, including H$\beta$ and [OIII] $\lambda\lambda 4959,5007$ Å. Following Cappellari et al. (2011), we determined a single best-fit template spectrum that was then adopted for all the bins, to avoid issues that might arise in the derived kinematics from imperfect velocity calibration of the various template stars in the library. We adopted a BIAS parameter of 0.4 and additive and multiplicative Legendre polynomials of degree 4 and 0, respectively. Because the data cube was already binned with a non-symmetric binning pattern, we were unable to carry out the kinematic fits using the two-sided point-symmetric option, so each bin was fit individually. Only the first two Gauss-Hermite moments, $V$ and $\sigma$, were included in the fit, with the other terms set to zero. We then masked several dozen AGN-contaminated spaxels within a few arcseconds from the center, retaining 447 bins in the final dataset, which is displayed in Figure 5.

The Schwarzschild (1979) orbit-superposition method is a standard approach for dynamical modelling of galaxies based on integrated-light stellar kinematics. In the last two decades, several independent implementations of this method have been widely used to determine black hole masses: the Nuker code (Gebhardt et al., 2000, 2003a; Thomas et al., 2004; Siopis et al., 2009), the Leiden code (van der Marel et al., 1998; Cretton et al., 1999; Cappellari et al., 2002), the MasMod code (Valluri et al., 2004, 2005), and the Heidelberg code (van den Bosch et al., 2008; van den Bosch & de Zeeuw, 2010).

In this paper, we use yet another recently developed code Forstand⁴⁴4The code is publicly available, for details see Vasiliev & Valluri (2019). (Vasiliev & Valluri, 2020). Although it can be applied to galaxies of any geometry, in the present case we limit ourselves to the axisymmetric case. We refer to that paper for a detailed description of the method, and here only briefly summarize the general workflow of orbit-superposition modelling and the particular aspects used in this study of NGC 4151.

1. 1.

   First, we determine the 3d luminosity density profile from the surface brightness profile. As explained in Section 3.3, we vary the assumed flattening $q$ and inclination $i$ independently, and each combination of these two parameters produces a different 3d density profile (although the spherically averaged profiles are all very similar). The mass density profile $\rho(\boldsymbol{r})$ is $\Upsilon$ times the luminosity density $j(\boldsymbol{r})$, where the mass-to-light ratio $\Upsilon$ is another free parameter that is varied at a later stage in the modelling pipeline. We denote the mass density profile with a fixed fiducial value of $\Upsilon_{0}$ as the “baseline” profile $\rho_{0}$.

2. 2.

   The total gravitational potential $\Phi(\boldsymbol{r})$ contains the contributions from the stars $\Phi_{\star}$ (which is related to $\rho_{0}$ by the Poisson equation), the black hole $\Phi_{\mathrm{BH}}\equiv-GM_{\mathrm{BH}}/r$, and optionally the dark matter halo $\Phi_{\mathrm{h}}$, for which we adopt a spherical NFW profile with a scale radius $r_{\mathrm{h}}$ and peak circular velocity $v_{\mathrm{h}}$. For simplicity, we fix $r_{\mathrm{h}}$ at $50\arcsec$ and vary only $v_{\mathrm{h}}$; since the inner part of the model is self-similar, it is sufficient to vary only the overall density normalization ($\rho\propto v_{\mathrm{h}}^{2}$).

3. 3.

   For each choice of model parameters $i$, $q$, $M_{\bullet}$, and $v_{\mathrm{h}}$, we construct an orbit library by integrating $N_{\mathrm{orb}}$ orbits for 100 dynamical times in the given potential and recording the LOSVDs of each orbit in each Voronoi bin (convolved with the PSF). The initial conditions for the orbits are sampled randomly: the positions are sampled from the stellar density profile, and the velocities – from a 3d Gaussian distribution with mean and dispersions given by the solution of an axisymmetric Jeans equation for the stellar density in the total potential. The random sampling ensures that all possible orbits supported by the given potential can enter the library, but the results (kinematic fit quality expressed in terms of $\chi^{2}$) do depend on the specific random realization, because some realizations might contain orbits that are more suitable for reproducing certain kinematic fluctuations caused by observational errors than other realizations. We stress that the scatter in $\chi^{2}$ values between different orbit libraries with the same model parameters is a general consequence of the highly flexible nature of the method, and the ability to consider different sets of initial conditions does not cause this effect but merely uncovers it. Nevertheless, this stochasticity is clearly undesirable when comparing $\chi^{2}$ values between models with different parameters, and we use two approaches to reduce its impact: (1) for each series of models with given $i$, $q$, and $v_{\mathrm{h}}$ but different $M_{\bullet}$, we use the same set of orbital initial conditions, and (2) we run series of models for several random realizations of the orbital initial conditions and consider the “ensemble average” $\chi^{2}$ contours to determine the range of best-fit model parameters.

4. 4.

   We then seek a solution for orbit weights that (1) reproduces the 3d density discretized on a $20\times 15$ cylindrical grid in the $R,z$ plane, and (2) minimizes the objective function $\mathcal{F}\equiv\mathcal{F}_{\mathrm{kin}}+\mathcal{F}_{\mathrm{reg}}$. The first term is the measure of fit quality for the kinematic constraints ($v$, $\sigma$, $h_{3}\dots h_{6}$). The second term is the regularization score, which penalizes large variations between orbit weights $\boldsymbol{w}$; in our case, $\mathcal{F}_{\mathrm{reg}}=\lambda\,N_{\mathrm{orb}}^{-1}\,\sum_{i=1}^{N_{\mathrm{orb}}}(w_{i}/\overline{w})^{2}$, where the mean orbit weight is $\overline{w}\equiv M_{\star}/N_{\mathrm{orb}}$, and the regularization coefficient $\lambda$ controls the tradeoff between smoothness (for large $\lambda$) and closeness of reproduction of kinematic constraints (for negligible $\lambda$). With a too small value of $\lambda$, there is a risk of overfitting (the model trying to reproduce noise in the data), which is undesirable because $\mathcal{F}_{\mathrm{kin}}$ of the best-fit solution exhibits large fluctuations between adjacent points in the parameter space. While a detailed analysis of model sensitivity and fidelity as a function of regularization coefficient is beyond the scope of this study, after some experiments we settle on a value $\lambda=100$ that produces reasonably smooth likelihood surfaces.

The solution for orbit weights in step 4, and the corresponding $\chi^{2}$ score, are obtained separately for each value of $\Upsilon$, but reusing the same orbit library constructed in step 3, only rescaling the velocity in the model by $\sqrt{\Upsilon/\Upsilon_{0}}$ (which corresponds to a rescaling of all masses by $\Upsilon/\Upsilon_{0}$). We start the search from a plausible initial guess for $\Upsilon$ and increase or decrease it by a factor $2^{1/16}\approx 1.04$, until reaching a difference between $\chi^{2}(\Upsilon)$ and $\chi^{2}_{\mathrm{min}}$ greater than 100. The whole process is then repeated from steps 1–3 for a different choice of other model parameters ($i$, $q$, $M_{\bullet}$, and $v_{\mathrm{h}}$). In total, we considered a few thousand models, which took $\sim 10^{3}$ CPU hours (the code is parallelized for multi-core architectures, so the wall-clock time is far shorter).

A technical detail worth mentioning is that when using Gauss–Hermite moments as observational constraints, one needs to decompose the LOSVDs of orbits in the model in the same basis, using the observed values $v$ and $\sigma$ as the parameters of the Gaussian–Hermite series for each bin, but expressing the measurement uncertainties $\delta v$, $\delta\sigma$ as uncertainties on the first two coefficients $\delta h_{1}$, $\delta h_{2}$, whose measured values are zero by construction. This makes the kinematic objective function $\mathcal{F}_{\mathrm{kin}}=\sum_{b=1}^{N_{\mathrm{bins}}}\sum_{c=1}^{6}\big{[}(h_{b,c}^{\mathrm{model}}-h_{b,c}^{\mathrm{data}})/\delta h_{b,c}\big{]}^{2}$ quadratic in orbit weights and amenable to efficient quadratic optimization solvers. After obtaining the best-fit solution for orbit weights, we then compute the final $\chi^{2}$ with respect to the original measured values $v,\sigma,h_{3}\dots h_{6}$ and their associated uncertainties, which is somewhat different from $\mathcal{F}_{\mathrm{kin}}$. Nevertheless, the shapes and locations of the minima are similar for both $\chi^{2}$ and $\mathcal{F}_{\mathrm{kin}}$ as functions of model parameters.

Given the relatively large number of free parameters ($i$, $q$, $v_{\mathrm{h}}$, $M_{\mathrm{BH}}$, and $\Upsilon$), we do not attempt to cover this parameter space exhaustively, but use a multi-stage strategy.

It is clear that $M_{\mathrm{BH}}$ is only constrained by the NIFS kinematics, since the SAURON data has lower spatial resolution, and moreover, we excise the central few arcsec of SAURON dataset because of AGN contamination. At the same time, the small spatial coverage of NIFS makes it insensitive to the dark matter halo normalization $v_{\mathrm{h}}$, which is merely a nuisance parameter in the present context. Since our primary goal is to determine $M_{\mathrm{BH}}$, we first consider a series of models fitted to NIFS kinematics alone and ignore the dark halo.

Figure 6 shows a four-dimensional grid of models in the parameter space of inclination $i$ (increasing from top to bottom rows), intrinsic axis ratio $q$ (increasing from left to right columns), and in each panel, black hole mass $M_{\mathrm{BH}}$ (abscissa) and mass-to-light ratio $\Upsilon$ (ordinate). The two leftmost panels in the lowest row show examples of models that use the cuspy Nuker density profile, and the best-fit $M_{\mathrm{BH}}$ remains near zero for all such models regardless of $i$ and $q$. The radius of influence $r_{\mathrm{infl}}\equiv G\,M_{\mathrm{BH}}/\sigma^{2}$ for a $2\times 10^{7}\,M_{\odot}$ black hole and $\sigma\simeq 100$ km s^-1 is $\sim 0\farcs 1$, at the limit of resolution of both kinematic and photometric data. Naturally, the stellar mass within this radius is comparable to $M_{\mathrm{BH}}$, but also varies by a factor of a few between the models with cuspy ($\sim 4\times 10^{7}\,M_{\odot}$) or cored ($\sim 1.5\times 10^{7}\,M_{\odot}$) profiles. We therefore conclude that a cuspy Nuker profile has too high stellar mass in the innermost region, which obviates the need for a black hole. As this is clearly in contradiction with observational evidence for a black hole as demonstrated by AGN activity, in the rest of the paper we consider only the models with the cored Zhao profile; we obtained very similar results with the cored MGE profile from Onken et al. (2014).

The lowest values of $\chi^{2}$ differ between panels, but the contours of $\Delta\chi^{2}\equiv\chi^{2}(M_{\mathrm{BH}},\Upsilon)-\chi^{2}_{\mathrm{min}}$ look similar in all panels, showing a “tilted valley” of acceptable models: higher $M_{\mathrm{BH}}$ values correspond to lower $\Upsilon$, such that the total gravitating mass within the region $\lesssim 0.5-0\farcs 6$ is approximately constant. In most panels, the marginalized 1d profiles of $\Delta\chi^{2}$ as a function of $M_{\mathrm{BH}}$, shown by green curves, have large, shallow minima in the range $0.5\times 10^{7}\,M_{\odot}\lesssim M_{\mathrm{BH}}\lesssim(3-4)\times 10^{7}\,M_{\odot}$.

Figure 7 shows a series of models constrained by both NIFS and SAURON kinematics, with each panel plotting the contours of $\Delta\chi^{2}$ as functions of $M_{\mathrm{BH}}$ and $\Upsilon$ for a given choice of other parameters. We fix the inclination to $i=30^{\circ}$ (the overall trends are similar for other values of $i$), and consider different choices of flattening $q$ (thickness increases from left to right) and dark matter halo normalization $v_{\mathrm{h}}$ (increases from top to bottom). We plot the contributions of NIFS and SAURON datasets to the total $\chi^{2}$ of each model separately by red and blue contours. It is clear that the NIFS contours are very similar to the ones in the third row of Figure 6, independently of $v_{\mathrm{h}}$, and they constrain a certain linear combination of $\Upsilon$ and $M_{\mathrm{BH}}$, as discussed earlier. By contrast, the SAURON contours are insensitive to $M_{\mathrm{BH}}$, but constrain a combination of $\Upsilon$ and $v_{\mathrm{h}}$. One could reasonably guess that a “universally acceptable” model must have compatible best-fit $\Upsilon$ from both datasets, which indeed minimizes the overall $\chi^{2}$. By adjusting the halo normalization $v_{\mathrm{h}}$, we can shift the location of best-fit SAURON $\Upsilon$ up or down, and it is natural to select such a value that maximizes the overlap with best-fit NIFS $\Upsilon$. Therefore, the dark matter halo normalization is not really an independent free parameter in our case, but rather should be determined from the consistency between the two datasets.

Comparing the series of models with different flattening $q$, we see that more disky models (smaller $q$) have higher best-fit $\Upsilon$ (especially for low inclination) and correspondingly lower $M_{\mathrm{BH}}$. This can be understood as follows: since we observe the galaxy at a nearly face-on orientation, the line-of-sight velocity dispersion $\sigma$ is determined by how far the stars travel in the vertical direction (i.e., thickness) and how large is the restoring force (i.e., the mass density). Indeed, for a thin isothermal disk with a surface mass density $\Sigma$ and scale height $h$, the vertical velocity dispersion is $\sigma=\sqrt{2\pi\,G\,\Sigma\,h}$ (Binney & Tremaine, 2008, problem 4.21). Therefore, when decreasing $q$ towards more disky models and hence decreasing $h$, we must simultaneously increase $\Upsilon$ and hence $\Sigma$ to keep the velocity dispersion at the observed value. On the other hand, when adding a dark matter halo, we increase the gravitating mass and hence reduce the stellar $\Upsilon$, but only in the outer parts (i.e., in the SAURON dataset), where stars are not overwhelmingly dominating the total potential. This suggests that the flattening $q$ is largely degenerate with $\Upsilon$ and $v_{\mathrm{h}}$, but nevertheless we find that the overall $\chi^{2}$ is lower for values of $q$ in the range $0.2-0.4$.

Finally, the assumed inclination does have a significant impact on the model properties. The mean rotational velocity of stars in the equatorial plane $\overline{v_{\phi}}(R)$ cannot exceed the circular velocity $v_{\mathrm{circ}}(R)\equiv\sqrt{R\,\partial\Phi/\partial R}$. In the case of a relatively cold disk, the difference (called asymmetric drift) is $v_{\mathrm{circ}}-\overline{v_{\phi}}\propto\sigma_{R}^{2}$ (Binney & Tremaine, 2008, equation 4.228). The projected rotational velocity is $\sim\overline{v_{\phi}}\,\sin i$, but for small inclination angles, the line-of-sight velocity dispersion is nearly independent of $i$, being determined by the mass density and thickness of the galaxy. Thus the observed line-of-sight velocity gradient sets the lower limit on the inclination angle $i_{\mathrm{min}}$, for which stars need to move on nearly circular (“cold”) orbits to produce the given rotational signal. For larger inclinations, stellar orbits have to be “warmer” (have higher eccentricity or even rotate in the opposite direction) in order not to exceed the projected rotational velocity.

Comparing different rows of Figure 6, we see that models with low inclination generally have higher $\chi^{2}$ even if the location of the minimum in the $M_{\mathrm{BH}}$–$\Upsilon$ plane is not strongly varying (except when both the inclination $i$ and the thickness $q$ are so low that the model is forced to have higher $\Upsilon$ to match the line-of-sight velocity dispersion). The tendency of orbit-superposition models to produce better fits at high inclination angles (edge-on orientations), even if the actual orientation is closer to face-on, has been noted in many studies (e.g., Section 5.3 in Thomas et al. 2007, or Lipka & Thomas 2021), and can be attributed to a greater flexibility in rearranging the orbits in the case of sub-maximal rotation. However, as discussed in Section 3.3, NGC 4151 appears fairly round at large radii, and certainly resembles a disk galaxy seen close to face-on rather than an intrinsically round galaxy seen edge-on. For this reason, we only consider models with $i\leq 30^{\circ}$, with an exception of two models with $i=40^{\circ}$ and $q\geq 0.5$. In any case, the range of acceptable values of $M_{\mathrm{BH}}$ does not strongly depend on the inclination, and we take $i=30^{\circ}$, $q=0.25$ for our fiducial series of models.

As mentioned in Section 4.3, point-symmetrized NIFS kinematic maps provide two independent variants of observational constraints and allow us to test the sensitivity of global model parameters to some asymmetries in the LOSVD. Indeed, the kinematic maps shown in Figure 4 are not symmetric with respect to reflection about the kinematic major axis (exchanging the upper and lower halves of the maps), although they are symmetric (for even moments) or antisymmetric (for odd moments) between the opposite quadrants (upper left and lower right quadrants are identical and represent one variant of data reduction, while upper right and lower left quadrants represent another, independent variant). In most cases we fitted both pairs of quadrants simultaneously, but in Figure 8 we show the results obtained from a series of models constrained by only one of the two variants of kinematic maps. The similar morphology of the $\chi^{2}$ contours supports the robustness of our results with respect to small variations in the kinematic maps.

Figures 9 and 10 show the NIFS and SAURON kinematic maps, respectively, for a fiducial model with $i=30^{\circ}$, $q=0.25$, $M_{\mathrm{BH}}=1.3\times 10^{7}\,M_{\odot}$, and $\Upsilon=2.5$, which has one of the lowest $\chi^{2}$ values. The maps actually look very similar for other nearby choices of parameters, even though the difference in $\chi^{2}$ may be rather significant (of order few tens). Variation of the black hole mass most noticeably affects $\sigma$ and $h_{4}$ maps in the innermost $\sim 0\farcs 2$.

Having considered the overall trends in the grid of models, we finally return to the main science question of this study – the measurement of the black hole mass. Summarizing the above discussion, we find that the addition of the SAURON kinematics does not tighten constraints on $M_{\mathrm{BH}}$, since it also brings another free parameter (dark matter normalization) that is largely degenerate with $\Upsilon$. The limited spatial extent of the SAURON data does not allow us to detect the spatial gradients of dynamical mass-to-light ratio associated with the gradual increase of the halo contribution with radius and hence to disentangle these two parameters. Since the halo properties are irrelevant for the present study, we focus primarily on the NIFS-only series of models (Figure 6), most of which have similar and rather broad ranges of $M_{\mathrm{BH}}$ masses with small $\Delta\chi^{2}$.

It is difficult to establish statistically strict constraints on $M_{\mathrm{BH}}$ for several reasons. First, the discrete nature of orbit-superposition models necessarily implies some noise in the value of $\chi^{2}_{\mathrm{min}}$. For each of the panels in that figure, we considered several random realizations of the orbit library, which exhibited variations of $\chi^{2}_{\mathrm{min}}\sim\mathcal{O}(10)$, and the locations of the minima were randomly scattered within regions roughly bounded by $\Delta\chi^{2}\lesssim 10$ (see, e.g., bottom row of Figure 2 in Vasiliev & Valluri 2020); in the plot, we show the contours averaged across these realizations. Second, the value of $\chi^{2}$ per constraint (reduced $\chi^{2}$) is $\sim 1.7$, indicating a moderately poor fit (or, more likely, underestimated measurement uncertainties). Third, there is little systematic study of statistical foundations of the Schwarzschild method, in particular, regarding the rigorous choice of confidence intervals on model parameters in terms of $\Delta\chi^{2}$. This choice cannot be guided by standard considerations applicable to the normally distributed errors (e.g., $\Delta\chi^{2}=1,4,\dots$ for $1\sigma$, $2\sigma$, etc.), since this does not take into account the highly non-parametric nature of the method: for each choice of “real” model parameters ($M_{\mathrm{BH}}$, $\Upsilon$, etc.), we consider the $\chi^{2}$ value produced by only one of many possible combinations of orbit weights (“hidden” parameters), instead of marginalizing over them (see Magorrian 2006 for a discussion).

Pending a rigorous statistical analysis, we adopt a simple qualitative prescription guided by our experiments with a large suite of models. Namely, we take the 1d $\Delta\chi^{2}$ profile, marginalized over $\Upsilon$ in each panel, and plot these profiles for several choices of $i$, $q$, and in the case of NIFS+SAURON models, $v_{\mathrm{halo}}$, to examine the overall range of acceptable $M_{\mathrm{BH}}$ values. Figure 11 shows $\sim 20$ such profiles, whose minima span the range $0.5\times 10^{7}\,M_{\odot}\lesssim M_{\mathrm{BH}}\lesssim 2\times 10^{7}\,M_{\odot}$; given the inherent noise in $\chi^{2}\gtrsim\mathcal{O}(10)$, we extend the range of acceptable values to $0.25\times 10^{7}\,M_{\odot}\lesssim M_{\mathrm{BH}}\lesssim 3\times 10^{7}\,M_{\odot}$. The corresponding mass-to-light ratios lie in the range $\Upsilon\simeq 2.5\pm 0.3$, consistent with the photometric estimates (Section 3.2).

The results presented in the previous section generally agree with those of Onken et al. (2014), which is unsurprising given the large error bars in both studies. Onken et al. (2014) found $M_{\mathrm{BH}}=(3.76\pm 1.15)\times 10^{7}\,M_{\odot}$ but assumed a distance of $D=13.9$ Mpc. When adjusted for the recently published Cepheid distance to NGC 4151 of $D=15.8\pm 0.4$ Mpc (Yuan et al., 2020), their best-fit mass becomes $M_{\mathrm{BH}}=(4.27\pm 1.31)\times 10^{7}\,M_{\odot}$. Our range of $M_{\mathrm{BH}}$ values is somewhat lower, which may be attributed to a somewhat lower velocity dispersion $\sigma$ in the central parts inferred from our reanalysis of the observational data. Their mass-to-light ratio $\Upsilon_{H}=0.34\pm 0.03$ is quoted for a different frequency band, and is generally consistent with our inferred $\Upsilon_{V}=2.5\pm 0.3$, as discussed in Section 3.2.

We also list in Table 3 all other $M_{\mathrm{BH}}$ determinations for NGC 4151 from other direct methods. As we did for the Onken et al. (2014) SD mass, we have adjusted the GD modeling mass to account for our adopted distance of 15.8 Mpc. We have also adjusted the RM masses so they have the same adopted value of $\langle{f}\rangle=4.3$ (Grier et al., 2013). The SD mass range we report here generally agrees with the two RM masses from Bentz et al. (2006) and De Rosa et al. (2018) and the H₂ GD modeling mass of Hicks & Malkan (2008).

Despite the improvements in the data reduction and analysis process, our uncertainties on $M_{\mathrm{BH}}$ are much larger than in previous studies of the same galaxy, primarily due to the fact that we have explored a wide range of parameters but have only reported results for marginalizing over $\Upsilon$. We acknowledge that the intrinsic shape and orientation of the galaxy are not well constrained, and instead of considering only one formally best-fit choice of these nuisance parameters, we accepted the $M_{\mathrm{BH}}$ values from a broad range of models. However, Figures 6, 7 demonstrate that even for a given choice of galactic geometry, there is a relatively large uncertainty on $M_{\mathrm{BH}}$ made possible by suitable rearrangement of orbits in our extensive orbit libraries. This confirms the expectations stemming from the analysis of mock data in Section 3.2 of Vasiliev & Valluri (2020), but stands at odds with more optimistic conclusions about the precision of $M_{\mathrm{BH}}$ measurements reached by some other studies, e.g., Neureiter et al. (2021). As these experiments did not use the same setup, it is clear that a more thorough investigation of this question is needed in the future.

The rather weak constraints on $M_{\mathrm{BH}}$ that we obtain in this study may seem overly pessimistic, but we recall that this galaxy is a rather difficult case from the observational perspective, since the AGN strongly dominates the light profile and limits the precision of kinematic measurements at small radii. For our estimated mass range $0.25\times 10^{7}\,M_{\odot}\lesssim M_{\mathrm{BH}}\lesssim 3\times 10^{7}\,M_{\odot}$, and our measured central value of $\sigma\sim 100$km s^-1 in the inner few pixels, the sphere-of-influence of the black hole is $1.1-11.2$pc or 0$\farcs 015-0\farcs$15. This implies that the radius of influence of the black hole is barely resolved by the innermost few pixels in the NIFS dataset. Moreover, the stellar luminosity profile in the same central region is also very difficult to constrain due to overwhelming AGN brightness, and as we have seen, for a plausible cuspy stellar profile consistent with observations, dynamical models prefer no black hole, which is clearly counterfactual. These considerations highlight the importance of spatially resolving the sphere of influence both photometrically and kinematically, especially in late type galaxies where black hole masses and stellar velocity dispersions are typically smaller than in elliptical galaxies.

Previous authors (e.g. Gültekin et al., 2009; Batcheldor, 2010; Gültekin et al., 2011) have focused on the estimated parameters of the $M_{\mathrm{BH}}$-$\sigma$ relation when $M_{\mathrm{BH}}$ measurements that do not resolve the black hole sphere-of-influence are included or excluded from the sample. Gültekin et al. (2011) found, for a sample of $M_{\mathrm{BH}}$ in early type galaxies with central velocity dispersions with median $\sigma\gtrsim 260$km s^-1, that excluding measurements which do not resolve the radius-of-influence of the black hole biases the estimated $M_{\mathrm{BH}}$-$\sigma$ relation. They concluded that instead of being excluded these measurements should be included as upper limits or with large error bars. For the elliptical galaxies and massive bulges in these studies, it is argued based on empirical evidence from repeated measurements, that resolving the radius-of-influence only helps to reduce the errors on the black hole mass estimate.

Using mock long-slit kinematic data Valluri et al. (2004) showed that if the radius-of-influence of a black hole is not resolved, nuclear kinematical data can be fitted without a black hole. As far as we are aware the work presented in this paper is the first modeling study in a late type galaxy that clearly demonstrates the drawbacks of not resolving the influence radius on the stellar dynamical $M_{\mathrm{BH}}$ measurement. Despite the fact that the AGN makes it difficult to resolve the influence radius of the black hole, its very presence necessitates a black hole. Yet we saw that because the central luminosity profile is not resolved below 0$\farcs$1, the assumption of a cuspy Nuker density profile (which is generally considered reasonable for a late type galaxy) would require no black hole to fit the kinematics. Therefore this study clearly illustrates the need for both photometric and kinematic data that resolve the influence radii of supermassive black holes, data that will become more readily available with JWST and the ELTs.

Efforts to model the RM observations presented by De Rosa et al. (2018) are currently underway and will remove the dependence on $\langle{f}\rangle$ and provide a fully-independent RM mass for comparison with the SD and GD-based masses. Such analyses rely on strong AGN variability during the monitoring campaign, as well as careful management of all noise sources and the observing cadence, and have only been possible for a handful of objects thus far. The observations presented by De Rosa et al. (2018) are of similar quality to previously successful RM modeling attempts (Pancoast et al., 2014; Grier et al., 2017a; Williams et al., 2018), and should therefore be sufficient for an accurate mass constraint.

Furthermore, NGC 4151 is the target of an Early Release Science program with JWST (ERS 1364, PI Bentz). Observations with the NIRSpec IFU will probe the nuclear stellar dynamics and will be directly compared with the observations and analysis presented in this work. NIRSpec is expected to provide some crucial advantages over AO-assisted ground-based observations with its stable and diffraction-limited PSF and significantly lower backgrounds. However, NIFS provides a higher spectral resolution, which may allow for more accurate measurements of the higher-order moments of the stellar absorption profiles. We will thus revisit the topic of the black hole mass in NGC 4151 from stellar dynamical modeling in the near future, once JWST has successfully launched.