BLR size in Realistic FRADO Model: The role of shielding effect

The effective size of Broad Line Region (BLR), so-called the BLR radius, in galaxies with active galactic nuclei (AGN) scales with the source luminosity. Therefore by determining this location either observationally through reverberation mapping or theoretically, one can use AGNs as an interesting laboratory to test cosmological models. In this article we focus on the theoretical side of BLR based on the Failed Radiatively Accelerated Dusty Outflow (FRADO) model. By simulating the dynamics of matter in BLR through a realistic model of radiation of accretion disk (AD) including the shielding effect, as well as incorporating the proper values of dust opacities, we investigate how the radial extension and geometrical height of the BLR depends on the Eddington ratio [and blackhole mass], and modeling of shielding effect. We show that assuming a range of Eddington ratios and shielding we are able to explain the measured time-delays in a sample of reverberation-measured AGNs.


INTRODUCTION
The delayed response of BLR to the variation of AGNs continuum is a probe to study the structure of BLR (Peterson et al., 1998). The main discovery is a linear log − log scaling law between the size of BLR and the AGN monochromatic luminosity, so-called the radius-luminosity (RL) relation. It has been observationally discovered (Kaspi et al., 2000;Peterson et al., 2004;Bentz et al., 2009Bentz et al., , 2013, and theoretically explained by Czerny and Hryniewicz (2011) who proposed that failed dust-driven wind (FRADO) is responsible for the formation of BLR. However, the departures from the law seen for some sources with high Eddington ratios (Du et al., 2015(Du et al., , 2018Grier et al., 2017), obviously requires the structure and extension of BLR to be studied in more detail, not only in observational side via reverberation-mapping (e.g. Panda et al., 2019), but also through theoretically-motivated models. Provided that the RL relation is concretely justified for all already-observed sources Rakshit et al., 2019), AGNs can be treated as astronomical candles to be used in cosmology (Watson et al., 2011;Haas et al., 2011).
In this article, we briefly report our pilot numerical study on the size of BLR based on the FRADO model. A concise background on the existing BLR models is addressed in section 2. In section 3 an introduction to FRADO model is provided, followed by section 4 on shielding effect. After a short description of our numerical setup in section 5, the method we adopted to find the BLR size is specified in section 6. We discuss our results and conclude the article in the last two sections.
The first model in the category of Dust-driven Winds is FRADO model (Czerny and Hryniewicz, 2011;Czerny et al., 2015Czerny et al., , 2017. There is only another existing model in this category for the formation of BLR in which the static atmosphere of disk is puffed up by irradiation in the dust-dominated region (Baskin and Laor, 2018). Differently from this static model, FRADO gives a dynamic view of the BLR. The advantage of the FRADO model is that the underlying physics is known, the approximate compensation of the inflow/outflow leads to approximately symmetric lines, and the model is likely to produce also single component line instead of disky profiles if the vertical motion is strong, as expected at higher Eddington ratios. However, the model applies only to Low Ionization Lines, and for High Ionization Lines a line-driven outflow is likely favored.

FRADO MODEL
FRADO model is one of the interesting theoretically-based models explaining the formation mechanism of the BLR in AGNs. It is based on the existence of dust at large radii where the AD atmosphere is cold enough (Dong et al., 2008). The dusty material, i.e. clumps composed of dust and gas, is lifted up by the AD local radiation flux. Once reaching considerably high altitudes above the disk surface, dust content of clump gets evaporated by the strong radiation coming from the central parts of an AGN and then follows a subsequent ballistic motion and forms a failed wind.
In the basic 1-D analytical form of FRADO model, only the vertical component of local flux of accretion disk at a given radius was considered. It directly gives the inner/outer radius of BLR. Asymmetry in line profiles can be seen due to the dust evaporation; line shape parameter FWHM/σ is consistent with data; and also the model predicts the dependence of virial factor on blackhole mass (Czerny et al., 2015(Czerny et al., , 2016(Czerny et al., , 2017. However, the analytical local flux of Shakura-Sunyaev (S-S) disk (Shakura and Sunyaev, 1973) that is used in this model is constant with height; radiation from other radii and also the radial motion of material are neglected, as well as the inner boundary condition in S-S AD.
In this paper, for the first time we calculate the realistic 3-D radiation pressure of AD in which the radiation coming from other radii is also contributing. The real dust opacities are also taken into account, and the inner boundary condition of S-S AD is included. The general form of the 3-D equation of motion, considering the forces acting on the clump to be the gravity of central blackhole and radiative force of AD, is thus given as below f (r, ψ, M BH ,ṁ, κ ext (λ), T sub , ρ, ϕ, λ, C) da dλ (1) Figure 1. Illustration of shielding models. The disk is shown in light grey in both models, while the visible area of the disk (not shielded) to the clump is recolored to dark grey. In Patch Model (a), the radius (r p ) of the patch which is not drawn in scale is linearly proportional to the height (H) of the clump from the AD surface. In Wall Model (b), the radius (r w ) of the cylindrical shell is linearly proportional to the height (h w ) of the wall while β = tan(θ) is constant. The visible area of the disk to the clump varies due to the changes in the size of patch and wall while the clump is moving. In order to keep the figures neat and clean, the conceptual location of BLR in , and also the torus [BLR out to R out ] are not drawn.
where G is the gravitational constant, M BH the blackhole mass, ψ the dust-to-gas mass ratio, κ ext the dust opacity as a function of wavelength,ṁ the dimensionless accretion rate, T sub the dust sublimation temperature, and r = RR + Hẑ is the vector showing the position of the clump in cylindrical coordinates in which the blackhole is at the origin as shown in the figure 1(a). The AD for simplicity is mapped onto the equatorial plane (z = 0) from which the clump's height is measured. The parameters of ρ and ϕ show the coordinates of infinitesimal surface areas of AD in polar coordinates, and C stands for some fixed factors and physical constants. The vector radiative force acting on a clump is calculated through a double integral evaluation over AD surface (visible area due to shielding effect) and an integral over the effective range of wavelengths for the adopted dust model of the clump.
The effect of the dust evaporation is also included; assuming dust instantly emits the absorbed energy in the form of blackbody radiation, we numerically calculate the irradiated energy absorbed by dust Q abs at each integration step. Once the Q abs for a dusty clump reaching a certain altitude above AD exceeds the level of energy emission corresponding to sublimation temperature of dust Q emit (T sub ), the clump loses its dust content. The radiation pressure force is then immediately switched off and the dustless cloud follows a ballistic motion, finally falling back onto the AD surface. The subsequent motion of ionized dustless gas driven by Thompson electron scattering is neglected.

SHIELDING EFFECT
It is well-known from all radiatively-driven wind models that efficient alleviation of the material above the disk is not possible if the wind launching region is irradiated (see e.g. Proga, 2007). Some shielding is necessary, and it can be due to either inner failed Compton-driven wind or cold material forming at the edge between the inner hot ADAF (Advection Dominated Accretion Flow) and a cold disk. This shielding effect protects dusty material from too intense radiation from the central parts of AD, including X-rays. Therefore it would make the local radiation pressure from the disk to be effective in launching the material upward avoiding too early dust evaporation in the cloud. In order to mimic this effect we model it in the following two configurations: In this model it is assumed that the radiation only comes from a small patch of the disk as shown in figure 1(a). The radius of patch is proportional to the height of the clump from the disk surface. Mathematically speaking, this dynamically variable radius can be expressed as a function of clump's height as r p = αH where α is a free constant parameter of the model. The location of patch passively moves with the clump and increases/decreases in size as the cloud goes up/comes down.

• Wall Model
In this model as shown in figure 1(b), as long as the clump is at rest on the AD surface, the radiation due to inner region of AD extending from inner radius of AD (R in ) up to half of the radial position of the clump (R/2) is shielded. The model can be illustrated as a circular wall barrier (cylindrical shell) surrounding the inner region of the AD whose radius and height are r w and h w = β r w , respectively; where the model constant parameter β = tan(θ). The range of variation for r w is min(r w ) = R in and max(r w ) = R out /2. Moreover, until the clump to be launched, the wall also obscures the radiation of almost half of the AD located behind the wall (viewed from the clump). As the clump flies up, due to an anti-correlation between the wall's radius and the clump's height, the wall shifts inwards while β is constant. As a result the clump sees a larger portion of AD irradiating. The moments when the clump is flying in an altitude H > βR above AD, it can be irradiated by the whole AD as if there is no shielding effect.
Obviously, the proposed models of shielding are geometrically different. The results from these models can provide us with useful information about how shielding effect really looks like, and how effective it is in the physics of BLR and its dynamics. We incorporate these shielding models into our computations by imposing proper dynamic upper-lower integral limits while evaluating the radiative force produced by AD.

NUMERICAL SETUP
In order to perform the simulation, we consider the following setup. The blackhole mass is set to be of 10 8 M (Panda et al., 2018) which is almost the mean value of the blackhole mass in quasar catalog of Shen et al. (2011). The underlying disk is the S-S disk model which radiates locally as a blackbody. We consider a range of AD radii from R in = 6R g to R out = 10 6 R g where R g is the gravitational radius of the blackhole defined as R g = GM BH /c 2 . The adopted dimensionless accretion rates (ṁ =Ṁ /Ṁ Edd whereṀ Edd is the Eddington accretion rate) are 0.01, 0.1, and 1. Actual values of wavelength-dependent dust opacity are incorporated into our computation. These values are computed using MCDRT (Szczerba et al., 1997) and KOSMA-τ PDR (Röllig et al., 2013) codes which give the mean extinction cross-section per dust mass (equivalent to opacity) for different dust models with a given distribution of dust grain size and different materials. We employ the MRN dust model (Mathis et al., 1977) which consists of silicate and graphite and has an identical grain size distribution for both dust populations. Silicate features are seen in the AGN spectra (e.g. Netzer et al., 2007); graphite is more problematic, we would rather expect amorphous grains at the basis of the UV spectra (see e.g. Czerny et al., 2004), but we use MRN for simplicity. Moreover, having assumed the gas and dust within the clump being strongly coupled we set the parameter ψ equal to 0.005 which is the mean value for Milky Way (Mathis et al., 1977). However, an estimate mean ratio of 0.008 is also adopted by others for a sample of AGNs (Shangguan et al., 2018). It is also assumed that the dust content of clumps sublimates at a temperature of 1500 K which represents approximately the mean value for the different grain species and sizes (see Baskin and Laor, 2018). The clumps in our model are launched with an azimuthal Keplerian velocity and zero vertical velocity from the surface of AD which has a thickness, D(R). The latter, which is function of radius, was computed separately for different initial values of AD size, accretion rate, black hole mass, and viscosity based on Rosseland mean opacity (Czerny et al., 2016;Różańska et al., 1999) without the effect of self-gravity of AD itself being included. This implicitly also means that the height of the clump (H) never vanishes since our mathematically thin disk is at z = 0. The Runge-Kutta method is then applied in order to integrate the equations of motion.

BLR SHAPE IN FRADO
In order to find the distribution of material above AD we need to determine the radial extension and height of BLR as follows.
The radial extension of BLR is by definition the range confined within BLR in to BLR out . Theoretically, BLR in is the onset of BLR where the dust is instantly sublimated upon departure from AD. The outer radius of BLR, i.e. BLR out , is where the dust survives irradiation by the whole disk even in spherically symmetric approximation of the radiation field. We determine this outer radius by finding the position where the sublimation height for a given radius is equal to the radius itself. Obviously, the shielding effect does not play a role at this large height. In order to find these two radial ends of BLR, for each different initial conditions we introduce a dense grid of (R, H) pairs ranging as R : [R in to R out ], and H : [D(R) to R out ]. We then numerically compute Q abs at each grid point. The geometrical location where Q abs = Q emit (T sub ) is a curved line called sublimation location, S(R), dividing the region above AD into two parts, above which dust can not survive irradiation. The crossing radius of S(R) and D(R) yields the BLR in ; and BLR out can be found by setting S(R) = R.
As for the geometrical height of BLR, we observe the trajectory of clumps being launched from a set of logarithmically spaced initial launching radii (R init ) ranging from BLR in to BLR out . Neglecting the relaxation due to probable encounters between clumps and loss/gain of momentum, we find the maximum heights (H peak ) the clump can attain per each initial launching radius within the set. However, in order to stress on the importance of the opening angle of the BLR clouds distribution, introducing the variable (P = H peak /R peak ) where R peak is the maximum height's corresponding radius, we find the maximum of P , i.e. P peak . Therefore, the shape of BLR is specified by [BLR in − BLR out ], and [D(R) − H peak (R)]. The R(P peak ) and H(P peak ) corresponding to the unique value of P peak represents the effective location of BLR clouds. This is based on the assumption that the the parameter R(P peak ) marks the region effectively exposed to the irradiation by the central parts of AD.

RESULTS & DISCUSSION
Having performed numerical computations for different initial conditions and for both models of shielding, we have obtained values for the size of BLR as collected in the table 1. Numerically computed values of the monochromatic luminosity of AD at 5100Å in erg/s depend only on the accretion rate and blackhole mass; the effects of the shielding and AD self-irradiation are not included (Czerny et al., 2003). Other columns of data listed in the table 1 are expressed in light-days, exceptṁ and P peak being dimensionless. The viewing angle for calculating the time-delay τ (equivalent to BLR radius assuming the light speed c propagation of the signals) is set to be i = 39.2 (Lawrence and Elvis, 2010).
The sublimation surface, S(R), obtained for all initial conditions and parameters of shielding models in table 1 is shown in figure 2. Upper and lower panels show the same plots but in linear and log scale, respectively; in the log plots the onset of BLR and disk thickness is clearly visible while the linear plots Frontiers Table 1. Output data for the size of BLR for blackhole mass of 10 8 M . Values of L 5100 are in ergs/s. Other data are in light-days exceptṁ and P peak being dimensionless. Patch model and Wall model are specified by their characteristic parameters of α and β, respectively.
Shielding Model (i = 39.2) m BLR in BLR out P peak R(P peak ) log(L 5100 ) log(τ in ) log ( (1) Reaching the peak with dust-content sublimated: dust-less failed wind (2) It never forms a failed wind neither dust-less nor dusty. All clumps just escape to torus (3) This row is just added for the purpose of comparison.
better illustrate the radial extension of BLR for different values of the accretion rate. Obviously, the Patch model provides a stronger shielding than Wall model, and the table is also sorted down based on the strength of shielding. However, it seems that the presented shielding models can behave similarly if certain values of α and β are chosen. For example, based on the figure 2 one can expect the two shielding models to show very similar patterns for α = 1/β. S(R) profiles in figure 2 for all values of α. and β show a convergence to the case of no-shielding at BLR out . This is related to the definition of the BLR out , with the requested condition set at large height of the cloud above the disk where the shielding becomes unimportant.
In addition the wall model interestingly shows a privilege of unique prediction of BLR in for all values of β that implies the Wall shielding does not work efficiently when the cloud height is small as it happens close to BLR in . However, testing with other values of β is required to firmly determine the model properties.
The two considered shielding models have significantly different geometrical properties, as shown in figure 1. In the Wall model a cloud is exposed to the radiation from a broad range of radii of AD while in the Patch model it is affected only by a local radiation. This may affect the cloud dynamics, and likely will  Figure 2. Representation of the sublimation location above AD equatorial plane for Eddington rates of 0.01(left panels), 0.1 (middle panels), 1 (right panels) in normal (upper panels) and logarithmic (lower panels) scales. Axes are normalized to BLR in from analytical FRADO model for the corresponding Eddington rate to be easily comparable. The area shaded in light-green covers the radial extension of BLR predicted by analytical FRADO model. Yellowish strip shows the inner part of torus. The onset of BLR and the thickness of AD are clearly visible in logarithmic panels, while in the linear ones the BLR radial extension for different accretion rates is easily comparable.
be reflected in the emission-lines shape which is beyond the scope of the current work. However, further studies including the cases with α = 1/β are required to highlight the differences between these shielding models properly, and to enable us to make a choice between these two.
As can be seen in figure 2, the ratio of (BLR out /BLR in ) increases with Eddington ratio for both analytical and numerical FRADO. However, the BLR size in numerical models with shielding is radially more extended than that of analytical ones which implies longer time-delays from the entire BLR.
According to the table 1, in case of no-shielding for accretion rates of 0.1 and 1, all launched clumps only escape to the torus and do not form a failed wind. These cases were omitted due to being out of interest in the current preliminary study. There are cases for accretion rates of 0.01 and 0.1 in Patch model in the table 1 with R peak ∼ BLR in yielding time-delays for BLR in , log(τ in ), being larger than that of R peak , log(τ peak ). They imply the presence of a very sharp cut-off in dynamics of material past the BLR in . These cases can not be considered as reliable. It is while for the accretion rate of 1, all shielding cases give a larger value for log(τ peak ). This means that for high accretion rate the presence of shielding is necessary. In other words, the higher the accretion rate is, the stronger shielding is required. A denser set of bins in the range of values of α and β in future studies will be more indicative.
As explained in section 6, since larger values of opening angle are of our interest, we chose the model showing the maximum P peak among the obtained values for all shielding models per the arbitrary accretion rate of 0.1 to perform a test with other blackhole masses. This maximum P peak corresponds to α = 2.0, however it does not necessarily mean the preference of Patch model over Wall model, or appropriateness of this α. We then repeated the simulation for this model for two other blackhole masses of 10 7 , and 10 9 M . Interestingly, while the Patch model with the adopted values of α, was too strong as mentioned above to be considered as reliable for accretion rates of 0.01 and 0.1, the consistent time-delay obtained for the blackhole mass of 10 9 M and accretion rate of 0.1 with α = 2.0 indicates the dependence of shielding on blackhole mass as well; i.e. the higher the blackhole mass is, the stronger shielding is required.
As for plotting the RL relation, since we only have the data for the blackhole mass of 10 8 M (except the case of α = 2.0), we take use of the slope of Bentz et al. (2013) relation. The RL relation based on reliable results for all adopted values of accretion rates are then shown in figure 3, and compared to observational data and analytical FRADO model (see also Naddaf et al. 2019 for RL relation based on τ in ). Time-delays are computed assuming the effective reflection of reprocessed radiation of continuum happens for the clouds located in farther side of AD, as in Czerny and Hryniewicz (2011), while the clouds in the closer side are almost self-obscured. As the whole material in the range between the onset of BLR and the maximum height can be irradiated by the continuum, the area between minimum/maximum time-delays corresponding to the minimum BLR in /maximum R(P peak ) obtained out of all reliable data in the table 1 for a given accretion rate are shown as shaded patterns in the figure 3. While the analytical FRADO only gives one line, the realistic 3-D approach covers almost the whole sample of observational data. The points corresponding to three blackhole masses for the test case of α = 2.0 with accretion rate of 0.1 are also depicted in the figure 3(B) which are located within the yellow strip, and follow the RL trend apparently. However testing with a wider-denser set of blackhole masses is also required in future studies to provide us with a full insight into the models.

CONCLUSIONS
In this preliminary work we theoretically studied the dynamics of matter in the BLR based on a realistic 3-D approach to FRADO model. We introduced two geometrical configurations for shielding effect, and based on the dynamics, we tried to predict the BLR size (radial extension and geometrical height) and its dependence on different Eddington ratios [and blackhole mass] to compare with observational data. The main findings can be concluded as below • The shielding must be included in the model, although low Eddington sources do not require it.
• The role of shielding must increase with the Eddington ratio of the source, as well as with blackhole mass.
• BLR size is predicted to be radially more extended compared to analytical FRADO which implies longer time-delays from the entire BLR.
• Future modelling of the emission line shapes based on the cloud dynamics will bring further constraints to the shielding geometry.

FUNDING
The project was partially supported by National Science Centre, Poland, grant No. 2017/26/A/ST9/00756 (Maestro 9), and by the Ministry of Science and Higher Education (MNiSW) grant DIR/WK/2018/12.