Interpretation of the resonant drag instability of dust settling in protoplanetary disc

TL;DR

This paper analyzes the resonant drag instability in protoplanetary discs, showing how dust settling and mode coupling lead to dust clumping, with growth rates depending on dust velocity and mixture parameters.

Contribution

It introduces a mode coupling framework between inertial waves and settling dust waves, explaining the instability mechanism in dust-gas mixtures.

Findings

Instability growth rate scales with the square root of the coupling term.

Higher dust settling velocities lead to faster dust clumping.

Mode coupling enables dust aggregation in protoplanetary discs.

Abstract

The recently discovered resonant drag instability of dust settling 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: the first mode is inertial wave having positive energy, while the second mode is a settling 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 settling of dust, which transports perturbations of dust density. The phase velocity of SDW is equal to the bulk settling velocity times the cosine of the angle formed by the wave vector and the rotation axis. In this way, the bulk settling of dust makes possible the coupling of SDW with the inertial wave and the onset of the instability. In accordance with the concept of the mode coupling, the instability growth rate is proportional to the square root of the dispersion equation coupling term, which itself contains the small mass fraction of dust in gas-dust mixture, the squared radial wavenumber of the modes, and the squared bulk settling velocity. Thus, the higher is the bulk settling velocity, the heavier clumps of dust can be aggregated by the instability of the same rate.

Equations

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

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

$$\frac{1}{{\rho }_b}\frac{\partial p_b}{\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 }_b}+\frac{\sigma }{{\rho }_b}\frac{\mathbf{V}}{\tau },$$

$$\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_b}{{\rho }_b}-\frac{\mathbf{V}}{\tau },$$

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

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

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

$${\mathbf{U}}_p=\mathbf{U}+\frac{{\rho }_b}{\hat{\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 }_b}{\hat{\rho }}\left\{\left(\mathbf{V}\nabla \left(\frac{\sigma }{\hat{\rho }}\right)\right)\mathbf{V}+2\frac{\sigma }{\hat{\rho }}\left(\mathbf{V}\nabla \right)\mathbf{V}\right\}=\frac{\nabla p_b}{{\rho }_b}-\frac{\nabla (p+p_b)}{\hat{\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 }_b}{\hat{\rho }}(\mathbf{V}\nabla )\left(\frac{{\rho }_b}{\hat{\rho }}\mathbf{V}\right)-$$

$$\frac{\sigma }{\hat{\rho }}(\mathbf{V}\nabla )\left(\frac{\sigma }{\hat{\rho }}\mathbf{V}\right)=\frac{\nabla (p+p_b)}{{\rho }_b}-\frac{\hat{\rho }}{{\rho }_b}\frac{\mathbf{V}}{\tau },$$

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

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

$$f\equiv \frac{\sigma }{{\rho }_b}<1.$$

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

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

$${\lambda }_{\ast }^{-1}\equiv \frac{g{\tau }^2}{\lambda }\ll 1,$$

$$V\sim g\tau ,U\sim \frac{g\tau }{\sqrt{{\lambda }_{\ast }^{-1}}},$$

$$V/U\sim \sqrt{{\lambda }_{\ast }^{-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_b}{{\rho }_b}-\frac{\nabla (p+p_b)}{\hat{\rho }},$$

$$\frac{\nabla (p+p_b)}{\hat{\rho }}=\frac{\mathbf{V}}{\tau }.$$

$$({\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}+(\mathbf{U}\nabla )\mathbf{u}=-\frac{\nabla p^{\prime }}{\hat{\rho }}+\frac{\nabla (p+p_b)}{\hat{\rho }}\frac{{\sigma }^{\prime }}{\hat{\rho }},$$

$$\frac{\nabla p^{\prime }}{\hat{\rho }}-\frac{\nabla (p+p_b)}{\hat{\rho }}\frac{{\sigma }^{\prime }}{\hat{\rho }}=\frac{\mathbf{v}}{\tau },$$

$$\nabla \cdot \left(\mathbf{u}-\frac{\sigma }{\hat{\rho }}\mathbf{v}-\frac{{\rho }_b}{\hat{\rho }}\frac{{\sigma }^{\prime }}{\hat{\rho }}\mathbf{V}\right)=0,$$

$$({\partial }_t-q{\Omega }_{0}x{\partial }_y){\sigma }^{\prime }+\nabla ({\sigma }^{\prime }\mathbf{U}+\hat{\rho }\mathbf{u})=0.$$

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

Full text

††pagerange: Interpretation of the resonant drag instability of dust settling 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. 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 a new vigorous instability found by Squire & Hopkins (2018a) who termed it as settling instability, SI hereafter. Along with the streaming instability, SI is another manifestation of the resonant drag instability 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 Squire & Hopkins (2018a).

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 stream of the settling dust, the electromagnetic wave in helical coil by inertial wave (IW) in gas environment, the space-charge wave by a settling 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) 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 general two-fluid perturbations is introduced for the particular case of dust settling to the midplane of protoplanetary disc. It is found that the energy of certain (neutral) mode of such perturbations akin to SDW is negative. SDW itself is a trivial negative energy mode, which 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 settling of dust that gives rise to a reservoir of negative energy in the form of SDW, from which the 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 relation, which is proportional to the mass fraction of dust in gas-dust mixture.

The study starts from general equations for local dynamics of gas-dust mixture in protoplanetary disc, step by step proceeding to a particular model of the perturbed gas-dust mixture with dust bulk settling to the disc midplane. Further, a variational principle for evolution of subsonic gas-dust perturbations is formulated and energy of perturbations is derived. These results are applied to modes of perturbations in order to demonstrate that SI is caused by the mode coupling between SDW and IW.

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}_{b}$. The latter describes the stationary circular 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_{b}$ and $\rho_{b}$ 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_{b}$, 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 $\sigma$. The dust particles stopping time $\tau=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 $\sigma=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\sigma$ 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 $\hat{\rho}\equiv\rho_{b}+\sigma$ is the total density of gas-dust mixture.

Since

[TABLE]

[TABLE]

one arrives at the new equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

2.1 Terminal velocity approximation

Let the characteristic time- and length-scales of gas-dust mixture dynamics be $t_{ev}$ and $\lambda$, 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]

Hereafter, 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 as terminal velocity approximation, 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, terminal velocity approximation is the case for 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 g\tau^{2}$.

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

[TABLE]

[TABLE]

This work deals with the analysis of the particular linear solution of the set of eqs. (18), (19), (13) and (14).

2.2 Equations for perturbations

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, $\sigma^{\prime}$. Those quantities obey the corresponding linear equations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The set of eqs. (20-23) 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. (20-23) become much more simple for the particular background model considered in the next Section.

3 Particular case of the dust vertical settling

In geometrically thin disc

[TABLE]

which may be considered as a constant value inside the small box. Note that in the last equality in eq. (24) 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}$. It is assumed that the box is located sufficiently far above the disc midplane for $g_{z}$ to dominate the radial effective gravity, see eq. (2). The corresponding restriction reads

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

on the local scale considered here.

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

[TABLE]

[TABLE]

[TABLE]

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

3.1 Equations for axisymmetric perturbations in the leading order in small dust fraction

Let additionally constrain the dust fraction $f\ll 1$ and consider axisymmetric perturbations only. Eqs. (30-32) yield

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Hence, in this the most simple case the small gas-dust perturbations of dust settling through the static gas environment to the disc midplane are described by eqs. (33-35, 37) 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 term of vertical balance in the RHS of eq. (35) 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. (36). The perturbed flow of dust is affected by gas through the term $\propto\nabla^{2}W$ in the RHS of eq. (36). In the limit $f\to 0$ the set of eqs. (33-37) splits into eqs. (33-35, 37), which determine the dynamics of solely the gas perturbations with ${\bf u}$ becoming identical to the velocity of gas, and eq. (36), 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 $\tau\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 part of this paper is focused on eqs. (33-37).

4 Variational principle for dynamics of small gas-dust perturbations

It is appropriate to reformulate eqs. (33-37) in terms of the variables

[TABLE]

where

[TABLE]

[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. (33-35) yields the following

[TABLE]

[TABLE]

[TABLE]

Note that eqs. (39-41) 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. (39-41) are the time derivatives of components of vorticity perturbation111In 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. (33-35) yields

[TABLE]

Taking additional derivatives over $z$ from equations (36), (42), (41) and employing eq. (37) one comes to the final set of equations for the new variables

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\kappa^{2}\equiv 2(2-q)\Omega_{0}^{2}$ is the epicyclic frequency squared and the term $\propto f\tau$ has been neglected in eq. (46).

The set of eqs. (43-46) can be derived from the requirement that action

[TABLE]

with $d^{3}{\bf x}\equiv dxdydz$ and the Lagrangian density

[TABLE]

be stationary with respect to arbitrary variations of $\chi_{i}$ in the small shearing box. 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, i.e. the dust back reaction on gas. In the absence of dust, ${\cal L}_{0}$ describes single-fluid vortical perturbations in rotating gas. The subscript ‘1’ denotes additional contribution to ${\cal L}$ responsible solely for the dynamics of dust in terms of the evolution of its density. The first term in ${\cal 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. If so, eqs. (43-46) are identical to the Euler-Lagrange equations

[TABLE]

where explicitly

[TABLE]

[TABLE]

Note that ${\cal L}$ is known up to an arbitrary (dimensional) constant factor since it describes the linear problem. The sign of ${\cal 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 the next Section.

4.1 Energy of gas-dust perturbations

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

[TABLE]

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

[TABLE]

Explicitly,

[TABLE]

[TABLE]

It is straightforward to check that eqs. (43-46) 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. Note that ${\cal E}_{0}$ is positive definite quantity222This is true in centrifugally stable flows only, while in general ${\cal E}_{1}$ can have any sign depending on the contribution of the cross term. Nevertheless, the net contribution of the cross term to $E$ can be small or even vanish for certain non-trivial perturbations. Then, owing to perturbations of dust density, ${\cal E}_{1}$ gives negative amount to the energy.

5 Modes of gas-dust perturbations

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

[TABLE]

That is, $\hat{\chi}_{i}$ are complex Fourier amplitudes of $\chi_{i}$. In eq. (54) ${\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 instability of dust settling to the disc midplane (SI).

Eqs. (43-46) yield

[TABLE]

[TABLE]

[TABLE]

[TABLE]

5.1 Energy of mode

In order to obtain expression for the energy density of mode, $\Re[\chi_{i}]$ is substituted into eq. (51) and subsequently averaged over the wavelength of mode. Eqs. (55-58) 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 frequency is known, i.e. there is a solution of a dispersion equation.

5.2 Dispersion equation

As it follows from eqs. (55-58), $\omega$ obeys equation

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

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

An accurate solution of eq. (60) for particular values of parameters is presented in Figs. 1, 2. 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 value of averaged energy density of mode, $\hat{\cal E}$, normalised by $|\varpi|^{2}$ is shown in Fig. 3.

As far as the dust fraction is negligible, $f=0$, there is no instability of gas-dust mixture. It is notable that 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. For the particular values of parameters taken in Figs. 1-3 this point corresponds to $k_{x}=6$. Further, inside the band of SI both modes acquire identical $\Re[\omega]$ (equivalently, phase velocities) and opposite $\Im[\omega]$, i.e. the solution of eq. (60) is given by complex conjugate pair of roots introducing damping and growing modes of gas-dust perturbations. At the same time, $\hat{\cal E}$ of such perturbations naturally vanishes, since this is a necessary condition for the energy of each mode to be conserved, see Friedman & Schutz (1978a). It is important that outside the band of SI, when $\omega$ is real, the modes have $\hat{\cal E}$ of different signs.

As far as $f=0$, eq. (60) splits into two separate dispersion equations

[TABLE]

and

[TABLE]

describing, respectively, two IW333Referred to as ‘epicyclic mode’ in Squire & Hopkins (2018a) propagating in opposite directions and one else mode, which is referred to as SDW. Indeed, eq. (64) can be obtained from eqs. (55-58) setting $\hat{\delta}=0$. Since the dust fraction is negligible, the remaining variables describe perturbations of gas only. The set of equations (55-58) 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, below the modes described by eq. (64) are referred as the IW444Note that according to eq. (21) such IW produce a non-zero perturbation of relative velocity.. On the contrary, eq. (65) follows from eqs. (55-58) 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 setting of dust with the velocity (28). The phase velocity of SDW is equal to settling velocity times the cosine of angle formed by the wave vector and the rotation axis.

It is important that one of the branches satisfying eq. (64) crosses the branch satisfying eq. (65) on the plane $(\omega,k)$ (or, alternatively, $(\omega,k_{x,z})$) as soon as the condition

[TABLE]

is true. Once (66) is satisfied, the phase velocities of IW and SDW are equal to each other, $\omega_{1}=\omega_{2}\equiv\omega_{c}=\omega_{p}$. Below eq. (66) is referred to as the condition of a mode crossing. Clearly, the IW participating in the mode crossing propagates downwards to the disc midplane. The second IW propagates upwards to the disc midplane and does not satisfy the condition (66) for any wavevector.

As one can see in Figs. 1,2, the non-negligible dust fraction causes the 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 Figs. 1, 2. 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, see Fig. 3. At the same time, the neutral waves coalesce as seen in Fig. 1 and give birth to a pair of growing and damping waves as seen in Fig. 2. Similar picture occurs in some applications from physics of hydrodynamical instability of shear flows, where it is regarded as a mode coupling, see Section 5.3.

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

[TABLE]

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

5.2.1 Approximate solution at the mode crossing

Eq. (60) is cubic with respect to $\omega$. Seeking the Cardano solution of eq. (60) at the mode crossing (66) one finds that the discriminant

[TABLE]

vanishes for $f=0$, and eq. (60) 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. (60) 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. (60) by the orders of $f^{1/2}$ at the mode crossing yields the following approximate real and imaginary parts of the frequency

[TABLE]

[TABLE]

It is straightforward to check that expressions (67) and (68) make $\hat{\cal E}$ as given by eq. (59) vanish for any small $f>0$. Also, eq. (68) is in accordance with eq. (5.13) of Squire & Hopkins (2018a) in the leading order in $f$ (their $\mu$). For small value of $f$ used in Fig. 2 the analytical increment/decrement (68) virtually coincides with an accurate solution of eq. (60).

5.3 The mode coupling

The picture revealed in Section 5.2 has been examined in the dynamics of waves in shear flows, see Fabrikant et al. (1998). Cairns (1979) suggested that the dispersion relation in a form of eq. (60) with small RHS can be interpreted as the coupling of modes having their own dispersion relations $D_{g}$ and $D_{p}$. The coupling causes the deviation of the mode frequencies by a small value $\Delta\equiv\omega-\omega_{c}$ at the mode crossing. Then,

[TABLE]

[TABLE]

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

[TABLE]

which immediately recovers eq. (68) in the leading order in $f^{1/2}$ after one employs the eqs. (61-63). Note that with the help of the mode crossing condition (66) the coupling term (63) can be reduced to the form

[TABLE]

so that the increment/decrement of SI becomes

[TABLE]

Eq. (70) can also be usefull in various estimations related to the problem. By virtue of eq. (66), the maximum growth rate of SI

[TABLE]

i.e. does not depend on the grain size. However, the larger the grains, the heavier the dust clumps, which can be aggregated by SI.

A more subtle thing is that the sign of the energy of wave is

[TABLE]

where $D(\omega)=0$ is a dispersion equation of wave, see Landahl (1962) and also Cairns (1979) for first discussion of this feature in hydrodynamics. In the next Section it is confirmed in application to the dynamics of gas-dust mixture. Since in the considered situation $\omega_{p}<0$, the coupling term is positive at the mode crossing $\epsilon|_{\omega_{c}}>0$, which implies that the instability occurs as soon as denominator in eq. (69) is negative, which in turn is true provided that the coupling waves have the energy of different signs. 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. The two cases described above, represent the pair of complex conjugate solutions of dispersion equation (60) discussed in Section 5.2.

5.4 Averaged variational principle

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

Let the modal solution have the form

[TABLE]

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

[TABLE]

According to the averaged variational principle, 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}\}$, which generates the set of eqs. (55-58) for modes of perturbations provided that the action similar to (47) be stationary with respect to arbitrary variations of $\tilde{\chi}_{i}$ and $\theta$. It reads

[TABLE]

It can be checked that eqs. (55-58) are equivalent to the Euler-Lagrange equations

[TABLE]

It is straightforward to check that equations (74) yield an alternative form of the dispersion equation (cf. eq. (60)), which is nothing but

[TABLE]

At the same time, the remaining Euler-Lagrange equation

[TABLE]

is satisfied identically since $\tilde{\chi}_{i}$ and $\omega$, $k_{x}$, $k_{z}$ are assumed to be constants.

Finally, the energy of neutral mode following from the averaged variational principal is explicitly

[TABLE]

Note that eq. (76) is identical to eq. (59) provided that $\Im[\omega]=0$. There is no cross term $\propto\tilde{\varpi}\tilde{\delta}$ in eq. (76) which leads to the conclusion that any perturbation of dust density gives negative amount to the energy of neutral mode with $\omega<0$. In the limit $f\to 0$ the first term in eq. (76) introduces the energy of IW, whereas the last term therein introduces the energy of SDW, which confirms the rule (71) employed for interpretation of SI, see Section 5.3. Clearly, in the case of dust settling to the disc midplane with the velocity (28) the mode crossing occurs at negative frequency $\omega_{c}=\omega_{p}<0$, therefore SDW has negative energy, while the energy of IW remains to be positive. As the dust fraction becomes non-negligible, whether the energy of the particular mode of gas-dust perturbations is either positive or negative is determined by the balance of terms in eq. (76) for the corresponding solution of eq. (60). Accurate computation demonstrates that the modes akin to IW and SDW have, respectively, positive and negative energies, see the branches (2) and (3) in Fig. 3 outside the band of SI.

6 Conclusions

The Lagrangian perturbation theory for local dynamics of gas-dust mixture with dust settling to protoplanetary disc midplane is constructed in order to show the existence of negative energy wave in such a flow. In general, the negative energy wave is akin to SDW, while it becomes identical to SDW in the absence of dust back reaction on gas. The dispersion equation for gas-dust perturbations with the account for dust back reaction on gas is reduced to the 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 IW, respectively. The dispersion curve of SDW crosses the dispersion curve of IW propagating downwards to the disc midplane at the point referred to as the mode crossing. It is revealed that SI occurs due to the energy transfer from (negative energy) SDW to (positive energy) IW in the vicinity of the mode crossing, where the phase velocities of both waves are equal to each other. The coupled modes are represented by the complex conjugate pair of frequencies, for which the analytical expressions are derived up to the next order in $f$, see eqs. (67) and (68). It is checked that the energy of the coupled modes vanishes.

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. The coupling term originates from the product of the last term in eq. (55) and the last term in eq. (58). In turn, the former comes from the last term in RHS of eq. (35), which is proportional to vertical 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. (58) comes from the first term in RHS of eq. (36), which exists as soon as the flow is rotational, see eq. (42). 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, see their discussion about the geostrophic balance. The models that introduce SI and the streaming instability differ from each other only by the direction of the effective gravity and the relative velocity in the background solution. Presumably, the streaming instability has the similar physical interpretation as given above for SI, however, the details are to be revealed in future work.

As one follows derivation of eq. (66) starting from the general equations (20-23), it becomes clear that the LHS of eq. (66) is nothing but the bulk settling velocity of dust given by eq. (28), 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. Clearly, this is because in the geometrically thin disc settling velocity is higher than the velocity of radial drift by factor $\sim(h/r_{0})^{-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 must be less sensitive to external turbulence inherent in astrophysical discs, see the recent numerical study of dust settling through turbulent gas by Lin (2019). On the contrary, enhanced effective viscosity produced by turbulence may cause an additional non-resonant growth of negative energy mode akin to SDW outside of the mode coupling.

At last, it is highly likely that the other resonant drag instabilities 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’).