Magnetohydrodynamic Simulations of Active Galactic Nucleus Disks and Jets
MMagnetohydrodynamicSimulations of Active GalacticNucleus Disks and Jets
Shane W. Davis and Alexander Tchekhovskoy Department of Astronomy, University of Virginia, Charlottesville, Virginia,22904, USA; email: [email protected] Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA),Department of Physics & Astronomy, Northwestern University, Evanston,Illinois, 60657; email: [email protected] .Rev. Astron. Astrophys. 2020.58:407–441This article’s doi:10.1146/annurev-astro-081817-051905Copyright © Keywords black hole physics, radiation transport, general relativity
Abstract
There is a broad consensus that accretion onto supermassive black holesand consequent jet formation power the observed emission from ac-tive galactic nuclei (AGNs). However, there has been less agreementabout how jets form in accretion flows, their possible relationship toblack hole spin, and how they interact with the surrounding medium.There have also been theoretical concerns about instabilities in stan-dard accretion disk models and lingering discrepancies with observa-tional constraints. Despite seemingly successful applications to X-raybinaries, the standard accretion disk model faces a growing list of ob-servational constraints that challenge its application to AGNs. Theo-retical exploration of these questions has become increasingly relianton numerical simulations owing to the dynamic nature of these flowsand the complex interplay between hydrodynamics, magnetic fields,radiation transfer, and curved spacetime. We conclude the following: • The advent of general relativistic magnetohydrodynamics(MHD) simulations has greatly improved our understanding ofjet production and its dependence on black hole spin. • Simulation results show both disks and jets are sensitive to themagnetic flux threading the accretion flow as well as possiblemisaglingment between the angular momentum of the accretionflow and the black hole spin. • Radiation MHD simulations are providing new insights intothe stability of luminous accretion flows and highlighting thepotential importance of radiation viscosity, UV opacity fromatoms, and spiral density waves in AGNs. a r X i v : . [ a s t r o - ph . H E ] J a n ontents
1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4082. CURRENT STATUS OF OBSERVATION AND THEORY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4102.1. Accretion Disk Theory in AGNs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4102.2. Observational Status: AGNs Versus XRBs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4122.3. Status of the Standard Jet Model Interpretation of AGNs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4173. OVERVIEW OF ACCRETION DISK AND JET SIMULATION METHODS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4203.1. General Relativistic MHD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4203.2. Radiation Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4224. MHD SIMULATIONS OF AGN ACCRETION DISKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4244.1. Thermal Instability in Radiation Pressure Dominated Accretion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4244.2. Global Simulations of AGN Accretions Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4255. SIMULATIONS OF AGN JETS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4285.1. The Disk-jet Connection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4285.2. Origin of Large-Scale Vertical Magnetic Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4295.3. Jet Acceleration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4315.4. Jet Emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4315.5. Tilted Disk Precession, Jets, Alignment, and Tearing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4325.6. Jet Interaction with Ambient Medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4336. SUMMARY AND OUTLOOK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435
1. INTRODUCTION
The early observations of quasars and the realization that they must be distant objects(Schmidt 1963) already indicated a need for objects much smaller than galaxies to beproducing immense energies to power their optical emission and radio lobes. Accretion
Quasar: a luminousactive galacticnucleus thatoutshines its hostgalaxy of interstellar gas onto a supermassive black hole was identified as the only viable sourcesfor such energy (Salpeter 1964; Lynden-Bell 1969).The first detailed model of the vertical and radial structure of disks with angular mo-mentum transport was provided by Shakura & Sunyaev (1973) with additional insightsfrom Lynden-Bell & Pringle (1974) and general relativistic generalization by Novikov &Thorne (1973). The key assumptions in the model were that the vertical scale height of thedisk was small compared with the radius everywhere and that the stress leading to angularmoment transport could be approximated as a constant fraction ( α ) of the total pressure inthe disk midplane. Even 45 years later, this seminal work remains the primary basis of ourunderstanding of accretion onto black holes in XRBs and AGN and we refer to this as the X-ray Binary (XRB):
X-ray emittingsource powered byaccretion fromcompanion star ontoa black hole orneutron star
Active galacticnucleus (AGN): compact emittingregion in the centerof a galaxy thoughtto be powered byaccretion onto asupermassive blackhole standard model of accretion disks. Although these steady disk models suffer from instabil-ities discussed below, they were successful in broadly explaining the spectral properties ofXRBs (see e.g. Remillard & McClintock 2006; Done, Gierli´nski & Kubota 2007) and theoptical/UV continuum of luminous AGNs (Shields 1978; Kolykhalov & Sunyaev 1984).Although this standard model is thought to apply in the luminous accretion regime,there are observational and theoretical arguments that this model breaks down when theaccretion rate is very low or high when compared with the Eddington (or critical accretion)rate, ˙ M E = 4 πGMm p ησ T c , (1)
408 Davis & Tchekhovskoy here M is the mass of the black hole, η is the radiative efficiency, and other values arefundamental constants with standard meanings. This is the accretion rate that wouldproduce an Eddington luminosity of emission ( L E (cid:39) M/M (cid:12) ) if the radiative efficiency
Eddingtonluminosity:
Theluminosity for whichoutward radiationpressure on theelectrons equals theinward gravitationalpull on the ions remains constant (typically assumed to be around 10%). At low accretion rates ( ˙ M (cid:46) .
01 ˙ M E ), flows are thought to be hotter and less radiatively efficient (Yuan & Narayan 2014).In these models, number densities are low enough that Coulomb collisions are less efficientat coupling electrons to ions and radiative cooling is weak, leading to flows dominated byradial advection of thermal energy (Narayan & Yi 1994) and/or outflows (Blandford &Begelman 1999). At high accretion rates ( ˙ M (cid:38) ˙ M E ), the assumption of a thin disk is nolonger self-consistent as the vertical scale height of the inner disk becomes comparable withradius. The flow is thought to again become less radiatively efficient as the photon diffusiontime becomes longer than the radial inflow time, allowing a larger fraction of the dissipatedenergy to be advected into the black hole. Models including these effects are called slimdisks (Abramowicz et al. 1988).All these accretion disk models rely on an anomalous viscosity or stress ( α prescription)to explain angular momentum transport because realistic physical viscosities are too lowby many orders of magnitude. Much of the theoretical work in the field for several decadesfocused on discerning the origin of this anomalous stress. After the identification of the Magnetorotationalinstability (MRI):
Instability driven bymagnetic fieldsinteracting withshear that providesa source of angularmomentumtransport role of the (MRI; Balbus & Hawley 1991, 1998) in driving magnetohydrodynamic (MHD)turbulence in disks, it is now widely believed that MRI turbulent stresses provide the dom-inant angular momentum transport mechanism. To date, this has not led to a significant
Magnetohydrodynamics(MHD):
Thecombined magneticand fluid dynamicsof conductingplasmas revision of the standard model because the Maxwell and Reynolds stresses provided by theMRI turbulence are broadly consistent with the α prescription ansatz (Balbus & Papaloizou1999). However, in addition to disk internal stresses, disk outflows can also remove the an-gular momentum from the accretion disk and affect the disk structure. This is particularlytrue in the case of MHD-driven outflows, where nonlocal stresses are dominant (Blandford& Payne 1982). This suggests that the processes of accretion and outflows may be tightlycoupled.Observations indicate that accretion systems produce various types of outflows. Ofparticular interest are collimated outflows, more commonly referred to as jets. Assisted Blandford-Znajek(BZ):
Mechanism bywhich magneticfields are twisted byblack hole spin toform an outflow with large-scale magnetic fields, jets can derive their power from the rotation of either theblack hole via the Blandford & Znajek (1977) effect or the rotation of the accretion diskvia the Blandford & Payne (1982) mechanism. The jets can cover many decades in radius,spanning from the event horizon of the black hole to the outskirts of the galaxies andbeyond. Modeling jets is challenging owing to their nonlinear nature that gets exacerbatedby the poorly understood disk-jet connection near their base and turbulent interactionswith the ambient medium as they propagate away from the center.Although work on the revision and development of analytic models continues, the dy-namic and three-dimensional nature of the problem of coupled radiation transfer and MHDin general relativistic spacetimes has proven to be a challenge for further analytic progressand simpler modeling efforts. The critical role played by dynamic magnetic fields in thelaunching of jets and evolution of accretion disks has led to an increased focus and relianceon numerical simulation modeling. In this review, we summarize recent developments andthe current status of numerical simulations of disks and jets. We review observational re-sults that challenge the existing models and motivate the incorporation of new or improvedtreatments of physical processes in state-of-the-art simulations. ••
Mechanism bywhich magneticfields are twisted byblack hole spin toform an outflow with large-scale magnetic fields, jets can derive their power from the rotation of either theblack hole via the Blandford & Znajek (1977) effect or the rotation of the accretion diskvia the Blandford & Payne (1982) mechanism. The jets can cover many decades in radius,spanning from the event horizon of the black hole to the outskirts of the galaxies andbeyond. Modeling jets is challenging owing to their nonlinear nature that gets exacerbatedby the poorly understood disk-jet connection near their base and turbulent interactionswith the ambient medium as they propagate away from the center.Although work on the revision and development of analytic models continues, the dy-namic and three-dimensional nature of the problem of coupled radiation transfer and MHDin general relativistic spacetimes has proven to be a challenge for further analytic progressand simpler modeling efforts. The critical role played by dynamic magnetic fields in thelaunching of jets and evolution of accretion disks has led to an increased focus and relianceon numerical simulation modeling. In this review, we summarize recent developments andthe current status of numerical simulations of disks and jets. We review observational re-sults that challenge the existing models and motivate the incorporation of new or improvedtreatments of physical processes in state-of-the-art simulations. •• Disks and Jets in AGNs 409 . CURRENT STATUS OF OBSERVATION AND THEORY2.1. Accretion Disk Theory in AGNs
There are numerous presentations of the standard disk model equations, so we do notattempt to reproduce them here, but in our opinion the seminal paper by Shakura &Sunyaev (1973) remains the best place to start. For a detailed discussion of how this modelapplies to AGNs, we recommend the textbook by Krolik (1999). Here, we briefly reviewthe key results relevant to this work. Broadly speaking, the accretion flow divides into tworegions: an inner region, in which radiation pressure is the primary support against gravity,and an outer region, in which gas pressure dominates. The transition radius is a functionof disk parameters and is given by r tr (cid:39) r g (cid:18) . η L . L E (cid:19) / (cid:18) α . M M (cid:12) (cid:19) / R / R R / z R / T (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r tr , (2)where r g = GM/c is the gravitational radius, and α is the famous parameter that relatesthe accretion stress τ rφ to the total midplane pressure via τ rφ = αP tot . The radial- andblack-hole-spin-dependent quantities R R , R z , and R T are defined in chapter 7 of Krolik(1999) and represent general relativistic correction factors. Here, P tot is the sum of gaspressure P g and radiation pressure P r because the standard model does not consider theimpact of magnetic pressure. One can derive different scalings of quantities of interest inthese inner and outer regions. For brevity, we focus on the scalings in the radiation-pressuredominated inner regions and assume electron scattering dominates the opacity. We find T eff = 6 . × K (cid:18) M (cid:12) M . η L . L E (cid:19) / R / R (cid:18) rr g (cid:19) − / (3) T mid = 9 × K (cid:18) M (cid:12) M . α (cid:19) / (cid:18) R z R T R R (cid:19) / (cid:18) rr g (cid:19) − / (4) ρ mid = 3 . × − g / cm M (cid:12) M . α (cid:18) η . . L E L (cid:19) R z R T R R (cid:18) rr g (cid:19) / , (5)where r is radius. The effective temperature T eff is an estimate of the surface temperatureand ρ mid and T mid are the midplane density and temperature, respectively. The ratio of h/r , where h is the disk scale height for the radiation dominated flow is hr = 0 .
75 0 . η L . L E R R R z r g r . (6)(See the sidebar titled Thermal Instability in the Standard Accretion Disk Model.)Although the standard accretion model approximately succeeds at explaining SEDs of SED:
Spectralenergy distribution
AGN, there are a number of observational and theoretical issues that have arisen. Shortlyafter the model was formulated, it was realized that the model is subject to inflow (orviscous) and thermal instabilities (Lightman & Eardley 1974; Shakura & Sunyaev 1976)when radiation pressure exceeds gas pressure in the disk midplane. Inserting the abovescalings into standard expression for the radiation and gas pressure yields the scaling, P r P g ∝ α / M / η − (cid:18) LL E (cid:19) (cid:18) rr g (cid:19) − / . (7)
410 Davis & Tchekhovskoy hermal Instability in the Standard Accretion Disk Model
The thermal instability of radiation pressure-dominated accretion disks (Shakura & Sunyaev 1976) can beunderstood via a simple argument about how the heating and cooling rates in an accretion flow depend ontemperature (Piran 1978). In an accretion flow, the cooling rate per unit mass Q − from a disk annulus atradius r is simply given by the local radiative flux F r divided by the mass surface density (mass per unitarea) Σ so that (accounting for both sides of the disk) Q − = 2 F r Σ . The rate at which energy is dissipated in a shear flow is related to the product of the accretion stresses τ rφ and the shear rd Ω /dr . Integrating over height and dividing by Σ yields a heating rate per unit mass Q + = 1Σ (cid:90) τ rφ (cid:18) − r d Ω dr (cid:19) dz (cid:39) hτ rφ ΩΣ . Using the standard model equations in the radiation pressure-dominated regime with electron scatteringopacity yields h ∝ F r ∝ P r / Σ and τ rφ ∝ P r . Since P r = aT / Q + ∝ T Σ , Q − ∝ T Σ . In deriving the standard model, we assume the disk is in thermal equilibrium with Q − = Q + . Now supposethat the temperature is perturbed upward. Because departures from thermal equilibrium happen on thethermal timescale while changes in the surface density happen on the much longer viscous timescale weassume Σ is constant. An increase in the temperature causes a large increase in the heating rate due tothe T dependence and a somewhat more modest change in the cooling rate due to the T dependence. So,as the temperature is perturbed upward, the heating rate goes up faster than the cooling rate, leading toa thermal runaway. The same argument applies in reverse to downward perturbations in temperature, andwe conclude that radiation pressure-dominated disks are thermally unstable.Hence, for the same Eddington ratio, the radiation pressure in at 10 M (cid:12) black hole isexpected to 100 times more dominated by radiation pressure than a black hole with 10 M (cid:12) .This extreme dominance of radiation pressure in AGNs could suggest that the thermaland inflow instabilities could have a larger impact in AGNs. It is generally thought thatsuch instabilities will lead to limit cycle behavior in which the disk oscillates between a gasdominated lower branch and an advection dominated upper branch (Nayakshin, Rappaport& Melia 2000; Janiuk, Czerny & Siemiginowska 2002). The relevant timescales are thedynamical time, thermal time, and inflow (viscous) time given by t dyn = 1Ω = 1 day M M (cid:12) (cid:18) r r g (cid:19) / (8) t th = 1 α Ω = 9 day M M (cid:12) . α (cid:18) r r g (cid:19) / (9) t in = 1 α Ω r h = 260 yr M M (cid:12) . α (cid:18) h . r (cid:19) (cid:18) r r g (cid:19) / , (10) ••
The thermal instability of radiation pressure-dominated accretion disks (Shakura & Sunyaev 1976) can beunderstood via a simple argument about how the heating and cooling rates in an accretion flow depend ontemperature (Piran 1978). In an accretion flow, the cooling rate per unit mass Q − from a disk annulus atradius r is simply given by the local radiative flux F r divided by the mass surface density (mass per unitarea) Σ so that (accounting for both sides of the disk) Q − = 2 F r Σ . The rate at which energy is dissipated in a shear flow is related to the product of the accretion stresses τ rφ and the shear rd Ω /dr . Integrating over height and dividing by Σ yields a heating rate per unit mass Q + = 1Σ (cid:90) τ rφ (cid:18) − r d Ω dr (cid:19) dz (cid:39) hτ rφ ΩΣ . Using the standard model equations in the radiation pressure-dominated regime with electron scatteringopacity yields h ∝ F r ∝ P r / Σ and τ rφ ∝ P r . Since P r = aT / Q + ∝ T Σ , Q − ∝ T Σ . In deriving the standard model, we assume the disk is in thermal equilibrium with Q − = Q + . Now supposethat the temperature is perturbed upward. Because departures from thermal equilibrium happen on thethermal timescale while changes in the surface density happen on the much longer viscous timescale weassume Σ is constant. An increase in the temperature causes a large increase in the heating rate due tothe T dependence and a somewhat more modest change in the cooling rate due to the T dependence. So,as the temperature is perturbed upward, the heating rate goes up faster than the cooling rate, leading toa thermal runaway. The same argument applies in reverse to downward perturbations in temperature, andwe conclude that radiation pressure-dominated disks are thermally unstable.Hence, for the same Eddington ratio, the radiation pressure in at 10 M (cid:12) black hole isexpected to 100 times more dominated by radiation pressure than a black hole with 10 M (cid:12) .This extreme dominance of radiation pressure in AGNs could suggest that the thermaland inflow instabilities could have a larger impact in AGNs. It is generally thought thatsuch instabilities will lead to limit cycle behavior in which the disk oscillates between a gasdominated lower branch and an advection dominated upper branch (Nayakshin, Rappaport& Melia 2000; Janiuk, Czerny & Siemiginowska 2002). The relevant timescales are thedynamical time, thermal time, and inflow (viscous) time given by t dyn = 1Ω = 1 day M M (cid:12) (cid:18) r r g (cid:19) / (8) t th = 1 α Ω = 9 day M M (cid:12) . α (cid:18) r r g (cid:19) / (9) t in = 1 α Ω r h = 260 yr M M (cid:12) . α (cid:18) h . r (cid:19) (cid:18) r r g (cid:19) / , (10) •• Disks and Jets in AGNs 411 here Ω = (
GM/r ) / is the Keplerian rotation rate and h/r is the aspect ratio of thedisk scale height to radius. Time-dependent models yield evolution along stable branchesoccurring on the inflow timescale, whereas evolution between branches occurs on the thermaltimescale. If the disk is as thin as commonly thought, the inflow timescale is very longcompared to observing timescales for most systems.It was shown that the thermal and inflow instabilities could be suppressed in alternativestress prescriptions (Sakimoto & Coroniti 1981) where the stress scaled not with the totalpressure but gas pressure only. However, as we discuss below, evidence from radiation MHDsimulations of accretion flows suggests the stress does generally scale with the total pressureand, thus, provide some evidence that thermal instability is present (Jiang, Stone & Davis2013; Mishra et al. 2016). On the observational side, there are also a number of potential discrepancies. Koratkar &Blaes (1999) provide a somewhat dated but still useful review of several of the problems,many of which remain unresolved. Some of the issues that have received the most attentionare the apparent discrepancies between the standard model and sizes of emission regionsdetermined by microlensing (e.g Morgan et al. 2010; Blackburne et al. 2011; Morgan et al.2018) and continuum reverberation mapping (see e.g. Edelson et al. 2015; Jiang et al. 2017;Homayouni et al. 2019). A summary of the microlensing size scale constraints from Morganet al. (2018) is shown in Figure 1. Their best-fit relation suggests that emitting regions areabout a factor of 3-4 greater than the standard model predicts, which is consistent with thetime delays inferred from continuum reverberation mapping.Variability is also a challenge for disk models. Although the XRB GRS 1915+109famously shows variability that looks somewhat like the limit cycle predictions (Nayak-shin, Rappaport & Melia 2000), XRBs show no generic evidence of such variability (Done,Gierli´nski & Kubota 2007). In AGNs, the predicted viscous timescales are expected to betoo long (in some cases, thousands of years) to probe observationally. Instead, the observedvariability poses a different challenge. Surveys and monitoring campaigns have identified alarge number of “changing look” AGNs, that show large-amplitude variations in the con-tinuum or emission lines on timescale of months to years (LaMassa et al. 2015; MacLeodet al. 2016), which are difficult to explain with “viscous” evolution timescales (Lawrence2018). The question becomes how does one obtain large-amplitude variability of accretiondisks on such short timescales?Others have looked at the spectral variations of the UV continuum of AGNs with massestimates and found them to be poorly fit by spectral models based on the α -disks (see e.g.,Kolykhalov & Sunyaev 1984; Laor & Netzer 1989; Sincell & Krolik 1997; Bonning et al. 2007;Davis, Woo & Blaes 2007). The discrepancy arises in the sense that the model spectra risetoo rapidly into the far- to extreme-UV, whereas the observed spectra tend to be flat ordeclining in λF λ above ∼
412 Davis & Tchekhovskoy igure 1: Discrepancies between microlensing measured sizes and predictions from thestandard model. The solid line gives the best fit relation and the shaded regions shows the1 σ uncertainties. For comparison, the black dot-dashed curve provides the estimate for thinaccretion disk models with the same mass (assuming L/L E = 1 / η = 0 . ••
412 Davis & Tchekhovskoy igure 1: Discrepancies between microlensing measured sizes and predictions from thestandard model. The solid line gives the best fit relation and the shaded regions shows the1 σ uncertainties. For comparison, the black dot-dashed curve provides the estimate for thinaccretion disk models with the same mass (assuming L/L E = 1 / η = 0 . •• Disks and Jets in AGNs 413 igure 2:
Left panel:
Comparison of observed broadband quasar spectral energy distribu-tions from PG 0052+251 (Shang et al. 2005) with relativistic multitemperature blackbodymodels chosen to match the reverberation-mapped mass and with accretion rates chosen tomatch the optical emission. The blue X-ray curves show the level and slope change of theX-ray emission. Note that there is a turnover in the data in the ultraviolet band around3 × Hz (1000˚A), whereas the model (red dashed curve) continues rising into the extremeultraviolet. This is characteristic of most of the quasars in the PG sample (Davis & Laor2011).
Right panel:
The best-fit model to the black hole X-ray binary source LMC X-3(Davis, Done & Blaes 2006). The
BeppoSAX data are shown in black, whereas the best-fitmodel is shown as a green curve. The model has two components, with the best-fit absorbedrelativistic accretion disk model shown in red. The unabsorbed relativistic accretion diskmodel is show in orange. The model provides a good fit to the data with only two freeparameters in the disk model.be emission from a geometrically thin, optically thick accretion disk, whereas the low/hardstates is dominated by harder X-ray continuum component that is generally thought to beproduced by inverse Compton scattering of photons in a hot corona above or interior to thedisk. Although XRBs can have significant variability on shorter timescales, this is mostlyassociated with the coronal component. The high/soft state tends to be characterized byrelatively low variability, particularly in disk-like component of emission.If one focuses solely on the seemingly disk-dominated high/soft state, then the predic-tions of the standard model seem to match the data well (Davis et al. 2005; Dunn et al.2011). A good example of this, LMC X-3, is shown in the right panel of Figure 2 withthe best-fit model being a fully relativistic spectrum derived from self-consistent radiationtransfer calculations through a vertically varying disk structure (Davis et al. 2005). There-fore, a key question is why is there this seeming discrepancy in how well standard diskmodels reproduce the properties of AGNs and black hole XRBs? One possibility is thatwe have overstated the agreement between the models and XRB observations. Althoughthere are compelling models for the origin of the coronae and the nature of the transitionsbetween the different spectral states (e.g., Esin, McClintock & Narayan 1997; Ohsuga et al.2009), there is no consensus and this remains an active area of research. Whether or notone accepts the premise that standard models work for XRBs and fail for AGNs, there are
414 Davis & Tchekhovskoy everal reasons to expect that accretion flows in AGNs and XRBs might differ. We focuson three possibilities.The first possibility is that the standard model predicts that disks around supermassiveblack holes should be more radiation pressure-dominated than those around stellar massblack holes. For accretion rates typically inferred in high/soft state XRBs, these flowsare expected to be radiation pressure-dominated, but possibly only by factors of ∼ ≥ T ∼ × K knownas the Fe opacity bump is thought to produce extreme variability and help drive outflowsin the envelopes of massive stars (Jiang et al. 2018). Looking at equations (3) and (4),we see that inner regions of near-Eddington accretion rate flows tend to be hotter thanthis in the inner-most radii of the disk, so electron scattering will dominate. However, asone moves out, there is always a broad range of radii that are in the relevant temperatureand density regime. For black hole masses of (cid:38) M (cid:12) this happens in the inner severalhundred gravitational radii, where the bulk of optical to UV radiation is produced. It isconceivable that dynamics driving extreme convection and outflows in massive stars arealso generically occurring in massive black holes.A third set of considerations is the impact of the different environments and feeding of • Disks and Jets in AGNs 415 T ( K ) − κ R ( c m / g ) Figure 3: Comparison of a sum of electron scattering opacity with Rosseland mean opacitiesfrom the OPAL project (black curves) and a Kramer’s-like estimate free-free opacity (redcurves). Line shapes correspond to densities of 10 − (solid), 10 − (dotted), and 10 − (dash-dotted) g / cm . The modest bump at ∼ × K is the Fe opacity bump.the two systems. AGNs are thought to be fed predominantly from gas in the large-scaleinterstellar medium (ISM) of their host galaxies, whereas XRBs are fed by Roche lobeoverflow or wind from a companion star. Outer regions of AGN disks are expected to beself-gravitating, and either system could be subject to H ionization instabilities (Lasota2001). AGN accretion flows are also expected to be coincident with nuclear star clusters oftheir host galaxies and may have non-trivial interaction with stars or stellar remnants (see,e.g., Miralda-Escud´e & Kollmeier 2005; McKernan et al. 2012).Can these differences in the outer regions of the disks lead to differences in the flow nearthe inner radii, where most of the emission is thought to originate? Because most of thekey timescales (inflow, thermal, dynamical) decrease with radius, one might expect the diskto lose information about its large-scale feeding and slowly adjust to a similar steady statemodel if other considerations (e.g. opacities, radiation pressure) do not drive differences inthe innermost regions. However, this argument may not hold if the large-scale geometry ofthe magnetic fields is influenced by conditions in the outer flow. Evidence from numericalsimulations suggests that the presence of a poloidal field may impact both the accretion
Poloidal: thepoloidal componentof magnetic field isthe sum of R and z components disk and jet formation, but the presence of poloidal field may be partially dependent onwhether such disks can efficiently accrete a large-scale poloidal field to their inner regions(Lubow, Papaloizou & Pringle 1994; Beckwith, Hawley & Krolik 2009; Zhu & Stone 2018) orgenerate it in situ (Beckwith, Hawley & Krolik 2008; McKinney, Tchekhovskoy & Blandford2012; Liska, Tchekhovskoy & Quataert 2020).We emphasize that radiation plays a key role, particularly in the first two of these con-siderations and may in the third as well. If one wants to understand the differences between
416 Davis & Tchekhovskoy
GNs and XRBs, then models and simulations in which radiation is treated explicitly arerequired. The physics involved also requires dynamical and multidimensional treatments.Such considerations strongly motivate general relativistic radiation MHD numerical simu-lation to study the evolution of such flows with realistic opacities, initial conditions, andboundary conditions appropriate for AGNs. Of these questions, the most difficult problemfor simulations to address is the effects of differences in large-scale feeding. The discrep-ancies in timescales between the largest and smallest radii are a significant impediment tonumerical simulations, which must resolve the shortest flow times in the inner regions butmust be run for many inflow times to achieve steady states at large radii. Nevertheless,simulations can examine what impact, if any, the magnetic field strength and geometryhave on the accretion flow because the presence and strength of the poloidal magnetic fieldcan be controlled through the initial conditions. Hence, we think that the key questionsdriving numerical simulation of luminous accretion flows in AGNs are the following: • What role does radiation and radiation pressure play in AGNs? Are radiationpressure-driven instabilities present, and how do they manifest in the flow structureand dynamics? • Are a significant fraction of observed AGNs accreting near enough to (or above) theEddington limit so that effects of advection and outflows become important? Howdo these effects modify the light curves and SEDs of accretion flows in this limit? • What role do dust and atomic opacities play in the dynamics and structure of AGNs?Are there significant outflows of gravitationally bound or unbound material launchedfrom the inner regions of AGNs even when the accretion rates are below the Eddingtonlimit for electron scattering opacity? • What role do magnetic fields play in AGNs? Does magnetic pressure support changethe structure of the disks? Is there any reason to believe that magnetic field geometriesdiffer between AGNs and XRBs?Many of these effects are equally important for understanding accretion in XRBs andpossibly ultraluminous X-ray sources, which may be the result of accretion above the Ed-dington limit. However, here we focus on the issues specific to AGNs.
It is generally agreed that for a black hole to launch jets it needs to be accreting. However,the ingredients necessary for jet formation are poorly understood. Is it the black hole spinor the rotation of the accretion disk that powers the jets? Complicating the situation, theability of black hole systems to produce jets depends on the state of the accretion flow:The same black hole can exist in jetted and jet-phobic states depending on how it is fed.Perhaps the best evidence for this comes from microquasars: XRBs show radio emissionthat typically varies during spectral state transitions (Fender, Belloni & Gallo 2004). Beforethe transition, the system is in a low/hard accretion state characterized by a nonthermalspectrum and accompanied by weak continuous jets. As the transition starts, the luminosityof the system goes up, and powerful transient jets emerge. As the luminosity declines, thejets disappear, revealing a near-thermal accretion flow. The power of transient jets appearsto correlate with the black hole spin (Narayan & McClintock 2012; Steiner, McClintock &Narayan 2013; but see Russell, Gallo & Fender 2013), suggesting that black hole rotationplays a role in powering these outflows. However, how this correlation emerges is unclear. ••
It is generally agreed that for a black hole to launch jets it needs to be accreting. However,the ingredients necessary for jet formation are poorly understood. Is it the black hole spinor the rotation of the accretion disk that powers the jets? Complicating the situation, theability of black hole systems to produce jets depends on the state of the accretion flow:The same black hole can exist in jetted and jet-phobic states depending on how it is fed.Perhaps the best evidence for this comes from microquasars: XRBs show radio emissionthat typically varies during spectral state transitions (Fender, Belloni & Gallo 2004). Beforethe transition, the system is in a low/hard accretion state characterized by a nonthermalspectrum and accompanied by weak continuous jets. As the transition starts, the luminosityof the system goes up, and powerful transient jets emerge. As the luminosity declines, thejets disappear, revealing a near-thermal accretion flow. The power of transient jets appearsto correlate with the black hole spin (Narayan & McClintock 2012; Steiner, McClintock &Narayan 2013; but see Russell, Gallo & Fender 2013), suggesting that black hole rotationplays a role in powering these outflows. However, how this correlation emerges is unclear. •• Disks and Jets in AGNs 417 f i e l d li n e t=0 t=t t=t t=t v °uid ¼ v ¯eld (a) (b) (c) (d) p m = B ' = ¼ B ' À B z © Figure 4: Illustration of jet formation by magnetic fields. (a) Consider a purely poloidal(i.e., toroidal field vanishes, B ϕ = 0) field line attached on one end to a stationary “ceiling”(which represents the ambient medium and is shown by a hashed horizontal line) and onthe other end to a perfectly conducting sphere (which represents the central black holeor neutron star and is shown by a gray filled circle) rotating at an angular frequency Ω.(b) After N rotations, at time t = t , the initially purely poloidal field line develops N toroidal loops. This magnetic spring pushes against the ceiling with an effective pressure p m ∼ B ϕ / π due to the toroidal field, B ϕ . As time goes on, more toroidal loops form, andthe toroidal field becomes stronger. (c) At some later time, t = t , the pressure becomes solarge that the magnetic spring, which was twisted by the rotation of the sphere, pushes awaythe ceiling and accelerates the plasma attached to it along the rotation axis, forming a jet.Asymptotically far from the center, the toroidal field is the dominant field component anddetermines the dynamics of the jet. (d) It is convenient to think of the jet as a collectionof toroidal field loops that slide down the poloidal field lines and accelerate along the jetunder the action of their own pressure gradient and tension (hoop stress). The rotation ofthe sphere continuously twists the poloidal field into new toroidal loops at a rate that, insteady state, balances the rate at which the loops move downstream. The power of the jet isdetermined by two parameters (equation 11): the rotational frequency of the central object,Ω, and the radial magnetic flux threading the object, Φ. Figure adapted from Tchekhovskoy(2015).To understand the disk-jet connection, let us first review the basics of jet formationwithout attempting to connect it to the disk. For simplicity, let us consider a perfectlyconducting spinning sphere, shown in Fig. 4(a), which is meant to represent the centralcompact object (a black hole, neutron star, white dwarf, or even normal star), and a per-fectly conducting “ceiling”, which is meant to represent the ambient medium. Suppose amagnetic field line connects the sphere to the ceiling. As the sphere rotates, it coils up thefield line into a magnetic spring (Fig. 4b), which pushes the ceiling away and acceleratesunder the action of its own pressure (Fig. 4c). One can view this acceleration in an alterna-tive way (Fig. 4d): the rotation of the central sphere continuously reprocesses the initiallyvertical field line into toroidal field loops that emanate from the sphere and slide up thejet. As they do so, they expand, their pressure drops, and the pressure gradient acceleratesthem away. The situation for black holes, qualitatively, is rather similar. Even though ablack hole does not have a physical surface, black hole rotation leads to the dragging ofthe inertial frames near the black hole. This causes the magnetic field lines to rotate in asurprisingly similar fashion to the perfectly conducting sphere.
418 Davis & Tchekhovskoy landford & Znajek (1977) showed, in the limit of slow rotation, that a spinning blackhole immersed into a large-scale vertical magnetic field would produce jets of power P BZ = k ( a/ r g ) Φ c, (11)where Φ BH is the magnetic flux threading the black hole event horizon, − ≤ a ≤ k is a dimensionless proportionality factor.Ignoring constant prefactors, we can obtain this expression from dimensional analysis bywriting that the jet power is the product of magnetic energy density, ∝ B , the cross-sectionof the base of the jet, ∝ r g , and the speed v ∼ c with which the energy flows through thejets, giving P ∝ a r g B c . Here, we introduce the a prefactor to account for the variationof jet power with black hole spin (it has to be an even power of spin because spin signchange, by symmetry, leaves the power unchanged). Switching from the field strength, B ,to the magnetic flux, Φ BH ∼ Br g , gives P ∝ ( a/r g ) Φ c , which has the same scaling asin eq. (11).With the advent of numerical simulations, it became clear that though equation (11)works well for small values of a , it underestimates the jet power for rapidly spinning blackholes (Komissarov 2001). Comparison against numerical solutions shows that replacing r g with the event horizon radius, r H = r g [1 + (1 − a ) / ], gives a more accurate expression, P jet = k ( a/r H ) Φ c × f ( a ) , (12)that serves well for most practical purposes at a (cid:46) .
95 with f ( a ) = 1 (Tchekhovskoy,Narayan & McKinney 2010). A higher-order correction, f ( a ) = 1 + 0 . ar g /r H ) − . ar g /r H ) , allows eq. (12) to maintain accuracy all the way up to a = 1 (Tchekhovskoy,Narayan & McKinney 2010; Pan & Yu 2015).Thus, there appears to be a clear relationship between the black hole spin and jet power.What, then, causes the jet power to change for any given black hole? Since in AGNs theblack hole spin does not change in our lifetime, changes in jet power can only come from thevariations in the black hole magnetic flux, Φ BH . To understand the disk-jet connection, wetherefore need to understand how the accretion physics determines the value of Φ BH . If wewere to remove the accretion disk, then the black hole would lose its magnetic flux: By theno hair theorem, the black hole can only have three hairs – mass, spin, and charge (Misner,Thorne & Wheeler 1973). Thus, the accretion disk holds the magnetic flux on the blackhole and prevents it from slipping away: Pressure of the magnetic flux on the black holemust be in some way balanced by the pressure of the accretion flow. Here, we can againturn to the dimensional analysis and write that the magnetic pressure on the black hole isproportional to Φ . This balances the ram pressure in the inflow, which should be roughlyproportional to the mass accretion rate, ˙ M . To characterize the relative strength of the two,we introduce a dimensionless black hole magnetic flux, φ BH = Φ BH / ( ˙ Mr g c ) / . Its value isset by the interaction with the disk. The nonlinearity and intrinsic time-variability of thedisk-jet interaction makes numerical simulations an attractive approach for constrainingthe allowed range of values of φ BH . As we discuss in Section 5, the allowed range of φ BH spans the range from zero to a maximum value around 50 for which the black hole magneticflux becomes as strong as the gravity that keeps the disk on an orbit around the black hole(Tchekhovskoy, Narayan & McKinney 2011). In this magnetically arrested disk (MAD) Magneticallyarrested disk(MAD): accretionstate wheremagnetic flux insidethe disk becomesstrong enough todisrupt the flow state, the magnetic flux is as strong as possible and weakly depends on the black hole spin, φ MAD ≈ − . ar g /r H ) h / . , where h = r × . h . is the half-thickness of the disk. • Disks and Jets in AGNs 419 his corresponds to the following maximum jet power: P MAD = η MAD ˙ Mc ≈ . h . a ˙ Mc , (13)where we used the fact that the jet efficiency is η MAD = kφ a ( r g /r H ) f ( a ) ≈ . a h . (Tchekhovskoy 2015). Thus, the jet power P jet takes on a value between 0 and P MAD depending on the strength of the black hole magnetic flux φ BH relative to its maximumpossible value, φ MAD :0 ≤ P jet = (cid:18) φ BH φ MAD (cid:19) P MAD ≤ P MAD = 1 . h . a ˙ Mc . (14)Now that we understand the jet power, there are still many questions left unanswered.These questions are driving detailed studies in various branches of astrophysics include thefollowing: • Where does the magnetic flux powering jets come from? Does it need to be draggedinward from the ISM or can it be generated inside the disk in situ? • How are transient jets launched during spectral state transitions in XRBs? Does ananalog of the state transitions and transient jets exist in AGNs? • How do jets accelerate to relativistic velocities? Is radiation pressure important inlaunching and accelerating the jets? • Do jets heat up the interstellar gas and affect galaxy evolution? Can jet feedbacklead to the M − σ relationship? • What makes jets shine? What can we learn from the observations of the black holeshadow with the Event Horizon Telescope? What can the jets tell us about thestrong-field gravity and general relativistic frame dragging that birthed them?We discuss several of these questions in Sections 3 and 5.
3. OVERVIEW OF ACCRETION DISK AND JET SIMULATION METHODS3.1. General Relativistic MHD
The numerical simulation of accretion flows onto black holes has a long history, with Wilson(1972) already considering the flow of fluid with nonzero angular momentum in the Kerrspacetime. Early work primarily focused on hydrodynamics models (e.g., Hawley, Smarr& Wilson 1984), but the realization that the MRI could provide the necessary angularmomentum transport led to the development and application of MHD methods (Balbus &Hawley 1998). With this realization it became standard for numerical simulations to solvesome version of the MHD equations, which includes conservation of mass, ∂ρ∂t + ∇ · ( ρ v ) = 0 , (15)conservation of momentum, ∂ ( ρ v ) ∂t + ∇ · ( ρ vv − BB + P ∗ ) = − ρ ∇ Φ , (16)conservation of energy, ∂E∂t + ∇ · [( E + P ∗ ) v − B ( B · v )] = − ρ v · ∇ Φ , (17)
420 Davis & Tchekhovskoy umerical Methods for Solving the MHD Equations
There are a significant number of different methods available for simulating MHD flows and an even largernumber of codes that implement them. Most methods for solving fluid dynamics equations fall into twocategories: Eulerian or Lagrangian. Eulerian methods tend to solve finite difference representations of theequations on a fixed mesh. Finite volume (Godunov) methods (LeVeque 2002) are particularly popularfor their shock-capturing capabilities because shocks commonly arise in astrophysical flows. Lagrangianmethods tend to utilize particle-based representations with fluid properties advected along with the particles,such as smoothed particle hydrodynamics (Liu & Liu 2003). Other examples include moving mesh codes thatcarry a time variable mesh along with the particles or fluid flow (e.g., Springel 2010; Duffell & MacFadyen2011).Both Eulerian and Lagrangian codes have their merits, but the majority of the simulations discussedin this review utilize mesh-based Eulerian schemes. This is partly due to their shock-capturing capabilitiesand their simplicity relative to moving mesh methods. Another major consideration is the need to preservethe divergence-free nature of the magnetic field: ∇ · B = 0. Although it should always hold in nature,this condition is not preserved by all integration methods. Several codes address this issue by allowing thedivergence to develop but then limiting it via a divergence cleaning scheme. Other integration schemes arespecifically designed to preserve the divergence-free condition to machine precision. A popular example isthe constrained transport scheme (Evans & Hawley 1988), which is most easily implemented on structuredEulerian mesh (c.f. Mocz, Vogelsberger & Hernquist 2014).and the induction equation, ∂ B ∂t − ∇ × ( v × B ) = 0 . (18)Here, ρ , B , v are density, magnetic field and flow velocity, P ∗ ≡ ( P g + B / I (with I theunit tensor), P g is the gas pressure, and the magnetic permeability µ = 1. The total gasenergy density is E = E g + 12 ρv + B , (19)and Φ is the gravitational potential. The precise form of these equations varies depending onapplication. Sometimes the energy equation is omitted and an isothermal equation of stateis adopted. Additional source terms may be added to the right hand sides of the equationsof momentum and energy, such as the effects of radiative cooling, heating, or radiationforces. This may require the solution of additional equations, such as the radiation transferequation or its moments.Early MHD simulations of accretion focused on shearing box simulations (Hawley, Gam-mie & Balbus 1995; Brandenburg et al. 1995; Stone et al. 1996), which adopt shearingboundary conditions and add the Coriolis force to mimic the effects of fluid rotation inan accretion disk but only simulate a small patch of the disk. Although these simulationsare still widely used to study the structure, dynamics and thermodynamics of the disk athigh resolution, concerns about convergence with resolution (Pessah, Chan & Psaltis 2007;Fromang & Papaloizou 2007; Ryan et al. 2017) and dependence on box size (e.g. Simon, ••
There are a significant number of different methods available for simulating MHD flows and an even largernumber of codes that implement them. Most methods for solving fluid dynamics equations fall into twocategories: Eulerian or Lagrangian. Eulerian methods tend to solve finite difference representations of theequations on a fixed mesh. Finite volume (Godunov) methods (LeVeque 2002) are particularly popularfor their shock-capturing capabilities because shocks commonly arise in astrophysical flows. Lagrangianmethods tend to utilize particle-based representations with fluid properties advected along with the particles,such as smoothed particle hydrodynamics (Liu & Liu 2003). Other examples include moving mesh codes thatcarry a time variable mesh along with the particles or fluid flow (e.g., Springel 2010; Duffell & MacFadyen2011).Both Eulerian and Lagrangian codes have their merits, but the majority of the simulations discussedin this review utilize mesh-based Eulerian schemes. This is partly due to their shock-capturing capabilitiesand their simplicity relative to moving mesh methods. Another major consideration is the need to preservethe divergence-free nature of the magnetic field: ∇ · B = 0. Although it should always hold in nature,this condition is not preserved by all integration methods. Several codes address this issue by allowing thedivergence to develop but then limiting it via a divergence cleaning scheme. Other integration schemes arespecifically designed to preserve the divergence-free condition to machine precision. A popular example isthe constrained transport scheme (Evans & Hawley 1988), which is most easily implemented on structuredEulerian mesh (c.f. Mocz, Vogelsberger & Hernquist 2014).and the induction equation, ∂ B ∂t − ∇ × ( v × B ) = 0 . (18)Here, ρ , B , v are density, magnetic field and flow velocity, P ∗ ≡ ( P g + B / I (with I theunit tensor), P g is the gas pressure, and the magnetic permeability µ = 1. The total gasenergy density is E = E g + 12 ρv + B , (19)and Φ is the gravitational potential. The precise form of these equations varies depending onapplication. Sometimes the energy equation is omitted and an isothermal equation of stateis adopted. Additional source terms may be added to the right hand sides of the equationsof momentum and energy, such as the effects of radiative cooling, heating, or radiationforces. This may require the solution of additional equations, such as the radiation transferequation or its moments.Early MHD simulations of accretion focused on shearing box simulations (Hawley, Gam-mie & Balbus 1995; Brandenburg et al. 1995; Stone et al. 1996), which adopt shearingboundary conditions and add the Coriolis force to mimic the effects of fluid rotation inan accretion disk but only simulate a small patch of the disk. Although these simulationsare still widely used to study the structure, dynamics and thermodynamics of the disk athigh resolution, concerns about convergence with resolution (Pessah, Chan & Psaltis 2007;Fromang & Papaloizou 2007; Ryan et al. 2017) and dependence on box size (e.g. Simon, •• Disks and Jets in AGNs 421 eckwith & Armitage 2012; Shi, Stone & Huang 2016) remain. (See the sidebar titledNumerical Methods for Solving the MHD Equations.)Global simulations of MHD black hole accretion flows were first performed with pseudo-Newtonian potentials (Hawley & Balbus 2002) but these were supplanted by general rel-ativistic MHD (GRMHD) simulations of accretion flows onto spinning and nonspinningblack holes (Gammie, McKinney & T´oth 2003; De Villiers, Hawley & Krolik 2003). InGRMHD, the covariant generalizations of the MHD equations are conservation of mass,stress-energy, and the relativistic Maxwell’s equations, which are evolved on a choice ofcoordinates (usually Boyer-Lindquist or Kerr-Schild) dictated by the Kerr spacetime.
GRMHD:
Generalrelativistic magneto-hydrodynamics
The first GRMHD simulations lacked realistic cooling and, therefore, best approximatelow-luminosity accretion flows. Nevertheless, they have led to a much better understandingof flow in the innermost regions near the black hole, including accretion in the plungingregion, driving of outflows, and launching of jets.In particular, it became clear that in the presence of large-scale vertical magnetic fluxjets are a typical outcome of radiatively inefficient, geometrically thick black hole accretiondisks with h/r ∼ . − h/r (cid:46) .
05 weresimulated via the inclusion of cooling terms that kept the disk thickness at a desired value(Shafee et al. 2008; Penna et al. 2010; Noble, Krolik & Hawley 2010). These simulations didnot produce jets, in agreement with the observations that show no detectable radio emission(and, hence jet activity) from thin disks in XRBs and the vast majority (90%) of quasars.The remaining 10% of quasars show powerful radio jets whose origins are still a mystery.Recent GRMHD simulations show that even very thin disks, h/r ∼ . cells) are requiredto achieve convergence in the value of the effective α -viscosity parameter in global GRMHDsimulations of weakly magnetized accretion disks (Porth et al. 2019). Global simulationsof strongly magnetized accretion disks are less sensitive to resolution (White, Stone &Quataert 2019). Adding and improving radiation transfer treatments in numerical simulations of accretionflows is an important focus of recent work. The relevant equation to solve is the radiativetransfer equation, 1 c ∂I ν ∂t + n · ∇ I ν = η ν − α ν I ν , (20)where I ν is the specific intensity, η ν is the emissivity, α ν = κ ν ρ is the extinction coefficient,and κ ν is the opacity. The form on the right-hand side is general if extinction includesscattering opacity contributions and the emissivity accounts for scattered radiation. Because I ν is a function of position, angle, and frequency, it is usually computationally expensivecompared to the standard MHD equations, which only depend on position. For this reason,some treatments focus on its angle and frequency-integrated moments corresponding toconservation of radiation energy, ∂E r ∂t + ∇ · F r = cρ (cid:0) κ P aT − κ E Er (cid:1) , (21)
422 Davis & Tchekhovskoy nd conservation of radiation momentum,1 c ∂ F r ∂t + ∇ · P r = − κ F ρc F r . (22)Here, E r , F r , and P r are the radiation energy density, flux, and pressure tensor, respectively.The opacities κ P , κ E , and κ F correspond to the Planck, energy, and flux mean opacities.In most applications, the Rosseland mean is used for κ F and the Planck mean is used for κ E .For the sake of brevity, we will summarize some of the most salient points and referthe reader to Mihalas & Mihalas (1984) for a comprehensive discussion of theses equations.Although these equations look simple, one should note that, as written, all variables are inthe Eulerian frame, but the simple isotropic forms for the emissivities and opacities onlyhold in the comoving frame. Hence, one usually must Lorentz transform the source terms,expand the right-hand sides in powers of v/c , or use a comoving frame approach whichintroduces additional acceleration terms. Compton scattering also generally introducesadditional terms. The source terms on the right hand side of equations (21) and (22)represent the transfer of energy and momentum from the radiation field to fluid. Hence,the negative of these terms are added to the right-hand side of equations (16) and (17),corresponding to net heating/cooling by the radiation field and the radiation force.Early efforts to include radiation predominantly focused on simulations utilizing flux-limited diffusion (FLD; Levermore & Pomraning 1981). In FLD, only equation (21) is Flux-limited diffusion(FLD): diffusion-likeapproximation toradiation transfer retained, and F r is computed from the gradient of E r using a diffusion approximation. Alimiter is utilized to keep | F r | ≤ cE r in the optically thin regime. This method has been usedin the context of shearing box (Turner 2004) and global accretion flow (Ohsuga et al. 2005)simulations. More recently, several groups have developed general relativistic radiationtransfer modules based on two moment methods with M1 closure (S¸adowski et al. 2013;McKinney et al. 2014; Mishra et al. 2016). In M1 closure schemes, equations (21) and (22) M1: approximateclosure scheme usedrelate radiationpressure tensor tothe radiation fluxand energy density or their relativistic generalization (corresponding to conservation of radiation stress-energy)are solved. Because P r appears in equation (22), it must be specified. This is generally doneby introducing a new tensor called the Eddington tensor f ≡ P r /E r which characterizes theangular distribution of the radiation field. In the M1 scheme, f is computed as a functionof E r and F r .Other approaches include direct solution of the angle-dependent transfer equation(Jiang, Stone & Davis 2014b) or the variable Eddington tensor (VET) method (Jiang,Stone & Davis 2012). In the direct solution approach, a frequency-averaged version of variable Eddingtontensor (VET): closure scheme usedcompute radiationpressure tensor froman approximatesolution of thetransfer equation. equation (20) is solved explicitly for a fixed number of angles in each zone. In the VETmethod, equations (21) and (22) are integrated directly, but a time-independent solutionof the transfer equations is used to compute the Eddington tensor. Becasue both the FLDand M1 methods make assumptions about the angular distribution of the radiation field byprescribing relations among E r , F r or P r , methods involving direct solution of the transferequation are generally thought to be more accurate. The trade-off is that they are also morecomputationally expensive and are algorithmically more complex to implement. For thisreason, incorporating general relativity in the radiation transfer equation is significantlymore challenging, so simulations performed with these methods have all been nonrelativis-tic MHD-based, whereas many of the simulations performed with M1 methods have beengeneral relativistic MHD-based. The formalism for the general relativistic transfer equationis well-developed (Cardall, Endeve & Mezzacappa 2013) and efforts to implement it in blackhole accretion simulations are ongoing. ••
422 Davis & Tchekhovskoy nd conservation of radiation momentum,1 c ∂ F r ∂t + ∇ · P r = − κ F ρc F r . (22)Here, E r , F r , and P r are the radiation energy density, flux, and pressure tensor, respectively.The opacities κ P , κ E , and κ F correspond to the Planck, energy, and flux mean opacities.In most applications, the Rosseland mean is used for κ F and the Planck mean is used for κ E .For the sake of brevity, we will summarize some of the most salient points and referthe reader to Mihalas & Mihalas (1984) for a comprehensive discussion of theses equations.Although these equations look simple, one should note that, as written, all variables are inthe Eulerian frame, but the simple isotropic forms for the emissivities and opacities onlyhold in the comoving frame. Hence, one usually must Lorentz transform the source terms,expand the right-hand sides in powers of v/c , or use a comoving frame approach whichintroduces additional acceleration terms. Compton scattering also generally introducesadditional terms. The source terms on the right hand side of equations (21) and (22)represent the transfer of energy and momentum from the radiation field to fluid. Hence,the negative of these terms are added to the right-hand side of equations (16) and (17),corresponding to net heating/cooling by the radiation field and the radiation force.Early efforts to include radiation predominantly focused on simulations utilizing flux-limited diffusion (FLD; Levermore & Pomraning 1981). In FLD, only equation (21) is Flux-limited diffusion(FLD): diffusion-likeapproximation toradiation transfer retained, and F r is computed from the gradient of E r using a diffusion approximation. Alimiter is utilized to keep | F r | ≤ cE r in the optically thin regime. This method has been usedin the context of shearing box (Turner 2004) and global accretion flow (Ohsuga et al. 2005)simulations. More recently, several groups have developed general relativistic radiationtransfer modules based on two moment methods with M1 closure (S¸adowski et al. 2013;McKinney et al. 2014; Mishra et al. 2016). In M1 closure schemes, equations (21) and (22) M1: approximateclosure scheme usedrelate radiationpressure tensor tothe radiation fluxand energy density or their relativistic generalization (corresponding to conservation of radiation stress-energy)are solved. Because P r appears in equation (22), it must be specified. This is generally doneby introducing a new tensor called the Eddington tensor f ≡ P r /E r which characterizes theangular distribution of the radiation field. In the M1 scheme, f is computed as a functionof E r and F r .Other approaches include direct solution of the angle-dependent transfer equation(Jiang, Stone & Davis 2014b) or the variable Eddington tensor (VET) method (Jiang,Stone & Davis 2012). In the direct solution approach, a frequency-averaged version of variable Eddingtontensor (VET): closure scheme usedcompute radiationpressure tensor froman approximatesolution of thetransfer equation. equation (20) is solved explicitly for a fixed number of angles in each zone. In the VETmethod, equations (21) and (22) are integrated directly, but a time-independent solutionof the transfer equations is used to compute the Eddington tensor. Becasue both the FLDand M1 methods make assumptions about the angular distribution of the radiation field byprescribing relations among E r , F r or P r , methods involving direct solution of the transferequation are generally thought to be more accurate. The trade-off is that they are also morecomputationally expensive and are algorithmically more complex to implement. For thisreason, incorporating general relativity in the radiation transfer equation is significantlymore challenging, so simulations performed with these methods have all been nonrelativis-tic MHD-based, whereas many of the simulations performed with M1 methods have beengeneral relativistic MHD-based. The formalism for the general relativistic transfer equationis well-developed (Cardall, Endeve & Mezzacappa 2013) and efforts to implement it in blackhole accretion simulations are ongoing. •• Disks and Jets in AGNs 423 . MHD SIMULATIONS OF AGN ACCRETION DISKS
To date, there have been relatively few simulations directly aimed at studying supermassiveblack hole accretion disks. Prior to the addition of radiation transfer in simulations, therewas usually no explicit dependence on the black hole mass in the simulations. This isbecause the bare MHD equations can be arbitrarily rescaled by an explicit choice of lengthscale or characteristic density. Choosing the characteristic length scale specified the mass,and choosing the density then specified the mass accretion rate. Hence, one could scalethe simulation to the supermassive black hole regime, but there was nothing to distinguishthe dynamics or thermodynamics of the rescaled simulations from that of stellar massblack holes. However, this rescaling freedom is lost when radiation is added, because thedependence of opacities and emissivities on temperature and density enforces an explicitchoice of length and density scales. In this sense, radiation transfer is fundamental to thedistinction between supermassive and stellar mass black holes. Although we have learneda great deal about accretion flows from pure MHD simulations, they cannot tell us aboutthe differences between AGN flows and XRB flows, so we choose to focus on radiativesimulations in this section of the review.The first sets of radiation MHD simulations of accretion flows utilized the FLD approxi-mation (Levermore & Pomraning 1981; Turner & Stone 2001) to perform local shearing boxsimulations (e.g. Turner et al. 2003; Turner 2004; Hirose, Krolik & Stone 2006). With theexception of Turner (2004), these simulations focused on patches of disk around ∼ M (cid:12) black holes, where the radiation-to-gas pressure ratio tended to be lower. At the same time,efforts were underway to use FLD to study the global structure of accretion flows. Initially,this took the form of two-dimensional, viscous hydrodynamic flows (Ohsuga et al. 2005),but was eventually generalized to MHD treatments (Ohsuga & Mineshige 2011). Thesecalculations were mainly aimed at studying the variation of the flow with accretion rateand also focused on 10 M (cid:12) black holes. These radiative MHD shearing box simulations were particularly interesting for testing thepredictions of thermal instability in radiation pressure-dominated accretion flows (Shakura& Sunyaev 1976). The initial results from the Zeus FLD simulations confirmed that gaspressure-dominated disks were stable (Hirose, Krolik & Stone 2006), as expected. Moreintriguingly, simulations found evidence for thermal stability even when radiation pressuredominates gas pressure (Turner 2004; Hirose, Krolik & Blaes 2009). However, these resultswere contradicted by Athena simulations utilizing the VET method (Jiang, Stone & Davis2013), which generally found runaway heating or collapse. Simulations using the AthenaFLD module also showed indications of thermal runaway, so the discrepancy was not purelya result of the radiation transfer algorithm. Instead, Jiang, Stone & Davis (2013) attributethe discrepancy to the radial width of the Zeus simulation box being too small, findingthat the simulations became more stable as they narrowed the radial dimensions of thesimulation domain.Despite reaching different conclusions about stability, these simulations all found thatthe disks could survive for many thermal timescales, which is indicative of a possible robust-ness against the standard thermal instability. Although the relation of the time-averagedstress in MHD simulations is roughly consistent with the α -disk assumption, the shearingboxes differed from the standard model assumptions in important ways. Magnetic pressure
424 Davis & Tchekhovskoy s not included in standard disk models, but in simulations it plays an important role insupporting the surface regions against gravity. Thermodynamics is more complicated inMHD runs, with a greater fraction of dissipation occurring near or above the photosphere,leading the density to drop faster with height than in standard disk models. And the timedependence of the relation between central pressure and stress is less direct, so that in-creases (decreases) in radiation pressure did not immediately result in increases (decreases)in the scale height or dissipation (Hirose, Krolik & Blaes 2009).Turner (2004) was the sole simulation among those above that was run for supermassiveblack hole conditions, but it only included free-free (bremsstrahlung) and electron-scatteringopacity. Jiang, Davis & Stone (2016) explored the role of opacity in such environments,using the Athena VET module to simulate several shearing boxes for conditions around a 5 × M (cid:12) black hole while including OPAL opacities. Despite much higher ratios of radiationto gas pressure compared with the lower-mass black holes discussed by Jiang, Stone & Davis(2013), the resulting simulations appeared to be thermally stable. In contrast, simulationsperformed with the same parameters, but only including free-free and electron-scatteringopacity, all showed runaway collapse on very short timescales, much faster than in the lowerblack hole mass runs. Jiang, Davis & Stone (2016) attribute the enhancement in stabilityto a combination of two effects. First, UV opacities introduce an anticorrelation in theoptical depth to the midplane and central temperature, opposite to the electron scattering-dominated assumption. Second, the opacity drives convection that isn’t present in theprevious simulations, which enhances the nonradiative fluxes carried in the simulations.However, the dependence of shearing box simulations on box size, combined with time-dependent global modeling (Grz¸edzielski, Janiuk & Czerny 2017) strongly motivates furtherstudy of UV opacity effects in global numerical simulations. Local (shearing box) simulations are useful for studying accretion physics at high resolution,but one ultimately requires global simulations to obtain a more complete picture of theaccretion flow. Furthermore, the apparent dependence of thermal instability on radial width(or aspect ratio) in shearing box simulations motivates global simulations for which thisartificial constraint is removed. As already noted, global simulations of black hole accretionhave advanced significantly over the past two decades, advancing from pseudo-NewtonianMHD (Hawley, Gammie & Balbus 1995) to full GRMHD (Gammie, McKinney & T´oth2003), sometimes employing ad hoc cooling function to maintain a thin disk (e.g., Pennaet al. 2010; Noble, Krolik & Hawley 2010). Without radiation transfer and realistic opacities,such simulations can be scaled to arbitrary masses. Therefore, these simulations apply toboth the solar mass and the supermassive black hole regime. However, in practice, suchsimulations have been primarily compared with XRBs (e.g. Kulkarni et al. 2011; Schnittman,Krolik & Noble 2013).For these reasons, AGN specific MHD numerical simulations are relatively recent andrelatively few in number. Nevertheless, there are several lines of related inquiry that areworth summarizing. Although their simulation was calibrated for solar mass black holes,Mishra et al. (2016) performed the first radiation GRMHD simulation of radiation pressure-dominated global disk. Consistent with the earlier shearing box simulations results, theyfound a runaway collapse of the disk, approximately on the thermal timescale. In contrast,S¸adowski (2016) performed global radiation GRMHD simulations in a similar regime that • Disks and Jets in AGNs 425 as stable to thermal and inflow instabilities.Thus, we have two sub-Eddington accretion rate simulations that provide different pre-dictions about the stability of accretion flows in this radiation pressure dominated regime.S¸adowski (2016) attribute the stability of their simulation to magnetic pressure from toroidalmagnetic fields, which is broadly consistent with predictions of earlier analytic models thatsuggested magnetic support could impact stability (Begelman & Pringle 2007; Begelman &Silk 2017). The difference between the S¸adowski (2016) and Mishra et al. (2016) results canthen be attributed to the magnetic field geometry. The required toroidal fields are producedby threading the disk with a large poloidal flux that then gets wound up by the shear inthe accretion flow; this is a picture that is broadly supported by other recent simulationresults (Mishra et al. 2020), including some in the AGN regime discussed below (Jiang et al.2019). These results suggest that magnetic support can stabilize accretion flows, as longas sufficiently large poloidal magnetic fields can be accreted into the inner regions of theaccretion disk.Radiation hydrodynamic simulations of line-driven AGN outflows (e.g., Proga, Stone &Kallman 2000) have also been used to model properties of broad absorption line quasars.Because these calculations treat the accretion disk and its emission as a fixed boundary con-dition, they cannot evolve the disk self-consistently, but it is notable that some simulationshave mass outflow rates that are a substantial fraction of the assumed inflow rate. Takentogether, these simulations suggest that the impact of atomic opacities and the enhancedradiation pressure in AGNs may lead to distinct differences from XRB accretion flows.Radiation GRMHD simulations of the radiatively inefficient accretion flow regime havealso advanced significantly in the past several years, motivated by the imaging of Sgr A*and M87 with the Event Horizon Telescope (Ryan et al. 2018; Event Horizon TelescopeCollaboration et al. 2019). Because such flows are relatively optically thin, they are wellsuited to Monte Carlo-based algorithms for radiation transfer, such as the BHLIGHT code(Ryan, Dolence & Gammie 2015). However, the standard implementation of these algo-rithms has not traditionally scaled well to the optically thick regime of luminous accretionflows, so these techniques have not been applied to luminous AGNs.To date, only a few radiation MHD simulations have been performed with AGN massesand UV appropriate opacities. These AGN-specific simulations have been performed forsuper-Eddington (Jiang, Stone & Davis 2019) and sub-Eddington (Jiang et al. 2019) regimesusing the Athena++ radiation MHD code with a pseudo-Newtonian potential. The simu-lations are initialized with a torus centered at 50-80 gravitational radii, seeded with a weakmagnetic field. The torus is unstable to the MRI, generating turbulence that allows anaccretion flow to self-consistently form at radii interior to the torus.The generic properties of the accretion flow are broadly similar to the two-dimensionalviscous radiation hydrodynamic simulations (Ohsuga et al. 2005; Kitaki et al. 2018) as wellas radiation GRMHD simulations performed for solar mass black holes (McKinney, Dai &Avara 2015; S¸adowski et al. 2015; S¸adowski & Narayan 2016) or super-Eddington accretionin tidal disruption events (TDEs; Dai et al. 2018). At the highest accretion rates, theaccretion disk is geometrically thick, with strong radiation pressure-dominated outflows atvelocities of 10-30% of the speed of light. However, in contrast to previous simulations andslim accretion disk models, the radiative efficiencies can remain relatively high, as much as ∼
426 Davis & Tchekhovskoy t higher accretion rates ( ∼ L Edd /c ), radiative efficiencies drop to ∼ ∼
20 gravitational radii and last for a duration of many thermaltimescales without any evidence for thermal or inflow instabilities. As shown by S¸adowski(2016), the stability is attributed to magnetic pressure support. Although magnetic pressureis lower than radiation pressure, it has a smaller scale height and provides the dominantsupport of the gas against tidal gravity. The radiation pressure scale height is large becauseradiation pressure within the disk is roughly constant, due to the dissipation primarilyoccurring in surface layers with relatively low optical depth. Such magnetic pressure supporthas been proposed to aid in explaining the size and timescale discrepancies discussed inSection 2.1 (Dexter & Begelman 2019).The large dissipation in surface layers is the result of strong radiative viscosity. In thesimulation with 7% of the Eddington accretion rate, the radiation viscosity is actually thedominant mechanism for angular momentum transport, and more mass flows inward in thesurface layers than in the midplane of the disk where MRI turbulence dominates. The largeradiation viscosity is the result of a combination of large photon mean free paths in thesesurface regions and large radiation energy density.As shown in Figure 5 the large surface dissipation leads to strong temperature inver-sions above the disk, reminiscent of models of AGN coronae. The simulations have bothtemperatures high-enough to potentially produce the hard X-ray coronae and more mod-erate temperatures needed to explain the presence of soft X-ray excesses observed in manyAGNs (Walter & Fink 1993; Kubota & Done 2018). Intriguingly, the fraction of emissionin optically thin regions is higher in the lower accretion rate simulations, leading to highertemperatures consistent with the observed relation that optical to X-ray spectral indexesget harder as luminosity decreases (Steffen et al. 2006). Although these results are promis-ing, we must keep in mind that all these effects are occurring in close vicinity to the blackhole, where relativistic effects are important, so we must ultimately test these conclusionswith general relativistic calculations. • Disks and Jets in AGNs 427
10 20 r sin θ/r g − − − − r c o s θ / r g AGN0 .
07 10 − T g / × K r sin θ/r g − − − − r c o s θ / r g AGN0 . − T g / × K Figure 5: Time and azimuthally averaged spatial distributions of gas temperature T g forthe two simulations runs onto 5 × black holes from Jiang et al. (2019). The left andright panels have accretion rates of 7% ( AGN0.07) and 20% (
AGN0.2 ), respectively. Thegas temperature is ≈ − × K in the optically thick part of the disk but rapidlyincreases to 10 − K in the optically thin coronal regions. Note how the corona is moreradially extended in the lower accretion rate simulation than in the higher accretion-ratecase.Although OPAL opacities are used in the Athena++ simulations of super-Eddingtonand sub-Eddington disks, they have relatively little impact on the resulting dynamics. Thisis because the inner regions of the accretion rate, where the flow reaches a steady state, aretoo hot for there to be significant enhancements in the opacity above electron scattering.Future simulations need to focus on larger radii to see the effects of such opacities.
5. SIMULATIONS OF AGN JETS5.1. The Disk-jet Connection
Numerical simulations are a powerful tool for quantifying the power of relativistic jets andunderstanding the disk-jet connection. A particularly useful way of doing so is throughmeasuring jet energy efficiency, or the ratio of jet-to-accretion power, η = P jet / ˙ Mc . Nu-merical simulations by different groups found vastly different values of η : For instance, for a = 0 .
99, the efficiency ranged from ∼
3% (McKinney 2005) to ∼
15% (Hawley & Kro-lik 2006). Perhaps not surprisingly, the simulations sample the large parameter space ofallowed jet powers, which can range from zero power (no jet) to some maximum power.In the context of magnetically powered jets, this would map into a range from zero large-scale magnetic flux to some maximum value that could be measured by the simulations(see eq. 14). To determine the maximum jet power, Tchekhovskoy, Narayan & McKinney
428 Davis & Tchekhovskoy , r g was much larger than the typical tori of ∼ r g . The large torus size translatedinto a great amount of vertical magnetic flux contained in the system, thereby flooding theblack hole with the magnetic flux and determining the maximum jet efficiency achievablefor a given spin. As the simulation started, the MRI led to the accumulation of gas andmagnetic flux on the black hole (Figure 6b). This resulted in black hole magnetic fluxincreasing (Figure 6f) to the point that it overcame the force of gravity that keeps the gason an orbit around the black hole: The flux became strong enough to escape from the blackhole by tearing through the disk (Figure 6c,d). Such periods of flux expulsion followedby periods of relatively quiet accretion and flux accumulation lead to oscillations in bothmagnetic flux and jet efficiency (Figure 6f,g). In this MAD accretion regime, first simulatedin the non-relativistic context (Igumenshchev, Narayan & Abramowicz 2003; Igumenshchev2008), the magnetic flux on the black hole is as strong as possible and is dynamically impor-tant. Its time-average value, (cid:104) φ BH (cid:105) ≈
50, translates due to eq. (14) into the time-average jetefficiency, (cid:104) η (cid:105) ≈ (cid:104) η (cid:105) > h/r ∼ . −
1. Radiatively efficient disks, such as those thought to power luminous quasars, aremuch thinner, with h/r ∼ . L/ . L E ). Unexpectedly, such thin disks led to the pro-duction of jets with efficiency η = 20 −
50% (Liska et al. 2019c), defying the theoreticalexpectations that geometrically thin disks are incapable of dragging large-scale poloidalmagnetic flux inward and powering the jets (Lubow, Papaloizou & Pringle 1994). In thesesimulations, the thin disk was formed by rapid cooling of a thick disk, which is similar tothe physics expected to take place during the hard-to-soft spectral state transition. Thus,these powerful jets might represent the transient jets seen in XRBs. They could also berelated to radio-loud quasars that make up 10% of all quasars. It is possible that radioquasar is a transient stage in the life of a quasar, and the radio quasars are undergoing aspectral state transition from low-luminosity AGNs to radio-loud quasars. The simulationsshow signs that these jets are indeed transient phenomena:Tthe magnetic flux strength onthe black hole appears to undergo steady decline, indicating that the large-scale magneticflux diffuses out of the black hole and out to larger radii through the disk, suggesting animpending shutoff of the jets. However, longer-duration simulations are needed to verifythis hypothesis.
An important question arises: Do black holes in nature receive enough large-scale magneticflux to produce such powerful jets? Indications are that such strong magnetic fluxes andpowerful jets are present in a wide range of systems ranging from AGNs (Zamaninasabet al. 2014; Ghisellini et al. 2014; Nemmen & Tchekhovskoy 2015) to tidal disruption events(TDEs; Tchekhovskoy et al. 2014) and gamma-ray bursts (Tchekhovskoy & Giannios 2015).
Tidal disruptionevent (TDE): transient emissionfrom the disruptionand partial accretionof a star thatwanders close to ablack hole
What is the origin of this magnetic flux? Whereas in AGNs the ISM contains a substantialamount of large-scale poloidal magnetic flux, sufficient to flood the black hole if the accretion ••
What is the origin of this magnetic flux? Whereas in AGNs the ISM contains a substantialamount of large-scale poloidal magnetic flux, sufficient to flood the black hole if the accretion •• Disks and Jets in AGNs 429 ˙ M c (e) (cid:0) B H (f) t [r g /c] (cid:1) [ % ] (g) (cid:2) (cid:2)
10 0 10 20 x [r g ] (cid:2) (cid:2) z [ r g ] (cid:2) (cid:2) y [ r g ] (a) t =0 (cid:2) (cid:2)
10 0 10 20 x [r g ] (cid:2) (cid:2) (cid:2) (cid:2) (b) t =5785 (cid:2) (cid:2)
10 0 10 20 x [r g ] (cid:2) (cid:2) (cid:2) (cid:2) (c) t =16485 (cid:2) (cid:2)
10 0 10 20 x [r g ] (cid:2) (cid:2) (cid:2) (cid:2) (d) t =27015 (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) (cid:2) Figure 6: Large-scale vertical magnetic field accumulates at the center and forms a dynami-cally important magnetic field that obstructs the accretion flow and leads to a magneticallyarrested disk (Tchekhovskoy, Narayan & McKinney 2011). Color shows the logarithm ofdensity (see color bar). (a) The horizontal and vertical slices through the initial condition, aradially extended (from r in = 15 r g out to r out ∼ × r g ) equilibrium hydrodynamic torusembedded with a weak (plasma β (cid:38) t (cid:39) r g , it becomes dynamically important and partially leaves theblack hole, leading to a drop in φ BH . Soon thereafter, the system settles into a quasi steadystate: The black hole overeats magnetic flux before shedding its excess into the accretiondisk, after which the cycle repeats again. (g) Jet energy efficiency undergoes oscillationsthat mirror those of φ BH . The time-average value of the efficiency shown with the dashedline, η (cid:39)
430 Davis & Tchekhovskoy n an accretion disk that feeds the black hole for months to years. In one remarkable TDE,the formation of the putative disk was accompanied by the launching of a powerful jet(Bloom et al. 2011; Burrows et al. 2011; Levan et al. 2011; Zauderer et al. 2011) that wouldrequire orders-of-magnitude-more large-scale flux on the black hole than that available in thedisrupted star (Tchekhovskoy et al. 2014). One possibility is that the flux could be providedby a preexisting accretion disk (Kelley, Tchekhovskoy & Narayan 2014). However, couldthe large-scale magnetic flux be produced in situ, by the turbulent dynamo in the accretiondisk? The problem of large-scale poloidal flux dynamo is long-standing and interdisciplinary:The causes the of 11-year cycle of magnetic flux polarity flips on the Sun and the muchlonger dynamo cycle in Earth’s core (Cowling 1981) are still unclear. Numerical simulationsindicate that the formation of powerful jets requires preexisting large-scale vertical magneticflux in the flow (Beckwith, Hawley & Krolik 2008; McKinney, Tchekhovskoy & Blandford2012). More recently, evidence points to the large-scale poloidal flux dynamo producinglarge-scale magnetic flux at least in some circumstances; in fact, in thick disks, h/r ∼ . h/r ∼ .
02, so far do not show evidence of thisprocess (Liska et al. 2019a).
Many AGN jets reach relativistic velocities, with typical values of bulk Lorentz factor of γ ∼
10. Relativistic motion is not limited to AGNs, but is also observed in XRBs, with γ ∼ few, and gamma-ray bursts, with γ (cid:38) θ ∝ r − . ± . , over more than five orders of magnitudein distance before levelling off at θ (cid:39) .
01 rad ∼ . γ ∼ γ ∼ θ − (Tchekhovskoy, McKinney & Narayan 2008). The physics behind this accelerationlaw is simple: The jet accelerates such that its opening angle θ does not exceed its beamingangle, γ − , so the jet maintains transverse causal contact, so that the interior of the jetavoids running into the edge of the jet. However, this acceleration cannot last indefinitely:Eventually, the jet converts most of its internal energy into bulk motion kinetic energy, andthe acceleration levels off at γ ∼ γ max . Recently, Chatterjee et al. (2019) simulated thepropagation of black hole jets over five orders of magnitude in distance and showed that thedevelopment of the single-power law jet shape extending for as many orders of magnitudeoccurs naturally if the jets are collimated by an extended accretion flow. Furthermore, thejet accelerates in a way similar to observations. It is broadly agreed upon that jet emission is due to a combination of synchrotron and inverseCompton processes. However, presently, there is no agreement on what accelerates theemitting electrons. In particular, we do not understand what accelerates the electrons that • Disks and Jets in AGNs 431 roduce the radio emission that we use as an indicator of jet activity. The problem is evenmore severe in the origin of high-energy gamma-ray (GeV and TeV) emission, where it isunclear whether the emission is coming from near the event horizon or from large distances.There are several candidate processes potentially responsible for accelerating the emittingelectrons and, possibly, positrons. They can be accelerated in magnetospheric gaps nearthe black hole event horizon (e.g., Hirotani & Okamoto 1998; Broderick & Tchekhovskoy2015; Hirotani & Pu 2016; Ptitsyna & Neronov 2016; Levinson & Segev 2017; Chen, Yuan &Yang 2018; Parfrey, Philippov & Cerutti 2019), in 3D magnetic kink instabilities and shocksthat can be triggered by jets running into the ISM (e.g., Bromberg & Tchekhovskoy 2016;Tchekhovskoy & Bromberg 2016; Barniol Duran, Tchekhovskoy & Giannios 2017), and ininterface instabilities (e.g., Chatterjee et al. 2019). Heating due to magnetized turbulencein the exterior to the jet is also possible, so what we are seeing could be the emission fromthe sheath that surrounds the jet (and not from the jet itself Mo´scibrodzka & Falcke e.g.,2013; Ressler et al. e.g., 2017).
A standard approach to modeling black hole accretion is to consider an accretion disklying in the equatorial plane of the black hole. However, typically the disk midplane istilted relative to that of the black hole, necessitating the consideration of tilted disks. Earlysimulations of tilted thick disks with h/r ∼ . h/r < α , warps propagate in the disks viscously. Bardeen & Petterson (1975) predicted Bardeen-Pettersoneffect (BP): anabrupt change in thetilt of the disk dueto the interaction ofviscosity withrelativisticprecession that this would lead to the inner parts of the accretion disk aligning with the black holemidplane and separating from the outer, misaligned part of the disk by a smooth warp (theBardeen-Petterson effect or BP). However, at the time it was impossible to take into accountthe magnetized turbulence responsible for the accretion. And, though such alignment was
Smoothed particlehydrodynamics(SPH): a method forsolving fluiddynamics byfollowing thetrajectories ofparticles rather thanintegratingequations on a mesh seen in non-elativistic smoothed particle hydrodynamics (SPH) simulations (e.g., Nelson &Papaloizou 2000), until recently general relativistic magnetized simulations of tilted disks,as thin as h/r = 0 .
08, have shown no sign of the BP alignment (Zhuravlev et al. 2014). Liskaet al. (2019c) simulated the thinnest disk, with h/r = 0 .
03, tilted by 10 degrees relative toa rapidly spinning black hole, and found that the inner regions of the disk, r (cid:46) r g alignedwith the black hole. This remarkable discovery is the first demonstration of BP alignmentin a general relativistic numerical simulation of a magnetized accretion disk. Importantly,the duration of these simulations is long enough for the inner regions of the accretion diskto achieve inflow equilibrium and quasi-steady state. The BP alignment was also seen in
432 Davis & Tchekhovskoy horter nonrelativistic simulations that approximated the effects of black hole rotation viaadditional torques in MHD equations (Hawley & Krolik 2019) and did not include theeffects of apsidal precession, which can qualitatively affect the alignment process (Nealonet al. 2016).Figure 7: Tilted thin accretion disks tear up into individual, independently precessingsubdisks. Whereas at small radii the jets align with the inner subdisks, at large radiithe jets align with the outer subdisk. As the disks precess, jets run into subdisks thatget in the way and can not only expel them, as seen in the movie of the simulation(https://youtu.be/mbnG5 UTTdk), but also lead to energy dissipation that modifies theradial emission profile of the disk. Figure adapated from Liska et al. (2020).Interestingly, a disk tilted by 65 degrees got torn up into several individually precessingsubdisks (Liska et al. 2020), as seen in Fig. 7. Note that tilted disks in SPH simulations alsoshow tearing but with an important difference: Tilted disks get torn into a large number ofthin rings (Nixon et al. 2012). This difference likely stems from the effects of magnetizedturbulence and large-scale magnetic fields that hold the disk together differently than hy-drodynamic viscosity. The observational manifestations of torn disks are far reaching. Theinteractions between adjacent disks can lead to additional dissipation and emission at largeradii. Jets and radiation from the inner disk interacting with outer subdisks can contributeadditional heating. Not all gas flows between the adjacent subdisks, and an interesting frac-tion blown away owing to the interaction with the jets and winds. These factors can leadto modifications in the emission profile, potentially reducing the tension with the observeddisk sizes (see Fig. 1). Disk tearing can also lead to a wide range of variability, for example,sub-disks precessing through our line of sight might act as absorbers, from time to timedimming the emission of the central source.
As the jets emerge from the black hole’s sphere of influence, they run into the ISM. Atthis point, their behavior qualitatively changes. In the case of the M87 jet, the jet stopscollimating (Section 5.3) and starts to decelerate (Mertens et al. 2016). Jet interaction ••
As the jets emerge from the black hole’s sphere of influence, they run into the ISM. Atthis point, their behavior qualitatively changes. In the case of the M87 jet, the jet stopscollimating (Section 5.3) and starts to decelerate (Mertens et al. 2016). Jet interaction •• Disks and Jets in AGNs 433 igure 8: (a) Low-power AGN jets (blue-green) succumb to global magnetic instabilities,stall within their host galaxies, and inflate quasi-spherical cavities (yellow). (b) High-powerjets maintain their stability, leave their host galaxies, and form strong backflows. Thus,magnetic instabilities can be the key to resolving a 40-year long puzzle on the cause ofFanaroff & Riley (1974) morphological dichotomy of AGNs (Tchekhovskoy & Bromberg2016).with the ISM is poorly understood. For instance, there is no agreement on the origin ofthe Fanaroff & Riley (1974) morphological dichotomy of AGN jets: Fanaroff-Riley typeI (FRI) jets appear to develop instabilities early on and often disrupt inside the galaxy,whereas FRII jets appear well collimated and stably propagate to outside of the galaxy.This can have important consequences for their parent galaxy: Whereas FRII jets leavethe galaxy unscathed and deposit their energy outside the galaxy, FRI jets inject theirenergy into the galaxy and can significantly affect the star formation and dynamics of gas.Thus, it is important to understand the stability properties of the jets. Whereas it is rathereasy to reproduce the stable FRII jet morphology in numerical hydrodynamic simulations(Clarke, Norman & Burns 1986), reproducing the FRI morphology turned out to be muchmore difficult. Some of the possibilities include Kelvin-Helmholtz (KH) instabilities in theshear layers (e.g., Kaiser & Alexander 1997; Perucho, Mart´ı & Hanasz 2005; Meliani &Keppens 2009; Perucho et al. 2010), and mass entrainment from stellar winds (Komissarov1994; Perucho et al. 2014; Wykes et al. 2015). Studies of magnetized jet stability weremore promising; however, it turns out that magnetized jet stability sensitively dependson the degree of azimuthal winding of the magnetic field (e.g., Guan, Li & Li 2014), afree parameter whose value is poorly understood. This introduces an uncertainty into thefactors that control jet stability. Tchekhovskoy & Bromberg (2016) carried out large-scalesimulations of magnetized relativistic jets interacting with the ISM that for the first timewere launched via the rotation at the base, as they are launched in nature. The advantageof this approach is that it organically determines the degree of azimuthal winding andtherefore leads to the same stability properties of the jets as in nature. Figure 8 showsthat a change by two orders of magnitude in jet power leads to a drastic change in jetmorphology: More powerful jets are stable, whereas less powerful jets are unstable, andthis is in qualitative agreement with the FRI/II dichotomy. Future work including magneticfields and jet precession, which naturally emerges as a result of tilted accretion, will help to
434 Davis & Tchekhovskoy efine the models of jet feedback so that they can be used as subgrid models in cosmologicalsimulations.
6. SUMMARY AND OUTLOOK
Our understanding of the central engines of AGNs has greatly advanced over past severaldecades. We have a fiducial model for how accretion proceeds in the general relativisticspacetime of a supermassive black hole and gained an understanding of the central roleplayed by magnetic fields in the process of jet formation and angular momentum trans-port. This general picture of accretion-powered emission is well supported by the existingobservational constraints and we owe many of these insights (in part) to the application ofstate-of-the-art numerical simulations of GRMHD flows. This is particularly true for ourunderstanding of jet production and jet-disk interactions. We are beginning to understandthe effects of magnetic field geometry on the dynamics of accretion flows, the ability of disksto transport large-scale magnetic fields, and the disk-jet connection. The jets appear to bemore resilient than previously thought and defy the standard expectations: They emergein systems (a) without any large-scale vertical magnetic flux, (b) with a large disk tilt, and(c) with extremely small disk thickness.However, many challenges remain with our theoretical understanding of accretion disks.Longstanding theoretical inconsistencies remain unresolved, and detailed comparisons be-tween observations and theory yield significant discrepancies. What many viewed as promis-ing early agreement between theory and observations has faltered in the face of new ob-servational constraints: far-UV SEDs, size-scale constraints, and variability studies. Forthese reasons, we have chosen to focus much of this review on the limitations of our currentunderstanding, with the hope of providing a framework to motivate future research. Wehave argued that global numerical simulations will be essential for solving some of the mostimportant theoretical issues, including the following: • It is possible that gas fed to the accretion flows on large scales may be misaligned withthe black hole spin. We must understand the warping of the accretion flow in suchcases, including the possible impact of disk tearing, and its impact on disk emissionand reprocessing. • Standard accretion disk models predict the radiation pressure-dominated central re-gions of accretion flows in AGNs to be violently unstable to thermal and inflowinstabilities, yet there is no clear evidence of instability in the generic variability ofthese sources. Local MHD simulations cannot address inflow instability and provideconflicting results on the presence of thermal instability. • Local simulations and simple models suggest that opacities due to atomic transitionscould strongly modify the structure of AGN accretion flows or even drive outflows.Such dynamics are already seen in simulations of massive stellar envelopes. • Magnetic pressure support seems to be a promising candidate for stabilizing accretionflows in the radiation pressure-dominated regime. Such support also leads to thickeraccretion flows and both may change the vertical structure of accretion flows andincrease the reprocessing of emission from the inner regions of the flow.All the above questions require global simulations with large dynamic range. They may alsorequire resolving relatively thin accretion flows, making them computationally expensive.Many of these questions also require increasingly complex treatments of radiation transfer ••
Our understanding of the central engines of AGNs has greatly advanced over past severaldecades. We have a fiducial model for how accretion proceeds in the general relativisticspacetime of a supermassive black hole and gained an understanding of the central roleplayed by magnetic fields in the process of jet formation and angular momentum trans-port. This general picture of accretion-powered emission is well supported by the existingobservational constraints and we owe many of these insights (in part) to the application ofstate-of-the-art numerical simulations of GRMHD flows. This is particularly true for ourunderstanding of jet production and jet-disk interactions. We are beginning to understandthe effects of magnetic field geometry on the dynamics of accretion flows, the ability of disksto transport large-scale magnetic fields, and the disk-jet connection. The jets appear to bemore resilient than previously thought and defy the standard expectations: They emergein systems (a) without any large-scale vertical magnetic flux, (b) with a large disk tilt, and(c) with extremely small disk thickness.However, many challenges remain with our theoretical understanding of accretion disks.Longstanding theoretical inconsistencies remain unresolved, and detailed comparisons be-tween observations and theory yield significant discrepancies. What many viewed as promis-ing early agreement between theory and observations has faltered in the face of new ob-servational constraints: far-UV SEDs, size-scale constraints, and variability studies. Forthese reasons, we have chosen to focus much of this review on the limitations of our currentunderstanding, with the hope of providing a framework to motivate future research. Wehave argued that global numerical simulations will be essential for solving some of the mostimportant theoretical issues, including the following: • It is possible that gas fed to the accretion flows on large scales may be misaligned withthe black hole spin. We must understand the warping of the accretion flow in suchcases, including the possible impact of disk tearing, and its impact on disk emissionand reprocessing. • Standard accretion disk models predict the radiation pressure-dominated central re-gions of accretion flows in AGNs to be violently unstable to thermal and inflowinstabilities, yet there is no clear evidence of instability in the generic variability ofthese sources. Local MHD simulations cannot address inflow instability and provideconflicting results on the presence of thermal instability. • Local simulations and simple models suggest that opacities due to atomic transitionscould strongly modify the structure of AGN accretion flows or even drive outflows.Such dynamics are already seen in simulations of massive stellar envelopes. • Magnetic pressure support seems to be a promising candidate for stabilizing accretionflows in the radiation pressure-dominated regime. Such support also leads to thickeraccretion flows and both may change the vertical structure of accretion flows andincrease the reprocessing of emission from the inner regions of the flow.All the above questions require global simulations with large dynamic range. They may alsorequire resolving relatively thin accretion flows, making them computationally expensive.Many of these questions also require increasingly complex treatments of radiation transfer •• Disks and Jets in AGNs 435 uilt on top of already sophisticated algorithms for evolving GRMHD. Progress will almostcertainly require utilizing efficiently optimized algorithms that can harness the capabilitiesof the next generations of computing infrastructure (exascale and beyond). Fortunately,there has already been substantial recent progress in developing the numerical tools andalgorithms needed for the future, and we are optimistic these questions can be addressedin the next decade.
FUTURE ISSUES
1. Although local simulations have advantages for studying accretion disks at high res-olution, questions about convergence and dependence on simulation domain sizesand magnetically powered outflows strongly motivate global numerical simulationsof accretion flows. The large dynamic range required by such simulations requiresthe development of algorithms and codes that model the required physics while effi-ciently utilizing the latest advances in high-performance computing infrastructure.2. Multiscale simulations that bridge the scales of the galaxy and that of the blackhole will be needed to self-consistently determine the mass supply near the blackhole and the feedback of the black hole on the galaxy.3. Realistic treatments of radiation transfer are essential for modeling the radiation-dominated regions of accretion flows to resolve longstanding questions about stabil-ity and to study the effects of atomic opacities. The ability of scaling simulationsto different black hole masses that is present in most previous nonradiative simula-tions is no longer possible when radiation transfer is included, requiring AGN- andXRB-specific simulations.4. Future simulations will need to address the evident discrepancies between observa-tions of AGNs and standard disk models. This includes size discrepancies, flatteningof spectra in the UV and absence of edges, the origin of AGN continuum variability,and lack of generic evidence of instability.
DISCLOSURE STATEMENT
The authors are not aware of any affiliations, memberships, funding, or financial holdingsthat might be perceived as affecting the objectivity of this review.
ACKNOWLEDGMENTS
We thank the anonymous referees for many helpful suggestions that improved the qualityof this review. We are grateful to our many colleagues and collaborators, too numerous tolist here, who have enlightened us over the years and shaped our views on AGN, accretiondisks and jets.
LITERATURE CITED
Abramowicz MA, Czerny B, Lasota JP, Szuszkiewicz E. 1988. ApJ 332:646–658Balbus SA, Hawley JF. 1991. ApJ 376:214–233Balbus SA, Hawley JF. 1998.
Reviews of Modern Physics
436 Davis & Tchekhovskoy albus SA, Papaloizou JCB. 1999. ApJ 521:650–658Bardeen JM, Petterson JA. 1975. ApJ 195:L65Barniol Duran R, Tchekhovskoy A, Giannios D. 2017. MNRAS 469:4957–4978Baron D, Stern J, Poznanski D, Netzer H. 2016. ApJ 832:8Beckwith K, Hawley JF, Krolik JH. 2008. ApJ 678:1180–1199Beckwith K, Hawley JF, Krolik JH. 2009. ApJ 707:428–445Begelman MC, Pringle JE. 2007. MNRAS 375:1070–1076Begelman MC, Silk J. 2017. MNRAS 464:2311–2317Beskin VS, Nokhrina EE. 2006. MNRAS 367:375–386Blackburne JA, Pooley D, Rappaport S, Schechter PL. 2011. ApJ 729:34Blandford RD, Begelman MC. 1999. MNRAS 303:L1–L5Blandford RD, Payne DG. 1982. MNRAS 199:883–903Blandford RD, Znajek RL. 1977. MNRAS 179:433–456Bloom JS, Giannios D, Metzger BD, Cenko SB, Perley DA, et al. 2011.
Science
Monthly Notices of the Royal Astronomical Society ••
Monthly Notices of the Royal Astronomical Society •• Disks and Jets in AGNs 437 uan X, Li H, Li S. 2014. ApJ 781:48Hawley JF, Balbus SA. 2002. ApJ 573:738–748Hawley JF, Gammie CF, Balbus SA. 1995. ApJ 440:742Hawley JF, Krolik JH. 2006. ApJ 641:103–116Hawley JF, Krolik JH. 2019. ApJ 878:149Hawley JF, Smarr LL, Wilson JR. 1984. ApJ 277:296–311Hirose S, Krolik JH, Blaes O. 2009. ApJ 691:16–31Hirose S, Krolik JH, Stone JM. 2006. ApJ 640:901–917Hirotani K, Okamoto I. 1998. ApJ 497:563–572Hirotani K, Pu HY. 2016. ApJ 818:50Homayouni Y, Trump JR, Grier CJ, Shen Y, Starkey DA, et al. 2019. ApJ 880:126Hubeny I, Blaes O, Krolik JH, Agol E. 2001. ApJ 559:680–702Iglesias CA, Rogers FJ. 1996. ApJ 464:943Igumenshchev IV. 2008. ApJ 677:317–326Igumenshchev IV, Narayan R, Abramowicz MA. 2003. ApJ 592:1042–1059Janiuk A, Czerny B, Siemiginowska A. 2002. ApJ 576:908–922Jiang YF, Blaes O, Stone JM, Davis SW. 2019. ApJ 885:144Jiang YF, Cantiello M, Bildsten L, Quataert E, Blaes O, Stone J. 2018. Nature 561:498–501Jiang YF, Davis SW, Stone JM. 2016. ApJ 827:10Jiang YF, Green PJ, Greene JE, Morganson E, Shen Y, et al. 2017. ApJ 836:186Jiang YF, Stone JM, Davis SW. 2012. ApJS 199:14Jiang YF, Stone JM, Davis SW. 2013. ApJ 778:65Jiang YF, Stone JM, Davis SW. 2014a. ApJ 796:106Jiang YF, Stone JM, Davis SW. 2014b. ApJ 784:169Jiang YF, Stone JM, Davis SW. 2019. ApJ 880:67Kaiser CR, Alexander P. 1997. MNRAS 286:215–222Kelley LZ, Tchekhovskoy A, Narayan R. 2014. MNRAS 445:3919–3938Kitaki T, Mineshige S, Ohsuga K, Kawashima T. 2018. PASJ 70:108Kolykhalov PI, Sunyaev RA. 1984.
Advances in Space Research
Active galactic nuclei : from the central black hole to the galactic environment
Kubota A, Done C. 2018. MNRAS 480:1247–1262Kulkarni AK, Penna RF, Shcherbakov RV, Steiner JF, Narayan R, et al. 2011. MNRAS 414:1183–1194LaMassa SM, Cales S, Moran EC, Myers AD, Richards GT, et al. 2015. ApJ 800:144Laor A, Davis SW. 2014. MNRAS 438:3024–3038Laor A, Netzer H. 1989. MNRAS 238:897–916Lasota JP. 2001. New A Rev. 45:449–508Lawrence A. 2012. MNRAS 423:451–463Lawrence A. 2018.
Nature Astronomy
Science
Finite-Volume Methods for Hyperbolic Problems . Cambridge University PressLevermore CD, Pomraning GC. 1981. ApJ 248:321–334Levinson A, Segev N. 2017. Phys. Rev. D 96:123006Li ZY, Chiueh T, Begelman MC. 1992. ApJ 394:459Lightman AP, Eardley DM. 1974. ApJ 187:L1Liska M, Chatterjee K, Tchekhovskoy Ae, Yoon D, van Eijnatten D, et al. 2019a. arXiv e-prints :arXiv:1912.10192Liska M, Hesp C, Tchekhovskoy A, Ingram A, van der Klis M, Markoff S. 2018. MNRAS 474:L81–
438 Davis & Tchekhovskoy arXiv e-prints
Liska M, Hesp C, Tchekhovskoy A, Ingram A, van der Klis M, et al. 2020. MNRASLiska M, Tchekhovskoy A, Ingram A, van der Klis M. 2019c. MNRASLiska M, Tchekhovskoy A, Quataert E. 2020. MNRAS 494:3656–3662Liu G, Liu M. 2003.
Smoothed Particle Hydrodynamics: A Meshfree Particle Method
Lubow SH, Papaloizou JCB, Pringle JE. 1994. MNRAS 267:235–240Lynden-Bell D. 1969. Nature 223:690–694Lynden-Bell D, Pringle JE. 1974. MNRAS 168:603–637MacLeod CL, Ross NP, Lawrence A, Goad M, Horne K, et al. 2016. MNRAS 457:389–404Masaki I. 1971. PASJ 23:425McKernan B, Ford KES, Lyra W, Perets HB. 2012. MNRAS 425:460–469McKinney JC. 2005. ApJ 630:L5–L8McKinney JC, Dai L, Avara MJ. 2015. MNRAS 454:L6–L10McKinney JC, Gammie CF. 2004. ApJ 611:977–995McKinney JC, Tchekhovskoy A, Blandford RD. 2012. MNRAS 423:3083–3117McKinney JC, Tchekhovskoy A, Blandford RD. 2013.
Science
Foundations of radiation hydrodynamics
Miralda-Escud´e J, Kollmeier JA. 2005. ApJ 619:30–40Mishra B, Begelman MC, Armitage PJ, Simon JB. 2020. MNRAS 492:1855–1868Mishra B, Fragile PC, Johnson LC, Klu´zniak W. 2016. MNRAS 463:3437–3448Misner CW, Thorne KS, Wheeler JA. 1973.
Gravitation
Mocz P, Vogelsberger M, Hernquist L. 2014. MNRAS 442:43–55Morgan CW, Hyer GE, Bonvin V, Mosquera AM, Cornachione M, et al. 2018. ApJ 869:106Morgan CW, Kochanek CS, Morgan ND, Falco EE. 2010. ApJ 712:1129–1136Mo´scibrodzka M, Falcke H. 2013. A&A 559:L3Murray N, Chiang J, Grossman SA, Voit GM. 1995. ApJ 451:498Nakamura M, Asada K. 2013. ApJ 775:118Narayan R, Igumenshchev IV, Abramowicz MA. 2003. PASJ 55:L69–L72Narayan R, McClintock JE. 2012. MNRAS 419:L69–L73Narayan R, Yi I. 1994. ApJ 428:L13–L16Nayakshin S, Rappaport S, Melia F. 2000. ApJ 535:798–814Nealon R, Nixon C, Price DJ, King A. 2016. MNRAS 455:L62–L66Nelson RP, Papaloizou JCB. 2000. MNRAS 315:570–586Nemmen RS, Tchekhovskoy A. 2015. MNRAS 449:316–327Netzer H, Trakhtenbrot B. 2014. MNRAS 438:672–679Nixon C, King A, Price D, Frank J. 2012. ApJ 757:L24Noble SC, Krolik JH, Hawley JF. 2010. ApJ 711:959–973Novikov ID, Thorne KS. 1973.
Astrophysics of black holes. In Black Holes (Les Astres Occlus) , eds.C Dewitt, BS DewittOhsuga K, Mineshige S. 2011. ApJ 736:2Ohsuga K, Mineshige S, Mori M, Kato Y. 2009. PASJ 61:L7–L11Ohsuga K, Mori M, Nakamoto T, Mineshige S. 2005. ApJ 628:368–381Pan Z, Yu C. 2015. ApJ 812:57Parfrey K, Philippov A, Cerutti B. 2019.
Phys. Rev. Lett. • Disks and Jets in AGNs 439 erucho M, Mart´ı JM, Hanasz M. 2005. A&A 443:863–881Perucho M, Mart´ı JM, Laing RA, Hardee PE. 2014. MNRAS 441:1488–1503Pessah ME, Chan Ck, Psaltis D. 2007. ApJ 668:L51–L54Piran T. 1978. ApJ 221:652–660Porth O, Chatterjee K, Narayan R, Gammie CF, Mizuno Y, et al. 2019. ApJS 243:26Proga D, Kallman TR. 2004. ApJ 616:688–695Proga D, Stone JM, Kallman TR. 2000. ApJ 543:686–696Ptitsyna K, Neronov A. 2016. A&A 593:A8Remillard RA, McClintock JE. 2006. ARA&A 44:49–92Ressler SM, Tchekhovskoy A, Quataert E, Gammie CF. 2017. MNRAS 467:3604–3619Richards GT, Hall PB, Vanden Berk DE, Strauss MA, Schneider DP, et al. 2003. AJ 126:1131–1147Russell DM, Gallo E, Fender RP. 2013. MNRAS 431:405–414Ryan BR, Dolence JC, Gammie CF. 2015. ApJ 807:31Ryan BR, Gammie CF, Fromang S, Kestener P. 2017. ApJ 840:6Ryan BR, Ressler SM, Dolence JC, Gammie C, Quataert E. 2018. ApJ 864:126Rybicki GB, Lightman AP. 1986.
Radiative Processes in Astrophysics
Sakimoto PJ, Coroniti FV. 1981. ApJ 247:19–31Salpeter EE. 1964. ApJ 140:796–800S¸adowski A. 2016. MNRAS 459:4397–4407S¸adowski A, Narayan R. 2016. MNRAS 456:3929–3947S¸adowski A, Narayan R, Tchekhovskoy A, Abarca D, Zhu Y, McKinney JC. 2015. MNRAS 447:49–71S¸adowski A, Narayan R, Tchekhovskoy A, Zhu Y. 2013. MNRAS 429:3533–3550Schmidt M. 1963. Nature 197:1040Schnittman JD, Krolik JH, Noble SC. 2013. ApJ 769:156Shafee R, McKinney JC, Narayan R, Tchekhovskoy A, Gammie CF, McClintock JE. 2008. ApJ687:L25Shakura NI, Sunyaev RA. 1973. A&A 24:337–355Shakura NI, Sunyaev RA. 1976. MNRAS 175:613–632Shang Z, Brotherton MS, Green RF, Kriss GA, Scott J, et al. 2005. ApJ 619:41–59Shi JM, Stone JM, Huang CX. 2016. MNRAS 456:2273–2289Shields GA. 1978. Nature 272:706–708Shull JM, Stevans M, Danforth CW. 2012. ApJ 752:162Simon JB, Beckwith K, Armitage PJ. 2012. MNRAS 422:2685–2700Sincell MW, Krolik JH. 1997. ApJ 476:605–619Slone O, Netzer H. 2012. MNRAS 426:656–664Springel V. 2010. MNRAS 401:791–851Steffen AT, Strateva I, Brandt WN, Alexander DM, Koekemoer AM, et al. 2006. AJ 131:2826–2842Steiner JF, McClintock JE, Narayan R. 2013. ApJ 762:104Stone JM, Hawley JF, Gammie CF, Balbus SA. 1996. ApJ 463:656Tchekhovskoy A. 2015.
Launching of Active Galactic Nuclei Jets . In
The Formation and Disruptionof Black Hole Jets , eds. I Contopoulos, D Gabuzda, N Kylafis, vol. 414 of
Astrophysics and SpaceScience Library
Tchekhovskoy A, Bromberg O. 2016. MNRAS 461:L46–L50Tchekhovskoy A, Giannios D. 2015. MNRAS 447:327–344Tchekhovskoy A, McKinney JC, Narayan R. 2008. MNRAS 388:551–572Tchekhovskoy A, Metzger BD, Giannios D, Kelley LZ. 2014. MNRAS 437:2744–2760Tchekhovskoy A, Narayan R, McKinney JC. 2010. ApJ 711:50–63Tchekhovskoy A, Narayan R, McKinney JC. 2011. MNRAS 418:L79–L83Tombesi F, Cappi M, Reeves JN, Palumbo GGC, Yaqoob T, et al. 2010. A&A 521:A57Turner NJ. 2004. ApJ 605:L45–L48
440 Davis & Tchekhovskoy urner NJ, Stone JM. 2001. ApJS 135:95–107Turner NJ, Stone JM, Krolik JH, Sano T. 2003. ApJ 593:992–1006Walter R, Fink HH. 1993. A&A 274:105Weymann RJ, Morris SL, Foltz CB, Hewett PC. 1991. ApJ 373:23–53White CJ, Stone JM, Quataert E. 2019. ApJ 874:168Wilson JR. 1972. ApJ 173:431Wykes S, Hardcastle MJ, Karakas AI, Vink JS. 2015. MNRAS 447:1001–1013Yuan F, Narayan R. 2014. ARA&A 52:529–588Zamaninasab M, Clausen-Brown E, Savolainen T, Tchekhovskoy A. 2014. Nature 510:126–128Zauderer BA, Berger E, Soderberg AM, Loeb A, Narayan R, et al. 2011. Nature 476:425–428Zhu Z, Stone JM. 2018. ApJ 857:34Zhuravlev VV, Ivanov PB, Fragile PC, Morales Teixeira D. 2014. ApJ 796:104 ••