Modeling the Reverberation Response of the Broad-line Region in Active Galactic Nuclei

The variable continuum emission of an active galactic nucleus (AGN) produces corresponding responses in the broad emission lines, which are modulated by light travel delays, and contain information on the physical properties, structure, and kinematics of the emitting gas region. The reverberation mapping technique, a time series analysis of the driving light curve and response, can recover some of this information, including the size and velocity field of the broad-line region (BLR). Here we introduce a new forward-modeling tool, the Broad Emission Line MApping Code, which simulates the velocity-resolved reverberation response of the BLR to any given input light curve by setting up a 3D ensemble of gas clouds for various specified geometries, velocity fields, and cloud properties. In this work, we present numerical approximations to the transfer function by simulating the velocity-resolved responses to a single continuum pulse for sets of models representing a spherical BLR with a radiatively driven outflow and a disklike BLR with Keplerian rotation. We explore how the structure, velocity field, and other BLR properties affect the transfer function. We calculate the response-weighted time delay (reverberation “lag”), which is considered to be a proxy for the luminosity-weighted radius of the BLR. We investigate the effects of anisotropic cloud emission and matter-bounded (completely ionized) clouds and find the response-weighted delay is only equivalent to the luminosity-weighted radius when clouds emit isotropically and are radiation-bounded (partially ionized). Otherwise, the luminosity-weighted radius can be overestimated by up to a factor of 2.


INTRODUCTION
Rapidly growing supermassive black holes (SMBHs) are fed by a surrounding accretion disk, which together compose the central engines of active galactic nuclei (AGN).Most, if not all, large galaxies host a SMBH and its mass, M • , is closely tied to the galaxy's evolution.Decades of M • estimates show a correlation to the host galaxy's stellar velocity dispersion (σ ⋆ ), the M • − σ ⋆ relation (McConnell & Ma 2013), implying co-evolution between the SMBH and the bulge.The power of an AGN is regulated by the SMBH accretion rate, which is limited by its Eddington luminosity, where L Edd ∝ M • and M • can be 10 6 − 10 10 solar masses ( M ⊙ ).Outflows and jets from AGN are thought to have a profound influence on galaxy evolution by suppressing star formation.Therefore, establishing how the SMBH mass function evolves with redshift and luminosity is critical to understanding galaxy formation and evolution (Di Matteo et al. 2005;Hopkins et al. 2008;Fabian 2012;Alexander & Hickox 2012).

The Broad Line Region
The AGN's central engine is surrounded by a geometrically and optically thick structure comprised of molecular gas and dust clouds, generally known as the dusty "torus" (e.g., Urry & Padovani 1995).Depending on grain composition and size, dust sublimates at temperatures between ∼ 1500 − 1800 K and cannot not survive the heating by the accretion disk's UV/optical continuum within the corresponding dust sublimation radius, R d ∼ 1 pc (assuming an AGN luminosity ∼ 10 45 erg s −1 ).Inside R d is the broad line region (BLR), a largely dust-free zone where dense gas clouds are photoionized by the continuum radiation (e.g., Baskin et al. 2014;Baskin & Laor 2018;Amorim et al. 2021a).Residing well within the SMBH's gravitational sphere of influence, the gas clouds emit spectral lines that are Doppler-broadened by a several-tens of thousand km s −1 (e.g., Netzer 2008).The BLR is therefore a unique diagnostic of the SMBH's mass and the physical processes that operate in its environment.
The mass of the black hole, M • , can be estimated from the velocity dispersion (∆v) and radius (R BLR ) of the BLR, where G is the gravitational constant (e.g., Bentz & Katz 2015).However, since AGN are too distant and the BLR is far too small in angular size to obtain spatially resolved spectra, the BLR's structure and dynamics are not understood in detail.Therefore, the virial factor, f , accounts for the unknown geometry and kinematics, and ∆v 2 R BLR /G is the virial mass, M •, vir (Homayouni et al. 2020;Pancoast et al. 2014;Netzer & Marziani 2010).Comparisons of M •, vir from the AGN that have been reverberation mapped to M • determined with the M • − σ * scaling relationship indicate f ≈ 4.5 (Woo et al. 2015;Batiste et al. 2017).∆v can be estimated from broad emission line (BEL) profile widths, but R BLR cannot be directly determined through imaging for most AGN (with a few exceptions discussed below).Instead, R BLR can be determined by the reverberation mapping technique, which utilizes the response of the broad lines to the AGN's continuum variability and is not fundamentally limited by distance.BLR reverberation mapping is currently the only practical means of probing the SMBH mass function over a large range in redshift and luminosity (Shen 2013, and references therein).
Arguably the most important result from BLR reverberation mapping is the tight relationship found between the measured lag (t ′ ) and the AGN luminosity, which implies a BLR radius -AGN luminosity relationship, t ′ ∝ R BLR ∝ L 0.5 AGN (e.g., Bentz et al. 2009Bentz et al. , 2013)).This provides the key for determining M • for very large samples, since both the AGN continuum flux and the BEL profile width can be measured from a single spectrum, yielding L AGN (hence R BLR ) and ∆v, respectively (e.g., McLure & Jarvis 2002;Vestergaard & Peterson 2006;Liu et al. 2012).
In a few nearby and luminous AGN, the BLR has been spatially resolved with near-infrared (NIR) interferometry.Very high resolution interferometry observations of 3C 273, IRAS 09149-6206, and NGC 3783 have been obtained with the GRAVITY instrument using the European Southern Observatory Very Large Telescope Interferometer (VLTI).These data suggest that all 3 AGN have a rotating, disk-like BLR.From the broad Paschen-α line of the quasar 3C 273, Sturm et al. (2018) inferred R BLR = 0.12 ± 0.03 pc and M • = 2.6 ± 1.1 × 10 8 M ⊙ .Amorim et al. (2020) spatially resolved the broad Brγ line in IRAS 09149-6206 and inferred R BLR = 0.075 pc and M • ∼ ×10 8 M ⊙ .More recently, Amorim et al. (2021b) used VLTI data to map broad Brγ in NGC 3783 and found a BLR with a radial distribution of clouds within R BLR ≈ 0.014 pc and then in Amorim et al. (2021a) reported M • = 2.54 +0.90  −0.72 × 10 7 M ⊙ .For all 3 objects, the BLR measurements obtained from GRAVITY are consistent with the R − L relation and reverberation mapping studies of those same AGN.

Reverberation Mapping
In the "point-source reverberation model", clouds respond to fluctuations in the continuum flux on time scales of hours, much smaller than the light crossing time of the whole BLR, which is on the order of days to weeks.The differences in time for a photon to travel between the accretion disk, a given cloud, and then (as a line photon) to an observer means there is a corresponding time lag between the observed continuum variations and the response of the BELs.From the perspective of a distant observer, the clouds that respond at the same delay time form parabolic isodelay surfaces, where r is the cloud's radial distance from the continuum source and Θ is the angle between the observer's line of sight (LOS) and the cloud's position vector.The observed variation of emission line luminosity with time, t, at LOS velocity, v || , is given by the convolution of the continuum luminosity, L c (t − t ′ ), with a transfer function, Ψ(t ′ , v || ), where L(t, v || ) is the velocity-resolved response of the emission line.The transfer function is the emission line's response to the continuum at t = t ′ and encodes information about the BLR's geometry and velocity field (Blandford & McKee 1982).Typically R BLR is recovered by cross-correlating the observed line light curve with the continuum light curve (L(t) and L c (t − t ′ ), respectively) to measure the response-weighted time lag, t ′ RW , which is assumed to be equivalent to the luminosity-weighted radius of the BLR; R LW ≈ ct ′ RW (Robinson & Perez 1990;Koratkar & Gaskell 1991;Pérez et al. 1992;Almeyda et al. 2020).The value of R LW inferred is then used in Equation 1 to find the SMBH mass; i.e., R BLR = R LW .However, due to the complex properties of the BLR, ct ′ RW is not necessarily an accurate estimate of R BLR .As part of this work, we compare the relation between the ct ′ RW and R LW values for our BLR models to understand how well the measured lag represents R LW .

Modeling Reverberation Responses
With sufficient data, the transfer function can be recovered using deconvolution techniques, such as with the Memecho code by Horne et al. (2004).In practice, however, deconvolving the transfer function accurately from an observed response is not easily done, as it requires very well sampled data over a sufficiently long observing period.Furthermore, even if the transfer function could be reliably obtained, deciphering the BLR's geometry and velocity parameters would not be a straightforward task, since they would have to be inferred from the shape of the transfer function, which would probably require modeling.Instead we can use computer simulations to predict L(t, v || ), given an input AGN light curve, geometry, and velocity field parameters.The simulated response light curves and BEL profiles for different BLR configurations can then be compared to observed AGN light curves and time-domain spectra.Numerical reverberation mapping studies of the velocity-resolved response and the transfer function have been presented by Welsh & Horne (1991); Pérez et al. (1992); Robinson (1995), and Peterson & Horne (2004).More complex reverberation mapping codes, such as the Code for AGN Reverberation and Modeling of Emission Lines (Caramel) by Pancoast et al. (2011) and, for the AGN torus, the TOrus Reverberation MApping Code (TORMAC) (Almeyda 2017), provide more detailed simulations.
TORMAC, developed by Almeyda et al. (2017) and further expanded by Almeyda et al. (2020), simulates the response from the dusty torus for select IR wavelengths given any input light curve.We have adapted and extended TORMAC for application to the BLR.This new code, the Broad Emission Line MApping Code (BELMAC), generates the velocity-resolved response for various BLR geometrical configurations, gas properties, cloud distributions, and velocity fields, given any input light curve.The overall goal is to develop BELMAC into a fast and flexible code that can be used to generate large model grids to fit reverberation mapping observations.However, it is useful to understand how various properties of the BLR affect the reverberation response (e.g., Welsh & Horne 1991;Pérez et al. 1992;Horne 1994).In this paper, we introduce the basic features of BELMAC and explore the effects of parameters controlling the BLR geometry, kinematics, and cloud properties on the velocity-resolved transfer function.As our main focus here is on the general properties of the BLR we simplify the calculation of the line emission by using hydrogen recombination theory.A narrow square-wave pulse is used as the as the driving light curve to produce a numerical approximation of the transfer function.
In forthcoming papers, we will describe the next version of BELMAC, which incorporates a grid of photoionization models for more sophisticated emissivity calculations.
BELMAC is outlined in Section 2 with a description of the geometry setup, velocity field, and the calculation of the cloud emission.In Section 3 we present our results for two example models, a spherical BLR with a radiative pressure driven radial outflow and a rotating thin disk.We discuss how our simple reprocessing models can help to understand the behavior of responses and timedomain spectra, the parameters that affect measured time lags, and how BELMAC compares to other reverberation mapping codes in Section 4. We also discuss how well our results represent various BELs.Finally, we summarize our conclusions in Section 5.

OUTLINE OF THE CODE: BELMAC
The BLR is represented as a 3D ensemble of discrete clouds, that are randomly distributed within a user-defined structure.BELMAC follows the same geometry set-up as TORMAC, but the clouds within the ensemble are dust-free, photoionized gas clouds and have velocities specified by various choices of velocity field.The general flow of BELMAC is illustrated in Figure 1 and described in detail in the following subsections.The user provides the AGN's bolometric luminosity, L AGN , spectral energy distribution (SED), and driving light curve.The geometry, velocity field, and distribution of the ensemble of BLR clouds are described by Y BLR , σ, i, C f , s, p, and M • .The cloud positions are randomly generated and each cloud's LOS velocity is calculated based on the type of velocity field specified, either Keplerian, radial flow, or turbulent.In all, there are 10 descriptive parameters for each distinctive BLR model.BELMAC can use any input light curve, but in order to obtain a numerical approximation of the transfer function, the input light curve is a single square-wave pulse to represent a δ-function.We refer to the approximate transfer function as the response function (Almeyda et al. 2017).For the duration of the pulse, ∼ 4 days, the ionizing luminosity increases by a factor of 2. The light-front corresponding to the light pulse propagates outward from the center of the system with time and the total luminosity of the BLR is computed by integrating over the cloud ensemble at each observer time step, taking into account light travel delays.

BLR Geometry and Cloud Structure
In common with many other studies (e.g., Armijos-Abendaño et al. 2022;Marconi et al. 2008;Robinson & Perez 1990;Davidson & Netzer 1979), we assume that the BLR gas is distributed in many small clouds.The ensemble of discrete, gas clouds are randomly distributed within an overall geometry of either a disk or a sphere.Cloud positions are defined by the spherical polar coordinates r, θ, ϕ, where r is the radial distance from the continuum source at the origin, θ is the polar angle, and ϕ is the azimuthal angle.In the models considered here, the BLR is assumed to have 'sharp' boundaries with the clouds uniformly distributed in elevation above the equatorial plane, 90 o − θ, and in ϕ.However, 'fuzzy' edges may also be arranged by drawing values of the elevation angle from a Gaussian distribution (Figure 2; Almeyda et al. 2017).The disk's angular width is set by the half-angle σ, where σ = 90 o is a spherical ensemble (Figure 2).The inclination, i, is the angle between the polar axis and the observer's LOS, therefore, for i = 0 o the disk is viewed face-on and for i = 90 o it is viewed edge-on.
The outer radius of the BLR is set by the dust sublimation radius, R d , where T sub is the dust sublimation temperature (Nenkova et al. 2008a).Here we take representative values, L AGN = 10 45 erg s −1 and T sub = 1500 K, then R d ≈ 10 18 cm (0.4 pc).For simplicity, we will treat the dust sublimation radius as a sharp boundary between the BLR and torus, rather than as a transition region, which would be more physically realistic (Baskin & Laor 2018).The inner radius of the BLR, R in , is set by the free parameter In all of the models presented here, we will use Y BLR = 20; similar to scaled sizes used in previous modeling (e.g.Du et al. 2015;Netzer 2020), giving R in ≈ 6 × 10 16 cm (0.02 pc), which corresponds to a light-crossing time ∼ 24 days.This is consistent with observed BLR time-delays that range from days to ∼a month (e.g., Fausnaugh et al. 2017) and the R − L relation (Bentz et al. 2013).The radial cloud distribution is given by a power-law with index p, such that the number of clouds between r and r + dr is, where A o is a normalization constant that is determined by the total number of clouds, (Almeyda et al. 2017).The BLR clouds are assumed to be spherical, with radii, R cl (r), specified by the integrated covering fraction of the BLR, where the cloud's cross-sectional area is A cl (r) = πR 2 cl (r).We assume the clouds are composed of pure hydrogen and have a constant mass, m cl (r) = m cl .For a pure hydrogen gas cloud, the ionized portion of the cloud has an electron density equal to the hydrogen density, which is the cloud's gas density n(r).The gas density is a power-law function of the cloud's radial distance from the continuum, n(r) = n o (r/R d ) s , where the values of the gas density at the outer radius of the BLR, n(R d ) = n o , and the power-law index, s, are free parameters.Therefore, the cloud radius R cl (r) ∝ r −s/3 and the cloud column density, will vary with distance from the ionizing source.Realistically, the mass of the clouds could change with distance depending on how the clouds are confined, such as by radiation pressure or magnetic fields (Netzer 2008).Table 1 summarizes the BELMAC geometry parameters with their respective descriptions, symbol, and values we will adopt as "standard".

Cloud Emission
To approximate the Hα luminosity of a cloud, we use hydrogen recombination theory (Osterbrock & Mathews 1986).For spherical clouds, d c (r) is the average path length through the cloud, which is proportional to R cl (r).The surface of a cloud facing the source is exposed to the ionizing flux and is entirely ionized.Assuming ionization equilibrium, the clouds are fully ionized within the Strömgren depth, where α B is the recombination coefficient for hydrogen (α B ≈ 2.6 × 10 −13 cm 3 s −1 ; Osterbrock & Mathews (1986)) and N s (r) is the Strömgren column density.The ionization state of a given cloud is described by the ionization parameter, the ratio of the local ionizing photon flux, Φ(r, t ′ ), to the cloud's gas density, where Q H (t ′ ) is the total ionizing photon luminosity.Since the continuum luminosity naturally undergoes temporal variability, Q H (t ′ ) and, hence, U (r, t ′ ) fluctuate with time.Q H (t ′ = 0) is obtained from any selected AGN's SED and L AGN , which sets U (R d , 0).To estimate Q H (0), for the models presented in this paper, we use the compiled SED created by Jin et al. (2012) who constructed an average SED using 17 type 1 AGN that have an average Eddington ratio Γ = 0.07 and , the cloud is fully ionized, or "matter-bounded", and has reached its maximum ability to emit hydrogen recombination line radiation.Otherwise a cloud is partially ionized and "radiationbounded."The cloud luminosity for radiation and matter-bounded clouds is, respectively, where α ef f Hα is the effective recombination coefficient for Hα (α ef f Hα ≈ 1.1 × 10 −13 cm 3 s −1 ).In the matter-bounded case for constant mass clouds, the cloud luminosity L cl (r) ∝ n(r) ∝ r s , but is independent of U (r, t ′ ) and therefore is not time-dependent.In the radiation-bounded case, the cloud luminosity does depend on U (r, t ′ ) and thus The cloud luminosities are summed at each time intervals in the observer's frame to calculate the total BLR luminosity as a function of time.
At the high densities of the BLR, some emission lines, such as the Balmer lines, may become optically thick (Osterbrock & Mathews 1986).It is also possible that dust grains may survive in the cooler, largely neutral gas beyond the ionization front (e.g., Baskin & Laor 2018), causing extinction of line radiation produced in the ionized zone.As a result, the line emission of a cloud may vary strongly between the illuminated and non-illuminated faces.To account for this effect, BELMAC includes an approximate treatment of anisotropic cloud emission (ACE), which is dependent on a cloud's location with respect to the continuum source and the observer's line of sight.Each cloud's luminosity is multiplied by its emission line fraction (elf ), which is the fraction of the illuminated-side an observer sees.However, if a cloud is matter-bounded (d s ≥ d c ) then we assume that it has isotropic cloud emision (ICE) and set elf = 1. Figure 3 shows from the observer's perspective how the elf changes the cloud's luminosity with position.
We summarize the parameters governing a cloud's emission and internal properties, along with their typical values in Table 1.

The Analytical Transfer Function
To understand how the basic features of the response function are influenced by the size of the BLR, the gas density, and radial cloud distribution, we discuss the analytical transfer function for a spherical shell.The radial cloud distribution and internal cloud properties determine the volume emissivity, ε V (r) ∝ L cl N (r)/r 2 , which will affect the shapes of the line profile and transfer function.We adopt Pérez et al. (1992)'s assumption that the cloud's emission efficiency behaves as a power-law with distance and approximate the volume emissivity as, where and the emissivity index η relates to our p and s indices (from Equations 5 and 7, respectively) as, for radiation and matter-bounded clouds, respectively.Recalling Equations 2 and 3, the analytical transfer function is given by, where V is the volume of the BLR.Solving Equation 14for a spherical shell BLR model we arrive at, and where τ = ct ′ /2R d is the delay in units of the light crossing time of the BLR, such that the line response is complete when τ = 1 (Robinson & Perez 1990;Pérez et al. 1992;Almeyda et al. 2017).In terms of τ , the duration of the square-wave pulse used in our models is 0.004.

Response Weighted Delay and Luminosity Weighted Radius
As mentioned in Section 1.2, the response-weighted time delay determined from reverberation mapping is considered a reasonable estimate of the radial extent of the BLR, ct ′ RW ≈ R LW = R BLR .However, it is unclear how accurate this assumption is due to the uncertain and complex cloud properties of the BLR.Robinson & Perez (1990) defined the dimensionless, response-weighted delay (RWD) as, to characterize the delay associated with the transfer function.This is related to the cross-correlation lag as t ′ RW = 2R d RWD/c in physical units.For a spherical shell containing isotropically emitting clouds, the RWD is exactly the dimensionless luminosity-weighted radius (LWR), (Robinson & Perez 1990;Almeyda et al. 2020).The analytical LWR and RWD/LWR for spherical shell BLRs with varying Y BLR , s, and p are shown in Figure 4. We will present and compare the LWR and RWD of various BLR models in Section 3.

BLR Cloud Dynamics
In order to compute the velocity-resolved response, we include three types of velocity fields; rotational, radial, and random turbulence.The broadest line widths observed are ≲ 0.1c (e.g., Assef et al. 2011;Fausnaugh et al. 2017), hence the dynamical time scale of the BLR is much greater than (≳ 10×) its light-crossing time.Therefore, we assume the clouds are static as they do not significantly move along their trajectories over timescales typically of interest for reverberation mapping.BELMAC is capable of combining radial, rotational, and random motions for a multi-component velocity field, however, for the purposes of this paper, we will present and discuss each type of motion separately.Furthermore, we will not explore the random, turbulent motion in the models presented here.

Rotational Velocity Field
In the rotating disk BLR, the clouds are assumed to follow circular, Keplerian motion, v Kep (r) = (GM • /r) 1/2 .The rotational velocity is scaled to the velocity at the inner radius, for a cloud at radius, r.The cloud orbits are randomly inclined relative to the disk's mid-plane and are constrained by the disk angular width, σ, which requires 90 o − σ ≤ |θ| ≤ 90 o + σ.Given a cloud with coordinates r, θ, and ϕ and a BLR inclination i, the velocity component along the LOS velocity to an observer is, v If we combine Equations 2 and 20 for an infinitesimally thin disk (σ → 0, θ = 90 o ), we find the response for a given orbital radius in the projected velocity and time delay space forms an ellipse.The ellipses are centered at (t ′ = r/c, v || Kep = 0), with semi-axes ±v rot sin i in velocity and r/c sin i in time delay (Pérez et al. 1992;Welsh & Horne 1991).

Radial Velocity Field
The radial motion is dependent on the primary outward driving mechanism, assumed here to be radiation pressure, and the inward pull of gravity.We further assume that the BLR intercloud medium is sufficiently tenuous that we can ignore the forces due to drag and the pressure gradient between the clouds.We follow Marconi et al. (2008) . Rotational and radial velocities with respect to distance scaled to R d .The escape velocity at R in (dashed gray line) is 6,578 km s −1 .The blue curves show the radial velocity given by Equation 24.The lighter the shade, the lower the s index.The solid, dashed, and dotted lines are p = 0, 1, and 2, respectively.The blue circles on the s = −1 and −2, p = 0, 1, and 2 curves mark the distance at which radiation pressure overcomes gravity and the clouds start accelerating outwards.Each of these points corresponds to the critical N col = 5.26 × 10 22 cm −2 .When s = 0, the outward, radiation pressure force is greater than gravity at all distances.The green curve is the rotational velocity field, v Kep , described by Equation 19, which does not depends on the s and p indices.
The equation of motion for a cloud of pure hydrogen gas at a distance r from the radiation source, with cross-sectional area A cl (r) and mass m cl is, where α f ≃ ∞ ν H L ν dν/L AGN , the fraction of the AGN's bolometric luminosity that is absorbed by the cloud.Although the absorption fraction will vary with U (r, t ′ ) and N col (r), in this paper we use the constant value α f = 0.5, following Netzer & Marziani (2010), throughout the BLR.Putting Equation 22 in terms of N col and recognizing that the Eddington luminosity is after simplifying by introducing the Eddington ratio, Γ = L AGN /L Edd , and the force multiplier, , where σ T is the Thomson cross-section.It can now be easily seen that for ΓF M > 1 the cloud is accelerated outwards and accelerated inwards when ΓF M < 1 (Marconi et al. 2008;Netzer 2008;Netzer & Marziani 2010).The radii at which ΓF M = 1 are when N col = 5.26 × 10 22 cm 2 are are indicated in Figure 5; for the s = −1 models at r = 2.55 × 10 17 cm, 7.12 × 10 17 cm, and 5.84 × 10 17 cm, and for s = −2, at r = 8.15 × 10 17 cm, 9.12 × 10 17 cm, and 9.61 × 10 17 cm for p = 0, 1, and 2, respectively.If the BLR is created by a wind from the accretion disk (e.g., Elvis 2000), the gas clouds at R in would have a velocity comparable to the escape velocity from the disk.Therefore, we assume the radial velocity at R in to be the local escape velocity.Integrating Equation 23we obtain, the radial velocity field.Figure 5 shows v Kep (r) and v rad (r), Equations 19 and 24, respectively, for several values of the gas density and cloud radial distribution power-law indices.
The component of the radial velocity along the observer's LOS is given by, Combining Equations 25 and 2 the delay-LOS velocity relation for a thin, spherical shell with radius r is a straight line in velocity and time delay space, where the slope, r c v rad , is positive for an outflow and negative for an inflow (Welsh & Horne 1991;Pérez et al. 1992).

RESULTS
Here we present the velocity-resolved reverberation response functions (2DRFs) of, firstly, a spherical shell BLR with a radiation pressure driven outflow and, secondly, a rotating disk BLR, for their respective parameters listed in Table 3.We also present the velocity-integrated response functions (1DRFs), time-averaged line profiles, and root-mean-square (RMS) line profiles.
To isolate the response amplitude, we subtract the initial BLR state (prior to the onset of the continuum pulse) from the response function at each τ .The 1DRF is normalized as L(τ where L o is the BLR luminosity in its initial state and L max is the peak luminosity of the BLR's response.We will present results from BLR models with both isotropic cloud emission (ICE) 10 a Calculated from the SED model by Jin et al. (2012) for M • = 10 8 M ⊙ and Γ = 0.1.See their Figure 14, panel 4-C.
and anisotropic cloud emission (ACE).The BLR models presented here include N tot = 10 5 clouds, but as the BLR probably contains ∼ 10 6 clouds (Arav et al. 1998;Dietrich et al. 1999), these can be considered a representative sub-sample.

Spherical BLR with Radial Outflow
In Figure 6 we show the 2DRF for an ICE BLR model that has a spherical geometry and a radiation pressure driven outflow, with gas density and cloud distributions n(r) ∝ r −1 and N (r) = constant (defined by s = −1 and p = 0), respectively.The left panel shows the parabolic isodelay surfaces corresponding to the continuum pulse at various delays (Equation 2) as it propagates through the BLR.The right panel shows the 2DRF, 1DRF, average line profile, and RMS line profile.At τ = 0, only clouds located along the LOS from the center of the BLR to the observer respond and these have the largest blueshifted velocities.The strongest changes in the line profile are shown by the RMS profile, which is peaked and skewed to the blueshifted side of the line (far-right sub-panel).The response remains entirely blueshifted until τ ≥ 1/2Y BLR = 0.025 (Equation 2 with r = R in , Θ = 90 o ), when the response-front crosses into the far-side hemisphere.The response covers the greatest spread in LOS velocity at the light-crossing time of the inner radius, τ = 1/Y BLR = 0.05 (bold curve in left panel), when it reaches the maximally redshifted LOS velocity.This also corresponds to the width of the flat-topped portion of the 1DRF (bottom-right sub-panel).The response then narrows in velocity space as τ increases and becomes entirely redshifted for τ ≥ 0.5; the instantaneous line profiles are now red-asymmetric.After crossing the inner cavity, the response amplitude decays as τ increases until the BLR light-crossing time is reached (τ = 1.0).It can be seen from Equation 14that the relatively rapid decay is a consequence of the steeply declining emissivity distribution, where ε V ∝ r −10/3 (Equation 12), in this case.Although changes in the line profile begin in the blue wing then propagate redward as τ increases, the time-averaged line profile is symmetric, since in this ICE model there is an equal amount of line emission from the redshifted and blueshifted hemispheres of the BLR.The 2DRFs of spherical shell BLR models, with p = 0, 1, 2 and s = 0, − 1, − 2 are shown in Figure 7 for both ICE and ACE cases.Their corresponding 1DRFs and delay-averaged line profiles are shown in the left and right panels of Figure 8, respectively.Considering first the ICE models with s = −1 and −2 (middle and bottom rows of Figures 7 and Figure 8), it can be seen that as p is increased from p = 0 to p = 2 the response narrows in velocity space and it's amplitude declines more gradually with increasing τ .For the larger values of p, the clouds with the highest velocities are located near R in (see Figure 5) and make-up a small fraction of the ensemble.These clouds, therefore, make a relatively small contribution to the response, which is dominated by clouds at larger radii, resulting in a response that is quite narrow in velocity-space.The more gradual decline in the tail (τ > 1/Y BLR ) of the 1DRF is due to the more gradual decrease in emissivity with r, where for s = −2, p = 0, 1, 2, ε V ∝ r −8/3 , r −5/3 , and r −2/3 .The 1DRF's of the ICE models in Figure 8 have a flat-topped portion at short delays, but the widths vary greatly between the s = 0 and s = −1 models.Given that n(R d ) = 10 9 cm −3 and that cloud sizes are constant when s = 0, a portion of the clouds are matter-bounded in the BLR's initial state, prior to the onset of the continuum pulse, and thus do not respond to the increase in the ionizing flux (see Equation 10).There is a transition from matter-to radiation-bounded clouds at Table 3. Statistics for matter-bounded clouds in the s = 0 models.The first columns is the cloud distribution index parameter p.The 2 nd − 4 th columns are, respectively, the percentage of matter-bounded clouds, the distance at which the clouds are radiation-bounded, and scaled radial depth, R d /R M BC , prior to continuum pulse initiates.The 5 th − 7 th columns are, respectively, the percentage of radiation-bounded clouds that become matter-bounded, the greatest distance the bound state transition occurred, and the apparent scaled radial BLR size after the response completes.The true inner radius is 6.17 × 10 16 cm and scaled size Y BLR = 20, for all BLR models presented here.

Matter-Bounded Clouds
, which is smaller than the true value, Y BLR .The transition radius varies with p and is listed in Table 3 along with the fraction of matter-bounded clouds during the initial state.Even a small matter-bounded cloud fraction corresponds to an R M BC that is several times larger than the true R in .For example, when p = 2, only ∼ 1% of clouds are matter-bounded and don't respond, but Y mbc BLR ≈ 4. In addition to the clouds initially emitting at full capacity, there are clouds with d s ≲ d c that are almost matter-bounded and become matter-bounded when the response-front reaches them.These clouds are at distances r > R M BC and cause a turn-over feature in the response.Over the course of the response, the fraction of radiation-bounded clouds that transition to matter-bounded for p = 0, 1, and 2 are 28%, 13%, and 4%, respectively.The clouds transitioning to the matter-bounded state further increases the apparent inner BLR radius to R apt M BC , which increases the width of the flat-topped portion of the 1DRF to 1/Y apt BLR .For example, Figure 9 shows the 2DRF, 1DRF, and averaged line profile for p = 0 and s = 0. Overplotted on the 1DRF are two analytical 1DRFs calculated using Equations 15 and 16 with Y mbc BLR = 2.60, determined from the transition radius R M BC , and Y apt BLR = 2.27, which was determined by R apt M BC .Although these values are similar, the analytical curve for Y apt BLR = 2.27 is clearly a better match to the 1DRF.The values of Y mbc BLR and Y apt BLR for the other BLR models are listed in Table 3.

Anisotropic Cloud Emission
Figures 7 and 8 also show the 2DRFs, 1DRFs and τ -averaged line profiles for the ACE models, where ε V also depends on elf and therefore becomes a function of θ and ϕ, as well as r.For ACE models in which all the clouds are radiation-bounded (as is the case for s = −1 and −2), the observer sees < 50% of the illuminated face for the clouds on the near-side of the BLR.Therefore, at short time delays, the blue-wings of the line profiles are suppressed and the 1DRF is no longer flat-topped.The response peak occurs at delays when the majority of clouds on the far-side are along the LOS (i.e., elf > 1) are responding, that is, when τ = 1/Y BLR for p = 0 and at later delays for the p = 1 and 2 models, where the clouds are more uniformly distributed.The τ -averaged line profiles of the ACE models are highly asymmetric and have half the amplitude of the corresponding ICE profiles.Clouds that are matter-bounded are assumed to emit isotropically and, therefore, have their elf values set to 1.The models that include matter-bounded clouds (which occurs for s = 0) therefore have 1DRFs that are more similar to their ICE counterparts.In the line profiles, the asymmetry between the red and blue wings is smallest in the p = 0 model, which is the model having the highest fraction of matter-bounded clouds, ∼ 34%.The largest asymmetry occurs in the s = 0, p = 2 model, which has the lowest fraction of matter-bounded clouds of ∼ 1%.

The Response Weighted Delay and Luminosity Weighted Radius
Figure 10 shows the variation with Y BLR of the LWR, RWD, and the ratio LWR/RWD for ICE and ACE BLRs, and for each combination of s and p.
For ICE models with only radiation-bounded clouds (s ̸ = 0), the LWR and RWD both decrease with increasing Y BLR , similar to the analytical relations shown in Figure 4, and the ratio LWR/RWD = 1.Additionally, since LWR is not dependent on elf , the LWR values for the ICE and ACE models are identical.However, the response peaks at longer delays for ACE models, as discussed in Section 3.1.2,and thus the RWD is greater than for the corresponding ICE models.Therefore, for ACE models with s ̸ = 0, LWR/RWD ≈ 0.8 for all p and Y BLR values.
For models with s = 0, the RWD increases with Y BLR due to the presence of non-responding, matter-bounded clouds.The largest increase in RWD occurs for s = 0, p = 0, when about ∼ 30% of the cloud ensemble is matter-bounded.As a result, LWR/RWD decreases as Y BLR increases, implying that the RWD overestimates the LWR and is dependent on the fraction of matter-bounded clouds and hence also Y apt BLR ).The overestimation of LWR due to the presence of matter-bounded clouds ranges

Keplerian Disk BLR
Here we present the results of models in which the BLR is configured as a Keplerian disk.As specific examples, we discuss the reverberation responses for the ICE and ACE cases for an edge-on (i = 90 o ), rotating thin disk (σ = 15 o ), where the gas density decreases as r −2 , U is constant (s = −2), cloud size increases with distance, and the cloud radial distribution is constant (p = 0, N (r) constant).In this case, all of the clouds are radiation-bounded.
Figure 11 displays the 2DRF, 1DRF, time-averaged, and RMS line profiles of the ICE model.For a specific orbital radius r, the projection of the isodelay surfaces and orbital velocities in the disk's plane is an ellipse in the 2DRF (Equation 21; Pérez et al. 1992).As can be seen in Figure 11, the characteristic rocket or bell-shaped 2DRF is created by the superposition of ellipses of varying r.Clouds at R in have the largest LOS velocity range and respond at a shorter time delay than the clouds located at larger radii.The major axis of the corresponding ellipse is therefore aligned with the velocity axis.Conversely, the ellipse associated with the clouds at R d , has the smallest velocity range and spans the largest range in time delay, and hence it's major axis is aligned with the time delay axis.Rather than the flat-top that was seen for the 1DRF of the spherical BLR, the disk 1DRF has a sharp peak at τ = 0.0, then a second, weaker peak in the response at the light-crossing  21.The system is entirely radiation-bounded clouds.This figure has an animated version that is available in the online version of this paper, which does not display the ellipses.In the animation, the 1DRF is plotted with τ from 0 to 1, along with the τ −varying line profile (gray).The changing line profile is a narrow, bright peak at τ = 0, then quickly becomes lower in amplitude and broader with double peaks at the extremities of the profile's wings.These low-amplitude peaks slowly move back towards the center of the profile until τ = 1.The τ −averaged line profile (black), RMS (blue), and, in the animation, the τ −varying line profiles (gray) are normalized to their respective maxima.
time of the inner radius, τ = 0.05, due to the decrease of clouds while the response-front crosses the inner cavity.In this model, ε V ∝ r −8/3 and as τ increases, more clouds at larger radii respond and, the tail of the 1DRF smoothly decreases until the light-crossing time of the BLR is reached (τ = 1).For an edge-on rotating disk, the τ −averaged line profile has distinctive, symmetrical redshifted and blueshifted peaks since the emission from the receding and approaching sides of the disk contributes equally.The response begins at v || rot ≈ 0 km s −1 and the amplitude is large enough  11, including an animated version that is available online, but for an ACE model.In the animation, the τ −varying line profile is a narrow, low-amplitude peak at τ = 0, then quickly becomes brighter in amplitude and broader with double peaks at the extremities of the profile's wings.These high-amplitude peaks slowly move back towards the center of the profile until τ = 1.
that the RMS profile has a central peak.The response then moves quickly and simultaneously to each wing, contributing to the weaker, symmetrical peaks in the RMS profile.
The effect of ACE decreases the response amplitude during short time delays because the observer sees very little of the illuminated faces of the responding clouds.The initial response amplitude is diminishes by almost half compared to the amplitude in the ICE 1DRF model in Figure 11.Instead, in the ACE model, the response peaks at τ = 0.05, when the observer receives the response from the clouds with elf > 0.5.For the same reason, although the response begins in the core of the line profile, the effect of ACE decreases the response amplitude and the RMS profile doesn't have a central peak.
Figure 13 shows the 1DRFs and line profiles for disk models inclined at i = 60 o with varying values of the s and p parameters for both ICE and ACE cases.We chose to present the i = 60 o BLR models, since an edge-on BLR would be obscured to an observer by the surrounding dusty torus.However, we still want to present a high inclined disk to show distinctly separated peaks in the line profiles.Other choices of intermediate inclinations produce qualitatively similar results.When s = −1 and −2, all of the clouds are radiation-bounded.For these cases of ICE models, ε V (r) decreases less rapidly as s decreases and p increases, therefore the response decays more gradually with τ .In the ACE models, the clouds on the far-side relative to an observer contribute more emission to the total response, so the peak response occurs at longer time delays as p increases.When clouds are concentrated near R in (i.e., p = 0), the peak is at τ = 1/Y BLR .
Both the ICE and ACE models have symmetrical line profiles with redshifted and blueshifted peaks, characteristic of rotating disks.Since the elf is symmetrical about the disk axis, the ICE and ACE line profiles have the same shape, but the ACE model profiles have half the amplitude.
When s = 0, a portion of the clouds are matter-bounded prior to the onset of the continuum pulse and hence do not respond.Consequently, both peaks in the 1DRF are broader and the second occurs at τ ≈ 1/Y apt BLR .The percentages of matter-bounded clouds, R M BC , and Y apt BLR for each model are listed in Table 4. Since elf = 1 (isotropic emission) when a cloud is matter-bounded, the ratio of the ICE and ACE line profiles amplitudes approaches 1 as the fraction of matter-bounded clouds increases (see Table 4).
The LWRs and RWDs for a disk with i = 60 o , shown in Figure 14, show the same general trends with Y BLR as the spherical BLR.In particular, the ICE models (left panel) and the corresponding spherical BLR models (left panel in Figure 10) show very similar behavior.However, the RWDs for the ACE models are longer for a disk than a sphere and therefore the LWR/RWD ratio has slightly  More details on the results of models with matter-bounded clouds is in Section 3.2.2.

Effects of Inclination
Cross-sectional diagrams of the disk-like BLR ensembles at i = 0 o , 45 o and 90 o with isodelay surfaces over plotted are shown in Figure 15.The effects of inclination on the 1DRFs and line profiles are shown in Figure 16 for both ICE and ACE models with s = −2 and varying p.When the disk is face-on to the observer, the 1DRF is zero until τ = 0.03 when the response-front reaches the clouds at the inner-radius.As i increases, the response begins at progressively shorter delays until, at i = 90 o − σ the response begins at τ = 0. Additionally for ICE models, the 1DRF peak becomes narrower and the tail overall extends to longer delays, as i increases.The ACE 1DRFs show much stronger, extended responses at larger inclinations than their ICE counterparts.
When the disk is face-on, there is a sudden decrease in the 1DRF at τ = 0.37.The inflection point corresponds to the delay, τ = 1 2 (1 − cos (90 o − σ)), when the isodelay surface intersects with the outermost edge of the disk surface facing the observer, shown in Figure 15 as a dashed line.The 1DRF peaks at this τ -inflection point when p = 1 and 2, as the largest number of clouds are responding at this delay.The response is complete at τ = 1 2 (1 − cos (90 o + σ)).The right panel of Figure 16 shows the line profiles for the same models, where the profile amplitudes have been normalized to the maximum amplitude for the ICE model.With respect to the line center, the peaks are located at ± sin i v min /v max , therefore the gap between them is zero at i = 0 o resulting in a single-peaked profile (Robinson 1995).Otherwise, when the disk is inclined, the line profile  shape has the characteristic, symmetrical peaks about the center of the profile.The trough between the peaks also gets deeper with increasing i and p.

Models with Matter-Bounded Clouds
Figure 17.The same as Figure 16, but with matter-bounded clouds in the ensemble (s = 0).
Figure 17 is the same as Figure 16, but for models with s = 0. Here, clouds within R M BC are matterbounded and do not respond.This creates an "apparent" inner cavity in the BLR that is larger than R in .Therefore, the response has a delayed start that is dependent on R apt M BC in addition to its dependence on i and σ, as discussed above.However, as R M BC < R d , this does not affect the clouds at the largest radii in the BLR and the response ends at the same delay, τ = 1 2 (1 − cos (90 o + i + σ)), regardless or whether or not matter-bounded clouds are present.The 1DRFs are double-peaked at i ≥ 30 o , where each peak represents the response-front encountering the near-side then far-side clouds at R M BC .The time-delay between the two peaks lengthens as i increases, since the second peak occurs at τ = (2Y mbc BLR ) −1 (1 − cos (90 o − i)) and also as p decreases due to the increase in the fraction of matter-bounded clouds.
In general, the ICE and ACE line profiles in the right panel of Figure 17 have a similar shape.However, the presence of matter-bounded clouds creates shoulder-like features in the line profiles of ACE BLR models when i ≥ 60 o .This is most clearly seen for the model with i = 90 o and p = 0, the symmetrical inflections are a result of matter-bounded clouds within R M BC emitting isotropically, then, as the isodelay surface propagates throughout the BLR, radiation-bounded clouds at larger radii respond and emit anisotropically.The width of these shoulders relates to R M BC and thus increases as p increases.The width of the shoulders broadens with increasing p and become more smoothed with decreasing i, therefore these features we do not clearly see for p > 1 and i ≤ 45 o .

DISCUSSION
We have presented a new BLR reverberation mapping forward-modeling code, BELMAC, that incorporates various velocity field models to calculate the velocity-resolved response.We assume the clouds are bound spheres, which follows similar modeling work such as Lawther et al. (2018); Goad et al. (2012); Marconi et al. (2009); Netzer (2008); Robinson & Perez (1990).The version of the code used here assumes hydrogen recombination to determine cloud luminosity, where the clouds can optionally emit isotropically (ICE) or anisotropically (ACE).We presented BLR cloud ensemble models with a radiation pressure driven outflow in a sphere and Keplerian orbits in a disk, respectively.
BELMAC does not provide a full hydrodynamic treatment of the gas flows and does not consider inter-cloud radiative transfer.However, representing the BLR as a cloud ensemble, offers the ability to quickly and efficiently explore multiple configurations of BLR geometries and velocity fields.It also allows for compatibility with photoionization model grids, which will be implemented in the forthcoming version.This representation may impose some limits on BELMAC's predictions.Ideally, reverberation mapping modeling codes would include a detailed treatment of hydrodynamics, incorporating radiative transfer and photoionization physics.However, these codes do not yet exist and would no doubt be computationally expensive.

The Radiation Pressure Driven Flow
The BLR clouds are photoionized and therefore subjected to radiation pressure forces due to scattering and absorption.Disregarding radiation pressure could lead to significant underestimations or overestimations of M • as determined by the virial method (see e.g., Marconi et al. 2008;Netzer & Marziani 2010).Although Keplerian motion was found to dominate the BLR velocity field in the AGN discussed in Section 1.2, there may also be weak radial flows.In fact, there is evidence for inflows in NGC 3783 (Amorim et al. 2021a).Furthermore, these AGN have low to moderate Eddington ratios (Γ < 0.2, based on the reported L AGN and M • ), but outflows are expected to be much stronger in high and super-Eddington AGN.The outflows may be associated with a disk wind, with the relative prominence of the disk and wind components scaling with AGN luminosity (Elitzur et al. 2014).This suggests that the dominant emission line producing structures in the BLR are dependent on the SMBH's accretion rate.
In our adopted model (Equation 23), there is an outflow when ΓF M > 1 where F M ∝ 1/N col , therefore the conditions for an outflow in our models depend on the Eddington ratio, Γ, and column density.Since Γ is dependent on L AGN and M • , it will be an output when BELMAC is applied to model real data, but here we have fixed the AGN values to Γ = 0.07 and M • = 10 8 M ⊙ .To achieve an outflow, the clouds must have N col < 5.26 × 10 22 cm −2 (see Equation 23).For models with s = 0, clouds have clouds a constant cross-sectional area and N col (r) ∼ 10 22 cm −2 , thus ΓF M > 1 throughout the entire BLR.When cloud cross-sectional areas vary with distance (s ̸ = 0), the clouds with N col (r) > 5.26 × 10 22 cm −2 are too massive for radiative forces to overcome gravity.As N col (r) decreases the F M increases, because the clouds expand to maintain constant mass.
The above density calculations are also dependent on the initially-defined gas density, n(R d ).If n(R d ) is large, e.g., n(R d ) = 10 10 cm −3 , ΓF M < 1 throughout the BLR and v rad (r) → GM • /r.Furthermore, in these models we have assumed that the mass of the clouds is conserved (i.e., m cl (r) = constant).In real BLRs, the relation between the sizes of clouds and their masses may vary, leading to a more complex gas density distribution.Varying the cloud mass would alter our adopted radial cloud motion model and could produce less-boxy line profiles.The last parameter to consider is α f , the fraction of L AGN absorbed by a cloud.A more realistic BLR model would let α f depend on the ionization parameter and N col , but here we have opted to keep that quantity constant (α f = 0.5; Equation 22).Future versions of BELMAC will include the radial and time dependence, α f (r, t ′ ).
The line profiles produced by the spherical shell-radial flow models tend to be boxy unlike the observed broad line profiles in most AGN, which typically resemble Gaussian or Lorentzian profile.This is a feature of the radial velocity field model adopted here.As can be seen in Figure 5, the radial velocity curves rapidly level off to a terminal velocity, within ∼ 0.1−0.3R d .As v(r) = constant in the outer regions, the dispersion in LOS velocity varies little with radius, resulting in flat topped profiles.Although, varying α f (r, t ′ ) could reduce the acceleration so that the terminal velocity is reduced at a larger radius.Here, the velocity range for the large majority of clouds in the ensemble is quite small, resulting in boxy line profiles that don't resemble most observed broad emission lines.For s = 0 and p = 0, the clouds experiencing the greatest acceleration are also matter-bounded clouds, which don't contribute to the reverberation response, but still contribute to the average line profile.Additionally, allowing the cloud mass to vary, such as varying cloud mass, could produce less boxy profiles.Finally, a BLR model with combined velocity fields, with non-radial components contributing to the LOS the velocity (such as turbulence and/or Keplerian motion) will tend to produce more Gaussian or Lorentz-like line profiles that would be more representative of observed BELs.

Representation of Real Broad Emission Line Responses
Utilizing simple recombination theory (Equation 10) for the cloud luminosity is a reasonably good approximation for Hα emission.Roughly speaking, the Hα flux emitted per unit area by a given cloud is proportional to the column density of ionized gas, N s (r).This in turn depends on U (r, t ′ ) and therefore responds to a changing ionizing radiation field.We can also make general statements about the emissivity distribution for other broad lines.For instance, the high-ionization emission lines, such as CIVλλ 1551, 1548 Å, are mainly emitted near R in where the ionizing radiation field is strongest (Osterbrock & Mathews 1986).The response of CIV would peak at a short time delay, then rapidly drop off with time delay.Therefore, we'd expect the response behavior of high-ionization lines like CIV to be somewhat similar to models with small values of s and p. Conversely, low-ionization lines, such as MgIIλλ 2795, 2803 Å, have an emissivity distribution that varies more slowly with radius and so can be approximated by specifying larger s and p values.The power-law approximation of ε(r) includes the s and p parameters and should provide a general guide to the dependence of a given line's reverberation response to U (r, t ′ ), n(r), and N col (r).However, the emissivity distributions for true emission lines in BLRs are certainly more complex functions than implied by Equation 12.For a more sophisticated calculation of the line emission fluxes, the next version of BELMAC will include photoionization model grids for selected BELs.
When the elf is used, Equation 11, it provides a reasonable approximation for anisotropic line emission from a distribution of clouds and the ACE models should be generally representative of emission lines that are optically thick.For example, ACE models can be utilized to approximate the response of Lyα that is emitted in the ionized part of the cloud and escapes almost entirely through the illuminated face, but has a very low probability of escaping through the non-illuminated face.ICE models represent lines that are optically thin, such as CIV, which can escape isotropically if the cloud contains no dust.Furthermore, we assume that a cloud which becomes matter-bounded also becomes effectively optically thin, allowing the line emission to escape isotropically.We could expect dust to survive in the partially ionized and neutral interiors of BLR clouds (Baskin & Laor 2018).Dust is not explicitly included in the current version of BELMAC, but ACE, modeled with the elf parameter, could also be used to approximate the effects of dust extinction in BLR clouds.

Inferring Black Hole Mass from the Reverberation Lag
As discussed in Section 1 an important motivation for BLR reverberation mapping is to determine R BLR and hence M • via the virial theorem (Equation 1).The BLR models presented Section 3 show that the RWD is not always equivalent to the LWR.Matter-bounded and anisotropically emitting clouds have a significant affect on the response-weighted lag, t ′ RW , which could lead to overestimates of the SMBH mass.If we assume the velocity dispersion is determined exactly and assume ct where M • is the true black hole mass and M •, vir is the mass inferred from reverberation mapping.Since we are assuming 2R d LWR = R BLR , if LWR/RWD < 1 then RWD is overestimating LWR and therefore also M • .This is the case for all of the models with matter-bounded clouds and/or ACE presented in Figures 10 and 14.
Figure 18 shows the LWRs and RWDs for the spherical and disk BLR models that include matterbounded clouds from Figures 10 and 14.The bottom panels show the RWD can be ∼ 1.1× LWR, for the smallest BLRs, to nearly 2× LWR for Y BLR ≳ 20.ACE also results in LWR/RWD < 1, even when 100% of the clouds are radiation-bounded.Figures 10 and 14 show that RWD ≈ 1.2−1.5×LWRfor ACE models with no matter-bounded clouds (s = −1, −2).As discussed above, ACE is a fair representation of optically thick emission lines and the possible extinction by dust embedded in the cooler regions of the clouds.Therefore, in the sample of BLR models explored here, R BLR is correctly determined from t RW only when all the clouds in the BLR ensemble are radiation-bounded and emitting isotropically.Given the conditions in real BLR, this scenario is unlikely to occur for Balmer lines, which tend to be optically thick (Osterbrock & Mathews 1986).This suggests the lag measured from cross-correlation does not relate to R BLR as straightforwardly as previously assumed.However, we will be revisiting these models with a version of BELMAC that implements photoionization models in forthcoming work.

Comparison to TORMAC
Although TORMAC was designed to simulate reverberation of the torus dust emission in AGN, it is worth discussing similarities and differences between TORMAC and BELMAC.
TORMAC includes physics potentially important for both the torus and the BLR; cloud occultation, cloud shadowing, and the effect of cloud orientation.A cloud is occulted and its emission is attenuated when another cloud is positioned along the LOS between the former cloud and the observer.The strength of the attenuation depends on the number of occulting clouds and their optical depths.Cloud shadowing occurs when a cloud is positioned within the shadow of another cloud intervening between it and the illuminating source.A cloud without direct exposure to the source is indirectly heated by the diffuse emission from neighboring clouds and emits isotropically (Almeyda et al. 2017(Almeyda et al. , 2020)).As a general point, we expect cloud occultation and shadowing to be less important for the BLR's gas clouds, even if the clouds themselves are optically thick for a given line, because the BLR as a whole is presumed to be optically thin (i.e., C f < 1).Therefore, these two effects are not included in the current version of BELMAC, but could be incorporated in later versions.Cloud emission anisotropy, on the other hand, is expected to be important in some situations and is incorporated in BELMAC via the elf parameter (Equation 11).This is based on the orientation-dependent fraction of the cloud's illuminated surface that is viewed by the observer.TORMAC includes a more sophisticated implementation of cloud orientation effects than BELMAC's elf parameter.The anisotropic emission of a cloud in TORMAC, adopted from the clumpy torus model by Nenkova et al. (2008a), is wavelength dependent.In this model, the spectrum of a cloud is approximated by averaging over the emission from many slabs of gas at varying orientation angles relative to the central illuminating source (Almeyda et al. 2017(Almeyda et al. , 2020)).
TORMAC also includes anisotropic illumination of the torus.Edge darkening of the AGN's accretion disk causes the illuminating radiation field to be anisotropic.As a result, the dust sublimation radius is not constant, but a function of polar angle and the dust can exist closer to the source near the midplane than towards the pole.For dust clouds arranged in a sphere and anisotropically illuminated, the cross-section of the inner dust free cavity would have a figure-8 shape (Almeyda et al. 2017;Baskin & Laor 2018;Kudoh et al. 2023).Although we do not explore anisotropic illumination here, it is included in BELMAC.For the BLR, this effect will result in both the ionization parameter (Equation 9) and the outer boundary, the dust sublimation radius, becoming functions of polar angle, θ.

Other Reverberation Mapping Modeling Codes
Mangham et al. ( 2019) categorizes 2 methods for deciphering the structure and kinematics of the BLR from velocity-resolved reverberation mapping data.The 'inverse' method constructs the 2DRF from data.For example, the Maximum Entropy Method (MEM) is a deconvolution technique utilized by the echo mapping code, Memecho, to recover the velocity-resolved transfer function from observations and fit reverberation mapping data (Horne 1994).However, it is still necessary to interpret the results to determine the geometry and dynamics.
The other method is forward modeling, which is the approach used by BELMAC, Caramel (Code for AGN Reverberation and Modeling of Emission Lines) by Pancoast et al. (2011), and BRAINS by Li et al. (2013).BELMAC and Caramel both create a 3D ensemble of clouds, however, Caramel utilizes a more prescriptive parameterization of BLR properties.For example, Caramel treats the clouds as particles that are representative of a gas density field and approximates the responding emission by setting the emissivity as a radial power-law function (Williams & Treu 2022).In contrast, BELMAC models the BLR as an ensemble of discrete clouds, whose line emission is calculated using nebular physics.We used hydrogen recombination theory for the models presented here, but a large grid of photoionization models will be used in the forthcoming version.Caramel includes a parameter estimation procedure to recover the best matching model for a time series dataset.This capability is currently not implemented in BELMAC, but in the future we will include a parameter estimation procedure.This will make BELMAC a valuable tool for BLR reverberation mapping modeling that will compliment existing modeling codes such as Caramel, while offering different capabilities.
The BLR Reverberation-mapping Analysis In AGNs with Nested Sampling (BRAINS) code, developed by Li et al. (2013), is based on Caramel and employs similar methodology.However, BRAINS allows the BLR emission to respond non-linearly and additionally, Li et al. (2018) added a 2-zone BLR geometry, where each zone may have a different cloud distribution, which could be used to model a disk-like BLR with a wind.BELMAC may also be used to model a 2-zone BLR, with different gas density and cloud distributions, velocity fields, and radial sizes, but we leave this discussion for a future paper.Lastly, BRAINS also models the velocity-resolved spectroastrometric response (Li & Wang 2023, and references therein).This is not currently a feature in BELMAC, but could be implemented with straightforward modifications.

CONCLUSION
We have presented a parameter exploration of velocity-resolved, reverberation response functions (numerical approximations of the transfer function) of the BLR using BELMAC.The main purpose of BELMAC is to serve as a versatile forward modeling tool for analyzing BLR reverberation mapping data.BELMAC is designed for flexibility, allowing various different combinations of geometry and dynamical models to be incorporated.Here, we focused on 2 distinct BLR models; a spherical shell with a radiation-pressure driven outflow and a thin disk with Keplerian motion.In the models presented in this paper, the broad Hα line emission is calculated using an analytical prescription based on recombination theory.The main results of this paper a summarized below.
1.The τ − v || relationship is linear for spherical radial outflows, and elliptical for rotating disks.
2. For models with isotropic cloud emission (ICE) and no matter-bounded clouds, the 1DRFs are flat-topped for spherical shell BLRs and double-peaked for disk-like BLRs, for delays less than the light-crossing time of the inner cavity (τ ≤ 1/Y BLR ).For τ > 1/Y BLR , the 1DRFs for both shells and disks decay as τ increases, but the shape of the decay depends on the gas density profile and how the clouds are distributed.Whether for a sphere or a disk, the 1DRFs become more extended as s decreases and p increases (ε V (r) decreases more slowly with r).
Additionally for the disk models, the time delay difference between the double-peaks will reduce with decreasing inclination.
3. For models with ICE and no matter-bounded clouds, the time-averaged line profiles are flattopped for spherical BLR models with radial outflows and double-peaked for rotating disk models.In spherical, radial outflow models, the response begins in the blue-wing of the line profile then propagates to the red-wing.Therefore, the RMS profiles tend to be skewed blueward, but the time-averaged line profiles are symmetric.In the rotating disk models, the response begins in the profile core, then propagates to the maximum velocity, ±v || Kep , in the wings, and finally moves symmetrically back towards the core.The averaged line profiles have symmetrical double-peaks, which have separations dependent on the disk's inclination.The RMS profiles also have symmetrical double-peaks, however, at large inclinations, the response amplitude at v || Kep ≈ 0 km s −1 and τ ≈ 0 is strong enough that the RMS spectra are triple-peaked.4. We assume that recombination theory provides a reasonable estimation of the broad Hα response.Considering lines other than the hydrogen Balmer Series, the models with ICE can be considered to represent emission lines that are likely to be optically thin, such as CIV, and likewise, the ACE models represent lines, such as Lyα, that are likely to be optically thick.
5. Both matter-bounded clouds and ACE can affect the response quite dramatically because these effects suppress the response of the inner clouds, or the near-side clouds, respectively.Both effects cause the 1DRF to become more extended in time delay.The models with matterbounded clouds respond as if the BLR has an effective inner radius that is greater than the true value of R in .In ACE models, the 1DRFs tend to peak at delays longer than 1/Y BLR .
In radial outflow models, ACE causes redward asymmetric line profiles.For rotating disks, the line and RMS profiles decrease in amplitude by half, but ACE doesn't affect the shapes.However, matter-bounded clouds emit isotropically and therefore partly mitigate the profile asymmetry and decreased amplitude.
6.The luminosity-weighted radius is equivalent to the response-weighted delay when the BLR is comprised entirely of radiation-bounded, isotropically emitting clouds.However, the RWD can significantly exceed the LWR when matter-bounded clouds are present or when clouds emit anisotropically.These effects in combination lead to over-estimations of R BLR by ∼ 10 − 100% implying that measured cross-correlation lags may also over-estimate R BLR in such circumstances.
We will use photoionization models grids to explore multi-emission line responses in a forthcoming paper.More complex combinations of geometries and velocity fields, such as a rotating disk with a biconical wind, will also be presented in a later paper.In the future, BELMAC will include a parameter estimation procedure for matching response models to BLR reverberation mapping data.BELMAC will compliment existing BLR reverberation modeling codes such as Caramel, Memecho, and BRAINS.This paper is based upon work supported by the National Science Foundation under AARG/Grant Number 2009508.

Figure 1 .
Figure1.A summary of BELMAC's operation.The parameters are explained in Section 2 and listed in Table1.The parameters in squares are fixed and provided by the user.The free parameters in thick-lined circles determine the BLR model set-up and parameters in thin-lined circles are calculated within the code.The final output of the code is the velocity-resolved response, L(t, v || ).

Figure 2 .
Figure 2. Geometrical parameters used for clouds arranged in a flared disk in BELMAC.R in and R d are the inner and outer radii, σ is half-angle width of the disk, and i is the inclination of the disk axis to the observer's line of sight.On the left is a sharp edged BLR configuration and a 'fuzzy' (Gaussian distribution in β) boundary is shown on the right.Adapted from Nenkova et al. (2008b).

Figure 3 .
Figure 3. Two geometrical BLR configurations explored in this paper, a sphere and thin disk (i = 45 o , σ = 15 o ).The color represents the elf parameter.The dark green means the observer views only the non-illuminated face of the cloud.The white means the observer sees only the fully illuminated side of the cloud.Intermediate shades indicate the fraction of the illuminated side see by the observer.

Figure 4 .
Figure 4.The analytical LWR (top row) and ratio RWD/LWR (bottom row) for a spherical shell for varying combinations of s and p parameters with Y BLR .
's and Netzer & Marziani (2010)'s prescription for clouds accelerated by radiation pressure due to Thomson scattering and absorption of ionizing radiation emanating from the accretion disk.

Figure 6 .
Figure 6.Left: Cross-section of a spherical BLR for Y BLR = 20 with isodelay surfaces over plotted for several values of the normalized time delay, τ .The arrows represent the velocity field, with the colors blue, red, and yellow indicate a blueshift, redshift, and no shift, respectively.An animation of this figure can be viewed in the online version of this paper.In the animation, the isoddelay surfaces appear with τ from 0 to 1, which is denoted in the upper left-hand corner.In the static figure, all the isodelay surfaces are shown and each isodelay curve is labeled with it's corresponding value τ along the border.The bold isodelay curve corresponds to the light-crossing time of R in , τ = 0.05.Right: Response function for a spherical BLR with a radial outflow, for a gas density index s = −1 and cloud distribution index p = 0.All of the clouds in this BLR model are radiation-bounded.The color-scale represents the normalized response amplitude.The lower panel shows the velocity averaged, 1DRF in black, with the corresponding analytical model, calculated using Equation 14 with η = −10/3, in dashed yellow.In the animated version, the 1DRF is plotted in-sync with the isodelay surfaces in the left figure.The dashed lines in the 2D and 1DRFs indicate the light-crossing time of R in , τ = 1/Y BLR = 0.05, and the RWD = 0.13.The right panel shows the time averaged line profile (black) and the RMS profile (blue).The animated version of this figure also plots the line profile changing with τ (gray), which is in-sync with the bottom panel and isodelay surfaces on the left.The profiles are each normalized to their respective maximum amplitude and binned to 100 km s −1 per bin.

Figure 7 .
Figure7.Velocity-resolved response functions (2DRFs) of spherical BLRs with radial outflow for varying s and p model parameters.The velocities are normalized to the maximum velocity of the respective BLR models.For instance, the response for s = 0, p = 0 is normalized to 17 × 10 3 km s −1 and for the s = −2, p = 2 model, to 7 × 10 3 km s −1 .Refer to Figure5for the maximum velocities in other cases.

Figure 8 .
Figure 8. Left: The velocity integrated, 1DRFs for the models shown in Figure 7 for ICE (blue) and ACE (red) models.The dashed blue and red lines indicate the RWDs for the ICE and ACE models, respectively.The values of Y apt BLR are given for each p when the ensemble includes matter-bounded clouds (see Section 3.1.1for details).Right: The delay time integrated, line profiles of Figure 7.The profiles of ACE models (red) are normalized to the peak of the ICE model profiles.For each model, the velocities are normalized to the maximum velocity of the broadest line, which is when s = 0 and p = 0.

Figure 9 .
Figure 9. Same as right panel in Figure 6, but with a gas density index s = 0 and therefore the ensemble includes matter-bounded clouds.The green and blue lines are the analytical 1DRFs for Y mbc BLR = 2.60 and Y apt BLR = 2.27, respectively.

Figure 10 .
Figure 10.The LWR and RWD of the spherical ICE (left) and ACE (right) BLR models of varying Y BLR for each combination of s and p parameters, with their ratio in the bottom panel.

Figure 11 .
Figure 11.Same as the right panel in Figure 6, but for a disk-like (σ = 15 o ) BLR with clouds moving in Keplerian orbits.The disk is viewed edge-on (i = 90 o ), s = −2, and p = 0.The dashed lines indicate the light-crossing time of R in and the measured RWD = 0.197.The white ellipses are defined by Equation 21.The system is entirely radiation-bounded clouds.This figure has an animated version that is available in the online version of this paper, which does not display the ellipses.In the animation, the 1DRF is plotted with τ from 0 to 1, along with the τ −varying line profile (gray).The changing line profile is a narrow, bright peak at τ = 0, then quickly becomes lower in amplitude and broader with double peaks at the extremities of the profile's wings.These low-amplitude peaks slowly move back towards the center of the profile until τ = 1.The τ −averaged line profile (black), RMS (blue), and, in the animation, the τ −varying line profiles (gray) are normalized to their respective maxima.

Figure 12 .
Figure12.The same as Figure11, including an animated version that is available online, but for an ACE model.In the animation, the τ −varying line profile is a narrow, low-amplitude peak at τ = 0, then quickly becomes brighter in amplitude and broader with double peaks at the extremities of the profile's wings.These high-amplitude peaks slowly move back towards the center of the profile until τ = 1.

Figure 13 .
Figure 13.The same as Figure 8, but for disk-like (σ = 15 o ) BLRs with clouds moving in Keplerian orbits and viewed from i = 60 o .

Figure 14 .
Figure 14.The same Figure 10, but a rotating disk with i = 60 o .

Figure 15 .
Figure 15.Isodelay surfaces over plotted for a disk-like, σ = 15 o , BLR viewed from left to right at, i = 0 o , i = 45 o , and i = 90 o .The isodelay surfaces are labelled with the corresponding delays, τ , across the top and down the right-hand side of each sub-figure.When s = 0, the light gray points are matter-bounded clouds, that don't respond.The light blue points represent clouds become matter-bounded when the continuum pulse reaches them and the gray clouds are always radiation-bounded.If s ̸ = 0, then all of the clouds in the ensemble are radiation-bounded.The dashed line in the i = 0 o panel is the isodelay surface corresponding to inflection point seen in the response in the left panel of Figure 16.The thicker, black line is the light-crossing time for R in .

Figure 16 .
Figure 16.1DTFs (left) and line profiles (right) for a disk-like (σ = 15 o ) BLR with clouds moving in Keplerian orbits and viewed from face on to edge-on, i = 0 o , 45 o , 90 o , for each row and increasing p = 0, 1, 2 for columns left to right, and s = −2.

Figure 18 .
Figure 18.The LWR, RWD, and LWR/RWD for (left) spherical and (right) disk-like (σ = 15 o , i = 60 o ) BLR models, but for only models containing matter-bounded clouds (s = 0).The blue and red lines represent ICE and ACE BLR models, respectively.

Table 1 .
BLR geometry and property parameters used to set-up for the ensemble and their standard model value.

Table 4 .
The same as Table3, but for a disk geometry with i = 60 o .