Formation of precessing jets by tilted black hole discs in 3D general relativistic MHD simulations

Abstract

Gas falling into a black hole (BH) from large distances is unaware of BH spin direction, and misalignment between the accretion disc and BH spin is expected to be common. However, the physics of tilted discs (e.g. angular momentum transport and jet formation) is poorly understood. Using our new GPU-accelerated code h-amr, we performed 3D general relativistic magnetohydrodynamic simulations of tilted thick accretion discs around rapidly spinning BHs, at the highest resolution to date. We explored the limit where disc thermal pressure dominates magnetic pressure, and showed for the first time that, for different magnetic field strengths on the BH, these flows launch magnetized relativistic jets propagating along the rotation axis of the tilted disc (rather than of the BH). If strong large-scale magnetic flux reaches the BH, it bends the inner few gravitational radii of the disc and jets into partial alignment with the BH spin. On longer time-scales, the simulated disc–jet system as a whole undergoes Lense–Thirring precession and approaches alignment, demonstrating for the first time that jets can be used as probes of disc precession. When the disc turbulence is well resolved, our isolated discs spread out, causing both the alignment and precession to slow down.

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

1 INTRODUCTION

The angular momentum of matter accreting on to a spinning black hole (BH) is expected to often be misaligned with the BH spin, in wide range of systems including X-ray binaries, active galactic nuclei (AGNs), tidal disruption events (TDEs), and BH–neutron star mergers, with observations indicating the presence of such discs in some systems (e.g. Hjellming & Rupen 1995; Caproni et al. 2006). Disc tilt has been invoked to explain quasi-periodic oscillations (QPOs) in BH systems (Stella & Vietri 1998; Ingram, Done & Fragile 2009) and variability in the jet orientation (Agudo et al. 2012; Stone & Loeb 2012; Lister et al. 2013, 2016; Tchekhovskoy et al. 2014). Without an accurate theoretical description of tilted accretion flows, 3D MHD simulations in general relativity (GR) are an excellent tool for understanding these important systems.

An effect of crucial importance for tilted accretion is Lense–Thirring (LT) precession (Thirring 1918). In GR, massive rotating objects distort nearby inertial frames (frame-dragging). Such twisting induces nodal precession in test particles on tilted orbits, depending on their distance to the central object (precession frequency ∝ 1/r³). We characterize warps in the disc in terms of the precession angle |$\mathcal {P}(r)$| with respect to the initial disc orientation and the tilt angle |$\mathcal {T}(r)$| with respect to the BH spin. MHD effects allow perturbations in |$\mathcal {T}(r)$| and |$\mathcal {P}(r)$| to travel radially, affecting the overall disc behaviour.

Disc warps can travel by means of pressure waves and viscosity (parametrized by the viscosity parameter α), which is generated by magnetic turbulence. In thin high-viscosity discs, with scaleheight H/R < α, pressure waves are damped and warps travel by viscous diffusion (Papaloizou & Pringle 1983). Here, we work in the opposite limit of thick low-viscosity discs, with H/R > α. In such discs, pressure waves can travel freely at about half the sound speed (Papaloizou & Lin 1995). In this wave-like limit, |$\mathcal {T}(r)$| oscillates in radius at r ≲ 20 r_g (Demianski & Ivanov 1997; Ivanov & Illarionov 1997; Lubow & Ogilvie 2000; Lubow, Ogilvie & Pringle 2002). This prediction was confirmed in simulations of tilted accretion discs using GR hydrodynamics (Mewes et al. 2016) and GRMHD (Fragile & Anninos 2005; Fragile et al. 2007). The latter simulations also showed LT precession of the disc around the BH at a constant rate (constant |${\rm d}\mathcal {P}/{\rm d}t$|⁠), with the disc behaving as a rigid body. Interestingly, the precession frequency was found to be consistent with observed Type-C QPOs.

On the other hand, jets can strongly affect the accretion flow and its orientation: when the magnetic flux on the BH is strong enough to obstruct the inner disc infall (Narayan, Igumenshchev & Abramowicz 2003), the associated jets extract large amounts of rotational energy from the BH and the disc (Tchekhovskoy, McKinney & Narayan 2012; Tchekhovskoy 2015) and compress the disc vertically (Tchekhovskoy, Narayan & McKinney 2011; McKinney, Tchekhovskoy & Blandford 2012; Tchekhovskoy & McKinney 2012). In this so-called magnetically arrested disc (MAD) regime, the jets can force the inner parts of tilted thick, radially extended accretion discs to align with the BH spin (McKinney, Tchekhovskoy & Blandford 2013; Polko & McKinney 2017). However, the behaviour of jets produced by tilted discs with smaller magnetic fluxes, in the so-called standard and normal evolution (Narayan et al. 2012) regime, and/or smaller radial extents remains completely unexplored.

Here, we study tilted thick disc–jet systems for a range of magnetic field strength and with disc size changing substantially over time, using first-principles GRMHD simulations. We describe our numerical method in Section 2 and our numerical set-up in Section 3. We present our results in Section 4 and conclude in Section 5.

2 H-AMR (HAMMER) CODE

We use a new massively parallel 3D GRMHD code h-amr (pronounced ‘hammer’) accelerated by Graphical Processing Units (GPUs). We developed h-amr based on a 2D serial open-source code harm2d (Gammie, McKinney & Tóth 2003; Noble et al. 2006). h-amr performs 10 times faster on a GPU than on a 16-core CPU. h-amr is parallelized via MPI with domain decomposition and scales well to thousands of GPUs, achieving weak scaling efficiency of 85 per cent on 4096 GPUs for a tile size of 100³ cells on the Blue Waters supercomputer (see Supporting Information [SI] Section 2.1).

h-amr features a staggered grid for constrained transport of magnetic field (Gardiner & Stone 2005), adaptive mesh-refinement (AMR, not used adaptively here), and a locally adaptive time-step (details to be described in future work). These specialized features substantially speed up the simulations, which we carry out on a spherical polar grid (see also SI Section 2.3). We use outflow boundary conditions in the r-direction, transmissive polar boundary conditions in the θ-direction and periodic boundary conditions in the ϕ-direction. The radial grid is uniform in log r and extends from just inside of the event horizon out to 10⁵r_g, where r_g = GM/c² is the gravitational radius, such that the outer boundary is causally disconnected over the simulation duration, ≳10⁵t_g, where t_g = r_g/c. If our grid were uniform in ϕ and θ, cells would become prohibitively small near the polar axis. To mitigate this problem, we adopt the approach of Tchekhovskoy et al. (2011) of stretching the polar cells in θ and, additionally, use two to four layers of static mesh-refinement to decrease the ϕ resolution, giving a speedup by a factor of 4–16 (see SI Sections 2.2– 2.3). Local adaptive time-stepping between the tiles gives an additional speedup by a factor of 2–3. In practice, these advanced features give a speedup by a factor of 10–30 (see SI Section 2.1). The high speed of the code allows us to study tilted discs at much higher resolution and over longer durations than was possible until now, as required to handle the large dynamic range necessary to study tilted accretion and jets in 3D.

3 NUMERICAL MODELS

Our simulations start with a hydrostatic torus (Fishbone & Moncrief 1976) for a BH spin a = 0.9375 and use a Kerr–Schild foliation. We use an ideal gas equation of state, p_g = (γ − 1)u_g, where p_g and u_g are the thermal pressure and energy density, and γ = 5/3. We place the torus inner edge at r_in = 12.5r_g and density maximum ρ_max = 1 at r_max = 25r_g. We tilt the torus relative to the BH spin by an angle, |$\mathcal {T}_{\rm init} = 30^{\circ }$|⁠. In contrast to earlier work in which the BH spin was tilted with respect to the grid (Fragile & Anninos 2005; Fragile et al. 2007, 2009; McKinney et al. 2013), we tilt the disc itself and leave the BH spin pointed along the polar axis, because an axisymmetric metric lowers the memory footprint.

We carried out five production models listed in Table 1, each differing in the amount of the initial magnetic flux and resolution. For our strongly magnetized models (denoted with ‘S-’), we insert in the torus a large magnetic field loop, described by the magnetic vector potential A_ϕ ∝ (ρ − 0.05)²r³, where ρ is the rest mass density. For our weakly magnetized models (denoted with W-), we use a smaller loop, A_ϕ ∝ (ρ − 0.05). We simulate the two models, S-R and W-R, using sufficiently high resolutions (resolved, denoted with ‘-R’; see Table 1) to resolve the magnetorotational instability (MRI; Balbus & Hawley 1991) that fuels magnetic turbulence. We simulate models S-HR and W-HR at twice the resolution in all three dimensions to check for convergence (highly resolved, denoted with ‘-HR’). Finally, we set up a fifth model, W-U to be physically identical to model W-R but with only half the resolution in all three dimensions (underresolved, denoted with -U; see Table 1), in order to investigate the effect of underresolving the MRI.

Table 1.

(Top panel) Parameters common to all models: BH spin a, the radii of the torus inner edge r_in and pressure maximum r_max, the initial tilt angle |$\mathcal {T}_{\rm init}$|⁠, and simulation duration t_sim. (Bottom panel) Short and full model name (i.e. S25A93 stands for size and BH spin, W for weakly magnetized, LR or HR for low or high resolution), the strength of the magnetic flux, the resolution N_{r, θ, ϕ} and the quality factor Q_{r, θ, ϕ} (the number of cells per MRI fastest growing wavelength in r-, θ- and ϕ-directions) at t = 5 × 10⁴ t_g. Short names of the models link to 3D animations in a.

Model	a	rin (rg)	rmax (rg)	|$\mathcal {T}_{\rm init}\ ({\rm deg})$|
All	0.9375	12.5	25	30
Short	Full model	B flux	Resolution,	Q factor	tsim
name	name	strength	Nr × Nθ × Nφ	Qr, Qθ, Qφ	[105tg]
S-R	S25A93	Strong	448 × 144 × 240	(9,9,32)	1.2
S-HR	S25A93HR	Strong	896 × 288 × 480	(30,28,80)	0.75
W-U	S25A93WLR	Weak	448 × 144 × 240	(2.5,3,16)	1.2
W-R	S25A93W	Weak	896 × 288 × 480	(20,18,53)	1.2
W-HR	S25A93WHR	Weak	1792 × 576 × 960	(59,55,138)	0.3

We set the numerical resolution, N_r × N_θ × N_ϕ (see Table 1), to attain a near-unity cell aspect ratio everywhere. We normalize the magnetic pressure, p_B, such that it is subdominant compared to p_g, by setting p_{g, max}/p_{B, max} = 100. Following the approach of Ressler et al. (2017), we enforce throughout the simulations that ρc² ≥ max (p_B/50, 10⁻⁶c²(r/r_g)⁻²) and u_g ≥ max (p_B/150, 10⁻⁷c²(r/r_g)^−2γ), approximating physical processes that mass-load relativistic jets at their base, where p_B ≫ ρc².

We note that using a spherical polar grid to simulate tilted disc–jet systems is challenging, because improper treatment of the transmissive polar boundary condition can lead to numerical artefacts. To demonstrate that the disc and jets pass freely through the polar coordinate singularity, we carried out test simulations using the parameters of model S-R (see SI Section 2.4). First, we verified that the two equivalent ways of simulating a tilted system – i.e. separately misaligning (i) the BH (a = 0.9375) or (ii) the disc by 30° relative to the grid – do indeed give consistent profiles for |$\mathcal {T}(r)$| and |$\mathcal {P}(r)$|⁠. Secondly, we verified that the simulation outcome of a BH-torus system is independent of its orientation relative to the grid. Namely, we simulated an aligned BH-torus system oriented along the polar axis and compared it to that oriented at 30° relative to the polar axis, for two values of the BH spin (a = 0 and a = 0.9375). Please see SI Section 2.4 for more detail.

4 RESULTS

In all five models (Table 1), magnetic turbulence develops and the gas reaches the BH as seen in Figs 1(a) and (b). Accretion occurs via two polar plunging streams, consistent with the findings of Fragile et al. (2007), who studied the dynamics of tilted accretion flows with weak magnetic flux at a resolution similar to our model W-U. As Figs 1(c) and (d) show, each of our models produces an energetically significant outflow of energy, with power |$P_{\rm outflow}\equiv \dot{M}c^2-\dot{E} \gtrsim (0.1{-}1)\dot{M}c^2$|⁠. Here, |$\dot{M}$| and |$\dot{E}$| are mass and total energy accretion rates on to the BH. For all five models, 3D animations are available (see SI or this YouTube playlist). Fig. 2 shows a volume rendering for model W-R after 1.4 × 10⁴ t_g and 10⁵ t_g. Model S-R looks similar, but its jet has a larger opening angle reflecting the stronger outflow power than in model W-R (see Fig. 1b). In all our models, we find that the jets are connected by magnetic field lines to the event horizon (see Fig. 2), and therefore appear to be powered by the extraction of BH rotational energy via the Blandford–Znajek (BZ; Blandford & Znajek 1977) mechanism. They reach relativistic Lorentz factors γ ≳ 10 at r ≲ 200r_g, similar to 2D jets (Beskin & Nokhrina 2006; McKinney 2006; Tchekhovskoy, McKinney & Narayan 2008).

Figure 1.

Time dependence of mass accretion rate |$\dot{M}$| in units of M_init/t_g [panels (a) and (b)], where M_init is the initial torus mass; outflow power P_outflow in units of M_initc²/t_g [panels (c) and (d)]; precession angle |$\mathcal {P}$| [panels (e) and (f)] and tilt angle |$\mathcal {T}$| [panels (g) and (h)], where both angles measured over the radial interval 50–150 r_g; and disc radial extent r_disc [panels (i) and (j)], which is the average radius, in the plane of the disc, weighted by rest mass density. The left and right columns show, respectively, the weak (W-) and strong (S-) magnetic field models, for the resolved runs (-R, thick dashed blue lines), highly resolved runs (-HR, solid thin green lines), and underresolved run (W-U, dotted red lines). All variables show approximate convergence between the resolved and highly resolved runs. Both differ greatly from the underresolved run (W-U).

Figure 2.

Volume rendering of density [log (ρ) in blue and green, indicating the disc] and jet magnetic field lines coloured by the scaled magnetic energy density [p_Br/ρc² > 0.5 in red and yellow, indicating the jets] for model W-R at 1.4 × 10⁴ t_g (left) and 10⁵ t_g (right). The magnetic field lines in the jets are shown with yellow–red lines, threading the entire volume rendering of the jet. We measure distances in units of r_g. The disc–jet system precesses as a whole around the BH spin vector, which is vertical in the figure. As the simulation progresses, the disc spreads outward, its density profile flattens, and the distance between the high- and low-density regions increases.

Comparison of the two panels in Fig. 2 shows that the disc orientation changes over time, i.e. the tilted disc precesses. This is also seen in Figs 1(e) and (f) through the increase in |$\mathcal {P}$|⁠. We measure |$\mathcal {P}$| and |$\mathcal {T}$| of the disc by calculating its angular momentum vector (Nelson & Papaloizou 2000; Fragile & Anninos 2005) as a function r. [We measure the location of the jet by isolating the highly magnetized region (using p_B > ρc²/2r) and weighing its position by the magnetic pressure on spherical shells (see SI Section 3).] In addition to precession, the disc also aligns with the BH spin, as seen in Figs 1(g) and (h) through the decrease in |$\mathcal {T}$|⁠. In the simplest case, model W-U, the disc precesses at a constant rate, |${\rm d}\mathcal {P}/{\rm d}t={\rm constant}$|⁠, as expected for a misaligned disc angular momentum vector under the action of a constant LT torque (see Fragile et al. 2007). LT precession also induces a twist in the innermost disc. That effect remains relatively constant over time (in agreement with Fragile et al. 2007), allowing for solid-body-like precession of the whole system.

Figs 1(e)–(h) show that in our models with well-resolved MRI, S-R and W-R, the alignment and precession slow down over time. Indeed, well-resolved MRI is crucial for capturing disc angular momentum transport and outward expansion (see Figs 1i and j) that brings the disc angular momentum out of reach of jet and LT torques (both are stronger near the BH) and impedes the alignment and precession. The approximate agreement between the -R and -HR models (with numerical resolutions in every dimension different by a factor of 2) indicates that our results are reasonably converged with the numerical resolution.

Figs 2 and 3(a) show that the jets precess together with the disc: the spatially averaged (over 50r_g ≤ r ≤ 150r_g) |$\mathcal {P}$| of the disc and jet closely track each other in time. Fig. 3(b) shows that the disc and jet not only precess together but also align together, with the tilt angle decreasing from ∼30° to ∼25°.

Figure 3.

Time-dependence of precession [panel (a)] and tilt [panel (b)] angles of the disc (blue) and the upper jet (red), as measured at 50 r_g ≤ r ≤ 150 r_g in model W-R, with 1σ error bars showing variation across this radius interval. The jet direction closely tracks the disc rotation axis. Jet wobbling by several degrees, evident in both panels, could give rise to jet-powered high-energy flares.

Fig. 4 shows the tilt |$\mathcal {T}$| (top row) and precession |$\mathcal {P}$| (bottom row) angles for models W-R (left column) and S-R (right column) versus r at t = 5 × 10⁴ t_g. The jets align with the disc far from the BH. Radial tilt oscillations lead to the peak in disc tilt angle at r ∼ 10r_g, as expected for a thick disc (Lubow & Ogilvie 2000; Lubow et al. 2002), in agreement with previous work (Fragile et al. 2007). Strong twisting close to the BH, indicated by the steep radial dependence of disc precession angle at r ≲ 10r_g, builds up early on in the simulation and remains constant thereafter (Fragile et al. 2007).

Figure 4.

Radial profiles of tilt |$\mathcal {T}$| and precession |$\mathcal {P}$| angle of an upper jet (red dotted lines; the other jet is similar), averaged over t = (4.5–5.5) × 10⁴t_g, closely follow those of the disc (blue curves) at r ≳ 3r_g (using spherical radius for both the jet and the disc, error bars indicate 1σ variation across this time interval). In both cases of weak (model W-R, left-hand panels) and strong (model S-R, right-hand panels) magnetic flux, jets propagate along the disc rotational axis, not the BH spin axis. At r ≲ 10r_g, radial tilt oscillations cause variations in the disc and jet tilt angles [panels (b) and (d)]. If strong magnetic flux is present, as in model S-R [panel (d)], it can bend the innermost ≲3r_g of the disc and jets into partial alignment with the BH spin.

Figs 4(a) and (b) show that in model S-R the disc and jets tend to align with the BH spin in the inner few r_g, while in the weak field case W-R the jet remains misaligned up to the event horizon. This result supports the notion that strong jets can align the inner parts of disc–jet systems with the BH spin (McKinney et al. 2013; Polko & McKinney 2017). However, the alignment can be very limited in distance, r ≲ 3r_g, and magnitude: the stronger magnetic flux in model S-R leads to an additional disc alignment relative to model W-R by only ∼4° at r ∼ 15–50r_g. Also, the LT torque acting on a twisted accretion disc can contribute to the alignment indirectly. When precessed material mixes due to viscous dissipation (see e.g. King et al. 2005; Sorathia, Krolik & Hawley 2013), the net result is alignment with the BH spin. Since jets and (MRI-driven) dissipation both depend on magnetic field strength, both could lead to differences in alignment between S-R and W-R.

5 DISCUSSION AND CONCLUSIONS

We carried out GRMHD simulations of tilted BH accretion discs at the highest resolutions to date (effective resolution ∼1 billion cells). We find that our discs, with an initial tilt of 30° relative to the BH of spin a = 0.9375, undergo LT solid-body-like precession (see also Fragile et al. 2007). Our simulations for the first time show that tilted precessing discs can launch relativistic jets. The jets (i) propagate along the disc rotation axis and (ii) precess together with the disc.

The amount of large-scale vertical magnetic flux significantly affects the orientation of the disc and jets: when the BH is saturated with the flux, as in our model S-R, it is able to warp the disc and jets into partial alignment with the BH spin in the inner few r_g (Fig. 4b). Differences in magnetic flux content can also explain differences with McKinney et al. (2013), who reported partial alignment of the jet with the BH spin up to ∼100 r_g. We have replicated in h-amr one of their models (A0.99N100T0.6). Their initial conditions contain much more magnetic flux in a 100-fold larger disc. As accretion drags the large-scale vertical magnetic flux inwards over time, not only the BH but also the inner accretion disc gets saturated with magnetic flux (i.e. reaches the MAD state), with average p_g/p_mag within 100 r_g decreasing to ∼1, much lower than 18 and 35 in our models S-R and W-R, respectively.

Type-C QPOs observed from BH XRBs have been interpreted as LT precession of the inner accretion flow (Stella & Vietri 1998; Ingram et al. 2009, 2016). If the jet base is X-ray bright (Markoff, Nowak & Wilms 2005), the precession can lead to quasi-periodic swings in the X-ray polarization angle in addition to those expected from the inner disc alone (Ingram et al. 2015). This can be tested with the Imaging X-ray Polarimetry Explorer due to launch in 2020. Observing longer term changes in large-scale jet orientation due to LT precession (e.g. Ekers et al. 1978; Bridle et al. 1979; Martí-Vidal et al. 2011; Kalamkar et al. 2016) can enable new tests of GR, BH accretion and jet physics (e.g. Stone & Loeb 2012). The precession of AGN jets can cause them to spread their power over a large area and heat the ambient gas instead of escaping out of the galaxy/cluster (Nawaz et al. 2016; Yang & Reynolds 2016). Together with the jet wobbling due to magnetic instabilities (Tchekhovskoy & Bromberg 2016), this can provide an explanation for the unexpectedly high temperature of the intracluster medium, known as the cooling flow problem.

In addition to solid-body-like precession, we find that tilted disc–jet systems align with the BH spin axis over longer, accretion time-scales (see also Foucart et al. 2011). Strong large-scale magnetic fields and the associated powerful jets accelerate this alignment (Figs 4a and b). Short-time-scale wobbles of several degrees in jet orientation superimposed on the smooth precession and alignment trends (Figs 3a and b) could boost jet emission in and out of our line of sight and result in jet-powered high-energy flares similar to those in the X-ray light curve of a TDE Swift J1644 (e.g. Bloom et al. 2011).

We found it important for the simulations to resolve the MRI throughout the disc (see also this 3D visualization of the effects of resolution on our weak-field models). Magnetic turbulence causes angular momentum transport and radial disc expansion (Figs 1i and j) that substantially slows down the precession (Figs 1e and f) and alignment (Figs 1g and h) by redistributing most of the angular momentum out of reach of the BH. Our thick discs are isolated, and external disc feeding from large radii (e.g. by a stellar debris stream in TDEs or a geometrically thin disc in AGNs and XRBs) can affect the disc expansion either by the addition of new gas or by applying external pressure. Also, geometrically thin discs (H/R < α) have larger viscosity, which can reduce disc expansion and affect precession and alignment. We use h-amr to study that limit in future work.

Acknowledgements

We thank Chris Fragile, Nick Stone and Asaf Pe'er for helpful comments and Mark Vanmoer for his help with visualization. This research was made possible by NSF PRAC award no. 1615281 at the Blue Waters sustained-petascale computing project and supported in part under grant no. NSF PHY-1125915. ML and MK were supported by the Netherlands Organisation for Scientific Research (NWO) Spinoza Prize, CH by the Amsterdam Science Talent Scholarship, AI by the NWO VENI grant (no. 639.041.437), AT by the TAC and NASA Einstein (grant no. PF3-140131) postdoctoral fellowships, and SM by the NWO VICI grant (no. 639.043.513).

SUPPORTING INFORMATION

Supplementary data are available at MNRASL online.

Simulation data is available upon request from AT at [e-mail]

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.