Linear Analysis of the Nonaxisymmetric Secular Gravitational Instability

TL;DR

This paper conducts a linear analysis of nonaxisymmetric secular gravitational instability in protoplanetary disks, revealing transient growth of perturbations that can influence planet formation, especially under certain dust-gas coupling conditions.

Contribution

It provides the first linear perturbation analysis of nonaxisymmetric SGI in protoplanetary disks using a two-fluid model, highlighting transient amplification effects.

Findings

Nonaxisymmetric perturbations grow significantly over finite times.

Dust amplification occurs even when gas remains stable.

Amplification depends on dust-gas coupling and azimuthal wavelength.

Abstract

In protoplanetary discs (PPDs) consisting of gas and dust particles, fluid instabilities induced by the drag force, including secular gravitational instability (SGI) can facilitate planet formation. Although SGI subject to the axisymmetric perturbations was originally studied in the absence of gas feedback and it then generalized using a two-fluid approach, the fate of the nonaxisymmetric SGI, in either case, is an unexplored problem. We present a linear perturbation analysis of the nonaxisymmetric SGI in a PPD by implementing a two-fluid model. We explore the growth of the local, nonaxisymmetric perturbations using a set of linearized perturbation equations in a sheared frame. The nonaxisymmetric perturbations display significant growth during a finite time interval even when the system is stable against the axisymmetric perturbations. Furthermore, the surface density perturbations do not show continuous growth but are temporally amplified. We also study cases where the dust component undergoes amplification whereas the gas component remains stable. The amplitude amplification, however, strongly depends on the model parameters. In the minimum mass solar nebula (MMSN), for instance, the dust fluid amplification at the radial distance 100 au occurs when the Stokes number is about unity. But the amplification factor reduces as the dust and gas coupling becomes weaker. Furthermore, perturbations with a larger azimuthal wavelength exhibit a larger amplification factor.

Tables

Table 1. Table 1: Our model parameters in the MMSN model with a solar mass host star at the radial distance 100 au.

$ϵ$	$\alpha$	$Q_g$	$Q_d$	$\mathrm{St}$	$a\left(\mu m\right)$
$0.01$	${10}^{-4}$	$17.7$	$6.66$	$10$	$55555$
$0.01$	${10}^{-4}$	$17.7$	$18.24$	$1$	$5555$
$0.01$	${10}^{-4}$	$17.7$	$17.71$	$0.1$	$555$
$0.01$	${10}^{-4}$	$17.7$	$17.70$	$0.01$	$55$
$0.01$	${10}^{-4}$	$17.7$	$17.70$	$0,001$	$5.5$
$0.01$	${10}^{-4}$	$17.7$	$17.70$	$0.0001$	$0.5$

Equations

$$\frac{\partial {\Sigma }_{\mathrm{g}}}{\partial t}+\nabla .({\Sigma }_{\mathrm{g}}{\mathbf{V}}_{\mathrm{g}})=0,$$

$${\Sigma }_{\mathrm{g}}\left(\frac{\partial {\mathbf{V}}_{\mathrm{g}}}{\partial t}+({\mathbf{V}}_{\mathrm{g}}.\nabla ){\mathbf{V}}_{\mathrm{g}}+2\mathbf{\Omega }\times {\mathbf{V}}_{\mathrm{g}}-{\Omega }^2\mathbf{r}\right)=-c_{\mathrm{s}}^2\nabla {\Sigma }_{\mathrm{g}}$$

$$-{\Sigma }_{\mathrm{g}}\nabla ({\psi }_{\mathrm{g}}+{\psi }_{\mathrm{d}})+\frac{{\Sigma }_{\mathrm{d}}({\mathbf{V}}_{\mathrm{d}}-{\mathbf{V}}_{\mathrm{g}})}{t_{\mathrm{stop}}}$$

$$+\frac{\partial }{\partial x_{\mathrm{k}}}[{\Sigma }_{\mathrm{g}}\nu (\frac{\partial v_{\mathrm{ig}}}{\partial x_{\mathrm{k}}}+\frac{\partial v_{\mathrm{kg}}}{\partial x_{\mathrm{i}}}-\frac{2}{3}{\delta }_{\mathrm{ik}}\frac{\partial v_{\lg }}{\partial x_{\mathrm{l}}})],$$

$$\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial t}+\nabla .({\Sigma }_{\mathrm{d}}{\mathbf{V}}_{\mathrm{d}})=D{\nabla }^2{\Sigma }_{\mathrm{d}},$$

$${\Sigma }_{\mathrm{d}}\left(\frac{\partial {\mathbf{V}}_{\mathrm{d}}}{\partial t}+({\mathbf{V}}_{\mathrm{d}}.\nabla ){\mathbf{V}}_{\mathrm{d}}+2\mathbf{\Omega }\times {\mathbf{V}}_{\mathrm{d}}-{\Omega }^2\mathbf{r}\right)=-c_{\mathrm{d}}^2\nabla {\Sigma }_{\mathrm{d}}$$

$$-{\Sigma }_{\mathrm{d}}\nabla ({\psi }_{\mathrm{g}}+{\psi }_{\mathrm{d}})+\frac{{\Sigma }_{\mathrm{d}}({\mathbf{V}}_{\mathrm{g}}-{\mathbf{V}}_{\mathrm{d}})}{t_{\mathrm{stop}}},$$

$${\nabla }^2({\psi }_{\mathrm{g}}+{\psi }_{\mathrm{d}})=4\pi G({\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}})\delta (z).$$

$$c_{\mathrm{d}}^2=\alpha c_{\mathrm{s}}^2[\frac{1+2(\mathrm{St}{)}^2+\frac{5}{4}(\mathrm{St}{)}^3}{(1+(\mathrm{St}{)}^2{)}^2}]$$

$$\xi =\alpha [\frac{1+(\mathrm{St})+4(\mathrm{St}{)}^2}{(1+(\mathrm{St}{)}^2{)}^2}].$$

$$\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial t}+2Ax\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial y}+{\Sigma }_{0g}(\frac{\partial (\delta v_{x\mathrm{g}})}{\partial x}+\frac{\partial (\delta v_{y\mathrm{g}})}{\partial y})=0,$$

$$\frac{\partial (\delta v_{x\mathrm{g}})}{\partial t}+2Ax\frac{\partial (\delta v_{x\mathrm{g}})}{\partial y}-2\Omega (\delta v_{y\mathrm{g}})=-\frac{\partial }{\partial x}(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})$$

$$-\frac{c_{\mathrm{s}}^2}{{\Sigma }_{0g}}\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial x}+\frac{ϵ(\delta v_{x\mathrm{d}}-\delta v_{x\mathrm{g}})}{t_{\mathrm{stop}}}$$

$$+\frac{4}{3}\nu \frac{{\partial }^2(\delta v_{x\mathrm{g}})}{\partial x^2}+\nu \frac{{\partial }^2(\delta v_{x\mathrm{g}})}{\partial y^2}+\frac{2A\nu }{{\Sigma }_{0g}}\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial y}+\nu \frac{{\partial }^2}{\partial y\partial x}(\delta v_{y\mathrm{g}}),$$

$$\frac{\partial (\delta v_{y\mathrm{g}})}{\partial t}+2Ax\frac{\partial (\delta v_{y\mathrm{g}})}{\partial y}+2B(\delta v_{x\mathrm{g}})=-\frac{\partial }{\partial y}(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})$$

$$-\frac{c_{\mathrm{s}}^2}{{\Sigma }_{0g}}\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial y}+\frac{ϵ(\delta v_{y\mathrm{d}}-\delta v_{y\mathrm{g}})}{t_{\mathrm{stop}}}$$

$$+\frac{4}{3}\nu \frac{{\partial }^2(\delta v_{y\mathrm{g}})}{\partial y^2}+\nu \frac{{\partial }^2(\delta v_{y\mathrm{g}})}{\partial x^2}+\frac{2A\nu }{{\Sigma }_{0g}}\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial x}+\nu \frac{{\partial }^2}{\partial x\partial y}(\delta v_{x\mathrm{g}}),$$

$$\frac{\partial (\delta {\Sigma }_{\mathrm{d}})}{\partial t}+2Ax\frac{\partial (\delta {\Sigma }_{\mathrm{d}})}{\partial y}+{\Sigma }_{0d}(\frac{\partial (\delta v_{x\mathrm{d}})}{\partial x}+\frac{\partial (\delta v_{y\mathrm{d}})}{\partial y})$$

$$=D(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2})\delta {\Sigma }_{\mathrm{d}},$$

$$\frac{\partial (\delta v_{x\mathrm{d}})}{\partial t}+2Ax\frac{\partial (\delta v_{x\mathrm{d}})}{\partial y}-2\Omega (\delta v_{y\mathrm{d}})=-\frac{\partial }{\partial x}(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})$$

$$-\frac{c_{\mathrm{d}}^2}{{\Sigma }_{0d}}\frac{\partial (\delta {\Sigma }_{\mathrm{d}})}{\partial x}+\frac{(\delta v_{x\mathrm{g}}-\delta v_{x\mathrm{d}})}{t_{\mathrm{stop}}},$$

$$\frac{\partial (\delta v_{y\mathrm{d}})}{\partial t}+2Ax\frac{\partial (\delta v_{y\mathrm{d}})}{\partial y}+2B(\delta v_{x\mathrm{d}})=-\frac{\partial }{\partial y}(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})$$

$$-\frac{c_{\mathrm{d}}^2}{{\Sigma }_{0d}}\frac{\partial (\delta {\Sigma }_{\mathrm{d}})}{\partial y}+\frac{(\delta v_{y\mathrm{g}}-\delta v_{y\mathrm{d}})}{t_{\mathrm{stop}}},$$

$$(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2})(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})=4\pi G(\delta {\Sigma }_{\mathrm{g}}+\delta {\Sigma }_{\mathrm{d}})\delta (z),$$

$$x^{\prime }=x,y^{\prime }=y-2Axt,z^{\prime }=z,t^{\prime }=t,$$

$$\frac{\partial }{\partial x}\equiv \frac{\partial }{\partial x^{\prime }}-2A{t}^{\prime }\frac{\partial }{\partial y^{\prime }},\frac{\partial }{\partial y}\equiv \frac{\partial }{\partial y^{\prime }},$$

$$\frac{\partial }{\partial t}\equiv \frac{\partial }{\partial t^{\prime }}-2A{x}^{\prime }\frac{\partial }{\partial y^{\prime }},\frac{\partial }{\partial z}\equiv \frac{\partial }{\partial z^{\prime }}.$$

$$\frac{\partial (\delta {\Sigma }_{\mathrm{g}})}{\partial t^{\prime }}+{\Sigma }_{0g}(i{k}_x-2i{k}_{y}A{t}^{\prime })\delta v_{x\mathrm{g}}+i{k}_y{\Sigma }_{0g}(\delta v_{y\mathrm{g}})=0,$$

$$\frac{\partial (\delta v_{x\mathrm{g}})}{\partial t^{\prime }}-2\Omega (\delta v_{y\mathrm{g}})=(i{k}_x-2i{k}_{y}A{t}^{\prime })[-(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})$$

$$-\frac{c_{\mathrm{s}}^2}{{\Sigma }_{0g}}(\delta {\Sigma }_{\mathrm{g}})]+\frac{ϵ(\delta v_{x\mathrm{d}}-\delta v_{x\mathrm{g}})}{t_{\mathrm{stop}}}-\frac{4}{3}\nu (k_x-2k_{y}A{t}^{\prime }{)}^2(\delta v_{x\mathrm{g}})$$

$$-\nu k_y^2(\delta v_{x\mathrm{g}})+\frac{2iA{k}_y\nu }{{\Sigma }_{0g}}(\delta {\Sigma }_{\mathrm{g}})-\nu k_y(k_x-2k_{y}A{t}^{\prime })(\delta v_{y\mathrm{g}}),$$

$$\frac{\partial (\delta v_{y\mathrm{g}})}{\partial t\prime }+2B(\delta v_{x\mathrm{g}})=i{k}_y[-(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})-\frac{c_{\mathrm{s}}^2}{{\Sigma }_{0g}}(\delta {\Sigma }_{\mathrm{g}})]$$

$$+\frac{ϵ(\delta v_{y\mathrm{d}}-\delta v_{y\mathrm{g}})}{t_{\mathrm{stop}}}-\frac{4}{3}\nu k_y^2(\delta v_{y\mathrm{g}})-\nu (k_x-2k_{y}A{t}^{\prime }{)}^2(\delta v_{y\mathrm{g}})$$

$$+\frac{2A\nu }{{\Sigma }_{0g}}(i{k}_x-2i{k}_{y}A{t}^{\prime })(\delta {\Sigma }_{\mathrm{g}})-\nu k_y(k_x-2k_{y}A{t}^{\prime })(\delta v_{x\mathrm{g}}),$$

$$\frac{\partial (\delta {\Sigma }_{\mathrm{d}})}{\partial t^{\prime }}+{\Sigma }_{0d}(i{k}_x-2i{k}_{y}A{t}^{\prime })\delta v_{x\mathrm{d}}+i{k}_y{\Sigma }_{0d}(\delta v_{y\mathrm{d}})$$

$$=D[-(k_x-2k_{y}A{t}^{\prime }{)}^2-k_y^2]\delta {\Sigma }_{\mathrm{d}},$$

$$\frac{\partial (\delta v_{x\mathrm{d}})}{\partial t^{\prime }}-2\Omega (\delta v_{y\mathrm{d}})=(i{k}_x-2i{k}_{y}A{t}^{\prime })[-(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})$$

$$-\frac{c_{\mathrm{d}}^2}{{\Sigma }_{0d}}(\delta {\Sigma }_{\mathrm{d}})]+\frac{(\delta v_{x\mathrm{g}}-\delta v_{x\mathrm{d}})}{t_{\mathrm{stop}}},$$

$$\frac{\partial (\delta v_{y\mathrm{d}})}{\partial t^{\prime }}+2B(\delta v_{x\mathrm{d}})=i{k}_y[-(\delta {\psi }_{\mathrm{g}}+\delta {\psi }_{\mathrm{d}})-\frac{c_{\mathrm{d}}^2}{{\Sigma }_{0d}}(\delta {\Sigma }_{\mathrm{d}})]$$

Full text

††pubyear: 2017††pagerange: Linear Analysis of the Nonaxisymmetric Secular Gravitational Instability–Linear Analysis of the Nonaxisymmetric Secular Gravitational Instability

1 Introduction

Understanding planet formation mechanisms in protoplanetary discs (PPDs) is still a controversial issue despite considerable achievements in recent years (Lissauer & Stevenson, 2007; Durisen et al., 2007; Helled et al., 2014). The most widely studied planet formation theories are the so-called core accretion model which is efficient in the inner region of a PPD (e.g., Mizuno, 1980; Stevenson, 1982; Pollack et al., 1996) and gravitational instability (GI) which operates mainly in the outer regions (e.g., Adams et al., 1989; Boss, 1997; Boley, 2009; Rafikov, 2009). However, alternative planet formation theories have been proposed during recent years, including inside-out planet formation scenario (Chatterjee & Tan, 2014) and tidal downsizing model (Nayakshin, 2010). In the core accretion model, when the mass of an already formed rocky core becomes about $10M_{\oplus}$, this planetary embryo is able to grow further through the accretion of its ambient gas and dust particles onto it (Cameron, 1973; Hayashi et al., 1977; Lissauer, 1993; D’Angelo et al., 2010).

Linear perturbation analysis and numerical simulations show that a gaseous disc with the surface density $\Sigma$, Keplerian angular velocity $\Omega$ and the sound speed $c_{\rm s}$ is gravitationally unstable subject to the axisymmetric perturbations if Toomre parameter, i.e., $Q\equiv c_{\rm s}\Omega/(\pi G\Sigma)$ becomes less than a critical value around unity (Toomre, 1964). This criterion has been successfully implemented in the star formation theories in galaxies (e.g., Collin & Zahn, 2008; Krumholz & Burkert, 2010; Romeo & Agertz, 2014; Goldbaum et al., 2016) and planet formation scenarios in the PPDs (e.g., Matzner & Levin, 2005; Kratter et al., 2008; Rafikov, 2005; Boley, 2009; Rafikov, 2009). Gammie (2001) suggested that Toomre diagnostic is only a necessary condition for disc fragmentation and it does not guarantee to have long-lived fragments in a turbulent self-gravitating disc (also see, Rice et al., 2003; Mejía et al., 2005; Boss, 2017). But recent numerical simulations of the self-gravitating discs have shown that small enough cooling time scale is not a sufficient condition for fragmentation (Tsukamoto et al., 2015). Takahashi et al. (2016) showed that formation of the spiral arms in a PPD and their fragmentation are the essential phases of a disc fragmentation. They proposed a revised condition for the fragmentation.

While the primary focus of most disc stability studies is the gas component and its evolution, about one percent of a PPD total mass is attributed to the dust particles with different sizes (Natta et al., 2007). Presence of dust particles provides valuable insights about dynamical and chemical structure of a PPD (Vasyunin et al., 2011; Akimkin et al., 2013; Woitke et al., 2016; Rab et al., 2017), its ionization level (Okuzumi et al., 2011; Akimkin, 2015; Ivlev et al., 2016; Rab et al., 2017) and radiative transfer through a disc (Akimkin et al., 2013; Rab et al., 2017). Formation of rocky planets or cores of gaseous giant planets is explained based on the collective instabilities associated with the dust component. In the early works on the clumping of dust particles either purely dusty discs have been considered (Safronov, 1972; Goldreich & Ward, 1973) or their implemented approximations were rather restrictive (Coradini et al., 1981; Sekiya, 1983; Noh et al., 1991).

More recent progress for understanding mechanisms of dust clumping in a PPD relies on the existence of a relative velocity between gas and dust particles and the associated drag force. A dust layer is generally believed to be formed at a PPD midplane due to the sedimentation of these particles. Their rotational velocity is Keplerian, whereas the gas component is rotating with a sub-Keplerian velocity because of the radial gradient of pressure. In a PPD as a mixture of gas and dust, therefore, these components may experience a relative velocity. Since the aerodynamically friction force is directly proportional to this relative velocity, motions of dust particles are affected by this force. The gas motion, however, is unaffected by the dust movement when the dust-to -gas density ratio is small.

A mechanism of the so-called streaming instability (SI) is triggered because of the dust movement through the gas (Youdin & Goodman, 2005; Youdin & Lithwick, 2007; Jacquet et al., 2011). However, its efficiency strongly depends on the dust-to-gas density ratio and the dust-gas coupling which is quantified in terms of the Stokes number, i.e., ${\rm St}=t_{\rm stop}\Omega$ where $t_{\rm stop}$ is the stopping time and $\Omega$ is Keplerian angular velocity. For relatively well-coupled particles (${\rm St}<1$), a mechanism known as secular gravitational instability (SGI) is driven by the drag force (Youdin, 2011; Shariff & Cuzzi, 2011; Michikoshi et al., 2012). This mechanism is a dissipative version of the classical Toomre analysis for the two-component discs (Youdin, 2011; Shariff & Cuzzi, 2011; Michikoshi et al., 2012; Takahashi & Inutsuka, 2014). The SGI has primarily been studied using a single fluid model where the dust dynamics is treated in a gaseous background with no backreaction of the dust. The onset of SGI, therefore, is found to be unconditional and it is triggered no matter how dusty layer is thin or thick (Youdin, 2011). This interesting feature of the SGI is lost when the gas feedback is included in a two-fluid model (Takahashi & Inutsuka, 2014). Numerical simulations of the SGI, however, are needed to address whether this mechanism leads to an appreciable enhancement of the dust surface density (Tominaga et al., 2018).

Recently observed multiple concentric ringlike structures in the PPDs are commonly interpreted as a sign of the newly born planets at their early formation phase (Andrews et al., 2011; Mayama et al., 2012; Yen et al., 2016; Loomis et al., 2017; Hendler et al., 2017; Dipierro et al., 2018; van Terwisga et al., 2018). But we note that a planet associated with these rings or gaps has not yet been observed directly. Non-planet-related scenarios, thereby, have also been proposed as mechanisms of the multiple ring formation (e.g., Okuzumi et al., 2016; Suriano et al., 2018). Takahashi & Inutsuka (2016) applied their two-fluid model of the SGI to explain HL Tau rings resulting from this instability. They studied growth time-scale and unstable wavelengths and proposed that SGI is able to create ring structures in the HL Tau disc. Latter & Rosca (2017), on the other hand, suggested that these ring-like structures probably unrelated to the SGI. In an alternative mechanism (Okuzumi et al., 2016), however, ring-like feature in the HL Tau disc is explained by incorporating sintering in a dust growth model.

Although the focus of the recent studies is to provide an explanation for the ring-like structures, some discs also exhibit complex non-axisymmetric patterns. Spiral arm-like structures, for instance, have been observed in the PPDs such as MWC 7588 (Grady et al., 2013; Benisty et al., 2015) and SAO 206462 (Muto et al., 2012; Garufi et al., 2013). Spiral arms that result from non-axisymmetric perturbations are observed in near infrared scattered light (Muto et al., 2012; Wagner et al., 2015) and also in sub-millimeter (Tobin et al., 2016).

These spiral features are commonly studied in terms of planet-disc interactions (Zhu et al., 2015; Lee, 2016; Bae & Zhu, 2017). In this scenario, spiral waves are generated due to the interaction of a planer with epicyclic oscillations in a disc. A set of these waves, which are in phase, create constructive interference and can create spiral arms (Ogilvie & Lubow, 2002; Zhu et al., 2015; Lee, 2016).

In an alternative explanation, however, the spiral arms are generated due to the non-axisymmetric development of the gravitational instability (Rice et al., 2003; Dong et al., 2015; Tomida et al., 2017; Juhász & Rosotti, 2018). Spiral arm formation in the galactic scale is a well-studied problem known as spiral density wave theory (Lin & Shu, 1964; Julian & Toomre, 1966; Toomre, 1981), and we think that a similar approach can be implemented in scenarios of GI-induced structures in PPDs. In this regard, SGI provides a natural route of ring-like patterns in a PPD. But prior linear studies of SGI are restricted to only axisymmetric perturbations. Therefore, the recent observed non-axisymmetric features in PPDs motivates us to investigate SGI subject to the non-axisymmetric perturbations.

In this paper, we generalize SGI to a case with nonaxisymmetric perturbations in a disc composed of the gas and dust particles. Our sheared two-fluid disc model and the main equations are presented in Section 2. We then present a linear analysis and a set of ordinary differential equations are obtained for the evolution of the perturbations in Section 3. These equations are solved numerically in Section 4. We conclude with a summary of our main results in Section 5.

2 Basic Equations

We persist with already studied SGI models and their general formulation (e.g., Takahashi & Inutsuka, 2014), however, our focus is to explore fate of the non-axisymmetric perturbations. A shearing sheet model (Goldreich & Lynden-Bell, 1965), as a representation of a razor thin disc small portion, is constructed where its center is at a fixed radial distance $r_{0}$. The $x$-axis is defined in the radial direction and the $y$-axis is oriented in the azimuthal direction. Thus, position of any point in this rotating plane is $(x,y)=(r-r_{\rm 0},r_{\rm 0}(\theta-\Omega t))$, where $\Omega=(GM/r_{0}^{3})^{1/2}$ is the Keplerian angular velocity at $r=r_{0}$. Here, $M$ is the star mass. In our analysis we use the relation $\Omega=-\Omega_{0}=B-A$, where $A$ and $B$ are the standard Oort constants and ${\bf\Omega}=\Omega{\bf z}$ is the Keplerian angular velocity.

Our basic equations, therefore, are the continuity, momentum, and the Poisson equations for the dust fluid and gas component:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The subscripts "g" and "d" stand for the gas and dust components. Thus, $\Sigma_{\rm g}$ and $\Sigma_{\rm d}$ are the gas and the dust surface densities. The velocities of the gas and the dust fluids are denoted by $\bf V_{\rm g}$ and $\bf V_{\rm d}$. Furthermore, $\psi_{\rm g}$ and $\psi_{\rm d}$ represent the gravitational potential associated with the gas and dust components. The gas sound speed and the velocity dispersion of the dust particles are $c_{\rm s}$ and $c_{\rm d}$ respectively. Here, $\nu$ is the kinematic viscosity which is written as $\nu=\alpha c_{\rm s}^{2}\Omega^{-1}$, where $\alpha$ is the dimensionless measure of the turbulent strength. Finally, $D$ is the diffusivity of the dust due to the gas turbulence.

In order to close the equations of the system, a relation between $c_{\rm s}$ and $c_{\rm d}$ and an equation for the diffusivity coefficient $D$ are needed. Youdin & Lithwick (2007) found the following relation:

[TABLE]

where ${\rm St}$ is the Stokes number, i.e. ${\rm St}=t_{\rm stop}\Omega_{0}$. Furthermore, the diffusivity of the dust $D$ is written as $D=\xi c_{\rm s}^{2}\Omega_{0}^{-1}$, where $\xi$ is the dimensionless diffusion coefficient and is defined as

[TABLE]

The Stokes number depends on the particle size, $a$, and the mean free path of the molecules, $\lambda$. When the size of particles is smaller than the mean free path of the molecules, i.e. $a<9\lambda/4$, which is known as Epstein regime, the Stokes number is written as ${\rm St}=({\rho_{\rm m}a}/{\rho_{g}v_{\rm th}})f_{\rm d}^{-1}\Omega$, where $\rho_{g}$ is the gas density and $\rho_{\rm m}=2$ g cm*-3* denotes the homogeneous material density of a dust particle. We also have $v_{\rm th}=(8/\pi)^{1/2}c_{\rm s}$ and $f_{d}=[1+({9\pi}/{128})({\left\|\Delta v\right\|}/{c_{\rm s}})^{2}]^{1/2}\simeq 1$, where the relative velocity between dust and gas, i.e. $\Delta v$, is much smaller than the sound speed at the disc midplane (e.g., Miyake et al., 2016).

3 Linear perturbations

The two-fluid equations admit an equilibrium configuration with a constant gas surface density $\Sigma_{\rm 0g}$ and a constant dust surface density $\Sigma_{\rm 0d}$. The dust-to-gas density ratio is defined via $\epsilon=\Sigma_{\rm 0d}/\Sigma_{\rm 0g}$ which is an input model parameter. The gas and dust components are undergoing Keplerian motion, i.e. ${\bf V}_{\rm 0g}={\bf V}_{\rm 0d}=2Ax{\bf j}$. Therefore, the linearised equations are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\delta(z)$ is Dirac delta function.

For exploring non-axisymmetric perturbations, it is more convenient to use a shearing sheet coordinates $(x^{\prime},y^{\prime},z^{\prime})$ introduced by Goldreich & Lynden-Bell (1965). We have

[TABLE]

where

[TABLE]

In this sheared system of coordinates, all disc quantities as perturbed as $\chi=\chi_{\rm 0}+\delta\chi$, where $\delta\chi$ is assumed to be small in comparison to the initial state. The disturbances, denoted by primes, are proportional to $\exp[i(k_{x}x^{\prime}+k_{y}y^{\prime})]$, where $k_{x}$ and $k_{y}$ are the radial and azimuthal wavenumbers respectively. Thus, the equations (8)-(14) become

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For perturbations with a non-zero $k_{y}$, it is more convenient to re-write equations (17)-(23) in terms of a new dimensionless time variable, i.e. $\tau=2At^{\prime}-k_{x}/k_{y}$ (Goldreich & Lynden-Bell, 1965). Therefore, we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Upon solving the Poisson equation (30) and using the approximation of the finite thickness of the disc (Vandervoort, 1970; Shu, 1984), we obtain the gravitational potential perturbation:

[TABLE]

where $H=c_{\rm s}/\Omega$ and $H_{\rm d}=\sqrt{(\alpha/{\rm St})H}$ are the gas disc scale height and the dust scale height respectively.

Equations (24)-(3) and (3) constitute main equations of the model to be solved subject to appropriate initial conditions. In the absence of drag force, mathematical forms of these equations are similar to Jog (1992) who studied GI in a two-component disc consisting of the gas and stars to mimic a galaxy. Our second component, however, is the dust fluid and its feedback on the gas is included via the drag force. Time evolution of the perturbations is studied by solving our main equations. For numerical purposes, however, it is better to re-write equations (24)-(3) and (3) in terms of dimensionless variables. In doing so, we introduce the following new variables:

[TABLE]

Here, parameters $Q_{\rm g}$ and $Q_{\rm d}$ stand for the Toomre parameter for the gas and dust components respectively. The epicyclic frequency is denoted by $\kappa$ where for a Keplerian disc it becomes $\kappa=\Omega$. Furthermore, the shear parameter is $\eta={2A}/{\Omega_{\rm 0}}$. Ratios of the gas and dust surface densities and their corresponding initial values are shown by $\theta_{\rm g}$ and $\theta_{\rm d}$ respectively. We also introduce a critical wavelength, i.e. $\lambda_{\rm crit}=4\pi^{2}G\Sigma_{\rm 0g}(1+\epsilon)/\kappa^{2}$, and the perturbation wavelength is written in terms of this critical wavelength.

Using introduced dimensionless variables, equations (24)-(3) and (3) are reduced to the following set of ordinary differential equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $W_{\rm g}=i[\delta u_{x\rm g}+\tau(\delta u_{y\rm g})]$ and $W_{\rm d}=i[\delta u_{x\rm d}+\tau(\delta u_{y\rm d})]$. Equations (3)-(3) are solved using Runge-Kutta method to determine evolution of the perturbations with time. In the next section, we present our solutions.

4 Numerical solutions

Although a wide range of the initial conditions can be implemented, we consider the following simple initial conditions. At the initial time $\tau_{\rm ini}$, the relevant quantities are

[TABLE]

When a perturbation oscillates with an amplitude less than unity, we consider it as a stable configuration. Note that stable axisymmetric perturbations also exhibit oscillatory behavior with time. The non-axisymmetric perturbations, however, are unstable once their amplitudes display significant growth only for a limited time. Unstable axisymmetric perturbations grow with an exponential profile, but growth of the nonaxisymmetric perturbations is not exponential and the concept of the instability corresponds to rapid transient amplification of the perturbations during a limited time. These transient amplifications have already been found in the linear stability analysis of the gaseous self-gravitating discs (Mamatsashvili et al., 2013), discs with gas and stars (Jog, 1992) and even nonaxisymmetric MRI (Balbus & Hawley, 1992). It is therefore quite normal that a mixture of gas and dust undergoes transient amplifications subject to the nonaxisymmetric perturbations. The fate of these transient patterns can be addressed in the non-linear regime, however, it is important to specify range of the model parameters for which the system becomes linearly unstable. We thereby consider initial states which are stable subject to the axisymmetric perturbations. But are they remain stable subject to the nonaxisymmetric perturbations? This is an important question that motivated us to perform stability analysis with a broad range of the model parameters.

We first investigate the evolution of the perturbations in Figure 1 for a fiducial set of the model parameters. We consider an initial configuration with the gas and dust large Toomre parameters to ensure stability of the system against to the axisymmetric perturbations. The model parameters, therefore, are assumed $Q_{\rm g}=15$, $Q_{\rm d}=11$, $\epsilon=0.01$, $X=3$, $\eta=1.5$ and $\alpha=10^{-4}$. The initial dimensionless time is $\tau_{\rm ini}=-10$. Note that we also verified that for other values of $\tau_{\rm ini}$ the results are qualitatively similar. We, however, note that the initial time $t=0$ corresponds to $\tau_{\rm ini}=-k_{x}/k_{y}$ that can be rewritten as $k_{x}=-\tau_{\rm ini}k_{y}$. Since the azimuthal wavenumber $k_{y}$ is a given model parameter in terms of the dimensionless parameter $X$, we can infer that perturbations with a larger $|\tau_{\rm ini}|$ have a larger radial wavenumber. Each panel of Figure 1 shows the evolution of the gas and dust perturbations, i.e. $\theta_{\rm g}$ and $\theta_{\rm d}$, as a function of $\tau$ for a given Stokes number, as labeled. Obviously, a given Stokes number is equivalent to a certain dust size. Black and red solid curves correspond to the gas and dust components respectively.

Figure 1 shows that the gas component is always stable irrespective of the Stokes number. Since the gas Toomre parameter is larger than its critical value, the gas component is stable subject to the perturbations and dust dynamics is unable to change this trend due to a small dust-to-gas density ratio. The response of the dust component, however, strongly depends on the adopted Stokes number that controls the magnitude of the drag force. For a large Stokes number, the dust component is also stable because the drag force is too weak to affect dust dynamics. But when the Stokes number is less than about 10, the dust component tends to be unstable subject to the nonaxisymmetric perturbations whereas the gas is still stable. We also find that the dust component is unstable for ${\rm St}=0.3$. If the Stokes number is adopted less than about 0.3, not only the dust component becomes stable but also its evolution is similar to the gas evolutionary trend. This behavior is understood in terms of the strong dust and gas coupling for the small Stokes numbers. Under this condition, the stability of the gas component is dictated to the dust component due to their strong coupling. Thus, for an intermediate range of the Stokes number $0.3\lesssim{\rm St}\lesssim 10$, while the gas component remains stable, the dust component undergoes transient growing patterns during a time interval almost independent of the gas component. But either the dust component undergoes a growing phase or just display oscillatory behaviour, the amplitudes of the perturbations decay at larger times.

In a realistic case, however, our model parameters can not be adopted independently as we did in Figure 1. These parameters depend on the gaseous disc model and the initial distribution of the dust particles. In agreement with most previous studies in this context, we consider the minimum mass solar nebula (MMSN; Hayashi, 1981) to represent our gaseous disc model. However, the entire structure of a PPD is unlikely to be described using the MMSN model and there are also alternative disc profiles (Nixon et al., 2018). Our stability analysis, therefore, is done in an MMSN model at a given radial distance. In this model, the surface density and sound speed are given as power-law functions of the radial distance (Hayashi, 1981):

[TABLE]

[TABLE]

In MMSN model with a solar mass host star, the gas Toomre parameter can therefore be expressed by

[TABLE]

Furthermore, the Toomre parameter associated to the dust component becomes

[TABLE]

The Stokes number therefore at the disc midplane becomes

[TABLE]

We can use the above relations for determining model parameters self-consistently in the MMSN model at a given radial location. In Table 1, we list our model parameters for the MMSN model at the radial distance 100 au.

Figure 2 displays the evolution of the gas surface density perturbation (black curve) and the dust surface density (red curve) in MMSN model with $M=1$ M*⊙* at the radial distance 100 au and for different Stokes numbers, as labeled. The gas Toomre parameter is $Q_{\rm g}=17.7$, whereas the Toomre parameter for the dust component becomes $Q_{\rm d}=18.24$, 18.43 and 17.89 corresponding to the Stokes numbers 1, 0.7 and 0.3, respectively. Note that for these Stokes numbers, dust particles are millimeter-sized in a range between 1.6 mm and 5.5 mm. We have chosen other model parameters as $\epsilon=0.01$, $X=2.78$, $\eta=1.5$ and $\alpha=10^{-4}$. We verified that both the gas and dust fluids are stable subject to the axisymmetric disturbances. But the response of the system to the nonaxisymmetric perturbations is quite different as we see in Figure 2. While the gas component remains stable subject to these perturbations irrespective of the adopted Stokes number, the dust component displays noticeable growth and the amplification factor reduces with decreasing the Stokes number from 1 to 0.3.

The dust abundance is a key parameter and its role is explore in Figure 3 for different values of $\epsilon$, as labeled. As in the previous figure, the model parameters correspond to MMSN model at the radial distance 100 au. We set ${\rm St}=0.3$, $Q_{\rm g}=17.7$, $\eta=1.5$ and $\alpha=10^{-4}$. For $\epsilon=0.1$, $0.01$ and $0.001$, we have $Q_{\rm d}=1.78$, 17.89 and 178.9 and the corresponding azimuthal dimensionless parameter becomes $X=2.5$, $2.78$ and $2.8$, respectively. As before, the integration is started at $\tau_{\rm ini}=-10$ and the initial condition (37) is implemented to determine evolution of the perturbations. The gas component exhibits an oscillatory behavior, however, its amplitude significantly decreases with increasing dust abundance. But the dust component undergoes a slightly higher growth rate when its abundance is increased.

In Figure 4, we investigate the role of the parameters $\eta$ and $\alpha$ in the evolution of the disturbances. Model parameters correspond to MMSN model at the radial distance 100 au and the runs are started at $\tau_{\rm ini}=-10$. In both panels, we adopt ${\rm St}=0.3$ and $\epsilon=0.01$. In the top panel, we have $\alpha=10^{-4}$ and different values of the shear parameter are considered, as labeled. For the adopted values of the shear parameter, i.e. $\eta=1$, 1.5 and 1.8, the associated gas and dust Toomre parameter become $Q_{\rm g}=249.9$, 17.7 and 2.8, and, $Q_{\rm d}=252.6$, 17.89 and 2.83, and, we also have $X=557.3$, 2.78 and 0.07, respectively. In the bottom panel, we assume that $Q_{\rm g}=17.7$, $X=2.78$ and $\eta=1.5$ and the dust Toomre parameter for $\alpha=10^{-3}$, $10^{-4}$ and $10^{-5}$ becomes $Q_{\rm d}=56.57$, $17.89$ and $5.65$, respectively. Note that both the gas and dust components are stable subject to the axisymmetric disturbances. In the top panel of Figure 4, role of $\eta$ that quantifies the differential rotation rate is shown for different values of this parameter. All model parameters are similar to the bottom panel and the viscosity coefficient is $\alpha=10^{-4}$. As the parameter $\eta$ decreases, the dust component tends to be more unstable. The bottom panel of Figure 4 shows that fate of the dust component against to the nonaxisymmetric perturbations strongly depends on the chosen viscosity coefficient. While for $\alpha=10^{-3}$, the dust fluid remains stable, the amplitude of the perturbations gradually increases as the viscosity coefficient tends to the smaller values.

In our stability analysis, we introduced the dimensionless parameter $X$ that measures azimuthal wavelength of the perturbations in terms of a critical wavelength $\lambda_{\rm crit}$. So far this parameter has been included as a fixed value. In Figure 5, we explore role of this parameter in evolution of the dust component in MMSN model at the radial distance 100 au for these model parameters: $Q_{\rm g}=17.7$, $Q_{\rm d}=21.26$, ${\rm St}=0.3$, $\epsilon=0.01$, $\eta=1.5$, $\alpha=10^{-4}$. As before, the starting time is $\tau_{\rm ini}=-10$. Each curve is marked with the corresponding value of $X$. Comparison of these curves show that the amplitude of the oscillations increases with increasing $X$.

In Figure 6, we consider typical cases that satisfy conditions of the axisymmetric stability (Takahashi & Inutsuka, 2014) but finite amplifications are found subject to the nonaxisymmetric perturbations. We set $X=2.78$, $\eta=1.5$, $\epsilon=0.01$ and $\alpha=10^{-4}$. In the top left hand panel of Figure 6, variations in $\theta_{\rm g}$ and $\theta_{\rm d}$ with $\tau$ are shown for ${\rm St}=0.7$, $Q_{\rm g}=7$ and $Q_{\rm d}=7.29$. The dust component is weakly coupled to the gas component and the corresponding perturbations exhibit a significant growth during a finite time period, whereas the gas component remains stable. In the top right hand panel, the model parameters are ${\rm St}=0.9$, $Q_{\rm g}=7$ and $Q_{\rm d}=7.26$. The dust component again displays a fairly strong amplification for this choice of the parameters. In the bottom panels, we explore stability of the cases with a smaller Toomre parameter. In the bottom left hand panel, we set ${\rm St}=0.7$, $Q_{\rm g}=5$ and $Q_{\rm d}=5.2$, whereas the bottom right hand panel corresponds to a case with the same gas and dust Toomre parameter but with a slightly larger Stokes number, i.e. ${\rm St}=0.9$. Since the Toomre parameter associated with each component is smaller, the amplitude of the perturbations are larger for both cases with Stokes numbers $\rm St=0.7$ and $\rm St=0.9$ . We find that amplitude of the nonaxisymmetric perturbations increases with time and then these oscillations are damping.

In all explored cases so far, we used a fixed initial dimensionless time $\tau_{\rm ini}=-10$ which then it corresponds to a given radial perturbation wavelength $\lambda_{x}$ if the azimuthal wavelength is treated as a given fixed value. This argument is based on an already introduced relation as follows $k_{x}=-\tau_{\rm ini}k_{y}$. We now investigate evolution of the perturbations with different initial dimensionless time $\tau_{\rm ini}$. This analysis, thereby, corresponds to evolution of the perturbations with different radial perturbation wavelength for a given fixed azimuthal wavelength. In the top panel of Figure 7, we exhibit evolution of the dust component for different values of the initial time $\tau_{\rm ini}$, as labeled. The other model parameters are ${\rm St}=0.7$, $\epsilon=0.01$, $Q_{\rm g}=17.7$, $Q_{\rm d}=18.43$, $X=2.78$, $\eta=1.5$ and $\alpha=10^{-4}$. Note that the implemented values $\tau_{\rm ini}=-1$, -5, -10 and -15 correspond to radial wavelength $\lambda_{x}=21.23$ au, 4.24 au, 2.12 au and 1.41 au, respectively. This plot shows that the amplitude of the perturbations gradually increases with time, however, its growth is suppressed after a certain period of time depending on the radial perturbation wavelength. In the middle panel of Figure 7, we consider cases with different initial time $\tau_{\rm ini}$ in an interval from $-1$ to $-20$ and the resulting maximum amplitude for the gas and dust components, i.e. $(\theta_{\rm g})_{\rm Max}$ and $(\theta_{\rm d})_{\rm Max}$ are shown as a function of the radial perturbation wavenumber $\lambda_{x}$. The adopted model parameters are $\rm St=0.7$, $Q_{\rm g}=17.7$, $Q_{\rm d}=18.43$, $\epsilon=0.01$, $X=2.78$, $\alpha=10^{-4}$ and for a Keplerian disc with $\eta=1.5$. We find that for $\lambda_{x}\leq 10$ au the dust component generally undergoes a larger $(\theta_{\rm d})_{\rm Max}$ in comparison to the gas component. In other words, the dust component is more unstable in comparison to the gas component for the short radial perturbation wavelength. In the bottom panel, we display growth time (i.e., the time when the amplitude becomes maximum) as a function of the radial perturbation wavelength for the explored cases. It shows that both dust and gas component evolves to their maximum amplitude during more or less the same time period. This time scale, however, is in an interval from 1000 yr (for long wavelengths) to 3000 yr (for short wavelengths). We find that evolution of the nonaxisymmetric growth time scale is relatively fast in comparison to the axisymmeric growth time scale.

5 Conclusions

We investigated evolution of the imposed nonaxisymmetric perturbations in a PPD by treating the system as a mixture of the coupled gas and dust particles. While response of the system subject to the axisymmetric perturbations may become unstable with an exponential growth, we find that nonaxisymmetric perturbations evolve with an oscillatory amplitude. But a strong amplification is found corresponding to the cases where are stable against to the axisymmetric perturbations. We can now summarize our main results:

• Growth of the nonaxisymmetric perturbations is not significant when the Stokes number is very large or very small. But there is always an intermediate range for the Stokes number where a high amplification in dust against to the nonaxisymmetric perturbations is found even when the gas component remains stable. However, the upper and lower limits of the Stokes number corresponding to the instability depend on the other model parameter including dust and gas Toomre parameters.

• In the MMSN model with a solar mass host star, we found that amplitude of the nonaxisymmetric perturbations at the radial distance 100 au increases when the Stokes number lies in a range between $10^{-4}$ and $10^{-2}$. The amplification factor, however, decreases with increasing the Stokes number.

• Perturbations with a larger azimuthal wavelength exhibit a relatively higher amplification.

• Turbulent coefficient has a stabilizing role in promoting nonaxisymmetric SGI. In the MMSN model, for example, we showed that growth of the perturbation is completely suppressed when the coefficient is $\alpha=10^{-3}$, whereas for smaller values of this coefficient, the amplification factor increases.

The final outcome of the implemented nonaxisymmetric perturbations cannot be addressed using the present linear analysis. Our study demonstrates that nonaxisymmetric perturbations may lead to dust transient amplification for a wider range of the model parameters in comparison to the axisymmetric perturbations. In the context of the nonaxisymmetric gravitational instability in galaxies, many authors found similar evolutionary trends in either purely gaseous discs or two-fluid (stars and gas) systems (e.g., Julian & Toomre, 1966; Jog, 1992; Fuchs, 2001; Michikoshi & Kokubo, 2016; Ghosh & Jog, 2018). The imposed nonaxisymmetric perturbations get amplified during an initial time interval due to the mutual interplay between various physical agents including shear and self-gravity of the disc, but eventually this amplification is suppressed because of the disc shear. Our model also predicts a similar trend for the evolution of the dust component even when the gas component remains stable. Although our basic equations are similar to previous studies relevant to gravitational instability in galaxies, there is a significant difference because each component (gas or dust) is permitted to exchange momentum via the drag force. In other words, dust and gas components are gravitationally coupled and they are subject to the momentum exchange by the drag force. We then found that the drag force is able to promote temporal growth of the dust component even in the cases where both components are stable subject to the axisymmetric perturbations. Furthermore, amplitude amplification of the dust component may persists in configurations where the gaseous disc responds via oscillations with non growing amplitudes.

The above mentioned findings thereby propose that nonaxisymmetric SGI will have a better chance to exist. If that is the case, does it mean the observed dust rings in PPDs are the final outcome of SGI? This important question can not be adequately addressed with a linear analysis because it is not clear if the resulting nonaxisymmetric perturbation evolves to spiral waves or collapsing fragments. But we found that growth time scale of the nonaxisymmetric SGI is very fast. It then implies that the resulting patterns are less axisymmetric and probably the observed dust rings are not caused by SGI. Further numerical simulations are needed to address this essential question.

Acknowledgements

We are grateful to referee for a constructive report that helped us to improve the manuscript. MS is also grateful to Henrik Latter for his constructive comments.