OPEN ACCESS

<figure><img src="image_1.png" /></figure>

CrossMark

# Magnetically Driven Turbulence in the Inner Regions of Protoplanetary Disks

David G. Rea<fnref n="1" />, Jacob B. Simon<fnref n="1" />, Daniel Carrera<fnref n="1" />, Geoffroy Lesur<fnref n="2" />, Wladimir Lyra<fnref n="3" />, Debanjan Sengupta<fnref n="3" />, Chao-Chin Yang (楊朝欽)<fnref n="4" />, and Andrew N. Youdin<fnref n="5" />

Received 2024 April 4; revised 2024 June 7; accepted 2024 June 11; published 2024 August 29

## Abstract

Given the important role turbulence plays in the settling and growth of dust grains in protoplanetary disks, it is crucial that we determine whether these disks are turbulent and to what extent. Protoplanetary disks are weakly ionized near the midplane, which has led to a paradigm in which largely laminar magnetic field structures prevail deeper in the disk, with angular momentum being transported via magnetically launched winds. Yet, there has been little exploration of the precise behavior of the gas within the bulk of the disk. We carry out 3D, local shearing box simulations that include all three low-ionization effects (ohmic diffusion, ambipolar diffusion, and the Hall effect) to probe the nature of magnetically driven gas dynamics 1–30 au from the central star. We find that gas turbulence can persist with a generous yet physically motivated ionization prescription (order unity Elsässer numbers). The gas velocity fluctuations range from 0.03 to 0.09 of the sound speed $c_s$ at the disk midplane to ~$c_s$ near the disk surface, and are dependent on the initial magnetic field strength. However, the turbulent velocities do not appear to be strongly dependent on the field polarity, and thus appear to be insensitive to the Hall effect. The midplane turbulence has the potential to drive dust grains to collision velocities exceeding their fragmentation limit, and likely reduces the efficacy of particle clumping in the midplane, though it remains to be seen if this level of turbulence persists in disks with lower ionization levels.

*Unified Astronomy Thesaurus concepts: Planet formation (1241); Magnetohydrodynamics (1964); Accretion (14); Protoplanetary disks (1300); Magnetohydrodynamical simulations (1966)*

## 1. Introduction

Whether gas motions in protoplanetary disks are turbulent or laminar has significant implications for numerous steps in the planet formation process. For instance, if present, strong turbulence can cause large collisional velocities between dust grains, hindering their growth via fragmentation or bouncing (Ormel & Cuzzi 2007; Güttler et al. 2010; Kothe 2016). Turbulence may concentrate particles in turbulent eddies and facilitate the first stages of planetesimal formation (Johansen et al. 2007; Cuzzi et al. 2008; Hartlep & Cuzzi 2020) though it is also possible for turbulence to prevent or slow down the formation of planetesimals driven by the streaming instability (e.g., Chen & Lin 2020; Gole et al. 2020; Umurhan et al. 2020; Lim et al. 2024). Turbulence may also play an important role in the migration of protoplanets, by generating stochastic net torques that result in a random walk migration, as opposed to Type I monotonic inward drift (Nelson & Papaloizou 2004; Johnson et al. 2006; Yang et al. 2009, 2012; Paardekooper et al. 2011; Stoll et al. 2017).

There are a number of mechanisms that may give rise to turbulence in disks, ranging from purely hydrodynamic (i.e., not requiring a magnetic field) mechanisms (see Lyra & Umurhan 2019) to those induced by magnetic fields. The most studied mechanism in the latter category is magnetorotational instability (MRI; Balbus & Hawley 1998). The MRI arises

naturally in sufficiently ionized Keplerian disks with a (relatively weak) magnetic field present and gives rise to correlated turbulent fluctuations in the gas velocity and magnetic field that produce an outward flux of angular momentum. However, most of the gas in protoplanetary disks is cold and poorly ionized, and as such, much of the disk is dominated by low-ionization physics (e.g., Gammie 1996; Wardle 2007): ohmic diffusion (where electrons are collisionally impeded by neutral species, i.e., low gas conductivity), the Hall effect (current induced by the drift velocity between electrons and ions) and ambipolar diffusion (drift between the bulk neutral gas and ions). At a given magnetic field strength, gas becomes more decoupled from the field as density increases, i.e., increasing frequency of collisions between electrons, ions, and neutral species. In order of increasing gas density, the neutrals decouple (ambipolar diffusion), then ions (Hall effect), and then electrons (ohmic dissipation). The importance of these effects has been studied extensively by both analytic (e.g., Gammie 1996; Balbus & Terquem 2001; Kunz & Balbus 2004; Desch 2004; Salmeron & Wardle 2005) and numerical works (e.g., Sano & Stone 2002; Fleming & Stone 2003; Lesur et al. 2014; Bai 2014, 2015; Simon et al. 2015b, 2018; Gressel et al. 2015; Bethune et al. 2017; Cui & Bai 2021).

Ohmic and ambipolar diffusion generally weaken MRI-driven turbulence. At large distances from the central star (~30–100 au, though the exact distance is sensitive to disk properties), a layered structure is thought to exist in which MRI-active surface regions ionized by far-ultraviolet (FUV) photons (Perez-Becker & Chiang 2011) surround a midplane

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Physics and Astronomy, Iowa State University, Ames, IA, 50010, USA; drea1@iastate.edu</footnote>

<footnote n="2">²Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³New Mexico State University, Department of Astronomy, P.O. Box 30001 MSC 4500, Las Cruces, NM 88003-8001, USA</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, The University of Alabama, Box 870324, Tuscaloosa, AL 35487-0324, USA</footnote>

<footnote n="5">⁵University of Arizona, Steward Observatory and LPL, Tucson, AZ 85719, USA</footnote>
region in which ambipolar diffusion dampens MRI-driven turbulence (Simon et al. 2013a, 2013b; Cui & Bai 2021). A similar structure has long been thought to be at work in the inner regions of disks but with ohmic diffusion quenching the midplane region (Gammie 1996; Fleming & Stone 2003). This classic dead zone paradigm, which preceded the layered model for the outer disk, has been modified significantly over recent years with the appreciation that the Hall effect fundamentally alters the gas dynamics at intermediate radii (1–30 au, Wardle 2007).

When the Hall effect is important, the orientation of the vertical magnetic field plays a role in determining the ensuing gas dynamics (Wardle & Ng 1999; Balbus & Terquem 2001; Kunz 2008; Bai 2015; Simon et al. 2015b). In particular, the Hall shear instability (HSI; Kunz 2008) occurs when electro-magnetic (whistler) waves couple to the background shear. The HSI is typically considered to be active only when the magnetic field points in the same direction as the gas rotation vector; more generally, it is strongest when the gas vorticity, wavevector, and magnetic field share a mutual axis. Considering only the vertical field **B** and the background Keplerian shear $\Omega$ (which typically shares an axis with the gas vorticity in disks), the coupling parameter can be expressed as $(\Omega \cdot \mathbf{B})$. The Hall effect can cause strong, persistent, and highly laminar magnetic torques, which act on the gas to transport angular momentum (Lesur et al. 2014; Bai 2015; Simon et al. 2015b) if the vertical component of the magnetic field is aligned with the angular momentum vector of the rotating disk gas $(\Omega \cdot \mathbf{B} > 0)$, and nonaxisymmetric modes of the HSI can produce intermittent bursts of magnetic torque (Simon et al. 2015b) even when the magnetic field is anti-aligned with the disk rotation $(\Omega \cdot \mathbf{B} < 0)$, though these bursts can be quenched by stronger ambipolar diffusion (Bai 2015; Simon et al. 2015b).

Another source of angular momentum transport in protoplanetary disks is magnetothermal winds (Bai et al. 2016; Bai 2017; Bethune et al. 2017); in this process, gas is heated and ionized in the upper disk layers. This ionization couples the gas to the magnetic field, which then launches material away from the disk along magnetic field lines extending beyond the disk surface, transporting angular momentum in the process as in the Blandford–Payne mechanism (Blandford & Payne 1982). The importance of this mechanism in angular momentum transport, the critical role that nonideal MHD effects play in regulating magnetic activity, and the (resulting) largely laminar magnetic fields within the disk (e.g., Lesur et al. 2014) have led to a new paradigm where accretion is assumed to be a largely laminar process, in contrast to the vigorous turbulence induced by the MRI.

Despite the dominance of largely laminar magnetic stresses in angular momentum transport, it is still not well understood how exactly the gas behaves in the inner disk (<30 au) where the Hall effect is dominant. Indeed, as shown by Simon et al. (2018), at ~30–100 au scales, the gas can be rendered turbulent for a large range of parameter values, even if magnetothermal winds dominate the angular momentum transport. This was found to be caused by zonal flows that create regions both stable and unstable to the MRI; the latter regions generate turbulence that then propagates throughout the domain (see also Cui & Bai 2021, which confirms the presence of outer disk turbulence with global simulations).

Here, we ask whether similar processes can be at work in the inner disk. While the results of Bai (2015) suggest (see their

Figure 6, for instance) that there may be some form of gas turbulence at these smaller radii, there are three questions that still need to be addressed in order to elucidate the nature of gas motions in the inner disk:

1. What is the strength of nonlaminar gas motions (if present) in the 1–30 au regions of protoplanetary disks?

2. Do these gas motions resemble disk turbulence as seen in earlier simulations such as in Simon et al. (2012) or Simon et al. (2015a)?

3. What is the origin of these gas motions?

In this paper, we will address the first two questions, with a focus on characterizing gas motions as a function of the relevant parameters (e.g., magnetic field strength and orientation). We carry out a series of 3D, local shearing box simulations that employ all three nonideal MHD terms; while such a setup cannot accurately capture magnetothermal winds (e.g., Fromang et al. 2013), our focus is on the gas motions themselves, particularly near the midplane where planet formation occurs.

This paper is organized as follows: In Section 2, we describe our simulations and methods in detail. In Section 3, we present the magnetic field and velocity structure from our fiducial simulation. In Section 4, we examine the remainder of our parameter space and investigate the possible sources of turbulent gas motions. In Section 5, we summarize our conclusions and comment on their implications for protoplanetary disk evolution and planetesimal formation.

## 2. Method

### 2.1. Shearing Box

We carry out our numerical simulations using Athena, a second-order accurate dimensionally unsplit Godunov code for solving the equations of MHD (Stone et al. 2008). We use the HLLD Reimann solver of Miyoshi & Kusano (2005) with the piecewise-parabolic method (Colella & Woodward 1984) for third-order spatial reconstruction and the corner-transport-upwind method (Colella 1990) time integrator. The constraint $\nabla \cdot \mathbf{B} = 0$ is maintained to machine precision via the constrained transport algorithm (Evans & Hawley 1988). Nonideal MHD terms are implemented using the methods of Simon & Hawley (2009) for ohmic resistivity, Bai & Stone (2011) for ambipolar diffusion, and Bai (2014) for the Hall effect.

We use the local shearing box approximation (e.g., Hawley et al. 1995; Stone & Gardiner 2010) for our computational domain in order to resolve turbulence on scales $\ll H$ (where $H$ is the gas scale height). The shearing box is a small patch of the disk at a fiducial radius $R_0$ and co-orbiting the central star with Keplerian angular velocity $\Omega = \Omega_0 \dot{e}_z$. The shearing box has Cartesian coordinates $(x, y, z) = (R - R_0, R_0\phi, z)$, where $R$ and $\phi$ correspond to the radial and azimuthal directions of a cylindrical coordinate system, respectively. To account for differential rotation, there is a linear background shear $v_0 = -q\Omega_0 x \dot{e}_y$, where $q$ is the shear parameter:

$$q = - \frac{d \ln \Omega}{d \ln R} \bigg|_{R=R_0} . \qquad (1)$$

The shear parameter $q = 3/2$ for a Keplerian disk.

The shearing box boundaries are periodic in $y$, shear periodic in $x$ (Hawley et al. 1995; Stone & Gardiner 2010), and a
modified outflow boundary in z (Simon et al. 2011b). Shearing box simulations of the MRI in vertically stratified disks exhibit significant mass outflow through the vertical domain boundaries (Suzuki & Inutsuka 2009; Simon et al. 2013a); we continuously replenish disk mass to compensate for mass loss as described in Simon et al. (2015b). Finally, we employ the orbital advection algorithm as described in Stone & Gardiner (2010) to subtract off the Keplerian shear and integrate this component separately. In units such that the magnetic permeability is unity, the equations of MHD are the continuity equation

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (2) $$

the momentum equation

$$ \frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \mathbf{v} - BB) + \nabla \left( P + \frac{1}{2} B^2 \right) = 2 \rho q \Omega_0^2 x - \rho \Omega_0^2 z - 2 \Omega_0 \hat{e}_z \times \rho \mathbf{v}, \qquad (3) $$

and the induction equation

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times (\mathbf{v} \times \mathbf{B}) = -\nabla \times \left[ \eta_O \mathbf{J} + \eta_H \frac{\mathbf{J} \times \mathbf{B}}{B} - \eta_A \frac{(\mathbf{J} \times \mathbf{B}) \times \mathbf{B}}{B^2} \right]. \qquad (4) $$

Here, $\rho$ is the mass density, $\mathbf{v}$ is the gas velocity, $\mathbf{B}$ is the magnetic field, and $\mathbf{J} = \nabla \times \mathbf{B}$ is the current density. The magnetic field is normalized by a factor $\sqrt{4\pi}$. From left to right, the right-hand side terms of the momentum equation represent radial tidal forces, vertical gravity, and the Coriolis force. We assume an isothermal equation of state with gas pressure $P = \rho c_s^2$, where $c_s$ is the isothermal sound speed $\sqrt{k_B T / \mu m_p}$, with Boltzmann constant $k_B$, disk temperature $T$, proton mass $m_p$, and mean molecular weight $\mu = 2.33$ (see Section 2.2 for specific values of these quantities). The parameters $\eta_O$, $\eta_H$, and $\eta_A$ are the ohmic, Hall, and ambipolar diffusion coefficients, respectively. When the only charged species are electrons and ions (i.e., in the absence of dust grains) and ignoring the effect of disk chemistry on ionization, these are (Balbus & Terquem 2001; Wardle 2007):

$$ \eta_O = \frac{c_s^2 m_e}{4\pi e^2} \frac{n}{n_e} \langle \sigma \mathbf{v} \rangle_e, \qquad (5) $$

$$ \eta_H = \frac{B_c}{\sqrt{4\pi} e n_e}, \qquad (6) $$

$$ \eta_A = \frac{B^2}{\gamma_i \rho \rho_i}, \qquad (7) $$

where $m_e$ is the electron mass, $n$ and $n_e$ are the number densities of the neutrals and the electrons, respectively, and $\rho_i$ is the ion mass density.

$$ \langle \sigma \mathbf{v} \rangle_e = 8.28 \times 10^{-9} \left( \frac{T}{100 \text{ K}} \right)^{1/2} \text{cm}^3 \text{s}^{-1} \qquad (8) $$

is the electron-neutral collision rate (Draine et al. 1983), $\gamma_i = \langle \sigma \mathbf{v} \rangle_i / (m_n + m_i)$, and

$$ \langle \sigma \mathbf{v} \rangle_i = 1.3 \times 10^{-9} \text{cm}^3 \text{s}^{-1} \qquad (9) $$

<table><caption>**Table 1**
Disk Model Parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>M<sub>star</sub></td><td>1.0 M<sub>☉</sub></td></tr><tr><td>M<sub>disk</sub></td><td>0.05 M<sub>☉</sub></td></tr><tr><td>R<sub>c</sub></td><td>100 au</td></tr><tr><td>γ</td><td>1.0</td></tr><tr><td>T<sub>X</sub></td><td>3.0 keV</td></tr><tr><td>L<sub>X</sub></td><td>10<sup>30</sup> erg s<sup>-1</sup></td></tr><tr><td>Σ<sub>FUV</sub></td><td>0.005 g cm<sup>-2</sup></td></tr></tbody></table>

is the ion-neutral collision rate (Bruce 2011). Here, $m_i$ and $m_n$ are the masses of the ions and neutrals, respectively.

## 2.2. Disk Model

We center our local simulations at specific radii in a model disk, which is described as follows. We consider an exponentially tapered surface density profile as in Lynden-Bell & Pringle (1974) and Hartmann et al. (1998),

$$ \Sigma(R) = \frac{M_{\text{disk}}(2 - \gamma)}{2\pi R_c^2} \times \left(\frac{R}{R_c}\right)^{-\gamma} \exp\left[-\left(\frac{R}{R_c}\right)^{2-\gamma}\right] \text{g cm}^{-2}, \qquad (10) $$

where $M_{\text{disk}}$ is the disk mass, $R_c$ is the characteristic radius of the disk, and $\gamma$ is the surface density power law index. The values of these model parameters are chosen to be representative of the disks observed by Andrews et al. (2009), and can be found in Table 1. The temperature profile is

$$ T(R) = 280 \left( \frac{R}{1 \text{ au}} \right)^{-1/2} \text{K}. \qquad (11) $$

Assuming the disk is vertically isothermal at each radius, vertical hydrostatic equilibrium requires a density profile

$$ \rho(x, y, z) = \rho(R_0) \exp\left(-\frac{z^2}{2H^2}\right) \qquad (12) $$

with scale height $H = c_s / \Omega_0$.

## 2.3. Ionization

We include contributions to the ionization rate from radioactive decay with $\zeta_{\text{RD}} = 10^{-19} \text{s}^{-1}$ (Umebayashi & Nakano 2009) and cosmic rays with $\zeta_{\text{CR}} = \zeta_{\text{CR},0} \exp(-\Sigma/96 \text{ g cm}^{-2}) \text{s}^{-1}$ (Umebayashi & Nakano 1981). The cosmic-ray flux ($\zeta_{\text{CR},0}$) is highly uncertain, but observations in the direction of $\zeta$ Persei suggest $\zeta_{\text{CR},0} \sim 10^{-16} \text{s}^{-1}$ (McCall et al. 2003), which we take to be our fiducial value in order to make contact with previous work (e.g., Lesur et al. 2014), though we note that some modeling predicts $\zeta_{\text{CR},0} \lesssim 10^{-19} \text{s}^{-1}$ in TW Hya (Cleeves et al. 2015). We include the contribution from stellar X-rays with a luminosity $L_X = 10^{30} \text{ erg s}^{-1}$ and photon energy $T_X = 3.0 \text{ keV}$ (Igea &
<figure><img src="image_3.png" /><figcaption>Figure 1. Initial ionization fraction $x_e$ as a function of height. The disk ionization can change in time as it is dependent on the density of the neutrals (Equation (15)). The ionization near the disk surface ($|z| \gtrsim 3H$) is dominated by FUV photons and by cosmic rays near the disk midplane.</figcaption></figure>

Glassgold 1999; Bai & Goodman 2009) as

$$ \zeta_{\text{XR}} = \left( \frac{L_X}{10^{29} \text{ erg s}^{-1}} \right) \left( \frac{R}{1 \text{ au}} \right)^{-2.2} \times \left( \zeta_1 [e^{-(N_{\text{H1}}/N_1)^a} + e^{-(N_{\text{H2}}/N_1)^a}] + \zeta_2 [e^{-(N_{\text{H1}}/N_2)^b} + e^{-(N_{\text{H2}}/N_2)^b}] \right), \quad (13) $$

where $N_{H1,2}$ is the column density of hydrogen integrating through the disk from above ($N_{H1}$) and from below ($N_{H2}$), and for $T_X = 3.0$ keV, $\zeta_1 = 6 \times 10^{-12} \text{ s}^{-1}$, $N_1 = 1.5 \times 10^{21} \text{ cm}^{-2}$, $a = 0.4$, $\zeta_2 = 1.0 \times 10^{-15} \text{ s}^{-1}$, $N_2 = 7.0 \times 10^{23} \text{ cm}^{-2}$, and $b = 0.65$. The X-ray ionization rate is only weakly sensitive to $T_X$ (Igea & Glassgold 1999).

We follow Bai (2014) and treat FUV ionization as another independent ionization source with

$$ \zeta_{\text{FUV}} = 10^{-6} \left(\frac{R}{1 \text{ au}}\right)^{-2} \exp\left(-\frac{\Sigma}{\Sigma_{\text{FUV}}}\right) \text{s}^{-1} \quad (14) $$

and $\Sigma_{FUV} = 0.005 \text{ g cm}^{-2}$. The ionization fraction $x_e$ is obtained by balancing these ionization sources with dissociative recombination in the absence of metals and dust (Gammie 1996; Fromang et al. 2002):

$$ x_e = g \frac{n_e}{n_n} = g \sqrt{\frac{\zeta_{\text{RD}} + \zeta_{\text{CR}} + \zeta_{\text{XR}} + \zeta_{\text{FUV}}}{n_n \alpha_{dr}}}, \quad (15) $$

where $\alpha_{dr} = 3 \times 10^{-6} (T/\text{K})^{-1/2} \text{ cm}^3 \text{ s}^{-1}$ is the dissociative recombination rate coefficient. Similar to Bai (2017), we multiply the ionization fraction by a term $g$, where

$$ g = \exp\left[\frac{0.3\Sigma_{\text{FUV}}}{\Sigma + 0.01\Sigma_{\text{FUV}}}\right]. \quad (16) $$

Without this term, ambipolar diffusion becomes increasingly important in the FUV-ionized layer, but this is not the case as shown by more sophisticated photochemistry calculations (Walsh et al. 2012). The resulting ionization fraction is shown in Figure 1 for various radii.

<figure><img src="image_4.png" /><figcaption>Figure 2. Initial ohmic (Λ), Hall (Ha), and ambipolar (Am) Elsässer numbers as a function of height, when radius $R = 5$ au and $\beta_0 = 10^4$. A smaller Elsässer number indicates a stronger nonideal MHD effect. Here, the Hall effect is the strongest nonideal effect within $|z| < 2H$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Initial ohmic, Hall, and ambipolar Elsässer numbers across radii 1–30 au. The range of midplane values from $|z| < 2H$ are shown; the vertical profiles at each radius are qualitatively similar to those shown in Figure 2. The Hall effect is the strongest nonideal MHD effect in our disk model for $|z| < 2H$ and at all radii considered.</figcaption></figure>

We quantify the diffusivity coefficients $\eta_O$, $\eta_H$, and $\eta_A$ (Equations (5)–(7)) via the dimensionless Elsässer numbers

$$ \Lambda = \frac{v_A^2}{\Omega_0 \eta_O}, \quad (17) $$

$$ \text{Ha} = \frac{v_A^2}{\Omega_0 \eta_H}, \quad (18) $$

$$ \text{Am} = \frac{v_A^2}{\Omega_0 \eta_A}, \quad (19) $$

where $v_A$ is the Alfvén speed

$$ v_A = \frac{|B|}{\sqrt{\rho}}, \quad (20) $$

The initial vertical profiles of $\Lambda$, Ha, and Am at 5 au are shown in Figure 2, and the initial values in the midplane region $|z| < 2H$ are shown for radii between 1 and 30 au in Figure 3.

The disk ionization prescription adopted by this work is generous. Because the cosmic-ray flux through protoplanetary
<table><caption>**Table 2**
Shearing Box Simulations</caption><thead><tr><th>Label<br>(1)</th><th>R<sub>au</sub><br>(2)</th><th>β<sub>0</sub><br>(3)</th><th>log<sub>10</sub>ζ<sub>CR,0</sub><br>(4)</th><th>Am<sub>0</sub><br>(5)</th><th>Ha<sub>0</sub><br>(6)</th><th>Λ<sub>0</sub><br>(7)</th><th>α&#773;<br>(8)</th><th>α&#773;<sub>|z|&lt;H</sub><br>(9)</th><th>δv/c<sub>s</sub><sup>2</sup><sub>|z|&lt;H</sub><br>(10)</th></tr></thead><tbody><tr><td>R1–b4p</td><td>1</td><td>+10<sup>4</sup></td><td>−16</td><td>0.27</td><td>0.0065</td><td>0.023</td><td>2.26 (−3)</td><td>7.86 (−3)</td><td>7.50 (−2)</td></tr><tr><td>R5–b3p</td><td>5</td><td>+10<sup>3</sup></td><td>−16</td><td>1.7</td><td>0.54</td><td>37</td><td>4.91 (−3)</td><td>6.43 (−3)</td><td>5.79 (−2)</td></tr><tr><td>R5–b4p</td><td>5</td><td>+10<sup>4</sup></td><td>−16</td><td>1.7</td><td>0.17</td><td>3.7</td><td>1.16 (−3)</td><td>1.91 (−3)</td><td>6.80 (−2)</td></tr><tr><td>R5–b4n</td><td>5</td><td>−10<sup>4</sup></td><td>−16</td><td>1.7</td><td>0.17</td><td>3.7</td><td>7.95 (−4)</td><td>5.21 (−4)</td><td>5.21 (−2)</td></tr><tr><td>R5–b4p–weak-cr</td><td>5</td><td>+10<sup>4</sup></td><td>−17</td><td>0.55</td><td>0.054</td><td>1.2</td><td>1.19 (−3)</td><td>1.57 (−3)</td><td>5.15 (−2)</td></tr><tr><td>R5–b5p</td><td>5</td><td>+10<sup>5</sup></td><td>−16</td><td>1.7</td><td>0.054</td><td>0.37</td><td>3.13 (−4)</td><td>6.09 (−4)</td><td>3.27 (−2)</td></tr><tr><td>R10–b4p</td><td>10</td><td>+10<sup>4</sup></td><td>−16</td><td>2.4</td><td>0.44</td><td>22</td><td>1.34 (−3)</td><td>2.06 (−3)</td><td>5.80 (−2)</td></tr><tr><td>R10–b4n</td><td>10</td><td>−10<sup>4</sup></td><td>−16</td><td>2.4</td><td>0.44</td><td>22</td><td>1.06 (−3)</td><td>9.22 (−4)</td><td>6.28 (−2)</td></tr><tr><td>R30–b4p</td><td>30</td><td>+10<sup>4</sup></td><td>−16</td><td>3.2</td><td>1.7</td><td>320</td><td>1.82 (−3)</td><td>2.37 (−3)</td><td>3.81 (−2)</td></tr><tr><td>R30–b4n</td><td>30</td><td>−10<sup>4</sup></td><td>−16</td><td>3.2</td><td>1.7</td><td>320</td><td>1.69 (−3)</td><td>2.17 (−3)</td><td>8.84 (−2)</td></tr></tbody></table>

**Note.** Table of parameters and key diagnostics for the shearing box simulations: Column (1) Simulation label. Column (2) Radial location of the shearing box in astronomical units. Column (3) Initial vertical magnetic field strength, see Equation (21). Column (4) Cosmic-ray ionization rate. Columns (5)–(7) Initial midplane values of the ambipolar, Hall, and ohmic Elsässer numbers, respectively. Column (8) Time-averaged $\alpha$ parameter, see Equation (22). Column (9) Same as column (8), but the volume average is computed within $z = \pm H$. Column (10) Time- and volume-averaged turbulent velocity as a fraction of the sound speed, computed within $z = \pm H$. Values written as $a$ (b) denote $a \times 10^b$.

disks is highly uncertain (disks may be shielded from cosmic-ray flux by, e.g., stellar winds Cleeves et al. 2013, 2015), or attenuated along sheared magnetic field lines (Fujii & Kimura 2022) we also examine a setup with cosmic-ray flux weaker by a factor of 10 (R5–b4p–weak-cr, see Table 2) as in Umebayashi & Nakano (1981).

## 2.4. Set of Simulations

We assume the shearing box domain is threaded by a vertical magnetic field $B = B_0\hat{e}_z$ whose initial strength is specified by the signed ratio between the gas and magnetic pressures

$$ \beta_0 \equiv \frac{B_0}{|B_0|} \frac{2\rho_0 c_{s,0}^2}{B_0^2}, \qquad (21) $$

where $\rho_0$ is the initial midplane density of the disk (and we define $P_0 = \rho_0 c_s^2$ to be the initial midplane pressure). A positive $\beta_0$ (+) means the initial vertical magnetic field is aligned with the disk angular momentum vector ($\Omega \cdot B_0 > 0$), and a negative $\beta_0$ (−) refers to the anti-aligned case ($\Omega \cdot B_0 < 0$).

We explore radii of 1, 5, 10, and 30 au. We take a fiducial $\beta_0 = +10^4$ (which corresponds to $|B_{z,0}| = 0.1, 0.01, 0.005$, and 0.001 G at 1, 5, 10, and 30 au, respectively) and explore a strong ($\beta_0 = +10^3$, $|B_{z,0}| = 0.04$ G) and weak ($\beta_0 = +10^5$, $|B_{z,0}| = 0.004$ G) magnetic field at 5 au. We also explore anti-aligned field orientations where $\beta_0 < 0$ at 5, 10, and 30 au, for a field strength $|\beta_0| = 10^4$. We set $c_s = \Omega_0 = \rho_0 = 1$ in dimensionless code units.

The dimensions for each simulation are $(L_x, L_y, L_z) = (4H, 8H, 12H)$ with a resolution of 32 grid zones per scale height $H$. This resolution is comparable to but slightly higher than previous 3D shearing box simulations that include the Hall effect alongside ohmic and ambipolar diffusion (e.g., Lesur et al. 2014; Bai 2015; Simon et al. 2015b, 2018). The full set of simulations can be found in Table 2.

## 3. Results

We have run local simulations centered at several radii between 1 and 10 au in our model disk; in this region, the Hall effect is the dominant nonideal process at the midplane

(Kunz & Balbus 2004; Wardle 2007, or see Figure 2 in this work). We also run two simulations at 30 au where the Hall effect is less dominant but still present (Simon et al. 2015b; Bai 2015).

### 3.1. Angular Momentum Transport

We examine the $R\phi$ component of the stress tensor, which is responsible for radial angular momentum transport. We quantify this in terms of the Shakura & Sunyaev (1973) $\alpha$ parameter.

$$ \alpha = \frac{\langle \rho v_x v_y - B_x B_y \rangle}{\langle \rho c_s^2 \rangle}, \qquad (22) $$

where angled brackets denote the volume average

$$ \langle Q \rangle = \frac{1}{L_x L_y L_z} \int Q(x, y, z) \, dx dy dz. \qquad (23) $$

When $|\beta_0| = 10^4$, we find $\bar{\alpha} \sim 10^{-3}$ (Figure 4; Table 2), where the overline denotes the time average

$$ \bar{Q} = \frac{1}{\tau} \int Q(t) \, dt. \qquad (24) $$

We average over times when $\alpha$ is in a statistically steady state, which corresponds to the last $\sim 140$ orbits for all but the run at 1 au; in that run, we average over the last $\sim 35$ orbits. Our simulations with a stronger (weaker) initial magnetic field exhibit distinctly larger (smaller) $\alpha$. The $\alpha$ stresses exhibit small amplitude fluctuations over time, typically by a factor of $\lesssim 2$, which may indicate some level of turbulent angular momentum transport in our simulations.

We compute the approximate corresponding mass accretion rate ignoring any contribution from disk winds and assuming a steady-state disk:

$$ \dot{M} = 3\pi\alpha c_s H \Sigma. \qquad (25) $$

At 5 au, our disk model suggests that an $\alpha = 10^{-3}$ drives a mass accretion rate of $\dot{M} \approx 5 \times 10^{-9} M_\odot \text{yr}^{-1}$. This value should be treated as a lower limit because we ignore the contribution from magnetic winds.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Shakura & Sunyaev (1973) α parameter over time for different field strengths and orientations (aligned with the disk rotation, +, and anti-aligned, —) at 5 au. The green, black, blue, and orange curves correspond to β₀ = +10³, +10⁴, -10⁴, and +10⁵, respectively. In all cases, the α stress exhibits small amplitude fluctuations over time. Furthermore, α decreases with increasing |β₀|, and is also lower for the anti-aligned run compared with the aligned field run of the same field strength throughout most of the evolution (though they appear to converge at later times, at least in the |β₀| = 10⁴ runs).</figcaption></figure>

## 3.2. Magnetic Field Structure and Evolution

We present the space-time diagram of the toroidal (azimuthal) magnetic field $\langle B_y \rangle_{xy}$ and the corresponding Maxwell stress $\langle -B_x B_y \rangle_{xy}$, where $\langle \cdot \rangle_{xy}$ denotes the horizontal average;

$$ \langle Q \rangle_{xy} = \frac{1}{L_x L_y} \int Q(x, y, z) dxdy \qquad (26) $$

for the runs with R = 5 au and β₀ = 10⁴ (Figure 5) and R = 5 au and β₀ = -10⁴ (Figure 6). In both the aligned and anti-aligned cases there is strong, large-scale Maxwell stress close to the disk surface. In the midplane of the disk, the Maxwell stress has small-scale vertical fluctuations in intensity. Additionally, the midplane Maxwell stress is extremely sensitive to the magnetic field polarity. In the aligned case, the Maxwell stress remains strong in both the midplane and surface regions of the disk, but in the anti-aligned case the midplane Maxwell stress (|z| ≲ 2H) is suppressed by a factor ~4.

The qualitative behavior of the toroidal field $\langle B_y \rangle_{xy}$ in our simulations also changes significantly between aligned and anti-aligned B_z,0. Both simulations exhibit thin sheets of current generated between regions of magnetic field with opposing signs, but here the similarities end. In the aligned case, multiple current sheets are present at the beginning of the simulation, though all are transient and last for fewer than ~50 orbits before the toroidal field organizes into an “even” configuration $[B_y(z) \sim B_y(-z)]$. In the anti-aligned case, there is typically only one current sheet present at any time, and it does not remain at the midplane at all times. The toroidal field is organized into an “odd” configuration,<fnref n="6" /> with the toroidal field changing sign across the current sheet. In our other simulations, we generally find multiple coexisting current sheets (R5–b3p, R5–b4p-weak-cr), and current sheets which are off-midplane and not

necessarily at a constant height over time (R5–b3p, R5–b4p-weak-cr, R5–b4n, all simulations with R_AU ≥ 10).

Our results are consistent with previous studies in that complex structures in $\langle B_y \rangle_{xy}$ are present regardless of field alignment or the inclusion of the Hall effect. For instance, Bai & Stone (2013a, 2013b), Bai (2013), and Lesur et al. (2014) found off-midplane current sheets. Bai (2015) often found off-midplane or multiple coexisting current sheets, and patterns where the sign of $\langle B_y \rangle_{xy}$ flips periodically for various radii and β₀. Lesur et al. (2014) found that all of their simulations that include the Hall effect eventually eject the current sheet and exhibit even symmetry and speculate that all of their simulations are eventually symmetric on long enough time-scales (≈10³ Ω⁻¹). However, our work and Bai (2015) find that the system starting from an even configuration can occasionally spontaneously generate a current sheet and reverse the sign of B_y (see Section 3.7); thus, current sheets could be present at later times in the disk evolution by spontaneous generation. Global disks do not possess the radial symmetry of the shearing box (i.e., the radial direction of the central star is known) and require an odd number of current sheets in or slightly vertically offset from the disk (see the simulations of, e.g., Bai 2017; Bethune et al. 2017; Hu et al. 2023) such that the toroidal magnetic field changes sign and the vertical gas outflow and magnetic field lines that leave the disk are bent radially away from the star both above and below the midplane.

The presence of a current sheet could have implications for the level of turbulence in our simulations. Although turbulence may be generated by magnetic reconnection via, e.g., the tearing instability (Galeev & Zelenyi 1976), there are often only 2–4 grid cells across the width of a given current sheet for the simulations in our work. Because numerical simulations possess inherent numerical diffusion in addition to physical diffusion (e.g., ohmic and ambipolar diffusion), it is not clear how strong a role physical diffusion processes have in affecting the current sheets in our simulations. With current sheet widths being close to the grid scale it seems likely that they are predominately regulated by numerical diffusion and that any turbulence produced is the result of numerical reconnection and not small-scale instabilities. In any case, we will evaluate the relationship between the turbulent velocity and the current sheets in Section 3.7, while keeping in mind its potential numerical origin.

We examine the toroidal and poloidal magnetic field, in the last snapshot of the fiducial simulation (R5–b4p; 5 au, β = +10⁴). Every other simulation is qualitatively similar. The magnetic field is dominated by the toroidal component (Figure 7). In a laminar velocity flow, fluid layers move smoothly past each other with little to no mixing; here, the magnetic field lines are organized and nearly parallel. Therefore, we characterize the magnetic field as laminar despite the field lines not representing a flow themselves. The poloidal field (Figure 8) has a single concentration of magnetic flux located at x ~ -0.5H that extends vertically for ~3–4H around the disk midplane. Magnetic organization of the vertical magnetic field is also found by Riols & Lesur (2019) with only ohmic and ambipolar diffusion. The poloidal field strength in this feature is comparable to that of the toroidal field, and appears alongside a weaker but more turbulent field with substantially tangled field lines; the convoluted poloidal magnetic field background may be an indicator of a turbulent flow despite a strong, laminar toroidal field.

<footnote n="6">⁶ This “odd” symmetry is only approximate considering how far the current sheet deviates from the midplane at times. Nonetheless, we maintain this nomenclature to be consistent with the literature.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Space-time diagram of $\langle B_y/\sqrt{P_0}\rangle_{xy}$ (top) and $\log_{10}\langle -B_x B_y/P_0\rangle_{xy}$ (bottom) for the run with $R = 5$ au and $\beta_0 = 10^4$. The toroidal field is characterized by an even symmetry ($B_y(z) \simeq B_x(-z)$) after 50 orbits. The Maxwell stress is large throughout the vertical extent of the disk. The disk surfaces ($|z| \gtrsim 3H$) exhibit a large-scale, mostly laminar Maxwell stress; this is distinct from the midplane region of the disk, which exhibits Maxwell stress with a small-scale structure.</figcaption></figure>

### 3.3. Velocity Structure

Having made contact with the literature via an analysis of the magnetic field structure and evolution, we now turn to the question of the gas motions themselves. As mentioned in Section 2.1, the velocity $v$ is the gas velocity relative to the background shear $v_0 = -q\Omega x \hat{v}_y$.

We follow Simon et al. (2015a) in extracting the velocity fluctuations from the total velocity field at each snapshot. To remove bulk horizontal flows, we subtract the horizontally averaged velocities at each height:

$$v' = v - \langle v \rangle_{xy}. \qquad (27)$$

We also remove any zonal flows—large-scale, axisymmetric radial variations in the azimuthal velocity—by subtracting the azimuthal average of $v'_y$:

$$v_y'' = v_y' - \langle v_y \rangle_y. \qquad (28)$$

The velocity fluctuation is then

$$\delta v = (\delta v_x, \delta v_y, \delta v_z) = (v_x', v_y'', v_z'). \qquad (29)$$

Although this method produces velocity fluctuations that appear disordered, it does not guarantee that the flow is turbulent; nonetheless, we refer to all three components of the velocity fluctuations as turbulent hereafter for simplicity. We examine the velocities in the last snapshot for the simulations at 5 au, $\beta_0 = \pm 10^4$. The y-component of the midplane velocities is shown in Figure 9. The total velocity field $v_y(x)$ (though still without the Keplerian shear flow) is dominated by a large-scale, zonal flow-like structure at a location that approximately matches that of the magnetic concentration in the poloidal field (Figure 8). The strength of this flow is sensitive to the polarity of the magnetic field, being stronger by a factor of $\sim 4$ in the aligned case (+). However, the strength of the small-scale fluctuations $\delta v(x)$ is not sensitive to the magnetic field polarity.

The flows that exist in the midplane (Figure 10) appear to be somewhat organized, with only a few apparent eddies. However, the velocity fluctuations in the vertical plane (Figure 11) do appear to be turbulent, with many eddies and stochastic flow. The level of the turbulence has a strong dependence on height in part due to the highly variable disk

ionization (as in, e.g., Bai 2015; Simon et al. 2018) and the strength of $\beta$ (e.g., see Figure 2 of Yang et al. 2018). Figure 12 shows the magnitude and components of the time- and horizontally rms-averaged turbulent velocities. Near the disk surface, the turbulent velocity magnitudes approach the sound speed $c_s$ in agreement with a number of other works that explored a large range of parameters and numerical setups (e.g., Fromang & Nelson 2009; Simon et al. 2011a, 2011b, 2015a; Flock et al. 2012; Yang et al. 2018); evidently, this increase in turbulent velocity up to approximately the sound speed is a robust feature of magnetized disks for a number of different parameters. We further quantify the midplane velocities by averaging $\langle \delta v^2 \rangle_{xy}^{1/2}$ over the innermost two scale heights, i.e., for $|z| < H$ (Table 2). At 5 au, the midplane turbulence is $\approx 0.068 c_s$ when $\beta_0 > 0$ and $\approx 0.052 c_s$ when $\beta_0 < 0$. In general, the strength of the midplane turbulence in our simulations differs by only a factor <2 between the corresponding $\beta_0 > 0$ and $\beta_0 < 0$ configurations. The magnitude of the components of $\langle \delta v \rangle_{xy}$ are each within a factor $\sim 3$ of the others. The vertical component of turbulent velocities was briefly examined by Bai (2015); those turbulent velocities near the midplane are $\sim 2-3$ times lower than that of the corresponding runs in our simulations, but Bai (2015) considered a cosmic-ray ionization rate of $10^{-17} \text{ s}^{-1}$, an order of magnitude lower than our fiducial value (see Section 3.8).

### 3.4. Turbulent Velocity Power Spectrum

Gas turbulence is often described as being composed of eddies having a range of spatial scales $l$ and wavenumber $k = 2\pi/l$. The energy spectrum is a measure of the energy contained at a wavenumber $k$. The Fourier transform for some flow quantity $Q$ is defined as

$$\tilde{Q}(k) = \int Q(x) e^{-ik \cdot x} dx^3. \qquad (30)$$

We examine the turbulent kinetic energy in Fourier space, defined as

$$E_{3D}(k) = \frac{1}{2} \sum_{i=x,y,z} |\sqrt{\rho_0} \delta v_i(k)|^2. \qquad (31)$$
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Space-time diagram of $(B_y/\sqrt{P_0})_{xy}$ (top) and $\log_{10}(-B_x B_y/P_0)_{xy}$ (bottom) for the run with $R = 5$ au and $\beta_0 = -10^4$. The toroidal field changes sign across an approximately horizontal (xy-plane) current sheet. Similar to the ($\beta_0 > 0$) case, the Maxwell stress is large near the disk surface; however, in the midplane the Maxwell stress is largely suppressed.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Slice in the xy-plane of the toroidal field $B_y(z=0)$ in the last snapshot of the 5 au, $\beta_0 = +10^4$ run. Streamlines are overlaid to show magnetic field lines. The field lines of $B_y$ are laminar, in that streamlines move in organized, parallel paths with minimal deviation.</figcaption></figure>

We then integrate over spherical shells of radius $k = |\mathbf{k}| = \sqrt{k_x^2 + k_y^2 + k_z^2}$ and width $dk$,

$$ 
\begin{aligned}
E(k)dk &= E_{3D}(\mathbf{k})dk^3 \\
E(k) &= \int_{k-dk/2}^{k+dk/2} E_{3D}(\mathbf{k}')4\pi k'^2 dk'
\end{aligned} 
\quad (32) 
$$

and normalize by the total energy $E_{\text{tot}} = \int E(k) dk$. Here, $E_{3D}(\mathbf{k})$ denotes the three-dimensional spectrum obtained via the

Fourier transform, and $E(k)$ is the one-dimensional energy spectrum obtained by integrating $E_{3D}(\mathbf{k})$. The resulting energy spectra are shown in Figure 13.

The spectra broadly resemble that of an energy cascade from small wavenumbers, where energy is introduced into the system at large scales, to large wavenumbers, where energy is dissipated by viscous and magnetic diffusion. However, the energy spectra noticeably lack a distinct inertial range at intermediate wavenumber (similar to, e.g., Fromang & Papaloizou 2007; Simon et al. 2009), which should follow a simple power law $E(k) \propto k^p$.

The precise form of the energy spectrum has important implications for planet formation, as turbulence is the dominant force driving collisions between grains (Ormel & Cuzzi 2007). The collision velocities between small grains are extremely sensitive to $p$, and can vary by up to an order of magnitude for an inertial range scaling $-5/3 < p < -3/2$ (Gong et al. 2021). Many previous numerical MHD simulations have reported either $p \sim -4/3$ for driven MHD turbulence (Lemaster & Stone 2009; Salvesen et al. 2014; Makwana et al. 2015; Grete et al. 2017; Gong et al. 2020) or $p \sim -1.5$ for turbulence driven by the MRI (Fromang 2010; Lesur & Longaretti 2011; Walker et al. 2016). Only a small number of the works referenced above that report $p = -4/3$ consider nonideal MHD (these being Lesur & Longaretti 2011; Salvesen et al. 2014; Walker et al. 2016) and these works only consider ohmic diffusion. Additionally, many of the above works also contain spectra with a small or nonexistent<fnref n="7" /> inertial range akin to the spectra in Figure 13.

Because the magnitude of turbulence has a strong dependence on the vertical distance from the midplane—low-Mach number turbulence located in the midplane region ($|z| \lesssim 3H$), and the relatively large-Mach number turbulence located near the disk surface ($|z| \gtrsim 3H$)—we quantify the effect of this height dependence on the energy spectra. We computed spectra for the $+\beta_0$ case at different heights above (and below) the midplane, averaged over the last $40 \Omega^{-1}$ of the simulation (Figure 14). Here, $E(k)$ is obtained by integrating the 2D spectra (at some height $z$) over annuli of radius $k = \sqrt{k_x^2 + k_y^2}$ and width $dk$. Despite the large variation in the strength of turbulence across height, each

<footnote n="7">⁷ Our preference for a definitive inertial range would be a clear power law scaling that lasts well over one decade in $k$, which is not typically seen in these previous works (but see, e.g., Salvesen et al. (2014)).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Slice in the xz-plane of the poloidal magnetic field $B_p(y=0)$ in the last snapshot of the 5 au, $\beta_0 = +10^4$ run. Streamlines are overlaid to show magnetic field lines. The poloidal field is much weaker than the toroidal field, except in a highly concentrated region of magnetic flux at $x \sim -0.5H$, where the strengths are comparable. Elsewhere in the shearing box, the poloidal magnetic field is weak but twisted and nonlaminar.</figcaption></figure>

spectrum has mostly similar scaling and no clear inertial range, though at large vertical distances from the midplane ($z \gtrsim 4H$), the slope is slightly shallower at small wavenumber. Since there remains no clear inertial range at these separate heights, we conclude that the lack of an inertial range in our energy spectra is not the result of computing the spectra for a mixture of low-Mach and high-Mach turbulence and is thus more likely to be a result of insufficient resolution. It remains unclear what resolution would be necessary to see a clear inertial range, though Sengupta & Umurhan (2023) suggest that the dissipation scale begins at approximately one-quarter the Nyquist frequency (8 grid cells)<fnref n="8" />

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Azimuthal (y) component of the midplane velocities for the simulations at 5 au, with $\beta_0 = +10^4$ (top) and $\beta_0 = -10^4$ (bottom). The total (though still shear-subtracted) velocity field (left panels) is dominated by large-scale flows, the strength of which is sensitive to the polarity of the magnetic field. However, the turbulent velocity (right panels) is largely insensitive to the polarity. Some regions are purposely saturated in their color in order to emphasize the small-scale turbulent modes.</figcaption></figure>

### 3.5. Turbulent Velocity Autocorrelation Function

We also examine the spatial autocorrelation function (ACF; Guan et al. 2009; Simon et al. 2012, 2013b; Bai & Stone 2013b) to examine the structure in the xy-plane,

$$ \text{ACF}(Q(\Delta\mathbf{x})) = \left[ \frac{\int Q(t, \mathbf{x}) Q(t, \mathbf{x} + \Delta\mathbf{x}) d^3\mathbf{x}}{\int Q(t, \mathbf{x})^2 d^3\mathbf{x}} \right], \quad (33) $$

where $Q$ is the quantity of interest. For vector quantities, we follow Guan et al. (2009) and Simon et al. (2012) and define $\text{ACF}(Q) = \sum_{i=x,y,z} \text{ACF}(Q_i)$. We consider the ACF at a vertical offset $\Delta z=0$. The ACFs (Figure 15) show two distinct components, consistent with previous studies: a strong, tilted centroid region, and a volume-filling background component. The centroid region denotes a strong, localized correlation in the field, while the background component corresponds to the velocity flow (and magnetic field) structures on the largest scale of the numerical domain (Simon et al. 2012). The ACFs of the total velocity and magnetic field contain strong background components; this indicates strong, large-scale structures in these fields. As expected, the ACF of the turbulent velocity does not contain this background component ($|\text{ACF}(\delta v)| < 0.2$ for most of the domain).

<footnote n="8">⁸ Sengupta & Umurhan (2023) examined hydrodynamic turbulence using the Pencil code with the sixth-order accurate hyperviscosity scheme of Lyra et al. (2017). This causes the velocity field to dissipate energy at the smallest spatial scales, preserving the power spectra at larger scales. It is not clear at what spatial scale energy dissipation becomes important in our simulations, which are third-order accurate in space.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The xy-slice of the turbulent velocity $\delta v(z = 0)$ in the last snapshot of the 5 au, $\beta_0 = +10^4$ run. Streamlines of $\delta v$ (Equation (29)) show deviations from the mean flow. The flows in the xy-plane do not appear to be strongly turbulent (i.e., the velocity structure appears somewhat ordered), possibly because fluctuations in this plane are stretched by the background shear.</figcaption></figure>

Because the autocorrelation traces out coherent, filamentary structure (one expects a large, positive correlation along such filaments), the centroid tilt angle from the azimuthal axis is related to the Maxwell and Reynolds stresses. In the case of the magnetic field,

$$ \frac{\langle B_x B_y \rangle}{\langle B^2 \rangle} = \sin 2\theta_B \\ \theta_{B, \text{tilt}} \approx \theta_B. \qquad (34) $$

Studies of MRI-driven turbulence in ideal MHD (Guan et al. 2009; Beckwith et al. 2011; Simon et al. 2012; Bai & Stone 2013b) typically find tilt angles of ~7° for the velocity autocorrelation and ~15° for the magnetic field autocorrelation. In Figure 15, $\theta_{v, \text{tilt}} \sim 6°$ for ACF($\nu$) and ACF($\delta\nu$), in good agreement with these previous works. However, we find $\theta_{B, \text{tilt}} \sim 3°$, which is significantly smaller than that in ideal MHD simulations of the MRI. The small tilt angle may be reflective of the degree to which the laminar toroidal field dominates over the other magnetic field components. For example, in the ideal MRI simulations of Simon et al. (2009) $\langle B_y^2 \rangle / \langle B_x^2 \rangle \sim 10$, while in our simulations $\langle B_y^2 \rangle / \langle B_x^2 \rangle$ is more typically ~100.

While the centroid is well contained by the simulation domain in the turbulent velocity autocorrelation, we note the presence of a tail along the centroid major axis that wraps around through the periodic boundary in azimuth into the other side of the domain. This implies that there may be flows that are correlated in azimuth on scales larger than the numerical

<figure><img src="image_13.png" /><figcaption>**Figure 11.** The xz-slice of the turbulent velocity $\delta v(y = 0)$ in the last snapshot of the 5 au, $\beta_0 = +10^4$ run. Streamlines of $\delta v$ (Equation (29)) are also shown. The gas appears to be turbulent in the xz-plane, with a large number of eddies and a disorganized flow.</figcaption></figure>

domain. Referring again to Figure 9, the turbulence consists of small length scales in the x-direction, but larger length scales (possibly larger than the azimuthal domain, particularly for the $\beta_0 > 0$ case) in the y-direction. While a larger azimuthal extent would better capture the tail of the centroid, the ACF ($\delta\nu$) along the major axis tail of the centroid is $\lesssim 0.3$ at the boundary of the domain; thus, the main centroid of strong correlation is still well contained within the simulation domain.

### 3.6. Vertical Diffusion

We follow Zhu et al. (2015) and Youdin & Lithwick (2007) to compute an estimate for the particle scale height $H_d$. We compute the vertical gas diffusion coefficient via

$$ D_{g,z} = \int_{0}^{\infty} \langle v_{g,z}(t) v_{g,z}(0) \rangle_{|z|<H} dt, \qquad (35) $$
<figure><img src="image_14.png" /><figcaption>Figure 12. Time-averaged and rms turbulent velocity profile normalized by the sound speed vs. vertical height. The turbulent velocity is large (~$c_s$) in the highly ionized upper layers of the disk. In the midplane, the turbulent velocities are weaker but still substantial, with $\delta v_{|z|<H} \sim 5-7 \times 10^{-2}$. The turbulent velocity components are consistently within a factor ~3 of each other.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Turbulent kinetic energy spectra at 5 au, for $\beta_0 = \pm 10^4$. Neither spectra exhibit a clear inertial range with a power law scaling $k^p$. The slopes are somewhat steeper than $p = -4/3$ or $p = -3/2$, which have been reported from previous studies of MHD turbulence (see the main text).</figcaption></figure>

where $v_{g,z}$ is the vertical velocity of the gas. We find $D_{g,z} \approx 1.0 \times 10^{-3}$, which is roughly consistent with previous numerical simulations with large resistivity (Gole et al. 2016).

<figure><img src="image_16.png" /><figcaption>Figure 14. Energy spectra for the $\beta_0 > 0$ case, computed at different heights above (solid) and below (dashed) the disk midplane, averaged for the last $40 \Omega^{-1}$ of the simulation. The spectra have similar scaling (and no inertial range) despite the large variation in turbulence across height.</figcaption></figure>

The particle scale height (neglecting particle mass, and thus momentum feedback), is (Dubrulle et al. 1995; Youdin & Lithwick 2007; Zhu et al. 2015)

$$H_d = \frac{H}{\sqrt{c_s^2 t_s D_{g,z}^{-1} + 1}}, \qquad (36)$$

where $t_s$ is the stopping time of the particle. The stopping time is the time over which a particle velocity is reduced by a factor $e$ by a constant drag force (e.g., gas headwind) and is proportional to the particle size in the Epstein drag regime (Weidenschilling 1977). Because the vertical velocity fluctuations at the disk midplane are similar for both the aligned and anti-aligned magnetic field cases, we only examine the vertical stirring for $\beta_0 > 0$; the results are shown in Figure 16.

Particles with dimensionless stopping time $t_s\Omega = 0.1$ can potentially be diffused to ~10% of the gas scale height $H$. Smaller particles with $t_s\Omega = 10^{-3}$ can potentially be diffused to ~70% $H$, a significant fraction of the gas scale height. Although these grains would not be well sedimented, a grain abundance of a few percent or more might still lead to strong particle concentration in the radial direction (Yang et al. 2018), though Lim et al. (2024) have shown that hydrodynamic turbulence reduces the efficacy of strong particle concentration overall. In addition, small grains have an outsized influence on the gas ionization fraction by providing an additional recombination channel for electrons (and ions, to a lesser extent; see also Section 4.3.1). Because the small grains might essentially permeate the gas, the behavior of the turbulence in the disk midplane may be qualitatively changed as the coupling between the magnetic field and the gas becomes weaker. However, as the current simulations do not contain particles, the values of $H_d$ in Figure 16 should be taken as upper limits.

### 3.7. Effects of Magnetic Field Strength and Current Sheet

We investigate the effect of different initial magnetic field strengths on the turbulent velocities at 5 au (Figure 17). For a weak initial magnetic field ($\beta_0 = +10^5$), $\bar{\alpha} \sim 3 \times 10^{-4}$. The toroidal field exhibits an initial transient current sheet, but afterward maintains an even symmetry for the duration of the simulation. The time-averaged turbulent velocity magnitude in the
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Time-averaged ACF at zero offset from the vertical ($\Delta z = 0$) of the total velocity field (left), turbulent velocity field (center), and magnetic field (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Particle scale height $H_d$ as a fraction of the gas scale height $H$ (from Equation (36)) vs. stopping time in units of $\Omega^{-1}$ given the vertical gas diffusion coefficient $D_{g,z} = 1.0 \times 10^{-3}$. Circles mark the particle scale height at each decade. Large particles are more settled, but smaller particles may be diffused to a significant fraction of the gas scale height.</figcaption></figure>

midplane region (which we define to be the innermost two scale heights, or $|z| < H$) is $\sim 3.3 \times 10^{-2} c_s$, smaller by a factor $\sim 2$ compared with the $\beta_0 = +10^4$ configuration. These time averages correspond to times when $\langle B_y \rangle_{xy}$ does not exhibit a current sheet.

For a strong initial magnetic field ($\beta_0 = +10^3$), $\bar{\alpha} \sim 5 \times 10^{-3}$, but the midplane turbulent velocity magnitudes are not significantly larger than the $\beta_0 = +10^4$ configuration. Looking more closely at the magnetic field geometry, we note the strong field case results in a toroidal field with even symmetry (no current sheet) after initial transients until $t \approx 80$ orbits, at which time two current sheets are spontaneously generated near the midplane (Figure 18). At least one current sheet lasts for the next $\sim 40$ orbits, after which it is ejected from the shearing box and the sign of $\langle B_y \rangle_{xy}$ changes. Similar spontaneous current sheet generation was seen by Bai (2015), although it is not understood why this occurs. In the range of times over which a current sheet exists, the turbulent velocities at the midplane average to $\sim 1.7 \times 10^{-1} c_s$, which is enhanced by a factor $\sim 2.5$ times greater than the $\beta_0 = 10^3$ configuration without a current sheet. This is a general behavior of our simulations; Figure 19 shows the differences in midplane turbulent velocities between times with and without a current

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Same as Figure 12, but varying the strength of the magnetic field (i.e., varying $\beta_0$). A stronger (weaker) magnetic field results in stronger (weaker) gas turbulence. In the strong field simulation ($\beta_0 = +10^3$) the presence of a current sheet is correlated with enhanced midplane gas turbulence. Reducing the cosmic-ray flux ($\zeta_{CR} = 10^{-17}$) does not significantly reduce the level of midplane gas turbulence.</figcaption></figure>

sheet for the simulations that allow for such a comparison (ignoring simulations with only initial, transient current sheets). We observe a general trend of stronger midplane turbulence being accompanied by a stronger magnetic field (smaller $|\beta_0|$). At times when the toroidal field exhibits a current sheet, the midplane turbulence strength is a factor of $\sim 1.5–3$ stronger than in times without a current sheet. One exception is the simulation with $R = 30$ au, $\beta_0 = -10^4$, where the midplane turbulence is $\sim 5\%–20\%$ larger at times without a current sheet.

Whether or not a current sheet manifests in a shearing box simulation appears to be stochastic (Bai & Stone 2013a; see their Section 4.4). Bai (2013) found that the current sheet requires some coupling of the magnetic field to the gas and cannot exist in locations of strong magnetic diffusion. Otherwise, the only magnetic field left is the initial, background net vertical field (which cannot be removed from the domain). Because numerical diffusion acts to decouple the magnetic field from the gas (albeit nonphysically), strong numerical diffusion may facilitate the ejection of current sheets over long enough timescales. Previous shearing box studies centered in the inner disk (Bai & Stone 2013a; Lesur et al. 2014; Bai 2015)
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Space-time diagrams of the toroidal field $\langle B_y \rangle_{xy}$ and Maxwell stress $\log_{10}\langle -B_x B_y \rangle_{xy}$ for the simulation at 5 au with $\beta_0 = +10^3$. In the toroidal field, current sheets are spontaneously generated at around $t = 80$ orbits, both of which are ejected from the simulation domain by 120 orbits. The times at which these current sheets exist are correlated with times of enhanced midplane turbulence.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Midplane turbulent velocities averaged over times with (filled) and without (unfilled) current sheets exhibited by $\langle B_y \rangle_{xy}$. Different shapes correspond to different initial configurations in $\beta_0$ and radius. In general, the midplane turbulence is enhanced by a factor $\sim 1.5-3$ when a current sheet is present.</figcaption></figure>

employed a lower numerical resolution than this work, and thus have more numerical diffusion. To test the effect of resolution on current sheet generation and longevity, we compared our standard resolution runs with those of a lower resolution (8 grid cells per $H$) and found that the lower resolution runs exhibit current sheets that are less complex and much shorter lived, always resulting in an even toroidal field geometry at later times. Whether or not shearing boxes can sustain a current sheet is unclear; Bai (2013) speculates that the shearing box framework may favor magnetic field geometry without a current sheet, and Lesur et al. (2014) found that the current sheets in their simulations were ultimately ejected on long timescales ($t > 1400 \Omega^{-1}$). However, we have also seen that current sheets can be generated spontaneously (also seen in Bai 2015); thus the strength of turbulence in real disk systems may from time to time be enhanced by current sheets.

## 3.8. Effects of Cosmic-ray Flux

As the cosmic-ray flux through the disk is highly uncertain (see Section 2.3 above), we also examine the effect of a decreased

cosmic-ray ionization rate ($\zeta_{CR,0} = 10^{-17}$; Umebayashi & Nakano 1981), as in, e.g., Bai & Stone (2011), Bai (2015), Simon et al. (2018; Figure 17). Magnetically driven turbulence should be weaker with a lower degree of coupling between the gas and magnetic field. Despite this, we find that the turbulence is fairly insensitive to the degree of cosmic-ray ionization—this comparison is done only for times when a current sheet is not present, to control for the current sheet’s influence on the turbulence. A possible explanation for this is that although cosmic rays are the dominant source of ionization at the midplane in our disk model, reducing the cosmic-ray flux by a factor of 10 only reduces the ionization fraction by a less drastic factor $\sim 3$ under the metal-and-grain-free dissociative recombination approximation of Gammie (1996) and Fromang et al. (2002). The effects of weaker coupling of the magnetic field to the gas motions may only become apparent with a lower cosmic-ray flux than considered here.

## 3.9. Effects of Radial Disk Location

We investigate the effects of radial disk location on the magnetic field and velocity structure at $\beta_0 = \pm 10^4$ (Figure 20). Between 1 and 30 au midplane values are between $\sim 0.04$ and $0.1 c_s$ regardless of the polarity of the initial magnetic field, and without a clear trend across the radii considered. At 1 au, the turbulent velocities exceed the sound speed at the disk surface.

Across this large range of radii, the magnetic diffusivities can change significantly (Figure 3). The Hall effect is the dominant nonideal term at the disk midplane throughout the entire range of 1–30 au considered here. However, ohmic diffusion is also important further inward. In particular, ohmic diffusion is stronger than ambipolar diffusion (and comparable to the Hall effect) at 1 au, comparable to ambipolar diffusion at the midplane at 5 au, and weaker than both ambipolar diffusion and the Hall effect for $R > 10$ au.

Ohmic diffusion is known to quench turbulence driven by the MRI (Gammie 1996; Fleming & Stone 2003). With such weak ohmic diffusion at 10 au and beyond, it is natural to ask whether we expect to see turbulence driven by the MRI in our simulations. Because the strength of turbulence is relatively unchanged from that at 1 au, where ohmic diffusion is strong (and in fact, where all three nonideal terms are strong), it seems unlikely that we see the onset of MRI turbulence in our
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Same as Figure 12, but varying the disk radial location. Solid lines and (+) denote +β₀, and dashed lines and (-) mark -β₀. Between 1 and 30 au, the gas turbulence appears relatively insensitive to radial location and initial magnetic field polarity.</figcaption></figure>

simulations. Although this is initially counterintuitive, we recall that ambipolar diffusion is also known to damp MRI-driven turbulence in the outer disk (Simon et al. 2013a, 2013b; Cui & Bai 2021), and ambipolar diffusion is strong at all radii considered in our disk model (Am₀ ≲ 3). Therefore, the similarity in the strength of turbulence between 1 and 30 au may be due to the consistently strong ambipolar diffusion.

The Hall effect is strong at all radii considered (Ha₀ ≲ 1.7). Because the Hall effect is not a diffusive effect (despite η_H commonly being cast as the “Hall diffusivity”), its influence on the gas turbulence is more difficult to characterize. In the following section, we will discuss how the Hall effect could be the physical origin of the turbulence.

# 4. Discussion

## 4.1. Origin of the Turbulence

The strength of the initial vertical magnetic field B_z,0 and the presence (or lack of) a current sheet are the two factors with the largest effect on the strength of the gas turbulence in our simulations (Figures 17 and 19). This indicates that the turbulence in our simulations is primarily generated by magnetic effects. Between 1 and 30 au (depending on the disk model), the Hall effect is the dominant nonideal magnetic effect (Figure 2); when the Hall effect is strong, as it is in our simulations, the orientation of the magnetic field is also important (Balbus & Terquem 2001; Kunz 2008; Lesur et al. 2014; Bai 2015; Simon et al. 2015b). The influence of the magnetic field polarity can be seen when comparing the + and -β₀ configurations above at 5 au (Section 3.2), as the toroidal field and Maxwell stress are qualitatively and quantitatively distinct between the two cases. This is evidence that the Hall effect has a strong presence in our simulations. This difference is also apparent at 10 au, but much less pronounced at 30 au, where the Hall effect is weaker; therefore, it may be difficult to obtain observational signals of Hall MHD at larger radii.

Previous works (e.g., Lesur et al. 2014) have shown that the HSI (Kunz 2008) produces strong, laminar toroidal magnetic fields (when β₀ > 0), which we also observe in the current set of simulations. In addition to this, our simulations exhibit large-scale, laminar azimuthal velocities (Figure 9), which are strong

only when β₀ > 0. This could indicate that the HSI is the origin of these large-scale flows (or that an anti-aligned magnetic field prevents the magnetic self-organization like that seen by Riols & Lesur (2019) for ohmic and ambipolar diffusion only). Despite this, there is a strong turbulent gas component that is insensitive to the initial magnetic field orientation. In addition to the large-scale, axisymmetric HSI, the presence of which is determined only by the sign of β₀, other “guises” of the HSI may be at work.

The nonaxisymmetric HSI (hereafter, naHSI; Kunz 2008; Simon et al. 2015b) occurs when whistler waves become physically decoupled from the epicyclic motion of radially perturbed bulk neutral fluid, and may be driven unstable if magnetic perturbations grow by Keplerian shear faster than they can be advected by Alfvén waves. The dispersion relation that described this behavior can be found in Simon et al. (2015b). Eventually stabilizing terms dominate over destabilizing influences; in other words, the naHSI is inherently a transient instability analogous to the nonaxisymmetric MRI (Balbus & Hawley 1992). Simon et al. (2015b) found that the naHSI was likely responsible for sporadic bursts of magnetic stress in their simulations, even while (Ω·B) < 0, which lends credence to the possibility that this version of the HSI may be present in our simulations regardless of magnetic field polarity.

In addition to the naHSI, a local version of the axisymmetric HSI may generally survive (on scales less than the domain size) even when β₀ < 0. Although differential Keplerian rotation Ω is the single largest source of shear in the disk, weaker shear exists between any adjacent fluid parcels with different vorticity. Similarly, the magnetic field will evolve beyond the initial vertical field to have an arbitrary orientation in three dimensions (as in Figure 8). Therefore, the presence of the HSI in a local region may not be determined only by the sign of Ω·B. Thus, the HSI may locally persist at arbitrary orientations in our simulations regardless of the background magnetic field polarity.

Beyond the possibility of the naHSI or a local version of the HSI, it is not clear if one or more secondary mechanisms are at work generating turbulence in our simulations. We will follow up on these possibilities and investigate the origin of the turbulence in a subsequent paper.

## 4.2. Impact on Particle Growth and Planetesimal Formation

Although strong ohmic diffusion like that found in the inner disk will suppress MHD turbulence driven by the MRI, the Hall effect may still generate turbulence, replacing the classical “dead zone” (as in Gammie 1996; Balbus & Terquem 2001; Kunz & Balbus 2004). In this case, even if disk accretion is primarily driven by magnetothermal winds, significant turbulence can persist near the disk midplane. The existence of this turbulence may have significant implications for the growth of dust grains in planet-forming disks. The midplane turbulence levels found in this work are well above the fragmentation limit of Güttler et al. (2010) for collisions between compact grains of equal sizes. More recently, Kothe (2016; see also a review by Simon et al. 2022) found that grain growth could not proceed between equally sized silicate or water ice aggregates for collisional speeds of 25–50 m s⁻¹ (3.8–7.5 × 10⁻² c_s at 5 au in our disk model). The strong velocities produced in our simulations may impede the collisional growth of grains; however, collisional grain growth might still occur at these speeds between grains or aggregates of different sizes or
porosity (e.g., mass transfer, sticking; see Güttler et al. 2010 for details).

The levels of turbulence found in this work may greatly reduce the efficacy of planetesimal formation by preventing the clumping of particles via, e.g., the streaming instability either by affecting the linear modes of the instability (Chen & Lin 2020; Umurhan et al. 2020) or preventing clumping in the nonlinear state (Gole et al. 2020; Lim et al. 2024). For instance, while Li & Youdin (2021) found the critical midplane particle-to-gas density threshold $\epsilon_{crit}$ for planetesimal formation via the streaming instability to be ~0.3–1 (depending on the particle size), Lim et al. (2024) found that driven turbulence imposes a stricter threshold of $\epsilon_{crit} \sim 1–3$ for turbulent velocities as low as 0.01 $c_s$—the turbulent velocities found in this work greatly exceed these levels. However, Yang et al. (2018) demonstrated that anisotropies in the turbulent flow can actually lead to enhanced radial concentration of particles despite a lack of sedimentation toward the midplane.

## 4.3. Uncertainties and Caveats

### 4.3.1. Dust Grains

Our setup uses a generous ionization prescription. While we might expect some disks to have (relatively) high levels of ionization, perhaps, e.g., during an FU Orionis-type stellar outburst, the precise ionization levels at the disk midplane are highly uncertain. While our results show that the strength of turbulence is relatively insensitive to the precise level of cosmic-ray flux, it is important to consider other ways in which the ionization can be impacted. Dust grains are the electron absorbers that have the single largest effect on ionization. Including a self-consistent treatment of dust grains or grain growth in our simulations is beyond the scope of this paper, but here we briefly consider the effect dust grains might have on our results.

Dust grains are a catalyst for recombination (Weisheit & Upham 1978) and “soak up” electrons and ions from the gas, reducing the gas conductivity by several orders of magnitude (Wardle & Ng 1999; Wardle 2007; Bai 2011a, 2011b) and by extension altering the strength of low-ionization MHD effects. In particular, the Hall effect will be amplified by reduced gas conductivity (Equation (6)): the drift velocity between electrons and ions (Hall current) will increase as the number density of neutrals, which induces a collisional drag on the ions, is increased (though this trend cannot continue indefinitely—as the number of neutrals increase, the gas will increasingly decouple from the magnetic field and approach a purely hydrodynamic system). Grain-free treatments (e.g., Gammie 1996; Fromang et al. 2002; this work) assume that dust has completely settled to the disk midplane and thus does not affect the dynamics or ionization of the bulk gas. However, in general, turbulence will diffuse dust grains away from the midplane (Fromang & Papaloizou 2006) and mix grains into the bulk gas.

Lesur et al. (2014) increased the diffusion coefficients by a factor of 1000 in two of their simulations to mimic the effects of micron-sized grains and found that the magnetic stress in the disk was still enhanced by the Hall effect, though they did not examine the gas velocities. Simon et al. (2015a) found that with 0.1 µm grains modeled by the chemical reaction network of Bai & Goodman (2009), Bai (2011a, 2014), turbulent velocities were nearly identical to the grain-free case (for $R = 10–100$ au).

We speculate that the grain abundance in that work was too low for the grains to noticeably affect the ionization fraction (Xu & Bai 2016), and as an extension, the gas turbulence.

### 4.3.2. Numerical Limitations

As with any study relying on local simulations, we must ask whether the observed behavior will appear in a global setup. While shearing boxes are unsuited for studies of large-scale structures (e.g., winds), turbulence is inherently a small-scale, local phenomenon. Likewise, global simulations that strive to capture both large and small-scale phenomena are often too computationally intensive. Early global simulations that included nonideal MHD suggested that between ~1 and 10 au, the disk is laminar (e.g., Gressel et al. 2015; Bethune et al. 2017) with accretion driven primarily by magnetically launched winds. However, these models employ much lower resolution (respectively, 24 and 16 grid cells per scale height at their most refined) and both use 2D (axisymmetric) domains (although Bethune et al. 2017 do consider a few 3D simulations). Fundamental differences between 2D and 3D turbulence and low resolution may have affected the results of these early global simulations. More recent 3D, nonideal MHD global simulations focusing on the outer disk (>30 au) have employed static mesh refinement (Bai 2017) to better resolve small-scale turbulence (if present) in the bulk of the disk. These simulations produced turbulent velocities that are in good agreement with shearing box simulations (Simon et al. 2018) centered in the outer disk and with comparable resolution. Thus, it is possible that the earlier global simulations show little turbulence due to the lower resolution (e.g., no static mesh refinement) not adequately resolving turbulent velocity fluctuations at small scales.

One limitation of our setup is the radial symmetry of the shearing box (see Bai & Stone 2013a; Lesur et al. 2014), which can result in vertical outflows that are radially bent in opposite directions at the top and bottom vertical boundary (see Figure 11) when there is no current sheet in the toroidal magnetic field. Real disks and global simulations do not have radial symmetry, and thus should always have a current sheet (i.e., magnetic fields and outflow are always bent away from the central star). Whether or not a current sheet develops in a shearing box is highly stochastic (Bai & Stone 2013a). In a global simulation, this radial symmetry is broken; thus, the toroidal magnetic field should always exhibit an overall sign change in the azimuthal field across the vertical direction, with poloidal field lines bent away from the central protostar (though this sign change will not necessarily be at the midplane either, e.g., Bai 2017; Hu et al. 2023). Because of the local box approximation of the simulations in this work, it is possible for current sheets to escape the simulations domain.

Lastly, our shearing box setup is isothermal. Such a setup cannot capture the turbulence generated by effects that rely on temperature gradients or gas cooling (e.g., Lyra & Umurhan 2019). However, our setup is particularly advantageous for identifying turbulence generated exclusively by nonideal MHD.

### 4.3.3. Consistency with Disk Observations

There is a shortage of observations of turbulence in the 1–10 au regions of protoplanetary disks with which to compare our simulations. The most direct constraints on turbulence
come from the turbulent broadening of molecular emission lines (e.g., Flaherty et al. 2017, 2018, 2020) but these studies are limited by spatial resolution to large disk radii $\gtrsim$30 au, and in most cases limited to the uppermost surface layer of the disk by the large optical depth of the molecular line being observed (e.g., CO 3–2 rotational line). While spatial resolutions of 5 au have been achieved for millimeter continuum emission (e.g., DSHARP; Andrews et al. 2018; Franceschi et al. 2023; Pizzati et al. 2023) to our knowledge, these observations have not been used to constrain turbulence directly via line-broadening. Continuum observations of dust ring width can be used to constrain turbulence (e.g., Rosotti et al. 2020; Jennings et al. 2022) but these attempts must still resolve structures on the order of the gas scale height, i.e., ~4 au at a radius of 40 au, which means that turbulence can only be constrained in this way for radii $\gtrsim$40 au. The MAPS observations resolved emission lines at ~20 au scales (Öberg et al. 2021; Sierra et al. 2021), thus still limiting constraints to the outer disk. Finally, Carr et al. (2004) used hot rovibrational H₂O emission to find strong turbulence near the disk surface, but at a small radius of 0.3 au. Indirect measurements of turbulence via modeling the spectral energy distributions of young stellar objects (e.g., Mulders & Dominik 2012) are thus far unable to distinguish between strong and weak turbulence, given degeneracies between the strength of turbulence and other parameters, such as the dust-to-gas ratio.

The simulations presented here demonstrate midplane turbulent velocities between ~3 and $9 \times 10^{-2}c_s$ across all tested values of $\beta_0$ and radii. Line-broadening observations of optically thick lines (thus constraining these measurements to well away from the disk midplane) place the upper limits of the turbulence in the outer disk of HD 163296 to be $\lesssim 0.05c_s$ (Flaherty et al. 2017; in this work, turbulent gas motions near the midplane were also quantifiable) and TW Hydra to be $\lesssim 0.08c_s$ (Flaherty et al. 2018). Flaherty et al. (2020) find strong turbulence (0.25–0.33$c_s$) around DM Tau, consistent with earlier work by Guilloteau et al. (2012), who found a range of turbulent velocities 0.12–0.69$c_s$. Furthermore, Flaherty et al. (2020) place upper limits of $<0.08c_s$ and $<0.12c_s$ around MWC 480 and V4046 Sgr, respectively. More recently, Paneque-Carreño et al. (2024) found turbulent velocities of $\approx 0.4–0.6c_s$ at relatively large heights above the midplane ($z/r = 0.2–0.3$, where $z(r)$ is the vertical (radial) coordinate of the disk) in IM Lup. The level of turbulence observed in the outer disk surface (where the gas is in the ideal MHD approximation) cannot be compared in a precise manner to the turbulence level found in this work (in the inner disk, where the Hall effect is dominant). However, the turbulence limits constrained by such observations (e.g., $\lesssim 0.03c_s$, $\lesssim 0.08c_s$) are much smaller than the turbulent velocities produced in our simulations near the disk surface, which approach sonic speed. Thus, if such turbulent velocities are present in physical disks, they are only present in the inner disk where the Hall effect is strong. Observations with spatial resolutions that can probe the innermost radial regions ($R \lesssim 5$ au) of disks are needed to verify that such turbulence is present.

## 5. Summary

We performed 3D local shearing box simulations representing the inner 1–30 au of a model protoplanetary disk including all three nonideal MHD effects (ohmic diffusion, the Hall effect, and ambipolar diffusion). These simulations have the

highest resolution per gas scale height among similar simulations to date, with 32 zones per *H* in every dimension. We explored how turbulence is affected by magnetic strength and polarity as well as radial location within the disk. Our key results are as follows:

1. In the presence of the Hall effect, weakly ionized disks can produce vigorous gas turbulence, with turbulent velocity fluctuations of at least $\sim 3–9 \times 10^{-2}c_s$ (16–75 m s⁻¹ in our disk model) at the midplane and $\sim c_s$ (0.4–1 km s⁻¹) near the disk surface.

2. This turbulence is sensitive to the initial magnetic field strength and the presence or absence of a current sheet associated with a flip in the horizontally averaged toroidal field. As such, the turbulence is largely magnetically driven in nature.

3. Despite the role of magnetic fields in generating turbulence, the strength of the turbulence is not sensitive to the polarity of the initial vertical field.

4. The magnetic field is dominated by a strong, laminar toroidal component, with a weaker, turbulent poloidal component. Often, but not always, a current sheet is present near the midplane.

5. The midplane Maxwell stress is extremely sensitive to the polarity of the initial vertical magnetic field.

Overall, these findings show that gas in the inner (e.g., 1–30 au) regions of protoplanetary disks can be turbulent. As mentioned previously, there are several caveats to bear in mind. The neglect of dust grains and the physically motivated but generous ionization prescription present in this work may not reflect the conditions present in real protoplanetary disks. In addition, the toroidal field current sheet, which we have shown can cause drastic variability in the levels of turbulence, may not be adequately modeled without both a global domain and an even higher resolution. While these are clear uncertainties to address in future works, our results demonstrate the importance of considering the detailed gas dynamics in modeling such disks even if the accretion flow is mostly laminar and dominated by magnetically driven winds.

## Acknowledgments

We thank Kevin Flaherty, Meredith Hughes, and Xiao Hu for useful discussions related to this work. We thank the anonymous referee whose comments improved the quality of this paper. We acknowledge support from NASA via the Theoretical and Computational Astrophysics Network (TCAN) program through grant 80NSSC21K0497. D.G.R. and J.B.S. acknowledge support from NASA under the Emerging Worlds program through grant 80NSSC20K0702. D.C. acknowledges support from NASA under the Emerging Worlds program through grant 80NSSC21K0037. G.L. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant Agreement No. 815559 (MHDiscs)). W.L. acknowledges support from the NASA Emerging Worlds program via grant 80NSSC22K1419 and by NSF via grant AST-2007422. C.C.Y. is grateful for the support from NASA via the Astrophysics Theory Program (grant #80NSSC21K0141 and #80NSSC24K0133) and the Emerging Worlds program (grant #80NSSC20K0347 and #80NSSC23K0653). The computations were performed using Stampede2 at the Texas Advanced Computing Center using XSEDE/ACCESS grant TG-AST120062.
*Software:* Athena (Stone et al. 2008), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007).

ORCID iDs

David G. Rea https://orcid.org/0000-0002-5000-2747

Jacob B. Simon https://orcid.org/0000-0002-3771-8054

Daniel Carrera https://orcid.org/0000-0001-6259-3575

Geoffroy Lesur https://orcid.org/0000-0002-8896-9435

Wladimir Lyra https://orcid.org/0000-0002-3768-7542

Debanjan Sengupta https://orcid.org/0000-0003-0801-3159

Chao-Chin Yang (楊朝欽) https://orcid.org/0000-0003-2589-5034

Andrew N. Youdin https://orcid.org/0000-0002-3644-8726

References

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Bai, X.-N. 2011a, ApJ, 739, 50

Bai, X.-N. 2011b, ApJ, 739, 51

Bai, X.-N. 2013, ApJ, 772, 96

Bai, X.-N. 2014, ApJ, 791, 137

Bai, X.-N. 2015, ApJ, 798, 84

Bai, X.-N. 2017, ApJ, 845, 75

Bai, X.-N., & Goodman, J. 2009, ApJ, 701, 737

Bai, X.-N., & Stone, J. M. 2011, ApJ, 736, 144

Bai, X.-N., & Stone, J. M. 2013a, ApJ, 769, 76

Bai, X.-N., & Stone, J. M. 2013b, ApJ, 767, 30

Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152

Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610

Balbus, S. A., & Hawley, J. F. 1998, RvMP, 70, 1

Balbus, S. A., & Terquem, C. 2001, ApJ, 552, 235

Beckwith, K., Armitage, P. J., & Simon, J. B. 2011, MNRAS, 416, 361

Bethune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Carr, J. S., Tokunaga, A. T., & Najita, J. 2004, ApJ, 603, 213

Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132

Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204

Colella, P. 1990, JCoPh, 87, 171

Colella, P., & Woodward, P. R. 1984, JCoPh, 54, 174

Cui, C., & Bai, X.-N. 2021, MNRAS, 507, 1106

Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432

Desch, S. J. 2004, ApJ, 608, 509

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton, NJ: Princeton Univ. Press)

Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485

Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icar, 114, 237

Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659

Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109

Flaherty, K. M., Hughes, A. M., Teague, R., et al. 2018, ApJ, 856, 117

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908

Flock, M., Henning, T., & Klahr, H. 2012, ApJ, 761, 95

Franceschi, R., Birnstiel, T., Henning, T., & Sharma, A. 2023, A&A, 671, A125

Fromang, S. 2010, A&A, 514, L5

Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71

Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597

Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751

Fromang, S., & Papaloizou, J. 2007, A&A, 476, 1113

Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18

Fujii, Y. I., & Kimura, S. S. 2022, ApJL, 937, L37

Galeev, A. A., & Zelenyi, L. M. 1976, ZhETF, 70, 2133

Gammie, C. F. 1996, ApJ, 457, 355

Gole, D., Simon, J. B., Lubow, S. H., & Armitage, P. J. 2016, ApJ, 826, 18

Gole, D. A., Simon, J. B., Li, R., Youdin, A. N., & Armitage, P. J. 2020, ApJ, 904, 132

Gong, M., Ivlev, A. V., Akimkin, V., & Caselli, P. 2021, ApJ, 917, 82

Gong, M., Ivlev, A. V., Zhao, B., & Caselli, P. 2020, ApJ, 891, 172

Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84

Grete, P., O'Shea, B. W., Beckwith, K., Schmidt, W., & Christlieb, A. 2017, PhPl, 24, 092311

Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010

Guilloteau, S., Dutrey, A., Wakelam, V., et al. 2012, A&A, 548, A70

Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hartlep, T., & Cuzzi, J. N. 2020, ApJ, 892, 120

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742

Hu, X., Li, Z.-Y., Wang, L., Zhu, Z., & Bae, J. 2023, MNRAS, 523, 4883

Hunter, J. D. 2007, CSE, 9, 90

Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848

Jennings, J., Booth, R. A., Tazzari, M., Clarke, C. J., & Rosotti, G. P. 2022, MNRAS, 509, 2780

Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Natur, 448, 1022

Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413

Kothe, S. 2016, PhD thesis, Technische Universität Braunschweig

Kunz, M. W. 2008, MNRAS, 385, 1494

Kunz, M. W., & Balbus, S. A. 2004, MNRAS, 348, 355

Lemaster, M. N., & Stone, J. M. 2009, ApJ, 691, 1092

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Lesur, G., & Longaretti, P.-Y. 2011, A&A, 528, A17

Li, R., & Youdin, A. N. 2021, ApJ, 919, 107

Lim, J., Simon, J. B., Li, R., et al. 2024, ApJ, 969, 130

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Lyra, W., McNally, C. P., Heinemann, T., & Masset, F. 2017, AJ, 154, 146

Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001

Makwana, K. D., Zhdankin, V., Li, H., Daughton, W., & Cattaneo, F. 2015, PhPl, 22, 042902

McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2003, Natur, 422, 500

Miyoshi, T., & Kusano, K. 2005, JCoPh, 208, 315

Mulders, G. D., & Dominik, C. 2012, A&A, 539, A9

Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413

Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293

Paneque-Carreño, T., Izquierdo, A. F., Teague, R., et al. 2024, A&A, 684, A174

Perez-Becker, D., & Chiang, E. 2011, ApJ, 735, 8

Pizzati, E., Rosotti, G. P., & Tabone, B. 2023, MNRAS, 524, 3184

Riols, A., & Lesur, G. 2019, A&A, 625, A108

Rosotti, G. P., Teague, R., Dullemond, C., Booth, R. A., & Clarke, C. J. 2020, MNRAS, 495, 173

Salmeron, R., & Wardle, M. 2005, MNRAS, 361, 45

Salvesen, G., Beckwith, K., Simon, J. B., O'Neill, S. M., & Begelman, M. C. 2014, MNRAS, 438, 1355

Sano, T., & Stone, J. M. 2002, ApJ, 570, 314

Sengupta, D., & Umurhan, O. M. 2023, ApJ, 942, 74

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14

Simon, J. B., Armitage, P. J., & Beckwith, K. 2011a, ApJ, 743, 17

Simon, J. B., Bai, X.-N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013a, ApJ, 775, 73

Simon, J. B., Bai, X.-N., Flaherty, K. M., & Hughes, A. M. 2018, ApJ, 865, 10

Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., & Beckwith, K. 2013b, ApJ, 764, 66

Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685

Simon, J. B., Blum, J., Birnstiel, T., & Nesvorný, D. 2022, arXiv:2212.04509

Simon, J. B., & Hawley, J. F. 2009, ApJ, 707, 833

Simon, J. B., Hawley, J. F., & Beckwith, K. 2009, ApJ, 690, 974

Simon, J. B., Hawley, J. F., & Beckwith, K. 2011b, ApJ, 730, 94

Simon, J. B., Hughes, A. M., Flaherty, K. M., Bai, X.-N., & Armitage, P. J. 2015a, ApJ, 808, 180

Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015b, MNRAS, 454, 1117

Stoll, M. H. R., Picogna, G., & Kley, W. 2017, A&A, 604, A28

Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137

Suzuki, T. K., & Inutsuka, S.-I. 2009, ApJL, 691, L49

Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617
Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69

Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4

Walker, J., Lesur, G., & Boldyrev, S. 2016, MNRAS, 457, L39

Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114

Wardle, M. 2007, Ap&SS, 311, 35

Wardle, M., & Ng, C. 1999, MNRAS, 303, 239

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Weisheit, J. C., & Upham, R. J. 1978, MNRAS, 184, 227

Xu, R., & Bai, X.-N. 2016, ApJ, 819, 68

Yang, C.-C., Mac Low, M.-M., & Johansen, A. 2018, ApJ, 868, 27

Yang, C.-C., Mac Low, M.-M., & Menou, K. 2009, ApJ, 707, 1233

Yang, C.-C., Mac Low, M.-M., & Menou, K. 2012, ApJ, 748, 79

Youdin, A. N., & Lithwick, Y. 2007, Icar, 192, 588

Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81