MHD instabilities in accretion disks and their implications in driving fast magnetic reconnection
Abstract
Magnetohydrodynamic instabilities play an important role in accretion disks systems. Besides the wellknown effects of the magnetorotational instability (MRI), the ParkerRayleighTaylor instability (PRTI) also arises as an important mechanism to help in the formation of the coronal region around an accretion disk and in the production of magnetic reconnection events similar to those occurring in the solar corona. In this work, we have performed threedimensional magnetohydrodynamical (3DMHD) shearingbox numerical simulations of accretion disks with an initial stratified density distribution and a strong azimuthal magnetic field with a ratio between the thermal and magnetic pressures of the order of unity. This study aimed at verifying the role of these instabilities on the formation of the coronal region and especially on the driving of turbulence and magnetic reconnection. Our simulations revealed an initial formation of largescale magnetic loops due to the PRTI followed by the development of a nearly steadystate turbulence driven by both instabilities. We have employed an algorithm to identify the presence of current sheets produced by the encounter of magnetic flux ropes of opposite polarity in the turbulent regions of both the corona and the disk. We computed the magnetic reconnection rates in these locations obtaining average reconnection velocities in Alfvén speed units of the order of in the accretion disk and in the coronal region (with mean peak values of order of ), which are consistent with the predictions of the theory of turbulenceinduced fast reconnection.
DRAFT
Luis H.S. Kadowaki
0000000269085634]Luis H.S. Kadowaki \altaffiliationFAPESP Fellowship
1 Introduction
Accretion phenomena are ubiquitous in Astrophysical systems at all scales (for reviews see, e.g., Pringle, 1981; Balbus & Hawley, 1998; Abramowicz & Fragile, 2013). In the absence of rotation, the accretion processes should be spherical (Bondi, 1952). However, in systems with significant rotation, accretion is through a disk and this process is strongly limited by the action of a centrifugal barrier (in a near Keplerian regime). Therefore, some process must act to transport angular momentum to outer regions of the disk and allow the accretion. These systems are known as accretion disks and are present around lowmass young stars, as well as in all sorts of stellar binary systems and in the nuclei of most active galaxies surrounding massive black holes. Despite their diversity, the structure and dynamics of these objects can be understood by a theory combining shear and turbulence. For decades, since the seminal work of Shakura & Sunyaev (1973) that parameterized the angular momentum transport in terms of an effective viscosity, there has been an extensive study in the literature on accretion disk physics, particularly in regard to the mechanisms that may drive an efficient () outward transport of the angular momentum in a differentiallyrotating disk. Since the early 1990’s, the magnetorotational instability (MRI) has been recognized as the potentially dominant process, far more efficient than any other transport mechanism (see, e.g., Chandrasekhar, 1960; Balbus & Hawley, 1991; Hawley et al., 1995; Balbus & Hawley, 1998). Besides the MRI, other magnetohydrodynamic (MHD) processes can play an important role in accretion disks, for instance, in the building and evolution of a hot, lowdensity magnetized corona around them. In particular, the ParkerRayleighTaylor instability (PRTI, see Parker, 1966) is an important mechanism to drive, under some conditions, the formation of magnetic loops arising from the disk out of an initially horizontal magnetic field (see also the role of the PRTI in the formation of magnetic loops in the solar corona in Parker, 1955; Isobe et al., 2005, and references therein).
The role of magnetic fields and their geometry in accretion disk systems, as well as their formation via dynamo action or the advection from the neighborhood, are also long standing problems. Besides extensive numerical studies of the MRI considering a zero net magnetic field flux which provides a low angular momentum transport (, see Brandenburg et al., 1995; Davis et al., 2010; Simon et al., 2011, 2012), the inclusion of a net vertical field flux (Hawley et al., 1995; Bai & Stone, 2013; Salvesen et al., 2016a, b) predicts accretion rates which are more consistent with those expected from the observations (, see King et al., 2007).
In recent years, motivated by the fact that strongly magnetized accretion disks could be found in AGNs and binary systems, Johansen & Levin (2008) performed local numerical simulations considering strong azimuthal fields (with no initial net vertical magnetic flux) and found that the PRTI is important in the vertical transport of the magnetic field from the disk to the outer parts of the system via magnetic buoyancy and also could keep strong magnetization inside the disk. On the other hand, Salvesen et al. (2016a, b) using more realistic boundary conditions verified that the disk cannot sustain a magneticpressure dominated disk in the case of outflow boundary conditions with no initial net magnetic vertical flux, and stressed the importance in adopting more realistic boundaries in order to have a better understanding of the fundamental physical processes that regulate the simulations.
Observations can also help to impose important constraints on the accretion phenomena. The complex emission features from disk systems indicate the existence of different regimes of accretion and frequently a hot, lowdensity magnetized disk corona is invoked in order to explain nonthermal highenergy emission in these systems, such as the well known Xray transitions observed in black hole binaries (BHBs, see, e.g., Fender et al., 2004; Belloni et al., 2005; Remillard & McClintock, 2006; Kylafis & Belloni, 2015). These transitions are characterized by a Low/Hard state attributed to the inverse Compton process in a geometrically thick optically thin accretion flow at a subEddington regime (see Esin et al., 1997, 1998, 2001; Narayan & McClintock, 2008), also known as ADAF (advectiondominated accretion flow, see Narayan & Yi, 1994, 1995; Abramowicz et al., 1995). These systems also show a High/Soft state generally attributed to the heating of a geometrically thin, optically thick accretion disk (Shakura & Sunyaev, 1973) at near Eddington regime and, between these two states, there is a transient one (of the order of a few days, Remillard & McClintock, 2006), whose origin is not fully understood yet, but could be explained by shocks in a jet (e.g., Romero et al., 2003; BoschRamon et al., 2005; Piano et al., 2012) or reconnection in the corona (e.g., de Gouveia Dal Pino & Lazarian, 2005; de Gouveia Dal Pino et al., 2010a, b; Kadowaki et al., 2015).
Magnetic reconnection in accretion disk/corona systems as a mean to explain flaring emission has been explored in several analytical and numerical works (see, e.g., de Gouveia Dal Pino & Lazarian, 2005; Igumenshchev, 2009; Soker, 2010; Uzdensky & Spitkovsky, 2014; Dexter et al., 2014; Huang et al., 2014; Kadowaki et al., 2015; Singh et al., 2015; Khiali et al., 2015). de Gouveia Dal Pino & Lazarian (2005), for instance, proposed an analytical model to explain the peculiar flaring state of the BHB (also named microquasar) GRS+, which was later extended to other few BHBs, AGN sources and protostars (see de Gouveia Dal Pino et al., 2010a, b). According to this model, fast reconnection between the magnetic field lines of the source’s magnetosphere and those of the accretion disk can produce plasmoids and accelerate particles to relativistic velocities via a firstorder Fermi process in the current sheets (see de Gouveia Dal Pino & Lazarian, 2005; Drake et al., 2006, 2010; Zenitani & Hoshino, 2008; Zenitani et al., 2009; Kowal et al., 2011, 2012; del Valle et al., 2016; Guo et al., 2016, and references therein), providing enough power to explain observed nonthermal flaring emission. More recently, employing a similar model, but with the fast reconnection driven by the background turbulence in the coronal plasma permeating the magnetic field lines, Kadowaki et al. (2015) and Singh et al. (2015) verified that this could explain the very highenergy emission of hundreds of lowluminosity (nonBlazar) AGNs and BHBs. A similar process has been proposed for energy extraction in the vicinity of Kerr black holes by Koide & Arai (2008). Magnetic reconnection of the field lines generated by buoyancy processes like the PRTI, has been also invoked to be an efficient process to heat the coronal region of accretion disks around black holes by Liu et al. (2002, 2003) and Huang et al. (2014).
The efficiency of the reconnection in such cases is a key ingredient to allow for highly variable, explosive emission, as in the case of the solar flares which show magnetic reconnection velocities in a range between , where is the Alfvén speed (see, e.g., Dere, 1996; Aschwanden et al., 2001; Su et al., 2013). Several mechanisms have been invoked to describe magnetic reconnection. Since the proposed SweetParker model (Sweet, 1958; Parker, 1957) that predicts very slow reconnection rates:
(1) 
where is the reconnection velocity of opposite magnetic fluxes, is the Lundquist number, the length of reconnection site, and the Ohmic resistivity, other models have been proposed in order to explain observed fast reconnection, like the Petschek’s Xpoint model. In this model, reconnection is forced to occur in single points, rather than in an entire flat current sheet (with a reconnection rate , see Petschek, 1964). This increases the reconnection rate, but Biskamp et al. (1997) showed numerically that this mechanism is unstable and makes the system evolve rapidly to a SweetParker configuration, unless the flow is collisionless and has localized variable resistivity. The main difficulties in explaining fast reconnection in collisional flows have been removed with the proposal of the model of Lazarian & Vishniac (1999) that considers the ubiquitous presence of turbulence in astrophysical environments to drive fast reconnection. The idea behind this model is that the wandering of the magnetic field lines in a turbulent flow allows for several patches to reconnect simultaneously making reconnection fast (Lazarian & Vishniac, 1999):
(2) 
where and are the velocity and size of the turbulent eddies at injection scale), and independent of the background resistivity, so that even nearly ideal MHD flows that have nearly negligible Ohmic resistivity, as most astrophysical flows, can undergo fast reconnection when there is turbulence. Extensive 3D numerical MHD testing of this model has been successfully performed in current sheets with embedded forced turbulence (e.g., Kowal et al., 2009, 2012; Takamoto et al., 2015), but it has never been thoroughly investigated in systems with natural driving sources of turbulence, like the MRI and the PRTI in accretion disks (see, however, tentative exploration of fast reconnection driven by currentkink instability turbulence in 3D MHD relativistic jet simulations by Singh et al., 2016).
Magnetic reconnection has also been detected in the coronal regions of global numerical simulations of accretion disks both in nonrelativistic systems (see, e.g., Romanova et al., 2002, 2011; Zanni & Ferreira, 2009, 2013), and general relativistic systems (e.g., McKinney et al., 2012; Dexter et al., 2014; Parfrey et al., 2015).
Motivated by the studies above that have highlighted the potential importance of turbulence and fast magnetic reconnection in magnetized accretion disk coronal flows to explain observed phenomena including particle acceleration and nonthermal flare emission features in compact sources, our main goal here is to explore these processes quantitatively in depth, which require very high resolution numerical simulations. For this reason, we have performed local 3DMHD numerical simulations employing a shearingbox approximation (Hawley et al., 1995) starting with initially strong horizontal magnetic fields. We could assess the role of the PRTI and MRI in the development of turbulence, large scale magnetic loops in the corona, and fast magnetic reconnection driven by turbulence. We have also used a modified algorithm based on the work of Zhdankin et al. (2013) and Kowal et al. (2009) to identify reconnection sites and evaluate statistically the magnetic reconnection rates in the accretion disk and corona.
The paper is organized as follows, in section 2 we describe the numerical method for the shearingbox approach (Hawley et al., 1995) and the initial and boundary conditions used in this work (see, e.g., Johansen & Levin, 2008; Davis et al., 2010; Bai & Stone, 2013); in section 3 we show the results obtained from the evolution of the averages of the physical quantities taken over the accretion disk and coronal regions, and compare with previous works (see, e.g., Miller & Stone, 2000; Johansen & Levin, 2008; Salvesen et al., 2016a, b); in section 4 we show the results of the identification of fast reconnection driven by turbulence in the disk and corona, as well as the statistical analysis of these sites and the corresponding reconnection rates. Finally, in section 5 we discuss these results and draw our conclusions.
2 Numerical method
The numerical solutions of our simulations are obtained from the magnetohydrodynamics equations (MHD) which describe the macroscopic behavior of a magnetized fluid. These equations, in the conservative and ideal form, are given by:
(3) 
(4) 
(5) 
and correspond to the equations of mass and momentum conservation, and induction, respectively. In these equations, is the density, is the thermal pressure of the gas, is the velocity field, is the identity matrix, and is the magnetic field. The equations have been solved in nondimensional form, thus without a factor . In addition, corresponds to the gravity acceleration and to the gravitational potential. In the present work, we adopt an isothermal equation of state with , where is the sound speed.
In this work, we use a shearingbox approach (see Hawley et al., 1995), useful to obtain the statistical properties of very small regions of accretion disks systems. Such approach consists in a linear shearing velocity in the azimuthal direction “” given by:
(6) 
where is the angular velocity at an arbitrary radius , and is the shear parameter for a Keplerian profile.
From the shearingbox approach (in the disk reference frame at ), the right side of the momentum equation (4) is rewritten as:
(7) 
where the source term “” (for ) corresponds to Coriolis forces, “” corresponds to the tidal gravity, and “” corresponds to the vertical gravity due to the vertical density stratification.
We have used the ATHENA code to obtain the numerical solution of the 3DMHD equations with the shearingbox approach (see Stone et al., 2008, 2010; Stone & Gardiner, 2010), for a Cartesian coordinate system. To compute the intercell fluxes of the computational grid, a HLLD Riemann solver has been employed (see Miyoshi & Kusano, 2005), while a secondorder RungeKutta scheme has been used to solve the equations in time. An orbital advection scheme has also been adopted (Stone & Gardiner, 2010), where the azimuthal component of the velocity () is split into an advection part () and another one involving only fluctuations () given by:
(8) 
This scheme improves the simulations, since the numerical integration of the advection part is not subject to the CourantFriedrichLewy (CFL) condition, making the time step () less restrictive (see, Johnson et al., 2008; Davis et al., 2010; Stone & Gardiner, 2010).
We have used shearingperiodic boundary conditions in the radial direction “” that reproduce the differential rotation through the azimuthal displacement of the radial boundaries (see Hawley et al., 1995). In the azimuthal direction “” we have used standard periodic boundary conditions. Finally, we have adopted outflow boundaries in the vertical direction “” 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, zerogradients are used for the velocity and magnetic fields, but the vertical velocity component is set zero () whether inflows are detected through the boundaries. Besides, the density is extrapolated assuming a vertical hydrostatic equilibrium.
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):
(9) 
where is the initial density in midplane of the accretion disk, 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 , and the thermal scale height as our unit of length (), that yields to . Assuming a constant , the magnetic field profile is given by (see also Johansen & Levin, 2008):
(10) 
where we have taken . It is important to emphasize that the criteria to trigger the PRTI will be slightly different compared with the wellknown work of Parker (1966), since we are dealing with a differential rotation system under a linear gravity. The presence of a differential rotation was studied analytically by Foglizzo & Tagger (1994, 1995), where they 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, more appropriate to the case of Keplerian disks and the present work. However, in all the works, an initial of the order of the unit is still required to trigger the PRTI.
We have assumed the initial velocity field as zero (), since an orbital advection scheme has also been adopted (see eq. 8 and Stone & Gardiner, 2010). However, a Gaussian perturbation with a amplitude of has been used in the azimuthal velocity component to trigger the linear phase of the PRTI and MRI (see Johansen & Levin, 2008). Finally, a density floor () has been applied to minimize the problems with the CFL condition (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 unit. 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^{1}^{1}1Different to the work of Gressel et al. (2011), our implementation does not conserve mass and for this reason, a density floor has been used. and allows us to use for the low and intermediate resolution and for the high resolution (see next section for more details).
3 Numerical results
We have performed a set of 3DMHD shearingbox simulations of accretion disks employing the ATHENA code (see Stone et al., 2008, 2010). In this system, the PRTI and MRI are triggered by initial perturbations, as described in section 2. We have used a resolution of (), (), and () cells per thermal height scale of the system (), considering a computational domain of size . We also considered three different values of (, , and ) to evaluate the evolution of both instabilities in the system. The parameters of the simulation are shown in Table 3 and each model name is composed by 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.
lccccccc
[2ex] Simulation & Computacional domain & Resolution & & & & & Vertical boundaries
[0.5ex]
[1.8ex]
\endfirstheadTable 0  Continued
[0.5ex]
[2ex]
Simulation & Computacional domain & Resolution & & & & & Vertical boundaries
[0.5ex]
[1.8ex]
\endheadContinue in the next page…
\endfoot\endlastfoot
R11b1 & & & & & & & Outflows
R21b1 & & & & & & & Outflows
R21b10 & & & & & & & Outflows
R21b100 & & & & & & & Outflows
R43b1 & & & & & & & Outflows
3.1 Diagnostics
In order to follow the evolution of the system and the instabilities, we will present here time and space average diagrams of relevant quantities. The space average of a given variable taken at the entire volume of the box is represented as , so that:
(11) 
Also, in order to track the evolution of at the vertical direction “”, we performed space averages at “” plane, at each height “” (horizontal averages), represented by the symbol , where:
(12) 
Finally, time averages are indicated by a index, so that time volume and horizontal averages are represented by and , respectively.
Another useful quantity is the angular momentum transport of the component (or in the cylindrical geometry) of the stress tensor:
(13) 
where is the Maxwell (magnetic) stress tensor, and is the Reynolds stress tensor. The total stress tensor is correlated with the angular momentum transport parameter (Shakura & Sunyaev, 1973) by the relation:
(14) 
We have also evaluated the contribution of the turbulent magnetic field component () in the evolution of the system. To determine its intensity we subtracted in each cell the horizontal average from the total magnetic field :
(15) 
To evaluate the power spectrum of the velocity field (), we performed a twodimensional Fourier transformation to obtain the and wavenumbers in each height of the domain:
(16) 
where and are the total number of cells in the radial and azimuthal directions, respectively. The spectrum has been obtained by integrating over circular shells with radius .
We adopt the R21b1 model (see table 3) as our reference one. We describe the main results of this simulation in the next section.
3.2 R21b1 model
Model R21b1 shows an initial configuration with an azimuthal magnetic field in a regime with to trigger the PRTI (see eq.10, section 2).
Figure 1 shows the initial profiles of the density distribution and magnetic field lines (upper diagram) and their evolution after (bottom diagram, where corresponds to the orbital period of the accretion disk), for a resolution of cells per scale height (, or cells, see Table 3). The configuration of the magnetic field lines at shows that the azimuthal component “” is stretched in the vertical direction “” due to the effects of magnetic buoyancy, producing loops during the exponential growth of the PRTI which attain larger amplitudes at the higher latitudes of the computational domain with a wavelength (in agreement with the results obtained by Johansen & Levin, 2008).
Figure 2 shows the evolution of the mean magnetic energy density (top diagram), and the Maxwell and Reynolds stress tensors and (bottom diagram) over orbital periods. The top diagram indicates an amplification of the and components during the first orbital periods. Estimating the exponential growth time from the linear perturbation solutions for the PRTI obtained by Kim et al. (1997) for the higher latitudes of the system () we find values () that are in agreement with our simulations, although Kim et al. (1997) have neglected differential rotation in their evaluation. After orbital periods, the total magnetic field decreases mainly due to expansion and upward transport through the outflow boundaries in the vertical direction “”, and a small dissipation due to turbulent magnetic reconnection probably also contributes to this decrease (see below). After , it achieves a nearly steadystate regime with a total magnetic field that is of its initial intensity. A similar behavior is found in the simulations of Johansen & Levin (2008), but these authors have obtained a smaller reduction in the field, probably due to the fact that they adopted different vertical boundary conditions that prevent the escape of azimuthal flux in that direction. Other simulation studies of stratified shearingboxes considering outflow vertical boundaries have also found that strong azimuthal fields cannot be maintained within the system due to the buoyancy effects upon them (see,e.g., Miller & Stone, 2000; Salvesen et al., 2016b, a). At , the decrease of the magnetic field intensity and consequent increase of induces the growth of the MRI (from the azimuthal and vertical magnetic field components) that will dominate the system evolution for the rest of the simulation.
The bottom diagram of Figure 2 shows that is dominated by the Maxwell stress tensor and follows the same trend of the magnetic energy density. The evolution of also indicates that there is a high accretion regime in the first orbital periods when it achieves a maximum value around (consistent with those expected from the observations, see King et al., 2007). We discuss in more details the behavior of in the following sections.
3.2.1 Accretion disk and corona evolution
In order to evaluate the evolution of the accretion disk and the corona separately, we have obtained volume averages considering three different regions in the computational domain. The region between and we assume to correspond to the accretion disk (see the black line in Figure 3), while the more rarefied regions above and below this zone we define as the upper corona, between and (red line of Figure 3) and the lower corona, between and (blue line of Figure 3), respectively. More specifically, we have considered the corona as the region where the initial density is smaller than . We have excluded the regions from this analysis in order to avoid boundary effects.
Figure 3 shows the evolution of the magnetic energy density, the Maxwell stress tensor, and the ratio between the thermal and magnetic pressures () in different layers of the domain for the same model of Figures 1 and 2. In the first orbital periods, the magnetic energy density in the disk decreases, whereas in the coronal region it increases due to the upward transport of the magnetic field from the central region to the corona. Like the average quantities evaluated for the entire system (Figure 2), after we see that the average magnetic energy density of the disk and the corona decreases. The variability seen in the coronal region corresponds to inversions of polarity produced by the action of a dynamo induced by the MRI, as we will see in Section 3.2.2 (see also, e.g., Davis et al., 2010; Shi et al., 2016).
The increase of with time changes significantly the behavior of the Maxwell stress tensor. Its maximum value is reached between and orbital periods and then decays over time until a saturated minimum value. The contribution of the accretion disk to around is important compared with the contributions of the coronal region. In the coronal region, the intensity of (red and blue lines of the bottom diagram of Figure 3) decreases until due to the growth of the magnetic field strength in this region induced by the PRTI. Thereafter, it saturates around an average value with a strong variability (between and ). This behavior demonstrates that even with an initial highly magnetized system (), the disk evolves to a gaspressure dominated regime, while the corona is magnetically dominated, as it should be.
3.2.2 Transition between PRTI and MRI
Once the thermal pressure becomes larger than the magnetic pressure (, see the last diagram of Figure 3), the MRI is expected to settle in and soon dominate the dynamics of the disk. The MRI evolves out of the azimuthal and also of the vertical component of the magnetic field produced by the PRTI (see also, e.g., Balbus & Hawley, 1992; Hawley et al., 1995; Foglizzo & Tagger, 1995; Johansen & Levin, 2008). Figure 4 highlights the transition between the PRTI and MRI regimes. The diagrams show the evolution of the magnetic field lines and the density distribution between 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 the PRTI (and to outside the computational domain) are observed. These diagrams also show that once the MRI is established, turbulence develops in the system which is a well known result (see, e.g., Hawley et al., 1995).
Figure 5 shows the evolution of the total magnetic energy density (black line) compared with the average (blue line) and fluctuating components (red line). In the first orbital periods, the energy density is dominated by the contribution of the mean field amplified by the loops formed thanks to the PRTI, but thereafter the turbulent component becomes relevant. This behavior coincides with the transition between the PRTI and MRI, which strengthens the conclusion that the turbulence generated by the MRI dominates the dynamics of the system after . The bottom diagram shows that after orbital periods, the ratio between the mean and fluctuating components remains around .
Figure 6 shows the evolution of the magnetic field components as a function of the height of the system. These diagrams are produced from the horizontal averages of each component () obtained at each height “”. The first and second diagrams of the Figure 6 show a change of regimes between and of the radial and azimuthal components of the magnetic field. Besides, during the regime where the PRTI dominates, there is an increase of the radial component with a peak of intensity in the middle of the disk, followed by an inversion of polarity in the coronal regions between and . In the higher latitudes of the corona, between and , there are smaller polarity inversions. The bottom diagram of Figure 6 shows the evolution of the vertical magnetic field component (), for which the horizontal averages result a nearly zero value (), though the volume average grows continuously with time during the PRTI regime (as indicate in Figure 2).
After , the diagrams of Figure 6 show that the radial and the azimuthal components of the magnetic field undergo polarity inversions on a time scale of approximately orbital periods. This is also a well known result consistent with previous studies of MRI in weak field regimes () and zero net flux (e.g., Brandenburg et al., 1995; Davis et al., 2010; Simon et al., 2011, 2012), which have produced similar butterfly diagrams to those we see in Figure 6. This pattern suggests the action of an alphaOmega 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^{2}^{2}2Note that since the mean vertical is null, it should not effectively participate in the amplification of the largescale components, except possibly at the beginning of the evolution of the dynamo process in the PRTI phase. in the azimuthal direction. The alpha effect allows for the amplification of the radial component “” from the azimuthal component “”, through the electromotive forces in the azimuthal direction , where and correspond to the fluctuations of velocity and magnetic field, respectively. We also note that in Figure 6, the polarity inversions of and components at the end of each half cycle (where the old field is replaced by a new one with 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, b; Guerrero & de Gouveia Dal Pino, 2008; Guerrero et al., 2009, 2016a, 2016b, and references therein).
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 the unit 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 dynamo process early in the PRTI regime.
3.3 Comparison between models with different initial values of
In this section, we show the evolution of the PRTI and MRI for different initial values of (, and ) in a computational domain of size and a resolution of (models R21b1, R21b10, and R21b100, respectively). The presence of a weak magnetic field is essential for the development of the MRI (even for azimuthal fields, see, e.g., Chandrasekhar, 1960; Balbus & Hawley, 1991; Hawley et al., 1995; Miller & Stone, 2000; Salvesen et al., 2016b, a). In contrast, the PRTI evolves in regimes where the magnetic pressure is of the order of the thermal pressure (; 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 7 show the evolution of the volume averages of the magnetic energy density (top diagram), the Maxwell stress tensor (second diagram), and (bottom diagram) for different initial values of ^{3}^{3}3Note that the “” index indicates the initial value of .. The reference model R21b1 with discussed in the previous sections is also shown for comparison (black line).
The diagrams of Figure 7 show clearly that the initial evolution of the magnetic field is different for models with distinct . Nevertheless, after orbital periods all systems evolve in a similar way achieving a nearly steady state situation. For R21b10 model (with ), 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 R21b100 model (with ) shows the standard evolution of the magnetic energy density driven by the MRI only, with an initial amplification until the saturation is achieved after orbital periods (e.g., Hawley et al., 1995; Miller & Stone, 2000).
The Maxwell stress tensor also shows convergence for all models after orbital periods (second diagram at the middle of Figure 7). The initial peak value achieved during the growth phase of the PRTI seen in R21b1 model gets smaller for increasing , as expected, and absent for .
The bottom diagram of Figure 7 shows the time evolution of the volumeaveraged value of in the disk (continuous line), upper (shortdashed line) and lower corona (longdashed line) 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 leads to a decrease in and then, after 20 orbital periods all 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 gaspressure dominated system surrounded by a much more magnetized corona (with ), which is in agreement with Miller & Stone (2000); Salvesen et al. (2016a). We note, however, that Johansen & Levin (2008) obtained a more magnetized disk () that was sustained for 2030 orbital periods, most likely due to their different adopted vertical boundary conditions which allowed for confinement of magnetic flux inside the domain. As stressed in Section 1, recently, Salvesen et al. (2016a) have investigated the influence of a different vertical boundary condition in shearingbox simulations of initially strongly magnetized accretion disks, and found that the disk cannot maintain its strongly magnetized initial state in the case of outflow boundary conditions, which is consistent with the present results. Finally, the differences found in all diagrams of Figure 7 show that the PRTI has an important role just in an early “transient” phase.
3.4 Numerical convergence
As 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, specially the modes (see Hawley et al., 1995). Besides, considering a zero net vertical flux in stratified shearingbox 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). On the other hand, in stratified shearingbox simulations with a strong azimuthal magnetic field, the size of the time step is strongly constrained by the CFL condition as a consequence of the highly magnetized corona that is built by the PRTI and magnetic buoyancy, independent of the initial , as seen in the previous section. In this work, a set of procedures have been taken to minimize this problem, such as the application of the Fargo scheme with a density floor (see, e.g., Davis et al., 2010; Bai & Stone, 2013), and a mass diffusion term in the mass conservation equation (see Gressel et al., 2011), as described in section 2. Other procedures can be used like the socalled Alfvén speed limiter adopted by Miller & Stone (2000) and Johansen & Levin (2008). 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 8 shows the time evolution of the magnetic energy density (top diagram) and the Maxwell stress tensor (bottom diagram) for the resolutions of (R11b1 model, black line), (R21b1 model, red line), and (R43b1 model, blue line). The highresolution simulation (R43b1) is numerically costly and we have evolved over a short time equivalent to orbital periods. During the first orbital periods (where the PRTI is dominant) until , all the resolutions show the same behavior for the magnetic energy density, indicating that the PRTI operates appropriately even in our lowerresolution 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 would not affect significantly the initial evolution of such instability. Nevertheless, the maximum value of for this model is slightly lower than for the intermediate and highresolution simulations (R21b1 and R43b1 models, respectively), indicating that resolutions lower than could affect the initial evolution of the system.
After orbital periods, when the MRI dominates, the differences in the evolution between the low and intermediate resolution models is significant. The lowresolution model continues to decrease the magnetic field without achieving saturation. Such 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 behavior is due to 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 highresolution models indicate a good agreement between them until the evolved time of the highresolution model R43b1 (around orbits), though it is not possible to determine yet whether it will achieve the same saturation as in R21b1 model.
4 The search for magnetic reconnection
As demonstrated, the PRTI is fundamental in the formation of magnetic loops in the coronal region of the accretion disk. During this process, the encounter and squeezing between loops may lead to magnetic reconnection. Indeed, magnetic reconnection is expected to occur whenever two magnetic fluxes of opposite polarity approach each other in the presence of finite magnetic resistivity. The ubiquitous microscopic Ohmic resistivity is enough to allow for magnetic reconnection, though in this case the rate at which the lines reconnect is expected to be very slow according to the SweetParker mechanism (Sweet, 1958; Parker, 1957). In the present analysis, we are dealing with ideal MHD simulations with no explicit resistivity. This in principle would prevent us from detecting magnetic reconnection. Nevertheless, in numerical MHD simulations magnetic reconnection can be excited because of the presence of numerical resistivity. This in turn, could make one to speculate that the identification of potential sites of magnetic reconnection in ideal MHD simulations would be essentially a numerical artifact. Nevertheless, according to the LazarianVishniac reconnection model (Lazarian & Vishniac, 1999; Kowal et al., 2009; Eyink et al., 2011), the presence of turbulence in real MHD flows is able to speed up the reconnection to rates nearly as large as the local Alfvén speed. This occurs because of the turbulent wandering of the magnetic field lines that allow them to encounter each other in several patches simultaneously making reconnection very fast. Even in subAlfvénic flows, where the magnetic fields are very strong, this fast reconnection process may be very efficient. We have seen here that the PRTI and the MRI are able to trigger turbulence in the flow and also the close encounter of magnetic field lines with opposite polarity in several regions and specially at the coronal loops. We will see below that these regions are loci of large increase of the current density in very narrow regions and are, therefore, potential sites of magnetic discontinuities or current sheets. The numerical resistivity here simply mimics the effects of the finite Ohmic resistivity, but the important physical real mechanism that may allow for the fast reconnection is the process described above due to the turbulence. This process has been extensively and successfully tested numerically in several studies involving ideal MHD simulations of turbulent flows (see e.g., the references above and the reviews of Lazarian et al., 2012, 2015; de Gouveia Dal Pino & Kowal, 2015).
4.1 Method
In order to identify magnetic reconnection sites in our disk/corona system, we have searched for current density peaks () and then evaluated the local magnetic reconnection rate () 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 , where 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 (, where corresponds to the index position of each cell) within a surrounding cubic subarray of the data (with a size of cells) and then we checked whether they are located between magnetic field lines of opposity polarity (in each component separately). In the present work, we are not interested to identify 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) have used an Xpoint 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 () for each cell of the last subsample (see Zhdankin et al., 2013):
(17) 
where corresponds to the secondorder 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 . We have assumed this direction as the thickness of the magnetic reconnection site. The eigenvectors of provide the three orthonormal vectors of the new coordinate system centered in the local reconnection site (see the left sketch of Figure 9).
The local magnetic reconnection rate has been evaluated in a similar way as in Kowal et al. (2009), where we averaged the inflow velocity (, the projection onto the direction) divided by the Alfvén speed at the edges of the reconnection site (see the second sketch of Figure 9):
(18) 
where the Alfvén speed is given by:
(19) 
and , , and 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 greater than half of the maximum local value given by (, see the second sketch of Figure 9). We have constrained the subsample to obtain for the most symmetric profiles, since it is expected a symmetry of both the magnetic and velocity field profiles around the magnetic reconnection site (as seen in Figure 9). As an example, Figure 10 shows the spatial distribution of the local maxima identified by the algorithm, in the coronal region (upper corona, ) at 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 10 show examples of rejected and accepted profiles of , and interpolated along the axis. This constraint reduces considerably the subsample size, but also reduces an over or underestimation of .
4.2 Magnetic field configuration
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’s flips of the magnetic field in each direction separately. These distributions can help to 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 11 shows such distributions (with bins of orbital period) for the model R21b1, where corresponds to the number of sign’s flip of the magnetic field component “” in the direction “”. The upper and lower diagrams correspond to the distribution taken in the coronal region and disk, respectively. Initially, between and orbital periods at the coronal region, the number of sign’s flips is dominated mainly by the component along the vertical direction “” () and the component along the radial direction “” (). 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 Simon et al., 2012; Hirose et al., 2006) and allowing for reconnection to occur in the “” and “” directions mainly. This behavior can be clearly seen in the first diagram of Figure 4. On the other hand, after , the number of sign’s flip of the component in the “” and “z” directions ( and , respectively) increases playing an important role in the reconnection too, which is compatible with the inversion of polarity of the azimuthal magnetic field during the MRI phase (as seen in the middle diagram of Figure 6). The sign’s flips in the azimuthal direction ( and ) are less relevant during all the evolution of the simulation, as expected, since the shear and stretching prevent local encounters of magnetic field lines of opposite polarity.
Similar behavior can be found in the disk region, except for the component in the “” direction (), for which the flips are less relevant than in the coronal region. During all the simulation, the sign’s flips are dominated by , and . This is not a surprise since the component is produced mainly by loops at the higher latitudes above and below the disk during the exponential growth of the PRTI (as described in section 3).
4.3 Magnetic reconnection rate
The distribution of was obtained for each time step of the simulations and Figures 12 depict the time evolution of the average value of for the upper and lower coronal regions together and the disk (continuous black line). These averages were calculated out of histograms computed over orbital periods. The diagrams also depict the standard deviation (blue shade) of each distribution.
In the corona, increases very fast in the early PRTI regime (at the first orbital periods) and then saturates (after orbital periods) to a mean value varying between and . In the disk, also increases fast in the PRTI regime, achieving a peak value () around orbital periods, then decreases to nearly constant value () after when the MRI sets in, following the same trend of the time evolution of the Maxwell stress tensor and magnetic energy (Figure 3).
These values of are in agreement with the predictions of the theory of fast magnetic reconnection driven by turbulence (Lazarian & Vishniac, 1999) and the numerical studies that have tested it (see, e.g., Kowal et al., 2009; Takamoto et al., 2015; Singh et al., 2015). They are also compatible with observations of fast magnetic reconnection in the solar corona associated to flares ( and up to , see, e.g., Dere, 1996; Aschwanden et al., 2001; Su et al., 2013). This indicates that the algorithm used here can to 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 the PRTI and MRI inducedturbulence to drive fast reconnection in accretion disks and their coronae.
Figure 13 shows the histograms of (diagrams on the left side) and the measure of the thickness of the reconnection regions (diagrams on the right side) for R21b1 model obtained between and orbital periods in the coronal region (upper diagrams) and disk (lower diagrams), where the reconnection rate, as the other quantities, achieves a nearly steadystate. In both cases, the distributions of 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 , but we have constrained the histogram to the range of between to which corresponds to of the total data (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 13 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 and a median of , whereas in the disk we have obtained an average of and a median of . The averages values of the thickness show that the reconnection sites occupy or cells (for the resolution of ), therefore numerical effects could affect the evaluation of . We have also performed tests with larger reconnection sites, but the averages of did not change significantly.
4.4 Correlation with turbulence
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 twodimensional power spectrum of the velocity field ^{4}^{4}4The power spectrum was obtained from the velocity field without the advection term “”, computed by the fargo scheme (see section 2). () for the and wavenumbers in different heights (“” direction) inside the computational domain. The power spectrum has been taken from the average in circular shells with radius in the Fourier space, as described in section 3.1 (see eq.16). We have also averaged in time with bins of orbital period to reduce the fluctuations of the spectra. Figure 14 shows the power spectrum compensated by a factor of in three different times, at , and , where the colors of each line correspond to the height ””. As we expected, the turbulence is not fully developed in the early stages of the simulation. At the power is still very small and the magnetic reconnection rate has values below (see Figure 12). On the other hand, after orbital periods the turbulent power spectrum is much larger and shows the typical cascading to small scales with a possible inertial range that approximates a Kolmogorov spectrum (dashed line).^{5}^{5}5We note that the behavior of the power spectrum for the MRI (and PRTI)driven turbulence is not fully understood yet. For instance, recently Walker et al. (2016), considering shearingbox simulations, have decomposed the energy power spectrum into a parallel (largescale shearaligned) and a perpendicular (smallscale fluctuation) component to the mean magnetic field. With this procedure, they obtained power law spectra close to and for the parallel and perpendicular components, respectively. A similar decomposition is out of the scope of the present work, but their results for the small scale, turbulent fluctuations are compatible with our results. Around this time, during the PRTI, the magnetic reconnection rate achieves its largest values between and in the disk and 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 in such regions (where the average value in the coronal region is higher than in the disk after orbital periods), indicating a correlation between the turbulence and .
4.5 Comparison between models
Finally, Figure 15 compares the time evolution of 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 12), the left diagrams show that the behavior of is similar to the one of the reference model R21b1, and consistent with the time evolution of these systems as discussed in Section 3.3. In the first orbital periods (before ), 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 orbital periods, regardless of the initial , converges to the average values discussed previously (between and ). At this stage, the volume averages of the magnetic energy density, the Maxwell stress tensor and also converge to a steadystate (as seen in Figure 3 and 7). The comparison of with different resolutions reveals a good convergence mainly in the disk region (bottom diagram in the right side of Figure 15), whereas in the coronal region (top diagram in the right side) the mean values of for the highresolution models are slightly smaller than for the lowresolution.
The right diagrams also show that the variability amplitude of decreases for the models with higher resolution. Above orbital periods, the model R43b1 shows values varying between and in the coronal region, and between and in the disk. On the other hand, the model R11b1 shows values between and in the coronal region, and between and in the disk. This is not a surprise, since it is harder to resolve the structures of the magnetic reconnection site in the lowest resolution. Table 4.5 summarizes the time average values of the magnetic reconnection rate obtained between and orbital periods, excepted for the highresolution simulation (model R43b1) whose time average was obtained between and orbital periods. Considering the standard deviations, all the models show values which are compatible, both for different values of and resolutions.
{longtable}lccc
[2ex] Simulation & & &
Name & Corona & Disk & Range
[0.5ex]
[1.8ex]
\endfirstheadTable 0  Continued
[0.5ex]
[2ex]
Simulation & & &
& Corona & Disk &
[0.5ex]
[1.8ex]
\endheadContinue in the next page…
\endfoot\endlastfootR11b1 & & & 20P100P
R21b1 & & & 20P100P
R21b10 & & & 20P100P
R21b100 & & & 20P100P
R43b1 & & & 20P58P
5 Discussion and conclusions
In this work, we have performed 3DMHD shearingbox numerical simulations (see Hawley et al., 1995) of accretion disks in order to capture with a resolution as high as possible, the long term dynamical evolution of the ParkerRayleighTaylor and magnetorotational instabilities (PRTI and MRI, respectively), and follow the formation of the disk corona, turbulence and magnetic reconnection. To allow for the growth of the PRTI and MRI, we have considered accretion disks with an initial strong azimuthal magnetic field, having , and .
As stressed in Section 1, this study is applicable to accretion phenomena in general, but is particularly relevant for accretion disks around stellar mass and supermassive black holes which are believed to sustain strong magnetic fields at least during certain accretion regimes. These could be produced by dynamo processes and/or by the transport of magnetic fields from a companion star or the surrounding medium (see, e.g., de Gouveia Dal Pino & Lazarian, 2005; Johansen & Levin, 2008; Wielgus et al., 2015; Salvesen et al., 2016a, b, and references therein).
Therefore, numerical studies on the formation of turbulent magnetized coronae and magnetic loops induced by these instabilities can provide important clues about magnetic reconnection events which can have an important role in the heating and particle acceleration in such regions and could explain the radio and gammaray emission observed in BHBs and active galactic nuclei (see de Gouveia Dal Pino & Lazarian, 2005; de Gouveia Dal Pino et al., 2010b; Kadowaki et al., 2015; Singh et al., 2015).
Our main results can be summarized as follows:

The PRTI, which dominates the early evolution of the system due to the small values of , leads to the formation of poloidal magnetic fields and loops which are transported from the midplane of the disk to the higher latitudes, allowing for the formation of a magnetized corona with complex structure and values around unit. The increase of in the disk, on the other hand, allows for the development of the MRI (which becomes dominant after orbital periods in all models), and both instabilities drive turbulence and a dynamo action in the system.

We have employed outflow boundaries conditions in the vertical direction which are more suitable to reproduce real accretion disk coronae. These boundaries allow for the transport of magnetic field to outside of the computational domain that causes the decay of the total magnetic field with time until it achieves a nearly steadystate regime, this due to the generation of new field flux by dynamo process mainly during the MRI regime.

Our systems end up with a gas pressuredominated disk (with ) and a magneticallydominated corona (whith ) which is consistent with the results found in Miller & Stone (2000); Salvesen et al. (2016a). This behavior is not in agreement with the results of Johansen & Levin (2008), who obtain a highly magnetized disk. This is probably due to their adopted vertical boundary conditions which prevent the escape of the azimuthal magnetic flux from the domain. In a more recent work, Salvesen et al. (2016a, b), considering similar boundary conditions (with initial zero net vertical magnetic flux) as in here, have also found that strong azimuthal fields cannot be maintained for long periods within the disk due to the magnetic buoyancy effects, unless the system has sufficient initial net vertical flux.

We have tested the numerical convergence of our results with different resolutions (, and ) and found a good agreement between the intermediate (our reference model, R21b1) and high (R43b1) resolution models which reach the steadystate regime at the same time when the MRI becomes selfsustained inside the system. The lowresolution model (R11b1), on the other hand, shows a continuous slow temporal decrease of the magnetic field after orbital periods without achieving a steadystate indicating that this model did not achieve appropriate resolution.

The values found for the Maxwell stress tensor in the steadystate regime () are smaller than those expected from observations. For instance, considering the viscous time scales in the most active (outburst) states of FU Orionis systems, these give a viscosity parameter of the order of (see Zhu et al., 2007). On the other hand, these values are compatible with earlier numerical studies of homogeneous and stratified systems that have imposed initial zero net vertical fields to trigger the MRI (see Hawley et al., 1995; Stone et al., 1996; Fromang & Stone, 2009; Davis et al., 2010), rather than an initial azimuthal magnetic field as in here. We obtained a large only during the PRTI phase, in the first orbital periods. Bai & Stone (2013) have demonstrated that a net vertical magnetic field applied to the simulations (generating stronger largescale fields) can increase the values of and explain the observational estimates discussed above, although they do not explain how these fields could arise naturally. In contrast, in our simulations, we find that the presence of the PRTI induces both the formation of poloidal fields and the increase of the Maxwell stress tensor to the observed values. We can speculate that this short period in which the PRTI grows and largescale poloidal fields develop could be related with flare events in accretions disk systems (as argued, e.g., in de Gouveia Dal Pino & Lazarian, 2005; de Gouveia Dal Pino et al., 2010b; Kadowaki et al., 2015; Singh et al., 2015).

Finally, the arising of magnetic loops due to the PRTI followed by the development of turbulence due both to the PRTI and MRI produce current density peaks in the coronal region and disk, indicating the presence of magnetic reconnection. To track this process, we employed a modified version of the algorithm developed by Zhdankin et al. (2013) to identify current sheets (with strong current density) produced by the encounter of magnetic field lines of opposite polarity in the turbulent regions of the computational domain. From this analysis, we evaluated the magnetic reconnection rates employing the method adopted by Kowal et al. (2009). Despite the high standard deviations derived from the method, we have found peak values for the reconnection rate () of the order of , and average values of the order of in the accretion disk and in the coronal region (for our reference model, R21b1), indicating the presence of fast magnetic reconnection events, as predicted by the theory of turbulenceinduced fast reconnection of Lazarian & Vishniac (1999). Regarding the histograms of , they do not resemble a normal distribution as they exhibit a tail at higher velocities. This is probably due to the limitations of the method employed which evaluates the reconnection rate only at two points at the edges of the magnetic reconnection site along the axis of the fastest decaying of the current density (obtained from the Hessian matrix). In future work, we intend to apply different methods to compare with the current results. Kowal et al. (2009), for instance, besides evaluating the magnetic reconnection rate with the same method used here, tried also another one considering the time derivative of the magnetic flux^{6}^{6}6Zhdankin et al. (2013) have also identified regions considering changes in the magnetic flux function, but from saddle points in a slice of the computational domain.. However, this is hard to apply in our simulations since the magnetic reconnection sites can move in space or change the direction with time (other methods to be considered include, e.g., Greco et al., 2008; Servidio et al., 2011).
Despite the limitations of the method, the observations of flares in the solar corona indicate magnetic reconnection rates in a range between (see, e.g., Dere, 1996; Aschwanden et al., 2001; Su et al., 2013) and strengthen our results. Numerical simulations of turbulent environments point to similar reconnection rates (see, e.g., Kowal et al., 2009; Takamoto et al., 2015; Singh et al., 2015), indicating that the algorithm used in the present work could be a useful tool for the identification of magnetic reconnection sites in numerical simulations. Furthermore, these results have important implications for the understanding of fast magnetic reconnection processes, flaring and nonthermal emission in accretion disks and coronae. In particular, as remarked in Section 1, recent observations of very rapidly variable, high energy emission associated to compact sources like Xray binaries and low luminosity AGNs have been interpreted as possibly due to fast reconnection in the coronal regions of these sources (e.g., de Gouveia Dal Pino et al., 2010a, b; Kadowaki et al., 2015; Singh et al., 2015; Kushwaha et al., 2017), so that our results offer some support to these studies (see also Asenjo & Comisso, 2017, about reconnection in Kerr spacetime). In forthcoming work, we plan to extend the present study exploring the formation of turbulent magnetized accretion disk coronae carrying out global simulations of these systems.
It is important to stress that the reconnection of the lines has been possible in our ideal MHD simulations because of the underlying numerical resistivity that mimics a small ohmic resistivity, while the presence of turbulence makes it fast (Lazarian & Vishniac, 1999). This is distinct from previous works that explored the effects of an explicit large resistivity () and viscosity (), but still keeping the Prandt number of the order of unit (, as in Fromang & Stone, 2009, using nonideal MHD simulations in shearingbox simulations).
Furthermore, we have dealt with isothermal shearingbox simulations, but an extension of the analysis to a nonisothermal approach could be interesting since thermal diffusivity may play an important role in the dynamics of the system, leading to the expansion of the disk and the development of convection (see Bodo et al., 2012, 2013) that may also have consequences on the formation of a hot magnetized corona.
References
 Abramowicz et al. (1995) Abramowicz, M. A., Chen, X., Kato, S., Lasota, J.P., & Regev, O. 1995, ApJ, 438, L37
 Abramowicz & Fragile (2013) Abramowicz, M. A., & Fragile, P. C. 2013, Living Reviews in Relativity, 16, 1
 Aschwanden et al. (2001) Aschwanden, M. J., Poland, A. I., & Rabin, D. M. 2001, ARA&A, 39, 175
 Asenjo & Comisso (2017) Asenjo, F. A., & Comisso, L. 2017, Physical Review Letters, 118, 055101
 Bai & Stone (2013) Bai, X.N., & Stone, J. M. 2013, ApJ, 767, 30
 Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
 Balbus & Hawley (1992) —. 1992, ApJ, 392, 662
 Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1
 Belloni et al. (2005) Belloni, T., Homan, J., Casella, P., et al. 2005, A&A, 440, 207
 Biskamp et al. (1997) Biskamp, D., Schwarz, E., & Drake, J. F. 1997, Physics of Plasmas, 4, 1002. http://scitation.aip.org/content/aip/journal/pop/4/4/10.1063/1.872211
 Bodo et al. (2012) Bodo, G., Cattaneo, F., Mignone, A., & Rossi, P. 2012, ApJ, 761, 116
 Bodo et al. (2013) —. 2013, ApJ, 771, L23
 Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195
 BoschRamon et al. (2005) BoschRamon, V., Aharonian, F. A., & Paredes, J. M. 2005, A&A, 432, 609
 Brandenburg et al. (1995) Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741
 Chandrasekhar (1960) Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
 Davis et al. (2010) Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52

de Gouveia Dal Pino & Kowal (2015)
de Gouveia Dal Pino, E. M., & Kowal, G. 2015, in Astrophysics and Space
Science Library, Vol. 407, Magnetic Fields in Diffuse Media, ed.
A. Lazarian, E. M. de Gouveia Dal Pino, & C. Melioli, 373
 de Gouveia Dal Pino et al. (2010a) de Gouveia Dal Pino, E. M., Kowal, G., Kadowaki, L. H. S., Piovezan, P., & Lazarian, A. 2010a, International Journal of Modern Physics D, 19, 729
 de Gouveia Dal Pino & Lazarian (2005) de Gouveia Dal Pino, E. M., & Lazarian, A. 2005, A&A, 441, 845
 de Gouveia Dal Pino et al. (2010b) de Gouveia Dal Pino, E. M., Piovezan, P. P., & Kadowaki, L. H. S. 2010b, A&A, 518, 5
 del Valle et al. (2016) del Valle, M. V., de Gouveia Dal Pino, E. M., & Kowal, G. 2016, MNRAS, 463, 4331
 Dere (1996) Dere, K. P. 1996, ApJ, 472, 864
 Dexter et al. (2014) Dexter, J., McKinney, J. C., Markoff, S., & Tchekhovskoy, A. 2014, MNRAS, 440, 2185
 Drake et al. (2010) Drake, J. F., Opher, M., Swisdak, M., & Chamoun, J. N. 2010, ApJ, 709, 963
 Drake et al. (2006) Drake, J. F., Swisdak, M., Che, H., & Shay, M. A. 2006, Nature, 443, 553
 Esin et al. (2001) Esin, A. A., McClintock, J. E., Drake, J. J., et al. 2001, ApJ, 555, 483
 Esin et al. (1997) Esin, A. A., McClintock, J. E., & Narayan, R. 1997, ApJ, 489, 865

Esin et al. (1998)
Esin, A. A., Narayan, R., Cui, W., Grove, J. E., & Zhang, S.N.
1998, ApJ, 505, 854
 Eyink et al. (2011) Eyink, G. L., Lazarian, A., & Vishniac, E. T. 2011, ApJ, 743, 51
 Fender et al. (2004) Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105
 Foglizzo & Tagger (1994) Foglizzo, T., & Tagger, M. 1994, A&A, 287, 297
 Foglizzo & Tagger (1995) —. 1995, A&A, 301, 293
 Fromang & Stone (2009) Fromang, S., & Stone, J. M. 2009, A&A, 507, 19
 Greco et al. (2008) Greco, A., Chuychai, P., Matthaeus, W. H., Servidio, S., & Dmitruk, P. 2008, Geophys. Res. Lett., 35, L19111
 Gressel et al. (2011) Gressel, O., Nelson, R. P., & Turner, N. J. 2011, MNRAS, 415, 3291
 Guerrero & de Gouveia Dal Pino (2007a) Guerrero, G., & de Gouveia Dal Pino, E. M. 2007a, A&A, 464, 341
 Guerrero & de Gouveia Dal Pino (2007b) —. 2007b, AN, 328, 1122
 Guerrero & de Gouveia Dal Pino (2008) —. 2008, A&A, 485, 267
 Guerrero et al. (2009) Guerrero, G., Dikpati, M., & de Gouveia Dal Pino, E. M. 2009, ApJ, 701, 725
 Guerrero et al. (2016a) Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016a, ApJ, 819, 104
 Guerrero et al. (2016b) —. 2016b, ApJ, 828, L3
 Guo et al. (2016) Guo, F., Li, H., Daughton, W., Li, X., & Liu, Y.H. 2016, Physics of Plasmas, 23, 055708
 Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
 Hirose et al. (2006) Hirose, S., Krolik, J. H., & Stone, J. M. 2006, ApJ, 640, 901
 Huang et al. (2014) Huang, C.Y., Wu, Q., & Wang, D.X. 2014, MNRAS, 440, 965
 Igumenshchev (2009) Igumenshchev, I. V. 2009, ApJ, 702, L72
 Isobe et al. (2005) Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2005, Nature, 434, 478
 Johansen & Levin (2008) Johansen, A., & Levin, Y. 2008, A&A, 490, 501
 Johnson et al. (2008) Johnson, B. M., Guan, X., & Gammie, C. F. 2008, ApJS, 177, 373
 Kadowaki et al. (2015) Kadowaki, L. H. S., de Gouveia Dal Pino, E. M., & Singh, C. B. 2015, ApJ, 802, 113
 Khiali et al. (2015) Khiali, B., de Gouveia Dal Pino, E. M., & Sol, H. 2015, ArXiv eprints, arXiv:1504.07592
 Kim et al. (1997) Kim, J., Hong, S. S., & Ryu, D. 1997, ApJ, 485, 228
 King et al. (2007) King, A. R., Pringle, J. E., & Livio, M. 2007, MNRAS, 376, 1740
 Koide & Arai (2008) Koide, S., & Arai, K. 2008, ApJ, 682, 1124
 Kowal et al. (2011) Kowal, G., de Gouveia Dal Pino, E. M., & Lazarian, A. 2011, ApJ, 735, 102
 Kowal et al. (2012) —. 2012, PhRvL, 108, 1102
 Kowal et al. (2009) Kowal, G., Lazarian, A., Vishniac, E. T., & OtmianowskaMazur, K. 2009, ApJ, 700, 63
 Kushwaha et al. (2017) Kushwaha, P., Sinha, A., Misra, R., Singh, K. P., & de Gouveia Dal Pino, E. M. 2017, ApJ, 849, 138
 Kylafis & Belloni (2015) Kylafis, N. D., & Belloni, T. M. 2015, A&A, 574, A133
 Lazarian et al. (2015) Lazarian, A., Eyink, G. L., Vishniac, E. T., & Kowal, G. 2015, in Astrophysics and Space Science Library, Vol. 407, Magnetic Fields in Diffuse Media, ed. A. Lazarian, E. M. de Gouveia Dal Pino, & C. Melioli, 311
 Lazarian & Vishniac (1999) Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700
 Lazarian et al. (2012) Lazarian, A., Vlahos, L., Kowal, G., et al. 2012, Space Sci. Rev., 173, 557
 Liu et al. (2003) Liu, B. F., Mineshige, S., & Ohsuga, K. 2003, ApJ, 587, 571
 Liu et al. (2002) Liu, B. F., Mineshige, S., & Shibata, K. 2002, ApJ, 572, 173
 McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
 Miller & Stone (2000) Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398
 Miyoshi & Kusano (2005) Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315. http://www.sciencedirect.com/science/article/pii/S0021999105001142
 Narayan & McClintock (2008) Narayan, R., & McClintock, J. E. 2008, New A Rev., 51, 733
 Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13
 Narayan & Yi (1995) —. 1995, ApJ, 452, 710
 Parfrey et al. (2015) Parfrey, K., Giannios, D., & Beloborodov, A. M. 2015, MNRAS, 446, L61
 Parker (1955) Parker, E. N. 1955, ApJ, 121, 491
 Parker (1957) —. 1957, J. Geophys. Res., 62, 509
 Parker (1966) —. 1966, ApJ, 145, 811
 Parker (1967) —. 1967, ApJ, 149, 535
 Petschek (1964) Petschek, H. E. 1964, NASA Special Publication, 50, 425
 Piano et al. (2012) Piano, G., Tavani, M., Vittorini, V., & et al.,. 2012, A&A, 545, A110
 Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137
 Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49

Romanova et al. (2002)
Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace,
R. V. E. 2002, ApJ, 578, 420
 Romanova et al. (2011) —. 2011, MNRAS, 416, 416
 Romero et al. (2003) Romero, G. E., Torres, D. F., Kaufman Bernadó, M. M., & Mirabel, I. F. 2003, A&A, 410, L1
 Ryan et al. (2017) Ryan, B. R., Gammie, C. F., Fromang, S., & Kestener, P. 2017, ApJ, 840, 6
 Salvesen et al. (2016a) Salvesen, G., Armitage, P. J., Simon, J. B., & Begelman, M. C. 2016a, MNRAS, 460, 3488
 Salvesen et al. (2016b) Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016b, MNRAS, 457, 857
 Servidio et al. (2011) Servidio, S., Dmitruk, P., Greco, A., et al. 2011, Nonlinear Processes in Geophysics, 18, 675
 Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
 Shi et al. (2016) Shi, J.M., Stone, J. M., & Huang, C. X. 2016, MNRAS, 456, 2273
 Simon et al. (2012) Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685
 Simon et al. (2011) Simon, J. B., Hawley, J. F., & Beckwith, K. 2011, ApJ, 730, 94
 Singh et al. (2015) Singh, C. B., de Gouveia Dal Pino, E. M., & Kadowaki, L. H. S. 2015, ApJ, 799, L20
 Singh et al. (2016) Singh, C. B., Mizuno, Y., & de Gouveia Dal Pino, E. M. 2016, ApJ, 824, 48
 Soker (2010) Soker, N. 2010, ApJ, 721, L189
 Stone & Gardiner (2010) Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142
 Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
 Stone et al. (2010) —. 2010, Athena: Gridbased code for astrophysical magnetohydrodynamics (MHD), Astrophysics Source Code Library, , , ascl:1010.014
 Stone et al. (1996) Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656
 Su et al. (2013) Su, Y., Veronig, A. M., Holman, G. D., et al. 2013, Nature Physics, 9, 489
 Sweet (1958) Sweet, P. A. 1958, in IAU Symposium, Vol. 6, Electromagnetic Phenomena in Cosmical Physics, ed. B. Lehnert, 123
 Takamoto et al. (2015) Takamoto, M., Inoue, T., & Lazarian, A. 2015, ApJ, 815, 16
 Uzdensky & Spitkovsky (2014) Uzdensky, D. A., & Spitkovsky, A. 2014, ApJ, 780, 3
 Walker et al. (2016) Walker, J., Lesur, G., & Boldyrev, S. 2016, MNRAS, 457, L39
 Wielgus et al. (2015) Wielgus, M., Fragile, P. C., Wang, Z., & Wilson, J. 2015, MNRAS, 447, 3593
 Yamada et al. (2010) Yamada, M., Kulsrud, R., & Ji, H. 2010, Reviews of Modern Physics, 82, 603
 Zanni & Ferreira (2009) Zanni, C., & Ferreira, J. 2009, A&A, 508, 1117
 Zanni & Ferreira (2013) —. 2013, A&A, 550, A99
 Zenitani et al. (2009) Zenitani, S., Hesse, M., & Klimas, A. 2009, ApJ, 696, 1385
 Zenitani & Hoshino (2008) Zenitani, S., & Hoshino, M. 2008, ApJ, 677, 530
 Zhdankin et al. (2013) Zhdankin, V., Uzdensky, D. A., Perez, J. C., & Boldyrev, S. 2013, The Astrophysical Journal, 771, 124. http://stacks.iop.org/0004637X/771/i=2/a=124
 Zhu et al. (2007) Zhu, Z., Hartmann, L., Calvet, N., et al. 2007, ApJ, 669, 483