Radiation driven outflow in active galactic nuclei: the feedback effects of scattered and reprocessed photons
Abstract
We perform timedependent, twodimensional, hydrodynamical, numerical simulations to study the dynamics of a slowly rotating accretion flow from subpc to pc scales under the irradiation from the central AGN. Compared to previous work, we improve the calculation of the radiative force due to Xrays. More importantly, in addition to radiative pressure and radiative heating/cooling directly from the central AGN, in the momentum equation we also include the force due to the scattered and reprocessed photons. We find that the accretion flow properties change significantly due to this “reradiation” effect. The inflow rate at the inner boundary is reduced, while the outflow rate at the outer boundary is enhanced by about one order of magnitude. This effect is more significant when the density at the outer boundary is higher. The properties of outflows such as velocity, momentum and energy fluxes, and the ratio of outflow rate and the accretion rate, are calculated. We find that the efficiency of transferring the radiation power into the kinetic power of outflow is typically , far below the value of which is assumed in some cosmological simulations. The effect of the temperature of the gas at the outer boundary () is investigated. When is high, the emitted luminosity of the accretion flow oscillates. This is because in this case the gas around the Bondi radius can be more easily heated to be above the virial temperature due to its high internal energy. Another question we hope to address by this work is the socalled “subEddington” puzzle. That is, observations show that the luminosity of almost all AGNs are subEddington, while theoretically the luminosity of an accretion flow can easily be superEddington. We find that even when the reradiation effect is included and outflow does become much stronger, the luminosity, while reduced, can still be superEddington. Other observational implications and some caveats of our calculations are discussed.
keywords:
accretion, accretion discs  hydrodynamics  methods: numerical  galaxies: active  galaxies: nuclei1 Introduction
Active galactic nuclei (AGNs) are thought to play an important role in the processes of galaxy formation and evolution. The most important evidence is perhaps the remarkable correlations between host galaxy properties and the mass of the supermassive black holes (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Graham et al. 2001; Kormendy et al. 2009) observed in the past decade. Many theoretical studies have also confirmed that AGN feedback can affect the ambient environment effectively (e.g., Silk & Rees 1998; Ciotti & Ostriker 1997, 2001, 2007; Ciotti, Ostriker & Proga 2009; King 2003; Merloni et al. 2004; Springel, Di Matteo & Hernquist 2005a; Sazonov et al. 2005; Di Matteo, Springel, & Hernquist 2005; Hopkins et al. 2005; Murray, Quataert & Tompson 2005; Somerville et al. 2008; Ostriker et al. 2010).
The AGN feedback can be in the forms of radiative (e.g., Ciotti & Ostriker 1997, 2001, 2007, Ciotti et al. 2009) or mechanical energy input (e.g., Springel, Di Matteo & Hernquist 2005b; Ostriker et al. 2010). In terms of the radiative feedback, both momentum feedback and energy feedback work, i.e., the ambient gas could be accelerated by radiation pressure and heated by radiative heating. Mechanical feedback is via wind/outflow. For example, in the context of hot accretion flow (e.g., advectiondominated accretion flow), Yuan, Xie & Ostriker (2009; see also references therein) studied the effect of global Compton scattering, namely the scattering between the photons produced at the innermost region of the ADAF and the gas at large radius. They found that when the luminosity is higher than , the heating is so strong that the flow will be heated to be above the virial temperature thus the accretion will be stopped. It is shown that this physical mechanism is able to explain the intermittent activity of some compact young radio sources (Yuan & Li 2011).
Provided the gas is moderately ionized and can interact with the UV continuum through many UV line transitions, radiation pressure force due to spectral lines (i.e., line force) has been shown to be very efficient at producing winds from accretion disks. This can occur in the small scale of the innermost accretion flow (e.g., Murray et al. 1995; Proga, Stone & Kallman 2000, hereafter PSK00) and much larger scale from pc to 10 pc (Proga 2007a; Proga, Ostriker & Kurosawa 2008; Kurosawa & Proga 2009, hereafter KP09). Some highly blueshifted broad absorption line features seen in the UV and optical spectra of AGN can be explained by this mechanism (e.g. Chartas, Brandt & Gallagher 2003; Crenshaw, Kraemer & George 2003; Proga & Kallman 2004; Hamann et al. 2008). However, not all AGN outflows can be explained by linedriven mechanism because of, e.g., the toolow luminosity, toohigh ionization state, or both (e.g., Chelouche & Netzer 2005; Kraemer et al. 2006; see Proga 2007b for a review). This has also been shown to be the case of black hole Xray binaries which are believed to be the scaleddown version of AGNs (e.g., Miller et al. 2006, 2008; Luketic et al. 2010).
Therefore, other mechanisms are required, such as thermal (e.g., Begelman, de Kool & Sikora 1991) and, especially, magnetic driving (e.g., Blandford & Payne 1982; Emmering, Blandford & Shlosman 1992). Most recently, based on global HD and MHD numerical simulations, Yuan, Bu & Wu (2012; see also Li, Ostriker & Sunyaev 2013) show that outflow must exist in hot accretion flows. Moreover, Yuan, Bu & Wu (2012) proposed a Blandford & Paynelike mechanism of producing winds. This mechanism should be quite universal since it is a result of magnetorotational instability (MRI), which is widely believed to exist in accretion flows and is the mechanism of transporting the angular momentum. Previous works on studying the production of MHD winds usually require the preexistence of large scale poloidal field and treat the disk as the boundary condition. Different from these works, Yuan, Bu & Wu (2012) simulate the disk and outflow selfconsistently and simultaneously, and the preexistence of the largescale poloidal magnetic field is not required. However, the simulation was done for hot accretion flows. Thus the next step it remains to be probed explicitly whether it can be equally effective for thin disks, although in principle it is expected to work.
Among the abovementioned linedriven works, KP09 studied the momentum and energy interaction between the radiation coming from the central AGN and the accretion flow in the region between pc and pc. They focus on the linedriven outflow. Compared with previous works, the luminosity of the central AGN is not fixed, but selfconsistently determined by evaluating the mass accretion rate at the inner boundary. They found that outflow can be driven efficiently from the inflow. When the temperature of the gas at the outer boundary is high, , a strong correlation between the mass outflow rate () at the outer boundary and the luminosity () exists, with . When and when , the correlation becomes flatter, . Another interesting result is, they found that when the density at the outer boundary is high, the accretion luminosity of the central AGN can well exceed the Eddington value. The superEddington accretion mainly proceeds in the equatorial plane and is possible because the radiation flux from the disc is significantly reduced in the equatorial direction due to the geometrical foreshortening effect (Fukue 2000; Watarai et al. 2005).
In almost all previous works on radiatively driven outflow including KP09, they properly consider the forces due to electron scattering and UV line processes due to the “firstorder” radiation coming from the central AGNs. However, the effect of the locally produced photons via scattering, the bremsstrahlung, and line radiation, are all neglected. Thus when radiation is absorbed and then reemitted the dynamical effects of the latter are ignored, essentially violating energy conservation. In principle these photons could play an important role. For example, they can produce an additional force in the vertical direction and thus puff the disk up. This may make the interaction between the radiation and the gas significantly stronger. In the present paper we want to investigate this “reradiation” effect.
Another motivation to consider this effect in the present work is that we hope to address the socalled following “subEddington” problem. Observations to a large sample of AGNs with 407 sources show that the distribution of their Eddington ratio () is well described by lognormal, with a peak at and a small dispersion of 0.3 dex (Kollmeier et al. 2006). In other words, almost all AGNs are radiating below . Later Steinhardt & Elvis (2010) extended this study to a much larger sample consisting of 62185 quasars from the Sloan Digital Sky Survey and they reached the similar conclusion (but see Kelly & Shen 2013 for a different opinion). On the other hand, it has long been known from the black hole accretion theory that the rotating accretion flow can well radiate above . Such an accretion flow is called slim disk or optically thick advectiondominated accretion flow (Abramowicz et al. 1988). This analytical result has later been confirmed by radiative hydrodynamical numerical simulations (Ohsuga et al. 2005).
It has been argued that the galacticscale fueling events are likely to have a broad mass distribution so this should not be the reason for the subEddington luminosity of AGNs. One possible way to solve this puzzle is that, the accretion rates are determined by the selfregulation of the radiative feedback. In other words, the strong radiation from the inner accretion flow may drive a strong outflow thus reducing the accretion rate to below the Eddington rate (). This should not be effective in the innermost region close to the black hole, since the interaction between radiation and accretion flow has been properly considered in Ohsuga et al. (2005). At larger scale, however, the simulations by KP09 still find superEddington accretion even though the radiative feedback is included in their simulations. In this paper, we include the reradiation effect to do the calculation again, to see whether the “subEddington” puzzle can be solved. It is vitally important that the effects of radiation feedback at the Bondi radius be correctly included, since it is in this region that the rate of inflow is determined.
In this paper, following the approach of KP09, we perform twodimensional hydrodynamical (HD) timedependent simulations of accretion flow irradiated by the angle dependent UV flux from an accretion disk and isotropic Xray flux from a corona. The models include radiation forces due to electron scattering and line processes from both the firstorder radiation and the reradiation, and radiative cooling/heating. Special attention is paid to the reradiation effect. In §2, we describe some basic aspects of our model, especially the radiative force due to reradiation. The numerical simulation method and the initial and boundary conditions are also described in this section. The results of the simulations are given in §3, and §4 provides some discussion. Finally, §5 is devoted to a summary.
2 Model
The basic scenario of the model, together with most of the formula for calculating the radiative cooling and heating, are the same with those presented in KP09 (see also PSK00). The readers are referred to that paper for details. However, we do make some changes compared to KP09. In addition to the inclusion of reradiation force mentioned in §1, we also improve, we believe, the calculation of the radiation force due to Xray.
2.1 Hydrodynamics of the accretion flow
Our aim is to study the dynamics of the accretion flow between an inner radius and outer radius . In addition to the gravitational force, this gas will inevitably be irradiated by the radiation emitted from the central AGN. The radiation will heat and push the gas by Compton scattering and line absorption, if the gas is not fully ionized. The dynamics of the accretion flow is described by the following set of equations.
(1) 
(2) 
(3) 
where , and are the mass density, internal energy density, gas pressure, and velocity, respectively; describes the radiative cooling and heating, simply as the net heating rate (see PSK00); is the gravitational acceleration of the central object; (see equation 12) is the total radiation force per unit mass. We adopt an adiabatic equation of state , and consider models with the adiabatic index , with, of course, allowance for heating and cooling via the term in equation (3). The implementations of and are described as follows.
2.2 The central AGN
The central AGN is described by a supermassive black hole with mass surrounded by an accretion flow. The nature of the accretion flow in a luminous AGN is still an unsolved problem. Usually authors assume that it is described by a standard thin disk. However, a thin disk alone can’t explain the origin of hard Xray emission widely detected in AGNs. For simplicity, we assume that the accretion flow is a combination of a standard thin disk (Shakura & Sunyaev 1973) plus a spherical corona. They are responsible for the optical/UV and Xray emissions respectively. The radiation from the corona is assumed to be isotropic.
In the pointsource approximation limit, the radiation flux from the central Xray corona can be written as
(4) 
where is the fraction of the total luminosity () in Xray. We adopt throughout the paper. is the Xray optical depth in the radial direction,
(5) 
where is the Xray opacity. We assume that the attenuation is dominated by Thomson scattering, i.e., . However, the scattering opacity must be overestimated. This is because unlike true absorption, scattering merely redirects the photons. The scattered Xray photons must be replenished by photons scattered from other propagation lines. So another extreme is to take . In the present paper, we also test this case for comparison (see §4.2). The real opacity should be a certain value between these limits. The radiation from the standard thin disk . Consequently, the disc radiation flux at a distance from the center can be written as
(6) 
where is the fraction of the total luminosity in the disc emission, is the UV optical depth and we assume
In our simulation, at each time step, the accretion luminosity will be calculated selfconsistently based on the mass inflow rate at the inner boundary . It is given by
(7) 
where is the mass inflow rate at and at a given time , is the radiative efficiency. We adopt in this paper. Following KP09, the mass accretion rate at a given time is actually timeaveraged over some time interval to reduce the fluctuation level. Obviously, a time lag is expected between the change of this accretion rate to the disk and the change of the luminosity from the central AGN . Since in the present paper the accretion flow we consider has small angular momentum, a small disk with a radius of will be formed within the inner boundary of our simulation. In this case, the time lag is roughly equal to the accretion timescale at the outer boundary of this small disk. If the small disk can be described by a standard thin disk, the lag time can be approximated as ; if it is a slim disk, the time lag should be . Here is the viscosity parameter, and we set in this paper. The black hole mass is assumed to be constant. The value of is not important for our aim.
One caveat here is that we assume the accretion rate onto the black hole is equal to the inflow rate at . This may not be true, however. Almost all global numerical simulations to hot accretion flows show that the accretion rate should decrease with decreasing radius (see Yuan, Wu & Bu 2012 for a review). The reason is identified to be because of continuous mass loss in outflow produced by an MHD process similar to the BlandfordPayne mechanism (Yuan, Bu & Wu 2012). This mechanism in principle also works in the case of a cold thin disk. Actually, the initial global radiation MHD numerical simulations to thin disk by Ohsuga et al. (2009; see also Ohsuga & Mineshige 2011) do show significant outflow, although the exact mechanism remains to be probed. Since the outflow in the context of a thin disk is still an open question, in the present paper we simply assume that the accretion rate is constant within and thus our computed luminosity is an upper bound on the true luminosity for given outer boundary conditions.
2.3 Heating/cooling and radiative force of the gas
The interaction between the radiation and the gas and the radiative processes we consider includes Compton heating/cooling, Xray photoionization and recombination, bremsstrahlung, and line cooling (cf. PSK00; see also Blondin 1994). The net cooling rate depends on the photoionization parameter and the characteristic temperature (or radiation temperature) of the radiation from the central AGN (). The gas photoionization state is determined by the photoionization parameter
(8) 
where is the number density of the local gas located at distance from the central AGN, is the mean molecular weight and is fixed to be 1. Only Xray radiation ionizes the gas. The sum of photoionization heatingrecombination cooling rate () and Compton heating/cooling rate () is
(9) 
with
(10) 
(11) 
Here (Sazonov et al. 2004) is the “characteristic temperature” of the Xray radiation and it is times larger than the Compton temperature (refer to PSK00 and Proga 2007a). For a lowluminosity AGN, will be much higher, K (Yuan, Xie & Ostriker 2009).
2.4 Radiative force
Model  Model  Reradiation  
Number  
(1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9) 
1  K5a  500  10  2  no  1.19(0.19)  15.31(9.12)  0.76(9.12) 
2  K6a  500  20  2  no  1.85(0.43)  25.61(24.29)  0.84(0.92) 
3  K6c  1250  20  2  no  1.72(0.14)  29.05(17.12)  0.98(0.57) 
4  K6d  2550  20  2  no  1.71(0.23)  31.17(13.22)  1.10(0.53) 
5  K7a  500  50  2  no  3.97(1.10)  21.54(27.46)  0.35(0.58) 
6  K7b  625  50  2  no  3.69(0.78)  20.80(17.05)  0.33(0.24) 
7  K7d  2500  50  2  no  3.14(1.03)  49.82(42.21)  1.03(0.95) 
8  K8a  500  100  2  no  5.74(0.70)  21.52(14.41)  0.22(0.13) 
9  K8d  2500  100  2  no  5.70(1.75)  34.24(50.46)  0.35(0.51) 
10  R5a  500  10  2  yes  0.77(0.04)  37.95(8.71)  2.88(0.70) 
11  R5b  625  10  2  yes  0.94(0.04)  31.53(5.07)  1.96(0.34) 
12  R5c  1250  10  2  yes  1.10(0.10)  24.28(10.67)  1.28(0.55) 
13  R5d  2500  10  2  yes  1.12(0.11)  26.88(8.05)  1.40(0.45) 
14  R6a  500  20  2  yes  1.19(0.16)  67.30(10.90)  3.35(0.78) 
15  R6b  625  20  2  yes  1.28(0.17)  65.14(10.73)  3.02(0.69) 
16  R6c  1250  20  2  yes  1.48(0.18)  39.10(12.13)  1.54(0.45) 
17  R6d  2500  20  2  yes  1.60(0.18)  44.04(11.45)  1.61(0.45) 
18  R7a  500  50  2  yes  2.02(0.02)  139.13(5.03)  3.99(0.14) 
19  R7b  625  50  2  yes  2.03(0.03)  137.86(4.83)  3.93(0.14) 
20  R7c  1250  50  2  yes  1.81(0.03)  127.04(11.80)  4.08(0.39) 
21  R7d  2500  50  2  yes  1.60(0.03)  194.29(1.92)  7.04(0.14) 
22  R8b  625  100  2  yes  3.14(0.27)  202.05(28.40)  3.75(0.59) 
23  R8c  1250  100  2  yes  4.56(1.09)  36.45(54.69)  0.49(0.82) 
24  R8d  2500  100  2  yes  4.14(1.20)  126.12(97.44)  2.14(1.91) 
25  R7cX0  1250  50  2  yes  1.76(0.22)  120.61(9.67)  4.03(0.58) 
26  R7clT  1250  50  0.5  yes  1.61(0.27)  28.59(11.05)  1.08(0.50) 
27  R7clTX0  1250  50  0.5  yes  1.19(0.01)  52.17(3.28)  2.54(0.16) 
28  R7chT  1250  50  6  yes  2.94(0.04)  245.72(10.33)  4.85(0.22) 
29  R7chTX0  1250  50  6  yes  2.26(0.46)  338.04(29.36)  8.98(1.85) 
30  R7chT2  1250  50  10  yes  2.02(0.97)  461.93(100.00)  15.43(5.07) 
31  R7chT2X0  1250  50  10  yes  2.74(0.51)  508.50(57.05)  11.18(3.12) 
Note: Values in brackets in Columns 7, 8 and 9 are the normalized standard deviations of the time series values. The symbol ‘X0’
in the model names means in those models, ‘lT’ means lower temperature , ‘hT’ means higher temperature .
Compared to KP09, in our present work we make two changes when calculating radiative forces. The radiative forces considered in KP09 are due to Compton scattering and line processes. The former is the only available radiative force if the gas is fully ionized. Obviously, this force is important only when the luminosity from the central AGN approaches the Eddington value. However, the force can be significantly enhanced if the gas is moderately ionized, since the opacity is much higher than scattering opacity due to the boundbound and boundfree transitions. So we should include the line force. For this force, following KP09, we assume that only optical/UV photons have contributions, but neglect the possible force line due to some metal line in the soft Xray band. Regarding the sources of photons, as we have emphasized in §1, in addition to the radiation from the central AGN, the “local” photons from the local radiative processes and scattered photons originally from the central AGNs are also included. Both of them contribute to the radiative force,
(12) 
The inclusion of the force due to the “reradiation” process () is the first change we make compared to KP09. It turns out that the inclusion of this effect does significantly change the results, as we will state in the following section.
For Xray photons, KP09 only consider the force due to Thomson scattering. However, if the gas is not fully ionized as in our present case, absorption processes such as photoionization should also contribute to the force. We characterize this force by
Based on the above considerations, the total radiative force due to the central AGN photons then can be written as:
(13) 
Here is the scattering coefficient for free electrons, and is the force multiplier (Castor et al. 1975) – the numerical factor which parameterizes by how much spectral lines increase the scattering coefficient (see PSK00 for details). We have calculated the dependence of the multiplier on the temperature and the ionization parameter . We found that the multiplier could be up to the order of , providing that the gas is moderately ionized. It drops rapidly when temperature or the ionization parameter is larger than some critical value (eg., , or ). Specifically, it is safe to ignore the line force ( of electron scattering force) when or , since then the gas is almost fully ionized.
Now we calculate the second term in equation (12), i.e, the reradiation force. Exact treatment of this force requires the full radiative transfer calculation which is beyond the scope of the present paper. Because the accretion flow is rotating (see below), the flow has a disklike shape. In this case, for simplicity we can adopt the assumption of planeparallel approximation. The “reradiation” photons would eventually escape from the accreting flow, and produce a net force in the vertical direction. Consider a rectangular box located within the disk, with the bottom of the box overlapped with the equatorial plane of the accretion disk while the upper side at height . “Reradiation” photons are steadily produced within this box and then escape. Given the planeparallel approximation and symmetry, the net radiative flux due to reradiation (and thus the force) should be only in the upper side of the box. Applying the Gauss theorem, we can easily calculate the vertical radiative force,
(14) 
where
(15) 
is the source term due to the firstorder scattered photons of the radiation from the central AGN; while
(16) 
and
(17)  
are the rate of bremsstrahlung and line cooling (cf. PSK00), respectively. The parameter is introduced to control line cooling, represents optically thin cooling, and represents optically thick cooling. We set throughout this paper. In general, the line cooling dominates over the other cooling processes. The formulae above only show the reradiation force caused by electron scattering while the line force is ignored. So the effect of “reradiation” calculated in this paper should be regarded as a lower limit in this sense. We did include the line force in some models and found that the change of the result is not very significant. This is perhaps because the reradiation force monotonically increases with height, and at a large height, the gas temperature becomes high so the line force is not important. One thing to note is that all our treatment is approximate; a full radiative transfer calculation is required in the future.
2.5 Simulation setup
We solve the dynamical equations (13) using the code ZEUSMP (Stone & Norman 1992a; Hayes et al. 2006). The simulations are performed in spherical polar coordinates assuming axial symmetry about the disc rotational axis () and in twodimensions. Our computational domain is defined to occupy the angular range and the radial range . In KP09, both and are fixed, i.e, they are same in all models. In the present work, we want to check the convergence for various values of , so we set various in different models but adopt , where is the innermost stable circular orbit (ISCO hereafter) of a Schwarzschild BH and . We set the BH mass, . Our standard numerical resolution in the direction consists of 144 zones with the zone size ratio, and consists of 64 zones in the direction with . Gridding in this manner ensures good spatial resolution close to the inner boundary.
For the initial condition, following KP09, we set the density and temperature of gas uniformly , i.e., and everywhere in the computational domain. The initial velocity of the gas is assigned to have the distribution described by the following latitudedependent angular momentum at the outer boundary.
(18) 
where is the “circularization radius” on the equatorial plane. In the present work, we set . The boundary conditions are specified in the following way. We apply an axisofsymmetry boundary condition at the pole (i.e., ) and reflecting boundary conditions at . For the inner and outer radial boundaries, we apply an outflow boundary condition (i.e., to extrapolate the flow beyond the boundary, we set values of variables in the ghost zones equal to the values in the corresponding active zones, see Stone & Norman (1992a) for more details). To represent steady conditions at the outer radial boundary, during the evolution of each model we apply the constraints that all hydrodynamical variables except radial velocity (always let it float) in the last zone in the radial direction to be equal to the initially chosen values, i.e., , , and , when (inflowing). We allow all hydrodynamical variables to float when (outflowing). This approach is to mimic the situation where there is always gas available for accretion.
So the free parameters in our model include the inner radius of the computational domain, the density and temperature of inflowing matter at the outer boundary. In this paper we focus mainly on the case of . But we do also consider other for comparison. Various and are adopted. The details of various models are listed in Table 1. The “yes/no” shown in Column 6 indicates whether reradiation is included in the simulations.
3 Results
Some results of various models are given in Table 1. Especially, Column 7 gives the timeaveraged luminosity calculated according to the mass accretion rate at the inner boundary. This may exceed the luminosity as seen at infinity due to absorption of radiation within the flow. Models 1, 2, 5 and 8 are identical to the models 5, 6, 7 and 8 in KP09 except that we use our improved term when calculating Xray force, i.e., the first term in equation (13). The last characters in the model name, i.e., “a”, “b”, “c”, and “d”, correspond to in units of , respectively. We define the mass inflow and outflow rates as below:
(19) 
(20) 
The timeaveraged mass outflow rate is given in Column 8. The Column 9 gives
(21) 
which is the ratio of mass outflow rate at the outer boundary and mass accretion rate at the inner boundary. It is important to note that the fraction of the material which is inflowing at the outer boundary which actually accretes on the black hole, defined as , is from the mass conservation ; so that in cases where the net accretion, is small.
3.1 Effects of reradiation
We take models 7 series (i.e., K7* and R7*) as our fiducial models and explore their properties in details. Fig. 1 compares the radial profiles of mass inflow and outflow rates, gas density and temperature between models K7a and R7a. The former does not consider reradiation while the latter does. Fig. 2 shows the snapshots of contours of gas density and temperature of the two models. From the two figures we can find the following changes after reradiation is taken into account.

Close to the equatorial plane the gas density decreases when reradiation is included, as shown by the bottomleft panel of Fig. 1. However, we can see from Fig. 2 that away from the equatorial plane the density becomes higher, i.e., the density scaleheight becomes larger. The reason is that the accretion flow expands upward due to the inclusion of the vertical radiative force from reradiation photons, then the vertically expanded inner edge of the accretion flow shields the radiation from the central AGN, which makes the density higher and the mass inflow rate larger.

From Fig. 2, we can see that the temperature becomes lower in most of the region. This is because the increase of density produces stronger radiative cooling. The decrease of temperature close to the equatorial plane is because the compression work induced by accretion becomes weaker due to the inclusion of the vertical force.

Perhaps the most important effect of including reradiation is that outflow becomes nearly one order of magnitude stronger, as clearly shown by the comparison between the top two panels of Fig. 1. This is easy to understand. As we state above, the inclusion of the vertical force due to reradiation photons makes the density scaleheight larger. Since the irradiation force from the central AGN , the gas will be more easily blown away if located at higher latitude.

The mass inflow rate at the inner boundary decreases by a factor of 2. This is obviously because of the enhanced mass loss via the outflow. It is interesting to note that the mass inflow rate at the outer boundary increases by a factor of 2. The reasons are twofold. One is that the radiation from the central AGN decreases because of the decrease of the inflow rate at the inner boundary, thus the outward radiation force becomes weaker. Another reason is that, as we stated above, the gas density in most of the region increases. This gas will more effectively shield the gas at the outer boundary from irradiation by the central AGN.

The flow morphology also changes with the inclusion of reradiation. Comparing the top and bottom panels of Fig. 2, we can see that in Model K7a, the inflow falls inward in a relatively small range of polar angle, i.e., . But in Model R7a, i.e, the reradiation effect is included, the inflow occurs in a relatively larger range, . In addition, it is interesting to note that above the inflow region, the motion of the gas is turbulent in Model K7a; while in Model R7a, the motion is less turbulent.
KP09 studied the correlation between the mass outflow rate at the outer boundary and the luminosity in units of Eddington luminosity , where is the Eddington luminosity. They found that when is not too high, , outflows exist only when . Specifically, when , as in our paper, a strong correlation is found for , i.e., with . The correlation index only becomes slightly smaller when the temperature at the outer boundary is higher. In the case of , becomes almost saturated when . KP09 explains this result as because the inflow is squeezed to a narrow range of when becomes high.
We have studied the effect of reradiation on the correlation. Fig. 3 shows the result. All the lines and dots except the four solid blue dots are taken from KP09. The four solid blue dots are the results when reradiation effect is included for . We can see that the correlation index roughly remain the same, and the correlation again becomes almost saturated when . The only change is the normalization: the outflow rate increases by almost one order of magnitude as we have stated above.
At last, as we have mentioned above, the reradiative force makes the density in most of the region larger. In this case, the attenuation of Xrays also becomes much more significant, thus the photoionization parameter decreases. As we have mentioned in section 2.4, the line force multiplier usually increases with decreasing . In other words, the line force should increase rather than decrease when we consider the reradiation effect.
3.2 Convergence with the inner boundaries
It is necessary to check whether the result converges with various inner boundaries . With this aim, we have implemented a series of runs (see Table 1) by varying the value of .
Fig. 4 shows the simulation results for the model R7 series. The top and middle panels show the radial profiles of inflow, outflow, and net rates for the four models with different , while the bottom panel shows the time variation of the luminosity of the models. The values of the outflow rate of the four models are almost same. Moreover, the radius where the outflow begins to become equal to the net rate is also roughly same, i.e., . For the inflow rate at the inner boundary, models R7a and R7b are almost same; while model R7c is only about 10% smaller than model R7a and model R7d is only about 12% smaller than model R7c. The deviations are within the acceptable range. Therefore, we think model R7 series are convergent as we vary the radius of the inner boundaries.
We also did the simulations with different for model R5, R6 and R8 series. The results are also listed in Table 1. For model R5 series, the largest discrepancy of timeaveraged luminosity among the four models is . For model R6 and R8 series, they are and , respectively. We found that the results are similar when reradiation is not included, as shown by Models K6a, K6c, and K6d. The change of the luminosity with is complicated: the luminosity can increase or decrease with increasing . Physically, when increases, the attenuation of Xray flux will become weaker, thus the corresponding radiation force becomes stronger. But on the other hand, when Xray radiation becomes weaker, the ionization will become weaker and correspondingly the line force stronger. Thus the final result depends on the competition of these two effects. But one relatively robust result is that the outflow mainly originates from . Another noteworthy result is that the emitted luminosities shown in Fig. 4 are all superEddington, although the reradiation effect has been taken into account. So the subEddington puzzle is not solved for these simulations.
Model  

R5c  24.28(10.67)  800  21.57(5.72)  174.16(48.80)  3.77(0.83) 
R6c  39.10(12.13)  1000  43.03(16.81)  547.17(459.39)  7.56(4.09) 
R7c  127.04(11.80)  700  98.00(5.52)  1098.91(37.47)  12.31(0.63) 
R8c  28.94(44.12)  600  23.58(64.40)  1110.45(10393.40)  91.55(1206.77) 
Note: The values in brackets are the normalized standard deviations of the time series values. For the Column , is
the timeaveraged mass fluxweighted radial velocity of outflow at (see panel c of Fig. 5). To avoid the influence of the outer
boundary, we actually calculate all the quantities at pc.
3.3 Dependence on gas density at the outer boundary
We can see from Table 1 that the significance of reradiative effects depends on the outer boundary density . For example, when , the timeaveraged luminosity decreases by a factor of 2 from model K7a to model R7a. The outflow rate at the outer boundary increases by a factor of about 7 from model K7a to model R7a. But when , the luminosity decreases by only , and the outflow rate increases by a factor of about 2.5 from model K5a to R5a. Usually we find that the smaller the gas density at the outer boundary is, the weaker the reradiation effect will be. This result is easy to understand. The thickness of the accretion disk is determined by the vertical component of the gravitational force, gradient of gas pressure, and the reradiation force. The former two forces are proportional to density while the reradiation force is proportional to the square of density (refer to equation 14). Therefore, the disk will become thicker and more exposed to the irradiation from the central AGN.
3.4 Properties of outflows
In this section, taking models R5c, R6c, K6c, R7c and R8c as examples, we study the properties of the outflows. This is a sequence of increasing outer density and larger and lager rates of nominal inflow. These properties could be used to compare with observations. Some main properties are listed in Table 2. From the left to right columns, they are the mass flux, radial velocity, momentum flux, kinetic energy flux, and thermal energy flux of outflow at the outer boundary, respectively.
Fig. 5 shows some additional properties of outflow for models R7c (the solid lines), R6c (the dashed lines) and K6c (the dotdashed lines). Panels (a) and (b) show the angular distribution of mass fluxes and densities at , respectively. For model R7c, we can see that most of the outflow occurs close to the disk surface, in the range of . In the range of it is the inflowing region and there is almost no outflow. Panel (c) shows the radial profiles of the mass flux weighted radial velocities of outflows. The typical value of is . We can see that in most of the radial range, the radial velocity does not keep increasing as we naively expect because of the acceleration by the radiation force. It even decreases somewhat at large radii. This is because of the continuous injection of new gas into the outflow. Panel (d) shows the angular distribution of at . It monotonically increases from close to the surface of inflowing region at to close to the axis for models R7c and R6c. The dotdashed line (model K6c) sharply decreases at highlatitude (small ) region. This is because some inflowing clumps often appear in that region which block the outflow (cf. the top panels of Fig. 2, model K6c has similar structure). Due to the same reason, of model K6c is very weak at highlatitude region (see the dotdashed lines in panel (a)). The black and blue lines in panel (e) show the radial profiles of the momentum fluxes of outflow () and radiation from the central AGN, respectively. In panel (f), the solid and dotted lines respectively represent the radial profiles of the kinetic () and thermal energy fluxes () of outflows. Their definitions are:
(22) 
(23) 
(24) 
We see that the momentum flux keeps increasing outwards, because more and more momentum of the radiation is transferred into the outflow gas. But we note that even at , the momentum flux of outflow is still several times smaller than that of radiation. It will be interesting to adopt a larger to see how the momentum flux of outflow can finally become closer to that of the radiation. From Panel (f) we see that most of the power of outflow is in the form of kinetic energy rather than thermal energy. The highest power is reached at , which is for model R7c. As a comparison, for this model the luminosity from the central AGN is . This is about three orders of magnitude higher than the kinetic power of outflow. Compared to the model of KP09, the efficiency of transferring the radiation into kinetic power of outflow is one order of magnitude higher. But it is still much lower than , which is often assumed in cosmological simulations (see discussions in Kurosawa, Proga & Nagamine 2009).
Comparing model R6c and K6c, we can see that model R6c has higher mass outflow rate at highlatitude region at and higher kinetic and thermal power at all radii. The outflow velocity is also slightly higher. Compared with model R6c, mass outflow rate of model R7c is higher most obviously at the moderatelatitude region at as mentioned earlier. This is due to the increase of density at that region (see panel (b)). Compared to model R6c, the luminosity of model R7c is higher hence the radiative force larger, so outflow stronger at small radii.
To better understand the overall energetics of the flow, we calculate the thermal energy of the input gas, the released gravitational energy, the energy absorbed from the center engine, the energy advected inside , the energy advected outside and the energy reemitted out per unit time for model R7c. We find that the first term is the smallest, . The second term is comparable to the fourth term, . The fifth term () is much smaller than the third and sixth terms (). In other words, the energy absorbed from the center is comparable to the energy reemitted out.
The value of (defined in equation 21) is given in the last column of Table 1. We can see that it becomes significantly larger when the reradiation effect is included, . This implies that outflow rate is larger than the accretion rate. It is interesting to note that this is also the case for hot accretion flow (Yuan, Wu & Bu 2012) and found by Li et al. (2013). Of course, in that case, the mechanism of producing outflow is completely different.
4 Discussion
4.1 Effect of changing the temperature at the outer boundary
So far we have fixed K. It will be interesting to see the effect of varying . With K, the Bondi radius is . If we adopt higher but maintain the outer boundary unchanged, the Bondi radius will be within our computational domain. This is an advantage since Bondi radius is the place where the inflow rate is determined (Li, Ostriker & Sunyaev 2013). We choose three values: K, K and K. The corresponding of the two latter cases is smaller than now. These models with new are listed in the end of Table 1 (Models 2631).
The top panel of Fig. 6 gives the emitted luminosity as a function of time for Models R7c (red line; K), R7clT (blue line; K); R7chT (green line; K), and R7chT2 (black line; K). It is interesting to note that for model R7chT2 which has the highest and smallest (), the light curve oscillates regularly. Similar behavior has been found before such as the feedback study in elliptical galaxies by Novak, Ostriker & Ciotti (2011). Such an oscillation is caused by radiative feedback. For model R7chT2, at some time the strong radiative heating heats the gas at a certain radius () above the local virial temperature. Then the accretion is partly suppressed since the gas becomes unbound. The accretion rate will thus begin to decrease, and subsequently the radiation from the central AGN will decrease. This makes the radiative heating weaker and the temperature of the gas will decrease until below the virial value. So the gas supply then recovers and the accretion becomes active again. The timescale of such a cycle should be determined by the accretion timescale at . From the top panel of Fig. 6, we see that the timescale is . This is roughly the freefall timescale at the outer boundary , which is also the accretion timescale for our low angular momentum accretion.
We would like to emphasize that for such an oscillation behavior to occur, it is not necessary to have . In fact, Yuan, Xie & Ostriker (2009) find the same oscillation result even when when they consider the global Compton heating effect in hot accretion flows. The reason why the other three models do not have such oscillation is perhaps as follows. One is that the luminosity of the active phase from the other three models is smaller (refer to Fig. 6), thus the radiative heating is weaker. Another reason is that the specific thermal energy of the gas in Model R7chT2 is the highest among the four models at the same radius, therefore it is the easiest to be heated to be above the virial temperature.
In this context we would like to mention that the peak luminosity found in Yuan, Xie & Ostriker (2009) is only , much smaller than that in Model R7chT2, which is superEddington. One of the main reasons for the discrepancy is that the emitted spectrum from the innermost accretion flow is different. In the current work we assume a typical quasar spectrum, i.e., and . Such a spectrum is likely produced by a standard thin disk plus a corona. However, in Yuan, Xie & Ostriker (2009), the accretion flow is hot (such as advectiondominated accretion flow) thus the emitted spectrum is quite different. The radiation temperature of the spectrum in that case is and . In this case the Compton heating will become much stronger (e.g., Proga 2007a). In addition, the angular distribution of radiation is also different. In the case of hot accretion flow, the radiation is nearly spherically symmetric. This again increases the heating close to the equatorial plane. We have run some test simulations and confirmed the above effects. And lastly, the effectiveness of Compton heating is proportional to density of the flow. The specific angular momentum of the accretion flow in Yuan et al. (2009) is much larger than that in the current work. Therefore, for the same density, the corresponding accretion rate and therefore the luminosity of Yuan et al. (2009) will be significantly smaller than in the current work.
4.2 Effect of Xray optical depth
So far we have treated the Thompson scattering opacity for the calculation of Xray optical depth as if it were absorption. However, as we have mentioned above, this is likely an overestimation (we neglect Xray absorption). In this section, we want to check the optical depth effect of the outward propagation of Xray from the central source by assuming . But in the calculation of the force we still keep the first term in equation (15) for comparison purpose. The last few models labeled with ‘X0’ in Table 1 are such models.
Fig. 7 compares the logarithmic density and temperature contours with (model R7c, top panels) and without (model R7cX0, bottom panels) Xray opacity. When , there will be no Xray attenuation. In this case, the ionization parameter increases and Xray heating is enhanced especially at large radii (cf. equations 1011). We can see from the bottomright panel of the figure that the gas temperature at low latitudes increases significantly. The gas density near the pole decreases greatly as the Xray pushes the gas away. Once the gas density becomes lower, the Xray heating becomes less effective, hence the temperature decreases. Note in table 1 the high values of (cf. equation 21) in the models having and . The outflow rates exceed the net accretion rates.
The top panel of Fig. 6 shows the light curves of various models which have been shown in the top panel of the same figure, except here . In the top panel, only the model with the highest (model R7chT2) oscillates. But here with the exception of model R7clTX0 which has the lowest , the light curves of all other three models oscillate. The discrepancy is because of the enhanced Xray heating. Comparing model R7chT2 in the top panel and model R7chT2X0 in the bottom panel, we can see that the timescale of the variability becomes smaller. This indicates that the virial radius becomes smaller. Another effect is that the amplitude of oscillation of model R7chT2X0 becomes smaller now. This is perhaps because a smaller fraction of the gas located at is heated to be unbound now.
4.3 Observational implications
We now briefly discuss the observational implications of our results. The first issue is the origin of outflow observed in AGNs. They have been widely observed in the form of absorption lines in a variety of AGNs, including radioquiet and radio loud, Seyfert and quasars, lowluminosity or highluminosity ones (see Crenshaw, Kraemer & George 2003 for a review). Especially, recently Tombesi et al. (2010, 2011, 2012) detected ultrafast outflows with velocities up to via Fe Kshell absorption lines. They found that the ultrafast outflows are located in the interval from the black hole on average. The mass outflow rates are constrained between and the average lower and upper limits on the mechanical power are and , respectively. These features are difficulty to explain with our simulations. In our models, the highest outflow velocity is about (see panel (d) of Fig. 5), and the outflows starts from (see Fig. 4), the mass outflow rates and outflow kinetic power at the outer boundary are and (see Table 2). However, based on these results we can not rule out the “linedriven” models as the origin of ultrafast outflow. This is because our simulations concentrate only on the large scale of the accretion flow, pc, the innermost accretion flow is not included and substantial radiative acceleration may occur within (e.g., PSK00). We did adopt different inner boundary and found that the properties of outflow, including the location of origin, remain largely unchanged. This is likely because the angular momentum of the accretion flow we consider is rather small so there is no real disk formed in our simulation domain. On the other hand, the discrepancies between our simulation results and observations indicate that at least the ultrafast outflow must come from the innermost region of the accretion disk.
5 Summary
In this paper we study the dynamics of accretion flow within pc from the central supermassive black hole irradiated by the central AGN. We focus on the outflow driven by the radiation via the Compton scattering and the line force. Compared to the previous works, the main improvement is that we include the reradiation force, which is the additional force produced by the scattered photons and the reprocessed photons. We find that the results change significantly. The accretion flow now becomes thicker and is thus more exposed to the radiation from the central source. Consequently, the outflow rate measured at the outer boundary increases by about one order of magnitude. We find that the correlation between the Eddington ratio and the outflow rate originally found in previous works still exists. The properties of the outflows are calculated and compared to observations, such as the ratio of the outflow rate and the accretion rate. We find that some observed outflows can not be explained by this model, but note that the inferred high ratios of occur naturally in our high Eddington ratio models. When the density of the gas at the outer boundary is sufficiently high, we still find superEddington accretion, even though we include the reradiation effect and the outflow becomes stronger. We also investigate the effect of different temperatures at the outer boundary (). We find that when is large ( K), the emitted luminosity from the accretion flow oscillates significantly because of the radiative heating feedback.
One of the main aims of the present research is the subEddington puzzle. Observations show that most AGNs are subEddington, while on the other hand both analytical theory and numerical simulations of accretion flows have clearly shown that superEddington accretion can exist. In this paper, we hope to investigate whether the feedback at larger scale can solve the puzzle, especially after the reradiation effect is taken into account. Unfortunately, we find that the superEddington accretion still exists (refer to Fig. 6) but at a reduced rate.
There are some caveats in our study. These issues can be improved in the future, which may be helpful to solve the subEddington puzzle. The first one is that we need to consider the probable large angular momentum of the accretion flow. When the angular momentum is large, we must include viscosity and a disk will be formed. In this case, such a disk itself will emit strong radiation, which may help to drive an outflow. Moreover, when the angular momentum is large, a strong outflow will be produced intrinsically by other mechanisms (Yuan, Bu & Wu 2012; Bu et al. 2013). For example, outflow can be produced by buoyancy or magnetocentrifugal force in a hydrodynamical or magnetodynamical accretion flow, respectively (Yuan, Bu & Wu 2012; see also Bai & Stone 2013). The second caveat is that in the current simulations, we only consider the radiation from within and we do not investigate the dynamics within . However, strong outflows have been shown to exist from the innermost accretion flow and they will interchange their momentum and energy with the gas. So the feedback of the outflow should also be included (e.g., Ostriker et al. 2010; Novak et al. 2011; Choi et al. 2012). All the abovementioned factors are expected to be able to make the outflow stronger and thus the inflow rate at and luminosity smaller. In addition, a simplified approach is adopted when we deal with the interaction between the radiation, from both the central force and reradiation, and the gas. Obviously a proper radiation transfer calculation is desired. And lastly, the opacity due to dust can also potentially greatly enhance the outflow (e.g., Novak et al. 2012; Roth et al. 2012; FaucherGiguère, Quataert & Hopkins 2013).
The values of found in our high inflow models are consistent with recent assessments of Moe et al. (2009), Dunn et al. (2010) and the references therein that winds carry out more mass in bright quasars than is being accreted on the central black holes.
Acknowledgments
FY thanks the hospitality of Princeton University where part of the work was done. We thank Daniel Proga and Ryuichi Kurosawa for useful discussions/comments, and the anonymous referee for constructive comments. This work was supported in part by the Natural Science Foundation of China (grants 11103059, 11121062, and 11133005), the National Basic Research Program of China (973 Program 2009CB824800), and the CAS/SAFEA International Partnership Program for Creative Research Teams. The simulations were carried out both at Shanghai Supercomputer Center and the Super Computing Platform of Shanghai Astronomical Observatory.
Footnotes
 pagerange: Radiation driven outflow in active galactic nuclei: the feedback effects of scattered and reprocessed photons–Radiation driven outflow in active galactic nuclei: the feedback effects of scattered and reprocessed photons
 We adopt this approximation since the formula for calculating the line force multiplier only applies to UV optically thin case (cf. the related discussion in KP09). For the same reason, the term does not appear in equation (13). We did some test calculations by including and found that the results change only slightly. But note that the gas can still get momentum and energy from the UV radiation although the radiation does not loss its momentum and energy. In this sense, the conservation of momentum and energy is not guaranteed.
 We note that this simplification will underestimate the force, since in principle we should use only Xray heating rate. We find that this approximation is feasible because the force is in general dominated by the line force.
References
 Abramowicz M. A., Czemy B., Lasota J. P., Szuszkiewicz E., 1988, ApJ, 332, 646
 Bai X. N., Stone J. M., 2013, ApJ, 767, 30
 Begelman M., de Kool M., Sikora M., et al., 1991, ApJ, 382, 416
 Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
 Blondin J. M., 1994, ApJ, 435, 756
 Bu D., Yuan F., Wu M., Cuadra J., 2013, MNRAS (arXiv:1303.2473), accepted
 Castor J. I., Abbott D. C., Klein R. I., 1975, ApJ, 195, 157
 Chartas G., Brandt W. N., Gallagher S. C., et al., 2003, ApJ, 595, 85
 Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 152
 Chelouche D., Netzer H., 2005, ApJ, 625, 95
 Ciotti L., Ostriker J. P., 1997, ApJ, 487, L105
 Ciotti L., Ostriker J. P., 2001, ApJ, 551, 131
 Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038
 Ciotti L., Ostriker J. P., Proga D., 2009, ApJ, 699, 89
 Crenshaw D. M., Kraemer S. B., George I. M., 2003, ARA&A, 41, 117
 Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
 Dunn J. P., Bautista M., Arav N., et al., 2010, ApJ, 709, 611
 Emmering R. T., Blandford R. D., Shlosman I., et al., 1992, ApJ, 385, 460
 FaucherGiguère C. A., Quataert E., Hopkins P. F., 2013, MNRAS (arXiv: 1301.3905), in press
 Fukue J., 2000, PASJ, 52, 829
 Ferrarese L., Merritt D., 2000, ApJ, 539, L9
 Gebhardt K., Bender R., Bower G., Dressler A., Faber S. M., Filippenko A. V., Green R., Grillmair C., Ho L. C., Kormendy J., Lauer T. R., Magorrian J., Pinkney J., Richstone D., Tremaine S., 2000, ApJ, 539, L13
 Graham A. W., Erwin P., Canon N., Trujillo I., 2001, ApJ, 563, L11
 Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., udDoula A., Mac Low M.M., 2006, ApJS, 165, 188
 Hamann F., Kaplan K. F., Hidalgo P. R., Prochaska J. X., HerbertFort S., 2008, MNRAS, 391, L39
 Hopkins P. F., Hernquist L., Cox T. J., et al., 2005, ApJ, 630, 705
 Kelly B. C., & Shen, Y. 2013, ApJ, 764, 45
 King A. R., 2003, ApJ, 596, L27
 Kollmeier J., Onken C. A., Kochanek C. S. et al., 2006, ApJ, 648, 128
 Kormendy J., Bender R., 2009, ApJ, 691, L142
 Kraemer S. B., Crenshaw D. M., Gabel J. R., et al., 2006, ApJS, 167, 161
 Kurosawa R., Proga D., 2009, MNRAS, 397, 1791 (KP09)
 Kurosawa R., Proga D., Nagamine K., 2009, ApJ, 707, 823
 Li J., Ostriker J. P., & Sunyaev R., 2013, ApJ, 767, 105
 Luketic S., Proga D., Kallman T. R., et al., 2010, ApJ, 719, 515
 Magorrian J., Tremaine S., Richstone D., et al., 1998, AJ, 115, 2285
 Moe M., Arav N., Bautista M. A., Korista K. T., 2009, ApJ, 706, 525
 Merloni A., Rudnick G., Di Matteo T., 2004, MNRAS, 354, L37
 Miller J. M., Raymond J., Fabian A. C., et al., 2006, Nature, 441, 953
 Miller J. M., Raymond J., Reynolds C. S., et al., 2008, ApJ, 680, 1359
 Murray N., Chiang J., Grossman S. A., Voit G. M., 1995, ApJ, 451, 498
 Murray N., Quataert E., Thompson T. A., et al., 2005, ApJ, 618, 569
 Novak G. S., Ostriker J. P., Ciotti L., 2011, ApJ, 737, 26
 Novak G. S., Ostriker J. P., Ciotti L., 2012, ApJ, 427, 2734
 Ohsuga K., Mori M., Nakamoto T., Mineshige S., 2005, ApJ, 628, 368
 Ohsuga K., Mineshige S., 2011, ApJ, 736, 2
 Ohsuga K., Mineshige S., Mori M., Kato Y., 2009, PASJ, 61, L7
 Ostriker J. P., Choi E., Ciotti L., Novak G. S., Proga D., 2010, ApJ, 722, 642
 Proga D., Stone J. M., Kallman T. R., 2000, ApJ, 543, 686 (PSK00)
 Proga D., Stone J. M., Kallman T. R., 2004, ApJ, 616, 688
 Proga D., 2007a, ApJ, 661, 693
 Proga D., 2007b, in Ho L. C., Wang J.M., eds, ASP COnf. Ser. Vol 373, The Central Engine of Active Galactic Nuclei. Astron. Soc. Pac., San Francisco, p. 267
 Proga D., Ostriker J. P., Kurosawa R., 2008, ApJ, 676, 101
 Roth N., Kasen D., Hopkins P. F., Quataert E., 2012, ApJ, 759, 36
 Sazonov S. Y., Ostriker J. P., Sunyaev R. A., 2004, MNRAS, 347, 144
 Sazonov S. Y., Ostriker J. P., Ciotti L., Sunyaev R. A., 2005, MNRAS, 358, 168
 Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
 Silk J., Rees M. J., 1998, A&A, 331, L1
 Springel V., Di Matteo T., Hernquist L., et al., 2005a, ApJ, 620, L79
 Springel V., Di Matteo T., Hernquist L., et al., 2005b, MNRAS, 361, 776
 Stone J. M., Norman M. L., 1992a, ApJS, 80, 753
 Stone J. M., Pringle J. E., Begelman M. C. 1999, MNRAS, 310, 1002
 Steinhardt C. L., Elvis M., 2010, MNRAS, 402, 2637
 Somerville R. S., Hopkins P. F., Cox T. J., et al., 2008, MNRAS, 391, 481
 Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Yaqoob T., Braito V., Dadina, M., 2010, A&A, 521, A57
 Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Braito V., Dadina M., 2011, ApJ, 742, 24
 Tombesi F., Cappi M., Reeves J. N., Braito V., 2012, MNRAS, 422, L1
 Watarai K. Y., Ohsuga K., Takahashi R., Fukue J., 2005, PASJ, 57, 513
 Yuan F., Bu D. F., Wu M. C., 2012, ApJ, 761, 130
 Yuan F., Li M., 2011, ApJ, 737, 23
 Yuan F., Xie F., Ostriker J. P., 2009, ApJ, 691, 98
 Yuan F., Wu M. C., Bu D. F., 2012, ApJ, 761, 129