AASTeX sample article

MHD instabilities in accretion disks and their implications in driving fast magnetic reconnection

[ Universidade de São Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Departamento de Astronomia
1226 Matão Street
São Paulo, 05508-090, Brasil
Elisabete M. de Gouveia Dal Pino Universidade de São Paulo, Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Departamento de Astronomia
1226 Matão Street
São Paulo, 05508-090, Brasil
James M. Stone Department of Astrophysical Sciences, Peyton Hall, Princeton University
Princeton, NJ 08544, USA
March 21, 2018
Abstract

Magnetohydrodynamic instabilities play an important role in accretion disks systems. Besides the well-known effects of the magnetorotational instability (MRI), the Parker-Rayleigh-Taylor 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 three-dimensional magnetohydrodynamical (3D-MHD) shearing-box 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 large-scale magnetic loops due to the PRTI followed by the development of a nearly steady-state 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 turbulence-induced fast reconnection.

accretion, accretion disks — magnetohydrodynamics (MHD) — instabilities — turbulence — magnetic reconnection
journal: ApJ\watermark

DRAFT

\correspondingauthor

Luis H.S. Kadowaki

0000-0002-6908-5634]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 low-mass 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 differentially-rotating disk. Since the early 1990’s, the magneto-rotational 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, low-density magnetized corona around them. In particular, the Parker-Rayleigh-Taylor 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 magnetic-pressure 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, low-density magnetized disk corona is invoked in order to explain non-thermal high-energy emission in these systems, such as the well known X-ray 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 sub-Eddington regime (see Esin et al., 1997, 1998, 2001; Narayan & McClintock, 2008), also known as ADAF (advection-dominated 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; Bosch-Ramon 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 first-order 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 non-thermal 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 high-energy emission of hundreds of low-luminosity (non-Blazar) 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 Sweet-Parker 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 X-point 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 Sweet-Parker 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 current-kink 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 non-relativistic 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 non-thermal 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 3D-MHD numerical simulations employing a shearing-box 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 shearing-box 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 shearing-box 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 shearing-box approach (in the disk reference frame at ), the right side of the momentum equation (4) is re-written 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 3D-MHD equations with the shearing-box 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 second-order Runge-Kutta 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 Courant-Friedrich-Lewy (CFL) condition, making the time step () less restrictive (see, Johnson et al., 2008; Davis et al., 2010; Stone & Gardiner, 2010).

We have used shearing-periodic 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, zero-gradients 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 well-known 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 condition111Different 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 3D-MHD shearing-box 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.

{longtable*}

lccccccc Simulation parameters.


[-2ex] Simulation & Computacional domain & Resolution & & & & & Vertical boundaries



[0.5ex]

[-1.8ex]

\endfirsthead

Table 0 - Continued

[0.5ex]


[-2ex]

Simulation & Computacional domain & Resolution & & & & & Vertical boundaries



[0.5ex]

[-1.8ex]

\endhead

Continue 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 two-dimensional 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 1: The top diagram shows the initial profile of the density distribution of the model R21b1, in a logarithmic scale (background colors). The scale color of the streamlines corresponds to the total intensity of the magnetic field. The bottom diagram shows the formation of magnetic loops in the higher latitudes of the system at due to the PRTI. The computational domain has a resolution of cells per thermal height scale of the system (), in a box of dimensions .

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 steady-state 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 shearing-boxes 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.

Figure 2: The upper diagram shows the time evolution of the magnetic energy density evaluated for the , and components. The solid line corresponds to the radial component “”, the dotted line to the azimuthal component “” and the dashed line to the vertical component “”. The bottom diagram shows the time evolution of the Shakura & Sunyaev (1973) viscosity parameter , where the solid line gives the contribution of the Maxwell stress tensor , and the dotted line the contribution of the Reynolds stress tensor . In both diagrams, volume averages have been obtained from the entire computational domain.

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).

Figure 3: Time evolution of the volume averages of the magnetic energy density, the Maxwell stress tensor and the ratio between the thermal and magnetic pressures () in different layers of the R21b1 model. The black, red and blue colors correspond to the averages obtained in the disk region (between and ), the upper corona (between and ), and the lower corona (between and ), respectively. The green line corresponds to the volume average obtained in the entire computational domain.

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 gas-pressure 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 4: Evolution of the magnetic field lines and density distribution in and for R21b1 model. The colors are the same as in the diagrams of Figure 1.

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 5: The upper diagram shows the time evolution of the total magnetic energy density (black line), the average componente (blue line) and the fluctuating component (red line). The bottom diagram shows the evolution of the ratio between the mean and the fluctuating components. After orbital periods, the fluctuations become dominant in the system.

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).

Figure 6: Time evolution and vertical profile of the horizontal averages of the , and components. It is possible to verify that after orbital periods, and present periodic variations with polarity inversions every orbital periods which are due to dynamo action triggered mainly by the MRI.

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 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 lines222Note that since the mean vertical is null, it should not effectively participate in the amplification of the large-scale 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 333Note 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).

Figure 7: The diagrams from top to bottom show the evolution of the total magnetic energy density, the Maxwell stress tensor, and , respectively, obtained in the disk and coronal regions. For all diagrams, the colors black, red and blue correspond to the volume averages of these quantities obtained for different initial values of . The black line also corresponds to the reference model (R21b1, see Figure 3).

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 volume-averaged value of in the disk (continuous line), upper (short-dashed line) and lower corona (long-dashed 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 gas-pressure 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 20-30 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 shearing-box 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 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). On the other hand, in stratified shearing-box 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 so-called 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 high-resolution 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 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 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 high-resolution simulations (R21b1 and R43b1 models, respectively), indicating that resolutions lower than could affect the initial evolution of the system.

Figure 8: Time evolution of the magnetic energy density (top diagram) and the Maxwell stress tensor (bottom diagram) for the resolutions of (black line), (red line), and (blue line).

After orbital periods, when the MRI dominates, the differences in the evolution between the low and intermediate resolution models is significant. The low-resolution 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 high-resolution models indicate a good agreement between them until the evolved time of the high-resolution 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 Sweet-Parker 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 Lazarian-Vishniac 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 sub-Alfvé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 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 () for each cell of the last subsample (see Zhdankin et al., 2013):

(17)

where 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 . 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).

Figure 9: The left sketch shows a new coordinate system centered in a local reconnection site. We have identified the orthonormal vector as the direction of the fastest decrease of . We have assumed this direction as the thickness of the magnetic reconnection site. The right sketch shows the details of the reconnection configuration, where the edges have been defined as the cells of .

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 .

Figure 10: Top: diagram of the coronal region above the disk (upper corona, ) at for R21b1 model. The colors in the background correspond to the current density magnitude. The white points correspond to the local maxima identified by the algorithm and the green points correspond to the magnetic reconnection sites with the most symmetric magnetic and velocity field profiles. Bottom: the diagrams show examples of rejected and accepted profiles of , and interpolated along the axis. Those evidencing substantial symmetry in velocity and magnetic field in the separatrix have been accepted.

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).

Figure 11: Distributions of the sign’s flip of the magnetic field components “” in the direction “” (). The counts have been organized in bins of orbital periods. The upper and lower histograms correspond to the distributions taken in the coronal region and disk, respectively.

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 induced-turbulence to drive fast reconnection in accretion disks and their coronae.

Figure 12: Time evolution of the magnetic reconnection rate () of the R21b1 model evaluated in the disk and corona separately. The continuous line corresponds to the mean value obtained from the histograms for each time step and the blue shade corresponds to the standard deviation.

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 steady-state. 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.

Figure 13: Histograms of the magnetic reconnection rates (left diagrams) and the thickness (right diagrams) of the reconnection sites in the coronal region (upper diagrams) and disk (lower diagrams). The distributions have been obtained between and orbital periods.

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 two-dimensional power spectrum of the velocity field 444The 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).555We 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 shearing-box simulations, have decomposed the energy power spectrum into a parallel (large-scale shear-aligned) and a perpendicular (small-scale 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 .

Figure 14: Compensated power spectrum of the velocity field “” taken at , and (from the top to the bottom). The colors indicate the height where the power spectrum have been evaluated. The dashed line corresponds to the Kolmogorov spectrum .

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 steady-state (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 high-resolution models are slightly smaller than for the low-resolution.

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 high-resolution 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 Time average of magnetic reconnection rate.

[-2ex] Simulation & & &

Name & Corona & Disk & Range

[0.5ex]

[-1.8ex]

\endfirsthead

Table 0 - Continued

[0.5ex]

[-2ex]

Simulation & & &

& Corona & Disk &

[0.5ex]

[-1.8ex]

\endhead

Continue in the next page…

\endfoot\endlastfoot

R11b1 & & & 20P-100P

R21b1 & & & 20P-100P

R21b10 & & & 20P-100P

R21b100 & & & 20P-100P

R43b1 & & & 20P-58P

Figure 15: The left diagrams show the time evolution of the magnetic reconnection rate for different initial values of (for the resolution). The black line corresponds to the reference model with (R21b1), whereas the red and blue lines correspond to the models with (R21b10) and (R21b100), respectively. The right diagrams show the time evolution of for different resolutions. In this case, the black line corresponds to the low resolution model (, R11b1), whereas the red and blue lines correspond to the intermediate (, R21b1) and high (, R43b1) resolution simulations, respectively. We have evaluated in the coronal region (upper diagram) and in the disk (lower diagram).

5 Discussion and conclusions

In this work, we have performed 3D-MHD shearing-box 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 Parker-Rayleigh-Taylor 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 gamma-ray 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 steady-state 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 pressure-dominated disk (with ) and a magnetically-dominated 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 steady-state regime at the same time when the MRI becomes self-sustained inside the system. The low-resolution model (R11b1), on the other hand, shows a continuous slow temporal decrease of the magnetic field after orbital periods without achieving a steady-state indicating that this model did not achieve appropriate resolution.

  • The values found for the Maxwell stress tensor in the steady-state 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 large-scale 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 large-scale 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 turbulence-induced 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 flux666Zhdankin 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 non-thermal 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 X-ray 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 non-ideal MHD simulations in shearing-box simulations).

Furthermore, we have dealt with isothermal shearing-box simulations, but an extension of the analysis to a non-isothermal 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.

Acknowledgments. The numerical simulations in this work have been performed in the Blue Gene/Q supercomputer supported by the Center for Research Computing (Rice University) and Superintendência de Tecnologia da Informação da Universidade de São Paulo (USP). This work has also made use of the computing facilities of the Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul), whose purchase was made possible by the Brazilian agency FAPESP (grant 2009/54006-4) and the INCT-A; and the cluster of the group of Plasmas and High-Energy Astrophysics (GAPAE), acquired also by the Brazilian agency FAPESP (grant 2013/10559-5). LHSK acknowledges support from the Brazilian agency FAPESP (postdoctoral grant 2016/12320-8) and CNPq (grant 142220/2013-2). EMGDP also acknowledges partial support from the Brazilian agencies FAPESP (grant 2013/10559-5) and CNPq (grant 306598/2009-4). Support from an international cooperation grant between Princeton University and the Universidade de São Paulo is gratefully acknowledged. Also, We would like to thank Kengo Tomida, Zhaohuan Zhu, Grzegorz Kowal and Ji-Ming Shi for useful discussions.








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
  • Bosch-Ramon et al. (2005) Bosch-Ramon, 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 e-prints, 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., & Otmianowska-Mazur, 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: Grid-based 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/0004-637X/771/i=2/a=124
  • Zhu et al. (2007) Zhu, Z., Hartmann, L., Calvet, N., et al. 2007, ApJ, 669, 483
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
244068
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description