11 Numerical Simulations

In simulating accretion disks around black holes, there are a number of challenging issues. First, there is quite a lot of physics involved: relativistic gravity, hydrodynamics, magnetic fields, and radiation being the most fundamental. Then there is the issue that accretion disks are inherently multi-dimensional objects. The computational expense of including extra dimensions in a numerical simulation is not a trivial matter. Simply going from one to two dimensions (still assuming axisymmetry for a disk) increases the computational expense by a very large factor ( or more). Going to three-dimensions and relaxing all symmetry requirements increases the computational expense yet again by a similar factor. Simulations of this size have only become feasible within the last decade and still only with a subset of the physics one is interested in and usually with a very limited time duration.

Another hindrance in simulating accretion disks is the very large range of scales that can be present. In terms of a grid based code, a disk with a scale height of requiring zones to resolve in the vertical direction at some radius , would require something of the order zones to cover each factor of that is treated in the radial direction. The azimuthal direction in a full three-dimensional simulation would require a comparable number of zones to what is used in the radial direction. Given that a very large calculation by today’s standards is zones, we can see that treating very thin disks () in a full three-dimensional simulation, even with a very modest number of zones in the vertical direction () would take all the resources one could muster. However, two-dimensional (axisymmetric) simulations of thin disks and three-dimensional simulations of thick, slim, or moderately thin disks have now become quite common.

The approach of numericists in many ways parallels that of theorists, so we structure this section much the same as the first part of this Living Review. That is, after a brief introduction to numerical techniques, we discuss how the various components of the matter description (à la Section 3) are implemented in numerical simulations. We then review a few special cases illustrating how analytic models (Sections 4 – 7) and numerical simulations can complement one another. We finish this numerical section with a few topics of special interest.

11.1 Numerical techniques

11.1.1 Computational fluid dynamics codes

There are many numerical codes available today that include relativistic hydrodynamics or MHD that are, or can be, used to simulate accretion disks. A partial list includes: Cosmos++ [18], ECHO [74], HARM [105], and RAISHIN [208] plus the codes developed by Koide et al. [152], De Villiers and Hawley [71], Komissarov [154], Antón et al. [19], and Anderson et al. [16]. This represents tremendous growth in the field as prior to 1998 there were only a couple such codes, largely derived from the early work of Wilson [318] and later developments by Hawley [126, 125]. Those early codes were only available to a handful of researchers. This was simply a reflection of the fact that the computational resources available to most researchers at that time were insufficient to make use of such codes, so there was little incentive to develop or acquire them. This has clearly changed.

Today there are primarily two types of numerical schemes being used to simulate relativistic accretion disks, differentiated by how they treat discontinuities (shocks) that might arise in the flow: the artificial viscosity scheme, still largely based upon formulations developed by Wilson [318]; and Godunov-type approaches, using exact or approximate Riemann solvers. Both are based on finite difference representations of the equations of general relativistic hydrodynamics, although the latter tend to also incorporate finite-volume representations. The artificial viscosity schemes have the advantages that they are more straightforward to implement and easier to extend to include additional physics. More importantly, they are computationally less expensive to run than Godunov schemes. Furthermore, recent work by Anninos and collaborators [17, 18] has shown that variants of this scheme can be made as accurate as Godunov schemes even for ultrarelativistic flows, which historically was one of the principal weaknesses of the artificial viscosity approach. Godunov schemes, on the other hand, appreciate the advantage that they are fully conservative, and therefore, potentially more accurate. They also require less tuning since there are no artificial viscosity parameters that need to be set for each problem. More thorough reviews of these two methods, with clear emphasis on the High-Resolution Shock-Capturing variant of the Godunov approach, are provided in the Living Reviews by Martí & Müller [182] and Font [97]. Other numerical schemes, such as smooth-particle hydrodynamics (SPH) and (pseudo-)spectral methods are less well developed for work on relativistic accretion disks.

11.1.2 Global vs. shearing-box simulations

Along with settling on a numerical scheme, a decision must also be made whether or not to try to treat the disk as a whole or to try to understand it in parts. The latter choice includes “shearing-box” simulations, the name coming from the type of boundary conditions one imposes on the domain to mimic the shear that would be present in a real disk [121]. The obvious advantage of treating the disk in parts is that you circumvent the previously noted problem of the large range in scales in the disk by simply ignoring the large scales. Instead one treats a rectangular volume generally no larger than a few vertical scale heights on a side and in some cases much smaller. In this way, for a moderate number of computational zones, one can get as good, or often much better, resolution over the region being simulated than can be achieved in global simulations. The obvious disadvantage is that all information about what is happening on larger scales is lost. By construction, no structures larger than the box can be captured and the periodic boundaries impose some artificial conditions on the simulations, such as no background gradients in pressure or density and no flux of material into or out of the box, meaning the surface density remains fixed. Perhaps most important in the context of a review on relativity is that the box is treated as a local patch within the disk, with a scale smaller than that of the curvature of the metric. Thus most general relativistic effects are not, and by construction need not be, included in shearing box simulations. Since this is a Living Review in Relativity, we will not dwell much on this class of simulations. Still, there are some significant highlights that should be mentioned.

By far the most extensive work using shearing box simulations has been aimed at better understanding the magneto-rotational instability (MRI). Starting from the earliest “proof-of-concept” simulations [121, 296], shearing box simulations have been used to demonstrate various properties of the saturated state of MRI turbulence [271, 272, 173, 283] and much about the vertical structure of MRI turbulent disks [46]. Shearing box simulations have also proven valuable in studying radiation-dominated disks. For example, Turner [307] studied the vertical structure of radiation dominated accretion disks and Hirose et al. [129] studied the gas-pressure dominated case, both including radiation using flux-limited diffusion. Hirose and collaborators subsequently used shearing box simulations to demonstrate that radiation-dominated disks are thermally stable [128], though they had long been held not to be [176, 280, 246]. As we said, though, there are limits imposed by the finite size of the shearing box. For instance, the viscous, Lightman–Eardley instability [176] can only be studied through global simulations.

11.2 Matter description in simulations

The minimum physics required for a global simulation of an accretion disk are gravity and hydrodynamics (assuming the disk is dense enough for the continuum approximation to hold). Since many disks have masses that are small compared to the mass of the central compact object, the self-gravity of the disk can often be ignored. Therefore, in the next three sections, gravity will simply mean that of the central black hole.

11.2.1 Hydrodynamics + gravity

The first researcher to develop and use numerical algorithms for simulating relativistic accretion flows was Wilson [316], who considered the spherical infall of material with a non-zero specific angular momentum toward a Kerr black hole using the full metric, although restricted to two spatial dimensions. Wilson was able to confirm the additional centrifugal support that the infalling material experienced due to the rotating black hole. For sufficiently high values of angular momentum, the material could not immediately accrete onto the black hole, instead forming a fat disk of the type described in Section 4.

In many ways, Wilson’s pioneering work was at least a decade ahead of its time and there was consequently a lull in activity until Hawley and Smarr collaborated with Wilson to revive his work [126, 125]. As an example of how analytic and numerical work can complement each other, it is worth noting that one of the test cases they used for their new two-dimensional relativistic hydrodynamics code was based on the analytic theory for relativistic thick disks, which had been worked out in the time since Wilson’s original simulations. Using the analytic theory, they constructed a series of disks with different (constant) specific angular momenta, from to . The case yielded a static solution as expected and confirmed the ability of their code to accurately evolve such a solution in multi-dimensions over a dynamical time. The cases showed greater time variability and illustrated the power of direct numerical simulations to extend our understanding of black hole accretion.

11.2.2 Magnetohydrodynamics + gravity

Magnetic fields can play many important roles in relativistic accretion disks, from providing local viscous stresses through turbulence that results from the magnetorotational instability (Section 8.2), to providing a mechanism for launching and confining jets (Section 10). Thus, the inclusion of magnetic fields in numerical simulations of relativistic accretion disks is important. Although the techniques for doing this were first described in [317], relatively little work was done in this area until fairly recently. Perhaps the two most important contributions so far from relativistic MHD simulations of accretion disks have been: 1) elucidating the behavior of MRI turbulent disks in the vicinity of black holes, and 2) exploring the many interrelations between magnetic fields and relativistic jets.

In at least one case, the inclusion of MHD + gravity is sufficient to adequately capture the dynamics of a real black hole accretion disk. Sgr A*, the black hole at the center of the Milky Way galaxy has such an anemically low luminosity [194] that numerical simulations that ignore radiation can safely be applied [76], as has been done by many authors [112, 227, 211, 75]. In most other cases, radiative processes must somehow be accounted for.

11.2.3 Radiation-Magnetohydrodynamics + Gravity

Probably the most glaring shortcoming of almost all numerical simulations of accretion disks and many other phenomena in astrophysics to date is the unrealistic treatment of radiation, which is most often simply ignored. This is not due to a lack of appreciation of its importance on the part of numericists, but simply a reflection of the fact that there are very few efficient ways to treat radiation computationally in multiple dimensions.

The optically-thin limit is an exception. In this case, the radiative cooling only enters the energy and momentum conservation equations as a source term

where is the cooling function of the gas (Section 3.4). The subsequent propagation of the radiation through the matter can be ignored or post-processed. Even so, including radiative cooling can be critically important to capturing the true dynamics of the disk [101, 76], as shown in Figure 14.

Figure 14: Pseudo-color plots of with contours of from four different general relativistic MHD simulations. The simulations all begin with the same initial conditions, but have different energy conservation and cooling treatments: The upper-left panel conserves internal energy and ignores cooling; the upper-right panel conserves internal energy and includes cooling; the lower-left panel conserves total energy and ignores cooling; the lower-right panel conserves total energy and includes cooling. The very different end states illustrate the importance of properly capturing thermodynamic processes. Image reproduced by permission from [101], copyright by AAS.

In the optically-thick limit, some progress can be made by employing the flux-limited diffusion approximation [174]. Global simulations of disks including a flux-limited diffusion treatment of radiation (but using pseudo-Newtonian gravity) are now being done by Ohsuga and collaborators [232, 231]. A few steps toward the goal of relativistic radiation MHD simulations of black hole accretion disks have also been taken in recent years. A method for treating optically thick accretion using a conservative, Godunov scheme was developed by Farris and collaborators [91]. The same basic method has now been used to examine both Bondi [100] and Bondi–Hoyle [326, 266] accretion. Simulations of accretion disks, though, must await the generalization of this method to treat radiation both in the optically thick and thin limits.

11.2.4 Evolving GRMHD

In most simulations of accretion disks around black holes, the self-gravity of the disk is ignored. In many cases this is justified as the mass of the disk is often much smaller than the mass of the black hole. This is also much simpler as it allows one to treat gravity as a background condition, either through a Newtonian potential or a relativistic metric (the so-called Cowling approximation in relativity). However, there are plausible astrophysical scenarios in which this approximation is not valid. Two of the more interesting are: 1) a tidally disrupted neutron star accreting onto a stellar-mass black hole; and 2) an overlying stellar envelope accreting onto a nascent black hole during the final dying moments of a massive star. Interestingly these scenarios are currently the most popular models of gamma-ray bursts [82, 319], which are the most powerful explosions in the Universe since the Big Bang and are observed at a rate of about one per day.

Treating a self-gravitating disk in general relativity requires a code that can simultaneously evolve the spacetime metric (by solving the Einstein field equations) and the matter, magnetic, and radiation fields. Codes capable of doing this have very recently become available, divided into the following “flavors”: self-gravity + hydro [24]; self-gravity + MHD [81, 286, 110, 88]; and self-gravity + radiation MHD  [91]. More detailed discussion of these methods and their application to problems in astrophysics is provided in the Living Review by Font [97].

The importance of evolving the gravitational field in the case of black hole accretion is well illustrated in the case of the runaway instability, which we described in Section 8.1.2. That instability has now been simulated directly, including evolving the spacetime metric [209, 210, 160]. These studies have confirmed earlier results that the runaway instability does not develop.

Another application of these types of codes that is relevant to our review is the reexamination of the Papaloizou–Pringle instability (discussed in Section 8.1.1). Through this instability, a massive, self-gravitating Polish doughnut orbiting around a black hole could become a strong emitter of large amplitude, quasiperiodic gravitational waves [149].

11.3 Polish doughnuts (thick) disks in simulations

As we already mentioned, it is computationally less expensive to simulate thick disks than thin, so it comes as no surprise that most of the numerical work for nearly four decades, starting with the work of Wilson [316], has been on thick disks. In this section we touch on a few highlights from this area of research

Some of this work has focused on confirming predictions from analytic theory. Igumenshchev and Beloborodov [137] used two-dimensional relativistic hydrodynamic simulations to confirm that: (i) The structure of the innermost disk region strongly depends on the black hole spin, and (ii) the mass accretion rate is proportional to the size of the energy gap between the inner edge of the disk and the cusp (see Figure 15), as predicted by Kozłowsk and collaborators [162].

Figure 15: Time-average mass accretion rate as a function of the energy gap for models with (circles), (squares), and (triangles). The bars show the variability of . The lines represent the predicted dependencies , where is the adiabatic index. Image reproduced by permission from [137], copyright by RAS.

Hawley [116, 117] studied the nonlinear evolution of the non-axisymmetric Papaloizou–Pringle instability (Section 8.1.1). He discovered that for radially slender tori, the principal mode of the PPI saturates with the formation of ellipsoidal overdensity regions he termed “planets” (see Figure 16 for an example). In cases where several wavenumbers were unstable, multiple planets would begin to form but then nonlinear mode-mode coupling would cause them to merge. For wide tori, particularly when the tori extended beyond their Roche limits, overflowed their cusps, and accreted into the hole, the PPI saturated at low amplitude or was prevented from growing altogether. The initial work only considered Schwarzschild black holes, but all of the results were later qualitatively confirmed for Kerr black holes [69].

Figure 16: Equatorial slice through hydrodynamic tori at saturation of the Papaloizou–Pringle instability showing formation of significant non-axisymmetric () overdensity clumps. The density contours are linearly spaced between and 0.0. This figure represents models A3p (left) and B3r (right) of [69]. Image reproduced by permission, copyright by AAS.

Another interesting note is that most of the global numerical simulations of MRI turbulent disks have used a Polish doughnut as the starting condition. This was first done by De Villiers and Hawley [70] (animations can be viewed at [68, 115]) and later by a number of other groups [105, 18]. As expected from linear analysis and earlier Newtonian simulations, the Maxwell stresses produced by the MRI-driven turbulence naturally force the final distribution of the angular momentum to be nearly Keplerian (see Figure 17), independent of the initial matter distribution.

Figure 17: Specific angular momentum as a function of radius at (thin line) and at orbits (thick line). The individual plots are labeled by model. In each case the Keplerian distribution for a test particle, , is shown as a dashed line. Image reproduced by permission from [72], copyright by AAS.

Although the magnetic field plays a crucial role through the action of the MRI, its amplitude saturates at a relatively weak level, always remaining subthermal in the disk (see Figure 18). On the other hand, simulations show that magnetic fields do tend to dominate regions outside the disk, including in the hot, magnetized corona that sandwiches the disk and in the evacuated, highly magnetized funnel region. The funnel region and the boundary between the funnel and the corona are where jets and outflows are generically observed in these simulations (Section 11.7).

Figure 18: Color contours of the ratio of azimuthally averaged magnetic to gas pressure, . The scale is logarithmic and is the same for all panels; the color maps saturate in the axial funnel. The body of the accretion disk is identified with overlaid density contours at , , , and of . The individual plots are labeled by model. In all cases, the magnetic pressure is low () in the disk, comparable to gas pressure () in the corona above and below the disk, and high () in the funnel region. Image reproduced by permission from [72], copyright by AAS.

Despite the increased complexity of the disk when MRI-driven turbulence is at play, many features and properties remain strikingly similar to what is revealed in hydrodynamic-only studies. For instance, a commonly seen feature in global, non-radiative MHD simulations of black hole accretion disks is an “inner torus” [120]. The inner torus is predominantly a hydrodynamic structure, as it is largely gas pressure supported. Generically, this inner torus consists of a sequence of equidensity and equipressure surfaces with a pressure maximum at , a cusp at , and a roughly parabolic, evacuated, though magnetized, funnel along the rotation axis. This structure is remarkably similar to the Polish doughnut model, as shown in Figure 19. Thus, many hydrodynamic torus results may retain relevance even in MRI turbulent disks [253].

Figure 19: On the left, equidensity contours calculated from an analytic Polish doughnut. On the right, equidensity contours from a numerical MHD simulation (model 90h from [99]). Note, though, that the contours on the left are linearly spaced, while those on the right are logarithmically spaced. Thus, the gradients represented on the left are shallower than those on the right. Image reproduced by permission from [253], copyright by ESO.

11.4 Novikov–Thorne (thin) disks in simulations

Early numerical simulations of black hole accretion disks nearly all focused on thick disks. This was mostly for computational convenience as thicker disks require fewer resources than thin ones. In recent years, though, the resources have become available to start testing thinner disks. This has enabled researchers to begin rigorously testing the Novikov–Thorne model (Section 5.3), especially the assumption that the internal torque of the disk vanishes at the ISCO. This is important as some researchers have argued that magnetic fields might nullify this hypothesis by maintaining stresses inside the ISCO [164, 104, 25].

Initial global MHD simulations did show some variations from the Novikov–Thorne model at the ISCO at the level of  10% [123, 226], while others [240] showed much smaller deviations  3%. Penna and collaborators [240] discussed possible reasons for the discrepancies, and suggested that the results are, in fact, in better agreement than they may have initially appeared. Their conclusion is that, in fact, the Novikov–Thorne model matches numerical results quite well (see Figure 20, plus animations can be viewed at [239]). Any remaining differences are small enough to have negligible effect on common applications of the Novikov–Thorne model, such as measuring black hole spin via spectral fitting [166, 327] (Section 12.1).

Figure 20: Left: Time-averaged rest mass density in the – plane for four GRMHD simulations with and various disk thicknesses. The dashed vertical line marks the ISCO. The disk opening angle, , and effective Shakura–Sunyaev viscosity, , are reported in each panel. The top three panels have and the inner edge of the disk is located outside the ISCO. The bottom panel has and the density increases monotonically down to the event horizon. Figure from [241]. Right: Various fluxes as functions of radius for a numerical Novikov–Thorne disk simulation. Top: Mass accretion rate. Second panel: Accreted specific angular momentum. Solid line is simulation data; dashed line gives Novikov–Thorne solution; dotted line is ISCO value. Note that the specific angular momentum does not drop significantly inside the ISCO, unlike for thick disks, such as in Figure 17. Third panel: The “nominal” efficiency, which is the total loss of specific energy from the fluid. Bottom panel: Specific magnetic flux. The near constancy of this quantity inside the ISCO is an indication that magnetic stresses are not significant in this region. Image reproduced by permission from [240].

We mention that since present day relativistic numerical simulations do not treat the radiation in geometrically thin, optically thick disks, the simulations of Novikov–Thorne disks have so far employed an ad hoc cooling prescription [278, 226]. This prescription conveniently makes the same assumption that the Novikov–Thorne model does: that all energy dissipated as heat in the disk is radiated away locally on roughly the orbital timescale. This is probably reasonable for an appropriate range of mass accretion rates, though it will be good to test this assumption with future global radiation-MHD simulations.

11.5 ADAFs in simulations

A lot of recent simulation work has been focused on exploring more ADAF-like flows under the action of realistic MRI turbulence [222, 322, 320]. To some extent, this is an extension of the Polish-doughnut simulations of the last decade, yet goes beyond it in at least two important respects: 1) the simulations cover a significantly larger spatial range (a few hundreds of versus a few tens); and 2) the simulations explore much longer temporal evolution (hundreds of thousands of versus tens of thousands). The result of this is that the simulations are able to explore steady-state accretion out to much larger radii ( as opposed to ). More can be expected from this work in the coming years.

11.6 Oscillations in simulations

Hydrodynamic simulations have also been used extensively to study the natural oscillation modes of relativistic disks, particularly as they might relate to QPOs (Section 12.4). Rezzolla and collaborators [263, 264] identified a mode they referred to as a -mode that occurred in a near 3:2 ratio with the radial epicyclic mode, which was subsequently confirmed through numerical simulations [325, 324] (animations can be viewed at [262]). The 3:2 ratio of these modes is important as the highest-frequency QPOs in black hole low-mass X-ray binaries are observed to occur in this ratio. Later, Blaes and collaborators [45] identified a different pair of modes, the vertical epicyclic and axisymmetric breathing mode, that have a near 3:2 ratio for a broader range of parameters. Similar to [325, 324], they also demonstrated that these modes could be identified in numerical simulations.

A number of MHD simulations have also been performed which claim to observe QPOs. Kato [145] performed a global MRI simulation in a pseudo-Newtonian potential and claimed to see QPOs in the power spectrum in a measure of the luminosity associated with radial infall. Machida and Matsumoto [179] have reported the formation of a one-armed spiral within the inner torus of a global MRI simulation, and suggested that this could be responsible for the low frequency QPO in the hard state. Schnittman, Krolik, and Hawley [274] have found tentative evidence for high frequency QPO features in light curves generated by coupling a general relativistic ray tracing code to a global GRMHD simulation with an assumed emission model. Henisey and collaborators [127] found similar tentative evidence in a sample of tilted disk simulations, possibly confirming earlier suggestions that disk tilt may be an important mechanism for driving high frequency QPOs in black hole accretion disks [143, 95]. Most recently Dolence and collaborators [78] report evidence for near-infrared and X-ray QPOs in numerical simulations of Sgr A*.

Aside from these few interesting examples, though, no robust QPO signal has generically emerged in MHD simulations [22, 261]. It is unclear what the implications of this are. It may be that some missing physics, such as radiation transport, plays a fundamental role in exciting QPOs. This is an important open problem in numerical simulations of black hole accretion disks.

11.7 Jets in simulations

In recent years several numerical simulations have demonstrated the generation of jets self-consistently from simulations of disks. A few researchers have claimed to produce jetted outflows from purely hydrodynamic interactions. Nobuta and Hanawa [228], for instance, were able to achieve jetted outflows driven by shock waves when infalling gas with high specific angular momentum collided with the centrifugal barrier. However, such simulations usually require very special starting conditions and the jets tend to be transient features of the flow; therefore, it is unlikely these are related to the well-collimated expansive jets observed in many black hole systems.¹³

Among MHD simulations, some distinction should be drawn between those that impose large scale magnetic fields that extend beyond the domain of the simulation and those that, at least initially, impose a self-contained magnetic field. In the first case, the disk is often treated merely as a boundary condition for the evolution of large-scale magnetic field. Shibata and Uchida [284, 285] used this technique to show that jets could be driven both by a pinching of the fields due to radial inflow through the disk and by the force caused when twisted fields unwind.

Since 1999, many MHD disk simulations, usually starting with locally confined, weak magnetic fields, have shown some sort of magnetically driven bipolar outflow [153, 146, 73, 124, 191, 189]. Since these jets are propagating at the Alfvén speed in a magnetically dominated region (), they proceed supersonically relative to the gas. However, some authors claim the outflow is only a transient state [146] and most do not achieve the very high Lorentz factors observed in astrophysical jets [73]. It is also unclear whether or not the jets in these simulations are connected with the Blandford–Znajek process [73] (Section 10), although see [188]. McKinney [189] and Komissarov [155] suggest that part of the problem is that the jets arising from MHD simulations of disks are often contaminated by the application of numerical floors on density and energy and through numerical diffusion, both of which can lead to anomalous mass loading of the jet. In their own carefully constructed simulations, they independently found that the magnetic structures exhibited by their simulated jets are, in fact, consistent with the expectations of the Blandford–Znajek model, as illustrated in Figure 21 (animations can be viewed at [187]).

Figure 21: On the left is a schematic diagram of the Blandford–Znajek mechanism [49] for an assumed parabolic field distribution. On the right is the result of a numerical simulation from [158] showing a very similar structure. Images reproduced by permission; copyright RAS.

Tchekhovskoy and collaborators [303] have used GRMHD simulations of black hole accretion disks plus jets to investigate possible explanations for the observed radio loud/quiet dichotomy in AGN.¹⁴ For a black hole surrounded by a thin disk, the Blandford–Znajek mechanism predicts the luminosity of the jet should scale as . This limits the range of power expected for realistic AGN spins to a factor of a few tens at most – too small to explain the observed differences. Thicker disks, on the other hand, produce jets whose luminosity can scale as or even  [303], providing sufficient range to perhaps explain the dichotomy.

There has also been some work recently trying to understand a different jet dichotomy – one that is observed in black hole low-mass X-ray binaries (LMXBs). These exhibit radio emission (associated with jets) whose properties change with the observed X-ray spectral and, to a less well determined extent, temporal properties of the accretion disk [92, 94]. Briefly, compact, steady jets are observed in the Low/Hard state, whereas jets appear absent in the High/Soft state. Fragile and collaborators [102] used numerical simulations to rule out disk scale height as the controlling factor, suggesting instead that perhaps the jets are intimately connected with the corona or failed MHD wind. Alternatively, it could be that magnetic field topology is the key factor [136, 190].

11.8 Highly magnetized accretion in simulations

Recently, work has begun to focus on highly magnetized disk configurations, for which some modes of the MRI may be suppressed. One motivation is that these might provide an alternate explanation for the Low/Hard state [35] (Section 12.3).

One way a highly magnetized state could come about is as a result of a thermal instability in an initially hot, thick, weakly-magnetized disk [179, 101]. Machida and collaborators [179] simulated this process for an optically thin disk, assuming bremsstrahlung cooling and pseudo-Newtonian gravity. They found that, indeed, the densest inner regions of the disk collapse down to a cool, thin, magnetically-supported structure. Fragile and Meier [101] extended these results by including more cooling processes and using a general relativistic MHD code.

Another way to achieve a highly magnetized state is to have the disk “drag” the magnetic field in from some distant region [42, 291]. If the field has a consistent net flux, then doing so will necessarily increase the strength of the field near the hole due simply to the smaller area through which the flux must thread. Recent numerical simulations have demonstrated this effect, first in the pseudo-Newtonian limit [135] (Figure 22), and more recently in relativistic simulations [304, 192]. These examples are interesting because they have led directly to a phenomenon known as a “magnetically arrested” accretion state [218], where the accumulation of magnetic field near the black hole is sufficient to temporarily halt the inflow of matter. This state has been shown to produce jet efficiencies  [304], where

is the mass accretion rate and is the energy flux, both taken at the black hole event horizon , and the angle brackets indicate a time-averaged quantity. This can only happen if at least some of the energy powering the jet is being extracted from the rotational energy of the black hole itself!

Figure 22: Top: Distributions of density in the meridional plane at different simulation times, showing a magnetically arrested state (left) and a non-arrested state (right). Bottom: Snapshot of magnetic field lines at the same simulation times. Image reproduced by permission from [135], copyright by AAS.