Origin of Weak Turbulence in the Outer Regions of Protoplanetary Disks

Jacob B. Simon; Xue-Ning Bai; Kevin M. Flaherty; A. Meredith Hughes

doi:10.3847/1538-4357/aad86d

1. Introduction

A long standing question in accretion disk theory is how exactly angular momentum is removed from orbiting gas, allowing it to accrete onto the central object. While this issue is universal to all accretion disks, understanding angular momentum transport in protoplanetary disks is particularly crucial to understanding planet formation. Indeed, whether angular momentum is transported via turbulence or some other process has significant implications for a range of planet formation stages, including the settling and growth of dust grains (Fromang & Papaloizou 2006; Youdin & Lithwick 2007; Birnstiel et al. 2010), the concentration of particles (Cuzzi et al. 2008; Johansen et al. 2009; Simon & Armitage 2014), and the migration of planetary bodies (Nelson & Papaloizou 2004; Baruteau et al. 2011; Lubow & Ida 2011; Paardekooper et al. 2011).

For many years, it has been theorized that turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1998) is responsible for angular momentum transport in such disks; angular momentum is redistributed radially via correlated turbulent fluctuations in the radial and azimuthal components of both magnetic fields and gas velocities. Recently, however, protoplanetary disk theory has undergone a paradigm shift. In particular, studies including nonideal magnetohydrodynamic (MHD) effects other than Ohmic resistivity (as conventionally considered, e.g., Gammie 1996) resulting from low ionization fractions have found that MRI-driven turbulence is largely reduced at the midplane in the outer disk (Simon et al. 2013a, 2013b), whereas it remains active in the upper layers due to ionization from FUV photons (e.g., Perez-Becker & Chiang 2011; Simon et al. 2013a). In the inner disk, MRI-driven turbulence can be either reduced significantly (Simon et al. 2015b) or quenched altogether (e.g., Bai & Stone 2013; Lesur et al. 2014), depending on the precise levels of ionization (Simon et al. 2015b). Furthermore, angular momentum can be additionally removed vertically through a magnetic wind (Salmeron et al. 2007; Suzuki & Inutsuka 2009; Bai 2017; Béthune et al. 2017), akin to the Blandford–Payne process (Blandford & Payne 1982) and radially through laminar torques (e.g., Lesur et al. 2014).

With the advent of ALMA, we can now directly test these ideas with high spatial and spectral resolution data. In our recent work, Flaherty et al. (2015, 2017, hereafter F15 and F17, respectively), we analyzed ALMA data using Markov chain Monte Carlo techniques to put a strong upper limit on turbulent velocities in HD 163296. We found that ${v}_{\mathrm{turb}}\lesssim 0.05{c}_{{\rm{s}}}$ throughout the entire disk column in the regions most easily resolved with ALMA (∼30–100 au and beyond), corresponding to a Shakura–Sunyaev α < 10⁻³. While these constraints are consistent with theoretical expectations from numerical simulations for the midplane region, where ambipolar diffusion damps MRI-driven turbulence (Simon et al. 2015a, hereafter S15), they are in fact inconsistent with predictions that ${v}_{\mathrm{turb}}\sim 0.1\mbox{--}1{c}_{{\rm{s}}}$ in the upper layers of the disk, where FUV ionization is sufficiently strong for the MRI to operate (e.g., Perez-Becker & Chiang 2011; Simon et al. 2013a). Furthermore, such a low α is inconsistent with the observed accretion rates (assuming a steady-state accretion flow driven by turbulence) onto the central star, $\dot{M}\sim {10}^{-7}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Mendigutía et al. 2013).

In addition, other observations and analyses have pointed toward weak turbulence in the outer regions of protoplanetary disks. We have recently found that TW Hya also exhibits weak turbulence at large radial distances, with ${v}_{\mathrm{turb}}\lesssim 0.08{c}_{{\rm{s}}}$ (Flaherty et al. 2018). Furthermore, modeling of the dust in the HL Tau system by Okuzumi et al. (2016) and Pinte et al. (2016) has suggested that turbulence in this system must be weak (Pinte et al. 2016 estimates α ∼ 10⁻⁴) in order for the rings to achieve their observed structure. If indeed turbulent velocities are as weak as these studies suggest, the inconsistencies between observations and theoretical predictions call into question the idea of turbulence-driven angular momentum transport. The question that then remains is: can magnetic winds drive accretion in a laminar fashion?

In this paper, we carry out a series of local (i.e., small corotating disk patch) simulations to explore accretion driven both by magnetically launched winds and turbulence, focusing in particular on the observational constraints from the HD 163296 disk. We find that, only in the limit of weak ionization and a weak vertical magnetic field threading the disk can we reproduce the small turbulent velocities as measured by F15 and F17. In particular, even when there is a substantial magnetic wind, the gas flow is not laminar, but instead quite turbulent. Such turbulence, if present in HD 163296, would have been detected in the data from F15 and F17. These results suggest that the outer disk is only very weakly accreting, if at all. One possible solution to the low-turbulence issue (as we discuss in detail below) is that ionizing radiation is being blocked from reaching the outer disk, where only a very weak magnetic field is present. Along these lines, we construct a preliminary model in which a large-scale magnetic field, mostly confined to small radial distances, can block ionizing radiation from the outer disk. We use this model to explain current observational constraints and make testable predictions for future observations, both of HD 163296 and other systems.

The outline of our paper is as follows. In Section 2, we describe the details of our calculations, including the algorithm and parameters chosen. In Section 3, we present our main results, and we follow up with an in-depth analysis of these results in Section 4. We then discuss the implications of our study and our preliminary model to explain the observations in Section 5. Finally, we wrap up with our main conclusions in Section 6.

2. Method

To carry out our investigation, we simulate the evolution of a small accretion disk patch using the local, shearing-box approximation (Hawley et al. 1995), including vertical stratification and all three nonideal effects (Ohmic and ambipolar diffusion, and the Hall effect). The local simulations are centered at 100 au in a model for the disk around HD 163296 (Rosenfeld et al. 2013), following the approach of S15. We employ the Athena code (Stone et al. 2008) with the implementation of ambipolar diffusion and the Hall effect described by Bai (2011, 2014), respectively. Furthermore, we use Athena in a configuration identical to that in S15. Following that work, we obtain the diffusivities using a look-up table calculated via a chemical network as described further in Bai & Stone (2013) and Bai (2013, 2014). We include cosmic-ray ionization at a rate ${\xi }_{\mathrm{cr}}={10}^{-17}\exp [-{\rm{\Sigma }}/(96\,{\rm{g}}\,{\mathrm{cm}}^{-2})]$ s⁻¹ (Umebayashi & Nakano 1981), though in two of our simulations, we set the prefactor of this quantity to 10⁻¹⁹ s⁻¹ (see the description below). We also include X-rays as in S15 with an X-ray luminosity L_X = 4 × 10²⁹ erg s⁻¹ and X-ray temperature corresponding to 1 keV, specifically appropriate to HD 163296 (Swartz et al. 2005; Günther & Schmitt 2009). Additionally, we include ²⁶Al decay at an ionization rate 10⁻¹⁹ s⁻¹.

In the surface layer, we further include the effect of far-UV ionization based on calculations from Perez-Becker & Chiang (2011). The FUV ionization is assumed to have a constant penetration depth of Σ_FUV = 0.01 g cm⁻². In these surface layers, we calculate the ionization rate by setting the ambipolar diffusion Elsässer number Am = γρ_i/Ω, where γ is the coefficient of momentum transfer for ion-neutral collisions, ρ_i is the mass density of ions, and Ω is the orbital frequency, and then calculating the diffusivity from Am. Compared with S15, we use a modified version of Am in the FUV-ionized layer more appropriate to the structure of HD 163296,

$\begin{eqnarray}&&{\mathrm{Am}}_{\mathrm{FUV}}=2.4\times {10}^{16}{\rho }_{0,\mathrm{mid}}\left(\displaystyle \frac{f}{{10}^{-5}}\right){\left(\displaystyle \frac{r}{1\mathrm{au}}\right)}^{3/2},\end{eqnarray} \tag{ 1 }$

where ρ_0,mid is the midplane gas density in cgs units, f is the ionization fraction, and r is the radial coordinate away from the star. In calculating the diffusivities, we assume no grains are present since S15 demonstrated that the inclusion of 0.1 μm grains made very little difference to the amplitude of turbulence. We further discuss the neglect of dust grains in Section 5.

The initial conditions are identical to the simulations of S15. The standard domain size is 4H × 8H × 8H in (x, y, z), where the vertical scale height $H=\sqrt{2}{c}_{{\rm{s}}}/{\rm{\Omega }}$ , with c_s and Ω being the sound speed and orbital frequency, respectively. We decrease the vertical extent of the box for some of our simulations in order to reduce the overhead associated with small time steps set by low-density gas. We employ a resolution of 36 grid zones per H in all calculations. All simulations assume an isothermal equation of state. We parameterize the net vertical flux threading the simulation volume in terms of ${\beta }_{0}\,\equiv \mathrm{sign}({B}_{z,0})2{P}_{0}/{B}_{z,0}^{2}$ , where P₀ is the initial midplane gas pressure and B_z,0 is the strength of the background vertical field. Note that due to the dependence of the Hall effect on the sign of this vertical field (e.g., Kunz & Lesur 2013; Lesur et al. 2014; Bai 2015; Simon et al. 2015b), β₀ includes this sign dependence. We run each calculation out so that a statistically steady state is reached.

In what follows, we vary β₀, the radial location of the shearing box R₀, and the strength of various ionization sources. These simulations are summarized in Table 1. The label of each simulation is R#-B#(p, n)-FXC, where the first # is the radial location of the shearing box in units of au, the second # is the exponent of β₀ (e.g., 3 for β₀ = 10³), (p, n) refers to a positive or negative value for β₀, and the inclusion of "F," "X," "C" denotes the inclusion of FUV, X-rays, and cosmic rays, respectively (in our simulations, ²⁶Al decay is always present as an additional ionization source). As mentioned above, some of the simulations were run in a shorter domain in order to ensure that the time step (which is controlled by the lowest density gas near the vertical boundaries) remains reasonable. For these calculations, we append the label with #H, where # is the total number of vertical scale heights spanned by the domain. Note that two additional runs were done with cosmic rays and ²⁶Al decay only, but with the cosmic-ray ionization rate set to 10⁻¹⁹ s⁻¹; these runs are labeled with "CWeak."

Table 1. Numerical Simulations

Label	R₀	β₀	FUV?	X-Ray?	CR?	${\left\langle \tfrac{\delta v}{{c}_{{\rm{s}}}}\right\rangle }_{{ \mathcal G },z\geqslant 2H}$	α
	(au)
R100-B3p-FXC-6H	100	10³	Yes	Yes	Yes	0.39	2.7 × 10⁻²
R100-B4p-FXC	100	10⁴	Yes	Yes	Yes	0.25	9.7 × 10⁻³
R100-B4n-FXC	100	−10⁴	Yes	Yes	Yes	0.45	7.4 × 10⁻³
R100-B5p-FXC	100	10⁵	Yes	Yes	Yes	0.73	2.7 × 10⁻³
R100-B6p-FXC	100	10⁶	Yes	Yes	Yes	0.61	8.4 × 10⁻⁴
R100-B7p-FXC	100	10⁷	Yes	Yes	Yes	0.38	2.3 × 10⁻⁴
R100-B3p-XC-5.25H	100	10³	No	Yes	Yes	0.21	1.1 × 10⁻²
R100-B4p-XC-6H	100	10⁴	No	Yes	Yes	0.30	2.9 × 10⁻³
R100-B4n-XC-6H	100	−10⁴	No	Yes	Yes	0.13	2.4 × 10⁻⁴
R100-B5p-XC-6H	100	10⁵	No	Yes	Yes	0.14	7.0 × 10⁻⁴
R100-B6p-XC-6H	100	10⁶	No	Yes	Yes	0.09	2.3 × 10⁻⁴
R100-B7p-XC-6.5H	100	10⁷	No	Yes	Yes	0.04	5.9 × 10⁻⁵
R100-B4p-X-5.5H	100	10⁴	No	Yes	No	0.29	1.7 × 10⁻³
R100-B4p-C-5.5H	100	10⁴	No	No	Yes	0.20	2.3 × 10⁻²
R100-B4p-CWeak-5.25H	100	10⁴	No	No	Very weak	0.08	1.8 × 10⁻³
R100-B5p-CWeak-5.25H	100	10⁴	No	No	Very weak	0.02	5.7 × 10⁻⁴
R200-B4p-XC-5.25H	200	10⁴	No	Yes	Yes	0.18	1.4 × 10⁻³

Download table as: ASCII Typeset image

3. Turbulent Velocities

Following the arguments in S15, we calculate the turbulent velocity, δv, by subtracting the shear flow, the laminar wind component, and the systematic increase or decrease in azimuthal velocity due to zonal flows and then adding the remaining velocity components in quadrature. In what follows, we time average the turbulent velocity over a period in which various flow quantities (e.g., Maxwell stress) are statistically constant in time. Figure 1 shows this velocity, normalized to ${c}_{{\rm{s}}}$ , versus z for a range of background magnetic field strengths, ranging from β₀ = 10⁷ to β₀ = 10³. The complete run list, along with a geometrically averaged (denoted by ${ \mathcal G }$ ) turbulent velocity (for $| z| \gt 2H$ ) is shown in Table 1 and is demonstrated pictorially in Figure 2. We also include the standard Shakura & Syunyaev (1973) α value in the table.

**Figure 1.** Horizontally and temporally averaged $\delta v/{c}_{{\rm{s}}}$ vs. z for magnetic field strengths characterized by β₀ = 10³ (red), 10⁴ (green), 10⁵ (blue), and 10⁷ (black), with all sources of ionization (solid lines), X-rays and cosmic rays only (dashed), and only cosmic rays at a reduced flux (dotted–dashed). The horizontal dotted line represents the approximate upper limit on the turbulent velocity $\delta v\lt 0.05{c}_{{\rm{s}}}$ from F15 and F17. The simulations that fall below this line (or are somewhat marginal) are for weak magnetic fields (β₀ > 10⁵) and reduced ionization.
Download figure:
Standard image High-resolution image

**Figure 1.** Horizontally and temporally averaged $\delta v/{c}_{{\rm{s}}}$ vs. z for magnetic field strengths characterized by β₀ = 10³ (red), 10⁴ (green), 10⁵ (blue), and 10⁷ (black), with all sources of ionization (solid lines), X-rays and cosmic rays only (dashed), and only cosmic rays at a reduced flux (dotted–dashed). The horizontal dotted line represents the approximate upper limit on the turbulent velocity $\delta v\lt 0.05{c}_{{\rm{s}}}$ from F15 and F17. The simulations that fall below this line (or are somewhat marginal) are for weak magnetic fields (β₀ > 10⁵) and reduced ionization.
Download figure:
Standard image High-resolution image

**Figure 2.** Spatial (geometric) and time average of $\delta v/{c}_{{\rm{s}}}$ for $| z| \geqslant 2H$ as a function of β₀ and ionization level. Numbers in blue are well above the upper limit of $0.05{c}_{{\rm{s}}}$ (from F17), whereas red numbers are below this limit. The purple numbers correspond to $\delta v\lt 0.1{c}_{{\rm{s}}}$ and may be considered marginal. The dashed lines suggest a region of parameter space in which a weaker magnetic field can be traded off for stronger ionization in order to maintain low turbulence values.
Download figure:
Standard image High-resolution image

It is clear from both the figures and the table, that for a large range in magnetic field strengths, β₀ = 10³–10⁷, there are significant turbulent velocities at large $| z|$ , well above the upper limit put on this velocity by F15 and F17, represented approximately by the dotted line in Figure 1. As β₀ decreases, the magnetic field structure tends toward being more laminar as shown by Simon et al. (2013a) and verified here. Thus, for different magnetic field configurations (largely laminar, largely turbulent, or some combination thereof), we still find velocities above observational constraints at large $| z|$ . For fields weaker than β₀ ∼ 10³–10⁴, the turbulence in the midplane is well below this upper limit and thus consistent with these observations. Even the simulations in which only one source of ionization (in addition to ²⁶Al decay) is included at the strength used previously (S15) show appreciable turbulence at large $| z|$ for all but the weakest vertical fields.

We have varied other parameters, including magnetic field polarity, R₀, and the strength of ionization sources. In nearly all cases, we find strong turbulence with δv ∼ 0.1–1 ${c}_{{\rm{s}}}$ at $| z| \sim 2\mbox{--}4H$ . Based on the data in Table 1 and Figure 2, there is a trade-off between field strength and ionization level; for higher ionization, a weaker background field is necessary to produce low turbulence values.

Note that for β₀ ≥ 10⁶, the MRI is largely under-resolved near the disk midplane. However, since the strength of MRI turbulence decreases with increasing β₀ (Bai & Stone 2011), the trend of the midplane δv values with β₀ suggests that δv will be quite small there, and will not likely strongly affect the turbulence in the upper layers.

Furthermore, due to a restrictively small time step, we could not run stratified simulations with β₀ < 10³ and Am ≪ 1. However, with a strong net field (e.g., β₀ ≤ 10²), one can conceive of two regimes. First, for Am sufficiently low to suppress/damp the MRI but still sufficiently large for there to be some coupling between gas and magnetic field, a strong field would drive a substantial wind, resulting in accretion rates well in excess of the observed rate (e.g., Simon et al. 2013a). Second, for extremely weak ionization (Am → 0), the field and gas are completely decoupled, and hence there is no accretion. An intermediate regime in which the coupling and field strength are at just the right level to give the right accretion rate may exist, but as we discuss in Section 4.3, even here, one expects to see turbulent gas velocities.

4. Origin of the Turbulence

Our results clearly show that even in the presence of a predominantly wind-driven accretion flow (i.e., the magnetic field structure appears relatively laminar), there are still significant turbulent velocities. This result is surprising in light of recent results demonstrating that for a sufficiently strong magnetic field, ambipolar diffusion quenches the turbulence (see Figure 16 in Bai & Stone 2011). To further understand these results, we now focus on determining the origin of this turbulence.

4.1. Local MRI Turbulence

We have examined the spatial structure of the magnetic field and ${\beta }_{z}\equiv 2P/{B}_{z}^{2}$ , where β_z is defined locally and not as the background value. Figure 3 shows β_z for a snapshot at orbit 30 in R100-B3p-XC-5.25H, for which β₀ = 10³. It appears that the magnetic flux has been redistributed into nearly axisymmetric structures. With this redistribution of flux, large regions of the box have β ≥ 10⁴, sufficiently large for the MRI to be active in those regions, based on Figure 16 of Bai & Stone (2011).⁸ That these regions are MRI active is supported by the right panel of the figure, in which a meridional slice of the domain is shown with magnetic field lines overplotted. In the regions of low β_z, the field lines are laminar, whereas they become quite tangled when β_z ≳ 10⁴.

**Figure 3.** Spatial distribution of the local ${\beta }_{z}\equiv 2P/{B}_{z}^{2}$ for a snapshot at orbit 30 in R100-B3p-XC-5.25H. The left panel shows slices along the midplane and y = 0 in the simulation domain. The right panel shows the y = 0 slice with magnetic field lines superimposed (red lines). There are large-scale, mostly axisymmetric, spatial variations in β_z. In large β_z regions, the magnetic field is tangled due to turbulence, whereas for β < 10³, the field is relatively laminar.
Download figure:
Standard image High-resolution image

Examining our remaining calculations, we find that for 10³ ≤ β₀ ≤ 10⁴ (with the exception of the run R100-B4p-CWeak-5.25H; see below) zonal flows are established near the midplane, which lead to regions of magnetic flux that are sufficiently weak to drive the MRI. In all of these cases, while the turbulence apparently originates near the midplane, velocity fluctuations are likely amplified as they propagate toward the lower density regions at large $| z|$ (Simon et al. 2011).

For simulations with 10⁵ ≤ β₀ ≤ 10⁷ (with the exception of the run R100-B5p-CWeak-5.25H; see below), the background magnetic field is sufficiently weak that there is MRI-driven turbulence throughout a large portion of the domain (again, based on the condition that weaker fields lead to turbulence; Bai & Stone 2011). It is worth pointing out again that for runs with β₀ ≥ 10⁵, the MRI is largely under-resolved near the midplane. Thus, higher resolution simulations will be required to determine the precise amplitude and structure of the turbulence in these cases. However, as we discussed above, for sufficiently weak magnetic field, the turbulent velocities induced by the MRI will likely be small enough that they will not appreciably change our results.

Beyond the MRI, other instabilities that could potentially generate the turbulence seen in our simulations are the Hall-shear instability (HSI; Kunz 2008) and the ambipolar diffusion shear instability (ADSI; Kunz 2008). We have run one additional simulation in which we neglect the Hall effect and set Am = 1 everywhere. We still see significant turbulence with a similar structure to that seen in the runs with all three nonideal terms. This result suggests that the HSI, while potentially present in our simulations, does not play the dominant role in driving turbulence in these simulations.

We cannot say for certain that the ADSI is absent in our simulations, as to our knowledge, no well-tested numerical study of the nonlinear ADSI has been done. However, as suggested by the work of Pandey & Wardle (2012), this instability can be considered as an extension of the MRI in the strong ambipolar diffusion dominated regime.

Taken together, these considerations strongly suggest that the MRI is largely responsible for generating turbulence in our simulations, even in the case of a relatively strong magnetic field, and that this turbulence originates from regions of relatively small vertical magnetic flux. Furthermore, in some of our simulations, a current sheet persists near the disk midplane, which is the result of large-scale toroidal fields of opposite polarities coming into contact. These current sheets generate velocity fluctuations, which likely become amplified as they move toward lower density regions (see Simon et al. 2011). In summary, even so-called laminar wind flows are not entirely laminar.

4.2. Zonal Flows

We have analyzed the axisymmetric structures seen in R100-B3p-XC-5.25H in order to better characterize and understand their properties. Figure 4 shows the spacetime diagram of the vertical magnetic energy averaged over all y and for $| z| \lt 0.5H$ . As the figure shows, after ∼15 orbits, very narrow bundles of vertical magnetic flux emerge and persist until the end of the simulation. The radial width of these structures (as measured at the base; see Figure 5) is approximately ∼0.5H–H, consistent with the zonal flows seen in Bai & Stone (2014), and especially Bai (2015), where the very narrow width of the magnetic flux bundles here roughly match that seen in that work. Bai & Stone (2014) explained the emergence of these zonal flows with a phenomenological model. Given the differences between the simulations in that work (e.g., unstratified, no Hall effect) and ours, it is beyond the scope of our work to match precisely the properties of the zonal flows we see here compared to these previous works. However, some generic properties should persist if these structures are indeed zonal flows, which we now test.

**Figure 4.** Spacetime diagram (*t, x*) of the vertical magnetic energy normalized to the initial midplane gas pressure for the simulation R100-B3p-XC-5.25H. The average of the magnetic energy was taken over all y and over $| z| \lt 0.5H$ . The emergence of very narrow bands of large magnetic flux emerge after ∼15 orbits. These structures have been seen in the simulations of Bai (2015) and coincide with the strong magnetic flux regions shown in Figure 3.
Download figure:
Standard image High-resolution image

Figure 5 plots the relative perturbation amplitude, ${\rm{\Delta }}Q/\langle Q\rangle \,=(Q-\langle Q\rangle )/\langle Q\rangle$ for some quantity Q that has already been averaged over all y, over $| z| \lt 0.5H$ , and in time from orbits 20 to 35, during which the radial movement of the magnetic flux bundles is small (see Figure 4). The angled brackets denote an average over x. Here, we plot this relative perturbation for the vertical magnetic energy, Maxwell stress, and the gas pressure. We also plot the perturbed (i.e., Keplerian-shear-subtracted) azimuthal velocity normalized by the sound speed $\delta {v}_{y}/{c}_{{\rm{s}}}$ , time-averaged over the same period of time as the other quantities. Note that we had to reduce the vertical magnetic energy and Maxwell stress by a factor of 10 for it to fit on the same plot. Consistent with previous zonal flow models (e.g., Johansen et al. 2009; Bai & Stone 2014), the magnetic flux and Maxwell stress are out of phase with the gas pressure by ∼π, and the azimuthal velocity is out of phase with the gas pressure by ∼π/2. Furthermore, in geostrophic balance, the gas pressure gradient balances with the Coriolis force,

$\begin{eqnarray}&&\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{2{\rm{\Omega }}}\displaystyle \frac{\partial \overline{\rho }}{\partial x}\approx {\rho }_{0}\overline{\delta {v}_{y}},\end{eqnarray} \tag{ 2 }$

where the over-bar represents the time and y, z average that is done in Figure 5 and ρ₀ is the midplane gas density. We have checked this balance by calculating both sides of Equation (2); we find that as a function of x both sides agree with each other to within a factor of order unity. Thus, the features that we observe in this simulation are indeed zonal flows.

Comparing the location of the magnetic flux bundles with that seen in the snapshot at 30 orbits (Figure 3), we see that they match exactly. Clearly, the turbulence seen in our simulations arises from the weakly magnetized, high-density regions of the zonal flows. This result is consistent with the fact that the MRI requires a weak magnetic field to become active, especially in the presence of strong ambipolar diffusion.

4.3. Caveat: The "CWeak" Simulations

Two of our simulations, R100-B4p-CWeak-5.25H and R100-B5p-CWeak-5.25H, which assume cosmic-ray ionization at a reduced rate (in addition to ²⁶Al decay), depict a very laminar magnetic field structure. Yet, there are clear indications that the gas flow is turbulent (see Figure 2). We have analyzed the velocity structure in these simulations and have discovered the presence of an oscillatory, vertically extended circulation pattern, nearly identical to that seen in Gole et al. (2016). While these flows are not turbulent per se, their radial structure is significantly smaller than what can be resolved with ALMA. Thus, these flows would manifest themselves as turbulent broadening in the observations (Flaherty et al. 2018, in preparation).

While the velocity fluctuations seen in these simulations are weak and below (or marginally above) the observational constraints (see Figure 2), there is a trend of increasing velocities with magnetic field strength. Thus, we anticipate that for stronger magnetic fields (β₀ ≤ 10³), the velocity fluctuations will increase.

Upon further examining the magnetic field structure in these simulations, we find that radial and toroidal fields are sitting at the midplane in a steady-state configuration. For such a steady-state configuration to exist, the continued shearing of radial field into toroidal field must be balanced by buoyancy and/or diffusion to vertically remove magnetic flux from the midplane region. Such departures from a stationary configuration will likely stir the gas to some extent, inducing fluid motions. In this particular case, the fluid motions induced are a natural oscillatory mode of the disk (Lubow & Pringle 1993; Gole et al. 2016). That we see stronger gas motions for stronger magnetic fields is consistent with the notion of the magnetic field stirring the gas. More work is required to fully understand the generation of these modes from the magnetic fields. However, our preliminary studies suggest that this is a physical effect, as is the increase in velocity fluctuations with magnetic field strength.

5. Discussion

In the following sections, we discuss the implications of our results for protoplanetary disk structure and evolution. We first put forth a preliminary model to explain current observations within the context of our simulations, followed by a discussion of support for the model, limitations, and predictions.

5.1. A Preliminary Model for Disk Evolution and Structure

While we have focused here on the HD163926 disk, for which there is significant evidence for low turbulence, there are other systems that display similarly low values (Flaherty et al. 2018, K. M. Flaherty et al. 2018, in preparation). These results, coupled with the numerical simulations presented here, call into question the notion of an outer disk that is strongly accreting, either through magnetically driven winds or turbulence. Of course, there is ample evidence that gas from the inner disk is accreting onto the star in many of these systems (for HD163926, see Mendigutía et al. 2013).

Here, we put forth a preliminary model, described pictorially in Figure 6, to explain current observations that show weak turbulence in the outer disk, but with continued accretion from the inner disk. A key component in this picture is the presence of a large-scale vertical field in the inner regions of the disk (≲30 au). This vertical field launches a wind (through a Blandford–Payne type process Blandford & Payne 1982) that blocks FUV photons (and possibly X-rays; see the discussion below) from reaching the outer disk. Furthermore, this magnetic field may also reduce the cosmic-ray flux in the outer disk through mechanisms such as magnetic mirroring and funneling (Cleeves et al. 2013, L. I. Cleeves 2018, private communication). A second component to this picture is that any large-scale vertical magnetic field present in the outer disk must be relatively weak (though, see the discussion below), having been advected inward, or (as supported by some recent global calculations; see Bai & Stone 2017) diffused to very large distances. In either case, most of the magnetic flux remains in the inner disk with less flux present on the scales of ∼100 au. This model thus satisfies the two conditions we have found necessary for significantly reduced turbulence: low ionization and weak vertical magnetic field in the outer disk.

5.2. Support for the Model, Limitations, and Caveats

There is both theoretical and observational support for a large-scale field in the inner disk. Recent global simulations of the inner disk region (Bai 2017) have found no evidence for outward transport of magnetic flux from the inner disk throughout the duration of the simulations, consistent with the maintenance of a strong magnetic flux in the inner disk once it has been developed. Furthermore, the presence of an inner disk field that launches a wind is consistent with observational constraints found in optical forbidden line studies (Simon et al. 2016), and by modeling the near-infrared emission properties of Herbig Ae stars (Bans & Königl 2012). Specific systems for which there is some evidence of a wind (though it is possible these winds are photoevaporative and not magnetic in origin) include V4046Sgr (Sacco et al. 2012), MWC480 (Fernandes et al. 2018), TW Hya (Pascucci et al. 2011), and HD 163296 (Klaassen et al. 2013). That the two systems currently shown to exhibit weak turbulence also show evidence of a wind, providing strong support for our model.

The degree to which ionizing radiation is prevented from reaching the outer disk is more uncertain, however. In particular, a substantial wind would be required to block X-rays altogether, and indeed, there is observational work (e.g., Pascucci et al. 2011, 2014) suggesting that X-rays can penetrate beyond 1 au in protoplanetary disks. It is worth mentioning, however, that while the absorption column for X-rays is quite small (of the order of 0.01 g cm⁻²; Igea & Glassgold 1999), and thus may be blocked by the wind, the scattering component can ionize gas to significantly deeper columns (∼1–10 g cm⁻² in a protoplanetary disk; Bai & Goodman 2009). How precisely this scattering component would penetrate through the wind and impinge upon the disk is not well known, and it may be the case that X-ray flux on the outer disk is significantly reduced as a result of scattering by the wind. Finally, for cosmic rays to be blocked, the magnetic field would likely need to be strongly inclined with respect to the disk rotation axis; if instead the field has a small inclination or is collimated, cosmic rays could reach the outer disk unimpeded.

However, we reiterate that even if no X-rays and cosmic rays are blocked, one can still achieve turbulent velocities consistent with observations so long as no FUV flux reaches the outer disk and a sufficiently weak magnetic field is present there (Figure 2). When the column required to block FUV is small (e.g., 0.01 g cm⁻² Perez-Becker & Chiang 2011), a modestly strong wind can prevent these photons from reaching the outer disk, as found in some of the global simulations of Bai (2017). There is observational support for a wind blocking FUV photons as well. For example, as argued in Bans & Königl (2012), a dusty FUV-shielding wind could explain the strong near-infrared excess seen in accreting systems around Herbig Ae stars. Indeed, given the large UV flux from the HD 163296 star (Meeus et al. 2012), it is not inconceivable that a substantial dust-filled wind will be required to block UV photons from reaching the outer disk.

While this preliminary picture is consistent with our observational and theoretical constraints, there remain several uncertainties with this picture and the simulations that inform it. First, as mentioned above, we have neglected the effect of dust on the ionization fraction. While S15 demonstrated that the inclusion of 0.1 μm grains with mass ratio 10⁻⁴ made very little difference to the amplitude of turbulence in the outer disk compared with the grain-free case,⁹ higher abundance of submicron grains could lower the ionization fraction and help suppress turbulence further. Indeed, such grains may very well be present in protoplanetary disks, at least in the upper layers (Bouwman et al. 2000, 2001; Sicilia-Aguilar et al. 2007; Zsom et al. 2011), and if they are sufficiently abundant, the ionization fraction could be driven to such a low value that no turbulence would ensue even in the presence of magnetic fields. However, the abundance levels of these small grains throughout the vertical extent of the disk are not well-constrained by observations. Additionally, calculations by Birnstiel et al. (2011) have demonstrated that in coagulation-fragmentation equilibrium, the number density of grains drops drastically below ∼0.1–1 μm size. These arguments suggest that these small grains may not play a large role in suppressing turbulence. However, the complete exclusion of very small grains remains an assumption that should be explicitly tested; this will be done in a future publication.

Furthermore, as with any study involving local simulations, there remains the question of whether the same behavior would appear in global simulations. Recent works by Bai (2017) and Béthune et al. (2017) claim their global simulations display mostly laminar fluid motions, which seem to contradict our finding of turbulent motions in local calculations. However, the vast majority of these global simulations were carried out in 2D (or focus on the inner disk, as in Bai 2017), and it is well-known that allowing for nonaxisymmetric behavior in fully 3D simulations can lead to fundamentally different behavior (Goodman & Xu 1994; Balbus & Hawley 1998; Simon et al. 2015b) compared to 2D. In addition, the one instance of a fully 3D global calculation focused on the outer regions of protoplanetary disks (see Béthune et al. 2017) was carried out at significantly lower effective resolution compared to our local simulations; in these global calculations, the number of grid zones per disk scale height was a factor of 2 (8) lower in radius and height (azimuth) compared to our local models, and the highly diffusive HLLE flux solver (e.g., O'Neill et al. 2012; Salvesen et al. 2014) was used instead of the HLLD solver used in the simulations presented here. Furthermore, the origin of turbulence, as described in Section 4, is a highly local phenomenon, occurring on scales ∼H as the result of zonal flows. Indeed, zonal flows of comparable scale to those seen here have been seen in global simulations themselves (Béthune et al. 2017). It seems likely that in full 3D and with sufficiently high resolution, patches of turbulence will be present in global simulations as well.

The presence of the oscillatory modes as seen in the "CWeak" simulations (see Section 4.3), while a physical effect, should also be verified within the context of global simulations. As pointed out by Gole et al. (2016), such modes would likely have somewhat different properties in global disks compared with local models. If these modes exhibited substantially smaller velocity amplitudes in global simulations, it is possible that the gas velocities in the "CWeak" simulations would fall below the observational limit regardless of background field strength, which may temper the claim that any large-scale magnetic field present in the outer disk should be weak. However, we suspect that the properties of the modes, particularly their amplitude, would not change significantly in global simulations since the physical mechanism driving the modes would be the same. In summary, while the numerical results presented here are encouraging, global simulations will be an important next step in testing this preliminary model, and we will pursue such calculations in future studies.

Finally, it is worth commenting on the applicability of our HD 163296-based simulations to TW Hya, the only other system currently shown to display weak turbulence (Flaherty et al. 2018). Compared to HD 163296, the TW Hya disk is lower in mass (Isella et al. 2007; Bergin et al. 2013; F15; Flaherty et al. 2018) and thus has lower gas densities at R ∼ 100 au, the scales most relevant to ALMA observations (F15; Flaherty et al. 2018). Since ambipolar diffusion is even more dominant, compared with other nonideal MHD effects, at lower densities (Kunz & Balbus 2004), we expect that the results we present here, which are extracted from ambipolar diffusion dominated simulations, are generally applicable to TW Hya as well as HD163926. In addition, since TW Hya is nearly face-on, the dominant turbulent velocity to be detected by observations would be the z component. We have examined the turbulent velocity components, namely xy versus z turbulent velocities, in all of our simulations. The vertical velocity component dominates over (by factors ranging from 2–6) or is very comparable to the xy velocities at heights ranging from 2H to 4H above the midplane. These results suggest that if present, turbulence would be observable in TW Hya as well.

5.3. Model Predictions

Despite the uncertainties discussed in the previous section, our model makes several powerful and testable predictions for disks that show weak turbulence. First, the outer disk (at scales of ∼100 au) should have weak ionization;¹⁰ at the very least, FUV ionization should be absent, and the degree of remaining ionization depends on how much magnetic field exists at these distances (see below). There is already evidence for weak ionization in the ratios of neutral and ionic molecules in the TW Hya disk (Cleeves et al. 2015, but see Mathews et al. 2013), and while recent work (Teague et al. 2016) has put large upper limits on turbulent velocities δv < 0.2– $0.4{c}_{{\rm{s}}}$ , a new analysis (Flaherty et al. 2018) has suggested that indeed turbulence in this system is quite weak.

Second, there should be a magnetically launched wind originating from small radii, sufficiently strong so as to largely reduce the FUV (and possibly X-ray) flux at the outer disk. As mentioned above, there are a number of systems that show some evidence of a wind.

Finally, any magnetic field present in the outer disk, and thus any wind launched from this region, must be weak. Conversely, if future observations find a disk with turbulence present, then it will be equally useful to search for high levels of ionization and the presence of magnetic field (e.g., through Zeeman observations) in the outer disk. Furthermore, these systems should show no sign of a wind launched from the inner disk capable of blocking ionizing radiation. Intriguingly, there is tentative evidence that the system DM Tau is lacking such a wind (Simon et al. 2016), which is consistent with preliminary work showing the presence of strong turbulence in this system (K. M. Flaherty et al. 2018, in preparation). These results, if verified, will strongly support our model.

Because of the highly nonuniform distribution of magnetic flux in our model, one further prediction is that accretion is likely very nonsteady. In the case of HD 163296, the (instantaneous) high accretion rate of ∼4.5 × 10⁻⁷ M_⊙ yr⁻¹ (Mendigutía et al. 2013) is only telling us about the accretion rate through the inner disk; the outer disk can maintain a low accretion rate, consistent with our model. More support for a nonsteady accretion rate can be acquired via a simple timescale argument. Given the disk mass, 0.09 M_⊙ (Isella et al. 2007), the instantaneous accretion rate would yield a lifetime of only ∼200,000 years, which is significantly shorter than the actual lifetime (∼3 Myr; Montesinos et al. 2009). However, the accretion rate of HD 163296 was found to be an order of magnitude lower ∼20 years ago (Mendigutía et al. 2013), which if more representative of the average accretion rate over the disk's lifetime, would be roughly consistent with the system's age, supporting the notion of a nonsteady disk. Even at this lower accretion rate, however, a significant amount of magnetic flux in the inner disk would be required to drive accretion. Ultimately, it appears inevitable that in our model, the inner disk will be depleted on a relatively short timescale. Not only does this predict a nonsteady accretion rate, but it may also be related to the formation of transition disks, as some models suggest (Wang & Goodman 2017).

Given the trade-off between field strength and ionization level, as shown in Figure 2, it is difficult to quantify the precise values of the field strength and ionization level required to generate weak turbulence consistent with observations. However, assuming that β₀ = 10⁷ in the outer disk, then to be consistent with observations, there should be no FUV flux in the outer disk (though X-rays and cosmic rays can still be present). For the HD 163296 disk, this field strength corresponds to ∼5–10 μG at 100 au, roughly consistent with estimates of Milky Way magnetic field strengths (Haverkorn 2015 and references therein). Further support for such a weak field comes from comparing the ionization levels in our model and the abundance of DCO+ from F17. The abundance of DCO+ puts a lower limit on the ionization fraction at the disk midplane of ∼10⁻¹¹. This is consistent with the ionization level in our simulations that include X-ray and cosmic-ray ionization (at the standard ionization rate of 10⁻¹⁷ s⁻¹). However, for the "weak CR" simulations, the ionization fraction is an order of magnitude lower than this, suggesting that in HD 163296, X-rays and cosmic rays are not blocked, whereas FUV photons are, and there is a weak vertical magnetic field at large distances from the star.

Of course, if the overall global magnetic field is stronger (β₀ < 10⁷), then one might expect that this stronger magnetic field would allow less ionizing flux from reaching the outer disk, possibly blocking X-ray photons in addition to FUV. This would be consistent with the simulation at β₀ = 10⁵ with both FUV and X-ray ionization removed. While potentially not applicable to HD 163296, as discussed above, this picture may apply to other disk systems. Ultimately, however, more sophisticated fully three-dimensional global calculations will be required to determine the magnetic field strength needed to block ionizing sources of radiation.

6. Conclusions

We have carried out a series of local, shearing-box simulations in order to further quantity the influence of magnetic fields on protoplanetary disks in the presence of low ionization. Our primary results are as follows:

1.
For a relatively strong background vertical field (β₀ = 10³), which drives a largely laminar wind, turbulence is still present at values that are inconsistent with observational constraints and arises from regions of zonal flows where the vertical magnetic field is sufficiently weak for the MRI to persist. Magnetic winds are not the solution to the weak turbulence issue.
2.
Only with a weak vertical magnetic field and a low ionization fraction are turbulent velocities consistent with observations. Furthermore, a trade-off exists; if the magnetic field strength (ionization fraction) is increased, then the ionization fraction (magnetic field strength) must be decreased to match observations.
3.
A preliminary picture to explain these observational and theoretical constraints is that a large-scale vertical magnetic field must be present in the inner disk to launch a strong wind that shields the outer disk from ionizing radiation, and the outer disk itself has a very weak (or no) large-scale field.
4.
Such a nonuniform distribution of magnetic flux implies that such disks are nonsteady, which, in the case of HD 163296, is supported by observational measurements of the accretion rate onto the central star.

While the details of this preliminary model need to be tested with fully 3D global MHD simulations, the model predicts that if indeed there is low turbulence inferred in the outer disk, then observations should constrain the ionization rate in the outer disk to be relatively small and should demonstrate the presence of a large-scale magnetically launched wind in the inner disk. With ALMA operations continuing well into the coming decade, we should have many opportunities to test these predictions.

We thank Phil Armitage, Ilse Cleeves, Ilaria Pascucci, Richard Teague, and Upasana Das for useful discussions related to this work. We also thank the anonymous referee whose comments improved the quality of this paper. J.B.S. acknowledges support from NASA grants NNX13AI58G and NNX16AB42G, and from NSF grant AST 1313021. The computations were performed on Stampede, Stampede 2, and Maverick through XSEDE grants TG-AST120062 and TG-AST140001. We thank the Kavli Institute for Theoretical Physics, supported in part by the National Science Foundation under Grant No. NSF PHY-1125915, for hospitality during the completion of the paper.

Software: Athena (Stone et al. 2008).

Origin of Weak Turbulence in the Outer Regions of Protoplanetary Disks

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Method

3. Turbulent Velocities