TRANSPORT OF LARGE-SCALE POLOIDAL FLUX IN BLACK HOLE ACCRETION

In this study, we carry out both axisymmetric and fully three-dimensional simulations around a Schwarzschild black hole. The initial condition, as in previous works (De Villiers et al. 2003, 2005; Hirose et al. 2004; Krolik et al. 2005; Hawley & Krolik 2006; Beckwith et al. 2008a, 2008b), consists of an isolated gas torus upon which a weak magnetic field is imposed. We use the torus model described by De Villiers et al. (2003) with an adiabatic index of Γ = 5/3. The initial angular momentum distribution is slightly sub-Keplerian, with a specific angular momentum at the inner edge of the torus (located at r = 30M) ℓ_in ≡ U_ϕ/U_t = 6.10. For this choice of ℓ_in, the pressure maximum is at r ≈ 40M, and the outer edge of the torus is located at r ≈ 60M. The angular period at the pressure maximum is Ω⁻¹ ≈ 250M. With this inner edge radius, the resulting accretion flow should have sufficient radial range to permit a detailed study of the advection of vertical flux. The choice of angular momentum parameter at the inner edge of the torus, ℓ_in, was made so that the initial torus (and hence the steady state accretion flow) would have a scale height similar to that of previous simulations (see, e.g., Beckwith et al. 2008b). The torus is initially surrounded by zones filled with zero-velocity gas set to a constant value of density and pressure. The ratios of these "vacuum" values to the torus maximum are 2 × 10⁻⁴ for the pressure, and 5 × 10⁻⁶ for the density. Because the code is relativistic, the time step is always limited by the time for light to cross the smallest cell; high Alfvén speeds in the low-density coronal region do not create any unique difficulty.

The initial magnetic field is homogeneous and vertical, filling the annular cylinder 35M ⩽ R ⩽ 55M (where R is the cylindrical radius). This configuration differs from the vertical field used in some previous studies (see, e.g., Koide et al. 1999; McKinney & Gammie 2004; De Villiers 2006) which fills all space, not just through the initial torus. By isolating the field in this manner, we know that any vertical flux that ends up crossing the black hole event horizon was brought in from field initially passing through the torus at large radius, rather than already being present at the horizon in the initial conditions. Also, because equilibrium field strengths scale strongly with radius, a global vertical field that is strong enough to be interesting at one radius is likely to be overwhelming at some radii while negligible at others. We prefer to see the global field strength emerge as an outcome of the evolution. In our configuration, the initial vertical field is subthermal within the torus and is unstable to the MRI. It should be noted, though, that this initial field configuration is not in dynamic equilibrium outside of the torus. Also, this is but one very specific realization of a torus embedded in large-scale field. Future experiments will need to explore a wider variety of initial configurations.

This initial field is computed from the curl of the four-vector potential, i.e.,

$$\begin{equation} F_{\alpha \beta } = \partial _\alpha A_\beta - \partial _\beta A_\alpha = A_{\beta, \alpha } - A_{\alpha, \beta }, \end{equation} \tag{ 1 }$$

with only A_ϕ ≠ 0. The vector potential corresponding to a uniform (Newtonian) vertical field geometry is proportional to the cylindrical radius, i.e.,

$$\begin{equation} A_{\phi } \propto r \sin \theta. \end{equation} \tag{ 2 }$$

The vector potential for the initial field is

$$\begin{equation} A_{\phi } = A_0 \left\lbrace \begin{array}{@{}lll} R_{\rm in} & \quad \mbox{$R < R_{\rm in}$}, \\ r \sin \theta & \mbox{$R_{\rm in} \le R \le R_{\rm out}$}, \\ R_{\rm out} & \quad \mbox{$R > R_{\rm out}$},\end{array} \right. \end{equation} \tag{ 3 }$$

where R_in and R_out are 35M and 55M, respectively. Outside of those cylindrical radii, the field is zero. A₀ was specified so that the ratio of the average thermal to average magnetic pressure inside the torus is β = 100. At the center of the torus β = 275, and β declines smoothly toward 1 near the edge of the torus. In the surrounding corona, β = 0.038. The initial mass, field, and β distribution is shown in Figure 1.

Zoom In Zoom Out Reset image size

The simulation we will focus on is three-dimensional with 256 × 256 × 64 zones in (r, θ, ϕ). The simulation is labeled "VD0m3d." Given the challenge of adequately resolving such models and of doing anything like a convergence study, we also carried out two-dimensional simulations at two different resolutions. In two dimensions, we used 256 × 256 zones and 512 × 512 zones in (r, θ). The two-dimensional simulations are labeled "VD0m2d" and "VD0h2d," where the "m" designates medium resolution and the "h" designates the higher resolution case.

For all models, the radial grid extends from an inner boundary at r_in = 2.104M to an outer boundary located at r_out = 500M. The radial grid was graded using a logarithmic distribution that concentrates zones near the inner boundary. An outflow condition is applied at both the inner and outer radial boundaries. In contrast to the zero net flux simulations that we have performed previously, magnetic fields pierce the outer boundary in the initial state. This net flux is kept fixed throughout the simulation, but otherwise the field at the boundary can evolve. When the fluid velocity normal to the boundary is directed outward, we perform a simple copy of fluid variables into the ghost zones that lie outside of the boundary. In the case where the fluid velocity is inward, the ghost zones are filled with a cold, low density, zero momentum fluid at the vacuum value so that only cold, very low density gas can enter the computational grid. To minimize the influence of the outer boundary, it is located as far away from the region of interest as is feasible. Any torques that arise from the interaction of the outer boundary with the large-scale vertical field should have minimal influence on the turbulent disk itself.

The θ-grid spans the range 0.01π ⩽ θ ⩽ 0.99π. This creates a conical cutout that surrounds the coordinate axis. Such a cutout will prevent flow through the coordinate axis, but it is computationally advantageous to avoid the coordinate singularity. The size of this θ-cutout is determined by two considerations. First, light travel times across narrow ϕ-zones near the axis set the time step, so the cutout should not be too small. In this simulation, Δt = 2.53 × 10⁻³M. The θ-cutout must not be too large; initial vertical field should not cross the θ-cutout before reaching the outer radial boundary at r = 500M. A reflecting boundary condition is used along the θ-cutout. The θ-zones were distributed logarithmically so as to concentrate zones around the equator.

The ϕ-grid spans the quarter plane, 0 ⩽ ϕ ⩽ π/2. The ϕ-grid is uniform and periodic boundary conditions are applied. The use of this restricted angular domain significantly reduces the computational requirements of the simulation. Schnittman et al. (2006) examined the variance in surface density as a function of ϕ and found the characteristic size of perturbations to be around 25°. Thus, while some global features may be lost by using a restricted domain, the character of the local MRI turbulence is captured.

In analyzing the evolution of the large-scale poloidal flux, we will make use of the poloidal flux function, which corresponds to the ϕ-component of the magnetic vector potential (van Ballegooijen 1989; Heyvaerts et al. 1996; Reynolds et al. 2006), which we must calculate from the simulation data. In GRMHD, we directly evolve the components of the Faraday tensor, F_μν, using the alternative form of the induction equation (for further details see De Villiers & Hawley 2003):

$$\begin{equation} \partial _\delta F_{\alpha \beta } + \partial _\alpha F_{\beta \delta } + \partial _\beta F_{\delta \alpha } = 0. \end{equation} \tag{ 4 }$$

The space–space components of F_μν are identified with the CT magnetic fields via $\mathcal {B}^{i} = [ijk] F_{jk}$:

$$\begin{equation} {\cal B}^r = F_{\phi \theta }; \;\; {\cal B}^\theta = F_{r \phi }; \;\; {\cal B}^\phi = F_{\theta r}. \end{equation} \tag{ 5 }$$

This identification allows us to write the induction equation in the familiar form

$$\begin{equation} \partial _t {\cal B}^{i} - \partial _j (V^{i} {\cal B}^j - {\cal B}^i V^{j}) = 0, \end{equation} \tag{ 6 }$$

where Vⁱ = Uⁱ/U^t. The tensor F_μν is related to the magnetic induction in the fluid frame, B^α, by the relation

$$\begin{equation} F_{\mu \nu } = \epsilon _{\alpha \beta \mu \nu }\,B^{\alpha }\,U^\beta, \end{equation} \tag{ 7 }$$

where $\epsilon _{\mu \nu \delta \gamma }=\sqrt{-g}\,[\mu \,\nu \,\delta \,\gamma ]$.

The Faraday tensor can also be written in terms of a vector potential:

$$\begin{equation} F_{\alpha \beta } = \partial _\alpha A_\beta - \partial _\beta A_\alpha = A_{\beta, \alpha } - A_{\alpha, \beta }. \end{equation} \tag{ 8 }$$

It is convenient to define an azimuthally averaged poloidal flux function, A_ϕ. This can be derived from the azimuthally averaged CT field data by noting that the total (spatial) derivative of A_ϕ is

$$\begin{equation} dA_\phi = \frac{\partial A_\phi }{\partial r} dr + \frac{\partial A_\phi }{\partial \theta } d\theta. \end{equation} \tag{ 9 }$$

From the definition of the Faraday tensor, we have

$$\begin{equation} F_{\phi \theta } = A_{\theta,\phi } - A_{\phi,\theta }; \;\; F_{r \phi } = A_{\phi,r} - A_{r,\phi }. \end{equation} \tag{ 10 }$$

After azimuthal averaging, this reduces to

$$\begin{equation} F_{\phi \theta } = - A_{\phi,\theta }; \;\; F_{r \phi } = A_{\phi,r}. \end{equation} \tag{ 11 }$$

The space–space components of F_μν are identified with the CT magnetic fields via $\mathcal {B}^{i} = [ijk] F_{jk}$ so that

$$\begin{equation} {\cal B}^r = F_{\phi \theta } = -A_{\phi,\theta }; \;\; {\cal B}^\theta = F_{r \phi } = A_{\phi,r}. \end{equation} \tag{ 12 }$$

We therefore have

$$\begin{equation} dA_\phi = F_{r \phi } dr - F_{\phi \theta } d\theta = {\cal B}^\theta dr - {\cal B}^r d\theta. \end{equation} \tag{ 13 }$$

To find the flux from this differential, we could in principle integrate over any poloidal surface. For the radial motion of vertical flux through the disk, the surface of greatest significance is the equatorial plane. Because this surface is interrupted by the event horizon, we must stretch it over half of the event horizon. We will therefore define a flux function that is the sum of two parts: the radial flux through the inner horizon boundary,

$$\begin{equation} \Psi (\theta)|_{r=r_{\rm in}} = \int ^{\theta }_{0} \, d\theta ^\prime \, {\cal B}^r (r_{\rm in},\theta ^\prime), \end{equation} \tag{ 14 }$$

and the vertical flux through the equator,

$$\begin{equation} \Phi (r)|_{\theta =\pi /2} = \int ^r_{r_{\rm in}} \, dr^\prime \, {\cal B}^\theta (r^\prime,\theta =\pi /2). \end{equation} \tag{ 15 }$$

Using these functions, we can determine how the accumulated flux through the top half of the event horizon depends on polar angle, Ψ(θ), and how much total vertical flux has been brought to within radius r. The full flux function along the total path is

$$\begin{equation} \mathcal {A} (r,\theta =\pi /2) = \Psi (\pi /2) - \Phi (r). \end{equation} \tag{ 16 }$$

This corresponds to the net flux piercing the surface covering the top hemisphere of the black hole horizon and the equatorial plane out to radius r. The full vector potential component A_ϕ(r, θ) is obtained by integrating the CT variables throughout the computational domain starting from r = r_in and θ = 0. $\mathcal {A} (r,\theta =\pi /2)$ is therefore identical to A_ϕ(r, θ = π/2). For the three-dimensional simulation, VD0m3d, ${\cal B}^r$ and ${\cal B}^\theta$are integrated over the ϕ-dimension prior to computing A_ϕ(r, θ). Given the periodic nature of the ϕ direction this is sufficient to compute the total flux piercing the volume bounded by the equatorial plane and the black hole horizon. We note that this procedure is consistent with the approach adopted in the calculations of van Ballegooijen (1989); Heyvaerts et al. (1996); and Reynolds et al. (2006).

The simulation begins with an isolated torus threaded with a column of vertical field. The first part of the simulation is marked by relatively rapid evolution of the coronal field and the surface layers of the torus. Figure 2 shows the evolution of the poloidal magnetic flux and the density distribution over the time period 1000M ⩽ t ⩽ 2500M. This process, along with the subsequent evolution of the flow is also shown in an animation included in the online edition of this work. We encourage the reader to refer to this material as it clarifies many of the issues that we shall discuss at length in the remainder of the text. As noted above, the initial state is not in equilibrium within the corona. There is differential rotation between the coronal gas and the torus, and the magnetic pressure within the vertical field column is not balanced by pressure within the initial atmosphere. The field expands radially, and stresses at the boundary of the disk begin to drive low density gas from the torus outward along field lines. In the inner half of the torus, low density surface layers move inward in two thin accretion streams above and below the equatorial plane, dragging poloidal field lines down to the black hole. This behavior was dubbed an "avalanche flow" when it was observed in the axisymmetric MHD accretion torus simulations of Matsumoto et al. (1996), and it is a common feature to accretion models that begin with vertical fields threading the disk (e.g., Koide et al. 1999). These surface features develop well before field deeper within the torus is significantly amplified by the MRI. The surface layers in the outer half of the torus also evolve, but here mass flows radially outward along the field lines, carrying the field along as it does so.

Zoom In Zoom Out Reset image size

By t = 2000M, net flux has begun to build up on the black hole horizon, well before any significant mass accretion has begun. This flux arrives there by a process we call the "coronal mechanism," which we discuss briefly here and will discuss again in greater detail in Section 3.2. A typical field line participating in the coronal mechanism may be traced from the outer boundary at high latitude more or less vertically down to the upper surface of the torus, where it bends sharply toward the hole near the boundary of the denser region. It goes in some distance before doubling back (forming a "hairpin" shape) and returning to the torus. There it passes through the equator, and turns inward to mirror its behavior on the other side of the equator. A short time later, the apex of each hairpin moves inward and closer to the equator, where, in the early stages of this simulation when there is little mass at small radii, it can encounter an oppositely oriented hairpin from the other hemisphere. Reconnection ensues, causing the magnetic flux distribution to change essentially instantaneously.

The four panels of Figure 2 show these structures in differing ways. In the first panel, an isolated hairpin field line can be clearly seen approaching very close to the horizon. By the time of the second panel, the bends in this field line have reconnected, attaching it to the event horizon and liberating a pair of closed loops. Also in the second panel, new hairpins can be seen forming along the top and bottom surfaces of the disk in the region 20M < r < 30M. These hairpins evolve further and become embroiled in the beginnings of disk turbulence (the channel modes) in the third and fourth panels.

Magnetic field enters the extremely low-density funnel around the rotation axis by an offshoot of the coronal mechanism just described. There is initially no field within the funnel. Once horizontal field lines have been brought close to the black hole, their radial, near-horizontal components become unstable to the formation of ballooning half-loops that rise upward into the funnel. These loops are the initial source of field for the funnel. Only the inner half of the rising loop enters the funnel and, as a result, the field direction changes sign abruptly at the funnel's outer edge along the centrifugal barrier.

The coronal process also occurs at radii exterior to the initial pressure maximum of the torus. There it acts to move net flux progressively outward. Above the disk surface, the field is stretched out radially, but subsequently bends back down toward the equator, where reconnection can occur. This transfers the net flux from field lines going through the torus (those lines now close in large loops) to vertical field lines at large radius.

Thus, in the early phase of the evolution, the global topology of the field is rapidly rearranged through coronal motions without any significant evolution within the main disk. Although these features result from the initial conditions, the process illustrates how (in principle) coherent large-scale poloidal flux can be rapidly moved over large radial distances.

Although low-density surface layers participate in the initial avalanche flow, the dense part of the torus is largely unaffected by these coronal motions. By t = 2500M, however, the MRI reaches a nonlinear amplitude throughout the torus and the characteristic structures associated with the "channel mode" have become visible within the torus. The next stage of evolution is illustrated in Figure 3, which shows snapshots every 500M in time from 3000 to 4500M. Gas in the torus moves radially both inward and outward in the characteristic channels as angular momentum is transported along field lines, which are themselves increasingly stretched into long radial filaments. The channel modes are unstable to nonaxisymmetric perturbations and rapidly breakup, resulting in a turbulent disk without prominent radial features. In the axisymmetric simulations, these relatively coherent large-scale radial flows continue to dominate throughout the evolution. This artificial persistence of coherent motion is a significant limitation of two-dimensional simulations with vertical field.

Zoom In Zoom Out Reset image size

Accretion into the black hole does not begin until 1650M, and the mass accretion rate at the inner edge of the disk reaches its long-term mean for the first time shortly after 4500M. The mass accretion rate then continues to vary by factors of ∼2 around this mean through the remainder of the simulation. The disk achieves an approximate inflow equilibrium within r ≃ 25M from time 10⁴M onward. Prior to 10⁴M, the mass in the inner disk (r ⩽ 25M) rises steadily, but after that time, it changes by about 10%, and without any prevailing trend. The mass accretion rate averaged from 10⁴M until the end of the simulation at 2 × 10⁴M is flat to within 10% for r ⩽ 25M. The total mass on the grid at the end of the simulation is 78% of its initial value; of this, 18.8% has gone into the black hole.

The divergence-free nature of the magnetic field means that the total flux is rigorously conserved except for losses at the boundaries of the computational domain. To measure the evolution of the total flux, we integrate along the black hole horizon from the axis to the midplane and then out along the equator to the outer boundary. Figure 4 shows this integral as a function of time. The figure shows that there is no flux lost through the outer radial boundary during the simulation, but considerable net flux is deposited on the horizon at the expense of flux through the equator. The majority of the flux brought to the horizon arrives by t ≃ 5000M, which implies that the net horizon field is created mostly by coronal motions, not through disk accretion. During the later stages of more or less steady state accretion, the horizon flux grows much more slowly. Relative to its magnitude at t = 5000M, the horizon flux increases by only 28% over the next 5000M and adds only another 11% during the last 10⁴M. By contrast, the total mass accreted onto the black hole grows by a factor of 9 from 5000M to 2 × 10⁴M, with a rate of increase that is approximately constant from ≃10⁴M onward.

Zoom In Zoom Out Reset image size

The appearance of the overall flow at the end of the simulation is given in Figure 5, which shows the azimuthally averaged density distribution and the β parameter, and in Figure 6, which shows the vector potential component A_ϕ at time t = 2 × 10⁴M. There are several distinct regions. First, there is a low density, strongly magnetized funnel in which the field lines are primarily radial. Second, there is a corona characterized by large-scale field lines that loop about creating islands and hairpins, with numerous current sheets separating closely spaced regions of oppositely directed field. Here there are regions of both relatively strong and relatively weak magnetization; β ranges from ∼0.1 to ∼10. A thin region of high β fluid running along the funnel wall marks a current sheet between field lines directed in opposite radial directions as well as the boundary of the region with strong toroidal field. Finally, the main disk body remains weakly magnetized, with β ranging from ∼10 to 100. The disk thickness H/R ranges from 0.1 to 0.15, using the vertical density moment (Noble et al. 2009) for the definition of H. As can be seen from Figures 3 and 6, substantial parts of the disk body (most of the disk inside r ≃ 40M) are actually disconnected from the large-scale field, their magnetic connections having been broken by the reconnection events that transferred field both to the horizon and to the outside of the torus. In the outer half of this region (20M ≲ r ≲ 40M), however, some field lines do pass through the equator and out through the disk and into the corona.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This initial magnetic field and torus configuration was also simulated in axisymmetry. While the initial evolutions of the two- and three-dimensional models are similar, once the MRI sets in within the torus the subsequent evolutions differ dramatically. In axisymmetry, the channel modes dominate throughout the simulation, creating coherent radial field lines that extend over large distances. These field structures strongly influence the accretion of both mass and magnetic flux. Figure 7 compares the flux through the horizon, the mass accretion rate, and the specific angular momentum accreted in simulations VD0m3d and VD0m2d. The axisymmetric simulation is characterized by very large fluctuations in mass accretion rate. Similar large fluctuations are seen in the value of the magnetic flux through the horizon. Note the contrast between the net specific angular momentum j_in carried into the hole by the accretion flow in the two- and three-dimensional simulations. For both plots, the dashed line indicates the specific angular momentum of the innermost stable circular orbit at 6M. The deep minima in the axisymmetric run correspond to strong torques created by radial field extending through the plunging region. These extended radial fields are associated with the coherent channel flows seen in axisymmetry, and they have an obvious strong effect on mass, angular momentum, and magnetic fluxes carried into the black hole. The three-dimensional simulation is qualitatively different. In three dimensions, the channel modes quickly breakup, yielding to genuine turbulence. The accretion rate into the hole fluctuates, but not nearly as strongly as in axisymmetry. Similarly, only in three dimensions can one speak sensibly about a well-defined net accreted angular momentum per unit rest mass, j_in. On the basis of these contrasts, we conclude that axisymmetric simulations, while useful for preliminary investigations of various initial configurations, are of only limited utility in investigating the long-term behavior of MHD turbulent disks.

Zoom In Zoom Out Reset image size

The dominant mechanism for bringing flux to the event horizon in simulation VD0m3d is one that operates primarily in the corona, and not in the disk accretion flow proper. Although most of the flux is brought in during the first 5000M in time, there is a slow increase in flux throughout the remainder of the simulation. The majority of this flux appears to be transported through the continuing operation of the coronal mechanism. The key aspect is that magnetic field is carried inward by infalling low-density fluid. Although this coronal mechanism operates outside the disk where most of the mass accretion happens, it, too, depends on angular momentum transport to proceed, and the basic physical mechanism for angular momentum transport is likewise the same: magnetic torques. Wherever the vertical field is perturbed so as to create a radial component, orbital shear acting on the radial field in turn creates toroidal field. Moreover, whenever the orbital frequency decreases outward, the relative signs of the radial and toroidal fields are always such as to make the r–ϕ element of the magnetic stress tensor, b^rb_ϕ/4π, negative. Such a stress transports angular momentum outward, tending to accentuate inflow, and therefore growth in the radial field component; this is, of course, the mechanism of the magnetorotational instability.

An example of the coronal mechanism in action is given in Figure 8. In this figure, color contours show the strength of the azimuthally averaged Maxwell stress per unit mass, overlaid with field lines derived from A_ϕ(r, θ). For clarity, we focus on those field lines that still pass through the equator in the accretion disk. The first frame, corresponding to t = 6000M, features a field that runs from the equator vertically upward through the disk and then down toward the hole before reversing and moving radially outward. The strong stress per unit mass in the vicinity of the hairpin carries angular momentum away from the associated matter, enabling it to move inward. In the next frame at t = 6200M, the hairpin moves in radially toward the black hole. The next two frames show a close-up view of the innermost part of the accretion flow at times t = 6440M and 6460M. By the time of the first one, the inner portion of the hairpin has reached and crossed the horizon. Only 20M later, the inside hairpin has reconnected, forming a closed field loop. Thus, the net result of the coronal mechanism at this stage is to attach those field lines forming the exterior of the hairpin to the black hole. These field lines then extend out into the corona, ultimately attaching to the disk at much larger radius.

Zoom In Zoom Out Reset image size

Although net local inflow of magnetic field is a necessary condition for flux accretion, it is not a sufficient one; the net flux changes only when reconnection occurs. Local motion alone is insufficient because magnetic flux is an essentially global concept, so local motion is not enough to determine net flux transport. This distinction is emphasized by contrasting the local magnetic inflow rate with global measures. We define the local rate by

$$\begin{equation} V_{\Phi }(r,\theta) \equiv \frac{1}{2\pi }\int \, d\phi V^r {\cal B}^\theta, \end{equation} \tag{ 17 }$$

averaged over the portion of the simulation during which the trapped flux grows slowly (5000M–2 × 10⁴M). This is shown in Figure 9, for which the radius was chosen to be 15M. As this figure shows, the local flux inflow speed is considerably greater just outside the funnel (|π − θ| ≃ 0.2π), where it is ≃(0.5–1) × 10⁻⁴, than in the midplane, where it is nearly 2 orders of magnitude smaller.

Zoom In Zoom Out Reset image size

This local rate of flux inflow is, however, much larger than the global rate of flux accumulation. The global rate is ∂A_ϕ/∂t, which can be inferred from Figure 4. For the purpose of this figure, the magnetic flux was separated into the flux piercing the midplane, Φ, and the flux through the upper hemisphere of the black hole, Ψ. During the early stages of the simulation (before 5000M), ∂Ψ/∂t ≃ 1 × 10⁻⁵, but this rate of increase diminishes at later times by at least an order of magnitude. Thus, the local inflow rate, V_Φ(θ ≃ 0.2π, 0.8π) ≃ 1 × 10⁻⁴ at r = 15M during the latter 3/4 of the simulation, is as much as 2 orders of magnitude greater than the net global rate. In other words, local field inflow alone does not suffice to create global flux inflow.

Likewise, reconnection on its own is also a necessary condition for flux accretion, but not a sufficient one. In order for flux accretion to occur, reconnection must result in a global rearrangement of the field topology. This fact is illustrated by the field structures displayed in Figure 8. In that figure, a narrow current sheet lies between the two sides of the hairpin, permitting reconnection to happen there. In order to transfer magnetic flux inward through the accretion flow and increase the net flux in the funnel, reconnection must preferentially destroy the side of the magnetic hairpin that lies closest to the midplane. This is possible because reconnection can also occur across the equator, something that the funnel field cannot do. Figure 10 illustrates the idea. Coronal hairpin field loops bring field down to the inner disk, returning along the disk back out to large radius where they eventually pass through the equator and connect to another hairpin in the opposite hemisphere. The innermost bend of these hairpins can attach to the horizon. When this happens, the field appears to be almost a closed loop, with one end of the loop inside the horizon (e.g., the two inner-disk field lines in the gray region of the diagram). If the two oppositely directed portions of one of these loops can meet at intermediate radius, the equatorial flux jumps inward. This inward jump drastically shortens the time until the entire field loop can accrete onto the black hole. When that happens, flux of one sign is removed from the horizon above the equator and flux of the other sign is removed from below the equator, leaving a net increase in flux piercing each horizon hemisphere.

Zoom In Zoom Out Reset image size

An example of a reconnection event within VD0m3d that leads to this sort of global flux rearrangement is shown in Figure 11. In this figure, color contours denote the radial component of the CT-field, ${\cal B}^r$, indicating the direction of the field line, which are again overlaid with white contours derived from the azimuthally averaged A_ϕ(r, θ). Initially (t = 14600M), field lines connect large radii to the black hole event horizon through the funnel region (deep blue region below midplane) and then proceed back out into the disk close to the midplane (yellow region below midplane). These field lines cross the midplane at large radius and proceed back to the black hole north of the midplane (light blue region north of midplane), before finally proceeding out to large radius through the funnel region (red region north of midplane). At time t = 14620M, the oppositely directly field lines that lie just above and below the midplane have reconnected at r ∼ 6M, and the resulting closed field loop (which is still attached to the black hole) is rapidly accreted. This process transfers flux from the outer disk to r ∼ 6M (where the reconnection event occurs) and then down to the black hole as the field loop is accreted, arriving at the black hole with the outer edge of the loop.

Zoom In Zoom Out Reset image size

The coronal mechanism resembles the flux advection mechanism proposed by Rothstein & Lovelace (2008) in the sense that most flux transport takes place outside of the disk body, and that flux in the corona moves inward because magnetic stresses send angular momentum off to infinity. They emphasize the contrast between the turbulent disk and non-turbulent flow outside of the disk. We likewise observe that large-scale, relatively laminar flows within the low-density corona easily move field about; flux-freezing is an excellent approximation there. On the other hand, there are also several points of contrast between our results and the model of Rothstein & Lovelace (2008). They argue that in the disk corona, defined as the region where β ∼ 1, the T^z_ϕ component of the magnetic stress tensor could transport angular momentum outward as effectively as the T^R_ϕ component to which attention is more often directed, allowing the disk surface to move inward. Although this is similar to what we see in our simulation during the initial avalanche phase, it is different from the dynamics prevailing at later times. During the accretion phase of the simulation, the coronal field lines develop MRI-like bends at high altitudes above the disk (≫H), so that the inflow occurs in the corona rather than along the disk surface. In addition, as shown by Figure 12, which contrasts the azimuthal averages of −b^rb_ϕ/ρ and −b^θb_ϕ/ρ, there is very little net angular momentum flux per unit mass in the polar angle direction except in the funnel (at latitudes not too far from the midplane, b^θb_ϕ ≃ b^zb_ϕ). Although the magnitude of T^θ_ϕ is often comparable to T^r_ϕ, it has very nearly the same probability of having either sign, so that its mean value is small. On the funnel edge, T^θ_ϕ is more likely to take angular momentum away from the disk but has magnitude on average somewhat smaller than T^r_ϕ, while in the outer corona, on average it brings angular momentum toward the midplane.

Zoom In Zoom Out Reset image size

The majority of the net flux on the black hole is brought there by the coronal mechanism during the initial 5000M of time. However, the flux on the horizon does increase over the last 1.5 × 10⁴M, and, while the coronal mechanism continues to operate during this interval, it is possible that part of the flux might be brought in through the accretion disk itself. As the right panel of Figure 9 makes plain, the local magnetic inflow rate V_Φ is roughly 2 orders of magnitude smaller in the midplane than it is at the funnel edge. Nonetheless, its magnitude is still consistent with the total flux inflow rate during the latter part of the simulation. It may therefore be possible for flux advection directly associated with the accretion flow in the disk body to account for some part of the flux accumulation during the quasi-steady accretion phase. In this section, we attempt to estimate the size of that contribution.

We begin by comparing the rate of magnetic flux accretion to the rate of mass accretion. If magnetic flux and mass moved inward in lock-step, the functions

$$\begin{equation} {\cal M}(r,t) \equiv \int _{0}^{t} \, dt^{\prime } \dot{M}(r,t^\prime)/M_0 \end{equation} \tag{ 18 }$$

and

$$\begin{equation} {\cal A}(r,t) \equiv A_\phi (r,t,\pi /2)/A_\phi (r_{\rm max},0,\pi /2) \end{equation} \tag{ 19 }$$

would be identical. These quantities are normalized to the initial mass, M₀, and the initial flux; r_max is the outer radius of the initial torus. In Figure 13, we show the time dependence of both these quantities at two fiducial radii, just outside the event horizon (left panel) and at r = 20M (right panel). The former, of course, tells about the flux attached directly to the black hole and the mass deposited there. The difference between the two plots is the mass and flux within the accretion disk between the horizon and r = 20M; here we will refer to this region as the "inner disk."

Zoom In Zoom Out Reset image size

In evaluating Figure 13, one must allow for an important distinction between the histories of mass and magnetic flux: the mass of the black hole must increase monotonically (and the mass of the inner disk, although not required to do so, generally does), but the net magnetic flux can increase or decrease, both by magnetic reconnection and by the radial motion of closed loops. Thus, the inflow of magnetic flux is highly intermittent, and the time derivative of flux within some radius can easily be negative. On short timescales, the ratio of the rate of flux accretion to mass accretion can vary by at least an order of magnitude, as well as fluctuate in sign. For this reason, a sensible comparison of the rate of magnetic flux inflow relative to that of mass can be made only in regard to long-term time averages. In addition, in order to evaluate the efficiency of this process in real disks, we must restrict these averages to times when the inner disk was in approximate inflow equilibrium: for this simulation, that means the latter half of the data, from t = 10⁴M to t = 2 × 10⁴M. Within this span, we smooth ${\cal A}(t)$ by taking a running average 1000M wide and then compose the ratio of the change in ${\cal A}$ to ${\cal M}$ from 10⁴M to 2 × 10⁴M using the smoothed data. We find that $\Delta {\cal A}/\Delta {\cal M}\simeq 0.47$ at the horizon, and ≃0.50 at r = 20M. In other words, relative to the initial amount embedded in the accretion flow, magnetic flux moves inward roughly half as fast as mass, and these rates are the same at both the horizon and r = 20M. If the mass is in a state of inflow equilibrium, the equality of these two ratios shows that the magnetic flux in the inner disk is also in equilibrium.

This ratio of inflow rates may be contrasted with the local ratio of magnetic flux to mass. Looking beyond the short-timescale fluctuations in the plots of ${\cal A}(t)$ in Figure 13, one sees that the average values of the flux at the horizon and the flux contained within r = 20M are very similar. In other words, the vertical magnetic flux through the inner disk, Φ(r = 20M), contributes very little to $\cal {A}$. Time-averaging over the last 10⁴M of data, we find that it is only ≃0.86% of the initial flux, and it appears that whatever flux there is is confined to the outer part of the inner disk. By contrast, averaged over the same period, 5% of the initial mass can be found in the inner disk. Thus, in the inner disk, the time-averaged flux/mass ratio is only 17% of the initial value.

To reconcile this flux/mass ratio with the accreted flux/mass ratio, there are only two possibilities: either the magnetic flux moves in at the same rate (or slower) than the mass and (at least) 2/3 of the flux reaching the horizon between 10⁴M and 2 × 10⁴M arrived there via a different route (most likely the coronal mechanism), or the magnetic flux on average moved inward three times faster than the mass. We view the former alternative as much more plausible. Note, however, any flux delivery from outside the inner disk must be balanced by flux losses from the inner disk to the horizon because the magnetic inflow rate at r = 20M is nearly the same as that at the horizon.

Next, we further explore the nature of the magnetic flux equilibrium in the inner disk, using two different views of this quantity. Figure 14 shows the distribution of the vertical flux function Φ(r) time averaged over the last 5000M of the simulation and normalized to the initial flux value. For comparison, we also plot in that figure both the initial flux distribution and the final mass distribution in the disk, both likewise normalized to the initial total. The radial derivative dΦ/dr shows the magnitude and sign of the vertical magnetic field piercing the equator. We see that, like the mass, vertical flux has spread away from its initial location, both in and out, but, as already mentioned, a much smaller part of the net magnetic flux than the mass resides in the inner disk.

Zoom In Zoom Out Reset image size

The time dependence of the magnetic flux in the inner disk over the last 10⁴M in time is illustrated in Figure 15. As this figure shows, not only is the magnetic flux contained within the inner disk a rather small part of the total, it fluctuates in time, frequently going negative. Little sign of any long-term trend can be seen, consistent with our contention that the inner disk magnetic flux has reached a statistical equilibrium. The frequent sign changes of the local net flux suggest that the poloidal magnetic field lines in the inner disk predominantly close within the inner disk.

Zoom In Zoom Out Reset image size

Moreover, the magnetic field corresponding to the net flux is a very small fraction of the total field within the disk. Figure 16 shows the azimuthally averaged value of ${\cal B}^\theta$ through the equator at the end time as a function of radius, along with the initial value. The MRI has generated considerable field of both signs throughout the disk. Local regions of the disk have a net vertical flux, but their contrasting signs lead to little contribution to the total disk flux. The global net flux is present only as a slight positive excess moving both inward and outward with time.

Zoom In Zoom Out Reset image size

To reach the state described by these figures, most of the large-scale magnetic connections to the disk matter must have been reconnected away. But this should come as no surprise because we have already seen that a significant part of the original flux has moved by reconnection from the disk to the horizon. The matter in the inner disk has little net flux because it has been largely bypassed through the coronal mechanism. Thus, while the MRI has created a large effective turbulent viscosity in the sense that considerable mass has accreted from the location of the initial torus, very little net flux has moved with it. In other words, within the disk, turbulence created by the MRI transports angular momentum rapidly and drives accretion, but very little flux moves with the accreting mass.

Lastly, we raise the question of whether in the very long-term the magnetic flux on the horizon would continue to grow. If it is regulated by a combination of the magnetic and gas pressure near r = 6M, and these are in an approximately time-steady state, one might expect that further flux accumulation on the horizon would have to stop. Figure 17 supports this idea. This figure shows that the magnetic pressure inside the funnel, the magnetic pressure in the inner accretion flow, and the gas pressure in the inner accretion flow are well correlated with one another in a temporal sense, even while all three change by more than two orders of magnitude over the course of the simulation. The time-dependent data show fluctuations that are so large, however, that the much slower accumulation of flux seen in Figure 13 is well below the noise, so our data on magnetic and gas pressures cannot answer any questions about long-term saturation of magnetic flux attached to the horizon.

Zoom In Zoom Out Reset image size

That the magnetic and gas pressures should be closely related is an old idea (Rees et al. 1982; MacDonald & Thorne 1982), but there is also a long history of controversy about whether the magnetic field pressure in the funnel should be closer in magnitude to the (poloidal) magnetic pressure or the gas pressure in the inner disk. For example, Ghosh & Abramowicz (1997) and Livio et al. (1999) argued that as the poloidal magnetic field strength in the disk is subthermal, the funnel field strength should be regulated by it rather than the gas pressure. In the simulation presented here, the total, i.e., poloidal plus toroidal, magnetic pressure in the disk is smaller than the gas pressure for the duration of the simulation (i.e., the field remains subthermal), and the magnetic pressure from the field in the funnel lies between these two candidate regulators, sometimes closer to one and sometimes closer to the other. Thus, the funnel field is consistently stronger than the poloidal field in the disk.

The relevance of the inner disk poloidal field to the strength of the funnel field seems to be rather limited based upon what we observe about the process carrying flux to the horizon. The field in the funnel is nearly radial, and its total intensity is determined by the poloidal flux that is delivered to the black hole over the course of the simulation. This flux delivery system is the coronal mechanism, and so it should not be surprising that the poloidal flux within the accretion disk itself plays at most a secondary role (e.g., by determining the rate of reconnection and accretion of flux loops). Instead, simple pressure matching in the inner regions of the accretion flow—be it gas, (total, mostly toroidal) magnetic, or radiation pressure—appears to be important in regulating how much flux may be attached to the black hole. We speculate that when the flux on the horizon has a total field intensity larger than would be consistent with the gas and magnetic pressure in the nearby accretion flow, the rate of reconnection along the funnel field boundary rises, so that the field approaching the black hole is entirely in closed loops and does not add to the net flux on the horizon.