The toy model

We follow the approach of HPB to study the dragging of an external magnetic field by an MHD turbulent accretion disk. As already noted, this results in a non-relativistic model but should be able to provide quantitative insights on the behavior of slowly rotating black holes (where the radius of marginal stability is in a rather low gravitational redshift regime). We consider a thin Keplerian accretion disk with geometric thickness $h\ll r$ extending down to the radius of marginal stability at $r=r_{\rm ms}$. We suppose that the disk has an effective magnetic diffusivity (due to reconnection in the MHD turbulence) of $\eta_*$ which is comparable to the effective viscosity $\nu_*$. In other words, the effective magnetic Prandtl number $P_{\rm m}=\nu_*/\eta_*$ is of order unity (see HPB for an explicit justification of $P_{\rm m}\sim 1$). Note that we employ the standard definition of ${P_{\rm m}}$ as the ratio of the viscosity to the magnetic diffusivity which is the reciprocal of that used in Lubow et al. (1994) and HPB.

We now suppose that an external uniform magnetic field with strength $B_0$ is present in the vertical direction (i.e., aligned with the normal to the accretion disk). We are interested in the dragging of this field by the accretion flow. Assuming the system remains axisymmetric at all times, and employing cylindrical polar coordinates $(r,z,\phi)$, the poloidal magnetic field structure is completely described by the flux function $A(r,z;t)$ via ${\bf B}_p=\nabla\times (A\hat{\bf\phi}/r)$. With such a definition, the magnetic flux threading a ring of radius $r$ at height $z$ from the disk plane is $2\pi A(r,z;t)$ and the magnetic field components are given by

$${\rm B}_r=-\frac{1}{r}\frac{\partial A}{\partial z},$$ (2)

$${\rm B}_{\rm z}=\frac{1}{r}\frac{\partial A}{\partial r}.$$ (3)

The flux function can be decomposed into three components,

$$A(r,z;t)=A_{\rm BH}(r,z;t)+a(r,z;t)+\frac{r^2B_0}{2}$$ (4)

where $A_{\rm BH}(r,z;t)$ is the flux function associated with cleaned black hole-threading field (generated by currents in the disk), the final term on the RHS is just the uniform imposed flux (generated by currents at infinity), and $a(r,z;t)$ accounts for all other (disk-threading) magnetic field structures (generated, in principle, by currents either in or out of the disk plane). For definiteness we suppose that, in the region exterior to the disk ($\vert z\vert>h$), the black hole-threading field has the form of a split monopole,

$$A_{\rm BH}(r,z;t)=A_*(t)\left(1-{\rm sgn}(z)\frac{z}{(z^2+r^2)^{1/2}}\right)\hspace{2cm}(\vert z\vert>h),$$ (5)

where $A_*(t)$ is $1/2\pi$ times the total hole-threading flux. To reiterate, this hole-threading field is generated by toroidal currents flowing in the disk ($\vert z\vert<h$) and is a vacuum solution to Maxwell's equations elsewhere. While the precise structure of the cleaned black hole field is unclear, the choice of the split monopole has support from recent General Relativistic MHD simulations (e.g., see Hirose et al., 2004; Komissarov 2005).

HPB showed that the time evolution of the flux function in the diffusive part of the disk plane, $A(r,0;t)$ ($r>r_{\rm ms}$), is given by

$$\frac{\partial A}{\partial t}+v_r\frac{\partial A}{\partial ... ...eft(\frac{1}{r}\frac{\partial A}{\partial r}\right) \right]=0,$$ (6)

where $v_r$ is the radial velocity of the accretion flow. The $(\partial A/\partial z)_{z=h}$ term represents the effect of magnetic tension due to the curvature of field lines across the disk plane, and thus depends on the structure of the magnetic field above and below the disk. For example, one could formulate steady-state MHD wind solutions which take the instantaneous value of $A(r,0;t)$ as a boundary condition. This would be a task of great complexity (note that the radial structure of the boundary condition would not admit self-similar wind solutions). Here, we make the following simplifications. Firstly, we assume that the magnetic field outside of the disk (hereafter referred to as the disk magnetosphere) is force-free, i.e., $(\nabla\times {\bf B})\times {\bf B}=0$. Secondly, we assume that the Alfvén speed in the disk magnetosphere is sufficiently high as to reduce the toroidal field to essentially zero (through the production of torsional Alfvén waves). Setting $B_\phi=0$, the field in the disk magnetosphere becomes potential ( $\nabla\times {\bf B}=0$) and the flux function obeys

$${\cal D}A=0$$ (7)

where ${\cal D}$ is the linear differential operator

$${\cal D}\equiv \frac{\partial}{\partial r}\left(\frac{1}{r}\... ...partial z}\left(\frac{1}{r}\frac{\partial}{\partial z}\right).$$ (8)

Noting that both the imposed uniform field and (exterior to the disk) the black hole-threading field $A_{\rm BH}(r,z;t)$ individually obey ${\cal D}A=0$, the structure of the disk magnetosphere is determined by solving the potential problem for $a(r,z;t)$, i.e., ${\cal D}a=0$.

At this point, a brief discussion of our $B_{\phi}=0$ assumption (which leads to the potential field condition) is in order. In the non-relativistic treatment here, we can always assume that the Alfvén speed in the disk magnetosphere is large enough such that any twist in the magnetic field is removed via a torsional Alfvén wave. However, our intent is to produce a toy-model for accretion onto a black hole so we should be wary of an assumption that so explicitly relies on non-relativistic physics. In a fully relativistic treatment, the force-free magnetosphere around a black hole accretion disk would be described by the Grad-Shafranov equation (BZ; MacDonald & Thorne 1982; Uzdensky 2004, 2005). It is found that the poloidal field structure depends on both the poloidal current distribution (which gives rise to toroidal fields) and the field line rotation (due to the fact that the field lines are frozen into an orbiting accretion disk, for example). In particular, the field structure is affected by the presence of an inner and outer light cylinder. Ultimately, a relativistic version of our model should study the flux dragging and the magnetization of the black hole including these physical effects. Here, we simply note that detailed studies of non-rotating (or slowly-rotating) black hole magnetospheres have shown that the field line rotation associated with a Keplerian accretion disk has only a small effect on the poloidal field as compared with the equivalent non-rotating configuration (MacDonald 1984; Uzdensky 2004). In this sense, Keplerian accretion disks are ``slow rotators'' (Uzdensky 2004).

For the rest of this paper, we explicitly consider the behavior of the magnetic field in the upper half of the $z$-plane, $z>0$ -- we suppose the system to be symmetric in the $z=0$ plane such that $B_z(r,z)=B_z(r,-z)$ and $B_r(r,z)=-B_r(r,-z)$. The tension term in eqn. 6 can be decomposed into

$$\left(\frac{\partial A}{\partial z}\right)_{z=h}=\left(\frac... ...right)_{z=h}+\left(\frac{\partial a}{\partial z}\right)_{z=h}.$$ (9)

The contribution from the hole-threading flux can be evaluated directly from eqn.5,

$$\left(\frac{\partial A_{\rm BH}}{\partial z}\right)_{z=h}\approx -\frac{A_*}{r},$$ (10)

where we have neglected a term which is smaller by a factor of $(h/r)^2$. The remaining contribution to eqn. 9 follows from the solution to the potential problem ${\cal D}a=0$ with boundary conditions $a(r=0,z;t)=0$ and $a(r,0;t)$ specified. As shown by HPB, this gives

$$\left(\frac{\partial A}{\partial z}\right)_{z=h}={\cal P}\in... ...x\frac{[a(x,0;t)-a(r,0;t)]}{\pi(r-x)^2}-\frac{a(r,0;t)}{\pi r}$$ (11)

where ``$\cal P$'' signifies the principal part of the integral. We can now write an explicit integro-differential equation for the time evolution of $a(r,0;t)$ in the diffusive part of the disk ($r>r_{\rm ms}$);

$$\frac{\partial a}{\partial t}+\frac{\partial A_*}{\partial t}+v_rrB_0+\left(v_r+\frac{\eta_*}{r}\right)\frac{\partial a}{\partial r} =$$

$$\eta_*[\frac{1}{h}{\cal P}\int_0^\infty dx\frac{[a(x,0;t)-a(r,0;t... ...- \frac{a(r,0;t)}{h\pi r}-\frac{A_*(t)}{hr} +\frac{\partial^2a}{\partial r^2}].$$ (12)

As part of our model, we must specify $h(r)$, $v(r)$ and $\eta_*(r)$. For definiteness, we define $h(r)$ by taking the ratio $h/r$ as a fixed parameter of our model (in principle, one could substitute a particular form for $h(r)$ resulting from a detailed disk model). To specify the radial velocity field, we follow Lubow et al. (1994) and split our disk into two zones which we dub an ``active'' and a ``dead'' zone. In the active zone ( $r_{\rm ms}<r<r_{\rm dead}$), we set $v_r=-\nu_*(1/r-1/r_{\rm dead})$ where $\nu_*=\alpha h^2(GM/r^3)^{1/2}$ (HRB), and $\eta_*=\nu_*/{P_{\rm m}}$. The magnetic Prandtl number ${P_{\rm m}}$ is a fixed and constant parameter of the active disk. Note that we have introduced the usual $\alpha$ of accretion disk theory (in contrast with HPB who implicitly employ $\alpha\sim 1$). In the dead zone ( $r_{\rm dead}<r<r_{\rm out}$), the diffusivity is still given by $\eta_*=\alpha h^2(GM/r^3)^{1/2}/{P_{\rm m}}$, but the velocity is set to zero. For computational necessities, we impose an outer cutoff on the system at $r=r_{\rm out}$. We assume that the disk beyond $r_{\rm out}$ is a perfect and static conductor. Hence the total magnetic flux threading a loop ( $r=r_{\rm out},z=0$) is constant and has the value $\pi r_{\rm out}^2B_0$. The inclusion of the dead-zone makes the evolution of the inner part of the system essentially independent of the position or exact nature of the $r=r_{\rm out}$ boundary. In particular, the dead zone acts as a reservoir of magnetic flux that can feed the actively accreting part of the disk -- only in the outermost parts of the active disk does the conservation of magnetic flux lead to a non-negligible magnetic pressure trying to ``suck'' magnetic flux out of the active disk. The physical nature of the dead zone will be discussed in Section 4.

Finally, we must specify boundary conditions on $a(r,0;t)$. The implementation of the inner boundary condition must capture the fact that the plunge region is extremely effective at sweeping in poloidal magnetic field that crosses within $r=r_{\rm ms}$. Consider a poloidal magnetic field line which is dragged towards the plunge region on the viscous timescale $t_{\rm visc}\approx (r_{\rm ms}/h_{\rm ms})^2(r_{\rm ms}^3/GM)^{1/2}\alpha^{-1}$. Once in the plunge region, the radial velocity of the disk material rapidly increases with no associated increase in the effective magnetic diffusivity (indeed, to the extent that the plunge region becomes a laminar rather than a turbulent flow, the effective magnetic diffusivity may well plummet to very small values). For the field strengths under consideration here (i.e., with an energy density much less that the kinetic energy density of the accretion flow) inward advection of the field line on a dynamical timescale $t_{\rm dyn}\approx (r_{\rm ms}^3/GM)^{1/2}$ will dominate all other processes. Since the characteristic evolution timescale of the system is $t_{\rm visc}\gg t_{\rm dyn}$, flux conservation gives that the vertical magnetic field in the $z=0$ plane in the plunge region compared with that in the disk just outside is

$$\frac{B_z({\rm plunge})}{B_z({\rm disk})}\approx \frac{t_{\r... ...n}}{t_{\rm visc}}\approx\alpha\left(\frac{h}{r}\right)^2\ll 1.$$ (13)

To a good approximation, we can say that the magnetic flux locally crossing the plunge region is zero. Thus, the only magnetic flux passing through a loop $(r<r_{\rm ms},z=0)$ is that which threads the black hole, i.e., $A(r\le r_{\rm ms},0;t)=A_*(t)$. To cancel the contribution from the externally imposed uniform field in this region, we must have

$$a(r,0;t)=-r^2B_0/2\hspace{2cm}(r<r_{\rm ms}).$$ (14)

Thus, the appropriate inner boundary condition for eqn 12 is $a(r=r_{\rm ms},0;t)=-r_{\rm ms}^2B_0/2$ and we must use eqn. 14 in the evaluation of the integral term of eqn. 12. The fact that $\partial a(r_{\rm ms},0;t)/\partial t=0$ allows us to use eqn. 12 to evaluate the rate of change of black hole-threading flux,

$$\frac{\partial A_*}{\partial t}= \eta_*(r_{\rm ms})\left[\fr... ...}{hr} +\frac{\partial^2a}{\partial r^2}+B_0\right]_{r=r_{ms}},$$ (15)

where we have used the continuity of $\partial a/\partial r$ across $r=r_{\rm ms}$ to combine the third and fourth terms on the left hand side of eqn. 12. We can justify this assumption of continuity as follows. Suppose that this derivative was discontinuous across $r=r_{\rm ms}$, resulting in a discontinuity in the strength of the vertical magnetic field. This would lead to a large magnetic pressure gradient and a very rapid rearrangement of material until continuity was achieved. We do note, however, that we expect a rather narrow transition zone just outside of $r=r_{\rm ms}$ where vertical magnetic field goes from zero to the value characteristic of the disk. We must spatially resolve this transition in our numerical model.

For the outer boundary condition, we set $a(r_{\rm out},z=0;t)=-A_*(t)$ for some $r_{\rm out}>r_{\rm dead}$. This amounts to bounding the entire system by a perfect and static conductor in the disk plane ($z=0$) for all $r>r_{\rm out}$, as discussed above.

With these assumptions, eqns. 12 and 15 completely describe the evolution of $a(r,0;t)$ and $A_*(t)$ from some initial state once we fix the magnetic Prandtl number ${P_{\rm m}}$, the disk thickness $h/r$, the characteristic radii of the problem ($r_{\rm ms}$, $r_{\rm dead}$, $r_{\rm out}$), the external field strength $B_0$, and the viscosity parameter $\alpha$. In fact $\alpha$ and $B_0$ are trivial parameters of the model, affecting only the scaling of the time coordinate and the absolute normalization of $a$, respectively. Furthermore, the inclusion of the dead-zone makes the evolution of the inner disk/field essentially independent of the location of the outer boundary $r=r_{\rm out}$. Hence, the non-trivial parameters describing this system are ${P_{\rm m}}$, $h/r$, and $r_{\rm dead}$. For our initial condition, we take

$$a(r,z=0,t=0)=\left\{ \begin{array}{ll} -r^2B_0/2 & (r<r_{\rm ms})\\ -r_{\rm ms}^2B_0/2 & (r\ge r_{\rm ms}) \end{array}\right.$$ (16)

This amounts to saying that the initial currents flowing in the disk are only those required to cancel the imposed uniform field in the plunge region.