Dust Dynamics in Hall-effected Protoplanetary Disks. I. Background Drift Hall Instability

Recent studies have shown that the large-scale gas dynamics of protoplanetary disks (PPDs) are controlled by nonideal magnetohydrodynamics (MHD), but how this influences dust dynamics is not fully understood. To this end, we investigate the stability of dusty, magnetized disks subject to the Hall effect, which applies to planet-forming regions of PPDs. We find a novel background drift Hall instability (BDHI) that may facilitate planetesimal formation in Hall-effected disk regions. Through a combination of linear analysis and nonlinear simulations, we demonstrate the viability and characteristics of BDHI. We find it can potentially dominate over the classical streaming instability (SI) and standard MHD instabilities at low dust-to-gas ratios and weak magnetic fields. We also identify magnetized versions of the classic SI, but these are usually subdominant. We highlight the complex interplay between magnetic fields and dust-gas dynamics in PPDs, underscoring the need to consider nonideal MHD like the Hall effect in the broader narrative of planet formation.


INTRODUCTION
Understanding the processes involved in planetesimal and planet formation is essential for unraveling the origins and evolution of our own Solar System and those around other stars (e.g., Ohashi et al. 2023;Lin et al. 2023).In the core accretion scenario, the formation of 1-100 km sized planetesimals -the building blocks of planets-is one of the most critical steps for planet formation (e.g., Chiang & Youdin 2010;Johansen et al. 2014;Nayakshin et al. 2022;Drążkowska et al. 2023).
Specifically, the growth of solid particles from micronsized grains via sticking is limited to pebbles ranging in size from millimeters to centimeters.Once the particles reach this size range, collisions often lead to either bouncing or fragmentation rather than further growth (see Blum & Wurm 2008 for a review).Moreover, the presence of gas drag can induce a swift inward migration of solid particles in the process (Whipple 1972;Weidenschilling 1977), creating a radial drift barrier (Birnstiel et al. 2010(Birnstiel et al. , 2012)).
One way to overcome these growth barriers is through the collective, self-gravitational collapse of a swarm of particles directly into planetesimals (Goldreich & Ward 1973;Youdin & Shu 2002).The particle swarm must reach a dust-to-gas mass density ratio, ϵ, significantly higher than unity (Shi & Chiang 2013).This contrasts the typical value of ϵ ∼ 0.01 found in the interstellar medium but can be reasonably expected in protoplanetary disks (PPDs, Testi et al. 2014), especially in their initial stages of evolution.Consequently, an additional mechanism is required to enhance the local dust-to-gas ratio.This mechanism may be facilitated through dust settling, radial drift, dust trapping by pressure bumps or vortices, and additional dust-gas instabilities (Chiang & Youdin 2010;Johansen et al. 2014).
On the other hand, magnetic fields are known to play a leading role in driving gas dynamics in PPDs (e.g., Williams & Cieza 2011;Crutcher 2012).The radial transport of angular momentum is mediated by magneto-rotational instability (MRI) induced turbulence (Balbus & Hawley 1991).The vertical transport of angular momentum is through magnetized disk winds (Bai & Stone 2013;Béthune et al. 2017;Bai 2017;Wang et al. 2019;Cui & Bai 2021).Moreover, non-ideal MHD effects -Ohmic resistivity, Hall effect, and ambipolar diffusion -have significant impact due to the extremely weak ionization in PPDs (Gammie 1996;Armitage 2011).It is, therefore, imperative to understand how planetesimals form in non-ideal MHD.
Dust-gas dynamics under the influence of magnetic fields have been investigated in the literature.Early works focused on the dust radial migration, vertical diffusion, and settling in MRI turbulent disks (Fromang & Nelson 2005;Fromang & Papaloizou 2006;Johansen & Klahr 2005;Johansen et al. 2006).More recent numerical simulations found dust clumping was common in ideal and non-ideal MHD-dominated disks, where the clumping is attributed to the SI (Johansen et al. 2007;Balsara et al. 2009;Tilley et al. 2010), weak radial diffusion in Ohmic resistive disks (Yang et al. 2018), or pressure bumps and zonal flows induced by the MRI turbulence (Johansen et al. 2009;Dittrich et al. 2013), ambipolar diffusion (Riols et al. 2020;Xu & Bai 2022), and the Hall effect (Krapp et al. 2018).
The Hall effect corresponds to the electron-ion drift.The effect of Hall drift on disk dynamics depends on the magnetic polarity that co-or counter-aligns with the angular momentum of the disk.The Hall effect suppresses linear MRI modes when co-aligned (Wardle 1999;Balbus & Terquem 2001), while diffusive MRI may operate when counter-aligned (Pandey & Wardle 2012;Latter & Kunz 2022).Local simulations demonstrate the nonlinear saturation of the MRI modified by Hall (Sano & Stone 2002a,b;Kunz & Lesur 2013;Bai 2015).Global simulation reveals its importance on the evolution of disk magnetic flux transport (Bai & Stone 2017).Analytical studies of MHD winds-driven accretion taken into account the Hall effect are also conducted (Wardle & Koenigl 1993;Königl et al. 2010;Salmeron et al. 2011).The Hall effect is applicable over the bulk disk from one to several tens of AU (Armitage 2015;Lesur 2021), i.e. the planet-forming regions of PPDs.
In this work, we study dust-gas interaction in a Halldominated disk.In addition to the classical SI, we are also motivated by recent work on magnetic RDIs related to MHD waves (Squire & Hopkins 2018a,b;Hopkins & Squire 2018).In the ideal-MHD limit, the fast magnetosonic, slow magnetosonic, and Alfvén wave can be in resonance with dust-gas drift (see Lin & Hsu 2022, for the 'Alfvén wave' RDI in PPDs).Incorporating the Hall effect enables two more waves, the electroncyclotron (whistler) wave and the ion-cyclotron wave (Lesur 2021), which can lead to new branches of RDIs.We indeed find such 'Hall RDIs', but these are usually sub-dominant.However, we also find a novel instability that is associated with the transport of magnetic fields by dust-induced gas flows.This non-RDI can dominate over other instabilities in dust-poor disks and could thus be important for planetesimal formation.
The paper is structured as follows.In § 2, we present the basic equations governing a dusty, magnetized PPD and specify the physical setups that are being considered.Additionally, we provide an overview of the linear problem and discuss selected analytic results in § 3.In § 4, we present numerical solutions to the linear stability problem.In § 5, we discuss the newly discovered instability mentioned above and present a toy model for it.In § 6, we present illustrative non-linear simulations of this instability.In § 7, we discuss our results in the context of PPDs and planetesimal formation.We summarize in § 8.

BASIC PHYSICS
We consider a three-dimensional (3D) disk composed of gas and dust orbiting a central star of mass M ⋆ .Cylindrical coordinates (r, ϕ, z) are centered on the star.The gas has density ρ g , velocity V , and pressure p.We assume an isothermal gas with p = C 2 s ρ g , where C s = H g Ω K is the strictly constant sound-speed, H g is gas pressure scale height, Ω K = GM * /r 3 is the Keplerian frequency, and G is the gravitational constant.The gas is threaded by a magnetic field B and we assume the Hall effect is the dominant non-ideal MHD term.
We consider a single species of dust grains that interacts with the gas through a drag force characterized by a stopping time τ s .For small grains with small τ s , dust grains are tightly -but not necessarily perfectly -coupled to the gas.In this limit, the dust grains can be modeled as a second, pressureless fluid with density ρ d and velocity W (Jacquet et al. 2011).
For the most part, we neglect gas viscosity, ohmic diffusion, and dust diffusion.However, we do include these dissipative effects in a few selected simulations for numerical stability.

Shearing Box Approximation
We consider the stability of the magnetized, dustygas on scales much smaller than the typical disk radius r.As such, we adopt the shearing box approximation (Goldreich & Lynden-Bell 1965) to model a local patch of the disk that rotates around the star at the Keplerian frequency Ω 0 = Ω K (r 0 ), where r 0 is the radius of the patch.Cartesian coordinates (x, y, z) in the box correspond to the (r, ϕ, z) directions in the global disk.
In the shearing box, Keplerian rotation appears as the linear shear flow U K = −(3/2)Ω 0 x ŷ.Defining the gas and dust velocities v and w relative to this shear flow, assuming axisymmetry (∂ y ≡ 0), and focusing on regions close to the disk midplane (z = 0) and neglecting the vertical component of stellar gravity, the equations that govern the magnetized gas and dust are: where (6) is a dimensionless measure of the radial pressure gradient from the background global disk, η H is the constant Hall diffusivity, B = B/|B| is the unit vector along the magnetic field, and µ 0 is the magnetic permeability.For clarity, hereafter we drop the subscript 0 on r 0 and Ω 0 .

Disk equilibrium and parameters
Our shearing box admits a steady-state equilibrium with constant densities, velocities, and a uniform vertical field B z0 .The velocities are given by where ∆ 2 = (1 + ϵ) 2 + St 2 , and v z = w z = 0. Here, is the Stokes number, which is assumed to be constant.Note that the presence of a constant vertical field has no impact on the equilibrium flow.
In smooth disks, η is of order h 2 g , where h g ≡ H g /r is the aspect-ratio.For convenience, we define the reduced pressure gradient which is the relevant parameter in linear theory when velocities are normalized by C s .The magnetic field strength B z0 is defined through the Alfvén speed itself characterized through the plasma beta parameter, which measures the ratio between thermal and magnetic pressure: To parameterize the Hall effect, we follow Latter & Kunz (2022) and define the Elsässer number which is used to set the Hall diffusivity η H .Note that Ha is positive (negative) when the equilibrium field is aligned (anti-aligned) with the disk rotation.In discussions, we also refer to the Hall length l H ,

Single-fluid models
Most of our analyses are based on the above explicit, two-fluid model equations.However, we also make use of a modified, single-fluid formulation that is more amenable to the numerical simulations that we present later.In this approach, the gas is approximated to be incompressible (so ρ g is a constant); we solve for the relative dust-gas velocity instead of solving for w directly; the Hall term is slightly augmented (which has no effect on linear stability); and dissipation is included for numerical stability.This single-fluid model is detailed in Appendix A. Furthermore, in the limit of tightly-coupled dust grains with St ≪ 1, it is possible to simplify the singlefluid model by approximating u with the 'terminal velocity approximation' (TVA).We describe the TVA in Appendix B.

LINEAR THEORY
We consider axisymmetric Eulerian perturbations to a variable f of the form where δf is a complex amplitude, σ is the complex growth rate, and k x,z are real wavenumbers.We take k x,z > 0 without loss of generality.Linearizing Eqs.1-5 gives: where W ≡ δρ g /ρ g , Q ≡ δρ d /ρ d , and δb = δB/ √ µ 0 ρ g .
See Lin & Hsu (2022) for a similar set of equations for the case with ohmic resistivity instead of the Hall effect.
The linearized system can be regarded as an eigenvalue problem which M is the 11-dimensional matrix representation of the right-hand-side of Eqs.19-29, and X = [W, Q, δv , δw , δb] T is the eigenvector of complex amplitudes.To solve Eq. 30, we use standard numerical methods in the eigvals package for PYTHON.

Dust-free Dispersion Relation
In the dust-free limit and approximating the gas to be incompressible (thus filtering out sound waves), the linearized equations yield the quartic dispersion relation with C 3 = C 1 = 0, see Lesur (2021).If ohmic resistivity and ambipolar diffusion are included, C 3 and C 1 become non-zero but are real.In this work, we only consider the Hall effect, in which case Eq. 31 is a bi-quadratic equation for σ 2 with coefficients C 2,0 given by where is the squared epicycle frequency with q ≡ −∂ ln Ω/∂ ln r being the local shear rate, and the Alfvén frequency ω A ≡ k z V A .Unless otherwise stated, we set q = 3/2 as appropriate for Keplerian flow, in which case κ = Ω.
According to Eq.31, a sufficient condition for instability is C 0 < 0, which implies at least one positive real root.To understand the MHD effects brought about by different negative Ha, Latter & Kunz (2022) considered three different regimes: Regime 1 corresponds to Ha < −1, i.e.Ha is sufficiently large and negative, we associate instability with the standard MRI, slightly modified by the Hall effect.Regime 2 corresponds to −1 ≤ Ha < −0.25, where instability is associated with the diffusive MRI (DMRI, Pandey & Wardle 2012).Regime 3 corresponds to −0.25 ≤ Ha < 0, i.e. small negative Ha, where MHD instabilities are suppressed.When Ha > 0, the situation becomes more straightforward: the MRI/HSI is deactivated when β z is sufficiently small.

Resonant Drag Instabilities
The presence of dust can drive drag-related instabilities when the mutual friction between dust and gas is accounted for.Squire & Hopkins (2018a,b) proposed the RDI theory as a robust framework for characterizing such instabilities.RDIs stem from the resonance between waves in the gas and the equilibrium relative drift between the gas and dust.This occurs when the condition ω gas Here, ω gas is the neutral frequency of a wave mode with wavevector k in the gas when no dust is present.
The resonance condition in our unstratified, axisymmetric shearing box can be expressed as: Several candidates for ω gas may be considered, depending on the physics involved in the pure gas.Once ω gas is specified, Eq. 34 can be used to find the relation between k x and k z of resonant or most unstable modes.

RDIs with inertial waves: classical SI
In a hydrodynamic disk with ϵ ≪ 1, the classic SI of Youdin & Goodman (2005) is an RDI with inertial waves (Squire & Hopkins 2018a).Inertial waves are most easily recovered by setting ω A → 0 and σ = −iω gas in Eqs.31-33.We find By utilizing Eq. 35 and writing u x = w x − v x (which is given explicitly by Eq.A8), the resonance condition stated in Eq. 34 can be rewritten as: That is, without dissipation, there is no minimum scale to the SI.

RDIs with Hall-MHD: No rotation
In the presence of an active magnetic field, one can consider a number of MHD waves for ω gas (Squire & Hopkins 2018a).Differential rotation (such as Keplerian shear) enables MHD instabilities such as the MRI, which complicates the application of the RDI recipe.For simplicity, we first neglect rotation, which is expected to be appropriate on small scales (large k x,z ).This limit also conveniently filters out inertial waves and allows us to focus on pure MHD effects.
Setting Ω = κ = 0 in Eq.31, we find In ideal MHD (l H → 0), we recover Alfvén waves with ω gas = ω A .Lin & Hsu (2022) indeed found RDIs associated with Alfvén waves on small scales, but they are unlikely relevant to PPDs as they only operate under nearly ideal conditions.The positive and negative solutions to Eq. 37 are known as whistlers (electron-cyclotron) modes, and ioncyclotron modes, respectively (Lesur 2021).Inserting Eq. 37 into the RDI condition Eq. 34, one can then solve for the resonant k z as a function of k x for these Hall RDIs.The explicit solutions are unwieldy and are thus relegated to Appendix C.
Here, we consider the strong Hall regime such that |l H k| ≫ 1.Then Inserting these into the RDI condition Eq. 34 and rearranging, we find for whistler waves: where in the second equality we utilized Eqs.15-16 to eliminate V A and l H in favor of Ha and β z .
For ion-cyclotron waves, we find which has the same functional form as for the classic SI, Eq. 36, since in both cases ω gas ∝ k z /k.Indeed, like the SI, Eq. 41 shows that For | Ha | = 0.5, this k x,c coincides with that of the classic SI, whence the two RDIs overlap.For | Ha | smaller (larger) than 0.5, the ion-cyclotron RDI occurs at smaller (larger) k x than the classical SI.

RDIs with Hall-MHD: Keplerian rotation
For Hall-only MHD, the general solution to Eq. 31 is given, in terms of the wave frequency, by Since RDI theory is based on neutral waves in the dustfree problem, we apply Eq. 43 to the RDI condition Eq. 34 to obtain the relation between k x and k z for resonant modes only when ω 2 gas > 0. This is met for most of the parameter space we explore, except when | Ha | is large and ω 2 gas < 0 for the MRI.

NUMERICAL RESULTS
In this section, we present numerical solutions to the eigenvalue problem specified by Eq. 30.In the first section, we vary Ha and fix other parameters to study the impact of the Hall effect on stability (Figure 1, §4.1).As alluded to above, we indeed find whistler and ioncyclotron RDIs, but they generally do not dominate over the classical SI or the MRI.However, we shall also find a non-RDI, novel "background drift Hall instability" (BDHI) that can dominate in certain regimes, which we further explore in §4.2 (Figures 2 and 3).We quote the dimensionless wavenumbers K x,z ≡ k x,z H g and normalize growth rates by Ω.
In Table 1, we summarize the parameters space and comment on key results of our calculations.With names of each mode reflecting one or more major model assumptions, e.g., mode "e02s01b4h-025" indicates its ϵ = 0.2, St = 0.1, β z = 10 4 and Ha = −0.25.Here we fixed the value of η = 0.05 in all of our calculations.The table is broken into two horizontal blocks separated by single solid lines for §4.1 and §4.2, respectively.

Hall effect on the Streaming Instability
The Hall Elsässer number Ha sets the direction and strength of the Hall effect with a given magnetic field intensity, here measured by the plasma beta parameter β z .We begin with β z = 10 4 , which is a relatively mild field strength (Cui & Bai 2021).Following Latter & Kunz (2022, see their Figure 4), we vary the Elsässer number Ha, including its sign, to adjust the possible impact of the Hall effect.For simplicity, here we fix grain size to St = 0.1 and consider a dust-poor disk with ϵ = 0.2, as considered in the 'AA' setup for the classic SI by Johansen & Youdin (2007).
Figure 1 shows growth rates of unstable modes as a function of wavenumbers.The absolute value of Ha increases from the left to right panels, signifying a weakening of the Hall effect.The upper and bottom lower rows display cases for Ha > 0 and Ha < 0, respectively, i.e. for fields aligned and anti-aligned with the disk rotation.The RDI condition for the classic SI, i.e., without MHD effects, is depicted by the solid line (Eq.36).We also plot the two solutions to Eq. 43 for RDIs associated with whistler waves (dashed) and ion-cyclotron modes (dotted) in the absence of rotation.
For cases with Ha > 0, we identify the MRI as the block of 'red-yellow' modes in the lower left corners.For Ha = 0.01, we also identify the classic SI modes associated with the solid curve; as well as the weakly growing, ion-cyclotron RDIs associated with the upper end of the dashed curve (vertical blue strip).As Ha increases, the ion-cyclotron modes shift to larger k x , as expected from Eq. 41; and by Ha = 2 their growth rates exceed the classic SI modes in the hydrodynamic limit.A narrow, horizontal band of whistler RDIs can also be identified along the dashed curve.However, regardless of the value of Ha, when Ha > 0, the system is dominated by the MRI so we do not expect of any the above dust-gas instabilities to be relevant.
For cases with Ha < 0, the situation becomes more interesting.When | Ha | ≪ 1 and the impact of the Hall effect is more pronounced (e.g.Ha = −0.01), the classic SI dominates as the system is close to the hydrodynamic limit.As | Ha | increases to a certain extent (Ha = −0.25)and the system moves slightly closer to Table 1.Selected unstable modes in Hall-dusty protoplanetary disks.List of parts modes adopted in the numerical calculation.The models in first part is presented in §4.1, the models in second are selected from our parameters space study (see in §4.2).The main parameters are listed in Columns 2 to 5: the dust to gas ratio ϵ, the Stokes number St, the plasma beta parameters βz for the vertical magnetic field, and the Elsässer number Ha to represent the Hall effect.The last column uses a comment to characterise main results of the calculations.the ideal limit, the classic SI weakens while the whistler RDI (the dashed curve) becomes pronounced.Furthermore, a new region of instability arises below it1 , with intensity on par with that of the classic SI.We will refer to it as the "Background Drift Hall Instability" (BDHI) for the reasons given below, where we further explore it.
As we increase | Ha | further (Ha = −0.5 and Ha = −2), the influence of MRI re-emerges and dominates the system.However, note that for Ha = −2 the ion-cyclotron RDI (dashed curve) attains a similar growth rate as the MRI, in which case one expects the two instabilities to operate simultaneously.Concluding this section, the majority of the new instability branches align well with the RDI theory.Moreover, at significantly negative values of Ha, MRI and ioncyclotron RDI exhibit comparable magnitudes.However, in cases with moderate and negative values of Ha, a novel instability, the BDHI, may prevail.

Parameter study
Here, we explore how the above instabilities depend on the disk parameters.In Figure 2 The classical SI manifests in the 'blue' region where Ha ≳ −0.3.We conclude this from the fact that the SI is fundamentally a non-magnetic phenomenon, so its growth should remain constant when varying magnetic parameters, and should equal to that for | Ha | → 0. This also simplifies the identification of parameter ranges where the BDHI dominates, namely the cyan region where β z ≳ 10 5 .
We next investigate the impact of dust properties.In Figure 3, we plot growth rates as a function of ϵ and St with fixed β z = 10 5 and Ha = −0.25.Similar to Figure 2, we also prepared a reference case for comparison.In the left panel of Figure 3, we set β z = ∞ as non-magnetic disk models.In such a disk, the source of instability should be the classical SI.Then, we vary ϵ from 0.01 to 3, and St from 10 − 3 to 10, depend on the right panel of Figure 2. Because BDHI could potentially assume a dominant position within this parameter space.In a dust-rich disk with ϵ ≳ 1 or with large grains of St ≳ 0.3, we find the classical SI is always the dominant instability with high growth rates (s max ∼ 0.1Ω).By contrast, in a dust-poor disk with ϵ ≲ 1 and 0.01 ≲ St ≲ 0.3, we find the BDHI dominates (orange-yellow).However, for tightly coupled grains with St ≲ 10 −2 at low dust-to-gas ratios, both the classical SI and BDHI are weak, making it indiscernible which is dominant.

BACKGROUND DRIFT HALL INSTABILITY
In §4, we identified a new branch of unstable modes -unrelated to any of the RDIs in §3.2 -that appear at low dust-to-gas ratios (ϵ < 1) and a weak field (β z ∼ 10 5 ) anti-aligned with the disk rotation with moderate Hall strength (−0.25 ≲ Ha < 0).See, for example, the dominant modes in the left panel of Figure 2. In this section, we reproduce this new instability in a toy model, which shows that it results from the advection of magnetic field perturbations relative to the dusty gas' center of mass.The fact that this new instability is not an RDI indicates that dust does not play an active role, specifically, that is through the drag forces in the linearized momentum equations.However, dust-gas drag also manifests through the equilibrium radial velocities.Here, we show that the new instability is related to the advection of magnetic fields by such an equilibrium flow.To this end, we repeat the standard calculation in Figure 1 (the second panel from the left in the bottom row, with Ha = −0.25)but artificially set the equilibrium gas radial velocity v x → 0 in the linearized induction equations Eq. 27-29.The result is shown in Figure 4 and we see that the new instability has disappeared.We, therefore, refer to it as the "background-drift Hall instability" (BDHI).

Toy model
As shown above, the BDHI arises from the advection of magnetic perturbations by the dust-induced background flow.This effect can be captured in the standard, pure gas incompressible Hall-MHD equations by adding a source term in the induction equation to transport magnetic field perturbations.The governing equations of such a toy model are: where Π is the total pressure, v ext is a velocity field prescribed below, and is the magnetic field deviation (not necessarily small) from the equilibrium, constant vertical field B z0 .In the toy induction equation Eq. 46, we augmented the Hall term slightly to facilitate our nonlinear simulations in §6 (see also Appendix A), but this does not affect the system's linear stability.
Our toy model should be viewed in conjunction with a corresponding dusty problem.We thus set v ext = v ext x, where v ext is the dust-induced, radial gas drift given by Eq. 7. Furthermore, in Eq. 45, we reduce the pressure and magnetic tension forces by a constant gas fraction f g .This reduction manifests in the dusty problem for tightly-coupled grains because dust-loading increases the gas' inertia (see Appendix B and Eq.B26).
One can regard the v in our toy model to represent the corresponding dusty-gas' center of mass velocity (which experiences no equilibrium drift).However, since the magnetic field is physically carried by the gas component only, there exists a drift between the magnetic field and the system's center of mass.This is the spirit of our toy model.

Linearized equations and dispersion relation
Linearizing the toy model about v = 0 and a constant vertical field gives: which may also be obtained from the linearized singlefluid formulation of the dusty problem (Eqs.A10-A13) by setting δϵ = St = 0 where they appear explicitly.These equations yield the dispersion relation with Without the imposed drift, v ext = C 3 = C 1 = 0, we recover the bi-quadratic dispersion relation for Hall-only MHD, see e.g.Lesur (2021) and §3.1.In the presence of the drift, the dispersion relation is a full quartic and must generally be solved numerically.Notice the Alfvén speed is reduced by dynamical coupling to dust.

An RDI-like interpretation
As the toy model is comprised only of a single fluid, the background drift Hall instability cannot be an RDI as described by Squire & Hopkins (2018a).However, it turns out that we can estimate the most unstable wavenumbers based on similar ideas, as follows.This is because our toy model also exhibits a drift: between the fluid and the magnetic field.
We consider the strong Hall regime with | Ha | ≪ 1.We first show that without the imposed magnetic drift, the system admits neutral waves.This is akin to the first step of the RDI recipe, where one identifies neutral waves in the corresponding pure gas system.
In this limit, C 3 = C 1 = 0, and Solving the dispersion relation Eq. 51 in the above limit, we find the gas wave frequencies ω gas = iσ satisfy ω 2 gas = k 2 z Ω 2 /k 2 , corresponding to inertial waves, or In the case of interest with Ha, η H < 0 and thus ω 2 gas > 0, these correspond to whistler waves modified by Keplerian rotation.Note that in the absence of rotation, Eq. 58 reduces to the positive solution of Eq. 38.
We now suppose magnetic perturbations with radial wavenumber k x are transported at a radial velocity v ext .We may then expect a resonance to occur if this transport velocity matches the phase velocity of neutral whistler waves, i.e. if ω gas /k x = v ext .This is akin to the RDI condition where the dust-gas drift velocity matches the phase speed of neutral waves in the gas.
The above RDI-like condition can be written in dimensionless form where θ ≡ K x /K z .Eq. 59 can be regarded as a quadratic equation for K 2 z as a function of K x .Defining We also require θ 2 > 3/ 2β z | Ha |ζ 2 x for real vertical wavenumbers.That is, modes must be sufficiently narrow in the radial direction compared to their vertical length scale.

Example solution
In Figure 5 we plot growth rates obtained from the dispersion relation of the toy model (Eq.51) and compare it with that obtained from the full model.Here, we set ϵ = 0.2, St = 0.1, η = 0.1; Ha = −0.25, and β = 10 5 , for which the BDHI dominates.We also over-plot the RDI-like condition given by Eq. 59.The toy model indeed reproduces the BDHI, although with slightly higher growth rates.The toy model also predicts significant growth for K x ≳ 10 3 ; while in the full model, this region has negligible growth.Nevertheless, the coincidence between the BDHI region with the curves supports an RDIlike interpretation of the instability.

DIRECT SIMULATIONS
In this section, we demonstrate the background drift Hall instability using direct numerical simulations.We use the dedalus spectral code (Burns et al. 2020) to evolve both the single-fluid model described in Appendix §A and the toy model described above.In both models, we solve for the magnetic vector potential A such that ∆B = ∇ × A, which ensures that ∇ • B = 0. Explicit dissipation is needed in spectral simulations for numerical stability.In the single-fluid model, we include viscosity and resistivity characterized by a constant kinematic viscosity coefficient ν 1 and constant resistivity coefficient χ, respectively.As this model also evolves the relative dust-gas drift u = w − v and the dust-to-gas ratio ϵ, we further include a viscous term ν 2 ∇ 2 u in the equation for u and dust diffusion D∇ 2 ϵ in the dust mass equation; where ν 2 and D are constant viscosity and diffusion coefficients, respectively.We solve for Q ϵ ≡ ϵ 0 ln (ϵ/ϵ 0 ), where ϵ 0 is the initial dust-to-gas ratio, to ensure that ϵ > 0. Our implementation of the single-fluid model in dedalus is described in Appendix A.5.
In the toy model, we only include gas viscosity (ν 1 ) and ohmic resistivity (χ) since this model does not involve u or ϵ.The dissipative toy model is described in Appendix D for ease of reference.More specifically, we evolve Eqs.D39 and D41 subject to We consider a periodic domain of size x ∈ [−L x /2, L x /2] and z ∈ [−L z /2, L z /2].We decompose all fields into N x × N z Fourier modes and evolve them using a second-order Runge-Kutta time integration.A standard dealiasing factor of 3/2 is employed.For simplicity, we set where the Re is the Reynolds number.We initialize the simulations in a steady state and then add either an eigenmode perturbation or random perturbations.

Code test
We first test our code implementation against linear theory.The linearized full and toy model equations with dissipation are given by Eqs.A10-A14 and Eqs.D42-D43 (subject to ∇ • δv = 0), respectively.We consider a system with an equilibrium ϵ = 0.2, St = 0.1, η = 0.1; Ha = −0.3, and β z = 10 5 .The Reynolds number Re = 2 × 10 5 .Recall the imposed radial flow in the toy model is set to v ext = 2ϵStηc s /∆ 2 ; which is the dust-induced, equilibrium gas radial flow in the full model.
In Figure 6, we solve the linearized equations numerically and show growth rates as a function of k x,z for the full model (left panel) and the toy model (middle panel).For comparison, we also plot growth rates in a full model with β z = 10 7 in the right panel, for which the background drift Hall instability is negligible and only the classical SI remains.
The toy model correctly isolates the background drift Hall instability from the corresponding full model, although the former slightly overestimates growth rates.The modes in the full model that extend to K x ∼ 1 and K z ∼ 10 are RDIs associated with Hall MHD.These are not captured by the toy model.The full model also captures the classical SI as a faint branch of modes below the BG-drift Hall instability.These classical SI modes survive in the β z = 10 7 case (right panel), wherein MHD effects are effectively absent.
For the dedalus runs, we consider the full and toy models with β z = 10 5 .We add a clean eigenmode solution from the linearized problem.We choose K x = 69.6 and K z = 34.3,which is approximately the most unstable mode in the full model and corresponds to the background drift Hall instability.The perturbation amplitude is normalized such that |δv y | = 10 −4 C s .We set the simulation domain to one wavelength in each direction, i.e.L x,z = λ x,z ≡ 2π/k x,z , and use a resolution of N x × N z = 64 × 64.We note that the Reynolds number associated with the box size Re box ≡ 4π 2 /K 2 Re ≃ 1300.
Figure 7 shows the evolution of the maximum velocity perturbation normalized by the Alfvén speed V A .We also plot in asterisks the expected evolution based on theoretical growth rates and find an excellent match with the simulation in the linear regime.Growth rates are also shown in the legend, which demonstrate a relative error of O(10 −3 ).

Nonlinear evolution
We now present simulations of the three models shown in Figure 6.Here, we set the box size to L x = L z = H g and use a resolution of N x × N z = 512 × 512.The range of relevant wavenumbers is delineated by the dotted boxes in Figure 7, which is obtained by assuming the maximum allowed wavelength is the box size and that a minimum of 10 collocation points are needed to resolve a wavelength.We destabilize the system by adding random perturbations to v y with a maximum amplitude of 10 −4 C s .
In Figure 8, we compare the evolution of the vertical velocity between the full model and the toy model, both with β z = 10 5 .Both models yield exponential growth with growth rates of 3.15 × 10 −2 Ω, which is somewhat smaller than the theoretical maximum growth rates of 3.36 × 10 −2 Ω (4.10 × 10 −2 Ω) in the full (toy) models.This is likely due to the reduced resolution compared to the code tests.The models saturate at similar amplitudes with |v z | ∼ 6V A . Figure 9 shows that both runs exhibit large-scale upwelling and downwellings, with small-scale substructures.The similarity between the two models indicates that the background drift Hall instability dominates the system up to the initial saturation.
However, the full model undergoes a secondary, albeit minor, growth around ∼ 330P , where |v z | reaches ∼ 10V A .We find this is associated with the development of small-scale, enhanced dust-density filaments surrounding low dust-density voids (see below).We suspect this is related to classical streaming-type instabilities operating.
We next compare the full models with β z = 10 5 and β z = 10 7 .Figure 10 shows the evolution of the maximum dust-to-gas ratio in the box and Figure 11 shows ϵ at the end of the simulations.For β z = 10 7 , where only the classic SI operates, we find the development of weak dust clumps around t ∼ 370P associated with the break up of initial SI filaments.However, these dust enhancements are not sustained and the system decays to a quasi-steady state with box-scale filaments and a negligible increase in ϵ.This case proceeds similarly to the Johansen et al. ( 2007)'s 'AA' run of the SI at low ϵ, which also exhibits transient growth.
On the other hand, the β z = 10 5 case evolves quite differently.The system first saturates to a state with max(ϵ) ∼ 0.3, which is sustained for ∼ 100 orbits.Since the toy model behaves similarly to the full model up to this point, it suggests that dust concentration is passively induced by the background drift Hall instability.However, unlike the β z = 10 7 case, dust then undergoes a secondary concentration and the system sustains max(ϵ) ∼ 0.5, or about twice the initial value, until the end of the simulation.The left panel of Figure 11 shows that the dust enhancements develop around small dust 'voids'.This again differs from the classical SI, where dust typically concentrates into filaments instead.

DISCUSSION
We find the classical SI dominates dust-rich systems with ϵ ≳ 1.Similarly, for Ha > 0, the MRI always dominates.We generally find Hall RDIs, for which Hall MHD and dust-gas interaction act in concert, to be negligible.The only exception is when Ha is large and negative, in which case the ion-cyclotron RDI co-dominates the system with the MRI, and they are distinctly associated with high and low wavenumber, respectively.Squire & Hopkins (2018b) also suggested that MHD RDIs are likely unimportant, at least in the disk midplane, be-     cause ohmic resistivity and ambipolar diffusion are both expected to dampen MHD waves (but see below).
On the other hand, the newly-discovered BDHI can dominate systems with weak fields (β ≳ 10 5 ), pebblesized grains with 10 −2 ≲ St ≲ 0.3, and low dust abun- dance (ϵ ≲ 1).The latter may be relevant to disk regions with particle stirring, as dust diffusion limits the midplane dust-to-gas ratio.Furthermore, dissipation, neglected in our main linear calculations, can suppress high-wavenumber classical SI and MRI modes, which could extend the relevance of the BDHI beyond the above regimes.
Consider a stratified disk with a vertically integrated dust-to-gas ratio Z = ϵH d /H g ≃ ϵ δ/St, where H d is the particle scale-height and δ = D/C s H g is the dimensionless diffusion coefficient (Lin 2021).For the BDHI-dominated disk presented in §6 with St = 0.1 and δ = 2 × 10 −5 (Figure 6, left panel), we find Z ≃ 0.003, which is smaller than the canonical value of 0.01.This suggests that the BDHI could facilitate planetesimal formation in dust-depleted regions of PPDs.

Application to protoplanetary disks
The new BDHI is found under a couple of conditions.First of all, the background magnetic field has a dominant vertical component.Secondly, the magnetic flux is anti-parallel to the disk rotation spin (Ha < 0).Third, the Hall effect needs to be of moderate strength (small | Ha |).We discuss the relevance of such conditions in the physical context of PPD formation under the effects of non-ideal MHD.
Several numerical and theoretical works (Lee et al. 2021;Zhao et al. 2021;Tsukamoto et al. 2015) discuss scenarios of non-ideal MHD disk formation under conditions where rotation and magnetic field are either parallel or anti-parallel.The actual behavior and the strength of the Hall effect depend on many factors such as the cosmic ionization rate, dust grain size distribu-tion, gas chemical composition, and temperature (Marchand et al. 2016).At PPD densities, it is reasonable to assume that the Hall effect is dominated by the electrons, giving a positive Hall coefficient η H .If the Hall effect is strong enough to have a dynamic effect on the formation process of the disk, then the spin-up of the disk is induced by the Hall drift regardless of the initial rotation of the disk-forming core.The resulting configuration is always anti-parallel, yielding a negative Ha independently of the initial conditions.Within the disk, there is a background radially outward drift of magnetic flux that is the result of the synergy between the Hall effect and the ambipolar diffusion, guaranteeing that B z will grow less rapidly with respect to the disk mass (Zhao et al. 2021).This also means that β z is reasonably high and will only increase with time.Typical values of β z ≳ 10 4 are easily reached.
Due to the vertical differential rotation, a strong azimuthal component of the magnetic field is generated through magnetic induction.As a consequence, the presumption of a magnetic configuration dominated by the vertical field should be valid only at the very early phase of PPD formation (Strong-B case in Lee et al. 2021, which does not last for long), or not too far away from the midplane, which is more physically relevant.Given that the BDHI proposed in this work happens at k z H g typically larger than order unity for β z = 10 4 , the development is limited to the midplane where the vertical field is dominant due to symmetry.As the PPD evolves and β z increases, the toroidal field is more developed and the region dominated by the vertical field becomes even narrower around the midplane.Meanwhile, the BDI develops at large wavenumbers and thus can still be a relevant mechanism for instability growth.
For a typical Minimum Mass Solar Nebula (MMSN) surface density profile (Hayashi 1981) and a temperature of a few hundred kelvin around a sub-Solar mass protostar, St = 0.1 corresponds roughly to millimeter-tocentimeter-sized grains within a few tens of AU.Recent observations suggest that dust grains may grow through coagulation up to millimeter sizes in such region (Liu et al. 2021).According to the giant planet instability models (e.g. de Sousa et al. 2020), this is where the giant planets were believed to have formed.The BDHI proposed in this study might be a viable channel to help the dust grains cross the centimeter barrier in such a regime.
Furthermore, the BDHI is weak when the dust-to-gas mass ratio ϵ is close to that of the typical ISM (0.01).The value ϵ = 0.2, with which the BDHI growth rate dominates, is mostly relevant for an early disk with moderate dust settling toward the midplane.As new evidence suggests that giant planets might have formed early and close-in (e.g.Manara et al. 2018;Morbidelli & Nesvorný 2020), the newly discovered BDHI could effectively be a plausible mechanism for early planetesimal formation during the class 0/I phases of PPD evolution.

Caveats and outlook
Our calculations are based on unstratified, local shearing box models with a single species of dust grains.The gas is magnetized but for the most part it is only subject to the Hall effect.How our results carry over to a more realistic setup should be addressed in future work.
In a global disk, the Hall effect can drive laminar gas accretion flows in the midplane (e.g.Bai 2017).This can lead to an 'azimuthal drift SI' (AdSI), even in the absence of a radial pressure gradient (Lin & Hsu 2022;Hsu & Lin 2022), which may facilitate planetesimal formation.We neglected this effect but it may be relevant to Hall-effected disks as these can naturally develop zonal flows that trap dust near pressure maxima (Béthune et al. 2017;Krapp et al. 2018).Future work should incorporate such a background accretion flow, for example by applying an appropriate torque onto the gas (McNally et al. 2017) and examining how it interacts with the BDHI.
Recent generalizations of the classical SI to a grain size distribution have shown that when the total dustto-gas ratio is less than unity, the instability can be significantly weakened compared to the single-size approximation adopted here (Krapp et al. 2019;Paardekooper et al. 2020;Zhu & Yang 2021).Since this is the regime where the BDHI is the most relevant, it will also be necessary to generalize the BDHI to polydisperse dust grains.
In a realistically stratified disk, the dust-rich midplane has nearly Keplerian rotation while the dust-poor, pressure-supported layers above and below rotate at a sub-Keplerian speed.Such a vertical gradient in the gas and dust rotation velocities, or vertical shear, can drive Kelvin-Helmholtz instabilities (Chiang 2008;Barranco 2009;Lee et al. 2010) and 'vertically shearing' SIs (Ishitsu et al. 2009;Lin 2021), which tend to disrupt the dust layer.Future work will also need to clarify how the BDHI operates in a stratified disk.
Finally, we have neglected ambipolar diffusion (AD), which is applicable to low gas density regions (Lesur 2021).Without an equilibrium azimuthal field, as considered here, AD behaves similarly to ohmic resistivity (Latter & Kunz 2022).Lin & Hsu (2022) did not find new phenomena when considering resistive MHD with dust; we therefore do not expect new instabilities when combining AD and dust dynamics in the case.However, in the presence of an equilibrium azimuthal field, AD can destabilize oblique modes with k x ̸ = 0 (Kunz & Balbus 2004).It remains to be seen how these AD shear instabilities are affected by dust; and whether or not AD and the background gas drift yield an analogous instability to the BDHI.

SUMMARY
In this paper, we examine the interplay between nonideal magnetohydrodynamics (MHD) and dust dynamics in PPDs.We focus on the Hall effect and dust-gas instabilities such as the streaming instability (SI) since they are directly relevant to planet-forming regions of PPDs.We consider vertical fields anti-parallel to the disk's rotation since otherwise, MHD instabilities dominate the system.
Our linear analyses confirm, in line with theoretical expectations, variations of the SI that are attributable to the Hall effect.These are associated with resonances between the equilibrium dust-gas relative drift and the phase velocity of whistler or ion-cyclotron waves in the magnetized gas.While these Hall-mediated 'resonant drag instabilities' (RDIs) do exist, they have smaller growth rates and span a narrower range of wavenumbers compared to the classical hydrodynamic SI or pure MHD instabilities.This implies that these Hall RDIs will unlikely replace or override the classical SI or MRI.
However, we also discover a novel instability attributable to the interaction of the Hall effect with the dust-induced equilibrium gas flow.While unrelated to the SI or RDIs, this new 'Background Drift Hall Instability' (BDHI) can dominate in low dust-to-gas ratio envi- ronments with weak magnetic fields.We reproduce the BDHI in a single-fluid toy model and briefly explore its nonlinear evolution using spectral simulations.We summarize all the instability modes involved in this study in Figure 12.As demonstrated in Figure 12, our analyses indicate that in anti-parallel Hall-dusty disks, the composition of instabilities arises from the classical SI, Hall-RDI (whistler and ion-cyclotron), and the newly discovered BDHI.
In a dust-poor disk, we find the BDHI drives much more significant dust enhancements than the classical SI.High-resolution simulations with larger Reynolds numbers than that considered in this work (Re ∼ 10 5 ) will be required to address whether or not the BDHI can drive dust densities to the point of gravitational collapse.If so, the BDHI would lower the super-solar metallicity threshold planetesimal formation as based on the classical SI.

AKNOWLEDGEMENT
This work was initiated at the 2022 ASIAA Summer Student Program.We thank Xue-Ning Bai and Zhaohuan Zhu for inspiring discussions and comments.YW gratefully acknowledges the support of the DUST-BUSTERS RISE project (grant agreement number 823823) for his secondment at the University of Arizona.YW further thanks ASIAA for his two-week visit to Taipei.M-KL is supported by the National Science and Technology Council (grants 111-2112-M-001-062-, 112-2112-M-001-064-, 111-2124-M-002-013-, 112-2124-M-002 -003-) and an Academia Sinica Career Development Award (AS-CDA110-M06).CC acknowledges funding from STFC grant ST/T00049X/1.Y-NL acknowledges the support of the National Science and Technology Council (111-2636-M-003-002) and the Ministry of Education Yushan Young Scholar Fellowship.

DATA AVAILABILITY
The data obtained in our simulations can be made available on reasonable request to the corresponding author.
The relative drift is then u x = u (0)  x St = −2ηrΩStf g , (B29) x St 2 = ηrΩf 2 g St 2 . (B30) The gas velocities can be read from the gas momentum equation in equilibrium, giving while u z = v z = 0.

B.2. Linearized momentum equation in the TVA
Linearizing Eq.B26 gives: where we utilized the results Notice also that the background radial gas flow does not appear in Eq.B33, as it cancels with the background radial dust-gas drift.That is the advection term on the LHS of the gas momentum equation (Eq.A2) does not enter the linear problem.However, the background radial gas flow still appears in the linearized induction equation (Eq.A3).This difference from the momentum equation drives a new instability described below.

C. RDI CONDITION IN THE ABSENCE OF ROTATION
The following are the specific forms of solutions for RDI condition when rotation is ignored ( §3.2.2, or green lines in Figure 12): . (C38) The " + " waves (eq.C37) are commonly referred to as whistler or electron-cyclotron modes, whereas those marked with " − " (eq.C38) are identified as ion-cyclotron modes.

D. TOY MODEL WITH DISSIPATION
For numerical stability, we add gas viscosity and ohmic resistivity terms to the momentum and induction equations in the toy model ( §5.2).We obtain: ) , we plot the maximum growth rates as a function of β z ∈ [10 2 , 10 6 ] and Ha ∈ [−0.5, −0.01].The left panel of Figure 2 is a dust-free reference case (ϵ = 0), in which only the MRI operates for sufficiently negative Ha (≲ −0.25).In the right panel of Figure 2, we set ϵ = 0.2 and St = 0.1.The figure shows a sharp transition in the growth rate at Ha ≃ −0.3 with the MRI pushed to more negative values of Ha but remains dominant in that regime.

Figure 1 .
Figure1.The growth rates of SI in an actively magnetized dusty-disk initialized with a vertical field in the limit of ideal MHD.The upper row is for Ha > 0, and the lower row is for Ha < 0. The black solid line is the RDI condition for the classic streaming instability without MHD effects (Eq.36).The violet dashed line and the violet dotted line are two solutions of Eq. 43.The above three lines are also marked in the upper right corner of the first panel in the upper row.The title is the parameters of dust in the model.The parameters of the magnetic field are marked in the lower right part of each panel.

Figure 2 .
Figure 2. The maximum growth rate as a function of βz (x-axis) and Ha (y-axis).For the left panel, we set ϵ = 0, i.e., dust-free and MRI-only disk, as a reference case.For the right panel, we fixed ϵ = 0.2 and St = 0.1 for each model, as shown in the title.The instability modules occupying the dominant position in each region are labeled in the corresponding regions in the right panel.

Figure 3 .
Figure 3.The maximum growth rate as a function of ϵ (x-axis) and St (y-axis).Similar to Figure 2, we also set a reference case in the left panel, with βz = ∞, i.e., non-magnetic disk.In the right panel, we fixed βz = 10 5 and Ha = −0.25 for each model, as shown in the title.Same as Figure 2, we labeled the SI occupying the dominant position in the corresponding regions in the right panel, and we marked the rough outline of BDHI using black solid lines.

Figure 4 .
Figure 4. Same as the case with Ha = −0.25 in Figure 1 but set the equilibrium gas radial velocity vx → 0 in Eq. 27-29.

Figure 6 .
Figure 6.Growth rates as a function of wavenumber in a dusty, magnetized disk with βz = 10 5 (left) and a corresponding toy model based on a modified, pure gas Hall-MHD equations (middle).The right panel is a full model with βz = 10 7 , for which MHD effects are negligible.The dotted box corresponds to the resolved wavenumbers in the dedalus simulation presented in §6.

Figure 7 .
Figure 7. dedalus simulations of the background-drift Hall instability based on the full, single-fluid model (blue solid) and the toy model (green solid).The simulations are initialized with a pure eigenmode with theoretical growth curves shown as crosses for the toy model (orange) and full model (red).

Figure 8 .
Figure 8. Evolution of the vertical velocity, normalized by the Alfén speed, in dedalus simulations of the full model (blue) and the toy model (orange).

Figure 9 .
Figure 9. Vertical gas velocities in the full model (left) and the toy model (right) at t = 300P for the simulations shown in Figure 8.

Figure 10 .
Figure 10.Evolution of the maximum dust-to-gas ratio in dedalus simulations of the full model with βz = 10 5 (blue; dominated by the background drift Hall instability) and βz = 10 7 (green; only the classical SI develops).

Figure 11 .
Figure 11.Dust-to-gas ratios at the end of the dedalus runs of the full model with βz = 10 5 (left) and βz = 10 7 (right).

Figure 12 .
Figure 12.Summary of unstable modes for model "e02s01b5h-025", relevant parameters are marked on the title.The black solid line is the classical SI, the violet dashed line and dotted line are the two resolutions of Hall-RDI, these there instabilities as shown in Figures1 and 4. The pink solid line is BDHI, as shown in Figure5.The green dashed line and dotted line are the two resolutions of Hall-RDI without Keplerian rotation.