The Broad Line Region and Black Hole Mass of NGC 4151

We present a reanalysis of reverberation-mapping data from 2005 for the Seyfert galaxy NGC 4151, supplemented with additional data from the literature to constrain the continuum variations over a significantly longer baseline than the original monitoring program. Modeling of the continuum light curve and the velocity-resolved variations across the H$\beta$ emission line constrains the geometry and kinematics of the broad line region (BLR). The BLR is well described by a very thick disk with similar opening angle ($\theta_o \approx 57^{\circ}$) and inclination angle ($\theta_i \approx 58^{\circ}$), suggesting that our sight line towards the innermost central engine skims just above the surface of the BLR. The inclination is consistent with constraints from geometric modeling of the narrow line region, and the similarity between the inclination and opening angles is intriguing given previous studies of NGC 4151 that suggest BLR gas has been observed temporarily eclipsing the X-ray source. The BLR kinematics are dominated by eccentric bound orbits, with $\sim10$% of the orbits preferring near-circular motions. With the BLR geometry and kinematics constrained, the models provide an independent and direct black hole mass measurement of $\log M_{\rm BH}/M_{\odot} = 7.22^{+0.11}_{-0.10}$ or $M_{\rm BH}=1.66^{+0.48}_{-0.34}\times10^7 M_{\odot}$, which is in good agreement with mass measurements from stellar dynamical modeling and gas dynamical modeling. NGC 4151 is one of the few nearby broad-lined Seyferts where the black hole mass may be measured via multiple independent techniques, and it provides an important test case for investigating potential systematics that could affect the black hole mass scales used in the local Universe and for high-redshift quasars.


INTRODUCTION
Reverberation mapping (Cackett et al. 2021) is one of only a few methods that are able to directly constrain the mass of a supermassive black hole through its gravitational effects on luminous tracers (stars or gas). Most methods that directly constrain black hole mass depend on spatial or angular resolution. The mass of Sgr A* in the Galactic Center has been determined from decades of monitoring the proper motions of individual stars (Ghez et al. 2000;Genzel et al. 2000;Ghez et al. 2008). Spatially resolved water maser clouds allowed accurate measurement of the mass of the central black hole bentz@astro.gsu.edu * Packard Fellow in the nucleus of NGC 4258 (Miyoshi et al. 1995). And dynamical modeling of spatially resolved stellar and gas kinematics have produced a collection of over 100 black hole mass measurements for galaxies in the nearby Universe (D 100 Mpc; see the review by Kormendy & Ho 2013). In a few special cases, VLT-GRAVITY has been able to push beyond these typical distance limitations. By combining 10-m class telescopes as an interferometer, GRAVITY has successfully probed photoionized gas at sub-pc angular resolutions in the nuclei of active galaxies IRAS 09149−6206 (z = 0.056; Gravity Collaboration et al. 2020) and 3C 273 (z = 0.158; Gravity Collaboration et al. 2018) and determined constraints on their black hole masses.
In contrast, reverberation mapping relies on time resolution, using light echoes to probe the physical arrange-ment and conditions of photoionized gas around an accreting supermassive black hole. With no angular resolution limit, reverberation mapping may be applied to active galaxies at any distance (e.g., Kaspi et al. 2007;Hoormann et al. 2019;Grier et al. 2019;Williams et al. 2021a,b). Results from reverberation mapping studies provide scaling relationships that allow quick estimation of large numbers of AGN black hole masses (e.g., Shen et al. 2011), permitting studies of black hole growth and evolution as a function of lookback time.
In effect, there are two black hole mass scales currently in use and assumed to be equivalent: one based on the results of stellar and gas dynamical modeling in mostly early type galaxies within D ≈ 100 Mpc, and one based on reverberation mapping results for active galaxies at larger distances. All black hole mass measurement techniques include inherent uncertainties and and potential systematic biases (cf. Graham et al. 2011;Peterson 2010;Kormendy & Ho 2013), and at the moment it is not clear that dynamical modeling and reverberation mapping give the same results for the same black holes. The Event Horizon Telescope results for Pōwehi, the nuclear black hole in M87, highlighted the need to compare direct black hole mass measurements: the mass derived from modeling of the interferometry data agreed with a previous measurement from stellar dynamics but disagreed with a gas dynamics measurement (Event Horizon Telescope Collaboration et al. 2019). But with bright AGNs being rare within the volume that allows for high angular resolution in galaxy nuclei, there are few opportunities to directly compare masses from reverberation mapping and stellar or gas dynamics.
NGC 4151 is one of the nearest broad-lined AGNs at z = 0.0033, and thus one of only a handful of AGNs where the black hole mass may be measured using multiple independent techniques. Its bright nuclear emission was first described by Campbell & Moore (1918) based on observations collected, in part, by Dr. Heber Curtis, and computations likely conducted by unpaid assistant 1 Miss Adelaide M. Hobe, who is credited with having carried out "the major part of the computations" in their study of bright line nebulae. Campbell & Moore noted that NGC 4151 had a spectrum resembling that of NGC 1068 as described by Fath (1909). Additional objects with similar properties were soon discovered, and Seyfert (1943) conducted the first detailed investigation of so-called extragalactic nebulae with bright 1 https://150w.berkeley.edu/celestial-observers-first-sixteen- berkeley-women-doctoral-graduates-astronomy-1913-1952 nuclear emission, including NGC 4151. We now know these objects as Seyfert galaxies. Variability on short timescales has come to be recognized as another typical characteristic of AGNs, including Seyfert galaxies. Bahcall et al. (1972) noted that flux variability from a central ionizing source would cause variations in surrounding photoionized gas in AGNs and some novae. This idea was then developed into a framework for mapping out the geometry and kinematics of the broad line region (BLR) in AGNs by Blandford & McKee (1982), a technique they dubbed "reverberation mapping".
As a prototypical Seyfert galaxy with strong historic variability (cf. the 110 year light curve presented by Oknyanskij et al. 2016) and observed rapid variability in its bright nuclear spectral lines (Cherepashchuk & Lyutyi 1973;Antonucci & Cohen 1983;Bochkarev 1984) NGC 4151 was a natural target for some of the first reverberation mapping studies (Peterson & Cota 1988;Clavel et al. 1990;Ulrich et al. 1991;Maoz et al. 1991). However, early spectroscopic monitoring programs were generally undersampled in the temporal domain because of an expectation from photoionization models that the BLR was an order of magnitude too large (Peterson et al. 1985). This size problem was eventually solved by replacing single zone photoionization models with models that included gas covering a range of temperatures and densities, such as the LOC Model (Baldwin et al. 1995).
Around the same time, it became clear that reverberation mapping could constrain the masses of the central black holes in these AGNs (Peterson & Wandel 1999. By combining the average time delay for a broad emission line with its Doppler-broadened width, the black hole mass could be determined modulo a scaling factor that included important details such as the inclination angle at which we view the system. As a stopgap, the use of a population-average scale factor f was introduced to bring black hole masses from reverberation mapping into broad agreement with the mass scale derived from stellar and gas dynamical modeling (Onken et al. 2004). However, from the beginning, reverberation mapping was understood to be able to provide all the information needed to recover the full geometry and kinematics of the BLR, thus precluding the use of f and allowing for an independent and direct measurement of M BH .
And so a series of intensive monitoring campaigns began in the early 2000s Denney et al. 2006) with the goals of improving the BLR measurements for objects that had previous measurements of poor quality or significantly undersampled data, and culminating in the acquisition of unambiguous velocity-resolved reverberation mapping data (e.g., Bentz et al. 2008Bentz et al. , 2009Denney et al. 2009;Grier et al. 2012).
Constraining the details of the BLR from velocityresolved reverberation data has been approached in two ways: either through forward modeling that explores the potential parameter space of BLR geometries and kinematics (e.g., Pancoast et al. 2011), or through the ill-posed inverse approach in which the time delay distribution as a function of velocity across the emission line (a velocity-delay map) is reconstructed directly from the data (e.g., Horne 1994;Skielboe et al. 2015;Anderson et al. 2021). The two approaches are complementary. Recovery of a velocity-delay map relies on a smaller set of core assumptions and is able to capture the full level of detail present in the data. However, the interpretation of a velocity-delay map is not straightforward and relies on comparison with models. Forward modeling, on the other hand, begins with a larger set of key assumptions to construct a fully self-consistent framework. The need to rely on some simplifying assumptions ensures that the models may not fully explore the level of detail available in the data. But the strength of forward modeling is that the results are relatively simple to interpret.
The superb data sets that have finally begun to be acquired by reverberation mapping programs have allowed the BLR structure and kinematics to be explored in detail for a modest number of AGNs using the forward modeling approach (Pancoast et al. 2014b;Grier et al. 2017;Williams et al. 2018Williams et al. , 2020Bentz et al. 2021;Villafaña et al. 2022a). While the exact details vary from AGN to AGN, these studies have generally found that the BLR, as probed by the Hβ emission line, is arranged in a thick disk geometry that we are viewing at low to moderate inclination. The kinematics are generally dominated by rotation that may also include a contribution from inflow or, in a few cases, some outflow. In the few instances where velocity-delay maps have been cleanly recovered (Bentz et al. 2010;Skielboe et al. 2015;Horne et al. 2021), the interpretations are in general agreement with the results from forward modeling.
Here, we reanalyze the spectroscopic monitoring data of NGC 4151 that were obtained in early 2005 ) as part of the push to achieve velocityresolved reverberation mapping data. Poor weather significantly shortened the duration of the program, and so the original goals of the program were scaled back to simply determining an accurate Hβ time delay from well-sampled, albeit short, light curves. In this work, we supplement the original observations with additional measurements from the literature, extending the temporal coverage of the continuum variations by an ex-tra ∼ 100 days. We model the continuum light curve and Hβ emission-line profiles with the phenomenological modeling code CARAMEL (Pancoast et al. , 2014a, and provide constraints on the BLR geometry and kinematics in NGC 4151 along with an independent and direct measurement of the black hole mass.

OBSERVATIONS AND DATA PREPARATION
The data analyzed in this work come from several sources. The observations that provide the Hβ spectra and the bulk of the continuum flux measurements were originally described by Bentz et al. (2006), and we provide a summary here. Long slit spectroscopy was collected on a ∼nightly basis between 2005 February 27 and 2005 April 10 at the MDM 1.3m McGraw-Hill Telescope. Two to four 1200 s spectra of NGC 4151 were collected each night for a total of 96 individual spectra. Typical CCD reductions were applied in IRAF and XVista.
For this work, we re-extracted the spectra from the calibrated 2D frames and trimmed the spectra while applying a common dispersion solution, so that each spectrum covered 4400 − 5700Å with a dispersion of 1.25Å pix −1 . In previous work with these data, the nightly spectra were averaged together at this point and a typical flux uncertainty of 2% was assumed. Here, we continued our analysis with each individual spectrum and we carried the error array along with the spectra through our analysis, capturing the higher signal-tonoise (S/N) in the emission lines compared to the continuum and the additional noise from night sky lines. A typical spectrum achieved S/N=100 per pixel in the continuum.
Variable seeing and slit placement over the two months of observations induced small variations in the wavelength solution, resolution, and flux calibration of the individual spectra. These variations were minimized with the van Groningen & Wanders (1992) scaling algorithm using the [O III] λ4959 emission line as an internal calibration source (the [O III] λ5007 line was saturated in several of our spectra and could not be used).
We then prepared the spectra for modeling by isolating the Hβ emission through spectral decomposition.  Figure 1. Example spectrum of NGC 4151 (black) with the model components that were subtracted from the spectrum (red) and the resulting isolated Hβ emission (blue). The vertical dotted lines show the limits of the Hβ emission that was used to constrain the dynamical models in this work. many weak emission lines in the spectrum from Fe II and other species were adequately modeled with single Gaussians. As our goal here is to simply isolate the Hβ emission, we do not attempt to interpret the model components beyond the goodness of fit that they provide to the observed spectra. Our process closely followed that of Bentz et al. (2021), where we first modeled the very high S/N mean spectrum, and then used the best-fit parameters for the mean spectrum as the initial parameters for the model of each individual spectrum. Once a good fit was identified, we subtracted the host-galaxy and power law continuum components, as well as the [O III], He II, and other faint emission lines that were blended with Hβ. Figure 1 displays an example spectrum in black, with the ULySS model components that were ultimately subtracted shown in red and the isolated Hβ spectrum in blue.
After isolating the Hβ emission in each spectrum, we compared all spectra that were collected on a single night of observations. In a few cases, weather conditions changed enough over the course of the observations (∼ 1.0 − 1.5 hours) that the Hβ profile in one spectrum deviated strongly from what was observed in the other spectra collected on that night. We discarded three of the 96 spectra at this stage, and then averaged together the remaining spectra on each night to create 32 nightly spectra. Finally, we cropped the nightly spectra outside of the range 4720 − 4950Å to focus solely on the Hβ emission. A small portion of the red wing of Hβ was excluded by this region because of its position underneath the core of the bright [O III] λ4959 emission line, where our spectral modeling process was sometimes unable to fully separate the two without leaving strong residuals behind.
Finally, we created the continuum light curve from our final sample of 93 re-extracted and scaled spectra by measuring the flux at rest-frame 5100Å. We supplemented this light curve with measurements from several additional sources to significantly extend the time baseline as well as to improve the temporal sampling, when possible, beyond that provided solely by the MDM observations. As described in the original analysis presented by Bentz et al. (2006), a few additional spectra were collected at the Crimean Astrophysical Observatory while the MDM observations were underway, and the continuum flux at 5100Å from the Crimean spectra provided an additional eight measurements. A further 24 measurements from V −band photometry that were collected as part of the MAGNUM project (Koshida et al. 2014) were included, extending the continuum light curve approximately 100 days before the start of the spectroscopic monitoring (although with a coarse temporal cadence). Finally, an additional three V −band measurements were added from photometry collected with the SARA telescopes (Roberts & Rumstay 2012).
All of these sets of measurements adopted different aperture sizes and are thus subject to different aperture losses and the inclusion of different amounts of host galaxy starlight. Furthermore, there are bandpass differences between the spectroscopic measurements and the broad-band photometry. To correct for these differences and calibrate each of the supplemental data sets to match the continuum fluxes measured from MDM spectroscopy, we follow the general procedure outlined by Peterson et al. (1991). We first identified measurements that were made close in time (∆t 1 day) to MDM measurements. For each supplemental data set, we fit a linear function to the close-in-time points from that observatory and from MDM, accounting for the uncertainties in each set of measurements. The linear fit was then used to scale the supplemental dataset so that it matched the MDM measurements, and after all of the supplemental data sets were appropriately scaled to match the MDM measurements, they were merged together. Figure 2 shows the four intercalibrated data sets with measurements from MDM spectroscopy in black, measurements from Crimean Astrophysical Observatory spectra in red, V −band photometry from the MAG-NUM project in blue, and V −band photometry from the SARA telescopes in green. The final continuum light curve was binned so that all measurements collected within ∆t = 0.5 days were averaged together, providing 51 measurements over 143 days, with near-daily sampling during the final 42 days. In comparison, the original continuum light curve analyzed by Bentz et al. (2006) included 37 measurements covering just 41 days. Continuum light curve for NGC 4151 with data sets color coded as follows: spectroscopy from MDM (black), spectroscopy from Crimean Astrophysical Observatory (red), V −band photometry from the MAGNUM project (blue), and V −band photometry from the SARA telescopes (green).
Modeling of the Hβ-emitting BLR was conducted with CARAMEL, a phenomenological modeling code that is described in detail by Pancoast et al. (2014a). CARAMEL explores the BLR geometry and kinematics through the reverberation response across the velocity-resolved profile of a broad emission line as a function of time. Here, we summarize the main components of the model.
The emissivity of the BLR in CARAMEL is represented as a large collection of massless point particles that surround a massive black hole and are distributed in position and velocity space. We note that these points should not be interpreted as physical objects, but rather as a Monte Carlo sampling of the line emissivity. Continuum flux that is incident on a point particle is processed instantaneously, and the distribution of time delays from the BLR depends on the spatial distribution of the point particles, while the velocity distribution of the point particles gives the broad emission-line wavelength profile.
Radial and angular distributions are used to parameterize the spatial distribution of particles. The radial positions of the particles are drawn from a gamma distribution that is flexible enough to represent a Gaussian (α > 1), an exponential (α = 1), or a cuspier profile (0 < α < 1). Experiments with different functional forms have shown that the results are insensitive to the choice of the gamma distribution. The Schwarzschild radius, R s = 2GM/c 2 , plus an additional possible minimum radius r min are used to shift the gamma distribution of particles away from the location of the black hole. A change of variables is performed to assist with interpretation of the modeling results, which are given in terms of (µ, β, where µ is the mean radius, β is the shape parameter, and F is r min in units of µ. The shifted gamma profile has a standard deviation given by σ r = µβ(1 − F ). An outer radius of r out = c∆t data /2 truncates the BLR, where ∆t data is the time difference between the first point in the modeled continuum light curve and the first point in the emission-line light curve. The truncation at r out assumes that the time baseline of the monitoring campaign is sufficiently long to track reverberation signals across the whole BLR. The particle angular distribution is arranged in a disk with a thickness that corresponds to an opening angle θ o , where θ o = 0 • is a thin disk and θ o = 90 • is a sphere. The disk inclination to the line of sight of an observer is set by θ i , where θ i = 0 • is face on and θ i = 90 • is edge on. The distribution of particles as a function of depth within the disk sets the line emission strength. For a single particle, the angle of displacement from the disk midplane is given by where U is a random number that is drawn from a uniform distribution between 0 and 1. The value of γ ranges from 1 to 5, with a value of 1 corresponding to particles that are distributed uniformly throughout the thickness of the disk, while a value of 5 corresponds to clustered particles along the face of the disk, or emission that is preferentially from the outer skin of the BLR. The asymmetry parameter ξ parameterizes the amount of obscuration along the midplane of the disk, where ξ → 0 causes the entire back half of the disk to be obscured and ξ = 1 has no midplane obscuration. Finally κ provides a weight to each particle where W is the fraction of continuum flux that is radiated back towards the observer as line flux and φ gives the angle between the observer's line of sight to the source and the particle's line of sight to the source. The value of κ ranges from −0.5 to 0.5. In the case of κ = −0.5, the particles preferentially emit back towards the ionizing source and an observer would see preferential emission from the far side of the disk. Whereas for κ = 0.5, the particles preferentially radiate away from the ionizing source and and observer would see preferential emission from the near side. The particles are distributed in velocity space through radial and tangential distributions. Some fraction of the particles, f ellip , have near-circular orbits within the Keplerian potential of the central black hole with mass M BH . The remaining particles (1 − f ellip ) are either inflowing (f flow < 0.5) or outflowing (f flow > 0.5). However, these orbits may be highly eccentric and generally bound, or they may be unbound, as determined by the parameter θ e . The possible values of the radial and tangential velocities define a plane, within which θ e describes the angle of the velocity components towards the circular velocity and away from the escape velocity. For θ e = 0 • , the orbits are drawn from a Gaussian distribution centered on the escape velocity. For θ e → 90 • , the inflowing or outflowing orbits approach the parameter space occupied by near-circular orbits. High values of θ e indicate that inflowing or outflowing orbits are very nearly circular, while θ e ≈ 45 • indicates that most of the inflowing or outflowing orbits are highly eccentric but still bound, and low values of θ e correspond to most particles being near the escape velocity and unbound.
Included in the line-of-sight component of the velocity vector for each point particle is a contribution from macroturbulence, given by v turb = N (0, σ turb )|v circ |, where v circ is the circular velocity and N (0, σ turb ) is a normal distribution centered on 0 and with standard deviation σ turb . Parameterization of the spatial and velocity distributions of the particles allows an emission-line profile to be calculated for each continuum flux measurement, assuming that the continuum flux tracks the ionizing flux from a central point source. Included in the modeled emission-line profiles is a nonvariable narrow emissionline component, as well as a smoothing parameter to account for the small remaining differences in spectral resolution that arise from variable seeing conditions throughout the monitoring campaign.
The continuum light curve must be interpolated in order to explore the range of possible time delays arising from the BLR and to properly compare the measured and the modeled emission line profiles. CARAMEL employs Gaussian processes to interpolate between continuum flux measurements as well as to extrapolate the continuum light curve beyond the start and end of the monitoring campaign, thus extending the range of time delays that may be probed with the models. The determination of the BLR model parameters includes the uncertainties on the Gaussian process model parameters, and so captures the effects of interpolating and extrapolating the continuum data within the quoted uncertainties.
For each model realization, 2000 individual point particles are used to represent the BLR. The continuum light curve is interpolated and emission-line profiles models are calculated for each epoch at which an emission-line measurement was acquired. A Gaussian likelihood function compares the modeled spectra with the measured spectra and adjusts the model parameters accordingly. CARAMEL utilizes a diffusive nested sampling code based on DNEST4 (Brewer & Foreman-Mackey 2018) to efficiently explore the parameter space of the models. DNEST4 allows for the use of a likelihood softening parameter, or statistical temperature T , which has the effect of increasing the measurement uncertainties. The likelihood softening parameter is able to account for underestimated measurement uncertainties and for the inability of the model to capture all of the complex and real details in the measurements. After completing 10,000 model runs, the value of T is determined in the post analysis by examining the distributions of the model parameters and choosing the largest value of T for which the distributions remain smooth and generally unimodal. To verify that convergence had been reached, we compared the values of the model parameters from the first half of the model runs to the values determined for the second half of the model runs, finding no significant difference between the parameters constrained from either half.

Model limitations
Exploring the parameter space of a BLR model requires quick and repeated calculations to compute many emission line time-series. Due to this constraint, the model must be simplified, and with these simplifications come limitations. First, the model excludes some physics such as radiation pressure and photoionization processes. Including these would require additional assumptions about the ionizing continuum, which is generally not observable, and the BLR gas properties. Neglecting radiation pressure is standard in reverberation mapping studies as gravity is assumed to be the dominant force affecting BLR kinematics, especially in low Eddington ratio sources like NGC 4151 (L bol /L Edd ∼ 0.01 − 0.1, Merritt 2022). Since we do not include photoionization processes, however, it is critical to understand that the models investigate the BLR emissivity distribution, and extreme care should be taken when using the model results to describe the BLR gas. For instance, a model with γ = 5 indicates that the observed emissivity is concen-   trated at the "skin" of the thick BLR disk. This could mean there is no gas in the inner parts of the disk or that the emission from the inner parts is obscured, but our models cannot distinguish between the two scenarios. Second, the model is parametrized in such a way that the maximum range of geometries can be described with the fewest number of free parameters. This requires the use of smooth, continuous functions to describe the BLR emissivity. Thus, short-timescale fluctuations in the emission line profile, corresponding to short sizescale fluctuations in the BLR, cannot be modeled.
Third, because photoionization processes are not included, the CARAMEL model does not fit the absolute flux scale of the emission line but rather re-scales the continuum fluctuations so that the two scales match. Modeling the absolute fluxes would again require knowledge of the ionizing continuum as well as the physical properties of the BLR gas. By re-scaling the fluxes in this way, the emission line profile shapes can be fit and provide constraints on the BLR without making additional assumptions.
Despite these limitations, proper interpretation of the model-as a Monte Carlo approximation of the BLR kinematics and emissivity field-provides significant information about properties of the BLR that are of the most interest, such as the size, orientation, overall structure and kinematics, as well as the black hole mass. Repeat modeling of the same AGN over multiple observing campaigns demonstrates that parameters expected to remain constant are robust (e.g., Arp 151, Pancoast et al. 2018). Modeling of multiple emission lines from the same AGN find ionization stratification consistent with theory (NGC 5548, Williams et al. 2020;NGC 3783, Bentz et al. 2021). Constraints on the BLR kinematics are consistent with those found from velocity-resolved RM and the maximum entropy method (e.g., Villafaña et al. 2022a), and inclination angles are consistent in the few cases in which independent measurements are available (Grier et al. 2017;this work). Tests of simulated data also show reassuring results, in that the key properties of the BLR and the black hole mass are robustly recovered even when the input model is significantly different than what is assumed in CARAMEL (Mangham et al. 2019).

RESULTS
The median and 68% confidence intervals for all of the Hβ BLR model parameters in NGC 4151 are listed in Table 1 √ T = 44.7. We note that the very high S/N in each epoch of spectroscopy considered here suggests that the statistical noise is negligible compared to modeling and systematic errors, thus requiring T >> 1. −0.63 lt-day, respectively. The radial width of the emission is σ r = 5.74 +1.64 −1.19 lt-day, and the radial distribution of the emission has a profile that is slightly more cuspy than an exponential (β = 1.18 +0.17 −0.14 ). The emission is distributed fairly uniformly through the disk, with a slight preference for stronger emission near the face of the disk (γ = 1.67 +0.98 −0.47 ) and almost complete obscuration along the midplane (ξ = 0.10 +0.13 −0.07 ). Most of the line emission is preferentially directed back towards the central illuminating source (κ = −0.36 +0.06 −0.06 ), with the observer seeing a strong response from the far side of the disk and weak response from the front. Figure 5 displays a representative geometric model, drawn from the posterior probability distribution for the Hβ emission-line response in NGC 4151.

DISCUSSION
Previous studies that model velocity-resolved reverberation mapping data have consistently found a preference for BLR geometries that resemble a moderatelyinclined thick disk Pancoast et al. 2014b;Grier et al. 2017;Williams et al. 2018;Bentz et al. 2021;Villafaña et al. 2022a). The opening angle and inclination angle that we find here suggest that NGC 4151 has the thickest and most highly inclined BLR disk among Seyfert 1s that have been studied with these methods (Villafaña et al. 2022b).
The inclination angle to our line of sight, θ i = 58.1 +8.4 −9.6 deg, is consistent with the value of θ i = 45 ± 5 • derived through geometric modeling of the narrow-line region as a bicone (Das et al. 2005;Fischer et al. 2013). And with an opening angle of θ o = 56.6 +15.8 −14.3 deg, the models suggest that, as observers, we are just able to peer over the edge of the BLR structure and view the innermost central engine. This interpretation is supported by observations of X-ray variability in NGC 4151 that seem to be caused by eclipsing material associated with the BLR traversing our line of sight (Puccetti et al. 2007;Wang et al. 2010).
In Table 2 we list all published direct black hole mass measurements for NGC 4151, including the black hole mass that we present here from modeling of velocityresolved RM data of log M BH /M = 7.22 +0.11 −0.10 or M BH = 1.66 +0.48 −0.34 × 10 7 M . Roberts et al. (2021) carried out a new stellar dynamical modeling analysis of the nuclear stellar kinematics presented by Onken et al. (2014) and determined that the models preferred a black hole mass of M BH = 0.25 − 3.0 × 10 7 M , which agrees very well with the value we present here. The mass based on gas dynamical modeling presented by Hicks & Malkan (2008) was adjusted to a galaxy distance of D = 15.8 ± 0.4 Mpc as determined from HST observations of Cepheid stars in NGC 4151 (Yuan et al. 2020), and also agrees with the mass we present here as well as the mass from stellar dynamical modeling.
Reverberation analyses of NGC 4151 that adopt a population-average scale factor of f = 4.82 ± 1.67 (Batiste et al. 2017) and are listed in Table 2 as 'Hβ RM', somewhat overestimate the black hole mass. Based on the inclination angles preferred by the best-fit BLR models, this effect is to be expected. Values of f range from 2.8 (Graham et al. 2011) to 5.5 (Onken et al. 2004 in the literature depending on the exact sample and the analysis methods employed, with most investigations settling on values of ∼ 4 − 5. Given the factor of 1/ sin θ i between the observed velocities along the line of sight and the true velocities, values of f ≈ 4 − 5 suggest that most Seyferts in the reverberation sample are viewed at inclinations of 25 • − 30 • . An inclination angle of ∼ 58 • , such as the models suggest for NGC 4151, would require a smaller than average scale factor to accurately calibrate the black hole mass. With the black hole mass constraints from the models presented here and the mean time delay and line width for Hβ that were reported by Bentz et al. (2006), we can infer a specific value of f = 1.8 +0. −9.6 deg. A similar analysis by Villafaña et al. (2022b) with an expanded sample, including the results we present here for NGC 4151, predicts the same individual scale factor but with somewhat smaller uncertainties, log f = −0.44 +0.73 −0.74 or f = 0.36 +1.57 −0.29 . Both studies agree with the value that we infer here. We note that if a scale factor of f = 1.8 is adopted for the reverberation analysis of De Rosa et al. (2018), the derived mass is somewhat lower than we find here, but agrees within 2σ. It is likely that the scale factor associated with a specific AGN may change as a function of time as the detailed geometry and kinematics in the BLR change on a dynamical timescale (a few years at the location of the Hβ-emitting BLR in Seyferts) and respond to large-scale changes in the ionizing flux. With 7 years passing between the 2005 observational campaign analyzed here and the 2012 program presented by De Rosa et al. (2018), BLR structural changes would not be unexpected. Modeling of the 2012 data is currently in progress and may provide interesting insights into the time evolution of the structure of the BLR in NGC 4151 (Robinson et al., in prep).
Finally, with the successful launch of JWST in December 2021, new observations of the nuclear stellar kinematics in NGC 4151 will soon be collected with NIRSpec as part of an Early Release Science program (ERS 1364, PI Bentz). NIRSpec is expected to provide some crucial advantages over AO-assisted ground-based observations, such as the NIFS observations analyzed by Roberts et al. (2021). With its stable and diffraction-limited PSF and significantly lower backgrounds, stellar kinematics mea-sured with NIRSpec may allow for tighter constraints to be placed on the black hole mass through stellar dynamical modeling, allowing further exploration of the consistency in masses derived from independent black hole mass measurement techniques.

SUMMARY
We have reanalyzed the 2005 monitoring data for NGC 4151, supplemented with additional measurements from the literature, and carried out an exploration of models for the full velocity-resolved BLR response to continuum variations. The modeling results find that the BLR is well represented by a very thick disk with an opening angle (θ o ≈ 57 • ) that is similar to the inclination angle (θ i ≈ 58 • ). The inclination angle is consistent with the value derived from bicone modeling of the narrow line region, and the similarity of the opening angle and inclination angle suggests that our line of sight to the innermost central engine is just barely free of obstruction from the BLR and dusty torus. This is an intriguing consideration since previous studies suggest that BLR gas has been observed to temporarily eclipse the central X-ray source along our sight line to NGC 4151. The kinematics of the BLR gas are found to be dominated by eccentric but bound orbits, with ∼ 10% of the orbits showing a preference for nearcircular motions. With the geometry and kinematics of the BLR constrained, the models provide an independent and direct black hole mass measurement of log M BH /M = 7.22 +0.11 −0.10 or M BH = 1.66 +0.48 −0.34 ×10 7 M , which is in excellent agreement with black hole masses determined from stellar dynamical modeling and gas dynamical modeling of NGC 4151.