On the nature of the resonant drag instability of dust streaming in protoplanetary disc

TL;DR

This paper investigates the resonant drag instability in protoplanetary discs, revealing how mode coupling between inertial waves and dust waves leads to dust density perturbations and growth rates dependent on dust mass fraction.

Contribution

It provides a detailed theoretical analysis of the mode coupling mechanism underlying RDI, including growth rate dependencies and the conditions for resonance and instability.

Findings

RDI growth rate scales as f^{1/2} for mode coupling.

High wavenumber RDI growth rate scales as f^{1/3} due to triple mode coupling.

Mode coupling occurs only with combined radial and vertical dust streaming.

Abstract

The recently discovered resonant drag instability (RDI) of dust streaming in protoplanetary disc is considered as the mode coupling of subsonic gas-dust mixture perturbations. This mode coupling is coalescence of two modes with nearly equal phase velocities: inertial wave (IW) having positive energy and a streaming dust wave (SDW) having negative energy as measured in the frame of gas environment being at rest in vertical hydrostatic equilibrium. SDW is a trivial mode produced by the bulk streaming of dust, which transports perturbations of dust density. In this way, settling combined with radial drift of the dust makes possible coupling of SDW with IW and the onset of the instability. In accordance with the concept of the mode coupling, RDI growth rate is proportional to the square root of the coupling term of the dispersion equation, which itself is proportional to mass fraction of dust, $f\ll 1$. This clarifies why RDI growth rate $\propto f^{1/2}$. When SDW has positive energy, its resonance with IW provides an avoided crossing instead of the mode coupling. In the high wavenumber limit RDI with unbounded growth rate $\propto f^{1/3}$ is explained by the triple mode coupling, which is coupling of SDW with two IW. It coexists with a new quasi-resonant instability accompanied by bonding of two oppositely propagating low-frequency IW. The mode coupling does not exist for dust streaming only radially in a disc. In this case RDI is provided by the obscured mechanism associated with the inertia of solids.

Tables

Table 1. Table 1: Formal types of the mode crossing for any g z > 0 subscript 𝑔 𝑧 0 g_{z}>0

type	$(k_x,Re\left[\omega \right])$ sector	group velocities	concavity
I			any value
II			any value
III
III′			${\partial }_{k_x}^2{\omega }_i>0$
IV			any value

Table 2. Table 2: Existence of the mode crossings of different types for any g z > 0 subscript 𝑔 𝑧 0 g_{z}>0 .

type	existence
I	$\forall g_x$: $k_z\le {\tilde{k}}_z$
II	:
III	: ,  :
III′	: ,  : ,
	:
IV	$\forall k_z$: $g_x>0$

Equations

$${\mathbf{U}}_0=-q{\Omega }_{0}x{\mathbf{e}}_y,$$

$$\frac{1}{{\rho }_g}\frac{\partial p_0}{\partial x}=-\frac{\partial \Phi }{\partial x}+{\Omega }_0^2(r_0+x)+2{\Omega }_0{U}_0,$$

$$\frac{1}{{\rho }_g}\frac{\partial p_0}{\partial z}=-\frac{\partial \Phi }{\partial z},$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y){\mathbf{U}}_g-2{\Omega }_0{U}_{g,y}{\mathbf{e}}_x+(2-q){\Omega }_0{U}_{g,x}{\mathbf{e}}_y+$$

$$({\mathbf{U}}_g\nabla ){\mathbf{U}}_g=-\frac{\nabla p}{{\rho }_g}+\frac{{\rho }_p}{{\rho }_g}\frac{\mathbf{V}}{t_s},$$

$$\nabla \cdot {\mathbf{U}}_g=0.$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y){\mathbf{U}}_p-2{\Omega }_0{U}_{p,y}{\mathbf{e}}_x+(2-q){\Omega }_0{U}_{p,x}{\mathbf{e}}_y+$$

$$({\mathbf{U}}_p\nabla ){\mathbf{U}}_p=\frac{\nabla p_0}{{\rho }_g}-\frac{\mathbf{V}}{t_s},$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y){\rho }_p+\nabla ({\rho }_p{\mathbf{U}}_p)=0.$$

$$\mathbf{U}\equiv \frac{{\rho }_g{\mathbf{U}}_g+{\rho }_p{\mathbf{U}}_p}{\rho },$$

$${\mathbf{U}}_g=\mathbf{U}-\frac{{\rho }_p}{\rho }\mathbf{V},$$

$${\mathbf{U}}_p=\mathbf{U}+\frac{{\rho }_g}{\rho }\mathbf{V},$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y)\mathbf{U}-2{\Omega }_0{U}_y{\mathbf{e}}_x+(2-q){\Omega }_0{U}_x{\mathbf{e}}_y+(\mathbf{U}\nabla )\mathbf{U}+$$

$$\frac{{\rho }_g}{\rho }\left\{\left(\mathbf{V}\nabla \left(\frac{{\rho }_p}{\rho }\right)\right)\mathbf{V}+2\frac{{\rho }_p}{\rho }\left(\mathbf{V}\nabla \right)\mathbf{V}\right\}=\frac{\nabla p_0}{{\rho }_g}-\frac{\nabla (p+p_0)}{\rho },$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y)\mathbf{V}-2{\Omega }_0{V}_y{\mathbf{e}}_x+(2-q){\Omega }_0{V}_x{\mathbf{e}}_y+$$

$$(\mathbf{U}\nabla )\mathbf{V}+(\mathbf{V}\nabla )\mathbf{U}+\frac{{\rho }_g}{\rho }(\mathbf{V}\nabla )\left(\frac{{\rho }_g}{\rho }\mathbf{V}\right)-$$

$$\frac{{\rho }_p}{\rho }(\mathbf{V}\nabla )\left(\frac{{\rho }_p}{\rho }\mathbf{V}\right)=\frac{\nabla (p+p_0)}{{\rho }_g}-\frac{\rho }{{\rho }_g}\frac{\mathbf{V}}{t_s},$$

$$\nabla \cdot \left(\mathbf{U}-\frac{{\rho }_p}{\rho }\mathbf{V}\right)=0,$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y){\rho }_p+\nabla (\rho \mathbf{U})=0.$$

$$f\equiv \frac{{\rho }_p}{{\rho }_g}<1.$$

$$g\equiv \left|\frac{\nabla (p+p_0)}{\rho }\right|,$$

$${\tau }_{\ast }\equiv t_s\max \{t_{ev}^{-1},{\Omega }_0\}\ll 1,$$

$${\lambda }^{-1}\equiv \frac{g{t}_s^2}{l_{ev}}\ll 1,$$

$$V\sim g{t}_s,U\sim \frac{g{t}_s}{\sqrt{{\lambda }^{-1}}},$$

$$V/U\sim \sqrt{{\lambda }^{-1}}.$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y)\mathbf{U}-2{\Omega }_0{U}_y{\mathbf{e}}_x+(2-q){\Omega }_0{U}_x{\mathbf{e}}_y+(\mathbf{U}\nabla )\mathbf{U}=$$

$$\frac{\nabla p_0}{{\rho }_g}-\frac{\nabla (p+p_0)}{\rho },$$

$$\frac{\nabla (p+p_0)}{\rho }=\frac{\mathbf{V}}{t_s}.$$

$$-\frac{{\partial }_z{p}_0}{{\rho }_g}\equiv g_z\approx {\Omega }_0^2{z}_0>0,$$

$$\eta \equiv \frac{{\Omega }_K-{\Omega }_0}{{\Omega }_K}\ll 1.$$

$$-\frac{{\partial }_x{p}_0}{{\rho }_g}\equiv g_x\approx 2\eta {\Omega }_0^2{r}_0>0.$$

$$\mathbf{U}=0,$$

$$\frac{\nabla (p+p_0)}{\rho }=\mathbf{g},$$

Full text

††pagerange: On the nature of the resonant drag instability of dust streaming in protoplanetary disc–LABEL:lastpage††pubyear: 2017

1 Introduction

Relative motion of gas and dust in protoplanetary discs drives the generic dynamical instability of gas-dust mixture due to interaction between two components via aerodynamic drag. Squire & Hopkins (2018b) recognised that instability arises each time there is synchronization between the streaming dust motion and the wave of any nature propagating in gas environment. They termed this type of instability as resonant drag instability (RDI). This gave a new look at the previously known streaming instability of Youdin & Goodman (2005), which arises due to the radial drift of solids in the midplane of protoplanetary disc. Currently, the streaming instability is a candidate for the solids growth mediator between pairwise sticking of particles and gravitational instability of dust-laden sub-disc via classic Safronov-Goldreich-Ward mechanism, see Safronov (1972) and Goldreich & Ward (1973), which ends with formation of planetesimals. This view on the role of the streaming instability was established in the work by Youdin & Johansen (2007) and Johansen & Youdin (2007) who investigated its non-linear regime, and subsequently by Johansen et al. (2007) who demonstrated that the resulting dust overdensities are gravitationally unstable. However, the subsequent simulations have shown that this scenario is sensitive to dust fraction: it requires the initial metalicity of disc be a few times higher than the solar metalicity, otherwise turbulence caused by the streaming instability leads to stirring rather than clumping of dust, see e.g. Johansen et al. (2009), Bai & Stone (2010). This constraint becomes more severe due to limitations for grain size.

Along with the known mechanisms of dust pile-up in long-living turbulent structures of disc such as zonal flows and anticyclonic vortices, for the review see e.g. Chiang & Youdin (2010), Johansen et al. (2014), Birnstiel et al. (2016), the local amount of dust in a disc can potentially be enhanced through the new vigorous instability found by Squire & Hopkins (2018a) (below SH) who termed it as settling instability, SI hereafter. Along with the streaming instability, SI is another manifestation of RDI caused by sedimentation of dust. It was revealed that, as compared with the streaming instability, SI growth rate is much higher for small grains with stopping time shorter than the Keplerian time, because its maximum value does not depend on grain size. Furthermore, SI operates at larger length-scales. These and other advantages of SI make it a promising candidate for essential link in the chain of physical processes leading to planetesimal formation, see the discussion in Section 9.2 of SH.

This work is focused on the nature of SI. The Lagrangian framework is employed to study the dynamics of small perturbations of gas-dust mixture in two-fluid approximation. Such a framework is well developed in application to perturbations in single-fluid dynamics. In particular, the Lagrangian describing the evolution of perturbations enables to introduce the energy of perturbations from the fundamental symmetry of the system with respect to translations in time. The energy of perturbations is conserved provided that the background flow is stationary, however, it is not necessarily a positive definite quantity in a moving fluid. The existence of negative energy perturbations in the flow is well established marker of instability. Indeed, any sink of energy in the system, such as e.g. flux of radiation at infinity or viscous damping, provides growth of amplitude of negative energy perturbations. Thereby, Friedman & Schutz (1978b) showed that such perturbations cause the secular instability of rotating stars with respect to gravitational radiation. The energy functional derived in terms of the Lagrangian displacement was referred to as the canonical energy, see Friedman & Schutz (1978a). The concept of the canonical energy has been used in various astrophysical applications including, for example, stability of relativistic stars, see e.g. Friedman & Morsink (1998), stability of magnetic stars, see e.g. Glampedakis & Andersson (2007), dynamic tides in stars, see e.g. Papaloizou & Ivanov (2010), Ogilvie (2013), and many others.

The physical sense of negative energy perturbations is easier to grasp noticing that the total mechanical energy of the corresponding perturbed flow being itself a positive definite quantity is less than the energy of the background flow. Probably, for the first time the concept of negative energy in application to space-charge waves in electron beams was discussed by L.J. Chu (1951), see the monograph by Briggs (1964) and the book by Pierce (2006). Kadomtsev et al. (1965) were the first to show the existence of negative energy electromagnetic waves in plasma, while the first analysis of such possibility for waves in hydrodynamical flows goes back to works by Landahl (1962) and Benjamin (1963). A bright exposition on the negative energy waves in context of fluid mechanics and various applications can be found, e.g., in the review by Ostrovskiĭ et al. (1986) and the monograph by Stepanyants & Fabrikant (1989).

One more effect, which provides the growth of amplitude of negative energy wave is coalescence, also termed as coupling, or resonant interaction with some other wave having positive energy. As soon as the former transfers the energy to the latter, the both waves experience an unbound growth. By the way, an example of such a coupling of space-charge waves with electromagnetic wave traveling along a helical coil surrounding an electron beam, which is described in the book by Pierce (2006) cited above, is quite similar to the subject of this work: roughly speaking, it is enough to replace the electron beam by the streaming of dust, the electromagnetic wave in helical coil by inertial wave (IW) in gas environment, the space-charge wave by a streaming dust wave (SDW) and finally, the spatial growth by the growth over time. A detailed description of the mode coupling and classification of its variants in application to the streaming of charged particles through a plasma can be found in Briggs (1964). As for hydrodynamic applications, a notable work was done by Cairns (1979) who discovered that the mode coupling is responsible for a classic Kelvin-Helmholtz instability. Two waves existing on a jump of fluid density in a uniform gravitational field acquire energy of different signs and couple with each other, as soon as the layers of different density are in a relative motion.

Another example, where the mode coupling is responsible for hydrodynamical instability, is the well-known Papaloizou-Pringle instability of uniform angular momentum rotating tori, see Papaloizou & Pringle (1984), Papaloizou & Pringle (1985), Papaloizou & Pringle (1987), Kojima (1986), Kojima (1989), Goldreich et al. (1986) and their citations. An exhaustive explanation of Papaloizou-Pringle instability belongs to Glatzel (1987a) and Glatzel (1987b). In a simplified two-dimensional model it was shown that small perturbations grow due to the coupling of surface gravity modes and sound modes attached to the inner and the outer boundaries of tori. In the subsequent paper Glatzel (1988) additionally investigated the mode coupling in supersonic plane shear flow.

In this work, the Lagrangian describing the linear dynamics of modes of two-fluid perturbations is introduced for the general situation of dust streaming both vertically and radially in protoplanetary disc. The problem is considered in the terminal velocity approximation (TVA). It is found that the energy of certain (neutral) mode of such perturbations akin to SDW can be negative. SDW itself is a trivial mode, which also can have negative energy and represents perturbations of dust density transported by the bulk stream of dust through the unperturbed gas environment in the absence of dust back reaction on gas. Physically, this is the bulk streaming of dust that gives rise to a reservoir of negative energy in the form of SDW, from which IW can draw the unlimited amount of energy. The latter occurs as soon as the phase velocity of SDW is sufficiently close to the phase velocity of IW. As the dust back reaction on gas is taken into account, both modes coalesce with each other giving rise to the pair of damping and growing coupled modes with strictly equal phase velocities and the vanished energy, thus, the band of SI appears. The strength of the mode coupling is determined by the coupling term in the dispersion equation, which is proportional to the mass fraction of dust in gas-dust mixture. The growth rate of SI takes maximum value at resonance between IW and SDW, where IW and SDW acquire equal phase velocities. It is proportional to the square root of the coupling term and, correspondingly, to the square root of the dust fraction. At the same time, it does not depend on the size of the particles, see also SH.

The analysis shows that generally, there are various types of resonances between IW and SDW, also called the mode crossings. Not all of them result in the mode coupling. If the energy of SDW is positive, the modes fall into the avoided crossing in the vicinity of resonance, which provides no SI. That is, the dust back reaction on gas breaks the dispersion curves describing SDW and IW and produces two separate branches of neutral modes. At the same time, there are resonances which exist in the high wavenumber limit. It turns out that the corresponding mode coupling gives rise to an unbounded growth rate of SI found by SH. For sufficiently high dust fraction the unbounded growth rate is proportional to its cube root, which is explained by the triple coupling of SDW and two low-frequency IW. In addition, the low-frequency IW bond to each other due to the dust back reaction on gas and give a quasi-resonant instability of a new type, which may be as strong as SI when the dust settling dominates the dust radial drift.

Finally, it is shown that there is no RDI within TVA when vertical gravity is negligible and dust streams only in the radial direction. The resonance between SDW and IW existing in this case provides the avoided crossing rather then the mode coupling. The known streaming instability studied since the work of Youdin & Goodman (2005) is found beyond TVA. As was shown by SH, its growth rate is proportional to the square root of the dust fraction, thus, it is also RDI. However, this sort of RDI arises with the account of terms describing the inertia of solids, which are of the next order in the particle stopping time. The physical mechanism responsible for streaming instability of Youdin & Goodman (2005), though resonant in nature, has nothing to do with the mode coupling.

The study starts from general equations for local dynamics of gas-dust mixture in protoplanetary disc and the subsequent review of TVA. The main part of the paper is focused on the particular model of the perturbed gas-dust mixture constructed in TVA for homogeneous background with dust streaming vertically and radially in a disc. In the last Section the streaming instability of dust drifting in radial direction is studied beyond TVA.

2 General equations for local dynamics of gas-dust mixture in a disc

In order to consider a small patch inside protoplanetary disc, the local Cartesian coordinates $x,y,z$ are introduced to represent, respectively, radial $r$, azimuthal $\varphi$ and vertical $z$ directions. The radial distance is measured ourwards with respect to the host star, while the vertical one is collinear with disc rotation axis. The centre of a reference frame is located at some point $r_{0},\varphi_{0},z_{0}>0$ above the disc midplane rotating with angular velocity $\Omega_{0}$ around the host star. Accordingly, there are $x\equiv r-r_{0},\,y\equiv r_{0}(\varphi-\varphi_{0}),\,z\to(z-z_{0})$ and it is assumed below that $\{x,y,z\}\ll h\ll r_{0}$, where $h$ is the disc scaleheight. This corresponds to the small shearing box approximation, see Goldreich & Lynden-Bell (1965) and the Appendix A of Umurhan & Regev (2004) for detailed derivation of fluid equations. In this work the same approximation is used to address the two-fluid dynamics of gas-dust mixture. Both gas and dust are considered as fluids with velocities, respectively, ${\bf U}_{g}$ and ${\bf U}_{p}$ measured with respect to some reference velocity ${\bf U}_{0}$. The latter describes the stationary circular shear motion of a single fluid in the gravitational field of the host star,

[TABLE]

obeying the following radial and vertical balance

[TABLE]

[TABLE]

where $\Phi$ is the Newtonian point-mass gravitational potential, $p_{0}$ and $\rho_{g}$ are pressure and density of the fluid, which can be obtained given the equation of state. In eq. (1) the shear rate $q$ is assumed to be constant normally close to the Keplerian value $q=3/2$ in protoplanetary discs. It can be shown, see the Appendix A of Umurhan & Regev (2004), that in the small shearing box approximation the Euler equation for gas reads

[TABLE]

where $U_{g,x}$ and $U_{g,y}$ are, respectively, $x$- and $y$-projections of ${\bf U}_{g}$, ${\bf e}_{x}$ and ${\bf e}_{y}$ are, respectively, $x$- and $y$-orts, while $p$ introduces addition to pressure, $p_{0}$, caused by the influence of dust. The last term in eq. (4) originates from aerodynamic drag of dust particles moving through the gas with the relative velocity ${\bf V}\equiv{\bf U}_{p}-{\bf U}_{g}$ and mass density $\rho_{p}$. The dust particles stopping time $t_{s}=const$ is a convenient quantity to parametrise aerodynamic drag, see Whipple (1972). On the scale much shorter than $h$ the low-frequency dynamics of gas is vortical, thus, eq. (4) is accompanied by the condition of the divergence-free motion

[TABLE]

Note that in the case of $\rho_{p}=0$ equations (4-5) follow from equations (A.28a-d) of Umurhan & Regev (2004) for unstratified fluid, i.e. provided that their $\rho^{\prime}$=0.

Similarly, the local dynamics of the pressureless fluid, which mimics the solids suspended in gas environment, is described by the following equation

[TABLE]

where $U_{p,x}$ and $U_{p,y}$ are, respectively, $x$- and $y$-projections of ${\bf U}_{p}$. The continuity equation for dust reads

[TABLE]

Again, eqs. (6-7) follow from equations (A.15-A.18) of Umurhan & Regev (2004) provided that their $P=0$, $\rho\to\rho_{p}$ and aerodynamic friction is added in its right-hand side (RHS). As it should be according to the Newton’s third law, the aerodynamic drag of dust acting on gas of unit volume equals to the aerodynamic friction of gas acting on dust of unit volume taken with the opposite sign.

In the case when dust is tightly coupled to the gas it is convenient to rewrite equations in terms of the center-of-mass velocity, see Youdin & Goodman (2005)

[TABLE]

where $\rho\equiv\rho_{g}+\rho_{p}$ is the total density of gas-dust mixture.

Since

[TABLE]

[TABLE]

one arrives at the new equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

2.1 Ordering of the terms entering the equations

Let the characteristic time- and length-scales of gas-dust mixture dynamics be $t_{ev}$ and $l_{ev}$, respectively, while the dust mass fraction

[TABLE]

Hereafter, the condition (15) implies that gas-dust mixture is not dominated by dust. If the absolute value of the specific pressure gradient, which governs the dynamics of gas-dust mixture, is

[TABLE]

then the restrictions

[TABLE]

[TABLE]

greatly simplify equations for gas-dust dynamics. More exactly, the terms in the left-hand side (LHS) of eq. (12) become small compared to the terms in its RHS, while the terms $\sim O(V^{2})$ become small compared to the rest terms in LHS of eq. (11). This can be justified as follows.

If the restrictions (16-17) are valid, the specific pressure gradient entering the RHS of eq. (12) is balanced by the leading term of this equation, which introduces aerodynamic drag. Dividing eq. (11) and eq. (12) by $g$ one finds that each of these equations consists of the dimensionless terms of various order in the small $\tau_{*}$ and $\lambda^{-1}$. So, the terms in the RHS of eq. (11) and eq. (12) as well as the gradient term $\sim U^{2}$ in the LHS of eq. (11) are of the zero order in both $\tau_{*}$ and $\lambda^{-1}$. The inertial terms $\sim U$ in the LHS of eq. (11) differ from unity by factor $\tau_{*}/\sqrt{\lambda^{-1}}$. Next, the inertial terms $\sim V$ in the LHS of eq. (12) are of the order of $\tau_{*}$, whereas the gradient terms $\sim UV$ in the LHS of eq. (12) and the gradient terms $\sim V^{2}$ in the both of these equations are, respectively, of the order of $\sqrt{\lambda^{-1}}$ and $\lambda^{-1}$. Additionally, the following order-of-magnitude relations are valid

[TABLE]

so that

[TABLE]

In order to study the instability of cooperative settling and radial drift of the dust in Section 3, all terms of the orders of $\tau_{*}$, $\sqrt{\lambda^{-1}}$ and $\lambda^{-1}$ in eqs. (11) and (12) are omitted, which corresponds to what is referred to as TVA, see Youdin & Goodman (2005). Physically, this means that the inertial forces acting on solids in the frame comoving with gas are small compared to drag force and effective gravity force measured in this frame. Note that the latter equals to the pressure gradient. As eqs. (16) and (17) imply, TVA is the case of sufficiently small solids with stopping time much shorter than the dynamical time-scale of the problem, as well as for length-scales of gas-dust mixture dynamics much longer than the solids stopping length $\sim gt_{s}^{2}$.

Later on, in Section 4 the full set of eqs. (11-14) is reconsidered in order to recover RDI in the case there is only the radial drift of the dust.

3 RDI within TVA

From now on, eqs. (11-12) read

[TABLE]

[TABLE]

This Section deals with the analysis of the particular linear solution of the set of eqs. (18), (19), (13) and (14). In what follows, where necessary, the usual dimensionless stopping time is used $\tau\equiv t_{s}\Omega_{0}$.

3.1 Stationary solution

In the reference frame rotating with $\Omega_{0}$, the local gravitational acceleration has both vertical and radial components. In the geometrically thin disc, the former is defined with the help of eq. (3),

[TABLE]

and may be considered as a constant value inside the small box. Note that in the last equality in eq. (20) the Keplerian value of angular velocity, $\Omega_{K}$, is replaced by $\Omega_{0}<\Omega_{K}$. This is justified by the small difference between $\Omega_{0}$ and $\Omega_{K}$, which is of the order of $(h/r_{0})^{2}$. This difference is usually parametrised by the dimensionless

[TABLE]

Eq. (2) gives the following relation between the radial (effective) gravitational acceleration and $\eta$:

[TABLE]

Depending on the height of the box above the disc midplane, $g_{x}$ and $g_{z}$ can be of any ratio with each other.

The most simple solution of eqs. (18), (19), (13) and (14) describing the dust settling is

[TABLE]

[TABLE]

[TABLE]

where ${\bf g}=-g_{x}{\bf e}_{x}-g_{z}{\bf e}_{z}$. Hence, the velocity of dust streaming through the gas environment is directed along the (effective) gravitational acceleration in the box.

Finally, eq. (13) along with eq. (24) implies that

[TABLE]

on the local scale considered here.

Eqs. (22-25) is the local variant of known Nakagawa et al. (1986) solution.

3.2 Equations for axisymmetric perturbations

The reduced set of equations describing the dynamics of gas-dust perturbations with the background (22-25) in the limit of small $f\ll\leavevmode\nobreak\ 1$ is derived in the Appendix A. The appropriate state variables are

[TABLE]

[TABLE]

$u_{z}$ and $\delta$, where $u_{x},u_{y},u_{z}$ are the components of the Eulerian perturbation of the centre-of-mass velocity and $\delta$ is the relative perturbation of the dust density, see the Appendix A for the details.

The basic equations take the form:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\kappa^{2}\equiv 2(2-q)\Omega_{0}^{2}$ is the epicyclic frequency squared.

Let the gas-dust perturbation have the form of the plane wave

[TABLE]

That is, $\hat{\chi}_{i}$ are complex Fourier amplitudes of $\chi_{i}$. In eq. (30) ${\bf kx}=k_{x}x+k_{z}z$, where $k_{x}$ and $k_{z}$ are wavenumbers, respectively, along radial and vertical directions in a disc. The frequency $\omega$ is in general a complex value, so that $\Im[\omega]>0$ implies the exponential growth of wave amplitude, i.e. the onset of SI.

Eqs. (26-29) yield

[TABLE]

[TABLE]

[TABLE]

[TABLE]

3.3 Energy of modes of gas-dust perturbations: derivation from variational principle

For class of neutral modes, when $\Im[\omega]=0$ and the amplitudes are constant, it is possible to formulate the variational principle describing the dynamics of small perturbations.

Let the modal solution have the form

[TABLE]

where tilded quantities are real and the phase $\theta\equiv-\omega t+{\bf kx}$. The tilded quantities in eq. (35) satisfy the set of eqs. (31-34), where the following replacement is made

[TABLE]

According to the averaged variational principle for modes, see Whitham (2011), there is a Lagrangian

[TABLE]

where $\tilde{\chi}_{i}\equiv\{\tilde{\varpi},\tilde{\phi},\tilde{u}_{z},\tilde{\delta}\}$ and $\partial_{i}\equiv\{\partial_{t},\partial_{x},\partial_{z}\}$, such that the corresponding Euler-Lagrange equations are equivalent to the set of eqs. (31-34) for modes of perturbations provided that the action

[TABLE]

with $d^{2}{\bf x}\equiv dxdz$ be stationary with respect to arbitrary variations of $\tilde{\chi}_{i}$ and $\theta$ in the small shearing box. The Lagrangian density is

[TABLE]

where

[TABLE]

and

[TABLE]

The subscript ‘0’ denotes the contribution to the full Lagrangian responsible for the dynamics of gas-dust perturbations with no account for the dust back reaction on gas. In the absence of dust, $L_{0}$ describes single-fluid vortical perturbations in rotating gas. The subscript ‘1’ denotes additional contribution to $L$ responsible solely for the dynamics of dust in terms of the evolution of its density. The first term in curly braces in $L_{1}$ is the cross-term, which describes both the action of gas on solids and the back reaction of solids on gas via the aerodynamic drag.

It can be checked that eqs. (31-34) correspond to the following Euler-Lagrange equations

[TABLE]

Also, it is straightforward to check that equations (40) yield an alternative form of the dispersion equation (cf. eq. (49) below), which is nothing but

[TABLE]

The remaining Euler-Lagrange equation is

[TABLE]

Note that ${L}$ is known up to an arbitrary (dimensional) constant factor since it describes the linear problem. The sign of ${L}$ is chosen so that the energy of IW propagating in rigidly rotating fluid (i.e. in the case $f=q=0$) be positive, see eq. (46) along with eq. (53) below.

A conserved quantity associated with the invariance of ${L}$ with respect to translations in time is the energy of neutral mode $E\equiv\int{\hat{E}}d^{2}{\bf x}$, where

[TABLE]

is the energy density of mode. Eq. (42) shows that $E$ changes due to the surface flux of

[TABLE]

Hence, this quantity is associated with the energy flux of mode, see Whitham (2011).

The energy density of mode incorporates the basic contribution, which is similar to the energy density of gas perturbations in rotating flow in the absence of dust, and the additional contribution related to perturbations of dust density

[TABLE]

Explicitly,

[TABLE]

where $k^{2}=k_{x}^{2}+k_{z}^{2}$, and

[TABLE]

For what follows, see Sections 3.7, 3.9 and 3.10, it is convenient to represent $\hat{E}$ in the form:

[TABLE]

In the limit $f\to 0$ $\hat{E}\to\hat{E}_{0}$, which is positive definite111This is true in centrifugally stable flows only. However, $\hat{E}_{1}$ has any sign depending on the signs of wavenumbers and frequency of mode, as well as on the sign of combination in the brackets of eq. (47). As the dust fraction becomes non-negligible, the balance of terms in eq. (45) for the corresponding mode frequency $\omega(k_{x},k_{z})$ determines whether $\hat{E}$ is either positive or negative quantity.

In the particular case of the dust settling, i.e. when $g_{x}\to 0$, it is possible to do a straightforward generalisation of the above variational principle onto perturbations of general form, see Appendix B. It can be shown that in the general case of dust streaming vertically and radially in a disc the variational principle of the form (131) should contain the higher derivatives of $\chi_{i}$. The development of such a formalism is beyond the scope of this study. Note that in the case $g_{x}\to 0$ eq. (48) becomes identical to eq. (138) provided that $\Im[\omega]=0$.

3.4 Dispersion equation

As it follows from eqs. (31-34), $\omega$ obeys an equation

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

with $\omega_{i}\equiv-(k_{z}/k)\kappa$ and $\omega_{p}\equiv-t_{s}g_{x}k_{x}-t_{s}g_{z}k_{z}$.

For $f\to 0$, eq. (49) splits into two separate dispersion equations

[TABLE]

and

[TABLE]

describing, respectively, two IW222Referred to as ‘epicyclic mode’ in SH propagating in opposite directions and one else mode, which is referred to as SDW. Indeed, eq. (53) can be obtained from eqs. (31-34) setting $\hat{\delta}=0$. Since the dust fraction is negligible, the remaining variables describe perturbations of gas only. The set of equations (31-34) in this particular limit is identical to those that describe IW in rigidly rotating fluid, see e.g. Landau & Lifshitz (1987), paragraph 14. In astrophysical context, one should take into account the more general case of shearing rotating flow and replace $2\Omega$, which is the epicyclic frequency for rigid rotation, by $\kappa$, see e.g. Balbus (2003), paragraph 3.2.3. For this reason, the modes described by eq. (53) are referred to as the IW below333Note that according to eq. (111) such IW produce a non-zero perturbation of relative velocity.. On the contrary, eq. (54) follows from eqs. (31-34) provided that $\hat{\varpi}=\hat{\phi}=\hat{u}_{z}=0$, which means that gas environment remains unperturbed, whereas arbitrary initial perturbations of dust are transported by the bulk streaming of dust with the velocity (24). The phase velocity of SDW is equal to the scalar product of the dust velocity and the normal to the wavefront.

Accordingly, as long as the dust fraction is negligible, $\hat{E}_{0}$ introduces the energy of IW, see eq. (46), whereas $f\hat{E}_{1}$ introduces the energy of SDW, see eq. (47).

3.5 The mode coupling

The form of the dispersion equation described in the previous Section has been examined in the dynamics of waves in shear flows, see Fabrikant et al. (1998). Cairns (1979) suggested that the solution of eq. (49) with small RHS can be interpreted as the coupling of modes having their own dispersion relations $D_{g}$ and $D_{p}$. The RHS of eq. (49) is referred to as coupling term.

The coupling term originates from the product of the last two terms in eq. (31) and the last term in eq. (34). In turn, the former comes from the last terms in RHS of eqs. (117) and (119), which are proportional to gravitational acceleration in a disc. Physically, this term is the deviation of the weight of solids acting on gas, which emerges due to perturbation of the dust density. At the same time, the last term in eq. (34) comes from the first term in RHS of eq. (120), which exists as soon as the flow is rotational, see eq. (126). Indeed, if it were not for the rotation, the pressure maxima would be absent in the perturbed gas environment. In the linear problem, the subsonic pressure maxima, which are necessary both for existence of IW and for dust clumping, arise due to the action of Coriolis force. The similar situation was noted by Jacquet et al. (2011) in context of the streaming instability of Youdin & Goodman (2005), see their discussion about the geostrophic balance. However, as it becomes clear from the study of SH, see also Section 4 of this paper, Jacquet et al. (2011) considered the non-resonant instability of gas-dust mixture different from RDI.

The coupling causes the deviation of the mode frequencies by a small value $\Delta\equiv\omega-\omega_{c}$ at the mode crossing. The mode crossing is the point in the phase space, where the phase velocities of modes are equal to each other, thus, eqs. (53-54) are jointly satisfied and

[TABLE]

Then,

[TABLE]

[TABLE]

Substitution of these approximations into eq. (49) yields the following solution

[TABLE]

This general result is employed in Sections 3.7-3.9 to show a simple way to reproduce the analytical estimations of SI growth rate obtained by SH.

A more subtle thing is that the sign of the expression under the square root in eq. (55) follows the sign of SDW energy taken at the mode crossing, see eq. (47). Therefore, as there is a mode crossing between SDW and IW, then SI occurs if and only if SDW has negative energy. Since the energy of IW is always positive, see eq. (46), this implies that the instability occurs as soon as the coupling waves have the energy of different signs. Following Cairns (1979), this allows for clear interpretation of SI as the resonant coalescence of waves in the gas-dust mixture, which provides an exchange with energy between the waves. The total energy of coupled mode formed by the coalesced waves remains constant, while the non-zero energy flux between the coalesced waves provides the growth of their amplitudes as soon as the energy flows from the negative energy wave to the positive energy wave. On the contrary, if the energy flows in the opposite direction, a coupled mode is damping. These cases represent the pair of complex conjugate solutions of the dispersion equation (49) approximately given by eq. (55).

The following Sections provide the details of the mode coupling during the dust settling or/and the dust radial drift. The frequencies and the energy of modes of gas-dust perturbations for finite non-zero $f>0$ are obtained from the accurate solution of eq. (49). Hereafter it is assumed that $k_{z}>0$. Eq. (49) indicates that the case $k_{z}<0$ is reduced to the previous one by the replacement $k_{x}\to-k_{x}$ and $\omega\to-\omega$.

3.6 Units

Wherever it is needed below, the frequencies are measured in units of $\Omega_{0}$ and the wavenumbers are measured in units of the inversed artificial length-scale, $L\lesssim h$. Accordingly, the gravitational acceleration is measured in units of $L\Omega_{0}^{2}$. In all particular calculations $q=3/2$ and $\tau=0.1$.

3.7 Particular case of the dust vertical settling

Here it is assumed that the box is located sufficiently far above the disc midplane so that $g_{z}$ solely controls the streaming of dust, cf. eqs. (20) and (21). The corresponding restriction reads

[TABLE]

An accurate solution of eq. (49) for particular values of parameters is presented in Fig. 1. The problem is symmetric with respect to $k_{x}\to-k_{x}$, so the case $k_{x}>0$ is considered only. Note that only two branches of $\omega$ with $\Re[\omega]<0$ are shown, while there exists a third one with $\Re[\omega]>0$. The corresponding profile of the energy density of mode, $\hat{E}$, normalised by $|\varpi|^{2}$ is shown on the bottom panel of Fig. 1.

As far as the dust fraction is negligible, $f=0$, there is no instability of gas-dust mixture. As expected, for $f>0$ SI sets on in the vicinity of the point, where the frequencies of two branches evaluated for $f=0$ are equal to each other. This is the mode crossing point. For the particular values of parameters taken in Fig. 1 it is located at $k_{x}=6$. The non-negligible dust fraction causes a small deviation of the branches corresponding to IW and SDW in such a way that in some band around the mode crossing the waves acquire identical phase velocities and increment/decrement, which takes maximum value approximately at the mode crossing, compare the curves marked by (1) with the curves marked by (2) and (3) in Fig. 1. Strictly speaking, the waves represented by curves (2) and (3) are the specific kind of waves propagating through the gas-dust mixture with dust settling to the disc midplane. Outside the band of instability these waves are akin to IW and SDW having, respectively, positive and negative energy densities. As one approaches the band of instability, the absolute value of the energy density of the both branches vanishes. At the same time, the neutral waves coalesce and give birth to a pair of growing and damping waves as seen in Fig. 1. The energy of such waves vanishes, since this is the only way for the energy to be conserved for alternating amplitudes of waves, see the Appendix B.2.1 and also Friedman & Schutz (1978a). Similar picture occurs in some applications from physics of hydrodynamical instability of shear flows, where it is regarded as a mode coupling discussed above in Section 3.5.

According to eqs. (53-54) with $g_{x}=0$ the mode crossing occurs as soon as the condition

[TABLE]

is true. Once (57) is satisfied, the phase velocities of IW and SDW are equal to each other $\omega_{c}=\omega_{i}=\omega_{p}$.

As one follows derivation of eq. (57) starting from the general equations (110-113), it becomes clear that the LHS of eq. (57) is nothing but the bulk settling velocity of dust given by eq. (24), also cf. eq. (4.4) of Squire & Hopkins (2018a). Therefore, the latter defines the characteristic wavelength of the mode coupling and the band of SI. As has been already recognised by Squire & Hopkins (2018a), the growing gas-dust perturbations of SI have larger length-scales than that of the streaming instability of Youdin & Goodman (2005). Clearly, this is because in the geometrically thin disc settling velocity can be higher than the velocity of radial drift by factor $\sim(h/r_{0})^{-1}$.

Equation (55) applied to the case $g_{x}\to 0$ yields

[TABLE]

where eq. (57) is taken into account. The estimate (58) is in accordance with eq. (5.13) of SH. For small value of $f$ used in Fig. 1 this analytical increment/decrement virtually coincides with an accurate solution of eq. (49). It is also equivalent to eq. (140) in the leading order in $f^{1/2}$.

Additionally, it should be noted that as far as the main restriction of solids tightly coupled to gas is valid, see eq. (16), the condition (57) fits well within the second restriction of the model, see eq. (17). Indeed, in the case of the Keplerian shear eq. (57) implies that

[TABLE]

which corresponds to the length-scale larger than the stopping length of solids settling to the disc midplane444cf. similar restriction obtained by Jacquet et al. (2011) who considered the streaming instability. $\sim g_{z}t_{s}^{2}$ by factor $\sim\tau^{-1}\gg 1$. The main restriction (16) is valid as far as $\Re[\omega]\sim\Omega_{0}$ and $\Im[\omega]\lesssim\Omega_{0}$.

3.8 Introducing the dust radial drift

As the box gets down to the disc midplane, the bulk streaming of dust is determined by combination of vertical and radial components of gravitational acceleration. The symmetry of the problem with respect to $k_{x}\to-k_{x}$ breaks and the general picture of the mode crossing (just ‘crossing’ below in this Section) between IW and SDW becomes more complex. It is shown in Fig. (2) on the plane of $(k_{x},-\Re[\omega])$. For the specific value of $k_{z}$ there are two certain branches of the oppositely propagating IW given by eq. (53), whereas the position of the straight line representing SDW, see eq. (54), is determined by the changeable $g_{x}$ and $g_{z}$. This line is horisontal for $g_{x}=0$, while it crosses the origin of coordinates, i.e. ($k_{x}=0,-\Re[\omega]=0$), for $g_{z}=0$, see the lines (1) and (4) in Fig. 2, respectively. In the former limiting case, there can be up to two crossings, whereas in the latter limiting case there are always two crossings, each one with its own branch of IW. In the general situation of $g_{x}>0$ and $g_{z}>0$, there can be from two and up to the four crossings, see the lines (2) and (3) in Fig. 2. Note that the lines of SDW shown in Fig. 2 intersect each other at the crossing considered for example in Fig. 1.

The condition of the crossing becomes

[TABLE]

where ‘+’ and ‘-’ are taken for $\omega_{c}<0$ and $\omega_{c}>0$, respectively.

In order to facilitate the understanding, it is appropriate to introduce the formal classification of the crossings, see table 1 and Fig. 2. The conditions for existence of the crossings of different types is collected in the table 2.

There are a few quantities that delimit the whole space of parameters with respect to the existence of crossings of various type. First, the value of $g_{x}=g_{x}^{{}^{\prime}}$, which restricts the existence of type-II crossings, is defined by the condition that the lines of SDW and IW do not touch each other for $g_{x}>g_{x}^{{}^{\prime}}$. It reads

[TABLE]

As soon as $g_{x}<g_{x}^{{}^{\prime}}$, there are two values of $k_{z}$, $k_{z}^{{}^{\prime\prime}}>k_{z}^{{}^{\prime}}$, corresponding to the positions of SDW line tangent to IW curve. It is straightforward to obtain them for given $g_{z}$ and $g_{x}<g_{x}^{{}^{\prime}}$, though the corresponding expressions are omitted here to save space. Note that in the limit $g_{x}\to g_{x}^{{}^{\prime}}$

[TABLE]

where

[TABLE]

For $k_{z}=\tilde{k}_{z}$ the crossing occurs at $k_{x}=0$, where the ‘negative’ branch of IW attains its minimum. Oppositely, as $g_{x}\to 0$ $k_{x}^{{}^{\prime}}\to 0$ and $k_{x}^{{}^{\prime\prime}}\to\tilde{k}_{z}$.

Next, for a higher $g_{x}^{{}^{\prime}}<g_{x}<g_{x}^{{}^{\prime\prime}}$, where

[TABLE]

type-III and type-III’ crossings co-exist in the range $k_{z}>\tilde{k}_{z}$ with the transition into each other at

[TABLE]

As $g_{x}\to g_{x}^{{}^{\prime\prime}}$, $\tilde{k}_{z}^{{}^{\prime}}\to\infty$ and type-III crossing disappears for all $k_{z}$. Also, note that in the limiting case $g_{x}=g_{x}^{{}^{\prime}}$

[TABLE]

A caveat should be done, that in this study of SI the distinction between type-III and type-III’ crossings looks completely formal. Still, it is shown here for the sake of rigorous analysis.

The formal picture presented above leads to the following summary:

• •

In the case of small radial drift, $0<g_{x}\ll g_{z}$, type-I (or type-III’) and type-II crossings are confined in the range $k_{z}\lesssim\tilde{k}_{z}$, while type-III crossing exists for unlimited $k_{z}$. For type-III crossing $|k_{x}|\gg k_{z}$.

• •

As radial drift and vertical settling are comparable to each other, $g_{x}\sim g_{z}$, type-II crossing disappears, while the whole range of $k_{z}$ is shared between type-I and type-III(’) crossings for $k_{z}\leq\tilde{k}_{z}$ and $k_{z}>\tilde{k}_{z}$, respectively.

• •

In the case of dominating radial drift, $g_{x}\gg g_{z}$, the whole range of $k_{z}$ is shared between type-I and type-III′ crossings, while the zone of type-III′ crossing shifts to high $k_{z}$ as $\tilde{k}_{z}$ increases with decreasing $g_{z}$. For type-III′ crossing $|k_{x}|\ll k_{z}$.

• •

Type-IV crossing exists in all mentioned cases.

Prediction of SI for the crossing of any type can be done using the rule revealed in Section 3.5. Examination of the energy density of SDW, see eq. (47), leads to the conclusion that type-I crossing produces SI for $k_{x}/k_{z}>g_{x}/g_{z}$ only. In the opposite case, the mode coupling is absent. Clearly, the greater the radial drift, the more severe becomes this restriction on SI, see Section 3.9 for details. On the contrary, type-II and type-III(’) crossings produce SI with no regards to the parameters, whereas, by the same reasoning, at the type-IV crossing the mode coupling is always absent, see Section 3.10.

According to the general formula (55), the approximate growth/damping rate at the crossing of any type producing SI is

[TABLE]

which must be evaluated under the condition (59). In the case of the crossing, which does not produce SI, eq. (65) provides an estimation for real correction to $\omega_{c}$. Eq. (65) is in accordance with eq. (5.12) of SH.

3.9 Avoided crossing of modes

In this Section type-I crossing is considered more closely. It can be seen that for $g_{x}>0$ the coupling term given by eq. (52) changes its sign at some $k_{x}>0$. For smaller $k_{x}$, the mode coupling ceases even for non-zero dust fraction and is replaced by the avoided crossing, see e.g. Stepanyants & Fabrikant (1989).

The condition $\epsilon=0$ along with the mode crossing condition (59) yields that for the given $g_{x}$ and $g_{z}$ the mode coupling ceases at

[TABLE]

[TABLE]

Accordingly, for $k_{z}<k_{z}^{(AC)}$ the type-I crossing causes the mode coupling, while for $k_{z}\geq k_{z}^{(AC)}$ it causes the avoided crossing. As $g_{x}$ becomes greater, $k_{z}^{\rm(AC)}$ steeply gets smaller than the characteristic $k_{z}=\tilde{k}_{z}$. Examination of eq. (65) for this case in the limit $g_{x}\gg g_{z}$, which corresponds to $k_{x}\gg k_{z}$ for type-I crossing, leads to conclusion that the growth rate of SI decreases like

[TABLE]

where $\Im[\omega]_{\rm stl}=\kappa\,(f/2)^{1/2}$ is the maximum growth rate given by eq. (58).

In the absence of the dust settling, i.e. for $g_{z}=0$, there is no coupling for type-I crossing, and one arrives at the following estimation for the real corrections to the frequencies of SDW and IW:

[TABLE]

which is in accordance with eq. (5.10) of SH to the leading order in $\tau$.

The transition of type-I crossing from the mode coupling to the avoided crossing as the radial drift of the dust becomes comparable to its settling is shown in Fig. 3. For several particular solutions of eq. (49) the components of gravitational acceleration are chosen in such a way that the mode crossing occurs at the same point as in Fig. 1. It is seen that as $g_{x}$ increases, the band of SI shrinks and the maximum growth rate becomes smaller. Next, for $g_{x}/g_{z}\approx 3/4$, which is determined by eqs. (66-67) for this mode crossing, the mode coupling ceases and the two separate curves of $\Re[\omega]$ emerge instead. Each of these curves consists of the parts, which can be associated with the modes akin to IW and SDW as demonstrated in Fig. 3. In other words, the branches representing the modes akin to IW and SDW become broken as the avoided crossing occurs. The corresponding analytical estimations done with the account of both the radial drift and the vertical settling are in a good agreement with an accurate result. Turning to the profiles of the energy of modes, one finds that as far as SI exists, they are similar those shown in Fig. 1. However, in the vicinity of the critical point defined by eqs. (66-67) there is another situation. Indeed, in the case marked as ‘3’ in Fig. 3 the mode coupling still occurs, though it is rather weak. Nevertheless, one of the corresponding branches of mode akin to SDW acquires the positive energy. Further, as the avoided crossing arises, one finds that the two separate curves of $\Re[\omega]$ marked as ‘4’ represent the positive energy modes.

3.10 Behaviour of mode crossings at negative radial wavenumber

All crossings except the type-I crossing take place at $k_{x}<0$. Eq. (65), or equivalently, the energy density of SDW (47) along with the rule noted in Section 3.5, indicate that type-II and type-III(′) crossings result in the mode coupling for any $f>0$, while type-IV crossing always provides the real correction to the frequency, thus, resulting in the avoided crossing. As follows from table 2, type-II and type-III crossings occur for the substantial settling only, as they exist for $g_{x}<g_{x}^{{}^{\prime}}$ and $g_{x}<g_{x}^{{}^{\prime\prime}}$, respectively. At the same time, type-III′ and type-IV crossings remain in the limit $g_{z}\to 0$. Additionally, unlike the type-II crossing, type-III,III′ and -IV crossings are not confined over $k_{z}$ and $|k_{x}|$.

In this Section the solution of eq. (49) is shown for the case when settling dominates the radial drift of the dust, see Fig. 4 for type-II crossing and Fig. 5 for type-III and type-IV crossings. The components of gravitational acceleration are chosen in such a way that the corresponding type-I crossing occurs at the same point as in Fig. 1. The mode coupling at the type-II crossing is similar to the case of dust settling, see Fig. 1, with the difference that band of instability is slightly wider at the type-II crossing for the same value of the dust fraction. In contrast, type-III crossing provides the band of instability in much greater range of $k_{x}$ with the growth rate by more than two times higher for the same value of the dust fraction. The bottom panel in Fig. 5 shows that in the vicinity of the mode coupling the energy of modes akin to IW and SDW still have the opposite signs, however, as soon as $\Re[\omega]$ of the mode akin to SDW becomes positive, its energy also acquires the positive values. Thus, the type-IV crossing of SDW which occurs with the other IW, provides the avoided crossing with no SI. As a result, the loop of the energy profile forms to the left from the mode coupling in Fig. 5. The analytical estimates made at type-III and type-IV crossings show less accordance with an accurate solution which may be caused by the proximity of the crossings to each other. It can be checked that for greater value of the dust fraction, or alternatively, for higher $k_{z}$ than shown in Fig. 5 the mode coupling expands beyond the type-IV crossing, while at the same moment the ordinary coupling of two modes is replaced by the novel variant of coupling between three modes, i.e. two IW and one SDW, see Section 3.11.

In the opposite situation when the radial drift of the dust dominates its settling, i.e. $g_{z}\ll g_{x}$, type-III′ crossing provides the mode coupling in the high-$k_{z}$ limit, see Fig. 6. This case can be considered analytically employing the limits $g_{x}\gg g_{z}$ and $k_{z}\gg\tilde{k}_{z}$. Then the growth rate given by eq. (65) is reduced to the following result

[TABLE]

Eq. (70) confirms that in this case the growth rate does not depend on the rate of radial drift, see Fig. 6. Instead, the band of the instability shrinks over the radial wavenumber. Eq. (70) can be compared with the growth rate in case of the dust settling, see eq. (58),

[TABLE]

which is formally even higher than the maximum value found in the absence of radial drift. However, as $g_{z}\to 0$, $\tilde{k}_{z}\to\infty$ and this branch of SI should be suppressed by viscosity in real disc.

3.11 Triple coupling of modes

Here the mode coupling provided by the type-III crossing is considered in the high-wavenumber limit $k_{z}>\tilde{k}_{z}$ in the case when settling dominates radial drift, see Fig. 7. As the dust fraction increases, the growth rate of SI behaves in a different way than it is expected from the ordinary mode coupling. In the high-$k$ limit the resonance condition (59) gives that $k_{x}\to-(g_{z}/g_{x})k_{z}$ and eq. (65) yields

[TABLE]

As it is seen on the bottom panel in Fig. 7, eq. (71) along with the more general eq. (65) provides satisfactory approximation of an accurate solution of eq. (49) for smaller $f=0.001$ only. As the dust fraction becomes higher, those equations overestimate actual growth rate evaluated at the type-III crossing. The reason for such discrepancy can be found out deriving the approximate solution to eq. (49) in a slightly different way than it is done in Section 3.5. Indeed, eq. (49) can be written in another form:

[TABLE]

Again, let the deviation from the type-III crossing caused by the coupling term be $\Delta=\omega-\omega_{i}=\omega-\omega_{p}$, where $\omega_{c}=\omega_{i}=\omega_{p}<0$. As far as

[TABLE]

$\Delta$ is determined by equation

[TABLE]

which is equivalent to eq. (55).

The condition (73) implies that the correction to $\omega_{c}$ with the account of coupling is much less than the difference between the frequencies of oppositely propagating IW. However, as $|k_{x}|,k_{z}\to\infty$, this difference tends to constant value $\sim\kappa g_{x}/g_{z}$, while the coupling term $|\epsilon|\sim O(k)\to\infty$ as far as the dust fraction is constant. In the marginal case

[TABLE]

Employing eq. (59) along with eq. (71), one obtains from eq. (75) the upper value of $f$ for the ordinary mode coupling

[TABLE]

As $f>f_{cr}$, it does not occur any more. In the opposite case $|\Delta|\gg|2\omega_{i}|$, $\Delta$ is determined by equation

[TABLE]

which yields the following result for the growth rate of the high-$k$ SI for $f>f_{cr}$:

[TABLE]

which recovers an equation (5.16) of SH. Fig. 7 demonstrates that imaginary part of eq. (78) is in much better agreement with an accurate solution rather than eq. (71). Note that for the case shown in fig. 7 $f_{cr}\approx 0.006$.

The transition from eq. (74) to eq. (77) implies than coupling of two modes is replaced by coupling of three modes, which are two IW and SDW. This new mechanism of instability, which can be referred to as the triple mode coupling, is responsible for different dependence of RDI growth rate on the dust fraction, i.e. $\Im[\omega]\sim O(f^{1/3})$ instead of the more common $\Im[\omega]\sim O(f^{1/2})$. Eq. (76) shows that the triple mode coupling is preferential for high wavenumbers in situations with the substantial settling of dust. Additionally, it occurs for particles of a greater size, as $f_{cr}\sim t_{s}^{-1}$. The weaker dependence of the growth rate on the small dust fraction can make the triple mode coupling the leading mechanism for dust clumping in certain regimes, which is advisable to check in the future work.

3.12 Bonding of inertial waves

It can be seen that in the case shown in fig. 7 the band of instability for $f=0.01$ extends unusually far towards $k_{x}\to-\infty$. It is shown below that this occurs due to the new instability of quasi-resonant nature.

The zone of lower negative $k_{x}$ is shown in Fig. 8 for the parameters used before in Fig. 7. As one can see, a band of new instability appears to the left from the band of SI produced by the type-III crossing. The top panel of fig. 8 shows that as this new instability sets on, the phase velocity of both modes akin to IW vanishes. In that sense it resembles the mode coupling, however, the other features are different:

i) the mode crossing between two IW is located formally at infinity, $k_{x}\to-\infty$;

ii) both IW have positive energy at $k_{x}\to-\infty$; also, it can be checked that, as $f>0$, both modes akin to IW have positive energy outside of the band of new instability.

By these reasons, it is preferable to use a different term for this mechanism of instability. Hereafter it is referred to as ’bonding of IW’.

Thus, as $k_{x}\to-\infty$, two IW propagating in the opposite directions become the low-frequency modes. The difference between their phase velocities becomes small, which makes possible the bonding of IW caused by the coupling term, i.e. induced physically by the drag force acting on gas. Analytical estimation for the growth rate of bonding instability of IW can be obtained along the following lines. Similar to the case with the triple mode coupling, the starting point is eq. (72). Now, the substitution is $\omega=\omega_{b}+\Delta$, where $\omega_{b}=0$ is not the frequency of the mode crossing as in the case of the mode coupling, but rather the average of the opposite frequencies of two IW taken in the case $f=0$. As far as $|\Delta|\ll|\omega_{p}|$, this yields

[TABLE]

Clearly, for the range of $k_{x}<0$ lower that the type-IV crossing both the coupling term $\epsilon$ and the frequency of SDW $\omega_{p}$ are positive, which implies that the bonding instability sets on as soon as

[TABLE]

The equality in the instability condition (80) yields the following approximate location of the bound of bonding instability along the radial wavenumbers:

[TABLE]

provided that there is a sufficient dust radial drift:

[TABLE]

Note that $|k_{x}|\gg k_{z}$ in this limit.

As far as the absolute value of $k_{x}/k_{z}$ is much higher than the corresponding absolute value of (81), the estimation of the growth rate is

[TABLE]

As for RDI, the growth rate of bonding instability $\propto f^{1/2}$, which is a feature of resonant mechanism. However, here it is referred to as ’quasi-resonant’ by the reasons mentioned above in this Section.

It was checked that eq. (83) gives a slight overestimation of an accurate value, nevertheless, the both values approach each other while $|k_{x}|\to\infty$. The frequency of corresponding gas-dust mode is imaginary, which means that bonding instability is provided by standing exponentially growing perturbations.

By virtue of eq. (81) one can put the upper limit on the growth rate of the bonding instability as

[TABLE]

which shows that it scales with the first power of the dust fraction, but still is comparable to SI growth rate in situation when dust settling dominates the radial drift of the dust. Eq. (84) implies that in geometrically thin disc the latter occurs as long as $h/r_{0}\lesssim f^{1/2}$, which is commonly the case. It is also important that, similar to SI, bonding instability does not depend on $t_{s}$.

4 RDI beyond TVA

The analysis performed above shows that there is no RDI555as well as the bonding instability in the absence of settling, $g_{z}\to 0$, within TVA. The solutions of eq. (49) become symmetric with respect to change $k_{x}\to-k_{x}$, so that both type-I and type-IV crossings provide the identical avoided crossing, see eq. (69). Hereafter it is assumed that $k_{x}>0$. Thus, in order to recover RDI obtained by SH in the case of the dust radial drift with no settling, the next order terms over the small $\tau_{*}$ and $\lambda^{-1}$ must be retained. Note that this concerns not only the terms $\sim UV$, but also the smallest terms $\sim V^{2}$ by the following reasoning. Since the condition of the mode crossing is now

[TABLE]

by the order of magnitude

[TABLE]

in the vicinity of resonance, which implies that $\lambda^{-1}\sim\tau_{*}$, i.e. the terms $\sim V^{2}$ in eq. (11-12) are comparable to the inertial terms in eq. (12).

4.1 Stationary solution

Given the vertical gravity be negligible, $g_{z}=0$, the general set of equations (11-14) have the following stationary solution

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the terms $\sim\tau^{3}$ and smaller have been omitted.

4.2 Reduced dispersion equation

Here the strategy is to rederive the dispersion equation from the general set of equations for axisymmetric perturbations relaxing TVA, though adding the corresponding higher-order terms in $\tau$ in the coupling term only. For this reason, it is referred to as the reduced dispersion equation. The higher-order corrections either to $D_{g}$ or to $D_{p}$ are not considered as those are not relevant to RDI, modifying the non-resonant dynamics of modes akin to IW or SDW and introducing (probably imaginary) corrections to the frequencies $\omega_{i}$ and $\omega_{p}$, which are $\sim f$, rather than $\sim f^{1/2}$ as for the RDI growth rate.

The set of full equations (11-14) yields the following equations for modes of linear gas-dust perturbations provided that the stationary solution is given by eqs. (86-89)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where it is assumed that $f\ll 1$, the terms $\sim f^{2}$ have been omitted and

[TABLE]

[TABLE]

The further derivation proceeds along the lines which lead to eqs. (31-34). That is, the coupling term is a product of combinations of terms containing $\hat{\delta}$ in eqs. (90-92) and the term containing $\hat{u}_{y}$ in eq. (97), after the latter has been rearranged with the help of eqs. (93-96).

First, the components of perturbation of the relative velocity can be represented in the following way

[TABLE]

[TABLE]

where it is assumed that

[TABLE]

and the higher order terms $\sim\tau^{3}$ have been omitted. Therefore, for purpose of this Section the term containing ${\bf v}$ in eq. (98) can be omitted hereafter. Eqs. (100-101) are substituted into eq. (96), so the second terms in square brackets of eqs. (100-101) yield together a single term $\sim O(f)$ which is also omitted in the new equation for $\hat{\delta}$.

At last, the set of eqs. (90-97) is rearranged as the following

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Eqs. (103-106) lead to the reduced dispersion equation

[TABLE]

where

[TABLE]

and it is assumed that $\omega_{p}$ is given by equation (102). Clearly, $\tilde{\epsilon}\to\epsilon$ in the leading order in $\tau$.

In turns out that in the vicinity of resonance (85) the solution of equation (108) is in a good agreement with the solution of the full dispersion equation derived from equations (90-97), see Fig. 1 of SH.

The approximate solution of eq. (49) at the mode crossing of IW and SDW given by eq. (55) with the replacement $\epsilon\to\tilde{\epsilon}|_{\omega_{c}}$ applied to equation (107) yields

[TABLE]

where the resonance condition (85) is used. Eq. (109) recovers the result of SH, see their equation (5.10).

Consideration of RDI of radially streaming dust reveals that the underlying mechanism goes beyond the concept of the mode coupling. In this situation, RDI is caused by a quasi-resonant process manifested in the growth/damping of the uncoupled modes both having positive energy and slightly different phase velocities at the mode crossing. Moreover, as the dust fraction increases, the growth rate increases along with the degree of the decoupling between IW and SDW, see eq. (109). At the same time, such a mechanism is generally weaker than the mode coupling as it emerges due to inertia of solids represented by the next-order terms in $\tau$ entering an augmented coupling term (108).

5 Conclusions

This study is focused on the nature of local instability of gas-dust mixture with dust streaming in protoplanetary disc. Previously, SH revealed that in the leading order in small $f$ axisymmetric perturbations grow due to RDI, which provides the growth rate $\propto f^{1/2}$. Here RDI is reconsidered in terms of perturbations of the center-of-mass velocity and the relative velocity of gas-dust mixture as well as perturbation of dust density. This allows for the analysis of RDI within TVA. The dynamics of modes of perturbations is considered in the general case of dust settling combined with dust radial drift see eqs. (31-34). The usefull framework of the Lagrangian perturbation theory for neutral modes of perturbations is applied to two-fluid model in order to show the existence of negative energy waves in such a flow from fundamental symmetry of the system with respect to translations in time, see eq. (48). Generally, the negative energy wave is akin to SDW, while it becomes identical to SDW in the absence of dust back reaction on gas, which is reproduced by the formal limit $f\to 0$. The energy of SDW is given by eq. (47). In the particular case of the dust settling without radial drift, there is a straightforward generalisation of variational principle onto arbitrary linear gas-dust perturbations, which makes it possible to show that the energy of growing/damping modes vanishes, see the Appendix B.

The dispersion equation (49) for gas-dust perturbations with the account for dust back reaction on gas is reduced to form typical for systems with the mode coupling. In the absence of the coupling term the dispersion equation splits into two independent equations describing one SDW and two oppositely propagating IW, respectively. As the dispersion curve of SDW crosses the dispersion curve of IW, the modes fall in resonance, and the corresponding point of phase space is referred to as the mode crossing. There can be up to four mode crossings for the given free parameters of the model, see Fig. 2. The formal classification of mode crossings is suggested in table 1. Their existence in phase space is outlined in table 2. Type-I crossing occurs for any combination of settling and radial drift, though in the limited range of wavenumbers. Type-II crossing occurs in the limited range of wavenumbers, but for the substantial settling only. Type-III crossing occurs for the substantial settling only, but at sufficiently high wavenumbers. Type-III′ crossing replaces type-III crossing as radial drift of the dust dominates its settling. At last, type-IV crossing occurs for all wavenumbers each time there is a radial drift of the dust. In the marginal case of dust streaming only vertically there are type-I,II crossings identical to each other, whereas, in the marginal case of dust streaming only radially there are type-I,IV crossings identical to each other.

The definite criterion for the existence of SI, i.e. RDI in the presence of the dust settling, at the given mode crossing is formulated in Section 3.5: the energy of SDW evaluated at the mode crossing must be negative. In the opposite case of the positive energy of SDW, the avoided crossing occurs, which gives no instability within TVA. The avoided crossing takes place for type-IV crossing and also for type-I crossing in the case of a sufficient radial drift of solids, see the conditions (66-67). In other situations the mode crossing provides SI. SI occurs due to the energy transfer from (negative energy) SDW to (positive energy) IW in the vicinity of the corresponding mode crossing. The coupled modes are represented by the complex conjugate pair of frequencies obtained from the accurate solution of the dispersion equation, see Figs. (1,3-6). The band of instability expands in the space of wavevectors as the coupling term in the dispersion equation, which is proportional to the dust mass fraction, becomes larger. SI attains the maximum growth rate approximately at the mode crossing. It is estimated analytically employing the usual mode coupling approach, see eq. (55). Eq. (55) also indicates that within TVA the correction to the mode frequency due to the dust back reaction on gas is generally independent on $t_{s}$. The latter applies to both the mode coupling and the avoided crossing.

The branch of SI having an unbounded asymptotics of growth rate at high wavenumbers is provided by mode coupling at type-III(′) crossing. The corresponding analytical estimate of maximum growth rate is given by eq. (71). Note that it is determined solely by the rate of vertical settling $\sim t_{s}g_{z}$ despite the presence of radial drift necessary for the very existence of the considered mode coupling. Also, this branch of SI takes place for wavenumbers of different signs, i.e. either for $k_{z}>0$ and $k_{x}<0$, or for $k_{z}<0$ and $k_{x}>0$. Further, as far as $k$ increases for constant dust fraction, the distance between type-III and type-IV crossings in phase space tends to a constant value, whereas the absolute value of the coupling term grows up. Eventually, the ordinary coupling between SDW and one of IW is replaced by the triple coupling between SDW and the both of IW, see Section 3.11 for details. As a result, the growth rate acquires a different dependence on the dust fraction, $\propto f^{1/2\to 1/3}$. The triple mode coupling yields eq. (78) for estimation of the maximum growth rate, which recovers the result of SH. Note that unlike eq. (71) it explicitly depends on both components of gravitational acceleration.

For the absolute value of radial wavenumber even higher than that of type-III,IV crossings, an additional branch of instability emerges, see Section 3.12 for details. For the sufficiently high dust fraction it is attached to the band of instability provided by the mode coupling at type-III crossing, see the curve marked ‘3’ on the bottom panel of Fig. 7. As is shown in Fig. 8, the new instability appears in the region where two low-frequency oppositely propagating IW coalesce with each other, while their phase velocity vanishes. As in the case of the mode coupling, the dust back reaction on gas, which introduces the coupling term of the dispersion equation, is responsible for this new process. However, the underlying mechanism of growth is different, since the coalesced modes have positive energy and there is no resonance between the modes in the limit $f\to 0$. Nevertheless, the low-frequency IW approach each other in phase space as $|k_{x}|\to\infty$, while the growth rate of new instability $\propto f^{1/2}$, see eq. (83). By these reasons, in this study it is referred to as quasi-resonant instability of bonding IW. The upper limit of the bonding instability estimated by eq. (84) confirms that this instability is important along with SI in the situation when the settling of the dust dominates its radial drift.

The short wavelengths typical for both the triple mode coupling and the bonding instability are susceptible to dissipation effects in a disc. The influence of effective viscosity of gas on these branches of instability should be considered in the future work in order to understand their relevance for structure formation in real protoplanetary discs. Besides, moderate dissipation may not necessarily lead to damping of high-wavenumber SI. In contrast, it may cause an additional non-resonant growth of negative energy mode akin to SDW. This issue will be addressed in the subsequent studies.

In the limit of the dust streaming only radially, $g_{x}\gg g_{z}$, or equivalently $g_{z}\to 0$, SI ceases to operate. The growth rate due to the mode coupling produced by type-I crossing vanishes according to eq. (68), while the mode coupling produced by type-III(′) crossing shifts towards $k\to\infty$ to be suppressed by any small amount of viscosity, see eqs. (70), (71) and (78). Strictly for $g_{z}=0$, there exist type-I and type-IV crossings identical to each other, which produce the avoided crossings with the account for dust back reaction on gas. Thus, there is no RDI within TVA when dust is subject only to radial drift. RDI is recovered by retaining the terms of the orders of $\tau_{*}$, $\sqrt{\lambda^{-1}}$ and $\lambda^{-1}$ from the full equations (11-12), which contribute to the coupling term of the dispersion equation, see eqs. (107-108). It turns out that in the vicinity of resonance between SDW and IW all those terms yield the imaginary correction $\sim O(\tau)$ to the coupling term. Therefore, the two separate branches of neutral modes representing the avoided crossing acquire small growing/damping. In this way, it becomes clear that the streaming instability of Youdin & Goodman (2005) is provided by resonant mechanism other than the mode coupling responsible for SI. This new mechanism remains obscured. Though, it should be related to the inertia of solids, because the latter is exactly what is neglected in TVA. Also, its relationship with the inertia of solids naturally explains the increase of the growth rate with the size of the particles.

The dispersion equation obtained beyond TVA, see eq. (107), is a reduced one, since the corrections $\sim O(f)$ to the dispersion relations of IW and SDW given, respectively, by $D_{g}$ and $D_{p}$, have been omitted there. However, as shows the analytical study of Jacquet et al. (2011), the modes can exhibit an intrinsic growth due to non-resonant mechanisms. This is especially probable for negative energy SDW, what remains to be checked further. The non-resonant growth may be important as $f\to 1$.

The relevance of the streaming instability to planetesimal formation has being extensively studied via local simulations of the non-linear dynamics of perturbations in dust-laden disc midplane, see the recent results by Yang et al. (2017). Although the numerical models usually incorporate the dust settling, SI has not been found so far, see the discussion by Squire & Hopkins (2018a). Besides, SI might be less sensitive to external turbulence inherent in astrophysical discs, see the recent numerical study of the dust settling through turbulent gas by Lin (2019). At last, it is highly likely that the other RDI revealed by Squire & Hopkins (2018a) are caused by coupling of SDW with waves of the other origin, such as sound waves and internal gravity waves ubiquitous in protoplanetary discs. These issues are relegated to future work.

Acknowledgments

The author acknowledges the support from the Program of development of M.V. Lomonosov Moscow State University (Leading Scientific School ’Physics of stars, relativistic objects and galaxies’).

Appendix A Equations for perturbations in TVA

The perturbed state of gas-dust mixture can be described by the small Eulerian perturbations of the centre-of-mass velocity, ${\bf u}$, the relative velocity, ${\bf v}$, the gas pressure, $p^{\prime}$, and the density of dust, $\rho_{p}^{\prime}$. Those quantities obey the corresponding linear equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The set of eqs. (110-113) describes the evolution of gas-dust perturbations for any stationary solution of eqs. (18), (19), (13) and (14) on the length-scales much shorter than the scaleheight of a thin protoplanetary disc. Eqs. (110-113) become much more simple for the particular background model considered in Section 3.1.

Once the background model is specified by eqs. (22-25), one arrives at the following equations for perturbations

[TABLE]

[TABLE]

[TABLE]

where $W\equiv p^{\prime}/\rho$ and the relative perturbation of the dust density $\delta\equiv\rho_{p}^{\prime}/\rho_{p}$. Note that eqs. (114-116) are valid for $f$ not necessarily of the small value.

A.1 Equations for axisymmetric perturbations in the leading order in small dust fraction

Let additionally the dust fraction be small $f\ll 1$ and perturbations have axial symmetry. Eqs. (114-116) yield

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Hence, in this the most simple case the small gas-dust perturbations of dust streaming through the static gas environment both vertically and radially are described by eqs. (117-119, 121) for perturbation of the centre-of-mass velocity, ${\bf u}$. These equations are identical to equations for single-fluid dynamics of vortical perturbations in a rotating plane shear flow, see e.g. equations (6-9) of Umurhan & Regev (2004) with $\partial_{y}=0$, except for the terms in the RHS of eqs. (117) and (119) proportional to perturbation of the dust density. At the same time, perturbation of the dust density is the only quantity that additionally describes the dynamics of dust by means of eq. (120). The perturbed flow of dust is affected by gas through the term $\propto\nabla^{2}W$ in the RHS of eq. (120). In the limit $f\to 0$ the set of eqs. (117-121) splits into eqs. (117-119, 121), which determine the dynamics of solely the gas perturbations with ${\bf u}$ becoming identical to the velocity of gas, and eq. (120), which separately determines the dynamics of dust perturbations. In the latter situation, solids move passively under the action of external gravity and aerodynamic drag. If additionally $t_{s}\to 0$, both the background relative velocity and the perturbation of relative velocity vanish, which results in the conservation of the density of dust frozen in the divergence-free motion of gas.

The main results of this paper come out from eqs. (117-121).

A.2 Equations for perturbations in the appropriate variables

It is appropriate to reformulate eqs. (117-121) in terms of the variables

[TABLE]

It is assumed that components of $\chi$ have non-zero derivatives over $t,x,z$, which are denoted below as $\partial_{k}\chi_{i}$. The usual Einstein’s rule of summation over the repeated upper and lower indices is assumed as well.

First, taking the curl of eqs. (117-119) yields the following

[TABLE]

[TABLE]

[TABLE]

Note that eqs. (123-125) are identical to equations describing IW in rigidly rotating fluid provided that $q=f=0$, see e.g. Landau & Lifshitz (1987), paragraph 14. The LHS of eqs. (123-125) are the time derivatives of components of vorticity perturbation666In this work vorticity is a curl of the center-of-mass velocity of gas-dust mixture rather than a curl of the velocity of gas.

Next, the divergence of eqs. (117-119) yields

[TABLE]

Taking additional derivatives over $z$ from equations (120), (126), (125) and employing eq. (121) one comes to the final set of equations for the new variables

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the term $\propto ft_{s}$ has been neglected in eq. (130).

Appendix B Variational principle for dynamics of general gas-dust perturbations in the case of the dust vertical settling

The set of eqs. (127-130) for $g_{x}=0$ can be derived from the requirement that action

[TABLE]

with the Lagrangian density

[TABLE]

be stationary with respect to arbitrary variations of $\chi_{i}$ in the small shearing box. Similar to the variational principle for modes, see eq. (36), the subscript ‘0’ denotes the contribution to the full Lagrangian responsible for the dynamics of gas-dust perturbations with no account for the drag force acting on gas, while the subscript ‘1’ denotes additional contribution to ${\cal L}$ responsible solely for the dynamics of dust. If so, eqs. (127-130) are identical to the Euler-Lagrange equations

[TABLE]

where explicitly

[TABLE]

[TABLE]

The first term in ${\cal L}_{1}$ is the generalized cross-term, cf. eq. (39). Like the Lagrangian (37), $\cal{L}$ is known up to an arbitrary constant factor chosen so that the energy of IW propagating in rigidly rotating fluid be positive.

B.1 Energy of general gas-dust perturbations

The energy of general perturbations associated with the invariance of ${\cal L}$ with respect to translations in time is $E\equiv\int{\cal E}d^{2}{\bf x}$, where

[TABLE]

is the energy density. Similar to the energy density of mode (45), it consists of the main part recovering the energy density of general gas perturbations in rotating flow in the absence of dust, and additional contribution related to general perturbations of dust density

[TABLE]

Explicitly,

[TABLE]

[TABLE]

It is straightforward to check that eqs. (127-130) guarantee that $\partial_{t}{\cal E}$ can be represented as the divergence of vector, which is identified with the energy flux vanishing along with perturbations $\chi_{i}$ as one goes to infinity. ${\cal E}_{0}$ is of course a positive definite quantity in centrifugally stable flows, while ${\cal E}_{1}$ can have any sign depending on the contribution of the cross term. However, the net contribution of the cross term to $E$ can be small or even vanish for certain non-trivial perturbations, i.e. for modes, see eq. (45). Then, owing to perturbations of dust density, ${\cal E}_{1}$ gives negative amount to the energy for the case $g_{x}=0$ considered here.

B.2 Energy of growing mode

In order to obtain expression for the energy density of mode, which is not necessarily neutral, $\Re[\chi_{i}]$ is substituted into eq. (135) and subsequently averaged over the wavelength of mode. Eqs. (31-34) provide the following result

[TABLE]

where $k^{2}\equiv k_{x}^{2}+k_{z}^{2}$.

In this way, one obtains the averaged energy density of the mode of gas-dust mixture perturbations for a particular wavevector ${\bf k}$, provided that its (complex) frequency is known, i.e. there is a solution of a dispersion equation.

B.2.1 Approximate solution at the mode crossing and the vanish of the mode energy

Eq. (49) is cubic with respect to $\omega$. Seeking the Cardano solution of eq. (49) for $g_{x}=0$ at the mode crossing (57) one finds that the discriminant

[TABLE]

vanishes for $f=0$, and eq. (49) acquires a double root. Then, any small non-zero $f>0$ makes $Q$ be a positive value, which implies that a double root of eq. (49) acquires the adjoint couple of imaginary parts, i.e. the instability sets on in some band around the mode crossing.

In the case of small non-zero $f>0$ the Taylor expansion of the Cardano solution of eq. (49) by the orders of $f^{1/2}$ at the mode crossing yields the following approximate real and imaginary parts of the frequency

[TABLE]

[TABLE]

Note that eq. (140) recovers eq. (58) in the leading order in $f^{1/2}$. It is straightforward to check that expressions (139) and (140) make $\hat{\cal E}$ as given by eq. (138) vanish for any small $f>0$.