Disc Tearing and Bardeen-Petterson Alignment in GRMHD Simulations of Highly Tilted Thin Accretion Discs

TL;DR

This study uses GRMHD simulations to explore how highly tilted thin accretion discs around rapidly spinning black holes undergo tearing and alignment, revealing complex precession behaviors and jet formation.

Contribution

First GRMHD simulations of very thin, highly tilted accretion discs around spinning black holes showing disc tearing and differential precession.

Findings

Inner disc aligns with black hole equator at small radii.

Disc tearing occurs due to frame-dragging overpowering viscosity.

Relativistic jets are produced with precession behaviors.

Abstract

Luminous active galactic nuclei (AGN) and X-Ray binaries (XRBs) often contain geometrically thin, radiatively cooled accretion discs. According to theory, these are -- in many cases -- initially highly misaligned with the black hole equator. In this work, we present the first general relativistic magnetohydrodynamic simulations of very thin (h/r~0.015-0.05) accretion discs around rapidly spinning (a~0.9) black holes and tilted by 45-65 degrees. We show that the inner regions of the discs with h/r<0.03 align with the black hole equator, though out to smaller radii than predicted by analytic work. The inner aligned and outer misaligned disc regions are separated by a sharp break in tilt angle accompanied by a sharp drop in density. We find that frame-dragging by the spinning black hole overpowers the disc viscosity, which is self-consistently produced by magnetized turbulence, tearing the disc apart and forming a rapidly precessing inner sub-disc surrounded by a slowly precessing outer sub-disc. We find that the system produces a pair of relativistic jets for all initial tilt values. At small distances the black hole launched jets precess rapidly together with the inner sub-disc, whereas at large distances they partially align with the outer sub-disc and precess more slowly. If the tearing radius can be modeled accurately in future work, emission model independent measurements of black hole spin based on precession-driven quasi-periodic oscillations may become possible.

Tables

Table 1. Table 1 : The tilt ( 𝒯 init subscript 𝒯 init \mathcal{T}_{\rm init} ), number of cells in r − limit-from 𝑟 r- , θ − limit-from 𝜃 \theta- and ϕ − limit-from italic-ϕ \phi- coordinates ( N r × N θ × N ϕ subscript 𝑁 𝑟 subscript 𝑁 𝜃 subscript 𝑁 italic-ϕ N_{r}\times N_{\theta}\times N_{\phi} ), disc thickness ( h / r ℎ 𝑟 h/r ) and time interval ( t i − t f subscript 𝑡 𝑖 subscript 𝑡 𝑓 t_{i}-t_{f} ) for each model.

Model	Full name	${𝒯}_{\mathrm{init}}$	$N_r\times N_{\theta }\times N_{\varphi }$	$h/r$	$t_i$-$t_f$ [${10}^4{t}_g$]
T45	T45HR03L	${45}^{\circ }$	2880$\times$864$\times$1200	$0.03$	$0$-$10.5$
T45-H	T45HR03H	${45}^{\circ }$	5760$\times$1728$\times$2400	$0.03$	$4.8$-$7.2$
UT45	T45HR015H	${45}^{\circ }$	5760$\times$1728$\times$2400	$0.015$	$4.8$-$7.0$
MT45	T45HR05L	${45}^{\circ }$	2880$\times$864$\times$1200	$0.05$	$0$-$4.5$
T65	T65HR03L	${65}^{\circ }$	2880$\times$864$\times$1200	$0.03$	$0$-$12$

Full text

††pagerange: Disc Tearing and Bardeen-Petterson Alignment in GRMHD Simulations of Highly Tilted Thin Accretion Discs–References††pubyear: 2015

1 Introduction

Because the gas supply of supermassive black holes (BHs) originates from far away, the angular momentum vector of the accreting gas will most likely be independent of the BH spin vector. In stellar-mass BHs, asymmetric supernova kicks can also lead to substantial misalignment between the two vectors. If the relative orientation of the two is random, the resulting accretion disc would half of the time make a 60-degree or greater angle relative to the BH equator. Tilted discs are expected in many luminous systems ranging from X-ray binaries (XRBs), active galactic nuclei (AGN), tidal disruption events (TDEs), and binary merger remnant discs (e.g. Hjellming & Rupen, 1995; Orosz et al., 2001; Greene et al., 2001; Caproni et al., 2006, 2007; Abbott et al., 2017). The physics of tilted accretion systems is of crucial importance for understanding the growth of most – if not all – supermassive BHs throughout cosmological time, and has profound implications for jet production, BH spin measurements, and energy release of BH accretion systems (e.g. Natarajan & Pringle 1998; Stella & Vietri 1998; Fiacconi et al. 2018).

General relativistic frame dragging by a rotating BH warps tilted accretion discs. The evolution of such warped discs depends sensitively on the ratio between their viscosity and their dimensionless scale height, or aspect ratio (Papaloizou & Pringle, 1983). In the analytic approach, the disc viscosity, typically parameterized through the $\alpha$-viscosity parameter, regulates the transport of angular momentum in the disc (Shakura & Sunyaev, 1973). The disc aspect ratio, $h/r$, the ratio of disc scale height $h$ to radius $r$, characterizes the internal pressure of the disc. When the disc is relatively thick, $h/r>\alpha$, such as in the hard state of XRBs and low-luminosity AGN, the warps are transmitted by pressure waves traveling at about half the speed of sound (Papaloizou & Lin, 1995). In this wave-like limit, analytic calculations (Ivanov & Illarionov, 1997; Lubow et al., 2002) and general-relativistic magnetohydrodynamic (GRMHD) simulations (Fragile et al., 2007; Morales Teixeira et al., 2014; Liska et al., 2018b, 2019a; White et al., 2019) have shown that the tilt of the disc oscillates as a function of radius within about 20 gravitational radii, such that material gets accreted at high inclination angles. When the disc is relatively thin, $h/r<\alpha$, warps are propagated through viscous diffusion, and radial tilt oscillations get damped by the dissipative effects of disc viscosity. This is thought to be the case in bright quasars and intermediate/soft state X-ray binaries. In this limit, which is the focus of this work, a so-called Bardeen-Petterson configuration is expected to emerge: the outer disc remains tilted whereas the inner disc aligns for $\theta\lesssim\pi/2$ (Bardeen & Petterson, 1975), or counter aligns for $\theta\gtrsim\pi/2$. Since, by Newton’s third law, the torque exerted by the BH on the inner disc is exactly equal and opposite to the torque exerted by the inner disc on the BH, Bardeen & Petterson (1975) alignment will not only align or counteralign the inner disc with the BH spin on shorter timescales, but will also torque the BH into alignment with the total angular momentum vector of the system on longer timescales. In cases where the disc contains most of the angular momentum, the Bardeen & Petterson (1975) effect substantially accelerates the alignment between disc and BH spin axes. This can possibly lead to much faster spin-up of supermassive BHs compared to direct accretion of angular momentum through the BH’s event horizon (e.g. Natarajan & Pringle 1998; King et al. 2005). However, as noted in King et al. (2005), in cases where the disc’s angular momentum is less than the black hole’s angular momentum and $\theta\gtrsim\pi/2$, the disc will be torqued into counteralignment with the black hole. In such cases rapid spin-down of the central black hole is expected. This may lead to a lower than expected spin for very massive SMBHs in a chaotic accretion scenario (King & Pringle, 2006).

However, across both disc radius and height, there tend to be large, non-linear, and anisotropic variations of the viscous stresses induced by the magneto-rotational instability (MRI, Balbus & Hawley 1991, 1998) in magnetized accretion discs, defying the simple $\alpha-$viscosity prescription (Penna et al., 2010; Sorathia et al., 2010; McKinney et al., 2012; Jiang et al., 2017). For example, recent GRMHD simulations (Liska et al., 2019b) of a very thin, $h/r=0.03$, magnetized accretion disc tilted by $10^{\circ}$ have shown that magnetic fields launch winds which counteract Bardeen & Petterson (1975) alignment in a strongly non-linear fashion, producing an aligned region that is much smaller than predicted for $\alpha-$discs (Kumar & Pringle, 1985; Nelson & Papaloizou, 2000; Lodato & Price, 2010; Nixon & King, 2012).

Crucially, when a thin *$\alpha-$*disc is tilted by $\mathcal{T}\gtrsim 45^{\circ}$, both simple 1D-evolution and smoothed particle hydrodynamics (SPH) simulations of warped discs suggest that frame-dragging by the spinning BH tears the disc apart into differentially precessing rings (Lodato & Pringle, 2006; Nixon et al., 2012b; Nealon et al., 2015). Just as in thin $\alpha-$discs, once a magnetized disc starts to tear, the viscosity may drop and encourage further tearing (Ogilvie, 1999; Nixon & King, 2012; Doǧan et al., 2018). This process has been observed in 1D calculations of warped disc evolution (Lodato & Pringle, 2006). However, since the tilt angle exceeds the disc’s angular thickness by more than an order of magnitude, the warp becomes highly non-linear – requiring a detailed treatment involving the 3D magnetized turbulence that is the glue that holds the disc together. GRMHD simulations make this treatment possible, and we use them in this paper to study whether and how tilted, thin, magnetized discs get torn apart. We pay special attention to the typical tearing radius and the physical prerequisites for tearing to occur. Forming such an understanding may pave the way for BH spin measurements based on quasi-periodic oscillations (QPOs) observed in XRB lightcurves (van der Klis, 1989). Such QPOs could be disc (Stella & Vietri, 1998; Ingram et al., 2009, 2016; Franchini et al., 2017; Motta et al., 2018) and/or jet (Kalamkar et al., 2016; Stevens & Uttley, 2016) precession.

In this work we present the first GRMHD simulations of highly tilted thin accretion discs in the diffusive limit of warp propagation ($h/r<\alpha$). In Section 2 we describe our code and initial conditions. We present our results in Sec. 3 and conclude in Sec. 4.

2 Numerical Method and Initial Conditions

For this work we use our recently developed GRMHD code H-AMR (Liska et al., 2018b; Chatterjee et al., 2019; Liska et al., 2019c). It evolves the GRMHD equations with a finite-volume–based method in modified Kerr-Schild coordinates (as in Gammie et al., 2003) and uses a constrained-transport scheme for magnetic field evolution (see Gardiner & Stone, 2005). Here we employ a logarithmic spherical polar grid with $3$ to $4$ levels of adaptive mesh refinement (AMR) and $4$ levels of local adaptive time-stepping, allowing us to focus the resolution on the regions of interest. In particular, magnetized turbulence in the disc needs to be resolved (Liska et al., 2018b), so we use rest-mass density $\rho$ as the refinement criterion in order to delineate the disc (as in Liska et al., 2019b). In this way, we achieve the following effective resolutions in spherical polar coordinates ($N_{r}\times N_{\theta}\times N_{\phi}$): $2880\times 864\times 1200$ in our low-resolution models, and, by doubling the resolution in every dimension, $5760\times 1728\times 2400$ in our-high resolution models (see Table 1). This resolves our thin discs by approximately $7$ to $14$ cells per scale height (see Sec. 3 for the MRI quality factors). We use outflow boundary conditions at the inner and outer radial boundaries, which we place inside the event horizon and at $r=10^{5}r_{g}$, respectively, where $r_{g}=GM/c^{2}$ is the gravitational radius. This way both boundaries are causally disconnected from the accretion system. Across the polar singularity in $\theta$ we use a transmissive boundary condition, which we implemented using a multi-faceted method that minimizes numerical dissipation in the polar region (for details, see Liska et al., 2018b).

All models, shown in Table 1, are initialized with a Kerr black hole (with $a=0.9375$) surrounded by a torus in hydrostatic equilibrium (in accordance with Fishbone & Moncrief, 1976) with its inner edge at $r_{\rm in}=12.5r_{g}$, its pressure maximum at $r_{\rm max}=25r_{g}$ (this results in torus outer edge located at $r_{\rm out}=200r_{g}$), and its density normalised by setting $\max\rho=1$. We use the equation of state of an ideal gas, $p_{g}=(\Gamma-1)u_{g}$, where $p_{g}$ and $u_{g}$ are thermal pressure and thermal energy density, and we use a polytropic index that corresponds to a non-relativistic monoatomic ideal gas, $\Gamma=5/3$. We seed the torus with a poloidal magnetic field defined by a covariant vector potential $A_{\phi}=(\rho-0.05)^{2}r^{3}$. We normalize the magnetic field strength by requiring that $\max p_{g}/\max p_{b}=30$, where $p_{b}$ is the magnetic pressure. Subsequently, we tilt the torus and magnetic field relative to the BH spin (and the grid) by an angle $\mathcal{T}_{\rm init}$ (see Liska et al. 2018b for details). Finally, we reduce the disc thickness to a target scale-height $h/r$ by cooling the gas at a rate slow enough to avoid disruption of the disc orbital dynamics. We do this by letting the internal energy decay throughout the simulation exponentially over time, with the time constant set by the orbital timescale (Noble et al., 2009). Since the disc needs some time to cool and reach the target scale height, we only include data after $t=10^{4}t_{g}$ into our analysis, where $t_{g}=r_{g}/c$. We initialize the high-resolution models T45-H and UT45 with a well-evolved state of model T45 at $t=4.8\times 10^{4}t_{g}$ (see Table 1 for details). In the case of T45-H, we, additionally, reduce the target thickness in the cooling function from $h/r=0.03$ to $h/r=0.015$.

3 Results

Figure 1 shows a colour map of the density for models UT45, T45-H, MT45 and T65 in their evolved state at $t\gtrsim 4\times 10^{4}r_{g}/c$. As can be seen from 3D animations (see this YouTube playlist), the inner disc always aligns with the BH spin in models T45-H and UT45, aligns most of the time in model T45, some of the time in model T65, and does not align at all in model MT45. This can be understood because the disc is thicker in model MT45, $h/r=0.05$, and may fall outside the diffusive warp propagation regime, $h/r<\alpha$, where the Bardeen & Petterson (1975) alignment is expected (e.g. Ivanov & Illarionov 1997).

For our thinnest discs, we can establish the dependence of disc structure on thickness. In particular, we observe that the Bardeen & Petterson (1975) alignment radius, $r_{\rm bp}$, increases from $r_{\rm bp}\sim 5r_{g}$ at $h/r=0.03$ (models T45 and T45-H) to $r_{\rm bp}\sim 10r_{g}$ at $h/r=0.015$ (model UT45). While this is consistent with the predicted analytic scaling of $r_{\rm bp}\sim(h/r)^{-8/7}$ (e.g. Kumar & Pringle 1985), the proportionality constant for our simulation is much smaller than in the analytic scaling. This discrepancy might be due to the torques applied on the disc by the large-scale magnetically-powered disc outflows that can counteract the Bardeen & Petterson (1975) alignment and reduce the value of $r_{\rm bp}$ (Liska et al., 2019b).

We note that there are important qualitative differences in the simulation results presented here and our previous work, which considered a disc with exactly the same initial conditions but tilted by a much smaller angle, $\mathcal{T}=10^{\circ}$ (Liska et al., 2019b). For this smaller value of tilt, the inner aligned and outer misaligned parts of the accretion disc were separated by a smooth warp (Liska et al., 2019b). However, for much larger tilt angles considered here, $\mathcal{T}=45{-}65$ degrees, the inner aligned and outer misaligned parts of the disc are separated by a discontinuity in tilt angle and sharp drop in density. This transition is referred to as a disc break (e.g., Lodato & Price 2010; Nixon & King 2012). As we discuss in Sec. 4, the development of a break will likely change the physical properties and observational signatures of a tilted disc and its precession.

The Bardeen & Petterson (1975) alignment is variable in time: as seen from the movies, as a general trend, an increase/decrease in BH mass accretion rate $\dot{M}_{\rm BH}$ (Fig. LABEL:fig:timeplotplota) tends to accompany the periods of misalignment/alignment of the inner disc with the BH in models T45 and T65. This might be due to rapid accretion giving the disc insufficient time to align with the BH spin axis (as postulated in e.g. Nealon et al. 2015; Liska et al. 2019b). Alternatively, the formation of a break may suppress accretion making the drop in mass accretion rate a consequence of the alignment instead of a cause.

Figures 1 and 3 show that at large values of tilt, $\mathcal{T}=65^{\circ}$, in model T65, the disc tears into $2$ or sometimes even $3$ differentially precessing sub-discs at radii ranging from $5r_{g}$ to $30r_{g}$. This is caused by the differential Lense & Thirring (1918) torques exceeding the viscous torques that hold the disc together. We observe that the disc also tears at lower tilt values, $\mathcal{T}=45^{\circ}$, in models T45 and UT45: however, this happens at later times, $t\gtrsim 10^{5}r_{g}/c$ for T45 and $t\gtrsim 7\times 10^{4}r_{g}/c$ for UT45, respectively. To reach the BH, the gas must pass between sub-discs via streamers, which are tentacle-like low density structures connecting adjacent sub-discs, as seen in Figs 1 and 3. Note that the BP-aligned part of the inner sub-disc in Fig 1 for T65 is fed partially by the misaligned outer part of the inner sub-disc and partially directly by the middle sub-disc, bypassing the misaligned part of the inner sub-disc. We leave quantification of mass transfer rate between the different sub-discs to future work.

As the phase difference in precession angle $\mathcal{P}$ between the inner and outer sub-disc(s) builds up (Fig. LABEL:fig:timeplotplotc), this contributes to effective cancellation of angular momentum between the two (Nixon et al. 2012a, b). This might explain a factor of $1.5{-}3$ enhancements in BH mass accretion rate (Fig. LABEL:fig:timeplotplota) at $t\sim 5\times 10^{4}r_{g}/c$ and $t\sim 15\times 10^{4}r_{g}/c$ in model T65, corresponding to the moments in time when $\mathcal{P}_{\rm inner}-\mathcal{P}_{\rm outer}=180+360n$, $n=0,1$ degrees, i.e., when the sub-discs are counter-rotating. Interestingly, the evolution of a sub-disc of radius $10r_{g}$ formed at $t\sim 4.5\times 10^{4}r_{g}/c$ (see the left panel of Fig. 3) differs from the medium-sized sub-disc of radius $20-30r_{g}$ present throughout the simulation. While the smaller sub-disc slowly accretes into the BH, the medium-sized sub-disc temporarily merges with the outer sub-disc around $t\sim 8\times 10^{4}r_{g}/c$ for a duration of $t\sim 3.0\times 10^{4}r_{g}/c$ (both in precession angle $\mathcal{P}$ and tilt angle $\mathcal{T}$, see Fig. LABEL:fig:timeplotplotc,d) before tearing again, but this time at a slightly larger radius of $r\sim 30r_{g}$. This merging may be driven by density fluctuations in the inner disc. Namely, during each precession cycle a significant portion of the inner disc mass falls into the BH (viscous timescale of inner disc equals $\sim 1-3\times 10^{4}r_{g}/c$), which reduces the differential Lense & Thirring (1918) precession rate between inner and outer disc, allowing them to merge and subsequently tear at a different radius.

To gain a better insight into the internal disc dynamics in models MT45, UT45 and T65 we analyze radial profiles of density, plasma $\beta=p_{g}/p_{B}$ and $\alpha-$viscosity, as shown in Fig. 4(a)-(d). Here, all vectors are calculated in a coordinate system ($r$, $\tilde{\theta}$, $\tilde{\phi}$) aligned with the local rotation axis of the disc (see Liska et al. 2019b for details). Since $\beta\gtrsim 1$ for $r>10r_{g}$ the disc remains (mostly) gas pressure dominated, except where the disc breaks or tears and the density drops. As expected for gas pressure dominated discs, we verified that the density weighted scale height, $(h/r)_{\rho}=\langle\tilde{\theta}-\langle\tilde{\theta}\rangle_{\rho}\rangle_{\rho}$, matches the thermal scale height, $(h/r)_{\rm thermal}=\langle c_{s}\rangle_{\rho}/\langle v_{k}\rangle_{\rho}$, where $v_{k}$ is the Keplerian 3-velocity, $c_{s}$ is the sound speed, and $\langle\dots\rangle_{q}$ indicates an angle-average weighted by the quantity $q$. Throughout the disc in models UT45 and T65 the sum of the Maxwell, $\alpha_{M}=b^{r}b^{\tilde{\phi}}/(p_{g}+p_{b})$, and Reynolds, $\alpha_{R}=\rho u^{r}u^{\tilde{\phi}}/(p_{g}+p_{b})$, stress contributions to the viscosity parameter remains larger than the disc scale height. This confirms that these two discs are in the $h/r<\alpha$ diffusive warp propagation regime (Papaloizou & Pringle, 1983). Here, $b^{\mu}$ and $u^{\mu}$ are the magnetic and velocity 4-vectors.

Surprisingly, the effective viscosity parameter, $\alpha_{\rm eff}=-v_{r}v_{k}/c_{s}^{2}$, which is a measure of the radial inflow speed, exceeds the sum of $\alpha_{R}$ and $\alpha_{M}$ by more than an order of magnitude. This indicates that the angular momentum transport in turbulent discs cannot be described by stresses induced through a local $\alpha-$viscosity. Most likely, large-scale magnetic and/or internal torques contribute to this discrepancy by transporting angular momentum outwards. The situation is similar to a model of a weakly tilted disk of Liska et al. (2019b, this model started out with a tilt angle of $\mathcal{T}=10^{\circ}$, but at late times the tilted reduced to $\mathcal{T}\sim 2^{\circ}$), but with several times smaller discrepancy between the effective viscosity parameter and combined Reynolds/Maxwell stresses than in this work. This suggests that tilt-related effects, such as viscous dissipation in warps (e.g. Lodato & Pringle 2006) and/or spiral shocks aligned with the line of nodes (Fragile & Blaes, 2008; White et al., 2019), are likely the more dominant dissipation mechanism at large tilt angles considered in this work ($\mathcal{T}\gtrsim 45^{\circ}$). Additionally, in our model T65, cancellation of angular momentum where the sub-discs intersect could explain the enhanced accretion rate during episodes of large misalignment between adjacent sub-discs (see Fig. LABEL:fig:timeplotplot). Observationally, the large, above unity, effective viscosity in the inner disc may cause accretion to proceed so fast that the electrons and ions do not have the time to equilibrate their temperatures, leading to a two-temperature plasma (Esin et al. 1997, see also Sec. 4), forming a rapidly precessing advection-dominated accretion flow that produces a hard spectrum (Narayan & Yi, 1994).

To verify numerical convergence, we compared the radial profiles of our model T45 and model T45-H, carried out a twice as high resolution, and found a good level of agreement. In addition, the number of cells per MRI wavelength, $Q_{r}\times Q_{\tilde{\theta}}\times Q_{\tilde{\phi}}$, saturates around $100\times 100\times 500$ for $r\lesssim 20r_{g}$ and $10\times 10\times 150$ for $r\gtrsim 20r_{g}$ in our low resolution models, indicating that the MRI turbulence is well-resolved in the inner and reasonably well-resolved in the outer regions of the accretion disc (Sorathia et al., 2010; Shiokawa et al., 2012).

All our models launch moderately strong Blandford & Znajek (1977) jets, which we define as magnetically-dominated regions with $p_{b}/\rho c^{2}>5$. Figure LABEL:fig:timeplotplot(b) shows that their energy outflow efficiency, or power measured in units of accretion power $\dot{M}_{\rm BH}c^{2}$, reaches $\eta_{\rm jet}\sim 1\%-10\%$. Interestingly, the jets appear to follow the orientation of the disc over a wide range of length scales. For instance, jets at small radii align with the inner disc, as seen in Figs 3 and 4(g-i). At larger radii, they get torqued via the corona, (roughly) defined as all gas not part of the jet and whose density is a factor $10^{3}$ smaller than that in the disc (see also Liska et al. 2019a), into (partial) alignment with the outer disc. This is also seen in Fig. 3. In fact, in model T65, the disc-jet interaction can be so strong that the jets running into outer sub-disc can push it into higher orbits (see this YouTube playlist). This may deprive the BH of its mass supply and quench the accretion at times beyond those simulated.

In nature, the outermost sub-disc would be essentially infinite in size, definitely much larger than in our simulations, and would have an essentially infinite precession period (Liska et al., 2018b). How does this affect potential quasi-periodic signals coming from such systems? Due to the lack of precession of the extremely large outermost sub-disc, all precessing sub-discs would be located at smaller radii. Thus, only the parts of the jet outside of the aligned region, $r\gtrsim r_{\rm bp}$, and inside of the outermost, non-precessing sub-disc would be expected to contribute to jet-driven QPOs (e.g., Kalamkar et al. 2016; Stevens & Uttley 2016).

4 Discussion and conclusion

In this work we have presented the first GRMHD simulations of highly tilted ($\mathcal{T}=45{-}65$ degrees), thin ($h/r=0.015{-}0.05$) accretion discs around rapidly spinning BHs ($a=0.9375$). We demonstrate for the first time that in the presence of realistic magnetized turbulence the inner parts of such discs can align with the BH spin axis, as predicted by Bardeen & Petterson (1975). However, we find that the alignment radius of $r_{\rm bp}\lesssim 5{-}10r_{g}$ (Fig. 1) is much smaller than predicted by analytic models. This discrepancy may be caused by a very large effective viscosity $\alpha_{\rm eff}$ giving the inner disc insufficient time to align (see Liska et al. 2019b for discussion).

The Bardeen & Petterson (1975) effect can torque the BH into alignment with the outer disc (as long as the disc’s angular momentum exceeds the BH’s angular momentum, see Natarajan & Pringle, 1998; King et al., 2005; Fiacconi et al., 2018) and lead to rapid BH spinup in such systems. Rapidly spinning black holes are known to launch powerful Blandford & Znajek (1977) jets in the presence of large scale poloidal or toroidal magnetic flux (e.g., Tchekhovskoy et al. 2011; McKinney et al. 2012; Liska et al. 2018a). Consistent with semi-analytic work (Nixon & King, 2012) and smoothed-particle hydrodynamics (SPH) simulations (Lodato & Price, 2010) the transition between the inner aligned and outer misaligned discs occurs over a very short distance and exhibits a sharp drop in density. As gas crosses this ‘break’, misaligned angular momentum cancels rapidly, possibly leading to enhanced dissipation of kinetic and magnetic energy causing non-thermal emission. We note that the absence of a break in recent MHD simulations of a thin ($h/r\sim 0.05$) disk tilted by $\mathcal{T}\sim 24^{\circ}$ (Hawley & Krolik, 2019) is not inconsistent with our work, as evidenced by the absence of a break in model MT45.

When the disc tilt is large, $\mathcal{T}\gtrsim 45^{\circ}$, we show for the first time that a magnetized thin $h/r=0.015{-}0.03$ disc can tear into multiple independently precessing sub-discs, as seen in Fig. 3. In future work we will investigate whether tearing can happen at smaller values of disc tilt for thinner discs and whether, therefore, even weakly misaligned discs in XRBs and AGN can be subject to disc tearing. Observationally, we expect disc tearing to lead to a wide range of interesting phenomena. Differential precession, as explained in Sec. 3, can lead to cancellation of angular momentum leading to a factor of few increase in the mass accretion rate. This may explain flaring in the hard-intermediate/ultra-luminous state of XRBs (Remillard & McClintock, 2006; McClintock & Remillard, 2006). For instance if a disc undergoes several tearing events in quick succession, its luminosity will increase. However, the inner disc density eventually drops, because the outer disc is unable able to keep up the supply of gas indefinitely, especially when a powerful jet injects energy and angular momentum into the outer disc (Sec. 3). This drop in density, together with a very short accretion time (due to unusually large effective viscosity, $\alpha_{\rm eff}\gtrsim 1$, in Fig. 4c,d), may lead to the decoupling of ions and electrons into a two-temperature plasma, reducing the ability of ions in the disc to cool and puffing up the inner thin disc into a hotter, and less radiatively efficient, thick advection-dominated accretion flow (ADAF, see Narayan & Yi 1994). In fact, for $\alpha_{\rm eff}>1$ a thick disc would be expected to always form when two-temperature thermodynamics effects are taken into account (see e.g. Esin et al. 1997; Ferreira et al. 2006; Marcel et al. 2018a, b; Liska et al. 2019b). Since the viscous torque is stronger for a larger disc thickness, the thicker disc would no longer undergo tearing. For instance, in agreement with Hawley & Krolik (2019), our thicker disc model MT45 with $h/r=0.05$, seen in the bottom-left panel of Fig. 1, does not show any signs of tearing. In the absence of disc tearing, the connection with the outer thin disc gets reestablished. The outer disc then feeds the inner one, and the BH mass accretion rate rises. This results in the density increase of the inner disc and its rapid cooling and collapse into a thin disc. The cycle then repeats. During each such cycle magnetic jets violently interact with the precessing sub-discs, making the tearing radius an interesting location for enhanced dissipation and (non-)thermal emission. Additionally, the streamers connecting torn sub-discs to each other can scatter and/or reradiate the emission from the central regions, substantially affecting the emergent spectrum and variability and making the discs appear larger than otherwise. This may resolve the puzzle of what makes AGN disc sizes exceed the predictions of an $\alpha-$disc model (e.g., Blackburne et al., 2011).

It has been suggested that BH spin could be measured based on precession induced Type-C QPOs (Stella & Vietri, 1998; Ingram et al., 2009; Franchini et al., 2017; Motta et al., 2018), whose frequency depends on BH spin and disc size. Such measurements would be able to independently verify the accuracy of the continuum fitting (e.g. McClintock et al. 2014) and iron-line methods (e.g. Reynolds & Fabian 2008) without making any assumptions about the disc’s emission near the innermost stable circular orbit (ISCO) or assuming that the system is aligned. Making self-consistent predictions for the tearing radius, based on e.g. disc thickness, tilt and magnetic field topology, requires clear theoretical understanding of the physics driving disc tearing. On a basic level, for a disc to tear, the differential Lense & Thirring (1918) torques need to exceed the viscous torques holding the disc together. For $\alpha$-discs, the effective torque counteracting breaking and tearing can be derived as function of warp amplitude (Ogilvie, 1999) making it possible to calculate criteria for disc breaking and tearing (Doǧan et al., 2018). However, disc tearing in GRMHD leads to a substantially different morphology compared to SPH simulations. More specifically, instead of tearing up into narrow ‘rings’ with $\Delta r\sim h$ as seen in SPH models in the $h/r<\alpha$ regime (Nixon et al., 2012b), our GRMHD models form radially extended sub-discs with $\Delta r\gg h$ (Fig. 3). We note that while SPH simulations in the thick disk, $h/r>\alpha$, regime form rings which are more radially extended (Nealon et al., 2015), their radial extent is still much smaller than presented in this work. This discrepancy suggests the disc tearing process may be more complicated when magnetized turbulence self-consistently determines the viscous-like coupling within the disc. We hypothesize that this might be caused by radial tension along magnetic field lines, which is neglected in the $\alpha$-disc approximation. These and other questions, such as how radiation pressure influences the behaviour of disc tearing into sub-discs, how a disc break influences the angular momentum transport, and what causes the apparent stochasticity in the disc tearing radius (see Sec. 3), will be addressed in future work.

5 Acknowledgments

This research was made possible by NSF PRAC awards no. 1615281 and OAC-1811605 as part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. ML and MK were supported by the NWO Spinoza Prize, AI by the Royal Society URF, CH by the NWO Research Talent grant (no. 406.18.535), SM by the NWO VICI grant (no. 639.043.513), and AT by the NSF grants 1815304, 1911080 and NASA grant 80NSSC18K0565. The simulation data presented in this work is available upon request to AT at atchekho@northwestern.edu.

6 Supporting Information

Additional Supporting Information may be found in the online version of this article: movie files. See our YouTube playlist for 3D visualizations of all models.