GLOBAL SIMULATIONS OF PROTOPLANETARY DISKS WITH OHMIC RESISTIVITY AND AMBIPOLAR DIFFUSION

The Ohmic and ambipolar diffusivities are evaluated dynamically for each grid-cell and are based on a precomputed look-up table, that has been derived self-consistently from a detailed chemical model accounting for grain-charging and gas-phase chemistry. The ionization rates that enter this equilibrium chemistry are governed by ionizing radiation entering the disk. The column densities that govern the attenuation of the external ionizing radiation are also computed on-the-fly, that is, they reflect the instantaneous state of the disc. In principle, this allows for a self-limiting of the emerging wind via shielding of ionizing radiation (Bans & Königl 2012).

In addition to some minor modifications to the grain charging prescription (as detailed below), we have improved the realism of our ionization model compared to our previous work. The main alterations concern the inclusion of two additional radiation sources, an un-scattered direct X-ray component, and hard UV photons. These have in common a shallow penetration depth but comparatively high flux levels, thus mainly affecting the ionization level of the disk's surface layers.

2.1.1. FUV Ionization Layer

2.1.2. Illumination by X-rays

A similar treatment is adopted for the unscattered X-ray component, for which we use the fit to the Monte-Carlo radiative-transfer results of Igea & Glassgold (1999) given by Equation (21) in Bai & Goodman (2009). Deviating from our previous local simulations (see Equation (1) in Gressel et al. 2011), our new X-ray illumination is

$$\begin{eqnarray}\begin{array}{rcl} {{\zeta }_{{\rm XR}}} & = & {{10}^{-15}}\;{{{\rm s}}^{-1}}\left[ {\rm exp} -{{\left( \frac{{{{\Sigma }}_{A}}}{{{{\Sigma }}_{{\rm sc}}}} \right)}^{\alpha }}+{\rm exp} -{{\left( \frac{{{{\Sigma }}_{B}}}{{{{\Sigma }}_{{\rm sc}}}} \right)}^{\alpha }}\ \right]\ {{r}^{-2}} \\ {} & {} & +6\times {{10}^{-12}}\;{{{\rm s}}^{-1}}\ \ {\rm exp} \;-{{\left( \frac{{{{\Sigma }}_{C}}}{{{{\Sigma }}_{{\rm ab}}}} \right)}^{\,\beta }}\ {{r}^{-2}} \\ \end{array}\end{eqnarray} \tag{ 4 }$$

where ${{{\Sigma }}_{A}}$ (${{{\Sigma }}_{B}}$) is the gas column to the top (bottom) of the disk surface, and ${{{\Sigma }}_{C}}$ is the column density along radial rays toward the star. The radial column contains a contribution from the inner disk (which is not part of the computational domain) and reaches in to an inner truncation radius of ${{r}_{{\rm tr}}}=5{{R}_{\odot }}$. Adopting values from Bai & Goodman (2009), the shape coefficients are ${{{\Sigma }}_{{\rm sc}}}=7\times {{10}^{23}}\;{\rm c}{{{\rm m}}^{-2}}$, and $\alpha =0.65$ for scattered X-rays, and ${{{\Sigma }}_{{\rm ab}}}=1.5\times {{10}^{21}}\;{\rm c}{{{\rm m}}^{-2}}$, and $\beta =0.4$, for absorption of the direct component, respectively. Both contributions additionally decay with the squared radius to account for the decrease in X-ray luminosity. To account for typical median values of ${{10}^{30}}\;{\rm erg}\;{{{\rm s}}^{-1}}$ in observational data of luminosities in young star clusters such as in the Orion Nebula (Garmire et al. 2000), we enhance the incident X-ray flux by a factor of $5\times$ compared to the above stated values. In view of future work, we envisage additional improvements to the X-ray ionization model by adopting the new results of Ercolano & Glassgold (2013), who consider the enhanced ionizing effects of X-rays due to the presence of heavier elements (assuming solar abundance). These alterations have been found to lead to enhanced ionizing fluxes at intermediate column densities compared to the original results of Igea & Glassgold (1999).

2.1.3. Modifications to the Disk Chemistry

Our chemical network is based on the one used by Ilgner & Nelson (2006), labeled model4, with grain surface reactions and a simplified gas-phase reaction set involving one representative metal and one molecular ion. Small changes in recent years include correcting the electron sticking probabilities for the grain charge (Bai 2011) and consistently treating the molecular ion such as HCO⁺ (with a mass of 29 protons) since this is the long-lived species in a series of several reactions. Further details on the reaction network and diffusivities can be found in Section 2.2 of Landry et al. (2013) as well as Section 4.2 of Mohanty et al. (2013).

The resulting initial ionization profile expressed in Elsasser numbers ${{{\Lambda }}_{{\rm O}/{\rm AD}}}\equiv v_{Az}^{2}\;{{({\Omega }\;{{\eta }_{{\rm O}/{\rm AD}}})}^{-1}}$ with ${{v}_{Az}}\equiv {{B}_{z}}/\sqrt{\rho }$ is plotted in Figure 1 at a radius of $1\;{\rm AU}$, together with the height-dependent plasma β parameter. The implications of this figure for our simulations are discussed in Section 3.1. The region with $|z|\gtrsim 5\;H$, where the plasma parameter becomes constant with height is caused by applying a lower limit on the gas density to avoid excessive Alfvén speeds in the disk halo. Following BS13, we chose this limit at 10⁻⁶ (in our global model, this value is with respect to the initial midplane value at any given cylindrical radius). We remark that the equilibrium density profile is artificially steep in our isothermal disk model with temperature constant on cylinders. In this respect, the density floor acts to mimic a disk with realistic radiation thermodynamics, where the irradiation heating of the upper disk layers naturally creates an increased pressure scale height and hence a much shallower density profile. Because ${{\eta }_{{\rm AD}}}$ (and hence ${{{\Lambda }}_{{\rm AD}}}$) depends on the magnetization of the plasma, the actual density profile in the disk halo has an influence on how well the gas is coupled to the magnetic field lines, which in turn affects the efficiency of the magneto-centrifugal wind launching mechanism. This caveat should be kept in mind when interpreting the mass loading of the wind, which we expect to be a function of the thermodynamics in real PPDs.

Zoom In Zoom Out Reset image size

We now describe the underlying protostellar disk model, which is largely identical to the one used in Gressel et al. (2013). The equilibrium disk model is based on a locally isothermal temperature, T, decreasing with cylindrical radius as $T(R)={{T}_{0}}\;{{(R/{{R}_{0}})}^{q}}$. For $q=-1$, such a profile leads to a constant opening angle throughout the disk as it is commonly used for global numerical simulations of PPDs, and which provides us with the reference disk model here. To study what effect disk flaring has on the absorption of direct FUV photons, we also produce a moderately flaring disk with $q=-3/4$. For this value, the resulting power-law index for the local opening angle, $h(R)\equiv H(R)/R$, is $\psi \equiv (q+1)/2=0.125$, which is in agreement with observational constraints for this value based on sub-mm observations (Andrews et al. 2009), and from multi-wavelength imaging (Pinte et al. 2008).⁴ The radial power-law exponents for the gas surface density, ${\Sigma }(R)$, are $-1/2$ for our fiducial model, and $-3/8$ for the flaring disk model, respectively.

Our temperature distribution is complemented by a power-law for the midplane density, ${{\rho }_{{\rm mid}}}(R)={{\rho }_{0}}\;{{(R/{{R}_{0}})}^{p}}$, with $p=-3/2$. The equilibrium disk solution is given by

$$\begin{eqnarray}&&\rho ({\boldsymbol{r}} )={{\rho }_{0}}{{\left( \frac{R}{{{R}_{0}}} \right)}^{p}}{\rm exp} \left( \frac{G\;{{M}_{\odot }}}{c_{{\rm s}}^{2}}\left[ \frac{1}{r}-\frac{1}{R} \right]\; \right),\ {\rm and}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\Omega }({\boldsymbol{r}} )={{{\Omega }}_{K}}(R){{\left[ (p+q){{\left( \frac{H}{R} \right)}^{2}}+(1+q)-\frac{qR}{r} \right]}^{\frac{1}{2}}},\end{eqnarray} \tag{ 6 }$$

which can be derived from the requirement of hydrostatic dynamic force balance in the vertical and radial directions, and where the Keplerian angular velocity ${{{\Omega }}_{K}}(R)=\sqrt{G\;{{M}_{\odot }}}\;{{R}^{-3/2}}$. In the above equations, the square of the isothermal sound speed results from the prescribed temperature profile as $c_{{\rm s}}^{2}=c_{{\rm s}0}^{2}\;{{(R/{{R}_{0}})}^{q}}$, with a parameter ${{c}_{{\rm s}0}}$ corresponding to a value of $H(R)\equiv {{c}_{{\rm s}}}{\Omega }_{K}^{-1}=0.05R$.

Guided by previous results of local 3D simulations that exhibit quasi-axisymmetric structure, we mainly focus on 2D-axisymmetric simulations. Our computational domain covers $r\in [0.5,5.5]\;{\rm AU}$. In the latitudinal direction, the grid covers eight pressure scale heights on each side of the disk midplane, that is, $\theta \in [\pi /2-\vartheta ,\pi /2+\vartheta ]$, with $\vartheta =8\;H/r$. In the case of the flaring disk, the coverage of scale heights varies from $8\;H$ at $r=1\;{\rm AU}$ to $6.5\;H$ at $r=5\;{\rm AU}$. The grid resolution is ${{N}_{r}}\times {{N}_{\theta }}\times {{N}_{\phi }}=512\times 192\times 64$, which means that we are able to resolve relevant MRI wave lengths. In the axisymmetric models we use ${{N}_{r}}\times {{N}_{\theta }}=1024\times 384$ cells, corresponding to $\geqslant 24$ grid points per pressure scale height, matching the resolution of the three-dimensional box simulations of BS13. For the two non-axisymmetric simulations, the azimuthal extent is restricted to $\phi \in [0,\varphi ]$, with $\varphi =\pi /2$ (a quarter wedge), or $\pi /4$ to limit the computational expense. We furthermore note that in the ideal-MHD case, where turbulence develops, the azimuthal domain size has been found to have an effect on the final saturated state (Beckwith et al. 2011; Flock et al. 2012; Sorathia et al. 2012), at least for the case of zero net flux.

Because our ionization model depends on real physical values of the gas column density, we have to choose a meaningful normalization factor ${{\rho }_{0}}$, which we fix according to a surface density of ${\Sigma }=150\;{\rm g}\;{\rm c}{{{\rm m}}^{-2}}$ at the location $R=5\;{\rm AU}$, compatible with the typical minimum-mass solar nebula (MMNS Hayashi 1981). In physical units, the disk temperature is $T=540\;{\rm K}$ at a radius of $R=1\;{\rm AU}$, and $T=108\;{\rm K}$ at $5\;{\rm AU}$, resembling typical expected conditions within the protosolar disk. While our values for ρ, T and Σ are similar to the MMSN values at $5\;{\rm AU}$, our adoption of different values for the radial power-law indices means that these values are different at $1\;{\rm AU}$ compared to those adopted by BS13.

Our model disk is initially threaded by a weak vertical magnetic field ${{B}_{z}}(R)$ corresponding to a midplane ${{\beta }_{p\,0}}\,\equiv 2\;p/{{B}^{2}}={{10}^{5}}$ independent of radius for our fiducial disk. To achieve this, the vertical magnetic flux has to vary as a power-law in radius, taking into account the radial distribution of the midplane gas pressure that itself depends on the density and temperature. In our disk model, the gas pressure, $p(R)$ decreases outward, implying that ${{B}_{z}}(R)$ falls off with radius, too. While such a configuration is generally expected when accounting for inward advection and outward diffusion of the embedded vertical flux (Okuzumi et al. 2014), our particular choice of keeping the relative field strength constant is of course arbitrary (but see Guilet & Ogilvie 2014). To preserve the solenoidal constraint to machine accuracy, we compute the poloidal field components from a suitably defined vector potential ${{A}_{\phi }}(r,\theta )$. The initial disk model is perturbed with random white-noise fluctuations in the three velocity components. The rms amplitude of the perturbation is chosen at one percent of the local sound speed. We furthermore add a white-noise perturbation in the magnetic field with an rms amplitude of a few $\mu {\rm G}$, which is on the sub-percent level compared with the net-vertical field of typically a few tens of mG.

To complete the description of our numerical setup, we need to specify BCs. For the fluid variables, we employ standard "outflow" conditions at the inner and outer radial domain edges, ${{r}_{{\rm in}}}$ and ${{r}_{{\rm out}}}$. This type is a combination of a zero-gradient condition in the case of $\;{\boldsymbol{\hat{n}}} \cdot {\boldsymbol{v}} ({{r}_{{\rm in}}},\theta )\gt 0\;$ (where ${\boldsymbol{\hat{n}}}$ denotes the outward-pointing normal vector), and reflective BCs in the opposite case, thus preventing inflow of material from outside the domain. At the upper and lower boundaries (that is, in the θ direction), we furthermore balance the ghost zone values such that a hydrostatic equilibrium is obtained. This helps to control artificial jumps of the fluid variables adjacent to the domain boundary as they are frequently encountered with finite-volume schemes.

In contrast to our earlier global simulations of layered accretion disks subject to Ohmic resistivity and containing an embedded gas-giant planet (Gressel et al. 2013), we here make a different choice for the vertical magnetic-field boundary condition. Whereas we previously applied magnetic BCs of the perfect conductor type (that is, forcing the normal field component to zero and allowing non-zero parallel field), we here use pseudo-vacuum conditions, conversely enforcing the perturbed parallel field to vanish at the surface and only allowing a perpendicular perturbed field. Note that we exclude the initial net-vertical field from the procedure such that only deviations from the background field are subject to the normal-field condition. The vertical flux threading the disk is thus preserved. Since the disk's upper layers are magnetically dominated, and the Lorentz force acts to align the flow lines with the magnetic field, letting the field lines cross the domain boundary is of course a prerequisite for launching an outflow. While we realize that forcing the radial and azimuthal components of the perturbed field to vanish at the boundary may unduly restrict the magnetic topology of the emerging wind solution, we postpone the study of less-restrictive but more cumbersome boundaries to future work.

Unlike in a radially periodic local shearing-box simulation, our global model is critically affected by the inner radial boundary condition. In a real PPD, we can expect that alkali metals will be thermally ionized in this region, leading to the development of the MRI on timescales short compared with the orbital timescale at the inner domain boundary of our model. This poses the question how to best interface the MRI-turbulent inner disk with the magnetically decoupled midplane of the outer disk that we model. It is likely that MRI-generated fields can efficiently leak into the outer disk via magnetic diffusion. In a first attempt to account for the MRI-active inner disk, we gradually taper the diffusivity coefficients to zero within the inner $0.5\;{\rm AU}$ of our disk model. Studying the inner edge of the dead zone will however require dedicated simulations (similar to the ones performed by Dzyurkevich et al. 2010) including this transition region.

The inputs to our models are very similar to those included in BS13, so it is instructive to compare the Elsasser number profiles in our model with their fiducial model as a means of understanding the similarities and differences between their results and ours. The fiducial model presented in BS13 is computed at $1\;{\rm AU}$, and has a dust-to-gas mass fraction of 10⁻⁴, so we can compare this with the dotted–dashed lines in Figure 1. The profiles shown there, and in Figure 1 of BS13 are similar, with ${{\beta }_{P}}$ decreasing below unity at disk heights $|z|\gtrsim 4.5\;H$, preventing MRI turbulence from developing at these high altitudes. In their initial condition and ours, ${{{\Lambda }}_{{\rm O}}}$ increases monotonically from the midplane upward. The two initial states also have similar profiles for ${{{\Lambda }}_{{\rm AD}}}$, displaying local maxima at intermediate disk heights ($z\sim \pm 2.5\;H$). The differences in the surface densities in the BS13 models and ours at $1\;{\rm AU}$, combined with our inclusion of a direct X-ray component in addition to the scattered component means that ${{{\Lambda }}_{{\rm AD}}}\sim 10$ at this local maximum in our model, whereas it only rises modestly above unity in BS13. We therefore expect that our model with a dust-to-gas mass fraction of 10⁻⁴ should contain narrow regions at intermediate heights that support the growth of MRI channel modes. It is noteworthy that BS13 also observed the development of the MRI at the beginning of their simulations, but report that the subsequent amplification of the field causes the AD to increase. Their disk then relaxes to completely laminar state in which angular momentum transport is driven by a magnetocentrifugal wind. The larger value of ${{{\Lambda }}_{{\rm AD}}}$ in our models may allow MRI turbulence to be sustained in these regions, or may instead lead to quenching of the MRI after it has reached a more nonlinear stage of development. In this paper, we define our fiducial model to be one in which the dust-to-gas mass fraction is 10⁻³, and the Elsasser numbers for this case are shown by the solid lines in Figure 1. The larger dust concentration reduces the gas-phase electron fraction and Elsasser numbers, and the local maximum in ${{{\Lambda }}_{{\rm AD}}}$ now peaks moderately above unity (but still attaining a larger peak value than the fiducial model in BS13).

We begin discussion of the simulation results by considering the 2D axisymmetric run O-b6, for which Ohmic resistivity is included and AD is neglected. To introduce our setup, and give the reader an impression of the general appearance of our PPD model, Figure 2 visualizes the magnetic field and poloidal flow vectors from the model. The color coding of the azimuthal field strength shows the MRI-turbulent surface layers with tangled poloidal field lines superimposed in white. In the outer disk, where the MRI has not fully developed yet, the parity of the horizontal field is odd with respect to reflection at the z = 0 line, consistent with the winding-up of vertical field, B_z, by the weak differential rotation ${{\partial }_{z}}\;{\Omega }$ of our disk model. While the inner domain boundary also shows an odd field symmetry, intermediate sections of the disk show a mixed parity, reflecting the presence of MRI modes with both even and odd vertical wave numbers. Quite remarkably, the MRI channels extend over several AU in the radial direction, although we remark that this may be overly pronounced in the axisymmetric case, where the growth of parasitic modes is likely to be artificially restricted (Pessah & Chan 2008). Because we use a net-vertical field, the linear MRI modes appear more pronounced than in the otherwise comparable 3D simulations of Dzyurkevich et al. (2010) without a net B_z field.

Zoom In Zoom Out Reset image size

The poloidal velocity field plotted as black arrows is mostly disordered in the turbulent regions but also shows some level of coherence in the upper layers, where a wind topology can be seen. While the wind is intermittent rather than steady and laminar, a general tendency for outflow is seen. This is consistent with similar observations of disk winds being launched from fully MRI-active accretion disks (see, e.g., Suzuki & Inutsuka 2009; Fromang et al. 2013). We quantify the mass loss via the vertical disk surfaces by evaluating

$$\begin{eqnarray}&&{{\dot{M}}_{{\rm wind}}}\equiv \left. 2\pi \int _{r={{r}_{{\rm i}}}}^{{{r}_{{\rm o}}}}\;\rho {{v}_{\theta }}r\;{\rm sin} \theta \;dr\; \right|_{\theta ={{\theta }_{1}}}^{{{\theta }_{2}}}\end{eqnarray} \tag{ 7 }$$

at the upper and lower disk surface, $\theta ={{\theta }_{1}},{{\theta }_{2}}$, respectively. For model O-b6, we find a net mass loss rate of ${{\dot{M}}_{{\rm wind}}}=1.47\pm 0.37\times {{10}^{-8}}\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}}$ (see also Table 2).

Table 2.  Summary of Simulation Results

	${{z}_{b}}$	${{z}_{A}}$	$T_{R\phi }^{{\rm Reyn}}$	$T_{R\phi }^{{\rm Maxw}}$	$T_{z\phi }^{{\rm Maxw}}$	${{\dot{M}}_{{\rm wind}}}$	${{\dot{M}}_{{\rm accr}}}$
	(H)	(H)	$({{10}^{-6}}{{p}_{0}})$	$({{10}^{-5}}{{p}_{0}})$	$({{10}^{-5}}{{p}_{0}})$	$({{10}^{-8}}\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}})$	$({{10}^{-8}}\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}})$
O-b6	⋯	7.60 ± 0.45	6.87 ± 14.4	7.44 ± 0.95	⋯	1.58 ± 0.59	0.14 ± 1.53
OA-b5	5.23 ± 0.07	7.31 ± 0.17	3.63 ± 0.19	2.22 ± 0.06	9.82 ± 0.08	0.79 ± 0.01	0.43 ± 0.01
OA-b6	7.22 ± 0.48	6.94 ± 0.37	−0.21 ± 0.18	0.79 ± 0.06	0.58 ± 0.06	0.19 ± 0.06	0.09 ± 0.02
OA-b7	7.31 ± 0.70	6.39 ± 0.16	0.07 ± 0.13	$\lt 0.01$	0.22 ± 0.01	0.03 ± 0.01	0.00 ± 0.03
OA-b5-d4	5.27 ± 0.07	7.33 ± 0.18	0.11 ± 0.30	2.88 ± 0.11	9.88 ± 0.12	0.76 ± 0.02	0.34 ± 0.04
OA-b5-flr	4.81 ± 0.03	6.90 ± 0.31	0.26 ± 0.21	1.78 ± 0.02	14.3 ± 0.02	1.57 ± 0.01	0.64 ± 0.02
OA-b5-flr-nx	4.78 ± 0.03	7.50 ± 0.30	2.28 ± 9.24	1.87 ± 0.04	13.0 ± 0.04	1.58 ± 0.01	0.64 ± 0.03
OA-b5-nx	5.10 ± 0.04	7.34 ± 0.13	0.94 ± 9.29	1.89 ± 0.06	7.84 ± 0.02	0.67 ± 0.01	0.30 ± 0.04

Note. The vertical position of the base of the wind, ${{z}_{b}}$, and the Alfvén point, ${{z}_{A}}$, are found independent on the radial location when measured in local scale heights, H. The viscous accretion stresses ${{T}_{R\phi }}$ are vertically integrated within $|z|\leqslant {{z}_{b}}$—note the different units for Reynolds and Maxwell stresses. The wind stress, ${{T}_{z\phi }}$, is measured at $z=\pm {{z}_{b}}$. All stresses depend weakly on radius; listed values are at $r=3\;{\rm AU}$.

Download table as:  ASCIITypeset image

We observe the formation of field arcs and material flowing downward along field lines toward the arc foot points. This can be seen in the lower disk half around $r=2.8$, and $r=4.3\;{\rm AU}$, for example, and illustrates the potential role of buoyancy instability in regulating the disk's magnetization. Amplification of the magnetic field through the growth of channel modes, and their break up through the action of parasites, apparently leads to the formation of these locally uprising field structures. As an alternative explanation, we remark that the dynamic evolution of such magnetic loops has previously been attributed to the effect of differential rotation onto the magnetic footpoints (Romanova et al. 1998).

In contrast, within the magnetically decoupled midplane layer, the magnetic field remains predominantly vertical, retaining the initial field configuration. Near the dead-zone edge, the effect of Ohmic diffusion gradually declines. There, the azimuthal magnetic field is allowed to change its value, resulting in localized current layers. This is very similar to the situation obtained in the disks that include AD, as discussed later.

We illustrate the vertical disk structure in Figure 3, where we plot, in the upper panel, the total Reynolds and Maxwell stresses relative to the initial midplane gas pressure, p₀. The profiles show the typical signature of a laminar dead zone and MRI-turbulent surface layers (Fleming & Stone 2003; Oishi et al. 2007; Gressel et al. 2011, 2012), where the mechanical and magnetic shear-stresses lead to transport of angular momentum. Because of the relatively weak net-vertical magnetic flux, the level of turbulence is modest, and the vertically integrated dimensionless Maxwell stress is $\simeq 8\times {{10}^{-5}}$. With the "viscous" mass accretion rate (i.e., that effected by internal stresses) approximated by

$$\begin{eqnarray}&&{{\dot{M}}_{{\rm visc}}}\approx {{\left. 3\pi \;{{{\Omega }}^{-1}}\int _{\theta ={{\theta }_{1}}}^{{{\theta }_{2}}}\left( T_{R\phi }^{{\rm Reyn}}+T_{R\phi }^{{\rm Maxw}} \right)r\;d\theta \; \right|}_{r={{r}_{{\rm ref}}}}},\end{eqnarray} \tag{ 8 }$$

we estimate a value of ${{\dot{M}}_{{\rm visc}}}\simeq 0.3\times {{10}^{-8}}\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}}$, at ${{r}_{{\rm ref}}}=2.5\;{\rm AU}$, if mass accretion were effected by turbulent transport of angular momentum alone. This can be compared to the actual mass accretion rate, which we can easily compute in our global model via

$$\begin{eqnarray}&&{{\dot{M}}_{{\rm accr}}}\equiv {{\left. 2\pi \int _{\theta ={{\theta }_{1}}}^{{{\theta }_{2}}}\;\rho {{v}_{r}}\;{{r}^{2}}{\rm sin} \theta \;d\theta \; \right|}_{r={{r}_{{\rm ref}}}}}.\end{eqnarray} \tag{ 9 }$$

Applying a radial average over $r\in [2,3]\;{\rm AU}$, and estimating fluctuations occurring within a time interval spanning $t\in [75,100]\;{\rm yr}$, we find ${{\dot{M}}_{{\rm accr}}}=(0.14\pm 1.53)\times {{10}^{-8}}{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}}$ (see also the last column of Table 2). Clearly, this diagnostic turns out to be subject to strong fluctuations for model O-b6—presumably due to the presence of strong radial motion within the MRI channel modes. Equation (9) will however prove useful when evaluating the radial drift of material in the localized current layers seen in the AD-dominated disk models. We note at this point that a simple energy conservation argument prevents the steady-state mass loss rate (to infinity) through a magneto-centrifugal wind being larger than the accretion rate through the disk. The large value of ${{\dot{M}}_{{\rm wind}}}$ reported above relative to ${{\dot{M}}_{{\rm accr}}}$ suggests that the wind loss rate is not yet converged. Indeed we note that BS13 report that increasing the vertical sizes of their shearing boxes leads to a reduction of the wind mass loss rates. A similar conclusion is reached by Fromang et al. (2013) for winds launched from turbulent disks. It appears that obtaining accurate and physically meaningful estimates for the mass fluxes through the upper boundaries of the disk will require simulations that are truly global in the vertical direction, such that the fast magnetosonic point is contained within the simulation domain (rather than outside of the boundary as it is at present for all of the models that we present in this paper, ensuring that the wind launching region is causally disconnected from the imposed boundary conditions).

Zoom In Zoom Out Reset image size

In the lower panel of Figure 3, the gas pressure, kinetic energy and magnetic pressure, are shown. Within the dead-zone layer, between $z=\pm 4\;H$, the disk is supported by gas pressure alone, which follows a Gaussian profile. Owing to the additional magnetic pressure support within the active layer, the effective scale height increases roughly twofold there. Note that unlike seen in the isothermal simulations of BS13, cf. their Figure 3, the magnetic pressure does not significantly exceed the value of the gas pressure in large parts of the disk corona. We attribute this difference to the dissipation heating (Hirose & Turner 2011) present in our simulations, which we note, however, may be overly pronounced in the axisymmetric models. The kinetic energy only becomes important very close to the vertical boundaries of our model. This is also reflected in the position of the Alfvén point, which is relatively poorly constrained and which we infer as $(7.60\pm 0.45)H$. Values for these vertical points are listed in Table 2 for all models. We report time- and space-averaged values obtained for $r\in [1,5]\;{\rm AU}$, and we have appropriately taken into account the flaring of the disk. Because of the turbulent character of the upper disk layers, we cannot obtain a good estimate for the location of the base of the wind in model O-b6.

We now discuss models that include both Ohmic resistivity and ambipolar diffusion. A key finding of BS13 is that the inclusion of AD inhibits linear growth of the MRI in the intermediate disk layers, that is, the regions that are not affected by Ohmic resistivity. Since the AD coefficient depends on the plasma parameter, one might expect that the effect of AD is less severe for high values of ${{\beta }_{P}}$. To test this, we have run simulations OA-b6, and OA-b7, with ${{\beta }_{p\,0}}={{10}^{6}}$, and 10⁷, respectively.

3.3.1. AD with Weak Vertical Fields

3.3.2. The Fiducial Model

In Figure 4, we show vertical profiles of the $R\phi$ stress components for our fiducial model, OA-b5, where owing to the stronger net vertical field the MRI is ultimately suppressed by AD and a laminar wind solution is established instead. Dashed vertical lines indicate the base of the wind, which we define according to Wardle & Königl (1993) as the vertical position, $z={{z}_{b}}$, where the rotation becomes super-Keplerian. Since we define our azimuthal velocity, ${{v}_{\phi }}$ as the deviation from the Keplerian background flow, this can be verified in Figure 4 by the fact that the Reynolds stress, $T_{R\phi }^{{\rm Reyn}}\equiv \rho {{v}_{R}}\;{{v}_{\phi }}$, vanishes at this point. The same is true for the vertical tensor component, $T_{z\phi }^{{\rm Reyn}}$, which makes it convenient to define the wind stress,

$$\begin{eqnarray}&&{{T}_{z\phi }}\equiv {{\left. T_{z\phi }^{{\rm Maxw}}\; \right|}_{z=+{{z}_{b}}}}-{{\left. T_{z\phi }^{{\rm Maxw}}\; \right|}_{z=-{{z}_{b}}}}\end{eqnarray} \tag{ 10 }$$

at this point. For the $R\phi$ component of the stress tensor, responsible for redistributing angular momentum radially, we furthermore obtain average values within the disk body $z\in [-{{z}_{b}},+{{z}_{b}}]$,

$$\begin{eqnarray}&&{{T}_{R\phi }}\equiv \frac{1}{2{{z}_{b}}}\int _{-{{z}_{b}}}^{+{{z}_{b}}}\left( T_{R\phi }^{{\rm Maxw}}+T_{R\phi }^{{\rm Reyn}} \right)dz,\end{eqnarray} \tag{ 11 }$$

which we list separately in Table 2, where values have been normalized to the midplane pressure, p₀. Already for a midplane plasma parameter as low as 10⁵, the wind stress exceeds the $R\phi$ (accretion) stress by a factor of four. BS13 also report that the wind stress exceeds the accretion stress; for their fiducial run, the discrepancy is bigger than an order of magnitude.

Zoom In Zoom Out Reset image size

In the lower panel of Figure 4, we see that in the absence of MRI turbulence, the hydrostatic balance extends further up in the disk and the gas pressure remains the dominant stabilizing force all the way up to the base of the wind. Even in the magnetically dominated wind layer, magnetic pressure gradients remain moderate and the scale height of the disk remains roughly constant up to ${{z}_{b}}$.

In Figure 5, we show a zoom-in of the inner part of the global field topology in our fiducial simulation OA-b5. The upper panel shows the resulting field configuration early on, when resistively modified MRI eigenmodes appearing inside $r=1.5\;{\rm AU}$ are still clearly discernible. Localized channel solutions form at $z\approx 3\;H$ near the interface between the magnetically decoupled region, characterized by ${{{\Lambda }}_{{\rm O}}}\ll 1$ and purely vertical field lines with ${{B}_{\phi }}\approx 0$, and the AD-dominated intermediate layer. In this small region, both Elsasser numbers are above unity (see Figure 1), allowing for linear growth of MRI channel modes. This transition layer is characterized by a zigzag-shaped abrupt change in the field lines (see Figure 8 below), associated with strong current sheets. At this point, the laminar wind solution is already well established in the FUV ionized surface layer, and has spread throughout the radial extent of the simulation domain.

Zoom In Zoom Out Reset image size

While the disk wind itself has already reached a steady state at this point in time, the strong field layers continue to evolve. In the lower panel of Figure 5, we accordingly show a later evolution stage, where the relic MRI channels have morphed into strong azimuthal field belts. These field belts can be understood as "undead zones," that is, a region that does not sustain MRI, but where diffusion can be balanced by the winding-up of pre-existing radial field via differential rotation (Turner & Sano 2008). It appears that the initial amplification of the magnetic field through the growth of channel modes causes the AD coefficient to increase until it quenches further development of the MRI. This notion is supported by a reference simulation, where the diffusion coefficients ${{\eta }_{{\rm O}}}$ and ${{\eta }_{{\rm AD}}}$ were held fixed at their initial value, and where the MRI continues to produce quasi-turbulent fields akin to the ones seen in Figure 2 for the Ohmic-only case. The reversal of the field direction seen at $r\simeq 1.5\;{\rm AU}$ is a relic of the local nonlinear development of the MRI modes early on (as may be seen in the upper panel at $r\simeq 0.75\;{\rm AU}$). This interface between the azimuthal field belts of opposite parity moves radially outward at a speed of $\approx 0.01\;{\rm AU}\;{\rm y}{{{\rm r}}^{-1}}$.

The emerging current layers are directly related to resistively modified vertically global MRI channel modes—see Latter et al. (2010), and Figure 8 in Section 3.5. With MRI channels potentially spanning large radial extents, as in the Ohmic run shown in Figure 2, such layers may be fed diffusively with magnetic field from the MRI-active inner disk. This is supported by our present models, where current layers appear to emerge from inner disk regions. What determines the exact shape of the surviving field configuration at late times is presently unclear. We speculate that the particular realization seen in the lower panel of Figure 5 is not necessarily a unique solution given the specific disk geometry and ionization model, but may to some extend depend on the early evolution history. However, simulations that were initiated with different random seeds, but were otherwise identical, only showed minor variations in the final field appearance. The configurations observed at late times are moreover found to remain quasi-stationary, at least during the evolution time of $100\;{\rm yr}$ that we have investigated.

Perhaps the most noteworthy result from run OA-b5 is the spontaneous emergence of the expected "hourglass" geometry for the magnetic field, allowing magnetic tension forces to accelerate the wind. For the adopted vertical field polarity, the field must bend radially outward near the upper and lower disk surfaces, and the azimuthal field should point in the positive direction in the lower hemisphere and reverse direction above the midplane, as observed in the upper panel of Figure 5.

3.4.1. Robustness of the Emerging Wind Geometry

3.4.2. Comparison with Local Models

Except for the reversal of the horizontal magnetic field components (which at this stage appears to occur rather gently instead of within a narrow current layer), the field structure observed in the upper panel of Figure 5 is similar to that described by BS13 and shown in their Figure 10. Near the midplane, where the dominance of the Ohmic resistivity retards the growth of currents, the field is dominated by the vertical component. Moving to higher altitudes, where the magnetic coupling increases, the field bends outwards because a radial velocity is generated through the azimuthal force balance between the Coriolis force and magnetic tension. We enter the FUV layer at disk altitudes $|z|\gtrsim 4.5\;H$, coinciding with the base of the wind and the wind itself, as illustrated by the velocity vectors.

At later times, the lower panel of Figure 5 shows that the varying polarity of the strong azimuthal field belts gives rise to two distinct field line geometries. In the inner region ($r\lesssim 1.5\;{\rm AU}$) the orientation of the field belts is the same as the one in the wind, and the vertical field lines bend smoothly at the transition into the wind base. Because some negative ${{B}_{\phi }}$ has diffused into the Ohmic-dominated midplane, there is a weak current layer forming at $z\simeq -2.15\;H$. In contrast, further out at $r\gtrsim 1.5\;{\rm AU}$, the horizontal-field belts are anti-aligned with the magnetic field direction of the wind. This can also be seen in the curvature of the field lines, which are concave (toward the star) in this region. To connect the different layers, the azimuthal and radial field components have to change their direction within two thin current layers located at $z\simeq \pm 3.2\;H$. To study these configurations in more detail, we plot separate vertical profiles for $r=1\;{\rm AU}$, and $r=2\;{\rm AU}$ in Figures 6 and 7, respectively. For easy reference, the radial positions are furthermore indicated in Figure 5 by means of dashed lines. The isothermal sound speeds, for the initial disk model, are ${{c}_{{\rm s}}}=1.5$, and $1.1\;{\rm km}\;{{{\rm s}}^{-1}}$ at the two respective radii.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The wind profile in the inner disk (see Figure 6), is qualitatively similar to the wind solutions previously found in the quasi-1D simulations with "even-z" symmetry, plotted in Figure 10 of BS13. The current layer is much weaker in our case, and the dip in v_r is barely visible in the left inset of Figure 6. This may be related to differences in the ionization model resulting in ${{{\Lambda }}_{{\rm AD}}}\gt 1$ in the region around $|z|\simeq 3\;H$, in our disk model. Differences in the particular diffusivity profiles also explain the pronounced peaks in ${{B}_{\phi }}$ in this "undead" transition layer, which are not seen in the simulations of BS13.

The profile at $r=2\;{\rm AU}$, with double field reversal in ${{B}_{\phi }}$, is illustrated in Figure 7. Because of the strong current layers that are emerging in this case, this configuration resembles more closely the situation seen in the local models (again see Figure 10 in BS13). Unlike in their simulation, we however observe two current layers, one on each side of the disk. These arise because the development of MRI channel modes, and their subsequent evolution into "undead" azimuthal field belts with opposite polarity to the background wind, forces the horizontal field to change direction three times. Strong current layers occur at altitudes where the field and gas couple more strongly, namely at $|z|\sim 3\;H$. Although the occurrence of strong azimuthal field belts is not precluded in local simulations (as their existence may depend on the details of the ionization model), it seems likely that the alternating field belts we observe will only be a robust feature of global simulations. Apart from their multiplicity, the current layers appear very similar to the ones presented in BS13. The radial velocity v_r shows a dip, while ${{v}_{\phi }}$ is modulated as a sine (the slight offset of ${{v}_{\phi }}$ in our simulation is due to the radial pressure support). The characteristic modulation (reflecting the abrupt kink in the field lines) is clearly seen in B_r and ${{B}_{\phi }}$, which agree very closely with the B_x, and B_y obtained by BS13. The emergence of such strong current layers naturally begs the question of their stability. Before we however discuss this issue in Section 3.6, we first want to take a look at the current sheets' likely origin as remnants of MRI eigenmodes.

To draw a closer connection between localized MRI channel modes and the current layers observed in our simulations, we compare the structure of the fastest growing MRI eigenmode with the initial evolution of our standard setup with ${{\beta }_{p\,0}}={{10}^{5}}$. The mode structure we find is qualitatively similar to that shown for a spatially constant Ohmic resistivity (see the appendix of Fromang et al. 2013), and resembles the eigenmodes shown in Figure 4 of Lesur et al. (2014), who performed Hall-MHD calculations that also included variable Ohmic resistivity. Dissipation modifies the MRI eigenmodes predominantly near the midplane, where the relative magnetization is lowest and the unstable MRI wavenumbers are the highest, explaining the similarity of the solutions for constant η, and spatially varying $\eta =\eta (z)$.

Here, instead of a stratified shearing box framework as previously employed by Latter et al. (2010), we used the framework introduced by McNally & Pessah (2014) to capture the full vertical structure of the baroclinic initial condition. The eigenmode analysis uses the linear equations described in McNally & Pessah (2014, Section 5) with the addition of terms for spatially varying Ohmic resistivity, and using the functional forms and coefficients appropriate for the cylindrical temperature structure of the disk initial condition (McNally & Pessah 2014, Section 3.3). The equations were discretized in real space with Chebyshev cardinal functions on a Gauss–Lobatto grid with 1000 points with the same vertical extent as the global simulation (±8 scale heights). Boundary conditions were enforced by the "boundary bordering" method (Boyd 2000) forcing the magnetic field perturbations to zero at the vertical edges of the domain. The analysis was performed at radial position $R=1\;{\rm AU}$ where the Ohmic resistivity profile, ${{\eta }_{{\rm O}}}(z)$, used was extracted from the simulation initial condition and interpolated to the grid used for the eigenmode calculation.

We moreover obtained eigenmode solutions that additionally included a fixed ${{\eta }_{{\rm AD}}}(z)$ profile, and these showed a nearly indistinguishable mode structure. Because of the extra nonlinearity introduced by the ${{\beta }_{P}}$ dependence of ${{\eta }_{{\rm AD}}}(z)$, we however restrict our comparison to the Ohmic-only case. In Figure 8, we compare the early stage of our fully nonlinear global simulation with the fastest growing eigenmode. The localized field reversals are also clearly seen in the upper panel of Figure 5, which shows the early evolution of model OA-b5 (however including AD). To isolate the linear MRI mode from the overlaid wind solution, we subtract respective field profiles from two distinct snapshots during the linear growth phase of the MRI. This is possible because the wind configuration emerges quickly, and then provides an approximately stationary background in which the MRI grows. Differencing two snapshots in time hence allows us to remove the wind part and extract the exponentially growing MRI eigenmode. While growing perturbations near the midplane are high-wavenumber and thus efficiently damped, the relative field strength in the disk halo is too strong to allow for the MRI modes to fit within the available domain. This only leaves an intermediate disk region near $|z|\simeq 3.5\;H$ for MRI channel modes. In the Ohmic+AD simulations, this corresponds to the position of the "undead" field belts (see lower panel of Figure 5) with adjacent current layers. In those simulations, the back-reaction of the growing eigenmodes on the ${{\eta }_{{\rm AD}}}(z)$ profile makes the MRI self-quenching, leaving the current layers as a remnant of the initial MRI channel mode.

Zoom In Zoom Out Reset image size

While we find the current layers to remain stable and long-lived in model OA-b5, some of the other simulations which we will discuss below also show signs of instability, leading to less regular fields.

3.6.1. Effect of Disk Flaring

An example of this is illustrated in Figure 9, where we show the vertical velocity, v_z (upper panel), and the azimuthal magnetic field, ${{B}_{\phi }}$ (lower panel), at the end of the flaring disk simulation OA-b5-flr. Both color scales have been restricted to highlight features in the current layer. In this model, the surface density interior to $1\;{\rm AU}$ is smaller than in OA-b5, and beyond $1\;{\rm AU}$ it is larger, altering the ionization structure sufficiently to put the current layer in a different regime. In particular, the increased ionization in the inner disk apparently allows the azimuthal field belts that form at intermediate heights to reach larger amplitudes. This, combined with the fact that these field belts are of opposite polarity relative to the background wind field, causes the two current layers located at $z\simeq \pm 3\;H$ to be stronger and to support narrow, radially directed flows, or "jets," similar to those in run OA-b5, as illustrated in Figure 7.

Zoom In Zoom Out Reset image size

Inside $r\simeq 1.4\;{\rm AU}$, the velocity near the current layer shows a distinct radial wavenumber, suggesting that the shear associated with the radial jets leads to growth of a Kelvin–Helmholtz (KH) instability. The pattern, however, remains more or less stationary, indicating that the instability is not growing significantly into the nonlinear stage. Although not obvious in Figure 9, wave-like perturbations to the magnetic field are also present in the current layer. The possibility of current-driven tearing modes being present also arises, but analysis of the shear rate associated with the radial jets suggests that the observed instability is driven by KH modes. The velocity associated with the shear is super-Alfvénic, and is therefore expected to stabilize the tearing modes that may in principle grow in a shearing background (Chen et al. 1997).

The wavelength of the velocity perturbations and magnetic field changes abruptly for $r\gtrsim 1.4\;{\rm AU}$, where now a characteristic vortex-pair pattern emerges. The exact reason for this abrupt transition is unclear, since most characteristic quantities, such as the Elsasser numbers ${{{\Lambda }}_{{\rm O}}}$, and ${{{\Lambda }}_{{\rm AD}}}$, only depend weakly and continuously on radius. The most likely explanation may be the radial variation of the azimuthal field (inside the "undead" belts), which drops in strength around this radial location. This is accompanied by a weakening and widening of the current layer, and a significant reduction in the inflow velocity of the radial accretion flow associated with the current layer, and a corresponding reduction in the shear relative to the background fluid. This suggests that the abrupt change in the nature of the disturbance observed at radius $\sim 1.8\;{\rm AU}$ is unlikely to be driven by a KH instability, but the weakening of the shear may allow a tearing-mode instability to develop instead. The appearance of disturbances that have the character of growing magnetic islands provides some support for this interpretation, but this is clearly a tentative explanation that will require further investigation in future work. The presence of a background wind, Keplerian shear, AD and a radially directed inflow associated with the current layer make an analysis of this problem somewhat complicated, and beyond the scope of this paper. We note, however, that the local studies on the development of parasitic modes in resistive MHD (Latter et al. 2009; Pessah 2010) may provide some insight, even though AD is neglected in those studies. Pessah (2010) has analyzed MRI parasitic modes in the presence of finite Ohmic resistivity and viscosity. In this paper, ${{{\Lambda }}_{{\rm O}}}$, is identified as the relevant dimensionless number (see their Figure 11), with tearing modes being the dominant parasitic modes for ${{{\Lambda }}_{{\rm O}}}\lesssim 1$, and KH modes dominating otherwise. Considering both ${{{\Lambda }}_{{\rm O}}}$, and ${{{\Lambda }}_{{\rm AD}}}$, we assess that our simulations are situated near this transition point.

3.6.2. Effect of Reduced Dust Fraction

Figure 10 shows the final appearance of the entire disk for model OA-b5-d4, which has a reduced dust fraction of 10⁻⁴ compared to a value of 10⁻³ in model OA-b5. On one hand, as the reader may check from Table 2, the two models are rather similar in their bulk properties, especially those related to the wind. This is unsurprising, since the reduced cross-section for recombination on grains affects the ionization fraction significantly only at large gas columns (see Figure 1), below the base of the wind. The grains are less important for the ionization state in the FUV layer, where the electrons are so abundant that many remain free once the grains charge to the Coulomb limit (Okuzumi 2009).

Zoom In Zoom Out Reset image size

In contrast to our fiducial model, the transition between the base of the wind and the Ohmic-dominated disk interior is very chaotic in model OA-b5-d4, whereas it is much more stable in model OA-b5. It appears that the larger values of ${{{\Lambda }}_{{\rm O}}}$ and ${{{\Lambda }}_{{\rm AD}}}$ in this model at heights $z\sim \pm 2.5\;H$ allow the channel modes to grow further into the nonlinear regime, such that parasitic modes are able to develop and cause the field belts to break up into regions with locally coherent field. The field amplification during this phase means AD strengthens, as in run OA-b5. Eventually the local evolution is driven by a combination of field winding and diffusion. The end result is the creation of local "undead" patches that evolve on timescales that are comparable to the run duration.

The different evolution of the horizontal field belts observed in runs OA-b5-flr and OA-b5-d4 compared with OA-b5 show that the evolution of the disk interior depends critically on the precise disk model that determines the Elsasser number profiles. In particular, the ionization state in the region between the Ohmic-dominated disk midplane and the AD-dominated intermediate disk layers depends sensitively on various model parameters. This suggests that a careful parameter study is required to obtain a reliable picture of how PPDs evolve over their life-times—especially since their properties change substantially over these time scales. Such a parameter study would appear to be most warranted when additional physics is included in the model, such as the Hall effect and more realistic thermodynamics. In a situation where the gas is allowed to be heated by resistive effects, the ionization balance may be shifted toward higher conductivity in the vicinity of the current layers, potentially leading to intermittent behavior (Faure et al. 2014; McNally et al. 2014).

A distinct instability emerges in the magnetically dominated corona of the laminar disk (see Figure 11). The instability is found to grow close to the smallest available wavenumber, and is only present in the highly resolved axisymmetric simulations. The perturbations are roughly aligned with the outward-pointing field lines of the laminar wind solution. Similar patterns have been observed in 2D-axisymmetric shearing box simulations with strong vertical fields performed by Moll (2012), who also finds field-aligned structures on the shortest length scales resolved on the grid. The instability also bears resemblance with the "striped wind" seen by Lesur et al. (2013). Similar to the instability seen by Moll (2012), we find a moderate modulation of the vertical mass flux by the instability pattern. We however do not find material falling down along the field lines as reported there for certain parameters.

Zoom In Zoom Out Reset image size

We identify three candidate mechanisms for creating such disturbances: (i) the vertical-shear instability investigated recently by Nelson et al. (2013), potentially modified by anisotropy from field-line tension, (ii) a KH-type instability related to the vertical gradient of the horizontal velocity near the base of the wind, and (iii) a buoyant interchange instability of the azimuthal field. It is important to note that the KH-type instability can create the elongated aspect ratio seen in Figure 11, when perturbations are carried-away by the overlaid laminar wind velocity—and we indeed find this to be the case in situations where the instability occurred before the wind was established. The conclusive hint that a KH instability causes the pattern in Figure 11 is a strong velocity gradient at the exact location of the transition into the FUV layer. This gradient can be clearly identified in Figure 7 (situated just above the wind base at ${{z}_{b}}$). The shear layer manifests itself as an additional steepening of the wind profiles, and can also be seen in the corresponding plot of BS13 (see their Figure 10).

The physical origin of this sharp velocity gradient can be traced back to a peak in the ambipolar part of the electromotive force in the induction equation. Because the coefficient ${{\eta }_{{\rm AD}}}(r,\theta )$ sharply drops at the position of the FUV interface, the term $\nabla \times ({{\eta }_{{\rm AD}}}\;[(\nabla \times {\boldsymbol{B}} )\times {\boldsymbol{\hat{b}}} \;]\times {\boldsymbol{\hat{b}}} )$ can become significant in the induction equation even for smoothly varying (or constant) ${\boldsymbol{B}}$. In the steady state, this curl of the ambipolar electromotive force needs to be balanced by the induction term, $\nabla \times ({\boldsymbol{v}} \times {\boldsymbol{B}} )$, requiring as a consequence a significant shear gradient in the horizontal velocity (again assuming a smooth ${\boldsymbol{B}} )$. Given its strong dynamical effect, this raises the question whether the coincidence of the wind base and the FUV interface is in fact accidental as claimed by BS13.

In our disk model, the FUV layer is imposed ad hoc, with a prescribed column density (Section 2.1.1). The resulting interface is smeared out to a prescribed width. Reflecting the predominant direct origin of the FUV photons, while scattering is expected to remain negligible, the tapering is done in the radial direction, and as a function of column density (rather than position). As a consequence, the resulting ${{\eta }_{{\rm AD}}}(r,\theta )$ may develop sharp steps on the numerical grid in the θ direction—in turn resulting in sharp velocity gradients. We accordingly ran a reference model, where the FUV transition was softer. In this model, we did not observe a strong velocity gradient near the wind base, and also found the instability to be largely suppressed. Also the flaring-disk model, OA-b5-flr, did not show any signs of the instability. This can potentially be explained by the FUV transition following a concave line in this case, producing a different aliasing of the tapering function on the spherical-polar mesh. We conclude that the instability seen in Figure 11 is likely caused by the particular numerical parametrization of the FUV layer in our model, and is entirely unrelated to the "striping" instability seen in Moll (2012) and Lesur et al. (2013). The instability may be be triggered by numerical effects at the grid scale, in a manner analogous to how physical secondary KH instabilities can be triggered by purely numerical effects as a shear layer thins to the grid scale (McNally et al. 2012). We remark that in the case of this wind instability too the mechanism behind the instability is genuinely physical, and may occur in real PPDs where there is a sharp local change in the ionization state.

The presence of secondary instability mandates investigating the non-axisymmetric evolution. This is done exemplarily for the two models OA-b5, and OA-b5-flr, whose 3D counterparts OA-b5-nx, and OA-b5-flr-nx are also listed in Table 2, showing only minor deviations from their respective 2D runs. Given the moderate magnetic Reynolds numbers in these regions, it is not per se obvious whether non-axisymmetric instabilities can transfer significant amounts of energy into modes with higher azimuthal wavenumber. To address this question, in Figure 12, we show vertically averaged azimuthal power spectra of the magnetic field. The color-coding indicates the logarithm of the power spectrum of ${{B}_{\theta }}$ as function of wavenumber, m, and radial position in the disk. In the upper panel, we show the result for model OA-b5-flr. With the exception of the KH features (between $R=1.5-2.0\;{\rm AU}$) also seen in Figure 9, almost all power resides in m = 0. Even in the regions of the secondary instabilities, the power in modes with $m\gt 0$ remains very small. The power seen adjacent to the inner radial domain boundary is caused by reducing the diffusion coefficients there, which is an attempt to mimic the MRI-active inner disk (see Section 2.3 for details).

Zoom In Zoom Out Reset image size

This can be contrasted to the case of the classic layered PPD with MRI-active surface layers (see model M1 in Gressel et al. 2013), which we plot in the lower panel of Figure 12. Even with the laminar midplane layer, a significant amount of power is found for high-m modes (note the different color bar in this figure). The lack of power near the radial boundaries is owing to a dissipative buffer region. Gressel et al. (2013) focused on the effects of a planet, which they inserted after the time shown here. We hence conclude that in the simulations that are dominated by axisymmetric configurations of AD-laminar winds, secondary instabilities do not lead to significant power in non-axisymmetric perturbations. This may potentially change in the presence of stronger current layers that are prone to nonlinear stage of parasitic modes as, for instance, can be seen for model OA-b5-d4 with a stronger dust depletion. When combined with stronger net-vertical fields and flaring, these dust-depleted disks may develop intermediate layers that are on the verge of becoming turbulent.