HALL EFFECT CONTROLLED GAS DYNAMICS IN PROTOPLANETARY DISKS. II. FULL 3D SIMULATIONS TOWARD THE OUTER DISK

The Hall effect is characterized by a physical length scale l_H. In the absence of charged grains, it reads as (Kunz & Lesur 2013)

$$\begin{equation} l_H\equiv \frac{v_A}{\omega _H}=\left (\frac{\rho }{\rho _i}\right) \left (\frac{v_A}{\omega _i}\right), \end{equation} \tag{ 1 }$$

where $v_A=B/\sqrt{4\pi \rho }$ is the Alfvén velocity, ω_i is the ion cyclotron frequency, ω_H = (ρ_i/ρ)ω_i is the Hall frequency, and ρ_i and ρ are the mass densities of the ions and the bulk of the gas, respectively, with ρ_i ≪ ρ for weakly ionized gas. Note that both v_A and ω_i are proportional to the magnetic field strength, so l_H is field-strength independent and is determined solely by the ionization fraction. In disks, it is natural to normalize l_H by the disk scale height H ≡ c_s/Ω. The associated Hall diffusivity η_H can be expressed as

$$\begin{equation} \eta _H=v_Al_H . \end{equation} \tag{ 2 }$$

At a fixed ionization fraction (∝ρ_i/ρ), we have η_H∝B/ρ.

AD is characterized by the frequency with which neutrals collide with ions γ_iρ_i, where γ_i is the coefficient of momentum transfer for ion-neutral collisions. In the disk, it is natural to normalize γ_iρ_i to the disk orbital frequency, by defining

$$\begin{equation} Am\equiv \frac{\gamma _i\rho _i}{\Omega }, \end{equation} \tag{ 3 }$$

which is the Elsasser number for AD. Generally, AD plays an important role in the gas dynamics when Am ≲ 10 (Bai & Stone 2011). The associated ambipolar diffusivity is given by

$$\begin{equation} \eta _A=v_A^2/\gamma _i\rho _i . \end{equation} \tag{ 4 }$$

At a fixed ionization fraction, we have η_A∝(B/ρ)².

The above definitions apply when electrons and ions are the main charged species, where the physics can be described most easily. Generalizations to include charged grains can be found in, e.g., Wardle (2007) and Bai (2011a), which are used in our vertically stratified simulations in subsequent sections.

Both η_H and η_A are inversely proportional to the ionization fraction, while different scalings of η_H and η_A with gas density indicate that AD becomes progressively more important toward low-density regions (e.g., the outer disk). To quantify this, we take the product of the two dimensionless numbers l_H/H and Am, which is independent of the ionization fraction:

$$\begin{equation} Am\cdot \left (\frac{l_H}{H}\right) =\frac{\gamma _i\rho }{\omega _i}\left (\frac{v_A}{c_s}\right) \propto \frac{\sqrt{\rho }}{c_s}. \end{equation} \tag{ 5 }$$

When adopting the minimum-mass solar nebula disk model (MMSN; Weidenschilling 1977; Hayashi 1981), we have at disk midplane ρ₀∝R^−11/4, c_s∝R^−1/4 and hence Am · (l_H/H)∝R^−9/8. More specifically, we find²

$$\begin{equation} Am\cdot \left (\frac{l_H}{H}\right)\approx 0.64\ \left (\frac{R}{10\,{\rm AU}}\right)^{-9/8}. \end{equation} \tag{ 6 }$$

In the outer region of PPDs, the value of Am is found to be of order unity for a wide range of disk radii (Bai 2011a, 2011b), and this formula provides a very useful relation in estimating the importance of AD and the Hall effect in PPDs. If we consider the Hall effect to be important when l_H/H ≳ 0.1, then the influence of the Hall effect extends to ∼50–60 AU.

For the MRI, the importance of the Hall effect and AD is characterized by their respective Elsasser numbers, defined as $v_A^2/\eta \Omega$, with η being η_H or η_A. With the AD Elsasser number introduced in Equation (3), the Hall Elsasser number can be written as (also see in Paper I)

$$\begin{equation} \chi \equiv \frac{\omega _H}{\Omega }. \end{equation} \tag{ 7 }$$

Note that χ depends on field strength (∝B) and also

$$\begin{equation} \frac{l_H}{H}=\frac{1}{\chi }\frac{v_A}{c_s}=\frac{1}{\chi }\sqrt{\frac{2}{\beta }} =\frac{X}{\sqrt{2\beta }}, \end{equation} \tag{ 8 }$$

where the plasma β = 8πP/B² is the ratio of gas to magnetic pressure, and X ≡ 2/χ is another commonly adopted dimensionless quantity in the literature (Sano & Stone 2002a, 2002b). The nonideal MHD term with an Elsasser number of order unity or less indicates that the term is physically significant.

A comparison between χ and Am reveals the relative importance between the Hall effect and AD, where a larger value indicates less importance. Using Equation (6), we find

$$\begin{equation} \frac{Am}{\chi }\approx 4.5\sqrt{\frac{\beta }{100}}\left (\frac{R}{10\, {\rm AU}}\right)^{-9/8}. \end{equation} \tag{ 9 }$$

Again, we see that Am and χ are likely of the same order for a wide range of disk radii given the expected magnetic field strength of β ≲ 100 (saturated β) in the outer disk.

In our definition, ω_H, l_H, and χ are all positive. On the other hand, the Hall term also depends on the polarity of the magnetic field relative to ${\boldsymbol{\Omega }}$. To distinguish the two cases, we always state explicitly in this paper the polarity of the background magnetic field B_z0 > 0 or B_z0 < 0 for fields aligned and anti-aligned with ${\boldsymbol{\Omega }}$.

The linear dispersion relations of the MRI for general axisymmetric perturbations in the Hall and AD regimes have been derived separately in Balbus & Terquem (2001) and in Kunz & Balbus (2004) and Desch (2004). The authors considered a general background field configuration ${\boldsymbol{B}}_0=B_{z0}{\boldsymbol{e}}_z+B_{\phi 0}{\boldsymbol{e}}_{\phi }$ and general axisymmetric perturbations of the form $\exp {({\rm i}{\boldsymbol{k}}\cdot {\boldsymbol{x}}+\sigma t)}$ with ${\boldsymbol{k}}=k_x{\boldsymbol{e}}_x+k_z{\boldsymbol{e}}_z$. The main results reveal that for linear MRI modes, the Hall term is coupled only to the vertical magnetic field, and the AD term is also coupled to the toroidal field. As a result, the presence of a background toroidal field has little effect on the Hall MRI, but it facilitates the MRI to operate in the AD-dominated regime with Am ≲ 1. A joint dispersion relation including all nonideal MHD terms was given by Pandey & Wardle (2012). They showed that while contributions from the Hall and AD terms are independent, the joint effect is that regimes stable to pure Hall MRI can be rendered unstable because of AD, a situation which again requires a net toroidal field and strong AD (Am ≲ 1).

Exploring the full parameter space of the MRI in the presence of the Hall and AD effects with different field orientations using nonlinear simulations is beyond the scope of this work. Here, we restrict ourselves to a pure vertical background field with either B_z0 > 0 or B_z0 < 0. This choice makes the dispersion relation much simpler, and the most unstable mode has a pure vertical wavenumber of k_z = k. For these modes, AD behaves the same way as Ohmic resistivity by replacing η_A with η_O (in the linear regime). This case also covers the most essential MRI physics relevant to PPDs: the Hall term is not directly coupled to the toroidal field, and for AD, the background toroidal field does not strongly affect the level of MRI turbulence for Am ≳ 1 (Bai & Stone 2011).

In reference to previous works (e.g., Wardle 1999), we show in Figure 1 the MRI growth rate for pure vertical modes k = k_z as a function of dimensionless wavenumber kv_A0/Ω and 1/χ₀, where subscript ₀ represents χ and v_A determined from the background field. Similarly, we use β₀ to denote the plasma β for the background field. The magnetic polarity is reflected from sgn(B_z0). We consider two cases with Am = 1 and Am = 100.

Zoom In Zoom Out Reset image size

For Am = 100 (very weak AD), the dispersion relation is well described by pure Hall MRI. For B_z0 > 0, the most unstable mode always has the maximum growth rate of 0.75 Ω⁻¹, and the most unstable wavelength λ_m shifts progressively to larger scales with $\lambda _m\propto \chi _0^{-1/2}$ as the Hall term strengthens (χ₀ → 0). Normalizing to disk scale height, we find

$$\begin{equation} \frac{\lambda _m}{H}\approx 4\pi \sqrt{\frac{l_H}{3H}}\left (\frac{2}{\beta _0}\right)^{1/4} \approx 0.5\sqrt{\frac{3l_H}{H}}\left (\frac{10^4}{\beta _0}\right)^{1/4}. \end{equation} \tag{ 10 }$$

For B_z0 < 0, unstable modes exist only when (1/χ₀) < 2, and the unstable wavenumber can extend virtually to infinity when (1/χ₀) > 1/2.

For Am = 1, we see that small-scale modes are strongly suppressed. For B_z0 < 0, the most unstable modes have wave numbers of k_mv_A/Ω ∼ 0.5. In the absence of the Hall effect (1/χ₀ = 0), λ_m is increased by a factor of ∼2 because of AD. For B_z0 > 0 and toward a stronger Hall term (1/χ₀ ≳ 5), λ_m is less affected by AD because it is shifted to larger scales, and the maximum growth rate is only slightly reduced.

Our unstratified shearing-box simulations mainly serve for calibrating and interpreting stratified simulation results. Therefore, we do not aim at a thorough parameter study but mainly focus on parameter regimes relevant to real PPDs. In this regard, we consider the following set of parameters.

• 1.  
  The Hall length l_H = 0.1H or 0.3H.
• 2.  
  Net vertical field strength, with β₀ = 10⁴ and 10⁵.
• 3.  
  Magnetic field polarity, B_z0 > 0 or B_z0 < 0.
• 4.  
  The value of Am = 1, occasionally 10 and 100.

Our simulations use a fixed box size of 4H × 4H × 2H in (x, y, z) dimensions. Note that our simulation box height is 2H rather than H typically used in unstratified shearing-box simulations, providing the potential to accommodate larger spatial structures while not being unrealistically tall for real disks. Our unstratified simulations can be performed with relatively high spatial resolution, 48 cells per H in the x−z plane (24 in the y dimension). We cannot afford the same resolution for our stratified runs in Sections 3–5. Therefore, we also conduct simulations with half the resolution to justify the use of lower resolution in our stratified simulations.

We have chosen the value of Am = 1 that is appropriate for the midplane region of the outer disk. From Equation (6), the Hall length of l_H ∼ 0.1 to 0.3H applies to the range of R ∼ 20 to 50 AU. Given β₀ = 10⁴ and 10⁵, the corresponding value of χ₀ ranges from 0.015 to 0.14.

For B_z0 < 0, and for this range of χ₀, there is no linearly unstable MRI mode. However, this does not necessarily relate to the nonlinear sustainability, given the relatively small value of l_H. Therefore, we first run simulations in the ideal MHD limit to time t = 60 Ω⁻¹, then we turn on nonideal MHD terms and evolve further to time t = 300 Ω⁻¹. In Figure 2, we show the time evolution of three runs with Am = 1, 10, and 100 at fixed l_H = 0.1 and β₀ = 10⁴. We see that for Am = 100, the MRI turbulence can be sustained but at a lower level, while for Am = 10 and 1, turbulence is suppressed. We have also tested with other values of β₀ and l_H, and we find that as long as Am = 1, no sustained MRI turbulence is possible. This implies that under this configuration the midplane region of the outer disk is likely the exact analog of the conventional "dead zone."

Zoom In Zoom Out Reset image size

For B_z0 > 0, the background field configuration is unstable to the MRI. We provide the list of runs and diagnostic quantities in Table 1. The runs are named in the form QxAyBz-Rw, where x = 10 l_H/H, y = Am, z = log₁₀(β₀), and w is the numerical resolution (24 or 48 per H). In all cases, we have fixed the value of Am = 1. We find that vigorous turbulence is quickly developed for all runs. Many of these runs show secular effects in their evolution (to be discussed later), so we run these simulations for a very long time to t = 1440 Ω⁻¹ and extract turbulence statistics by performing time and volume averages after t = 1120 Ω⁻¹ (denoted by the overline). Major diagnostic quantities include kinetic energy density $E_k=\overline{\rho v^2/2}$, magnetic energy density $E_M=\overline{B^2}/2$, Maxwell stress $\alpha ^{\rm Max}\equiv -\overline{B_xB_y}$, and Reynolds stress $\alpha ^{\rm Rey}\equiv \overline{\rho v_xv_y}$ (the normalization factor $\rho _0c_s^2$ is omitted because it is one in code units). The total Shakura–Sunyaev α is α^Max + α^Rey. Another useful diagnostic is α^mag ≡ α^Max/E_M (e.g., Hawley et al. 2011; Sorathia et al. 2012), which is considered a useful indicator for numerical convergence.

First, we find that for a relatively large l_H = 0.3H and a relatively strong field β₀ = 10⁴ a strong zonal field is gradually built up on relatively long timescales (∼10³ Ω⁻¹). This results from the concentration of the vertical magnetic flux pertaining to the Hall effect, as studied in detail in Kunz & Lesur (2013). In Figure 3, we show the final snapshot of our run Q3A1B4-R48 at time t = 1440 Ω⁻¹, which clearly shows the zonal field structure. On the other hand, we find that the zonal field coexists with vigorous turbulence and gives an α value of ∼10⁻³. The presence of vigorous turbulence, rather than remaining in the low-transport state, is largely due to relatively strong AD with Am = 1: additional magnetic dissipation acts against the buildup of magnetic flux (Kunz & Lesur 2013). We do not observe prominent zonal field structures in runs with smaller l_H and weaker magnetic fields.

Zoom In Zoom Out Reset image size

In the meantime, we find that in essentially all of our unstratified simulations, the density variation also shows significant zonal structure, leading to strong zonal flows to balance the pressure gradient of the zonal density variation (Johansen et al. 2009). Such density variation is not captured in Kunz & Lesur (2013) because of their usage of incompressible code. The density structure from our run Q3A1B4-R48 is shown in the bottom panel of Figure 3, where the density variation reaches ∼50%. As a result, the kinetic energy displayed in Table 1 is largely dominated by the kinetic energy associated with the zonal flow (v_y ∼ 0.2–0.3c_s). Other runs develop weaker zonal density variations and weaker zonal flows as well. These structures also develop over long timescales of ∼10³ Ω⁻¹ and show secular variations. A full discussion on such zonal flows is beyond the scope of this paper, but phenomenologically we observe from our unstratified simulations that a stronger zonal flow is launched when l_H is larger and when the background field is stronger.³

In all of our simulations, sustained MRI turbulence at the level of α ∼ 10⁻³ is obtained. A stronger background vertical field leads to stronger turbulence, and a larger l_H also leads to stronger turbulence until the zonal field configuration is developed, where the turbulence level is reduced. We caution that for the parameters considered here, the most unstable MRI mode is not well resolved. For the best-resolved case (run Q3A1B4-R48), we find from Equation (10) that the most unstable wavelength amounts to about 13 cells. We do not expect our simulations to show unambiguous convergence on the value of α (and in fact the value of α is also affected by the development of the zonal flows, which show long timescale variations). Nevertheless, by looking at the value of α_mag, we find that low- and high-resolution simulations give consistent values in all cases except for run Q1A1B5. Moreover, by inspecting snapshots in runs with different resolutions, we find that their evolutionary behaviors are qualitatively similar in all cases. This gives us confidence that the 24 cells per H adopted in our stratified runs is adequate to capture the essential properties of the MRI in the Hall–AD regime.

In sum, our unstratified simulations of the MRI in the presence of both the Hall effect and AD indicate that under conditions appropriate for the outer region of PPDs (Am ∼ 1), the MRI cannot be self-sustained in the midplane if B_z0 < 0, while for B_z0 > 0, the self-sustained turbulence exists at the level of α ∼ 10⁻³. When B_z0 > 0, we observe zonal fields develop when the Hall term and background field are relatively strong, and we find that zonal flows develop in all cases.

The global evolution of the system is largely controlled by magnetic fields, so we first discuss the overall evolution of a large-scale toroidal field from our simulations as a standard diagnostic. Starting from runs with β₀ = 10⁵, R30b5H+, R30b5H0, and R30b5H−, we show in Figure 4 the time evolution of the horizontally averaged B_y for all three runs. Because the initial conditions for these simulations are identical (except for a negative sign of B_z0 for run R30b5H−), these runs initially proceed in a similar way. The Hall and AD terms become progressively more important as (midplane) magnetic fields become stronger, and the three runs then evolve differently. All three cases show a prominent level of dynamo activities emanating from the surface layer, where the sign of mean B_y alternates over time. The alternation behavior is quite irregular and to some extent similar to ideal MHD simulations with a modestly strong vertical magnetic flux (β₀ ≳ 10³, Bai & Stone 2013a). It contrasts with the conventional MRI dynamo (zero net vertical magnetic field in ideal MHD) that exhibits very periodic cycles of about 10 orbits (e.g., Davis et al. 2010; Shi et al. 2010).

Zoom In Zoom Out Reset image size

We next discuss simulations with β₀ = 10⁴, with three runs R30b4H+, R30b4H0, and R30b4H−. Similar to the weaker field case, all three runs develop vigorous turbulence mainly in the surface layer because of FUV ionization (see next subsection). The time evolution of horizontally averaged B_y for the three runs is shown in Figure 5. We see that the MRI dynamo is suppressed in all cases, and the mean toroidal field is predominantly one sign. This is generally a consequence of a stronger net vertical field, which is an analog of the ideal MHD case (Bai & Stone 2013a). While the system is turbulent, the toroidal field is always the dominant field component, and when the dynamo is suppressed, this field component is dominated by the mean field. However, by viewing individual simulation snapshots, localized patches possessing opposite signs for the toroidal field do exist in runs R30b4H0 and R30b4H−. In the latter case, the region with opposite B_y gradually grows and eventually leads to the reversal of the mean toroidal field in the disk (bottom panel of the figure). We have continued this run further and found that the mean B_y will reverse again after another ∼300 Ω⁻¹, and this cycle is likely to continue. Similarly, the positive B_y region started to dominate the upper half of the disk near the end of our run R30b4H0.

Zoom In Zoom Out Reset image size

The secular evolution of the mean field discussed above exists in all of our simulations to a certain extent, and this is partly related to the limitations of the shearing-box approach: with an imposed net vertical field that presumably connects to infinity, the mean field in the disk should be in causal contact with the field beyond, but the causal connection is truncated with the prescribed outflow boundary condition. Because most activities in disks are magnetically driven, the secular evolution of the mean fields also makes the level of turbulence in disks time variable. For example, in run R30b4H−, the midplane region exhibits stronger turbulent activities around time t = 480–600 Ω⁻¹ with turbulent velocity about a factor of three higher than some other periods. Therefore, readers should bear in mind the potential uncertainties that are due to such variabilities.

Based on the time evolution of the mean field, we extract useful diagnostic quantities and average them in time from t = 480 Ω⁻¹ onward for simulations with β₀ = 10⁵, and from t = 360 Ω⁻¹ onward for simulations with β₀ = 10⁴. The mean vertical profiles of these quantities are shown in Figures 6 and 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The relative importance of various nonideal MHD effects can be best viewed from the top left panel of Figure 6 and the left panel of Figure 7. They show the profiles of the Elsasser numbers (based on the Hall-free run in each case, but the runs with the Hall term generally give almost the same profiles). Clearly, Ohmic resistivity is completely negligible with Λ ≫ 100 at all heights. With β₀ = 10⁵, both the Hall effect and AD are important within z ∼ ±2–2.5H with χ and Am being around one, and the range of influence of AD extends higher from the midplane than the Hall effect. The Hall effect is less important relative to AD with stronger net flux β₀ = 10⁴ because the resulting total field is stronger. Beyond z ∼ 2.5H, FUV ionization catches up and all nonideal MHD effects are greatly reduced. Beyond z = ±3H, the gas behaves essentially in the ideal MHD regime with Am > 100.

Vigorous MRI turbulence takes place beyond about z ∼ ±2.5H thanks to FUV ionization. As a result, the profile of the Maxwell stress $T_{R\phi }^{\rm Max}=-B_xB_y$ peaks at around z = ±3H, as shown in the top right panel of Figure 6 and the middle panel of Figure 7. Beyond z ∼ ±3H, the Maxwell stress drops because disk density drops, and it enters the magnetically dominated corona (plasma β < 1). All three runs at a given β₀ show very similar properties in this region because the gas behaves in the ideal MHD regime. Runs with β₀ = 10⁴ have a systematically higher Maxwell stress than the corresponding β₀ = 10⁵ runs by a factor of three to four as a result of a stronger background field.

The midplane region is where three simulations at fixed β₀ are expected to differ because of the Hall effect. The most prominent difference lies in the Maxwell stress. The runs with B_z0 > 0 give the highest stress that peaks at the midplane. This is related to the Hall-shear instability (Kunz 2008), which operates only when B_z0 > 0 and is responsible for generating stronger horizontal magnetic fields and hence Maxwell stress in the inner disk (Lesur et al. 2014; Paper I). Here, the effect is much less prominent than in the inner disk studied in Paper I and Lesur et al. (2014) because the Hall effect is only modestly significant (χ ∼ 1). The runs with B_z0 < 0 give the lowest midplane Maxwell stress, while the Maxwell stress from R30b5H0 (without the Hall term) lies in between. This is again consistent with the expectation from Paper I that the horizontal magnetic field tends to be reduced for negative B_z0.

As discussed in Section 2, for B_z0 > 0, the midplane region is unstable to the MRI, and the level of MRI turbulence is expected to be stronger than in the Hall-free case. For B_z0 < 0, self-sustained MRI turbulence is not expected because of the Hall effect. To characterize the level of turbulence, we consider the vertical component of the rms velocity, which is shown in the bottom left panel of Figure 6 and the right panel of Figure 7 for the two sets of runs. They are computed based on the turbulent kinetic energy at each height.⁴ In the same way, we define δv_z to be the rms vertical velocity fluctuation within z = ±2H for all of our runs, which is included in Table 2.

We see that the turbulent rms vertical velocity reaches ∼0.3–0.8c_s at the disk surface (z ∼ ±4H) for all of these runs, while it is reduced by more than one order of magnitude to ∼0.01–0.03c_s around the disk midplane. For β₀ = 10⁵, the run with B_z0 > 0 gives a higher midplane turbulent velocity, the run with B_z0 < 0 gives the lowest, and the Hall-free run lies in between, which is consistent with our expectation. Nevertheless, the difference is within a factor of two, so the role of the Hall effect in the midplane turbulent activities is only modest. While we caution that the level of turbulence in the B_z0 > 0 case may be underestimated because of the lack of numerical resolution, the overall scenario is similar to the Hall-free case: turbulence in the midplane region is substantially reduced largely because of AD, consistent with the earlier stratified AD simulations of (Simon et al. 2013a, where the midplane region was termed an "ambipolar-damping" zone). In the case of B_z0 < 0, because the MRI cannot be self-sustained at the disk midplane, the midplane turbulent motion is largely induced by the strong MRI turbulence from the disk surface layer. This is a direct analog of the conventional "Ohmic dead zone" of the inner disk (e.g., Fleming & Stone 2003; Turner & Sano 2008; Oishi & Mac Low 2009).⁵

For β₀ = 10⁴, we find that the level of midplane turbulence in all three runs are very similar (modulo some secular variations not reflected in the time-averaged plots), despite the marked difference in Maxwell stress. We have checked that for B_z0 > 0, the midplane Maxwell stress is dominated by contributions from the large-scale field ($-\overline{B}_x\overline{B}_y$), while for B_z0 < 0, the midplane Maxwell stress is almost entirely due to the turbulent field. The turbulent contributions of the midplane Maxwell stress from the two runs R30b4H+ and R30b4H− are in fact similar. We have also checked that for β₀ = 10⁵ the midplane Maxwell stress is always dominated by turbulent stress. The low level of turbulence in run R30b4H+ may be considered to be a consequence of the strong mean toroidal field ($\overline{B_y}$), which dominates the magnetic field strength and tends to suppress turbulent motions (but see also Section 4.4).

Overall, based on the six simulations with different strengths and polarities of the net vertical field, it is clear that the Maxwell stress profile (and hence the radial transport of angular momentum) is layered. Moreover, it appears that δv_z ≈ 0.01–0.02c_s is a good proxy for the level of turbulence in the midplane region of the outer disks, with much stronger turbulence in the FUV ionization layer at the disk surface.

Outflow is always launched in shearing-box simulations in the presence of a net vertical magnetic flux (e.g., Suzuki & Inutsuka 2009). Although this outflow may serve as a wind-launching mechanism, the kinematics of the outflow is not well characterized in shearing-box simulations because the rate of mass outflow does not converge with simulation box height (Fromang et al. 2013) and there are also symmetry issues (Bai & Stone 2013a). Therefore, we do not aim at fully characterizing the outflow properties but simply provide some basic diagnostics for reference. We calculate the rate of mass outflow leaving the simulation box $\dot{M}_{\rm out}$. This is computed by time averaging the sum of vertical mass flux at the two vertical boundaries. We also calculate the zϕ component of the Maxwell stress tensor $T_{z\phi }^{\rm Max}=-B_zB_\phi$, which determines the rate of wind-driven angular momentum transport (if the outflow is eventually incorporated into a global magnetocentrifugal wind). If the flow is laminar, T_zϕ can be conveniently evaluated at the base of the wind, where the toroidal velocity transitions from sub-Keplerian to super-Keplerian (Bai & Stone 2013b; Bai 2013). Because most of our simulation runs are highly turbulent at the disk surface, and there are ambiguities in defining the base of the wind (and whether the outflow can become a global wind at all; Bai & Stone 2013a), we simply provide a reference value of time-averaged $|T_{z\phi }^{\rm Max}|$ evaluated at z = ±4.5H in Table 2.

The value of the Shakura–Sunyaev α for a stratified disk can be written as

$$\begin{equation} \alpha =\frac{\int T_{R\phi }dz}{c_s^2\int \rho dz}, \end{equation} \tag{ 12 }$$

where T_Rϕ has contributions from both the Maxwell stress (− B_xB_y) and Reynolds stress (ρv_xv_y), leading to α^Max and α^Rey in Table 2. From the lower right panel of Figure 6, we see that the vertical profile of the Reynolds stress is generally a factor of several smaller than the Maxwell stress. Because of uncertainties in characterizing the outflow from shearing-box simulations, we truncate the vertical integral at z = ±4.5H in the equation above. For the six runs, the values of α are found to be around 1.5–2 × 10⁻³ for β₀ = 10⁵ and 5–6 × 10⁻³ for β₀ = 10⁴. As discussed in the previous subsection, a stronger net vertical field yields larger stress. These α values are also comparable to results from Simon et al. (2013a) with FUV penetration depth 0.01g cm⁻², but they have incorporated more realistic chemistry.

In a steady state, the total accretion rate driven by radial transport of angular momentum (given by α) and the putative wind-driven accretion (given by T_zϕ) can be approximately written as (e.g., Bai 2013)

$$\begin{eqnarray} \dot{M}&\approx& \frac{2\pi }{\Omega }\alpha c_s^2\Sigma +\frac{8\pi }{\Omega }R|T_{z\phi }|,\nonumber\\ \dot{M}_{-8}&\approx& 0.82\left (\frac{\alpha }{10^{-3}}\right)R_{\rm AU}^{-1/2} +4.1\left (\frac{|T_{z\phi }|}{10^{-4}\rho c_s^2}\right)R_{\rm AU}^{-3/4}, \end{eqnarray} \tag{ 13 }$$

where R_AU is the radius measure in AU, and we have assumed an MMSN disk model in the second equation, with $\dot{M}_{-8}$ being the accretion rate measured in 10⁻⁸ M_☉ yr⁻¹.

Using the values from Table 2 with R = 30 AU, we find that based on radial transport alone, the resulting accretion rate is about 0.24–0.33 × 10⁻⁸ M_☉ yr⁻¹ for the three runs with β₀ = 10⁵ studied here, which is somewhat smaller than desired. If there were contributions from disk wind, the estimated wind-driven accretion rate is about 0.7  ×  10⁻⁸ M_☉ yr⁻¹. The sum of the two contributions just matches the desired rate of 10⁻⁸ M_☉ yr⁻¹. For β₀ = 10⁴, the accretion rate resulting from radial angular momentum transport gives ∼0.72–0.91 × 10⁻⁸ M_☉ yr⁻¹, with a potential contribution from the wind to give ∼5  ×  10⁻⁸ M_☉ yr⁻¹.

For our 30 AU simulations, we find in using Equation (8) and from the Elsasser number plots in Figures 6 and 7 that l_H ≈ 0.2H around the disk midplane, which is about the threshold value to trigger the zonal field configuration in the unstratified case, as discussed in Kunz & Lesur (2013). In Section 2 we showed in Figure 3 that a strong zonal field and zonal flow are present in unstratified simulations when β₀ = 10⁴ and B_z0 > 0. To check whether our stratified simulations reveal similar behaviors, we show in Figure 8 the time evolution of mean gas density ρ and Bz for runs 30AUb4H± and 30AUb5H±, averaged in the y and z dimensions, within the disk region −2H ⩽ z ⩽ 2H.

Zoom In Zoom Out Reset image size

We find that, strikingly, for all runs, the vertical magnetic flux is concentrated into thin (axisymmetric) shells, while in regions outside these shells, the net vertical flux is close to zero. In the meantime, there are very prominent radial density variations that are characteristic of strong zonal flow. There is a clearly secular evolution of the vertical magnetic flux distribution and zonal flows; this is also related to the secular behaviors discussed in Section 4.1. At first glance, these features appear to be consistent with those shown in Figure 3 from our unstratified simulations. However, there are distinct differences. In particular, both the B_z0 > 0 and B_z0 < 0 cases show such zonal fields, while from unstratified simulations a zonal field is expected only from the B_z0 > 0 case. Also, the width of the zonal field is very small (<0.5H), while from unstratified simulations the width is generally wider than H.

In fact, we find that the concentration of magnetic flux is a generic behavior in shearing-box simulations with a net vertical magnetic flux. Not only simulations with the Hall effect, but also our Hall-free simulations at 30 AU, together with many simulations at other disk radii, show this behavior to some level. We also find that the concentration is less prominent when the net vertical field is weaker, as one compares the top and bottom panels in Figure 8. Accompanying the magnetic flux concentration is the strong zonal flow, the density of which varies across the domain up to ∼30%. An enhanced zonal flow in the presence of a net vertical magnetic flux was reported in Simon & Armitage (2014) based on stratified shearing-box simulations in the AD-dominated outer disk. Such zonal flows also exist in our earlier simulations including both Ohmic resistivity and AD closer in (at 10–20 AU; Bai 2013), and we have verified that, in general, there is only one single "wavelength" of the density/pressure variation across the radial domain, regardless of the radial domain size (Bai 2013, unpublished).

In our follow-up study (Bai & Stone 2014), we show that in the presence of a net vertical magnetic field, the magnetic flux concentration is a natural consequence of the MRI itself. Regions with a strong magnetic flux yield a larger Maxwell stress, which drives enhanced zonal flows with flux-concentrated regions corresponding to density minima. AD further enhances the concentration to yield sharp vertical flux profiles. From Figure 8, we see that the location where magnetic flux concentrates significantly correlates with the density minimum, especially in runs R30b4H+ and R30b5+. The correlation is less pronounced in runs with B_z0 < 0, which reflects the fact that MRI is not operating in the midplane. Moreover, because the surface layer is fully MRI turbulent, the magnetic flux distribution is also likely affected by the surface MRI in the ideal MHD regime.

In sum, the zonal field and zonal flow observed in our stratified simulations are not due to the Hall effect as reported in unstratified simulations, but are correlated phenomena generic to the MRI with a net vertical magnetic flux. Although the saturation of the zonal flow is artificially affected by the simulation box size, its association with the magnetic flux concentration is very likely a physical phenomenon. Our local simulations here serve as a first study of the PPD gas dynamics including all nonideal MHD effects, and more works, especially global simulations, are needed to quantify their properties.

For run R5b5H+, we see from the top panel of Figure 9 that the system is able to achieve a largely laminar state as desired. More interestingly, the toroidal field changes sign almost exactly at the disk midplane, automatically maintaining the reflection symmetry (more specifically, even-z symmetry; see Figure 9 of Bai & Stone 2013b). Achieving this field geometry is essential for physically launching a magnetocentrifugal wind. Our full 3D simulation therefore supports the procedure adopted in Paper I where the reflection symmetry across the midplane was enforced.

Checking the time-averaged vertical profiles of the main diagnostic quantities, we find that the result is almost identical with Figure 9 of Paper I (with a slight difference because our box extends to z = 6 H rather than 8 H). For this solution, the horizontal magnetic field near the midplane is strongly amplified by the Hall shear instability, and the flip of this horizontal field creates a strong current density at the midplane. This contrasts with the study by Bai (2013), where without the Hall term, the strong current layer was found to be located offset from the midplane at z_SC ≈ 1.3H in this particular case (see his Table 2 for run S-R5-b5). It appears that with the Hall term included, the horizontal magnetic field tends to flip right across the midplane rather than from the upper layers.

In Figure 10, we further show the vertical profiles of time-averaged Maxwell stress and vertical turbulent velocities. For run R5b5H+, the Maxwell stress profiles peak close to the disk midplane at a rather high level close to $10^{-2}\rho _0c_s^2$ as a result of the Hall shear instability. The dip at the midplane is due to the flip of the horizontal field, all in agreement with the results in Paper I. For the turbulent velocity profiles, however, we find that an appreciable level of turbulence is present in this run. The turbulent velocity is again on the order of 0.01c_s around the midplane and increases toward the surface layer at a level very similar to that in the outer disk studied in the previous section. Since we expect the system to be stable to the MRI, the turbulence mainly originates from elsewhere: at the midplane, we find that the strong current layer tends to exhibit small amplitude corrugation from time to time, resembling the tearing modes in reconnecting with the current sheet. Such corrugating motion is likely the source of most random velocities that propagate toward disk surface layers and become amplified because of the rapid density drop.

Zoom In Zoom Out Reset image size

In sum, for B_z0 > 0, our 3D simulation with the full box well reproduces the quasi-1D simulations with enforced reflection symmetry in Paper I: we expect that accretion is mainly driven by the magnetocentrifugal wind, with a significant contribution from the radial transport of angular momentum via the large-scale Maxwell stress and magnetic braking (see Table 2 of Paper I). The wind-driven accretion flow mostly proceeds in the strong current layer where the toroidal magnetic field flips (Bai & Stone 2013b), and here it takes place exactly at disk midplane. Our 3D simulation further reveals the presence of turbulence, which largely originates from the midplane region where the relatively strong large-scale horizontal magnetic fields flip. The level of turbulence is similar to that in the outer disk. We also comment that because the system is stable to the MRI, magnetic flux concentration into thin shells is not observed in this simulation.

For run R5b5H−, the system is expected to be unstable to the MRI in a narrow range of disk height at about |z| ∼ 2–3H. This can roughly be identified from the left panel of Figure 9 in Paper I, where the Hall Elsasser number χ passes one at around z = 2.5H, and the plasma β there is still much larger than unity (based on the Hall-free run in dashed lines). A detailed explanation of the onset of the instability is given in Section 5.2 of Paper I, but, in brief, it is related to the fact that for B_z0 < 0 the Hall term makes the most unstable MRI wavelength shift to a shorter wavelength when the Elsasser number χ₀ is of order unity, allowing the unstable modes (that are otherwise too long) to fit into the disk. Using full 3D simulations, we see from the middle panel of Figure 9 that the large-scale toroidal magnetic field flips in a highly periodic manner, and the origin of the periodic flips directly connects to the unstable region. Interestingly, the toroidal fields in the upper and lower halves always have opposite signs, and the midplane horizontal field is very weak (and goes through zero). We have also found that the overall mean field evolution can be almost exactly reproduced from our quasi-1D simulation of Paper I. An outflow is launched whose mass outflow rate is smaller than but the same order of magnitude as the rate from our run R5b5H+ (see Table 2). Therefore, at a given time, the magnetic field configuration can be considered physical for a magnetocentrifugal wind. However, because the toroidal (and hence radial) field constantly changes sign, the wind keeps oscillating between the radial inward and outward directions. This is inconsistent with global wind geometry and reflects the limitation of the local shearing-box framework (Bai & Stone 2013a). While periodic field flips are likely physical phenomena that are inherent with the onset of the MRI, global simulations are necessary to determine the fate of the outflow.

The onset of the MRI also leads to some level of turbulence, as seen from the bottom panel of Figure 10. Away from the region where the MRI operates, turbulent motion largely results from a passive response to the MRI activities, and the midplane has the weakest level of turbulent motion. Despite different origins, the level of turbulence is comparable to run R5b5H+, especially at the surface.

The fact that the mean toroidal field periodically changes sign makes it ambiguous in estimating the role of disk wind in transporting angular momentum (the net wind-driven accretion rate would be zero considering the periodic flips). Here we set it aside and look at the radial transport of angular momentum from the Maxwell stress, as shown in the top panel of Figure 10. We see that the Maxwell stress peaks at about |z| ∼ 4H, but at a relatively low level. We estimate the total α to be only about 4.5 × 10⁻⁴, corresponding to an accretion rate of ∼1.6 × 10⁻⁹ M_☉ yr⁻¹ using Equation (13). This is about an order of magnitude smaller than the expected level of 10⁻⁸ M_☉ yr⁻¹.

We further performed run R5b4H− with a stronger net vertical field β₀ = 10⁴. Based on Paper I, we expect the system to be stable to the MRI and to develop a laminar magnetocentrifugal wind. This is again confirmed using full 3D simulations, with the general wind properties almost identical to the one obtained in Paper I. In particular, our full 3D run automatically obeys the reflection symmetry across the midplane, confirming that solutions with enforced symmetry in Paper I are generally physical. Note that the toroidal field is close to zero near the midplane as a result of the Hall term. The level of random motion in our run R5B4H− is systematically weaker than all other runs, confirming its intrinsically laminar nature. One can obtain from Table 2 the Maxwell stress and the wind stress to derive the accretion rate resulting from the radial transport and wind, or directly look at Table 2 of Paper I for more accurate estimates. We see that the radial transport is completely negligible compared with the wind-driven accretion rate. The latter gives ∼10⁻⁷ M_☉ yr⁻¹, which is an order of magnitude more than sufficient.

In sum, it appears that for B_z0 < 0, while results from our shearing-box simulations are likely robust, they also raise puzzling issues regarding the mechanism to transport angular momentum. For a relatively weak net vertical field (β₀ ∼ 10⁵), MRI sets in, leading to a periodically oscillating outflow. Based on shearing-box simulations, we are unable to tell whether the outflow drives angular momentum transport, but radial transport of angular momentum by Maxwell stress appears to be too inefficient. For a relatively strong net vertical field (β₀ ∼ 10⁴), the system unambiguously launches the magnetocentrifugal wind with a physical wind geometry, but it drives a more rapid accretion than typically observed. At this point, it is unclear how the system can achieve an accretion rate at the desired level of ∼10⁻⁸ M_☉ yr⁻¹, an issue that can only be clarified from global simulations.

At 15 AU, our quasi-1D simulations suggest a laminar configuration for B_z0 > 0 with β₀ = 10⁴, while MRI should set in otherwise. In Figure 11, we show the space–time plot of the horizontally averaged toroidal field. In Figure 12, we further show the time-averaged profiles of Maxwell stress and vertical turbulent velocity for all four runs, where the time averages are taken from time t = 420 Ω⁻¹ onward. We see that for all four runs the system eventually settles into a state where the large-scale toroidal field remains one sign across the entire disk, so the symmetry of the outflow would be undesirable for a global wind. Nevertheless, we again set aside the issue of symmetry and focus on other properties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For B_z0 > 0, and comparing runs R15b5H+ with R15b4H+, it is counterintuitive to notice from both figures that a stronger mean toroidal magnetic field is generated when the net vertical field is weaker (R15b5H+), leading to stronger Maxwell stress around disk midplane. A look at all of the simulation data reveals that for run R15b4H+ almost the entire vertical magnetic flux is concentrated in a single thin shell, while the rest of the radial zones have a net vertical flux that is effectively zero. As a result, magnetic field amplification by the Hall shear instability is suppressed for the bulk of the disk. A strong zonal flow is also formed with a high density contrast of 30%, and the shell of magnetic flux coincides with the density minimum. The highly nonuniform distribution of magnetic flux also makes the flow properties in this run deviate from the wind solution in Paper I (see his Table 2). On the other hand, for run R15b5H+, the magnetic flux distribution is much more uniform, leading to effective growth of the horizontal magnetic field that is due to Hall shear instability, producing stronger Maxwell stress at the disk midplane. At this point, it is again unclear how realistic the level of magnetic flux concentration is, and the results shown here should be treated with caution.

For B_z0 < 0, we see that for both runs with β₀ = 10⁴ and 10⁵, the initial evolution of mean B_y closely resembles our run R5b5H− with quasi-periodic flips. This is again because the MRI sets in a thin layer where the Hall Elsasser number transitions through order unity. Later on, a field of one sign takes over and dominates the entire disk. There are also MRI activities in the FUV layer, though the level is weaker than in their 30 AU counterpart (e.g., seen from the peak Maxwell stress). To some extent, this location represents a transition between the 5AU and 30 AU cases: in the former case, the MRI is triggered mainly in the Hall-dominated layer, and in the latter case, the MRI is active mainly in the FUV layer. As usual, the midplane horizontal magnetic field is suppressed due to the Hall effect, and most of the Maxwell stress originates from the FUV layer.

From the values of α^Max and $T_{z\phi }^{\rm Max}$ listed in Table 2 and using Equation (13), we see that the net vertical magnetic flux has to be at least β₀ = 10⁴ in order for the accretion rate to reach levels comparable to 10⁻⁸ M_☉ yr⁻¹. On the other hand, if a magnetocentrifugal wind is operating, the level of $T_{z\phi }^{\rm Max}$ from a weak net vertical field with β₀ = 10⁵ is sufficient to drive the accretion rate above the desired level. Overall, the turbulent velocity is smallest at midplane either because of weak MRI turbulence (B_z0 > 0 with weak field) or induced random motion from MRI activities in the disk surface (B_z0 < 0), similar to the 30 AU case.

At 60 AU, the relative importance of the Hall effect is reduced by a factor of about two compared with the 30 AU case (see Equation (9)), and it is only marginally important at the disk midplane. AD is the dominant effect in most regions of the disk. Also, given the approximately constant penetration column density of the FUV ionization, it penetrates deeper geometrically (in terms of disk scale height) at the more tenuous outer disk. In Figure 13, we show the space–time plot of the horizontally averaged toroidal field. In Figure 14, we further show the time-averaged profiles of Maxwell stress and vertical turbulent velocity for all four runs, where the time averages are taken from time t = 300 Ω⁻¹ onward. The general evolution of the systems are in many ways similar to our runs at 30 AU, where the MRI drives vigorous turbulence in the surface FUV layer, with the midplane region only weakly turbulent. Here we mainly focus on the differences and the overall trend toward larger disk radii.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At β₀ = 10⁵, dynamo activities constantly flip the mean toroidal field similar to but apparently more regularly than in the 30 AU case for both magnetic polarities. For β₀ = 10⁴, the dynamo is suppressed, and the entire disk is dominated by a mean toroidal field with a single sign. When B_z0 < 0, we do not observe the mean field changing sign as for the 30 AU counterpart shown in Figure 5. In fact, the toroidal field in the entire disk has the same sign throughout the saturated state of the simulation, so we do not expect this sign flip to occur. We speculate that the flip we observed at 30 AU is associated with the relatively strong Hall effect at the disk midplane, but it is unlikely to occur toward the outer disk as the Hall effect becomes less dominant.

At 60 AU, the contrast in Maxwell stress between the B_z0 > 0 and B_z0 < 0 cases at disk midplane is still very evident. The level of turbulence is found to be higher for runs with a weaker net vertical field β₀ = 10⁵, which may be because in runs with β₀ = 10⁴ the turbulent motion is limited by the relatively strong large-scale toroidal field, but it may also be due to the strong concentration of magnetic flux into thin shells where a large fraction of the simulation domain has effectively zero net vertical flux.

Deeper penetration of FUV ionization allows the MRI to be fully active over thicker surface layers. With a larger gas density in the active layer (in code units), the turbulent magnetic fields from the MRI become stronger (at fixed β₀) than their 30 AU counterparts, giving larger values of α^Max. Again, we find that for the Maxwell stress alone to drive the accretion rate of ∼10⁻⁸ M_☉ yr⁻¹, the net vertical flux needs to be β₀ ∼ 10⁴ or stronger. The magnetocentrifugal wind, if operating in the outer disk, would drive an accretion with rate ∼0.4–4  ×  10⁻⁸ M_☉ yr⁻¹ for β₀ = 10⁵–10⁴.

In this work, we have studied the gas dynamics of PPDs, focusing on regions toward the outer disk (from 5 to 60 AU), taking into account all nonideal MHD effects in a self-consistent manner. In these regions, the Hall effect generally dominates near the disk midplane, AD plays an important role over a more extended region across the disk height, and the very surface layer behaves in the ideal MHD regime because of FUV ionization. In the presence of the Hall effect, the gas dynamics depends on the polarity of the net vertical magnetic field (B_z0) threading the disk relative to the rotation axis (along $\hat{z}$). Because the Hall effect becomes progressively less important relative to AD with increasing disk radius, we estimate based on the MMSN disk model that the Hall effect controlled polarity dependence extends to about 60 AU.

We first conducted unstratified MRI simulations that included both the Hall effect and AD. We find that at the conditions expected in the outer region of PPDs (the midplane plasma β₀ for the net vertical field being 10⁴–10⁵) the MRI leads to turbulence when B_z0 > 0, but the MRI cannot be self-sustained when B_z0 < 0. For B_z0 > 0, the level of MRI turbulence is of the order α ∼ 10⁻³ (with AD Elsasser number Am = 1). We confirm that the strong zonal field configuration of Kunz & Lesur (2013) can be achieved with a sufficiently strong Hall effect, and we find that in the meantime it leads to strong zonal flows. In addition, a numerical resolution of 24 cells per H = c_s/Ω is in general adequate to resolve the bulk properties of the MRI turbulence.

We then focus on self-consistent stratified MRI simulations at various disk radii from 5 AU to 60 AU, with the main results summarized as follows.

At a relatively small disk radius (∼5 AU), and for B_z0 > 0, we confirm and justify the results from Paper I that the system launches a strong magnetocentrifugal wind and is able to achieve a physical wind geometry, with the horizontal magnetic field flipped exactly at the disk midplane. Although Maxwell stress is enhanced because of the Hall shear instability, accretion is largely driven by the wind and proceeds primarily through the midplane strong current layer. In addition, the midplane region is weakly turbulent, which likely results from the flip of the relatively strong horizontal magnetic field. The turbulent motion gets amplified toward the disk surface as gas density drops.

For B_z0 < 0, our full 3D simulations confirm the results from Paper I that the system is unstable to the MRI in a thin Hall-dominated layer when the net vertical field is relatively weak (β₀ = 10⁵). This results in periodic flips of the large-scale horizontal magnetic field over time, accompanied by a radially oscillating disk outflow or wind. The fate of the outflow or wind and whether it drives accretion are uncertain based on shearing-box simulations. Radial transport of angular momentum by Maxwell stress is found to be too inefficient by one order of magnitude. A stable magnetocentrifugal wind with a physical wind geometry can be achieved with a stronger net vertical field (β₀ = 10⁴), which very efficiently drives accretion with $\dot{M}\gtrsim 10^{-7}\, M_{\odot }$ yr⁻¹. It is uncertain whether and how the system can achieve the typically observed rate of 10⁻⁸ M_☉ yr⁻¹.

At a relatively large disk radius (∼30 AU), we find that the Hall effect mainly affects the Maxwell stress at disk midplane, with B_z0 > 0 (B_z0 < 0) giving an enhanced (reduced) stress that is similar to that found at the inner disk (Paper I; Lesur et al. 2014). Nevertheless, the strongest Maxwell stress results from vigorous MRI turbulence in the surface layer because of FUV ionization (Perez-Becker & Chiang 2011; Simon et al. 2013a). While self-sustained MRI is expected at disk midplane when B_z0 > 0 but not when B_z0 < 0, the level of turbulence in all cases appears very similar, with a vertical turbulent velocity of the order δv_z ∼ 0.01–0.03c_s. The turbulent motion in the latter case is largely induced from stronger turbulence in the surface layer, analogous to the conventional "Ohmic dead zone" picture (e.g., Fleming & Stone 2003). Overall, the gas dynamics in the outer regions of PPDs shows a clear layered structure consisting of a highly turbulent surface FUV ionization layer and a weakly turbulent midplane region that is due to damping through a combination of AD, the Hall effect, and large-scale magnetic field structure.

We find that for a relatively weak net vertical field (β₀ = 10⁵), the MRI dynamo leads to repeated flips of the large-scale toroidal field with very irregular cycles. Dynamo activities tend to be suppressed for stronger fields (β₀ = 10⁴). Our simulations also show secular behavior in the evolution of the mean toroidal field, especially in simulations at 30 AU. This is to a certain extent related to the limitations of the shearing box because the net vertical magnetic flux ought to be connected to infinity but gets truncated by the vertical boundary condition without reaching all of the critical points (e.g., Fromang et al. 2013).

We also find that most of our simulations show a strong concentration of vertical magnetic flux into a thin axisymmetric shell at a certain radial location, and the rest of the regions have a net vertical flux that is close to zero. The concentration is generally stronger in simulations with a stronger net vertical field (β₀ = 10⁴) and toward outer disk radii (≳ 15 AU). Accompanying the magnetic flux concentration is an enhanced density variation (zonal flow) across the radial domain, with most flux concentrated in low-density regions. The concentration differs from the zonal field from unstratified Hall-MRI simulations (Kunz & Lesur 2013), but it is a generic property of the MRI turbulence in the presence of a net vertical magnetic flux, which is studied in detail in our followup paper (Bai & Stone 2014). Global simulations are required to determine the saturation amplitude of the concentration and transport of magnetic flux across the disk.

Although all of our simulations launch disk outflows, it is uncertain whether such outflows (at ≳ 15 AU) can be incorporated into a global magnetocentrifugal wind because of the MRI dynamo and symmetry issues (Bai & Stone 2013a), but if they do, the level of net vertical flux of β₀ = 10⁵ and stronger is generally sufficient to drive accretion at the desired level of 10⁻⁸ M_☉ yr⁻¹. On the other hand, to purely rely on radial transport of angular momentum by Maxwell and Reynolds stresses, the level of the net vertical field must be β₀ = 10⁴ or stronger, assuming an MMSN disk model. This level of field translates to physical field strength according to

$$\begin{equation} B=18.6\beta ^{-1/2}R_{\rm AU}^{-13/8}\ {\rm G}. \end{equation} \tag{ 14 }$$

For reference, we find for β₀ = 10⁴ that B_z0 ∼ 0.7 mG at 30 AU.

Combining the results from this paper and Paper I, we see that the Hall effect has a major influence on the disk dynamics in the inner region of PPDs (≲ 15 AU), where the wind properties, stability to the MRI, and the amplification or reduction of the midplane horizontal magnetic field show a strong dependence on the polarity of the net vertical magnetic field. The Hall effect also affects the stability to the MRI in the midplane region of the outer disk, although the level of midplane turbulence appears insensitive to the Hall effect. Overall, it is likely that wind-driven accretion dominates the inner disk, and accretion can be largely driven by the MRI through the surface FUV layer in the outer disk. This picture was already outlined in the discussion of Bai (2013), which incorporated numerical simulation results without the Hall effect (Bai & Stone 2013b; Simon et al. 2013a). On the other hand, the detailed behavior in the inner disk region, as well as the transition from the largely laminar inner disk to the MRI turbulent outer disk, is expected to have a strong polarity dependence because of the Hall effect, as summarized in the previous subsection and also in Paper I and Lesur et al. (2014).

Several observational consequences are expected based on our current simulation results. First, the fact that the inner disk launches a magnetocentrifugal wind can be detectable through gas tracers. In fact, signatures of low-velocity disk outflow have been routinely inferred from blue-shifted emission line profiles, such as from CO, O i, and Ne ii lines (e.g., Hartigan et al. 1995; Pascucci & Sterzik 2009; Pontoppidan et al. 2011; Herczeg et al. 2011; Sacco et al. 2012; Rigliaco et al. 2013). While conventionally interpreted as signatures of photoevaporation (e.g., Gorti et al. 2009; Owen et al. 2010), a magnetocentrifugal wind is likely to produce similar signatures because it possesses low velocities near the launching point before getting strongly accelerated and diluted. In reality, both mechanisms are likely to cooperate to launch the disk wind. We note that a purely photoevaporative wind is likely to be angular-momentum conserving because the radial driving force does not exert any torque on the outflow, while a strongly magnetized magnetocentrifugal wind is more likely to be angular-velocity conserving near the base of the wind where the gas is forced to move along supra-thermal magnetic fields anchored to the disk (e.g., Spruit 1996). Searching for distinguishable signatures between the two scenarios would be important for understanding the nature of the observed disk outflows.

Second, we expect the level of turbulence in the outer disk to be layered, where the level of turbulence is expected to be of the order δv_z ∼ 10⁻²c_s at the midplane and increases to near sonic level toward the disk surface. An empirical constraint on the level of turbulence in the outer regions of PPDs has already been reported based on the turbulent line width of the CO (3–2) transition (Hughes et al. 2011). This line is optically thick and probes the disk surface layer with a line width constrained to be ≲ 10%–40% of sound speed, consistent with a fully turbulent surface layer. With superb sensitivity and resolution, ALMA is expected to constrain the variations in turbulence level at different disk heights using different line tracers. Together with theoretical efforts in modeling turbulent line width (e.g., Simon et al. 2011, 2013b), we expect ALMA data to provide direct evidence of a layered structure in the outer PPDs.

Third, the weakly turbulent outer disk with a toroidal-dominated field configuration may lead to grain alignment and dust polarization (Cho & Lazarian 2007). We have found that in the outer disk (≳ 30 AU) the net vertical field needs to be β₀ ∼ 10⁴ or stronger for Maxwell stress to drive an accretion rate of 10⁻⁸ M_☉ yr⁻¹. For this level of net vertical field, we see that the MRI dynamo is suppressed, and the entire field is dominated by a large-scale toroidal magnetic field, whose strength at the disk midplane corresponds to a plasma β ∼ 10–20 (e.g., see Figure 13). Using Equation (14), we find that the midplane toroidal field can be at least ∼3–8 mG at 60–100 AU. Based on Equation (1) of Hughes et al. (2009), and using the MMSN disk model at midplane with grain size of 10–100 μm and dust aspect ratio s = 3, we find that the critical strength for grain alignment to occur is ∼1–40 mG at 60–100 AU. Although there are large theoretical uncertainties, we see that the match is marginal, and the field strength in the outer disk can either be just enough to promote grain alignment, or a little too weak to align the grains. Several observational attempts to search for dust polarization in Class II disks have failed (Hughes et al. 2009, 2013). Very recently, however, the successful detection of dust polarization toward younger sources has been reported, with an inferred field configuration resembling a large-scale toroidal field (Rao et al. 2014; Stephens et al. 2014). This might indicate that the disk magnetic field fades over time. Again, future dust polarization observations by ALMA will likely provide better constraints on the geometry, strength, and evolution of disk magnetic fields.

From this work, together with Paper I, we have explored the main parameter space in the gas dynamics of PPDs using local shearing-box simulations. There are other unexplored parameters and uncertainties, including the abundance and size distribution of grains, where small grains may reduce the importance of the Hall effect and AD and hence promote the MRI (Bai 2011a, 2011b). Also, the cosmic-ray ionization rate may be reduced and modulated by the stellar wind (Cleeves et al. 2013; it is less affected by disk wind), the X-ray luminosity can be highly variable because of stellar flares (Wolk et al. 2005; Ilgner & Nelson 2006), and FUV photons may be shielded by the dust in the disk wind from the inner disk (Panoglou et al. 2012; Bans & Königl 2012). It is likely that grain abundance and FUV ionization are more sensitive parameters (Bai & Stone 2013b; Simon et al. 2013a; Paper I), and X-ray ionization is less sensitive but also important (Bai 2011a; Paper I).

Probably the largest uncertainties in our work come from the use of a local shearing-box framework, and there are several outstanding issues related to net vertical magnetic flux. With a net vertical flux, it is well known that the properties of the disk outflow are not well characterized in shearing-box simulations, largely because the vertical gravitational potential is ever increasing in the local approximation (Fromang et al. 2013; Bai & Stone 2013b). Issues related to the symmetry and fate of the outflow are notorious (Bai & Stone 2013a, 2013b). Moreover, the evolution of a large-scale magnetic field can be affected by the vertical outflow boundary condition. In addition, magnetic flux concentration and zonal flows are again not well characterized in shearing-box simulations (Bai & Stone 2014). Global disk simulations with vertical stratification and net vertical magnetic flux have recently been carried out (Suzuki & Inutsuka 2014), yet many of these issues remain not quite addressed because of limited domain size in the θ dimension. In the future, it is crucial to perform global simulations with a sufficiently large vertical domain to accommodate the disk outflow or wind and the fine resolution in the disk to resolve the disk microphysics.