Black hole magnetosphere with small scale flux tubes--II. Stability and dynamics

TL;DR

This study uses 3D force-free simulations to explore how small-scale magnetic flux tubes near black holes influence X-ray emission, revealing that confinement prevents jet formation and leads to localized energy dissipation.

Contribution

It demonstrates how magnetic flux tube dynamics under confinement affect energy release and jet formation near black holes, providing insights into X-ray emission mechanisms.

Findings

Flux tubes become kink unstable and dissipate energy close to the black hole.

Strong confinement prevents jet formation despite flux tube inflation.

Energy dissipation may heat plasma and produce X-ray emission.

Abstract

In some Seyfert Galaxies, the hard X-rays that produce fluorescent emission lines are thought to be generated in a hot corona that is compact and located at only a few gravitational radii above the supermassive black hole. We consider the possibility that this X-ray source may be powered by small scale magnetic flux tubes attached to the accretion disk near the black hole. We use three dimensional, time dependent force-free simulations in a simplified setting to study the dynamics of such flux tubes as they get continuously twisted by the central compact star/black hole. We find that, the dynamical evolution of the flux tubes connecting the central compact object and the accretion disk is strongly influenced by the confinement of the surrounding field. Although differential rotation between the central object and the disk tends to inflate the flux tubes, strong confinement from surrounding field quenches the formation of a jet-like outflow, as the inflated flux tube becomes kink unstable and dissipates most of the extracted rotational energy relatively close to the central object. Such a process may be able to heat up the plasma and produce strong X-ray emission. We estimate the energy dissipation rate and discuss its astrophysical implications.

Equations

$${\mathbf{B}}_H\equiv 4\pi {\mathcal{J}}_H\times \mathbf{n},$$

$${\mathbf{E}}_H\equiv {\alpha }_H{\mathbf{E}}_{\parallel }=R_H{\mathcal{J}}_H,$$

$$R_H\equiv 4\pi \approx 377\text{ohms}.$$

$${\mathbf{B}}_H={\mathbf{E}}_H\times \mathbf{n}.$$

$${\sigma }_H\equiv E_n/(4\pi ),$$

$$\partial {\sigma }_H/\partial t{+}^{(2)}\nabla \cdot {\mathcal{J}}_H+({\alpha }_H{j}_n)=0,$$

$${\mathcal{J}}^{\prime }=\mathbf{n}\times {\mathbf{B}}_{\parallel }^{\prime }/(4\pi ).$$

$${\mathbf{E}}_{\parallel }^{\prime }=R{\mathcal{J}}^{\prime },$$

$${\mathbf{B}}_{\parallel }^{\prime }=\frac{1}{\eta }{\mathbf{E}}_{\parallel }^{\prime }\times \mathbf{n}.$$

$${\mathbf{E}}_{\parallel }^{\prime }=\gamma ({\mathbf{E}}_{\parallel }+\mathbf{v}\times {\mathbf{B}}_{\perp })-(\gamma -1)({\mathbf{E}}_{\parallel }\cdot \hat{v})\hat{v},$$

$${\mathbf{B}}_{\parallel }^{\prime }=\gamma ({\mathbf{B}}_{\parallel }-\mathbf{v}\times {\mathbf{E}}_{\perp })-(\gamma -1)({\mathbf{B}}_{\parallel }\cdot \hat{v})\hat{v},$$

$$\frac{\partial \mathbf{E}}{\partial t}$$

$$\frac{\partial \mathbf{B}}{\partial t}$$

$$\mathbf{J}$$

$$E_r$$

$$E_{\varphi }$$

$$E_r=-\frac{{\Omega }_{d}r}{c}{B}_z.$$

$$\frac{\partial \mathbf{E}}{\partial t}$$

$$\frac{\partial \mathbf{B}}{\partial t}$$

$$\lambda (\mathbf{r})={\lambda }^{\ast }(\mathbf{r})=\frac{K_{\mathrm{abs}}}{\Delta t}{\left(\frac{d}{d_{\max }}\right)}^3.$$

$$i_{\varphi }=\{\begin{array}{ll}{i}_0\cos \left(\frac{2\pi }{r_0}(r-r_2)\right){\left(\frac{r_2}{r}\right)}^{\alpha },& r_2\le r\le r_2+\frac{3}{4}{r}_0,\\ 0,& \text{otherwise}.\end{array}$$

$$t_{\mathrm{kink}}^{\prime }=\frac{r}{v_A},$$

$$t_{\mathrm{kink}}=\gamma \frac{r}{v_A},$$

$$z=k(r-r_{\mathrm{LC}}),$$

$$\gamma \approx r/r_{\mathrm{LC}}.$$

$$t_{\mathrm{kink}}\approx \frac{r^2}{r_{\mathrm{LC}}{v}_A}\approx \frac{r^2}{r_{\mathrm{LC}}c}.$$

$$t_{\mathrm{flow}}\approx \frac{z}{v\cos \theta }\approx \frac{\sqrt{z^2+(r-r_{\mathrm{LC}}{)}^2}}{c}.$$

$$\frac{t_{\mathrm{kink}}}{t_{\mathrm{flow}}}\approx \frac{(h/k+1{)}^2}{h\sqrt{1+1/k^2}}.$$

$$\mathbf{E}=-(\mathbf{\omega }\times \mathbf{r})\times \mathbf{B}/c,$$

$$(\mathbf{E}\times \mathbf{B}{)}_{\varphi }=\frac{\omega r\sin \theta }{c}{B}_p^2,$$

$$\omega =\frac{(\mathbf{E}\times \mathbf{B}{)}_{\varphi }}{B_p^{2}r\sin \theta }.$$

$$E_r=-\eta \frac{B_{\varphi }}{\gamma }-\frac{\Omega r}{c}{B}_z.$$

$$E_r=-\frac{\omega r}{c}{B}_z,$$

Full text

††pubyear: 2019††pagerange: Black hole magnetosphere with small scale flux tubes–II. Stability and dynamics–B

1 Introduction

In some Seyfert Galaxies there is a strong yet variable continuum X-ray source that irradiates the accretion disk and produces fluorescent emission lines. It is usually believed that the primary X-rays come from a hot corona located on either side of the accretion disk, where energetic electrons up-scatter optical/UV emission from the disk into X-ray range. Typical X-ray luminosity from the corona can be comparable to the disk thermal luminosity, indicating quite high X-ray efficiency. Recent observations have made good progress in constraining the geometry of the corona. Some pieces of evidence, including X-ray reverberation mapping (e.g., Uttley et al., 2014; Kara et al., 2016), disk emissivity profile modeling (Wilkins & Fabian, 2011; Wilkins & Gallo, 2015) and micro-lensing (Morgan et al., 2008; Chartas et al., 2009; Mosquera et al., 2013; Reis & Miller, 2013), suggest that the corona is quite compact and located at a few ($\lesssim 20$) gravitational radii above the black hole, which seems consistent with the so-called “lamppost” geometry. However, the formation and heating mechanism of the corona is still unknown.

Another hint is that these Seyfert Galaxies usually do not create powerful jets. Meanwhile, their accretion rates are in the intermediate regime (the luminosity $L$ is in the range $0.01L_{\rm Edd}\lesssim L\lesssim L_{\rm Edd}$), so the accretion disk is likely to be geometrically thin (e.g., Yuan & Narayan, 2014). Given these circumstances, we consider the possibility that the X-ray corona may be powered by small scale magnetic flux tubes on the accretion disk around the black hole (Yuan et al., 2019, henceforth Paper I). These flux tubes may have foot points located at different radii on the disk or one foot point on the disk and the other on the black hole. This kind of field configurations may arise if the large scale net flux is weak so that the inner disk is dominated by small coherence length scale structures produced by magneto-rotational instability (MRI) in the disk; or, similar to the solar surface, flux tubes may emerge from hot spots on the disk due to magnetic buoyancy. Differential rotation causes the flux tubes to twist and inflate (e.g., Achterberg, 1996a, b; Uzdensky & Goodman, 2008); in particular, they may get tangled up on or near the axis. This may lead to strong dissipation as the field relaxes/untangles, and most of the extracted rotational energy is converted into particle kinetic energy, then radiated away as X-rays. In such a way, the black hole may directly power the X-ray corona, instead of launching a powerful jet.

In Paper I, we set up simplified models where the magnetic field attached to the disk and the black hole is still axisymmetric, but has relatively small coherence length scales along the radial direction—on the order of a few gravitational radii (similar to Parfrey et al., 2015; Uzdensky, 2005). We found that, under the axisymmetric constraint, there exist steady configurations with disk-hole linking “closed” flux tubes, consistent with Uzdensky (2005). Furthermore, the extent of the closed zone is determined by the pressure balance between the spin induced twist in the closed zone and the poloidal confinement from both the closed zone itself and the overlying external field from the disk. The maximum extent decreases with black hole spin; for a typical confinement and black hole spin, the closed zone does not exceed a height of tens of gravitational radii. In a real dynamic situation, the equilibrium may be unstable, or the pressure balance may not be easily maintained. The flux tubes may become highly dynamic and produce significant dissipation. It is the goal of this paper to follow the time-dependent evolution of the flux tubes.

The paper is organized as follows. We describe our simplified setup and numerical method in §2. The simulation results are presented in §3 and §4. We then discuss the implications and conclusions in §5.

2 Numerical approach

2.1 The Membrane paradigm

The boundary conditions for the electromagnetic fields at the event horizon of a Kerr black hole can be interpreted in a way such that the event horizon behaves like a resistive membrane (e.g., Thorne et al., 1986). In particular, define the tangential magnetic field at the horizon as $\mathbf{B}_{H}=\alpha_{H}\mathbf{B}_{\parallel}$, where $\mathbf{B}_{\parallel}$ is the tangential magnetic field as measured by a fiducial observer (FIDO, taken to be ZAMO) close to the event horizon, and $\alpha$ is the lapse function there (though the magnitude of $\mathbf{B}_{\parallel}$ diverges at the event horizon, $|\alpha_{H}\mathbf{B}_{\parallel}|$ remains finite), then one can define a surface current $\boldsymbol{\mathcal{J}}_{H}$ from Ampere’s law by requiring that the tangential magnetic field terminates at the event horizon:

[TABLE]

where $\mathbf{n}$ is a unit vector normal to the event horizon. The surface current is driven by the tangential component of the horizon’s renormalized electric field and satisfies an Ohm’s law

[TABLE]

where the surface resistivity

[TABLE]

The Ohm’s law and Ampere’s law are equivalent to the radiative boundary condition at the event horizon

[TABLE]

If one also defines the surface charge density $\sigma_{H}$ by Gauss’s law,

[TABLE]

where $E_{n}$ is the normal component of the electric field as measured by a FIDO at the event horizon, then one can obtain the charge conservation relation

[TABLE]

where ${}^{(2)}\nabla\cdot\boldsymbol{\mathcal{J}}_{H}$ is the 2-dimensional divergence of the 2-dimensional surface current; $j_{n}$ is the normal component of the volume current leaving the horizon as measured by a FIDO.

In our flat spacetime toy model, we take into account the electromagnetic property of the event horizon using a similar membrane formulation. Suppose the event horizon can be described as a resistive membrane, rotating at an angular velocity $\boldsymbol{\Omega}$. A FIDO that is sitting on the membrane rotates with the membrane at the same velocity $\mathbf{v}=\boldsymbol{\Omega}\times\mathbf{r}$ and sees a tangential magnetic field $\mathbf{B}_{\parallel}^{\prime}$ terminating at the membrane (we use prime to denote quantities measured by the rotating FIDO), which requires a surface current density

[TABLE]

The surface current needs to be driven by a tangential electric field

[TABLE]

where $R$ is the surface resistivity. To mimic a black hole we take $R=4\pi$. For a general Ohmic resistivity $R$, the above relations give (in what follows we write $R/4\pi\equiv\eta$)

[TABLE]

The FIDO-measured fields are related to the lab frame values through a Lorentz transformation

[TABLE]

where $\gamma=1/\sqrt{1-v^{2}/c^{2}}$ is the Lorentz factor of the FIDO.

In Appendix A we show a comparison between steady state configurations obtained using the above membrane boundary conditions in flat spacetime and those for exact Kerr spacetime. The similarity lends further support that the membrane analogy captures the electromagnetic properties of the black hole reasonably well.

2.2 Simulation setup

We consider the force-free limit of magnetohydrodynamics, which is appropriate when the magnetization is sufficiently high. In flat spacetime, the system of equations can be written as (e.g., Gruzinov, 1999; Blandford, 2002)

[TABLE]

We use the three dimensional, Cartesian, special relativistic force-free code developed by Spitkovsky (2006) to follow the evolution of the system. To mimic the black hole, we put a disk-like membrane with radius $r_{1}$ and angular velocity $\Omega$ at the center of the $z=0$ plane, on which we apply the boundary condition described in §2.1. More specifically, we specify $\boldsymbol{E}_{\parallel}$ using equations (9–11):

[TABLE]

Here $r$ denotes the radial direction on the disk-like membrane, and $\phi$ indicates the azimuthal direction. The “black hole” (or rather, “black disk”) is surrounded by a thin accretion disk, modeled by a perfectly conducting surface extending over the region $r>r_{2}$ on the $z=0$ plane. The inner boundary of the disk $r_{2}$ may or may not correspond to the innermost stable circular orbit (ISCO). In what follows we show two classes of examples: one has $r_{2}\gtrapprox r_{1}$ to mimic a rapidly spinning black hole with a disk extending close to the event horizon; the other has $r_{2}=2r_{1}$ to represent a disk further away from the hole (e.g., truncated disk). In both cases the angular velocity of the central compact object satisfies $\Omega r_{1}/c=0.9$. We have experimented with different $\Omega$ values and found the results qualitatively similar. Suppose the angular velocity of the disk is $\Omega_{d}(r)$, then the boundary condition on the accretion disk is

[TABLE]

To minimize the effect of Cartesian grid stair-stepping at the boundary of the disk and the membrane, we apply smoothing of the fields over a few grid cells111Near the edge of the membrane, we apply the boundary condition as $E_{\parallel}=sE_{\parallel,\rm membrane}+(1-s)E_{\parallel,\rm local}$, where $s=[1-\tanh((r-r_{1})/\delta)]/2$ and the smoothing length scale $\delta$ is 1 to 2 grid cells.. For the outer boundary condition of the simulation domain, we use an absorbing layer that damps all electromagnetic disturbances. The absorbing layer satisfies the modified Maxwell equations that contain an electric and magnetic conductivity term (e.g., Cerutti et al., 2015)

[TABLE]

where the conductivity varies with distance $d$ to the boundary of the physical domain

[TABLE]

We start with a vacuum potential magnetic field generated by a surface current distribution $i_{\phi}$ in the disk, then turn on the rotation of the “black hole” and the accretion disk, and watch how the field evolves. In what follows, for simplicity, we set the accretion disk to be non-rotating: $\Omega_{d}(r)=0$, which is a rough approximation when the central object is rotating with a rate $\Omega$ much larger than every radius of the disk.222For a real black hole, if the inner boundary of the disk is at the ISCO, then the black hole rotates faster than every radius of the disk when its spin $a>0.3594$, or equivalently, the angular velocity of the event horizon $\Omega_{H}>0.093c/r_{g}$. We will also briefly discuss the effects of the disk rotation in section 3.4.

The initial disk current distribution is

[TABLE]

The current distribution has two oppositely directed components (see the first row of Figure 1); $r_{0}$ and the exponent $\alpha$ control the separation and relative strength of the two components, respectively. The resulting field geometry has two dipole-like loops (the second row of Figure 1). As $\alpha$ increases, the outer current loop becomes weaker, so does the outer magnetic field. As a result, the inner magnetic field loop becomes larger in size.

In our toy model the magnetic flux is fixed in the accretion disk; this could be a reasonable assumption if the overall flux changes on the accretion time scale, which is much longer than the orbital time scale relevant for our simulation here.

3 Simulation results: the dynamics

3.1 The case of $\eta=0$: a perfectly conducting star

When we set $\eta=0$, equation (15) and (16) indicate that the central membrane behaves like a perfectly conducting star. Any field line that threads the star will have its foot point corotating with the star. If the other end is attached to the disk, differential rotation between the disk and the star will build up more and more twist along the field line. Normally under axisymmetric constraint, the magnetic field will be inflated by the increasing toroidal twist, and eventually open up to infinity (e.g., Uzdensky et al., 2002). This is a direct consequence of the relativistic Ferraro’s Law of isorotation (e.g., Blandford & Znajek, 1977): if a steady state can be reached in axisymmetry, the angular velocity along the field line must be a constant. However, in a real 3D case, the situation can be different.

3.1.1 Different behaviors as $\alpha$ changes

In the first series of experiments, we fix the inner boundary of the ‘accretion disk’ at $r_{2}$ and the current loop separation at $r_{0}$, but change the exponent $\alpha$ to vary the relative strength of the two current loops. We find qualitatively different behaviors when $\alpha$ changes.

For sufficiently large $\alpha$, the outer current loop is relatively weak compared to the inner one. As we turn on the rotation of the star, a torsional Alfvén wave is emitted from the star and propagates out along the field lines threading the star. However, the other end of the field line is attached to the disk and is nonrotating. Therefore these field lines get twisted up, building up a magnetic tower (e.g., Lynden-Bell, 1996) which eventually breaks out from the overlying field. The inner field loop becomes open as a result, and there is a current layer separating the two polarities of magnetic field (Figure 2). This current layer marks the ‘jet wall’: within the jet wall, the field is rotating with the angular velocity of the ‘star’ and there is a Poynting flux going out to infinity; outside the jet wall, the field is attached to the accretion disk and is non-rotating.

When we decrease $\alpha$, the opening angle of the ‘jet’ decreases, as can be seen clearly in the middle column of Figure 2. For small enough opening angle, we start to see the development of kink instability in the jet. The kink mode growth rate increases as the opening angle decreases. For very small $\alpha$, namely, when the outer current loop is much stronger than the inner one, we do not see the jet break out: the kink mode develops very early and completely destroys the outflow. Figure 3 shows 3D renderings of field lines at three different time slices for the simulation with $\alpha=0.2$, where the inner flux loop first opens up, then gets destroyed by the kink instability. When the kink is fully developed, we see that the field lines from the star connect to one side of the accretion disk, forming a closed structure. As the star rotates, the closed structure precesses—its foot points on the disk appear to rotate around the star, in an opposite direction, with an apparent angular velocity that is about half of the angular velocity of the star. The field lines in and around this closed structure are actually continuously changing partners, or equivalently, reconnecting: if a particular point on the star and a particular point on the disk is connected by a field line at some point, they can only maintain the connectivity for a short amount of time; the rotation of the star will increase the toroidal twist in this connecting field line, which eventually breaks—the disk-end becomes an open field line, while the star-end gets connected to an open field line from the opposite side of the disk. The star-end remains closed, but it is now connected to a different point on the disk. As a result, the field lines on the accretion disk around the closed zone alternate between closed and open configurations, with continuous reconnection happening at the separatrix layer that lies between the open and closed field lines. A movie corresponding to this simulation is available online333Direct link: https://youtu.be/htSQFD64JGI.

3.1.2 The criterion for kink instability

In these configurations, when the field lines attached to the central membrane open up, the toroidal field becomes dominant toward the edge of the jet funnel (on the boundary the poloidal field reverses and goes through zero but the toroidal field is nonzero), as shown in Figure 4. Jet-like outflows dominated by toroidal field are prone to kink instability, which grows on Alfvén wave crossing time scales in the rest frame of the flow (e.g. Lyubarskij, 1992; Appl et al., 2000; Porth & Komissarov, 2015; Bromberg & Tchekhovskoy, 2016). If in the lab frame the kink growth time scale is shorter than the flow time scale, the jet can become kink unstable. Equivalently speaking, when the jet material stays in causal contact on the flow time scale, kink will be able to develop. In the above examples, we see that for different $\alpha$ the jet opening angle under axisymmetry constraint is different, and this affects the growth of the kink mode. We can obtain a rough criterion for the kink instability as follows.

The kink growth time scale in the rest frame of the flow is approximately

[TABLE]

where $r$ is the radius of the jet and $v_{A}$ is the Alfvén velocity. In the lab frame, we have

[TABLE]

where $\gamma$ is the Lorentz factor of the flow. From the simulation results we see that the jet boundary can be described approximately as

[TABLE]

where $r_{\rm LC}$ is the radius of the light cylinder of the star, and $k$ is a constant characterizing the opening angle of the jet. The Lorentz factor of the flow can be estimated as (e.g., Lyubarskij, 1992)

[TABLE]

This is exact for Michel monopole solution (Michel, 1973), and is also a good approximation for our case, as can be seen in Figure 5. Now the kink growth time scale is

[TABLE]

On the other hand, the flow time scale is

[TABLE]

Writing $z=hr_{\rm LC}$, we get the ratio

[TABLE]

Figure 6 shows the ratio $t_{\rm kink}/t_{\rm flow}$ as a function of $h$ for different opening angles. When $k$ gets larger than 4, namely the opening angle becomes smaller than $\sim 14^{\circ}$, $t_{\rm kink}/t_{\rm flow}$ can become less than 1 and the kink can develop. Our $\alpha=0.4$ case is close to marginal instability and agrees well with the above criterion.

3.2 Effect of changing the membrane resistivity $\eta$

Changing the membrane resistivity $\eta$ can also change the behavior of the system qualitatively. In Figure 7, we show a series of experiments where the current distribution in the disk is kept the same, while the membrane resistivity varies. We still set $r_{2}=1.08r_{1}$, and the exponent $\alpha$ is 0.9. For zero resistivity case, the inner field loop quickly opens up, and remains a stable steady state thereafter. The configuration is axisymmetric, and we can directly calculate the angular velocity $\omega$ of the field line in the following way: since

[TABLE]

the $\phi$ component of $\boldsymbol{E}\times\boldsymbol{B}$ is thus

[TABLE]

where the subscript p denotes poloidal component. So we have

[TABLE]

The lower panels of Figure 7 show the angular velocity $\omega$ of the field line. Indeed, when $\eta=0$, $\omega=\Omega$. As we increase the resistivity of the membrane, the angular velocity of the field line decreases. This is because the field lines can now slip on the membrane, and the resistivity allows them to rotate at a slower rate than the membrane itself. The toroidal field also decreases, meaning that the field lines are not twisted up as efficiently.

Close to the center of our disk-like membrane, the magnetic field is nearly uniform and points along z direction. We can use the membrane boundary condition to determine the angular velocity of these field lines. Equation (15) suggests that on the membrane, we have

[TABLE]

Right outside the membrane, the plasma is perfectly conducting so we have

[TABLE]

and the outgoing Alfvén wave along the vertical uniform field has

[TABLE]

Using the three equations, we obtain the angular velocity of the field line to be

[TABLE]

when $\Omega r\ll c$, $\gamma\approx 1$. In the simulations as we increase the resistivity $\eta$, the field line angular velocity decreases, agreeing with this analytical expression.

For sufficiently large resistivity, the inner magnetic flux loop no longer opens up. At $t=0$ when we turn on the rotation of the membrane, the membrane boundary condition still results in a torsional Alfvén wave being emitted from the foot points of the field lines on the black hole. The Alfvén wave propagates along the flux tubes, reaches the accretion disk, and gets reflected. The interaction between the outgoing wave and the reflected wave tries to establish an equilibrium configuration where the field lines linking the membrane and the accretion disk all have the same angular velocity as their foot points on the disk, namely, $\omega=\Omega_{d}=0$. Despite their zero angular velocity, the field lines do develop a toroidal component as they are dragged by the rotation of the membrane and have to slide on the membrane. From Figure 7 we can see that the vertical extent of the membrane-disk link decreases as the resistivity increases, because there is less twist in the magnetic field, or equivalently, less toroidal pressure pushing on the external field, when resistivity increases. This asymptotes to a vacuum potential field when we set $\eta$ to be very large. The existence and behavior of the closed zone with disk-hole linking magnetic flux agree with our results in Paper I.

We can also view this from a slightly different angle. The closed field lines lie in the sub-Alfvénic regime, as such an equilibrium can only be established when signals—Alfvén waves—can propagate back and forth along the field lines. Indeed, the initial Alfvén pulse has to be reflected at some point (e.g., the foot point on the disk) in cases where a closed equilibrium develops. We find from our simulations that the necessary condition for this is that the reflection point lies within the effective light surface as determined by the angular velocity from Equation (35) if the field lines were to open up. In Figure 7, as the resistivity $\eta$ increases, the effective light surface is moved to larger radii, and when it is further away than the reflection point of the closed flux bundle, the field lines no longer open up.

Around the transition between the open final state and the closed final state, the system is prone to the kink instability, as shown by the middle column in Figure 7. Although the opening angle of the outflow does not seem to change with the resistivity $\eta$, increasing $\eta$ does reduce the outflow velocity and puts the light surface at larger radius. When a large portion of the jet base lies within the light surfaces, namely within causal contact, kink instability will start to develop. Figure 8 shows another set of examples where we put the disk current further away from the central membrane, with $r_{2}=2r_{1}$ and $\alpha=1.0$. The configuration is prone to kink in all different resistivity cases, but the kink grows fastest near the open/closed transition.

3.3 The case of a “black hole” ($\eta=1$)

In what follows, we consider in more detail the case corresponding to a black hole: $\eta\equiv R/4\pi=1$. Depending on the current distribution in the disk, there are a few different regimes.

When the inner current loop is sufficiently weak compared to the outer one (small $\alpha$ in our setup), we get a final state that is an axisymmetric, disk-hole linking, closed configuration, with zero angular velocity (since $\Omega_{d}=0$) as we have seen in the last column of Figure 7. This is very similar to what we found in Paper I: the closed region remains sub-Alfvénic and appears to be stable.

On the other hand, when the inner current loop is too strong (large $\alpha$ in our setup), the resulting zero-angular-velocity configuration will develop a toroidal magnetic field that is too strong to be confined by the pressure provided by the field from the outer disk. The inner closed loop will be able to break out and open up eventually, reaching an angular velocity that is about half the angular velocity of the black hole, as required by the membrane boundary condition.

In the intermediate $\alpha$ range, we see the development of the $m=1$ kink mode. Figure 9 shows a few snapshots of the kink mode from a particular simulation where the disk current distribution has a positive loop $j_{1}$ located at $r=2.0r_{1}$ and a negative loop $j_{2}$ located at $r=5.33r_{1}$, with $j_{2}=0.7j_{1}$. The development of the kink is very similar to the case corresponding to the first row of Figure 2.

3.4 Effects of the accretion disk rotation

So far we have set the accretion disk at rest, with $\Omega_{d}(r)=0$. In reality, the accretion disk is orbiting around the black hole with Keplerian angular velocity, and the shear motion could affect the evolution of the magnetic field on the accretion disk. In figure 10 we show such an example. We take the case corresponding to the last row in Figure 2, and now allow the accretion disk to rotate with angular velocity $\Omega_{d}(r)=f\Omega(r/r_{2})^{-3/2}$, where we set $f<1$ to model a disk rotating slower than the black hole. Different $f$ values show qualitatively similar results. The rotation of the disk causes the field loop on the accretion disk to open up, launching Poynting flux/magnetized wind that flows out from every radius of the disk. If an axisymmetric constraint were imposed, at large radii the field would asymptote to a monopole configuration with alternating polarities. The interesting fact here is that the disk wind now imposes more pressure on the central “jet” column: the opening angle of the column is reduced and it is easier for the kink instability to set in. This again lends support to our basic conclusion that strong confinement from the disk field/disk wind makes the jet more prone to kink instability when the field coherence length scale is small.

4 On the energy extraction and dissipation

We can directly measure the energy extraction from the central rotating membrane by integrating the Poynting flux over the surface of the membrane. Figure 11 shows four examples of power output as a function of time. The total power is normalized to the fiducial value $P_{0}=\Omega^{2}\Phi^{2}/c$, where $\Phi$ is the magnetic flux threading the membrane at $t=0$. We can see that the power output from the membrane is indeed on the order of $P_{0}$. The actual prefactor depends on the resistivity of the membrane and the way the magnetic flux is distributed on the membrane. Blue line and red line both correspond to final steady, open configurations, but when the membrane has significant resistivity (red line), the power output becomes much smaller. For the case of the yellow line ($\eta=0$), the increase around $t/T=10$ corresponds to the fully developed kink mode: although the field lines are fixed to the perfectly conducting central object and rotate at the same angular velocity, the kink mode changes the current distribution on the field lines, thus changing the Poynting flux as well. In the case of the cyan line (black hole, $\eta=1$), around $t/T\approx 20$ the kink instability has fully developed; we also see a slight increase in the total power output from the membrane.444For a stable, closed configuration similar to the last column of Figure 7, based on our results from Paper I, the energy extraction rate $P\propto\omega I$, where $\omega$ is the angular velocity of the field line and $I$ is the poloidal current that flows along the field lines threading the black hole, we should have $P=0$ as $\omega=0$. In reality we see small amounts of energy extraction, part of which gets dissipated volumetrically in the closed region, while the majority gets dissipated along the field lines connecting to the edge of the membrane where a boundary current layer lies.

The energy extracted from the central membrane may have three outcomes: (1) flowing along the open field lines toward infinity; (2) flowing along closed field lines to the disk; (3) dissipated in the magnetosphere. Since our accretion disk is perfectly conducting and non-rotating, which acts like a reflecting wall, there is no Poynting flux entering the disk, namely, component (2) should be zero. Component (1) becomes negligible when most of the field lines are closed, or open flux bundle gets destroyed by the kink mode (dotted lines in Figure 11). As a result, in these cases most of the extracted rotational energy gets dissipated in the magnetosphere. In what follows we would like to get a sense of where the dissipation tends to take place in the magnetosphere. When force-free condition is satisfied, dissipation does not happen because $\boldsymbol{E}\cdot\boldsymbol{J}=0$ is ensured by the force-free condition. However, as we have seen, the system may evolve toward a state where the force-free condition gets violated in certain regions, and dissipation has to happen there. In our current simulation scheme, the dissipation is numerical in origin. Force-free gets violated when $\boldsymbol{E}\cdot\boldsymbol{B}\neq 0$ or $E>B$: we then directly reduce the corresponding component of $\boldsymbol{E}$ to restore it back to force-free. This can also be regarded as efficient radiation that removes the electromagnetic energy from the system.

Based on the energy conservation law, we have

[TABLE]

Writing the energy density of the electromagnetic field as $\mathcal{E}=(B^{2}+E^{2})/(8\pi)$, and the Poynting vector as $\boldsymbol{P}=\boldsymbol{E}\times\boldsymbol{B}c/(4\pi)$, we get

[TABLE]

For ideal force-free, $\dot{\mathcal{E}}_{\rm diss}=0$, but it will be nonzero in regions where force-free condition is violated. We can directly calculate $\dot{\mathcal{E}}_{\rm diss}$ using the left hand side of Equation (37) to see where the dissipation happens.

Figure 12 shows an example for the case where the membrane is perfectly conducting and the system eventually reaches a steady, axisymmetric open configuration, corresponding to the last row of Figure 2. In this case, $\dot{\mathcal{E}}_{\rm diss}=\nabla\cdot\boldsymbol{P}$. By comparing the plot of $\dot{\mathcal{E}}_{\rm diss}$ with the plot of current density, we see that the current sheet at the jet wall is the main site for dissipation.

Figure 13 shows another example where the evolution is highly dynamic and asymmetric. This case corresponds to the second row of Figure 2 where the kink instability fully destroys the outflow. We can see again that the dissipation sites coincide with places where the current is strong. In the first column, a transient state with field lines attached to the central star fully open, there is a current layer separating the field lines from the star and those from the disk, which have different polarities. The current distribution and dissipation are similar to those in Figure 12. When the kink has fully developed, we can see a strong current layer wrapping around the closed flux bundle. This is where the separatrix lies, and where the field lines reconnect and alternate between closed/open configurations. Indeed these appear to be strong dissipation sites. However, in order to accurately account for the energy dissipation rate, we need to include more physical dissipation mechanisms. This will be studied in a future work.

5 Discussion and Conclusions

In this paper, we use a simple setup in special relativistic force-free framework to study the dynamics of small scale flux tubes connecting the central black hole/star and the surrounding thin accretion disk. For a perfectly conducting star, the field that anchors on the star will always extract significant amount of energy from the star, but whether an outflow can be formed depends on the competition between the pressure of the twisted field in the magnetic tower and the external field that confines it. Strong confinement from the surrounding field could suppress the formation of an outflow and lead to the dissipation of most of the extracted energy near the central star. Weak confinement, on the other hand, facilitates the formation of a relatively stable outflow of Poynting flux with a wide opening angle as the magnetic twist builds up. For a black hole, since it behaves like a resistive membrane and the field lines can slip on it, there is an additional regime where small scale flux tubes can be well confined and little energy is extracted or dissipated. In the intermediate confinement regime, the black hole twists up the field and forms a magnetic tower, which gets disrupted by $m=1$ (kink) instability. In this particularly interesting regime, energy extraction increases but most of it is dissipated in the closed region disturbed by the instability. This demonstrates that, in principle, with the field configuration near the black hole being filled by small scale flux tubes, it is possible that the rotational energy is extracted from the black hole but gets dissipated relatively nearby, instead of forming an outflow in the form of a Blandford-Znajek jet. This may be relevant for the heating of the compact X-ray emitting corona in Seyfert Galaxies.

We have seen that in the regime where the kink instability disrupts the outflow, a significant fraction of the extracted rotational energy $P\sim P_{0}=\Omega^{2}\Phi^{2}/c$ gets dissipated in the magnetosphere, especially at the coherent current sheets formed by the kink instability. These current sheets also rotate around the central compact object, with a rate that is determined by the kink time scale and typically some fraction of the angular velocity of the black hole. This could lead to quasi-periodic signatures in the observed radiation, especially if the radiation mechanism at the current sheet is to some extent beamed. However, in an actual accretion disk the flux tubes may be much more irregular and variable, leading to more stochasticity in the observed signal.

The main dissipation mechanism is likely to be reconnection at these current sheets. Since the environment around the black hole is filled by soft optical/UV photons from the accretion disk, particles accelerated in the reconnection process will be quickly cooled through inverse Compton scatterings and may also initiate pair cascades on the soft photons. The reconnection outflow may be strongly Compton dragged, and the reconnection proceeds in a regime as described by Beloborodov (2017) where the spectrum is dominated by bulk Comptonization from the reconnection outflow and has a high energy cutoff near the pair production threshold.

To get an estimation of the dissipation rate in realistic systems, let us consider a typical Seyfert Galaxy, MCG-6-30-15, which has a supermassive black hole of mass $M=2\times 10^{6}M_{\odot}$, bolometric luminosity $L\sim 8\times 10^{43}\,\rm{erg}\,\rm{s}^{-1}\sim 0.4L_{\rm Edd}$, and X-ray luminosity $L_{X}\sim 2\times 10^{43}\,\rm{erg}\,\rm{s}^{-1}\sim 0.1L_{\rm Edd}$(Reynolds et al., 1997). Measuring $B$ in terms of the characteristic field $B_{\rm Edd}\equiv\sqrt{8\pi p_{\rm Edd}}=\sqrt{2L_{\rm Edd}/3cr_{g}^{2}}=\sqrt{m_{p}c^{2}/r_{e}^{2}r_{g}}=3.6\times 10^{5}M_{6}^{-1/2}$ G (Phinney, 1983), we get

[TABLE]

If the magnetic field near the black hole reaches a value $B\sim B_{\rm Edd}$, then the amount of power dissipated due to the kink instability will be enough to account for all the X-ray emission.

The kind of relatively strong magnetic field with small coherence length scales may be a generic outcome of MRI in the accretion disk, or associated with hot spots on the accretion disk as the flux tubes emerge buoyantly similar to solar magnetosphere. As we have seen from GRMHD simulations, when the accretion disk is threaded by a strong net magnetic flux, it tends to become magnetically arrested and form a powerful jet (McKinney & Blandford, 2009; Avara et al., 2016); however, in situations where the net magnetic flux is small, the disk magnetic field may be dominated by the locally amplified components. MRI and the subsequent buoyancy or gravitational instabilities can produce fields with alternating polarities (e.g., Davis et al., 2010; Zhu & Stone, 2018); they may also be quite asymmetric. As these structures are fed to the black hole, we may see continuous flux tube inflation and kink instability happening on a much more stochastic basis. We have seen from the simulations that the kink instability tends to happen when the confinement from the surrounding field is strong; this confinement can be produced by magnetic field from currents in the outer disk, or a disk wind that exerts pressure towards the central axis region.

We have so far approximated the black hole as a resistive membrane in a flat spacetime model. Whether the other GR effects are important will be tested using a time-dependent, GR force-free code in the future. Meanwhile, we used an ad hoc way to deal with the current sheet where force-free condition is violated. When the non-force-free regions have negligible volume, it can be shown that the overall dynamics does not depend sensitively on the dissipation mechanisms in these boundary layers, and force-free evolution gives a reasonable zero-th order picture (East et al., 2015; Zrake & East, 2016; Nalewajko et al., 2016) 555Force-free evolution becomes problematic if a large volume violates the force-free condition (Lyutikov et al., 2017; Li et al., 2018).. However, to get a better estimation of the energy dissipation rate, we need to use more realistic prescriptions for the dissipation in the force-free code. Another limitation of our approach is that we use the force-free electrodynamics, neglecting the inertia of the plasma attached to the coronal field lines. In regions where the mass loading is small such that the Alfvén speed is close to c and much larger than the sound speed, force-free is a good approximation. We can get a rough estimation of the plasma magnetization for actual astrophysical systems like MCG-6-30-15 (see also similar discussion in Yuan et al., 2019). The magnetic field near the black hole may be on the order of $B\sim 0.3B_{\rm Edd}\sim 10^{5}$ G; the density of the electrons $n_{e}$ can be estimated from the Compton y parameter $y=4kT\rm{Max}(\tau_{\rm es},\tau_{\rm es}^{2})/(mc^{2})$>1, where $\tau_{\rm es}=n_{e}\sigma_{T}r$ is the Thomson scattering optical depth, $T$ is the temperature of the electrons, and $r$ is the size of the region. Since $\tau_{\rm es}\gtrsim 1$, we get

[TABLE]

So the electron magnetization is

[TABLE]

If proton number density is the same as the electron number density, the total magnetization would be order 1. But the corona emitting region may well be electron/positron pair dominated (Zdziarski et al., 1996; Beloborodov, 2017), meaning that the total magnetization remains high. In such situations, the force-free approximation gives a reasonable description. However, another complication is that the mass loading could be quite large for flux tubes that are directed inward at small enough angles to the plane of the disk for which gravitational attraction will overcome centrifugal force and allow disk gas to be unstable to flow along the magnetic field lines (Blandford & Payne, 1982; Yuan et al., 2019). This could change the characteristic wave speeds along the magnetic field, affecting the criterion as to which parts of the flux tubes are more likely to break up and reconnect when the flow becomes supersonic, a process that could enhance the amount of dissipation on the axis near the black hole. A MHD approach is needed to capture such dynamics. Finally, here we do not include the accretion disk motion and field evolution self-consistently; a better approach would be to model the disk using full MHD and allow all scales of magnetic structures to form and evolve self-consistently. These will be investigated in future works.

Acknowledgements

We thank Omer Bromberg, Alexander Yuran Chen, Dmitri Uzdensky for helpful discussions, and the anonymous referee for useful suggestions. YY acknowledges support from the Lyman Spitzer, Jr. Postdoctoral Fellowship awarded by the Department of Astrophysical Sciences at Princeton University. AS is supported by NASA ATP grant 80NSSC18K1099 and by Simons Foundation (grant 267233). DRW is supported by NASA through Einstein Postdoctoral Fellowship grant number PF6-170160, awarded by the Chandra X-ray Center, operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. The simulations presented in this work used computational resources supported by the PICSciE-OIT High Performance Computing Center at Princeton University.

Appendix A Steady state tests of the membrane boundary condition

To test the membrane boundary condition, we first apply it to a relaxation scheme that solves the steady-state, Grad-Shafranov equation in flat spacetime (we have taken $c=1$):

[TABLE]

Here $\Psi$ is the flux function; $\omega(\Psi)$ is the angular velocity of the field line and $I(\Psi)$ is the enclosed poloidal current.

In what follows, we show one example where the membrane is spherical with a radius $r_{*}=1$ and rotating with angular velocity $\Omega=0.5$; the disk extends from $r_{\rm in}=2r_{*}$ and has an angular velocity profile $\Omega_{d}=r^{-3/2}$; the magnetic flux on the disk is fixed to be $\Psi(r>r_{\rm in},\theta=\pi/2)=(r-r_{0})/r_{\rm in}$, where $r_{0}=2r_{\rm in}$ is the location of the separatrix field line. In the steady state solution, the angular velocity of the field line is determined by the disk angular velocity at the foot point, while the current $I(\Psi)$ is given by the boundary condition (9). Specifically, in the slow rotation limit, we need

[TABLE]

Since outside the membrane,

[TABLE]

Equations (42)(43) then give

[TABLE]

Figure 14 shows an example of the steady state solution when we take the resistivity $R=4\pi$. As a comparison, Figure 15 shows an example of the exact solution in Kerr spacetime. Both are obtained using the method described in Yuan et al. (2019). The qualitative similarity suggests that as far as the electromagnetic property of the black hole is concerned, it indeed behaves like a resistive membrane.

Appendix B The case of a dipole field

When the current distribution in the accretion disk is only a single unidirectional loop, the magnetic field threading the central object will always open up, no matter how large the resistivity of the membrane is, how far away the current loop is from the central object, and how fast the central object is spinning (as long as $\Omega>0$).

We can first understand this under axisymmetry constraint, using an argument similar to Uzdensky (2005). If the field lines attached to the central object were to remain closed, the small bundle of field lines near the polar region would connect to the disk at large radius $r\to\infty$. We will show that the steady state Grad-Shafranov equation (41) cannot be satisfied there. Firstly, these field lines have $\omega\to 0$ following their foot points on the disk (our examples in the main paper has a non-rotating disk; however, this argument is still true for a Keplerian disk). Secondly, it is reasonable to assume that near the polar region, $\Psi\propto\theta^{\gamma}$, then boundary condition (48) suggests that

[TABLE]

So the right hand side of equation (41) is

[TABLE]

Now the left hand side of equation (41), without the $\omega$ term, is a diffusion-like operator and roughly scales as $\Psi/r^{2}$. As a result, near the foot points of the polar field lines on the disk, the left hand side is $\sim 1/r^{2}$ times the right hand side and the equation cannot be satisfied for very large $r$. We can also see that for the equilibrium to exist, we need the foot point distance to not be much larger than $r\sim\eta/\Omega$.

In our time-dependent force-free simulations, we indeed get steady, axisymmetric final solutions for these cases. Figures 16 and 17 show two examples. As the dipole field opens up, there exists a current sheet separating the outgoing and ingoing magnetic fluxes; the field attached to the central object is rotating while that attached to the accretion disk is not rotating (as we set $\Omega_{d}=0$). At large distances, the current sheet lies along a polar angle around $\theta=60^{\circ}$. We can understand this opening angle as follows. Since the initial current distribution is isolated, the rotation induced outgoing flux solution is asymptotically split-monopole-like. Due to magnetic flux conservation, the amount of outgoing magnetic flux should be equal to the ingoing flux. Therefore, the separation lies at $\theta\approx 60^{\circ}$.