Bardeen–Petterson alignment, jets, and magnetic truncation in GRMHD simulations of tilted thin accretion discs

ABSTRACT

Prevalent around luminous accreting black holes, thin discs are challenging to resolve in numerical simulations. When the disc and black hole angular momentum vectors are misaligned, the challenge becomes extreme, requiring adaptive meshes to follow the disc proper as it moves through the computational grid. With our new high-performance general relativistic magnetohydrodynamic (GRMHD) code H-AMR, we have simulated the thinnest accretion disc to date, of aspect ratio H/R ≈ 0.03 ≈ 1.7°, around a rapidly spinning (a ≈ 0.9) black hole, using a cooling function. Initially tilted at 10°, the disc warps inside ∼5 r_g into alignment with the black hole, where r_g is the gravitational radius. This is the first demonstration of Bardeen–Petterson alignment in MHD with viscosity self-consistently generated by magnetized turbulence. The disc develops a low-density high-viscosity (α_eff ∼ 1.0) magnetic-pressure–dominated inner region at r ≲ 25r_g that rapidly empties itself into the black hole. This inner region may, in reality, due to thermal decoupling of ions and electrons, evaporate into a radiatively inefficient accretion flow if, as we propose, the cooling time exceeds the accretion time set by the order unity effective viscosity. We furthermore find the unexpected result that even our very thin disc can sustain large-scale vertical magnetic flux on the black hole, which launches powerful relativistic jets that carry |$20\!-\!50{{\ \rm per\ cent}}$| of the accretion power along the angular momentum vector of the outer tilted disc, providing a potential explanation for the origin of jets in radio-loud quasars.

accretion, accretion discs, black hole physics, MHD, methods: numerical, galaxies: jets

1 INTRODUCTION

Black holes (BHs) in X-ray binaries (XRB) and possibly active galactic nuclei (AGN) cycle during their lifetimes through different accretion states characterized by the total luminosity and spectral hardness. It is widely believed that the accretion luminosity expressed as a fraction of the Eddington limit L_Edd is an important factor in determining the BH’s accretion state (Esin, McClintock & Narayan 1997; McClintock & Remillard 2006; Remillard & McClintock 2006). At low luminosities (L ≲ 0.01L_Edd), in the low-hard state, a promising model is the advection-dominated accretion flow (ADAF, Ichimaru 1977; Narayan & Yi 1994, 1995a,b). In an ADAF, the disc surface density is so low that the plasma can decouple into a two-temperature electron–ion plasma (Shapiro, Lightman & Eardley 1976). Since the ions are unable to cool on the accretion time, most of the dissipated energy is advected inwards in the disc or expelled in outflows leading to a low radiative efficiency. It is also known through general relativistic magnetohydrodynamic (GRMHD) simulations that these thick accretion discs can sustain and advect inwards large-scale poloidal magnetic flux (e.g. De Villiers, Hawley & Krolik 2003; McKinney 2006; Beckwith, Hawley & Krolik 2008; McKinney & Blandford 2009; Tchekhovskoy, Narayan & McKinney 2011; McKinney, Tchekhovskoy & Blandford 2012; Tchekhovskoy & McKinney 2012) which launches powerful jets when it reaches the central BH (Blandford & Znajek 1977). The large-scale height makes ADAFs numerically easy to study since they do not require high resolutions and have short viscous times. They have been studied extensively in GRMHD and are relatively well understood.

However, the observed emission from X-ray binaries in the high-soft state, and from high-Eddington fraction AGN, is incompatible with the ADAF solution. Their thermal emission spectrum requires a geometrically thin, optically thick accretion disc (Novikov & Thorne 1973; Shakura & Sunyaev 1973). Though X-ray binaries and AGN only spend a rather short amount of time in the thin disc regime, most of the BH growth and feedback may still occur there since the accretion rate is several orders of magnitude higher. Indeed, a typical bright quasar (or XRB) radiates at |${\sim }10{{\ \rm per\ cent}}$| of the Eddington rate, |$\dot{M}_{\rm Edd} = 10L_{\rm Edd}/c^2$|⁠, and its disc thickness is extremely small, |$H/R\simeq 0.01\times (\dot{M}/0.1\dot{M}_{\rm Edd})^{0.9}$| (see e.g. Fig. 3 in Piran, Sa̧dowski & Tchekhovskoy 2015). Thus, it is crucial to study these systems to understand the growth and feedback of supermassive BHs.

Modelling of relativistic iron line profiles suggests that most bright local AGN harbour rapidly spinning supermassive BHs, with spin parameter a > 0.8 (Reynolds 2014). Because the infalling material is unaware of the orientation of the BH spin, its angular momentum vector is expected to be misaligned with respect to the black hole spin vector, resulting in a tilted accretion disc. In fact, there exist several such candidates for both X-ray binaries (Hjellming & Rupen 1995; Greene, Bailyn & Orosz 2001; Maccarone 2002) and AGN (Caproni, Abraham & Mosquera Cuesta 2006; Caproni et al. 2007). Analytic theory predicts that the inner parts of a thin disc would align with the BH midplane due to the Bardeen and Petterson effect (hereafter BP, Bardeen & Petterson 1975; Papaloizou & Pringle 1983; Kumar & Pringle 1985; Pringle 1992; Ogilvie 1999) out to a Bardeen–Petterson radius, r_BP. The alignment can have profound consequences for the growth and feedback of supermassive BHs (Rees 1978; Scheuer & Feiler 1996; Natarajan & Pringle 1998) since a larger r_BP implies a larger alignment torque on the BH, forcing it into alignment with the outer disc on very short time-scales. Aligned accretion subsequently leads to rapid BH spinup by the disc, even if most of the infalling material is initially misaligned (and even counter-aligned in some cases, see King et al. 2005).

The BP effect is caused by the interplay between general relativistic Lense–Thirring precession (Lense & Thirring 1918) and viscosity. BP alignment has been observed in pioneering smoothed particle hydrodynamics (SPH) simulations (Nelson & Papaloizou 2000; Lodato & Pringle 2007; Lodato & Price 2010). However, because these simulations are non-relativistic and hydrodynamic, they are unable to include the full effects of GR and treat the anisotropy of magnetized turbulence accurately. In fact, GRMHD is a powerful way to model the non-linear nature of the magnetorotational instability (MRI) driven turbulence (Balbus & Hawley 1991) responsible for the viscosity whilst including the full effects of GR. However, GRMHD work at a moderate thickness of H/R = 0.08 did not find Bardeen–Petterson alignment (Morales Teixeira et al. 2014; Zhuravlev et al. 2014). Moreover, there is growing evidence that the interaction between the disc, magnetized corona, and the jets should be taken into account. Jets can torque the inner accretion disc into alignment (McKinney, Tchekhovskoy & Blandford 2013) before they align with the outer accretion flow (Liska et al. 2018). The corona, which we define as the hot bloated magnetic pressure supported flow surrounding the thin disc, on the other hand, is not expected to align since thick flows cannot exhibit Bardeen–Petterson alignment (Ivanov & Illarionov 1997; Papaloizou & Pringle 1983). Thus, full GRMHD simulations, which describe the entire thin-disc–corona–jet system, are uniquely positioned to address the more than 40-year-old fundamental problem whether BP alignment occurs in the inner parts of thin tilted discs.

Indications are that jets in systems that contain a thin disc are rare: only 10 per cent of quasars, or luminous AGN, are observed to produce relativistic jets and the associated radio emission (e.g. Sikora, Stawarz & Lasota 2007), while there are no convincing observations in soft-state X-ray binaries (though see Rushton et al. 2012). It is crucial to understand what factors are responsible for the formation and destruction of jets since they can be the dominant feedback mode in galaxy clusters (see e.g. Fabian 2012). Early theoretical work suggested that thin accretion discs are not expected to have jets since the poloidal magnetic flux may diffuse out before it can advect inwards (Lubow, Papaloizou & Pringle 1994). However, the non-uniform vertical structure of accretion discs and their turbulence may aid in the inward advection of poloidal magnetic flux (Rothstein & Lovelace 2008; Guilet & Ogilvie 2012, 2013).

The small vertical extent of thin discs makes them very difficult to study numerically, with the computational cost scaling as (H/R)⁻⁵ per accretion time. Because of this high cost, numerical work studying the physics of such discs has been mostly limited to shearing box simulations and semi-analytical studies. In fact, there are no 3D GRMHD simulations available for thin discs of aspect ratios H/R < 0.05, and the thinnest discs so far simulated in 3D GRMHD (0.05 < H/R < 0.1) were aligned, which enabled the vertical wavelength of the MRI to be resolved by a grid focused on the equatorial plane, leading to cells compressed in the θ-direction and elongated in the r- and ϕ-directions (e.g. Shafee et al. 2008; Noble, Krolik & Hawley 2009, 2010; Penna et al. 2010; Morales Teixeira et al. 2014; Avara, McKinney & Reynolds 2016; Morales Teixeira, Avara & McKinney 2017). Studying tilted discs, whose orbital motion does not conform to the main directions of the grid, is more difficult than aligned ones because one can assume neither axisymmetry nor use elongated cells to speed up the simulations.

In this work, we present the thinnest global GRMHD accretion disc simulations to date and study Bardeen–Petterson alignment and jet launching. We describe our setup in Section 2, present the results and discussion in Sections 3 and 4, and conclude in Section 5.

2 NUMERICAL MODELS

We use for this work our state-of-the-art GPU-accelerated GRMHD code H-AMR (Liska et al. 2018), which finds its heritage in the HARM2D code (Gammie, McKinney & Tóth 2003; Noble et al. 2006) but has been significantly expanded to take advantage of vectorization and SIMD instructions, achieving 1.2 × 10⁶ zone-cycles per second on a single Intel Skylake 3.3GHz CPU core, and include advanced features discussed below. We have developed a CUDA version of the code that reaches 10⁸ zone-cycles/s on an NVIDIA Tesla V100 GPU. It also features a staggered grid for constrained transport of magnetic fields (Gardiner & Stone 2005), uses an Harten-Lax-van Leer (HLL) Riemann solver (Harten 1983) and advanced features such as adaptive mesh refinement (AMR) and local adaptive time-stepping (LAT) that bring down the cost of the simulation described in this work by two extra orders of magnitude (in comparison to using a uniform grid and a global, fixed time-step).

Our thin disc model considers a spinning BH of spin parameter a = 0.9375. We insert an initial Fishbone & Moncrief torus (Fishbone & Moncrief 1976), which is maintained in hydrostatic equilibrium by the BH’s vertical component of gravity counterbalancing the pressure forces in the disc. This torus has an inner radius r_in = 12.5r_g, where r_g = GM_BH/c² is the gravitational radius, with the pressure maximum at r_max = 25r_g. We use an ideal gas law equation of state, p_g = (Γ − 1)u_g, with u_g the gas internal energy density, p_g the gas pressure, and Γ = 5/3 the adiabatic index. We insert a poloidal magnetic field in the torus described by a covariant vector potential A_ϕ = (ρ − 0.05)²r³, with ρ the gas density of the torus. The magnetic field is subsequently normalized by setting |$\bar{\beta }={\max p_{\rm g}}/{\max p_{\rm b}}=30$|⁠, where p_b is the magnetic pressure and both maxima are taken over the torus separately. The disc is tilted by |$\mathcal {T}_0 = 10^\circ$| with respect to the BH equator (see Liska et al. 2018 for details). Since the code is scale-free, we set the initial gas density maximum to ρ_max = 1. To maintain the desired value of disc thickness, (H/R)_target = 0.03, we cool the disc towards its target temperature on the Keplerian time-scale using a prescribed source term (Noble et al. 2009). We disable this cooling function when p_b/(ρc²) ≳ 5 at r ≳ 10r_g to avoid cooling the jets.

In this work, we, for the first time, use the full AMR capability of H-AMR, in order to ensure sufficient resolution within the thin disc. For the refinement criterion, we use a density cutoff equal to ∼4 per cent of the maximum disc density. The lower refinement levels are set such that jumps in spatial resolution are limited to a factor 2. To avoid noise from sporadic refinement and derefinement, we only derefine at a two times lower density than set as the refinement criterion. Typically, we use three levels of AMR and attain a speedup by a factor 32–60 in comparison to an equivalent uniform grid. By evolving lower AMR levels or parts of the grid further from the black hole, which have larger cell sizes, at a larger time-step, LAT gives an additional speed-up of factor ∼5 while reducing inversion errors in relativistic regions (Chatterjee et al. 2019). The effective resolution of 2880 × 864 × 1200 cells in r-, θ-, and ϕ-directions, respectively, resolves the target disc thickness, (H/R)_target = 0.03, by approximately 8 cells in all three dimensions, and the base grid of 720 × 216 × 300 cells guarantees that the jets and corona are also sufficiently resolved. As we will see in Section 3.5, the fastest growing MRI wavelength is resolved by ≳10 cells over most of the disc.

Operating in Kerr–Schild spherical polar coordinates, we place the inner r-boundary just inside the event horizon and the outer r-boundary at 10⁵r_g, so that the flow is unaffected by the boundaries. We use transmissive polar boundary conditions in the θ-direction (Liska et al. 2018), and periodic boundary conditions in the ϕ-direction.

In BH powered jets, the gas either drains off the field lines into the BH or gets flung away along the field lines into the jet. This leads to a runaway drop in density around the jet’s stagnation surface, at which the outflow velocity vanishes. To avoid the development of vacuum regions and the breakdown of ideal MHD, we replenish the density in the regions where the density drops too low. For this, we follow the approach of Ressler et al. (2017) and approximate physical processes that mass-load relativistic jets at their base by applying a density floor of ρ_floorc² = p_b/10 throughout the jet. This adds a small amount of density on the field lines and does not noticeably affect the energetics of the jets.

3 RESULTS

We start our analysis at t = 10⁴r_g/c when the disc has cooled to its equilibrium thickness, reaching a highly turbulent state, as seen in Fig. 1. We measure throughout this work vector quantities for the disc, corona, and jet in tilted spherical polar coordinates |$r,\tilde{\theta },\tilde{\phi }$| that are aligned with the disc’s angular momentum at each radius. We define the jet–corona boundary at p_b = 1.5ρc² and the corona–disc boundary at ρ = 1.2; this is approximately 10⁻² times the maximum density in the disc after the initial cooling completes.

Figure 1.

Vertical slice at t = 18, 091r_g/c through density of our tilted disc with an aspect ratio H/R = 0.03, the thinnest to date in a GRMHD simulation, shows that the inner parts of the disc (r ≲ 5r_g) align with the rapidly spinning BH’s equator (horizontal in the figure; see the inset for a zoom-in). The warp is accompanied by a ∼180° sweep in the precession angle (see Fig. 2b), creating a ‘negative’ tilt around r ∼ 8r_g. Thick black curve shows that the disc tilt |$\mathcal {T}$| smoothly decreases towards small radii before flattening at 0° for r ≲ 5r_g (see also Fig. 2a). This is the first demonstration of Bardeen & Petterson (1975) alignment in a GRMHD simulation, i.e. in GR and in the presence of large- and small-scale turbulent magnetic stresses. Magenta lines show the corona–jet boundary, p_b = 1.5ρc², and cyan lines the disc–corona boundary, ρ = 1.2.

3.1 Bardeen–Petterson alignment

Fig. 1 shows that the inner parts of the accretion disc (r ≲ 5r_g) align with the BH equator. This is the first demonstration of BP alignment in GRMHD. This alignment happens relatively soon, already ∼(500–1000)r_g/c after accretion starts, and persists throughout the simulation. The outer parts remain tilted, and the disc develops a smooth warp in between.

We can see this more quantitatively by introducing the tilt, |$\mathcal {T}$|⁠, and precession, |$\mathcal {P}$|⁠, angles. We define them as polar and azimuthal angles between the disc’s angular momentum and the BH spin vector (e.g. Fragile & Anninos 2005; Fragile et al. 2007). While this definition works well for the disc and the corona, it produces large fluctuations when applied to jets since they can undergo strong kinks and pinches, which violently change the jet angular momentum direction (by more than 10°) on very short time and length scales. For this reason, we adopt a more robust way of calculating the tilt and precession angles for the jets, as described in Appendix B of Liska et al. (2018). Namely, instead of considering the angular momentum vector orientation, we calculate the jet tilt and precession angles at each radius from the jet’s average position, (X, Y, Z), weighted by the magnetic pressure p_b:

\begin{eqnarray*} \epsilon =p_{\rm b}\times \operatorname{H}\left({p_{\rm b}}-1.5 {\rho c^2}\right), \end{eqnarray*}

(1)

where |$\operatorname{H}$| is the Heaviside step function that zeroes out the weight outside of the magnetized jet body.

Fig. 2(a) shows that the inner parts of the accretion disc, at r ≲ 5r_g, are well-aligned with the BH equator: The tilt angle |$\mathcal {T}^{\rm inner}_{\rm disc}\ll 1^\circ$| of the inner disc (r ≲ 5r_g) is much smaller than the tilt angle |$\mathcal {T}^{\rm outer}_{\rm disc}\gtrsim 5^\circ$| of the outer disc (r ≳ 10r_g). The alignment radius remains steady at r_BP ≈ 5r_g during the course of our simulation. Note that |$\mathcal {T}^{\rm outer}_{\rm disc}\simeq 5^\circ$| is smaller than the initial tilt of the disc, |$\mathcal {T}_0 = 10^\circ$|⁠. This is because the disc as a whole undergoes global alignment, which is distinctly different than the BP effect (Section 3.3).

Figure 2.

Radial profiles of quantities in our thin disc system, averaged over the time interval, 2 × 10⁴ < t < 2.25 × 10⁴r_g/c. The shaded regions represent the time variability over this interval. (a) Tilt angles of the disc, corona, and jet. Whereas the outer disc is tilted by |$\mathcal {T}\approx 5^\circ$|⁠, the inner disc at r ≲ 5r_g is aligned with the BH equator, indicating the presence of Bardeen & Petterson (1975) alignment. While the disc aligns rapidly near the BH, the corona and jet show a single radial tilt oscillation peaking at r ∼ 10r_g, which is characteristic behaviour of thick accretion discs that do not show Bardeen–Petterson alignment (e.g. Fragile et al. 2007; Liska et al. 2018). That the corona and the jets at smaller radii align less readily with the BH than the disc, suggests that the alignment is driven by the disc dynamics rather than that of the corona or jet. (b) The precession angle for the disc, corona, and jet are roughly consistent with each other at large radii, suggesting they become co-aligned. As the disc aligns at small radii, r ≲ 5r_g, the precession angles become ill-defined, and we do not show them there. (c) The effective viscosity and the sum of the Maxwell and Reynolds stresses exceed the disc’s scale height (h/r ∼ 0.03), suggesting our disc is in the viscosity-dominated warp propagation regime (see Sections 3.5 and 4.1) (d) The half-opening angle of our jets (red curve) exceeds the opening angle of the jets in a thick H/R ∼ 0.3 disc model (blue dotted curve, Liska et al. 2018), suggesting that a smaller disc thickness in this work provides less pressure support for the jet and causes its opening angle to widen.

Since the corona is relatively thick, spanning an opening angle of around 30–70° (e.g. Fig. 1), it is not expected to show BP alignment (Section 1). For such thick structures as the corona, the analytic theory predicts that BP alignment is suppressed and accretion instead occurs in a misaligned fashion through radial tilt oscillations (e.g. Ivanov & Illarionov 1997; Lubow, Ogilvie & Pringle 2002). Indeed, Fig. 2(a) shows a single radial tilt oscillation in the corona and jet peaking around r ∼ 10r_g. However, the corona still manages to align within r ≲ 2.5r_g, possibly due to torque from the disc. Since the angular momentum of the corona is negligible, the disc can easily affect its alignment.

While the inner jets are relatively closely aligned with the BH, the outer jets are torqued into misalignment by the corona. This is consistent with the previous work (Liska et al. 2018), which in the context of thick H/R ∼ 0.3 discs found that the outer disc–corona system is responsible for reorienting and collimating the jets. Similarly, the precession angle of the disc, corona, and jet are closely related at large radii, as seen in Fig. 2(b). At small radii, as the system becomes aligned, the precession angle becomes ill-defined, and we do not show it.

Since the disc shows much (∼2×) better alignment with the BH than the jet and corona, it is unlikely that the jet can be responsible for torquing the inner disc into (partial) alignment as has been demonstrated for thicker discs (e.g. McKinney et al. 2013). If this were the case, one would expect the jet tilt to be smaller than the disc tilt and, since the jet actually has to transmit its torque through the corona to the disc, one would expect better alignment for the corona as well.

3.2 Global precession and alignment

In addition to exhibiting the BP effect, which depends on space and is independent of time and in which the inner part of the disc aligns with the BH, the entire accretion system additionally undergoes global alignment that depends on time and is independent of radius. Figs 3(a) and (b) shows the time evolution of the average tilt, |$\mathcal {T}$|⁠, and precession, |$\mathcal {P}$|⁠, angles for the disc, jet, and corona. The disc and corona are angular momentum averaged, the jet is magnetic pressure averaged (Section 3.1). Similar to thick H/R ∼ 0.3 discs (Liska et al. 2018, 2019a), our thin disc aligns as a whole with the BH spin. Sorathia, Krolik & Hawley (2013) proposed that alignment (BP alignment or global alignment) may be caused by the turbulent mixing between disc annuli with different precession angles (Fig. 2b shows that the precession angle decreases as a function of radius), which leads to cancellation of misaligned angular momentum and thus produces net alignment.

Figure 3.

The tilt angle (a) and precession angle (b), as functions of time for the entire disc (blue), coronal wind within r < 250r_g (black) and jet (between 10r_g < r < 100r_g, red). The shaded regions show the standard deviations within 10r_g < r < 40r_g. The whole system precesses and over time aligns with the BH. The analytic estimate for the precession rate (green) is only valid at later times. (c) The jet half-opening angle at r = 200r_g decreases as function of time, possibly due to a fall in jet power (Figs 4a and b).

As we discussed above, this is not BP alignment for two reasons. First, BP alignment is characterized by a steady state solution that is established on a time-scale shorter than an accretion time of the disc (our disc lost |${\sim } 20 {{\ \rm per\ cent}}$| of its initial mass by the end of this simulation): otherwise, most of the disc’s mass would accrete misaligned. Second, if the disc is fed externally, such as by a larger thin disc or by fallback material in a tidal disruption event, such global alignment may disappear as the inner disc is replenished by the gas carrying the misaligned angular momentum on the accretion time of the inner disc. Indeed, Liska et al. (2018) showed that global alignment becomes slower as the disc size becomes larger due to viscous spreading.

The precession period of around 10⁵r_g/c is consistent with a Type-C quasi-periodic oscillation (QPO) frequency of |$0.2\rm {Hz}$| for a 10M_⊙ BH (e.g. Ingram et al. 2016). Although Type-C QPOs are not observed for the clean thin discs thought to be present in the soft state, we note that the precessing corona in our simulation could give rise to the QPOs typically observed in the Comptonised radiation from X-ray binaries during the transition from hard to soft state. Note that the precession angle of the corona is smaller than that of the disc, as seen in Fig. 3(b). Therefore, the precession of the corona lags that of the disc. Thus, the hard variations might lag the soft variations (see Liska et al. 2019a). Note that, due to global alignment, which causes the outer disc to align with the black hole to within 2° by 6 × 10⁴r_g/c (Fig. 3a), the precession cannot be sustained for more than a single period.

We can analytically estimate the precession period P_LT by calculating the total perpendicular angular momentum L_⊥ and Lense–Thirring torque τ_⊥ directly from the stress energy tensor |$T^{\mu }_{\nu }$| and test-particle LT precession rate Ω_LT:

\begin{eqnarray*} L_{\bot}=\int\!\!\!\int\!\!\!\int _{r_{\rm in}}^{r_{\rm max}}T^{t}_{\tilde{\phi }}\sin \mathcal {T}\, dV, \end{eqnarray*}

(2)

\begin{eqnarray*} \Omega _{\rm LT}=\frac{1}{r^{3/2}+|a|}\left(1-\sqrt{1-4\frac{|a|}{r^{3/2}}+3\frac{a^2}{r^2}}\right), \end{eqnarray*}

(3)

\begin{eqnarray*} \tau _{\bot}=\int\!\!\!\int\!\!\!\int _{r_{\rm in}}^{r_{\rm max}}\Omega _{LT} \times T^{t}_{\tilde{\phi }}\sin \mathcal {T}\, dV, \end{eqnarray*}

(4)

\begin{eqnarray*} P_{\rm LT}=2 \pi \times \frac{L_{\bot }}{\tau _{\bot }}. \end{eqnarray*}

(5)

Fig. 3(b) shows that the above analytic expression for the precession period approximately agrees with the simulation at later times if we assume r_in = 1.5 × r_isco. However, at earlier times, this overestimates the precession rate. To remedy this, r_in would need to increase in time, however it is unclear why this would happen.

3.3 Outflow and radiative efficiency

Fig. 4(a) shows that even after the disc has approached the equilibrium thickness at t ∼ 10⁴r_g/c, both the mass |$\dot{M}$| and energy |$\dot{E}$| accretion rates at the event horizon, r_H = r_g[1 + (1 − a²)^1/2] = 1.35r_g, show several peaks. They stabilize to within a factor of a few at t ≳ 4 × 10⁴r_g/c, suggesting that by then the disc has reached a quasi-steady state. That |$\dot{M} c^2\gt \dot{E}$| shows that the BH accretes more energy than it ejects out. To quantify this, we introduce dimensionless energy outflow efficiencies for the relativistic jets,

\begin{eqnarray*} \eta _{\rm jet}=-\frac{\dot{E}_{jet}}{\langle \dot{M}c^2 \rangle _{t}}, \end{eqnarray*}

(6)

where |$\dot{E}_{jet}$| is the energy accretion rate in the jets at r = 20r_g (chosen to be larger than r_H for robustness of identification of the jets), and the disc winds,

\begin{eqnarray*} \eta _{\rm wind}=\frac{\langle \dot{M}c^2 \rangle _{t}-\dot{\langle E \rangle _{t}}+\dot{E}_{jet}}{\langle \dot{M}c^2\rangle _{t}}. \end{eqnarray*}

(7)

In the above, e.g. |$\langle \dot{M} \rangle _{t}$| is the running time-average of the mass accretion rate over an interval of ±500r_g/c. Note that, from these equations, any of the rest mass energy of disc material not accreted on to the BH horizon, and not ejected in the jets, is by definition lost in a wind. Fig. 4(b) shows that the jet efficiency resides in the range |$\eta _{\rm jet} = 20\!-\!50{{\ \rm per\ cent}}$|⁠, implying that as much as |$50{{\ \rm per\ cent}}$| of the mass–energy accreted by the BH can be extracted by large-scale magnetic flux from the black hole spin energy and carried out in the form of Poynting-flux dominated jets. Fig. 5 shows that the field lines responsible for the jet launching are anchored in the BH: they extract the energy via the Blandford & Znajek (1977) effect. That the jets are so efficient is a surprising result: the standard expectation is that thin discs are incapable of holding on to large-scale poloidal magnetic flux and hence are not expected to have powerful jets (Lubow et al. 1994).

Figure 4.

Time-dependence of various quantities in our simulation. (a) Mass |$\dot{M}$| and energy |$\dot{E}$| accretion rates, which initially oscillate but stabilize later on. (b) The total outflow efficiency (black) and its constituents: jets (blue), coronal wind (red), and radiative (green) contributions. The total efficiency of |$60\!-\!80{{\ \rm per\ cent}}$| is unprecedentedly high for such a thin disc as ours and exceeds the standard thin disc efficiency, |$17.9{{\ \rm per\ cent}}$| (Novikov & Thorne 1973), by a factor of 3. Particularly surprising is the high jet efficiency, |$\eta _{\rm jet}\sim 20{{\ \rm per\ cent}}$|⁠, indicating that thin discs are capable of producing powerful relativistic jets. (c) The magnetic flux on the BH (blue) and in the disc (blue). While the disc magnetic flux, Φ_disc remains roughly constant at late times, t ≳ 4 × 10⁴r_g/c, the flux on the BH decreases, suggesting that the BH might be leaking its flux into the disc and that the jet formation by our thin disc might be a transient phenomenon reminiscent of transient jets in XRBs (see Section 4).

Figure 5.

Transverse slices of logarithm of density ρ and the logarithm of the magnitude of the proper velocity γv, at t ≈ 6 × 10⁴r_g/c. The field lines (black) that connect to the BH launch relativistic jets (demarcated by magenta lines at p_b = 1.5ρc²), while field lines threading the disc give rise to a non-relativistic coronal wind (disc–corona boundary is shown with the cyan line at ρ = 1.2, see Section 3). The twin low-density polar jets accelerate to γv ∼ 5c, whereas the winds are limited to at most γv ≲ c, or v ≲ 0.7c, which occurs near the edges of the jets.

In addition to the jets, the disc also launches a sub-relativistic disc wind, which is magnetic pressure dominated and thus similar to a corona, with energy efficiency |$\eta _{wind}\sim 20{{\ \rm per\ cent}}$|⁠. The field lines threading this outflow are anchored in the disc, suggesting it may be magneto-centrifugally driven by the Blandford & Payne (1982) mechanism. However, the precise nature of the launching mechanism is not entirely clear, since time-dependence, (magnetic) pressure gradient forces, and buoyancy forces, as well as energy extraction from the BH, may also play significant roles in launching this outflow. The outflow has a radial velocity of around 0.02 − 0.2c between 10r_g < r < 500r_g. This compares favorably to AGN ultra-fast outflows, whose velocities lie in the range of 0.1–0.4c (e.g. Tombesi et al. 2010, 2011), but is on the (very) high end of winds detected in the soft state of XRBs (e.g. Ponti et al. 2012, 2016; Miller et al. 2016).

To characterize the radial distribution of the magnetic flux, we define a 1D poloidal magnetic flux function, Φ(B, r), as the maximum of the poloidal magnetic flux at each radius of magnetic field B,

\begin{eqnarray*} \Phi (B, r)=\max _{\tilde{\theta }_{\rm m}}\int\!\!\!\int _{0}^{\tilde{\theta }_{\rm m}} B(r,\tilde{\theta },\tilde{\phi })dA_{\tilde{\theta }\tilde{\phi }}, \end{eqnarray*}

(8)

where |$dA_{\tilde{\theta }\tilde{\phi }} = \sqrt{-g}d\tilde{\theta }d\tilde{\phi }$| is the area element and g is the metric determinant. This gives us the magnetic flux on the BH, Φ_BH = 0.5Φ(|B_r|, r_H), and the the magnetic flux content of the disc, Φ_disc = max _rΦ(B_r, r). Fig. 4(c) shows that at early times Φ_BH ≃ Φ_disc, implying that most of the positively-oriented magnetic flux resides on the BH and very little in the disc. Over time, both initially decline roughly in the same proportion, but at t ≳ 4 × 10⁴r_g/c the value of Φ_disc flattens out whereas Φ_BH continues to decline. This suggests that the magnetic flux diffuses out of the BH and into the disc. The stability of Φ_BH on short time-scales suggests that the disc is not in the magnetically arrested disc regime (MAD, Igumenshchev, Narayan & Abramowicz 2003; Narayan, Igumenshchev & Abramowicz 2003). However, it is not yet known whether the magnetic flux expulsions and associated magnetic flux variations, as characteristic of thick MADs (Tchekhovskoy et al. 2011), are also present at the disc thickness considered here. Assuming that the thinner the disc the easier it is to saturate the BH with magnetic flux (Tchekhovskoy et al. 2014), it is likely that our disc is close to the saturation of the MAD state, since the dimensionless magnetic flux (⁠|$\Phi _{\rm BH}/(\dot{M}r_{g}^2c)^{1/2}\approx 25$|⁠) is only a factor of 2 smaller than the saturation value for thick discs (Tchekhovskoy et al. 2011; McKinney et al. 2012).

So far, we have considered mechanical outflows. To compute the full efficiency, we also need to account for radiative losses in the disc. A rough estimate for the radiative efficiency is obtained by integrating the cooling rate, |$-\dot{U}$|⁠, used to keep our disc thin (at the equilibrium thickness), from the photon orbit at r_photon = 1.43r_g to a large enough radius beyond which the disc luminosity is negligible, e.g. r_max = 100r_g, over volume:

\begin{eqnarray*} \eta _{\rm rad}=\frac{-\int\!\!\!\int\!\!\!\int _{r_{\rm photon}}^{{\rm r_{max}}} \dot{U}dV}{\langle \dot{M} c^2 \rangle _{t}}, \end{eqnarray*}

(9)

where |$dV = \sqrt{-g}dr d\tilde{\theta }d\tilde{\phi }$| is the volume element. This gives a radiative efficiency consistent with the Novikov & Thorne (1973) value of |${\sim } 18{{\ \rm per\ cent}}$| (Fig. 4b), which is encouraging for the continuum fitting method used to determine BH spin (see e.g. McClintock, Narayan & Steiner 2014). However, as also seen in Fig. 4(b), the total efficiency η_tot = η_jet + η_wind + η_rad of our disc reaches 60–80 per cent.

3.4 Jet geometry

A major surprise of this work is the finding of powerful jets, even in our thin disc accretion system. How do the jets from thin discs compare to the more familiar jets from thick discs? To carry out the comparison, we estimate the cross-sectional area A of the jet, defined by p_b > 1.5ρc², at each radius r:

\begin{eqnarray*} A=\int\!\!\!\int \operatorname{H}(p_b-1.5\rho c^2) dA_{\tilde{\theta }\tilde{\phi }}, \end{eqnarray*}

(10)

Assuming the jet cross-section is circular, we obtain its effective half-opening angle as

\begin{eqnarray*} \Delta \theta =\sqrt{\frac{A}{\pi r^2}}. \end{eqnarray*}

(11)

Fig. 2(d) compares the radial dependence of our jet opening angle to that of a jet from a thick H/R ∼ 0.3 disc, as found in non-radiative GRMHD simulations described in Liska et al. (2018). The jet opening angle in the present work is much wider, presumably due to lack of the support in collimation pressure from the much thinner and cooler ambient medium consisting out of the disc and corona. This wider opening angle associated with thin discs may lead to differences in the radiative flux (Fragile, Wilson & Rodriguez 2012), possibly explaining the absence of any observations of soft state jets in XRBs. Fig. 3(c) shows that Δθ decreases over time, possibly due to a decrease in jet power, as implied by the decreasing mass accretion rate and approximately constant jet efficiency, as seen in Figs 4(a) and (b). If the power of the jet is lower, its pressure is also lower, so the ambient medium would compress the jet into a smaller opening angle.

3.5 Radial structure

Fig. 6 shows a vertical slice through a late-time state of the system. The accretion disc exhibits a sharp transition in density around |x| ∼ 25r_g: at smaller radii, the accretion disc is lower-density and contains a larger magnetic flux (larger number of magnetic field lines) than at larger radii. We can see the same more quantitatively in Fig. 7, which shows disc radial profiles at late time. The disc reaches a quasi-steady state with a low-density, low plasma-β inner disc truncated at r ∼ 25r_g coupled to a high density, high plasma-β outer disc at larger radii. Here, β is calculated by taking the ratio of the density weighted average gas pressure p_g and the density weighted average magnetic pressure p_b. Note the presence of a density bump before the ISCO due to mass accumulation, as is characteristic for thin discs (Penna et al. 2010). To understand how this density and magnetic structure affects the disc dynamics, we compute the Maxwell viscosity, α_M, Reynolds viscosity, α_R, and effective viscosity, α_eff, parameters as follows,

\begin{eqnarray*} \alpha _M=\frac{\langle b_{r}b_{\tilde{\phi }}\rangle _{\rho }}{\langle p_{\rm b}+p_{\rm g}\rangle _{\rho }}, \end{eqnarray*}

(12)

\begin{eqnarray*} \alpha _R=\frac{\langle (\rho c^2+u_{\rm g}+p_{\rm g})\delta u_{r} \delta u_{\tilde{\phi }}\rangle _{\rho }}{\langle p_{\rm b}+p_{\rm g}\rangle _{\rho }} \end{eqnarray*}

(13)

\begin{eqnarray*} \alpha _{\rm eff}=\frac{\langle u_{\tilde{\phi }}u_{r}\rangle _{\rho }}{\langle c_{i}^{2}\rangle _{\rho }}, \end{eqnarray*}

(14)

\begin{eqnarray*} c_{i}=\sqrt{\frac{p_{\rm g}}{\rho c^2}}. \end{eqnarray*}

(15)

Here, c_i is the isothermal sound speed and δu = u − 〈u〉_ρ is the velocity deviation from the mean. The density-weighted average of a given quantity |$\rm x$| is denoted by |$\langle \rm x \rangle _{\rho }$|⁠. The averaging is performed within ∼10 scaleheights of the disc’s midplane, more specifically for ρ > 0.01.

Figure 6.

Late-time vertical slice through density in our GRMHD simulation of a thin disc (compare to Fig. 1 at earlier times). The disc reaches a steady state where poloidal magnetic flux (black lines) gets trapped in the inner disc. The density maximum in the disc at |x| ∼ 25r_g marks the transition from the inner low-density, high-magnetization disc to the high-density, low-magnetization outer disc (see also Fig. 7a). In Sec 4.2 we argue that the high effective viscosity α_eff ∼ 1 in the inner disc may cause radiative cooling to become inefficient and the flow to transition into an ADAF.

Figure 7.

Radial profiles averaged over 6 × 10⁴ < t < 6.25 × 10⁴r_g/c. Radii within the ISCO (r_ISCO ≈ 2.04r_g) are shaded grey. (a) The accretion disc transitions from an inner low-β low-density flow to an outer high-β high-density flow at r ≈ 25r_g. (b) The radial dependence of the viscous stresses. Since the effective viscosity exceeds the Reynolds and Maxwell stresses by one order of magnitude, α_eff ≫ α_{M, R}, it is likely that large-scale torques transport angular momentum outwards. The dotted lines denote negative values of α_eff. The disc is outflowing beyond the stagnation surface at r ≈ 23r_g. (c) Within 25r_g the density scale height (H/R)_ρ increases above the thermal scale height (H/R)_thermal due to magnetic pressure support in the inner disc. (d) We maintain more than 10 cells per MRI wavelength, Q, in all 3 dimensions, |$r,\tilde{\theta },\tilde{\phi }$| over most of the disc.

Fig. 7(b) shows that the Reynolds stress exceeds the Maxwell stress by a factor of ∼2 in the inner disc while in the outer disc the Maxwell stress dominates. The effective viscosity exceeds the sum of Reynolds and Maxwell viscosity contribution by an order of magnitude and becomes negative in the corona. This might be indicative of large-scale magnetic torques/winds removing the angular momentum from the inner disc and is not surprising since we define the effective viscosity based on a 1D stationary α-disc model that does not take into account any time-variability and magnetic effects. In fact, deviations of (GR)MHD results from the α-viscosity description are rather typical (Avara et al. 2016; Morales Teixeira et al. 2014; McKinney et al. 2012; Penna et al. 2010; Sorathia, Reynolds & Armitage 2010) and indicative of the intrinsic limitations of the α-disc (Shakura & Sunyaev 1973), especially at greater field strengths. It is an interesting question whether the reported discrepancy (King, Pringle & Livio 2007) between the values for the α-viscosity found in local MHD simulations (derived from Maxwell and Reynolds stresses) and those constrained from observations (derived from the effective viscosity set by accretion timescale) might be due to the limitations of the α-disc description rather than a mismatch between the simulations and observations.

Interestingly, while the Maxwell and Reynolds stresses remain positive throughout the disc, the effective viscosity (in other words u_r) becomes negative for r > 23r_g, causing the disc to spread out (Fig. 7b). Viscous spreading is characteristic of all finite-size accretion discs and is caused by an outwards flux of conserved angular momentum as mass moves inwards. However, in contrast to the thick H/R ∼ 0.3 discs considered in our previous work (Liska et al. 2018), where viscous spreading caused the precession to stall, it is unlikely to significantly affect the internal disc dynamics in this work since the timescale for viscous spreading is more than an order of magnitude longer than the simulation runtime.

Due to the global alignment of the disc over time, the tilt angle evolves from the initial 10° to ∼2° at late times, and therefore our simulation sweeps through a wide range of tilt angles. Throughout the evolution, we do not see any evidence of significant tilt-related effects on the internal disc dynamics, such as standing shocks/sharp entropy gradients aligned with the lines of nodes (Fragile & Blaes 2008).

We define the thermal and density disc scale heights, respectively, as

\begin{eqnarray*} (H/R)_{\rm thermal}=\frac{\langle c_i \rangle _{\rho ^2}}{\langle v_{\tilde{\phi }} \rangle _{\rho ^2}}, \end{eqnarray*}

(16)

\begin{eqnarray*} H/R)_{\rm \rho}=\langle |\tilde{\theta }-\tilde{\theta }_{avg}|\rangle _{\rho } \end{eqnarray*}

(17)

Here, c_i is the isothermal sound speed and v is 3-velocity, and |$\tilde{\theta }_{avg}$| is the average θ-position of the disc’s midplane in tilted coordinates. Fig. 7(c) shows that while the thermal scale height remains approximately constant for r > 3r_g, (H/R)_thermal ≈ 0.03, the density scale height develops a bump, (H/R)_ρ ∼ 0.08 at r ∼ 10r_g, due to the excess magnetic pressure support in a strongly magnetized, β ≪ 1, disc. Around the ISCO the scaleheight increases, because the disc becomes super Keplerian, while the cooling function assumes a Keplerian disc.

To quantify the degree to which our simulation resolves the magnetized turbulence in the disc, we compute the quality factors |$Q_{r,\tilde{\theta },\tilde{\phi }}$|⁠, which give the number of cells per MRI wavelength in each of the three directions:

\begin{eqnarray*} Q_{r,\tilde{\theta },\tilde{\phi }}=\frac{2\pi }{\Delta _{r,\tilde{\theta },\tilde{\phi }}}\frac{\langle v^{\rm A}_{r,\tilde{\theta },\tilde{\phi }} \rangle _w}{\langle \Omega \rangle _w}, \end{eqnarray*}

(18)

\begin{eqnarray*} v^{A}_{r,\tilde{\theta },\tilde{\phi }}=\sqrt{\frac{b_{r,\tilde{\theta },\tilde{\phi }}^2}{(\rho c^2+u_{\rm g}+p_{\rm g}+b^2)}}, \end{eqnarray*}

(19)

\begin{eqnarray*} \Omega =\frac{v_{\tilde{\phi }}}{r}, \end{eqnarray*}

(20)

where |$\Delta _{r,\tilde{\theta },\tilde{\phi }}$| is the size of the cell. The Alfven speed, v_A, and angular frequency, Ω, are volume-averaged with |$w=\sqrt{b^2\rho }$| as the weight since a large fraction of the mass in our thin discs resides in equatorial current sheets where the magnetic field vanishes. Fig. 7(d) shows that the whole disc is well resolved with |$Q_{r,\tilde{\theta }}\gt 30$|⁠, |$Q_{\tilde{\phi }} \gt 200$| in the inner disc and |$Q_{r,\tilde{\theta }}\gt 10$|⁠, |$Q_{\tilde{\phi }}\gt 70$| over most of the outer disc. This satisfies the numerical convergence criteria for MRI turbulence, Q_θ > 10 and Q_ϕ > 20 (see e.g. Hawley, Guan & Krolik 2011).

4 DISCUSSION

4.1 Bardeen–Petterson alignment

We found that the inner ∼5r_g of a thin disc, H/R ≃ 0.03, initially tilted by 10° relative to the central spinning BH, undergoes alignment with the BH equator (Section 3.1). This is the first demonstration of the Bardeen & Petterson (1975) effect in a GRMHD simulation, in the presence of non-local and anisotropic turbulent MHD stresses.

This confirmation of the BP effect has profound consequences for the growth and spin evolution of supermassive BHs (SMBHs), since BP alignment is a crucial ingredient that has been assumed to take place for misaligned accretion episodes (e.g. Volonteri et al. 2005; King & Pringle 2006; Fanidakis et al. 2011). Because the alignment radius acts as lever arm helping to torque the BH, the BP effect can torque the BH and align its spin vector with the outer misaligned accretion flow on much shorter timescales than otherwise (Rees 1978; Scheuer & Feiler 1996; Natarajan & Pringle 1998). This rapid reorientation of BH spin has the potential to create the right conditions for rapid BH spin-up. If BH spin reorientation occurs on a shorter time-scale than the time-scale of a single accretion episode in a chaotic accretion scenario (in which the direction of the supplied gas angular momentum randomly changes between different accretion episodes, see e.g. Volonteri et al. 2005), then supermassive BHs can be efficiently spun up (e.g. Natarajan & Pringle 1998). In the opposite case, the accretion-supplied angular momenta would tend to cancel out, and the central BHs would be on average spun down.

Due to the high cost of GRMHD simulations, it is appealing to use them to calibrate both computationally cheaper SPH simulations and analytic theory. Our simulations are in the viscosity-dominated, H/R < α regime (Papaloizou & Pringle 1983) since both our effective viscosity parameter, α_eff ≃ 1, and local Maxwell plus Reynolds stress related viscosity parameter, α_M + α_R ≲ 0.1, exceed H/R = 0.03 (see Fig. 7). However, we did not find Bardeen–Petterson alignment in any of our simulations featuring thicker discs with <H/R = 0.1 (Liska et al. 2019a). Because all of these simulations have α_{M, R} ≲ H/R < α_eff, it appears likely that values of local Maxwell and Reynolds stresses determine the transition from the wave- (H/R > α) to the viscosity- (H/R < α) dominated regime, where BP alignment is expected.

However, the alignment radius we find substantially differs from both the analytical theory and SPH simulations. Specifically, Kumar & Pringle (1985) found analytically using the corrected Bardeen & Petterson (1975) equations in Papaloizou & Pringle (1983) a radius r_BP ≈ 300 r_g, while Nelson & Papaloizou (2000) find r_BP ≈ 100 r_g for α = 0.05 and h/r = 0.03. SPH simulations of Nelson & Papaloizou (2000) find r_BP ≈ 30r_g for α = 0.1 and h/r = 0.03. Other more recent SPH work considers a smaller disc thickness h/r ≈ 0.013 and smaller viscosity α ≈ 0.03 (Lodato & Price 2010), which makes comparison to our work difficult. However, assuming a scaling relation, r_BP ∼ [(h/r)⁻²/α]^4/7 (Kumar & Pringle 1985), these results are in similar disagreement.

One contributing factor to this disagreement with SPH simulations may be the ∼2 times larger density scale height, caused by the buildup of magnetic pressure, between 5r_g < r < 20r_g in our simulations (see Section 3.5). We indeed see that r_BP (Fig. 2a) coincides with the radius where the disc becomes thicker (Fig. 7c). Naively, analytically one would only predict (Kumar & Pringle 1985) an ≈0.4 times smaller alignment radius for a similar increase in disc thickness, insufficient to account for the full extent of the discrepancy. However since H/R ≈ α_R + α_M the disc may not be in the H/R < α viscosity dominated regime required for BP alignment. Another contributing factor might be the presence of large-scale magnetic torques in the system, which can affect the alignment radius in at least two ways. First, these torques might induce coupling between the inner aligned disc and outer misaligned disc-corona-jet system moving r_BP inwards. Second, as discussed in Section 3.5, large-scale magnetic torques can remove angular momentum from the disc and increase the inflow velocity. Because BP alignment is expected to occur more rapidly when the radial inflow velocity is smaller (Bardeen & Petterson 1975), this might reduce the efficiency of BP effect (Nealon, Price & Nixon 2015).

To test if an external torque can explain the smaller-than-predicted BP alignment radius, we calculate the ratio between the external magnetic torque dragging misaligned angular momentum inwards and the LT torque. Note that since α_M + α_R ≪ α_eff (Fig. 2c) we can safely neglect internal stresses and assume that accretion is driven solely by an external magnetic torque. This viscous torque is given by (e.g. Shakura & Sunyaev 1973),

\begin{eqnarray*} T_{\rm mag}=\alpha _{\rm eff}(h/r)^2 \times v_{\rm K} \times T^r_{\tilde{\phi }} \times \sin (\mathcal {T}), \end{eqnarray*}

(21)

while the LT torque is given by,

\begin{eqnarray*} T_{\rm LT}=T^r_{\tilde{\phi }} \times \sin (\mathcal {T}) \times \Omega _{\rm LT}, \end{eqnarray*}

(22)

Assuming α_eff ∼ 1 and using equation (5) we conclude that these two torques are equal around r ∼ 15r_g. This still overestimates r_BP by a factor ∼3, but suggests that large scale torques may indeed contribute to the discrepancy between our work and SPH simulations. Note that in this very crude calculation we neglected that the LT torque acts perpendicular to the magnetic torque. Namely, T_LT ≳ T_mag does not guarantee BP alignment, since the disc may keep precessing as a rigid body without aligning (Liska et al. 2019b). For BP alignment, misaligned angular momentum also needs to mix azimuthally such that net alignment is produced (Sorathia et al. 2013). This mixing may take place on timescales (much) longer than the viscous time and thus explain this remaining factor ∼3 discrepancy.

A smaller r_BP could have a significant effect on the predictions of SMBH growth models because smaller values of r_BP lead to less rapid alignment between BH and outer disc and perhaps consequently less rapid spin-up. For the same reason, our result implies that initially misaligned X-ray binary systems will take even longer to align than previously predicted (King & Nixon 2016), indicating that there could be many misaligned X-ray binaries today, as implied by the Lense–Thirring precession QPO model of Ingram, Done & Fragile (2009). Future work should study the effect of BP alignment and disc warp on the measured values of BH spin (e.g. McClintock et al. 2014).

The inclusion of a gas pressure dominated equation of state with adiabatic index Γ = 5/3 as in this work is only applicable in the outer accretion disc of X-ray binaries (see e.g. Zhu & Narayan 2013), while the inner part could be radiation pressure dominated with Γ = 4/3. We find that changing the adiabatic index to Γ = 4/3 for a thin disc tilted by 45° behaves qualitatively similar to Γ = 5/3 (Liska et al. 2019b). This is not unexpected since the most dominant effect of a softer equation of state is that the disc becomes thinner for a given specific internal energy u_g/(ρc²). However, the cooling function (Noble et al. 2009) automatically adapts to the equation of state in order to maintain the desired thermal scale height and thus the absence of any strong dependence on the adiabatic indices is not unexpected.

4.2 Disc evaporation?

The transition from a low-viscosity, high-density outer disc into a high-viscosity, low-density inner disc (See Section 3.5 and Figs 6/7) might provide clues into a long-standing puzzle in accretion physics: How do cool thin discs (Shakura & Sunyaev 1973) transition into hot thick radiatively-inefficient accretion flows (Narayan & Yi 1994) near the black hole? This work shows that magnetically driven winds can lower the disc density and increase the inflow speed with respect to the outer disc. Subsequently, the ions and electrons may become weakly coupled and the cooling time-scale may become limited by the timescale for Coulomb collisions and other plasma processes to equilibrate the temperature of the hot non-radiative ions with the radiatively cooled synchrotron emitting electrons (Shapiro et al. 1976). This may prevent the inner disc from cooling and can conceivably lead to a radiatively inefficient thick accretion flow at a radius r ≲ 25r_g for this setup. Indeed, the finite time-scale for electron–ion coupling implies that one would generally expect such a disc to form for |$\dot{M}\lt \alpha ^2\dot{M}_{\rm Edd}$| (e.g. Esin et al. 1997). This disc-evaporation mechanism, through the elevated α-viscosity in the inner disc, is attractive in that it does not require conduction of heat from the corona to the disc (e.g. Meyer & Meyer-Hofmeister 1994; Liu et al. 1999; Czerny et al. 2000; Qian, Liu & Wu 2007).

How does the high-viscosity, low-β inner disc, seen in Fig. 7(a), form? In our simulation, it may have formed due to the rapid cooling of an initial torus threaded with poloidal magnetic flux (Sikora & Begelman 2013; Begelman & Armitage 2014): the cooling causes the thermal pressure to decrease, but – due to vertical magnetic flux conservation – the magnetic flux stays about the same. This causes the disc to become more strongly magnetized and plasma β to drop. Shearing box simulations seeded with strong vertical magnetic flux appear to develop a similarly highly magnetized accretion state with strong outflows (Salvesen et al. 2016; Bai & Stone 2013). If this scenario is indeed the case, it would require the presence of large scale magnetic flux in the accretion disc prior to the disc becoming thin, which would limit the applicability of this simulation to the intermediate states for X-ray binaries.

Future work will investigate the effect of the different initial magnetic field geometries (e.g. Liska, Tchekhovskoy & Quataert 2018), exploring if large scale poloidal magnetic flux is indeed a necessary ingredient for the high α-viscosity inner disc. It will also include electron–ion coupling, and on-the-fly radiation transfer, to accurately model the cooling of the disc.

4.3 Jet launching

This work shows that thin discs down to at least H/R = 0.03 can efficiently launch relativistic Blandford & Znajek (1977) jets of substantial power, carrying out |${\gtrsim }20{{\ \rm per\ cent}}$| of the accretion power, over time-scales comparable to the accretion time (Section 3.4). This suggests that even such thin discs as considered in this work are capable of retaining for their accretion time large-scale poloidal (vertical) magnetic flux on the BH, a necessary ingredient for launching relativistic jets (Blandford & Znajek 1977). This is particularly interesting given that simple analytical arguments suggest that thin discs should lose their large-scale magnetic flux to outward diffusion (Lubow et al. 1994). It is possible that the large scale external torques may overcome this problem by dragging flux inwards before it has time to diffuse out (see also Guilet & Ogilvie 2012, 2013).

How can we reconcile the formation of powerful jets from thin discs with observations? There are no observations that have convincingly detected jets from thin discs in the soft state of X-ray binaries (though see Rushton et al. 2012), however, about |$10{{\ \rm per\ cent}}$| of quasars are radio loud and form jets (Sikora et al. 2007). Because our simulated jets have wide opening angles, Δθ ∼ 20° (Fig. 2c), they might become less optically thick and more difficult to detect (see also Russell et al. 2011; Fragile et al. 2012). Another possible explanation is that our simulations do not apply to the soft state of XRBs (Section 4.4).

4.4 A transitional disc?

An interesting possibility is that our simulations apply to transitional discs, in the middle of the hard-to-soft state transition (e.g. Fender, Belloni & Gallo 2004). In fact, we set up our simulations in a very similar way: the initial thick torus rapidly cools down to the target thickness, H/R = 0.03, which is much smaller than the initial thickness, H/R ∼ 0.3. Since thick accretion discs may be able to generate and advect large scale poloidal magnetic flux through large-scale dynamo action (e.g. Liska et al. 2018), they are expected to retain a substantial amount of it after their collapse into a thin disc. This can lead both to a highly viscous inner disc that evaporates into an ADAF (Section 4.2) and sustains a strong jet (Section 4.3, see also Ferreira & Pelletier 1993; Ferreira et al. 2006; Sikora & Begelman 2013; Begelman & Armitage 2014). Indeed XRBs in the hard-to-soft state transition are known to produce jets (e.g. Fender et al. 2004), while radio-loud quasars may contain such transitional discs (Tchekhovskoy 2015). Spectral modelling of two-temperature magnetically truncated discs has proven successful in explaining both emission in X-Ray and radio during XRB state transitions (Marcel et al. 2018a,b).

However, as proposed in Lubow et al. (1994), this flux may slowly diffuse out and cause the jet to shut down. In addition, if this large scale poloidal magnetic field indeed leads to evaporation of the inner disc into an ADAF (Section 4.2), the truncation radius between the inner thick and outer thin disc will move inwards. This is consistent with observational evidence of the truncation radius moving in during the evolution towards the soft state (Esin et al. 1997; Done, Gierliński & Kubota 2007; Ingram & Done 2011). Fig. 4(c) indeed shows signs of magnetic flux diffusing out of the BH: the flux in the disc (Φ_disc) stays roughly constant while the flux on the BH (Φ_BH) drops. However, the drop is small and appears to be leveling off. Several mechanisms have been suggested that can prevent the poloidal magnetic flux from diffusing out in thin discs (Rothstein & Lovelace 2008; Guilet & Ogilvie 2012, 2013). Future simulations spanning much longer runtimes can probe if thin discs are able to retain poloidal magnetic flux for a more extended time period, or are always transitional.

5 CONCLUSIONS

In this work, we have performed the thinnest disc GRMHD simulations to date. We started with an H/R = 0.03 accretion disc tilted by 10° relative to a rapidly spinning a = 0.9375 BH. Using three AMR levels, we carried out GRMHD simulations at sufficiently high effective resolution, 2880 × 864 × 1200, which for the first time resolved the MRI turbulence in a thin disc in all three dimensions with near-cubical cells (of order unity aspect ratio). Our results can be summarized in three key points.

First, we have confirmed for the first time that the inner parts of tilted thin discs can align with the BH equatorial plane as theorized 40 years ago by Bardeen & Petterson (1975), even when the full effects of GR, anisotropic MRI turbulence and torquing of the disc by magnetized corona and jets are included. The disc aligns with the BH within the BP radius, r_BP ≃ 5r_g, whose value is expected to increase for thinner discs (e.g. Kumar & Pringle 1985). The development of a BP configuration can have profound consequences for the evolution of BH spins in AGN, as the large lever arm of r_BP out to which the disc is aligned can torque the BH into alignment with the outer, tilted disc on a much shorter timescale than without the BP effect (e.g. Scheuer & Feiler 1996; Natarajan & Pringle 1998).

Second, we have shown that an accretion disc can develop an inner low-density, high-viscosity disc coupled to an outer high-density, low-viscosity disc at r ≲ 25r_g. We suggested that the order unity viscosity of the inner disc we find might lead to it evaporating into a radiatively inefficient accretion flow when the electron-ion coupling time exceeds the accretion time (e.g. Esin et al. 1997). This high viscosity may be caused by the presence in the initial conditions of large-scale poloidal magnetic flux, which removes the angular momentum through large-scale outflows. Large-scale poloidal magnetic flux may be present in thin discs during hard-to-soft state transitions (e.g. Sikora & Begelman 2013; Begelman & Armitage 2014), as discussed in Section 4.4.

Third, we have shown that BH accretion systems with thin discs, if initially threaded with large-scale poloidal magnetic flux, can launch powerful Blandford & Znajek (1977) jets on the viscous timescale, with their power reaching |$20\!-\!50{{\ \rm per\ cent}}$| of the accretion power. This challenges the standard paradigm that thin discs in the soft state cannot advect inwards poloidal magnetic flux needed to launch jets (Lubow et al. 1994) and is seemingly in tension with the lack of any clear detection of jets in X-ray binaries. However the morphology of our jets, specifically their twice as large opening angle as of those produced by thick discs (e.g. McKinney 2006; Liska et al. 2018 Chatterjee et al. 2019), may make them more optically thin and thus more difficult to detect (see also Fragile et al. 2012). Another possibility is that our simulations describe transitional discs in the hard-to-soft state transition which are known to produce powerful jets (e.g. Fender et al. 2004) and, like our simulations, may naturally harbor large scale poloidal magnetic flux (Sikora & Begelman 2013), which is required to produce powerful jets (Blandford & Znajek 1977). This flux may eventually diffuse out (e.g. Begelman & Armitage 2014) causing the jets to shut down. Outwards flux diffusion might indeed be present in our simulation (Section 4.4).

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article: movie file (link).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank Chris Fragile and Cole Miller for useful suggestions. AI thanks James Matthews for useful discussions. This research was made possible by NSF PRAC award no. 1615281 and OAC-1811605 at the Blue Waters sustained-petascale computing project and supported in part under grant no. NSF PHY-1125915. Simulation data is available by writing to [e-mail]. ML and MK were supported by the Netherlands Organisation for Scientific Research (NWO) Spinoza Prize, AI by the Royal Society URF, AT by Northwestern University, NSF grant 1815304 and NASA grant 80NSSC18K0565, and the TAC and NASA Einstein (grant no. PF3-140131) postdoctoral fellowships.