Efficient Generation of Jets from Magnetically Arrested Accretion on a Rapidly Spinning Black Hole
We describe global, 3D, time-dependent, non-radiative, general-relativistic, magnetohydrodynamic simulations of accreting black holes (BHs). The simulations are designed to transport a large amount of magnetic flux to the center, more than the accreting gas can force into the BH. The excess magnetic flux remains outside the BH, impedes accretion, and leads to a magnetically arrested disc. We find powerful outflows. For a BH with spin parameter , the efficiency with which the accretion system generates outflowing energy in jets and winds is %. For , we find %, which means that more energy flows out of the BH than flows in. The only way this can happen is by extracting spin energy from the BH. Thus the simulation represents an unambiguous demonstration, within an astrophysically plausible scenario, of the extraction of net energy from a spinning BH via the Penrose-Blandford-Znajek mechanism. We suggest that magnetically arrested accretion might explain observations of active galactic nuclei with apparent .
keywords:black hole physics — (magnetohydrodynamics) MHD — accretion, accretion discs — galaxies: jets — gamma-rays: bursts — methods: numerical
Relativistic jets are a common feature of accreting black holes (BHs). They are found in both stellar-mass BHs (Remillard & McClintock, 2006) and supermassive BHs in active galactic nuclei (AGN, Tremaine et al. 2002). Jets can be very powerful, with their energy output sometimes exceeding the Eddington limit of the BH. This suggests an efficient mechanism for their production.
In seminal work, Penrose (1969) showed that a spinning BH has free energy that is, in principle, available to be tapped. This has led to the popular idea that the energy source behind relativistic jets is the rotational energy of the accreting BH. Blandford & Znajek (1977, hereafter BZ77) came up with an astrophysical scenario in which this could be achieved. In their picture, magnetic field lines are kept confined around the BH by an accretion disc. The rotation of space-time near the BH twists these lines into helical magnetic springs which expand under their own pressure and accelerate any attached plasma. In the process, energy is extracted from the spinning BH and is transported out along the magnetic field, making a relativistic jet. The BZ mechanism is a promising idea since magnetic fields are common in astrophysical accretion discs and so the requirements for this mechanism are easily met.
where is a numerical constant whose value depends on the magnetic field geometry (it is 0.053 for a split monopole geometry and 0.044 for a parabolic geometry), is the angular frequency of the BH horizon, is the magnetic flux threading one hemisphere of the BH horizon (the integral is over all , at the BH horizon, and the factor of converts it to one hemisphere), is an area element in the plane, and is the determinant of the metric. Here is the dimensionless BH spin parameter (sometimes also called ), is the radius of the horizon, is the gravitational radius of the BH, and is the BH mass. A simpler version of equation (1) with was originally derived by BZ77 in the limit . Tchekhovskoy et al. (2010) showed that the modified form written here, with , is accurate even for large spins up to , while for yet larger spins, they gave a more accurate 6th order approximation, .
Using equation (1), let us define the efficiency with which the BH generates jet power, , as the ratio of the time-average electromagnetic power that flows out of the BH, , to the time-average rate at which rest-mass energy flows into the BH, ,
where is the dimensionless magnetic flux threading the BH and is a time-average. Thus the efficiency with which a spinning BH can generate jet power depends on BH spin via the angular frequency and on the dimensionless magnetic flux . The strength of is very uncertain.
It is generally agreed that is non-zero, since magnetic flux is transported to the accreting BH by turbulent accretion. However, the key elements of this process are not agreed upon (Lubow et al., 1994; Spruit & Uzdensky, 2005; Rothstein & Lovelace, 2008; Beckwith et al., 2009; Cao, 2011). This leads to a large uncertainty in the value of . Recent time-dependent general relativistic magnetohydrodynamic (GRMHD) numerical simulations have found a rather low efficiency, , even when the central BH is nearly maximally spinning (McKinney, 2005; De Villiers et al., 2005; Hawley & Krolik, 2006; Barkov & Baushev, 2011). With such a modest efficiency it is not clear that we are seeing energy extraction from the BH. The jet power could easily come from the accretion disc (see Ghosh & Abramowicz, 1997; Livio et al., 1999).
Observationally, there are indications that some AGN in the universe may have extremely efficient jets that require (Rawlings & Saunders, 1991; Ghisellini et al., 2010; Fernandes et al., 2011; McNamara et al., 2011; Punsly, 2011). A non-spinning BH usually has , and might under special circumstances have tens of percent (e.g., Narayan et al. 2003). However, a non-spinning BH can never give , since this requires the system to produce more energy than the entire rest mass energy supplied by accretion. Values of are possible only by extracting energy from the spin of the BH. Thus, taken at face value, any robust observation of in an AGN implies that the Penrose/BZ process must be operating. This raises the following important question: Is it possible to show via a numerical simulation, an astrophysically plausible BH accretion scenario that gives jet efficiency ? To our knowledge, this has not been demonstrated with a GRMHD simulation.
Here we describe numerical simulations in which we arrange our setup such that the accreting BH receives as much large-scale magnetic flux as can be pushed into the BH by accretion. In nonradiative MHD, the limiting flux is proportional to . In fact, we supply more flux than this, so some of the flux remains outside the BH where it impedes the accreting gas, leading to a “magnetically arrested disc” (MAD, Narayan et al. 2003, see also Bisnovatyi-Kogan & Ruzmaikin 1974, 1976; Igumenshchev et al. 2003). The goal of the present simulations is to maximize and to make the jet efficiency as large as possible. As we show below, we do obtain larger efficiencies than reported in previous numerical experiments. Most interestingly, we find efficiencies greater than for a rapidly spinning BH. These experiments are the first demonstration of net energy extraction from spinning BHs via the Penrose/BZ mechanism in an astrophysically plausible setting. In §2, we discuss our numerical method, the physics of MAD accretion and our problem setup, and in §3 we discuss the results and conclude.
2 Numerical Method and Problem Setup
We have carried out time-dependent simulations of BH accretion for two values of BH spin (see Table LABEL:tab1). We use the GRMHD code HARM (Gammie et al., 2003; McKinney & Gammie, 2004) with recent improvements (McKinney, 2006; Tchekhovskoy et al., 2007, 2009; McKinney & Blandford, 2009). The numerical method conserves mass, angular momentum and energy to machine precision. We neglect radiative losses, so the simulations correspond to an ADAF mode of accretion (Narayan & McClintock, 2008).
The simulations are carried out in spherical polar coordinates modified to concentrate resolution in the collimating polar jet and in the equatorial disc. We use Kerr-Schild horizon-penetrating coordinates. We place the inner radial boundary inside the BH outer horizon (but outside the inner horizon), which ensures that no signals can propagate through the horizon to the BH exterior. The outer radial boundary is at , which exceeds the light travel distance for the duration of the simulation. Thus both boundaries are causally disconnected. We use a logarithmically spaced radial grid, for (see Table LABEL:tab1 for values of ), where we choose so that there are grid cells between the inner radial boundary and the BH horizon. For , the radial grid becomes progressively sparser, , with a smooth transition at . At the poles, we use the standard reflecting boundary conditions, while in the azimuthal direction, we use periodic boundary conditions. In order to prevent the extent of cells near the poles from limiting the time step, we smoothly deform the grid a few cells away from the pole so as to make it almost cylindrical near the BH horizon; this speeds up the simulations by a factor .
Numerical MHD schemes cannot handle vacuum. Therefore, whenever the fluid-frame rest-mass energy density, , falls below a density floor , where is the magnetic pressure in the fluid frame, or when the internal energy density, , falls below , we add mass or internal energy in the frame of a local zero angular momentum observer so as to make or (McKinney & Blandford, 2009). The factor sets the maximum possible Lorentz factor of the jet outflow. To investigate the effect of this factor on jet efficiency, we have tried two values, , . There is little difference in the results. In any case, we track the amount of mass and internal energy added in each cell during the course of the simulation and we eliminate this contribution when calculating mass and energy fluxes.
Model A0.99f (Table LABEL:tab1) uses a resolution of along -, -, and -, respectively, and a full azimuthal wedge, . This setup results in a cell aspect ratio in the equatorial region, . To check convergence with numerical resolution, at , well after the model reached steady-state, we dynamically increased the number of cells in the azimuthal direction by a factor of . We refer to this higher-resolution simulation as model A0.99fh and to A0.99f and A0.99fh combined as model A0.99fc. We also ran model A0.99 with a smaller azimuthal wedge, . We find that the time-averaged jet efficiencies of the four A0.99xx models agree to within statistical measurement uncertainty (Table LABEL:tab1), indicating that our results are converged with respect to azimuthal resolution and wedge size.
Our fiducial model A0.99fc starts with a rapidly spinning BH () at the center of an equilibrium hydrodynamic torus (Chakrabarti, 1985; De Villiers & Hawley, 2003). The inner edge of the torus is at and the pressure maximum is at (see Fig. 1a). At the initial torus has an aspect ratio and fluid frame density (in arbitrary units). The torus is seeded with a weak large-scale poloidal magnetic field (plasma ). This configuration is unstable to the magnetorotational instability (MRI, Balbus & Hawley, 1991) which drives MHD turbulence and causes gas to accrete. The torus serves as a reservoir of mass and magnetic field for the accretion flow.
Equation (1) shows that the BZ power is directly proportional to the square of the magnetic flux at the BH horizon, which is determined by the large-scale poloidal magnetic flux supplied to the BH by the accretion flow. The latter depends on the initial field configuration in the torus. Usually, the initial field is chosen to follow isodensity contours of the torus, e.g., the magnetic flux function is taken as , where the constant factor is tuned to achieve the desired minimum value of in the torus, e.g., . The resulting poloidal magnetic field loop is centered at and contains a relatively small amount of magnetic flux. If we wish to have an efficient jet, we need a torus with more magnetic flux, so that some of the flux remains outside the BH and leads to a MAD state of accretion (Igumenshchev et al., 2003; Narayan et al., 2003). We achieve this in several steps. We consider a magnetic flux function, , and normalize the magnitude of the magnetic field at each point independently such that we have everywhere in the torus. Using this field, we take the initial magnetic flux function as and tune such that . This gives a poloidal field loop centered at . The loop has a much larger spatial size and more magnetic flux than the usual initial field loop configuration considered in other studies. We maintain a nearly uniform radial distribution of for which lets us resolve the fastest growing MRI wavelength with more than cells over a wide range of radii.
Note that the above technique of starting with a large amount of poloidal magnetic flux in the torus is just a convenient trick to achieve a MAD state of accretion within the short time available in a numerical simulation. In our simulations magnetic flux is rapidly advected to the center from a distance , whereas in nature we expect magnetic flux to be advected from a distant external medium at and to grow on the corresponding accretion time. The latter time is much too long to be currently simulated on a computer, hence the need to speed up the process in the simulations. Note that field advection is likely to be more efficient in a thick, ADAF-like disc rather than in a thin disc (e.g., Lubow et al., 1994; Spruit & Uzdensky, 2005; Rothstein & Lovelace, 2008; Cao, 2011), hence we consider the present simulations to be a reasonable proxy for real ADAFs in nature. In any case, the key point is that, regardless of how a given system achieves a MAD state of accretion — whether it is by slow advection of field from large distances in a real system or through rapid advection of magnetic flux from short distances in our simulations — once the system has reached this state we expect its properties to be largely insensitive to its prior history. We thus believe the results obtained here are relevant to astrophysical objects with MADs.
Figure 1 shows results from the fiducial simulation A0.99fc for a rapidly spinning BH with . The top two panels (a) on the left show the initial torus, with purely poloidal magnetic field. The succeeding panels show how the accretion flow evolves. With increasing time, the MRI leads to MHD turbulence in the torus which causes the magnetized gas to accrete on the BH. In the process, magnetic flux is brought to the center and accumulates around the BH in an ordered bipolar configuration. These field lines are twisted into a helical shape as a result of space-time dragging by the spinning BH and they carry away energy along twin jets.
The rate of accretion of rest mass, , and rest mass energy, , at radius are given by
where is the radial contravariant component of the -velocity, and the integral is over all , at fixed . The negative sign means that the flux is defined to be positive when rest mass flows into the BH. Figure 1(e) shows the rest mass energy flux into the BH, , as a function of time (the flux has been corrected for density floors, see §2). Until a time , the MRI is slowly building up inside the torus and there is no significant accretion. After this time, steadily grows until it saturates at . Beyond this time, the accretion rate remains more or less steady at approximately code units until the end of the simulation at . The fluctuations seen in are characteristic of turbulent accretion via the MRI.
Figure 1(f) shows the time evolution of the dimensionless magnetic flux at the BH horizon. Since the accreting gas continuously brings in new flux, continues to grow even after saturates. However, there is a limit to how much flux the accretion disc can push into the BH. Hence, at , the flux on the BH saturates and after that remains roughly constant at a value around . The corresponding dimensionless magnetization parameter (Gammie et al. 1999; see Penna et al. 2010 for definition) is (much greater than ), indicating that the flow near the BH is highly magnetized. Panel (b) shows that magnetic fields near the BH are so strong that they compress the inner accretion disc vertically and decrease its thickness. The accreting gas, of course, continues to bring even more flux, but this additional flux remains outside the BH. Panels (c) and (d) show what happens to the excess flux. Even as the gas drags the magnetic field in, field bundles erupt outward (Igumenshchev, 2008), leaving the time-average flux on the BH constant. For instance, two flux bundles are seen at in Figure 1(c) which originate in earlier eruption events. Other bundles are similarly seen in Figure 1(d). During each eruption, the mass accretion rate is partially suppressed, causing a dip in (Fig. 1e); there is also a corresponding temporary dip in (Fig. 1f). Note that, unlike in 2D (axisymmetric) simulations (e.g., Proga & Begelman, 2003), there is never a complete halt to the accretion (Igumenshchev et al., 2003) and even during flux eruptions, accretion proceeds via spiral-like structures, as seen in Figure 1(d).
In analogy with , let us define the rate of inward flow of total energy (as measured at infinity) as follows,
where is the stress-energy tensor.
Figure 2 shows plots of and vs for the two simulations, A0.5 and A0.99fc. The fluxes have been averaged over time intervals and , respectively, to reduce the effect of fluctuations due to flux eruptions. The time intervals have been chosen to represent quasi-steady magnetically-arrested accretion. The calculated fluxes are very nearly constant out to , indicating that both simulations have achieved steady state in their inner regions. Consider first the results for model A0.5 with (Fig. 2a). We find units and units. The difference between these two fluxes represents the energy returned by the accretion flow to the external universe. Our simulations are non-radiative, with no energy lost via radiation. Hence the energy outflow is entirely in the form of jets and winds.
We define energy outflow efficiency as the energy return rate to infinity divided by the time-average rest mass accretion rate:
For model A0.5, the efficiency we obtain, , is much larger than the maximum efficiencies seen in earlier simulations for this spin. The key difference is that, in our simulation, we maximized the magnetic flux around the BH. This enables the system to produce a substantially more efficient outflow.
In the more extreme model A0.99f with (Fig. 2b), we find units and units. The net energy flux in this simulation is out of the BH, not into the BH, i.e., the outward energy flux via the Penrose/BZ mechanism overwhelms the entire mass energy flux flowing into the BH. Correspondingly, the efficiency is greater than : . Since the system steadily transports net total energy out to infinity, the gravitational mass of the BH decreases with time. Where does the energy come from? Not from the irreducible mass of the BH, which cannot decrease in classical GR. The energy comes from the free energy associated with the spin of the BH. The BZ effect, which has efficiency (eq. 2 with from Fig. 1f and ), accounts for most of the extracted energy.
Since greater than % efficiency has been a long-sought goal, we ran model A0.99fc for an unusually long time (). There is no indication that the large efficiency is a temporary fluctuation (see Fig. 1g). As a further check, we calculated efficiencies for each of the runs, A0.99, A0.99f, A0.99fh (Table LABEL:tab1), to estimate the uncertainty in . We conclude that and that an outflow efficiency % is achievable with a fairly realistic accretion scenario. We note, however, that by changing the initial setup, e.g., the geometry of the initial torus and the topology of the magnetic field, it might be possible to obtain even larger values of . This is an area for future investigation.
Our outflows are in the form of twin collimated relativistic jets along the poles and less-collimated sub-relativistic winds (Lovelace, 1976; Blandford & Payne, 1982). The former are mostly confined to streamlines that connect to the BH, while the latter emerge mostly from the inner regions of the accretion flow. The bulk of the outflow power is in the relativistic component. The energy outflow efficiency shows considerable fluctuations with time (Fig. 1g), reaching values as large as % for prolonged periods of time, with a long-term average value, . This may explain sources with very efficient jets (McNamara et al., 2011; Fernandes et al., 2011; Punsly, 2011). The quasi-periodic nature of the fluctuations in suggests magnetically-arrested accretion as a possible mechanism to produce low-frequency QPOs in accreting stellar-mass BHs (Remillard & McClintock, 2006) and variability in AGN (Ghisellini et al., 2010) and GRB outflows (Proga & Zhang, 2006). Additional studies are necessary to ensure the convergence of variability properties with numerical resolution.
We conclude that rapidly spinning BHs embedded in magnetically-arrested accretion flows can produce efficient outflows with . Such flows could be relevant for understanding astrophysical systems with extremely efficient jets. The fiducial model A0.99fc presented here, which is designed to mimic magnetically arrested systems in nature, has a net energy flux away from the BH and demonstrates that net extraction of energy out of an accreting BH is viable via the Penrose/BZ effect.
We thank S. Balbus, K. Beckwith, A. Benson, V. Beskin, I. Contopoulos, L. Foschini, D. Giannios, J. Goodman, I. Igumenshchev, B. Metzger, R. Penna, R. Rafikov, A. Spitkovsky, J. Stone, F. Tombesi, D. Uzdensky for discussions. We thank the anonymous referee for useful suggestions. AT was supported by a Princeton Center for Theoretical Science Fellowship. AT and RN were supported in part by NSF grant AST-1041590 and NASA grant NNX11AE16G. We acknowledge support by the NSF through TeraGrid resources provided by NICS Kraken, where simulations were carried out, NICS Nautilus, where data were analyzed, and NCSA MSS, where data were backed up, under grant numbers TG-AST100040 (AT), TG-AST080026N (RN) and TG-AST080025N (JCM).
- Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
- Barkov & Baushev (2011) Barkov M. V., Baushev A. N., 2011, New Ast., 16, 46
- Beckwith et al. (2009) Beckwith K., Hawley J. F., Krolik J. H., 2009, ApJ, 707, 428
- Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan G. S., Ruzmaikin A. A., 1974, Ap&SS, 28, 45
- Bisnovatyi-Kogan & Ruzmaikin (1976) Bisnovatyi-Kogan G. S., Ruzmaikin A. A., 1976, Ap&SS, 42, 401
- Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
- Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
- Cao (2011) Cao X., 2011, ApJ, 737, 94
- Chakrabarti (1985) Chakrabarti S. K., 1985, ApJ, 288, 1
- De Villiers & Hawley (2003) De Villiers J.-P., Hawley J. F., 2003, ApJ, 589, 458
- De Villiers et al. (2005) De Villiers J.-P., Hawley J. F., Krolik J. H., Hirose S., 2005, ApJ, 620, 878
- Fernandes et al. (2011) Fernandes C. A. C., et al., 2011, MNRAS, 411, 1909
- Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444
- Gammie et al. (1999) Gammie C. F., Narayan R., Blandford R., 1999, ApJ, 516, 177
- Ghisellini et al. (2010) Ghisellini G., et al., 2010, MNRAS, 402, 497
- Ghosh & Abramowicz (1997) Ghosh P., Abramowicz M. A., 1997, MNRAS, 292, 887
- Hawley & Krolik (2006) Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103
- Igumenshchev et al. (2003) Igumenshchev I., Narayan R., Abramowicz M., 2003, ApJ, 592, 1042
- Igumenshchev (2008) Igumenshchev I. V., 2008, ApJ, 677, 317
- Livio et al. (1999) Livio M., Ogilvie G. I., Pringle J. E., 1999, ApJ, 512, 100
- Lovelace (1976) Lovelace R. V. E., 1976, Nature, 262, 649
- Lubow et al. (1994) Lubow S. H., Papaloizou J. C. B., Pringle J. E., 1994, MNRAS, 267, 235
- McKinney (2005) McKinney J. C., 2005, ApJ, 630, L5
- McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
- McKinney & Blandford (2009) McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126
- McKinney & Gammie (2004) McKinney J. C., Gammie C. F., 2004, ApJ, 611, 977
- McNamara et al. (2011) McNamara B. R., Rohanizadegan M., Nulsen P. E. J., 2011, ApJ, 727, 39
- Narayan et al. (2003) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69
- Narayan & McClintock (2008) Narayan R., McClintock J. E., 2008, New Astronomy Review, 51, 733
- Penna et al. (2010) Penna R. F., et al., 2010, MNRAS, 408, 752
- Penrose (1969) Penrose R., 1969, Nuovo Cimento Rivista Serie, 1, 252
- Proga & Begelman (2003) Proga D., Begelman M. C., 2003, ApJ, 592, 767
- Proga & Zhang (2006) Proga D., Zhang B., 2006, MNRAS, 370, L61
- Punsly (2011) Punsly B., 2011, ApJ, 728, L17
- Rawlings & Saunders (1991) Rawlings S., Saunders R., 1991, Nature, 349, 138
- Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
- Rothstein & Lovelace (2008) Rothstein D. M., Lovelace R. V. E., 2008, ApJ, 677, 1221
- Spruit & Uzdensky (2005) Spruit H. C., Uzdensky D. A., 2005, ApJ, 629, 960
- Tchekhovskoy et al. (2007) Tchekhovskoy A., McKinney J. C., Narayan R., 2007, MNRAS, 379, 469
- Tchekhovskoy et al. (2009) Tchekhovskoy A., McKinney J. C., Narayan R., 2009, ApJ, 699, 1789
- Tchekhovskoy et al. (2010) Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50
- Tremaine et al. (2002) Tremaine S., et al., 2002, ApJ, 574, 740