GEOMETRIC AND DYNAMICAL MODELS OF REVERBERATION MAPPING DATA

The energy emitted by active galactic nuclei (AGNs) is argued to be the result of matter accreting onto supermassive black holes at the center of galaxies (Lynden-Bell & Rees 1971). However, the details of the geometry and kinematics of the region around the accretion disk are not well understood. In the standard model of AGNs, the region around the accretion disk is called the broad-line region (BLR) due to the broad emission lines from rapidly moving clouds of material near the black hole (Antonucci 1993; Urry & Padovani 1995). Models for the BLR attempt to explain the many categories of AGNs by the observer's viewing angle and the covering fraction of inflowing or outflowing material (see, e.g., Murray et al. 1995; Elvis 2000), since the degree to which the BLR geometry and kinematics vary between individual systems is also unknown.

The mass of the central black hole is a fundamental parameter in galaxy evolution, as suggested by a relation between black hole mass and stellar velocity dispersion of the host galaxy bulge, the M_BH–σ_⋆ relation (see, e.g., Bennert et al. 2011). While the black hole masses of very nearby galaxies can be measured by observing the orbits of stars or gas, a different approach is needed for more distant galaxies because the gravitational sphere of influence of the black hole cannot be spatially resolved (e.g., Ferrarese & Ford 2005). In active galaxies, the small size of the BLR, estimated to be around ∼10¹⁴–10¹⁶ m (Wandel et al. 1999; Kaspi et al. 2000; Bentz et al. 2006), inhibits direct imaging of the accretion disk and orbiting BLR clouds.

Reverberation mapping provides a method to determine the black hole mass, along with the geometry and the kinematics of the BLR (Blandford & McKee 1982; Peterson 1993; Peterson et al. 2004). Without relying on spatially resolving the gravitational sphere of influence of the black hole, reverberation mapping provides a powerful tool for studying black holes over a range of redshifts and black hole masses (e.g., Peterson et al. 2004; Woo et al. 2007; Bentz et al. 2009; Denney et al. 2009). The method relies on the large time variability of AGN luminosity, spanning timescales of days to years (e.g., Webb & Malkan 2000). Reverberation mapping data consist of a time series of frequent measurements of the intensity of the continuum and broad-line emission. The line emission strength is assumed to be proportional to the continuum emission strength, but with a time lag due to the light travel time from the central ionizing source to the BLR material. The time lag between line and continuum flux contains information about the size of the BLR, while the shape of the spectral line encodes the velocity information. An estimate of the black hole mass can be calculated assuming the BLR clouds orbit in the Keplerian potential of the black hole with velocities determined by the width of the spectral line and at a radius given by the average lag between the line and continuum fluxes.

A weakness of this traditional reverberation mapping method is that the relation between velocity and position observables of the clouds and the black hole mass depends on an unknown dimensionless proportionality constant that depends on the geometry and kinematics of the BLR. In practice, this so-called virial factor is estimated based upon some external criteria, such as the average factor that makes the M_BH–σ_⋆ relation consistent between different black hole mass estimators and between samples of active and inactive galaxies (Onken et al. 2004; Collin et al. 2006; Woo et al. 2010; Greene et al. 2010; Graham et al. 2011). Ideally, we would like to infer directly the morphology and kinematics of the BLR, and the black hole mass, including its uncertainty (for a discussion of potential systematic errors in reverberation mapping see Krolik 2001). The models for the structure and kinematics may include a net inflow or outflow of BLR clouds, among other physically motivated models (e.g., Murray et al. 1995; Marconi et al. 2008). Developing such a method is the goal of this paper.

The data required for reverberation mapping encode the geometry and kinematic information to some degree, depending upon the quality, in the form of the transfer function (or response function) that maps the continuum emission onto the line emission. The average lag used to estimate black hole mass is the first moment of the transfer function. Previous analysis involved estimating the transfer function and then interpreting the transfer function in relation to a model of the BLR (see Krolik & Done 1995; Done & Krolik 1996; Horne et al. 2003; Bentz et al. 2010). This is necessary because the transfer function is a function of time lag, not position within the BLR, so its interpretation requires a physical model. Estimating the transfer function requires inverting a linear integral equation, and while the method of Krolik & Done (1995) and Done & Krolik (1996) uses regularized linear inversion, thus allowing for uncertainty estimation, other inversion methods such as "maximum entropy" do not allow for straightforward uncertainty estimation or model selection (e.g., Horne et al. 2003; Bentz et al. 2010).

Our method of analyzing reverberation mapping data simplifies the process of obtaining a transfer function and then interpreting the result using different models. We compare reverberation mapping data directly with models of the BLR, obtaining uncertainty estimates as well as allowing for model selection. Once we have found models and model parameters that fit the data, we can easily compute the transfer function and average time lag. Our goal is to constrain the geometry and kinematics of the BLR and provide an internally consistent factor for the black hole mass. We note that the traditionally determined average time lag is exactly equivalent to a model where the BLR is a face-on ring of a given radius (response = δ-function) or a spherical shell (response = step function). This implicit assumption drives the inference on the average lag and its result, as we will show in this paper.

An important part of our method for directly modeling reverberation mapping data requires that we predict the AGN continuum light curve between the observations. Recent work has found that AGN continuum light curves are well-modeled by a damped random walk: Kelly et al. (2009) used a continuous time stochastic process; Kozłowski et al. (2010) used the formalism of Press et al. (1992) and Rybicki & Press (1992, 1994). This model for AGN variability was applied to ∼900 AGNs by MacLeod et al. (2010) in order to correlate variability with other parameters of AGNs. Zu et al. (2010) were then able to use this model for AGN variability to improve the standard analysis of reverberation mapping data, including a better understanding of the uncertainties involved. They model the continuum light curve using Gaussian Processes to recover the transfer function, assumed to be a top-hat. As with Kozłowski et al. (2010), they use an exponential covariance matrix to relate the continuum flux at different points in the time series. We also use Gaussian Processes to model the continuum light curve, as well as a slightly more general exponential covariance matrix. Our method improves upon the approach of Zu et al. (2010) by modeling the reverberation mapping data directly in terms of a geometric and dynamical model, rather than recovering the transfer function.

Bottorff et al. (1997) have also modeled reverberation mapping directly in an attempt to understand the BLR dynamics in the well-studied AGN NGC 5548. They expand upon the hydromagnetically driven outflow model of Emmering et al. (1992) and use one set of parameter values to compare their model with NGC 5548. While the specific models presented here are clearly not as sophisticated from a physical point of view, our method improves upon that approach by finding the best fit parameter values of our simple models and believable estimates of their uncertainties.

We consider two types of reverberation mapping data sets: velocity-unresolved, where there is a time series of the continuum flux and a time series of the integrated line flux, and velocity-resolved, where the data consist of a continuum flux time series, and a series of entire line spectra as a function of time.

The paper is organized as follows. In Section 2, we define and describe the physical problem. In Section 3, we outline our methods in the formalism of Bayesian probability theory and describe the algorithms we use to compare reverberation mapping data to mock data created from a model of the BLR. In Section 4, we test our method using simple models of the BLR and show that we are able to recover the parameter values of our test systems. Finally, in Section 5, we summarize our conclusions. Flux units throughout the paper are arbitrary, but computed consistently within our method.

Throughout this paper, we assume a simple model for the BLR, described as follows. The AGN is defined to be at (0, 0, 0), and the observer is at (d, 0, 0). We model the distribution of BLR gas by defining the gas density profile ρ(x, y, z), assumed to be normalized such that

$$\begin{equation} \int _V \rho (x, y, z) \, dV = 1, \end{equation} \tag{ 1 }$$

where dV = dx dy dz and V is all of space. We assume that the gas absorbs the ionizing radiation, but is not self-shielding, so that gas at larger radii is still illuminated. It should be noted that our approach is fully general and can support more complex models of the optical properties of the BLR, as well as its geometry and dynamics.

2.1. Velocity-unresolved Reverberation Mapping

If the continuum flux varies with time according to f_cont(t), then the total line flux as a function of time is given by

$$\begin{equation} f_{\rm line}(t) = A \int _V \ f_{\rm cont}(t - l(x, y, z)) \rho (x, y, z) \, dV, \end{equation} \tag{ 2 }$$

where l(x, y, z) is the lag, or time delay, associated with BLR gas at position (x, y, z), and A is a response coefficient. The lag l for each position is simply the excess light travel time from taking a path starting at (0, 0, 0) that travels to some gas at (x, y, z), where the light is absorbed and re-emitted as line emission, and that finally travels to the observer, relative to a direct path straight from the AGN to the observer:

$$\begin{eqnarray} l(x,y,z) &=& (\sqrt{x^2 + y^2 + z^2}+ \sqrt{(x - d)^2 + y^2 + z^2} - d)/c.\nonumber\\ \end{eqnarray} \tag{ 3 }$$

For any case of interest, $d \gg \sqrt{x^2 + y^2 + z^2}$, and therefore this is well approximated by

$$\begin{equation} l(x,y,z) \approx (\sqrt{x^2 + y^2 + z^2} - x)/c, \end{equation} \tag{ 4 }$$

which is the formula adopted throughout this paper. See Figure 1 for an illustration of this model.

Zoom In Zoom Out Reset image size

Note that Equation (2) is a special case of the general equation

$$\begin{equation} f_{\rm line}(t) = A \int \Psi (\tau) f_{\rm cont}(t - \tau) \, d\tau, \end{equation} \tag{ 5 }$$

where Ψ(t) is the so-called transfer function, which gives the response of the line flux to a delta-function pulse in the continuum flux.² Thus, for any particular system, if we can infer the density of BLR clouds throughout space, we can automatically deduce the corresponding transfer function:

$$\begin{equation} \Psi (\tau) = \int _V \delta {\left(\tau - l(x, y, z)\right)}\rho (x, y, z) \, dV. \end{equation} \tag{ 6 }$$

The meaning of this equation is that each location in space contributes to the transfer function at the value of the location's lag, with the size of the contribution being proportional to the amount of gas at that location.

2.2. Velocity-resolved Reverberation Mapping

Now suppose that the BLR gas is in motion, such that the system can be described by a time-invariant distribution function g defined over the phase space of a single particle:

$$\begin{equation} g(x, y, z, v_x, v_y, v_z) = \rho (x, y, z)g(v_x, v_y, v_z | x, y, z). \end{equation} \tag{ 7 }$$

The motion of the gas along the line of sight is assumed to affect the wavelength of re-emitted light, but its distribution function is assumed to be time invariant and therefore does not vary during the observing campaign. Then the emission-line profile at time t will be a function of the line-of-sight velocity, v_los:

$$\begin{eqnarray} f_{\rm line}(v_{\rm los},t) &=& A \int _{v_y, v_z} \int _V \ f_{\rm cont}(t - l(x, y, z)) \nonumber\\ &&\times g(x, y, z, v_x, v_y, v_z) \, dx \, dy \, dz \, dv_y \, dv_z,\qquad \end{eqnarray} \tag{ 8 }$$

where v_los is in the x-direction. This is the velocity-resolved equivalent of Equation (2).

Our method for constraining the geometry and kinematics of the BLR is an application of Bayesian Inference (Sivia & Skilling 2006). In general, to infer parameters θ from data D, we begin by assigning a prior probability distribution p(θ) describing our initial uncertainty about the parameters. Sampling distributions p(D|θ) are also assigned to describe our uncertainty about how the data are related to the parameters. Once specific data D = D* are obtained, our updated state of knowledge about the parameters is described by the posterior distribution, given by Bayes' rule:

$$\begin{equation} p(\theta | D=D^*, I) \propto p(\theta |I) p(D|\theta, I)|_{D=D^*}. \end{equation} \tag{ 9 }$$

Here, I is any background information we have about the problem. In complex problems, where θ consists of a large number of parameters, Monte Carlo methods are used to produce random samples from the posterior distribution for θ. Methods such as Nested Sampling (Brewer et al. 2010) can also provide the normalization constant for the posterior, known as the evidence, which is the key quantity for comparing the entire model with an alternative (Sivia & Skilling 2006).

In our method, the parameters θ to be inferred are those describing the spatial profile of the BLR gas, and the continuous continuum flux, f_cont(t). Since it is impossible to represent a continuum in a computer, we instead infer f_cont(t) evaluated at 500 time points, covering a time interval larger than the continuum data. The continuum modeling technique is described in detail in the following section.

Throughout this paper, both the continuum flux and line flux time series are considered part of the data set D:

$$\begin{equation} D = \left\lbrace \mathbf {y}_{\rm line}, \mathbf {y}_{\rm continuum} \right\rbrace. \end{equation} \tag{ 10 }$$

The prior information consists of the times at which the line flux and continuum flux are measured, t, and the error bars on the line flux and continuum flux measurements, σ:

$$\begin{equation} I = \left\lbrace \left(\mathbf {t}, {\bm \sigma }\right)_{\rm line}, \left(\mathbf {t}, {\bm \sigma }\right)_{\rm continuum}\right\rbrace. \end{equation} \tag{ 11 }$$

The likelihood function is chosen to be Gaussian, centered around the model-predicted line flux time series:

$$\begin{equation} p(D|\theta) = \prod _{i=1}^n \frac{\exp \big[-\frac{1}{2}\left(\frac{y_{i, \rm line} - m_i(\theta)}{(\kappa \sigma _i)}\right)^2\big]}{(\kappa \sigma _i)\sqrt{2\pi }}, \end{equation} \tag{ 12 }$$

where κ is a "noise boost" parameter to account for the presence of unknown systematic effects not included in the reported error bars, such as those due to flux calibration, wavelength calibration, and continuum subtraction.

Once the posterior distribution is obtained, many different algorithms are available for exploring it and computing summaries such as marginal distributions for parameters. We have implemented our model with two methods, the first is Metropolis–Hastings, a Markov Chain Monte Carlo (MCMC) algorithm, which provides samples from the posterior probability distribution function (PDF) for the model parameters. The second is Diffusive Nested Sampling (Brewer et al. 2010), which provides samples from the posterior PDF and an estimate of the evidence value for the model. Although the evidence calculation makes the second algorithm significantly slower than the first, Diffusive Nested Sampling is much faster than alternative MCMC-based implementations of Nested Sampling (Brewer et al. 2010). The results presented here to test the method use the MCMC algorithm, while the Diffusive Nested Sampling algorithm is used to apply the method to real reverberation mapping data (B. J. Brewer et al. 2011, in preparation).

3.1. Continuum Interpolation

In order to create a mock line flux time series to compare with the data, it is necessary to interpolate between the continuum flux data points. Linear interpolation is the simplest approach, but it does not provide an estimate of the uncertainty in the interpolation, suggesting that we know precisely the value of the continuum f(t) at all times between the measured data points. If we want to obtain reliable uncertainties in our results, we should acknowledge the uncertainty introduced by the interpolation process.

To account for this, we consider the entire continuum function f_cont(t) to be an unknown parameter to be inferred from the data. The prior distribution for f_cont(t) is a Gaussian Process (MacKay 2003; Rasmussen & Williams 2006), which is a convenient class of probability distributions over function space. Given a mean function μ(t) and a covariance function C(t₁, t₂), the probability distribution for the function value f at any finite set of times is a multivariate Gaussian:

$$\begin{equation} p(\mathbf {f}|\mu, C) = \frac{1}{\sqrt{(2\pi)^n {\rm det\, }\mathbf {C}}} \exp \left(-\frac{1}{2}(\mathbf {f} - {\bm\mu })^T\mathbf {C}^{-1}(\mathbf {f} - {\bm\mu })\right), \end{equation} \tag{ 13 }$$

where μ is a vector of means at the relevant time points and C is the covariance matrix, obtained by evaluating the covariance function at the relevant times. In the reverberation mapping problem, f_cont(t) is constrained by two data sets: the continuum measurements and the line measurements. We parameterize the covariance function and mean function with four hyperparameters: μ (the long-term mean), σ (the long-term standard deviation), τ (typical timescale of variations), and α (a smoothness parameter between 1 and 2), such that the mean function is a constant μ(t) = μ and the covariance function is

$$\begin{equation} C(t_1, t_2) = \sigma ^2 \exp \left[-\left(\frac{|t_2 - t_1|}{\tau }\right)^\alpha \right]. \end{equation} \tag{ 14 }$$

The posterior distribution function for f(t) given some continuum data (but not the line data) is shown in Figure 2. Note that outside the areas where we have data, the uncertainty gets large, but in areas where the data are well sampled, the uncertainty in the interpolation is small. We keep track of f(t) at 500 times, both slightly preceding and following the data. Further interpolation between these 500 points is linear. Five hundred continuum parameters is sufficient to render the distance between continuum flux points much smaller than the maximum monitoring cadence, allowing us to resort to linear interpolation only on scales not probed by the data. We change the 500 parameters in the same way as the model parameters, with every new proposal for the continuum function related to the one before. The function f(t) can be parameterized by 500 variables with standard normal priors, which are converted to f(t) values by multiplication with the Cholesky decomposition of C. We note that our Gaussian Process method for interpolation, in the special case α = 1, is equivalent to the method of Zu et al. (2010), apart from computational details. α = 1 has also been used in detailed studies of quasar variability (e.g., MacLeod et al. 2010).

Zoom In Zoom Out Reset image size

3.2. Creating Mock and Simulated Data

In order to illustrate our method, we have developed simple models of the BLR region geometry and dynamics. As this method is fully general, it is also possible to implement more complex models within the framework described so far. We showcase these simple models by creating simulated data with known true parameter values in our models. This allows us to test our code as well as to explore the accuracy and precision of the results obtainable by this method of reverberation mapping analysis. Such tests on simulated data also allow us to ascertain the data quality needed to perform inferences regarding increasingly complicated model parameters. We showcase both geometry-only and geometry plus kinematics models, where the latter are the same as the first with the addition of velocity information given to the BLR gas. We show transfer functions for the geometry models and velocity-resolved transfer functions for the kinematics models.

To ensure that in our method the true parameter values are recovered, we save instances of each model and use them as simulated reverberation mapping data, adding noise and varying the time series characteristics to match reverberation mapping campaigns of varying quality. A simulated data set consists of line flux and continuum flux measurements. Given a BLR model, the continuous continuum light curve is all that we need to create, since the mock line flux measurements can be obtained from the model and continuous continuum light curve. We create continuous continuum light curves by using the hyperparameters of the Gaussian Processes continuum interpolation. The hyperparameters contain information about the timescales and levels of variability in an AGN continuum time series. We use values for the hyperparameters from interpolation of the Lick AGN Monitoring Project (LAMP; Walsh et al. 2009; Bentz et al. 2009) continuum time series of Arp 151, one of the most variable AGNs in the LAMP sample. The values used for the hyperparameters were μ = 75 (arbitrary units), σ = 30 (same units as μ), τ = 6 × 10⁶ s, and α = 1.5 (dimensionless).

4.1. Geometry Model: Ring/Disk/Shell

4.1.1. Model Definition

We use a flexible geometry model of the BLR gas density to test our method when only integrated line flux measurements are used instead of the full spectral shape. The model is that of a spherical shell centered on the central engine with parameters allowing partial, axisymmetric illumination of the shell and varying inclination of the resulting ring/disk. Examples of possible configurations, ranging from a complete shell to a thin ring/disk, are shown in Figure 3. The parameters of the model are the mean radius of the disk, r₀, the thickness of the disk in the radial direction, σ_r, the illumination angle of the shell, and the inclination of the shell. The illumination angle is defined so that values approaching 0 define an increasingly thin ring/disk and a value of π/2 defines a spherical shell. The inclination angle is defined so that values approaching 0 define a face-on ring/disk and a value of π/2 is an edge-on ring/disk. We use a normal distribution to define the radial thickness of the shell, so that r₀ and σ_r are the average and 1σ width of a normal distribution. The normal distribution is created in the x, y, and z Cartesian coordinates.

Zoom In Zoom Out Reset image size

It is important to set appropriate prior probability distributions for each model parameter. For parameters where we know the order of magnitude of the parameter value we use a flat prior in the parameter. Examples of parameters with flat priors in the parameter include the inclination angle and the illumination angle, which may only vary between 0 and π/2. For parameters where we do not know the order of magnitude of the parameter value we need a prior that treats many orders of magnitude equally, so we use a flat prior in the log of the parameter. Examples of parameters with flat priors in the log of the parameter include r₀ and σ_r. These choices of prior probability express complete ignorance in the value of a parameter within some reasonable range, but it is necessary to make the distinction between whether or not the order of magnitude of a parameter value is known. In the cases considered in the remainder of this paper, the posterior is much narrower than the prior, and therefore the inference is dominated by the likelihood, i.e., the information contained in the data.

The underlying spherical symmetry of these models and the angular dependence of the ring/disk model allow us to use spherical coordinates. In order to sample the gas phase-space density at a finite number of points, we use a grid in log(r), ϕ, and cos(θ). Using equal steps in cos(θ) instead of θ means that the volume of each grid point depends only on the radius, r. The density is then multiplied by the volume of the grid point to find the total mass of gas in each grid point. The emissivity of each grid point also depends on the radius r because the continuum ionizing radiation flux falls off as r⁻², requiring more gas mass at larger radii to have the same line flux contribution as less gas mass at smaller radii. In general, the illumination parameter allows us to model any axisymmetric ionizing flux.

We test our method to recover the BLR model parameters by creating simulated data, the true parameter values of which are given in Table 1. The continuous continuum function is obtained using the hyperparameters from the Gaussian Processes interpolation of Arp 151 reverberation mapping data, as described in Section 3.1, and evaluated at 120 consecutive "observations" one day apart. The line flux time series for each model are generated using this continuous continuum function and a given set of model parameters. The line flux time series contain 60 "observations" one day apart, starting 60 days after the start of the continuum flux "observations." These simulated data are meant to represent excellent reverberation mapping data, with an observation campaign of similar length to recent campaigns (see, e.g., Bentz et al. 2009), but without gaps due to difficult weather conditions. Additional noise has also been added to the simulated data. Most simulated data sets have line flux errors of 1.5%, which represents very favorable observing conditions, but we have also tested simulated data with errors of 5% to reflect the current typical error of reverberation mapping line flux measurements.

Table 1. Simulated Geometry Data True Parameter Values

Data	Model	Uncertainty (%)a	ro	σr	Inclination Angle	Illumination Angle
			(1014 m)	(1014 m)	(rad)	(rad)
1	Inclined disk	1.5	5	0.3 rmean = 1.5	0.79	0.22
2	Inclined disk	5	5	1.5	0.79	0.22
3	Edge-on disk	1.5	5	1.5	π/2	0.22
4	Face-on disk	1.5	5	1.5	0.0	0.22
5	Shell	1.5	5	1.5	...	π/2

Notes. Each simulated data set consists of 60 line emission data points and the same 120 continuum emission data points, where the line emission time series start half-way through the continuum time series. ^aLine flux uncertainty of the simulated data set.

Download table as:  ASCIITypeset image

4.1.2. Testing the Geometry Model

The first test is whether we can recover the parameter values of the simulated data using the MCMC algorithm described in Section 3. Since our one flexible geometry model encompasses a number of different geometries, such as a shell, thin or thick ring or disk, we do not have to consider model selection at this point. We test the many possible geometries of this model by creating five simulated data sets, whose true parameter values are given in Table 1. The simulated data sets include an inclined disk with line flux errors of 1.5% and 5% and an edge-on disk, a face-on disk, and a shell with line flux errors of 1.5%. The MCMC algorithm is typically run for 150,000 iterations and all parameter values are recovered to within two standard deviations of the posterior probability distributions of the parameters, with 10/13 recovered to within one standard deviation. This is as expected, since we should find the true parameter value to lie within 1σ about ∼68% of the time and to lie within 2σ about ∼95% of the time. The mean and standard deviation of the posterior distributions are given in Table 2, with the exception of many of the angular parameters, where the quoted 1σ uncertainty does not adequately describe the posterior distribution. Part of the reason for the standard deviation of the angular parameters not describing the posterior is due to the uneven step size in θ, so that values of the illumination angle close to π/2 and values of the inclination angle close to 0 rad have much poorer angular resolution. This might lead to an angular parameter being quoted as having a mean of 1.22 rad and a 1σ uncertainty of 0.29, as for the illumination angle of the shell model simulated data, but while this uncertainty may seem large, it corresponds to an uncertainty of only 1–2 grid points in θ. The posterior distributions for the face-on disk and shell simulated data are shown in Figures 4 and 5. Select joint probability distributions between the inclination and illumination angles are also shown in Figure 6 in order to show the degeneracies between different models. In particular, for the shell model, the inclination is not constrained unless the illumination angle is small, or rather, unless the sphere of BLR gas is not entirely illuminated.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Simulated Geometry Data Recovered Parameter Values

Data	Model	ro	σr	Inclination Angle	Illumination Angle
		(1014 m)	(1014 m)	(rad)	(rad)
1	Inclined disk	4.53 ± 0.47	2.93 ± 0.78	...	...
2	Inclined disk	4.81 ± 1.10	2.24 ± 1.56	...	...
3	Edge-on disk	4.76 ± 0.65	1.83 ± 0.60	...	...
4	Face-on disk	4.93 ± 0.67	3.10 ± 1.31	0.23 ± 0.13	0.29 ± 0.10
5	Shell	4.48 ± 0.52	1.85 ± 0.45	...	1.22 ± 0.29

Notes. Results for 150,000 iterations using an MCMC algorithm. See Section 4 for a discussion of why average values for the angular parameters are not quoted for most simulated geometry data sets.

Download table as:  ASCIITypeset image

The posterior PDFs for the five simulated data sets show that the edge-on disk, face-on disk, and shell geometries allow for excellent recovery of the parameter values with estimates of the uncertainty. For the two inclined disk simulated data sets, there is some degeneracy in the angular parameters, leading to large uncertainties in their average values. The MCMC algorithm finds a more likely geometry configuration than the true configuration for the inclined disk data sets, although the true configuration is still a valid possibility with posterior local maxima at the true parameter values. With the increased simulated line flux error from 1.5% to 5%, however, it becomes increasingly difficult to recover the angular parameters, and only the mean radius is recovered with a small enough uncertainty as to be useful in describing the BLR. This emphasizes the importance of obtaining high quality line flux data in reverberation mapping campaigns.

The time series and transfer functions for the face-on disk and shell MCMC geometry model tests are shown in Figures 7 and 8, respectively. The time series show the simulated data overlaid with mock data created with parameters sampled randomly from the posterior probability distributions. The fit of the mock data to the simulated data is excellent for all five models. The variety in the shape of the simulated data time series, all well-fit by their respective models, shows that the MCMC algorithm for model parameter value recovery is robust for a wide range of models. The transfer functions also show a variety of shapes. For a thin shell geometry, thinner than the shell of simulated data set 5, our resulting transfer function agrees with the analytic form of a top-hat function (see Peterson 1993).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Dynamical Model

4.2.1. Model Definition

In order to constrain the kinematics of the BLR and the mass of the central black hole, we must model the velocity distribution of the BLR gas in the context of a dynamical model. For simplicity of illustration and speed of computation, we consider here a cylindrically symmetric model where the BLR gas is considered to be made of test particles in bound orbits within the spherical Keplerian potential of the black hole. We parameterize the model in terms of energy and angular momentum, constants of the BLR gas motion, so we are guaranteed velocity and geometry distributions that do not evolve in time, and are therefore stationary during the monitoring campaign. In future papers, we will generalize the model to include unbound orbits to describe inflows and outflows, and also other physical mechanisms, such as radiation pressure or winds (Marconi et al. 2008; Netzer & Marziani 2010).

The model is illustrated in Figure 9. For any choice of angular momentum L, energy E, and black hole mass M_BH, the motion of the BLR test particles is then described by the standard conservation equation resulting in elliptical orbits in the plane perpendicular to the angular momentum. Given our cylindrical symmetry, we will consider families of angular momenta obtained by rotation around the z-axis and defined by the polar angle θ₀ (see Figure 9). The spatial density of the BLR is then given by

$$\begin{equation} P(r, \theta, \phi | E, L, \theta _0) \propto \frac{1}{v} \times \frac{1}{| \sqrt{\sin ^2\theta _0 - \cos ^2\theta }|}, \end{equation} \tag{ 15 }$$

where the angular term comes from integrating over the uniform distribution of azimuthal angle ϕ₀ of the angular momentum vector, and v is the total magnitude of the velocity vector:

$$\begin{equation} v = \sqrt{2E + \frac{2GM_{{\rm BH}}}{r}}. \end{equation} \tag{ 16 }$$

Owing to the symmetry of our model, we can consider only θ₀ < π/2 (i.e., L_z>0), obtaining the following limits on the allowed θ coordinate for the BLR:

$$\begin{eqnarray} &&\frac{\pi }{2}-\theta _0 < \theta < \frac{\pi }{2}+\theta _0. \end{eqnarray} \tag{ 17 }$$

As θ₀ approaches zero, the model represents a thin disk, while as θ₀ approaches π/2 the model covers the whole sphere. Conservation of energy and angular momentum limits the radial coordinate to the range:

$$\begin{equation} r>-\frac{GM_{{\rm BH}}}{2E} - \frac{1}{2} \sqrt{\left(\frac{GM_{{\rm BH}}}{E}\right)^2 + \frac{2L^2}{E}}, \end{equation} \tag{ 18 }$$

$$\begin{equation} r< -\frac{GM_{{\rm BH}}}{2E} + \frac{1}{2} \sqrt{\left(\frac{GM_{{\rm BH}}}{E}\right)^2 + \frac{2L^2}{E}}. \end{equation} \tag{ 19 }$$

Finally, E and L are connected by the usual condition

$$\begin{equation} |L| \le \frac{GM_{{\rm BH}}}{\sqrt{-2E}}. \end{equation} \tag{ 20 }$$

For every allowed value of r, θ, and ϕ, the component of the velocity vector along the line of sight can be computed in the standard manner, resulting in two solutions per position, in general (outbound and inbound; four if one considers L_z < 0 as well). More complex geometries and kinematics can be obtained by superpositions of multiple sets of E, L, and θ₀ values within the same potential given by M_BH. However, this further increases the dimensionality of parameter space and computing time. Therefore, in this example we will only use one such set.

Zoom In Zoom Out Reset image size

We apply prior probability distributions to the model parameters as described for the geometry model. The priors for the extra parameters in the dynamical model not part of the geometry model are as follows. The parameter θ₀ has a flat prior in the parameter ranging from 0 to π/2. The parameters M_BH, E, and L have flat priors in the log of the parameter.

In addition, in order to impose a BLR gas geometry, we model the distribution of illuminated gas, as the product of the spatial distribution given by the dynamical model with that imposed by one of our geometrical models, representing in this case the illumination function. This results in a broad range of geometries, giving the model a considerable flexibility (for example, in the future one could think of an anisotropic illumination function to model dust obscuration). The procedure is illustrated in Figure 10. Note that the radial distribution of the illuminated gas is not Gaussian anymore, as in the ring/disk/shell geometry model. The mean radius then is not the r₀ parameter of the geometry model, but must be computed numerically for each set of geometric and dynamical parameters. Similarly, the mean width is no longer σ_r and must be computed numerically.

Zoom In Zoom Out Reset image size

A model spectrum at a given time is obtained by summing all the line of sight velocities, weighted by the spatial density of illuminated gas multiplied by the continuum flux at an epoch corresponding to the appropriate time lag. In order to compare with real data, the model spectrum is then convolved with a Gaussian to represent instrumental broadening. Since we do not expect real data to match our model perfectly, we introduce a relatively large uncertainty in the form of the spectral line by adding Gaussian noise with a variance of σ²(F) = α F + β, where α = 0.00018 and β = 0.025. This model for the variance assumes both a dependence on spectral line flux F through the α parameter and a dependence on the continuum uncertainty through the β parameter. The units of α are flux and the units of β are flux². The specific values of α and β are related to the arbitrary flux units of our simulated spectra and result in a signal to noise of ∼4. Conservatively, this signal-to-noise ratio is lower than typically achieved in state-of-the-art spectral monitoring campaigns. Examples of synthetic spectra at a resolution of FWHM =13.1 Å, or ∼800 km s⁻¹ at the wavelength of Hβ, are shown in Figure 11. The face-on disk systems (top and middle panels of Figure 11) have velocity bins of ∼120 km s⁻¹, while the sphere system (bottom panel) has velocity bins of ∼20 km s⁻¹. Notice how the line shapes are clearly different even for models with the same black hole mass. This is a clear illustration of the power of velocity resolved reverberation mapping as a diagnostic of the BLR geometry as well as kinematics.

Zoom In Zoom Out Reset image size

4.2.2. Testing the Dynamical Model

We test our dynamical model by creating simulated data sets consisting of time series of the continuum flux and of the line profiles of a broad line. The line profiles of the simulated data sets are shown in Figure 11. The kinematic parameters E and L are initially chosen to satisfy nearly circular orbits of the BLR gas at the mean radius given by the illumination function. A disk of broad line emitting material can be constrained by either the illumination function or the value of θ₀.

The first simulated data set is a thin disk viewed nearly face-on, with dynamics imposed by a single value of energy and angular momentum. The thin disk is constrained by the value of θ₀, while the illumination function describes the whole sphere being illuminated. This means that all allowed orbits lie in the disk and that the rest of the sphere does not contain broad line emitting gas. The second simulated data set is also a thin disk viewed nearly face-on with a single value of energy and angular momentum, but for this case the illumination function constrains the disk. We choose a value of θ₀ close to π/2 so that orbits are allowed in the entire sphere. The third simulated data set is a fully illuminated sphere with orbits that are also allowed in the entire sphere, so again θ₀ is close to π/2. This is still an axisymmetric configuration, as the BLR gas density imposed by the kinematics depends upon the θ coordinate. The true parameter values of the three simulated data sets used to test the kinematics model are shown in Table 3.

Table 3. Simulated Dynamics Data True Parameter Values

Data	Model	〈S/N〉a	MBH	Mean Radius	Mean Width	Inclination Angle	Illumination Angle	θo
			(107 M☉)	(1014 m)	(1013 m)	(rad)	(rad)	(radians)
1	Face-on disk	4.6	1	1.132	4.484	0.1	π/2	0.3
2	Face-on disk	4.4	1	1.195	4.022	0.1	0.3	π/2
3	Sphere	3.5	1	1.195	4.022	0.1	π/2	π/2

Notes. ^aAverage signal to noise of the line flux profile. Each simulated data set consists of 60 line emission profiles and the same 120 continuum emission data points, where the line emission time series start half-way through the continuum time series. The simulated line emission profiles are created by taking the true model and adding Gaussian noise with a variance of v = α × Flux + β, where α = 0.00018 and β = 0.025.

Download table as:  ASCIITypeset image

We test each of the three simulated data sets assuming only one set of kinematic parameters. The parameter values inferred using our method are shown in Table 4, while the full posterior PDFs are shown for all parameters of interest for the first simulated data set in Figures 12 and 13. The posterior PDFs for the black hole mass, average radius of BLR gas mass, and average width of the BLR gas mass are also shown for the second and third simulated data sets in Figures 14 and 15. They show that the black hole mass, average radius, and average width of the BLR are well determined for all three simulated data sets. The angular parameters are also well determined when physically possible. For example, for the first dynamics simulated data set of a face-on disk with orbits confined to the disk, the inclination angle and θ₀ are determined to within one or two grid points, while the illumination angle is only constrained to be ≳0.3 rad. The illumination angle cannot be determined more accurately because the BLR gas emission only comes from the disk, so as long as the entire disk is illuminated the spectrum is not sensitive to further changes in the illumination angle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Simulated Dynamics Data Recovered Parameter Values

Data	Model	MBH	Mean Radius	Mean Width	Inclination Angle	Illumination Angle	θo
		(107 M☉)	(1014 m)	(1013 m)	(rad)	(rad)	(rad)
1	Face-on disk	0.95 ± 0.05	1.04 ± 0.04	4.27 ± 0.05	0.11 ± 0.07	...	0.31 ± 0.02
2	Face-on disk	1.10 ± 0.13	1.12 ± 0.04	3.80 ± 0.03	0.12 ± 0.06	0.31 ± 0.02	1.40 ± 0.13
3	Sphere	1.00 ± 0.04	1.21 ± 0.08	3.95 ± 0.04	0.12 ± 0.07	1.46 ± 0.06	1.50 ± 0.05

Notes. Results for 470 × 10³ (data 1), 330 × 10³ (data 2), and 110 × 10³ (data 3) iterations using an MCMC algorithm. See Section 4 for a discussion of why average values for the angular parameters are not quoted for the illumination angle of data 1.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Finally, we also compute the velocity-resolved transfer functions for the three simulated data sets, shown in Figure 16. As expected, the transfer functions for the face-on disk configurations show little response at very small lags, while the sphere configuration shows the highest intensity of response at small lags. The transfer functions for the face-on disk configurations are similar, but clearly lead to different line profiles, again illustrating the power of modeling the full data set rather than just trying to model the transfer function.

Zoom In Zoom Out Reset image size

We introduce and test a new method for analyzing reverberation mapping data of AGNs by directly modeling the BLR. We illustrate our method by creating simple geometry and dynamical models of the BLR. Using a model of the BLR geometry to reproduce the integrated line flux time series from reverberation mapping data allows us to estimate the average radius of the BLR, as well as the mean width, illumination function, and inclination angle to the line of sight. Models of the BLR that include geometry and dynamical information allow us to additionally estimate the black hole mass and obtain an estimate of the extent to which the BLR gas orbits are confined to a disk or the whole sphere.

Our method of analysis provides several advantages over previous methods. First, previous methods rely upon cross-correlation to obtain a mean radius for the BLR and a virial relation with unknown virial coefficient to obtain an estimate of the black hole mass. Our method estimates the black hole mass self-consistently, without the need for a virial coefficient. Second, work modeling reverberation mapping data has previously focused on modeling the velocity-resolved or unresolved transfer function. However, the implications for the geometry and kinematics of the BLR are not clear for such analysis, as the transfer function is a function of the lag between the continuum and line emission. Instead of modeling the transfer function and then interpreting the transfer function in terms of a geometrical or dynamical model of the BLR, we focus on modeling the BLR directly. This allows us to extract more information and thus constrain the models more tightly. Finally, our fast method provides estimates of the uncertainty in the model parameter values and can be used with numerical algorithms such as Nested Sampling that allow for model selection. Our main results can be summarized as follows.

• 1.  
  We create simulated data sets using the geometry model with known true parameter values and find that we can recover these values with uncertainties that depend upon the random uncertainty of the reverberation mapping data. We can recover the mean radius of the BLR to within ∼0.1 dex and the mean width of the BLR to within ∼0.2 dex for simulated data with an integrated line flux uncertainty of 1.5%. We can also place constraints on the inclination and illumination with uncertainties of ∼0.2  rad for simulated data with face-on and spherical geometry configurations and 1.5% integrated line flux uncertainty. Current integrated line flux uncertainties of about ∼5% are on the edge of what would allow for successful recovery of more than just a mean radius for the BLR.
• 2.  
  We create simulated data sets using the dynamical model that consist of time series of a broad-line profile and we compare them to mock spectra made using our model. Despite the larger number of free parameters in our dynamical model, we find that we can recover all the parameters physically possible because the line profile is a stronger constraint on the model than the integrated line flux. We can recover the black hole mass and the mean radius of the BLR to within ∼0.05 dex, for simulated data with a line profile signal-to-noise ratio of ∼4 per spectral pixel. We can also recover the mean width of the BLR to within ∼0.1 dex and the inclination angle and illumination angle to within ∼2 grid spacings over which the BLR density is defined.

The small random uncertainties obtained in our tests of the simple geometry and dynamical models are partly due to the inherent assumption that our simulated data are drawn directly from the set of possible model configurations. In order to simulate the expected systematic error in applying simple models to complicated real BLR systems, we have added substantial Gaussian noise to instances of the model in order to create our simulated data sets. The time series of line profiles, in the case of the dynamical model, is very constraining, and leads to the reduced random uncertainty in the mean radius and mean width of the BLR by a factor of two for the dynamical model, as compared to the geometry model. When applying the method to real data we expect larger uncertainties, owing to modeling errors. The uncertainties quoted here should therefore be considered as lower limits to the overall precision of the method for data of comparable quality. This emphasizes the importance of good quality data and increasingly more realistic models, for recovering detailed information about the BLR from reverberation mapping data.

While we have created and tested both simple geometry and dynamical models, our method is more general, allowing for use of any geometry or dynamical model that can be simply parameterized. We plan to expand our library of models to include inflowing or outflowing BLR gas, which may be needed to explain some of the line profile asymmetries of current reverberation mapping data.

We thank the referee for helpful comments. We thank our friends and collaborators in the LAMP 2008 project for many insightful conversations. We are grateful to Chris Kochanek and Vardha Nicola Bennert for helpful suggestions on the manuscript. We acknowledge support by the NSF through CAREER award NSF-0642621, and by the Packard Foundation through a Packard Fellowship. A.P. also acknowledges support by the NSF through the Graduate Research Fellowship Program.