1 Introduction

Driving Disk Winds and Heating up Hot Coronae by MRI Turbulence


We investigate the formation of hot coronae and vertical outflows in accretion disks by magneto-rotational turbulence. We perform local three-dimensional (3D) MHD simulations with the vertical stratification by explicitly solving an energy equation with various effective ratios of specific heats, . Initially imposed weak vertical magnetic fields are effectively amplified by magnetorotational instability (MRI) and winding due to the differential rotation. In the isothermal case (), the disk winds are mainly driven by the Poynting flux associated with the MHD turbulence and show quasi-periodic intermittency. On the other hand, in the non-isothermal cases with , the regions above 1-2 scale heights from the midplane are effectively heated up to form coronae with the temperature of 50 times of the initial value, which are connected to the cooler midplane region through the pressure-balanced transition regions. As a result, the disk winds are mainly driven by the gas pressure with exhibiting more time-steady nature, although the nondimensional time-averaged mass loss rates are similar to that of the isothermal case. Sound-like waves are confined in the cool midplane region in these cases, and the amplitude of the density fluctuations is larger than that of the isothermal case.

Subject headings:
accretion, accretion disks — MHD — stars: winds, outflows — planetary systems: protoplanetary disks — turbulence


1. Introduction

In accretion disks magnetohydrodynamical (MHD) turbulence is believed to work as effective viscosity and play an essential role in the outward transport of the angular momentum and the radial motion of the gas (e.g. Lynden-Bell & Pringle, 1974; Balbus & Hawley, 1998). Magnetorotational instability (MRI hereafter Velikhov, 1959; Chandrasekhar, 1960; Balbus & Hawley, 1991) is a promising source of such turbulent viscosity in accretion disks.

Properties of the MRI in accretion disks have been widely studied. MHD simulations in the local shearing coordinates (Goldreich & Lynden-Bell, 1965) have been extensively and intensively performed without vertical density stratification (e.g., Hawley et al., 1995; Matsumoto & Tajima, 1995; Sano et al., 2004), and with the vertical stratification due to the gravity by a central object (e.g., Stone et al., 1996; Miller & Stone, 2000; Davis et al., 2010; Shi et al., 2010; Okuzumi & Hirose, 2011; Sai et al., 2013).

Accretion disks threaded with global magnetic fields are supposed to drive disk winds, which was first suggested by Blandford & Payne (1982) as an origin of jets from black hole accretion disks. This picture of magnetocentrifugally driven disk winds was extended to protostellar jets (Pudritz & Norman, 1983). Outflows and jets from various types of accretion disks have been observed, e.g., from young stars (Ohashi et al., 1997; Coffey et al., 2004), and from active galactic nuclei (Tombesi et al., 2010).

These two pictures – MRI-driven turbulence in accretion disks and magnetically driven winds from disk surfaces involving coherent field lines – are supposed to have a close link. For instance, turbulence in disks could drive the vertical motions of the gas by the MHD turbulent pressure. Based on these considerations, Suzuki & Inutsuka (2009) and Suzuki et al. (2010) recently proposed the onset of vertical outflows from MRl-turbulent accretion disks. They performed MHD simulations in local shearing boxes with the vertical stratification and net vertical magnetic fields by taking a special care of the outgoing boundary conditions (Suzuki & Inutsuka, 2006) at the vertical surfaces of the simulation boxes. They found that the disk winds are driven from the upper and lower boundaries by the Poynting flux associated with the MHD turbulence. This process can contribute to the mass loading to the basal regions of the global disk winds introduced above. Such vertical outflows are also observed in 3D global simulations (Flock et al., 2011; Suzuki & Inutsuka, 2013).

Later on, such MRI-turbulent driven disk winds have been further studied from various aspects. Bai & Stone (2013) examined properties of the disk winds with stronger net vertical magnetic fields. The connection between the upflows in local simulations and the global disk winds/jets is investigated (Lesur et al., 2013; Bai & Stone, 2013). Fromang et al. (2013, see also Suzuki & Inutsuka 2010) pointed out the mass flux of the disk winds depends on the simulation box size, which shows that we need great cares to handle disk winds by local shearing boxes.

Although there are limitations in the shearing box approximation, it is still useful technic to study the basic properties of MRI-turbulent driven vertical outflows. Most of previous simulations adopt an isothermal equation of state to mainly focus on the dynamics of the disk winds, apart from a limited number of works that considered detailed heating and cooling processes for black hole accretion disks (Turner, 2004; Hirose et al., 2006) and protoplanetary disks (Hirose & Turner, 2011). So far there has been no systematic studies done for the MRI-turbulent driven winds with different ratios of specific heats, even within a framework of the shearing box approximation, which is the main focus of the present paper.

2. Simulation Setups

In Suzuki & Inutsuka (2009) & Suzuki et al. (2010) we performed 3D MHD simulations in local stratified shearing boxes (Stone et al., 1996) by solving the ideal MHD equations with an isothermal equation of state. In this paper, we extend these works by explicitly solving an energy equation in a Lagrangian form,


where is the Keplerian rotation frequency, is the specific energy per mass which is related to the gas pressure, , the density, , and the effective ratio of specific heats, , as


and the other variables in Equation (1) have the conventional meanings. The terms involving in Equation (1) originate from the central object; is the potential due to the vertical component of the gravity and denotes the tidal potential.

In the energy equation above we do not explicitly consider external cooling and heating processes, e.g., radiation cooling/heating, thermal conduction, and etc. Instead, we study their effect phenomenologically by assuming different but spatially uniform from 1 to in different cases. In our simulations, the gas is heated up mainly by the dissipation of magnetic energy, which we discuss later in this section. Taking (isothermal condition) indicates that we (implicitly) assume that the temperature is kept constant by an unspecified cooling that can balance the magnetic heating. On the other hand, larger corresponds to suppressing cooling; and correspond to the adiabatic conditions for diatomic and monoatomic gases, respectively.

Apart from explicitly solving the energy equation, we basically follow Suzuki & Inutsuka (2009) when performing the numerical simulations. We adopt a second-order Godunov-CMoCCT scheme (Sano et al., 1999), in which we solve the nonlinear Riemann problems with the magnetic pressure at the cell boundaries for the compressive waves and adopt the consistent method of characteristics (CMoC) for the evolution of magnetic fields (Clarke, 1996). At the top and bottom boundaries, we prescribe the outgoing boundary condition by using the seven MHD characteristics (Suzuki & Inutsuka, 2006).

We fix the ratio of the box sizes of the , , and axes to 1:4:8, which are respectively resolved by 32, 64, and 256 spatially uniform grids. Each cell is elongated in the direction with being twice as large as and . The lengths in the simulations are normalized by the initial pressure scale heights, , as


where is the initial sound speed, and we also use the initial temperature, , which has the dimension of . Hereafter subscripts ‘0’ represent the initial state. In this paper, stands for an “isothermal” sound speed, ; the usual sound speed, , is expressed as . We initially set up the hydrostatic density structure with the constant temperature, ,


where is the initial density at the midplane. In order to perform the simulations stably, we adopt a floor value, throughout the simulations. In the non-isothermal cases, the simulation boxes are larger than that of the isothermal case to treat the extended coronae as shown later. In the cases with , the initial densities near the surface regions, , are smaller than , hence is used for these regions. In the simulations, we use the unit of , , and ; accordingly, and . We initially impose the weak net vertical magnetic field, , with the plasma value, at the midplane. We start the simulations with giving small random perturbations with 0.5% of as seeds of the MRI. The simulations are kept running until 200 rotations ().

In the simulations with , the gas is heated by the dissipation of the magnetic energy, where the heating is done by the numerical effect in the sub-grid scale since we do not explicitly take into account resistivity. Physically, we expect that the heating is due to cascading turbulence at small scales. As a disk is heated up, the sound speed increases, and accordingly the pressure scale height also increases. Therefore, to study the vertical structure, the coordinate in units of the scale height,


is a key quantity, where the subscripts, , of indicate the average over a horizontal () plane. While for the isothermal situation, generally as a result of the heating for .

1.0 0.9,3.6,7.2 7.2 140–200 0.907
1.01 1.0,4.0,8.0 7.1 120–200 0.895
1.03 1.35,5.4,10.8 7.1 100–200 0.896
1.1 4.0,16.0,32.0 7.3 150–200 0.951
1.2 4.0,16.0,32.0 7.1 130–200 0.906
1.3 4.0,16.0,32.0 7.3 130–200 0.937
1.4 4.0,16.0,32.0 7.4 120–200 0.934
5/3 4.0,16.0,32.0 7.2 120–200 0.890

Note. – The initial box size (2nd column) and the final vertical box size measured in the final scale height (3rd column; Equation 6) averaged over are compared. ,, and in the 2nd column are the sizes of , , and components of the simulation box. The 4th column shows the final mass, , normalized by the initial mass, , remained in the box.

Table 1Simulation runs with different ratios of specific heats, .

The mass flux of the disk winds depends on the vertical box size in units of the scale height (Suzuki et al., 2010; Fromang et al., 2013). In order to compare the properties of the disk winds with different , it is desirable to adopt the same vertical box size in units of the final scale height after the steady-state conditions are achieved,


where and indicate the locations of the bottom and top boundaries of the simulation box, and we take the average over (Table 1) after the magnetic fields are amplified to the saturated state, which we describe later in this section. We cannot estimate in advance, hence we perform simulations with different box sizes and pick up one that gives the desirable value for each . For the isothermal case (), we fix the vertical box size , indicating and for the box sizes in the and directions. For the other cases with , we tune the initial box sizes to give the final vertical box sizes after the saturation in a range of .

Figure 1.— Evolution of the box-averaged of the three cases with (black solid), 1.03 (red dashed), and 1.0 (blue dotted). The vertical axis, , indicates time in units of one rotation.
Figure 2.— Time-distance diagrams of of the cases with (top), 1.03 (middle) and 1.0 (lower). The horizontal axis () denotes time in units of rotation. The vertical axis shows ; on the left shown is distance in units of the initial scale height, , and on the right shown is distance in units of the final scale height, .

In order to check the saturation of the magnetic fields, we monitor the box-averaged values (Shakura & Sunyaev, 1973), which is the sum of the Reynolds and Maxwell stresses normalized by the gas pressure,


where is the difference of toroidal velocity from the Keplerian shear flow, . In Figure 1 we show the results of the cases with (black solid line), 1.03 (red dashed line), and 1.0 (blue dotted line). We also display the time evolution of the vertical structure of horizontally averaged in diagrams in Figure 2 to see more details. In the case with , the saturation is observed after , and we use rotations (1 rotation) to estimate the vertical box size in the final scale height (Equation 6). Although in the case with the box averaged almost saturates after (Figure 1), we use a more conservative , because at the midplane is growing in (middle panel of Figure 2). In the isothermal case (), we use rotations, because becomes saturated later than in the other two cases. On the right of each panel of Figure 2, we show in units of the final scale height, . The interval between the ticks with varies with height in the non-isothermal cases because the temperature is higher near the surfaces as shown in 3.

The initial and final box sizes and ’s are summarized in Table 1. For the cases with , the vertical box size measured in the final scale height shrinks more than a factor of 4 because of the heating, which we examine in §3. The derived is also used to examine the time-averaged vertical structures of each case in §3.

The horizontal box sizes measured in the scale height of the non-isothermal cases vary with height because of the variation of the temperature. Because of the higher temperature near the surfaces (§3), the horizontal box sizes in units of the final scale height, and , are smaller in the surface regions. of the non-isothermal cases. Although these horizontal sizes are insufficient to quantitatively discuss the saturation of the magnetic fields (e.g., Guan et al., 2009), we anticipate that it is still meaningful to compare the vertical disk structures with different .

3. Results

Figure 3.— Snapshot of the simulation with at 180 rotation time (). White lines indicate magnetic field lines, colors on the back show normalized temperature, , (left legend), color contours denote iso-density surfaces (right legend), and arrows indicate velocity fields.

Continued from the previous section we examine the time evolutions of the box averaged in Figure 1. After the quasi-steady saturated states are achieved, the three cases show different behavior. The isothermal case () show large fluctuations of from 0.015 to 0.04. On the other hand, the case with exhibits much milder behavior with kept in a range of 0.01-0.015. The case with shows intermediate behavior; although the increase of the is faster with showing large fluctuations in the earlier time, , it settles down to a softer state later, between those for and 1.0. Interestingly enough, the similar trends are observed in the disk winds, which we discuss the details in §3.3.

Figure 3 exhibits a snapshot structure of the case with at 180 rotation time (). The turbulent magnetic field, mostly dominated by the toroidal component, is amplified by MRI and winding due to the differential rotation. The temperature contour on the back shows that the temperatures in increase up to more than 50 times of the initial value, while the temperature in the midplane does not increase so much. -shaped field lines, which are typical for Parker (magnetic buoyancy) instability (Parker, 1955; Nishikori et al., 2006) are seen in the surface regions, and the vertical outflows are observed from both the upper and lower surfaces, similarly to those seen in the isothermal simulations (Suzuki & Inutsuka, 2009; Suzuki et al., 2010).

3.1. Coronae

Figure 4.— Time-distance diagrams of temperatures, , of the cases with (upper) and 1.03 (lower). The labels for the horizontal and vertical axes are the same as in Figure 2.
Figure 5.— Time-distance diagrams of densities, , of the cases with (top), 1.03 (middle), and 1.0 (bottom). The labels for the horizontal and vertical axes are the same as in Figure 2.

Figures 4 and 5 present diagrams of the temperatures and the densities of the three cases with , & 1.0, whereas we do not show the temperature of the isothermal () case. As shown in Figure 4 as well as Figure 3, the gas in the surface regions of the non-isothermal cases is heated up. The heating is done by the dissipation of the magnetic energy which is amplified by the MRI and the winding due to the differential rotation. Since the simulations do not explicitly take into account resistivity, the dissipation of the magnetic field occurs numerically in the sub-grid scales, whereas we expect that this actually takes place by cascading of small-scale turbulence.

The two cases of Figure 4 show an initial temperature rise in , because the initial densities in the upper regions, , are larger than the hydrostatic value (Equation 4); the gas initially flows down from there toward the midplane, which causes the initial heating. The entire region is eventually relaxing to the quasi-steady state, which we analyze from now. In both the cases, the regions of are most effectively heated up to in the case with and in the case with . Because the heating is done by the dissipation of the magnetic field which is amplified by the MRI and the winding, the increase of the temperature follows the increase of shown in Figure 2. This is more easily seen in the case with ; the temperature increases in , which is delayed compared to the increase of .

The top panel of Figure 4 illustrates that the hot regions with are clearly separated from the cool midplane region with by the transition regions located at ; from now we call these hot regions coronae. On the other hand, in the small (1.03) case, the temperatures of the surface regions are not so high with in the upper regions.

Figure 6.— Comparison of the time and horizontally averaged vertical structures of the cases with (solid), 1.1 (dot-dashed), 1.03 (dashed), and 1.0 (dotted). The horizontal axis is vertical distance in units of the initial scale height for the left panels and in units of the final scale heights, , for the right panels. From top to bottom, the temperatures normalized by the initial value, the densities, and the vertical velocities normalized by the initial constant sound speed, , (left panels) or normalized by the final local sound speeds, , (right panels) are compared. The gray thick lines in the middle panels indicate the initial density profile.

For more quantitative inspection we display the time and horizontally averaged vertical structures of the four cases, , & 1.4 in Figure 6, where the time averages are taken over in Table 1 after the quasi-steady saturation is achieved (Figure 1). The top panel illustrates that the coronae with also form in the case with . We have found that the formation of the coronae with is universal in the cases with . In these cases, with the large increases of the temperature (around for and for ) the density (middle panels) rapidly decreases to keep the nearly pressure-balanced structure. We call this pressure balanced structure with temperature rise a transition region.

The existence of the hot coronae above the cool midplane region is quite different from some results of recent stratified shearing box simulations. For instance, Bodo et al. (2012) performed the simulations with fixing the temperatures to the initial value () and the velocities to zero () at the top and bottom boundaries. They showed that the maximum temperature is obtained at the midplane and the temperature monotonically decreases to the surfaces. This shows that the boundary condition at the boundaries play a significant role in the vertical temperature and velocity structures.

Coupled with the formation of the hot coronae () or the warm regions (), the gas is lifted up from the midplane region to the upper regions as shown in Figure 5. In the case with (top panel), it is clearly seen that the coronae are gradually filled with denser material in by the evaporation from the midplane region, and at later times the evaporated gas is mainly supported by the gas pressure of the hot coronae (§3.3). The time averaged density structures of the cases and 1.4 (middle panels of Figure 6) quantitatively show the large uplifted gas compared with the initial hydrostatic profile with the constant (gray thick lines). The supply of the gas to the coronal regions is also seen in the case with (middle panel of Figure 5) although the amount is smaller than those of larger cases. Even in the isothermal case (bottom panel of Figure 5) quasi-periodic uplifting motions of the gas are observed, which are by the Poynting flux associated with the breakups of channel flows (Suzuki & Inutsuka, 2009). The uplifted gas is finally connected to disk winds as will be discussed in §3.3.

3.2. Saturation of Magnetic Field

Figure 7.— Comparison of the time and horizontally averaged vertical structures of the (top), (middle), and (bottom) components of magnetic energies, , normalized by the initial gas pressure at the midplane, (left), and quality factors, , for MRI (right) of the cases with (solid), 1.1 (dot-dashed), 1.03 (dashed), and 1.0 (dotted). The horizontal axis, , is measured in the final scale heights.
Figure 8.— Comparison of the time and horizontally averaged vertical structures of the plasma values (upper) and the values (bottom) of the cases with (solid), 1.1 (dot-dashed), 1.03 (dashed), and 1.0 (dotted). The horizontal axis, , is measured in the final scale heights.

Magnetic field plays a central role in the formation of the coronae because its dissipation leads to the heating. The left panels of Figure 7 compare the time and horizontally averaged magnetic energies of the four cases with , 1.1, 1.03, and 1.0. The toroidal () component dominates as expected because of the winding by the differential rotation. In the surface regions, , these four cases show similar profiles one another. On the other hand, in the midplane region larger cases give lower magnetic energies of the and components in particular. only weakly depends on because the strength of is mostly controlled by the winding. As a result, in the midplane region of larger cases more or less coherent toroidal magnetic field dominates (Figure 3).

The lower levels of and are mainly because of the insufficient resolution. The right panels of Figure 7 display the time and horizontally averaged -th components () of quality factors with respect to MRI (Noble et al., 2010; Hawley et al., 2011),


where is Alfvén velocity along with an –th direction and is the grid size of an –th component. measures the number of mesh points which resolves the most unstable wavelength of MRI. According to Sano et al. (2004), is a necessary condition for a vertical magnetic field to get a linear growth rate close to the analytic prediction from MRI. and of the cases with and 1.1 at the midplane is small . In these cases the initial scale height, , is resolved by 8 mesh points. Although at later time one scale height, , is eventually resolved by a larger number of grids, , it is still insufficient to capture fine turbulent structure by the MRI at the midplane (Simon et al., 2009; Hawley et al., 2011, 2013; Parkin & Bicknell, 2013a, b). Therefore, the magnetic field strength at the midplane in the large is supposed to be underestimated and not to develop to the physically saturated state. On the other hand, in the small cases and the surface regions of the large cases we can safely discuss the saturation of the magnetic fields.

Figure 8 displays the time and horizontally averaged plasma values,

(upper panel) and values,

(lower panel). is larger and is smaller for smaller . In the surface regions, this is explained by the larger gas pressure due to the hotter coronae for larger since the magnetic field strengths are similar among different cases (left panels of Figure 7). The large and small values at the midplane region in the cases with and 1.1 are mainly because of the weaker magnetic field strength owing to the insufficient resolution (Figure 7).

In the surface regions, , the values of the different cases are almost spatially constant (Figure 8), which indicates that the magnetic energy decreases with increasing (left panels of Figure 7) in the same manner as the decrease of the gas pressure, whereas the level of depends on reflecting the coronal gas pressure in the numerator of . In other words, the magnetic energy dissipates in these surface regions; the magnetic energy becomes relatively important with an elevating altitude from the midplane with decreasing density, and at the locations, , with times (Figure 6) and , i.e., nearly equipartition, the magnetic energy starts to be gradually converted to other forms. This is the source of the coronal heating and the driving of the disk winds. Although in our simulations the dissipation of the magnetic energy is due to numerical dissipation at the sub-grid scales, in reality the MHD turbulence is supposed to dissipate through the energy cascading to smaller scales.

3.3. Vertical Outflows

Figure 9.— Comparison of the time and horizontally averaged vertical structures of the gas pressure (dashed), the magnetic pressure (dotted), and the magnetic tension (solid) terms of the cases with (left), 1.03 (middle), and 1.0 (right). See text (Equations 9 & 10 and their explanations) for details.

The bottom panels of Figure 6 show that vertical outflows are accelerated in the surface regions. The left bottom panel () shows that in larger cases faster outflows are driven at higher altitudes. This is because the material is more extended to higher altitudes and the winds are effectively accelerated by the gas pressure of the hotter gas. If we plot the same quantity in the diagram (right bottom panel), the profiles of the shown four models resemble one another; the vertical velocities reach the local sound speeds at where the densities decrease to times . By these outflows, the total mass in the simulation box of each case gradually decreases with time. Typically 5-10 % of the initial mass is lost until the end of the simulations at rotations (Table 1).

We inspect how the driving mechanism of these disk winds is different for different . By rearranging Equation (1), we can derive an Eulerian form of the energy equation:


where and . One can see from the last equality, the energy flux can be divided into the two parts; the first bracket indicates the energy advected with and the second bracket indicates the work acting on gas; the second part exactly appears in the spatial derivative term in the Lagrangian form (Equation 1). Here, we examine this work component in the simulations. The component of the second part can be explicitly written as


where and . The first, second, and third terms indicate the work by gas pressure, magnetic pressure, and magnetic tension, respectively. The last term is canceled out by a term from the advection component. We compare the first – third terms of the three cases with (left panel), 1.03 (middle panel), and 1.0 (right panel) in Figure 9.

If a line in Figure 9 decreases with increasing in the region, the force acts on gas to drive a vertical upflow, and vice versa in the region. In the isothermal () case, the magnetic pressure (dotted line) and tension (solid line) comparably contribute to driving the vertical outflows, while the contribution from the gas pressure is quite small. The “injection regions” of the magnetic tension form around as pointed out by Suzuki & Inutsuka (2009); from these regions, the Poynting fluxes associated with the tension are injecting toward both the midplane and surface directions.

On the other hand, in the case with , the gas pressure largely dominates the magnetic components in driving the disk winds. In this case, the magnetic energy in low-altitude regions once dissipates to heat up the gas. The gas pressure, which increases owing to the magnetic heating, finally contributes to driving the vertical outflows. This is in contrast to the isothermal case, in which the magnetic forces directly drive the vertical outflows. The behavior of the case with lies between the two cases. While the largest contribution is from the magnetic pressure, the magnetic tension and the gas pressure also play a significant role.

Inspecting all the cases with different , we can conclude the following results on the driving mechanisms of the disk winds: While in the isothermal and small cases the vertical outflows are mainly driven by the Poynting flux, in the large cases the gas pressure dominantly drives the vertical outflows. However, we should note that, even in the large cases the magnetic fields play an important role, because the gas pressure is maintained by the dissipation of the magnetic energy which is amplified by the MRI and the winding due to the differential rotation. Since the disk winds are driven from the surface regions where the MRI is well-captured, increasing the numerical resolution does not change the obtained properties of the disk winds so much, although the magnetic field at the midplane of the large cases would be affected because the numerical resolution there is insufficient at the moment (Figure 7).

Magnetic buoyancy (Parker, 1955) is also considered to operate in the surface regions and contribute to the outflows (Nishikori et al., 2006; Machida et al., 2013). As discussed in Figure 3, –shaped field lines, typical for the Parker instability, are sometimes observed. Compared to the isothermal (see also Suzuki et al., 2010) and small cases, however, the contribution from the magnetic buoyancy is smaller in the large cases. In general, the Parker instability sets in for magnetically dominated (small ) condition with strong stratification (short pressure scale height). In the large cases the time and horizontally averaged is not so small (Figure 8) and the scale height in the surface regions is not short, which tend to suppress the Parker instability compared to the small cases.

Figure 10.— Time-distance diagrams of the normalized mass flux, of the cases with (upper) and 1.0 (lower). The labels for the horizontal and vertical axes are the same as in Figure 2.

The difference of the driving mechanisms influences the time-dependency of the disk winds. Figure 10 compares the diagrams of the mass flux of the cases with (upper panel) and 1.0 (lower panel). The two cases show quite different appearances. The isothermal case (lower panel) shows a clearer on-off nature of the disk winds. As reported in Suzuki & Inutsuka (2009), it is related to the quasi-periodic breakups of channel-mode flows with 5-10 rotation times. In the case with , the mass flux increases to the saturated value in , after the gas is supplied to the coronal regions (top panel of Figure 5). After that exhibits smoother structure in the coronal regions, , with rather quasi-steady vertical outflows. The midplane region () with (Figure 6) seems to be separated from the coronal regions by the transition regions at . The quasi-steady nature of the disk winds in this case is related to the fact that the vertical outflows are mainly driven by the gas pressure. The coronal temperatures are more or less uniformly distributed in , and the force by the gas pressure gradient are more time-steady, which is in contrast to the strong intermittency of the Poynting flux-driven outflows in the isothermal case.

3.4. Dependence on

Figure 11.— Dependences of time averaged three quantities on . Shown are the sum of the nondimensionalized mass flux of the vertical outflows, , from the upper and lower surfaces (eq.11; top), the peak temperature, the sum of the box-averaged Maxwell and Reynolds stresses, , (middle), and (bottom). In the shaded region at the left most location the results of the isothermal case () are plotted.

We summarize typical time-averaged quantities of the simulations with different . The top panel of Figure 11 shows the dependence of the mass flux of the disk winds. The shown quantity is


where the subscripts, ‘top’ and ‘bot’, indicates the top and bottom boundaries of the simulation box, and the subscript ‘mid’ stands for the midplane. The variables in the brackets on the right hand side are also time and horizontally averaged. The meaning of Equation (11) would be clear, the sum of the mass fluxes from the upper and lower simulation boundaries, which is further normalized by the time and horizontally averaged at the midplane. Note that the sign for the mass flux from the bottom surface is to pick up the outflowing direction.

The data points show that the nondimensional mass flux, , seems to be almost independent from and is distributed within a factor of 2, which is consistent with the trend obtained from the time-averaged vertical structure (bottom right panel of Figure 6)4. This indicates that the time-averaged dose not depend on the properties of the vertical outflows, either Poynting flux-driven winds with strong intermittency (smaller ) or more time-steady gas pressure-driven wind (larger ), which is a little surprising. We suppose that the main reason of the insensitive to is that the original source of the vertical outflows is the magnetic energy. As discussed in §3.2 the magnetic energy gradually dissipates to be converted to other forms of energies at the locations, , which is almost independent from , with times (Figure 6) and (nearly equipartition; Figure 8) .

In the small regime the magnetic pressure and tension directly accelerate the Poynting flux-driven vertical outflows, while in the large regime the magnetic energy is firstly transferred to the internal energy, the coronal heating in other words, and then the disk winds are driven by the gas pressure of the hot coronae. The locations (at ) of the energy conversion regulate the final mass flux, , which is insensitive to .

The relative comparison of among different cases of the simulations is meaningful since the vertical box sizes are tuned to give 7.1–7.4 in units of the final scale heights (Table 1). However, we should cautiously note that the absolute values of should be taken with cares because they depend on the vertical box sizes (Suzuki et al., 2010; Fromang et al., 2013); a larger vertical box would give a smaller .

The middle panel of Figure 11 presents the box- and time-averaged values. We do not find any monotonic trend of with and the values are typically . However, we cannot proceed detailed saturation arguments (e.g. Simon et al., 2009; Hawley et al., 2011; Parkin & Bicknell, 2013a), since the numerical resolution is not sufficient particularly in the midplane region of the large cases (right panels of Figure 7).

The bottom panel of Figure 11 compares the maximum temperatures, , normalized by the initial value, , of the different cases. is derived from the time and horizontally averaged vertical structure. The derived is located in a region of in the non-isothermal cases, whereas temperature becomes locally and transiently higher than . One can see a clear increasing trend with , because larger simply corresponds to smaller net cooling (cooling - heating). Moreover, jumps up from to , which corresponds to the change of the regime from the Poynting flux-driven winds to the gas pressure-driven winds.

3.5. Wave Phenomena

Figure 12.— Comparison of the time and horizontally averaged vertical structures of (dotted) and (solid) for the cases with (left), 1.03 (middle), and 1.0 (right). They correspond to the net energy fluxes associated with Alfvénic and sound-like waves, respectively. See text (Equations 12 & 13 and their explanations) for details.

In Suzuki & Inutsuka (2009) we examined the vertical energy flux associated with wave-like activities and found Alfvénic and sound-like waves propagating to both upward and downward directions. Following the procedure in Suzuki & Inutsuka (2009), we inspect the vertical structures of two quantities, and . is the Poynting flux of magnetic tension, as discussed in Equation (9), and can be rewritten as


where are Elsässer variables, which correspond to the amplitudes of Alfvén waves propagating to the -directions. Thus, corresponds to the net Poynting flux associated with propagating Alfvénic disturbances to the direction. is also rewritten as


where denote the amplitudes of sound waves5 propagating to the -directions. Here, note that the sound speed is expressed as since in this paper we define as isothermal sound speed.

Figure 12 compares these quantities of the three cases with (left panel), 1.03 (middle panel), and 1.0 (right panel). As already discussed in Suzuki & Inutsuka (2009), in the isothermal case sound-like waves propagating to the midplane are observed. The peak values are located at the injection regions, for the Alfvénic waves (dotted line); Alfvénic waves are injected from these regions mostly associated with the breakups of channel flows.

The case with exhibits similar structures except for (solid line) in the midplane region, ; the direction of the sound-like waves is upward (to both ) from the midplane, which is opposite to that in the isothermal case. A speculative explanation is that the downward magnetic tension forces from the injection regions, which drive downward sound-like waves, are not relatively sufficient because the upward forces by the gas pressure are comparably important in this case (middle panel of Figure 9).

The case with shows very different behavior. The energy flux of Alfvénic waves is mostly dominated by that of sound-like waves. In particular, it is nearly zero in the midplane region between the transition regions at , which separate the cool midplane from the above hot coronae. The midplane region is protected from the magnetic perturbations in the upper coronae because the Alfvénic disturbances are reflected at the transition regions (see §3.6). The direction of the sound-like waves in this region is to the midplane. The large jumps are also seen in the sound-like waves at the transition regions; the sound speed also changes abruptly there owing to the change of the temperature, which causes the reflection of the sound-like waves. Thus, the sound-like waves are confined in the cool midplane region.

Figure 13.— Comparison of the time and horizontally averaged nondimensional density fluctuations, , of cases with (solid), 1.03 (dashed), and 1.0 (dotted).

Sound waves are associated with density perturbations. Therefore, we expect that the different structures of will give different density perturbations. In Figure 13 we compares the vertical structures of the nondimensional density perturbations, which are calculated as


The cases with (dashed line) and 1.0 (dotted line) exhibit similar vertical structures; larger in the surface regions decrease to at the midplane. On the other hand, the case with shows complicated structure. The jumps at coincide with the transition regions between the upper hot coronae and the lower cool midplane. In the midplane region, is larger than those obtained in the other two cases, probably because of the confinement of the sound-like waves in the midplane region (Figure 12). However in the coronal regions is much smaller than those in the other cases.

3.6. Time Evolution of

Figure 14.— Same as Figure 10 but for the toroidal magnetic fields.

The evolution of toroidal magnetic fields, , is widely discussed in terms of “dynamo” activity in disks by MRI (Brandenburg et al., 1995; Nishikori et al., 2006; Shi et al., 2010). In Figure 14, we compare the diagrams of the cases with (upper panel) and 1.0 (lower panel). The isothermal case (lower panel) exhibits usual quasi-periodic changing of the sign of , as seen by previous works (e.g. Davis et al., 2010). The case with also shows quasi-periodic oscillations of the , but the period is shorter than that in the isothermal case and the amplitude (contrast between red and blue regions) is smaller. The tendency is similar for other large cases. This indicates that the toroidal magnetic fields change the sign before the magnetic fields are amplified to the level as strong as that obtained in the isothermal case. This is probably because in the large case the magnetic fields are more subject to the gas motion because of the high condition (Figure 8).

The difference of the time-dependencies is consistent with the time evolutions of as shown in Figure 1. For larger , shows smaller fluctuations with time, because the dynamics are largely controlled by the gas pressures, which distribute more uniformly in a more time-steady manner, rather than by the magnetic fields, which intermittently pile up associated with channel flows.

Another interesting feature of the non-isothermal case with is that the butterfly pattern triggered in the coronal regions is very vague in the midplane region. The left panel of Figure 9 also shows that the Poynting flux is negligibly small in this region. They imply that the miplane region seems to be protected from the magnetic perturbations in the coronal regions. The transition region which separates the cool midplane and the hot corona accompanies the large density difference to satisfy the pressure balance structure (Figure 6). This inevitably leads to the large jump of the Alfvén speed, , across the transition region. Then, Alfvénic perturbations arising from magnetic tension suffer reflection and the downward Poynting flux from the coronal regions are mostly reflected back and cannot penetrate the transition regions to the midplane. Reflection of Alfvén waves is widely discussed on the sun because it is very efficient at the transition region between the cool chromosphere and the hot corona (e.g., Suzuki & Inutsuka, 2005; Matsumoto & Suzuki, 2012). A signature of reflected Alfvén waves is actually observed on the solar surface by the HINODE satellite (Fujimura & Tsuneta, 2009). These are quite similar to what we observe in the present simulations with relatively large .

4. Summary and Discussion

Using the 3D MHD simulations in the local stratified shearing boxes with weak net vertical magnetic fields, we have studied the formation of the hot coronae and the disk winds induced by the MRI turbulence in the accretion disks. Taking into account the effect of the cooling phenomenologically with the effective ratio of the specific heats, , we have inspected how the basic properties are affected by the cooling. Since the simulations are performed in the nondimensional form without any physical scale, the simulations are applicable to various objects. The amplifications of the magnetic fields are observed in all the simulation runs with different from 1 to 5/3, which give the time- and box-averaged , whereas we should take these values with cares because the numerical resolution is not sufficient in the midplane region.

The properties of the coronae and the disk winds are not so affected by the numerical resolution because the simulations well capture the MRI there, and we have found that the results are classified into the two regimes. In the small regime, the temperatures in the surface regions are not high because the effect of the heating is weak owing to the small . The vertical outflows are directly driven by the Poynting flux associated with the amplified turbulent magnetic fields, they are more time-dependent, involved with the intermittent breakups of large-scale channel flows.

In the large regime, the hot coronae form by the dissipation of the magnetic energy and the temperatures are times of the initial values with the time-averaged peak temperature slowly increasing with . The vertical outflows are mainly driven by the gas pressure of the hot coronae. Because the spatial distribution of the gas pressure is more uniform than that of the magnetic energy, the disk winds stream out in a more time-steady manner than in the small regime. The hot coronae are connected to the cool midplane through the sharp transition regions. Across the transition regions, both the sound and Alfvén speeds change abruptly because of the jumps of the temperature and the density. The transition regions work as the walls against wave activity. Sound-like waves are confined in the cool midplane region with giving the larger amplitudes of the density perturbations. The midplane region is also protected from the magnetic perturbations in the upper coronae.

Although the driving mechanisms and the time-dependencies of the vertical outflows are different for the small and large regimes, the time-averaged nondimensional mass fluxes, , are similar each other. This is because in both the regimes the origin that drives the vertical outflows is the magnetic energy that is amplified by the MRI. The magnetic energy is gradually converted at the locations with , directly to the kinetic energy of the disk winds in the small regime, or firstly to the thermal energy of the hot coronae that is finally transferred to the disk winds in the large regime. The location of the energy conversion, which is insensitive to , controls the final . We should cautiously note that the derived mass flux, , depends on the vertical box size (Suzuki et al., 2010; Fromang et al., 2013), although the comparisons of different cases with the same vertical box size make sense.

Our treatment of spatially uniform is too much simplified. In reality, should be non-uniform. Deeper regions near the midplane tend to be more optically thick, which gives larger , while in surface regions is smaller. In more elaborated models with applications to specific objects, adequate cooling and heating processes should be included with radiative transfer (e.g., Hirose et al., 2006).

The authors thank Prof. Shu-ichiro Inutsuka for many fruitful discussion. The authors also appreciate constructive comments by an anonymous referee. This work was supported in part by Grants-in-Aid for Scientific Research from the MEXT of Japan, 22864006. Numerical simulations in this work were carried out with SR16000 at the Yukawa Institute Computer Facility and with the Cray XT4 and XC30 operated in CfCA, National Astrophysical Observatory of Japan.


  1. affiliation: Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan
  2. affiliation: Nihon Unisys, Ltd., Toyosu, Koto-ku, Tokyo, 135-8560, Japan
  3. affiliation: Department of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan
  4. Furthermore, since the normalization, , only weakly depends on (see top and middle panels of Figure 6), the time-averaged mass flux itself is almost independent from
  5. Strictly speaking, these are magnetosonic waves, namely the fast mode in the high plasma, and the slow mode that propagates along in the low plasma. Note also that the signs are opposite for and , reflecting the transverse and longitudinal characters.


  1. Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30
  2. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
  3. —. 1998, Reviews of Modern Physics, 70, 1
  4. Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  5. Bodo, G., Cattaneo, F., Mignone, A., & Rossi, P. 2012, ApJ, 761, 116
  6. Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741
  7. Chandrasekhar, S. 1960, Proceedings of the National Academy of Science, 46, 253
  8. Clarke, D. A. 1996, ApJ, 457, 291
  9. Coffey, D., Bacciotti, F., Woitas, J., Ray, T. P., & Eislöffel, J. 2004, ApJ, 604, 758
  10. Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52
  11. Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122
  12. Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71
  13. Fujimura, D., & Tsuneta, S. 2009, ApJ, 702, 1443
  14. Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125
  15. Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010
  16. Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  17. Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84
  18. Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102
  19. Hirose, S., Krolik, J. H., & Stone, J. M. 2006, ApJ, 640, 901
  20. Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30
  21. Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61
  22. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
  23. Machida, M., Nakamura, K. E., Kudoh, T., et al. 2013, ApJ, 764, 81
  24. Matsumoto, R., & Tajima, T. 1995, ApJ, 445, 767
  25. Matsumoto, T., & Suzuki, T. K. 2012, ApJ, 749, 8
  26. Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398
  27. Nishikori, H., Machida, M., & Matsumoto, R. 2006, ApJ, 641, 862
  28. Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959
  29. Ohashi, N., Hayashi, M., Ho, P. T. P., & Momose, M. 1997, ApJ, 475, 211
  30. Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65
  31. Parker, E. N. 1955, ApJ, 122, 293
  32. Parkin, E. R., & Bicknell, G. V. 2013a, ApJ, 763, 99
  33. —. 2013b, MNRAS
  34. Pudritz, R. E., & Norman, C. A. 1983, ApJ, 274, 677
  35. Sai, K., Katoh, Y., Terada, N., & Ono, T. 2013, ApJ, 767, 165
  36. Sano, T., Inutsuka, S., & Miyama, S. M. 1999, in Astrophysics and Space Science Library, Vol. 240, Numerical Astrophysics, ed. S. M. Miyama, K. Tomisaka, & T. Hanawa, 383
  37. Sano, T., Inutsuka, S.-i., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321
  38. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  39. Shi, J., Krolik, J. H., & Hirose, S. 2010, ApJ, 708, 1716
  40. Simon, J. B., Hawley, J. F., & Beckwith, K. 2009, ApJ, 690, 974
  41. Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656
  42. Suzuki, T. K., & Inutsuka, S.-i. 2005, ApJ, 632, L49
  43. —. 2006, Journal of Geophysical Research (Space Physics), 111, 6101
  44. —. 2009, ApJ, 691, L49
  45. —. 2013, submitted to ApJ
  46. Suzuki, T. K., Muto, T., & Inutsuka, S.-i. 2010, ApJ, 718, 1289
  47. Tombesi, F., Sambruna, R. M., Reeves, J. N., et al. 2010, ApJ, 719, 700
  48. Turner, N. J. 2004, ApJ, 605, L45
  49. Velikhov, E. P. 1959, Zh. Eksp. Teor. Fiz., 36, 1398
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
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

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 description