Parsec-scale accretion and winds irradiated by a quasar
We present numerical simulations of properties of a parsec-scale torus exposed to illumination by the central black hole in an active galaxy (AGN). Our physical model allows to investigate the balance between the formation of winds and accretion simultaneously. Radiation-driven winds are allowed by taking into account radiation pressure due to UV and IR radiation along with X-ray heating and dust sublimation. Accretion is allowed through angular momentum transport and the solution of the equations of radiation hydrodynamics. Our methods adopt flux-limited diffusion radiation-hydrodynamics for the dusty, infrared pressure driven part of the flow, along with X-ray heating and cooling. Angular momentum transport in the accreting part of the flow is modeled using effective viscosity. Our results demonstrate that radiation pressure on dust can play an important role in shaping AGN obscuration. For example, when the luminosity illuminating the torus exceeds , where is the Eddington luminosity, we find no episodes of sustained disk accretion because radiation pressure does not allow a disk to form. Despite the absence of the disk accretion, the flow of gas to smaller radii still proceeds at a rate through the capturing of the gas from the hot evaporative flow, thus providing a mechanism to deliver gas from a radiation-pressure dominated torus to the inner accretion disk. As increases, larger radiation input leads to larger torus aspect ratios and increased obscuration of the central black hole. We also find the important role of the X-ray heated gas in shaping of the obscuring torus.
Active galactic nuclei (AGN) owe their radiative output to accretion of gas which is supplied by the host galaxy. Many AGNs show evidence for obscuration of the radiation from the central black hole, and many show evidence for outflows or winds from the central regions. These phenomena likely arise in the region 1 pc from the center, where there is a transition between the predominantly cold and dusty galaxy interstellar medium and the viscous accretion flow onto the black hole. In this region, the dust from the galaxy can reprocess the strong radiation field from the center into the infrared (IR), and the relatively large IR opacity can efficiently couple the gas dynamics to the radiation field. Numerical models for dusty accretion flows require treatment of the coupled radiation and gas dynamics in addition to the effects of radiation on the dust and gas: dust destruction, radiative heating and cooling. In previous papers we have presented algorithms for treating these processes using simplified assumptions about the origin of the IR radiation field and also assuming a source of gas on the domain boundary. In this paper we present models which self-consistently treat the reprocessing of the radiation from the center into IR, and which have no source of gas outside the computational domain. These models demonstrate that geometrically thick obscuration, outflows, and inflow to the center can all be produced in a radiation-pressure dominated region near 1 pc.
The specific goal of our models include:
Our previous work supports the idea that obscuration in luminous AGN can originate from dusty IR-driven winds. We explore this idea further without artificially separating between accretion disk and wind.
We calculate mass-loss rate and accretion rate. These are time-dependent and can help to understand the fate of episodic accretion events in luminous AGNs.
Further developing a physical model of the pc-scale AGN yields such properties of the wind as covering factors, column densities, outflow rates, kinetic luminosities as functions of observing angle, and their dependence on accretion rate and luminosity. These can provide input for simulations on the larger galactic scale and for models for feedback.
The plan of this paper is as follows: A discussion of previous work by ourselves and others is presented in Section 2. Also included are simple estimates demonstrating the potential importance of IR radiation radiation pressure; this section can be skipped by the expert reader. In Section 3 equations of radiation hydrodynamics used to describe the problem are reviewed together with adopted approximations, along with the adopted approximations for angular momentum transport and interaction of the radiation with the matter. Numerical methods and approximations used to solve the system of radiation hydrodynamic equations are summarized in Section 4 and the results are presented in Section 5.
In the AGN region close to the black hole accreting material likely forms an accretion disk in which angular momentum is transferred via the instability associated with the distortion of magnetic fields in a conductive shearing flow (Balbus & Hawley, 1991). This occurs within cm, where , is the mass of the black hole, and subscripts indicate scaled quantities, i.e. . The inner disk region is responsible for most of the energy generation in AGN but in order to reach the inner disk, gas must cross many decades in radius from regions associated with the host galaxy interstellar medium. In order to do so, the gas must lose most of its angular momentum; the details of how this occurs on large scales are likely to be very different from the behavior in the inner disk itself. This is due to both the properties of the cooler gas associated with the host galaxy, notably the effects of dust, and also to the effect of preheating by the strong radiation from the inner disk. In this paper, we present models for preheated dusty accretion flows and examine the behavior of the material as a function of the intensity of the illumination from the inner disk.
Radiative preheating can heat accreting material to temperatures greater than the local virial temperature. For gas pressure dominated material, this temperature is K, where is the distance in parsecs. If so, the heated material can be levitated into a quasi-static geometrically thick structure, or can flow outward in a wind. Such structures have implications for AGN unification, i.e. the hypothesis that the properties of both type I and type II AGNs can be explained if the two are intrinsically similar but are viewed from different angles (Rowan-Robinson, 1977; Antonucci, 1984; Antonucci & Miller, 1985). The crucial feature is an optically- and geometrically-thick obscuring structure which prevents a direct view of the inner disk and black hole from many directions (Antonucci & Miller, 1985). Outflows also can couple the properties of the black hole, i.e. its luminosity, with the properties of the host galaxy, since the outflows can deposit significant energy on galaxy scales. Also outflows can limit or regulate the net accretion rate onto the black hole. The importance of radiation as a regulating agent has been explored in by, e.g. Novak et al. (2012); Proga et al. (2008); Proga (2007) and its connection to feedback has been demonstrated by (e.g. Silk & Rees, 1998; King, 2005; Ostriker et al., 2010; Ciotti & Ostriker, 2007).
Winds can play a crucial role in the observational appearance of AGNs. Winds are produced in luminous AGN in the form of UV broad absorption lines (BALs)(Trump et al., 2006), UV warm absorbers (Crenshaw, 1997), and X-ray warm absorbers (Reeves et al., 2003; Blustin et al., 2003; Gibson et al., 2009; Gofford et al., 2011; Page et al., 2011). Evidence for outflows is also seen in the optical spectra from distant obscured quasars (Greene et al., 2012). These flows can be important to the AGN mass budget, and their rich spectra can provide detailed information about the outflowing gas: its ionization, temperature, speed, geometrical distribution, and elemental abundances. (Murray et al., 2005)
In dynamical models light from the central engine is obscured by the flows which are either inflowing (accretion flows) or outflowing (winds). Dynamical models for obscuration include: (i) Quasi-stationary models in which the obscuration is a torus supported by rotation in the radial direction. and by turbulent motion of molecular clouds in the vertical direction. Collisions between clouds would likely raise the temperature of the clouds to a fraction of the virial temperature K, which is too high to account for observed molecular and dust emission. Magnetic fields have been suggested as a mechanism to provide elasticity to the clouds and avoid dissipation (Krolik & Begelman, 1988; Beckert & Duschl, 2004; Nenkova et al., 2008). (ii) Various mechanisms, including irradiation (i.e. Foulkes et al., 2010), can potentially produce global warps in a locally geometrically thin accretion disk (Phinney, 1989; Sanders et al., 1989; Pringle, 1992). Such disks can potentially provide obscuration for a significant fraction of lines of sight. (iii) Geometrically thick flows can be associated with a global magnetic field which allows a magnetohydrodynamic (MHD) wind (Konigl & Kartje, 1994; Elitzur & Shlosman, 2006) or directly supports a quasi-static torus (Lovelace et al., 1998), or a combination of radiative and magnetic driving (Emmering et al., 1992; Everett, 2005; Keating et al., 2012; Konigl & Kartje, 1994).
Dust can play a key role in the properties of accreting or outflowing material near 1 pc. Dust survival in the presence of irradiation from the black hole and inner disk depends on the equilibrium temperature determined by balancing blackbody cooling with the radiative heating by radiation from the inner disk. This predicts that dust can survive outside of the sublimation radius where the dust sublimation temperature is K, , and is the radiative luminosity of the inner disk. Evidence for dust in AGN obscuration comes from observations of several nearby AGN where the nucleus is resolved by interferometric techniques. For example, in NGC 1068 such observations reveal a multi-component, multi-temperature dusty medium: an inner, relatively small ( pc) and hotter ( K) component embedded in a larger ( pc), cooler ( K) component (Jaffe et al., 2004; Raban et al., 2009). In the Circinus galaxy the observed elongated, 0.4 pc diameter component is interpreted as a disk-like structure seen almost edge on. The disk-like structure is coincident with that inferred from the VLBI maps of maser emission (Greenhill et al., 2003), being embedded in a much larger rounded component. This is interpreted as a geometrically thick torus with temperature K (Tristram et al., 2007). Dust affects the dynamics via its effects on radiation. Dust opacity, , in the infrared (IR) is typically (Semenov et al., 2003), where is the electron scattering opacity. The critical luminosity at which IR radiation becomes dynamically important is obtained by equating the dust radiation pressure per H atom with the gravitational force: , where Eddington luminosity, is the corresponding quantity for material dominated by the electron scattering opacity:
The high opacity of dust can effectively trap IR radiation, if it is generated within an optically thick medium. In the torus region, the energy density of trapped radiation can exceed the thermal energy density, and we can define a dusty virial temperature by equating the internal radiation energy density with the gravitational energy: , where is the gas number density. When the temperature of trapped IR exceeds this value, the material can be levitated, prevented from accreting, or it can flow out in a wind. This can produce a static torus (Krolik, 2007; Shi & Krolik, 2008) or an outflow if the radiation temperature in the torus exceeds .
AGN fueling depends on the action of angular momentum loss processes, and no single mechanism is likely to be dominant over the range of length scales from galactic (kpc) down to the black hole event horizon. On galactic scales, angular momentum can be transferred by long-range non-axisymmetric gravitational torques acting in the stellar and gas system(e.g. Lynden-Bell & Kalnajs, 1972). A particular implementation of this effect has been proposed by Shlosman et al. (1989, 1990), who suggested that from a typical galactic scale of 10 kpc down to 10 pc, gas accretes via a ”bar within bar” instability: successive bar-type instabilities in a gaseous, self-gravitating disk. Numerical simulations (e.g. Hopkins et al., 2012) give support for the idea that torques generated by non-axisymmetric self-gravitating instability are indeed responsible for the transport of the angular momentum from smaller to larger radii.
In a self-gravitating environment of the outskirts of an AGN, a similar effect can operate due to non-axisymmetric self-gravitating instabilities in a dense gas system. If the gas cools faster than the dynamical time-scale , where is the angular frequency, then the disc fragments into individual self-gravitating clumps. On the other hand, if these two time-scales are comparable, the disk settles into a gravito-turbulent state (i.e. Paczynski, 1978; Gammie, 2001). The non-linear development of self-gravitating instabilities provide dissipation and stresses for the angular momentum transport in a gravito-turbulent state . At radii pc we would generally expect magnetic field to be less important and self-gravity to be responsible for the angular momentum transport and accretion (Wada, 2012; Wada & Norman, 2001). At smaller radii the role of self-gravity diminishes, and magnetic stresses become responsible for angular momentum transport and accretion (Bisnovatyi-Kogan & Ruzmaikin, 1974; Narayan et al., 2003; Igumenshchev, 2008), including the regime where both mechanisms can be important (Fromang et al., 2004). The thick torus morphology of (Wada, 2012) was found by Schartmann et al. (2014) to be in a reasonably good agreement of the observations.
The stability of a geometrically, thin self-gravitating disk can be described by the Toomre parameter (Toomre, 1964) : where is the disk surface mass density, is the sound speed and is the epicyclic frequency, and is the specific angular momentum. It is usually assumed that a quasi-stationary, self-gravitating accretion disk hovers on the border of instability maintaining marginal stability and providing (Gammie, 2001; Goodman, 2003; Rafikov, 2009).
The dynamical model of AGN obscuration in which the torus consists of a radiation driven dusty outflow was explored in Dorodnitsyn et al. (2011, 2012); Dorodnitsyn & Kallman (2012) (hereafter Papers I-III, respectively). Numerical 2.5D (3D axially symmetric) radiation hydrodynamics simulations showed that the reprocessing of the UV and X-ray radiation into IR and the corresponding IR radiation pressure on dust grains can drive a strong wind. In Papers I and II we calculated the structure of the radiation-driven wind assuming that a certain fraction of the incident UV and X-ray radiation are converted into IR. Thus we did not follow the conversion of X-rays and UV into the IR. Instead, the IR radiation field was assumed to enter through the boundary of the computational domain. Flux-limited radiation-hydrodynamics was invoked to solve for the radiation-driven wind, originating from an accretion disk. The disk at the equator was treated as a boundary condition, serving as a source of matter for the wind. In Paper III we calculated the conversion of external X-rays into IR, preserving the equatorial symmetry and the disk boundary conditions. Mass loss rates were found to be in the range . This can be compared with the Eddington accretion rate, calculated taking into account electron scattering: , which suggests that it is difficult or impossible to sustain accretion when the gas is affected by such intense radiation pressure. Accretion with such intense winds should be intermittent with strongly varying accretion rate.
3.1 Equations of Gas Dynamics and Radiation Transport
Our physical model describes the time-evolution of a three-dimensional distribution of gas and dust in the gravitational field of a supermassive BH, adopting radiation hydrodynamics in axial symmetry. Evolution of dusty gas exposed to IR radiation is described by equations of radiation hydrodynamics. To first order in , these can be written (e.g. Mihalas & Mihalas, 1984):
where is the convective derivative, is the mass density and is the number density of the gas; is the velocity, , and are gas pressure and gas energy density, is the gas temperature. We assume a polytropic equation of state: for the gas, and ; throughout this work we assume .
Radiation input from X-ray and UV illumination is taken into account by the term in the equation (3), where and are respectively heating and cooling functions. IR radiation field is described by the radiation energy density, ; radiation flux, ; and by the radiation pressure tensor, ; opacities: , are the Planck mean and energy mean absorption opacities. Opacities as well as other dependent variables are evaluated in the co-moving frame of reference; other radiation-related notation include the Planck function , where is the Stefan-Boltzmann constant and is the speed of light; denotes the contraction , where . The validity of equations (1)-(4) in different regimes of the radiation-matter interaction is discussed in Paper II.
The Radiation pressure combines pressure of continuum IR radiation on dust, and the UV pressure in spectral lines and the radiation pressure due to Thomson scattering on free electrons. Gravitational forces are assumed to be only due to a supermassive black hole, . Self-gravity of the gas is not taken into account. Angular momentum transport is allowed by an anomalous viscous stress-tensor, , and viscous dissipation taken into account via a dissipation function .
Following previous studies of gaseous, self-gravitating disks (Lin & Pringle, 1987) we describe the mechanism of angular momentum transport through the introduction of the anomalous viscous stress tensor. The anomalous stress tensor Cartesian coordinates reads
where , and is the effective kinematic viscosity.
The viscous force is , where the summation over a set of similarly indexed terms is implied. Notice that in case of a geometrically thin equatorial part of the accretion disc, the viscous force simplifies to . The viscous heating function, that appears in (3) reads
The effective viscosity,
depends on the critical spatial scale, which defines the size of the unstable perturbation. The latter can be stabilized by the shear, giving , or the perturbation is stabilized by the effective pressure, i.e. the velocity dispersion, giving the Jeans scale , where is the local dispersion velocity, Lin & Pringle (1987) argue that ; here, for simplicity we assume that . When performing numerical simulations (see further in the text) we experimented with and did not find any noticeable difference in the results.
3.2 Photoionization equilibrium and X-ray heating and cooling
The radiation field in our calculation is divided into the incident continuum from the inner disk and black hole, and the diffuse IR produced by thermal emission within the accretion flow. The incident continuum is treated using a single stream approximation. We trace X-rays from a central point source and assume attenuation according to , where the optical depth with respect to X-ray absorption is , and is the X-ray opacity.
We assume that the ionization state of the gas can be parametrized by ionization parameter, :
where is the column density in ; is the local X-ray flux; X-ray luminosity of the nucleus: , where is the total accretion luminosity. This relation also serves as a definition of the model parameter where is the Eddington luminosity, and is the mass of a BH in .
Energy deposition by external X-rays is calculated making use of heating and cooling functions. These include Compton heating and cooling:
photo-ionization heating-recombination cooling:
(and ), bremsstrahlung and line cooling:
We take the parameter representing optically thin line cooling.
Heating and cooling rates , and are calculated making use of the XSTAR photo-ionization code (Kallman & Bautista, 2001) assuming ionizing continuum with (energy) power law index of =1. Results were approximated to a reasonable () accuracy by the analytical formulae, originally derived by Blondin (1994) for a 10 keV bremsstrahlung spectrum with K. Formulae (9)-(11) are the slightly improved version by Dorodnitsyn et al. (2008) which incorporate newer atomic data.
3.3 Cold gas and dust
Interaction of X-ray radiation with the cold molecular and dusty gas depends on various physical and chemical processes (Maloney et al., 1996; Krolik & Lepp, 1989). In general, the proper incorporation of this material calls for consideration of all the radiative, chemical and dust network of processes. To couple this to the radiation hydrodynamics we adopt a simplified treatment which nevertheless preserves the qualitative behavior of gas and dust in the presence of an intense radiation field. We separate the interaction into different regimes depending on whether the gas is highly ionized, primarily atomic, or largely molecular. To decide between different regimes we check if an effective ionization parameter meets certain threshold requirements. Thus, following (Maloney et al., 1996), cooling of the cold gas depends on the value of the effective ionization parameter, which is a switch that determines the cooling channel and whether the gas is assumed to be primarily in atomic or molecular state. For details of how we treat this problem, refer to Paper III.
Dust and gas coupling
Throughout this work we treat dust and gas as one fluid. This is a simplification that implies both momentum and positional coupling. Momentum coupling means that the momentum of the accelerated dust grain is effectively shared with the gas molecules and ions; the positional coupling means that the grain doesn’t travel far alone with respect to the gas. On the other hand, momentum coupling and no positioning coupling are the usual assumptions in the calculations of dusty radiation-driven winds, such as winds from AGB stars (Gilman, 1972; Lamers & Cassinelli, 1999). The conditions, most notably, the characteristic length scales, in our wind are different from those of a wind from a star. In the following and throughout our calculations, we adopt the following fiducial values for the parameters of the dust: the grain radius, , the density of the dust grain , and the dust-to-gas mass ratio . Above the sublimation temperature, the dust is assumed to be destroyed. The assumption of dust and gas coupling is justified for the conditions relevant to our simulations, as can be shown via simple estimates.
The analysis, given by Krumholz & Thompson (2013) for the dust grain collisional stopping length, gives pc (see also Murray et al., 2005). This length scale is a factor smaller than the smallest length scale in our simulations. Sticking of impinging electrons and ions and interaction with intense X-ray and UV radiation makes the grains to possess a charge which further couples them to gas (Draine & Salpeter, 1979). If there is no significant magnetic field, then the interaction between grains and ions is governed by the grain stopping time which can be estimated taking into account the fact that dust grains are much heavier than ions. If so Spitzer (1978) gives yrs where is the thermal velocity of gas. Comparing it with the sound crossing time we have: .
A different estimate comes from that AGN torus region should have small scale magnetic fields (i.e. see the discussion in Krolik & Begelman (1988)). On the microphysical level, the gas at the torus region should be turbulent, with small-scale magnetic field of the order of the equipartition level. In this case the grain stopping time is defined by , where is the gyroradius. The dust grain that possesses a charge, will gyrate in the magnetic field with the gyration radius: . The small-scale equipartition magnetic field thus gives pc for . Thus the gyroradius is small compared with the size of the region, and compared with our numerical resolution. Dust particles must move with the small scale field, and hence with the entire plasma, on length scales .
These estimates show that and the dust-grain mixture can be treated as a single fluid. The assumption of the positional coupling can break at high altitude, i.e. close to the torus funnel where the density is very small. However, the intense X-ray radiation of AGN makes an important difference between the torus winds and a dusty stellar wind. Our numerical results show, a posteriori, that no dynamically significant amounts of dust survive in the region of low density because of the strong X-ray radiation there. Though it would be interesting to consider the effects of dust separation at high altitudes this is beyond the scope of the current work.
Dust and molecular phase
If the ionization parameter is small the gas is primarily molecular. When , we assume radiative equilibrium between molecular gas and dust, i.e. . The equilibrium dust temperature is found assuming that dust can directly reprocess X-rays to IR (see Paper III). If and we assume that the exchange between the radiation temperature, and is determined from equations (3)-(4). The gas temperature can contribute to IR and vice versa. In equation (4) we assume that , where is the infrared flux. The details of how there radiation pressure is calculated are given in Section 3.4. For we assume that the dust is destroyed. The opacity of the gas, in this case, is taken to be that of Thomson.
3.4 IR radiation pressure
Radiation pressure consists of two components: i) pressure of the continuum IR radiation on dust, ii) Another component of the radiation pressure - the radiation pressure in lines, is due to interaction of the UV radiation with spectral lines (Section 3.5). The total radiation pressure:
The continuum radiation pressure: where is the radiation pressure due to Thomson scattering and is the local radiation flux.
In order to make the problem tractable in the framework of radiation hydrodynamics we make several simplifications: We assume that the UV flux contributes only to the radiation pressure in spectral lines (see next paragraph). Similarly, X-ray flux contributes to X-ray heating; Both fluxes are attenuated through a simple, prescription. In all our models we assume .
The IR radiation pressure is given by:
where the total flux mean opacity, consists of absorption opacity and the Thomson scattering opacity. For simplicity, we do not differentiate between , and .
If optical depth is large , the diffusion coefficient is: . where is the photon mean free path, and is the dust opacity per unit mass. Notice, that in the diffusion regime cancels out from the formula for .
That is if the optical depth of the gas-dust mixture in the IR , radiation pressure where is the the radiation temperature. When radiation pressure scales as .
The diffusion approximation tacitly assumes that optical depth . Most of the torus where IR pressure is important indeed has . When the mean free path, , and , and . In a free-streaming limit however, there should be . That is, when optical depth becomes small, or when , the diffusion approximation is no longer applicable. Thus, at it must be modified to obtain the correct limiting behavior in the radiation free-streaming limit.
Thus, in order to describe correctly regions of small , we adopt the flux-limited diffusion approximation (Alme & Wilson, 1974; Minerbo, 1978; Levermore & Pomraning, 1981). In this approximation is replaced by , where is the flux limiter. The flux limiter we adopted in Paper I,II and in the current work is that of Levermore & Pomraning (1981):
where . One can see, that if , then , and ; and if and .
3.5 Radiation pressure in spectral lines
Where the gas is hot and partially ionized, we take into account the pressure of UV radiation in spectral lines:
where is the optical depth parameter, is the line strength parameter (different from the accretion efficiency used earlier), is the Thomson cross-section, is a proton’s mass, and ; a parameter was introduced to limit the effect of very strong lines (Owocki et al., 1988; Stevens & Kallman, 1990). Local UV flux is calculated from the total bolometric luminosity adopting a simple attenuation model and assuming the nucleus radiates in UV. The parameter in (17) reads (Stevens & Kallman, 1990):
where is found from
As it follows from the results of our numerical simulations, the large opacity of dusty gas with respect to UV makes component of the acceleration negligible in the dusty component of the flow.
4 Numerical Methods
The system of radiation-hydrodynamics equations (1)-(4) is solved with our radiation-hydrodynamics code. The hydrodynamic part of this system is solved with the methods implemented in the ZEUS-MP (Hayes et al., 2006) version of the original ZEUS code (Stone & Norman, 1992). The radiation-hydrodynamic part is designed and tested in Papers I-III and conforms with methods and structure of the original ZEUS code (Stone & Norman, 1992).
We calculate the dynamics throughout the entire volume occupied by the accreting and outflowing gas. In this way, we address a variety of possible ways in which an accretion inflow and the wind can co-exist. The addition of this new capability to treat accretion requires to compute the term in the momentum equation, (2), and the viscous heating term in the energy equation, (3). During the time-step update this is done adopting an operator-splitting technique.
We solve the time-dependent system of equations (1)-(4) in cylindrical coordinates , adopting a uniform cylindrical grid, with the . The grid extends from approximately the dust sublimation radius, pc to pc in the radial, and from pc to pc in the vertical direction. Three components of the velocity, , , are taken into account, assuming azimuthal symmetry, , i.e. adopting the approximation.
Boundary conditions (BC) for , , at the boundaries of the computational domain are set to be outflowing. The radiative BC are free-streaming at all boundaries. Radiation BC are implemented in a way that preserves second order accuracy, required for proper interoperability with other numerical methods in our original diffusion solver.
Radiation hydrodynamics is a highly non-linear problem, and obtaining the numerical solution of the system (1)-(4) is challenging. For example, gas temperature depends on X-ray heating which depends on density. Temperature gradients create strong local radiation pressure, providing a feedback loop to dynamics (see Paper II for the description of the radiation time-scales).
The solution of coupled equations, (3)-(4) for and requires the solution of a general fourth order algebraic equation (Paper II). This is very sensitive to the range of parameters of this equation such as time step, etc. To make things worse, a strong non-linearity affects many parts of our system: X-ray heating-matter interaction; coupling between thus obtained gas temperature and radiation energy equation; finally, in the momentum equation there is a strong dependence on the gradients of the radiation energy density, and again, implicitly on the rate of energy deposition from X-ray heating. To tackle this problem(s) we often have to split the time-step and then separately sub-cycle the radiation hydrodynamics part of equations in time until the desired time-step is achieved. A similar line of arguments but perhaps less restrictive applies to coupling between the dynamics and the effects of viscosity.
4.1 Initial conditions
Simulations start from a rotating, polytropic, i.e. , constant angular momentum torus that is imbedded in low-density gas. An analytical solution that describes such a torus does exist, (Papaloizou & Pringle, 1984) but does not take into account radiation pressure. We calculate the initial distribution of the density from an approximate formula which in the non-radiative case reads:
where is the distance to the locus of the density maximum, ; the parameter measures the distortion of the torus’s polar cross-section. Notice that in the radiation-dominated case . To calculate this initial distribution we do the following: we calculate from (19) assuming a polytropic torus with ; substitute in (19) total pressure , where is the universal gas constant. Thus, we assume that the gas and radiation are initially in equilibrium and calculating the distribution of ; in the rest of the simulations we take .
4.2 Input Parameters
Radiation energy input into the flow is controlled by the major input parameter . We present several models for which illustrate the range of behaviors spanned by interesting values of . Local radiation flux depends on the attenuation and thus on the density scale, of the initial density distribution, . All models have which corresponds to the initial torus mass, . The boundaries of the computational domain are located at and and remain fixed for all models. Other parameters are , , and . The evolution of the radiative flow was followed for dynamical times, .
5.1 General effect of radiation on a torus
Radiation influences gas dynamics in our models in several ways:
X-ray irradiation of the torus creates an overheated thin layer - a ’skin’ where the ionization parameter and X-rays heat the gas to a temperature comparable to the local gas virial temperature, so that it expands into the inner cone-shaped throat or funnel at small radius. Eventually, this hot gas completely fills the throat of the torus. Some of this gas is accreted i.e. crosses the left boundary of the computational domain, and some flows to larger radii as a thermal wind. Confinement by the cold and dense torus creates a funnel flow - a fast outflow of a hot and diluted gas within an approximately conical region surrounding the disk axis.
UV and X-rays are absorbed in the interior of the torus, and reradiated as thermal IR at the local equilibrium temperature. This may increase the internal torus radiation temperature to or above the local radiation virial temperature. Then the torus expands and the outflow commences.
Radiation deposited in the inner disk also exerts a net bulk radiation pressure on the torus, pushing it outward.
Viscous torques redistribute angular momentum within the torus allowing the gas to move inward, i.e. to accrete. This process is enhanced by radiative heating of the disk interior.
Back-reaction from the hot gas on the torus. Hot gas created in step (i) pushes the torus further outward ablating and compressing it in the vertical direction.
These processes are illustrated in the following results. We present four models, which differ by their values of the Eddington ratio . Dynamical and obscuring properties of this models are described in Sections 5.6, 5.7, and 5.9.
5.2 Model with
At low Eddington ratio the effects of internal radiation pressure are small compared with the other radiation influences described above. Figure 1 shows the density map for the model with . This model is characterized by two distinct components: a hot evaporative outflow, and a colder thin accretion disk. Radiation heats the gas in the inner torus ”funnel” through Compton and photo-ionization heating, raising the temperature to K. This occurs within a Compton heating time , where is the Compton energy transfer per scattering. This is fast compared with the dynamical time associated with gravitational motion. Thermal expansion times are also rapid compared with the dynamical time so that after yrs, hot gas evaporated from the torus fills approximately half the domain, This can be seen at the earliest time-step shown in Figure 1 and corresponds to step (i) above. Temperatures within this flow range from to a fraction of the Compton temperature, K. An outflow occurs mostly in a form of evaporative flow with K. This gas leaves the computational domain as a wind. For example, the maximum velocities at the and directions: , at pc and at pc.
The influence of infrared radiation pressure on the flow depends on the value of via the dust effective temperature. At low , the infrared radiation input into the interior of the disk is not strong enough to generate any significant radiation pressure on dust. This can be seen from Figure 1 at times greater than yrs: the cold core of the thin disk is not disrupted rapidly. Consequently, there is no infrared-driven wind formed. Inside the torus, viscous torques redistribute angular momentum of the gas allowing the gas to move inward, i.e. to accrete. The important characteristic of this model is that it shows a significant period of disk accretion. This is apparent at time yrs; a tail of higher angular momentum material moving out at large radii is a signature of accretion. This structure is also compressed by wind material.
At later times, yrs evaporation from the disk significantly depletes the disk mass. This reduces the net inflow, so that evaporation from the inner disk edge ultimately stops the accretion of disk material. The disk is truncated and hot, much more compact but still geometrically thick torus which appears roughly at the same place where the initial torus was located. The mass of this torus is much less than the initial mass .
5.3 Model with
At larger Eddington ratio radiation pressure plays a stronger dynamical role. This is illustrated in Figure 2 which shows a sequence of time snapshots from the evolution of the distribution of the density for .
Comparison of Figure 2 with the previous figure shows that the increased illumination decreases the timescales for all radiation-related processes. The evaporative flow is established earlier and is stronger, as is the later evaporation of the disk. There is no net inflow of the cold thin disk core; rather there is net outflow of this material starting from the earliest times shown. This is due in part to the effects of infrared radiation pressure, which disperses the disk and prevents the establishment of a viscous-dominated inflow, and also to evaporation owing to the penetration of X-rays into the disk interior. Figure 2 shows that, at early times, the dispersal of the disk is strongly affected by radial flow in the disk plane. This is due to IR radiation pressure which tends to act in the direction opposite to maximum temperature gradient. Owing to the presence of warm or hot (K) gas both above and below the disk plane, the maximum gradient is in the radial direction.
The maximum velocity is observed in the hot flow where it approaches in the radial, and in the vertical direction (Figure 2, upper left). At later times, in all three other time-steps, the velocities are in the range and . Following Paper III we calculate the average bulk velocity of the flow , where is the total or partial volume occupied by the flow. The analysis indicates that on average the most massive wind is relatively slow: for the first snapshot shown, for the second, and and for the rest. These velocities are close to quasi-stationary values, and in accord with the evolution of morphology. This shows similar shapes for the torus at the last three time snapshots.
The gas temperature in the =0.05 model is shown in Figure 3, at the time 6.8 yrs corresponding to the first panel in Figure 2. Also shown are velocity streamlines and the region where dust survives. One can see that the cold gas is mostly accumulated within an equatorial disk. The streamlines illustrate the radial nature of the IR-driven part of the flow, originating in the torus interior. The effects of wind-compression of the disk at large radii 2 pc are also apparent. A shock-like structure at 1.8 pc is associated with the radial flow encountering the boundary of the initial torus. The apparently geometrically thick part of such disk at pc is mostly a relic left from the original torus that has been compressed by the hot gas. In its thinnest part, the equatorial disk, shown in Figure 3 and Figure 2 is under-resolved, spanning only grid cells in height.
The dust distribution for the configuration in Figure 3 coincides with the distribution of the cold gas shown in blue color. As we noted before, in our simulations we do not have dust as a separate component, but rather take into account the enhanced opacity due to the presence of dust. That is done whenever . The corresponding region where dust survives is shown in Figure 3 inside the dashed contour.
At late times the disk at distances greater than 1.2 pc is dispersed into a hot cloud which has a small net speed, in either the vertical or the radial direction. This gas resembles a failed wind, in which cooling or weak illumination prevent the gas from attaining a temperature capable of escaping the system.
5.4 Models with : evaporative wind
The behavior of models with larger luminosity resemble that of the model except that step (ii) of the outline, i.e. the effect of internal IR from absorbed X-rays, is much more rapid and apparent at short times. In Papers II and III we found that when luminosity approaches radiation pressure can disrupt or reverse the accretion flow. If so, a key assumption used in those papers, the existence of a thin accretion disk situated at the equator which is unaffected by preheating, may be invalid; here we can test whether similar results are found without requiring a thin disk boundary condition. Figure 4 shows the evolution of the torus illuminated by from the onset of the simulations (upper left) towards a complete disruption of the torus into filaments and clouds (lower right). The torus is continuously evaporating, replenishing the warm gas. The maximum velocity is in the range for the radial velocity, and in for the vertical one. This approximately equals to , where is the escape velocity at 1pc.
Except for the first time snapshot shown, . The phase of an intense radiation-driven wind lasts for about yrs. Compared with the results for , the strong radial flow in the disk plane is established much earlier and transports material more rapidly. At this stage, . As the wind becomes progressively more thermally-driven, the mass-loading declines and so does velocity, to An equatorial disk-like outflow persists for approximately another yrs. This outflow is clearly seen in Figure 5, where polar distributions of and are shown. In this figure, an accretion flow can also be recognized as a region of negative radial velocity. In this disk where , both radiation pressure and gas pressure are important, and the majority of the radial radiation-driven flow is formed. The left panel of Figure 5 also shows a fast thermal wind.
At later times ( yrs), the part of the wind that has been shielded behind the dense torus evolves into a disrupted network of over-dense filaments. This further evolves into a clumpy structure with total mass .The results of this simulation are in accord with results from Paper I: at the external illumination significantly affects the flow to the extent that the interior of the viscous disk is affected. Radiation pressure on dust can either prevent accretion or it can reduce the accretion rate by squeezing the inflow into a self-shielding, thin equatorial disk.
5.5 Models with inflow-outflow accretion.
Figure 6 shows a series of time snapshots from a simulation for . Time steps shown are the same as in the three previous models. The strong illumination in this model removes much of the complexity seen in the previous models. At yr a complex tail is seen in the plots, at large radii outside the torus. This is likely a result of thermal instability. The disk behind the denser torus has a layered, sandwich-like structure: cooler and denser layers are intermittent with less dense and hotter ones. The dynamical structure is similar to a previous model, depicted in Figure 5. At late times ( yrs) there is no evidence for moderate temperature (K) gas outside the torus; the radiation is strong enough to drive all evaporated material to the Compton temperature and thus to escape. At these times the torus possesses a fairly constant evaporative wind. As expected, this configuration evaporates faster compared to less luminous cases: within a few dynamical times the initial torus is stripped of most of the gas and the only densest core of the original torus is left evaporating at the rate (Section 5.7).
5.6 Mass-Accretion rate
The results in the previous sections show that, when radiation heating is strong (i.e. ), gas that was initially in the disk is dispersed into a hot wind and pushed outward by radiation pressure at a rate which is faster than the viscosity can diffuse gas inward toward the black hole. If so, viscous accretion stops. Figure 7 shows time-dependent properties of the accretion rate, , i.e. the mass flux through the inner boundary of our computational domain.
The model with the lowest radiation input, i.e., demonstrates a period when viscous accretion disk occurs at the equator (c.f. Figure 1, lower left panel), which corresponds to a rise of the accretion rate by four orders of magnitude within a very short period of time near 2 yrs. The disk thickness can be estimated, assuming radiation pressure support: , which is in relatively good agreement with Figure 1. Other studies of large scale accretion including irradiation by the central AGN found similar behavior: despite strong illumination accretion proceeds in the equatorial region because the radiation flux is significantly reduced due to the geometrical foreshortening effect and because the equatorial flow is optically thick (Kurosawa & Proga, 2009). Until accretion proceeds at the level of , where is the critical ”Eddington” accretion rate, for simplicity taken to be that of a non-rotating BH: . After peaking at at the accretion rate declines and eventually stops. This is due to truncation of the accretion disk at . That is, due to rapid accretion the disk depletes and the reduced density cannot effectively support screening of its inner parts from overheating and evaporation. After some time, hot gas fills the throat of the torus and accretion resumes at much lower level through the capturing of the gas from the wind.
Accretion is interrupted for a period of , after which it resumes at the level of . This second episode of accretion is happening via a different mechanism of accretion: through capturing of the gas from the wind. Note that what we call the is not the actual rate at which the gas is crossing the innermost stable orbit of the BH. We only register the gas when it is crossing the boundary of the computational domain. The potential energy of the accreted matter will be released in a form of radiation - something that we don’t address here as it happens much closer to the BH. This will correspondingly increase .
At , the character of the flow changes from net accretion to net outflow. That is, our models with do not show disk accretion existing for significant amounts of time. It is possible that more dense models (i.e. higher values of ) will allow for the disk to form even at levels of greater than . At this level of radiative input, radiation heating and radiation pressure hampers the formation of the accretion flow. In a model it takes for the initial configuration to adjust and to achieve quasi-stationary accretion. All models except the one with accrete at less than the Eddington rate, . The model with has the weakest evaporative wind and thus the lowest accretion rate: , while the model with has the highest: .
5.7 Mass-loss rate of the wind
The wind mass-loss rate, is calculated at the upper, bottom and at the outer cylindrical boundary. Figure 8 (left) shows the time-dependent properties of this quantity. If the entire torus is expelled within dynamical times, with no material accreted, then on average the wind mass loss rate satisfies . Mass-loss rates observed in the simulations are typically much smaller: In all runs the core of the initial torus that contains most of mass survives complete evaporation (c.f. Figure 6).
From comparing Figure 7 with Figure 8 (left), the wind mass-loss rate mirrors the accretion rate. In the model as the inner edge of an accretion disk progresses towards the radiation source, the mass loading from the evaporation increases and produces an increase in . Similar behavior can be noticed for other models.
Models with higher have correspondingly higher : For after peaking at , settles at the level of . Models with higher values of peak at gradually reducing to lower values at later times. The peak is associated with the excretion disk and the corresponding strong radial outward motion of the gas (Figure 5, right panel). Without the outer boundary, the mass-loss rate would be for , gradually increasing to for . After an initial period of approximately the mass loss rate settles to a quasi-stationary value.
As mentioned above, at radiation input is not strong enough to prevent the formation of an accretion disk (c.f. Figure 1). The accretion flow brings new material in proximity to the inner boundary of the computational domain where it can be evaporated by radiative heating. This is a likely explanation of a quicker recovery of after yr in a model.
5.8 Energy of the wind
To assess the dynamical importance of the wind, we calculate the kinetic energy fluxes at the outer boundaries of the computational domain. The resultant kinetic luminosity of the wind is , where the integration is done over the upper and lower boundary that is excluding the excretion disk. It is this relatively dilute gas that can be potentially detected as X-ray or UV absorbers. During initial period of relaxation, thermal output of the wind is , where is the thermal output and is the kinetic output. After yrs for all models. Within one orbital period the energy output rises to and stays remarkably constant at this this level until the end of the simulation: . However, at all times it remains much less than the total luminosity: . Comparing the two panels of Figure we see that the kinetic output of the model is actually the smallest. Since the average velocity of the flow is not large, km/s, the dramatic rise of mass loss rate in this model that is happening at is not mirrored in the kinetic energy output, c.f. Figure 8 (right). The energy that the wind carries away is quite low compared to the luminosity of the accretion disk. This is typical for stellar winds, driven by the radiation pressure (Lamers & Cassinelli, 1999), and was also found in simulations of outflows from accretion flows (Kurosawa et al., 2009).
5.9 Obscuring properties
One of key predictions from our models is the dependance the obscuring properties of the flow on the inclination angle. The apparent torus opening angle is directly relevant to the appearance of an AGN as type I or type II object. Our previous results suggest that at high radiation input not only radiation pressure but also pressure of the hot wind are dynamically important. The torus opening angle is determined by a balance between the gradient of the pressure of gas and radiation from the torus with pressure of the hot gas and radiation pressure from outside the torus. That is, as a low density, hot gas fills the throat the torus is squeezed towards the equatorial plane. This, together with the external radiation pressure tends to make the torus aspect ratio, smaller.
As the column density increases towards the equator we are interested in the critical angle, at which the wind becomes opaque. That is when , where is the Thomson optical depth, and is measured along a given line of sight, and the viewing angle, is measured from the vertical axis. Figure 9, shows the distribution of models with . That is for a given inclination angle, we calculate a total number of model time steps which are approximately opaque. Figure 10 show column densities of models at a given inclination. Colors correspond to different . Each point of a particular color indicates a model at some evolutionary time step. Model time span corresponding to the entire span of the time evolution of our simulations, are shown. The total number of models shown, .
When the episode of disk accretion, leads to a significant depletion of the gas, and thus this set of models has a lower proportion of obscuring models at any inclination. As increases, no more episodes of disk accretion are observed and larger radiation input leads to larger aspect ratios. The distributions of the column densities are broader at higher radiation input, reflecting the importance of radiation pressure. The results are that a model with has the smallest and a model with has the largest aspect ratio of the torus. The region of the torus which is truly Thomson thick subtends a polar angle at most (for ); lower column densities extend to .
In our previous work we adopted a pre-existing thin accretion disk as a boundary condition. This allowed us to study quasi-stationary accretion disk wind solutions but did not allow to address the question of whether accretion can be stopped by radiation pressure and heating effects. Together with these earlier findings, our current results suggest that, qualitatively, obscuration can proceed in accretion or outflow mode, depending on the radiative input. The first scenario comprises a global accretion flow that is slowed, puffed up and kept geometrically thick by IR radiation pressure. The second describes a radiation-driven obscuring outflow. There is also an intermediate mode in which equatorial accretion disk surrounded by slowly in-, or outflowing envelope. The quantitative boundaries between these scenarios depend on details of heating/cooling, dust formation and destruction.
Some of results presented here may be sensitive to the assumptions regarding microphysics or computational simplifications. A simple analytical model developed in Paper I of this series predicts that as the temperature of the dusty gas reaches K. the radiation pressure in the optically thick regime becomes larger than gravity and outflow begins. This is confirmed in the current studies. For example, the dynamics of the torus depends primarily on the energy density of infrared radiation and on the ratio of the dust opacity to the opacity of the fully ionized gas, , making the radiation pressure launching of the wind depend on how dust is treated. In this paper we adopt relatively simple assumptions about the dust formation and destruction. In later work we plan to examine in more detail the effects of the treatment of dust formation on the hydrodynamics.
It is found that gas accretion is observed throughout the full length of the simulations. However, the details are somewhat surprising: the main path for accretion at is the accretion from capturing of the thermally driven, evaporative flow. Accretion and outflows are both possible at . In Paper III we already cast some doubt on the applicability of the accretion disk as a boundary condition due to the finding that the wind mass loss exceed accretion rate needed to produce the required radiative input. Our current results suggest that pc-scale disk accretion is likely to be intermittent.
To calculate the X-ray flux we adopt a simple attenuation model neglecting other radiation transfer effects. Effectively, we adopt a single scattering approximation when calculating local X-ray heating of the gas. It has been shown (Roth et al., 2012) that multiple scattering of X-rays can be important (e.g. Sim et al., 2010; Higginbottom et al., 2014). However, no previous work has self-consistently included multiple-scattering effects with the dynamics. In this work, we have chosen to focus on the dynamics while neglecting multiple scattering.
A lack of the detailed radiation transfer can potentially have an important influence on the results. This is anticipated to improve in the future work. An additional approximation that we employ is that the frequency integral characterizing the coupling between the radiation field and the dust opacity is constant. As a result of reprocessing and expansion the IR radiation can be shifted to longer wavelengths in the regions of the flow at the largest radii. Though it would be interesting to investigate this effect but this is beyond the scope of this paper.
Other limitations are associated with adopted boundary and initial conditions. One approach is to allow matter to constantly enter the domain at the boundary (see e.g., Proga & Begelman, 2003; Mościbrodzka & Proga, 2008, 2013). However this introduces multiple degrees of freedom such as the distribution of matter at the boundary, initial angular momentum of gas, etc. In this work we start from a bound torus that is close to dynamical equilibrium and study how it reacts to external illumination. The limitation of this approach is that we are fundamentally limited by the initial reservoir of gas. We do not allow matter to enter, and all of the initial torus will be eventually depleted.
The exact characteristics of the accretion flow depend on the assumption of the effective viscosity. Limitations of numerical resolution make it unfeasible to capture the effects of angular momentum transport in the self-gravitating dusty gas self-consistently. We have neglected self-gravity and assumed axisymmetry. Our results suggest that at parsec scales the gas is self-gravitating and can develop strong over-dense filaments and clumps on dynamical time-scales. Long range gravitational torques from the filaments would lead to angular momentum transport. Heating of these filaments can provide a negative feedback and make the whole structure hover on the border of collapsing into stars, and maintain a Toomre parameter . In the denser part of the accretion disk, such as close to the equator, gravito-thermal instability may generate self- gravitating turbulence that provide an effective viscosity. We will explore these processes in future work.
We have shown that conversion of external UV and X-ray into IR radiation becomes important at luminosities . At characteristic luminosity the effect of radiation heating is profound: both direct X-ray heating and the radiation pressure from IR radiation that is converted from X-rays are altering the dynamics of the torus.
The temperature of the hot gas reaches K: such gas escapes in the form of a fast thermally-driven wind with characteristic velocity of ; the temperature in the optically thick dusty gas approaches K, and the radiation pressure of infrared photons on dust grains is capable of driving a massive wind with velocities of . The typical morphology consists of two kinds of winds: a thermal wind located closer to the BH and the IR-driven wind situated further away. An interface between these two winds typically includes a thin ”conversion” layer where X-rays are converted into infrared.
We find the possibility of three types of gas flows illuminated by intense radiation of AGN: i) A global accretion flow that is slowed, puffed up and kept geometrically thick by IR radiation pressure; ii) An equatorial accretion disk surrounded by slowly in-, or outflowing envelope; iii) A radiation-driven obscuring outflow. At late stages of evolution, the interaction of radiation with the outflowing dense gas creates a complex structure consisting of compression waves, filaments and clouds. We hypothesize that strongly compressed gas and the web of shocked filaments and clumps that is found in our simulations can be possible sites for maser emission. More work is required to study the details of this regime.
Large amounts of hot, photo-ionized, and dilute gas are routinely observed in our simulations. Ram pressure of this component is important in shaping of the cold flow. Both radiation pressure and the dynamical pressure of the hot flow contributes to the ”receding torus” effect at high luminosities. The spectroscopic detection of this gas can be an important task for future missions such as Astro-H and Athena.
- affiliation: X-ray Astrophysics Laboratory, Code 662, NASA Goddard Space Flight Center, Greenbelt, MD, 20771, USA
- affiliation: Department of Astronomy/CRESST, University of Maryland, College Park, MD 20742, USA
- affiliation: Space Research Institute, Profsoyuznaya st., 84/32, 117997, Moscow, Russia
- affiliation: X-ray Astrophysics Laboratory, Code 662, NASA Goddard Space Flight Center, Greenbelt, MD, 20771, USA
- affiliation: Department of Physics and Astronomy, University of Nevada, Las Vegas, NV 89154, USA
- Alme, M. L., & Wilson, J. R. 1974, ApJ, 194, 147
- Antonucci, R. R. J. 1984, ApJ, 278, 499
- Antonucci, R. R. J., & Miller, J. S. 1985, ApJ, 297, 621
- Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
- Beckert, T., & Duschl, W. J. 2004, A&A, 426, 445
- Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45
- Blondin, J. M. 1994, ApJ, 435, 756
- Blustin, A. J., et al. 2003, A&A, 403, 481
- Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157
- Ciotti, L., & Ostriker, J. P. 2007, ApJ, 665, 1038
- Crenshaw, D. M. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 128, Mass Ejection from Active Galactic Nuclei, ed. N. Arav, I. Shlosman, & R. J. Weymann, 121
- Dorodnitsyn, A., Bisnovatyi-Kogan, G. S., & Kallman, T. 2011, ApJ, 741, 29
- Dorodnitsyn, A., & Kallman, T. 2012, ApJ, 761, 70
- Dorodnitsyn, A., Kallman, T., & Bisnovatyi-Kogan, G. S. 2012, ApJ, 747, 8
- Dorodnitsyn, A., Kallman, T., & Proga, D. 2008, ApJ, 687, 97
- Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 77
- Elitzur, M., & Shlosman, I. 2006, ApJ, 648, L101
- Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992, ApJ, 385, 460
- Everett, J. E. 2005, ApJ, 631, 689
- Foulkes, S. B., Haswell, C. A., & Murray, J. R. 2010, MNRAS, 401, 1275
- Fromang, S., Balbus, S. A., Terquem, C., & De Villiers, J.-P. 2004, ApJ, 616, 364
- Gammie, C. F. 2001, ApJ, 590, 174
- Gibson, R. R., et al. 2009, ApJ, 692, 758
- Gilman, R. C. 1972, ApJ, 178, 423
- Gofford, J., et al. 2011, MNRAS, 414, 3307
- Goodman, J. 2003, MNRAS, 339, 937
- Greene, J. E., Zakamska, N. L., & Smith, P. S. 2012, ApJ, 746, 86
- Greenhill, L. J., et al. 2003, ApJ, 590, 162
- Hayes, J. C., Norman, M. L., Fiedler, R. A., Bordner, J. O., Li, P. S., Clark, S. E., ud-Doula, A., & Mac Low, M. 2006, ApJS, 165, 188
- Higginbottom, N., Proga, D., Knigge, C., Long, K. S., Matthews, J. H., & Sim, S. A. 2014, ApJ, 789, 19
- Hopkins, P. F., Hayward, C. C., Narayanan, D., & Hernquist, L. 2012, MNRAS, 420, 320
- Igumenshchev, I. V. 2008, ApJ, 677, 317
- Jaffe, W., et al. 2004, Nature, 429, 47
- Kallman, T. R., & Bautista, M. A. 2001, ApJS, 133, 221
- Keating, S. K., Everett, J. E., Gallagher, S. C., & Deo, R. P. 2012, ApJ, 749, 32
- King, A. 2005, ApJ, 635, L121
- Konigl, A., & Kartje, J. F. 1994, ApJ, 434, 446
- Krolik, J. H. 2007, ApJ, 661, 52
- Krolik, J. H., & Begelman, M. C. 1988, ApJ, 329, 702
- Krolik, J. H., & Lepp, S. 1989, ApJ, 347, 179
- Krumholz, M. R., & Thompson, T. A. 2013, MNRAS, 434, 2329
- Kurosawa, R., & Proga, D. 2009, ApJ, 693, 1929
- Kurosawa, R., Proga, D., & Nagamine, K. 2009, ApJ, 707, 823
- Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999, Introduction to Stellar Winds
- Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
- Lin, D. N. C., & Pringle, J. E. 1987, MNRAS, 225, 607
- Lovelace, R. V. E., Romanova, M. M., & Biermann, P. L. 1998, A&A, 338, 856
- Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1
- Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561
- Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics, ed. Mihalas, D. & Mihalas, B. W.
- Minerbo, G. N. 1978, J. Quant. Spec. Radiat. Transf., 20, 541
- Mościbrodzka, M., & Proga, D. 2008, ApJ, 679, 626
- —. 2013, ApJ, 767, 156
- Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569
- Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69
- Nenkova, M., Sirocky, M. M., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 147
- Novak, G. S., Ostriker, J. P., & Ciotti, L. 2012, MNRAS, 427, 2734
- Ostriker, E. C., McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975
- Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914
- Paczynski, B. 1978, Acta Astron., 28, 91
- Page, M. J., Carrera, F. J., Stevens, J. A., Ebrero, J., & Blustin, A. J. 2011, MNRAS, 416, 2792
- Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721
- Phinney, E. S. 1989, in NATO ASIC Proc. 290: Theory of Accretion Disks, ed. F. Meyer, 457–+
- Pringle, J. E. 1992, MNRAS, 258, 811
- Proga, D. 2007, ApJ, 661, 693
- Proga, D., & Begelman, M. C. 2003, ApJ, 582, 69
- Proga, D., Ostriker, J. P., & Kurosawa, R. 2008, ApJ, 676, 101
- Raban, D., Jaffe, W., Röttgering, H., Meisenheimer, K., & Tristram, K. R. W. 2009, MNRAS, 394, 1325
- Rafikov, R. R. 2009, ApJ, 704, 281
- Reeves, J. N., O’Brien, P. T., & Ward, M. J. 2003, ApJ, 593, L65
- Roth, N., Kasen, D., Hopkins, P. F., & Quataert, E. 2012, ArXiv e-prints
- Rowan-Robinson, M. 1977, ApJ, 213, 635
- Sanders, D. B., Phinney, E. S., Neugebauer, G., Soifer, B. T., & Matthews, K. 1989, ApJ, 347, 29
- Schartmann, M., Wada, K., Prieto, M. A., Burkert, A., & Tristram, K. R. W. 2014, MNRAS, 445, 3878
- Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
- Shi, J., & Krolik, J. H. 2008, ApJ, 679, 1018
- Shlosman, I., Begelman, M. C., & Frank, J. 1990, Nature, 345, 679
- Shlosman, I., Frank, J., & Begelman, M. C. 1989, Nature, 338, 45
- Silk, J., & Rees, M. J. 1998, A&A, 331, L1
- Sim, S. A., Proga, D., Miller, L., Long, K. S., & Turner, T. J. 2010, MNRAS, 408, 1396
- Spitzer, L. 1978, Physical processes in the interstellar medium
- Stevens, I. R., & Kallman, T. R. 1990, ApJ, 365, 321
- Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753
- Tassoul, J. 1978, Theory of rotating stars, ed. Tassoul, J.-L.
- Toomre, A. 1964, ApJ, 139, 1217
- Tristram, K. R. W., et al. 2007, A&A, 474, 837
- Trump, J. R., et al. 2006, ApJS, 165, 1
- Wada, K. 2012, ApJ, 758, 66
- Wada, K., & Norman, C. A. 2001, ApJ, 547, 172