MOMENTUM TRANSPORT FROM CURRENT-DRIVEN RECONNECTION IN ASTROPHYSICAL DISKS

As matter accretes from astrophysical disks onto a central body (such as protostars, neutron stars, and black holes), angular momentum is rapidly transported outward. This redistribution of angular momentum is far too rapid to be explained by collisional viscosity. It is believed that turbulence initiated and sustained by magnetohydrodynamic (MHD) plasma instabilities can be responsible for the enhancement of effective viscosity in accretion disks (Shakura & Sunyaev 1973; Pringle 1981). Disks appear to be hydrodynamically linearly stable according to the Rayleigh criterion. However, it is not clear whether through a nonlinear process a hydrodynamically stable flow can become turbulent, a topic that has been debated vigorously. Recently, it has been shown experimentally that hydrodynamic turbulence in astrophysically relevant flows at Reynolds numbers up to several million cannot provide the required effective viscosity (Ji et al. 2006).

The ineffectiveness of hydrodynamic turbulence in transporting momentum has motivated the search for other instabilities that could yield enhanced viscosity. A candidate for such instability is the magnetorotational instability (MRI). The MRI, originally studied by Velikhov (1959) and Chandrasekkar (1961), was analyzed for astrophysical disks by Balbus and Hawley in the early 1990s (Balbus & Hawley 1991). They showed that a differentially rotating disk in the presence of a weak magnetic field is MHD unstable. Thus, under appropriate conditions, the MRI can lead to turbulent momentum transport in disks (see, e.g., Balbus & Hawley 1998, and references therein). The momentum transport arises from Maxwell or Reynolds fluid stresses. In this paper, we investigate whether those stresses and momentum transport can be provided by a current-driven MHD instability, rather than through the flow-driven MRI. In particular, we consider in depth a tearing instability, which is a resistive MHD instability that causes magnetic reconnection and can persist in a strong magnetic field (in regions in which the MRI could be stabilized by the strong field).

Magnetic fields have long been considered to be prevalent in some disk regions as well as to play a key role in jets emanating from disk regions. Observational evidence is in hand for the existence of strong magnetic fields in disks, although the observations are limited. In recent years, direct measurements have shown strong evidence of a magnetic field in protostellar disks. Existence of a significant azimuthal magnetic field of about 1 kG has been confirmed by observations around FU Orionis in the innermost regions of its accretion disk (Donati et al. 2005). The complex topology of magnetic fields around stars during their formation has also been confirmed through observation and extrapolation of the reconstructed surface magnetic field. Recent observation and analysis of magnetic fields on the surface of classical T Tauri stars (cTTSs) suggest a strong octopolar field (≈1.2 kG) and a smaller dipolar field (≈0.35 kG) (Donati et al. 2007).

There is some indirect evidence of a strong magnetic field in the protosolar nebula based on the magnetization of primitive meteorites. It was suggested that an MHD nebular dynamo can generate a magnetic field with an intensity of up to 10 G at the distance of a few astronomical units from the Sun (Levy 1978), and it was shown that this magnetic field can produce flares in a nebular corona (Levy & Araki 1989), which may account for the magnetized meteorites.

The source of large-scale magnetic fields of cTTS could be the combination of fossil fields (from interstellar media) and fields generated through dynamo action (Donati et al. 2007). However, the source of large-scale magnetic fields in disks is still unknown. Large-scale magnetic fields in disks may be supplied externally and may also be generated through a kinematic alpha–omega dynamo (Reyes-Ruiz & Stepinski 1999; Torkelsson & Brandenburg 1994; Rudiger et al. 1995). MRI can also produce a turbulent small-scale MHD dynamo (Brandenburg et al. 1995; Hawley et al. 2001); however, the generation of a large-scale magnetic field by MRI is not known to occur. Based on the observational evidence of a magnetic field and the modeling of a self-generated large-scale magnetic field (MHD dynamo) in disks, we examine current-driven instabilities as a candidate for angular momentum transport in an MRI stable magnetized disk with large currents.

In this paper, we explore the possibility of angular momentum transport in disks by current-driven instabilities in situations where MRI is stable (due to a strong magnetic field) and cannot provide the required effective viscosity. We recognize three classes of situations in which current-driven instabilities may be important in transporting angular momentum in astrophysical disks: (1) in the inner disk region around a young magnetized protostar, (2) in the innermost and active regions of weakly ionized protostellar and protoplanetary disks, and (3) in the upper and lower surfaces layers of a protoplanetary disk and magnetized coronas. Below we further discuss these situations.

The strong and complex large-scale magnetic fields detected in protostellar disks can provide the free energy (current) for current-driven instability and magnetic reconnection in these systems. Disk–star interaction and accreting matter to the stars along the complex magnetic field lines have been studied through MHD simulations by Long et al. (2008, and reference therein). However, with open and complex field lines attached to the disk observed at radii only up to several stellar radii and with disk differential rotation, magnetic reconnection may be important in the inner region of a disk close to a young protostar (such as a cTTS) and cause radial angular momentum transport.

Current-driven instabilities may also be important in protoplanetary and protostellar disks (e.g., protosolar nebulas, disks around cTTSs). However, in weakly ionized protoplanetary disks, the magnetic field might not be well coupled to the gas throughout the whole disk. Therefore, magnetic reconnection can only play a role in some parts of these disks which are thermally ionized (the innermost region r < 0.1 AU) or nonthermally ionized by cosmic rays (in the outer active region). Other possible candidates for momentum transport in disks are dynamo modes and flow-driven MRI. The notion of angular momentum transport through magnetic stresses from the oscillating dynamo modes has been examined by Stepinski & Levy (1990) in a protostellar nebula in the absence of flow-driven MRI turbulence. It was estimated that a magnetic field strength of 10²–10³ G is needed for transport of angular momentum from a protostar into a surrounding disk. Torkelsson & Brandenburg (1994) also found stronger magnetic torques from quadrupolar dynamo solutions. The effectiveness of MRI in transporting angular momentum transport in weakly ionized disks is limited and strongly depends on the ionization factor. It has been shown that in the so-called dead zone between the active layers in protoplanetary disks, MRI cannot operate (Gammie 1996; Fleming & Stone 2003; Zhu et al. 2010). However, including nonideal effects such as the Hall effect and ambipolar diffusion has been shown to change the stability and the saturation of MRI in weakly ionized protoplanetary disks (Wardle 1999, 2007; Balbus & Terquem 2001).

In strongly magnetized coronas where MRI is stable, the role of current-driven instabilities can become significant. Magnetic activity at the surface of the disks, possibly originated from the internal dynamos, may produce coronas, which can transport angular momentum vertically and radially. Using local simulations, Miller & Stone (2000) showed the generation of strongly magnetized corona (force-free) through an MRI-driven turbulence originated in the core of a disk with initially weak magnetic field. They also investigated angular momentum transport from MRI in a vertically extended local domain. However, the possibility of angular momentum transport in the strongly magnetized, MRI stable force-free corona was not studied. Theoretical models have also been developed to study the possibility of angular momentum transport by the coronal magnetic field (Goodman 2003; Uzdensky & Goodman 2008). The current-driven instabilities can also become important in corona, upper and lower layers of protoplanetary disks, with significant currents, where MRI neither operates in the core region (dead zone) nor is unstable in the coronas.

In this paper, we specify the equilibrium magnetic field (and corresponding current density) in disk geometry, which is chosen to be sufficiently strong that the flow-driven MRI is stable. We then study current-driven instabilities. Unlike flow-driven MRI and pressure-driven interchange instabilities, which can be treated locally, current-driven instabilities are driven due to equilibrium gradients with scale length comparable to the global scale lengths. It is shown that gradients in global equilibrium parallel current provide the driving source for both ideal and resistive current-driven (tearing) instabilities. Ideal current-driven instabilities are global modes (with global velocity perturbations) with fast Alfvénic growth rates which are independent of resistivity. However, for tearing modes, resistivity is important in the vicinity of the radius at which the wavenumber parallel to the magnetic field vanishes (i.e., k_|| = ${\bm k}$ · ${\bm B}/B$ = 0). Magnetic reconnection (or tearing) occurs at such radial locations. Although the reconnection occurs locally in radius, the tearing modes are global in their radial extent. Therefore, tearing modes are usually treated using asymptotic matching where solutions outside the reconnecting layer (outer layer or ideal η = 0 solutions) are matched with the inner layer (reconnecting layer η ≠ 0) solutions.

We investigate the physics of momentum transport through three sets of calculations of increasing completeness. First, we perform linear stability analysis to establish conditions for instability and to investigate radial structure of the modes. Second, through quasilinear analytic theory, we calculate the Maxwell stress in the outer ideal region and show its dependence on the azimuthal flow shear and instability growth rate. We find that in the presence of azimuthal Keplerian flow the radial and azimuthal magnetic fluctuations are in phase and cause a nonzero Maxwell stress for a single (a single spatial Fourier component) current-driven mode and thus momentum transport. We also examine the structure of the stress for a single mode in the quasilinear regime. We show that the stress and the resulting momentum transport are localized to the vicinity of the reconnection layer. However, an ideal current-driven instability has a global Maxwell stress, and therefore produces more global transport.

Third, we compute the simultaneous nonlinear evolution of multiple modes, each corresponding to different reconnection layers (separated radially). With the inclusion of nonlinear mode coupling, the Maxwell stress becomes global, extending over the full plasma cross-section. We perform nonlinear computations which are (1) weakly nonlinear driven and (2) strongly nonlinear driven. In the weakly nonlinear driven case, many current-driven modes are linearly unstable while the effect of nonlinear coupling to stable modes is weak. In the strongly nonlinear case, however, many stable modes are nonlinearly driven by linearly unstable modes. It is shown that in the weakly nonlinear driven case, the nonlinear structure of Maxwell stress is global mainly due to the growth of an ideal current-driven instability. However, in the strongly nonlinear case, the global structure of the Maxwell stress is mainly due to the nonlinear coupling of many tearing modes. We also estimate the effectiveness of the current-driven instabilities for angular momentum transport using the standard Shakura–Sunyaev α model for both the weakly nonlinear driven and strongly nonlinear driven cases. We find that Maxwell stress is much larger for the strongly nonlinear case compared to the weakly nonlinear case and α_SS is about 10⁻³ and 10⁻² for the weakly nonlinear driven and strongly nonlinear driven cases, respectively.

The paper is organized as follows. The characteristics of momentum transport by current-driven instabilities and the model used are described in Sections 2 and 3, respectively. Single-mode calculations are presented in Section 4. In Section 4.1, equilibrium and linear stability analysis are described. Analytical quasilinear calculation of Maxwell stress and the structure of a single mode in the quasilinear regime are presented in Section 4.2. We present full nonlinear multiple-mode computations in Section 5. Multiple computations for a weakly nonlinear driven and a strongly nonlinear driven case are presented in Sections 5.1 and 5.2, respectively. We summarize in Section 6.

The momentum transport caused by current-driven instability has the same origin in fluctuation-induced stresses as does the well-studied flow-driven MRI. In that sense, the mechanism for transport is the same, although the underlying MHD instability is different in its energy source and spatial structure. The change in flow is given by the MHD momentum equation

$$\begin{equation} \rho \frac{\partial \textbf V }{ \partial t } = - \rho \textbf V\cdot \nabla \textbf V + \textbf J \times \textbf B. \end{equation} \tag{ 1 }$$

For a rotating, cylindrical plasma we decompose the flow and magnetic fields into spatially mean and fluctuating quantities, where mean values (denoted by "〈〉") are averaged over axial and azimuthal directions and fluctuating values subtract out the mean. If we then average the above equation over the axial and azimuthal directions we find

$$\begin{equation} \rho \frac{\partial \langle \textbf V\rangle }{ \partial t } = - \rho \langle \textbf V \cdot \nabla \textbf V\rangle + \langle \textbf J \times \textbf B\rangle, \end{equation} \tag{ 2 }$$

which describes the time evolution of the mean flow which depends upon radius and time. We see that the mean flow can evolve from fluctuating flow (arising from a Reynolds stress, the first term on the right-hand side (RHS)) and fluctuating magnetic field (the Maxwell stress or Lorentz force of the second term on the RHS). These terms, quadratic in the fluctuations, can be nonzero upon spatial averaging. MHD current-driven instabilities produce spatially fluctuating quantities that, through the Reynolds and Maxwell stress, alter the radial profile of the axial and azimuthal mean flow. These forces can decrease the flow at one radius, while increasing it another, essentially transporting momentum radially. In this sense, the mechanism is identical to momentum transport from the well-studied MRI. Crucial differences lie in the energy source, spatial structure, and dynamics of the underlying instabilities.

The MRI is a global ideal instability (occurs at zero resistivity) with the stresses acting globally in radius. The current-driven tearing instability also has an amplitude that is global in radius, but the forces on the RHS of Equation (2) act locally in radius for a single instability. For a single tearing instability, magnetic reconnection occurs at one radial location—that for which the wavelength parallel to the equilibrium magnetic field becomes infinite (Biskamp 1993). That is, the instability amplitude is constant along the helical magnetic field at this radius, although it varies in the azimuthal and axial directions. Only in the vicinity of this radius is electrical resistivity important. The resistive layer of narrow radial extent, within which reconnection occurs, functions as a boundary layer between two ideal regions. Within this region amplitudes of velocity and current density can become large, leading to locally strong forces that alter flow.

An additional difference between the MRI and tearing instability is that momentum transport in the MRI can be understood as the instability acting to diminish its energy source (flow gradient). For the tearing instability, the analog is the reduction of the current density, which is a strong effect of the instability. Reduction of the flow gradient can be viewed as a parasitic effect of the tearing mode. Interestingly, the MRI also alters the magnetic field (Ebrahimi et al. 2009), so both instabilities alter the flow and field. We should note that an ideal current-driven instability, despite having the same energy source as tearing instability, exhibits different characteristics. A major difference is that for an ideal current-driven instability (with no reconnecting layer in the plasma volume), the fluctuation-induced forces are global as will be shown below.

The tearing mode, being a resistive instability, does not lend itself to simple analytic calculation of instability-induced transport. However, in Section 4.2, an analytic calculation of the transport is provided for the ideal region. We find that in the presence of azimuthal flow, the radial and azimuthal magnetic fluctuations are in phase and result in a nonzero Maxwell stress. The resulting stress (and the direction of angular momentum transport) depends on the global equilibrium as well as the global radial mode structure. We therefore employ quasilinear computations to show the structure of stresses and outward momentum transport. It will be shown that the structures of the stresses for ideal current-driven and tearing modes are very different, therefore they can affect the momentum transport differently. With multiple tearing modes present, nonlinear three-wave interactions alter the structure of the modes and the forces broaden radially.

Throughout this work, we employ the MHD equations in doubly periodic (r, ϕ, z) cylindrical geometry for both the analytic and computational studies. All variables are decomposed as $f(r,\phi,z,t)=\sum _{(m,k)} \widehat{f}_{m,k}(r,t) e^{i(m \phi + kz)}=\langle f(r,t)\rangle +\widetilde{f}(r,\phi \,z,t)$, where 〈f〉 is the mean (m = k = 0) component, and $\widetilde{f}$ is the fluctuating component (i.e., all other Fourier components with m ≠ 0 and k ≠ 0). Note that the mean component is a function of radius. We consider an azimuthal equilibrium flow $\textbf V_0 = V_{\phi }(r) \hat{\phi }$ in a current-carrying disk configuration plasma. (Here, V_ϕ(r) is the equilibrium value of the mean azimuthal flow 〈V_ϕ〉.) In order to excite the current-driven instabilities, both the vertical magnetic field and the azimuthal magnetic field, $\textbf B = B_z(r) \hat{z} + B_{\phi }(r) \hat{\phi }$, are imposed (Figure 1). The single fluid MHD equations are

$$\begin{eqnarray} \frac{\partial \textbf A }{ \partial t } &=& -\textbf {E} = S\textbf V\times \textbf B - \eta \textbf J,\\ \rho \frac{\partial \textbf V }{ \partial t } &=& -S \rho \textbf V \cdot \nabla \textbf V + S\textbf J \times \textbf B +P_m \nabla ^2 \textbf V \nonumber\\ &&-S \frac{\beta _0}{2}\nabla P -S \rho \nabla \Phi,\\ \frac{\partial P }{ \partial t } &=& -S\nabla \cdot (P \textbf V) - S (\Gamma -1) P \nabla \cdot \textbf V,\\ \frac{\partial \rho }{ \partial t } &=& -S\nabla \cdot (\rho \textbf V),\\ \textbf B &=& \nabla \times \textbf A,\\ \textbf J &=& \nabla \times \textbf B, \end{eqnarray} \tag{ 3 }$$

where the variables, ρ, P, V, B, J, Γ, and Φ are the density, pressure, velocity, magnetic field, current, ratio of the specific heats, and gravitational potential, respectively. Time and radius are normalized to the resistive diffusion time τ_R = 4πr²₂/c²η₀ and the outer radius r₂, making the normalized outer radius unity. Velocity V and magnetic field B are normalized to the Alfvén velocity V_A, and the magnetic field on axis B₀, respectively. The parameter $S = \frac{\tau _R}{\tau _A}$ is the Lundquist number (where τ_A = r₂/V_A), and $P_m= \frac{\tau _R}{\tau _{{\rm vis}}}$ measures the ratio of characteristic viscosity ν₀ to resistivity η₀ (the magnetic Prandtl number), where τ_vis = 4πr²₂/c²ν₀ is the viscous diffusion time. The factor β₀ = 8πP₀/B²₀ is the beta normalized to the axis value. The resistivity profile η is uniform. The boundary conditions in the radial direction are as appropriate to dissipative MHD with a perfectly conducting boundary: the tangential electric field, the normal component of the magnetic field, and the normal component of the velocity vanish, and the tangential component of the velocity is the rotational velocity of the wall. The azimuthal (ϕ) and axial (z) directions are periodic.

Zoom In Zoom Out Reset image size

We pose an initial value problem that consists of the equilibrium plus a perturbation of the form $\tilde{f}(r, \phi,z,t)=\tilde{f}_{m,k} (r,t) \mathrm{exp}(im\phi +ikz)$. Equations (1)–(6) are then integrated forward in time using the DEBS code. The DEBS code uses a finite difference method with a staggered grid for radial discretization and pseudospectral method for azimuthal and vertical coordinates. In this decomposition, each mode satisfies a separate equation of the form $\partial \tilde{f}_{m,k} / \partial t = L_{m,k} \tilde{f}_{m,k} +\sum _{(m^{^{\prime }},k^{^{\prime }})} N_{m,k,m^{^{\prime }},k^{^{\prime }}}$, where L_{m, k} is a linear operator that depends on $\tilde{f}_{0,0}(r,t)$, and $N_{m,k,m^{^{\prime }},k^{^{\prime }}}$ is a nonlinear term that represents the coupling of the mode (m, k) to all other modes (m', k'). (This latter term is evaluated pseudospectrally.) The time advance is a combination of the leapfrog and semi-implicit methods (Schnack et al. 1987).

We present three types of initial value computations: linear, nonlinear single mode, and fully nonlinear computations. In linear computations, the initial conditions consist of an equilibrium 〈f(r)〉 plus a single-mode $\tilde{f}_{m,k}(r,0) \mathrm{exp}(im\phi +ikz)$ perturbation. Only the mode (m, k) is then evolved; in particular, the equilibrium (the m = 0, k = 0 mode, or 〈f〉) is not evolved, and remains fixed in time. In the nonlinear single-mode computations, however, the m = 0, k = 0 component (the background) is allowed to evolve self-consistently. The evolution of the background profile $\tilde{f}_{0,0}$ can affect the evolution of the mode (m, k) and cause the mode to saturate. In a fully nonlinear computation, all modes are initialized with small random amplitude and are evolved in time, including the full nonlinear term ($N_{m,k,m^{^{\prime }},k^{^{\prime }}}$).

The plasma rotates azimuthally with a mean Keplerian flow 〈V_ϕ(r)〉∝r^−1/2. The initial (i.e., at t = 0) radial equilibrium force balance (Equation (2)) is satisfied by $\frac{\beta _0}{2}\nabla p + \rho \nabla \Phi = \rho V_{\phi }^2/r$, where ∇Φ = GM/r², and a magnetic force-free condition $\textbf {J}\times \textbf {B} =0$. The initial pressure and density profiles are assumed to be radially uniform. Pressure and density are evolved; however, they remain fairly uniform during the computations.

We consider a cylindrical disk-shaped plasma with aspect ratio L/(r₂ − r₁) (where L is the vertical height, and r₁ and r₂ are the inner and outer radii, respectively; Figure 1). The inner and outer radial boundaries are perfectly conducting, concentric cylinders that can rotate independently at specified rates. Periodic boundary conditions are used in the vertical and azimuthal directions.

The aspect ratio (L/(r₂ − r₁)) used in the nonlinear computations is 1.3. The nonlinear computations are performed in a thick-disk approximation where vertical and radial distances are of the same order. The range of parameters used in the computations is β = 1–10, S = 10⁴, P_m = 0.1–20.

Here, we present linear and quasilinear single-mode computations for current-driven instabilities. We first choose equilibria which are unstable for current-driven instabilities and then perform linear computations (Section 4.1). To study whether a single current-driven instability transports momentum, we examine the Maxwell stress in the ideal MHD region (outer region) through quasilinear analytical calculations (Section 4.2). We also examine the structure of quasilinear stresses for a single mode (with specific m and n) with the equilibrium described in Section 4.1.

Unlike flow-driven MRI and pressure-driven interchange instabilities, which can be treated locally, current-driven instabilities have global characteristics. Gradients in global equilibrium parallel current provide the driving source for the ideal and resistive current-driven instabilities. Ideal current-driven instabilities, the so-called kink modes, are helical long wavelength structures with global (broad radial extent) velocity perturbations. These ideal modes do not scale with resistivity and have fast Alfvénic growth rates. In an ideal plasma, the restoring force is infinite and fluid is frozen to the field. However, in the presence of resistivity, the field lines can break up and at locations where the parallel wavenumber is zero (the so-called resistive or reconnecting layer), field lines can reconnect. Both structure and growth rate of a resistive current-driven mode (tearing mode) are affected by resistivity. The solutions around the reconnecting layer are solved by including the resistivity (the so-called inner layer solutions) and are asymptotically matched with solutions outside the reconnecting layer (the so-called outer layer solutions). We will show that the structures of the stresses for ideal current-driven and tearing modes are very different, therefore they can affect the momentum transport differently.

4.1. Linear Computations

We first choose equilibrium profiles that are unstable for current-driven tearing instabilities. We start with a force-free plasma $\textbf {J}\times \textbf {B} =0$. Current flows parallel to the magnetic field line J_|| = λ(r)B, where J_|| is the parallel current. Figure 2(a) shows a typical equilibrium λ(r) = J_||/B profile used in the computations at t = 0. The equilibrium vertical and azimuthal magnetic field components can be obtained from the force balance equation $\nabla \times \textbf {B} = \lambda (r) \textbf {B}$, and are shown in Figure 2(a). The magnitude of the magnetic fields is large enough to make this equilibrium MRI stable.

Zoom In Zoom Out Reset image size

As described in the previous section, perturbations are assumed in the form of $\mathrm{exp}(i k \textbf r) \sim \mathrm{exp} (im\phi + i 2 \pi n z/L)$ in cylindrical geometry, where m and n are the azimuthal and axial mode numbers, respectively (k_z = 2πn/L). At locations where parallel wavenumber is zero $\textbf {k} \cdot \textbf {B} = m B_{\phi }/r + 2 \pi n B_z/L=0$ (the so-called resonant or reconnecting surfaces), resistivity becomes important and reconnection can occur. From this resonant condition, a field line winding number q = −m/n = 2πrB_z/LB_ϕ can be defined. The profile of axial winding number q is shown in Figure 2. Here, conventionally positive q corresponds to a negative axial mode (n) and vice versa. It will be shown that, depending on the equilibrium profile (λ), current-driven modes with −5 ⩽ q ⩽ 5 can be linearly unstable.

To verify that the equilibrium is tearing mode unstable, we first perform linear single-mode computations (single m and n) in the absence of azimuthal flow. The number of radial mesh points used is 250. The modal structure and the growth rates of the tearing modes depend on the equilibrium profile, i.e., parallel current profile λ(r). Therefore, we perform four sets of linear computations with four different equilibria. The four equilibria are distinguished by the maximum value of the equilibrium λ = J_||/B profile. We find that both resistive current-driven modes (tearing modes) and ideal current-driven modes with azimuthal mode numbers m = 0, 1, 2, 3 are linearly unstable for sets of equilibria with λ_max = 10.6 (shown in Figure 2), λ_max = 14.2, and λ_max = 17.2. For equilibria with λ_max = 9, only m = 0, 1, 2 tearing modes are linearly unstable. Below, we further discuss the mode structure and the properties of tearing modes and ideal current-driven modes.

In the presence of resistivity, current-driven instabilities (tearing) can become linearly unstable. The radial structures of radial magnetic and velocity eigenfunctions for m = 1 tearing mode with axial mode number n = −1 are shown in Figure 3(a) (for λ_max = 10.6) and Figure 3(b) (for λ_max = 17.2). As can be seen, the reconnecting component of magnetic field ($\widetilde{B}_r$) is nonzero around the reconnecting surface (r = 0.6 for (1, −1) in Figure 3(a)). Moreover, the radial fluctuating velocity is concentrated and changes sign around the reconnecting surface. One of the characteristics of tearing instability is the jump in the logarithmic derivative of $\widetilde{B}_r$ across the resistive layer ($\Delta ^{\prime } = (\widetilde{B}_r^{\prime }|_{\rm{rs}}^+ - \widetilde{B}_r^{\prime }|_{\rm{rs}}^-)/\widetilde{B}_r|_{\rm{rs}}$), where rs denotes reconnecting surface. For the tearing mode to be unstable, this jump should be positive (Furth et al. 1963). The growth rate of the tearing mode scales with resistivity and Δ' as γ_tearing∝η^3/5Δ'^4/5 (or ∝S^−3/5Δ'^4/5). As is seen in Figure 3, Δ' is positive for m = 1 tearing modes, and is larger in Figure 3(b) with larger free energy (current gradient λ_max), which results in larger growth rate. We also note that for larger λ_max, the (1, −1) mode becomes a double tearing mode (reconnects at two locations), and the radial velocity perturbation changes sign twice in radius (absolute values of the eigenfunctions are shown in Figure 3(b)).

Zoom In Zoom Out Reset image size

The growth rates of m = 0–3 tearing modes for four sets of λ profiles with λ_max = 9, 10.6, 14.2, and 17.2 are given in Table 1. For smaller λ_max, the tearing mode growth rates are smaller, and the growth rates increase with λ_max approaching the ideally unstable limit (approaching Alfvénic growth rates). Modes with azimuthal mode numbers 0–3 and axial mode numbers 1, −1 have tearing mode structure with a corresponding reconnecting surface. However, the mode with axial mode number zero (1, 0) is ideally unstable for equilibria with λ_max = 10.6, 14.2, and 17.2. The radial structure of this mode is shown in Figure 4 which has a kink-like characteristic. Unlike the tearing modes, ideal current-driven modes do not necessarily have to have a reconnecting surface (with $\textbf {k} \cdot \textbf {B}$) within the plasma and the fluctuating velocity components can be very global without changing sign (kink-like). The growth rates of the ideal current-driven modes also do not scale with resistivity. It has been shown that the current-driven modes are unstable for the equilibria chosen (Figure 2), and therefore we can study the effect of both tearing and ideal current-driven modes on the momentum transport through nonlinear computations (Section 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Growth Rates (γτ_A) of Tearing Modes (m, n), for Four λ_max

λmax	(1, 0)	(1, 1)	(1, −1)	(0, 1)	(2, 1)	(2, −1)	(3, 1)	(3, −1)
9	Stable	0.030	0.031	0.0362	0.0155	0.0165	Stable	Stable
10.63	0.058	0.059	0.045	0.057	0.048	0.0445	0.029	Stable
14.2	0.23	0.148	0.101	0.132	0.141	0.0557	0.116	0.0105
17.2	0.3	0.31	0.21	0.275	0.3	0.14	0.255	0.082

Download table as:  ASCIITypeset image

4.2. Analytical Quasilinear Calculations

In order to investigate momentum transport from current-driven instabilities, we obtain MHD stresses using quasilinear calculations of a single mode with an initial azimuthal Keplerian flow. The question is whether a single current-driven mode can transport momentum and affect the azimuthal flow profile. To obtain more insight into this question, we first analytically examine the ideal MHD equations with azimuthal flow in cylindrical geometry. We then construct the Maxwell stress $\langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle = 2 {\rm Re}(r^2 \widetilde{B}_{\phi }^*\widetilde{B}_r)$ from the linearized solutions in the outer (ideal) region. Through analytical quasilinear calculations, we aim to identify whether the Maxwell stress is nonzero for a current-driven instability and therefore momentum can be transported. We do not intend to analytically solve the complete sets of solutions (inner and outer solutions), and we only present the simplified outer solutions. To construct the radial structure of the quasilinear stresses, we use the eigenfunctions from the linearized computations.

The linearized incompressible ideal MHD equations in the presence of mean azimuthal flow are

$$\begin{eqnarray} &&\rho \left(\frac{\partial \widetilde{\textbf {V}}}{\partial t}+ (\textbf V_0 \cdot \nabla) \widetilde{\textbf {V}} + (\widetilde{\textbf {V}} \cdot \nabla) \textbf V_0\right) \nonumber\\ &&\quad= -\nabla (\textbf B_0 \cdot \widetilde{\textbf {B}}) +(\textbf B_0 \cdot \nabla) \widetilde{\textbf {B}}+(\widetilde{\textbf {B}} \cdot \nabla) \textbf B_0, \end{eqnarray} \tag{ 9 }$$

$$\begin{equation} \frac{\partial \widetilde{\textbf {B}}}{\partial t}= \nabla \times (\widetilde{\textbf {V}} \times \textbf B_0) + \nabla \times (\textbf V_0 \times \widetilde{\textbf {B}}). \end{equation} \tag{ 10 }$$

We assume perturbation of the form Q(r, ϕ, z, t) = Q(r) × exp(iωt + i(mϕ + kz)), and the equilibrium magnetic field and flow are $\textbf {B} = B_{\phi }(r) \textbf {e}_{\phi } + B_z(r) \textbf {e}_z$ and $\textbf {V}_0 = V_{\phi }(r) \textbf {e}_{\phi }$, respectively. The linearized equations can be combined and be presented in the form of one ordinary differential equation for $\widetilde{V}_r$,

$$\begin{equation} (r \widetilde{V}_r)^{\prime \prime } + a_1 (r \widetilde{V}_r)^{\prime } + a_2 (r \widetilde{V}_r)=0, \end{equation} \tag{ 11 }$$

where

$$\begin{eqnarray} a_1&= \frac{1}{A} \frac{dA}{dr} - \frac{m G}{r \bar{\omega } (1-M^2)} + \frac{2 m \rho \bar{\omega } V_{\phi }}{r^2 F^2 (1-M^2)}\nonumber\\ a_2&=\frac{1}{A} \bigg[\frac{d}{dr} \bigg(\frac{m G A }{r \bar{\omega } (1-M^2)}\bigg) - \frac{d}{dr} \bigg(\frac{2 m B_{\phi } F}{r \bar{\omega } (m^2 + k^2 r^2)}\bigg) \nonumber\\ &+ \frac{2 B_{\phi }}{r \bar{\omega }} \bigg(\frac{B_{\phi }}{r}\bigg)^{\prime } + \frac{4 k^2 B_{\phi }^2}{r \bar{\omega } (1-M^2) (m^2 + k^2 r^2) } -\frac{F^2 (1-M^2)}{r \bar{\omega }}\nonumber\\ &- \bigg(\frac{2 B_{\phi } m^2 F}{r^2 \bar{\omega }^2(m^2 + k^2 r^2)}\bigg) \bigg(V_{\phi }^{\prime } -\frac{V_{\phi }}{r}\bigg)\nonumber\\ & - \bigg(\frac{4 k^2 F B_{\phi } M^2}{\bar{\omega }^2 (m^2 + k^2 r^2)(1-M^2)}\bigg)\frac{V_{\phi }}{r}\nonumber\\ & - \bigg(\frac{4 \rho k^2B_{\phi }}{F(m^2 + k^2 r^2) (1-M^2)}\bigg)\frac{V_{\phi }}{r}\nonumber\\ & - \frac{2 \rho k^2 V_{\phi } G}{\bar{\omega } (m^2 + k^2 r^2) (1-M^2)} \bigg] \end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray*} A= \frac{r F^2 (1-M^2)}{\bar{\omega } (m^2+k^2 r^2)},\quad && G = (1-M^2)V_{\phi }^{\prime } - (1+M^2)V_{\phi }/r \end{eqnarray*}$$

$$\begin{eqnarray*} \bar{\omega } = \omega +m V_{\phi }/r,\quad && M = \sqrt{\rho } \frac{\bar{\omega }}{F} \end{eqnarray*}$$

$$\begin{eqnarray*} F= \frac{m B_{\phi }}{r}+ k B_z,\quad && k_{\perp }B = - m B_z/r+ k B_{\phi }, \end{eqnarray*}$$

where ω = ω_r + iγ, and ω_r and γ are the real frequency and the growth rate, respectively. The equation for parallel magnetic field perturbations $\frac{\widetilde{\bf{B}} \cdot \bf{B}}{B} = \widetilde{B}_{||}$ in terms of radial velocity perturbation is obtained as

$$\begin{eqnarray} \widetilde{B}_{||}&=&i \frac{A}{B} \left[\frac{\partial }{\partial r}(r \widetilde{V}_r)\right.\nonumber\\ &&\left.-\left(\frac{2 m B_{\phi }}{(1-M^2)F r^2}+\frac{m G}{r \bar{\omega } (1-M^2)}\right) (r \widetilde{V}_r)\right]. \end{eqnarray} \tag{ 13 }$$

To obtain the Maxwell stress term, we write the parallel magnetic field perturbation (Equation (13)) in terms of the radial magnetic field perturbation $\widetilde{B}_r = F \widetilde{V}_r/\bar{\omega }$ and the parallel magnetic field perturbation without mean flow $\widetilde{B}_{||}^{(0)}$,

$$\begin{eqnarray} \widetilde{B}_{||}&=& (1-M^2) \widetilde{B}_{||}^{(0)} - i M^2 \left(\frac{2 m B_{\phi }}{r B (m^2+k^2 r^2)}\right) (r \widetilde{B}_r)\nonumber\\ && +i M^2\left(\frac{2 m F V_{\phi }}{r \bar{\omega } B (m^2+k^2 r^2)}\right)(r \widetilde{B}_r), \end{eqnarray} \tag{ 14 }$$

where $\widetilde{B}_{||}^{(0)}$$= i [F r (r \widetilde{B}_r)^{\prime } - r (k_{\perp }B) \lambda (r \widetilde{B}_r)] /(B (m^2 + k_z^2 r^2))$ is the parallel magnetic field in the absence of mean flow obtained from Equation (13), and also known from the Newcomb equation (Newcomb 1960), and λ = J_||/B. In the absence of mean flow, it can be seen that $\widetilde{B}_{||}^{(0)}$ is purely imaginary and $\widetilde{B}_r$ and $\widetilde{B}_{||}^{(0)}$ are out of phase, and therefore parallel and azimuthal Maxwell stresses vanish for the ideal region. However, with the azimuthal mean flow the Maxwell stress $\langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle = 2 {\rm Re}(r^2 \widetilde{B}_{\phi }^*\widetilde{B}_r)$ from a current-driven mode is nonzero and is obtained as

$$\begin{eqnarray} \langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle &=& 2 \gamma \left(2 k \rho r^2 \frac{(\omega _r +m V_{\phi }/r)}{k_{\perp } F^2} | \widetilde{B}_{||}^{(0)}| \widetilde{B}_r\right)\nonumber \\ &&+ 2 \gamma \left(4 k \rho r^2\frac{(\omega _r +m V_{\phi }/r)}{ k_{\perp } (m^2 + k^2 r^2)F^2 B} m B_{\phi } \widetilde{B}_r^2 \right)\nonumber\\ &&- 2 \gamma \left(2 k \rho r^2\frac{m V_{\phi }}{ k_{\perp } (m^2 + k^2 r^2)F B} \widetilde{B}_r^2 \right), \end{eqnarray} \tag{ 15 }$$

where

$$\begin{equation} \widetilde{B}_{\phi }=\frac{k \widetilde{B}_{||}}{k_{\perp }} - \frac{i B_{z}}{(Brk_{\perp })}\frac{\partial }{\partial r}(r \widetilde{B}_r). \end{equation} \tag{ 16 }$$

The azimuthal Maxwell stress (Equation (15)) has been calculated using Equations (14) and (16) assuming $\widetilde{B}_r$ is purely real. With mean azimuthal flow, Equation (15) shows that the joint effect of mode growth and mean flow produces a nonzero stress term in the outer region. We note that for an ideal current-driven mode without a reconnecting surface (F ≠ 0 everywhere), the stress term can be global. The sign and the structure of the resulting stress depend on the global equilibrium as well as the global radial mode structure and do not follow from Equation (15) alone. However, for a special case of an ideal current-driven mode m = 1, n = 0 (k = n/R = 0), discussed in Section 4.1, we can examine the sign of the Maxwell stress using the local WKB approximation. For this mode, the Maxwell stress is simplified as $\langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle = 2 {\rm Re}\lbrace [\frac{i r^2}{m}\frac{\partial }{\partial r}(r \widetilde{B}_r)]^*\widetilde{B}_{r}\rbrace$, where $\widetilde{B}_r$ is complex, and we have used Equation (16) for the toroidal magnetic perturbation $\widetilde{B}_{\phi }$ (or equivalently from $\nabla \cdot \textbf B =0$). Using a local approximation, a negative Maxwell stress $\langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle = -2 r^3 \frac{k_r}{m} \widetilde{B}_r^2$ is obtained, where k_r is the local radial wavenumber. The negative Maxwell stress results in a bidirectional Lorentz force and causes outward momentum transport. Although WKB is not a valid approximation for the current-driven instability and does not provide the structure of the transport, it demonstrates an outward transport, which is consistent with the global solutions given below.

To obtain the radial structure of stresses, we employ quasilinear computations. Using the linear eigenfunctions for $\widetilde{B}_r$ (as shown in Figure 4), the radial structure of the Maxwell stress term for an ideal current-driven mode, m = 1, n = 0, is shown in Figure 5. Similarly, the quasilinear Reynolds stress term, $\langle r^2 \widetilde{V}_{\phi } \widetilde{V}_r\rangle$, can also be constructed from the linearized solutions for the m = 1, n = 0 mode (as shown in Figure 5 by the dashed-dotted line). We note that the total stress term, $\langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle - \langle r^2 \widetilde{V}_{\phi } \widetilde{V}_r\rangle$, is negative everywhere, which causes a transport of angular momentum outward. The total azimuthal fluid force consisting of Lorentz and inertia terms $\langle \widetilde{\textbf {J}} \times \widetilde{\textbf {B}}\rangle _{\phi } - \rho \langle \widetilde{\textbf {V}} \cdot \nabla \widetilde{\textbf {V}}\rangle _{\phi } = \frac{1}{r^2}[\frac{d}{d r}\langle r^2 \widetilde{B}_{\phi } \widetilde{B}_r\rangle - \rho \frac{d}{d r}\langle r^2 \widetilde{V}_{\phi } \widetilde{V}_r\rangle ]$, is bidirectional and causes transport of momentum outward, as shown in Figure 5.

For a resistive current-driven mode (a tearing mode) since the growth rate scales as S^−3/5 or η^3/5, where S is the Lundquist number, the Maxwell stress term is small in the outer region (η = 0) and the main contribution arises from the inner layer solution (resistive region). Equation (15) only presents the outer solution and is singular around the reconnection layer (F = 0). To obtain the radial structure of stresses for a tearing mode with a reconnecting surface $\textbf k \cdot \textbf B= F =0$, the inner layer equations need to be solved (Ebrahimi et al. 2008). We therefore perform nonlinear single tearing mode computations with two sets of equilibria (λ = 10 and λ = 17). The radial structure of the quasilinear Maxwell and Reynolds stress terms during the linear phase for two cases is shown in Figure 6. The mode structure for the case with λ_max = 10 is a single-mode tearing structure. As discussed, the main contribution of MHD stresses for the tearing mode arises in the inner region, which is also confirmed by the computations (Figure 6(a)). As is seen in Figure 6(a), the MHD stresses are localized around the reconnecting surface (r ≈ 0.6), and are small in the outer ideal region. Because the computation is in the visco-resistive regime, the Maxwell stress is much larger than the Reynolds stress. The localization of the stresses leads to a localized flattening of the azimuthal flow within the resistive layer. It is also seen in Figure 6(a) that the azimuthal Lorentz force ($\langle \widetilde{J} \times \widetilde{B}\rangle _{\phi }$) is bidirectional and transports momentum outward. As the free energy for the instability increases, for the larger current (λ_max = 17) case, the double tearing mode becomes unstable and has two reconnecting surfaces at two radii (r ≈ 0.5, r ≈ 0.8). The radial structure of the Maxwell stress shown in Figure 6(b) is broader than the single tearing mode case.

Zoom In Zoom Out Reset image size

It should be mentioned that a more generalized set of ideal compressible equations including equilibrium flows in a cylindrical current-carrying plasma (Bondeson et al. 1987) and the gravitational force (Keppens et al. 2002; Blokland et al. 2005) has been obtained. They used the Frieman & Rotenberg (1960) formalism to obtain first-order differential equations for displacement ξ_r and Π (total kinetic and magnetic pressure perturbations). Using these generalized equations (Equations (14)–(20) in Blokland et al. 2005), after some algebra we have also obtained the parallel magnetic field perturbation ($\widetilde{\textbf {B}} \cdot \textbf {B}/B$) which reduces to Equation (13) in the incompressible limit. We therefore find that in the incompressible regime by ignoring the pressure perturbation (and with uniform equilibrium pressure), but including the gravitational force perturbation, the resulting Maxwell stress term remains the same as Equation (15).

We also note that in a current-free plasma with a uniform axial magnetic field and a differential rotation, Equations (11) and (12) reduce to a differential equation for a global axisymmetric m = 0 MRI (Velikhov 1959; Chandrasekkar 1961),

$$\begin{equation} \frac{\partial ^2 \widetilde{V}_{r}}{\partial r^2}+ \frac{1}{r}\frac{\partial \widetilde{V}_{r}}{\partial r} - \left[\frac{1}{r^2} + k^2 + \frac{2r \Omega ^{\prime }\Omega k^2}{(\omega _A^2-\gamma ^2)} - \frac{4\Omega ^2 \gamma ^2 k^2}{(\omega _A^2-\gamma ^2)^2}\right]\widetilde{V}_{r}=0, \end{equation} \tag{ 17 }$$

where Ω = V_ϕ(r)/r, $\textbf {B}=B_0 \textbf {e}_z$, and ω²_A = k²B₀²/ρ. It can be shown that the stresses from an axisymmetric MRI can be global and depend on flow shear (Ebrahimi et al. 2009).

Here, we investigate full nonlinear dynamics of momentum transport from current-driven instabilities in a disk geometry using multiple tearing mode computations. In the nonlinear computations, with multiple modes (in both azimuthal and vertical directions), the additional effect of nonlinear mode coupling is revealed. Moreover, the transfer of energy from fluctuations to the mean field can occur during the nonlinear saturation. We perform full nonlinear computations when the current-driven modes with different azimuthal and axial mode numbers are included in the computations. Current-driven modes can nonlinearly interact and cause nonlinear growth. Three sets of nonlinear computations are performed in a disk geometry. The first two multiple-mode computations are weakly nonlinear driven and discussed in Section 5.1. A multiple-mode computation, which is strongly nonlinear driven, is presented in Section 5.2. In the weakly nonlinear driven case, many current-driven modes are linearly unstable while the effect of nonlinear coupling to stable modes is weak. In the strongly nonlinear case, however, many stable modes are nonlinearly driven by linearly unstable modes.

5.1. Weakly Nonlinear Driven

Here, two sets of computations are performed in a disk geometry with similar force-free equilibrium shown in Figure 2 with (1) λ_max = 9 and (2) λ_max = 14.2. Both computations start with a Keplerian flow profile with an on-axis amplitude of V_ϕ/V_A = 0.8 (with P_m = 1, S = 10⁴, β₀ = 10, and radial, azimuthal, and axial resolutions n_r = 220, 0 < m < 21, and −43 < n < 43, respectively).

The computations start with a current-carrying equilibrium, and the free energy from the parallel current causes the current-driven instabilities to grow. The radial magnetic energy for different current-driven tearing modes is shown in Figure 7. Tearing modes which are linearly unstable (with the growth rates given in Table 1 for equilibria with λ_max = 9 and λ_max = 14.2), are also driven linearly in the nonlinear computations. Moreover, other modes with higher azimuthal mode numbers are driven nonlinearly. As shown in Figure 7(a), modes with tearing parity m = 0, 1, 2 (both n = 1 and n = −1) linearly start to grow and saturate around t = 0.015τ_R. For the tearing mode m = 1, n = 1, there is a second nonlinear growth before saturation. Modes with higher azimuthal mode numbers are also shown, and it is seen that the m = 3 mode linearly grows with small growth rate. This mode is linearly stable in the absence of flow and here with the azimuthal flow becomes linearly unstable. For the case with equilibrium λ_max = 9, the current gradient (free energy) is not large enough to make the ideal current-driven modes (including non-resonant) unstable. Thus for this equilibrium, only resonant resistive tearing modes are unstable.

Zoom In Zoom Out Reset image size

When the free energy, i.e., current gradient, increases both resistive and ideal current-driven modes become unstable. Figure 7(b) shows the radial magnetic energies for several current-driven modes with the equilibrium λ_max = 14.2. As seen, tearing modes m = 0, 1, 2, and 3 linearly grow as expected from linear stability analysis (Table 1). These tearing modes also have a nonlinear growth around t = 0.002τ_R. In addition to linearly unstable modes, the m = 4 tearing mode which is linearly stable, nonlinearly starts to grow around t = 0.002τ_R. For this equilibrium, as shown in Figure 7(b), the ideal current-driven mode m = 1, n = 0 is linearly unstable and saturates at about two order of magnitude higher amplitude.

As current-driven modes grow and nonlinearly saturate through modifying the source of instability, i.e., the current gradient, they also modify the mean azimuthal flow profile. Figure 8 shows the modification of the azimuthal flow profile during the nonlinear evolution for both cases λ_max = 9 and λ_max = 14.2. For the case with equilibrium λ_max = 9, the modification of flow profile is mainly due to the resistive tearing modes. The azimuthal flow profiles at two times t₁ = 0.012, before saturation, and t₂ = 0.023, after the saturation of tearing modes, are shown in Figure 8(a) with λ_max = 9. At t = t₁, the flow is modified around r = 0.6 due to the m = 1 tearing mode. As discussed before, the azimuthal flow is modified through the Maxwell and Reynolds stresses. The structure of the fluid stresses is shown in Figure 9. The Maxwell stress at t = t₁ is localized around r = 0.6 and causes the modification of the flow and transport of momentum outward. After the nonlinear saturation, the Maxwell stress from all the tearing modes becomes broader (Figure 9(b)), and thus causes a broader flow modification as shown at t = t₂ (Figure 8(a)). The structure of the total nonlinear Maxwell stress (from multiple tearing computation) is broader than the structure for a single tearing mode (Figure 9(a)). The structure is more concentrated in the plasma core (around r  =  0.5 to r  =  0.8), where the current gradient is large.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The momentum transport from current-driven modes becomes stronger when ideal modes are also present. The nonlinear evolution of the azimuthal flow profile for the case with larger current gradient λ_max = 14.2 is shown in Figure 8(b). The structure of Maxwell and Reynolds stress terms is also shown in Figure 10. As seen, the Maxwell stress is concentrated in the core region at t = t₁ before the nonlinear saturation (Figure 10(a)) and the flow modification also occurs around the region where the stresses are peaked from r = 0.4 to r = 0.7. Around the nonlinear saturation at t = t₂ and after the saturation at t = t₃, it is shown that fluid stresses become more global (Figures 10(b) and (c)) and cause a global modification of the azimuthal flow profile and momentum transport outward (Figure 8(b)). The global modification is due to both resistive tearing and ideal current-driven modes.

Zoom In Zoom Out Reset image size

It is also worth estimating the effectiveness of the current-driven instabilities for angular momentum transport using the standard Shakura–Sunyaev α model. In this model, the effective viscosity is parameterized by the quantity α_SS which is related to the Reynolds and Maxwell stress terms by $\langle \rho \widetilde{V}_r \widetilde{V}_{\phi }\rangle - \langle \widetilde{B}_r \widetilde{B}_{\phi }\rangle = \alpha _{\rm{SS}} \langle P\rangle$, where P is the average pressure. The normalized alpha used for our simulations is α_SS/β₀/2. For the two weakly nonlinear driven cases presented here, λ_max = 9 and λ_max = 14.2, the time-averaged Shakura–Sunyaev α, 〈α_SS〉, during nonlinear saturation are 〈α_SS〉 ≈ 5 × 10⁻⁶ and 〈α_SS〉 ≈ 4 × 10⁻⁴, respectively. We have also calculated the time-averaged Shakura–Sunyaev α for a weakly nonlinear case with λ_max = 14.2 and β₀ = 1, which is 〈α_SS〉 ≈ 2.5 × 10⁻³.

5.2. Strongly Nonlinear Driven

We further investigate the effect of nonlinear mode coupling in a different equilibrium setting. The computations are performed in a force-free equilibrium for which the components of mean magnetic fields are shown in Figure 11 (solid lines). Full nonlinear computation starts with a Keplerian flow profile with an on-axis amplitude of V_ϕ/V_A = 0.8 (with P_m = 5, S = 10⁴, β = 10, and radial, azimuthal, and axial resolutions n_r = 250, 0 < m < 11, and −43 < n < 43, respectively). The free energy, parallel current, is mainly concentrated in the inner-half plasma region and causes a current-driven instability with mode number m = 1, n = −1 to become linearly unstable. The radial magnetic energy for different current-driven tearing modes is shown in Figure 12. Almost all the modes are stable up to around t/τ_R = 0.004 and plasma is in a quasi-single-mode state (with m = 1, n = −1). Two other tearing modes, m = 1, n = −2 and m = 2, n = −1 also grow linearly with small growth rates. As shown in Figure 12(a), around t/τ_R = 0.004, due to nonlinear mode coupling the axisymmetric m = 0 mode becomes nonlinearly unstable and saturates at a large amplitude comparable to the amplitude of the initial m = 1 linearly driven mode. Other non-axisymmetric modes are also driven nonlinearly and a turbulent state is formed. The magnetic spectrum for m = 0–2 during the nonlinear state is shown in Figure 12(b); as can be seen, the spectrum shows a broad range of magnetic fluctuations. The transition from a quasi-single-mode to a multiple-mode state occurs due to strong mode–mode coupling and transfer of energy from m = 1 modes to m = 0 modes. In the weakly nonlinear driven case (Section 5.1), the equilibrium was linearly unstable for most of the current-driven modes, and the nonlinearly driven modes (m = 3 and m = 4 modes for the cases with λ_max = 9. and λ_max = 14.2, respectively) would not grow to large amplitudes. Here, most of the tearing modes are driven nonlinearly and saturate at large amplitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Tearing fluctuations during the nonlinear state (t/τ_R > 0.005) can affect the mean profiles through the fluctuation-induced convolution terms. Mean magnetic fields during the nonlinear saturation, which are affected by the magnetic fluctuations, are shown in Figure 11 (the dashed lines). As can be seen, the vertical magnetic field changes sign and toroidal flux is redistributed by the fluctuation-induced term, $\langle \widetilde{\bm V} \times \widetilde{\bm B}\rangle$, the so-called dynamo term. The tearing fluctuations also affect the flow profile and cause momentum transport. The radial structures of Maxwell and Reynolds stresses during the two states, quasi-single mode and multiple modes, are shown in Figure 13. During the quasi-single-mode state (t/τ_R < 0.004), the Maxwell stress transports momentum outward, but it is localized in the inner region around the reconnection layer. However, during the multiple-mode state (t/τ_R > 0.005) the structure of total stresses is very global and all the tearing modes contribute to the momentum transport. It can also be seen that the Maxwell stress is much stronger (about an order of magnitude) for the strongly nonlinear driven case compared to the weakly nonlinear driven case (Figure 10).

Zoom In Zoom Out Reset image size

We have also performed strongly nonlinear driven cases with the same force-free equilibrium shown in Figure 11 but with P_m = 20, β = 10, and P_m = 20, β = 1. The structure of stresses after nonlinear saturation (after the growth and saturation of the nonlinearly driven mode, m = 0) is shown in Figure 14. As seen, Maxwell stress is broad due to nonlinear mode coupling. We have calculated the time-averaged Shakura–Sunyaev α during the nonlinear saturation, which is 〈α_SS〉 ≈ 4.3 × 10⁻³ and 〈α_SS〉 ≈ 3.3 × 10⁻² for the two cases with β = 10 and β = 1, respectively.

Zoom In Zoom Out Reset image size

We have shown that tearing instabilities can transport angular momentum effectively in thick disks. A single tearing mode transports momentum locally, in the vicinity of the radial reconnection location. However, with multiple modes simultaneously unstable and nonlinearly coupled to each other, the resultant Maxwell stress and momentum transport can span the entire plasma cross-section. It is also shown that an ideal current-driven instability without a reconnecting surface has a global Maxwell stress and therefore can produce more global momentum transport. However, the global Maxwell stress from many tearing modes causes much stronger momentum transport compared to an ideal current-driven instability with a global stress. These conclusions are drawn from a hierarchy of calculations, including linear stability, quasilinear analytic calculations, and computation of the full nonlinear evolution of multiple modes.

These instabilities, although present in rotating plasmas, are driven magnetically by the plasma current. We have selected for study disk equilibria for which the magnetic field is sufficiently strong that the plasma is stable to the flow-driven MRI. Determination of the relevance of tearing instabilities to astrophysical disks requires extension of these studies to a wider range of conditions, such as to thin disks and different magnetic field strength. As the disk becomes thinner, the relevant tearing modes (with smaller wavelengths) may become more stable. The effectiveness of tearing modes in transporting momentum in disks also depends on the ratio of magnetic energy to the flow energy. Although in the linear regime, ideal current-driven instabilities are independent of resistivity, nonlinear momentum transport from both tearing and ideal current-driven instabilities in the low magnetic diffusivity regime could be an interesting topic for a future work.

The authors are grateful for useful discussions with Dalton Schnack, Ellen Zweibel, and Amitava Bhattacharjee.