MHD Instabilities in Accretion Disks and Their Implications in Driving Fast Magnetic Reconnection

As initial conditions, we have used a density profile obtained from a magnetostatic equilibrium in the vertical direction since we have adopted a strong azimuthal magnetic field to trigger the PRTI. The density profile is given by (see Johansen & Levin 2008)

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}\exp [-{z}^{2}/2{H}_{\beta }^{2}],\end{eqnarray} \tag{ 8 }$$

where ${\rho }_{0}=1$ is the initial density in the midplane of the accretion disk, ${H}_{\beta }=\sqrt{1+{\beta }_{0}^{-1}}{c}_{s}/{{\rm{\Omega }}}_{0}$ is the scale height of the gas, and β₀ is the initial ratio between the thermal and magnetic pressures. For consistency with the work of Johansen & Levin (2008), we have adopted ${{\rm{\Omega }}}_{0}=1$ and the thermal scale height as our unit of length ($H={c}_{s}/{{\rm{\Omega }}}_{0}=1$), which yields ${c}_{s}=1$. Assuming a constant β₀, the magnetic field profile is given by (see also Johansen & Levin 2008)

$$\begin{eqnarray}&&{B}_{y}=\sqrt{{\rho }_{0}\exp [-{z}^{2}/2{H}_{\beta }^{2}]{{\rm{\Omega }}}_{0}^{2}/{\beta }_{0}},\end{eqnarray} \tag{ 9 }$$

where we have taken ${B}_{x}={B}_{z}=0$. It is important to emphasize that the criteria for triggering the PRTI will be slightly different compared with the well-known work of Parker (1966), since we are dealing with a differential rotation system under a linear gravity. The presence of differential rotation was studied analytically by Foglizzo & Tagger (1994, 1995), who obtained new instability conditions and studied the correlations between the PRTI and MRI. Kim et al. (1997), on the other hand, studied the instability criteria under a linear gravity, which is more appropriate for the case of Keplerian disks and the present work. However, in all works, an initial β₀ of the order of unity is still required to trigger the PRTI.

We have assumed the initial velocity field to be zero (${\boldsymbol{u}}=0$), since an orbital advection scheme has also been adopted (see Equation (7) and Stone & Gardiner 2010). However, a Gaussian perturbation with an amplitude of $\delta u\sim {10}^{-3}$ has been used in the azimuthal velocity component to trigger the linear phase of the PRTI and MRI (see Johansen & Levin 2008).

We have used shearing-periodic boundary conditions in the radial direction x that reproduce the differential rotation through the azimuthal displacement of the radial boundaries (see Hawley et al. 1995). We have used standard periodic boundary conditions in the azimuthal direction y and outflow boundaries in the vertical direction z, which are more appropriate to the study and evolution of the coronal region in the surroundings of accretion disks (see Bai & Stone 2013). In outflow boundary conditions, zero gradients are used for the velocity and magnetic fields, but the vertical velocity component is set to zero (${v}_{z}=0$) whenever inflows are detected through the boundaries. Besides, the density is extrapolated assuming a vertical hydrostatic equilibrium.

Finally, a density floor (${\rho }_{\min }$) has been applied to minimize the numerical limitations due to the high "physical" Alfvén speed (see, e.g., Bai & Stone 2013) since the rarefied corona is subject to the action of a strong magnetic field (that is transported by the magnetic buoyancy) with the ratio between the thermal and magnetic pressure of the order of unity. We notice that this floor could change the physics of the problem by breaking down the initial magnetostatic equilibrium. However, we verified along the simulation that the horizontally averaged densities are larger than the adopted density floor (as in Bai & Stone 2013), so that this procedure has a minimum impact on the evolution of the system. We have also used a mass diffusion term in the mass conservation equation (see, e.g., Gressel et al. 2011) to alleviate the density floor condition⁴ and allow us to use ${\rho }_{\min }={10}^{-6}$ for the low resolution; ${\rho }_{\min }={10}^{-6}$, 10⁻⁷, or 10⁻⁸ for the intermediate resolution; and ${\rho }_{\min }={10}^{-5}$ for the high resolution. For two of the models with intermediate resolution (see the next subsection for more details), due to the weak magnetic field intensity applied as an initial condition and the magnetostatic equilibrium adjustment, the density profile in the vertical direction has been found to decay faster than in the other models for floor values of 10⁻⁶ or 10⁻⁵, thus affecting strongly the initial evolution of the system, as mentioned above. In these cases, we adopted floor values of 10⁻⁷ and 10⁻⁸, which are lower than the horizontally averaged density at the high altitudes in these models.

In this work, we have used resolutions of ∼11 ($128\,\times 128\times 128$), ∼21 ($256\times 256\times 256$), and ∼43 ($512\times 512\,\times 512$) cells per thermal height scale of the system ($H={c}_{s}{{\rm{\Omega }}}^{-1}$), considering a computational domain of size $12H\times 12H\times 12H$. We also considered three different values of β₀ (1, 10, and 100) to evaluate the evolution of the MRI, PRTI, and turbulence in the system. The parameters of the simulations are shown in Table 1, and each model name is composed of the resolution (R11, R21, and R43) plus the initial value of β₀ (b1, b10, and b100). The diagnostics used to analyze these models are presented in the next section.

Table 1.  Simulation Parameters

Simulation	Computational Domain	Resolution	β0	Ω	q	${\rho }_{\min }$	Vertical Boundaries
R11b1	$12H\times 12H\times 12H$	$128\times 128\times 128$	1.0	1.0	1.5	10−6	Outflows
R21b1	$12H\times 12H\times 12H$	$256\times 256\times 256$	1.0	1.0	1.5	10−6	Outflows
R21b10	$12H\times 12H\times 12H$	$256\times 256\times 256$	10	1.0	1.5	10−7	Outflows
R21b100	$12H\times 12H\times 12H$	$256\times 256\times 256$	102	1.0	1.5	10−8	Outflows
R43b1	$12H\times 12H\times 12H$	$512\times 512\times 512$	1.0	1.0	1.5	10−5	Outflows

Download table as:  ASCIITypeset image

In order to follow the evolution of the system and the instabilities, we will present here the time- and space-average diagrams of the relevant quantities. The space average of a given variable f taken at the entire volume of the box is represented as $\langle f\rangle$, so that

$$\begin{eqnarray}&&\langle f\rangle =\displaystyle \frac{\int {fdxdydz}}{\int {dxdydz}}.\end{eqnarray} \tag{ 10 }$$

Also, in order to track the evolution of f at the vertical direction z, we performed space averages at the xy plane, at each height z (horizontal averages), represented by the symbol ${\langle f\rangle }_{{xy}}$, where

$$\begin{eqnarray}&&{ \langle f \rangle }_{{xy}}=\displaystyle \frac{\int {fdxdy}}{\int {dxdy}}.\end{eqnarray} \tag{ 11 }$$

Finally, the time averages are indicated by a t index, so that time volume and horizontal averages are represented by ${\langle f\rangle }_{t}$ and ${\langle f\rangle }_{{xyt}}$, respectively.

In this work, we have evaluated the angular momentum transport parameter α (Shakura & Sunyaev 1973) through the relation

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{{T}_{{xy}}}{P}=\displaystyle \frac{{T}_{{xy}}^{{\rm{Max}}}+{T}_{{xy}}^{{\rm{Rey}}}}{\rho {c}_{s}^{2}},\end{eqnarray} \tag{ 12 }$$

where ${T}_{{xy}}^{\mathrm{Max}}=-{B}_{x}{B}_{y}$ is the Maxwell (magnetic) stress tensor and ${T}_{{xy}}^{\mathrm{Rey}}=\rho {u}_{x}{u}_{y}$ is the Reynolds stress tensor. We have also evaluated the ratio between the thermal and magnetic pressures, given by

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\rho {c}_{s}^{2}}{{B}^{2}/2}.\end{eqnarray} \tag{ 13 }$$

The volume and horizontal averages of Equations (12) and (13) are represented by an upper bar (i.e., $\langle \overline{\alpha }\rangle$, ${\langle \overline{\alpha }\rangle }_{{xy}}$, $\langle \overline{\beta }\rangle$, and ${\langle \overline{\beta }\rangle }_{{xy}}$) to indicate the ratio of the averages of two different quantities, for example:

$$\begin{eqnarray}&&\langle \overline{\alpha }\rangle =\displaystyle \frac{\langle {T}_{{xy}}\rangle }{\langle \rho {c}_{s}^{2}\rangle }\,\,\mathrm{and}\,\,{\langle \overline{\alpha }\rangle }_{{xy}}=\displaystyle \frac{{\langle {T}_{{xy}}\rangle }_{{xy}}}{{\langle \rho {c}_{s}^{2}\rangle }_{{xy}}}.\end{eqnarray} \tag{ 14 }$$

To evaluate the power spectrum of the velocity field ($| \widetilde{{\boldsymbol{u}}}({k}_{x},{k}_{y},z){| }^{2}$), we performed a two-dimensional Fourier transformation to obtain the k_x and k_y wavenumbers in each height z of the domain:

$$\begin{eqnarray}\widetilde{{\boldsymbol{u}}}({k}_{x},{k}_{y},z) & = & \displaystyle \frac{1}{\sqrt{{N}_{x}{N}_{y}}}\\ & & \times \ \displaystyle \sum _{n,m=0}^{{N}_{x}-1,{N}_{y}-1}{\boldsymbol{u}}({x}_{n},{y}_{m},z){e}^{-2\pi i\left(\displaystyle \frac{{k}_{x}{x}_{n}}{{N}_{x}}+\displaystyle \frac{{k}_{y}{y}_{m}}{{N}_{y}}\right)},\end{eqnarray} \tag{ 15 }$$

where N_x and N_y are the total number of cells in the radial and azimuthal directions, respectively. The spectrum $P(k,z)$ has been obtained by integrating $| \widetilde{{\boldsymbol{u}}}({k}_{x},{k}_{y},z){| }^{2}$ over annular areas between k − dk and k, where $k=\sqrt{(}{k}_{x}^{2}+{k}_{y}^{2})$.

We adopt the R21b1 model (see Table 1) as our reference model. We describe the main results of this simulation in the next section.

Figure 1 shows the time evolution of the magnetic field (streamlines) and the density distribution (background colors in a logarithmic scale) at t = 0, $1.5P$, $4P$, and $8P$ (where $P=2\pi {{\rm{\Omega }}}_{0}^{-1}$ corresponds to one orbital period of the accretion disk), for the R21b1 model (see Table 1). The configuration of the magnetic field lines at t = 1.5P shows that the azimuthal component B_y is stretched in the vertical direction z due to the effects of magnetic buoyancy, producing loops during the exponential growth of the PRTI (which are in agreement with the results obtained by Johansen & Levin 2008). The development of the turbulence, and the decrease of the intensity of the magnetic field and its transport from the midplane of the disk to the coronal region due to PRTI (and to outside the computational domain) are observed at t = 4P and 8P.

Zoom In Zoom Out Reset image size

Figure 2 shows the evolution of the volume-averaged magnetic energy density $\langle {B}^{2}\rangle /2$ (top diagram) and the α-parameter (bottom diagram) over 100 orbital periods. The top diagram indicates an amplification of the B_x and B_z components during the first five orbital periods. Estimating the exponential growth time from the linear perturbation solutions for the PRTI obtained by Kim et al. (1997), for the higher altitudes of the system (3H–6H), we find values ($\sim 1.5P-4.5P$) that are in agreement with our simulations, although Kim et al. (1997) ignored differential rotation in their evaluation. After five orbital periods, the total magnetic field decreases due to expansion and upward transport through the outflow boundaries in the vertical direction z. A small dissipation due to turbulent magnetic reconnection probably also contributes to this decrease (see more details in Section 4). At $t\sim 30P$, the decrease of the magnetic field intensity and consequent increase of β induces the growth of MRI (from the azimuthal and vertical magnetic field components), which will dominate the system evolution for the rest of the simulation.

Zoom In Zoom Out Reset image size

The bottom diagram of Figure 2 shows that α is dominated by the Maxwell stress tensor and follows the same trend as the magnetic energy density, as we expected. The diagram also indicates that there is a high accretion regime in the first 10 orbital periods when it achieves a maximum value around $\langle \overline{\alpha }\rangle \sim 0.3$ (consistent with those expected from the observations; see King et al. 2007).

3.1.1. Accretion Disk and Corona Evolution

In order to evaluate the evolution of the accretion disk and the corona separately, we have obtained the horizontal averages of the B_x and B_y magnetic field components,⁵ and the β parameter as a function of the height of the system.The first and second diagrams in Figure 3 show a change of regimes between t = 5P and t = 10P of the radial and azimuthal components of the magnetic field. Besides, during the regime where PRTI dominates, there is an increase of the radial component with a peak of the intensity in the middle of the disk, followed by an inversion of polarity in the coronal regions between $| z| =2H$ and $4H$. In the higher altitudes of the corona, between $| z| =4H$ and $6H$, there are smaller polarity inversions.

Zoom In Zoom Out Reset image size

After t = 10P, the first and second diagrams of Figure 3 show that the radial and azimuthal components of the magnetic field undergo polarity inversions on a timescale of approximately 10 orbital periods. This is also a well-known result consistent with previous studies of the MRI in weak-field regimes (β ≫ 1) and zero net flux (e.g., Brandenburg et al. 1995; Davis et al. 2010; Simon et al. 2011, 2012; Shi et al. 2016), which have produced similar butterfly diagrams to those we see in Figure 3. This pattern suggests the action of an alpha–Omega dynamo process triggered both by differential rotation (Omega effect) and the turbulent cyclonic motions driven by the PRTI and MRI (alpha effect). The first one stretches the radial and vertical magnetic field lines⁶ in the azimuthal direction. The alpha effect allows for the amplification of the radial component B_x from the azimuthal component B_y, through the electromotive forces in the azimuthal direction ${({\boldsymbol{u}}\times {\boldsymbol{b}})}_{y}$, where ${\boldsymbol{u}}$ and ${\boldsymbol{b}}$ correspond to the fluctuations of the velocity and magnetic field, respectively. We also note that in Figure 3, the polarity inversions of the B_x and B_y components at the end of each half cycle (where the old field is replaced by a new one with the opposite sign) are in phase opposition, as expected in a shear dynamo process (for details, see, e.g., Brandenburg et al. 1995; Guerrero & de Gouveia Dal Pino 2007a, 2007b, 2008; Guerrero et al. 2009, 2016a, 2016b, and references therein).

Finally, the third diagram of Figure 3 shows the time evolution of $\langle \overline{\beta }{\rangle }_{{xy}}$. In the coronal region, the intensity of $\langle \overline{\beta }{\rangle }_{{xy}}$ decreases due to the growth of the magnetic field strength in this region induced by the PRTI. Thereafter, it saturates around an average value of $\langle \overline{\beta }{\rangle }_{{xy}}\sim 0.8$ with a strong variability (between ∼0.5 and 1.0). Besides, even with the initial highly magnetized system (${\beta }_{0}=1$), the disk evolves to a gas-pressure-dominated regime with a horizontal-average value around 10. Once the thermal pressure becomes larger than the magnetic pressure ($\beta \gt 1$), the MRI is expected to settle in and soon dominate the dynamics of the disk. The MRI evolves out of the azimuthal and vertical components of the magnetic field produced by the PRTI (see also, e.g., Balbus & Hawley 1992; Foglizzo & Tagger 1995; Hawley et al. 1995; Johansen & Levin 2008).

These results, although out of the main scope of this work, are compatible with the notion that the MRI starts to dominate the dynamics of the system when β grows to values greater than unity inside the disk. This inhibits the PRTI and at the same time gives rise to dynamo amplification of the azimuthal field and some amplification of the radial field which started to grow by the dynamo process early in the PRTI regime.

In this section, we show the evolution of the PRTI and MRI for different initial values of β₀ (1, 10 and 100) in a computational domain of size $12H\times 12H\times 12H$ and a resolution of 21H⁻¹ (models R21b1, R21b10, and R21b100, respectively). The presence of a weak vertical magnetic field is essential for the development of the MRI (see, e.g., Chandrasekhar 1960; Balbus & Hawley 1991; Hawley et al. 1995; Miller & Stone 2000; Salvesen et al. 2016a, 2016b). In contrast, the PRTI evolves in regimes where the magnetic pressure is of the order of the thermal pressure (${\beta }_{0}\sim 1;$ see, e.g., Parker 1966). It is expected, therefore, that for systems with initially weak magnetic fields (or large β), the PRTI should not dominate the initial evolution of the system.

The diagrams of Figure 4 show the time evolution of $\langle {B}^{2}\rangle /2$ (top diagram) and $\langle \overline{\alpha }{\rangle }_{\mathrm{mag}}=\langle -{B}_{x}{B}_{y}\rangle /\langle \rho {c}_{s}^{2}\rangle$ (bottom diagram). The reference model R21b1 with ${\beta }_{0}=1$ discussed in the previous sections is also shown for comparison (black line). The top diagram shows clearly that the initial evolution of the magnetic field is different for models with distinct β₀. Nevertheless, after 20 orbital periods, all systems evolve in a similar way, achieving a nearly steady-state situation. For the R21b10 model (with ${\beta }_{0}=10$), the initial amplification of the magnetic field and its subsequent decrease at the steady-state phase indicates that the magnetic buoyancy still plays an important role in the transport of the magnetic field from the disk to outside the system through the vertical boundaries, while the R21b100 model (with ${\beta }_{0}=100$) shows the standard evolution of the magnetic energy density driven by the MRI only, with an initial amplification until saturation is achieved after 20 orbital periods (e.g., Hawley et al. 1995; Miller & Stone 2000).

Zoom In Zoom Out Reset image size

The α-parameter also shows convergence for all models after 20 orbital periods (bottom diagram of Figure 4). The initial peak value achieved during the growth phase of the PRTI seen in the R21b1 model gets smaller for increasing β₀, as expected, and absent for ${\beta }_{0}=100$.

The diagrams of Figure 5 show the time evolution of $\langle \overline{\beta }{\rangle }_{{xy}}$ for the models discussed in this section. Consistent with the other diagrams, the initial growth of the magnetic field due to the PRTI in the models with initial ${\beta }_{0}\lt 100$ leads to a decrease in β and then, after 20 orbital periods, all of the models converge to the same steady-state values, both in the disk and corona, when the MRI sets in. At this stage, regardless of the initial β₀, the disk becomes a gas-pressure-dominated system surrounded by a much more magnetized corona (with $\beta \simeq 1$), which is in agreement with Miller & Stone (2000) and Salvesen et al. (2016a). Finally, the differences found in all diagrams show that the PRTI has an important role just in an early "transient" phase.

Zoom In Zoom Out Reset image size

As is well known, the evolution of the MRI triggered by an initial azimuthal magnetic field is highly dependent on the resolution of the computational domain, since it is not possible to resolve all growing wavenumbers, especially the k_z modes (see Hawley et al. 1995). Besides, considering a zero net vertical flux in stratified shearing-box simulations, Ryan et al. (2017) have found that the α-parameter is resolution dependent, which is in contrast to the results of, e.g., Davis et al. (2010). In this section, we will discuss the role of the resolution in the evolution of the PRTI and MRI and check whether the reference model (R21b1) converges numerically to a stationary state.

Figure 6 shows the time evolution of $\langle {B}^{2}\rangle /2$ (top diagram) and $\langle \overline{\alpha }{\rangle }_{\mathrm{mag}}$ (bottom diagram) for the resolutions of 11H⁻¹ (R11b1 model, black line), 21H⁻¹ (R21b1 model, red line), and 43H⁻¹ (R43b1 model, blue line). The high-resolution simulation (R43b1) is numerically costly, and we have evolved over a short time equivalent to ∼58 orbital periods. During the first five orbital periods (where the PRTI is dominant) until t = 20P, all of the resolutions show the same behavior for the magnetic energy density, indicating that the PRTI operates appropriately even in our lower-resolution simulation (R11b1). It is known that the growth rate of the PRTI is favored by large wavenumbers (see, e.g., Parker 1966, 1967; Kim et al. 1997), so we might expect that the resolution of 11H⁻¹ would not affect significantly the initial evolution of such an instability. Nevertheless, the maximum value of $\langle \overline{\alpha }{\rangle }_{\mathrm{mag}}$ for this model is slightly lower than for the intermediate and high-resolution simulations (R21b1 and R43b1 models, respectively), indicating that resolutions lower than 11H⁻¹ could affect the initial evolution of the system.

Zoom In Zoom Out Reset image size

After 20 orbital periods, when the MRI dominates, the differences in the evolution between the low- and intermediate-resolution models are significant. The magnetic field in the low-resolution model continues to decrease without achieving saturation. Such a decrease indicates that the MRI in this case is unable to regenerate the magnetic field at the same rate that magnetic flux is transported throughout the vertical boundaries. As mentioned, such a behavior is due to the sensitivity of the MRI to the resolution when triggered by an initial azimuthal magnetic field (see Hawley et al. 1995). The comparison between the intermediate- and high-resolution models indicate a good agreement between them until the evolved time of the high-resolution model R43b1 (around 58 orbits), though it is not possible to determine yet whether it will achieve the same saturation as in the R21b1 model.

In order to identify magnetic reconnection sites in our disk/corona system, we have searched for current density peaks (${\boldsymbol{J}}={\rm{\nabla }}\times {\boldsymbol{B}}$) and then evaluated the local magnetic reconnection rate (V_rec) in the surroundings of each peak. To this aim, we have adapted the algorithm of Zhdankin et al. (2013) and extended it to a 3D analysis. The algorithm is well described in Zhdankin et al. (2013), and we will summarize the most important steps in this section. First, we have selected a sample of cells with a current density value greater than ${j}_{5}=5\langle | {\boldsymbol{J}}| \rangle$, where $\langle | {\boldsymbol{J}}| \rangle$ is the volume average of the total current density taken over the disk and the corona separately. From this first sample, we have selected those cells that have local maxima ($| {\boldsymbol{J}}{| }_{{ijk}}^{\max }$, where ijk corresponds to the index position of each cell) within a surrounding cubic subarray of the data (with a size of $7\times 7\times 7$ cells), and then we checked whether they are located between magnetic field lines of opposite polarity (in each component separately). In the present work, we are not interested in identifying null points since the magnetic reconnection topology is complex in a 3D regime (see, e.g., Yamada et al. 2010). Besides, this step is quite different from the original algorithm because Zhdankin et al. (2013) used an X-point model (see Petschek 1964) to represent the magnetic reconnection sites, and then they identified such events as saddle points in the magnetic flux function.

We have assumed that the cells of the last subsample are possible sites of reconnection. However, the topology of these sites is complex (as mentioned above), and they are not necessarily aligned with one of the axes of the Cartesian coordinate system. So, we adopted a new coordinate system centered in the local reconnection site obtained from the eigenvalues and eigenvectors of the current density 3D Hessian matrix (H_ijk) for each cell of the last subsample (see Zhdankin et al. 2013):

$$\begin{eqnarray}{H}_{{ijk}}=\left[\begin{array}{ccc}{\partial }_{{xx}}| {\boldsymbol{J}}{| }_{{ijk}} & {\partial }_{{xy}}| {\boldsymbol{J}}{| }_{{ijk}} & {\partial }_{{xz}}| {\boldsymbol{J}}{| }_{{ijk}}\\ {\partial }_{{yx}}| {\boldsymbol{J}}{| }_{{ijk}} & {\partial }_{{yy}}| {\boldsymbol{J}}{| }_{{ijk}} & {\partial }_{{yz}}| {\boldsymbol{J}}{| }_{{ijk}}\\ {\partial }_{{zx}}| {\boldsymbol{J}}{| }_{{ijk}} & {\partial }_{{zy}}| {\boldsymbol{J}}{| }_{{ijk}} & {\partial }_{{zz}}| {\boldsymbol{J}}{| }_{{ijk}}\end{array}\right],\end{eqnarray} \tag{ 16 }$$

where H_ijk corresponds to the second-order partial derivatives of the current density magnitude. The highest eigenvalue indicates that the associated eigenvector corresponds to the direction of the fastest decrease (or the highest variance) of $| {\boldsymbol{J}}|$. We have assumed this direction as the thickness of the magnetic reconnection site. The eigenvectors of H_ijk provide the three orthonormal vectors of the new coordinate system centered on the local reconnection site (see the upper sketch of Figure 7).

Zoom In Zoom Out Reset image size

The local magnetic reconnection rate has been evaluated in a similar way as in Kowal et al. (2009), where we averaged the inflow velocity (${V}_{\mathrm{in}}={V}_{\hat{{{\boldsymbol{e}}}_{1}}}$, the projection onto the $\hat{{{\boldsymbol{e}}}_{1}}$ direction) divided by the Alfvén speed at the edges of the reconnection site (see the bottom sketch of Figure 7):

$$\begin{eqnarray}&&\left.\left\langle \displaystyle \frac{{V}_{\mathrm{in}}}{{V}_{{\rm{A}}}}\right\rangle =\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{V}_{\hat{{{\boldsymbol{e}}}_{1}}}}{{V}_{{\rm{A}}}}\right|}_{\mathrm{lower}}-{\left.\displaystyle \frac{{V}_{\hat{{{\boldsymbol{e}}}_{1}}}}{{V}_{{\rm{A}}}}\right|}_{\mathrm{upper}}\right),\end{eqnarray} \tag{ 17 }$$

where the Alfvén speed is given by

$$\begin{eqnarray}&&{V}_{{\rm{A}}}=\sqrt{\displaystyle \frac{{B}_{\hat{{{\boldsymbol{e}}}_{1}}}^{2}+{B}_{\hat{{{\boldsymbol{e}}}_{2}}}^{2}+{B}_{\hat{{{\boldsymbol{e}}}_{3}}}^{2}}{\rho }},\end{eqnarray} \tag{ 18 }$$

and ${B}_{\hat{{{\boldsymbol{e}}}_{1}}}$, ${B}_{\hat{{{\boldsymbol{e}}}_{2}}}$, and ${B}_{\hat{{{\boldsymbol{e}}}_{3}}}$ correspond to the projection of the magnetic field onto the three eigenvectors of the Hessian matrix. As in Zhdankin et al. (2013), we have identified the edges and the cells belonging to a given local magnetic reconnection site considering only those with current density $| J|$ greater than half of the maximum local value given by $| {\boldsymbol{J}}{| }_{{ijk}}^{\max }$ ($| J| \gt 1/2| {\boldsymbol{J}}{| }_{{ijk}}^{\max }$, see the bottom sketch of Figure 7). We have constrained the subsample to obtain V_rec for the most symmetric profiles, since a symmetry of both the magnetic and velocity field profiles is expected around the magnetic reconnection site (as seen in Figure 7). As an example, Figure 8 shows the spatial distribution of the local maxima identified by the algorithm, in the coronal region (upper corona, $3H\lt z\lt 6H$), at t = 5P for the R21b1 model. The white points correspond to the total subsample and the green points to the constrained subsample, as described above. The diagrams of Figure 8 show examples of rejected and accepted profiles of ${V}_{\hat{{{\boldsymbol{e}}}_{1}}}$, ${B}_{\hat{{{\boldsymbol{e}}}_{2}}}$, and ${B}_{\hat{{{\boldsymbol{e}}}_{3}}}$ interpolated along the ${\hat{e}}_{1}$ axis. This constraint reduces considerably the subsample size, but also reduces the standard deviation of V_rec (see a detailed discussion in Sections 4.3 and 4.5).

Zoom In Zoom Out Reset image size

Finally, despite the challenge to represent the real topology of the reconnection site in our 3D-MHD simulations, as mentioned above, Figure 9 shows an example of one of the accepted regions using the line integral convolution (LIC) method combined with 2D maps of the current density (top diagram) and the magnetic field intensity (bottom diagram) at t = 5P. The diagrams have been obtained by interpolating the data in the surroundings of the reconnection site and corresponding to areas of 21 cells². In this example, the ${\hat{e}}_{3}$ axis of the reconnection region (see Figure 7) is approximately aligned to the y-axis of the Cartesian coordinate system, allowing us to visualize the topology in the xz plane. The magnetic field magnitude was evaluated using the B_x and B_z components. The bottom diagram of Figure 9 shows at least three possible reconnection events (white square and circles) at regions where the magnetic lines converge where the magnetic intensity decreases. However, the algorithm identified only one of these regions (white square) that corresponds to the local maxima of $| {\boldsymbol{J}}|$ (see the top diagram of Figure 9). It is beyond the scope of this work to provide a detailed analysis of all possible reconnection sites, but we can speculate that the two other regions (circles) correspond to already reconnected sites where substantial magnetic energy has already dissipated.

Zoom In Zoom Out Reset image size

As mentioned above, the algorithm checks whether each sample of cells is located between magnetic field lines of opposite polarity. We have also analyzed the distribution in time of the sign flips of the magnetic field in each direction separately (x, y, and z) at spatial scales corresponding to the size of a cell. These distributions can help identify what magnetic field components are reconnecting inside the system as a function of time, e.g., during the PRTI and MRI phases. As an example, Figure 10 shows such distributions (with bins of 10 orbital periods) for the model R21b1, where ${{\rm{\Delta }}}_{j}({B}_{i}{)}_{-}^{+}$ corresponds to the number of sign flips of the magnetic field component i in the direction j. The upper and lower diagrams correspond to the distribution taken in the coronal region and disk, respectively. Initially, between 0 and 10 orbital periods at the coronal region, the number of sign flips is dominated mainly by the B_x component along the vertical direction z (${{\rm{\Delta }}}_{z}({B}_{x}{)}_{-}^{+}$) and the B_z component along the radial direction x (${{\rm{\Delta }}}_{x}({B}_{z}{)}_{-}^{+}$). This behavior indicates that reconnections of the azimuthal field are not relevant during this initial phase of the PRTI at the corona, showing that the field lines initially twist, producing flux tubes in the azimuthal direction (see also Hirose et al. 2006; Simon et al. 2012) and allowing for reconnection to occur in the x and z directions mainly. This behavior can be clearly seen in the diagrams of Figure 1. On the other hand, after t = 10P, the number of sign flips of the B_y component in the x and z directions (${{\rm{\Delta }}}_{x}({B}_{y}{)}_{-}^{+}$ and ${{\rm{\Delta }}}_{z}({B}_{y}{)}_{-}^{+}$, respectively) increases, playing an important role in the reconnection, too, which is compatible with the inversion of the polarity of the azimuthal magnetic field during the MRI phase (as seen in the middle diagram of Figure 3). The sign flips in the azimuthal direction (${{\rm{\Delta }}}_{y}({B}_{x}{)}_{-}^{+}$ and ${{\rm{\Delta }}}_{y}({B}_{z}{)}_{-}^{+}$) are less relevant over the entirety of the evolution of the simulation, as expected, since the shear and stretching prevent local encounters of magnetic field lines of opposite polarity.

Zoom In Zoom Out Reset image size

A similar behavior can be found in the disk region, except for the B_z component in the x direction (${{\rm{\Delta }}}_{x}({B}_{z}{)}_{-}^{+}$), for which the flips are less relevant than in the coronal region. During the entire simulation, the sign flips are dominated by ${{\rm{\Delta }}}_{z}({B}_{y}{)}_{-}^{+}$, ${{\rm{\Delta }}}_{x}({B}_{y}{)}_{-}^{+}$, and ${{\rm{\Delta }}}_{z}({B}_{x}{)}_{-}^{+}$. This is not a surprise since the B_z component is produced mainly by loops at higher altitudes above and below the disk during the exponential growth of the PRTI (as described in Section 3).

The distribution of $\langle {V}_{\mathrm{rec}}\rangle$ was obtained for each time step of the simulations, and Figure 11 depicts the time evolution of the average value of V_rec for the upper and lower coronal regions together and the disk (continuous black line). These averages were calculated out of 800 histograms computed over 100 orbital periods. The diagrams also depict the standard deviation (blue shade) of each distribution.

Zoom In Zoom Out Reset image size

In the corona, $\langle {V}_{\mathrm{rec}}\rangle$ increases very fast in the early PRTI regime (during the first three orbital periods) and then saturates (after 10 orbital periods) to a mean value varying between ∼0.11 and 0.23. In the disk, $\langle {V}_{\mathrm{rec}}\rangle$ also increases fast in the PRTI regime, achieving a peak value (∼0.2) around ∼6 orbital periods, then decreases to nearly constant value (0.1–0.18) after t = 20P when the MRI sets in, following the same trend of the time evolution of the α-parameter and magnetic energy (Figure 4).

These values of $\langle {V}_{\mathrm{rec}}\rangle$ are in agreement with the predictions of the theory of fast magnetic reconnection driven by turbulence (Lazarian & Vishniac 1999) and numerical studies that have tested it (see, e.g., Kowal et al. 2009; Singh et al. 2015; Takamoto et al. 2015). They are also compatible with observations of fast magnetic reconnection in the solar corona associated with flares (∼0.001–0.1 and up to 0.5; see, e.g., Dere 1996; Aschwanden et al. 2001; Su et al. 2013). This indicates that the algorithm used here can be an efficient tool for the identification of magnetic reconnection sites in numerical MHD simulations of systems involving turbulence. Furthermore, it has demonstrated the ability of PRTI- and MRI-induced turbulence to drive fast reconnection in accretion disks and their coronae.

Figure 12 shows the histograms of V_rec (diagrams on the left side) and the measure of the thickness of the reconnection regions (diagrams on the right side) for the R21b1 model obtained between 20 and 100 orbital periods in the coronal region (upper diagrams) and disk (lower diagrams), where the reconnection rate, as the other quantities, achieves a nearly steady state. In both cases, the distributions of V_rec and the thickness do not resemble a normal distribution, showing a long tail on the right (a skewed distribution). We should note that the algorithm has identified in the tail very few events with magnetic reconnection rates larger than 1.0, but we have constrained the histogram to the range of V_rec between 0.0 to 0.5, which corresponds to 99.9% of the total data (a similar procedure was applied for the histograms of the thickness). This may reflect the limitations of the method to calculate the magnetic reconnection rate since the inflow velocities at the upper and lower edges in the reconnection site are not perfectly symmetric. Figure 12 also depicts for the skewed distributions both the median (evaluated from the total data) and the mean values with the standard deviations. For the coronal region, we have obtained an average reconnection rate of the order of 0.17 ± 0.10 and a median of ∼0.16, whereas in the disk we have obtained an average of 0.13 ± 0.09 and a median of ∼0.11. The average values of the thickness show that the reconnection sites occupy two or three cells (for the resolution of 21H⁻¹); therefore, numerical effects could affect the evaluation of V_rec. We have also performed tests with larger reconnection sites, but the averages of V_rec did not change significantly.

Zoom In Zoom Out Reset image size

As mentioned in Section 4.1, we have constrained our sample in the analysis above considering only symmetric profiles of the velocity and magnetic fields around the reconnection sites. Here, we repeated the same analysis using the whole sample (imposing no restriction) and another one including both nonsymmetric and symmetric profiles. In both cases, the mean values and standard deviations are higher than the ones shown previously (in Figure 11), as we expected since V_rec has been evaluated through Equation (17), where we averaged the ratio between the inflow velocity and Alfvén speed at the top and bottom of the reconnection site. When considering the whole sample, we have obtained ${V}_{\mathrm{rec}}^{\mathrm{corona}}\,=0.24\pm 0.24$ and ${V}_{\mathrm{rec}}^{\mathrm{disk}}=0.20\pm 0.25$, while for the sample including both symmetric and nonsymmetric profiles, we obtained ${V}_{\mathrm{rec}}^{\mathrm{corona}}=0.17\pm 0.13$ and ${V}_{\mathrm{rec}}^{\mathrm{disk}}=0.14\pm 0.15$. With regard to the thickness of the reconnection regions in these cases, the average values have not changed with respect to the previous results (Figure 12). This is not a surprise since the criteria to evaluate the thickness are independent of the symmetry of the velocity and magnetic profiles.

In order to verify more quantitatively the correlation between the magnetic reconnection rate and the turbulence, we have performed a Fourier analysis and evaluated the two-dimensional power spectrum of the velocity field⁷ ($P(k,z)=| \widetilde{{\boldsymbol{u}}}({k}_{x},{k}_{y},z){| }^{2}$) for the k_x and k_y wavenumbers in different heights (z direction) inside the computational domain. The power spectrum has been taken from the average in annular areas between k − dk and k (where $k=\sqrt{(}{k}_{x}^{2}+{k}_{y}^{2})$) in the k_x−k_y Fourier space, as described in Section 2.4 (see Equation (15)). We have also averaged in time with bins of one orbital period to reduce the fluctuations of the spectra. Such decomposition has been chosen in order to verify the turbulence level in the disk and coronal regions separately, since the diagrams of V_rec (see Figure 11) have shown significant differences between these two regions after t = 20P.

Figure 13 shows the power spectrum for the low-, intermediate-, and high-resolution simulations (R11b1, R21b1, and R43b1 models, respectively) compensated by a factor of ${k}^{5/3}$ in three different times, at t = 1P, 5P, and 50P, where the colors of each line correspond to the height z. As we expected, the turbulence is not fully developed in the early stages of the simulations. On the other hand, after five orbital periods, the turbulent power spectrum is much larger and shows the typical cascading to small scales with a slope (${k}^{\nu }$) between $-1.9\lt \nu \lt -2.4$ in the coronal region and disk (for all the resolutions). At t = 50P, the slope in the middle plane of the disk (z = 0) is $\nu \sim -1.7$ (i.e., Kolmogorov-like) and increases as a function of the altitude, reaching a value of $\nu \sim -2.5$ at the coronal region, which is closer to a $-8/3$ power law, typical of a 2D turbulent spectrum distribution (Equation (15)), suggesting a change of regimes as we go from the disk to the corona, i.e., from large to small values of β. We should remember that even at these evolved stages, the corona keeps traces of large-scale coherent magnetic loops embedded in the turbulent flow, developed during the PRTI regime. No significant differences are seen when considering different resolutions.⁸

Zoom In Zoom Out Reset image size

At t = 1P, for the R21b1 model (second diagram of Figure 13), the power is still very small, and the magnetic reconnection rate has values below 0.1 (see Figure 11). Around t = 5P, during the PRTI, the magnetic reconnection rate achieves its largest values between 0.1 in the disk and 0.2 in the coronal region. The power in the middle of the disk is smaller than in the coronal region and is consistent with the behavior of $\langle {V}_{\mathrm{rec}}\rangle$ in such regions (where the average value in the coronal region is higher than in the disk after 20 orbital periods). Similar behavior is found at t = 50P, indicating a correlation between the turbulence and V_rec.

Figure 13 also shows that the low-resolution model R11b1, despite the different behavior in Figure 6 that indicates a decrease of the volume-averaged magnetic energy density with time, has a velocity field power spectrum similar to the higher resolution models R21b1 and R43b1. This similarity between models of different resolutions will be further stressed in the next subsection.

Figure 14 compares the time evolution of $\langle {V}_{\mathrm{rec}}\rangle$ for simulations with different values of β₀ (models R21b1, R21b10, and R21b100, left diagrams) and resolutions (models R11b1, R21b1 and R43b1, right diagrams). Despite the high standard deviations (see, e.g., Figure 11), the left diagrams show that the behavior of $\langle {V}_{\mathrm{rec}}\rangle$ is similar to the one of the reference model R21b1 and consistent with the time evolution of these systems as discussed in Section 3.2. In the first orbital periods (before 20P), in the PRTI regime, the models with different β₀ show an increase of the reconnection rate, but with a significant delay for those with higher β₀. Above 20 orbital periods, regardless of the initial β₀, $\langle {V}_{\mathrm{rec}}\rangle$ converges to the average values discussed previously (between 0.11 and 0.23). At this stage, the volume averages of the magnetic energy density and α also converge to a steady state (as seen in Figure 4). The comparison of $\langle {V}_{\mathrm{rec}}\rangle$ with different resolutions reveals a good convergence mainly in the disk region (bottom diagram in the right side of Figure 14), whereas in the coronal region (top diagram in the right side), the mean values of V_rec for the high-resolution models are slightly smaller than for the low-resolution ones.

Zoom In Zoom Out Reset image size

Table 2 summarizes the time average values of the magnetic reconnection rate $\langle {V}_{\mathrm{rec}}{\rangle }_{t}$ , the thickness $\langle \delta {\rangle }_{t}$, and the number $\langle N{\rangle }_{t}$ of reconnection sites identified by the algorithm obtained between 20 and 100 orbital periods, except for the high-resolution simulation (model R43b1), whose time average was obtained between 20 and 58 orbital periods. Considering the standard deviations, all of the models show $\langle {V}_{\mathrm{rec}}{\rangle }_{t}$ values that are compatible both for different values of β and resolutions.

Table 2.  Time Average of V_rec, the Thickness (δ), and the Number of Identified Reconnection Sites (N)

Simulation	$\langle {V}_{\mathrm{rec}}{\rangle }_{t}\pm \sigma$		$\langle \delta {\rangle }_{t}\pm \sigma$		$\langle N{\rangle }_{t}\pm \sigma$		${\rm{\Delta }}T$
Name	Corona	Disk	Corona	Disk	Corona	Disk	Range
R11b1	0.20 ± 0.10	0.12 ± 0.07	0.17 ± 0.04	0.15 ± 0.02	25 ± 9	10 ± 5	20P−100P
R21b1	0.17 ± 0.10	0.13 ± 0.09	0.10 ± 0.03	0.08 ± 0.01	134 ± 34	107 ± 16	20P−100P
R21b1a	0.24 ± 0.24	0.20 ± 0.25	0.10 ± 0.03	0.08 ± 0.01	559 ± 74	402 ± 42	20P−100P
R21b1b	0.17 ± 0.13	0.14 ± 0.15	0.10 ± 0.03	0.08 ± 0.01	216 ± 44	163 ± 21	20P−100P
R21b10	0.17 ± 0.10	0.12 ± 0.09	0.09 ± 0.03	0.08 ± 0.01	137 ± 34	106 ± 15	20P−100P
R21b100	0.17 ± 0.10	0.12 ± 0.08	0.10 ± 0.03	0.08 ± 0.01	135 ± 33	105 ± 15	20P−100P
R43b1	0.15 ± 0.09	0.12 ± 0.08	0.05 ± 0.01	0.04 ± 0.01	726 ± 103	675 ± 50	20P−58P

Notes.

^aSample with the whole data (without restrictions; see Section 4.3) . ^bSample with symmetric and nonsymmetric profiles (see Section 4.3).

Download table as:  ASCIITypeset image

The averaged thickness $\langle \delta {\rangle }_{t}$ for the high-resolution model is significantly smaller than that of the intermediate- and low-resolution models, as we expect, since the structures of the magnetic reconnection sites are better resolved for the R43b1 model. Besides, $\langle {V}_{\mathrm{rec}}{\rangle }_{t}$ is not strongly affected, neither by the thickness nor by the numerical resistivity⁹ for different resolutions. This is not a surprise since, according to the turbulence-induced fast reconnection theory (see Lazarian & Vishniac 1999; Eyink et al. 2011, 2013), the presence of turbulence speeds up the reconnection independently of the Ohmic resistivity (here mimicked by the numerical resistivity, as stressed before; see also Equation (2)).

Finally, as the resolution increases, the number of identified reconnection sites increases as a consequence of the better-resolved magnetic structures, which improve the statistical analysis of V_rec. For this reason, the right diagrams in Figure 14 show that the amplitude variability of $\langle {V}_{\mathrm{rec}}\rangle$ decreases for the models with higher resolution. Above 20 orbital periods, the model R43b1 shows $\langle {V}_{\mathrm{rec}}\rangle$ values varying between ∼0.13 and 0.17 in the coronal region, and between ∼0.11 and 0.13 in the disk. On the other hand, the model R11b1 shows $\langle {V}_{\mathrm{rec}}\rangle$ values between ∼0.10 and 0.30 in the coronal region, and between ∼0.0 and 0.26 in the disk. Despite the increase in the number of identified reconnection sites, the distributions in time seen in Figure 10 do not change significantly for different resolutions.