Velocity-resolved Reverberation Mapping of NGC 3227

We describe the results of a new reverberation mapping program focused on the nearby Seyfert galaxy NGC 3227. Photometric and spectroscopic monitoring was carried out from 2022 December to 2023 June with the Las Cumbres Observatory network of telescopes. We detected time delays in several optical broad emission lines, with Hβ having the longest delay at τcent=4.0−0.9+0.9 days and He ii having the shortest delay with τcent=0.9−0.8+1.1 days. We also detect velocity-resolved behavior of the Hβ emission line, with different line-of-sight velocities corresponding to different observed time delays. Combining the integrated Hβ time delay with the width of the variable component of the emission line and a standard scale factor suggests a black hole mass of MBH=1.1−0.3+0.2×107 M ⊙. Modeling of the full velocity-resolved response of the Hβ emission line with the phenomenological code CARAMEL finds a similar mass of MBH=1.2−0.7+1.5×107 M ⊙ and suggests that the Hβ-emitting broad-line region (BLR) may be represented by a biconical or flared disk structure that we are viewing at an inclination angle of θ i ≈ 33° and with gas motions that are dominated by rotation. The new photoionization-based BLR modeling tool BELMAC finds general agreement with the observations when assuming the best-fit CARAMEL results; however, BELMAC prefers a thick-disk geometry and kinematics that are equally composed of rotation and inflow. Both codes infer a radially extended and flattened BLR that is not outflowing.


INTRODUCTION
The prevalence of supermassive black holes with masses of 10 6 M BH /M ⊙ 10 10 in the centers of massive galaxies was one of the surprise discoveries enabled by the Hubble Space Telescope (e.g., Ford et al. 1994;Harms et al. 1994;Kormendy & Richstone 1995;van der Marel et al. 1997).More surprising still was the finding that the masses of these black holes scale with observable properties of their host galaxies (e.g., bentz@astro.gsu.edu* Packard Fellow Ferrarese & Merritt 2000;Gebhardt et al. 2000).After three decades of study, we now understand that these black holes generally go through at least one period of active accretion during their lifetimes, and that the accretion process releases huge amounts of energy across the electromagnetic spectrum, energy that is deposited into the gas within the host galaxy and circumgalactic environment and thus seems to play a role in modulating the growth of the galaxy and of the black hole itself (see the reviews by, e.g., Fabian 2012;Heckman & Best 2014;Fan et al. 2023).
The mass of a black hole is one of its fundamental properties, and it sets the limit for the amount of feedback power that may be released during an accretion event.Thus it is an important measurement to constrain as we attempt to further understand the detailed physics involved in black hole feeding and feedback.Sagittarius A * is the nearest supermassive black hole, and its mass has been precisely determined through monitoring the orbits of individual stars via their proper motions over decade timescales (Ghez et al. 2000;Genzel et al. 2000;Ghez et al. 2008).Unfortunately, all other galaxies are so distant that we cannot achieve the necessary spatial resolution in their nuclei to apply the same technique.
A variety of techniques have been developed to determine black hole masses in other galaxies.These include dynamical modeling of the bulk motions of stars or gas (see the review by Kormendy & Ho 2013) and kinematic studies of masing gas clouds (e.g., Miyoshi et al. 1995).High spatial resolution is a critical component of these techniques, which generally limits their application to galaxies with D 100 Mpc.
Black hole masses at cosmological distances are most often determined through reverberation mapping (RM; see the recent review by Cackett et al. 2021), a technique that makes use of the photoionized gas deep in the potential well of an actively accreting black hole.RM utilizes spectrophotometric monitoring with high temporal sampling to track changes in the continuum emission (expected to arise from the accretion disk) and the "echoes" of those changes in the broad emission-line fluxes.The time delay between the two gives the typical light travel time across the broad line region, and thus a physical size.With reverberation mapping, spatial resolution is replaced by temporal resolution, and black holes at significant distances become available for study.
One effect of having these varied techniques and their disparate limitations is that we have two black hole mass scales currently in use: one determined from nearby, mostly early-type inactive galaxies based on stellar and gas dynamical modeling, and another based on reverberation mapping of more distant active galaxies that tend to be mostly later morphological types.Because broad-lined active galaxies with D 100 Mpc are rare, there are very few galaxies that meet the requirements for multiple techniques that would allow us to check that these two mass scales are in agreement.
NGC 3227 is a nearby (z = 0.00377) spiral galaxy that is interacting with its elliptical companion, NGC 3226.Both galaxies in the pair show evidence for nuclear activity, a fact that was first remarked upon by Curtis (1918) who noted that the pair of galaxies contained "stellar nuclei".NGC 3226 is classified as a LINER (Véron-Cetty & Véron 2006), while NGC 3227 was one of the original sample of 12 Seyfert (1943) galaxies.The broad hydrogen lines in the nucleus of NGC 3227, compared to the widths of the nebular lines, was noted by Dibai & Pronik (1968) and Rubin & Ford (1968), though other early studies have sometimes classified it as a Type II Seyfert (e.g., Khachikian & Weedman 1974).
Similarly, early studies were divided regarding the potential detection of nuclear variability in NGC 3227 (e.g., Netzer 1974;Liutyi 1977).While Tennant & Mushotzky (1983) detected X-ray variability on timescales of ∼ 1 day in NGC 3227, Salamanca et al. (1994) were the first to systematically monitor the AGN with optical spectroscopy on short timescales and demonstrate clear variability in the continuum and broad emission lines.With only 26 measurements over the course of 6 months, however, their light curves were severely undersampled when it came to trying to constrain any time delays between them.A similar campaign with comparable time sampling described by Winge et al. (1995) was also plagued by the same difficulties.A re-analysis of both datasets by Peterson et al. (2004) using a set of best practices that had been developed for reverberation mapping studies concluded that any potential time delays were formally consistent with zero because of the coarse temporal sampling.Denney et al. (2010) thus undertook a new monitoring program with the intent of determining the first accurate and precise BLR time delay measurement for NGC 3227, finding that the Hβ broad line emission lagged the continuum variations by 3.8 ± 0.8 days and showed evidence for different time delays as a function of velocity across the emission line profile.Further campaigns by De Rosa et al. (2018) and Brotherton et al. (2020) sought to improve on these results by more fully exploring the velocity-resolved time delays to determine the BLR geometry and kinematics, in the hopes of deriving a more direct constraint on the black hole mass.Preliminary modeling of the De Rosa et al. (2018) data sets suggest that the low amplitude of variability in the light curves may hinder any firm conclusions from being drawn (Robinson 2022), and while Brotherton et al. (2020) provide a velocity-delay map for Hβ based on their monitoring data, the lack of details in the map limited any potential for interpretation.
We thus undertook a new monitoring program focused on NGC 3227, which we present in this work.Strong variability coupled with high temporal sampling allowed us to measure broad emission-line time delays and to further explore the velocity-resolved response of Hβ.These promising characteristics also enabled a more thorough exploration of the constraints that may be placed on the structure and kinematics of the BLR in this nearby AGN, based on two independent modeling codes that have been developed specifically for use with reverberation mapping data.

OBSERVATIONS
NGC 3227 is located in the direction of the constellation Leo at α = 10h23m30.6sand δ = +19 • 51 ′ 54 ′′ with a redshift of z = 0.00377.A surface brightness fluctuations (SBF) distance to the companion galaxy NGC 3226 puts the pair at a distance of D = 23.7 ± 2.6 Mpc (Tonry et al. 2001 with the adjustments of Blakeslee et al. 2001).Rather than rely on aperture photometry, which may include large amounts of host-galaxy starlight that varies with the seeing and also dampens the true brightness variations of the central AGN, we instead used the image subtraction algorithms of Alard & Lupton (1998); Alard (2000).In general, image subtraction works by scaling individual images to match a reference image and then subtracting the two, removing all sources of constant flux and leaving behind only the variable flux.Image subtraction methods are commonly applied in studies of extragalactic transients and other time domain phenomena, including microlensing (Udalski et al. 2008(Udalski et al. , 2015)), supernovae (Riess et al. 2001;Miknaitis et al. 2007;Melinder et al. 2008), tidal disruption events (Holoien et al. 2016;Brown et al. 2018), as well as AGN variability (Fausnaugh et al. 2016;Grier et al. 2017;Bentz et al. 2021a).The inner 10 ′ × 10 ′ region of the reference image, centered on the nucleus of NGC 3227 and oriented with North up and East to the left.
We began by registering all of the images to a common grid using the algorithm of Siverd et al. (2012).We then selected the highest quality images (lowest background and best seeing) from the subset that were obtained at SSO, and stacked them to build a reference image (see 1).The reference image was convolved to match, and then subtracted from, each of the 376 individual images.Finally, aperture photometry with a radius of 9 pixels was performed at the location of the AGN in all of the subtracted images, giving us a light curve of the variable AGN counts.
While the LCO network employs identical telescopes and detectors across sites, we have previously found slight offsets between the measurements obtained at different LCO locations (Bentz et al. 2021a(Bentz et al. , 2023)).Accordingly, we compared the light curves obtained at each location to the light curve obtained at SSO and again found that, when working in counts, small adjustments were necessary to bring them all into agreement.For a pair of light curves, we identified the set of measurements that were taken close together in time (usually within 0.5 − 1 days) and determined the slope and intercept of the relationship between the measurements from one observatory versus the other, taking into account the errors in both measurements.Good agreement would be indicated by a slope of 1.0 and intercept of 0.0, with some scatter given the small and random variations of the AGN on very short timescales.In all comparisons, the slope and intercept were found to deviate somewhat from these expected values, demonstrating a difference in the counts that are registered on the detectors at different observatory locations, likely due to factors like slight variations in the overall sensitivities of the CCDs or differences in atmospheric transparency.We applied a small additive and multiplicative correction factor to each light curve to bring them into agreement with the SSO light curve, which was the basis for the reference image.
We then calibrated the light curves using the AGN flux in the reference image.To separate the AGN from its bright and spatially extended host galaxy, we employed the two-dimensional surface brightness modeling package Galfit (Peng et al. 2002(Peng et al. , 2010)).We first constructed a PSF model for the reference image by fitting multiple Gaussians to a small cutout centered on an isolated field star, with the background sky brightness modeled as a gradient in x and y.Four Gaussians were found to be sufficient to capture the detailed shape of the PSF.We then modeled the reference image using this PSF to simulate the AGN and several field stars, a gradient in x and y for the background sky component, and multiple Sérsic profiles for the galaxies NGC 3227 and NGC 3226.Our choice of starting parameters for the Sérsic profiles was guided by previous surface brightness models of NGC 3227 using highresolution Hubble Space Telescope images (Bentz et al. 2009;Bentz & Manne-Nicholas 2018).Once we had a good model that allowed us to cleanly separate the AGN flux from the surrounding host galaxy, we determined the final calibration of the reference frame by setting the zeropoint such that the recovered V −band magnitudes for several field stars matched the values tabulated by the AAVSO Photometric All Sky Survey catalog (Henden & Munari 2014) DR10.The calibrated V −band magnitude for the AGN in the reference image was then combined with the variable flux measured from image subtraction to determine the calibrated V −band light curve.
We next examined the residual counts around several non-varying field stars in the subtracted images and compared these measurements to the uncertainties reported on the aperture photometry at the location of the AGN.Previous work by Zebrun et al. (2001); Hartman et al. (2005) and others has demonstrated that the uncertainties in image subtraction measurements are often underestimated, and we found that to hold true in the case of our analysis of NGC 3227.We followed the basic procedure outlined by Hartman et al. (2004) to determine that a multiplicative factor of ∼ 9 needed to be applied to provide more realistic uncertainties.
Finally, we examined the effectiveness of binning the light curve for decreasing some of the noise that may re- sult from images that were obtained under marginal observing conditions.We found that a bin size of 0.5 days was a good compromise between retaining much of the temporal sampling while also decreasing some of the noise.In Figure 2, we show the final calibrated and binned V −band light curve for the central AGN in NGC 3227 during this monitoring campaign.The bottom panel shows the scatter in the residuals for a nonvarying field star of similar magnitude for comparison with the AGN variability.The light curve is tabulated in Table 2, with flux densities based on the absolute spectrophotometric calibration of Vega provided by Colina et al. (1996).

Spectroscopy
Spectroscopic monitoring was carried out on the LCO 2m telescope network between UT dates 2023 February 02 and 2023 June 01 under programs LCO-2023A-004 and ANU-2023A-002.The FLOYDS crossdispersed spectrographs provided wavelength coverage of 5400 Å−1.0 µm with a dispersion of 3.51 Å/pix in the first order, and wavelength coverage of 3200 − 5700 Å in the second order with a dispersion of 1.74 Å/pix.Spectra were obtained through a 6 ′′ slit oriented N-S along a position angle of 0 • .Observations were scheduled every ∼ 24 hr, and a typical visit consisted of a 900 s science exposure together with an arc lamp and a flat field taken through the same slit.
A total of 51 spectra were acquired over the course of the program, with 42 acquired by Faulkes Telescope South and 9 by Faulkes Telescope North.Observations were downloaded from the LCO archive after being processed by the LCO pipeline, which splits the 2D spectra into separate orders, rectifies each order, and applies typical CCD reductions including biases, darks, and flats.The pipeline also provides initial wavelength and flux calibrations based on stored wavelength and sensitivity functions.Focusing on the 2nd order spectra, which have a higher spectral resolution and also include the [O III] λλ4959, 5007 Å doublet (see discussion below), we began with the wavelength-and flux-calibrated 2D spectra and carried out cosmic ray cleaning, spectral extraction through a 10 pix extraction width, and improved wavelength calibration using the arcs taken at the same pointing as each spectrum.We then adopted a common dispersion of 1.6 Å for all spectra and cropped them to have the same starting and ending wavelengths.
Next, we applied the spectral scaling method of van Groningen & Wanders (1992) to account for slight variations in the spectra caused by imprecise flux calibration, varying weather and seeing conditions, and slight differences in telescope pointing.The algorithm applies small shifts in wavelength and flux and a small amount of smoothing to minimize the differences in each spectrum relative to a reference spectrum, where the reference spectrum is generally created from a subset of the best-quality spectra.We applied the algorithm focusing on the spectral region around the [O III] λλ4959, 5007 Å emission lines, as they are known to not vary on the short timescales that are covered by a single observing season (Peterson et al. 2013).This allows the [O III] doublet to be used as a set of internal comparison sources and also explains why most reverberation mapping experiments have focused on the Hβ emission line: even though Hα is brighter, there are no similar strong and relatively unblended narrow lines nearby to use for correction of the slight differences in the spectra from night to night.We adopted an absolute flux calibration of F (5007) = 7.8 × 10 −13 erg s −1 cm −2 , which is based on the values determined for NGC 3227 over the last decade through a wide slit on photometric nights, as reported by De Rosa et al. ( 2018) and Brotherton et al. (2020).
Figure 3 displays the mean and root-mean-square (rms) of all the final calibrated and scaled spectra.The rms spectrum highlights the variable spectral components, which include the AGN continuum and broad emission lines, while non-varying components like the narrow lines are eliminated.

Line Width Measurements
For each of the broad emission lines visible in the spectra of NGC 3227, we report the line widths in both the mean spectrum and the rms spectrum.Line widths were determined by setting a local linear continuum under each line and then measuring the distribution of flux directly from the data above the continuum within a specified wavelength window.Line widths were quantified in two ways: the full width at half-maximum flux (FWHM) and the velocity dispersion of the line profile (σ line ).The uncertainties in the line widths were determined by selecting a random subset of spectra, building a new mean and rms spectrum from this subset, recording the FWHM and σ line , and repeating the process 1000 times to build up distributions of the measurements.In Table 1, we report the median and inner 68% confidence interval as the line width and its associated uncertainties, for both FWHM and σ line as measured in the mean spectrum and the rms spectrum.
We note that all reported line widths have been corrected for instrumental resolution by assuming that ∆λ 2 obs ≈ ∆λ 2 true + ∆λ 2 disp (Peterson et al. 2004).The width of the [O III] λ5007 Å emission line was reported to be FWHM = 485 km s −1 when measured through a narrow slit at high spectral resolution by Whittle (1992), and we take this measurement to represent ∆λ true .The measured width of [O III] λ5007 Å in our spectra (∆λ obs ) is FWHM = 16.58Å.We thus adopted a resolution correction of ∆λ disp = 14.47 Å or 866 km s −1 .Note-Heliocentric Julian dates are provided as HJD−2400000 (days).V −band flux densities have units of 10 −15 erg s −1 cm −2 Å−1 while emission-line fluxes have units of 10 −13 erg s −1 cm −2 .Table 2 is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.Light curves of the broad emission lines Hβ, He II, Hγ, and Hδ were generated by fitting a local linear continuum under each line and integrating the flux above this continuum for each calibrated and scaled spectrum.While simplistic, this method accounts for all of the flux in a line without having to rely on the discretion of the user in choosing an appropriate set of model parameters to properly fit asymmetric or complicated line shapes.The non-varying narrow component of each line is included as a simple flux offset.Table 2 tabulates the light curves.Table 3 lists basic information and variability statistics about each light curve; for each spectral feature listed in column (1), we give the number of measurements in column (2), and the average and median temporal sampling, respectively, in columns (3) and (4).Column (5) lists the mean flux and standard deviation, while column (6) lists the mean fractional error.Column (7) gives the noise-corrected fractional variation, which is computed as

Emission-Line Light Curves and Time Delays
where σ 2 is the variance of the fluxes, δ 2 is their meansquare uncertainty, and F is the mean flux.Column (8) lists the ratio of the maximum to the minimum flux.
The typical time delay between each emission-line light curve and the continuum light curve was then assessed.Throughout this work, we adopt the final V −band light curve from photometry as the continuum light curve because it has several advantages over the continuum as measured from the spectra.In particular, the use of broad-band photometry, even on a small telescope, allows a high S/N to be reached in a short amount of observing time.The wide field of view of the detectors available on small telescopes provide ample field stars that are observed simultaneously with the science target, providing a more accurate calibration than can generally be achieved for ground-based long-slit spectroscopy.The abundance of small optical telescopes around the globe allows a significantly higher temporal sampling to be achieved, and the image subtraction methods employed in our analysis produce a high-fidelity light curve without the damping effect of nonvariable flux components.Finally, the V band was chosen because of its widespread availability and because it includes a mostly line-free portion of the spectrum for local AGNs like NGC 3227.
Time delays were measured using the interpolated cross-correlation function (ICCF) method of Gaskell & Sparke (1986); Gaskell & Peterson (1987) with the modifications of White & Peterson (1994).The cross-correlation function is determined twice with the ICCF method, first by interpolating one light curve and then the other.The two resultant CCFs are averaged together to provide the final CCF, and we display these in Figure 4 along with the continuum and emission-line light curves.
The peak of the CCF (r max ) occurs at the time shift for which the two light curves are most highly correlated.This time delay is often reported as τ peak .A slight variation is to report τ cent , which is the time delay associated with the centroid of the points around τ peak above some threshold (typically 0.8r max ).We prefer τ cent because it is less susceptible to noise than is τ peak (cf.Peterson et al. 2004).
The uncertainties on τ peak and τ cent were derived following the flux randomization and random subset sampling (FR/RSS) method of Peterson et al. (1998Peterson et al. ( , 2004)).The FR/RSS method accounts for uncertainties in the time delay arising from noise in the flux measurements as well as the inclusion (or exclusion) of each individual data point.A large number of FR/RSS realizations (N = 1000) build up distributions of τ peak and τ cent measurements.From these distributions, we report the medians and inner 68% confidence intervals as the time delay values and their associated uncertainties.Time delays in the observed frame are tabulated in Table 1.
We also measured the emission-line time delays using JAVELIN (Zu et al. 2011).JAVELIN assumes a damped random walk for the behavior of the driving continuum light curve, and optimizes the parameters for a top-hat transfer function to predict the responding emission-line light curve(s).Even though JAVELIN is capable of fitting multiple emission-line light curves at the same time, we were unsuccessful in fitting them all simultaneously, so we instead fit each one separately.The best-fit time delays are listed in Table 1, and while they agree well with τ peak , they generally do not agree with τ cent .Differences between τ peak and τ cent are caused by the shape of the CCF, which is in turn related to the transfer function.We speculate that in this case, JAVELIN's assumption of a top hat may not be a good match to the actual transfer function of NGC 3227.
Using the same methods as described above, we also divided the Hβ broad emission line into 5 separate bins, with similar flux in all bins, and generated a light curve for each one to assess the Hβ gas response as a function of line-of-sight velocity.The time delay between the continuum and the light curve associated with each bin was then determined through cross correlation following the same methods outlined above.As we show in Figure 5, there are clear differences in the time delay associated with each bin across the broad emissionline profile, with the longest time delays near the line center and the shortest time delay associated with the blue wing.This velocity-resolved behavior is qualitatively similar to what has been observed in previous velocity-resolved analyses of NGC 3227 (Denney et al. 2010;De Rosa et al. 2018;Brotherton et al. 2020).

CARAMEL Modeling
Using the full velocity-resolved response of an emission line to constrain the geometry and kinematics of the BLR is usually approached in two ways: through the inverse approach, in which deconvolution of the emission-line response produces a velocitydelay map (e.g., Horne et al. 2004;Skielboe et al. 2015;Anderson et al. 2021), or through forward modeling, in which a self-consistent set of models are used to predict the emission-line response and then compared to the observations (e.g., Pancoast et al. 2014;Rosborough et al. 2023).The first approach has the advantage of relying on a relatively small set of assumptions and has the ability to make use of all the details in the data, but short baselines, irregular sampling, and noisy data complicate the deconvolution process, and the results are difficult to interpret without extensive modeling.The second approach is relatively easy to interpret, but is fundamentally limited by a larger set of core assumptions and the level of flexibility and completeness in the models.
CARAMEL is a phenomenological modeling code that was developed to explore the BLR using reverberation data sets.It is described in detail by Pancoast et al. (2014) and we have previously used it in the study of other nearby AGNs (Bentz et al. 2021b(Bentz et al. , 2022(Bentz et al. , 2023)).By exploring the parameter space of geometric and kinematic models that may represent the BLR, CARAMEL is able to provide a direct and primary constraint on the black hole mass (without resorting to the use of a scale factor) and can also provide insight into various characteristics of interest, such as the inclination of the BLR to our line of sight.In the few cases where independent measurements are available, they have been to found to be in general agreement with CARAMEL results, such as the mass of the black hole and the BLR inclination in NGC 4151 (Bentz et al. 2022), and the general structural parameters of the BLR in NGC 3783 (Gravity Collaboration et al. 2021;Bentz et al. 2021b).
We followed the approach described in our previous work (Bentz et al. 2021b(Bentz et al. , 2022(Bentz et al. , 2023) ) and first carried out a decomposition of each spectrum using the ULySS package (Koleva et al. 2011).We modeled all of the continuum and emission-line features in the spectra and then subtracted the continuum components and any emission lines close in wavelength to Hβ, leaving an isolated Hβ feature corresponding to each observational epoch.These isolated observed Hβ profiles were then fed into CARAMEL for comparison with the line profiles created from each model, along with the driving continuum light curve.Note-Tabulated values are the median and 68% confidence intervals.In Figure 6 we display the posterior probability distribution function for the model parameters preferred by CARAMEL, and we tabulate the median and 68% confidence interval for each parameter in Table 4.The median values of these parameters indicate that the Hβemitting region of the BLR may be represented by the surface or "skin" of a biconical or flared disk structure with its axis of symmetry inclined to our line of sight at an angle of θ i = 33.2+13.5 −9.1 degrees, with preferential radiation back towards the central source (κ = −0.28+0.19  −0.15 ), and with gas motions that are dominated by rotation (f ellip = 0.71 +0.14  −0.37 ).We include an illustrative example of this geometry in Figure 7, and note that strong inward anisotropy of hydrogen emission from the BLR is expected due to the high densities and correspondingly large line optical depths in the gas (Ferland et al. 1992).(Brotherton et al. 2020).The time delay of Hβ is expected to change as a function of the central AGN luminosity.The relationship for these two quantities follows a power-law shape with a slope of ≈ 0.5 for a sample of local Seyferts (the R BLR −L relationship; Bentz et al. 2006Bentz et al. , 2009Bentz et al. , 2013)), and each individual AGN is expected to have its own R BLR − L relationship as its luminosity varies (e.g., Bentz et al. 2007;Kilerci Eser et al. 2015), so we would not expect to always measure the same time delay for a specific emission line in a single AGN.
On the other hand, the two Hβ time delay measurements reported by De Rosa et al. ( 2018) lie below the R BLR −L relationship, while the Hβ time delay reported by Brotherton et al. (2020) lies above it, even after ensuring uniformity in the assumed distances, Galactic extinction, and method of starlight correction when deriving the AGN luminosity.The longer time delay reported by Brotherton et al. (2020), when compared to the measurements presented here and by Denney et al. (2010), seems to tentatively agree with the expectation of a steeper slope for the R BLR −L relationship of a single AGN when the optical luminosity is used as a proxy for the ionizing luminosity (Kilerci Eser et al. 2015).However, the AGN luminosity was quite similar during the monitoring program of Brotherton et al. (2020) and the programs described by De Rosa et al. ( 2018), yet the time delays that were derived are quite different.We note that the amplitude of variability relative to the noise was significantly lower in the light curves presented by De Rosa et al. (2018).If the measurements of τ cent and their larger uncertainties are preferred over those of τ jav , then the measurements of De Rosa et al. ( 2018) are within the typical scatter around the R BLR − L relationship, and they agree with the time delay measurement reported by Brotherton et al. (2020) at the 3σ level.

Black Hole Mass
Because of its proximity, NGC 3227 is one of the few Seyfert galaxies that has previous dynamical mass measurements.We note that dynamical masses are linearly dependent on the adopted distance to the galaxy, so any comparison of masses from different techniques must keep this detail in mind.NGC 3227 is near enough that its redshift of z = 0.00377 cannot be used to accurately predict its distance.Instead, throughout this work we have adopted the SBF-based distance to NGC 3226 as the distance to NGC 3227 (given their clear ongoing interaction).The SBF distance of D = 23.7 ± 2.6 Mpc is based on the analysis of Tonry et al. (2001), with slight adjustments by Blakeslee et al. (2001).In addition, a careful assessment of the predicted distance based on the Tully & Fisher (1977) relationship by Robinson et al. (2021) finds D = 24.3 ± 4.9 Mpc for NGC 3227, which is in good agreement.
Previous mass constraints from stellar dynamical modeling (Davies et al. 2006) and from gas dynamical modeling (Hicks & Malkan 2008) used different distances from our adopted value and from each other, both relying on various measurements of the redshift instead.Adjusting the dynamical masses to account for their assumed distances relative to our adopted distance gives a stellar dynamical mass of M BH = (1.9 ± 0.9) × 10 7 M ⊙ and a gas dynamical mass of M BH = (3.1 +1.5 −0.6 )×10 7 M ⊙ .While these values agree with each other, they have generally not agreed with reverberation-based masses, being a factor of 4 − 5 larger (Denney et al. 2010;De Rosa et al. 2018).
For most reverberation analyses, black hole mass is calculated from the time delay (τ ) and line width (V ) as: where c is the speed of light, and G is the gravitational constant.An extensive analysis of reverberation datasets and analysis procedures by Peterson et al. (2004) provided a list of best practices, particularly recommending the use of τ cent for the time delay and σ line (rms) for the line width.The scale factor f encodes all the geometric and kinematic details of the BLR, including the inclination angle to our line of sight.Since these details are generally unknown for individual AGNs, a population-average scale factor, f , is usually adopted so that the M BH − σ ⋆ relationship for AGNs with masses from reverberation is brought into general agreement with the relationship for mostly quiescent galaxies with black hole masses from dynamical modeling (e.g., Onken et al. 2004).For the combination of τ cent and σ line (rms), the value of f has been found to range from 2.8 (Graham et al. 2011) to 5.5 (Onken et al. 2004), depending on the exact details and the sample that is included.We generally prefer to adopt f = 4.8 from Batiste et al. (2017) because of their careful treatment of the effects of galaxy morphology, including bars, on the determination of σ ⋆ .While in practice, the adoption of f minimizes any bias in the full sample of reverberation masses, it also means that the reverberation mass for any particular AGN may be over/underestimated by a factor of a few.If we assume that the variance in f values between different AGNs is mainly due to their random inclinations to our line of sight (as suggested by the analysis of Williams et al. 2018), then f ≈ 5 implies a typical inclination of 27 • .For an AGN that is observed at a lower (more face-on) inclination, f ≈ 5 would significantly underestimate the correction needed to convert the line-of-sight velocity component to the true velocity when determining M BH .And in the case of NGC 3227, a factor of ∼ 4 discrepancy between the reverberation mass and the dynamical masses could be resolved if the AGN were found to be viewed at an angle of ∼ 15 • .
Using the measurements for Hβ listed in Table 1 indicates M BH = 1.07 +0.23  −0.25 × 10 7 M ⊙ .This is somewhat larger than the masses that were determined by Denney et al. (2010) andDe Rosa et al. (2018), but the 1σ range formally agrees with the low end of the uncertainties on the much higher mass reported by Brotherton et al. 2020.The time delay for Hβ found in this work and the analysis of Denney et al. (2010) are nearly identical, however, somewhat perplexingly, while the FWHM measurements of Hβ agree within the uncertainties, σ line is ∼ 300 km s −1 (22%) larger in our measurements leading to a larger derived mass.The rms line profile presented by Denney et al. (2010) shows a clear double-peaked shape without strong wings, while the rms line profile in our observations has strong wings and the same blue peak but without the corresponding red peak, potentially indicating a significant change in the properties of the reverberating gas.
The Hβ mass from our time delay and line width measurements is in good agreement with the mass from CARAMEL modeling of log M/M ⊙ = 7.09 +0.35  −0.34 or M BH = 1.23 +1.52  −0.67 × 10 7 M ⊙ , and in this case, the good agreement is likely because the inclination angles preferred by CARAMEL (θ i = 33.2+13.5 −9.1 ) are similar to the population average value implied by f ≈ 5. Interestingly, this inclination angle agrees with the results of a reanalysis of the geometry of the spatially resolved narrow line region (Falcone et al. 2023), while an earlier analysis suggested a much smaller value of θ i ≈ 15 • (Fischer et al. 2013).
Additionally, while the black hole mass from this reverberation program agrees with the dynamical masses, there are still good reasons to revisit the dynamical constraints.Gas dynamical modeling is known to be susceptible to biases when the gas exhibits noncircular motions, and this worry should be taken seriously in the presence of an AGN (e.g., Verdoes Kleijn et al. 2006;Jeter et al. 2019).While stellar dynamical modeling is usually held to be more accurate, the stellar dynamical modeling analysis of Davies et al. (2006) relied on fairly shallow but high spatial resolution observations of the nucleus of NGC 3227, and the S/N of the observations only allowed V and σ to be constrained for the stellar kinematics.Best practices in stellar dynamical modeling require that the kinematic constraints include not only V and σ, but also additional higher order terms to reduce the degeneracy between velocity anisotropy of the stellar distribution and M BH (Binney & Mamon 1982).Because of this potential degeneracy, deeper observations that provide more detailed information regarding the nuclear stellar kinematics are needed to accurately assess the dynamical constraints on the black hole mass in NGC 3227.

BLR Structure and Kinematics
Modeling of the velocity-resolved Hβ response in NGC 3227 with CARAMEL suggests that the BLR may be represented by the surface of a biconical structure or flared disk, oriented at ∼ 33 • to our line of sight and with gas motions that are dominated by rotation around a black hole with M BH ≈ 1.2 × 10 7 M ⊙ .While there are independent measurements of the inclination and black hole mass with which we can compare the CARAMEL constraints, as we have described above, there are few options for independent tests of the other BLR properties determined by CARAMEL.
One promising new tool that is currently being developed is BELMAC (Broad Emission-Line Mapping Code; Rosborough et al. 2023), an extension of TORMAC (Almeyda et al. 2017(Almeyda et al. , 2020) ) that has been adapted for application to the BLR.BELMAC simulates the velocityresolved response of an emission line to an input driving light curve.Using photoionization grids and a 3D ensemble of gas clouds, BELMAC simulates the observed response of the gas clouds to continuum variations for various user-specified geometries, velocity fields, and cloud properties.
While BELMAC is not yet able to independently explore and optimize the models to determine a best fit to the observed behavior of an emission line (though that capability is currently under development), we were able to use it in its current form to carry out two important tests: (1) when given the results of the CARAMEL models, could BELMAC confirm that the models provided a good fit to the data?, and (2) what types of geometric and kinematic parameters allow BELMAC to most closely replicate the observed integrated Hβ light curve and mean and rms line profiles?
We began by calculating the photoionization grids using CLOUDY v17.03 (Ferland et al. 2017) assuming a bolometric luminosity for NGC 3227 of log L bol = 43.16erg s −1 (where L bol = 9λL 5100 ; Kaspi et al. 2000;Merritt 2022) and an SED based on the average of 17 local, unobscured AGNs with Eddington ratios < 0.1 from Jin et al. (2012).The BELMAC model results shown in Figure 8 were then created using the median values of r min , r mean , θ i , and M BH listed in Table 4. Al-  4. though γ and ξ are not BELMAC parameters, they imply a biconical geometry, as shown in Figure 7, which is a BELMAC geometry configuration.To obtain a better match to the Hβ data, θ o was used as a guide rather than a fixed parameter.Assuming a bicone, θ o is approximately the separation angle between each cone and the midplane.Therefore, the angular distribution of clouds ranged from 5 − 55 • , where 0 • is the midplane, for each cone.Similarly, the emission distribution described by CARAMEL's β and σ r parameters provide estimates for the cloud distribution, gas density distribution, and covering fraction.In BELMAC, the number density of clouds with radius and gas density of the clouds with radius are, respectively, N (r) ∝ r p and n H (r) ∝ r −s , where p and s are the free parameters.The cross-sectional area of the clouds is primarily determined from the covering fraction.To approximately match CARAMEL's emission distribution, we adopted p = −2, s = −1.45,with a gas density at the outer radius of log n H = 9.3 cm −3 , and a covering fraction of 0.4.Lastly, rotational and radial motions were included in the BELMAC model, where the radial cloud motion is influenced by radiation pressure and gravity.Therefore, using the CARAMEL parameters listed in Table 4 to guide the above described BELMAC parameters, the radial motion is entirely an inflow.In general, the model is able to reproduce the integrated light curve fairly well, and does a reasonable job matching the mean and rms line profiles, though there is some mismatch in the wings of the mean line profile and the blue peak of the rms profile.In comparison, Figure 9 shows the best BELMAC model results to the data.Since there is currently no parameter optimization capability, the parameters were manually adjusted until a reasonable match to the data was found.However, the same r min , r mean , and M BH determined by CARAMEL were used, as well as the same photoionization grids.The best model found is a disk with a half-angular width 30 • (analogous to θ o ) and inclined to our line of sight at 35 • .About 53% of the clouds in the model have bound elliptical orbits and 47% are unbound hyperbolic orbits, of which 98% are inflowing.The geometry is described by p = −1.7,s = −1.3,where the gas density at the outer radius is log n H = 9.2 cm −3 , and the covering fraction is 0.33.Other than the overall shape of this BELMAC BLR model, which is a thick disk rather than a bicone, and a preference for a stronger inflow component, the parameters are in general agreement with those inferred by CARAMEL, lending additional credence to the derived inclination angle and black hole mass.Furthermore, while some of the details appear to be slightly different between the preferred models from the two codes, we note that it is not straightforward to directly compare them since CARAMEL models the emission itself while BELMAC includes a full distribution of gas clouds, including those that may be shielded or only weakly emitting.Nevertheless, both agree that the Hβemitting BLR is extended in radius and of similar predicted size, that the geometry is flattened rather than spherical, and that the kinematics are not dominated by outflows (as might be expected for a disk wind) but instead are dominated by rotation and/or inflow.

SUMMARY
We have carried out a new optical reverberation mapping campaign focused on the nearby broad-lined Seyfert NGC 3227.Using the LCO global network of telescopes, we employed broad-band photometry to track the continuum variations coupled with spectroscopy to track the responses of the broad emission lines.Time delays relative to the continuum variations were found for the Hβ, He II, Hγ, and Hδ broad emis-sion lines.We also detected velocity-resolved time delays across the broad Hβ emission-line profile.
Modeling of the velocity-resolved Hβ response with CARAMEL finds a black hole mass of log M/M ⊙ = 7.09 +0.35  −0.34 or M BH = 1.2 +1.5 −0.7 × 10 7 M ⊙ , which agrees well with the mass determined from the Hβ time delay and line width when atypical scaling factor is assumed, M BH = 1.1 +0.2 −0.3 × 10 7 M ⊙ .This good agreement is likely due to the intermediate inclination that is preferred by the models, θ i ≈ 33 • , which is similar to the populationaverage inclination angle that is implied by typical RM scale factors of f = 4 − 5.The models also suggest that the Hβ-emitting BLR may be represented by the inner surface of a biconical or flared disk structure, and that the gas motions are dominated by rotation.
Finally, using the photoionization-based modeling code BELMAC, we confirmed that the CARAMEL results provide a reasonable fit to the Hβ light curve and mean and rms line profiles, though there is some mismatch in the line wings.To correct for this mismatch, BELMAC prefers a thick disk geometry, rather than a bicone, and roughly equal weight between rotation and inflow.While they disagree about these details, both codes agree on several general properties, including that the Hβ-emitting BLR is extended in radius, flattened rather than spherical, and that the kinematics are dominated by rotation and/or inflow rather than outflow.Future improvements to the BELMAC code will allow a more thorough exploration of reasonable BLR parameter combinations for this AGN and others, and a more detailed comparison with the modeling results from CARAMEL.
MCB gratefully acknowledges support from the NSF through grant AST-2009230. MM
Figure 1.The inner 10 ′ × 10 ′ region of the reference image, centered on the nucleus of NGC 3227 and oriented with North up and East to the left.

Figure 2 .
Figure 2. The final calibrated and binned V −band light curve for NGC 3227.For comparison with the AGN variations and uncertainties, the bottom panel shows the scatter in residuals for a non-varying field star of similar brightness.
Figure 3.Mean (top) and rms (bottom) spectra of NGC 3227 based on the observed-frame spectra obtained throughout the monitoring program.The strong broad emission lines are labeled for ease of identification.

Figure 4 .
Figure 4. Light curves (left) for the continuum and broad emission lines throughout the course of the monitoring program.Cross-correlation functions (right) are for each light curve relative to the continuum, and in the case of the continuum is the auto-correlation function.

Figure 6 .
Figure 6.Histograms displaying the posterior distributions of the BLR model parameters for Hβ.

Figure 7 .
Figure 7. Representative geometric model of the Hβ emission in the BLR of NGC 3227.The left panel is oriented edge on, with an Earth-based observer at large +x values, and the right panel shows the Earth-based observer's view.The transparency of each point is related to the strength of the response of the gas, with more opaque points having a stronger response to continuum fluctuations.
Radius vs. Luminosity Previous measurements of the Hβ time delay in NGC 3227 range from ∼ 2 days (De Rosa et al. 2018) to ∼ 7 days

Figure 8 .
Figure8.Hβ and model light curves (left) as black points and blue curve, respectively.Mean and rms line profile in black and grey, respectively, and the model mean and rms line profiles in blue and light blue, respectively, for broad Hβ.The narrow Hβ line profile is in yellow and the overall fit is in red (right).The light curve and line profile models were created using BELMAC with parameters interpreted from Table4.

Figure 9 .
Figure 9.The same as Figure 8, but with model parameters adjusted to best match the Hβ light curve, line profile, and rms profile.

Table 1 .
Emission-Line Time Lags and Widths

Table 2 .
Note-Reported time lags are in the observer's frame, while, per convention, line widths are measured in the rest frame.Continuum and Emission-Line Light Curves

Table 3 .
Light Curve Statistics

Table 4 .
Broad-line region model parameter values was supported by National Science Foundation grant No. 2050829 to Georgia State University under the Research Experiences for Undergraduates program.SR was supported by NSF grant AST-2009508.CAO was supported by the Australian Research Council (ARC) through Discovery Project DP190100252.MV was supported by the NSF through AST-2009122.TT gratefully acknowledges supprt by NSF through grant AST-1907208.