EDP Sciences
Free Access
Issue
A&A
Volume 566, June 2014
Article Number A56
Number of page(s) 17
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201423660
Published online 11 June 2014

© ESO, 2014

1. Introduction

The formation and evolution of protoplanetary discs is a long-standing problem in accretion theory. Observations indicate that T Tauri stars accrete mass from their surrounding discs at a rate of 10-910-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 (so-called α-disc models; Shakura & Sunyaev 1973), or transported vertically by large-scale 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 angular-momentum 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 contender1 to explain accretion in many types of astrophysical systems. However, protoplanetary discs are cold (T ≲ 300 K at 1 au), dense (n ≳ 1013 cm-3), and thus poorly ionised. Their evolution is well-described by a set of magnetohydrodynamic (MHD) equations that include three non-ideal 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 non-ideal 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 layered-accretion model has been studied extensively with increasingly complex analytical models (Sano et al. 2000; Fromang et al. 2002; Ilgner & Nelson 2006a,b) and three-dimensional 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 non-linear 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 Hall-dominated magnetorotational turbulence, whereby the vertical magnetic field self-organises into coherent, long-lived axisymmetric (“zonal”) structures and the rate of angular-momentum transport becomes vanishingly small. While interesting in itself, this result was based on a crude representation of a protoplanetary disc – an incompressible, dust-free and unstratified shearing box without ion-neutral 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 three-dimensional 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 non-ideal 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 finite-volume code Pluto. A simplified ionisation model is used to determine vertical diffusivity profiles for all three non-ideal 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 non-ideal MHD in the stratified shearing-box approximation. In Sects. 2.2 and 2.3, we introduce our numerical method and define the diagnostics used to quantify the efficiency of angular-momentum 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.13.2) and then with all three non-ideal MHD effects treated self-consistently (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 shearing-box model

We study the evolution of a poorly ionised protoplanetary disc in the shearing-box approximation (Hawley et al. 1995). A local patch of the disc centred at R0 with an extent ~H (the vertical scale height) rotates with the disc at a constant angular speed Ω = ΩK(R0). The local frame is defined by (ex,ey,ez), where ex is aligned with the radial direction, ey with the azimuthal direction, and ez with the vertical direction. The computational domain has a size (Lx,Ly,Lz) with ; we assume the disc to be geometrically thin, which implies HR0. This allows us to neglect curvature terms – the so-called Hill approximation (Hill 1878). Taking the flow to be locally isothermal, the equations of non-ideal MHD that govern the evolution of the mass density ρ, the velocity u, and the magnetic field B are, respectively, where cs is the isothermal sound speed, Ω = Ω ez is the local angular velocity, g = 2qΩ2xex − Ω2zez is the local gravitational acceleration with q = −dlnΩ/dlnR (=3/2 for a Keplerian disc), eb = B/ | B |, and J = ×B. We include all three non-ideal 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 non-linear.

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 Hcs and Bz0 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Ωxey. 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 polarity-dependent way (Wardle 1999; Balbus & Terquem 2001).

The boundary conditions are shearing-periodic in the x direction and periodic in the y direction. The vertical boundary conditions are:

  • vertical hydrostatic equilibrium for ρ;

  • zero vertical gradient for vx and vy;

  • outflow for vz; and

  • vertical field for B.

We have checked that the final condition can be replaced by a zero-current condition on the boundary without having any significant impact on our results.

2.2. Numerical method

We use a modified version of the finite-volume 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 second-order 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. Face-centered 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 second-order-accurate Runge-Kutta scheme2.

We choose our units such that Ω = cs = ρ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 non-ideal MHD effects at R0 = 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 1-OHA-5 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 vx 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 angular-momentum transport and the production of outflows. These rely on two averaging procedures: The amount of transport is determined by the Reynolds stress Rij = ρvivj and the Maxwell stress Mij = −BiBj. This allows us to define an effective (radial) transport parameter, (8)Outflows are characterised by the mass-loss rate, (9)and the (dimensionless) surface magnetic stress (10)which is evaluated at the top (zt = 4.5H) and bottom (zb = −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).

Table 1

List of runs discussed in this paper.

These diagnostics can be directly related to large-scale quantities, such as the accretion rate acc. Assuming a 1 M central object, smooth surface density and temperature profiles, and neglecting pre-factors 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 R0 = 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 ζ:

The ionisation fraction xe is obtained by balancing these ionisation sources with dissociative recombination in a metal- and dust-free environment (Gammie 1996; Fromang et al. 2002): (12)where αdr = 3 × 10-6T−1/2 cm2 s-1 is the dissociative recombination rate coefficient for molecular ions.

thumbnail Fig. 1

Ionisation fraction xe 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 xFUV ~ 10-510-4 for a penetration depth ~10-2 g cm-2 (Perez-Becker & 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 electron-neutral collision rate (Draine et al. 1983), ρi is the ion mass density, γi = ⟨ σvi/ (mn + mi), and is the ion-neutral collision rate (Draine 2011); mn and mi 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 X-rays to temperatures up to ~8000 K (Aresu et al. 2011). Our simple isothermal model thus breaks down for column densities 10-310-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.

thumbnail 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 = 105.

Open with DEXTER

thumbnail 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 magnetic-field 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 magnetic-field strength, is the Hall lengthscale: (18)When H ≳ 0.2H, KL13 found that unstratified simulations of the Hall-MRI 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 zonal-field regime.

The non-ideal 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 ΩH2, 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 ΩBz0 > 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 MRI-unstable, 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 1-OH-5, 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, By(z), at t = 10 Ω-1 with that predicted by linear theory, where the latter is calculated with a pseudo-spectral 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 fastest-growing 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 By(z) from the numerical simulation into symmetric (σ = 1) and antisymmetric (σ = −1) parts, each of which can then be compared to the two fastest-growing 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 buoyancy-driven 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.

Table 2

Fastest-growing eigenmodes in run 1-OH-5 with growth rate γ (in units of Ω) and symmetry label σ (see Eq. (19)).

thumbnail Fig. 4

Comparison of By(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 = 105. Top: n = 1 eigenmode and symmetric component of By(z). Bottom: n = 2 eigenmode and antisymmetric component of By(z).

Open with DEXTER

3.2. Non-linear evolution of stratified discs subject to Ohmic and Hall diffusion

3.2.1. Fiducial runs

In Fig. 5, we present space–time diagrams of By for runs 1-O-5 (top) and 1-OH-5 (bottom). The Ohmic run (1-O-5) 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 (1-OH-5), this butterfly pattern is replaced by a strong quasi-steady 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. 6-top). 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 Hall-dominated region (−3Hz ≲ 3H; Fig. 6-bottom). This stress is not a “turbulent” stress in the usual sense but rather is a consequence of the large-scale 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 1-O-5.

thumbnail Fig. 5

Space–time diagram of the horizontally averaged azimuthal magnetic field, By, in the Ohmic (1-O-5; top) and Ohmic-Hall (1-OH-5; bottom) runs.

Open with DEXTER

The presence of a vertical magnetic flux is known to trigger outflows in shearing-box 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: mass-loss rates increase by more than two orders of magnitude. However, the outflows do not have any well-defined orientation: they can be directed toward either positive or negative x, which is a phenomenon previously observed in ideal-MHD simulations presented by Fromang et al. (2013) and Bai & Stone (2013a). This is one consequence of the neglect of curvature terms in the shearing-box approximation. Because of this asymmetry, we find no average stress exerted at the disc surface by the outflow.

thumbnail Fig. 6

Space–time diagram of the logarithm of the horizontally averaged Maxwell stress, log ⟨ −BxBy, in the Ohmic (1-O-5; top) and Ohmic-Hall (I-OH-5; bottom) runs.

Open with DEXTER

The very large mass-loss rate found in the Ohmic-Hall 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 1-O-5 and 1-OH-5.

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 Hall-amplified azimuthal magnetic field dominates the total pressure in the disc midplane.

thumbnail Fig. 7

Horizontally averaged density profile (top) and vertical velocity (bottom) versus height z, averaged over ≈140 orbits, from runs 1-O-5 (solid line) and 1-OH-5 (dashed line).

Open with DEXTER

3.2.2. Origin of the Hall-driven azimuthal field

thumbnail Fig. 8

Contributions in run 1-OH-5 to the right-hand 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 Hall-MHD 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 horizontalmagnetic-field 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 1-OH-5 is shown in Fig. 8 with each of these terms identified4. We focus on the inner disc region | z | < 3H, where Hall and Ohmic diffusion dominate the dynamics.

thumbnail Fig. 9

Space–time evolution of the logarithm of the horizontally-averaged magnetic stress, log ⟨ Mxy, in the Ohmic (1-O-5; top), Ohmic-ambipolar (1-OA-5; middle), and Ohmic-ambipolar-Hall (1-OHA-5; bottom) runs.

Open with DEXTER

We first describe the contributions of the right-hand side of the x-component of the induction equation (Fig. 8-top). 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 Bx < 0 in the disc midplane), and the outflow term, which carries Bx out of the midplane. This quasi-equilibrium is different than that found in ideal-MHD simulations, in which the driving term is due to stretching. On the other hand, the y-component of the induction equation (Fig. 8- bottom) is a resistive MHD quasi-equilibrium, 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 self-consistent picture of run 1-OH-5. We first note that Moreover, since the fluctuations δBz are small compared to the mean vertical field Bz0, we have zByBz ⟩ ≃ Bz0zBy and zBxBz ⟩ ≃ Bz0zBx. Retaining only the driving terms in Eq. (18) and neglecting the outflow, we find These reduced equations for Bx and By describe a non-local version of the “Hall-shear” instability (cf. Kunz 2008, Eq. (46)): Keplerian shear generates an azimuthal magnetic-field component from a radial one, and the Hall effect conservatively reorients this azimuthal field back into the radial direction. When qΩBz0 > 0, this feedback loop can lead to growth. Note that the Hall-shear instability is not a version of the MRI and does not rely on the Coriolis force (though it can be obtained from the Hall-MRI 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 vd is the drift velocity of the electrons relative to the ions. Hence, the term on the right-hand 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 non-linear saturated state results from a balance between this Hall-shear 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. Non-linear 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 (1-O-5), Ohmic-ambipolar (1-OA-5), and Ohmic-ambipolar-Hall (1-OAH-5) runs at R0 = 1 au is presented in Fig. 9. The first two simulations are quantitatively similar to runs O-b5 and OA-b5 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.

thumbnail Fig. 10

Horizontally averaged Maxwell (top) and Reynolds (bottom) stress versus height z, averaged over 700 Ω-1, from runs 1-O-5 (dash-dotted), 1-OA-5 (dashed), and 1-OAH-5 (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 Mxy = 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 1-OAH-5 (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 quasi-equilibrium state, producing mass-loss rates a factor of 20 larger than those obtained in the Ohmic-ambipolar simulations.

thumbnail Fig. 11

Horizontally averaged turbulent plasma β versus height z, averaged over 700 Ω-1, from runs 1-O-5 (dash-dotted), 1-OA-5 (dashed), and 1-OAH-5 (solid).

Open with DEXTER

3.3.2. Parameter study

Field polarity:

it is well known that the linear Hall-MRI is sensitive to the field polarity (Wardle 1999; Balbus & Terquem 2001). To test whether or not this result carries over to the non-linear regime, we have performed a simulation identical to run 1-OHA-5 but with a reversed mean field: (run 1-OHA-mB). In this case, we find that the level of transport, the surface stress, and the mass-loss rate to be smaller than in the run without the Hall effect (run 1-OA-5). 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 non-linearly 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.

thumbnail Fig. 12

Horizontally and temporally averaged Maxwell (top) and Reynolds (bottom) stress versus height z for runs 1-O-3 (dash-dotted), 1-OA-3 (dashed), and 1-OAH-3 (solid).

Open with DEXTER

Zero vertical-net-flux configuration:

the zero vertical-net-flux 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 non-ideal effects (run 1-OHA-znf). This simulation is identical to run 1-OHA-5 except that the initial magnetic-field configuration is given by Bz = Bz0sin(2πx/Lx), where Bz0 is the initial field strength from run 1-OHA-5. This initial condition rapidly collapses into a quiet state with virtually no transport and no outflow (see Table 1).

Mean field strength:

the runs 1-xxx-5 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 = 103 (runs 1-O-3, 1-OA-3, 1-OAH-3). 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 ideal-MHD 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 1-OHA-5; 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 = 103. 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, R0 (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 = 105, varying R0 from 1 au to 10 au. The resulting Maxwell stress profiles are presented in Fig. 135.

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 mass-loss 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 1-OH-5. 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.

thumbnail Fig. 13

Horizontally averaged Maxwell stress versus height z, averaged over 700 Ω-1, from runs 1-OAH-5 (dash-dotted), 5-OAH-5 (dashed), and 10-OAH-5 (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 ambipolar-diffusion–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 off-midplane 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 long-lived σ = −1 solutions when only Ohmic and ambipolar diffusion are included (see run 1-OA-5-e in Fig. 14). These solutions start to develop a strong off-midplane 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.

thumbnail Fig. 14

Space–time evolution of the horizontally averaged azimuthal magnetic field, By, as a function of time in run 1-OA-5-e. 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 self-organisation in Hall-dominated magnetorotational turbulence

KL13 have shown that the saturated state favoured by Hall-dominated 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 (By ⟩ = 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 By ⟩ ~ 200 ⟨ Bz. As pointed out by KL13, this magnetic configuration (By ⟩ ≫ ⟨ Bz) 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, micron-sized 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μm-sized 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 H2 are drastically reduced. The diffusivities associated to the three non-ideal 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 non-ideal 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 angular-momentum transport observed in our Hall-dominated 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 gz is the vertical component of gravity, vz is the vertical gas velocity, and σvg is the rate of momentum exchange of grains with neutral gas molecules. We approximate σvg ~ a2cs and mg ~ a3ρS, where ρS is the grain material density and a is the grain radius. Using gz = Ω2z, we find to be the equilibrium height reached by the grains when they are dragged by the outflow. The continuity equation implies that ρvz is conserved along z, defining the mass-outflow rate. Assuming ρS = 1 g cm-3 and using Eq. (9) to quantify the mass-loss 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 mass-loss 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 R0 = 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 10-OA-gr and 10-OHA-gr). Even in this worst-case 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 non-linear behaviour of poorly ionised, vertically stratified protoplanetary discs. For the first time, all three relevant non-ideal MHD effects – Ohmic dissipation, ambipolar diffusion, and the Hall effect – are self-consistently included. To accomplish this feat, we have implemented an original formulation of ambipolar diffusion and the Hall effect in the finite-volume 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 angular-momentum 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 MRI-driven outflows are also modified by the Hall effect. The mass-loss 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 mass-loss rates) are affected by the vertical boundary conditions in the shearing-box approximation (Fromang et al. 2013; Lesur et al. 2013). Shearing-box outflows should be treated with caution.

Surprisingly, our stratified simulations do not produce the axisymmetric (“zonal”) magnetic-field 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 build-up of azimuthal flux is not possible in unstratified shearing boxes. It is suggestive that KL13 did not find the zonal-field route to saturation in unstratified shearing boxes with By ⟩ ≫ ⟨ Bz. 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 well-known 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 non-linear 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 long-term evolution of the large-scale magnetic field.

The relatively large magnetic stresses generated by the Hall effect are not associated with turbulent fluctuations. Rather, they are associated with large-scale amplification of an azimuthal magnetic field via the shearing of a Hall-induced radial field (cf. Kunz 2008). Indeed, the Reynolds stresses are always very small in the midplane (with Rxy ⟩ ≲ 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 large-scale stress. One consequence of our result is that the Hall-dominated 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 non-ideal 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 time-dependent, 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 angular-momentum transport in protoplanetary discs.


1

Hydrodynamical processes are an interesting alternative, despite serious theoretical difficulties; see Turner et al. (2014b) for a review.

2

Super-time-stepping schemes are of no use in our case since the Hall effect is directly incorporated into the Riemann solver and causes one of the most severe CFL conditions on the time step.

3

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.

4

The Hall term is artificially large in Fig. 8-top. 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 post-processing horizontally averaged quantities.

5

The Reynolds stress is negligible in these runs.

6

In a realistic plasma, the whistler wave speed is limited by finite Larmor radius effects that are excluded in the MHD approximation.

Acknowledgments

G.L. is indebted to Wing-Fai 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 PCIG09-GA-2011-294110. Support for M.W.K. was provided by NASA through Einstein Postdoctoral Fellowship Award Number PF1-120084, issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. Support for S.F. was provided by the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/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.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) 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

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 finite-volume 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 finite-volume 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 Hall-MHD dictate (A.2)where is the total pressure and xH = η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 ill-defined 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 Hall-MHD 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 cell-averaged variables V from the conserved variables U.

  • 2.

    Reconstruct the primitive variables at the cell edges using a “well-chosen” second-order, total-variation-diminishing (TVD) spatial-reconstruction scheme (see below). This defines a left state (VL) and a right state (VR) at each cell face.

  • 3.

    Compute the left and right conserved variables UR /L from VR /L.

  • 4.

    Compute the face-centered J from the cell-averaged B using finite-difference formulae. Care is taken at this stage to apply the proper boundary conditions; shearing-sheet 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 face-centered current: FL/R = F(UL / R,J).

  • 6.

    Compute the Godunov flux F using the whistler-modified 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 edge-centered 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 SL is the smallest algebraic signal speed for the left state and SR is the largest algebraic signal speed for the right state. Since we are solving the Hall-MHD equations, SR /L includes both the fast magnetosonic speed and the whistler wave speed. Therefore, we choose where v is the flow speed, cf is the fast magnetosonic speed, and cw is the whistler wave speed. As whistler waves are dispersive, we choose cw 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 B0 plus small sinusoidal perturbations δB perpendicular to B0. 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 ωHvA/H.

In these numerical calculations, the Nyquist frequency is such that kH = 80. We find very good agreement for the whistler branch up to half of the Nyquist frequency. The non-propagating ion-cyclotron 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.

thumbnail 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

Table A.1

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 Hall-MHD 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 higher-order numerical diffusivity. To verify this, we present the numerical damping rate γ/ωH for a whistler wave at kH = 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 second-order convergence as the resolution per wavelength is increased. This is not the case for the minmod limiter, where only first-order 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 Hall-MHD. 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 ~cwU(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 higher-order time-integration 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 higher-order time-integration scheme is not required in Pluto. We use an explicit second-order accurate Runge-Kutta scheme unless otherwise stated.

Appendix A.2.3: Unstratified shearing box

thumbnail Fig. A.2

Snapshot of Bz in our unstratified Hall-MRI run at t = 30. The flow exhibits a zonal-field structure similar to that discovered by KL13.

Open with DEXTER

KL13 demonstrated that the Hall-dominated MRI with By ⟩ ≲ ⟨ Bz saturates by producing large-scale axisymmetric structures in the magnetic field (“zonal fields”). To further test our Hall-MHD integration scheme and to control our numerical diffusivity, we have tried to reproduce this non-linear 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, Bz0 = 0.03, and η = ηA = 0. The resulting saturated state is shown in Fig. A.2 and exhibits a zonal-field 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 zonal-field 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 grid-scale instabilities sometimes arise if the shearing-sheet 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 small-scale 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: Bz0 = 0.025 and By0 = 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 fastest-growing mode to have (kx,kz) = 2π × (−0.75,1). To check that this conforms to the expectation from linear theory, we show the theoretical growth rate of the ambipolar-dominated MRI as a function of (kx,kz) in Fig. B.2. We find that the mode seen in the simulation corresponds to the fastest-growing 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.

thumbnail Fig. B.1

Snapshot of Bx in our ambipolar-MRI 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

thumbnail Fig. B.2

Growth rate of the ambipolar-MRI for Bz0 = 0.025 and By0 = 0.1 as a function of the horizontal and vertical wavenumbers, kx and kz. The largest growth rate available in the simulation obtains for (kx,kz) = 2π × (−0.75,1) (white cross).

Open with DEXTER

All Tables

Table 1

List of runs discussed in this paper.

Table 2

Fastest-growing eigenmodes in run 1-OH-5 with growth rate γ (in units of Ω) and symmetry label σ (see Eq. (19)).

Table A.1

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

thumbnail Fig. 1

Ionisation fraction xe versus height z at 1, 5, and 10 au.

Open with DEXTER
In the text
thumbnail 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 = 105.

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 4

Comparison of By(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 = 105. Top: n = 1 eigenmode and symmetric component of By(z). Bottom: n = 2 eigenmode and antisymmetric component of By(z).

Open with DEXTER
In the text
thumbnail Fig. 5

Space–time diagram of the horizontally averaged azimuthal magnetic field, By, in the Ohmic (1-O-5; top) and Ohmic-Hall (1-OH-5; bottom) runs.

Open with DEXTER
In the text
thumbnail Fig. 6

Space–time diagram of the logarithm of the horizontally averaged Maxwell stress, log ⟨ −BxBy, in the Ohmic (1-O-5; top) and Ohmic-Hall (I-OH-5; bottom) runs.

Open with DEXTER
In the text
thumbnail Fig. 7

Horizontally averaged density profile (top) and vertical velocity (bottom) versus height z, averaged over ≈140 orbits, from runs 1-O-5 (solid line) and 1-OH-5 (dashed line).

Open with DEXTER
In the text
thumbnail Fig. 8

Contributions in run 1-OH-5 to the right-hand side of the horizontally averaged induction equation deduced from Eqs. (20a; top) and (20b; bottom).

Open with DEXTER
In the text
thumbnail Fig. 9

Space–time evolution of the logarithm of the horizontally-averaged magnetic stress, log ⟨ Mxy, in the Ohmic (1-O-5; top), Ohmic-ambipolar (1-OA-5; middle), and Ohmic-ambipolar-Hall (1-OHA-5; bottom) runs.

Open with DEXTER
In the text
thumbnail Fig. 10

Horizontally averaged Maxwell (top) and Reynolds (bottom) stress versus height z, averaged over 700 Ω-1, from runs 1-O-5 (dash-dotted), 1-OA-5 (dashed), and 1-OAH-5 (solid).

Open with DEXTER
In the text
thumbnail Fig. 11

Horizontally averaged turbulent plasma β versus height z, averaged over 700 Ω-1, from runs 1-O-5 (dash-dotted), 1-OA-5 (dashed), and 1-OAH-5 (solid).

Open with DEXTER
In the text
thumbnail Fig. 12

Horizontally and temporally averaged Maxwell (top) and Reynolds (bottom) stress versus height z for runs 1-O-3 (dash-dotted), 1-OA-3 (dashed), and 1-OAH-3 (solid).

Open with DEXTER
In the text
thumbnail Fig. 13

Horizontally averaged Maxwell stress versus height z, averaged over 700 Ω-1, from runs 1-OAH-5 (dash-dotted), 5-OAH-5 (dashed), and 10-OAH-5 (plain line).

Open with DEXTER
In the text
thumbnail Fig. 14

Space–time evolution of the horizontally averaged azimuthal magnetic field, By, as a function of time in run 1-OA-5-e. 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
thumbnail 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
thumbnail Fig. A.2

Snapshot of Bz in our unstratified Hall-MRI run at t = 30. The flow exhibits a zonal-field structure similar to that discovered by KL13.

Open with DEXTER
In the text
thumbnail Fig. B.1

Snapshot of Bx in our ambipolar-MRI 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
thumbnail Fig. B.2

Growth rate of the ambipolar-MRI for Bz0 = 0.025 and By0 = 0.1 as a function of the horizontal and vertical wavenumbers, kx and kz. The largest growth rate available in the simulation obtains for (kx,kz) = 2π × (−0.75,1) (white cross).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.