Dust Segregation in Hall-dominated Turbulent Protoplanetary Disks

We show below the continuity and momentum equations for the gas, and one of the dust species (we consider three dust species in total). In addition we present the induction equation in a full non-ideal regime:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{\rho }_{{\rm{g}}}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{g}}}{{\boldsymbol{v}}}_{{\rm{g}}}) & = & 0,\\ {\partial }_{t}{\rho }_{{\rm{d}}}+{\rm{\nabla }}\cdot ({\rho }_{{\rm{d}}}{{\boldsymbol{v}}}_{{\rm{d}}}) & = & 0,\\ {\rho }_{{\rm{g}}}({\partial }_{t}{{\boldsymbol{v}}}_{{\rm{g}}}+{{\boldsymbol{v}}}_{{\rm{g}}}\cdot {\rm{\nabla }}{{\boldsymbol{v}}}_{{\rm{g}}}) & = & -{\rm{\nabla }}P-{\rho }_{{\rm{g}}}{\rm{\nabla }}{\rm{\Phi }}+{\boldsymbol{J}}\times {\boldsymbol{B}}-{{\boldsymbol{F}}}_{\mathrm{drag}},\\ {\rho }_{{\rm{d}}}({\partial }_{t}{{\boldsymbol{v}}}_{{\rm{d}}}+{{\boldsymbol{v}}}_{{\rm{d}}}\cdot {\rm{\nabla }}{{\boldsymbol{v}}}_{{\rm{d}}}) & = & -{\rho }_{{\rm{d}}}{\rm{\nabla }}{\rm{\Phi }}+{{\boldsymbol{F}}}_{\mathrm{drag}},\\ {\partial }_{t}{\boldsymbol{B}} & = & {\rm{\nabla }}\times ({{\boldsymbol{v}}}_{{\rm{g}}}\times {\boldsymbol{B}})-{\rm{\nabla }}\times { \mathcal E },\end{array}\end{eqnarray} \tag{ 1 }$$

where ${\rho }_{{\rm{g}}}$ and ${{\boldsymbol{v}}}_{{\rm{g}}}$ are the density and velocity of the gas, whereas ${\rho }_{{\rm{d}}}$ and ${{\boldsymbol{v}}}_{{\rm{d}}}$ correspond to the dust density and velocity, respectively. The term ${{\boldsymbol{F}}}_{\mathrm{drag}}$ denotes the drag force between gas and dust, specified in Equation (3). The non-ideal part of the electromotive force is given by

$$\begin{eqnarray}&&{ \mathcal E }={\eta }_{{\rm{O}}}{\boldsymbol{J}}\,+\,{\eta }_{{\rm{H}}}{\boldsymbol{J}}\times {\widehat{{\boldsymbol{e}}}}_{{\rm{B}}}\,-\,{\eta }_{{\rm{A}}}{\boldsymbol{J}}\times {\widehat{{\boldsymbol{e}}}}_{{\rm{B}}}\times {\widehat{{\boldsymbol{e}}}}_{{\rm{B}}},\end{eqnarray} \tag{ 2 }$$

with ${\widehat{{\boldsymbol{e}}}}_{{\rm{B}}}$ the unit vector along ${\boldsymbol{B}}$, and where the current density is ${\boldsymbol{J}}\equiv {\mu }_{0}^{-1}\,{\rm{\nabla }}\times {\boldsymbol{B}}$.

We furthermore adopt a locally isothermal equation of state where the disk temperature determines the sound speed ${c}_{{\rm{s}}}$, given by ${c}_{{\rm{s}}}^{2}(r)=P/{\rho }_{{\rm{g}}}$, with P being the gas pressure. The scalar variable ${\rm{\Phi }}\equiv -{{GM}}_{\odot }/r$ is the cylindrical gravitational potential of a central solar mass star.

For the drag force acting on the dust, we consider the Epstein regime as introduced by Whipple (1972) and Weidenschilling (1977), that is,

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{drag}}={{\rm{\Omega }}}_{{\rm{K}}}\,{\mathrm{St}}^{-1}\,({{\boldsymbol{v}}}_{{\rm{g}}}-{{\boldsymbol{v}}}_{{\rm{d}}}),\end{eqnarray} \tag{ 3 }$$

where ${{\rm{\Omega }}}_{{\rm{K}}}=\sqrt{{M}_{\odot }G/{r}^{3}}$ is the Keplerian angular frequency and $\mathrm{St}\equiv {{\rm{\Omega }}}_{{\rm{K}}}\,{t}_{\mathrm{stop}}$ is the so-called Stokes number. In the context of the present work, St is considered to be constant, which implies a radial dependence of the stopping time, ${t}_{\mathrm{stop}}$. The Epstein regime is valid as long as the mean free path of gas molecules is roughly larger than the particle radius. This is generally true for grain sizes smaller than a few centimeters, when adopting the disk model of Section 2.2. When including the drag force, the integration time of our simulations is short compared with the drift timescale. Hence, we can safely neglect the back-reaction of the dust onto the gas.

Neglecting charges carried by grains, the Ohmic (${\eta }_{{\rm{O}}}$), ambipolar (${\eta }_{{\rm{A}}}$), and Hall (${\eta }_{{\rm{H}}}$) diffusion coefficients can be expressed in terms of the species masses, m, number densities, n, and collision frequencies with neutrals, γ, as (Bai 2011)

$$\begin{eqnarray}\begin{array}{rcl}{\eta }_{{\rm{O}}} & = & \displaystyle \frac{{c}^{2}{\gamma }_{e}{m}_{e}{\rho }_{{\rm{g}}}}{4\pi {e}^{2}{n}_{e}}\sim \displaystyle \frac{1}{{x}_{e}},\\ {\eta }_{{\rm{H}}} & = & \displaystyle \frac{{cB}}{4\pi {{en}}_{e}}\sim \displaystyle \frac{1}{{x}_{e}}\displaystyle \frac{B}{{\rho }_{{\rm{g}}}},\\ {\eta }_{{\rm{A}}} & = & \displaystyle \frac{{B}^{2}}{4\pi {\gamma }_{i}{\rho }_{{\rm{g}}}{\rho }_{i}}\sim \displaystyle \frac{1}{{x}_{e}}{\left(\displaystyle \frac{B}{{\rho }_{{\rm{g}}}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 4 }$$

where ${x}_{e}={n}_{e}/n$ is the ionization fraction. It is common to quantify these terms using the Elsässer numbers

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{\rm{O}}}\equiv \displaystyle \frac{{v}_{A}^{2}}{{\eta }_{{\rm{O}}}{\rm{\Omega }}},\qquad {{\rm{\Lambda }}}_{{\rm{A}}}\equiv \displaystyle \frac{{v}_{A}^{2}}{{\eta }_{{\rm{A}}}{\rm{\Omega }}},\qquad {{\rm{\Lambda }}}_{{\rm{H}}}\equiv \displaystyle \frac{{v}_{A}^{2}}{{\eta }_{{\rm{H}}}{\rm{\Omega }}},\end{eqnarray} \tag{ 5 }$$

which characterize the coupling between the magnetic field and the neutral fluid at a length scale proportional to the characteristic wavelength of the MRI unstable modes (e.g., Wardle & Salmeron 2012).

Moreover, the relative importance of the Hall term can, in general, be defined via the so-called Hall length, ${l}_{{\rm{H}}}\equiv {\eta }_{{\rm{H}}}{v}_{{\rm{A}}}^{-1}$, (Pandey & Wardle 2008; Kunz & Lesur 2013). Adopting an incompressible non-stratified shearing box disk model, Kunz & Lesur (2013) found that when ${l}_{{\rm{H}}}/H\gtrsim 0.2$ the MRI saturates to a new regime where the self-organization mechanism start to be operative. The typical length H can be adopted as the hydrostatic scale height of a vertical stratified disk. Thus, ${l}_{{\rm{H}}}/H$ provides a dimensionless control parameter for this new saturation regime of the MRI.

In Section 5, we use the following disk model. The density and the sound speed are

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}=\displaystyle \frac{{\rm{\Sigma }}(r)}{\sqrt{2\pi }H(r)},\quad {c}_{{\rm{s}}}={c}_{{\rm{s}}0}\,{r}_{\mathrm{au}}^{-p},\end{eqnarray} \tag{ 6 }$$

with a surface density, ${\rm{\Sigma }}(r)$, and the hydrostatic scale height, H(r), defined as

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}{r}_{\mathrm{au}}^{-\sigma }\quad \mathrm{and}\quad H(r)={c}_{{\rm{s}}}/{{\rm{\Omega }}}_{{\rm{K}}},\end{eqnarray} \tag{ 7 }$$

respectively. We initialize the vertical and azimuthal magnetic field components as

$$\begin{eqnarray}&&{B}_{z0}=\sqrt{\displaystyle \frac{2{\mu }_{0}{\rho }_{{\rm{g}}}{c}_{{\rm{s}}}^{2}}{{\beta }_{z}}}\,\quad \mathrm{and}\,\quad {B}_{\phi 0}=\sqrt{\displaystyle \frac{2{\mu }_{0}{\rho }_{{\rm{g}}}{c}_{{\rm{s}}}^{2}}{{\beta }_{\phi }}},\end{eqnarray} \tag{ 8 }$$

respectively, where our units are such that ${M}_{\odot }=1$, ${\mu }_{0}=1$. The radial distance, r, is given in au. We typically quote elapsed time in terms of the orbital period, $t\equiv 2\pi {{\rm{\Omega }}}_{{\rm{K}}}^{-1}({r}_{0})$ at the inner radius ${r}_{0}=1$ of the disk.

The initial plasma-β parameter is set to be constant everywhere, ${{\rm{\Sigma }}}_{0}=1.1\times {10}^{-4}{M}_{\odot }/{r}_{0}^{2}\simeq 980\,{\rm{g}}\,{\mathrm{cm}}^{-2}$ and ${c}_{{\rm{s}}0}\,=0.05\,{r}_{0}{{\rm{\Omega }}}_{{\rm{K}}}({r}_{0})$. We moreover choose $\sigma =1$ as the power-law index for the surface density profile, and $p=1/2$ for the sound speed, yielding a constant aspect ratio of ${h}_{0}\,\equiv H/r=0.05$.

For the purpose of complementing the basic cylindrical model with a radial structure, we assume a radially increasing ionization fraction. In doing so, we stress the importance of maintaining consistency between the coefficients in terms of their dependence on x_e as shown in Equation (4).

In order to compare our model with previous results, we will make use of the Hall length, ${l}_{{\rm{H}}}$, which in this simplified model is normalized by the initial density profile. We assume

$$\begin{eqnarray}&&\displaystyle \frac{{\eta }_{{\rm{H}}}}{| {\boldsymbol{B}}| }={q}_{{\rm{H}}}\displaystyle \frac{{h}_{0}}{\sqrt{{\rho }_{0}}}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{1+w}=\displaystyle \frac{{l}_{{\rm{H}}}}{\sqrt{{\rho }_{{\rm{g}}}}},\end{eqnarray} \tag{ 9 }$$

where ${\rho }_{0}={{\rm{\Sigma }}}_{0}/\sqrt{2\pi }H({r}_{0})$ and ${q}_{{\rm{H}}}={l}_{{\rm{H}}}/H$ at $r={r}_{0}$, that is,

$$\begin{eqnarray}&&{q}_{{\rm{H}}}=\displaystyle \frac{{l}_{{\rm{H}}}({r}_{0})}{H({r}_{0})}.\end{eqnarray} \tag{ 10 }$$

This defines a Hall length proportional to the aspect ratio of the disk, and a diffusion coefficient, ${\eta }_{{\rm{H}}}$, which is a function of the initial density profile. In all our models, we use a vertical domain ${L}_{z}=4H$, then ${q}_{{\rm{H}}}/4={l}_{{\rm{H}}}({r}_{0})/{L}_{z}$. With the ionization fraction an increasing power law of radius, that is, ${x}_{e}\sim {r}^{u}$, with $u\gt 0$, we are looking for a limit on the w exponent to establish a radial profile for the Hall diffusion.

Setting w = 0.5 and using Equation (10) from Kunz & Lesur (2013), we estimate an ionization fraction of

$$\begin{eqnarray}&&{x}_{e}\simeq 3.9\times {10}^{-12}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{1/2},\end{eqnarray} \tag{ 11 }$$

yielding a maximum of ${x}_{e}\sim 8.7\times {10}^{-12}$ at radius r = 5. We remark that we do not use a dedicated ionization model to compute x_e, but attempt to qualitatively describe radial profiles in the disk inner regions (that is, up to two scale heights from the midplane). The value obtained for the ionization fraction is in agreement with those presented by others (Bai & Goodman 2009; Lesur et al. 2014; Keith & Wardle 2015; Béthune et al. 2017), and is at least one order of magnitude below the profile captured by Rodgers-Lee et al. (2016), implying that we are in a stronger Hall regime.

The Hall scheme implemented here is robust in terms of stability when the Hall effect does not strongly dominate the dynamics, i.e., ${\rm{H}}/{\rm{I}}\,\lesssim 100$. As pointed out before (e.g., Huba 2003; Bai 2014), in the case where the Hall term dominates, whistler waves can introduce spurious oscillations, affecting the stability of the numerical scheme.

We have attempted to remedy this issue by applying compact spectral filters (à la Lele 1992) on the edge-centered electric field within the CT step, but the approach has not been able to fully remove the grid-level noise. In the context of solar physics, Vögler et al. (2005) have used a hyperdiffusivity scheme to improve the stability of their ideal-MHD scheme. Adopting a similar strategy, we incorporate an artificial resistivity term that efficiently smooths strong gradients, removing cell-level noise without affecting the large-scale dynamics, as described below.

We note that the artificial resistivity is only required in the Hall-dominated regime. Empirically, we have determined that it is sufficient to apply the artificial resistivity only near a maximum (or minimum) of ${\boldsymbol{B}}$ to obtain a stable evolution, resulting in a parameterization (at the grid location ${{\boldsymbol{x}}}_{i}$) of

$$\begin{eqnarray}&&{\eta }_{\mathrm{art}}={\eta }_{\mathrm{art}}^{0}{\left(\displaystyle \frac{{\displaystyle \sum }_{s}\max (\delta {{\boldsymbol{B}}}_{x}^{(s)},\delta {{\boldsymbol{B}}}_{y}^{(s)},\delta {{\boldsymbol{B}}}_{z}^{(s)})}{| {\boldsymbol{B}}{| }_{\mathrm{ref}}}\right)}^{{q}_{{\rm{a}}}},\end{eqnarray} \tag{ 12 }$$

where ${\eta }_{\mathrm{art}}^{0}$ has the dimension of diffusivity, and s refers to the spatial direction (i.e., x, y or z). For instance, if we take s = y and consider a cell with its center at the indexed position $({\phi }_{i},{r}_{j},{z}_{k})$, we have

$$\begin{eqnarray}\begin{array}{rcl}\delta {{\boldsymbol{B}}}_{x,i+1/2,j,k}^{(y)} & = & (| {B}_{x,i+1/2,j,k}-{B}_{x,i+1/2,j-1,k}| \\ & & +\,| {B}_{x,i+1/2,j+1,k}-{B}_{x,i+1/2,j,k}| \,)\,f(a),\end{array}\end{eqnarray} \tag{ 13 }$$

with similar expressions for $\delta {{\boldsymbol{B}}}_{x}^{(x)}$ and $\delta {{\boldsymbol{B}}}_{x}^{(z)}$, respectively, but of course taking the difference along index i and k. Because we apply the diffusion only over a local maximum (minimum), we have

$$\begin{eqnarray}f(a)\equiv \left\{\begin{array}{ll}0 & \mathrm{if}\ \ a\gt 0\\ \displaystyle \frac{1}{2} & \mathrm{if}\ \ a\leqslant 0\end{array}\right.,\end{eqnarray} \tag{ 14 }$$

with $a\equiv ({B}_{x,i+1/2,j+1,k}-{B}_{x,i+1/2,j,k})$ $\times \,({B}_{x,i+1/2,j,k}-{B}_{x,i+1/2,j-1,k})$, in the case of Equation (13), and correspondingly for the other components.

We find that ${{\eta }_{\mathrm{art}}}_{0}\simeq {10}^{-5}$ is typically sufficient to improve the stability of our numerical simulations. A weighting of ${q}_{{\rm{a}}}=2$ helps to increase the contrast between regions where the artificial diffusion is needed or not, making the artificial dissipation more localized. A seven-point stencil with $| {\boldsymbol{B}}{| }_{\mathrm{ref}}\equiv (| {\boldsymbol{B}}{| }_{i,j,k}+| {\boldsymbol{B}}{| }_{i\pm 1,j,k}+| {\boldsymbol{B}}{| }_{i,j\pm 1,k}+| {\boldsymbol{B}}{| }_{i,j,k\pm 1})/7$ proves to be convenient for the normalization. The artificial resistivity term is also included when evaluating the CFL condition, i.e., we apply and effective Ohmic diffusion which is the sum of the physical and artificial terms, $\eta ={\eta }_{{\rm{O}}}+{\eta }_{\mathrm{art}}$.

The level of dissipation introduced by the artificial resistivity in our simulations can be measured via the magnetic Reynolds number Rm, which we define as

$$\begin{eqnarray}&&{\mathrm{Rm}}_{\mathrm{art}}\equiv \delta {v}_{\mathrm{rms}}\,L/{\eta }_{\mathrm{art}},\end{eqnarray} \tag{ 15 }$$

that is, taking the root mean square of the velocity perturbations, $\delta {v}_{\mathrm{rms}}$, as the characteristic speed. We then compute the value of ${\mathrm{Rm}}_{\mathrm{art}}$ setting L as the vertical cell size $L\,={\rm{\Delta }}z\lesssim r{\rm{\Delta }}\phi \lt {\rm{\Delta }}r$, where $r{\rm{\Delta }}\phi$ and ${\rm{\Delta }}r$ represent the azimuthal and radial cell sizes, respectively. This choice for the length scale returns the minimum possible value for ${\mathrm{Rm}}_{\mathrm{art}}$, providing a robust upper limit to the artificial dissipation.

In general, the artificial resistivity introduces values of ${\mathrm{Rm}}_{\mathrm{art}}$ which are larger than 100 for more than 80% of the active domain, $\geqslant 10$ for 90%, and never smaller than unity. In all of our runs presented in Section 4.2, we find that around 5% of the active domain has a magnetic Reynolds number in the range $1\lt {\mathrm{Rm}}_{\mathrm{art}}\lt 5$. This value drops below 3%, when zonal flows develop. In the case of the runs of Sections 5 and 6, we include a background resistivity profile corresponding to a range of $10\lt \mathrm{Rm}\lt 100$. Introducing artificial resistivity in these models affects around 6% of the active domain with $1\lt {\mathrm{Rm}}_{\mathrm{art}}\lt 5$. We understand that such levels of dissipation have a minimal impact on the dynamics and do not affect the results that we are going to present.

An extra complication is posed by the fact that, the permissible Hall–MHD time step might be significantly smaller than that for ideal MHD. For that reason, we use a sub-cycling scheme that helps to speed up the integration. Both in FARGO3D and Nirvana-iii, the sub-cycling integrates the magnetic field using the HDS scheme⁵ with the Hall–MHD time step, producing an intermediate solution B*. Using B*, we integrate the rest of the induction equation (advection + diffusion) with the time-step constraint from the CFL condition, but without considering the Hall effect. In the strong Hall regime, we find that the time step can be up to a hundred times smaller than the orbital advection time step. Note that we include the FARGO/orbital advection scheme (Masset 2000; Stone & Gardiner 2010) in both codes, reducing the computational integration time and the numerical diffusion. As sub-cycling can potentially affect the propagation of high-frequency waves, we have performed additional convergence tests for oblique waves and the linear growth of the MRI (both local and global modes) with the sub-cycling enabled; all these tests are described in Appendix B. For the oblique wave convergence test, we carried out runs with five and 11 sub-steps, respectively, and we generally obtained a slightly different error convergence rate, that is, up to 20% smaller than that recovered without the sub-cycling. In contrast, the linear growth of the MRI does not appear to be affected by the sub-cycling in our tests. The HDS sub-cycling procedure typically decreases the computational time by a factor between three and five, depending on the adopted physical parameters.

When the sub-cycling scheme is applied for the HDS update, the time step is dominated by the Ohmic resistivity. We then use a super-time-stepping (STS) scheme based on the solution of Alexiades et al. (1996) (both in FARGO3D and Nirvana-iii) and second-order Runge–Kutta–Legendre (RKL2, Meyer et al. 2012) in the case of Nirvana-iii. We do not include a full test of the implementation here, but we have performed runs with and without the scheme for a few specific models, recovering nearly identical results. In Nirvana-iii, we apply Strang splitting for the RKL2 step to maintain second-order accuracy in time. In FARGO3D, we update the magnetic field using the STS scheme after we apply the Hall scheme and before we solve the ideal-MHD induction equation, and the STS scheme follows the implementation described in Simon et al. (2013). The speed-up is within a factor of four in the case of FARGO3D for the more demanding simulations, where Ohmic diffusion induces a time step ten times smaller than the orbital advection step.

The numerical models presented here are all built to be periodic in the vertical and azimuthal directions. We use reflecting radial boundaries for both the velocity and magnetic field. The azimuthal velocity is extrapolated to the Keplerian profile, and we set the vertical velocity such that ${\partial }_{r}{v}_{z}=0$. Moreover, reflecting boundaries are also applied to the azimuthal (${{ \mathcal E }}_{\phi }$) and vertical (${{ \mathcal E }}_{z}$) components of the electromotive force, whereas for the radial component we have ${\partial }_{r}{{ \mathcal E }}_{r}=0$. In combination, these boundary conditions conserve the magnetic and mass flux to machine precision while simultaneously preserving ${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$. On the downside, reflecting boundary conditions induce an accumulation of mass and vertical magnetic flux close to the radial boundaries when the MRI is fully developed. To minimize the amount of magnetic stresses at the radial boundaries, we define buffer zones (with a width of $0.2\,{r}_{0}$) where we apply Ohmic diffusion with a magnitude of ${\eta }_{{\rm{O}}}=3\times {10}^{-4}$. Inside these buffer zones we further define a region of size $0.05\,{r}_{0}$ where we restore the density to the initial density profile ρ. In FARGO3D, following de Val-Borro et al. (2006), the density is modified as $\rho \to (\rho \,\tau +{\bar{\rho }}_{0}{\rm{\Delta }}t)/({\rm{\Delta }}t+\tau )$, where $\tau \propto R(r){{\rm{\Omega }}}^{-1}$, with R(r) a parabolic function that goes to one at the domain boundary and zero at the interior boundary of the wave-damping zone. In Nirvana-iii, we follow a similar procedure using the error function instead.

While the described procedure allows us to reduce the accumulation of mass at the inner boundary—which, if left unchecked, would create undesired large-scale vortices—it has the downside of generating a systematic drainage of the mass density. In a typical simulation, we measure mass losses between 2% and 10% of the initial mass contained in the domain. The amount of mass lost increases when the saturation of the Hall MRI is reached. After that, the mass within the domain remains roughly constant, and we thus do not expect the overall findings to be severely affected.

We discuss below the resolution requirements for the linear growth of the MRI when combining Hall–MHD with Ohmic diffusion. Despite the fact that our simulations are performed in a global cylindrical mesh, we make the simplified assumption that, for every radius, we can compute the local value of the maximum growth rate, ${\gamma }_{{\rm{m}}}$, for the MRI. We furthermore assume a disk with a Keplerian rotation profile, and define ${k}_{{\rm{m}}}$ as the wavenumber with the maximum growth rate.

From linear theory, we obtain the following expression for the fastest growing wavenumber ${k}_{{\rm{m}}}$ (Wardle & Salmeron 2012; Mohandas & Pessah 2017):

$$\begin{eqnarray}&&{k}_{{\rm{m}}}^{2}=\displaystyle \frac{-4{\gamma }_{{\rm{m}}}^{2}({\gamma }_{{\rm{m}}}^{2}+{{\rm{\Omega }}}^{2})}{2{v}_{A}^{2}(2{\gamma }_{{\rm{m}}}^{2}-3{{\rm{\Omega }}}^{2})-({\gamma }_{{\rm{m}}}^{2}+{{\rm{\Omega }}}^{2})[3{\rm{\Omega }}{\eta }_{{\rm{H}}}-4{\gamma }_{{\rm{m}}}{\eta }_{{\rm{O}}}]},\end{eqnarray} \tag{ 16 }$$

and the corresponding maximum growth rate, ${\gamma }_{{\rm{m}}}$, in terms of the Hall diffusivity as

$$\begin{eqnarray}&&{\eta }_{{\rm{H}}}=\displaystyle \frac{24\,{\rm{\Omega }}\,{\eta }_{{\rm{O}}}{\gamma }_{{\rm{m}}}}{9{{\rm{\Omega }}}^{2}-16{\gamma }_{{\rm{m}}}^{2}}\,-\,\displaystyle \frac{2{\rm{\Omega }}{v}_{A}^{2}}{{\gamma }_{{\rm{m}}}^{2}+{{\rm{\Omega }}}^{2}}.\end{eqnarray} \tag{ 17 }$$

If we only consider the Hall effect, Equations (16) and (17) can be combined to yield

$$\begin{eqnarray}&&{k}_{{\rm{m}}}^{2}=\displaystyle \frac{2\,{\rm{\Omega }}}{{\eta }_{{\rm{H}}}}\equiv \displaystyle \frac{2\sqrt{\beta }}{{l}_{{\rm{H}}}H}.\end{eqnarray} \tag{ 18 }$$

This result was obtained in the incompressible limit, but because of the low compressibility of our disk model, we can still establish the connection between the pressure scale height, H, and the wavenumber, ${k}_{{\rm{m}}}$, via Equation (18).

We guarantee a minimum resolution for the local linear instability that in general is around seven cells at r₀. Hawley et al. (1995) suggested a minimum of five cells for codes that use finite difference schemes, such as FARGO3D and eight cells for shock-capturing codes as Nirvana-iii. While the effective resolution increases with radius, the maximum growth rate becomes smaller. For instance, in the models F3D-bz1.4-bp0, and F3D-bz5.3-bp0 (see Section 2.2), the minimum resolution is eight cells at $r={r}_{0}$ for the mode growing at a rate ${\gamma }_{{\rm{m}}}\simeq 0.68$. The resolution is maximal around r = 4.5, where ${\lambda }_{{\rm{m}}}/{\rm{\Delta }}z\simeq 32$, and modes are reasonably resolved down to growth rates of ${\gamma }_{{\rm{m}}}\lesssim 0.06$.

It is well known that the dynamics of the Hall effect in a differentially rotating disk depend on the relative orientation between the angular velocity and the background vertical magnetic field (Wardle & Salmeron 2012; Kunz & Lesur 2013; Bai 2014). For a given range of the strength of the Hall effect, the resulting level of Maxwell stress can be amplified if both quantities are aligned, i.e., ${\boldsymbol{B}}\cdot {\boldsymbol{\Omega }}\gt 0$. This behavior can be traced back to the linear regime, where the additional rotation of the MRI channel solution caused by the Hall current has a destabilizing effect compared to the ideal case (Wardle & Salmeron 2012).

For the following discussion, we adopt the setup described in the Appendix B of Béthune et al. (2016), which was developed with the aim of comparing the results with the previous work by O'Keeffe & Downes (2014). The initial conditions are described in Table 1 (top row) and we run the setup using three different numerical schemes, that is,

• (i)  
  the unsplit HLL-based solver of Tóth et al.,
• (ii)  
  the split HLLD + HDS solver in Nirvana-iii,
• (iii)  
  the MOC + HDS implementation in FARGO3D.

For none of the runs did we make use of the artificial resistivity. As can be checked from Equation (18), the fastest growing Hall-modified MRI mode, with ${\lambda }_{{\rm{m}}}\equiv 2\pi {({H}^{2}{l}_{{\rm{H}}}^{2}/4\beta )}^{1/4}$ is resolved with approximately three cells at r = 4, which is not to the optimal resolution suggested in Section 3.5.

In Figure 1, we show the evolution of the dimensionless stress, α, with time, and we compute the time average, $\langle \alpha {\rangle }_{{\rm{T}}}$, of this quantity between t = 40 and t = 80. We define α as the sum of the $(r,\phi )$-component of the Reynolds and Maxwell stresses, normalized by the gas pressure. We denote these two contributions by ${\alpha }_{{\rm{R}}}$ and ${\alpha }_{{\rm{M}}}$, respectively. We normalize the stress as $\alpha \equiv {N}_{r}^{-1}{\sum }_{j}{\alpha }_{j}$, where

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{j} & \equiv & {\alpha }_{{\rm{R}}}+{\alpha }_{{\rm{M}}},\\ & = & \displaystyle \frac{{\displaystyle \sum }_{i,k}{(\rho {v}_{r}\delta {v}_{\phi })}_{i,j,k}}{{c}_{{\rm{s}}j}^{2}{\displaystyle \sum }_{i,k}{\rho }_{i,j,k}}-\displaystyle \frac{{\displaystyle \sum }_{i,k}{({B}_{r}{B}_{\phi })}_{i,j,k}}{{c}_{{\rm{s}}j}^{2}{\displaystyle \sum }_{i,k}{\rho }_{i,j,k}},\end{array}\end{eqnarray} \tag{ 19 }$$

with indices $(i,j,k)$ as introduced in Section 3.1, and where $\delta {v}_{\phi }$ is the deviation from the mean azimuthal velocity, and N_r denotes the number of cells in the radial direction. Since the quantities in Equation (19) have, in general, different staggering, they need to be interpolated to the r_j coordinate.

Zoom In Zoom Out Reset image size

The results are presented in Table 2, where we can see that the α-values obtained with Nirvana-iii using the HLL solver and FARGO3D only display a slight discrepancy in the average magnetic energy for the aligned configuration, but otherwise agree well.

Table 2.  Measurements of the Mean α-stress and Magnetic Energy

Hall-MHD	$\langle \alpha {\rangle }_{{\rm{T}}+}$	$\langle \alpha {\rangle }_{{\rm{T}}-}$	$\langle {E}_{B}{\rangle }_{{\rm{T}}+}$	$\langle {E}_{B}{\rangle }_{{\rm{T}}-}$
F-3D (HDS)	0.17±0.02	0.066±0.004	0.32±0.03	0.12±0.01
N-iii (HLL)	0.15±0.02	0.067±0.005	0.22±0.02	0.12±0.01
N-iii (HDS)	0.23±0.02	0.12±0.01	0.30±0.04	0.17±0.02
Ideal MHD	$\langle \alpha \rangle$	$\langle {E}_{B}\rangle$
F-3D (HDS)	0.11±0.01	0.19±0.02
N-iii (HDS)	0.16±0.01	0.25±0.03

Note. Averages between t = 40–80 for the different schemes. The magnetic energy, E_B, is normalized to the initial pressure, i.e., ${E}_{B}\equiv {B}^{2}/(2{\rho }_{0}\,{c}_{{\rm{s}}}^{2})$. Brackets $\langle .{\rangle }_{{\rm{T}}}$ indicate time averages and the sign ${\rm{T}}\pm$ refers to the aligned (+) and anti-aligned (−) configuration.

Download table as:  ASCIITypeset image

The $\langle \alpha {\rangle }_{{\rm{T}}}$ obtained with the MOC+HDS (for FARGO3D) and with HLL solver (for Nirvana-iii) are in agreement with those reported by Béthune et al. (2016) but are larger than those reported by O'Keeffe & Downes (2014), where $\langle \alpha {\rangle }_{{\rm{T}}-}\sim 0.0092$ and $\langle \alpha {\rangle }_{{\rm{T}}+}\sim 0.075$ were obtained.

We moreover observe a difference between the timescales in which the MRI is found to saturate. We obtain the saturation level after t = 20 for $\langle \alpha {\rangle }_{{\rm{T}}+}$ and t = 30 for $\langle \alpha {\rangle }_{{\rm{T}}-}$. In contrast, O'Keeffe & Downes (2014) obtained that the stress saturates between 10 and 15 inner orbits. These differences might be related to the use of a different sound speed. In addition, because of the single-fluid approach, our Hall diffusion coefficient only evolves with the magnetic field, whereas in the multifluid approach presented by O'Keeffe & Downes (2014), this coefficient is also updated with the densities and velocities of the charged species. In view of this, the models are likely not directly comparable.

Regarding the magnetic energy evolution, we recover a similar trend as seen by Béthune et al. (2016). The magnetic energy has the same saturation timescale as the α-parameter, and remains stable around a constant mean value between the inner orbits 30 and 90. This is another important difference with respect to O'Keeffe & Downes (2014) where they found a secular growth as a function of time.

The $\langle \alpha {\rangle }_{{\rm{T}}}$ values obtained with Nirvana-iii using the HLLD solver are within 30% larger compared with those from the other numerical schemes. We conjecture that this is simply a consequence of the higher intrinsic accuracy of the HLLD scheme,⁶ since the discrepancy persists for measurements of α in ideal-MHD simulations. Using the same initial conditions as in Table 2, but setting ${l}_{{\rm{H}}}=0$, we found $\langle \alpha {\rangle }_{{\rm{T}}+}=0.16$ with N-iii, and $\langle \alpha {\rangle }_{{\rm{T}}+}=0.11$ for F-3D.

After we have recovered the known dichotomy in the weak Hall regime with our Hall–MHD implementations, we now move on to a regime that involves a much stronger Hall effect. Note that, in the following sections, we exclusively focus on the aligned case, i.e., where ${\boldsymbol{B}}\cdot {\boldsymbol{\Omega }}\gt 0$. This is motivated by our interest in the strong Hall regime, that is, where $| {\eta }_{{\rm{H}}}| \gt 2{v}_{{\rm{A}}}^{2}/{\rm{\Omega }}$. When considering the anti-aligned configuration in this region of parameter space, the Hall effect will always suppress the MRI (Wardle & Salmeron 2012) irrespective of the level of dissipation. Note, however, that for a weaker Hall effect, there exists a limited region of MRI growth in the anti-aligned configuration.

Béthune et al. (2016) performed an extensive study of self-organization by the Hall effect using a cylindrical disk model assuming globally constant initial conditions. They concluded that a strong Hall effect is able to decrease the turbulent transport, in favor of producing large-scale azimuthal zonal flows accumulating vertical magnetic flux. Depending on the plasma-β parameter and the Hall length, these zonal flows can either evolve in the form of axisymmetric rings or vortices. Their exploration of parameter space suggests that the general dynamics is not dramatically altered by the inclusion of Ohmic diffusion or AD, nor by a non-zero toroidal magnetic flux. Béthune et al. (2016) moreover discuss the possibility that the vortices generated might eventually behave as dust traps via establishing regions of super-Keplerian rotation. Their simulations, however, do not include dust directly, and we addressed this topic in Section 6 in the context of radially varying models.

In the following, we will begin by adopting the model B3L6 of Béthune et al. (2016) as a starting reference (see the second line in Table 3). With these initial conditions, the vertical wavenumber ${\lambda }_{{\rm{m}}}$, yielding maximum linear growth for the Hall–MRI, is resolved with 22 cells at r = 2 (also see Section 3.5). The detailed numerical procedure for the setup is described in Section 3.4, above.

Table 3.  Overview of Simulation Runs

	${q}_{{\rm{H}}}$	${\beta }_{\phi }$	${\beta }_{z}$	Bϕ	Bz
				(mG)	(mG)
F3D-bz1.4-bp0	1	⋯	104	0	48.1
F3D-bz1.4-bp0-2	2	⋯	104	0	48.1
F3D-bz1.4-bp0-4	4	⋯	104	0	48.1
F3D-bz5.3-bp0	1	⋯	5 × 103	0	68.0
F3D-bz5.3-bp0-2	2	⋯	5 × 103	0	68.0
F3D-bz5.3-bp0-4	4	⋯	5 × 103	0	68.0
F3D-bz5.3-bp5.3-2	2	5 × 103	5 × 103	68.0	68.0
F3D-bz5.3-bp50-2	2	50	5 × 103	680.	68.0
F3D-bz5.3-bp5.3-4	4	5 × 103	5 × 103	68.0	68.0
F3D-bz5.3-bp50-4	4	50	5 × 103	680.	68.0
F3D-bz1.3-bp0	1	⋯	103	0	152.1
F3D-bz5.3-bp0-fd	1	⋯	5 × 103	0	68.0
F3D-bz5.3-bp50-2-fd	2	50	5 × 103	680.	68.0

Note. Values of ${B}_{\phi }$ and B_z are given at $r=1\,\mathrm{au}$. The models with "f" have an azimuthal domain of 2π and those denoted with a "d" include dust. Model labels state the parameter values used, such as "F3D-bz1.4-bp0-2," which for instance translates to: ${\beta }_{z}={10}^{4}$, ${\beta }_{\phi }=\infty$, ${q}_{{\rm{H}}}=2$.

Download table as:  ASCIITypeset image

In Figure 2, we show the vertical magnetic field at t = 50, that is, before turning on the Hall effect, and at t = 300, when the fluctuations have undergone an extended phase of self-organization. In our case, four bands of vertical magnetic flux are obtained, which is in good agreement with the model B3L6 described by Béthune et al. (2016).

Zoom In Zoom Out Reset image size

4.2.1. Confinement of Magnetic Flux

In order to compare the stress level and the amount of flux confined in the zonal flows, we compute the radial profile $\langle {B}_{z}\rangle$ by taking vertical and azimuthal averages, normalized by the initial value B_z0. We follow the same procedure for the absolute value of the Maxwell stress, $\langle {B}_{r}{B}_{\phi }\rangle$. These quantities are plotted in Figure 3, with the aim of showing how the vertical zonal flow is confined between regions of enhanced Maxwell stress. The values correspond to the time average between the $t=160\mbox{--}260$.

Zoom In Zoom Out Reset image size

This mechanism of confinement was first explained by Kunz & Lesur (2013) in the context of local shearing box simulations and was studied in a cylindrical domain by Béthune et al. (2016). By means of a mean-field theory, they showed that the Hall effect introduces a component proportional to the turbulent Maxwell stress in the evolution equation of the vertical magnetic field. This contribution of the Maxwell stress can be interpreted as an extra dissipation (of either sign) added to the Ohmic diffusion—in other words, it can favor the accumulation of magnetic flux in regions of positive curvature of the mean Maxwell stress. Once the zonal flow is formed, the turbulent field in between starts to flow into the regions of confinement. This generates a region of low stress and vanishing vertical flux between the azimuthal bands. On the other hand, turbulent fluctuations may survive in regions located close to the zonal flows, which can also be seen in Figure 3.

After around 300 orbits, the zonal flows remain quasi-stationary, and three separate zonal flows are clearly distinguishable, with a fourth seeming to appear close the outer boundary at $r\sim 4$. Figure 3 also shows that, on average, the Maxwell stress is larger near the radial buffer zones. Reflecting boundary conditions adopted in the radial direction (enabling conservation of the magnetic flux to machine precision) bring about the accumulation of magnetic flux close to the buffer zones, and this may eventually favor the concentration of Maxwell stress as well.

4.2.2. Effect on the Azimuthal Velocity

We first discuss the runs with fixed ${\beta }_{z}$ and varying ${q}_{{\rm{H}}}$, and then turn to results using three different initial ${\beta }_{z}$, while keeping ${q}_{{\rm{H}}}$ fixed.

The parameters used for models F3D-bz1.4-bp0, F3D-bz1.4-bp0-2 and F3D-bz1.4-bp0-4 are listed in Table 3, along with the other models. For this first fiducial set of simulations, we fix ${\beta }_{z}={10}^{4}$ and we set ${q}_{{\rm{H}}}=1$, 2, and 4, respectively. Zonal flows are recovered in all these configurations. The number of zonal flows that we observe increases as we increase ${q}_{{\rm{H}}}$ from one band to four bands.

We consider two more models where we change the plasma-β parameter to ${\beta }_{z}=5\times {10}^{3}$ and ${\beta }_{z}={10}^{3}$, but we fix ${q}_{{\rm{H}}}=1$. We recover one zonal flow of vertical magnetic flux independently of the initial ${\beta }_{z}$, and in addition to the zonal flow, large-scale vortices show up when we decrease the ${\beta }_{z}$ parameter.

In Figure 5, we plot (in a clockwise sense), the vertical average of ${B}_{z}/{B}_{0}$ together with the the vertical component of the vorticity ${\omega }_{z}\equiv {({\rm{\nabla }}\times {{\boldsymbol{v}}}_{{\rm{g}}})}_{z}$, the Maxwell stress, and the deviation from the Keplerian rotation profile, ${v}_{\phi }/{v}_{{\rm{K}}}-1$. The plasma-β parameter increases from ${\beta }_{z}={10}^{3}$ (top), followed by 5 × 10³ (middle) and finally 10⁴ (bottom panel). For the case of ${\beta }_{z}={10}^{4}$, a concentration of substantial Maxwell stresses is located in the region between the zonal flow and the inner radius. Adjacent to the zonal flow, the Maxwell stress decreases considerably. For the other two cases, it is possible to recognize regions with a lower amount of stress across the vortices. The Maxwell stress decreases as well between the zonal flows, and turbulent fluctuations persist near the regions with enhanced vertical magnetic field.

Zoom In Zoom Out Reset image size

In order to confirm the ability of these models to affect the dust evolution, we show, in the same figure, the vertical average of the deviation from Keplerian rotation, ${v}_{\phi }/{v}_{{\rm{K}}}-1$. In contrast to the setup described in section Section 4.2, the initial conditions here impose an equilibrium rotation profile with sub-Keplerian velocity. This implies that stronger local pressure gradients are need to locally reach super-rotation. However, for all the models, we find super-Keplerian regions with a velocity deviation of 1%, i.e., ${v}_{\phi }/{v}_{{\rm{K}}}-1\simeq 0.01$. These regions are stable for $t\gtrsim 100$.

5.1.1. Generalized Flow Vorticity

In the absence of dissipative terms—that is, viscosity and Ohmic diffusion—the magneto-vorticity ${{\boldsymbol{\omega }}}_{{\rm{M}}}$ satisfies a local conservation equation (Polygiannakis & Moussas 2001; Pandey & Wardle 2008; Kunz & Lesur 2013) if the fluid is incompressible. This quantity is defined as

$$\begin{eqnarray}&&{{\boldsymbol{\omega }}}_{{\rm{M}}}\equiv {\rm{\nabla }}\times {{\boldsymbol{v}}}_{{\rm{g}}}+{\omega }_{{\rm{H}}}\,{\widehat{{\boldsymbol{e}}}}_{{\rm{B}}},\end{eqnarray} \tag{ 21 }$$

where ${\omega }_{{\rm{H}}}\equiv | B| {n}_{e}e/\rho$ is the Hall frequency. In the limit of small Mach numbers, ${{\boldsymbol{\omega }}}_{{\rm{M}}}$ can still be regarded as an approximately conserved quantity.

The local conservation of ${{\boldsymbol{\omega }}}_{{\rm{M}}}$ implies that a concentration of vertical flux has to be anti-correlated with a low-vorticity region. Despite the fact that we include Ohmic diffusion, and the flow is moderately compressible, we do expect to find regions where ${\omega }_{z}$ decreases simultaneously when concentrations of the magnetic flux appear. To test this hypothesis, we compute ${\omega }_{z}$ using the vertical average of the perturbed velocity field. We plot the vorticity in the bottom right sector of the disks shown in Figure 5. By comparing the top left and bottom right sectors, respectively, one can see that negative vorticity regions coincide with the places where the vertical magnetic flux is accumulated—a trend that becomes more prominent for stronger fields, i.e., as ${\beta }_{z}$ decreases.

When including a weak azimuthal net flux, that is, ${\beta }_{\phi }\sim {\beta }_{z}$, Béthune et al. (2016) found that the global picture of self-organization remains intact, despite the presence of a slightly enhanced turbulent activity. We study the effect of a net azimuthal flux on the dynamics by first including an initial ${\beta }_{\phi }$ equal to ${\beta }_{z}=5\times {10}^{3}$ (see row 7 in Table 3), and then increasing the former by a factor of a hundred, that is, ${\beta }_{\phi }=50$ and ${\beta }_{z}=5\times {10}^{3}$ (see row 8 in Table 3).

In Figure 6, we show the vertical magnetic field after 200 orbits for ${\beta }_{z}=5\times {10}^{3}$, ${q}_{{\rm{H}}}=2$ and different values of ${\beta }_{\phi }$. In panels (a) and (b) of Figure 7, we plot the mean radial profile, between t = 100–200, of the vertical magnetic field and the Maxwell stress. When ${\beta }_{\phi }\simeq {\beta }_{z}$, we recover two zonal flows with adjacent regions of turbulent fluctuations. This is in agreement with the field morphology found with ${\beta }_{\phi }=\infty$, where we obtained three zonal flows.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The disk dynamics, however, evolve differently when ${\beta }_{\phi }=50$ is used. Panel (b) of Figure 7 shows a prominent region of Maxwell stress in the outer disk, where large-scale turbulent perturbations dominate the dynamics. The self-organization mechanism is not clearly distinguishable despite the fact that two adjacent zonal flows are recognized. These ring-like structures are radially thinner than the zonal flows obtained with a weak azimuthal field, but are also stable in time. Furthermore, they are correlated with super-Keplerian velocity regions, so again we obtain azimuthally large-scale regions where the radial drift of the dust is slowed down. Figure 6 conveys a clear trend where the inclusion of a net azimuthal field prevents the self-organization, but azimuthally large-scale flux concentrations are still possible to obtain.

We have inspected the Reynolds and Maxwell components associated with the α-parameter introduced in Section 4—see Equation (19). The contribution of the Reynolds stress is below 1% of the total stress, except in the regions where the zonal flows are sustained. In these regions, the Maxwell stress decreases and the local perturbations of the magnetic field are comparable to the velocity turbulent fluctuations. When no azimuthal flux is considered, the Maxwell contribution to the α-parameter is, on average, ∼10⁻⁴. Over the zonal flows it drops to ∼10⁻⁶, whereas the Reynolds contribution has a radially constant mean of ∼10⁻⁶. When ${\beta }_{\phi }=50$, the ratio between the Maxwell and Reynolds contributions is comparable to the case with zero net azimuthal flux, but the amplitude of the perturbations is higher by a factor of ∼100 for both.

The decreased appearance of self-organized structures as a consequence of including an azimuthal field, ${{\boldsymbol{B}}}_{\phi }$, much stronger than ${\boldsymbol{B}}$_z is in agreement with the work of Lesur et al. (2014). They included vertical stratification in a shearing box model, resulting in a strong azimuthal field generated by the Hall effect. Furthermore, they found that the mean azimuthal flux can be 200 times larger than the mean vertical flux, even if the initial azimuthal flux is negligible. A similar result was obtained by Béthune et al. (2017) using a spherical global disk configuration with a magnetized wind in a regime where ${\beta }_{z}={10}^{2}$. The authors showed that AD favors the accumulation of vertical magnetic field into zonal flows. Despite the fact that the Hall effect is negligible (compared with the AD), it can act against the self-organization if the wind can drive a magnetic stress in regions of strong field.

Even though we do not include vertical stratification, we observe that the inclusion of a strong azimuthal flux can alter the dynamics, inhibiting the organization of the zonal flows between channels of strong Maxwell stress and creating a more turbulent flow. However, a sufficient strong Hall diffusion (${q}_{{\rm{H}}}\gtrsim 2$) leads to azimuthally large-scale structures of the vertical magnetic field that generate super-Keplerian velocity regions which are stable for more than 100 orbits. This potentially implies that, even without a clear appearance of self-organized features, it is possible to find flow regions where the dust drift slows down or even halts.

5.2.1. Spectral Energy Distributions

In order to establish a more quantitative comparison between these cases, we compute the spectral energy distribution of the vertical and azimuthal field between radius $r=2$ and 4 for models F3D-bz5.3-bp0-2, -bp5.3-2, and -bp50-2.

In Figure 8, we plot the spectral energy distribution, $m\,| { \mathcal F }({B}_{i}){| }^{2}$, for the azimuthal and vertical field as a function of the azimuthal mode number, m, where the smallest mode is m = 4. For the runs with ${\beta }_{\phi }=5\times {10}^{3}$, the maximum of the distribution is located at large-scale modes, but the peak is clearly reached inside the azimuthal domain $\pi /2$ with an almost flat distribution at lower m. Increasing the azimuthal net flux to ${\beta }_{\phi }=50$ shows a different distribution. The maximum is reached around $m\sim 4$ for the components B_z and ${B}_{\phi }$. This implies that a full 2π domain might be needed in order to correctly capture the maximum of the energy distribution. The energy grows to larger length scales, which can be recognized in Figure 6, where large-scale field perturbations dominate the vertical field structure.

Zoom In Zoom Out Reset image size

The 1D spectrum in Figure 8 allows us to distinguish the trends of the different spectral energy distributions, but it might be affected by the specific choice of the spatial domain. To confirm the differences and obtain a clearer picture, we show in Figure 9 the spectral energy distribution for B_z as a function of radius. To this end, we divide the radial domain into 80 uniformly spaced bins and compute the azimuthal spectrum for each annulus. The bottom panel shows the case with ${\beta }_{\phi }=5\times {10}^{3}$, where it is possible to recognize the location of the two zonal flows. In these regions, the energy is mainly concentrated in low-m modes, and drops as we move outward in the disk. Perturbations are restricted to modes $m\lesssim 40$. The upper panel shows the spectral distribution for ${\beta }_{\phi }=50$, where the energy has a more homogeneous distribution and extends to modes with $m\gtrsim 100$.

Zoom In Zoom Out Reset image size

As a baseline nonlinear test, we perform the shock tube test problem described in O'Sullivan & Downes (2006). The shock propagates under the combined effect of ambipolar and Ohmic diffusion, along with the Hall effect. Because the original problem is formulated for three-fluid MHD, in order to obtain the correct diffusivity coefficients, we have to solve a multifluid system of equations (see O'Sullivan & Downes 2006). Using the multifluid version of FARGO3D, we update the densities of the two charged species via a continuity equation and we solve their velocities, ${{\boldsymbol{v}}}_{{\rm{i}}}$, assuming a steady state in the momentum equation. This yields

$$\begin{eqnarray}&&{ \mathcal E }+({{\boldsymbol{v}}}_{{\rm{i}}}-{{\boldsymbol{v}}}_{{\rm{g}}})\times {\boldsymbol{B}}-\displaystyle \frac{{\rho }_{{\rm{g}}}\,{K}_{i}}{{\alpha }_{i}}({{\boldsymbol{v}}}_{{\rm{i}}}-{{\boldsymbol{v}}}_{{\rm{g}}})=0,\end{eqnarray} \tag{ 27 }$$

where ${{\boldsymbol{v}}}_{{\rm{g}}}$ denotes the velocity of the neutrals, and ${\alpha }_{i}$ is the charge-to-mass ratio of the species. Assuming constant coefficients, K_i, the velocity of each of the charged species can be obtained via

$$\begin{eqnarray}&&({{\boldsymbol{v}}}_{{\rm{i}}}-{{\boldsymbol{v}}}_{{\rm{g}}})(| {\boldsymbol{B}}{| }^{2}+{\xi }^{2})=\displaystyle \frac{({ \mathcal E }\cdot {\boldsymbol{B}})}{\xi }{\boldsymbol{B}}+\xi \,{ \mathcal E }+{ \mathcal E }\times {\boldsymbol{B}},\end{eqnarray} \tag{ 28 }$$

where $\xi ={\rho }_{{\rm{g}}}\,{K}_{i}/{\alpha }_{i}$ and ${ \mathcal E }$ was defined in Equation (2) in Section 2.1. In our case, the diffusion coefficients are ${\eta }_{{\rm{A}}}\equiv {r}_{{\rm{A}}}$, ${\eta }_{{\rm{H}}}\equiv {r}_{{\rm{H}}}$ and ${\eta }_{{\rm{O}}}\equiv {r}_{{\rm{O}}}$, where ${r}_{{\rm{A}}}$, ${r}_{{\rm{H}}}$ and ${r}_{{\rm{O}}}$ are computed using Equations (10) and (11) in O'Sullivan & Downes (2006). The initial conditions are such that the left and right states are given by

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{{\rm{L}}} & = & 1.7942,\,\,{v}_{z{\rm{L}}}=-0.9759,\\ {v}_{y{\rm{L}}} & = & -0.6561,\,\,{B}_{y{\rm{L}}}=1.74885,\\ {\rho }_{1{\rm{L}}} & = & 8.9712\times {10}^{-8},\,\,{\rho }_{2{\rm{L}}}=1.7942\times {10}^{-3},\end{array}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{{\rm{R}}} & = & 1.0,\,\,{v}_{z{\rm{R}}}=-1.751,\\ {v}_{y{\rm{R}}} & = & 0.0,\,\,{B}_{y{\rm{R}}}=0.6,\\ {\rho }_{1{\rm{R}}} & = & 5\times {10}^{-8},\,{\rho }_{2{\rm{R}}}=1\times {10}^{-3},\end{array}\end{eqnarray} \tag{ 30 }$$

respectively. The sound speed is ${c}_{{\rm{s}}}=0.1$ and B_z = 1.0 for both states, while ${v}_{x}={B}_{x}=0.0$. The coefficients for the charged species are

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }_{1} & = & -2\times {10}^{9},\,\,{\alpha }_{2}=1\times {10}^{5},\\ {K}_{1} & = & 4\times {10}^{2},\,\,{K}_{2}=2.5\times {10}^{6}.\end{array}\end{eqnarray} \tag{ 31 }$$

The shock propagates along the z direction in a domain of unit size, using 500 and 1000 cells, respectively. The initial discontinuity is located at z = 0.25, and the boundary conditions are fixed at the initial state. In Figure 12, we show the numerical solution at t = 2.7. The L1 error for a grid resolution of ${\rm{\Delta }}z=2\times {10}^{-2}$ is ${e}_{1}=3.66\times {10}^{-2}$ and with a resolution of ${\rm{\Delta }}z=1\times {10}^{-3}$, it is ${e}_{1}=9.8\times {10}^{-3}$, giving the scaling ${e}_{1}\sim {\rm{\Delta }}{z}^{1.9}$ which is close to the expected second-order convergence.

Zoom In Zoom Out Reset image size

In this section, we consider the linearized Hall–MHD equations of an incompressible fluid with only the Lorentz force in the momentum equation. Let us assume a background magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}{\hat{x}}_{1}$ subject to the fluctuations $\delta {\boldsymbol{B}}=(0,\delta {B}_{2},\delta {B}_{3})$ and a perturbed velocity field $\delta {\boldsymbol{v}}=(0,\delta {v}_{2},\delta {v}_{3})$. All perturbations have the form $\delta f=\delta {f}_{0}\exp [i\,(\omega t-{\boldsymbol{k}}\cdot {\boldsymbol{x}})]$ along the direction x₁, that is, ${\boldsymbol{k}}=\pm k{\hat{x}}_{1}$. It can be shown that

$$\begin{eqnarray}&&\delta {v}_{\mathrm{2,3}}=-\displaystyle \frac{{\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}}{{\mu }_{0}{\rho }_{0}\omega }\,\delta {B}_{\mathrm{2,3}},\end{eqnarray} \tag{ 32 }$$

where $\delta {B}_{\mathrm{2,3}}$ is the non-trivial solution of the system $A\,\delta B=0$, with A being the matrix

$$\begin{eqnarray}A=\left[\begin{array}{cc}i(-{\omega }^{2}+{v}_{A}^{2}{k}^{2}) & -\omega \,{\eta }_{{\rm{H}}}{k}^{2}\\ \omega \,{\eta }_{{\rm{H}}}{k}^{2} & i(-{\omega }^{2}+{v}_{A}^{2}{k}^{2})\end{array}\right],\end{eqnarray} \tag{ 33 }$$

where ${v}_{A}={B}_{0}/\sqrt{{\mu }_{0}{\rho }_{0}}$ is the Alfvén speed, and ${\eta }_{{\rm{H}}}$ is considered constant and normalized by B₀. The system has four independent solutions, which correspond to two circularly polarized waves propagating along $\pm {x}_{1}$. For non-zero ${\eta }_{{\rm{H}}}$, these waves are commonly know as the whistler mode and the ion-cyclotron mode, respectively. Characteristically, these are dispersive waves, that is, the phase velocity, $\omega /k$, depends on the wavenumber. More precisely,

$$\begin{eqnarray}&&\displaystyle \frac{\omega }{k}=\displaystyle \frac{{\eta }_{{\rm{H}}}}{2}({\boldsymbol{k}}\cdot {\widehat{{\boldsymbol{e}}}}_{{\rm{B}}})\pm {\left[\displaystyle \frac{1}{4}{(k{\eta }_{{\rm{H}}})}^{2}+{v}_{A}^{2}\right]}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 34 }$$

with ${\widehat{{\boldsymbol{e}}}}_{{\rm{B}}}$ the unit vector in the direction of B₀.

To test the newly developed Hall–MHD modules, we perform a convergence study of an oblique wave propagation following a similar ansatz as the one presented in Gardiner & Stone (2008) for circularly polarized waves in the ideal-MHD regime. We define a Cartesian periodic box with dimensions $(3,1.5,1.5)$, resolved by $(2N,N,N)$ grid cells in x, y, and z, respectively. The propagation is along the oblique coordinate $x1={k}_{x}x+{k}_{y}y+{k}_{z}z$, with

$$\begin{eqnarray}\begin{array}{rcl}{k}_{x} & = & 2\pi \cos (\alpha )\cos (\beta ),\,\,{k}_{y}=2\pi \cos (\alpha )\sin (\beta ),\,\,\mathrm{and}\\ {k}_{z} & = & 2\pi \cos (\alpha ).\end{array}\end{eqnarray} \tag{ 35 }$$

We moreover set $\alpha =2/3$ and $\beta =2/\sqrt{5}$, in order to fit one wavelength $\lambda =1$ inside the box. The initial velocity field is ${\boldsymbol{v}}=(0,{10}^{-6}\cos (2\pi {x}_{1}),{10}^{-6}\sin (2\pi {x}_{1}))$ and the initial magnetic field is computed via a vector potential in such a way that ${\boldsymbol{B}}$ satisfies Equation (32) with a background field ${{\boldsymbol{B}}}_{0}=1.0{\hat{x}}_{1}$. The gas is treated as isothermal where ${P}_{0}={c}_{{\rm{s}}}^{2}{\rho }_{0}$, with an initial density ${\rho }_{0}=1.0$ and gas pressure P₀ = 1.0. The Hall diffusion is set to ${\eta }_{{\rm{H}}}=0.1$.

We compute the L1 error for the centered components of the magnetic field as defined in Equation (70) of Gardiner & Stone (2008), and we plot the results for different resolutions in Figure 13. For the whistler mode Nirvana-iii returns an error ${e}_{1}\sim {{\rm{\Delta }}}_{z}^{2.3}$ and ${e}_{1}\sim {{\rm{\Delta }}}_{z}^{1.94}$, with the HLLD and HLL solvers, respectively. This scaling is different for the ion cyclotron mode, where ${e}_{1}\sim {{\rm{\Delta }}}_{z}^{1.87}$ with HLLD solver and ${e}_{1}\sim {{\rm{\Delta }}}_{z}^{2.1}$ with the HLL scheme. In the case of FARGO3D, the error goes as ${e}_{1}\sim {\rm{\Delta }}{z}^{1.94}$ for the whistler mode and ${e}_{1}\sim {\rm{\Delta }}{z}^{1.87}$ for the ion cyclotron mode. All the results are reasonably close to the expected second-order convergence of the implemented schemes.

Zoom In Zoom Out Reset image size

We now study the growth rate of linear MRI modes under the Hall effect. We use an axisymmetric, 2D cylindrical domain with a grid of 64 × 64 cells, where the radial domain is fixed to $r\in [0.8,1.2]$. We adopt a strategy where we initialize simulations such that the vertical box size matches the wavelength of a given individual mode. In this way, we run one simulation for each mode in a box with ${L}_{z}=\lambda \equiv 2\pi /{k}_{z}$, maintaining the same effective resolution for all the modes.

We use a similar initial condition as Sano & Stone (2002), where the plasma-β parameter is set to $\beta =800$, ${\rho }_{0}=1.0$, ${c}_{{\rm{s}}}=0.1$, and we apply a Keplerian background velocity field. We then add perturbations to the azimuthal and radial velocities of the form ${v}_{0}\cos ({kz})$ and ${v}_{0}\sin ({kz})$, respectively, where ${v}_{0}={10}^{-6}{c}_{{\rm{s}}}$. The Hall diffusivity coefficient, ${\eta }_{{\rm{H}}}$, is defined in such a way that the Elsässer numbers are ${{\rm{\Lambda }}}_{{\rm{H}}}=1$ when ${\boldsymbol{B}}\cdot {\boldsymbol{\Omega }}\gt 0$ and ${{\rm{\Lambda }}}_{{\rm{H}}}=-1$ when ${\boldsymbol{B}}\cdot {\boldsymbol{\Omega }}\lt 0$. We integrate for two orbits at r = 1, and measure the growth rate between $t=1.5-1.9$ orbits at the same radius. By means of an exponential fit to the maximum value of the radial magnetic field at r = 1, we obtain the numerical solutions shown in the left panel of Figure 14, which agree with the analytic solution to within 3% for the critical mode and 0.02% for the maximum growth rate. We attribute deviations to the discrepancy between the local Cartesian approximation and the cylindrical computational setup.

Zoom In Zoom Out Reset image size

The eigenvectors that we show are the solutions of the linear Hall-MRI equation in a cylindrical coordinate system assuming axisymmetric perturbations. The perturbed MHD system of equations is linear but not fully algebraic. Thus, to compute the semi-analytic solutions we use a spectral code that expands the radial velocity and magnetic field using Chebyshev polynomials. The numerical setup is initialized using random white noise in the velocity field of the order ${10}^{-8}{c}_{{\rm{s}}}$ in a global isothermal disk. We adopt a box of size $({L}_{r},{L}_{z})=[1,1]$ with $256\times 256$ cells, with periodic boundaries in the vertical and azimuthal directions. The radial boundary conditions we apply are ${v}_{r}={b}_{r}={b}_{\phi }={\partial }_{r}{b}_{z}={\partial }_{r}{v}_{z}=0$. The disk is assumed to rotate with a Keplerian profile, and the initial conditions are

$$\begin{eqnarray}\begin{array}{rcl}\beta & = & 31250,\quad \rho =1.0,\quad {c}_{{\rm{s}}}=0.25,\\ {\eta }_{{\rm{O}}} & = & 0.003,\quad \mathrm{and}\quad {l}_{{\rm{H}}}=1.0.\end{array}\end{eqnarray} \tag{ 36 }$$

In the right panel of Figure 14, we plot the eigenvectors obtained after 37 orbits. We also compute the growth rate of the modes via a linear fit of $\mathrm{log}({B}_{r}(t))$. We obtained a growth rate $\gamma =0.0989$ with F-3D, $\gamma =0.0989$ ($\gamma =0.0988$) with N-iii HLLD (HLL). All the values are in excellent agreement with the expected analytic value of ${\gamma }_{0}=0.0989$ (see also Béthune et al. 2016).