Thanatology in protoplanetary discs
The combined influence of Ohmic, Hall, and ambipolar diffusion on dead zones^{⋆}
^{1}
Univ. Grenoble Alpes, IPAG,
38000
Grenoble,
France
email:
geoffroy.lesur@ujfgrenoble.fr
^{2}
CNRS, IPAG, 38000
Grenoble,
France
^{3}
Department of Astrophysical Sciences, 4 Ivy Lane, Peyton Hall,
Princeton University, Princeton
NJ
08544,
USA
^{4}
Laboratoire AIM, CEA/DSM–CNRS–Université Paris 7, Irfu/Service
d’Astrophysique, CEASaclay, 91191
GifsurYvette,
France
Received:
17
February
2014
Accepted:
8
April
2014
Protoplanetary discs are poorly ionised due to their low temperatures and high column densities and are therefore subject to three “nonideal” magnetohydrodynamic (MHD) effects: Ohmic dissipation, ambipolar diffusion, and the Hall effect. The existence of magnetically driven turbulence in these discs has been a central question since the discovery of the magnetorotational instability (MRI). Early models considered Ohmic diffusion only and led to a scenario of layered accretion, in which a magnetically “dead” zone in the disc midplane is embedded within magnetically “active” surface layers at distances of about 1–10 au from the central protostellar object. Recent work has suggested that a combination of Ohmic dissipation and ambipolar diffusion can render both the midplane and surface layers of the disc inactive and that torques due to magnetically driven outflows are required to explain the observed accretion rates. We reassess this picture by performing threedimensional numerical simulations that include all three nonideal MHD effects for the first time. We find that the Hall effect can generically “revive” dead zones by producing a dominant azimuthal magnetic field and a largescale Maxwell stress throughout the midplane, provided that the angular velocity and magnetic field satisfy Ω·B > 0. The attendant large magnetic pressure modifies the vertical density profile and substantially increases the disc scale height beyond its hydrostatic value. Outflows are produced but are not necessary to explain accretion rates ≲ 10^{7} M_{⊙} yr^{1}. The flow in the disc midplane is essentially laminar, suggesting that dust sedimentation may be efficient. These results demonstrate that if the MRI is relevant for driving mass accretion in protoplanetary discs, one must include the Hall effect to obtain even qualitatively correct results.
Key words: accretion, accretion disks / instabilities / magnetohydrodynamics (MHD) / protoplanetary disks / stars: formation
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2014
1. Introduction
The formation and evolution of protoplanetary discs is a longstanding problem in accretion theory. Observations indicate that T Tauri stars accrete mass from their surrounding discs at a rate of 10^{9}–10^{7} M_{⊙} yr^{1} on timescales of about 1–10 Myr (Hartmann et al. 1998). However, the physical mechanism driving this accretion has remained elusive for decades. Such a mechanism is not only necessary for understanding how the central protostellar object accretes mass at the observed rates but also represents a crucial component in any theory for the formation of planets. To effectively accrete, angular momentum must be removed from the infalling gas and transferred to other fluid elements within the disc and/or in the surrounding medium. It can either be transported radially, in which case the transport is usually modelled as a viscous process (socalled αdisc models; Shakura & Sunyaev 1973), or transported vertically by largescale magnetically driven outflows or winds (Blandford & Payne 1982).
The rediscovery of the magnetorotational instability (MRI; Velikhov 1959; Chandrasekhar 1960) by Balbus & Hawley (1991, 1998) in the context of astrophysical discs vastly improved our understanding of enhanced angularmomentum transport in differentially rotating plasmas. This instability produces vigorous turbulence and efficiently transports angular momentum outwards in discs with an effective viscosity α ≳ 10^{3} (Hawley et al. 1995). The MRI has since become the main contender^{1} to explain accretion in many types of astrophysical systems. However, protoplanetary discs are cold (T ≲ 300 K at 1 au), dense (n ≳ 10^{13} cm^{3}), and thus poorly ionised. Their evolution is welldescribed by a set of magnetohydrodynamic (MHD) equations that include three nonideal effects (Norman & Heyvaerts 1985; Wardle & Ng 1999): Ohmic dissipation, ambipolar diffusion, and the Hall effect. In a magnetised plasma composed of neutral molecules, ions, and electrons, Ohmic dissipation is caused by collisions between electrons and neutrals, ambipolar diffusion by collisions between ions and neutrals, and the Hall effect by the velocity difference (“drift”) between electrons and ions (or, more generally, between positively and negatively charged species). These nonideal effects tend to decouple the gas from the magnetic field, casting doubt upon the relevance of the MRI to protoplanetary discs (Blaes & Balbus 1994; Wardle 1999; Kunz & Balbus 2004; Desch 2004).
The impact of Ohmic dissipation on stratified protoplanetary discs was first studied by Gammie (1996) and Sano & Miyama (1999) and led to a model of layered accretion at distances of about ~1–10 au from the central protostellar object: a quiescent magnetically “dead” midplane that is surrounded by turbulent magnetically “active” surface layers located at a few disc scale heights. This layeredaccretion model has been studied extensively with increasingly complex analytical models (Sano et al. 2000; Fromang et al. 2002; Ilgner & Nelson 2006a,b) and threedimensional numerical simulations (Fleming et al. 2000; Fleming & Stone 2003; Turner et al. 2007; Ilgner & Nelson 2008). Depending upon the chemical network and the presence of well–mixed dust grains, these models predict a vertically integrated turbulent transport that is roughly consistent with the observed accretion rates. However, recent work by Bai & Stone (2013b) and Simon et al. (2013a) has shown that ambipolar diffusion can render the surface layers magnetically inactive. Magnetically driven outflows, which naturally occur in stratified shearing boxes with a mean vertical field (Suzuki & Inutsuka 2009; Moll 2012; Ogilvie 2012; Lesur et al. 2013; Fromang et al. 2013; Bai & Stone 2013a), are then needed to produce accretion rates comparable to those observed.
Although the influence of the Hall effect on the linear properties of the MRI is well understood (Wardle 1999; Balbus & Terquem 2001), its impact on the nonlinear evolution and, in particular on the dynamics of putative dead zones, has remained mostly speculative (Sano & Stone 2002; Salmeron & Wardle 2003, 2005; Wardle 2007; Salmeron & Wardle 2008; Wardle & Salmeron 2012). In a recent publication, Kunz & Lesur (2013, hereafter KL13) discovered a new saturation mechanism for Halldominated magnetorotational turbulence, whereby the vertical magnetic field selforganises into coherent, longlived axisymmetric (“zonal”) structures and the rate of angularmomentum transport becomes vanishingly small. While interesting in itself, this result was based on a crude representation of a protoplanetary disc – an incompressible, dustfree and unstratified shearing box without ionneutral drift. It is natural to ask whether this result applies to stratified models that include Ohmic dissipation and ambipolar diffusion.
This is precisely the question we address in this paper. We perform threedimensional numerical simulations of the MRI in a compressible, vertically stratified shearing box, taking Ohmic dissipation, ambipolar diffusion, and the Hall effect into account. This marks the first time that all three nonideal magnetohydrodynamic (MHD) effects have been accounted for in a numerical simulation of the MRI. To accomplish this feat, we have implemented an original formulation of ambipolar diffusion and the Hall effect in the finitevolume code Pluto. A simplified ionisation model is used to determine vertical diffusivity profiles for all three nonideal MHD effects at distances of 1, 5, and 10 au from the central protostellar object.
The paper is organised as follows. In Sect. 2.1, we present the equations of nonideal MHD in the stratified shearingbox approximation. In Sects. 2.2 and 2.3, we introduce our numerical method and define the diagnostics used to quantify the efficiency of angularmomentum transport and the production of outflows. We defer to Appendices A and B for the technical details and tests of our numerical implementations of the Hall effect and ambipolar diffusion. We then present the model used to determine the ionisation rate and the diffusivity tensor (Sect. 2.4). The results are presented, first with models that neglect ambipolar diffusion (Sects. 3.1–3.2) and then with all three nonideal MHD effects treated selfconsistently (Sect. 3.3). Since our ionisation model does not account for the chemical and dynamical effects of dust grains, we briefly address their impact on our results in Sect. 4. Finally, we provide a summary of our results and comment on their implications for the evolution of protoplanetary discs and the formation of planetesimals in Sect. 5.
2. Method of solution
2.1. The stratified shearingbox model
We study the evolution of a poorly ionised protoplanetary disc in the shearingbox approximation (Hawley et al. 1995). A local patch of the disc centred at R_{0} with an extent ~H (the vertical scale height) rotates with the disc at a constant angular speed Ω = Ω_{K}(R_{0}). The local frame is defined by (e_{x},e_{y},e_{z}), where e_{x} is aligned with the radial direction, e_{y} with the azimuthal direction, and e_{z} with the vertical direction. The computational domain has a size (L_{x},L_{y},L_{z}) with ; we assume the disc to be geometrically thin, which implies H ≪ R_{0}. This allows us to neglect curvature terms – the socalled Hill approximation (Hill 1878). Taking the flow to be locally isothermal, the equations of nonideal MHD that govern the evolution of the mass density ρ, the velocity u, and the magnetic field B are, respectively, where c_{s} is the isothermal sound speed, Ω = Ω e_{z} is the local angular velocity, g = 2qΩ^{2}xe_{x} − Ω^{2}ze_{z} is the local gravitational acceleration with q = −dlnΩ/dlnR (=3/2 for a Keplerian disc), e_{b} = B/  B , and J = ∇×B. We include all three nonideal MHD effects: Ohmic dissipation, the Hall effect, and ambipolar diffusion, which are characterised by their respective diffusivity coefficients η_{O}, η_{H}, and η_{A}. These coefficients are evaluated for a simple ionisation model in Sect. 2.4. Note that η_{H} and η_{A} are functions of B =  B , making the induction equation nonlinear.
It should be pointed out that the use of the isothermal shearing box approximation automatically removes large scale radial gradients and vertical thermal buoyancy effects.
The above equations admit a simple solution corresponding to an unperturbed Keplerian disc: (4)where H ≡ c_{s}/Ω and B_{z0} is a constant. We take this solution as our initial equilibrium state from which the MRI may develop. We often make use here of the velocity peculiar to the equilibrium shear flow, v = u + qΩxe_{y}. The strength of the mean vertical magnetic field is quantified by a modified plasma β_{0}, (5)This quantity equals the usual plasma β in the midplane times a sign function that describes the polarity of the magnetic field threading the disc. The Hall effect impacts the disc dynamics in a polaritydependent way (Wardle 1999; Balbus & Terquem 2001).
The boundary conditions are shearingperiodic in the x direction and periodic in the y direction. The vertical boundary conditions are:

vertical hydrostatic equilibrium for ρ;

zero vertical gradient for v_{x} and v_{y};

outflow for v_{z}; and

vertical field for B.
We have checked that the final condition can be replaced by a zerocurrent condition on the boundary without having any significant impact on our results.
2.2. Numerical method
We use a modified version of the finitevolume code Pluto (Mignone et al. 2007) to numerically integrate Eqs. (1)–(3). Ohmic dissipation is treated using the resistivity module that is included in the publicly available version of the code. We have introduced original implementations of ambipolar diffusion (Appendix B) and the Hall effect (Appendix A), which will be included in a future release of the code. Our version of Pluto implements a standard Godunov method using secondorder accurate spatial reconstruction with a monotonized central flux limiter, a modified HLL Riemann solver (described in Appendix A), and the constrained transport method (CT) of Evans & Hawley (1988), which maintains ∇·B = 0 to machine precision. Facecentered electromotive forces (EMFs) computed by the Riemann solver are interpolated to cell corners using arithmetic averaging. Flock et al. (2010) demonstrated that such averaging can affect MRI growth rates when combined with HLLD or Roe Riemann solvers. However, since we use an HLL solver, we do not anticipate any numerical artefacts from the arithmetic averaging. Indeed, we show that we recover the linear properties of the MRI in Sect. 3.1. The equations, including diffusion terms, are advanced in time explicitly using a secondorderaccurate RungeKutta scheme^{2}.
We choose our units such that Ω = c_{s} = ρ_{0} = 1, which implies H = 1. Unless otherwise stated, the total integration time is set to τ = 1000 Ω^{1}, which corresponds to ≈160 local orbital periods. Time averages are computed from t = 200 Ω^{1} to τ to avoid transients resulting from the initial conditions. We consider a box of size 4 × 8 × 12 and resolution 64 × 64 × 192. While this resolution (16 points per H) is poorer than that used in most contemporary ideal MHD simulations of the MRI, we find it to be more than enough to capture the relevant physical processes in our simulations. This is because the large magnetic diffusivities involved in our calculations render much of the flow laminar. Besides, at this resolution one run which includes all of the nonideal MHD effects at R_{0} = 1 au requires 70 000 core hours on the best available CPUs. Doubling this resolution would lead to an increase in the CPU time by a factor of 32, making a systematic exploration of the parameter space that is impractical, given our computational resources. To test convergence, we repeated run 1OHA5 with a doubled resolution for τ = 100 Ω^{1}, finding a quantitative difference of <10% in the transport diagnostics between 80 Ω^{1} and 100 Ω^{1} (see Table 1). This increases our confidence in the results. Unless otherwise stated, the initial conditions are the equilibrium (4) to which we add white noise on v_{x} with a typical rms amplitude of 0.057.
The evolution of the MRI in our simulations generically produces outflows, and our box continuously loses mass (and horizontal magnetic flux). To mimic a steady state with a global accretion flow replenishing the disc, we use a mass renormalisation procedure: at each time step, we multiply ρ in each cell by a constant factor such that the total mass in the box is conserved (Ogilvie 2012). This procedure implies that the total momentum in the box is not conserved, since the velocity field is kept constant during the renormalisation.
2.3. Diagnostics
We use several diagnostics to quantify the efficiency of angularmomentum transport and the production of outflows. These rely on two averaging procedures: The amount of transport is determined by the Reynolds stress R_{ij} = ρv_{i}v_{j} and the Maxwell stress M_{ij} = −B_{i}B_{j}. This allows us to define an effective (radial) transport parameter, (8)Outflows are characterised by the massloss rate, (9)and the (dimensionless) surface magnetic stress (10)which is evaluated at the top (z_{t} = 4.5H) and bottom (z_{b} = −4.5H) of the disc surface. These (arbitrary) locations are chosen to be sufficiently far from the vertical boundaries and yet far enough from the disc midplane to always be in the magnetically dominated region of the outflow (i.e., where β < 1).
List of runs discussed in this paper.
These diagnostics can be directly related to largescale quantities, such as the accretion rate Ṁ_{acc}. Assuming a 1 M_{⊙} central object, smooth surface density and temperature profiles, and neglecting prefactors in Fromang et al. (2013), where ϵ ≡ H/R. The above expression is made up of two contributions to the accretion rate. The first term in the square brackets leads to the traditional viscous αdisc accretion rate of Shakura & Sunyaev (1973), while the second term accounts for the stress from a magnetically driven outflow. For the latter, we have taken the stress to be opposite on both sides of the disc. We show later that this is not necessarily supported by our simulations; the contribution of this term should thus be considered an upper bound.
2.4. Ionisation profiles
Our disc model is the minimum mass solar nebula (MMSN; Hayashi 1981), whose column density and temperature profiles, are sampled at R_{0} = 1, 5, and 10 au. To compute the ionisation profiles, we assume hydrostatic balance in the vertical direction and consider the following contributions to the ionisation rate ζ:

Xray ionisation due to 3 keV photons (Igea & Glassgold 1999; Bai & Goodman 2009, see their Eq. (21));

cosmicray ionisation with ζ_{cr} = ζ_{0}exp(−Σ /96 g cm^{2}) s^{1} (e.g. Umebayashi & Nakano 1981) and ζ_{0} = 10^{16} s^{1}, corresponding to the cosmicray ionisation rate that is observed in the direction of ζ Persei (McCall et al. 2003) ;

radioactive decay with ζ_{rad} = 10^{19} s^{1} (Umebayashi & Nakano 2009).
The ionisation fraction x_{e} is obtained by balancing these ionisation sources with dissociative recombination in a metal and dustfree environment (Gammie 1996; Fromang et al. 2002): (12)where α_{dr} = 3 × 10^{6}T^{−1/2} cm^{2} s^{1} is the dissociative recombination rate coefficient for molecular ions.
Fig. 1
Ionisation fraction x_{e} versus height z at 1, 5, and 10 au. 

Open with DEXTER 
The final term in Eq. (12) represents the contribution from FUV radiation, which is known to almost fully ionise carbon and sulfur and lead to ionisation fractions x_{FUV} ~ 10^{5}–10^{4} for a penetration depth ~10^{2} g cm^{2} (PerezBecker & Chiang 2011). We use (13)as a rough estimate for the FUV ionisation fraction, which captures the essence of the FUV ionisation front. A similar approach has been used by Bai & Stone (2013b), replacing the exponential in Eq. (13) by a step function. The resulting ionisation profiles at R = 1, 5, and 10 au are given in Fig. 1. For simplicity, we assume that these profiles are fixed throughout the course of our simulations. This simplification is not likely to be satisfied in actual protoplanetary disks, since changes in the density profile, as well as turbulent mixing in the disc, could lead to local changes in the ionisation fraction. A more sophisticated treatment, in which the chemistry evolves alongside the dynamics, will be employed in a future publication.
These ionisation profiles, along with a choice of β_{0}, determine the diffusivity coefficients for Ohmic, Hall, and ambipolar diffusion (in the absence of dust grains; Balbus & Terquem 2001; Wardle 2007): where is the electronneutral collision rate (Draine et al. 1983), ρ_{i} is the ion mass density, γ_{i} = ⟨ σv ⟩ _{i}/ (m_{n} + m_{i}), and is the ionneutral collision rate (Draine 2011); m_{n} and m_{i} are the average masses of the neutrals and ions, respectively. Due to our normalisation of the magnetic field in Eqs. (2) and (3), factors of and 4π appear in our expressions for the Hall and ambipolar diffusivities. Note that the upper layers of protoplanetary discs at 1 au are likely to be strongly ionised and heated by Xrays to temperatures up to ~8000 K (Aresu et al. 2011). Our simple isothermal model thus breaks down for column densities ≲10^{3}–10^{2} g cm^{2}. In response, we arbitrarily multiply our diffusivities by a constant factor of exp(−Σ/0.01 g cm^{2}) to mimic a hot and fully ionised gas in the upper layers.
Fig. 2
Ohmic (Λ; Eq. (15)), Hall (Ha; Eq. (16)), and ambipolar (Am; Eq. (17)) Elsasser numbers in the initial state versus height z at 1 au (top) and 10 au (bottom), assuming a constant vertical magnetic field with β_{0} = 10^{5}. 

Open with DEXTER 
Fig. 3
Hall lengthscale ℓ_{H} (Eq. (18)) versus height z at 1, 5, and 10 au. The dotted line denotes the threshold above which saturation via zonal magnetic fields occurs in unstratified simulations with net vertical magnetic flux (KL13). 

Open with DEXTER 
Introducing the Alfvén speed , the diffusivities (12) can be cast in terms of dimensionless Elsasser numbers: Example profiles of these Elsasser numbers are given in Fig. 2; note that they evolve throughout the simulation only via changes in density and magneticfield strength. These profiles are similar to those used by Bai & Stone (2013b)^{3}. Another useful way of quantifying Hall diffusion, which is independent of the magneticfield strength, is the Hall lengthscale: (18)When ℓ_{H} ≳ 0.2H, KL13 found that unstratified simulations of the HallMRI with net vertical magnetic flux saturate by the production of strong axisymmetric (“zonal”) magnetic fields (see also Appendix A). In Fig. 3, we present ℓ_{H} as a function of z for our ionisation profiles, demonstrating that most of the disc midplane is in the putative zonalfield regime.
The nonideal MHD terms introduce second derivatives of the magnetic field into the induction equation. This often yields a tight constraint on the maximum allowed stable timestep in our simulations. To circumvent this problem, we follow Bai & Stone (2013b) by introducing a capping procedure on the magnetic diffusivities η_{O,H,A}. Whenever one of these diffusivities becomes larger than η_{cap} = 10 ΩH^{2}, we automatically set its value to η_{cap}. We have briefly explored the impact of the cap value on our results: increasing η_{cap} for η_{O} and η_{A} does not quantitatively impact our results at 1 au. However, increasing η_{cap} for η_{H} by a factor of 4 tends to increase α significantly (by ~50% at 1 au). Therefore, our stress values for ΩB_{z0} > 0 should be understood as lower bounds when the Hall effect is included.
3. Results
3.1. Linear evolution of stratified discs subject to Ohmic and Hall diffusion
We first concentrate on the influence of the Hall effect on the linear evolution of an MRIunstable, stratified disc. We set η_{A} = 0 to isolate the effects of Hall diffusion but include Ohmic dissipation. We use a setup that is identical to run 1OH5, except that we seed the instability with very small perturbations to the velocity field (RMS = 10^{6}) so as to produce a clean linear phase. We compare the vertical profile of the resulting azimuthal magnetic field, B_{y}(z), at t = 10 Ω^{1} with that predicted by linear theory, where the latter is calculated with a pseudospectral method similar to that used by Latter et al. (2010). Each linear eigenmode is characterised by its growth rate γ and its symmetry with respect to the midplane, σ = ±1: (19)The fastestgrowing eigenmodes, their calculated growth rates, and their symmetries are listed in Table 2. For the ionisation profiles, we consider that the n = 1 and n = 2 eigenmodes have very similar growth rates, making the identification of each individual mode difficult. To isolate each mode, we decompose B_{y}(z) from the numerical simulation into symmetric (σ = 1) and antisymmetric (σ = −1) parts, each of which can then be compared to the two fastestgrowing eigenmodes obtained by linear theory (Fig. 4). We find that the most unstable modes are accurately captured by our implementation of the Hall effect at this resolution. Moreover, the measured growth rate γ_{num} = 0.726 Ω in our numerical simulation matches the theoretical value for the fastest growing mode to less than one percent.
The above calculation assumed an isothermal fluid, for which buoyancydriven modes are absent. We refer the reader to Urpin & Rüdiger (2005) for a local, linear calculation that couples vertical buoyancy to the Hall effect.
Fig. 4
Comparison of B_{y}(z) from the simulation (red circles) and from the linear calculation (black line), which include Ohmic and Hall diffusion computed at 1 au with β_{0} = 10^{5}. Top: n = 1 eigenmode and symmetric component of B_{y}(z). Bottom: n = 2 eigenmode and antisymmetric component of B_{y}(z). 

Open with DEXTER 
3.2. Nonlinear evolution of stratified discs subject to Ohmic and Hall diffusion
3.2.1. Fiducial runs
In Fig. 5, we present space–time diagrams of ⟨ B_{y} ⟩ for runs 1O5 (top) and 1OH5 (bottom). The Ohmic run (1O5) presents an alternating pattern in the upper layers commonly referred to as a “butterfly diagram”, which was first described in the context of magnetorotational turbulence by Brandenburg et al. (1995) and Stone et al. (1996). When the Hall effect is included (1OH5), this butterfly pattern is replaced by a strong quasisteady azimuthal magnetic field that fills the disc midplane. This strong azimuthal field is accompanied by a weaker radial field (not shown), which feeds the azimuthal component via shear.
The Ohmic run exhibits the usual layered accretion model (Gammie 1996; Fleming & Stone 2003) with a very weak Maxwell stress in the disc midplane and a turbulent layer in the ionised atmosphere (Fig. 6top). When the Hall effect is included, we find that the magnetic stress is increased by almost two orders of magnitude (from α = 2.5 × 10^{3} to α = 4.5 × 10^{1}) and is clearly located in the Halldominated region (−3H ≲ z ≲ 3H; Fig. 6bottom). This stress is not a “turbulent” stress in the usual sense but rather is a consequence of the largescale magnetic structure shown in Fig. 5. In the upper atmosphere ( z  > 4), we find a more classical turbulent layer qualitatively similar to the active layer in run 1O5.
Fig. 5
Space–time diagram of the horizontally averaged azimuthal magnetic field, ⟨ B_{y} ⟩, in the Ohmic (1O5; top) and OhmicHall (1OH5; bottom) runs. 

Open with DEXTER 
The presence of a vertical magnetic flux is known to trigger outflows in shearingbox simulations (Suzuki & Inutsuka 2009; Lesur et al. 2013; Fromang et al. 2013; Bai & Stone 2013a); our simulations exhibit outflows as well. We find that these outflows are sensitive to the presence of the Hall effect: massloss rates increase by more than two orders of magnitude. However, the outflows do not have any welldefined orientation: they can be directed toward either positive or negative x, which is a phenomenon previously observed in idealMHD simulations presented by Fromang et al. (2013) and Bai & Stone (2013a). This is one consequence of the neglect of curvature terms in the shearingbox approximation. Because of this asymmetry, we find no average stress exerted at the disc surface by the outflow.
Fig. 6
Space–time diagram of the logarithm of the horizontally averaged Maxwell stress, log ⟨ −B_{x}B_{y} ⟩, in the Ohmic (1O5; top) and OhmicHall (IOH5; bottom) runs. 

Open with DEXTER 
The very large massloss rate found in the OhmicHall run can be explained by the modification of the vertical equilibrium, which is no longer hydrostatic. As revealed by Fig. 7, the pressure associated with the strong azimuthal magnetic field puffs up the disc, thereby increasing the amount of mass in the disc atmosphere where the outflow originates (near z = 4.5H). On the other hand, the vertical velocity is not significantly affected by the Hall term. This explains the increase in the rate of mass loss by two orders of magnitude between runs 1O5 and 1OH5.
We note that a similar thickening of the disc atmosphere has been found by Hirose & Turner (2011) using resistive MHD simulations including radiative transfer. This thickening has been used to explain the infrared excess in Herbig stars, which cannot be explained by a disc in vertical hydrostatic equilibrium (Turner et al. 2014a). Our model presents an alternative scenario, in which the strong magnetic pressure driven by the Hallamplified azimuthal magnetic field dominates the total pressure in the disc midplane.
Fig. 7
Horizontally averaged density profile (top) and vertical velocity (bottom) versus height z, averaged over ≈140 orbits, from runs 1O5 (solid line) and 1OH5 (dashed line). 

Open with DEXTER 
3.2.2. Origin of the Halldriven azimuthal field
Fig. 8
Contributions in run 1OH5 to the righthand side of the horizontally averaged induction equation deduced from Eqs. (20a; top) and (20b; bottom). 

Open with DEXTER 
The presence of a very strong mean azimuthal magnetic field in the disc midplane of our HallMHD simulations is surprising. To understand its origin, we decompose the horizontal components of the horizontally averaged induction equation (Eq. (3)) into several pieces, of which each is identified with a physical process: We have identified five types of terms:

outflow:
describes the vertical transport of horizontal flux tubes.It measures the rate at which the outflow evacuates horizontalmagneticfield lines from the disc midplane;

quantifies the horizontal bending of vertical field lines due to horizontal motions;

shear:
measures the creation of azimuthal field from radial field by the mean Keplerian shear;

Hall:
denotes the contributions from the Hall effect; and

Ohmic:
denotes the contributions from Ohmic dissipation.
The saturated state of run 1OH5 is shown in Fig. 8 with each of these terms identified^{4}. We focus on the inner disc region  z  < 3H, where Hall and Ohmic diffusion dominate the dynamics.
Fig. 9
Space–time evolution of the logarithm of the horizontallyaveraged magnetic stress, log ⟨ M_{xy} ⟩, in the Ohmic (1O5; top), Ohmicambipolar (1OA5; middle), and OhmicambipolarHall (1OHA5; bottom) runs. 

Open with DEXTER 
We first describe the contributions of the righthand side of the xcomponent of the induction equation (Fig. 8top). We find that the stretching term is totally negligible and that the equilibrium is dominated by a balance between the Hall term, which amplifies the radial field (note that B_{x} < 0 in the disc midplane), and the outflow term, which carries B_{x} out of the midplane. This quasiequilibrium is different than that found in idealMHD simulations, in which the driving term is due to stretching. On the other hand, the ycomponent of the induction equation (Fig. 8 bottom) is a resistive MHD quasiequilibrium, in which the shear production of the azimuthal field from the radial field is balanced by the outflow term in the atmosphere and the Ohmic term in the disc midplane. The Hall terms are completely negligible here.
This analysis allows us to present a selfconsistent picture of run 1OH5. We first note that Moreover, since the fluctuations δB_{z} are small compared to the mean vertical field B_{z0}, we have ∂_{z} ⟨ B_{y}B_{z} ⟩ ≃ B_{z0}∂_{z} ⟨ B_{y} ⟩ and ∂_{z} ⟨ B_{x}B_{z} ⟩ ≃ B_{z0}∂_{z} ⟨ B_{x} ⟩. Retaining only the driving terms in Eq. (18) and neglecting the outflow, we find These reduced equations for ⟨ B_{x} ⟩ and ⟨ B_{y} ⟩ describe a nonlocal version of the “Hallshear” instability (cf. Kunz 2008, Eq. (46)): Keplerian shear generates an azimuthal magneticfield component from a radial one, and the Hall effect conservatively reorients this azimuthal field back into the radial direction. When qΩB_{z0} > 0, this feedback loop can lead to growth. Note that the Hallshear instability is not a version of the MRI and does not rely on the Coriolis force (though it can be obtained from the HallMRI dispersion relation by taking the highly diffusive limit – see Rüdiger & Kitchatinov 2005 and Wardle & Salmeron 2012). It is, however, sensitive to the polarity of the vertical magnetic field. This is discussed more extensively in Sect. 3.3.2.
It is possible to understand the physical content of Eq. (20) by observing that , where v_{d} is the drift velocity of the electrons relative to the ions. Hence, the term on the righthand side of (21a) is simply the magnetic stretching term due to electron motion. In protoplanetary discs, the number of free electrons can be very small, so that the drift speed can be large for a given current density. As a result, the usual stretching term of ideal MHD caused by the bulk motion of the gas is replaced by a stretching term due to the electron motion. The electron drift produces a radial field from the vertical mean field (Eq. (21a)), which is then sheared by the Keplerian flow to produce an azimuthal field (Eq. (21b)). It is this azimuthal field that sustains the drift through Ampère’s law.
In our simulations, the nonlinear saturated state results from a balance between this Hallshear instability, an outflow which transports the magnetic energy away from z ~ H, and Ohmic diffusion which diffuses the field from the midplane to the base of the outflow (see Fig. 8).
3.3. Nonlinear evolution of stratified discs subject to Ohmic, Hall, and ambipolar diffusion
3.3.1. Fiducial runs
With an understanding of how a strong Hall effect impacts resistive protoplanetary discs, we now turn to simulations that also include ambipolar diffusion. A comparison of the magnetic stresses measured in the Ohmic (1O5), Ohmicambipolar (1OA5), and OhmicambipolarHall (1OAH5) runs at R_{0} = 1 au is presented in Fig. 9. The first two simulations are quantitatively similar to runs Ob5 and OAb5 of Bai & Stone (2013b), and we recover their main statistical properties. However, the inclusion of Hall diffusion dramatically changes the picture. Similar to the behaviour discussed in Sect. 3.2, we find a large vertical band of strong Maxwell stress in the disc midplane, which substantially increases the transport by almost two orders of magnitude.
Fig. 10
Horizontally averaged Maxwell (top) and Reynolds (bottom) stress versus height z, averaged over 700 Ω^{1}, from runs 1O5 (dashdotted), 1OA5 (dashed), and 1OAH5 (solid). 

Open with DEXTER 
The resulting stress profiles averaged over the entire simulation from t = 300 Ω^{1} are presented in Fig. 10. The addition of ambipolar diffusion tends to weaken the activity in the surface layers (cf. Bai & Stone 2013b). The addition of the Hall effect does not alter this picture, and we find that the surface activity is even more reduced in the OAH run. However, magnetic activity is greatly enhanced in the midplane, where stresses as high as M_{xy} = 3 × 10^{2} are found. The Reynolds stress is not significantly affected by the Hall effect. This implies that the flow is not “turbulent” in the usual sense.
The presence of a strong Hall effect in the midplane also affects the vertical structure of the disc. We find a turbulent plasma beta in the disc midplane in run 1OAH5 (Fig. 11). The hydrostatic equilibrium is significantly affected by the magnetic structure and the disc scale height increases, as in Sect. 3.2. Outflows are also affected in this new quasiequilibrium state, producing massloss rates a factor of 20 larger than those obtained in the Ohmicambipolar simulations.
Fig. 11
Horizontally averaged turbulent plasma β versus height z, averaged over 700 Ω^{1}, from runs 1O5 (dashdotted), 1OA5 (dashed), and 1OAH5 (solid). 

Open with DEXTER 
3.3.2. Parameter study
Field polarity:
it is well known that the linear HallMRI is sensitive to the field polarity (Wardle 1999; Balbus & Terquem 2001). To test whether or not this result carries over to the nonlinear regime, we have performed a simulation identical to run 1OHA5 but with a reversed mean field: (run 1OHAmB). In this case, we find that the level of transport, the surface stress, and the massloss rate to be smaller than in the run without the Hall effect (run 1OA5). Therefore, when , the Hall effect weakens the efficiency of the MRI and the strength of the outflow. This sensitivity to the field polarity was already pointed out in Sect. 3.2.2 as a natural consequence of the fact that the radial drift of electrons drives the nonlinearly saturated state. A similar trend was found by Sano & Stone (2002) in unstratified shearing boxes, though the Hall effect was much weaker in their case and the result less extreme.
Fig. 12
Horizontally and temporally averaged Maxwell (top) and Reynolds (bottom) stress versus height z for runs 1O3 (dashdotted), 1OA3 (dashed), and 1OAH3 (solid). 

Open with DEXTER 
Zero verticalnetflux configuration:
the zero verticalnetflux configuration has been largely explored in the past, with both unstratified (Hawley et al. 1996; Fromang & Papaloizou 2007) and stratified simulations (Brandenburg et al. 1995; Stone et al. 1996; Fleming & Stone 2003; Simon et al. 2012, 2013b). We have only explored one such situation including all three nonideal effects (run 1OHAznf). This simulation is identical to run 1OHA5 except that the initial magneticfield configuration is given by B_{z} = B_{z0}sin(2πx/L_{x}), where B_{z0} is the initial field strength from run 1OHA5. This initial condition rapidly collapses into a quiet state with virtually no transport and no outflow (see Table 1).
Mean field strength:
the runs 1xxx5 have an initially weak mean vertical magnetic field in the disc midplane. To test how this initial condition might influence the subsequent evolution, we have run a set of simulations with β_{0} = 10^{3} (runs 1O3, 1OA3, 1OAH3). The resulting profiles are presented in Fig. 12, which may be directly compared to the weak field case (Fig. 10). An initially stronger mean field significantly increases the transport in the active surface layers of the Ohmic run, as anticipated using results from idealMHD simulations (cf. Hawley et al. 1995). A strong Maxwell stress also appears in the disc midplane, which is most likely due to horizontal field lines diffusing down from the active layers (Turner & Sano 2008). Taken as a whole, these three simulations exhibit Maxwell stresses in the midplane that are roughly one order of magnitude larger than those found in run 1OHA5; this scaling is qualitatively similar to the scaling found in ideal MHD simulations (Hawley et al. 1995). Note, however, that α (the vertically integrated stress) does not increase that steeply (see Table 1); the stress is less vertically distributed for β_{0} = 10^{3}. Fitting our results with a simple power law, we obtain .
The outflow is also affected by the mean field strength. We recover approximately the scalings and found by Bai & Stone (2013b) in runs that include only Ohmic and ambipolar diffusion. When the Hall effect is included, the scaling for still holds but the scaling for Ṁ_{outflow} is shallower with .
Distance from the central protostellar object:
the ionisation profile depends strongly upon the radial location in the disc, R_{0} (Sect. 2.4). At larger distances, the gas density is lower, and not only do cosmic rays penetrate deeper into the disc interior, but also the recombination time is longer. Both effects lead to an increase in the ionisation fraction. For this reason, the saturation properties of the MRI vary with radius. We have performed three simulations with β_{0} = 10^{5}, varying R_{0} from 1 au to 10 au. The resulting Maxwell stress profiles are presented in Fig. 13^{5}.
We find that the Maxwell stress in the disc midplane decreases as one moves to larger radii, since Ha increases and approaches Am. Note that the outflow surface stress and the massloss rate decrease in a similar fashion (see Table 1). By contrast, the Maxwell stress tends to increase in the surface layers, contributing roughly to 1/5 of the vertically integrated stress at 10 au. We find that this component of the stress at 5 and 10 au is due to turbulence in the surface layer, in a way that is qualitatively similar to that found in run 1OH5. Overall, this trend indicates that, as we move outward, we approach the point at which the stress contribution due to the surface layers dominates that of the vertically integrated disc. This corresponds to the classic layered accretion picture, which is a limit that corresponds to “region III” described by Bai (2013) and studied extensively by Simon et al. (2013a,b) in situations where ambipolar diffusion suppresses the MRI in the disc midplane.
Fig. 13
Horizontally averaged Maxwell stress versus height z, averaged over 700 Ω^{1}, from runs 1OAH5 (dashdotted), 5OAH5 (dashed), and 10OAH5 (plain line). 

Open with DEXTER 
Symmetries:
the saturated states of all the OAH runs exhibit the same σ = 1 symmetry (see Eq. (19)). However, this symmetry implies that the outflow points inwards on one side of the disc and outwards on the other side. A more physical outflow solution has σ = −1. It has been shown by Bai & Stone (2013b) that ambipolardiffusion–dominated simulations can exhibit solutions having outflows with the proper symmetry at  z  ≳ 3.5H if particular initial conditions are chosen. These solutions exhibit a strong offmidplane current layer, which was found to be stable by Bai & Stone (2013b). We have tried to find similar solutions including the Hall effect but have failed: we have only found σ = 1 solutions when Hall diffusion is included. However, we are able to find longlived σ = −1 solutions when only Ohmic and ambipolar diffusion are included (see run 1OA5e in Fig. 14). These solutions start to develop a strong offmidplane current layer (around t ~ 1000 Ω^{1}), which is very similar to the one presented by Bai & Stone (2013b). However, we find that this layer is eventually ejected on longer timescales, and the system ultimately relaxes into a σ = 1 symmetry. This casts doubt upon the stability and even the applicability of these outflow solutions in actual protoplanetary discs.
Fig. 14
Space–time evolution of the horizontally averaged azimuthal magnetic field, ⟨ B_{y} ⟩, as a function of time in run 1OA5e. The run starts with odd symmetry for the azimuthal field. At t ≃ 1000 Ω^{1} a current layer forms and is ejected at t = 1410 Ω^{1}, leaving an azimuthal field with even symmetry. 

Open with DEXTER 
3.3.3. Magnetic selforganisation in Halldominated magnetorotational turbulence
KL13 have shown that the saturated state favoured by Halldominated magnetorotational turbulence in unstratified shearing boxes with net vertical magnetic flux is characterised by a strong axisymmetric (“zonal”) magnetic field and a vanishingly low level of turbulent transport. Remarkably, none of the simulations presented in this paper exhibit this behaviour. This disparity is not caused by the different numerical algorithms employed. Indeed, we were able to reproduce their results with our fiducial 16 points per H (see Sect. A.2.3). Instead, the difference is due to the strong azimuthal field that is naturally produced in our stratified simulations. In unstratified simulations, the net magnetic flux is conserved and, if there is initially no net azimuthal flux, none will be generated (⟨ B_{y} ⟩ = 0). In contrast, the outflow boundary conditions imposed at the top and bottom of the stratified shearing box allows a net azimuthal field to develop: azimuthal flux of opposite polarity is ejected during the generation of the outflow. As a result, stratified simulations with the Hall effect can produce a (very large) mean azimuthal field relative to the mean vertical field, typically with ⟨ B_{y} ⟩ ~ 200 ⟨ B_{z} ⟩. As pointed out by KL13, this magnetic configuration (⟨ B_{y} ⟩ ≫ ⟨ B_{z} ⟩) does not saturate via the production of zonal magnetic fields.
4. Influence of dust grains
Dust grains comprise ~1% by mass of protostellar cores and, by extension, protoplanetary discs (Hayashi 1981). In the latter, micronsized grains can dramatically increase the rate of magnetic diffusion, mainly because of their propensity to soak up free charges at high densities. This tends to decouple the gas from the magnetic field (Semenov et al. 2004; Wardle 2007; Bai 2011). Since these models usually involve complex chemical networks and various grain distributions, it is difficult to obtain a good physical intuition for the effect grains can have on the physics of accretion disc turbulence. To clarify the situation, we discuss the qualitative effects of dust grains on the chemistry and dynamics of a protoplanetary disc. Our discussion focuses on the disc midplane, where the ionisation rate is very low (typically ~).
4.1. Diffusivity tensor
Dust grains preferentially capture free electrons, since electrons have less inertia than the ions. This process occurs quickly, so that grains are usually the dominant negative charge carriers when they are well mixed with the gas. Since grains carry negative charges, they also tend to increase the effective recombination rate with ions, acting as a catalyst for recombination. This decreases the total ionisation fraction by one or two orders of magnitude, depending upon the grain size. A typical situation is presented in Bai (2011, fig. 1), where we clearly see that the presence of 0.1μmsized grains decreases the ionisation fraction and makes singly charged grains the dominant charge carriers.
In addition to the modification of the ionisation equilibrium, dust grains also have an impact on the gas dynamics. From a plasma point of view, the presence of charged grains indicates that some charge carriers are much heavier than electrons and ions. This suggests that the average mobility of charge carriers and their coupling time with neutral H_{2} are drastically reduced. The diffusivities associated to the three nonideal MHD effects (η_{O}, η_{H}, η_{A}) are therefore significantly altered.
To compute the diffusivities in such a plasma, one should follow the derivation of the complete diffusivity tensor of Wardle (2007), which takes into account all charged species. This calculation is beyond the scope of this paper, but we can rely on the work of Salmeron & Wardle (2008) to investigate the impact of dust grains on protoplanetary disc midplanes. Consider their Fig. 1, which presents diffusivity profiles for various grain contents. Without grains, we observe that the dominant nonideal effect in the midplane is Hall at 5 and 10 au, which matches our own diffusivity profiles (see Sect. 2.4). When 1 grains are introduced, the respective ratios of Ohmic, Hall, and ambipolar diffusion are largely unaffected in the midplane, although each of these diffusivities are increased by roughly three orders of magnitude. The most dramatic modification comes from grains: at 10 au, those authors found that ambipolar diffusion dominates in the disc midplane, followed by the Hall effect at about a factor of 10 smaller. At 5 au, all three diffusivities are comparable in the midplane, although their absolute values are increased by five to six orders of magnitude.
4.2. Vertical distribution
Despite the presence of a strong magnetic torque and the consequent large rates of angularmomentum transport observed in our Halldominated simulations, the flow appears to be predominantly laminar. One might then expect any population of large dust grains to slowly settle into the disc midplane and thereby affect the diffusivity tensor. However, the driven outflows may lift up these grains and, in doing so, mimic the role traditionally ascribed to turbulent stirring. It is therefore natural to ask whether or not grains are expected to sediment under the conditions found in our simulations.
To this end, let us consider the vertical equilibrium for dust grains obtained by balancing vertical gravity with gas drag: where g_{z} is the vertical component of gravity, v_{z} is the vertical gas velocity, and ⟨ σv ⟩ _{g} is the rate of momentum exchange of grains with neutral gas molecules. We approximate ⟨ σv ⟩ _{g} ~ a^{2}c_{s} and m_{g} ~ a^{3}ρ_{S}, where ρ_{S} is the grain material density and a is the grain radius. Using g_{z} = Ω^{2}z, we find to be the equilibrium height reached by the grains when they are dragged by the outflow. The continuity equation implies that ρv_{z} is conserved along z, defining the massoutflow rate. Assuming ρ_{S} = 1 g cm^{3} and using Eq. (9) to quantify the massloss rate, we obtain (22)This estimate suggests that dust grains with sizes up to 1 mm (30 μm) can be lifted up at 1 au (10 au) for typical outflow massloss rates found in our simulations (i.e. Ṁ_{outflow} ≳ 10^{4}).
4.3. Impact on this work
Overall, it appears that grains should not have too much of an impact on our results at 5 and 10 au: the Hall effect still dominates in the midplane, and the increase in the diffusivities by three orders of magnitude leads to Elsasser numbers close to the midplane values we adopted at R_{0} = 1 au. Hence, the role played by the Hall effect is likely to be qualitatively comparable to what we presented in Sect. 3.3 though quantitatively different. To test this conjecture, we performed two simulations at 10 au in which the diffusivities are artificially increased by a factor of 1000 throughout the vertical extent of the disc (runs 10OAgr and 10OHAgr). Even in this worstcase scenario, the Hall effect still significantly enhances the magnetic stress in the disc (see α values in Table 1).
For 0.1 μm grains, the impact is less obvious. The Hall effect can still be a major player at 5 au, but the dramatic increase of the diffusivities brings us into new territory, for which we have no numerical simulations to guide our intuition. This regime is quite difficult to study numerically, since the diffusivities are much larger than the cap values we used in this work (see Sect. 2.4). Increasing the cap values to the predicted diffusivities would lead to extremely small time steps, making explicit numerical simulations, such as ours, impractical at this time.
5. Summary
In this paper, we have explored the linear and nonlinear behaviour of poorly ionised, vertically stratified protoplanetary discs. For the first time, all three relevant nonideal MHD effects – Ohmic dissipation, ambipolar diffusion, and the Hall effect – are selfconsistently included. To accomplish this feat, we have implemented an original formulation of ambipolar diffusion and the Hall effect in the finitevolume code Pluto.
Our results demonstrate that none of these effects can be safely neglected at distances ~1–10 au from the central protostellar object. En route to this conclusion, we have confirmed previous work on the effects of Ohmic dissipation and ambipolar diffusion on magnetorotational turbulence in stratified discs: Ohmic diffusion quenches the MRI in the disc midplane, and ambipolar diffusion suppresses turbulence in the surface layers, making the traditional viscous disc model insufficient to explain the observed accretion rates (α ≈ 6 × 10^{4} at 1 au). Our principal finding is that the Hall effect alters this picture dramatically. A strong azimuthal magnetic field is produced in the disc midplane, which generates a vertical magnetic pressure gradient strong enough to substantially increase the vertical scale height of the disc. This azimuthal field correlates with a weak radial field to produce very efficient angularmomentum transport with α ~ 10^{2}. These values are compatible with accretion rates Ṁ ≲ 10^{7} M_{⊙} yr^{1} (Eq. (11a)) without a need to invoke external stresses generated by disc outflows (as in Bai & Stone 2013b).
The MRIdriven outflows are also modified by the Hall effect. The massloss rates increase by a factor ≲20 over those found in runs with only Ohmic dissipation and ambipolar diffusion. The surface stresses also increase by a factor of ~4 at 1 au. However, we caution that these numbers may be inappropriate for global disc models, since outflows (and in particular the massloss rates) are affected by the vertical boundary conditions in the shearingbox approximation (Fromang et al. 2013; Lesur et al. 2013). Shearingbox outflows should be treated with caution.
Surprisingly, our stratified simulations do not produce the axisymmetric (“zonal”) magneticfield structures and the attendant steep reduction in turbulent transport found by KL13 in unstratified simulations with a net vertical magnetic flux. This is despite our simulations being in the right regime (ℓ_{H} ≳ 0.2H). We speculate that this difference is due to the strong azimuthal magnetic field produced in the disc midplane for the stratified case, which “hides” the vertical field from the dynamics. This is related to our choice of vertical boundary conditions, which permits the evacuation of magnetic flux tubes from the computational domain and the production of a net azimuthal field (despite none being present initially). This buildup of azimuthal flux is not possible in unstratified shearing boxes. It is suggestive that KL13 did not find the zonalfield route to saturation in unstratified shearing boxes with ⟨ B_{y} ⟩ ≫ ⟨ B_{z} ⟩. We therefore conjecture that both the global topology of the magnetic field threading the accretion disc and the radial and vertical ionisation profiles may play a role in determining which saturated state is favoured.
It is wellknown that the influence of the Hall effect on the linear stability of the disc depends upon the field polarity Ω·B (Wardle 1999; Balbus & Terquem 2001), which is a result we have found to hold true in the nonlinear regime as well (see also Sano & Stone 2002). Configurations with Ω·B > 0 show enhanced transport, while configurations with Ω·B < 0 show reduced transport (viz. α ~ 10^{4} and no significant outflow). This suggests that one could have disconnected regions of the disc with different field polarities, which are either actively accreting or relatively quiet and evolve according to the global topology and the longterm evolution of the largescale magnetic field.
The relatively large magnetic stresses generated by the Hall effect are not associated with turbulent fluctuations. Rather, they are associated with largescale amplification of an azimuthal magnetic field via the shearing of a Hallinduced radial field (cf. Kunz 2008). Indeed, the Reynolds stresses are always very small in the midplane (with ⟨ R_{xy} ⟩ ≲ 10^{6}). This effect is distinct from the process found by Turner & Sano (2008), in which the magnetic field Ohmically diffuses from the upper active layers down into the midplane, where it exerts a largescale stress. One consequence of our result is that the Halldominated MRI “turbulence” may be unable to stir up dust grains that would otherwise slowly settle into the disc midplane under the action of vertical gravity. This could be a problem, since observations have implied that subμm grains are present in the upper atmospheres of protoplanetary discs (Pinte et al. 2008). On the other hand, such small grains could have been lifted up by a weak outflow (cf. Sect. 4.2).
Finally, using Salmeron & Wardle (2008) results, we have shown that sized grains should not have much of an impact on our conclusions: the Hall effect remains the dominant nonideal term in disc midplanes at 5 and 10 au with diffusivity coefficients comparable to that at 1 au without grains. The 0.1 sized grains have a more severe impact on the diffusivities, though Hall remains an important player at 5 au. We have not explored this regime numerically because of the cost of our explicit numerical scheme. We note, however, that there is no reason a priori to neglect the Hall effect at these distances.
While a more sophisticated treatment of the chemistry, in which dust grains are included and species abundances are timedependent, is necessary to make accurate quantitative predictions, the results presented here provide unequivocal evidence that Ohmic dissipation, ambipolar diffusion, and the Hall effect must all be taken into account to obtain a realistic description of angularmomentum transport in protoplanetary discs.
Hydrodynamical processes are an interesting alternative, despite serious theoretical difficulties; see Turner et al. (2014b) for a review.
Note that Bai & Stone (2013b) included a numerical cap to limit diffusion coefficients in their Fig. 1 for  z  < 2, which we do not include in our plot.
The Hall term is artificially large in Fig. 8top. This is because of the capping procedure used in the code (Sect. 2.4): highly magnetised cells see their Hall diffusivities artificially decreased, which is an effect we can only approximate when postprocessing horizontally averaged quantities.
Acknowledgments
G.L. is indebted to WingFai Thi for his many suggestions and clarifications concerning the radiative and chemical processes relevant to protoplanetary discs. G.L. also thanks Andrea Mignone and Gábor Tóth for fruitful discussions regarding the numerical implementation of the Hall effect during the Astronum 2013 conference. Support for G.L. was provided by the European Community via contract PCIG09GA2011294110. Support for M.W.K. was provided by NASA through Einstein Postdoctoral Fellowship Award Number PF1120084, issued by the Chandra Xray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS803060. Support for S.F. was provided by the European Research Council under the European Union’s Seventh Framework Programme (FP7/20072013)/ERC Grant agreement n° 258729. Most of the the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujfgrenoble.fr), which is supported by the RhôneAlpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This work was granted access to the HPC resources of IDRIS under allocation x2014042231 made by GENCI (Grand Equipement National de Calcul Intensif).
References
 Aresu, G., Kamp, I., Meijerink, R., et al. 2011, A&A, 526, A163 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Bai, X.N. 2011, ApJ, 739, 50 [NASA ADS] [CrossRef] (In the text)
 Bai, X.N. 2013, ApJ, 772, 96 [NASA ADS] [CrossRef] (In the text)
 Bai, X.N., & Goodman, J. 2009, ApJ, 701, 737 [NASA ADS] [CrossRef] (In the text)
 Bai, X.N., & Stone, J. M. 2013a, ApJ, 767, 30 [NASA ADS] [CrossRef] (In the text)
 Bai, X.N., & Stone, J. M. 2013b, ApJ, 769, 76 [NASA ADS] [CrossRef] (In the text)
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] (In the text)
 Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [NASA ADS] [CrossRef] (In the text)
 Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235 [NASA ADS] [CrossRef] (In the text)
 Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Blaes, O. M., & Balbus, S. A. 1994, ApJ, 421, 163 [NASA ADS] [CrossRef] (In the text)
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] (In the text)
 Brandenburg, A., Nordlund, Å, Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] (In the text)
 Chandrasekhar, S. 1960, Proc. Nat. Acad. Sci., 46, 253 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Desch, S. J. 2004, ApJ, 608, 509 [NASA ADS] [CrossRef] (In the text)
 Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press) (In the text)
 Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485 [NASA ADS] [CrossRef] (In the text)
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [NASA ADS] [CrossRef] (In the text)
 Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 [NASA ADS] [CrossRef] (In the text)
 Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464 [NASA ADS] [CrossRef] (In the text)
 Flock, M., Dzyurkevich, N., Klahr, H., & Mignone, A. 2010, A&A, 516, A26 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18 [NASA ADS] [CrossRef] (In the text)
 Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] (In the text)
 Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] (In the text)
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] (In the text)
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 464, 690 [NASA ADS] [CrossRef] (In the text)
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] (In the text)
 Hill, G. W. 1878, Am. J. Math., 1, 5 [CrossRef] (In the text)
 Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30 [NASA ADS] [CrossRef] (In the text)
 Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848 [NASA ADS] [CrossRef] (In the text)
 Ilgner, M., & Nelson, R. P. 2006a, A&A, 445, 205 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Ilgner, M., & Nelson, R. P. 2006b, A&A, 445, 223 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Ilgner, M., & Nelson, R. P. 2008, A&A, 483, 815 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Kunz, M. W. 2008, MNRAS, 385, 1494 [NASA ADS] [CrossRef] (In the text)
 Kunz, M. W., & Balbus, S. A. 2004, MNRAS, 348, 355 [NASA ADS] [CrossRef] (In the text)
 Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 [NASA ADS] [CrossRef] (In the text)
 Latter, H. N., Fromang, S., & Gressel, O. 2010, MNRAS, 406, 848 [NASA ADS] (In the text)
 Lesur, G., Ferreira, J., & Ogilvie, G. I. 2013, A&A, 550, A61 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2003, Nature, 422, 500 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] (In the text)
 Moll, R. 2012, A&A, 548, A76 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Norman, C., & Heyvaerts, J. 1985, A&A, 147, 247 [NASA ADS] (In the text)
 Ogilvie, G. I. 2012, MNRAS, 423, 1318 [NASA ADS] [CrossRef] (In the text)
 PerezBecker, D., & Chiang, E. 2011, ApJ, 735, 8 [NASA ADS] [CrossRef] (In the text)
 Pinte, C., Padgett, D. L., Ménard, F., et al. 2008, A&A, 489, 633 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Rüdiger, G., & Kitchatinov, L. L. 2005, A&A, 434, 629 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Salmeron, R., & Wardle, M. 2003, MNRAS, 345, 992 [NASA ADS] [CrossRef] (In the text)
 Salmeron, R., & Wardle, M. 2005, MNRAS, 361, 45 [NASA ADS] [CrossRef] (In the text)
 Salmeron, R., & Wardle, M. 2008, MNRAS, 388, 1223 [NASA ADS] (In the text)
 Sano, T., & Miyama, S. M. 1999, ApJ, 515, 776 [NASA ADS] [CrossRef] (In the text)
 Sano, T., & Stone, J. M. 2002, ApJ, 577, 534 [NASA ADS] [CrossRef] (In the text)
 Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef] (In the text)
 Semenov, D., Wiebe, D., & Henning, T. 2004, A&A, 417, 93 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
 Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685 [NASA ADS] [CrossRef] (In the text)
 Simon, J. B., Bai, X.N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ, 775, 73 [NASA ADS] [CrossRef] (In the text)
 Simon, J. B., Bai, X.N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ, 764, 66 [NASA ADS] [CrossRef] (In the text)
 Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656 [NASA ADS] [CrossRef] (In the text)
 Suzuki, T. K., & Inutsuka, S.i. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] (In the text)
 Tóth, G., Ma, Y., & Gombosi, T. I. 2008, J. Comput. Phys., 227, 6967 [NASA ADS] [CrossRef] (In the text)
 Turner, N. J., & Sano, T. 2008, ApJ, 679, L131 [NASA ADS] [CrossRef] (In the text)
 Turner, N. J., Sano, T., & Dziourkevitch, N. 2007, ApJ, 659, 729 [NASA ADS] [CrossRef] (In the text)
 Turner, N. J., Benisty, M., Dullemond, C. P., & Hirose, S. 2014a, ApJ, 780, 42 [NASA ADS] [CrossRef] (In the text)
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014b, in Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, & T. Henning (University of Arizona Press) (In the text)
 Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617 [NASA ADS] (In the text)
 Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69 [NASA ADS] [CrossRef] (In the text)
 Urpin, V., & Rüdiger, G. 2005, A&A, 437, 23 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Velikhov, E. P. 1959, Sov. Phys.JETP, 36, 1398 (In the text)
 Wardle, M. 1999, MNRAS, 307, 849 [NASA ADS] [CrossRef] (In the text)
 Wardle, M. 2007, Ap&SS, 311, 35 [NASA ADS] [CrossRef] (In the text)
 Wardle, M., & Ng, C. 1999, MNRAS, 303, 239 [NASA ADS] [CrossRef] (In the text)
 Wardle, M., & Salmeron, R. 2012, MNRAS, 422, 2737 [NASA ADS] [CrossRef] (In the text)
Online material
Appendix A: Numerical implementation of the Hall effect
In this appendix, we provide details of the modifications made to the Pluto code (v4.0) to include the Hall effect. While these modifications are quite general to any conservative finitevolume numerical scheme, the interested reader would benefit by familiarising themselves with both Mignone et al. (2007) and the Pluto user manual.
Appendix A.1: Conservative finitevolume scheme
The implementation of Hall effect in Pluto has proven to be a difficult task. Our naive attempts to include the Hall effect as a source term in a way that is analogous to that used for Ohmic diffusion have failed. Since the Hall effect is in essence a dispersive term rather than a true diffusive term, we have instead incorporated the Hall effect into the heart of the conservative integration scheme.
We begin by writing Eqs. (1)–(3) in conservative form: (A.1)where U is a vector of conserved quantities, F is the conservative flux function, and S is the source term function. The equations of isothermal HallMHD dictate (A.2)where is the total pressure and x_{H} = η_{H}/B is the control parameter for the Hall effect.
In the above formulation, we have explicitly used J = ∇×B in the flux function, which implies that the flux depends upon the conserved variables and their derivatives. This has dramatic consequences on the mathematical nature of the Riemann problem since it is illdefined in this formulation (to wit, if B is discontinuous, then J is not defined). Another way to see this is to note that the linearised flux function, which defines the characteristic speeds of the system, depends upon the wavelength of the perturbation (due to the presence of whistler waves). If a discontinuity appears in the flow, one obtains infinitely fast characteristic speeds, which are, of course, unphysical.^{6} Because of this mathematical difficulty, we have not found any sensible way to derive accurate Riemann solvers, such as the Roe or HLLD solvers. Instead, we follow Tóth et al. (2008) and implement an HLL solver, which approximates the dynamics of HallMHD by assuming that J is an external parameter (i.e., not related to B). In doing so, we circumvent the difficulties exposed above at the expense of large numerical diffusivities.
The algorithm we designed to integrate Eqs. (A.1) and (A.2) using Pluto is as follows:

1.
Compute the primitive cellaveraged variables V from the conserved variables U.

2.
Reconstruct the primitive variables at the cell edges using a “wellchosen” secondorder, totalvariationdiminishing (TVD) spatialreconstruction scheme (see below). This defines a left state (V_{L}) and a right state (V_{R}) at each cell face.

3.
Compute the left and right conserved variables U_{R /L} from V_{R /L}.

4.
Compute the facecentered J from the cellaveraged B using finitedifference formulae. Care is taken at this stage to apply the proper boundary conditions; shearingsheet boundary conditions have to be applied explicitly to the currents to avoid spurious oscillations at radial boundaries due to interpolation errors.

5.
Compute the left and right fluxes using the left and right states and the facecentered current: F_{L/R} = F(U_{L / R},J).

6.
Compute the Godunov flux F^{∗} using the whistlermodified HLL solver (described below).

7.
Evolve the conservative variables in time according to the equations of motion (see Mignone et al. 2007 for details). If constrained transport is used (Evans & Hawley 1988; Balsara & Spicer 1999), then the Godunov flux is used to compute the edgecentered electromotive forces (EMFs), which are then used to evolve the magnetic field.
To compute the Godunov flux F^{∗}, we follow the HLL scheme by using the approximation where S_{L} is the smallest algebraic signal speed for the left state and S_{R} is the largest algebraic signal speed for the right state. Since we are solving the HallMHD equations, S_{R /L} includes both the fast magnetosonic speed and the whistler wave speed. Therefore, we choose where v is the flow speed, c_{f} is the fast magnetosonic speed, and c_{w} is the whistler wave speed. As whistler waves are dispersive, we choose c_{w} to be equal to the whistler speed at the grid scale: where Δx is the grid spacing in the direction under consideration.
Appendix A.2: Numerical tests
Appendix A.2.1: Linear dispersion relation
We test our integration scheme by first considering a simple configuration with a uniform mean magnetic field B_{0} plus small sinusoidal perturbations δB_{⊥} perpendicular to B_{0}. We then compute the frequency of the oscillation obtained in the code and compare it to the theoretical dispersion relation (see KL13 for details). The results are presented in Fig. A.1, where we have quantified the Hall effect by the Hall lengthscale ℓ_{H} and the Hall frequency ω_{H} ≡ v_{A}/ℓ_{H}.
In these numerical calculations, the Nyquist frequency is such that kℓ_{H} = 80. We find very good agreement for the whistler branch up to half of the Nyquist frequency. The nonpropagating ioncyclotron wave is not obtained at large k because of the large numerical dissipation: the dissipation rate of this wave is faster than its oscillation period.
Fig. A.1
Dispersion relation for whistler waves. Black line: analytical prediction; circles: eigenfrequencies measured in Pluto using our implementation of the Hall effect. 

Open with DEXTER 
Numerical damping rate γ/ω_{H} of a whistler wave as a function of the number of points per wavelength n using the monotonized centered (MC), Van Leer (VL), and minmod (MM) slope limiters.
Appendix A.2.2: Convergence and numerical dissipation
The question of convergence in HallMHD was raised by Tóth et al. (2008). Since the whistler wave speed ∝(dx)^{1} and since numerical dissipation is roughly proportional to the wave speed, the convergence of the code at second order is not guaranteed in the presence of fast whistler waves. As was demonstrated by Tóth et al. (2008), this problem can be solved using a symmetric slope limiter in the reconstruction scheme, which naturally leads to a higherorder numerical diffusivity. To verify this, we present the numerical damping rate γ/ω_{H} for a whistler wave at kℓ_{H} = 20 as a function of the number of points per wavelength n in Table A.1. We find that both the monotonized centered and Van Leer slope limiters show secondorder convergence as the resolution per wavelength is increased. This is not the case for the minmod limiter, where only firstorder convergence is obtained. These results support the Tóth et al. (2008) argument and demonstrate that great care is needed when choosing the slope limiter in HallMHD. Unless otherwise stated, we always use the monotonized centered slope limiter.
Tóth et al. (2008) also demonstrated that this algorithm leads to numerical dissipation ~c_{w}U^{(4)}(Δx)^{3}, where U^{(n)} stands for the nth derivative with respect to x. This implies that

numerical dissipation damps whistler waves at the grid scale in one oscillation period; this explains why it is not possible to measure the numerical eigenfrequency in Sect. A.2.1 at the Nyquist frequency;

at a given resolution, the damping rate decreases as k^{4} (This has been verified by our numerical implementation.). Hence, numerical dissipation decreases very rapidly as one goes to larger scales.
Finally, KL13 have shown that higherorder timeintegration schemes are one way to guarantee the stable propagation of whistler waves if one desires spectral accuracy without any numerical dissipation. Because of the numerical dissipation inherent to our algorithm, which is a natural consequence of the Godunov flux as defined above, the use of a higherorder timeintegration scheme is not required in Pluto. We use an explicit secondorder accurate RungeKutta scheme unless otherwise stated.
Appendix A.2.3: Unstratified shearing box
Fig. A.2
Snapshot of B_{z} in our unstratified HallMRI run at t = 30. The flow exhibits a zonalfield structure similar to that discovered by KL13. 

Open with DEXTER 
KL13 demonstrated that the Halldominated MRI with ⟨ B_{y} ⟩ ≲ ⟨ B_{z} ⟩ saturates by producing largescale axisymmetric structures in the magnetic field (“zonal fields”). To further test our HallMHD integration scheme and to control our numerical diffusivity, we have tried to reproduce this nonlinear behaviour in Pluto. We use very similar parameters to the run ZB1H1 of KL13. We use a 4 × 4 × 1 box with a resolution 64 × 64 × 16. We set ℓ_{H} = 0.55, B_{z0} = 0.03, and η = η_{A} = 0. The resulting saturated state is shown in Fig. A.2 and exhibits a zonalfield structure similar to the one presented by KL13. This demonstrates that our scheme, despite being relatively diffusive due to the use of an HLL solver, manages to capture zonalfield structures with a resolution of 16 points per H.
Appendix B: Ambipolar diffusion
Ambipolar diffusion is implemented as a source term in a way analogous to that used for Ohmic resistivity. However, we noticed that gridscale instabilities sometimes arise if the shearingsheet boundary conditions are not enforced on the current density (a similar thing was observed with the Hall effect). Therefore, we use the current density computed in the conservative scheme for the Hall effect to compute the EMFs associated with ambipolar diffusion and Ohmic dissipation. This prevents smallscale instabilities from occurring at the radial boundaries.
We have tested our ambipolar diffusion module against the predictions of Kunz & Balbus (2004). An isothermal shearing box of size 4 × 4 × 1 and resolution 128 × 32 × 32 is threaded by a mean magnetic field with both vertical and horizontal components, whose magnitudes are characterised by their respective Alfvén speeds: B_{z0} = 0.025 and B_{y0} = 0.1. We introduce ambipolar diffusion by setting Am = 1. A typical snapshot of the linear growth phase is shown in Fig. B.1. We find the fastestgrowing mode to have (k_{x},k_{z}) = 2π × (−0.75,1). To check that this conforms to the expectation from linear theory, we show the theoretical growth rate of the ambipolardominated MRI as a function of (k_{x},k_{z}) in Fig. B.2. We find that the mode seen in the simulation corresponds to the fastestgrowing mode from linear theory that is available in the simulation. The theoretical growth rate for this mode is γ = 0.171; we find numerically γ = 0.17. This demonstrates the accuracy of our implementation of ambipolar diffusion in Pluto.
Fig. B.1
Snapshot of B_{x} in our ambipolarMRI test simulation at t = 100, just before saturation. The oblique nature of MRI channel mode, as predicted by Kunz & Balbus (2004), is evident. 

Open with DEXTER 
Fig. B.2
Growth rate of the ambipolarMRI for B_{z0} = 0.025 and B_{y0} = 0.1 as a function of the horizontal and vertical wavenumbers, k_{x} and k_{z}. The largest growth rate available in the simulation obtains for (k_{x},k_{z}) = 2π × (−0.75,1) (white cross). 

Open with DEXTER 
All Tables
Numerical damping rate γ/ω_{H} of a whistler wave as a function of the number of points per wavelength n using the monotonized centered (MC), Van Leer (VL), and minmod (MM) slope limiters.
All Figures
Fig. 1
Ionisation fraction x_{e} versus height z at 1, 5, and 10 au. 

Open with DEXTER  
In the text 
Fig. 2
Ohmic (Λ; Eq. (15)), Hall (Ha; Eq. (16)), and ambipolar (Am; Eq. (17)) Elsasser numbers in the initial state versus height z at 1 au (top) and 10 au (bottom), assuming a constant vertical magnetic field with β_{0} = 10^{5}. 

Open with DEXTER  
In the text 
Fig. 3
Hall lengthscale ℓ_{H} (Eq. (18)) versus height z at 1, 5, and 10 au. The dotted line denotes the threshold above which saturation via zonal magnetic fields occurs in unstratified simulations with net vertical magnetic flux (KL13). 

Open with DEXTER  
In the text 
Fig. 4
Comparison of B_{y}(z) from the simulation (red circles) and from the linear calculation (black line), which include Ohmic and Hall diffusion computed at 1 au with β_{0} = 10^{5}. Top: n = 1 eigenmode and symmetric component of B_{y}(z). Bottom: n = 2 eigenmode and antisymmetric component of B_{y}(z). 

Open with DEXTER  
In the text 
Fig. 5
Space–time diagram of the horizontally averaged azimuthal magnetic field, ⟨ B_{y} ⟩, in the Ohmic (1O5; top) and OhmicHall (1OH5; bottom) runs. 

Open with DEXTER  
In the text 
Fig. 6
Space–time diagram of the logarithm of the horizontally averaged Maxwell stress, log ⟨ −B_{x}B_{y} ⟩, in the Ohmic (1O5; top) and OhmicHall (IOH5; bottom) runs. 

Open with DEXTER  
In the text 
Fig. 7
Horizontally averaged density profile (top) and vertical velocity (bottom) versus height z, averaged over ≈140 orbits, from runs 1O5 (solid line) and 1OH5 (dashed line). 

Open with DEXTER  
In the text 
Fig. 8
Contributions in run 1OH5 to the righthand side of the horizontally averaged induction equation deduced from Eqs. (20a; top) and (20b; bottom). 

Open with DEXTER  
In the text 
Fig. 9
Space–time evolution of the logarithm of the horizontallyaveraged magnetic stress, log ⟨ M_{xy} ⟩, in the Ohmic (1O5; top), Ohmicambipolar (1OA5; middle), and OhmicambipolarHall (1OHA5; bottom) runs. 

Open with DEXTER  
In the text 
Fig. 10
Horizontally averaged Maxwell (top) and Reynolds (bottom) stress versus height z, averaged over 700 Ω^{1}, from runs 1O5 (dashdotted), 1OA5 (dashed), and 1OAH5 (solid). 

Open with DEXTER  
In the text 
Fig. 11
Horizontally averaged turbulent plasma β versus height z, averaged over 700 Ω^{1}, from runs 1O5 (dashdotted), 1OA5 (dashed), and 1OAH5 (solid). 

Open with DEXTER  
In the text 
Fig. 12
Horizontally and temporally averaged Maxwell (top) and Reynolds (bottom) stress versus height z for runs 1O3 (dashdotted), 1OA3 (dashed), and 1OAH3 (solid). 

Open with DEXTER  
In the text 
Fig. 13
Horizontally averaged Maxwell stress versus height z, averaged over 700 Ω^{1}, from runs 1OAH5 (dashdotted), 5OAH5 (dashed), and 10OAH5 (plain line). 

Open with DEXTER  
In the text 
Fig. 14
Space–time evolution of the horizontally averaged azimuthal magnetic field, ⟨ B_{y} ⟩, as a function of time in run 1OA5e. The run starts with odd symmetry for the azimuthal field. At t ≃ 1000 Ω^{1} a current layer forms and is ejected at t = 1410 Ω^{1}, leaving an azimuthal field with even symmetry. 

Open with DEXTER  
In the text 
Fig. A.1
Dispersion relation for whistler waves. Black line: analytical prediction; circles: eigenfrequencies measured in Pluto using our implementation of the Hall effect. 

Open with DEXTER  
In the text 
Fig. A.2
Snapshot of B_{z} in our unstratified HallMRI run at t = 30. The flow exhibits a zonalfield structure similar to that discovered by KL13. 

Open with DEXTER  
In the text 
Fig. B.1
Snapshot of B_{x} in our ambipolarMRI test simulation at t = 100, just before saturation. The oblique nature of MRI channel mode, as predicted by Kunz & Balbus (2004), is evident. 

Open with DEXTER  
In the text 
Fig. B.2
Growth rate of the ambipolarMRI for B_{z0} = 0.025 and B_{y0} = 0.1 as a function of the horizontal and vertical wavenumbers, k_{x} and k_{z}. The largest growth rate available in the simulation obtains for (k_{x},k_{z}) = 2π × (−0.75,1) (white cross). 

Open with DEXTER  
In the text 